伪随机序列生成原理_weierstrass判别法一致收敛[通俗易懂]

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 = 0
    if aa1 == aa2:  # 如果 aa1 和 aa2 相同的话，得分 +1
        score += 1
    else:
        def changeAA(aa):
            B_list = ['A','V','L','I','P','W','F','M'] # Hydrophobic
            J_list = ['D','E'] # Acidic
            O_list = ['K','R','H'] # Basic
            U_list = ['Q','N'] # Neutral
            X_list = ['G','S','T','Y','C'] # Polar
            if 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 aa
        if changeAA(aa1) == changeAA(aa2):  # 如果 aa1 和 aa2 不相同但是属于同一类的话，得分 +0.5
            score += 0.5
        else:  # 如果 aa1 和 aa2 不相同而且也不属于同一类的话，得分 +0
            score += 0
    return 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, motif
    for 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.id
        proteinseq =  str(recordR.seq)
        for recordM in SeqIO.parse(fileMotif,'fasta'):
            motifid = recordM.id
            motifseq = 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：相似性分值；

今天的文章伪随机序列生成原理_weierstrass判别法一致收敛[通俗易懂]分享到此就结束了，感谢您的阅读。