如何用BioPython通过PSI BLAST生成标准化PSSM矩阵？

我有一个包含8000条长序列的.fasta文件，想知道能不能用Python的BioPython包结合PSI-BLAST生成PSSM矩阵？目前我用的代码如下：

for fasta in files:
    alignment = AlignIO.read(fasta, "fasta")
    summary_align = AlignInfo.SummaryInfo(alignment)
    consensus = summary_align.dumb_consensus()
    my_pssm = summary_align.pos_specific_score_matrix(consensus, chars_to_ignore = ['N', '-'])
    file_pssm = fasta+"pssm"
    with open(file_pssm) as f:
         f.write(my_pssm)

但这段代码生成的矩阵只有0和1，我需要的是经过标准化的实际PSSM评分值，有没有更好的实现方法？

你当前代码用的AlignInfo.SummaryInfo.pos_specific_score_matrix方法，本质是生成氨基酸存在/缺失的二进制标记矩阵，并非基于同源序列频率、背景概率计算的标准化PSSM评分。要生成PSI-BLAST风格的实用PSSM，需要结合NCBI BLAST+工具调用PSI-BLAST，再处理输出结果，具体步骤如下：

1. 调用PSI-BLAST生成原始PSSM #

BioPython本身不直接运行PSI-BLAST，但可以通过subprocess调用本地安装的NCBI BLAST+工具，批量生成包含PSSM的输出文件。

示例代码：

import subprocess
import os

def run_psiblast(fasta_file, db_path, output_pssm, iterations=3):
    # 构建PSI-BLAST命令，可根据需求调整参数
    cmd = [
        "psiblast",
        "-query", fasta_file,
        "-db", db_path,
        "-num_iterations", str(iterations),
        "-out_pssm", output_pssm,
        "-evalue", "0.001",
        "-num_threads", "4"  # 多线程加速处理
    ]
    # 执行命令
    subprocess.run(cmd, check=True)

# 批量处理序列文件（建议单序列/小批次处理，避免内存过载）
for fasta in files:
    base_name = os.path.splitext(fasta)[0]
    output_pssm = f"{base_name}.pssm"
    run_psiblast(fasta, "你的本地BLAST数据库路径", output_pssm)

2. 解析并标准化PSSM #

PSI-BLAST生成的原始PSSM包含计数数据，需要转换为基于背景频率的log-odds评分（即标准化PSSM）。

示例解析与标准化代码：

import math

# 标准氨基酸背景频率（可根据研究需求调整）
background_freq = {
    'A': 0.078, 'R': 0.051, 'N': 0.041, 'D': 0.052,
    'C': 0.024, 'Q': 0.034, 'E': 0.062, 'G': 0.072,
    'H': 0.026, 'I': 0.059, 'L': 0.091, 'K': 0.056,
    'M': 0.025, 'F': 0.041, 'P': 0.046, 'S': 0.068,
    'T': 0.059, 'W': 0.013, 'Y': 0.033, 'V': 0.073
}

def parse_and_standardize_pssm(pssm_file):
    standardized_pssm = []
    with open(pssm_file, 'r') as f:
        lines = f.readlines()
        # 跳过注释行，读取有效评分数据
        for line in lines[3:-6]:
            parts = line.strip().split()
            if len(parts) < 22:
                continue
            # 提取原始氨基酸计数
            counts = list(map(int, parts[2:22]))
            total = sum(counts)
            # 计算log-odds评分（加1平滑避免除零）
            log_odds_row = []
            for idx, count in enumerate(counts):
                aa = list(background_freq.keys())[idx]
                observed_freq = (count + 1) / (total + 20)
                log_odds = math.log2(observed_freq / background_freq[aa])
                log_odds_row.append(round(log_odds, 2))
            standardized_pssm.append(log_odds_row)
    return standardized_pssm

# 使用示例
pssm = parse_and_standardize_pssm("example.pssm")
# 保存标准化后的PSSM
with open("standardized_pssm.txt", 'w') as f:
    for row in pssm:
        f.write('\t'.join(map(str, row)) + '\n')

3. 关于原始代码的问题说明 #

你用的pos_specific_score_matrix方法仅做了"该位置是否出现某氨基酸"的二值标记，没有结合同源序列的进化信息、背景频率计算，因此生成的矩阵只有0和1，不具备PSI-BLAST PSSM的生物学意义。真正的PSSM需要通过多轮PSI-BLAST迭代搜索同源序列，统计频率后转换为log-odds评分，才能体现序列位置的保守性。

内容的提问来源于stack exchange，提问作者Noorulain Islam