biopython之成對序列比對

锆石先生發表於2024-06-26
Python

成對序列比對

Bio.Align.PairwiseAligner

請注意，Bio.pairwise2在1.80版中被棄用。作為替代，請考慮使用Bio.Align.PairwiseAligner。

成對序列比對是透過最佳化兩個序列之間的相似性得分來將它們彼此比對的過程。Bio.Align模組包含PairwiseAligner類，用於使用Needleman-Wunsch、Smith-Waterman、Gotoh(three-state)和Waterman-Smith-Beyer全域性和區域性成對比對演算法進行全域性和區域性比對，並具有許多更改比對引數的選項。

說明：
http://biopython.org/DIST/docs/tutorial/Tutorial.html#sec127

def count_mismatches(s1, s2):
    # return sum(1 for a, b in zip(s1, s2) if a != b and (a != '-' and b != '-'))
    return sum(1 for a, b in zip(s1, s2) if a != b)

# 匯入模組
from Bio import Align

# To generate pairwise alignments, first create a PairwiseAligner object:
aligner = Align.PairwiseAligner(match_score=1.0, mismatch_score=-1.0, gap_score=-1.0, extend_gap_score=-1.0)
# 設定為全域性比對，區域性比對則是 aligner.mode = "local"
aligner.mode = "global"
seqA = 'TCAGTCG'
seqB = 'TCAGTCG'
alignments = aligner.align(seqA, seqB)
for alignment in alignments:
    print(alignment)
    print(alignment.score) #得分
    print(alignment.query)
    print(alignment.target)
    # print(alignment.aligned)
    # print(format_alignment(alignment))
    mismatches = count_mismatches(alignment.query, alignment.target) #計算mismatch的個數
    print(mismatches)

"""輸出：
TCAGTCG
|||||||
TCAGTCG

7.0
TCAGTCG
TCAGTCG
0
"""

Bio.pairwise2

請注意，Bio.pairwise2在1.80版中被棄用。作為替代，請考慮使用Bio.Align.PairwiseAligner。

說明：
https://biopython.org/docs/1.79/api/Bio.pairwise2.html
http://biopython.org/DIST/docs/tutorial/Tutorial.html#sec176

此模組中的對齊函式的名稱遵循約定

  • global + XX
  • local + XX

XX是2個字元的程式碼，表示它所採用的引數。

第一個字元表示匹配（和不匹配）的引數，第二個字元表示間隙罰分的引數。
匹配引數是：

  • x 沒有引數。 相同的字元得分為1，否則為0。
  • m 匹配分數是相同字元的分數，否則不匹配得分。
  • d 字典返回任意一對字元的分數。
  • c 回撥函式返回分數。

gap 懲罰引數是：

  • x 無gap罰款。
  • s 同樣的開放和延長兩個序列的gap罰分。
  • d 序列具有不同的開放和延長間隙罰分。
  • c 回撥函式返回間隙罰分。

示例1：

from Bio import pairwise2
from Bio.pairwise2 import format_alignment

seqA='TCAGTCG'
seqB='TCAGTCG'
alignments = pairwise2.align.globalxx(seqA, seqB)
for alignment in alignments:
    print(alignment)
    print(format_alignment(*alignment))
    print(alignment[2]) #比對得分
    # 第一個是query序列，第二個是target序列，第三列是得分，後面兩個是比對的起始和終止位置

示例2：
做全域性比對。 相同的字元給出1個點，每個不相同的字元扣除1個點，開啟間隙時扣除0.5分，並且在延長時扣除0.1分。

from Bio import pairwise2
from Bio.pairwise2 import format_alignment

seqA='TCAGTCG'
seqB='TCAGTCG'
alignments = pairwise2.align.globalms(seqA, seqB, match=1, mismatch=-1, open=-.5, extend=-.1)
for alignment in alignments:
    print(alignment)
    print(format_alignment(*alignment))
    print(alignment[2]) #比對得分
    # 第一個是query序列，第二個是target序列，第三列是得分，後面兩個是比對的起始和終止位置

示例3：
對齊函式還可以使用已包含在Biopython中的已知矩陣。

from Bio import pairwise2
from Bio.Align import substitution_matrices
matrix = substitution_matrices.load("BLOSUM62")
for a in pairwise2.align.globaldx("KEVLA", "EVL", matrix):
    print(format_alignment(*a))

示例4：
計算錯配數。

from Bio import pairwise2

def count_mismatches(s1, s2):
    # return sum(1 for a, b in zip(s1, s2) if a != b and (a != '-' and b != '-'))
    return sum(1 for a, b in zip(s1, s2) if a != b)

bc = "TCAGTCG"
seq1 = "TCAGTCG"
alignments = pairwise2.align.globalxx(bc, seq1)
for alignment in alignments:
    # 計算錯配數
    mismatches = count_mismatches(
        alignment[0], alignment[1]
    )
    print(alignment)
    print(format_alignment(*alignment))
    print(alignment[2])
    print(mismatches)

示例5：
獲得兩兩序列的相似比（這裡用的是全域性比對）

def get_similarity_pct(query, target):
    # 全域性比對，相同的殘基就給1分，不同和gap不扣分
    alignments = pairwise2.align.globalxx(query, target)
    seq_length = min(len(query), len(target))
    matches = alignments[0][2]
    percent_match = (matches / seq_length) * 100
    return percent_match

相關文章