Python3 提取CDS

根據(jù)物種基因組和注釋文件，可以編寫腳本，提取特殊的序列結(jié)構(gòu)，完成個(gè)性化的分析。本文利用Python3提取擬南芥的最長編碼序列，并轉(zhuǎn)化為氨基酸序列輸出到文件中。

基因結(jié)構(gòu)

參考鏈接：WIKIVERSITY

  • 真核生物和原核生物具有不同的基因結(jié)構(gòu)，首先我們認(rèn)識(shí)一下真核生物的基因結(jié)構(gòu)。接下來我們根據(jù)基因結(jié)構(gòu)，提取分析不同結(jié)構(gòu)。


    真核生物的基因結(jié)構(gòu)
簡并堿基

個(gè)別物種的參考基因組中存在特殊的簡并堿基，即在基因組某個(gè)位置上存在多種堿基類型，一般可以用單個(gè)堿基表示。

  • 類型
R:  A/G
Y:  C/T
M:  A/C
K:  G/T
S:  G/C
W:  A/T
H:  A/T/C
B:  G/T/C
V:  G/A/C
D:  G/A/T
N:  A/T/C/G
  • 簡并堿基的原因
    基因組測序樣品為多倍體，在同一位點(diǎn)檢測到多種堿基信號(hào)；測序樣品為混合樣品，在同一位點(diǎn)檢測到多種堿基信號(hào)；測序錯(cuò)誤；
  • 如何處理
    在擬南芥的基因組中存在的簡并堿基，在翻譯的時(shí)候統(tǒng)一翻譯成氨基酸X（參考NCBI的注釋文件）。
提取基因組基因最長蛋白編碼序列腳本
#------------------------------------------------------------------
#------------------------------------------------------------------
#作者：香波地海
#時(shí)間：2019/09/12
#腳本名稱：TAIR10.1_sub_longest_cds.py
#腳本功能：根據(jù)擬南芥注釋文件的features，提取最長編碼蛋白序列
#運(yùn)行模式：python Script.py genome.fasta genome.gff
#返回結(jié)果格式：
#   >geneid1
#   ATGTTTGGGAAACCCTGCGATGCTACGCT
#   >geneid2
#   ATGCCCGTAGCTAGCGATCGTAGCTAGCTAGCT
#---------------------------------------------------------------------------------------------------

#---------------------------------------------------------------------------------------------------
#舉例說明提取過程
#---------------------------------------------------------------------------------------------------

# cds = {}
#如果不是第一個(gè)gene信息的話，則處理上個(gè)基因里面的多個(gè)編碼蛋白序列，計(jì)算長度提取最長的編碼序列；

# cds['gene1'] = {}
# cds['gene1']['pro1'] = [('chr1','+',1,100)]
# cds['gene1']['pro1'] = [('chr1','+',1,100),('chr1','+',150,200),('chr1','+',250,300),('chr1','+',350,400),('chr1','+',450,500),]

# cds['gene1']['pro2'] = [('chr1','+',1,100)]
# cds['gene1']['pro2'] = [('chr1','+',1,100),('chr1','+',120,200),('chr1','+',220,300),('chr1','+',320,400),('chr1','+',420,500),]

# cds['gene2'] = {}
# cds['gene2']['pro3'] = [('chr1','+',1,100)]
# cds['gene2']['pro3'] = [('chr1','+',1,100),('chr1','+',110,200),('chr1','+',210,300),('chr1','+',310,400),('chr1','+',410,500),]

# cds['gene2']['pro4'] = [('chr1','+',1,100)]
# cds['gene2']['pro4'] = [('chr1','+',1,100),('chr1','+',100,200),('chr1','+',200,300),('chr1','+',300,400),('chr1','+',400,500),]

#處理最后一個(gè)基因的多個(gè)編碼蛋白序列，提取最長的蛋白編碼序列

#---------------------------------------------------------------------------------------------------
#---------------------------------------------------------------------------------------------------
import sys
import re
from Bio import SeqIO


