20230425_BLAST软件升级：新增可精确定位错配碱基位置功能

BLAST:

BLAST（Basic Local Alignment Search Tool）是一种常用的序列比对工具，用于比较两个或多个生物序列之间的相似性。BLAST算法是基于比对的方式，通过寻找两个序列之间的局部比对来确定它们的相似性，因此BLAST主要适用于寻找同源序列（homologous sequences）。

BLAST算法的流程如下：

• 预处理：生成一个参考序列数据库，并将输入序列进行格式转换以适应比对算法的需要。
• 查询：将查询序列与参考序列数据库中的序列进行比对。
• 局部比对：BLAST使用一种称为Seed-and-Extend的方法来进行局部比对，即先找到两个序列之间的“种子”序列（seed sequence），然后根据这个种子序列向两端扩展以获得更长的比对序列。
• 扩展比对：BLAST使用基于动态规划的Smith-Waterman算法来扩展种子序列的比对结果，得到最终的比对结果。
• 结果输出：BLAST将比对结果按照相似性得分（score）从高到低排序，并提供详细的比对报告和统计信息。

BLAST算法有多个不同的实现，包括BLASTN（用于比对核酸序列）、BLASTP（用于比对蛋白质序列）等，可以根据需要选择适当的实现。BLAST算法还可以通过调整参数来提高比对的准确性和速度。

BLAST软件的使用：

BLAST的使用需要先建立一个序列数据库，然后将查询序列与参考序列数据库中的序列进行比对。

#BLAST建库：BLAST建库的过程可以通过将多个序列文件合并，然后使用makeblastdb命令将其转换为BLAST所需的格式。makeblastdb命令需要指定输入文件名、输出数据库名以及数据库类型。例如，要建立一个核酸序列数据库，可以执行以下命令：

makeblastdb -in sequences.fasta -dbtype nucl -out sequences_db

#其中，-in指定输入文件名，-dbtype指定数据库类型（nucl表示核酸序列，prot表示蛋白质序列），-out指定输出数据库名。

#BLAST比对：BLAST比对可以使用blastn、blastp、tblastx等命令进行。其中，blastn用于核酸序列比对，blastp用于蛋白质序列比对，tblastx用于DNA和蛋白质序列比对。这些命令需要指定查询序列、参考数据库以及输出文件名。例如，要比对一个核酸序列，可以执行以下命令：

blastn -query query.fasta -db sequences_db -out output.txt

#其中，-query指定查询序列文件名，-db指定参考数据库名，-out指定输出文件名。

BLAST输出结果解读：

以下是一个假设的BLAST表格输出结果的例子：
query1  subject1  100.00  200  2  0  1  200  1  200  1e-100  400  
query2  subject1  90.00   190  1  1  1  200  1  190  1e-80   380   


假设我们使用BLAST进行了一个核酸序列比对，选择了表格输出格式（-outfmt 6或7）。表格输出格式是一种简洁的输出方式，每一行代表一次比对结果，各列之间以制表符（tab）分隔。以下是各列的具体含义：
qseqid：查询序列（query sequence）的标识符，通常是序列的名称或ID。
sseqid：匹配序列（subject sequence）的标识符，通常是序列的名称或ID。
pident：序列标识度（percentage of identical matches），表示查询序列和匹配序列在比对区域内的相同碱基对所占的百分比。
length：比对区域的长度，即比对上的碱基对数目。
mismatch：不匹配的碱基对数目。
gapopen：比对区域内的间隙（gap）开放次数。
qstart：比对区域在查询序列上的起始位置。
qend：比对区域在查询序列上的终止位置。
sstart：比对区域在匹配序列上的起始位置。
send：比对区域在匹配序列上的终止位置。
evalue：预期值（E-value），表示在随机情况下，可以找到具有相同或更高得分的比对次数的期望值。较小的E-value意味着比对结果更具有生物学意义。
bitscore：比对得分（bit score），反映了查询序列和匹配序列的相似性。较高的bit score意味着比对结果更具有生物学意义。

上面例子中，我们有两个比对结果。第一行表示query1和subject1的比对结果：标识度为100%，比对长度为200，有2个不匹配碱基，查询序列和匹配序列的比对区域分别为1-200和1-200，E-value为1e-100，bit score为400。第二行表示query1和subject2的比对结果，具体的解释方法与第一行类似。

BLAST软件功能升级：

