[已解决]python方法计算文本中二肽在蛋白序列中出现的次数

Dawnstar

本人刚接触python 两周不到，但着急处理一个问题。
想请教给位大神：我现在又一批文本文档（具体不少于2万），文档内容是表示每个物种中的不同蛋白序列，如下所示:

>YP_009440948.1 NADH dehydrogenase subunit 6 (mitochondrion) [Absidia glauca]
MNAILLDLLAFGSVLSGILVITSRNPIISVLFLIAVFVNVACYLILLGINFIGLTYLIIYVGAIAILFLFVVMMLNIKLVELQDSAENYSNPYPLAFVLGTLFVSGLGLSNSNISKIDLPSIFDSINLFSFKSNKLETLFVSHSNWDNVFVSLDQINSVGQVLYTSHALFLVIASMILLLAMVGPIVLCLKPTKRLS
>YP_009440949.1 GIY-YIG endonuclease (mitochondrion) [Absidia glauca]
MKNNSFVQTVLTDNGWTQEESLVSIHPLSSNDTQYHSFTFKSTPVKVYHNCEINAQLILDEIRDKFGIYLWLNTVNGIMYVGSAKDLSKRLINYWTPFKSVSQCIIEMNINRNIIYK
>YP_009440950.1 NADH dehydrogenase subunit 1 (mitochondrion) [Absidia glauca]
MLLSLIEVLIVIVPLLLSVAFMTIAERKAMGSMQRRLGPNRVGYYGLLQPVADALKLFVKESVLPAHSNKALFLLAPVISLIVSLVSWGVMPFGSGLTLADLSLGMLYLLAVSSLGVYGVIFAGWAANSKYAFLGSLRSTAQMVSYEVVMGLIILTVVLLVGSLNLTEIIQSQISIWYIIPLLPLSLMFLISAIAETNRAPFDLPEAESELVAGFFTEHSSVPFVMFFLGEYASIILMSSLVSILFLGGYLVPFVSFENPTFVSFEGLSLGLKTSLILFIYIWVRASFPRLRYDQLMSFTWTGMLPLALGFIILVPCILVAFEIA
>YP_009440951.1 GIY-YIG endonuclease (mitochondrion) [Absidia glauca]
MLNNKFYYYGSSKDLGTRLKYHYYVTPKDSNKFGLFLKTVGWDYFSVTIVELCDSKDLAERETWYLQKYRPLLNTLFEVGEWPGVKFHSESTKTLISKTLTGKTHSEETKLKMSQSHQGEKNIFFNKSLPKATLDAAALVNSNLVWVYNAETKTLLKESPISSKRQTAKILGISYNSVVKYLDTDKSFKGFLMYSKEKAPV
>YP_009440952.1 ATP synthase F0 subunit 8 (mitochondrion) [Absidia glauca]
MPQLVPFYFLNQVSFAFLLLMVLLYVVSKYILPNILLVQSARMFLASK

我现在想计算每个文本文档中的两个氨基酸如（LL）在整个物种中出现的总次数（PS:每个肽键记为一次重复，如--LLLL--这个多肽序列，应该记为3个），想请问一下，我这程序应该怎样写呢？

谢谢各位大神！

最佳答案

月排行榜 / 总排行榜

DT_Nelson

2019-7-21 19:02:13

Dawnstar 发表于 2019-7-21 17:01
请问，content = re.sub('\n', '', content) #去除换行符，这个的意思是不是要把所有的蛋白序列变为一个 ...

哦，原来是这样吗，那么把那句话删去吧。

1. import re
   
2. 
   
3. filepath = "a.txt" #你的文档路径，包含所有氨基酸的次序
   
4. with open(filepath, 'r') as f:
   
5.         content = f.read() #读取内容
   
6. 
   
7. content = re.sub(r'>YP.+?]\n', '', content) #去除无关信息，只保留氨基酸序列
   
8. 
   
9. total = set(content) ^ {'\n'} #保存除'\n'（换行符）外所有的氨基酸字符
   
10. 
    
11. # 初始化字典，并完成两个不同氨基酸的查找
    
