如何利用GFF注释与基因组序列提取Pun1基因CDS？结果不符排查

问题描述 #

使用Biopython和BCBio.GFF从辣椒（Capsicum annuum）基因组文件GCF_002878395.1_UCD10Xv1.1_genomic.fna和注释文件genomic.gff中提取Pun1基因（GeneID:107859694）的CDS时，得到的序列与NCBI参考CDS（NM_001324769.1:37-1359）不匹配。

参考CDS序列 #

NM_001324769.1:37-1359 LOC107859694 [organism=Capsicum annuum] [GeneID=107859694] [region=cds]
ATGGCTTTTGCATTACCATCATCACTTGTTTCAGTTTGTAACAAATCTTTTATCAAACCTTCCTCTCTCA
CCCCCTCTACACTTAGATTTCACAAGCTATCTTTCATCGATCAATCTTTAAGTAATATGTATATCCCTTG
TGCATTTTTTTACCCTAAAGTACAACAAAGACTAGAAGACTCCAAAAATTCTGATGAGCTTTCCCATATA
GCCCACTTGCTACAAACATCTCTATCACAAACTCTAGTCTCTTACTATCCTTATGCTGGAAAGTTGAAGG
ACAATGCTACTGTTGACTGTAACGATATGGGAGCTGAGTTCTTGAGTGTTCGAATAAAATGTTCCATGTC
TGAAATTCTTGATCATCCTCATGCATCTCTTGCAGAGAGCATAGTTTTGCCCAAGGATTTGCCTTGGGCG
AATAATTGTGAAGGTGGTAATTTGCTTGTAGTTCAAGTAAGTAAGTTTGATTGTGGGGGAATAGCCATCA
GTGTATGCTTTTCGCACAAGATTGGTGATGGTTGCTCTCTGCTTAATTTCCTTAATGATTGGTCTAGCGT
TACTCGTGATCGTACGACAACAACTTTAGTTCCATCTCCTAGATTTGTAGGAGATTCAGTCTTCTCTACA
CAAAAATATGGTTCTCTCATTACGCCACAAATTTTGTCCGATCTCAACCAGTGCGTACAGAAAAGACTCA
TTTTTCCTACAGATAAGTTAGATGCACTTCGAGCTAAGGTGGCAGAAGAATCAGGAGTAAAAAATCCAAC
AAGGGCTGAAGTTGTTAGCGCTCTTCTTTTCAAATGTGCAACAAAGGCATCATCATCAATGCTACCATCA
AAGTTGGTTCACTTCTTAAACATACGTACTATGATCAAACCTCGTCTACCACGAAATGCCATTGGAAATC
TCTCGTCTATTTTCTCCATAGAAGCAACTAACATGCAGGACATGGAGTTGCCAACGTTGGTTCGTAATTT
AAGGAAGGAAGTTGAGGTGGCATACAAGAAAGACCAAGTCGAACAAAATGAACTGATCCTAGAAGTAGTA
GAATCAATGAGAGAAGGGAAACTGCCATTTGAAAATATGGATGGCTATAAGAATGTGTATACTTGCAGCA
ATCTTTGCAAATATCCATACTACACTGTAGATTTTGGATGGGGAAGACCTGAAAGGGTGTGTCTAGGAAA
TGGTCCCTCCAAGAATGCCTTCTTCTTGAAAGATTACAAAGCTGGGCAAGGCGTGGAGGCGCGGGTGATG
TTGCACAAGCAACAAATGTCTGAATTTGAACGCAATGAGGAACTCCTTGAGTTCATTGCCTAA

错误输出序列 #

TTAGGCAATGAACTCAAGGAGTTCCTCATTGCGTTCAAATTCAGACATTTGTTGCTTGTG
CAACATCACCCGCGCCTCCACGCCTTGCCCAGCTTTGTAATCTTTCAAGAAGAAGGCATT
CTTGGAGGGACCATTTCCTAGACACACCCTTTCAGGTCTTCCCCATCCAAAATCTACAGT
GTAGTATGGATATTTGCAAAGATTGCTGCAAGTATACACATTCTTATAGCCATCCATATT
TTCAAATGGCAGTTTCCCTTCTCTCATTGATTCTACTACTTCTAGGATCAGTTCATTTTG
TTCGACTTGGTCTTTCTTGTATGCCACCTCAACTTCCTTCCTTAAATTACGAACCAACGT
TGGCAACTCCATGTCCTGCATGTTAGTTGCTTCTATGGAGAAAATAGACGAGAGATTTCC
AATGGCATTTCGTGGTAGACGAGGTTTGATCATAGTACGTATGTTTAAGAAGTGAACCAA
CTTTGATGGTAGCATTGATGATGATGCCTTTGTTGCACATTTGAAAAGAAGAGCGCTAAC
AACTTCAGCCCTTGTTGGATTTTTTACTCCTGATTCTTCTGCCACTTAGCTCGAAGTGCA
TCTAACTTATCTGTAGGAAAAATGAGTCTTTTCcgtgcggta