#---------------------------------------------------------------------------------------------------
#寫一個(gè)函數(shù)readgff(),讀取基因組注釋文件，返回gene、mrna、exon、cds、longest_cds的字典
#---------------------------------------------------------------------------------------------------
def readgff(gff_file):
    '''

    ''' 
    g_id = ''
    gene = {}
    mrna = {}
    exon = {}
    cds = {}
    longest_cds = {}

    with open(gff_file) as gff:
        for line in gff:
            if line.startswith('#'):
                continue
            else:
                line = line.strip()
                (chrs,ref,features,start,end,dot,strand,start_point,note)= line.split('\t')

                if features == 'gene' :
                #-----------------------------------------------------------------------------------
                #在讀取下一個(gè)基因之前，篩選上一個(gè)基因的轉(zhuǎn)錄本，構(gòu)建每個(gè)基因?qū)?yīng)的最長轉(zhuǎn)錄本字典
                #-----------------------------------------------------------------------------------
                    if g_id != '':
                        if cds.get(g_id,None):
                            temp = {}
                            for pro, pos in cds[g_id].items():
                                sum = 0
                                for i in pos:
                                    sum += int(i[-1]) - int(i[-2]) + 1
                                temp[pro] = sum
                            longest_pro = max(temp,key=temp.get)
                            longest_cds[g_id] = cds[g_id][longest_pro]
                    #-------------------------------------------------------------------------------
                    #如果g_id為空值，說明這是注釋的第一個(gè)基因，那么讀取基因ig，構(gòu)建新的基因字典；如果g_id不是空值，
                    #那么找到上一個(gè)基因的最長轉(zhuǎn)錄本之后，記錄和存儲(chǔ)下一個(gè)基因的各種信息
                    #-------------------------------------------------------------------------------
                    g_id = note.split(';')[0].replace('ID=','')
                    gene[g_id] = (chrs,strand,start,end)
                    cds[g_id] = {}
                    exon[g_id] = {}
                    mrna[g_id] = {}
    
                if features == 'mRNA':
                    r_id = note.split(';')[0].replace('ID=','')
                    mrna[g_id][r_id] = (chrs,strand,start,end)

                if features == 'exon':
                    r_id = note.split(';')[1].replace('Parent=','')
                    if exon[g_id].get(r_id,None):
                        exon[g_id][r_id].append((chrs,strand,start,end))
                    else:
                        exon[g_id][r_id] = [(chrs,strand,start,end)]

                if features == 'CDS':   
                    c_id = note.split(';')[0].replace('ID=cds-','')
                    if cds[g_id].get(c_id,None):
                        cds[g_id][c_id].append((chrs,strand,start,end))
                    else:
                        cds[g_id][c_id] = [(chrs,strand,start,end)]
                else:
                    continue
    #-----------------------------------------------------------------------------------------------
    #分析最后一個(gè)基因的最長轉(zhuǎn)錄本
    #-----------------------------------------------------------------------------------------------
    if cds.get(g_id,None):
        temp = {}
        for pro, pos in cds[g_id].items():
            sum = 0
            for i in pos:
                sum += int(i[-1]) - int(i[-2]) + 1
            temp[pro] = sum
        longest_pro = max(temp,key=temp.get)
        longest_cds[g_id] = cds[g_id][longest_pro]
    return gene,mrna,exon,cds,longest_cds

#---------------------------------------------------------------------------------------------------
#定義一個(gè)函數(shù)用于根據(jù)基因位置列表，提取并拼接得到orf的函數(shù)sequence(postion)
# longest_cds['gene2'] = [('chr1','+',1,100),('chr1','+',100,200),('chr1','+',200,300),
#                                                       ('chr1','+',300,400),('chr1','+',400,500),]
#---------------------------------------------------------------------------------------------------

def sequence(position):
    '''
    根據(jù)函數(shù)參數(shù)元組position提供的位置參數(shù)，提取基因的堿基序列，"+"正鏈的按照起始位置切片提取，"-"負(fù)鏈的轉(zhuǎn)換成反向互補(bǔ)序列提取，
    最后輸出dna字符串
    5'----ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAGCAT----3'
    3'----TACCGGTAACATTACCCGGCGACTTTCCCACGGGCTATCGTA----5'

    '''
    hash = {'A':'T','T':'A','C':'G','G':'C','N':'N'}
    scaffold,strand,start,end = position
    if strand == '+':
        dna = Fa[str(scaffold)][int(start)-1:int(end)].upper()
    if strand == '-':
        dna1 = Fa[str(scaffold)][int(start)-1:int(end)][::-1].upper()
        dna = ''.join([hash[i] for i in dna1])
    return dna