12. result = {i+j:content.count(i+j) for i in total for j in total if i!=j}
    
13. 
    
14. #连续相同氨基酸(如'LLLL')的查找
    
15. for i in total:
    
16.         result[2*i] = 0
    
17.         match = re.findall(fr'(({i})\2+)', content)
    
18.         for a, b in match:
    
19.                 result[2*i] += len(a)-1
    
20. 
    
21. lis = sorted(result, key=lambda x:result[x], reverse=True)
    
22. for i in lis:
    
23.         print(i, result[i], sep=':')

复制代码

跳转到最佳答案楼层

我是一个椭圆

1. a='''>YP_009440948.1 NADH dehydrogenase subunit 6 (mitochondrion) [Absidia glauca]
   
2. MNAILLDLLAFGSVLSGILVITSRNPIISVLFLIAVFVNVACYLILLGINFIGLTYLIIYVGAIAILFLFVVMMLNIKLVELQDSAENYSNPYPLAFVLGTLFVSGLGLSNSNISKIDLPSIFDSINLFSFKSNKLETLFVSHSNWDNVFVSLDQINSVGQVLYTSHALFLVIASMILLLAMVGPIVLCLKPTKRLS
   
3. >YP_009440949.1 GIY-YIG endonuclease (mitochondrion) [Absidia glauca]
   
4. MKNNSFVQTVLTDNGWTQEESLVSIHPLSSNDTQYHSFTFKSTPVKVYHNCEINAQLILDEIRDKFGIYLWLNTVNGIMYVGSAKDLSKRLINYWTPFKSVSQCIIEMNINRNIIYK
   
5. >YP_009440950.1 NADH dehydrogenase subunit 1 (mitochondrion) [Absidia glauca]
   
6. MLLSLIEVLIVIVPLLLSVAFMTIAERKAMGSMQRRLGPNRVGYYGLLQPVADALKLFVKESVLPAHSNKALFLLAPVISLIVSLVSWGVMPFGSGLTLADLSLGMLYLLAVSSLGVYGVIFAGWAANSKYAFLGSLRSTAQMVSYEVVMGLIILTVVLLVGSLNLTEIIQSQISIWYIIPLLPLSLMFLISAIAETNRAPFDLPEAESELVAGFFTEHSSVPFVMFFLGEYASIILMSSLVSILFLGGYLVPFVSFENPTFVSFEGLSLGLKTSLILFIYIWVRASFPRLRYDQLMSFTWTGMLPLALGFIILVPCILVAFEIA
   
7. >YP_009440951.1 GIY-YIG endonuclease (mitochondrion) [Absidia glauca]
   
8. MLNNKFYYYGSSKDLGTRLKYHYYVTPKDSNKFGLFLKTVGWDYFSVTIVELCDSKDLAERETWYLQKYRPLLNTLFEVGEWPGVKFHSESTKTLISKTLTGKTHSEETKLKMSQSHQGEKNIFFNKSLPKATLDAAALVNSNLVWVYNAETKTLLKESPISSKRQTAKILGISYNSVVKYLDTDKSFKGFLMYSKEKAPV
   
9. >YP_009440952.1 ATP synthase F0 subunit 8 (mitochondrion) [Absidia glauca]
   
10. MPQLVPFYFLNQVSFAFLLLMVLLYVVSKYILPNILLVQSARMFLASK'''
    
11. length=len(a)
    
12. i=0
    
13. count=0
    
14. sum=0
    
15. while i<length:
    
16.     if a[i]=='L':
    
17.         j=i+1
    
18.         while j<length:
    
19.             if a[j]=='L':
    
20.                 sum=sum+1
    
21.                 count=count+1
    
22.             else:
    
23.                 i = i + count
    
24.                 count=0
    
25.                 break
    
26.             j = j + 1
    
27.     i=i+1
    
28. print(sum)
    
29. 
    
30. 
    
31. 
    

复制代码

计算结果是19个肽链,用word查找也是19个,不过这个程序我自己现在都有些懵逼

DT_Nelson

emm,能不能说清楚你到底要找什么

DT_Nelson

如果只是找LL，正则表达式可以解决

