[Python] How can I get CDS using annotation.gff and gene sequence.fna?

This is the CDS for Pun1 from https://www.ncbi.nlm.nih.gov/datasets/gene/id/107859694/products/

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
​

My code:

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"))

cds_record_generator = GFF.parse(
"ncbi_dataset/capsicum annuum/data/genomic.gff", base_dict=genome, limit_info=dict(gff_id=["NC_061112.1"]))

cds_count = 0
gene_count = 0
mrna_count = 0

cds_pieces = []
strand = 0

for record in cds_record_generator:
for feature in record.features:
if feature.type == "gene":
gene_count += 1
gene_id = feature.qualifiers.get("Dbxref", [""])[0]
if "GeneID:107859694" in gene_id:
print("Pun1 gene found")
print(f"location: {feature.location}")
strand = int(feature.location.strand)
try:
with open("pun1_sequence.txt", "w") as f:
f.write(
str(feature.extract(record).seq)
)
f.close()
except Exception as e:
print(e)
print("-" * 50)

for sub in feature.sub_features:
if sub.type == "mRNA":
mrna_count += 1
for subsub in sub.sub_features:
if subsub.type == "CDS":
cds_count += 1
phase = int(
subsub.qualifiers.get('phase', [0])[0])
seq_piece = subsub.extract(record).seq
if phase:
print("phase: ", phase)
seq_piece = seq_piece[phase:]
cds_pieces.append((
int(subsub.location.start),
int(subsub.location.end),
int(subsub.location.strand),
seq_piece
))


print(f"strand: {strand}")

if strand == 1:
cds_pieces.sort(key=lambda x: x[0])
else:
cds_pieces.sort(key=lambda x: x[0], reverse=True)

cds = Seq(''.join(str(piece[3]) for piece in cds_pieces))

record = SeqRecord(
cds,
id="Pun1",
name="Pun1",
description="Pun1 gene from Capsicum annuum"
)

SeqIO.write(record, "Pun1_CDS.fasta", "fasta")

print("gene count: ", gene_count)
print("mrna count: ", mrna_count)
print("cds count: ", cds_count)


Using my code, I get the following:

TTAGGCAATGAACTCAAGGAGTTCCTCATTGCGTTCAAATTCAGACATTTGTTGCTTGTG CAACATCACCCGCGCCTCCACGCCTTGCCCAGCTTTGTAATCTTTCAAGAAGAAGGCATT CTTGGAGGGACCATTTCCTAGACACACCCTTTCAGGTCTTCCCCATCCAAAATCTACAGT GTAGTATGGATATTTGCAAAGATTGCTGCAAGTATACACATTCTTATAGCCATCCATATT TTCAAATGGCAGTTTCCCTTCTCTCATTGATTCTACTACTTCTAGGATCAGTTCATTTTG TTCGACTTGGTCTTTCTTGTATGCCACCTCAACTTCCTTCCTTAAATTACGAACCAACGT TGGCAACTCCATGTCCTGCATGTTAGTTGCTTCTATGGAGAAAATAGACGAGAGATTTCC AATGGCATTTCGTGGTAGACGAGGTTTGATCATAGTACGTATGTTTAAGAAGTGAACCAA CTTTGATGGTAGCATTGATGATGATGCCTTTGTTGCACATTTGAAAAGAAGAGCGCTAAC AACTTCAGCCCTTGTTGGATTTTTTACTCCTGATTCTTCTGCCACTTAGCTCGAAGTGCA TCTAACTTATCTGTAGGAAAAATGAGTCTTTTCcgtgcggta
​

But they do not match. Why?

I want to extract the CDS from whole chromosome sequence and annotation of Camsicum annuum for a gene Pun1. how can I do that?

Continue reading...