上述BLAST输出的结果只包含了查询序列与参考序列之间的mismatch错配碱基数量。然而，这并不足以提供足够的信息来理解生物体内发生的过程，因为无法确定这些错配碱基在参考序列上的具体位置。因此，需要对BLAST进行功能升级。

(base) root@dell-server:/home/dell/111111/test/newtest# cat blastn-identify-mismatch-location.py
import os
import sys
import subprocess
from Bio import SearchIO
import csv
#python 6.py query.fasta  reference.fasta
#awk '{if ($4=="150" && $5 == "1") print }' query_reference_with_mismatch-blast.out

# 定义输入文件和数据库名称
query_file = sys.argv[1]
reference_file = sys.argv[2]
blast_db = sys.argv[2].strip().split(".fasta")[0]
blast_output_file = sys.argv[1].strip().split(".fasta")[0]  + "_" +  sys.argv[2].strip().split(".fasta")[0] + "-blast.out"
blast_output_with_mismatch_file = sys.argv[1].strip().split(".fasta")[0]  + "_"  + sys.argv[2].strip().split(".fasta")[0] + "_with_mismatch-blast.out"

# 使用 makeblastdb 建立数据库
makeblastdb_cmd = f"makeblastdb -in {reference_file} -dbtype nucl -out {blast_db}"
subprocess.run(makeblastdb_cmd, shell=True, check=True)
# 使用 blastn 进行比对
blastn_cmd = f"blastn -query {query_file} -db {blast_db} -out {blast_output_file} -outfmt '6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore qseq sseq'"
subprocess.run(blastn_cmd, shell=True, check=True)
"""
# 读取并解析 BLASTn 输出文件
with open(blast_output_file, 'r') as output:
    csv_reader = csv.reader(output, delimiter='\t')
    for row in csv_reader:
        query_id, subject_id, pident, length, mismatch, gapopen, qstart, qend, sstart, send, evalue, bitscore, query_seq, subject_seq = row
        qstart, qend, sstart, send = map(int, (qstart, qend, sstart, send))
        # 找到 mismatch 的位置
        for i, (q_base, s_base) in enumerate(zip(query_seq, subject_seq)):
            if q_base != s_base:
                mismatch_position_in_query = qstart + i
                mismatch_position_in_subject = sstart + i
                print(f"Query: {query_id}, Subject: {subject_id}, Mismatch: {q_base} -> {s_base}, Query pos: {mismatch_position_in_query}, Subject pos: {mismatch_position_in_subject}")
"""
# 读取并解析 BLASTn 输出文件
with open(blast_output_file, 'r') as output, open(blast_output_with_mismatch_file, 'w') as output_with_mismatch:
    csv_reader = csv.reader(output, delimiter='\t')
    csv_writer = csv.writer(output_with_mismatch, delimiter='\t')
    for row in csv_reader:
        query_id, subject_id, pident, length, mismatch, gapopen, qstart, qend, sstart, send, evalue, bitscore, query_seq, subject_seq = row
        qstart, qend, sstart, send = map(int, (qstart, qend, sstart, send))
        # 找到 mismatch 的位置
        mismatch_info = []
        for i, (q_base, s_base) in enumerate(zip(query_seq, subject_seq)):
            if q_base != s_base:
                mismatch_position_in_query = qstart + i
                mismatch_position_in_subject = sstart + i
                mismatch_info.append(f"{q_base}{mismatch_position_in_subject}{s_base}")
        # 将 mismatch 信息作为附加列添加到原始输出文件中的每一行
        new_row = row + [';'.join(mismatch_info)]
        csv_writer.writerow(new_row)

os.system('awk -F \'\t\' -v OFS=\'\t\' \'{print $1,$2,$3,$4,$5,$6,$7,$8,$9,$10,$11,$12,$15}\'  %s | column -t -s $\'\t\' | sort -n -k 9   >  tmp.txt  && mv tmp.txt %s'%(blast_output_with_mismatch_file,blast_output_with_mismatch_file))
os.system('rm -rf *.nsq  *.nin *.nhr  ')

升级后的BLAST软件输出的新结果：

# 执行python  blastn-identify-mismatch-location.py   query.fasta  reference.fasta命令后输出的新结果：

query1  subject1  100.00  200  2  0  1  200  1  200  1e-100  400  A236T,C269G
query2  subject1  90.00   190  1  1  1  200  1  190  1e-80   380    G566C

上述比对结果显示，query1与subject1之间存在两个mismatch碱基，分别位于第236和269个碱基处；而query2与subject1之间存在一个mismatch碱基，位于第566个碱基处。