DT_Nelson

1. import re
   
2. 
   
3. filepath = "a.txt" #你的文档路径，包含所有氨基酸的次序
   
4. with open(filepath, 'r') as f:
   
5.         content = f.read()
   
6. 
   
7. 
   
8. content = re.sub(r'>YP.+?]\n', '', content)
   
9. content = re.sub('\n', '', content)
   
10. 
    
11. temp = 0
    
12. total = set(content)
    
13. result = {2*i:0 for i in total}
    
14. 
    
15. for i in total:
    
16.         temp = 0
    
17.         match = re.findall(fr'((?P<t>{i})(?P=t)+)', content)
    
18.         for a, b in match:
    
19.                 temp += len(a)-1
    
20.         result[2*i] += temp
    
21. 
    
22. print(result)

复制代码

DT_Nelson

DT_Nelson 发表于 2019-7-21 14:34

用re模块，可以找出所有二肽的次数，不知道是不是你想要的？

Dawnstar

DT_Nelson 发表于 2019-7-21 13:53
emm,能不能说清楚你到底要找什么

我最终是想找所有二肽组合（也就是20*20=400种二肽在分别在相应物种中出现的次数，然后计算每一种二肽出现的次数占400种二肽出现的总次数的比例），大神不知道我说明白没有。开始只是连一种二肽的次数编程我都写不出。

Dawnstar

我是一个椭圆 发表于 2019-7-21 13:14
计算结果是19个肽链,用word查找也是19个,不过这个程序我自己现在都有些懵逼

其实我没太看懂这个代码是计算的什么（真是惭愧），还劳请大神注释一下。
其实不止想计算LL这个二肽，我是想计算20种氨基酸（20*20）400种二肽在每个物种中出现的总次数，再计算每个二肽在物种张所有蛋白中出现的次数占400种氨基酸出现的总次数的比例。
想了很久，觉得应该：
1、先将每个物种文本中的'>'行删除再计算每条蛋白中指定二肽的次数，再将每个蛋白中指定二肽次数相加，
2、’指定二肽‘ 为一个400个字符串组成的列表，再用迭代的方式将400种二肽在真个物种中出现的次数相加，
3、再计算单个二肽出现次数/400种二肽出现的总次数；
4、如果可以，是否可以将400种二肽作为一个表格的横行，物种名（有大约二万个物种）作为纵列，将计算的单个二肽出现次数占总次数的比例填充在表格中，这样就可以一目了然的看到每个物种每种二肽所占比例了。

DT_Nelson

Dawnstar 发表于 2019-7-21 15:21
我最终是想找所有二肽组合（也就是20*20=400种二肽在分别在相应物种中出现的次数，然后计算每一种二肽出 ...

我明白了

1. import re
   
2. 
   
3. filepath = "a.txt" #你的文档路径，包含所有氨基酸的次序
   
4. with open(filepath, 'r') as f:
   
5.         content = f.read() #读取内容
   
6. 
   
7. content = re.sub(r'>YP.+?]\n', '', content) #去除无关信息，只保留氨基酸序列
   
8. content = re.sub('\n', '', content) #去除换行符
   
9. 
   
10. total = set(content)
    
11. 
    
12. # 初始化字典，并完成两个不同氨基酸的查找
    
13. result = {i+j:content.count(i+j) for i in total for j in total if i!=j}
    
14. 
    
15. #连续相同氨基酸(如'LLLL')的查找
    
16. for i in total:
    
17.         result[2*i] = 0
    
18.         match = re.findall(fr'(({i})\2+)', content)
    
19.         for a, b in match:
    
20.                 result[2*i] += len(a)-1
    
21. 
    
22. lis = sorted(result, key=lambda x:result[x], reverse=True)
    
23. for i in lis:
    
24.         print(i, result[i], sep=':')  

复制代码

DT_Nelson

Dawnstar 发表于 2019-7-21 15:46
其实我没太看懂这个代码是计算的什么（真是惭愧），还劳请大神注释一下。
其实不止想计算LL这个二肽，我 ...

emm,不好意思哈，我只能帮你把出现次数算出来，至于绘制表格，可能还得你自己想点别的办法