#---------------------------------------------------------------------------------------------------
#定義一個(gè)函數(shù)dna_to_aa()，將核苷酸序列轉(zhuǎn)變成氨基酸序列。輸入為dna序列，輸出為aa序列。
#---------------------------------------------------------------------------------------------------

def dna_to_aa(dna_seq):
    '''
    定義密碼子與氨基酸對(duì)應(yīng)的字典，其中終止密碼子用*表示

    '''
    aa_dict = {
    'TTT':'F', 'TTC':'F', 'TTA':'L', 'TTG':'L', 'CTT':'L', 'CTC':'L', 'CTA':'L', 'CTG':'L', 'ATT':'I', 'ATC':'I','ATA':'I', 
    'ATG':'M', 'GTT':'V', 'GTC':'V', 'GTA':'V', 'GTG':'V', 'TCT':'S', 'TCC':'S', 'TCA':'S', 'TCG':'S', 'CCT':'P', 'CCC':'P',
    'CCA':'P', 'CCG':'P', 'ACT':'T', 'ACC':'T', 'ACA':'T', 'ACG':'T', 'GCT':'A', 'GCC':'A', 'GCA':'A', 'GCG':'A','TAT':'Y',
    'TAC':'Y',  'CAT':'H', 'CAC':'H', 'CAA':'Q', 'CAG':'Q', 'AAT':'N', 'AAC':'N', 'AAA':'K', 'AAG':'K', 'GAT':'D', 'GAC':'D', 
    'GAA':'E', 'GAG':'E', 'TGT':'C', 'TGC':'C', 'TGG':'W', 'CGT':'R', 'CGC':'R','CGA':'R','CGG':'R', 'AGT':'S', 'AGC':'S', 
    'AGA':'R', 'AGG':'R', 'GGT':'G', 'GGC':'G', 'GGA':'G', 'GGG':'G', 'TGA':'*', 'TAA':'*', 'TAG':'*',
    }

    
    #定一個(gè)空列表，用來存儲(chǔ)翻譯完成的氨基酸序列
    aa_seq = []
    for i in range(0,len(dna_seq)-2,3):
        codons = dna_seq[i:i+3]
        if aa_dict.get(codons):
            aa_seq.append(aa_dict[codons])
        else:
            aa_seq.append('@')
    aa = "".join(aa_seq)
    return aa
#---------------------------------------------------------------------------------------------------
#---------------------------------------------------------------------------------------------------


try:
    genome = sys.argv[1]
    gff = sys.argv[2]
    gene,mrna,exon,cds,longest_cds = readgff(gff)
    Fa = {rec.id:rec.seq for rec in SeqIO.parse(genome,"fasta")}
    print('There are {} genes in the genome.'.format(len(gene)))
    print('There are {} genes transcripted mRNA.'.format(len(mrna)))
    print('There are {} genes transcripted cds.'.format(len(exon)))
    print('There are {} longest proteins produced.'.format(len(longest_cds)))
#---------------------------------------------------------------------------------------------------
#---------------------------------------------------------------------------------------------------
#輸出最長轉(zhuǎn)錄本
# longest_cds['gene2'] = [('chr1','+',1,100),('chr1','+',100,200),('chr1','+',200,300),('chr1','+',300,400),('chr1','+',400,500),]
except IndexError:
    print('please input an gff and fasta file!')
finally:
    #---------------------------------------------------------------------------------------------------
    #---------------------------------------------------------------------------------------------------
    with open("out1.fa",'w') as W:
        for name,positions in longest_cds.items():
            seq = ''
            for p in positions:
                seq += sequence(p)
            aa = dna_to_aa(seq)
            W.write(">" + name + "\n")
            W.write(aa + "\n")

    print('The programs of {} have completed!'.format(sys.argv[0]))

#------------------------------------------------------------------#------------------------------------------------------------------
#結(jié)束了!
#------------------------------------------------------------------#------------------------------------------------------------------
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容