從細菌GFF文件提取CDS序列並轉換為氨基酸序列

最近在上生物信息學原理，打算記錄一些課上的作業。第一次作業：如題。

基本思路：

　　　　　　1.從GFF中讀取CDS的起始終止位置以及正負鏈信息。GFF格式見 http://blog.sina.com.cn/s/blog_8a4f556e0102yd3l.html.

　　　　　　2.利用起始/終止位置等信息從FNA文件中提取CDS序列。FNA格式見 http://boyun.sh.cn/bio/?p=1192.

　　　　　　3.利用CDS序列及密碼子表得到FAA文件並輸出。

注意：最需要注意的一點是：當GFF中CDS位於負鏈時，需要進行鹼基互補配對，即反向互補(5'到3'配3'到5')。

下面給出python代碼。python3.6

轉載請保留出處

從細菌GFF文件提取CDS序列並轉換為氨基酸序列

 

#bioinformatics homework
import re  3 class CDS2AA():  4     pa = re.compile(r'\s+')  5     Pa = re.compile(r'[TCAG]TG')                 # 細菌起始密碼子NTG
    def __init__(self,fna,gff):  7         self.fna = fna  8         self.gff = gff  9     def N2M(self,sequence): 10         hash = {'A': 'T', 'T': 'A', 'C': 'G', 'G': 'C'} 11         sequence = ''.join([hash[i] for i in sequence])     #正負鏈轉換
        return sequence[::-1] 13     def Get_CDS_index(self,line):               #獲取CDS信息
        line = self.pa.split(line) 15         CDS = (line[0], line[3], line[4], line[6], line[7])  #這里字符串分割有的文件是會出問題的，所以要看文件格式而定
        return CDS 17     def Seq2AA(self,sequence,hash): 18         AA = hash[sequence[:3]] 19         if self.Pa.match(sequence[:3]): 20             AA = 'M'                                 #起始密碼子
        for i in range(3, len(sequence) - 3, 3): 22             AA += hash[sequence[i:i + 3]] 23         return AA 24     def CDS2AA(self): 25         fn = open(self.fna, 'r') 26         gf = open(self.gff,'r') 27         r = open('AA_sequence.txt', 'w') 28         w = open('CDS.txt', 'w') 29         hash_AA = {}  # AA hash,sequence2AA
        with open('AA.txt', 'r') as f:                         #AA.txt 為密碼子表
            for line in f: 32                 line = line.strip() 33                 if line: 34                     line = self.pa.split(line) 35                     hash_AA[line[0]] = line[1]      #AA hash
        hash_CDS = {}  # CDS hash,CDS2sequence
        for line in fn: 38             line = line.strip() 39             if line.startswith('>'): 40                 A = self.pa.split(line)[0].replace('>', '') 41                 hash_CDS[A] = ''
            else: 43                 hash_CDS[A] += line 44  fn.close() 45         for line in gf: 46             line = line.strip() 47             if 'CDS' in line: 48                 i = self.Get_CDS_index(line) 49                 sequence = hash_CDS[i[0]][int(i[1]) - 1:int(i[2])] 50                 sequence = sequence[int(i[4]):]                         # i[4] 表示密碼子開始位置
                if i[3] == '-': 52                     sequence = self.N2M(sequence) 53                 #w.write(i[0] + '\n' + sequence + '\n')
                #后面是一堆正則，主要是對序列做注釋的，看文件格式而定
                
                s1 = self.pa.split(line) 57                 p1 = re.compile(r'ID=(.*?);.*?Dbxref=(.*?);.*?Name=(.*?);.*?gbkey=(.*?);.*?product=(.*?);.*?protein_id=(.*?);') 58                 p2 = re.compile(r'.*?product=(.*?);.*?protein_id=(.*?);') 59                 m1 = re.findall(p1,line) 60                 m2 = re.findall(p2,line) 61                 s = '>'+s1[0]+'_'+m1[0][0]+'\tName='+m1[0][2]+'\tdbxref='+m1[0][1]+'\tprotein='+m1[0][4]+'\tprotein_id='+m1[0][5]+'\tgbkey='+m1[0][3] 62                 w.write(s + '\n' + sequence + '\n') 63                 p = '>' + s1[0]+'\tproduct:'+m2[0][0]+'\tproduct_id:'+m2[0][1] 64                 AA = self.Seq2AA(sequence, hash_AA) 65                 r.write(p + '\n' + AA + '\n') 66  w.close() 67  r.close() 68 
if __name__=='__main__': 70     fna = 'GCA_000160075.2_ASM16007v2_genomic.fna'
    gff = 'GCA_000160075.2_ASM16007v2_genomic.gff'
    m = CDS2AA(fna,gff) 73     m.CDS2AA()

 

 

 

出現的一些問題我會慢慢完善。后面的有意思作業題目會陸續上傳。