Dawnstar

DT_Nelson 发表于 2019-7-21 16:01
emm,不好意思哈，我只能帮你把出现次数算出来，至于绘制表格，可能还得你自己想点别的办法

已经非常感谢了，我再好好理解一下，膜拜

Dawnstar

DT_Nelson 发表于 2019-7-21 15:51
我明白了

请问，content = re.sub('\n', '', content) #去除换行符，这个的意思是不是要把所有的蛋白序列变为一个字符串呢？
如果是的话，可能这样计算结果不是很精确：因为把所有蛋白序列变为一个字符串以后，前一个蛋白序列最后一个氨基酸与下一个蛋白第一个氨基酸就会形成一个新的肽键，这样在计算二肽的时候就会多记一次，这样算出的结果与实际应该会有一些偏差。

DT_Nelson

Dawnstar 发表于 2019-7-21 17:01
请问，content = re.sub('\n', '', content) #去除换行符，这个的意思是不是要把所有的蛋白序列变为一个 ...

哦，原来是这样吗，那么把那句话删去吧。

1. import re
   
2. 
   
3. filepath = "a.txt" #你的文档路径，包含所有氨基酸的次序
   
4. with open(filepath, 'r') as f:
   
5.         content = f.read() #读取内容
   
6. 
   
7. content = re.sub(r'>YP.+?]\n', '', content) #去除无关信息，只保留氨基酸序列
   
8. 
   
9. total = set(content) ^ {'\n'} #保存除'\n'（换行符）外所有的氨基酸字符
   
10. 
    
11. # 初始化字典，并完成两个不同氨基酸的查找
    
12. result = {i+j:content.count(i+j) for i in total for j in total if i!=j}
    
13. 
    
14. #连续相同氨基酸(如'LLLL')的查找
    
15. for i in total:
    
16.         result[2*i] = 0
    
17.         match = re.findall(fr'(({i})\2+)', content)
    
18.         for a, b in match:
    
19.                 result[2*i] += len(a)-1
    
20. 
    
21. lis = sorted(result, key=lambda x:result[x], reverse=True)
    
22. for i in lis:
    
23.         print(i, result[i], sep=':')

复制代码

Dawnstar

DT_Nelson 发表于 2019-7-21 19:02
哦，原来是这样吗，那么把那句话删去吧。

灰常感谢大神的帮助。

Dawnstar

通过请教给位大神，终于在今日写出了能正常跑的程序。
如下:
import os

x = ''  
y = ''
aa = ['G','A','S','T','C','V','L','I','M','P','F','Y','W','D','E','N','Q','H','K','R']
dipeptide = []
for x in aa:
    for y in aa:
        z = x + y  #得到所有的二肽
        dipeptide.append(z) # 将所有二肽整合到一个列表中        
path = 'G:\\Protein probability\\viral_split\\'  #文件夹目录
files = os.listdir(path)  #遍历文件夹下的所有文件名称
text = 'text'
file1 = open (path + text, 'w')
print(file = file1,end="\t")
for each in dipeptide:
    print(each.strip(),file = file1,end="\t")  #打印二肽到文件
print(file = file1)   
for file in files:
    if '.fa'in os.path.splitext(file)[1]:  #获取所有含‘.fa’的文件
        fa_path = path + file
        content = open(fa_path , 'r')  #读文档内的内容
        seq1 = []  
        dimertime = []
        j = {}      #将二肽建一个字典
        for s in dipeptide:
            j[s] = 0
        for seq in content:                    
            if '>' in seq:
                del seq
            else:
                seq = seq.strip() #strip去掉末尾空格和换行符
                seq1.append(seq)
        for b in seq1:
            c=list(b)
            for each_char1_index in range(len(c)-1):
                dimer = c[each_char1_index]+c[each_char1_index+1]
                if dimer in dipeptide:
                    j[dimer] +=1   
        print(file,file = file1,end="\t")   #打印文件名到文件     
        for a in dipeptide:
            print(j[a],file = file1, end = "\t")
        print(file = file1)
        
file1.close()

在此本小白向给位大神致敬！