1. 目的：

比较motif和一条长的氨基酸序列的相似性，其中motif长度比较短(4-30左右)，另一条序列长度较长(150-1500左右)。

2. 思路：

将整个motif看作一个移动单元，从长氨基酸序列的N端开始，向C端移动。在此之前规定motif和长序列对应位置上的氨基酸相似性评分规则，具体如下：
假设aa1是motif上每个位置的氨基酸，aa2是长序列上每个位置的氨基酸。
1. 当 aa1 == aa2 时，分值 +1；
2. 当 aa1 != aa2，但是aa1和aa2属于同一类氨基酸时，分值 +0.5；
3. 当 aa1 != aa2，而且二者不属于同一类氨基酸时，分值 +0；
其中氨基酸的具体分类如下：

疏水氨基酸：['A','V','L','I','P','W','F','M']
酸性氨基酸：['D','E']
碱性氨基酸：['K','R','H']
中性氨基酸：['Q','N']
极性氨基酸：['G','S','T','Y','C']

评分矩阵如下所示：

当该motif在长序列上移动时，每次都会产生一组分数，将这组分数取平均分作为该motif和长序列上对应位置的序列之间的相似性分数，通过设置阈值，来对相似的序列进行筛选。(如下图所示：)

3. Python3脚本：

该脚本中默认的分数阈值是0.75，(threshold=0.75)。

from Bio import SeqIO'''设置打分函数 get_score()'''
def get_score(aa1,aa2): # aa1 和 aa2是motif和长序列上对应位置的氨基酸score = 0if aa1 == aa2:  # 如果 aa1 和 aa2 相同的话，得分 +1score += 1else:def changeAA(aa):B_list = ['A','V','L','I','P','W','F','M'] # HydrophobicJ_list = ['D','E'] # AcidicO_list = ['K','R','H'] # BasicU_list = ['Q','N'] # NeutralX_list = ['G','S','T','Y','C'] # Polarif aa in B_list:aa = 'B'elif aa in J_list:aa = 'J'elif aa in O_list:aa = 'O'elif aa in U_list:aa = 'U'elif aa in X_list:aa = 'X'return aaif changeAA(aa1) == changeAA(aa2):  # 如果 aa1 和 aa2 不相同但是属于同一类的话，得分 +0.5score += 0.5else:  # 如果 aa1 和 aa2 不相同而且也不属于同一类的话，得分 +0score += 0return score### 2. 以 motif 作为移动单元，在 protein 序列上移动，计算每个位置单元的平均分数(Sum(Score_MovingUnit)/len(motif))，理论上每一个 protein 会产生 (len(proteinSeq) - len(motif) + 1)个分数值，从中选取最大分值(分值最大为1分)对应的一个或多个，并返回对应的起始和终止位置。
'''序列比对函数 blast()'''
def blast(motif,protein,threshold_defult=0.75): # 阈值默认 0.75，如果 motif 有2个，那么至少满足其中一个是完全一样，另一个是同一类别。allmeanScore = []position = []posScore = []posSeq = []# 默认情况下 len(motif) < protein，即以短的那条序列为 motif 序列if len(motif) > len(protein):motif, protein = protein, motiffor ir in range(len(protein)-len(motif)+1):score = []for im in range(len(motif)):score.append(get_score(protein[ir+im],motif[im]))meanScore = sum(score)/len(motif)allmeanScore.append(meanScore)for index, value in enumerate(allmeanScore):if value >= threshold_defult:position.append([index+1,index+len(motif)])posScore.append(value)posSeq.append([protein[index:index+len(motif)]])return position, posScore, posSeq'''文件的输入和输出函数 file_in_out()'''
def file_in_out(fileMotif,fileProtein,outFile,threshold=0.75):outF = open(outFile,'w')outF.write('%s\t%s\t%s\t%s\t%s\t%s\t%s\n' % ('Protein1 ID','Motif ID','Start site in Protein ID','End site in Protein ID','Matched sequence in Protein ID','Motif Sequence','Score')) for recordR in SeqIO.parse(fileProtein,'fasta'):proteinid = recordR.idproteinseq =  str(recordR.seq)for recordM in SeqIO.parse(fileMotif,'fasta'):motifid = recordM.idmotifseq = str(recordM.seq)result = blast(motifseq,proteinseq,threshold_defult=threshold)for i in range(len((result[1]))):outF.write('%s\t%s\t%s\t%s\t%s\t%s\t%.4f\n' % (proteinid,motifid,result[0][i][0],result[0][i][1],str(result[2][i])[2:-2],motifseq,result[1][i]))outF.close()### 提供自己的目标输入文件filemotif.fasta和fileprotein.fasta，输出文件result.csv
filemotif = 'motif.fasta'
fileprotein = 'protein.fasta'
outfile = 'result.csv'
file_in_out(fileMotif=filemotif,fileProtein=fileprotein,outFile=outfile)##测试时可以将这两条序列放到两个 fasta文件中，file_in_out()的作为输入文件。
#Seqmotif = 'RRAR'
#Seqprotein = 'APAAGRPATTPPAPPPEEAASPAPPASPSPPGPDGDDAASPDNSTDVRAALRLAQAAGENSRFFVCPPPSGATVVRLAPARPCPEYGLGRNYTEGIGVIYKENIAPYTFKAYIYKNVIVTTTWAGSTYAAITNQYTDRVPVGMGEITDLVDKKWRCLSKAEYLRSGRKVVAFDRDDDPWEAPLKPARLSAPGVRGWHTTDDVYTALGSAGLYRTGTSVNCIVEEVEARSVYPYDSFALSTGDIIYMSPFYGLREGAHREHTSYSPERFQQIEGYYKRDMATGRRLKEPVSRNFLRTQHVTVAWDWVPKRKNVCSLAKWREADEMLRDESRGNFRFTARSLSATFVSDSHTFALQNVPLSDCVIEEAEAAVERVYRERYNGTHVLSGSLETYLARGGFVVAFRPMLSNELAKLYLQELARSNGTLEGLFAAAAPKPGPRRARRAAPSAPGGPGAANGPAGDGDAGGRVTTVSSAEFAALQFTYDHIQDHVNTMFSRLATSWCLLQNKERALWAEAAKLNPSAAASAALDRRAAARMLGDAMAVTYCHELGEGRVFIENSMRAPGGVCYSRPPVSFAFGNESEPVEGQLGEDNELLPGRELVEPCTANHKRYFRFGADYVYYENYAYVRRVPLAELEVISTFVDLNLTVLEDREFLPLEVYTRAELADTGLLDYSEIQRRNQLHELRFYDIDRVVKTDGNMAIMRGLANFFQGLGAVGQAVGTVVLGAAGAALSTVSGIASFIANPFGALATGLLVLAGLVAAFLAYRYISRLRSNPMKALYPITTRALKDDARGATAPGEEEEEFDAAKLEQAREMIKYMSLVSAVERQEHKAKKSNKGGPLLATRLTQLALRRRAPPEYQQLPMADVGGA'