错误原因分析 #

1. 负链处理缺失：Pun1基因位于负链（strand=-1），提取的基因组序列是负链方向的原始序列，未进行反向互补转换，导致最终序列与正链的参考CDS完全反向。
2. CDS拼接逻辑不严谨：负链的CDS片段在基因组上的位置是从大到小排列的，虽然代码做了降序排序，但未对序列进行反向互补，无法得到正确的正链CDS序列。
3. 转录本选择模糊：注释文件中该基因可能存在多个mRNA转录本，代码未指定对应参考序列的转录本（NM_001324769.1），可能误选其他转录本的CDS片段。

正确提取代码 #

from Bio import SeqIO
from BCBio import GFF
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord

# 加载基因组序列
genome = SeqIO.to_dict(SeqIO.parse(
    "ncbi_dataset/capsicum annuum/data/GCF_002878395.1_UCD10Xv1.1_genomic.fna", "fasta"))

# 解析GFF注释，限制到目标染色体
record_generator = GFF.parse(
    "ncbi_dataset/capsicum annuum/data/genomic.gff", base_dict=genome, limit_info=dict(gff_id=["NC_061112.1"]))

cds_pieces = []
target_mrna_id = "NM_001324769.1"  # 指定参考转录本ID

for record in record_generator:
    for feature in record.features:
        if feature.type == "gene" and "GeneID:107859694" in feature.qualifiers.get("Dbxref", [""])[0]:
            print("找到Pun1基因")
            # 遍历基因下的所有mRNA，匹配目标转录本
            for mrna_feature in feature.sub_features:
                if mrna_feature.type == "mRNA" and target_mrna_id in mrna_feature.qualifiers.get("Dbxref", [""])[0]:
                    print(f"匹配到目标转录本: {target_mrna_id}")
                    strand = mrna_feature.location.strand
                    # 收集该转录本下的所有CDS片段
                    for cds_feature in mrna_feature.sub_features:
                        if cds_feature.type == "CDS":
                            phase = int(cds_feature.qualifiers.get("phase", [0])[0])
                            # 提取CDS序列
                            cds_seq = cds_feature.extract(record).seq
                            # 处理相位：如果相位不为0，截断前端对应长度的碱基
                            if phase != 0:
                                cds_seq = cds_seq[phase:]
                            # 存储CDS的起始位置和序列，用于排序
                            cds_pieces.append((int(cds_feature.location.start), cds_seq))
                    break
            break

# 根据链方向排序CDS片段
if strand == 1:
    # 正链：按起始位置升序排列
    cds_pieces.sort(key=lambda x: x[0])
else:
    # 负链：按起始位置降序排列
    cds_pieces.sort(key=lambda x: x[0], reverse=True)

# 拼接CDS片段
cds_concat = Seq(''.join(str(p[1]) for p in cds_pieces))

# 负链情况下，对拼接后的序列做反向互补
if strand == -1:
    cds_concat = cds_concat.reverse_complement()

# 生成SeqRecord并保存
cds_record = SeqRecord(
    cds_concat,
    id="Pun1_CDS",
    name="Pun1",
    description=f"Pun1 CDS from Capsicum annuum (transcript: {target_mrna_id})"
)

SeqIO.write(cds_record, "Pun1_CDS_correct.fasta", "fasta")
print("正确的CDS序列已保存到 Pun1_CDS_correct.fasta")

关键优化点说明 #

• 指定目标转录本：通过target_mrna_id匹配参考序列对应的mRNA，避免提取错误转录本的CDS。
• 负链序列转换：明确判断链方向，对负链拼接后的序列执行reverse_complement()，得到正链的CDS序列。
• 严谨的相位处理：保留CDS相位的截断逻辑，确保每个CDS片段的阅读框正确。
• 精简逻辑：移除不必要的计数和冗余代码，聚焦于目标CDS的提取。

内容的提问来源于stack exchange，提问作者alkyl official