CDS序列的全面解析

一、CDS序列的基本概念

在分子遺傳學中，CDS序列是代表編碼蛋白質的部分，是完整的編碼序列。CDS是Coding Sequence的縮寫，即編碼序列。CDS通常是從起始密碼子（AUG）到終止密碼子（UAA，UAG，UGA）之間的所有核苷酸序列。因此，CDS可以唯一地描述一個蛋白質的氨基酸序列，是進行遺傳變異和系統發育分析的關鍵基礎。

二、CDS序列在基因注釋中的重要性

在基因組學中，基因注釋是指將基因組序列轉化為有意義的功能元件。其中，CDS序列是基因注釋中的一個重要部分。通過對CDS進行分析，可以確定基因的起始位點和終止位點，從而確定氨基酸序列。此外，CDS序列還可以用於預測蛋白質結構和功能，以及比較不同物種之間的CDS序列差異。

三、CDS序列的提取方法

1. NCBI在線工具提取

NCBI提供了一種在線工具，可以將給定基因的CDS序列提取出來。用戶只需要提供該基因在NCBI數據庫中的ID或AC（Accession Number）即可。下面是一個Python程序，利用Biopython庫實現了同樣的功能：

from Bio import SeqIO
from Bio.Seq import Seq
from Bio.Alphabet import generic_dna

# 讀取基因組序列
genome_seq = SeqIO.read("genome_seq.fasta", "fasta")

# 讀取基因注釋信息
gene_annotation = open("gene_annotation.gff")

# 遍歷每一個基因，提取CDS序列
for record in SeqIO.parse(gene_annotation, "gff"):
    if record.type == "gene":
        gene_seq = genome_seq[record.start:record.end]
        cds_seq = Seq("", generic_dna)
        for feature in record.features:
            if feature.type == "CDS":
                cds_seq += feature.extract(gene_seq)

        # output CDS seq here

2. 已知基因組和基因注釋信息的情況下，直接提取

在已知基因組和基因注釋信息的情況下，可以利用基因注釋文件中的CDS起始位點和終止位點，直接從基因組序列中提取CDS序列。下面是一個Python程序：

from Bio import SeqIO
from Bio.Seq import Seq
from Bio.Alphabet import generic_dna

# 讀取基因組序列
genome_seq = SeqIO.read("genome_seq.fasta", "fasta")

# 讀取基因注釋文件
gene_annotation = open("gene_annotation.gff")

# 遍歷每一個基因，提取CDS序列
for line in gene_annotation:
    if line.startswith("#"):
        continue
    fields = line.strip().split("\t")
    if fields[2] == "CDS":
        gene_seq = genome_seq[ int(fields[3])-1 : int(fields[4]) ]
        if fields[6] == "-":
            gene_seq = gene_seq.reverse_complement()
        cds_seq = Seq("", generic_dna)
        for i in range(7, len(fields)):
            if fields[i].startswith("Parent="):
                parent = fields[i][7:]
                for line2 in gene_annotation:
                    if line2.startswith("#"):
                        continue
                    fields2 = line2.strip().split("\t")
                    if (fields2[2] == "mRNA" or fields2[2] == "transcript") and fields2[8].find(parent) != -1:
                        transcript_seq = genome_seq[ int(fields2[3])-1 : int(fields2[4]) ]
                        if fields2[6] == "-":
                            transcript_seq = transcript_seq.reverse_complement()
                        cds_seq += gene_seq[len(gene_seq) - ( int(fields[4]) - int(fields2[4]) ) : len(gene_seq) - ( int(fields[3]) - int(fields2[4]) )]
                        break
        # output CDS seq here

四、CDS序列的應用舉例

1. 基於CDS序列預測蛋白質結構與功能

CDS序列是預測蛋白質結構和功能的重要基礎。通過CDS序列，可以預測蛋白質的氨基酸序列，從而獲取蛋白質的結構和功能信息。下面是一個Python程序，利用Phyre2等工具預測蛋白質結構：

import requests

# 獲取CDS序列
cds_seq = "ATGATGAACAGGACTG..."
 
# 使用Phyre2預測蛋白質結構
response = requests.post('https://www.sbg.bio.ic.ac.uk/phyre2/webscripts/phyre2_server.cgi', data={
    'seq': cds_seq,
    'output': 'json',
    'email': 'example@example.com'
})
 
# 獲取結果
result = response.json()

2. 比較不同物種的CDS序列差異

比較不同物種的CDS序列差異可以幫助研究人員了解不同物種之間的遺傳差異。下面是一個Python程序，比較兩個物種的CDS序列：

from Bio import SeqIO
from Bio.Seq import Seq
from Bio.Alphabet import generic_dna
from Bio.Blast import NCBIWWW
from Bio.Blast import NCBIXML

# 讀取基因組序列和基因注釋文件
genome_seq1 = SeqIO.read("genome1.fasta", "fasta")
genome_seq2 = SeqIO.read("genome2.fasta", "fasta")
gene_annotation1 = open("gene_annotation1.gff")
gene_annotation2 = open("gene_annotation2.gff")

# 遍歷每一個基因，提取CDS序列，並用BLAST比對
for line in gene_annotation1:
    if line.startswith("#"):
        continue
    fields = line.strip().split("\t")
    if fields[2] == "CDS":
        gene_seq1 = genome_seq1[ int(fields[3])-1 : int(fields[4]) ]
        if fields[6] == "-":
            gene_seq1 = gene_seq1.reverse_complement()
        cds_seq1 = Seq("", generic_dna)
        for i in range(7, len(fields)):
            if fields[i].startswith("Parent="):
                parent = fields[i][7:]
                for line2 in gene_annotation1:
                    if line2.startswith("#"):
                        continue
                    fields2 = line2.strip().split("\t")
                    if (fields2[2] == "mRNA" or fields2[2] == "transcript") and fields2[8].find(parent) != -1:
                        transcript_seq1 = genome_seq1[ int(fields2[3])-1 : int(fields2[4]) ]
                        if fields2[6] == "-":
                            transcript_seq1 = transcript_seq1.reverse_complement()
                        cds_seq1 += gene_seq1[len(gene_seq1) - ( int(fields[4]) - int(fields2[4]) ) : len(gene_seq1) - ( int(fields[3]) - int(fields2[4]) )]
                        break
        result_handle = NCBIWWW.qblast("blastn", "nt", str(cds_seq1))
        blast_record = NCBIXML.read(result_handle)
        for alignment in blast_record.alignments:
            for hsp in alignment.hsps:
                if hsp.expect < 0.001:
                    print("****Alignment****")
                    print("sequence:", alignment.title)
                    print("length:", alignment.length)
                    print("e-value:", hsp.expect)
                    print(hsp.query[0:75] + "...")
                    print(hsp.match[0:75] + "...")
                    print(hsp.sbjct[0:75] + "...")
                    print("\n\n\n")

        for line in gene_annotation2:
            ...

五、總結

本文從CDS序列的基本概念開始介紹了CDS序列在基因注釋中的重要性，並詳細講解了CDS序列的提取方法。此外，本文還結合實際案例，展示了CDS序列在預測蛋白質結構與功能，比較不同物種之間的遺傳差異方面的應用，具有很好的指導意義。