蛋白质序列编码

Question

我正在开发一个 python 程序来计算突变残基的数字编码和一组字符串（蛋白质序列）的位置，存储在 fasta 格式文件中，每个蛋白质序列用逗号分隔。我正在尝试找到发生突变的位置和序列。

我的fasta文件如下：

MTAQDDSYSDGKGDYNTIYLGAVFQLN,MTAQDDSYSDGRGDYNTIYLGAVFQLN,MTSQEDSYSDGKGNYNTIMPGAVFQLN,MTAQDDSYSDGRGDYNTIMPGAVFQLN,MKAQDDSYSDGRGNYNTIYLGAVFQLQ,MKSQEDSYSDGRGDYNTIYLGAVFQLN,MTAQDDSYSDGRGDYNTIYPGAVFQLN,MTAQEDSYSDGRGEYNTIYLGAVFQLQ,MTAQDDSYSDGKGDYNTIMLGAVFQLN,MTAQDDSYSDGRGEYNTIYLGAVFQLN

示例：
下图（基于另一组fasta文件）将解释这背后的算法。在该图中，第一个框表示输入文件序列的对齐。最后一个框代表输出文件。如何在 Python 中使用我的 fasta 文件执行此操作？

示例输入文件：

MTAQDD,MTAQDD,MTSQED,MTAQDD,MKAQHD


        positions  1  2  3  4  5  6                  1  2  3  4  5  6    
protein sequence1  M  T  A  Q  D  D                     T  A     D
protein sequence2  M  T  A  Q  D  D                     T  A     D    
protein sequence3  M  T  S  Q  E  D                     T  S     E    
protein sequence4  M  T  A  Q  D  D                     T  A     D    
protein sequence5  M  K  A  Q  H  D                     K  A     H

     PROTEIN SEQUENCE ALIGNMENT                   DISCARD NON-VARIABLE REGION    

        positions  2  2  3  3  5  5  5    
protein sequence1  T     A     D       
protein sequence2  T     A     D       
protein sequence3  T        S     E    
protein sequence4  T     A     D       
protein sequence5     K  A           H

变异的残基被拆分为单独的列

输出文件应该是这样的：

position+residue   2T  2K  3A  3S  5D  5E  5H    
       sequence1   1   0   1   0   1   0   0    
       sequence2   1   0   1   0   1   0   0    
       sequence3   1   0   0   1   0   1   0    
       sequence4   1   0   1   0   1   0   0    
       sequence5   0   1   1   0   0   0   1

    (RESIDUES ARE CODED 1 IF PRESENT, 0 IF ABSENT)

以下是我尝试过的两种方法：

ls= 'MTAQDDSYSDGKGDYNTIYLGAVFQLN,MTAQDDSYSDGRGDYNTIYLGAVFQLN,MTSQEDSYSDGKGNYNTIMPGAVFQLN,MTAQDDSYSDGRGDYNTIMPGAVFQLN,MKAQDDSYSDGRGNYNTIYLGAVFQLQ,MKSQEDSYSDGRGDYNTIYLGAVFQLN,MTAQDDSYSDGRGDYNTIYPGAVFQLN,MTAQEDSYSDGRGEYNTIYLGAVFQLQ,MTAQDDSYSDGKGDYNTIMLGAVFQLN,MTAQDDSYSDGRGEYNTIYLGAVFQLN'.split(',')
pos = [set(enumerate(x, 1)) for x in ls]
a=set().union(*pos)
alle = sorted(set().union(*pos))
print '\t'.join(str(x) + y for x, y in alle)
for p in pos:
    print '\t'.join('1' if key in p else '0' for key in alle)

（这里我得到了突变和非突变残基的列，但我只想要突变残基的列）

from pandas import *
data = 'MTAQDDSYSDGKGDYNTIYLGAVFQLN,MTAQDDSYSDGRGDYNTIYLGAVFQLN,MTSQEDSYSDGKGNYNTIMPGAVFQLN,MTAQDDSYSDGRGDYNTIMPGAVFQLN,MKAQDDSYSDGRGNYNTIYLGAVFQLQ,MKSQEDSYSDGRGDYNTIYLGAVFQLN,MTAQDDSYSDGRGDYNTIYPGAVFQLN,MTAQEDSYSDGRGEYNTIYLGAVFQLQ,MTAQDDSYSDGKGDYNTIMLGAVFQLN,MTAQDDSYSDGRGEYNTIYLGAVFQLN'  
df = DataFrame([list(row) for row in data.split(',')])
df = DataFrame({str(col+1)+val:(df[col]==val).apply(int) for col in df.columns for val in set(df[col])})
print df.select(lambda x: not df[x].all(), axis = 1)

（这里是给出输出，但不是按顺序输出，即先是 2K，然后是 2T，然后是 3A。）

我应该怎么做？

Answer 1

get_dummies 函数为您提供最大的帮助：

In [11]: s
Out[11]: 
0    T
1    T
2    T
3    T
4    K
Name: 1

In [12]: pd.get_dummies(s, prefix=s.name, prefix_sep='')
Out[12]: 
   1K  1T
0   0   1
1   0   1
2   0   1
3   0   1
4   1   0

以及那些具有不同值的列：

In [21]: (df.ix[0] != df).any()
Out[21]: 
0    False
1     True
2     True
3    False
4     True
5    False

把这些放在一起：

In [31]: I = df.columns[(df.ix[0] != df).any()]

In [32]: J = [pd.get_dummies(df[i], prefix=df[i].name, prefix_sep='') for i in I]

In [33]: df[[]].join(J)
Out[33]: 
   1K  1T  2A  2S  4D  4E  4H
0   0   1   1   0   1   0   0
1   0   1   1   0   1   0   0
2   0   1   0   1   0   1   0
3   0   1   1   0   1   0   0
4   1   0   1   0   0   0   1

注意：我按如下方式创建了初始 DataFrame，但是根据您的情况，这可能会更有效地完成：

df = pd.DataFrame(map(list, 'MTAQDD,MTAQDD,MTSQED,MTAQDD,MKAQHD'.split(',')))