4. 结果：

部分结果如下所示：

每列含义：
Protein1 ID：长序列的序列ID；
Motif ID：motif的序列ID；
Start site in Protein1 ID：与motif相似的这段序列在长序列中的起始位置；
End site in Protein1 ID：与motif相似的这段序列在长序列中的终止位置；
Matched sequence in Protein1 ID：长序列中与motif相似的这段序列；
Motif Sequence：该motif的序列；
Score：相似性分值；

比较motif和一条长序列的相似性相关推荐

两个文件比对_Edlib：方便快速的长序列比对软件包
序列比对(sequence alignment)又称序列联配, 为确定两个或多个序列之间的相似性(similarity) 或同源性(homology) ,将序列按照一定规律进行排列的操作.序列相似性和 ...
evaluate函数使用无效_使用Keras和Pytorch处理RNN变长序列输入的方法总结
最近在使用Keras和Pytorch处理时间序列数据,在变长数据的输入处理上踩了很多坑.一般的通用做法都需要先将一个batch中的所有序列padding到同一长度,然后需要在网络训练时屏蔽掉paddi ...
【重磅综述】长序列数据分析相关资源哪里找？一文读懂长序列测序数据分析的机遇与挑战！...
简介标题:长序列测序数据分析的机遇与挑战杂志:GenomeBiology 影响因子:10.806 发表时间:2020年05月08日 ...
【使用DTW将不等长序列对齐】
使用DTW将不等长序列对齐 DTW(dynamic time warping)动态时间规整,可以用来计算两个不等长时间序列的"距离",进而能比较两个时间序列的相似度.并且在这个过程 ...
pytorch中如何处理RNN输入变长序列padding
一.为什么RNN需要处理变长输入假设我们有情感分析的例子,对每句话进行一个感情级别的分类,主体流程大概是下图所示: 思路比较简单,但是当我们进行batch个训练数据一起计算的时候,我们会遇到多个训练 ...
「谷歌大脑」提出通过对长序列进行摘要提取，AI可自动生成「维基百科」
原文来源:arXiv 作者:Peter J. Liu.Mohammad Saleh.Etienne Pot.Ben Goodrich.Ryan Sepassi.Łukasz Kaiser.Noam S ...
【 MATLAB 】使用 MATLAB 作图讨论有限长序列的 N 点 DFT（强烈推荐）（含MATLAB脚本）
这篇博文本来是和上篇博文一起写的:[ MATLAB ]离散傅里叶级数(DFS)与DFT.DTFT及 z变换之间的关系但是这篇博文我最初设计的是使用MATLAB脚本和图像来讨论的,而上篇博文全是公式, ...
补零对有限长序列频谱及DFT的影响
2019独角兽企业重金招聘Python工程师标准>>> 补零对有限长序列频谱及DFT的影响 https://wenku.baidu.com/view/b304a4ce5fbfc77d ...

比较motif和一条长序列的相似性

1. 目的：

2. 思路：

3. Python3脚本：

4. 结果：

比较motif和一条长序列的相似性相关推荐

最新文章

热门文章