序列分析与比对 ================ 序列分析与比对是生物信息学最核心的技术领域之一。从一条DNA序列的碱基组成分析，到两条序列之间的相似性比对，再到多条序列的进化关系推断，序列分析贯穿了生物信息学研究的各个环节。当你获得一条新的基因序列时，首先需要了解它的基本特征，如GC含量、开放阅读框、重复序列等；然后需要通过序列比对找到数据库中与之相似的已知序列，推断其功能；如果有多条同源序列，还需要进行多序列比对和进化分析，揭示它们之间的进化关系。本章将系统介绍序列特征解析、双序列比对、多序列比对、序列基序发现、分子进化分析等内容，帮助读者掌握序列分析的核心方法与工具。 序列特征解析 ------------ 序列特征解析是对单条序列进行基本性质分析的过程，这是序列分析的第一步。在进行任何深入的比对或功能分析之前，了解序列的基本特征是非常重要的，因为这些特征会影响后续分析的策略选择和结果解读。 碱基组成与GC含量 ~~~~~~~~~~~~~~~~ 碱基组成是描述核酸序列基本特征的最基本指标。在DNA序列中，四种碱基腺嘌呤（A）、胸腺嘧啶（T）、鸟嘌呤（G）和胞嘧啶（C）的比例分布因物种和基因组区域而异。GC含量是指序列中G和C碱基所占的比例，通常用百分比表示。不同物种的基因组GC含量差异很大，例如人类基因组整体GC含量约为41%，而某些细菌基因组的GC含量可以超过70%。在同一基因组内，不同区域的GC含量也有显著差异，基因密集区域往往具有较高的GC含量，形成所谓的GC富集区（isochores）。当你需要判断一段序列的来源物种、预测基因位置、或者设计PCR引物时，GC含量都是一个重要的参考指标。GC含量还会影响DNA的物理性质，高GC含量的DNA片段具有更高的熔解温度，这在引物设计和杂交实验中需要特别考虑。 AT/GC偏倚是对GC含量的进一步细化，它描述了序列中AT碱基对与GC碱基对之间的偏好性。在某些基因组区域，如复制起始点附近，常常存在明显的AT偏倚。AT/GC偏倚可以用GC偏斜度（GC skew）来量化，即（G-C）/（G+C），这个指标在预测复制起始点和链不对称性方面有重要应用。 .. code-block:: python from Bio.SeqUtils import gc_fraction from Bio.Seq import Seq seq = Seq("ATGGAGGAGCCGCAGTCAGATCCTAGCGTCGAGCCCCCTCTGAGTCAGGAAACATTTTCAGACCTATGGAAACTACTTCCTGAAAACAACGTTCTGTCCCCCTTGCCGTCCCAAGCAATGGATGATTTGATGCTGTCCCCGGACGATATTGAACAATGGTTCACTGAAGACCCAGGTCCAGATGA") gc_content = gc_fraction(seq) * 100 print(f"GC含量: {gc_content:.2f}%") g_count = seq.count('G') c_count = seq.count('C') a_count = seq.count('A') t_count = seq.count('T') gc_skew = (g_count - c_count) / (g_count + c_count) at_skew = (a_count - t_count) / (a_count + t_count) print(f"GC偏斜度: {gc_skew:.4f}") print(f"AT偏斜度: {at_skew:.4f}") 序列长度分布 ~~~~~~~~~~~~ 序列长度分布是描述序列数据集基本特征的另一个重要指标。在转录组测序数据中，不同基因的转录本长度差异很大，从几百个碱基到几万个碱基不等。了解序列长度分布有助于选择合适的分析方法和参数，例如在RNA-seq分析中，较长的转录本可能产生更多的读段，需要在差异表达分析中进行长度校正。在基因组组装中，contig和scaffold的长度分布（特别是N50值）是评估组装质量的关键指标。当你评估一个基因组组装结果的质量时，N50值是一个重要的参考，它表示将所有contig按长度从大到小排列后，累计长度达到总长度50%时的contig长度。 .. code-block:: python def calculate_n50(lengths): sorted_lengths = sorted(lengths, reverse=True) total = sum(sorted_lengths) cumsum = 0 for length in sorted_lengths: cumsum += length if cumsum >= total / 2: return length return 0 contig_lengths = [15000, 12000, 8000, 5000, 3000, 2000, 1500, 1000, 500, 200] n50 = calculate_n50(contig_lengths) print(f"N50值: {n50}") print(f"总长度: {sum(contig_lengths)}") print(f"最长contig: {max(contig_lengths)}") print(f"最短contig: {min(contig_lengths)}") 重复序列识别 ~~~~~~~~~~~~ 重复序列是基因组中广泛存在的一类序列，根据其排列方式可以分为串联重复和散在重复两大类。串联重复是指重复单元首尾相连排列在一起，包括微卫星序列（也称简单序列重复SSR，重复单元1-6个碱基）和小卫星序列（重复单元7-100个碱基）。微卫星序列在法医学个体识别、遗传连锁分析和群体遗传学研究中有着广泛应用，因为不同个体之间微卫星的重复次数存在高度多态性。散在重复是指重复单元散布在基因组各处，主要包括转座子（transposon）和逆转录转座子（retrotransposon）。在人类基因组中，散在重复序列约占基因组总长度的45%，其中LINE（长散在核元件）和SINE（短散在核元件）是最主要的类型。重复序列的识别对于基因组注释、引物设计和序列比对都非常重要，因为重复区域可能导致比对结果出现假阳性，在进行BLAST搜索之前通常需要先屏蔽重复序列。 简单序列与低复杂度区域 ~~~~~~~~~~~~~~~~~~~~~~~~ 简单序列是指由单一碱基或少数几种碱基重复组成的序列，包括均聚物（如poly-A、poly-T）和二核苷酸重复（如ATATAT）。低复杂度区域是指碱基组成严重偏离随机分布的序列区域，这些区域的序列复杂度很低，在序列比对中容易产生虚假的高分匹配。当你进行BLAST搜索时，如果不先过滤低复杂度区域，可能会得到大量假阳性的比对结果。NCBI BLAST默认使用DUST程序（针对核酸序列）和SEG程序（针对蛋白质序列）来过滤低复杂度区域。低复杂度区域的检测通常基于序列的信息熵或K2统计量，当某个区域的复杂度低于设定阈值时，该区域会被标记为低复杂度区域并在比对中被屏蔽。 CpG岛预测 ~~~~~~~~~~ CpG岛是基因组中GC含量较高且CpG二核苷酸出现频率高于预期的区域。在哺乳动物基因组中，CpG二核苷酸是DNA甲基化的主要位点，大部分CpG位点都被甲基化，导致CpG二核苷酸在基因组中的出现频率远低于随机预期。然而，在基因启动子区域，CpG位点通常保持未甲基化状态，形成了所谓的CpG岛。CpG岛的定义通常采用Gardiner-Garden和Frommer的标准：长度至少200bp，GC含量大于50%，观测CpG与预期CpG的比值大于0.6。当你需要预测基因的启动子位置、分析基因的表观遗传调控时，CpG岛的识别是一个重要的步骤。人类基因组中约70%的基因启动子与CpG岛关联，因此CpG岛预测是基因预测和功能注释的重要辅助手段。 开放阅读框预测 ~~~~~~~~~~~~~~ 开放阅读框（ORF）是DNA序列中从起始密码子（ATG）到终止密码子（TAA、TAG或TGA）之间的连续编码区域。ORF预测是基因预测的基础步骤，在分析一条新的DNA序列时，寻找ORF是判断该序列是否可能编码蛋白质的第一步。原核生物的基因结构相对简单，ORF通常直接对应一个蛋白质编码基因。真核生物的基因含有内含子，ORF预测更加复杂，需要考虑剪接位点的影响。ORF预测需要同时考虑六种阅读框（正链三种、负链三种），因为基因可能位于DNA的任意一条链上。当你从转录组数据中获得了新的转录本序列，想要判断它是否编码蛋白质以及编码多长的蛋白质时，ORF预测是必不可少的步骤。 .. code-block:: python from Bio.Seq import Seq dna_seq = Seq("ATGGAGGAGCCGCAGTCAGATCCTAGCGTCGAGCCCCCTCTGAGTCAGGAAACATTTTCAGACCTATGGAAACTACTTCCTGAAAACAACGTTCTGTCCCCCTTGCCGTCCCAAGCAATGGATGATTTGATGCTGTCCCCGGACGATATTGAACAATGGTTCACTGAAGACCCAGGTCCAGATGATAGCTGA") stop_codons = ['TAA', 'TAG', 'TGA'] for strand, seq in [('+', dna_seq), ('-', dna_seq.reverse_complement())]: for frame in range(3): orf_start = None for pos in range(frame, len(seq) - 2, 3): codon = str(seq[pos:pos+3]) if codon == 'ATG' and orf_start is None: orf_start = pos elif codon in stop_codons and orf_start is not None: orf_len = pos + 3 - orf_start if orf_len >= 150: protein = seq[orf_start:pos+3].translate() print(f"链{strand} 阅读框{frame}: ORF位置 {orf_start}-{pos+3}, 长度{orf_len}bp, 蛋白质长度{len(protein)-1}aa") orf_start = None 信号肽预测 ~~~~~~~~~~ 信号肽是分泌蛋白和膜蛋白N端的一段短肽序列，负责引导新合成的蛋白质穿越细胞膜进入分泌途径。信号肽通常长度为15-30个氨基酸，包含三个区域：带正电荷的n区、疏水的h区和极性的c区。信号肽在蛋白质完成跨膜转运后通常会被信号肽酶切除。当你分析一个新预测的蛋白质序列时，信号肽的预测可以帮助判断该蛋白质是分泌蛋白还是胞内蛋白，这对于理解蛋白质的功能和亚细胞定位非常重要。SignalP是最常用的信号肽预测工具，它基于神经网络和隐马尔可夫模型进行预测，最新版本还支持区分信号肽和跨膜螺旋。 跨膜结构域预测 ~~~~~~~~~~~~~~ 跨膜结构域是膜蛋白中穿越脂质双分子层的疏水区域，通常由20-25个疏水氨基酸组成α螺旋。跨膜蛋白在细胞信号传导、物质运输和能量转换等过程中发挥关键作用，人类基因组中约30%的基因编码跨膜蛋白。当你分析一个未知功能的蛋白质序列时，跨膜结构域的预测可以帮助判断它是否为膜蛋白以及跨膜次数，这对于蛋白质功能分类和药物靶点筛选非常重要。TMHMM和Phobius是常用的跨膜结构域预测工具，其中Phobius还能同时区分信号肽和跨膜螺旋，避免将信号肽误判为跨膜区域。 亚细胞定位预测 ~~~~~~~~~~~~~~ 亚细胞定位预测是根据蛋白质序列特征预测其在细胞内的定位位置，如细胞核、线粒体、叶绿体、内质网、细胞膜等。蛋白质的亚细胞定位与其功能密切相关，核定位的蛋白质通常参与转录调控，线粒体定位的蛋白质通常参与能量代谢。当你通过实验或计算预测获得了一个新蛋白质序列时，亚细胞定位预测可以为功能注释提供重要线索。亚细胞定位的预测依据包括信号肽、跨膜结构域、核定位信号、线粒体导肽等序列特征。常用的预测工具有TargetP、WoLF PSORT、DeepLoc等，其中DeepLoc基于深度学习方法，预测准确率较高。 等电点与分子量计算 ~~~~~~~~~~~~~~~~~~ 等电点（pI）是蛋白质净电荷为零时的pH值，分子量是蛋白质中所有氨基酸残基的质量之和。等电点和分子量是蛋白质最基本的物理化学参数，当你需要设计蛋白质纯化方案（如等电聚焦电泳、离子交换层析）或进行蛋白质质谱鉴定时，这两个参数是必不可少的。等电点的计算基于蛋白质中各可电离基团的pKa值，包括N端氨基、C端羧基和侧链可电离基团。分子量的计算则相对简单，将所有氨基酸残基的分子量相加即可。在实际应用中，等电点和分子量常用于2D凝胶电泳中蛋白质点的识别和质谱数据的匹配。 .. code-block:: python from Bio.SeqUtils.ProtParam import ProteinAnalysis protein_seq = "MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGP" analyzed_seq = ProteinAnalysis(protein_seq) mw = analyzed_seq.molecular_weight() pi = analyzed_seq.isoelectric_point() print(f"蛋白质长度: {len(protein_seq)} 氨基酸") print(f"分子量: {mw:.2f} Da ({mw/1000:.2f} kDa)") print(f"等电点(pI): {pi:.2f}") aa_percent = analyzed_seq.get_amino_acids_percent() print(f"疏水残基比例: {aa_percent['A']+aa_percent['V']+aa_percent['L']+aa_percent['I']+aa_percent['F']+aa_percent['W']:.2%}") 疏水性图谱 ~~~~~~~~~~ 疏水性图谱是沿蛋白质序列逐位计算疏水性得分并绘制的曲线图，最常用的疏水性标度是Kyte-Doolittle标度。在Kyte-Doolittle标度中，每个氨基酸被赋予一个疏水性值，正值表示疏水，负值表示亲水。疏水性图谱在蛋白质结构预测中有重要应用，连续的疏水区域通常对应蛋白质的跨膜螺旋或核心疏水区域，而亲水区域通常位于蛋白质表面。当你需要预测蛋白质的跨膜区域、分析蛋白质的折叠模式、或者判断某段序列是否可能形成疏水核心时，疏水性图谱是一个直观且有效的工具。通常使用窗口滑动法计算平均疏水性值，窗口大小一般选择7-19个残基，跨膜预测通常使用19个残基的窗口。 .. code-block:: python from Bio.SeqUtils.ProtParam import ProteinAnalysis protein_seq = "MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGP" analyzed_seq = ProteinAnalysis(protein_seq) kyte_doolittle = { 'I': 4.5, 'V': 4.2, 'L': 3.8, 'F': 2.8, 'C': 2.5, 'M': 1.9, 'A': 1.8, 'G': -0.4, 'T': -0.7, 'S': -0.8, 'W': -0.9, 'Y': -1.3, 'P': -1.6, 'H': -3.2, 'E': -3.5, 'Q': -3.5, 'D': -3.5, 'N': -3.5, 'K': -3.9, 'R': -4.5 } window_size = 9 hydrophobicity = [] for i in range(len(protein_seq) - window_size + 1): window = protein_seq[i:i+window_size] avg_score = sum(kyte_doolittle[aa] for aa in window) / window_size hydrophobicity.append(avg_score) for i, score in enumerate(hydrophobicity): region = "疏水区" if score > 1.0 else ("亲水区" if score < -1.0 else "中性区") print(f"位置 {i+1}-{i+window_size}: {score:.2f} ({region})") 抗原表位预测 ~~~~~~~~~~~~ 抗原表位是抗原分子中被抗体或T细胞受体识别的特定区域，分为线性表位（连续氨基酸序列）和构象表位（由蛋白质三维折叠形成的非连续区域）。抗原表位预测在疫苗设计、抗体开发和免疫诊断等领域有重要应用。当你需要设计针对某个病原体蛋白的疫苗、开发诊断用抗体、或者研究蛋白质的免疫原性时，抗原表位预测可以帮助筛选潜在的免疫反应区域。线性表位的预测主要基于氨基酸的理化性质，如亲水性、柔韧性、表面可及性和抗原性指数。构象表位的预测则需要蛋白质的三维结构信息，基于结构特征如表面可及性、空间邻近性和保守性进行预测。BepiPred和Discotope分别是线性表位和构象表位的常用预测工具。 RNA二级结构特征 ~~~~~~~~~~~~~~~ RNA二级结构是RNA分子通过碱基配对形成的局部结构，包括茎环结构（stem-loop）、发夹结构（hairpin）、内部环（internal loop）、凸起（bulge）和假结（pseudoknot）等类型。茎环结构是最基本的RNA二级结构单元，由一段互补配对的茎区和顶端的环区组成。假结是一种更复杂的结构，当环区的碱基与茎区外的碱基配对时形成，假结在核糖体RNA和核酶中经常出现。RNA二级结构对RNA的功能至关重要，tRNA的三叶草结构、rRNA的复杂折叠和核酶的催化活性都依赖于特定的二级结构。当你研究非编码RNA的功能、设计siRNA或核酶时，RNA二级结构预测是必不可少的步骤。RNAfold是ViennaRNA软件包中的核心工具，基于最小自由能算法预测RNA二级结构。 DNA弯曲度与曲率 ~~~~~~~~~~~~~~~ DNA弯曲度与曲率描述了DNA分子偏离直线构象的程度。某些DNA序列具有内在的弯曲性，这种弯曲性影响DNA与蛋白质的结合和染色质的高级结构。DNA弯曲通常与连续的A-tract（连续的腺嘌呤残基）相关，当A-tract以约10.5bp的周期出现时，会产生明显的DNA弯曲。当你研究转录因子与DNA的结合、分析核小体定位偏好、或者预测染色质可及性时，DNA弯曲度是一个需要考虑的因素。DNA弯曲度的预测可以基于序列特征，特别是A-tract的位置和周期性。 序列复杂度与模糊碱基 ~~~~~~~~~~~~~~~~~~~~ 序列复杂度衡量序列的信息含量，低复杂度区域的信息含量很低，通常由少数几种碱基重复组成。序列复杂度的计算可以基于信息熵或语言学方法。在序列比对中，低复杂度区域容易产生虚假匹配，因此需要在比对前进行过滤。 模糊碱基（也称为IUPAC模糊码）用于表示碱基位置的不确定性。例如，N表示任意碱基，R表示嘌呤（A或G），Y表示嘧啶（C或T），S表示强氢键碱基（G或C），W表示弱氢键碱基（A或T）等。当你处理测序数据中的低质量位点、设计简并引物、或者表示多态性位点时，模糊碱基是常用的表示方式。在序列比对中，模糊碱基需要特殊处理，因为它们代表多种可能的碱基，不能简单地按照精确匹配来处理。 序列转换 ~~~~~~~~ 序列转换是对核酸序列进行基本操作的过程，包括互补、反向和反向互补三种操作。互补是将每种碱基替换为其互补碱基（A-T，G-C），反是将序列的顺序颠倒，反向互补则是先互补再反向。当你设计PCR引物时，需要获取模板链的反向互补链序列；当你在基因组浏览器中查看负链基因时，看到的也是反向互补序列。反向互补操作在生物信息学分析中使用极为频繁，几乎所有的序列分析工具都内置了序列转换功能。 .. code-block:: python from Bio.Seq import Seq seq = Seq("ATGGAGGAGCCGCAGTCAGATCCTAGC") print(f"原始序列: {seq}") print(f"互补序列: {seq.complement()}") print(f"反向序列: {seq[::-1]}") print(f"反向互补: {seq.reverse_complement()}") iupac_codes = { 'R': 'AG', 'Y': 'CT', 'S': 'GC', 'W': 'AT', 'K': 'GT', 'M': 'AC', 'B': 'CGT', 'D': 'AGT', 'H': 'ACT', 'V': 'ACG', 'N': 'ACGT' } ambiguous_seq = Seq("ATGGRATGC") for i, base in enumerate(ambiguous_seq): if base in iupac_codes: print(f"位置{i+1}: {base} -> {iupac_codes[base]}") 双序列比对 ---------- 双序列比对是生物信息学中最基本的分析方法之一，它通过将两条序列进行对齐排列，找出它们之间的相似性和差异性。双序列比对的结果可以帮助判断两条序列是否具有同源性、推断序列的功能、或者分析序列之间的进化距离。当你获得一条新的基因序列并想知道它的功能时，最直接的方法就是将它与数据库中已知功能的序列进行比对。 算法基础 ~~~~~~~~ 双序列比对的算法基础是动态规划方法，它通过构建一个得分矩阵来寻找两条序列之间的最优比对路径。理解这些算法的原理有助于选择合适的比对方法和参数，以及正确解读比对结果。 Needleman-Wunsch算法是全局比对算法，它对两条序列的整个长度进行比对，适用于长度相近且整体相似的序列。该算法通过填充一个二维得分矩阵，从矩阵的右下角回溯到左上角，找到全局最优比对路径。当你需要比较两个完整的蛋白质序列或基因序列时，全局比对是合适的选择。 Smith-Waterman算法是局部比对算法，它只比对两条序列中最相似的片段，适用于长度差异较大或只有局部区域相似的序列。与全局比对不同，局部比对的得分矩阵中不允许出现负值，回溯从矩阵中的最大值开始，直到遇到零值为止。当你需要在一个长基因组序列中搜索某个基因的片段，或者比较两个只有部分区域相似的蛋白质时，局部比对更为合适。 .. code-block:: python from Bio import pairwise2 from Bio.pairwise2 import format_alignment seq1 = "ATGGAGGAGCCGCAGT" seq2 = "ATGGAGCCGCAGTCAG" global_alignments = pairwise2.align.globalms(seq1, seq2, 2, -1, -0.5, -0.1) print("全局比对 (Needleman-Wunsch):") for alignment in global_alignments[:2]: print(format_alignment(*alignment)) local_alignments = pairwise2.align.localms(seq1, seq2, 2, -1, -0.5, -0.1) print("局部比对 (Smith-Waterman):") for alignment in local_alignments[:2]: print(format_alignment(*alignment)) 点矩阵法（Dot plot）是一种直观的序列比较方法，将两条序列分别作为二维矩阵的行和列，在对应碱基相同的格点上标记点。相似区域会在图上形成对角线，而插入缺失会导致对角线的中断。点矩阵法虽然计算简单，但能直观展示序列之间的相似性模式，特别适合发现序列中的重复区域、反向互补区域和重排事件。 编辑距离（Levenshtein距离）是将一条序列转换为另一条序列所需的最少编辑操作次数，编辑操作包括替换、插入和删除。编辑距离在序列相似性度量中有广泛应用，当你需要快速计算两条短序列之间的差异时，编辑距离是一个简单有效的指标。 打分矩阵 ~~~~~~~~ 打分矩阵（也称替换矩阵）是序列比对中评估碱基或氨基酸替换得分的核心工具。打分矩阵的设计基于生物学观察：某些替换在进化中比其他替换更常出现，因此应该给予更高的得分。选择合适的打分矩阵对比对结果有重要影响。 PAM（Point Accepted Mutation）系列矩阵是基于观察到的氨基酸替换频率构建的。PAM1矩阵表示每100个氨基酸残基中平均发生1次可接受的突变，PAM250矩阵对应更远的进化距离。PAM矩阵的构建过程是先从高度相似的序列对中统计替换频率，然后通过矩阵自乘外推到更大的进化距离。当你比对进化距离较近的序列时，应该选择PAM值较小的矩阵；比对进化距离较远的序列时，应该选择PAM值较大的矩阵。 BLOSUM（BLOcks SUbstitution Matrix）系列矩阵是基于蛋白质家族中保守模块的氨基酸替换频率构建的。BLOSUM62是最常用的氨基酸替换矩阵，它基于序列同一性不超过62%的蛋白质模块构建。BLOSUM90适用于比对高度相似的序列，BLOSUM45适用于比对远缘序列。与PAM矩阵相比，BLOSUM矩阵直接基于实际观察到的替换数据，不依赖外推，因此在大多数实际应用中表现更好。 DNA替代矩阵相对简单，最基本的是一致矩阵（匹配得1分，不匹配得0分或负分）。更精细的DNA替代矩阵会考虑转换（嘌呤之间或嘧啶之间的替换，即A-G或C-T）和颠换（嘌呤与嘧啶之间的替换）的区别，因为转换在进化中比颠换更常出现，应该给予较轻的罚分。 空位罚分 ~~~~~~~~ 空位罚分模型决定了对比对中插入缺失（gap）的惩罚方式。线性空位罚分对每个空位给予固定的罚分，不考虑空位的长度，这种模型过于简单，不符合生物学实际。仿射空位罚分是最常用的模型，它将空位罚分分为开罚分（gap open penalty，打开一个新空位的罚分）和延伸罚分（gap extend penalty，延伸已有空位的罚分），开罚分通常远大于延伸罚分，这反映了生物学观察：一个连续的插入缺失比多个分散的插入缺失更可能发生。凸函数空位罚分中，空位罚分随空位长度增加的速率逐渐降低，这种模型在理论上更精确，但计算复杂度更高，实际应用较少。当你设置BLAST或比对的参数时，空位罚分的选择直接影响比对结果，较小的空位罚分会产生更多的空位，较大的空位罚分则倾向于产生更连续的匹配区域。 启发式比对 ~~~~~~~~~~ 动态规划算法虽然能保证找到最优比对，但计算复杂度为O(mn)，对于长序列和大规模数据库搜索来说太慢。启发式比对算法通过牺牲一定的灵敏度来大幅提高搜索速度，是目前数据库序列搜索的主流方法。 BLAST家族是最广泛使用的启发式比对工具。blastn用于核酸序列对核酸序列的比对，blastp用于蛋白质序列对蛋白质序列的比对，blastx将核酸序列翻译成六种阅读框的蛋白质序列后与蛋白质数据库比对，tblastn将蛋白质序列与核酸数据库的六框翻译进行比对，tblastx将查询核酸序列和数据库核酸序列都进行六框翻译后比对。当你从转录组测序中获得了一条新的转录本序列，想要找到与之相似的已知蛋白质时，blastx是最合适的选择，因为它能利用蛋白质水平的保守性来检测远缘同源关系。 FASTA是另一种启发式比对工具，包括fasta（蛋白质对蛋白质）、tfasta（蛋白质对核酸翻译）、fastf和fasts等变体。FASTA的灵敏度通常略高于BLAST，但速度稍慢，在需要更高灵敏度的搜索中可以考虑使用。 词长（k-mer）与种子扩展是BLAST算法的核心思想。BLAST首先在查询序列中找到所有长度为w的词（word），然后在数据库中搜索精确匹配这些词的序列片段作为种子，再向两侧扩展种子形成更长的比对。词长的选择影响搜索的灵敏度和速度，较短的词长提高灵敏度但降低速度，较长的词长提高速度但可能遗漏远缘同源序列。 比对统计显著性用E值（Expect value）和比特分（Bit score）来衡量。E值表示在随机情况下期望得到的得分不低于当前比对得分的命中数，E值越小表示比对结果越显著。比特分是标准化的比对得分，不受打分矩阵和空位罚分设置的影响，便于不同比对之间的比较。当你解读BLAST结果时，E值是最重要的统计指标，通常E值小于0.05或0.01被认为具有统计显著性。低复杂度区域过滤（DUST用于核酸、SEG用于蛋白质）和邻接重复屏蔽是BLAST搜索前的预处理步骤，用于避免低复杂度区域和重复序列产生虚假的高分匹配。 比对工具与参数 ~~~~~~~~~~~~~~ BLAST+是NCBI提供的命令行版BLAST工具包，比网页版功能更强大、更适合批量处理。makeblastdb用于构建BLAST数据库，blastn和blastp分别用于核酸和蛋白质序列搜索。当你需要对大量序列进行批量BLAST搜索、使用自定义数据库、或者将BLAST集成到分析流程中时，命令行版BLAST是首选工具。 .. code-block:: bash # 构建BLAST数据库 makeblastdb -in custom_sequences.fasta -dbtype nucl -out custom_db # blastn搜索 blastn -query query.fasta -db custom_db -evalue 1e-5 -outfmt 6 -num_threads 4 -out blast_results.tsv # blastp搜索 blastp -query protein_query.fasta -db nr -evalue 1e-10 -matrix BLOSUM62 -outfmt "6 qseqid sseqid pident length evalue bitscore" -out blastp_results.tsv # 短序列比对（如引物搜索） blastn -query primers.fasta -db custom_db -task blastn-short -evalue 1000 -word_size 7 -out short_blast.tsv BLAST的关键参数包括-evalue（E值阈值）、-word_size（词长）、-gapopen（空位开罚分）、-gapextend（空位延伸罚分）、-matrix（替换矩阵）和-num_threads（线程数）。E值阈值决定了报告的比对结果的最小显著性水平，默认值为10，在严格的同源性分析中通常设置为1e-5或更小。词长影响搜索策略，blastn默认词长为11，对于短序列搜索可以减小词长以提高灵敏度。空位罚分的设置需要根据序列的预期差异程度进行调整，较大的空位开罚分适合比对高度相似的序列，较小的空位开罚分适合比对远缘序列。 BLAST的输出格式通过-outfmt参数控制，格式0-5为传统的文本格式，格式6为制表符分隔格式（便于程序处理），格式7为带注释的制表格式，格式10为逗号分隔格式，格式11为ASN.1格式，格式12-14为XML和JSON格式。当你需要将BLAST结果导入Excel或用脚本处理时，制表格式（格式6）是最方便的选择。 短序列比对（如引物、gRNA）需要特殊的参数设置。由于短序列的信息量有限，需要降低E值阈值、减小词长、并使用适合短序列的打分矩阵，才能获得可靠的比对结果。 多序列比对 ---------- 多序列比对（MSA）是将三条或更多序列同时进行比对的过程，它在双序列比对的基础上进一步揭示了序列之间的保守性和变异性模式。多序列比对是构建系统发育树、发现功能位点和进行蛋白质结构预测的基础。当你获得了多个物种的同源基因序列，想要分析它们之间的进化关系和功能保守性时，多序列比对是必不可少的步骤。 算法与方法 ~~~~~~~~~~ 多序列比对的计算复杂度随序列数量呈指数增长，精确的动态规划方法只适用于少量短序列，因此实际应用中采用各种近似策略。 渐进比对（Progressive alignment）是最经典的多序列比对策略，Clustal系列工具采用这种方法。渐进比对首先计算所有序列两两之间的距离，构建一棵引导树，然后按照引导树的拓扑顺序逐步将序列加入比对。渐进比对的优点是速度快，缺点是一旦早期比对出现错误，后续比对无法纠正，存在"一旦错位，永远错位"的问题。 迭代比对方法通过多轮迭代来优化比对结果，克服了渐进比对的错误累积问题。MUSCLE采用两阶段迭代策略，第一阶段快速构建初始比对，第二阶段通过树重分割和重新比对来优化。MAFFT提供了多种策略，fft-ns-i采用快速傅里叶变换加速两两比对，l-ins-i考虑局部一致性进行迭代优化，e-ins-i允许长插入缺失的迭代优化。PRANK的特点是在比对过程中考虑插入和删除的进化事件，避免将插入错误地比对为替换，这对于进化分析特别重要。 一致性比对方法通过整合多个双序列比对的结果来构建多序列比对。T-Coffee首先生成所有序列对的双序列比对和多个结构比对，然后计算一致性权重，最终基于一致性信息构建多序列比对。MSAProbs结合了隐马尔可夫模型和渐进比对方法，利用后验概率来指导比对决策，在准确率方面表现突出。 基于隐马尔可夫模型（HMM）的方法将多序列比对建模为一个概率模型，通过训练HMM的参数来优化比对。HMMER是这类方法的代表，它构建一个配置文件HMM来表示序列家族的比对特征，然后用于搜索数据库中的同源序列。SAM（Sequence Alignment and Modeling）系统也基于HMM方法，提供了完整的序列比对和建模工具。 基于图论的方法如POA（Partial Order Alignment）将有向无环图作为比对的数据结构，能够高效处理大量序列的比对。基于分治策略的方法将多序列比对问题分解为更小的子问题分别求解，然后合并结果。基于质心估计的方法如CentroidAlign通过估计比对的后验概率来选择最可靠的比对列。联合比对的稀疏性方法如UPP和SEPP利用稀疏性来处理超大规模序列集的比对问题。 代表性工具 ~~~~~~~~~~ Clustal Omega是Clustal系列的最新版本，采用mBed算法进行序列聚类和HHalign进行轮廓比对，速度比ClustalW快数百倍，适合处理大规模序列集。ClustalW和ClustalX是早期版本，ClustalW是命令行工具，ClustalX提供图形界面，虽然速度较慢但仍然广泛使用。 .. code-block:: bash # 使用Clustal Omega进行多序列比对 clustalo -i input_sequences.fasta -o aligned.fasta --outfmt fasta --force # 使用MAFFT进行多序列比对（快速策略） mafft --auto input_sequences.fasta > aligned.fasta # 使用MAFFT高精度策略 mafft --linsi input_sequences.fasta > aligned_linsi.fasta # 使用MUSCLE进行多序列比对 muscle -in input_sequences.fasta -out aligned_muscle.fasta # 使用trimAl清洗比对结果 trimal -in aligned.fasta -out trimmed.fasta -automated1 MUSCLE以速度和准确率的平衡著称，适合中等规模序列集的快速比对。MAFFT提供了多种策略选项，fft-ns-i适合快速比对大量序列，l-ins-i适合需要高准确率的比对，e-ins-i适合包含长插入缺失的序列。当你需要比对大量序列且对速度有要求时，MAFFT的fft-ns-i策略是首选；当比对准确率更重要时，l-ins-i策略更为合适。 T-Coffee利用一致性方法整合多个信息源，准确率较高但速度较慢。MSAProbs基于概率模型，在蛋白质序列比对中准确率表现优异。PRANK考虑插入删除的进化模型，特别适合进化分析中的比对。Kalign采用基于Wu-Manber算法的快速比对策略，适合大规模序列集。FAMSA是专为超大规模序列集设计的快速比对工具。PSI-Coffee是T-Coffee的扩展版本，整合了结构信息来提高比对质量。多线程和GPU加速版本如ClustalOmega的OpenMP版本和GPU-MAFFT可以显著缩短比对时间。 比对后处理 ~~~~~~~~~~ 比对后处理是多序列比对分析的重要环节，包括比对质量优化、格式转换、可视化和质量评估等步骤。 比对清洗用于去除比对中不可靠的区域，保留高质量的比对列。Gblocks通过严格的统计标准筛选保守的比对列，trimAl提供了多种自动和手动清洗策略，BMGE基于熵值和相似性评分来选择比对列。当你构建系统发育树时，使用清洗后的比对可以减少噪声对树构建的影响，提高进化推断的可靠性。 比对格式转换是日常分析中的常见需求，常用的比对格式包括FASTA（最简单通用）、Clustal（带序列编号）、Phylip（系统发育分析的标准格式）、Nexus（包含树和比对信息）和Stockholm（Pfam数据库使用的格式）。当你需要将比对结果输入不同的分析软件时，格式转换是必不可少的步骤。 比对编辑工具如Jalview、AliView和SeaView提供了图形化的比对查看和编辑功能，可以手动调整比对中的错误。比对区域选择包括保守区选择（用于设计通用引物或寻找功能位点）和可变区选择（用于物种鉴定或多样性分析）。 间隙处理策略包括完全删除（删除含间隙的列）、成对删除（只在计算两两距离时删除间隙位置）和编码为缺失（将间隙视为缺失数据）。比对质量评估的指标包括统计保守性（每列中优势残基的频率）、列熵（每列的信息熵）和Sum-of-pairs得分（所有序列对比对列贡献的总和）。 比对可视化中最常用的是序列徽标（Sequence Logo），它用字母的高度表示每个位置上各残基的频率和信息量。WebLogo是最经典的在线序列徽标生成工具，ggseqlogo是R语言中的序列徽标包，可以方便地整合到R分析流程中。一致性序列是从多序列比对中提取的共有序列，每个位置取出现频率最高的残基。 .. code-block:: python from Bio import AlignIO from Bio.Align import AlignInfo import numpy as np alignment = AlignIO.read("aligned_sequences.fasta", "fasta") summary = AlignInfo.SummaryInfo(alignment) consensus = summary.dumb_consensus(threshold=0.5) print(f"序列数量: {len(alignment)}") print(f"比对长度: {alignment.get_alignment_length()}") print(f"一致性序列: {consensus}") for i in range(alignment.get_alignment_length()): column = alignment[:, i] freq = {base: column.count(base) / len(column) for base in set(column)} info_content = AlignInfo.SummaryInfo(alignment).information_content( start=i, end=i+1, chars_to_ignore=['-', 'N'] ) 序列基序发现 ------------ 序列基序是DNA或蛋白质序列中反复出现的短序列模式，通常与特定的生物学功能相关。在DNA序列中，基序通常是转录因子结合位点；在蛋白质序列中，基序可能是功能位点或结构标记。当你分析启动子区域的序列特征、寻找转录因子结合位点、或者发现蛋白质家族的保守功能区域时，序列基序发现是关键的分析步骤。 基序表示模型 ~~~~~~~~~~~~ 共有序列是最简单的基序表示方式，它列出基序每个位置上最常出现的碱基或氨基酸。例如，TATA框的共有序列可以表示为TATAWAWR（其中W=A或T，R=A或G）。共有序列虽然直观，但丢失了各位置上碱基频率的详细信息。 位置频率矩阵（PFM）记录了基序每个位置上各碱基出现的次数，是最原始的频率数据。位置权重矩阵（PWM）是在PFM的基础上，将频率转换为权重得分，通常使用对数似然比来计算。位置特异性得分矩阵（PSSM）与PWM类似，有时特指经过背景频率校正和log转换后的矩阵。当你需要判断一条序列是否包含某个转录因子的结合位点时，PWM/PSSM可以计算该序列与基序的匹配得分，得分越高表示匹配越好。 .. code-block:: python import numpy as np pfm = np.array([ [10, 0, 0, 0], [0, 0, 10, 0], [0, 0, 10, 0], [10, 0, 0, 0], [0, 5, 5, 0], [0, 0, 0, 10], [0, 0, 0, 10], [5, 0, 5, 0], ]) background = np.array([0.25, 0.25, 0.25, 0.25]) pseudocount = 0.5 total = pfm.sum(axis=1, keepdims=True) ppm = (pfm + pseudocount) / (total + 4 * pseudocount) pwm = np.log2(ppm / background) bases = ['A', 'C', 'G', 'T'] test_seq = "TATATATA" score = sum(pwm[i, bases.index(base)] for i, base in enumerate(test_seq)) print(f"序列 {test_seq} 的PWM得分: {score:.2f}") 序列徽标是基序的图形化表示，每个位置上字母的高度与该位置的信息量成正比，字母的堆叠顺序从下到上按频率递增。信息量可以用相对熵（Kullback-Leibler散度）来衡量，它反映了基序各位置上的碱基分布与背景分布之间的差异程度。基序的简化字母表示使用IUPAC模糊码来概括每个位置的碱基偏好性，这种表示方式在文献中经常出现。 基序发现算法 ~~~~~~~~~~~~ 基序发现算法可以分为枚举法和概率方法两大类。 枚举法穷举所有可能的短序列模式，统计其在输入序列中出现的频率，判断是否显著富集。DREME是枚举法的代表工具，它专门用于发现短而精确的DNA基序。FIMO（Find Individual Motif Occurrences）不是发现新基序的工具，而是用已知的PWM在序列中扫描匹配位点的工具。枚举法的优点是结果确定性强，缺点是对于较长的基序计算量很大。 概率方法通过迭代优化来发现基序。MEME（Multiple EM for Motif Extraction）使用期望最大化（EM）算法，从随机初始基序出发，反复迭代更新基序模型和匹配位置，直到收敛。Gibbs采样是一种随机采样方法，通过随机选择基序位置并迭代更新来搜索最优基序。MEME适合发现序列中普遍存在的基序，Gibbs采样适合发现可能只存在于部分序列中的基序。 比较基因组学方法利用跨物种的序列保守性来发现功能基序。PhyloCon通过比较多个物种的同源序列来发现保守的调控元件。PhastCons基于系统发育隐马尔可夫模型来检测保守的非编码元件。当你有多个物种的同源序列数据时，比较基因组学方法可以发现单物种分析难以检测的弱信号。 基于ChIP-seq数据的基序发现是当前最常用的转录因子结合位点发现方法。HOMER的findMotifs.pl脚本可以从ChIP-seq峰序列中发现富集的基序，同时进行已知基序的匹配和从头基序发现。MEME-ChIP是MEME Suite中专门为ChIP-seq数据设计的流程，整合了MEME、DREME、CentriMo等多个工具。基于深度学习的基序发现方法如DeepBind和DanQ利用神经网络从大规模数据中学习序列到结合亲和力的映射关系，在预测精度上取得了显著提升。 基序富集分析用于判断某个基序在一组序列中是否显著富集。AME（Analysis of Motif Enrichment）是MEME Suite中的基序富集分析工具，GSEA的基序版本可以分析基序在排序基因列表中的富集趋势。 代表性工具与数据库 ~~~~~~~~~~~~~~~~~~ MEME Suite是最完整的基序分析工具包，包括MEME（从头基序发现）、DREME（短基序发现）、FIMO（基序扫描）、MAST（基序搜索）、AME（基序富集分析）和Tomtom（基序比较）等工具。HOMER除了基序发现外，还提供了ChIP-seq峰注释和基因本体分析功能。RSAT（Regulatory Sequence Analysis Tools）是一套完整的调控序列分析工具包，包括序列检索、基序发现、基序扫描等功能。Gibbs Motif Sampler是最早的基序发现工具之一，MDscan专门用于ChIP-seq数据的基序发现，XXmotif整合了位置先验信息来提高发现精度，ProSampler用于蛋白质基序发现，GLAM2支持基序的环状排列。 JASPAR是免费开放的真核转录因子结合谱数据库，提供了转录因子PWM模型和基因组范围的结合位点预测。TRANSFAC是商业化的转录因子数据库，学术版提供有限的功能。UniPROBE存储体外蛋白质结合谱数据。HOCOMOCO提供人类和小鼠转录因子的高质量结合模型。CIS-BP涵盖原核和真核生物的转录因子结合数据。当你需要将发现的基序与已知转录因子进行匹配时，Tomtom工具可以将你的基序与JASPAR等数据库中的基序进行比较，找到最相似的已知基序。FIMO和MAST用于在序列中扫描已知的基序匹配位点。 分子进化分析 ------------ 分子进化分析利用分子序列数据来研究生物的进化关系和进化过程。通过比较不同物种的同源序列，可以推断它们之间的系统发育关系、估计进化速率、检测自然选择的作用，甚至重建祖先序列。当你需要确定不同物种中同源基因的进化关系、检测基因是否受到正选择、或者估计物种分化的时间时，分子进化分析是核心方法。 进化模型 ~~~~~~~~ 进化模型描述了分子序列在进化过程中发生替换的数学规律，是系统发育树构建和进化假设检验的基础。选择合适的进化模型对于获得可靠的进化推断至关重要。 核苷酸替代模型描述了DNA序列中碱基替换的模式。JC69（Jukes-Cantor）是最简单的模型，假设所有碱基之间的替换速率相等。K80（Kimura两参数）模型区分了转换和颠换，赋予转换更高的速率。F81（Felsenstein）模型允许不同的碱基频率。HKY85模型结合了转换/颠换差异和碱基频率差异。GTR（通用时间可逆）模型是最一般的时间可逆模型，允许六种碱基对之间有不同的替换速率。在实际分析中，通常还需要考虑位点间替换速率的异质性，通过添加伽马分布（+Γ）来建模速率变异，以及添加不变位点比例（+I）来建模完全不发生替换的位点。当你使用ModelTest或jModelTest选择最优模型时，这些扩展模型都会被考虑。 氨基酸替代模型描述了蛋白质序列中氨基酸替换的模式。Poisson模型假设所有替换等概率发生，是最简单的模型。Dayhoff模型基于早期观察的氨基酸替换频率。JTT模型基于更大的数据集更新了替换频率。WAG和LG模型分别针对球蛋白和更广泛的蛋白质数据集构建。mtREV模型专门用于线粒体编码的蛋白质。当你进行蛋白质水平的系统发育分析时，LG模型通常是较好的默认选择。 密码子替代模型描述了编码序列中密码子的替换模式，能够区分同义替换和非同义替换。Goldman-Yang模型和Muse-Gaut模型是两种主要的密码子模型框架。dN/dS比率（ω）是密码子模型中的核心参数，它表示非同义替换速率与同义替换速率的比值。ω<1表示纯化选择（purifying selection），ω=1表示中性进化，ω>1表示正选择（positive selection）。当你检测基因是否受到正选择时，dN/dS分析是最常用的方法。 系统发育树构建 ~~~~~~~~~~~~~~ 系统发育树表示了分类单元之间的进化关系，树上的分支模式反映了物种或基因的分化历史。系统发育树的构建方法可以分为距离法、最大简约法、最大似然法和贝叶斯法。 距离法首先计算所有分类单元两两之间的进化距离，然后基于距离矩阵构建树。邻接法（NJ）是最常用的距离法，它逐步选择距离最近的分类单元对进行合并，直到所有分类单元都加入树中。NJ法速度快，适合快速构建初始树或处理大数据集。最小进化法（ME）寻找使树的总枝长最小的拓扑结构。非加权对组法（UPGMA）假设分子钟恒定，产生超度量树，当进化速率在不同谱系间差异较大时可能产生错误的拓扑结构。距离矩阵的计算需要选择合适的进化模型，如p-距离、Kimura距离或基于特定替换模型的距离。 .. code-block:: python from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor from Bio import AlignIO from Bio import Phylo alignment = AlignIO.read("aligned_sequences.fasta", "fasta") calculator = DistanceCalculator('identity') dm = calculator.get_distance(alignment) constructor = DistanceTreeConstructor() nj_tree = constructor.nj(dm) upgma_tree = constructor.upgma(dm) Phylo.draw_ascii(nj_tree) 最大简约法（MP）寻找需要最少进化步骤（替换事件）的树拓扑结构。Fitch算法用于计算无权简约得分，加权简约允许不同类型的替换有不同的权重。最大简约法的优点是不需要显式的进化模型，缺点是在长枝吸引（long-branch attraction）的情况下可能产生错误的拓扑，即进化速率快的两个远缘分类单元可能被错误地聚在一起。 最大似然法（ML）在给定进化模型和树拓扑的条件下，计算观测数据的似然值，选择似然值最大的树。ML法的树搜索策略包括最近邻交换（NNI）、子树修剪重接（SPR）和树二分重接（TBR），从简单到复杂搜索更广的树空间。分支支持度评估包括非参数自举（bootstrap，通常1000次重复）、SH检验和近似似然比检验（aLRT）。ML法是目前最广泛使用的建树方法，在准确率和计算效率之间取得了较好的平衡。 贝叶斯法通过马尔可夫链蒙特卡洛（MCMC）采样来估计树的后验概率分布。MrBayes是最常用的贝叶斯建树工具，BEAST可以构建时间标度的系统发育树（分子钟树），同时估计分化时间。贝叶斯法的优点是可以直接获得每个分支的后验概率作为支持度，缺点是MCMC运行时间长，且需要判断链是否收敛。 合并树（consensus tree）用于整合多棵树的信息，严格合并树只包含在所有输入树中都出现的分支，多数合并树包含在超过50%输入树中出现的分支。树可视化工具包括FigTree（轻量级桌面工具）、ggtree（R语言包，功能强大）、iTOL（在线工具，支持丰富注释）和Dendroscope（支持大树可视化）。当你需要将系统发育树用于论文发表时，ggtree和iTOL可以生成高质量的注释图。 进化分析工具 ~~~~~~~~~~~~ MEGA是最流行的图形化进化分析软件，集成了序列比对、模型选择、树构建和树编辑等功能，适合初学者使用。IQ-TREE是高效的最大似然法建树工具，内置了ModelFinder进行自动模型选择，支持超快速自举分析。RAxML是大规模最大似然法建树的标准工具，特别适合处理大数据集。FastTree使用近似最大似然法，速度比RAxML快数倍，适合超大规模数据集的快速分析。PhyML是另一个常用的ML建树工具，提供了多种树搜索选项。 .. code-block:: bash # IQ-TREE：自动模型选择+建树+自举 iqtree2 -s aligned_sequences.fasta -m MFP -bb 1000 -alrt 1000 -nt AUTO # RAxML：最大似然法建树 raxmlHPC -s aligned_sequences.fasta -n output -m GTRGAMMA -p 12345 -# 20 # FastTree：快速近似ML建树 FastTree -gtr -nt aligned_sequences.fasta > tree.nwk # MEGA命令行（模型选择） megacc -a model_selection.mao -d aligned_sequences.fasta PAUP*是经典的系统发育分析软件，支持多种建树方法。MrBayes和BEAST/BEAUti是贝叶斯建树的主要工具。PAML（Phylogenetic Analysis by Maximum Likelihood）专注于密码子水平的进化分析，其CODEML程序是检测正选择的标准工具。HyPhy提供了灵活的选择压力检测方法，包括MEME、FEL、SLAC、BUSTED等方法。PHYLIP是最早的系统发育分析软件包之一。DendroPy和ETE Toolkit是Python库，可以在脚本中程序化地处理系统发育分析。 进化假设检验 ~~~~~~~~~~~~ 进化假设检验是分子进化分析的核心应用之一，通过统计检验来判断进化假设是否得到数据支持。 分子钟检验判断不同谱系的进化速率是否恒定。似然比检验比较有分子钟约束和无分子钟约束的模型，如果两个模型的似然值差异显著，则拒绝分子钟假设。相对速率检验通过外群来比较两个内群分类单元的进化速率。 正选择检测是进化分析中最受关注的应用。dN/dS>1表明基因受到正选择，即非同义突变被自然选择保留。分支模型允许不同分支上有不同的ω值，可以检测特定谱系上的正选择。位点模型允许不同密码子位点上有不同的ω值，M7（beta分布）vs M8（beta分布+ω>1的位点类别）的似然比检验是检测位点水平正选择的标准方法。分支-位点模型同时考虑分支和位点的差异，可以检测特定分支上特定位点的正选择，BEB（Bayes Empirical Bayes）方法可以识别受到正选择的具体位点。 .. code-block:: bash # PAML CODEML：位点模型检测正选择 # codeml配置文件 (codeml.ctl) # seqfile = aligned_codons.phy # treefile = tree.nwk # outfile = codeml_out.txt # runmode = 0 # model = 0 # NSsites = 7 8 codeml codeml.ctl # HyPhy：MEME检测偶发性正选择 hyphy meme --alignment aligned_codons.fasta --tree tree.nwk # HyPhy：BUSTED检测基因水平正选择 hyphy busted --alignment aligned_codons.fasta --tree tree.nwk # HyPhy：FEL检测位点水平选择 hyphy fel --alignment aligned_codons.fasta --tree tree.nwk RELAX方法检测选择压力的放松或加强，而不是正选择本身。MEME（Mixed Effects Model of Evolution）检测个别位点上的偶发性正选择。FEL（Fixed Effects Likelihood）使用固定效应模型估计每个位点的dN和dS。SLAC（Single-Loss Ancestor Counting）基于祖先序列重建的计数方法估计选择压力。BUSTED检测基因水平的偶发性正选择。共进化分析检测不同位点之间的进化相关性，方法包括CD检测和互信息方法。重组检测识别序列中的重组事件，RDP、GARD和Phi test是常用的重组检测方法。趋同进化检测判断不同谱系是否独立进化出相似的分子特征。 进化历史重建 ~~~~~~~~~~~~ 祖先序列重建（ASR）是根据现存物种的序列和系统发育树推断祖先节点上的序列。边际重建为每个位点独立推断最可能的祖先状态，联合重建考虑位点之间的相关性推断整体最可能的祖先序列。当你需要研究蛋白质的进化历史、重建古代基因的功能、或者设计具有特定性质的蛋白质时，祖先序列重建是一个有力的工具。不确定性可视化可以展示祖先序列推断的可信度。 基因重复与丢失是基因家族进化的重要事件，Notung工具可以协调基因树与物种树之间的不一致，识别基因重复和丢失事件。水平基因转移（HGT）检测判断基因是否通过非垂直遗传的方式在不同物种之间转移。物种树与基因树冲突（ILS，Incomplete Lineage Sorting / Deep Coalescence）是另一个导致基因树与物种树不一致的原因，它源于群体遗传学效应。多序列比对与树的一致性评估是确保进化分析可靠性的重要步骤。 序列比较的特殊应用 ------------------ 序列比较技术除了基本的同源性检测外，还有许多重要的特殊应用场景，这些应用将序列比对技术扩展到了更广泛的生物学问题。 直系同源与旁系同源判定 ~~~~~~~~~~~~~~~~~~~~~~~~ 直系同源（ortholog）和旁系同源（paralog）是同源基因的两种类型。直系同源基因是由物种分化事件产生的同源基因，通常保留相似的功能；旁系同源基因是由基因重复事件产生的同源基因，功能可能发生分化。当你需要将一个物种的基因功能注释转移到另一个物种时，找到直系同源基因是最可靠的方式，因为直系同源基因更可能保留相同的功能。OrthoFinder是目前最常用的直系同源基因鉴定工具，它同时推断直系同源群和物种树。InParanoid用于鉴定种内旁系同源和种间直系同源基因。eggNOG数据库提供了预计算的直系同源群和功能注释。 基因家族聚类是将多个物种的同源基因按照进化关系进行分组。OrthoMCL基于MCL图聚类算法进行基因家族聚类，OrthoFinder在聚类的同时还推断基因树的进化关系。当你分析一个新物种的基因组时，基因家族聚类可以帮助识别基因的扩张和收缩，这些事件可能与物种的适应性进化相关。 全基因组比对与共线性分析 ~~~~~~~~~~~~~~~~~~~~~~~~~~ 全基因组比对是比较两个或多个完整基因组序列的过程，用于发现基因组结构变异、保守区域和重排事件。MUMmer是快速的全基因组比对工具，适合比较两个近缘物种的完整基因组。LASTZ用于远缘物种之间的全基因组比对。minimap2最初为长读长比对设计，也支持全基因组比对。当你需要比较两个品种的基因组差异、寻找结构变异、或者分析基因组重排时，全基因组比对是必要的步骤。 .. code-block:: bash # MUMmer全基因组比对 nucmer --maxmatch -p genome_comparison reference.fasta query.fasta delta-filter -m genome_comparison.delta > filtered.delta show-coords -rcl filtered.delta > coords.txt # minimap2全基因组比对 minimap2 -x asm5 reference.fasta query.fasta > alignment.paf # MCScanX共线性分析 MCScanX gene_collinearity 共线性分析检测基因组之间保守的基因排列顺序，是基因组进化研究的重要方法。MCScanX是最常用的共线性分析工具，可以检测共线性区域并分析基因重复事件。SynMap是CoGe平台上的在线共线性分析工具。基因组重排检测可以发现倒位、易位、融合和裂解等大规模基因组结构变异。 系统发育足迹分析与增强子识别 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 系统发育足迹分析（phylogenetic footprinting）通过比较多个物种的非编码序列来发现保守的调控元件。由于功能重要的调控元件在进化中受到选择压力的约束，它们在近缘物种之间会表现出高于背景的序列保守性。增强子和保守非编码元件（CNE）的识别是系统发育足迹分析的重要应用。phastCons和phyloP是基于多物种比对数据检测保守元件的工具，phastCons基于隐马尔可夫模型识别保守区域，phyloP对每个位点进行保守性检验。当你需要预测基因的调控元件、发现非编码区的功能区域时，系统发育足迹分析是一种有效的方法。 宏条形码与短读长比对 ~~~~~~~~~~~~~~~~~~~~ 宏条形码序列比对是微生物组研究中的关键步骤。VSEARCH和usearch是常用的宏条形码序列处理工具，支持序列去重、聚类和比对功能。VSEARCH是usearch的开源替代品，提供了与usearch兼容的功能接口。 短读长比对到参考基因组是基因组重测序分析的核心步骤。BWA（Burrows-Wheeler Aligner）是最广泛使用的短读长比对工具，支持多种比对模式。Bowtie2采用FM-index进行快速比对，适合中等长度的读段。HISAT2是专门为RNA-seq数据设计的剪接感知比对工具，能够将跨越外显子-内含子边界的读段正确比对到基因组上。这些工具虽然主要用于基因组比对而非传统的序列比对，但它们的核心技术——基于Burrows-Wheeler变换的索引和比对——与序列比对的理论基础密切相关。 .. code-block:: bash # BWA比对流程 bwa index reference.fasta bwa mem -t 8 reference.fasta reads_R1.fastq reads_R2.fastq > aligned.sam # Bowtie2比对 bowtie2-build reference.fasta ref_index bowtie2 -x ref_index -1 reads_R1.fastq -2 reads_R2.fastq -S aligned.sam # HISAT2剪接感知比对 hisat2-build reference.fasta ref_index hisat2 -x ref_index -1 reads_R1.fastq -2 reads_R2.fastq -S rna_aligned.sam 比对性能与统计 -------------- 理解比对性能和统计原理对于正确解读比对结果和选择合适的比对方法至关重要。比对结果不仅仅是序列的排列，更包含了丰富的统计信息，这些信息帮助我们判断比对结果的可靠性和生物学意义。 比对得分标准化 ~~~~~~~~~~~~~~ 原始比对得分依赖于所使用的打分矩阵和空位罚分参数，不同参数设置下的得分不能直接比较。比特分（Bit score）通过标准化消除了打分参数的影响，其计算公式为S' = (λS - lnK) / ln2，其中S是原始得分，λ和K是Karlin-Altschul统计量。比特分使得不同搜索之间的结果可以相互比较，当你需要比较不同BLAST搜索的结果时，应该使用比特分而不是原始得分。 统计显著性估计 ~~~~~~~~~~~~~~ BLAST比对的统计显著性基于极值分布（Gumbel分布）理论。Karlin-Altschul统计理论证明，随机序列比对得分服从极值分布，E值的计算公式为E = Kmn·e^(-λS)，其中m和n分别是查询序列和数据库的长度，S是比对得分。E值表示在随机情况下期望得到的得分不低于当前得分的命中数。P值与E值的关系为P = 1 - e^(-E)，当E值很小时，P值约等于E值。当你解读BLAST结果时，E值是最重要的统计指标，E值越小表示比对结果越不可能由随机因素产生。 .. code-block:: python import numpy as np from scipy.stats import gumbel_r m = 1000 n = 10000000 K = 0.13 lambda_param = 0.318 raw_score = 50 E_value = K * m * n * np.exp(-lambda_param * raw_score) bit_score = (lambda_param * raw_score - np.log(K)) / np.log(2) P_value = 1 - np.exp(-E_value) print(f"查询序列长度: {m}") print(f"数据库长度: {n:,}") print(f"原始得分: {raw_score}") print(f"E值: {E_value:.2e}") print(f"比特分: {bit_score:.2f}") print(f"P值: {P_value:.2e}") raw_scores = np.linspace(20, 80, 100) E_values = K * m * n * np.exp(-lambda_param * raw_scores) significant = E_values < 0.05 print(f"E值<0.05所需最低原始得分: {raw_scores[significant][0]:.1f}") 多重检验校正是当进行多次比对时需要考虑的统计问题。如果你对1000条序列分别进行BLAST搜索，即使每次搜索的E值阈值为0.05，也可能得到50个假阳性结果。Bonferroni校正是最保守的校正方法，将显著性阈值除以检验次数。FDR（False Discovery Rate）校正控制假发现率，比Bonferroni校正更宽松，在大规模分析中更为常用。 比对敏感性与特异性 ~~~~~~~~~~~~~~~~~~ 比对敏感性（sensitivity）是指真实同源序列被正确检出的比例，特异性（specificity）是指检出的同源序列中真正同源的比例。敏感性和特异性之间存在权衡关系，提高敏感性通常会降低特异性，反之亦然。当你调整BLAST的E值阈值时，降低阈值会提高特异性但降低敏感性，提高阈值则相反。在实际应用中，需要根据研究目的选择合适的平衡点：在需要全面检索同源序列的探索性研究中，可以适当放宽阈值；在需要高置信度结果的确认性研究中，应该使用严格的阈值。 比对评估基准 ~~~~~~~~~~~~ 比对评估基准数据集用于客观评估比对算法的性能。BAliBASE是蛋白质多序列比对的标准基准数据集，包含多种结构参考比对。HomFam基于Homstrad结构比对数据库构建。SABmark提供超家族和扭曲家族两个难度级别的比对基准。PREFAB结合了结构比对和序列比对信息。当你需要比较不同比对工具的性能或选择最适合你数据的工具时，这些基准数据集提供了客观的评估依据。 比对评分函数用于量化比对结果的质量。Sum-of-pairs得分计算比对中所有序列对的匹配得分之和。Column Score计算正确比对列的比例。TC Score（Total Column Score）计算完全正确的比对列的比例。比对一致性评估通过比较计算比对与参考比对的一致性来衡量比对质量。 长读长比对是近年来随着第三代测序技术发展而出现的新需求。长读长数据（如PacBio和Oxford Nanopore）的错误率较高（5-15%），且读长可达数十kb，传统的短读长比对工具不适用。minimap2是目前最常用的长读长比对工具，它采用最小化子（minimizer）索引和动态规划比对相结合的策略，能够高效处理高错误率的长读长数据。lra是另一个长读长比对工具，专门针对长读长的错误模式进行了优化。