基于多序列比对结果的一致性序列鉴定

from Bio import AlignIO  # 导入Biopython中的AlignIO模块用于读取比对文件
import sys  # 导入sys模块用于处理命令行参数

def find_consensus_regions(alignment, min_length=100):
    """
    查找一致序列区域。
    
    参数:
        alignment (MultipleSeqAlignment): 多序列比对对象
        min_length (int): 最小的一致性区域长度，默认为10个碱基
        
    返回:
        list of tuples: 包含开始和结束位置的一致性区域列表
    """
    consensus_regions = []  # 存储找到的一致性区域
    current_region_start = None  # 当前一致性区域的起始位置
    current_stretch = 0  # 当前连续相同字符的数量
    
    # 遍历比对中的每一列
    for i in range(alignment.get_alignment_length()):
        column = alignment[:, i]  # 获取当前列的所有序列字符
        if all(base == column[0] for base in column):  # 检查该列中所有字符是否相同
            if current_region_start is None:
                current_region_start = i  # 如果没有开始位置，则设置当前列为起始位置
            current_stretch += 1  # 增加当前连续相同字符的数量
        else:
            if current_stretch >= min_length:
                # 如果当前连续相同字符数量达到或超过最小长度，则保存该区域
                consensus_regions.append((current_region_start, i - 1))
            current_region_start = None  # 重置起始位置
            current_stretch = 0  # 重置连续相同字符数量
    
    # 检查最后一段连续相同字符是否满足条件
    if current_stretch >= min_length:
        consensus_regions.append((current_region_start, i))

    return consensus_regions

def main():
    """
    主函数，用于处理命令行参数并调用查找一致序列区域的函数。
    """
    if len(sys.argv) != 2:
        print("用法: python script.py <alignment_file>")
        sys.exit(1)

    alignment_file = sys.argv[1]  # 获取命令行提供的比对文件路径
    alignment = AlignIO.read(alignment_file, "fasta")  # 读取FASTA格式的多序列比对文件
    
    consensus_regions = find_consensus_regions(alignment)  # 查找一致序列区域
    
    with open('consensus_sequences.fasta', 'w') as output_file:
        for start, end in consensus_regions:
            # 提取一致性区域的序列，并写入输出文件
            sequence = alignment[:, start:end+1]._records[0].seq
            output_file.write(f">Consensus_{start}_{end}\n{sequence}\n")

if __name__ == "__main__":
    main()  # 如果脚本作为主程序运行，则调用main函数