基因组学分析 ============ 基因组学分析是生物信息学中最核心的应用领域之一，涵盖了从测序数据获取到基因组解读的完整流程。当你拿到一份测序数据时，首先需要了解这些数据是如何产生的、质量如何，然后需要将它们比对到参考基因组上或者从头组装出基因组序列，接下来需要对基因组进行注释以识别其中的基因和功能元件，最后还需要检测变异、比较不同基因组之间的差异，揭示进化规律和功能差异。本章将系统介绍高通量测序技术、测序数据预处理、序列比对、基因组组装、基因组注释、变异检测、比较基因组学和泛基因组学等内容，帮助读者掌握基因组学分析的完整技术体系。 高通量测序技术 -------------- 高通量测序技术，也称为下一代测序（Next-Generation Sequencing，NGS）技术，是当代基因组学研究的基石。理解测序技术的基本原理和特点，对于正确选择测序策略、设计实验方案以及解读分析结果都至关重要。不同的测序平台在读长、通量、错误率和成本方面各有优劣，选择合适的测序平台是实验设计的第一步。 第一代测序（Sanger测序） ~~~~~~~~~~~~~~~~~~~~~~~~ Sanger测序，也称为双脱氧链终止法，是由Frederick Sanger于1977年发明的DNA测序方法，是人类基因组计划所使用的主要测序技术。Sanger测序的基本原理是在DNA合成反应中加入四种双脱氧核苷酸（ddNTP），每种ddNTP带有不同的荧光标记。当ddNTP被掺入到正在延伸的DNA链中时，由于它缺少3'羟基，DNA链的延伸就会终止。通过电泳分离不同长度的DNA片段，并根据荧光信号读取每个位置的碱基，就可以得到DNA序列。Sanger测序的读长通常为700-1000bp，准确率非常高，错误率低于0.1%，因此至今仍然是验证性测序的金标准。当你需要验证一个克隆的序列、确认PCR产物的序列、或者对少量样本进行Sanger测序时，这种方法仍然是最可靠的选择。然而，Sanger测序的通量很低，每次反应只能读取一条序列，成本也相对较高，不适合大规模基因组测序项目。 第二代测序（Illumina测序） ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 第二代测序技术的代表是Illumina平台，它采用桥式PCR扩增和边合成边测序（Sequencing by Synthesis，SBS）的策略。在Illumina测序中，DNA文库片段首先与流动池（flow cell）表面的接头序列互补结合，然后通过桥式PCR在原位扩增形成簇（cluster），每个簇包含数千个相同的DNA分子拷贝。测序时，加入带有荧光标记的可逆终止子dNTP，每个循环只有一种碱基被掺入，通过拍摄荧光图像确定每个簇当前掺入的碱基，然后去除荧光基团和终止基团，进入下一个循环。Illumina测序的读长通常为75-300bp，通量极高，单次运行可以产生数十亿条读段，错误率约为0.1-1%，主要以替换错误为主。Illumina平台是目前应用最广泛的测序平台，当你需要进行全基因组测序、转录组测序、表观基因组测序等大规模测序项目时，Illumina通常是首选平台。它的主要优势在于高通量和低单碱基成本，主要局限在于读长较短，对于高度重复区域和结构变异的检测能力有限。 第三代测序（单分子实时测序） ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 第三代测序技术的代表是PacBio（Pacific Biosciences）的单分子实时测序（Single Molecule Real-Time，SMRT）技术。PacBio测序不需要对DNA分子进行扩增，而是直接观察单个DNA分子的合成过程。在SMRT测序中，DNA聚合酶被固定在零模波导孔（Zero-Mode Waveguide，ZMW）底部，当带有荧光标记的dNTP被掺入到正在延伸的DNA链中时，荧光信号被实时检测。PacBio测序的读长非常长，平均读长可以达到10-25kb，最长可超过100kb，这使得它在基因组组装、结构变异检测和全长转录本测序方面具有独特优势。早期PacBio测序的主要问题是单次读取的错误率较高（约10-15%），以插入缺失错误为主。然而，通过循环一致性测序（Circular Consensus Sequencing，CCS）技术，即对同一条环状DNA分子进行多次测序并生成一致性序列，PacBio HiFi读段的准确率可以达到99.9%（Q30以上），同时保持10-25kb的读长。当你需要获得高质量的基因组组装、检测复杂结构变异、或者获取全长转录本序列时，PacBio HiFi测序是目前最佳的选择之一。 第四代测序（纳米孔测序） ~~~~~~~~~~~~~~~~~~~~~~~~~ 第四代测序技术的代表是Oxford Nanopore Technologies（ONT）的纳米孔测序。纳米孔测序的原理是将DNA分子驱动通过一个蛋白质纳米孔，当不同的碱基通过纳米孔时，会引起离子电流的特异性变化，通过检测这些电流变化就可以识别通过纳米孔的碱基序列。纳米孔测序的读长极长，平均读长可以达到10-100kb，最长可超过2Mb，是目前所有测序技术中读长最长的。纳米孔测序的另一个显著优势是设备便携，MinION测序仪只有U盘大小，可以通过USB接口连接笔记本电脑使用，这使得现场测序和实时分析成为可能。当你需要在野外进行病原体快速鉴定、在临床环境中进行实时基因组监测、或者需要超长读长来解析基因组中最复杂的重复区域时，纳米孔测序是理想的选择。纳米孔测序的错误率约为5-10%，以插入缺失和替换错误为主，但通过最新的化学试剂和碱基识别算法的改进，准确率正在持续提升。此外，纳米孔测序还可以直接检测DNA修饰（如5mC、6mA），无需亚硫酸盐处理，这为表观遗传学研究提供了新的工具。 测序平台对比 ~~~~~~~~~~~~ 不同的测序平台在读长、通量、错误率、成本和运行时间等方面存在显著差异，选择合适的平台需要根据具体的研究目标和预算来决定。Sanger测序读长中等（700-1000bp），准确率最高（>99.9%），但通量最低，单碱基成本最高，运行时间约数小时到一天。Illumina测序读长最短（75-300bp），准确率高（>99%），通量最高，单碱基成本最低，运行时间从数小时到数天不等。PacBio HiFi测序读长长（10-25kb），准确率高（>99.9%），通量中等，单碱基成本中等，运行时间约一天。Oxford Nanopore测序读长最长（10kb-2Mb+），准确率较低（90-99%，最新化学试剂可达Q20+），通量中等，单碱基成本中等，运行时间从数分钟到数天（可实时读取）。在实际研究中，往往需要根据具体需求选择平台，例如基因组组装优先考虑长读长平台，大规模群体重测序优先考虑Illumina平台，而现场快速检测则考虑纳米孔平台。近年来，多种平台联合使用的混合策略也越来越普遍，例如用长读长数据辅助组装、用短读长数据校正错误。 .. code-block:: python platforms = { "Sanger": {"读长": "700-1000bp", "准确率": ">99.9%", "通量": "低", "单碱基成本": "高", "运行时间": "数小时-1天"}, "Illumina": {"读长": "75-300bp", "准确率": ">99%", "通量": "极高", "单碱基成本": "低", "运行时间": "数小时-数天"}, "PacBio HiFi": {"读长": "10-25kb", "准确率": ">99.9%", "通量": "中等", "单碱基成本": "中等", "运行时间": "约1天"}, "ONT": {"读长": "10kb-2Mb+", "准确率": "90-99%(Q20+)", "通量": "中等", "单碱基成本": "中等", "运行时间": "分钟-数天"} } print(f"{'平台':<15} {'读长':<15} {'准确率':<15} {'通量':<8} {'成本':<8}") print("-" * 65) for name, info in platforms.items(): print(f"{name:<15} {info['读长']:<15} {info['准确率']:<15} {info['通量']:<8} {info['单碱基成本']:<8}") 双端测序与单端测序 ~~~~~~~~~~~~~~~~~~ 双端测序（Paired-End Sequencing）和单端测序（Single-End Sequencing）是Illumina平台的两种基本测序模式。单端测序只从DNA片段的一端进行测序，产生一条读段。双端测序则从DNA片段的两端分别进行测序，产生两条读段，这两条读段之间有一段已知大致长度的间隔（插入片段大小，insert size）。双端测序相比单端测序提供了更多的信息：两条读段的配对关系可以帮助确定读段在基因组上的相对位置，提高比对的准确性；当一条读段无法唯一比对到基因组时，其配对读段的位置信息可以帮助确定正确的比对位置；双端读段之间的距离和方向信息可以用于检测结构变异，如缺失、插入、倒位等。当你进行基因组重测序、转录组分析或结构变异检测时，双端测序通常是更好的选择。单端测序的成本较低，适用于一些对比对精度要求不高的应用场景，如小RNA测序、ChIP-seq等。在选择测序模式时，需要权衡信息量和成本，大多数基因组学应用推荐使用双端测序。 测序深度与覆盖度 ~~~~~~~~~~~~~~~~ 测序深度（Sequencing Depth）和覆盖度（Coverage）是评估测序数据量的两个关键指标。测序深度是指基因组中每个碱基被测序读段覆盖的平均次数，通常用X表示，例如30X表示基因组中每个碱基平均被30条读段覆盖。覆盖度是指基因组中被至少一条读段覆盖的碱基占基因组总碱基数的比例，通常用百分比表示。测序深度和覆盖度是不同的概念，30X的测序深度并不意味着100%的覆盖度，因为测序读段在基因组上的分布是不均匀的，某些区域可能被多次覆盖，而另一些区域可能完全没有被覆盖到。不同应用对测序深度的要求不同：人类全基因组测序通常需要30X以上的深度，外显子组测序通常需要50-100X，而检测体细胞突变可能需要数百甚至上千X的深度。当你设计测序实验时，需要根据研究目标和样本类型确定合适的测序深度，深度过低可能导致变异漏检，深度过高则增加成本而不一定带来等比例的信息增益。覆盖度不足的区域通常出现在GC含量极端的区域、重复序列区域和端粒附近，这些区域的测序覆盖是基因组学分析中的一个持续挑战。 .. code-block:: python genome_size = 3_000_000_000 read_length = 150 total_reads = 600_000_000 sequencing_depth = (total_reads * read_length) / genome_size print(f"基因组大小: {genome_size/1e6:.0f} Mb") print(f"读段总数: {total_reads/1e6:.0f} M") print(f"读段长度: {read_length} bp") print(f"测序深度: {sequencing_depth:.1f}X") import math coverage_probability = 1 - math.exp(-sequencing_depth) print(f"理论覆盖度: {coverage_probability*100:.2f}%") 文库构建类型 ~~~~~~~~~~~~ 不同的研究目的需要构建不同类型的测序文库，文库类型决定了测序数据所代表的基因组区域和分子信息。全基因组测序（Whole Genome Sequencing，WGS）文库覆盖整个基因组，适用于基因组组装、变异检测等需要全面基因组信息的分析。全外显子组测序（Whole Exome Sequencing，WES）文库通过探针捕获技术富集基因组中的外显子区域，仅对编码蛋白质的区域进行测序，成本远低于WGS，适用于大规模群体中蛋白质编码区变异的筛查。RNA-seq文库来自反转录的cDNA，用于基因表达量分析和转录本结构鉴定。ChIP-seq文库来自染色质免疫沉淀后的DNA片段，用于检测蛋白质与DNA的相互作用，如转录因子结合位点和组蛋白修饰位置。ATAC-seq文库来自转座酶可接近的染色质区域，用于绘制开放染色质图谱。Hi-C文库捕获空间上邻近但线性距离可能很远的DNA片段，用于研究染色质三维结构。Methyl-seq文库经过亚硫酸盐处理或酶切富集，用于检测DNA甲基化。RAD-seq（Restriction-site Associated DNA Sequencing）通过限制性内切酶消化基因组DNA，仅对酶切位点附近的片段进行测序，适用于群体遗传学和系统发育研究中的SNP发现。当你需要根据研究目标选择合适的文库类型时，需要综合考虑研究问题、样本数量、预算和所需的数据分辨率。 捕获测序与扩增子测序 ~~~~~~~~~~~~~~~~~~~~ 捕获测序和扩增子测序是两种常用的靶向测序策略，它们都用于对基因组的特定区域进行深度测序，但实现方式不同。捕获测序（也称为杂交捕获）通过设计特异性探针与目标区域互补杂交，将目标DNA片段从文库中富集出来。捕获测序可以同时靶向数百个基因或数Mb的基因组区域，灵活性高，适用于中等规模到大规模的靶向测序项目。当你需要对大量样本的外显子组或特定基因集进行测序时，捕获测序是常用的选择。扩增子测序通过PCR扩增目标区域，然后对扩增产物进行测序。扩增子测序的设计更加简单，成本更低，适合对小规模目标区域（如几个到几十个基因）进行超高深度测序。在临床基因检测中，扩增子测序常用于已知致病基因的热点区域筛查。捕获测序的优势在于目标区域范围灵活、覆盖均匀性较好，但成本较高、实验周期较长。扩增子测序的优势在于成本低、操作简单、深度高，但目标区域受限于引物设计，且PCR偏好性可能导致覆盖不均匀。 多重测序与Barcode/Index ~~~~~~~~~~~~~~~~~~~~~~~~ 多重测序（Multiplex Sequencing）是在同一个测序通道中同时测序多个样本的技术，通过给每个样本添加唯一的条形码序列（barcode或index）来区分不同样本的读段。在文库构建过程中，每个样本的DNA片段两端会加上带有特定barcode序列的接头，测序完成后根据barcode序列将读段分配回各自的样本，这个过程称为demultiplexing。多重测序大大降低了每个样本的测序成本，提高了测序仪的利用率，是大规模群体测序项目的基础。当你需要对数十到数百个样本进行测序时，多重测序是必需的技术手段。barcode的设计需要考虑碱基平衡性，即每个位置的四种碱基比例应尽量均衡，以确保测序过程中荧光信号的正常检测。常见的barcode设计策略包括Illumina的i7和i5双index系统，以及更灵活的UMI（Unique Molecular Identifier）系统。UMI是在PCR扩增前给每个原始DNA分子添加的唯一分子标签，通过UMI可以识别和去除PCR重复，提高变异检测的准确性，在液体活检和低频突变检测中尤为重要。 测序质量控制指标 ~~~~~~~~~~~~~~~~ 测序质量控制是评估测序数据是否合格的关键步骤，常用的质量指标包括Q30、簇密度和碱基平衡。Q30是指Phred质量分数大于30的碱基所占的比例，Phred质量分数Q表示碱基识别的错误概率P，关系为Q=-10log10(P)，Q30对应的错误概率为0.1%，即99.9%的准确率。Illumina平台通常要求Q30比例达到75%以上才算合格。簇密度是指流动池上每个mm2的簇数量，过高或过低的簇密度都会影响测序质量，簇密度过高会导致簇之间的信号串扰，簇密度过低则浪费测序通量。碱基平衡是指每个测序循环中四种碱基的比例分布，良好的碱基平衡对于Illumina测序的图像分析至关重要，因为不平衡的碱基分布会导致信号校准困难。当你收到一批测序数据时，首先应该检查这些质量控制指标，判断数据是否满足分析要求，不合格的数据可能需要重新测序或采取额外的质量控制措施。 .. code-block:: python import math def phred_to_error_prob(q): return 10 ** (-q / 10) for q in [10, 20, 30, 40]: error_prob = phred_to_error_prob(q) accuracy = (1 - error_prob) * 100 print(f"Q{q}: 错误概率 = {error_prob:.4f} ({error_prob*100:.2f}%), 准确率 = {accuracy:.2f}%") total_bases = 150_000_000 q30_bases = 120_000_000 q30_percentage = (q30_bases / total_bases) * 100 print(f"\nQ30比例: {q30_percentage:.1f}% ({'合格' if q30_percentage >= 75 else '不合格'})") 测序数据预处理 -------------- 测序数据预处理是基因组学分析流程中承上启下的关键环节。原始测序数据通常包含接头序列、低质量碱基、PCR重复等各种技术噪音，这些噪音会严重影响后续比对和变异检测的准确性。数据预处理的目标是在保留尽可能多的真实生物学信号的前提下，去除技术噪音，生成高质量的分析数据。一个典型的预处理流程包括质量评估、质量控制、去除接头、修剪低质量碱基、去除短读段、去除PCR重复以及污染检测等步骤。 原始数据质量评估 ~~~~~~~~~~~~~~~~ 原始数据质量评估是预处理的第一步，也是最基础的质控环节。FastQC是目前最常用的测序数据质量评估工具，它可以对FASTQ文件进行全面的质控分析，生成包含多个质量指标的报告。FastQC报告的主要模块包括：每个位置的碱基质量分布（Per Base Sequence Quality）、每个读段的平均质量分布（Per Sequence Quality Score）、碱基含量分布（Per Base Sequence Content）、GC含量分布（Per Sequence GC Content）、读段长度分布（Sequence Length Distribution）、重复读段比例（Sequence Duplication Levels）、过度表达的序列（Overrepresented Sequences）以及接头含量（Adapter Content）等。当你拿到一批新的测序数据时，首先运行FastQC查看质量报告是一个好习惯，这可以帮助你了解数据的整体质量状况，发现潜在的问题，并据此制定合适的质控策略。当你需要同时评估多个样本的质量时，MultiQC可以将多个FastQC报告汇总成一个综合报告，方便批量查看和比较。 .. code-block:: bash # FastQC质量评估 fastqc -t 8 -o qc_reports/ sample_R1.fastq.gz sample_R2.fastq.gz # MultiQC汇总多个FastQC报告 multiqc -o multiqc_report/ qc_reports/ # 查看FastQC结果摘要 ls qc_reports/ # 输出: sample_R1_fastqc.html sample_R1_fastqc.zip # sample_R2_fastqc.html sample_R2_fastqc.zip 质量控制与修剪 ~~~~~~~~~~~~~~ 质量控制与修剪是去除测序数据中技术噪音的核心步骤。Trimmomatic是一个功能全面的质控工具，它可以执行滑动窗口质量修剪、去除接头序列、裁剪读段两端低质量碱基以及过滤短读段等多种操作。cutadapt专门用于去除接头序列和引物序列，支持多种接头类型和错误容忍设置。fastp是一个集成了质量评估、质量修剪、接头去除、重复过滤等多种功能的快速质控工具，它的速度比Trimmomatic快数倍，同时还能生成详细的质控报告。fastx_toolkit是一组用于处理FASTQ/FASTA文件的小工具集，包括质量格式转换、质量修剪、序列过滤等功能。当你需要对大规模测序数据进行质控时，fastp因其速度优势是很好的选择；当需要精细控制质控参数时，Trimmomatic提供了更灵活的配置选项。质控参数的选择需要根据数据质量和分析目标来确定，过于严格的质控可能丢失有效数据，过于宽松的质控则可能保留低质量数据影响后续分析。 .. code-block:: bash # Trimmomatic质控 trimmomatic PE -threads 8 \ sample_R1.fastq.gz sample_R2.fastq.gz \ sample_R1_paired.fq.gz sample_R1_unpaired.fq.gz \ sample_R2_paired.fq.gz sample_R2_unpaired.fq.gz \ ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \ SLIDINGWINDOW:4:20 \ LEADING:20 TRAILING:20 \ MINLEN:50 # fastp一键质控（更快） fastp -i sample_R1.fastq.gz -I sample_R2.fastq.gz \ -o sample_R1_clean.fq.gz -O sample_R2_clean.fq.gz \ --detect_adapter_for_pe \ --cut_front --cut_tail \ --cut_window_size 4 --cut_mean_quality 20 \ --length_required 50 \ --thread 8 \ --html fastp_report.html 去除接头 ~~~~~~~~ 接头序列是文库构建过程中添加到DNA片段两端的已知序列，用于在流动池上的扩增和测序引物结合。当测序读段比DNA插入片段长时，测序会读穿插入片段进入接头区域，产生所谓的接头污染。接头污染在短插入片段文库中尤为常见，如果不去除接头序列，这些序列会干扰后续的比对和分析。AdapterRemoval是一个专门用于去除接头序列和合并重叠双端读段的工具，它在去除接头的同时还可以将重叠的双端读段合并为一条更长的读段，这对于古代DNA等短片段样本特别有用。Skewer是另一个快速的接头去除工具，支持部分接头匹配和模糊匹配，在处理接头污染时更加灵活。当你发现FastQC报告中接头含量过高时，需要使用这些工具去除接头序列。去除接头后，建议再次运行FastQC确认接头已被有效去除。 .. code-block:: bash # AdapterRemoval去除接头并合并重叠双端读段 AdapterRemoval --file1 sample_R1.fq.gz --file2 sample_R2.fq.gz \ --adapter1 AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \ --adapter2 AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT \ --collapse \ --minlength 50 \ --output1 sample_R1_trimmed.fq.gz \ --output2 sample_R2_trimmed.fq.gz \ --outputcollapsed sample_collapsed.fq.gz # Skewer快速接头去除 skewer -x TruSeq3-PE.fa -t 8 -m pe \ sample_R1.fq.gz sample_R2.fq.gz \ -o sample_trimmed 修剪低质量碱基 ~~~~~~~~~~~~~~ 测序过程中碱基识别的准确率随着测序循环数的增加而逐渐降低，这导致读段末端通常含有较多低质量碱基。Phred质量分数是衡量碱基识别可靠性的标准指标，Q20表示错误概率1%，Q30表示错误概率0.1%。滑动窗口法是修剪低质量碱基的常用策略，它从读段的5'端或3'端开始，以固定窗口大小计算窗口内的平均质量分数，当平均质量低于设定阈值时，从该位置开始截断读段。这种方法比简单截断固定长度更加智能，因为它根据实际质量状况决定修剪位置。当你处理Illumina测序数据时，通常需要修剪3'端的低质量碱基，因为Illumina测序的质量衰减主要发生在读段的3'端。修剪阈值通常设置为Q20或Q15，具体取决于数据的整体质量和后续分析的精度要求。 .. code-block:: bash # cutadapt修剪低质量碱基 cutadapt -q 20,20 \ --minimum-length 50 \ -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \ -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT \ -o sample_R1_trimmed.fq.gz \ -p sample_R2_trimmed.fq.gz \ sample_R1.fq.gz sample_R2.fq.gz # 使用Python解析Phred质量分数 .. code-block:: python def parse_fastq_quality(quality_string, offset=33): scores = [ord(c) - offset for c in quality_string] return scores quality_string = "IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII" scores = parse_fastq_quality(quality_string) print(f"碱基数: {len(scores)}") print(f"平均质量分数: {sum(scores)/len(scores):.1f}") print(f"最低质量分数: {min(scores)}") 去除短读长 ~~~~~~~~~~ 经过接头去除和质量修剪后，一些读段可能变得非常短，这些短读段在比对时容易产生非特异性匹配，增加假阳性比对的风险，同时也会浪费计算资源。因此，质控流程中通常会设置一个最小读长阈值，将短于该阈值的读段过滤掉。最小读长阈值的选择需要根据具体应用来确定，对于基因组比对，通常设置为30-50bp，因为过短的读段很难在基因组上找到唯一比对位置。对于转录组分析，考虑到某些外显子可能很短，可以适当放宽阈值。当你设置过滤阈值时，需要在保留足够数据量和保证比对准确性之间取得平衡，过高的阈值会丢弃过多数据，过低的阈值则可能引入低质量比对。 去除重复读段 ~~~~~~~~~~~~ PCR扩增是文库构建过程中的必要步骤，但扩增会导致同一原始DNA分子产生多个相同的拷贝，这些重复读段（PCR duplicates）在变异检测中会造成等位基因频率的偏差，需要被去除。Picard MarkDuplicates是最常用的去重工具，它通过比对读段的比对位置和方向来识别PCR重复，保留每个重复组中质量最高的一条读段。对于全基因组测序数据，去重是GATK最佳实践流程的标准步骤。对于RNA-seq数据，由于基因表达量的差异，高表达基因自然会产生大量相同的读段，这些不是PCR重复而是生物学重复，因此RNA-seq分析通常不去除重复读段。当你进行DNA测序数据的变异检测时，去除PCR重复是必要的质控步骤，但需要注意，去重会降低有效测序深度，因此实验设计时需要考虑到去重后的深度是否仍然满足分析需求。 .. code-block:: bash # Picard MarkDuplicates去除PCR重复 picard MarkDuplicates \ I=aligned_sorted.bam \ O=aligned_dedup.bam \ M=duplication_metrics.txt \ REMOVE_DUPLICATES=true # 查看去重统计信息 cat duplication_metrics.txt # samtools去重（更快） samtools markdup -r aligned_sorted.bam aligned_dedup.bam # 使用UMI去重（更精确） umi_tools dedup \ -I aligned_sorted.bam \ --output-stats=dedup_stats \ -S aligned_umi_dedup.bam 污染检测及去除 ~~~~~~~~~~~~~~ 测序数据中可能混入非目标物种的DNA序列，这些污染序列可能来自实验过程中的交叉污染、环境微生物污染或宿主DNA污染。Kraken是一个基于k-mer精确匹配的分类工具，可以快速将读段分类到不同的物种，用于检测测序数据中的微生物污染。FastQ_Screen通过将读段比对到多个参考基因组来检测污染来源，它可以告诉你数据中有多少比例的读段可以比对到人、小鼠、大肠杆菌等常见污染源基因组。DeconSeq通过比对去除已知污染物种的序列。当你处理微生物组测序数据时，宿主DNA去除是一个关键步骤，例如在人类微生物组研究中，可能需要去除90%以上的人类DNA读段。污染检测应该在质控流程的早期进行，以便及时发现数据质量问题并采取相应措施。 .. code-block:: bash # Kraken2污染检测 kraken2 --db kraken2_db \ --paired sample_R1.fq.gz sample_R2.fq.gz \ --report kraken_report.txt \ --output kraken_output.txt # 查看污染物种分布 kraken2-report --db kraken2_db kraken_report.txt | head -20 # FastQ_Screen检测污染来源 fastq_screen --conf fastq_screen.conf \ --paired sample_R1.fq.gz sample_R2.fq.gz \ --outdir screen_results/ 测序错误的统计与校正 ~~~~~~~~~~~~~~~~~~~~ 测序错误校正是对原始测序数据进行错误修复的高级质控步骤。与简单的质量修剪不同，错误校正试图利用读段之间的重叠信息来纠正碱基识别错误，从而在保留读段长度的同时提高碱基准确率。BayesHammer是SPAdes组装器中集成的错误校正模块，它基于贝叶斯统计模型对k-mer频谱进行分析，区分真实k-mer和错误k-mer，然后对读段中的错误碱基进行校正。Lighter是一个无需参考基因组的快速错误校正工具，它通过采样k-mer来估计k-mer频谱，然后进行校正。BLESS利用k-mer频谱的分布特征来识别和纠正错误。当你进行基因组从头组装时，错误校正可以显著提高组装质量，特别是对于高错误率的长读长数据，错误校正是组装前的必要步骤。对于Illumina短读长数据，错误校正对组装质量的提升相对有限，因为Illumina本身的错误率已经较低。 过滤宿主来源序列 ~~~~~~~~~~~~~~~~ 在微生物组测序、病原体检测和宏基因组学研究中，测序数据中往往含有大量宿主DNA序列，这些宿主序列需要被去除以减少对微生物序列分析的干扰。BMTagger是一个常用的宿主序列过滤工具，它通过将读段比对到宿主参考基因组来识别和去除宿主来源的读段。Bowtie2也可以用于宿主序列过滤，通过将读段先比对到宿主基因组，然后提取未比对的读段作为非宿主序列。当你进行人类微生物组研究时，口腔、肠道等样本中人类DNA的比例可能高达90%以上，有效去除宿主序列对于提高微生物序列的覆盖深度至关重要。在临床病原体检测中，宿主序列去除同样重要，因为临床样本中宿主DNA通常占绝对多数。过滤宿主序列时需要注意选择合适的比对参数，避免将微生物序列错误地比对到宿主基因组上而被误去除。 质量转换格式 ~~~~~~~~~~~~ 在基因组学分析中，不同的分析步骤可能需要不同的文件格式，格式转换是日常操作中常见的需求。FASTQ是原始测序数据的标准格式，包含序列信息和碱基质量信息。SAM（Sequence Alignment/Map）是比对结果的标准文本格式，BAM是SAM的二进制压缩格式，CRAM是比BAM更高效的压缩格式。uBAM（unmapped BAM）是一种特殊的BAM格式，保留了原始测序数据的信息但使用BAM格式存储，常用于GATK流程中。当你使用GATK进行变异检测时，通常需要将FASTQ转换为uBAM格式，然后进行比对和后处理。samtools是格式转换的主要工具，使用view子命令可以在不同格式之间转换。FASTQ到SAM/BAM的转换通常通过比对工具（如BWA、Bowtie2）完成，而SAM/BAM到FASTQ的转换可以使用samtools fastq命令。理解这些格式之间的转换关系对于构建完整的分析流程非常重要。 高通量测序数据比对到参考基因组 ------------------------------ 将测序读段比对到参考基因组是基因组学分析中最基础的步骤之一。比对的目的是确定每条读段在参考基因组上的来源位置，为后续的变异检测、基因表达定量等分析提供位置信息。比对算法的效率和准确性直接影响整个分析流程的质量，理解比对算法的基本原理有助于正确选择比对工具和优化比对参数。 比对算法基础 ~~~~~~~~~~~~ 现代短读长比对算法的核心数据结构是Burrows-Wheeler变换（BWT）和FM-index。Burrows-Wheeler变换是一种对字符串进行可逆变换的方法，变换后的字符串具有更好的压缩性和可搜索性。FM-index是基于BWT构建的全文索引结构，它利用BWT的特性和辅助数据结构实现了高效的精确匹配搜索。FM-index的搜索算法基于向后搜索（backward search），从模式串的末尾开始逐字符向前搜索，每一步利用字符出现次数表（C数组）和_occurrence_数组来缩小搜索范围。FM-index的空间效率极高，人类基因组的FM-index仅需约3GB内存，使得在普通服务器上进行全基因组比对成为可能。哈希表是另一种常用的索引结构，它将参考基因组的所有k-mer（固定长度的子串）及其位置存储在哈希表中，查询时通过k-mer直接查找位置。哈希表的查询速度非常快，但空间占用较大。种子扩展策略结合了哈希表和BWT的优点，首先通过种子（通常是k-mer）快速定位候选比对位置，然后通过动态规划算法在种子周围进行扩展，寻找最优比对。当你理解了这些算法原理后，就能更好地理解不同比对工具的特点和适用场景。 短读长比对工具 ~~~~~~~~~~~~~~ BWA-MEM是最常用的短读长比对工具，它基于BWT/FM-index实现了高效的种子扩展比对算法。BWA-MEM自动选择最优的种子长度，支持嵌合体比对（chimeric alignment），对长读段（70bp-1Mbp）都有良好的表现，是基因组重测序和变异检测的标准比对工具。Bowtie2是另一个广泛使用的短读长比对工具，它同样基于FM-index，但采用了不同的种子策略，支持间隙比对和双端读段比对，在灵敏度和速度之间提供了灵活的调节选项。STAR是专为RNA-seq数据设计的剪接感知比对工具，它使用后缀数组索引，可以高效地识别跨越外显子-内含子边界的剪接读段，是RNA-seq分析的标准比对工具。HISAT2是另一个剪接感知比对工具，它使用分层索引（FM-index + 局部索引）来处理大规模基因组，内存占用远低于STAR，同时保持较高的比对速度和灵敏度。Minimap2虽然主要用于长读长比对，但也支持短读长比对模式。当你进行DNA测序数据的变异检测时，BWA-MEM通常是首选；进行RNA-seq分析时，STAR或HISAT2是更好的选择。 .. code-block:: bash # BWA-MEM比对（DNA测序标准流程） bwa index reference.fasta bwa mem -t 16 -M reference.fasta \ sample_R1_clean.fq.gz sample_R2_clean.fq.gz \ | samtools view -bS - > aligned.bam # Bowtie2比对 bowtie2-build reference.fasta ref_index bowtie2 -x ref_index -1 sample_R1_clean.fq.gz -2 sample_R2_clean.fq.gz \ -S aligned.sam --threads 16 # STAR RNA-seq剪接感知比对 STAR --runMode genomeGenerate \ --genomeDir star_index \ --genomeFastaFiles reference.fasta \ --sjdbGTFfile annotation.gtf STAR --genomeDir star_index \ --readFilesIn sample_R1.fq.gz sample_R2.fq.gz \ --readFilesCommand zcat \ --outSAMtype BAM SortedByCoordinate \ --runThreadN 16 # HISAT2剪接感知比对 hisat2-build reference.fasta ref_index hisat2 -x ref_index -1 sample_R1.fq.gz -2 sample_R2.fq.gz \ -S aligned.sam -p 16 长读长比对工具 ~~~~~~~~~~~~~~ 长读长比对需要处理高错误率和长读段的特殊挑战，传统的短读长比对算法不能直接适用。Minimap2是目前最流行的长读长比对工具，它使用最小化子（minimizer）索引策略，可以快速将PacBio和Nanopore读段比对到参考基因组，同时也支持剪接感知比对和基因组间比对。NGMLR是专门为长读长结构变异检测设计的比对工具，它考虑了长读长中常见的连续插入缺失错误，在结构变异断点附近的比对质量优于Minimap2。GraphMap是一个高灵敏度的Nanopore读段比对工具，使用基于图的比对策略。lra是一个新兴的长读长比对工具，支持PacBio和Nanopore数据，在比对速度和准确性方面都有不错的表现。BLASR是PacBio早期开发的比对工具，现在已逐渐被Minimap2取代。当你使用长读长数据进行结构变异检测时，Minimap2和NGMLR是主要选择；进行全长转录本比对时，Minimap2的剪接感知模式可以处理内含子跨越。 .. code-block:: bash # Minimap2长读长比对 minimap2 -ax map-pb -t 16 reference.fasta pacbio_reads.fasta \ | samtools sort -o aligned_pb.bam minimap2 -ax map-ont -t 16 reference.fasta nanopore_reads.fasta \ | samtools sort -o aligned_ont.bam # Minimap2剪接感知比对（全长转录本） minimap2 -ax splice -t 16 reference.fasta isoform_reads.fasta \ | samtools sort -o aligned_isoform.bam # NGMLR长读长比对（适合结构变异检测） ngmlr -t 16 -r reference.fasta -q nanopore_reads.fasta \ -o aligned_ngmlr.bam 比对参数优化 ~~~~~~~~~~~~ 比对工具的参数设置对比对结果有显著影响，合理的参数优化可以提高比对的准确性和效率。错配罚分（mismatch penalty）控制比对中允许的碱基不匹配数量，较低的错配罚分允许更多的错配，提高灵敏度但可能增加假阳性比对。开放空位罚分（gap open penalty）和延伸空位罚分（gap extension penalty）控制比对中插入缺失的惩罚，较高的空位罚分会减少间隙的数量，适合错误率较低的数据。对于Illumina短读长数据，通常使用默认参数即可获得良好的比对结果，因为Illumina数据的错误率较低且以替换错误为主。对于PacBio和Nanopore长读长数据，需要降低错配罚分和空位罚分以容忍较高的测序错误率。当你比对不同来源的数据时，需要根据数据的错误特征调整比对参数，例如Nanopore数据的插入缺失错误较多，需要降低空位罚分；而RNA-seq数据需要启用剪接感知模式。 比对后处理 ~~~~~~~~~~ 比对生成的原始SAM文件通常需要经过一系列后处理步骤才能用于下游分析。排序（sorting）是将比对记录按照参考基因组上的位置排列，这是大多数下游工具的要求。索引（indexing）是为排序后的BAM文件创建索引，使得可以快速随机访问特定基因组区域的比对记录。统计量计算包括总比对率、唯一比对率、双端读段比对率等指标，用于评估比对质量。samtools是比对后处理的核心工具集，samtools view用于格式转换（SAM/BAM/CRAM之间），samtools sort用于排序，samtools index用于创建索引，samtools flagstat用于快速统计比对信息，samtools stats用于生成详细的比对统计报告。当你完成比对后，首先应该运行samtools flagstat查看基本的比对统计信息，确认比对率是否在预期范围内，然后进行排序和索引，为后续分析做准备。 .. code-block:: bash # 比对后处理完整流程 # 1. 排序 samtools sort -@ 16 -o aligned_sorted.bam aligned.bam # 2. 索引 samtools index aligned_sorted.bam # 3. 快速统计 samtools flagstat aligned_sorted.bam # 4. 详细统计 samtools stats aligned_sorted.bam > samtools_stats.txt # 5. 格式转换 samtools view -b aligned_sorted.bam > output.bam samtools view -C -T reference.fasta aligned_sorted.bam > output.cram 比对质量评估 ~~~~~~~~~~~~ 比对质量评估是确保比对结果可靠性的重要步骤。Qualimap是一个综合性的比对质量评估工具，它可以计算覆盖度分布、插入片段大小分布、GC偏倚、比对质量分布等多种指标，并生成可视化的报告。samtools stats提供详细的比对统计信息，包括读段数量、比对率、插入片段大小、碱基分布等。samtools flagstat快速汇总比对的总体情况，如总读段数、比对读段数、正确配对的读段数等。samtools idxstats报告每个参考序列的读段数量和长度信息，可以检测染色体级别的覆盖偏差。当你评估比对质量时，需要关注几个关键指标：总体比对率应该达到70-90%以上（具体取决于样本类型），覆盖度应该相对均匀（GC偏倚不应过大），插入片段大小分布应该与文库构建的预期一致。如果比对率过低，可能是样本质量问题或参考基因组选择不当；如果覆盖度严重不均匀，可能需要考虑GC偏倚校正。 .. code-block:: bash # Qualimap比对质量评估 qualimap bamqc -bam aligned_sorted.bam -outdir qualimap_report # samtools详细统计 samtools flagstat aligned_sorted.bam samtools idxstats aligned_sorted.bam samtools coverage aligned_sorted.bam # 计算覆盖度分布 samtools depth -a aligned_sorted.bam | \ awk '{sum+=$3; count++} END {print "平均覆盖度:", sum/count}' 比对结果过滤 ~~~~~~~~~~~~ 比对结果中可能包含低质量的比对记录，这些记录会干扰后续的变异检测和其他分析。常见的过滤策略包括：去除未比对的读段（flag 4），去除低质量比对（MAPQ低于设定阈值，通常设为20或30），去除PCR重复（使用Picard MarkDuplicates），去除次级比对（flag 256）和补充比对（flag 2048）。MAPQ（Mapping Quality）是比对质量的量化指标，表示比对位置正确的概率的对数值，MAPQ=0表示比对位置不可靠，MAPQ=60表示比对位置非常可靠。当你进行变异检测时，通常需要过滤MAPQ<20的比对记录，因为低质量比对可能来自重复区域或测序错误。对于RNA-seq数据，过滤策略需要更加谨慎，因为一些基因家族成员之间的序列高度相似，过于严格的过滤可能丢失真实的比对。 .. code-block:: bash # 过滤未比对和低质量比对 samtools view -b -q 20 -F 4 -F 256 aligned_sorted.bam > aligned_filtered.bam # 过滤参数说明 # -q 20: MAPQ >= 20 # -F 4: 去除未比对读段 # -F 256: 去除次级比对 # -F 2048: 去除补充比对 # 使用sambamba过滤（更快） sambamba view -h -t 8 -f bam \ -F "mapping_quality >= 20 and not (unmapped or secondary_alignment or supplementary)" \ aligned_sorted.bam > aligned_filtered.bam 比对可视化 ~~~~~~~~~~ 比对可视化是核查比对结果和验证变异的重要手段。IGV（Integrative Genomics Viewer）是最常用的本地比对可视化工具，它可以加载BAM文件和参考基因组，以图形化方式展示读段在基因组上的比对情况，支持查看单个位点的碱基组成、插入缺失、剪接位点等信息。当你需要确认一个变异是否真实存在时，在IGV中查看该位点的比对情况是最直接的方法。UCSC Genome Browser和Ensembl Browser是基于Web的基因组浏览器，可以将比对数据上传到服务器上进行在线查看。JBrowse是一个可以本地部署的基因组浏览器，适合在服务器上查看大规模比对数据。可视化检查在变异验证中尤为重要，特别是对于结构变异和低频变异，仅依赖计算工具的输出可能不够可靠，人工查看比对模式可以提供额外的信心。 多比对读段处理 ~~~~~~~~~~~~~~ 多比对读段（multi-mapping reads）是指可以比对到参考基因组多个位置的读段，这种情况通常发生在重复序列区域和基因家族成员之间。在SAM/BAM格式中，多比对读段通过标志位和标签来区分：初级比对（primary alignment，flag 0或16）是比对工具认为最可能的比对位置，次级比对（secondary alignment，flag 256）是其他可能的比对位置，补充比对（supplementary alignment，flag 2048）是嵌合体比对中非主要部分的比对。XA标签列出了所有比对位置，HI标签指示比对的层级。当你处理RNA-seq数据时，多比对读段的处理尤其重要，因为同源基因之间的读段分配直接影响表达量估计。常见的处理策略包括：丢弃所有多比对读段（简单但丢失信息）、按比对质量分数分配（按比例分配给各比对位置）、或使用EM算法估计最可能的来源（如RSEM和Salmon的方法）。 嵌合体比对检测 ~~~~~~~~~~~~~~ 嵌合体比对（chimeric alignment）是指一条读段的一部分比对到基因组的一个位置，另一部分比对到距离较远的另一个位置。嵌合体比对可能由多种原因产生：融合基因（fusion gene）导致一条读段跨越两个不同基因的连接处、结构变异（如缺失、倒位）导致读段的一部分比对到异常位置、剪接读段在非剪接感知比对中的表现、以及测序错误或文库构建问题。BWA-MEM通过supplementary alignment标记嵌合体比对，STAR在RNA-seq比对中专门报告嵌合体比对用于融合基因检测。当你进行肿瘤转录组分析寻找融合基因时，嵌合体比对是主要的信号来源，但需要仔细过滤假阳性，因为基因组结构变异和反式剪接也会产生嵌合体比对信号。 比对格式及压缩与索引 ~~~~~~~~~~~~~~~~~~~~ SAM/BAM/CRAM是比对数据的三种主要格式。SAM是文本格式，人类可读但占用空间大。BAM是SAM的二进制压缩格式，占用空间约为SAM的1/3到1/5，是日常分析中最常用的格式。CRAM是比BAM更高效的压缩格式，通过利用参考基因组序列进一步压缩数据，占用空间约为BAM的1/2到1/3。SAM/BAM格式的核心字段包括：QNAME（读段名称）、FLAG（比对标志位）、RNAME（参考序列名称）、POS（比对起始位置）、MAPQ（比对质量）、CIGAR（比对操作字符串）、RNEXT（配对读段参考序列）、PNEXT（配对读段位置）、TLEN（模板长度）、SEQ（序列）、QUAL（质量）。CIGAR字符串使用M（匹配/错配）、I（插入）、D（缺失）、N（跳过，如内含子）、S（软裁剪）、H（硬裁剪）等操作符描述比对细节。BAM文件的索引（BAI文件）使用B-tree结构存储每个基因组区域的偏移量，使得可以快速随机访问特定区域的比对记录。当你处理大规模比对数据时，使用CRAM格式可以显著节省存储空间，但需要注意CRAM文件依赖参考基因组序列，迁移数据时需要同时保存参考基因组。 基因组组装 ---------- 基因组组装是从测序读段重建原始基因组序列的过程，是基因组学研究中最重要的计算挑战之一。当没有参考基因组可供比对时（如新物种的基因组研究），从头组装（de novo assembly）是获得基因组序列的唯一途径。即使有参考基因组可用，从头组装也能发现参考基因组中不存在的结构变异和新序列。基因组组装的质量直接影响后续所有分析的可靠性，因此理解组装算法的原理和评估组装质量的方法至关重要。 组装基础概念 ~~~~~~~~~~~~ 基因组组装涉及一系列基本概念，理解这些概念是掌握组装算法的前提。读长（read length）是指单条测序读段的碱基数，读长越长，包含的信息量越大，组装越容易。K-mer是指长度为k的DNA子串，在德布鲁因图组装中，k-mer是构建图的基本单元。覆盖度（coverage）是指每个碱基被测序读段覆盖的平均次数，高覆盖度提供了更多的冗余信息用于纠错和解决组装歧义。组装图是组装算法的核心数据结构，主要有两种类型：德布鲁因图（de Bruijn graph）和重叠-布局-一致性图（Overlap-Layout-Consensus graph，OLC图）。德布鲁因图将每条读段分解为一系列k-mer，然后以k-mer之间的重叠关系构建有向图，适用于短读长组装。OLC图首先计算所有读段之间的两两重叠（overlap），然后确定读段的排列顺序（layout），最后生成一致性序列（consensus），适用于长读长组装。读长重叠（overlap）是指两条读段之间存在一段相同的序列，重叠的长度和相似度是判断重叠是否可靠的标准。Contig是由读段拼接而成的连续序列片段，是组装的基本输出单位。Scaffold是将contig按照顺序和方向排列，并用N填充contig之间的间隙（gap）形成的更大结构。N50是评估组装连续性的核心指标，它表示将所有contig按长度从大到小排列后，累计长度达到总长度50%时的contig长度。N90和L50是类似的指标，N90表示累计长度达到90%时的contig长度，L50表示达到50%总长度所需的contig数量。组装大小与基因组大小估计可以通过k-mer频谱分析来实现，k-mer频谱的峰值位置对应测序深度，峰值的宽度反映覆盖度的变异，而杂合位点和重复序列会在频谱中产生额外的峰。GenomeScope是一个基于k-mer频谱估计基因组大小、杂合度和重复序列比例的工具。组装复杂性的度量主要考虑两个因素：杂合度和重复序列比例，高杂合度会导致等位基因被组装成独立的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 = [50000, 35000, 28000, 20000, 15000, 12000, 8000, 5000, 3000, 1000] n50 = calculate_n50(contig_lengths) n90 = calculate_n90(contig_lengths) total_length = sum(contig_lengths) print(f"总长度: {total_length:,} bp") print(f"N50: {n50:,} bp") print(f"Contig数量: {len(contig_lengths)}") def calculate_n90(lengths): sorted_lengths = sorted(lengths, reverse=True) total = sum(sorted_lengths) cumsum = 0 for length in sorted_lengths: cumsum += length if cumsum >= total * 0.9: return length return 0 print(f"N90: {calculate_n90(contig_lengths):,} bp") 短读长组装 ~~~~~~~~~~ 短读长组装主要基于德布鲁因图算法。德布鲁因图组装的基本流程是：首先将所有读段分解为k-mer，然后以k-1 mer为节点、k-mer为边构建德布鲁因图，接着在图上进行路径遍历寻找欧拉路径，最终将路径转换为contig序列。SPAdes是目前最常用的短读长组装器，它采用多k-mer策略，先使用较小的k-mer构建图以捕获低覆盖度区域，再使用较大的k-mer解决重复区域，最后将不同k-mer的组装结果进行整合。SPAdes特别适合细菌基因组的组装，也支持宏基因组和转录组数据的组装。ABySS是一个分布式组装器，可以利用集群的多个计算节点并行处理大规模数据，适合大型基因组的组装。SOAPdenovo2是华大基因开发的短读长组装器，专门针对大型基因组优化，在人类基因组组装中表现良好。Velvet是早期的德布鲁因图组装器，虽然现在使用较少，但它的算法原理对理解德布鲁因图组装非常有帮助。最优k-mer选择是短读长组装的关键参数，KmerGenie可以通过扫描不同k-mer值下的独特k-mer数量来推荐最优k-mer大小。读长纠错是组装前的重要步骤，BayesHammer和Musket是两种常用的纠错工具，纠错可以减少德布鲁因图中的错误分支，提高组装质量。脱泡（bubble）去除是德布鲁因图简化的重要操作，气泡通常由测序错误或杂合位点产生，去除气泡可以简化图结构。Scaffolding是将contig连接为scaffold的过程，SSPACE、Opera和SOAPdenovo的scaffolding模块利用双端读段和mate-pair读段的配对信息来确定contig之间的顺序和方向。Gap闭合是填补scaffold中N间隙的步骤，GapFiller、GapCloser和Sealer通过将读段比对到间隙两端的contig上，延伸序列来闭合间隙。当你组装一个细菌基因组时，SPAdes通常是首选工具；对于大型真核生物基因组，短读长组装的连续性通常不够理想，需要结合长读长数据。 .. code-block:: bash # SPAdes细菌基因组组装 spades.py -1 sample_R1_clean.fq.gz -2 sample_R2_clean.fq.gz \ -o spades_output --threads 16 --memory 64 # ABySS分布式组装 export NP=16 abyss-pe k=64 name=abyss_assembly \ in='sample_R1.fq.gz sample_R2.fq.gz' # KmerGenie最优k-mer选择 kmergenie sample_R1.fq.gz,sample_R2.fq.gz -o kmergenie_report # GapCloser闭合间隙 GapCloser -a scaffold.fasta -b config.txt -o gapclosed.fasta 长读长组装 ~~~~~~~~~~ 长读长组装主要基于重叠-布局-一致性（OLC）策略，因为长读长之间的重叠关系更加可靠，不需要将读段分解为k-mer。Canu是PacBio和Nanopore数据的标准组装器，它的工作流程包括三个阶段：纠错（利用读段之间的重叠对原始读段进行纠错）、修剪（去除读段中不可靠的区域）和组装（构建重叠图并寻找路径）。Flye是一个通用的长读长组装器，它使用重复图（repeat graph）来处理重复序列，即使原始读段错误率较高也能产生较好的组装结果。Miniasm是一个极简的长读长组装器，它不做纠错，直接利用原始读段之间的重叠构建组装图，速度非常快但组装错误率较高。Hifiasm是专门为PacBio HiFi数据设计的组装器，它利用HiFi读段的高准确率实现了分相组装（phased assembly），可以分别组装出两套单倍型序列，是目前PacBio HiFi数据组装的最佳选择。HiCanu是Canu的HiFi模式，针对HiFi数据的低错误率进行了优化。Verkko是另一个HiFi组装器，它结合了HiFi读段和Oxford Nanopore超长读段，可以产生端到端的染色体级别组装。对于Nanopore数据，Shasta是一个快速的Nanopore组装器，Flye的纳米孔模式也表现良好，NECAT则针对Nanopore数据的错误特征进行了优化。长读长自纠错是长读长组装的重要步骤，Canu的纠错模块通过多序列比对生成一致性序列来纠错，Racon则通过快速一致性计算来纠错。三重叠群一致性纠错（Triple Consensus Correction）是在组装后进一步提高序列准确率的步骤，Racon通过将读段比对回组装结果进行多轮纠错，Medaka使用神经网络对Nanopore组装结果进行纠错，Homopolish利用不同平台数据的互补性进行纠错。当你组装一个新物种的基因组时，如果有PacBio HiFi数据，Hifiasm是首选；如果只有Nanopore数据，Flye是可靠的选择。 .. code-block:: bash # Hifiasm PacBio HiFi组装（分相组装） hifiasm -o asm -t 32 -l0 pacbio_hifi_reads.fasta.gz # 提取单倍型组装结果 awk '/^S/{print ">"$2;print $3}' asm.bp.hap1.p_ctg.gfa > hap1.fasta awk '/^S/{print ">"$2;print $3}' asm.bp.hap2.p_ctg.gfa > hap2.fasta # Flye Nanopore组装 flye --nano-raw nanopore_reads.fasta.gz \ --genome-size 3g --threads 32 \ --out-dir flye_output # Canu长读长组装 canu -p asm -d canu_output \ genomeSize=3g \ -pacbio-raw pacbio_reads.fasta.gz # Racon纠错 minimap2 -x map-ont reference.fasta ont_reads.fasta | \ sort -k6,6 -k7,7n | \ racon ont_reads.fasta - assembly.fasta > corrected.fasta # Medaka Nanopore组装纠错 medaka_consensus -i ont_reads.fasta -d assembly.fasta \ -o medaka_output -t 16 -m r941_min_high_g360 混合组装 ~~~~~~~~ 混合组装策略结合了短读长和长读长的优势，短读长提供高准确率，长读长提供长距离连续性。DBG2OLC是一种混合组装策略，它首先用短读长构建德布鲁因图获得contig，然后将长读长与contig进行比对，利用长读长作为骨架将contig连接起来。MaSuRCA是另一个混合组装器，它将短读长拼接为超读段（super-read），然后与长读长一起进行OLC组装。SPAdes的hybrid模式可以利用长读长数据来改善短读长组装的连续性。桥接组装策略利用长读长来连接短读长组装产生的contig，长读长可以跨越contig之间的重复区域，提供连接信息。10x Genomics链接读长（Linked Reads）技术通过微滴分区给同一DNA分子的短读段添加相同的barcode，Supernova组装器利用这些barcode信息来辅助scaffold构建。光学图谱（Optical Map）技术如Bionano可以提供超长距离的限制性酶切位点信息，Bionano Solve可以将光学图谱与序列组装进行比对，辅助scaffold构建和错误检测。Hi-C辅助染色体挂载是获得染色体级别组装的关键技术，Hi-C数据反映了基因组三维空间中远距离的相互作用，3D-DNA、YaHS和SALSA等工具利用Hi-C数据的接触矩阵将scaffold排列到染色体上。当你需要将contig级别的组装提升到染色体级别时，Hi-C辅助挂载是目前最经济有效的方法。 .. code-block:: bash # Hi-C辅助染色体挂载（3D-DNA） # 第一步：比对Hi-C数据 bwa mem -5SP reference.fasta hic_R1.fq.gz hic_R2.fq.gz | \ pairtools parse -c reference.chrom.sizes -o pairs.pairs pairtools sort -o pairs.sorted.pairs pairs.pairs pairtools dedup --mark-dups -o pairs.dedup.pairs pairs.sorted.pairs # 第二步：3D-DNA挂载 run-3d-dna.sh reference.fasta pairs.dedup.pairs # YaHS Hi-C挂载（更快） yahs reference.fasta pairs.dedup.pairs # Bionano光学图谱辅助 bionano_solve scaffold.fasta bionano.cmap output_dir 组装评估 ~~~~~~~~ 组装完成后需要对组装质量进行全面评估，以确定组装结果是否满足下游分析的要求。QUAST是最常用的组装评估工具，它可以计算N50、总长度、最大contig长度等基本指标，还可以通过与参考基因组比对来评估组装的准确性和完整性。BUSCO是评估基因组组装基因空间完整性的标准工具，它通过搜索近缘物种的单拷贝直系同源基因来评估组装中包含的完整基因比例，BUSCO结果通常用C（完整）、S（单拷贝）、D（重复）、F（片段化）和M（缺失）五个指标表示。CEGMA是BUSCO的前身，使用核心真核基因来评估完整性，现在已被BUSCO取代。基于读长回比对的评估方法将原始测序读段比对回组装结果，统计覆盖度和覆盖均匀性，covN50是基于覆盖度的N50指标，可以反映组装的覆盖完整性。基于k-mer谱的评估方法比较原始读段的k-mer频谱与组装结果的k-mer频谱，Merqury是这一方法的代表工具，它可以计算组装的k-mer一致性和错误率，不需要参考基因组。评估指标还包括错误组装数量、结构变异准确性等，QUAST可以检测错误组装（misassembly），包括重排错误（relocation）、倒位错误（inversion）和转位错误（translocation）。组装完整性评估中，BUSCO的单拷贝直系基因比例是最常用的指标，完整单拷贝基因比例达到90%以上通常被认为是较好的组装。当你评估组装质量时，需要综合考虑连续性（N50）、完整性（BUSCO）和准确性（QUAST、Merqury）三个维度，单一指标不能全面反映组装质量。 .. code-block:: bash # QUAST组装质量评估 quast.py assembly.fasta -r reference.fasta -g annotation.gtf -o quast_report # BUSCO基因完整性评估 busco -i assembly.fasta -l insecta_odb10 -m genome -o busco_result -c 16 # Merqury k-mer一致性评估 merqury.sh reads.meryl assembly.meryl merqury_output # 计算组装统计信息 .. code-block:: python from collections import namedtuple def assembly_stats(contig_lengths): total = sum(contig_lengths) sorted_lengths = sorted(contig_lengths, reverse=True) num_contigs = len(sorted_lengths) max_len = sorted_lengths[0] min_len = sorted_lengths[-1] avg_len = total / num_contigs def nx(x): cumsum = 0 for l in sorted_lengths: cumsum += l if cumsum >= total * x: return l return 0 print(f"总长度: {total:,} bp") print(f"Contig数量: {num_contigs}") print(f"最长: {max_len:,} bp") print(f"最短: {min_len:,} bp") print(f"平均: {avg_len:,.0f} bp") print(f"N50: {nx(0.5):,} bp") print(f"N90: {nx(0.9):,} bp") assembly_stats([50000, 35000, 28000, 20000, 15000, 12000, 8000, 5000, 3000, 1000]) 基因组注释 ---------- 基因组注释是在基因组序列上识别和标记生物学功能元件的过程，是基因组学研究中从序列到功能的关键桥梁。一个未经注释的基因组只是一串ATCG字母，只有通过注释才能了解其中包含哪些基因、这些基因编码什么蛋白质、基因组中有哪些重复序列和调控元件等信息。基因组注释通常分为三个层次：重复序列注释、基因结构注释和功能注释，这三个层次逐步深入，从识别基因组的基本组成到赋予其生物学意义。 重复序列注释 ~~~~~~~~~~~~ 重复序列是基因组中含量最丰富的功能元件，在人类基因组中约占50%以上。重复序列注释是基因组注释的第一步，因为重复序列会干扰基因预测，需要先将其识别和屏蔽，然后再进行基因预测。重复序列可以分为几大类：转座元件（Transposable Elements，TEs）是能够从基因组一个位置移动到另一个位置的DNA序列，是重复序列中最主要的类型；简单重复（Simple Repeats）是由1-6个碱基的短序列串联重复组成，也称为微卫星；串联重复（Tandem Repeats）是较长的重复单元首尾相连排列，包括小卫星；散在重复（Interspersed Repeats）是分散在基因组各处的重复序列，主要包括LINE、SINE和LTR反转录转座子；低复杂度区域（Low Complexity）是碱基组成严重偏离随机分布的序列。 从头重复识别（de novo repeat identification）是在没有已知重复序列数据库的情况下，通过序列之间的相似性比较来发现新的重复序列家族。RepeatModeler是目前最常用的从头重复序列建模工具，它整合了RECON和RepeatScout两种算法，可以自动发现和分类重复序列家族。LTRharvest专门用于识别LTR反转录转座子，它通过搜索LTR对的特征来发现完整的LTR元件。MITE Hunter专门用于识别MITE（微型反向重复转座元件）。基于数据库的注释是将基因组序列与已知的重复序列数据库进行比对，RepeatMasker是最常用的重复序列屏蔽工具，它使用RMBlast或CrossMatch搜索引擎将基因组序列与RepBase和DFAM数据库进行比对，然后屏蔽匹配的重复区域。RepBase是人工审校的重复序列数据库，主要包含真核生物的重复序列；DFAM是基于HMM的转座元件数据库，覆盖面更广。重复序列家族分类与命名遵循Wicker等人的分类系统，将转座元件分为I类（反转录转座子，通过RNA中间体转座）和II类（DNA转座子，通过DNA中间体转座），每个类别进一步分为多个超家族和家族。重复序列掩码（masking）是将注释出的重复序列替换为N或转换为小写字母，以便后续基因预测时跳过这些区域。串联重复注释使用TRF（Tandem Repeats Finder）来识别基因组中的串联重复序列，TRF通过检测序列的周期性模式来发现串联重复。LTR反转录转座子完整鉴定需要同时检测5'和3'LTR以及它们之间的内部序列，LTR_FINDER和LTR_retriever是两个互补的工具，LTR_retriever整合多个LTR预测工具的结果，可以高精度地鉴定完整的LTR反转录转座子。当你注释一个新物种的基因组时，重复序列注释是必不可少的步骤，特别是对于植物基因组，重复序列比例可能高达80%以上。 .. code-block:: bash # RepeatModeler从头重复序列建模 BuildDatabase -name my_species -engine ncbi assembly.fasta RepeatModeler -database my_species -pa 16 -LTRStruct # RepeatMasker重复序列屏蔽 RepeatMasker -pa 16 -lib consensi.fa.classified \ -gff -dir repeatmasker_output assembly.fasta # TRF串联重复注释 trf assembly.fasta 2 7 7 80 10 50 2000 -d -h # LTR_retriever完整LTR鉴定 LTR_FINDER -seq assembly.fasta -harvestout -time 600 ltr_retriever -genome assembly.fasta -inharvest LTR_FINDER.harvest.out 基因结构注释 ~~~~~~~~~~~~ 基因结构注释是识别基因组中基因的位置、结构和边界的过程，包括确定基因的起始和终止位置、外显子-内含子边界、转录起始和终止位点等信息。基因结构注释的方法可以分为三大类：从头预测、同源预测和基于转录组证据的预测，这三类方法各有优缺点，通常需要整合多种证据来获得可靠的基因模型。 从头预测方法仅依赖基因组序列本身的统计特征来预测基因，不需要外部证据。隐马尔可夫模型（Hidden Markov Model，HMM）是经典的从头基因预测方法，Augustus和GlimmerHMM是两个基于HMM的基因预测工具，它们通过训练物种特异的HMM参数来识别编码区域和剪接位点。基于深度学习的方法近年来取得了显著进展，GeneMark-ETP结合了深度学习和自训练策略，Helixer使用深度神经网络直接从基因组序列预测基因结构。基于内容感应的方法如SNAP和Genie通过分析序列的编码特征来预测基因。BRAKER是一个集成了GeneMark和Augustus的自动化基因预测流程，它可以利用RNA-seq数据来训练和优化预测参数。 同源预测方法利用已知物种的基因信息来预测目标物种的基因。基于蛋白质比对的方法将已知蛋白质序列比对到基因组上，推断基因结构，Exonerate和GeneWise是两个常用的蛋白质比对基因预测工具，GeneWise特别适合处理移码和测序错误。基于EST/RNA-seq比对的方法将转录组数据比对到基因组上，直接提供基因结构的证据，PASA是一个整合转录组比对结果来构建基因模型的工具，TransDecoder从转录本序列中预测开放阅读框。基于已知物种同源的方法如GeMoMa利用近缘物种的基因注释来预测目标物种的基因。 证据整合是将多种预测方法的结果合并为一致的基因模型，EvidenceModeler（EVM）是最常用的证据整合工具，它为不同来源的预测结果分配权重，通过加权投票确定最终的基因模型。MAKER是一个集成了重复序列屏蔽、基因预测和证据整合的自动化注释流程。GATE是另一个证据整合框架。基因模型过滤与修正是整合后的必要步骤，需要去除不完整的基因模型、修正异常的基因结构。非编码RNA注释包括tRNA（使用tRNAscan-SE）、rRNA和其他非编码RNA（使用Infernal搜索Rfam数据库），mikado是一个整合多种转录本注释的工具。假基因鉴定使用Pseudopipe和PseudoFinder来识别失去功能的基因拷贝。UTR和剪接位点预测确定基因的非翻译区和剪接位点位置。当你进行基因组注释时，整合多种证据来源是获得高质量注释的关键，单一方法的预测结果往往不够可靠。 .. code-block:: bash # BRAKER自动化基因预测流程 braker.pl --genome=assembly.fasta \ --bam=rna_seq_aligned.bam \ --species=my_species \ --cores=16 # Augustus基因预测（需先训练） augustus --species=human assembly.masked.fasta > augustus.gff # PASA转录组证据整合 Launch_PASA_pipeline.pl -c alignAssembly.config -C -R \ --genome assembly.fasta \ --transcripts transcriptome.fasta # EvidenceModeler证据整合 evm.py --weights weights.txt \ --genome assembly.fasta \ --ab_initio augustus.gff \ --protein_evidence protein_align.gff \ --transcript_evidence pasa.gff # tRNAscan-SE tRNA注释 tRNAscan-SE -B -o tRNA.out assembly.fasta 功能注释 ~~~~~~~~ 功能注释是赋予基因结构注释产物生物学功能的过程，它回答的是基因做什么的问题。功能注释通常从序列比对开始，通过将预测的蛋白质序列与已知功能的蛋白质数据库进行比对，推断其可能的功能。 序列比对到UniProt/Swiss-Prot/TrEMBL是最基本的功能注释步骤。Swiss-Prot是人工审校的高质量蛋白质数据库，注释信息准确但覆盖面有限；TrEMBL是自动注释的蛋白质数据库，覆盖面广但注释质量参差不齐。当你将预测的蛋白质序列与Swiss-Prot进行BLAST比对时，如果找到高相似度的匹配，就可以将匹配蛋白的功能转移给目标蛋白。基因本体（Gene Ontology，GO）注释为基因提供标准化的功能描述，GO分为三个本体：分子功能（Molecular Function）、生物过程（Biological Process）和细胞组分（Cellular Component）。InterProScan通过搜索多个蛋白质特征数据库来为蛋白质分配GO术语，eggNOG-mapper基于直系同源群来注释功能，Blast2GO通过BLAST搜索结果来分配GO术语。KEGG通路映射将基因映射到代谢和信号通路中，KAAS是KEGG官方的注释服务器，KofamKOALA使用HMM谱来分配KO编号，GhostKOALA适用于大规模数据集的KEGG注释。蛋白质结构域注释识别蛋白质中的功能模块，Pfam是基于HMM的蛋白质家族数据库，InterPro整合了多个特征数据库，SMART专注于信号和调控结构域，PROSITE基于模式识别。酶学委员会编号（EC编号）注释将基因映射到酶学分类系统，PRIAM和EFICAz是两个专门的EC编号注释工具。亚细胞定位预测确定蛋白质在细胞中的位置，TargetP预测信号肽和转运肽，SignalP专门预测信号肽，TMHMM预测跨膜螺旋，WoLF PSORT综合预测亚细胞定位。基因符号自动分配根据同源性将标准基因符号分配给预测的基因，AnnotationHub是R/Bioconductor中的注释资源。注释文件输出格式主要有三种：GFF3（General Feature Format version 3）是最通用的基因组注释格式，GTF（Gene Transfer Format）主要用于基因模型描述，GBK（GenBank格式）整合了序列和注释信息。当你完成功能注释后，这些注释信息将成为后续功能基因组学分析的基础。 .. code-block:: bash # InterProScan功能注释 interproscan.sh -i proteins.fasta -f tsv,xml,gff3 \ -dp -goterms -iprlookup -pa -t p -cpu 16 # eggNOG-mapper功能注释 emapper.py -i proteins.fasta --output eggnog_result \ --cpu 16 --data_dir eggnog_db # BLAST比对Swiss-Prot blastp -query proteins.fasta -db swissprot \ -outfmt 6 -evalue 1e-5 -max_target_seqs 1 -num_threads 16 \ -out blast_swissprot.tsv # SignalP信号肽预测 signalp -fasta proteins.fasta -org euk -format short -gff3 # TMHMM跨膜区域预测 tmhmm --short < proteins.fasta > tmhmm_result.txt 注释质量评估与工作流 ~~~~~~~~~~~~~~~~~~~~ 注释质量评估是确保注释结果可靠性的重要环节。BUSCO不仅可以评估组装质量，也可以评估注释基因集的完整性，通过检查注释基因中包含的完整单拷贝直系同源基因比例来评估注释的完整性。注释的组装一致性检查使用GffCompare来比较注释与转录组比对结果的一致性，检测遗漏的基因和错误的外显子边界。注释的可视化使用JBrowse、Apollo和IGV等工具，JBrowse适合在Web上展示基因组注释，Apollo支持多人协作编辑注释，IGV适合本地查看注释细节。注释更新与版本控制是长期维护基因组注释的重要实践，随着新证据的出现，注释需要不断更新和改进。大规模注释工作流将上述步骤自动化，Galaxy提供了图形化的注释流程界面，GeneMark-ETP+是GeneMark的增强版自动化流程，BRAKER流程整合了GeneMark和Augustus的自动化注释。当你需要注释一个新基因组时，使用成熟的自动化工作流可以大大提高效率，但仍需人工检查关键基因的注释质量。 变异检测 -------- 变异检测是从测序数据中识别基因组变异的过程，是医学基因组学和群体遗传学的核心技术。从单核苷酸变异到大规模结构变异，不同类型的变异需要不同的检测策略和工具。变异检测的准确性直接影响下游的功能解读和临床决策，因此理解变异检测的原理和最佳实践至关重要。 变异类型与表示 ~~~~~~~~~~~~~~ 基因组变异按照大小和类型可以分为多种类别。单核苷酸变异（Single Nucleotide Variant，SNV）是单个碱基的改变，当SNV在群体中的频率超过1%时称为单核苷酸多态性（SNP）。短插入缺失（Indel）是1-50bp的插入或缺失，是仅次于SNV的最常见变异类型。结构变异（Structural Variant，SV）是大于50bp的基因组改变，包括缺失（deletion）、重复（duplication）、倒位（inversion）、易位（translocation）和插入（insertion）等类型。拷贝数变异（Copy Number Variant，CNV）是基因组区域拷贝数的改变，可以是缺失型（拷贝数减少）或重复型（拷贝数增加）。线粒体变异是线粒体DNA上的变异，由于线粒体DNA的高拷贝数和异质性，线粒体变异检测有其特殊性。复杂变异包括多核苷酸多态性（MNP，连续多个碱基同时改变）和混合型变异。 变异的VCF/BCF格式是存储变异信息的标准格式。VCF（Variant Call Format）是文本格式，BCF是VCF的二进制压缩格式。VCF文件的核心字段包括：CHROM（染色体名称）、POS（变异位置，1-based）、ID（变异标识符，如rs编号）、REF（参考碱基）、ALT（变异碱基）、QUAL（变异质量分数）、FILTER（过滤标记）、INFO（附加信息）和FORMAT（样本格式）。INFO字段包含丰富的变异注释信息，如AC（等位基因计数）、AF（等位基因频率）、AN（等位基因总数）、DP（测序深度）等。FORMAT字段定义了每个样本的基因型信息格式，如GT（基因型）、AD（等位基因深度）、DP（样本深度）、GQ（基因型质量）等。变异标准化是将变异表示为最简形式的过程，vt和bcftools norm可以将左对齐和拆分多等位基因位点，确保变异表示的一致性。当你处理变异检测数据时，理解VCF格式的各个字段含义是正确解读变异信息的基础。 .. code-block:: bash # bcftools变异标准化 bcftools norm -f reference.fasta -Oz -o normalized.vcf.gz raw_variants.vcf.gz # 拆分多等位基因位点 bcftools norm -m -any -Oz -o split_variants.vcf.gz normalized.vcf.gz # 解析VCF文件 .. code-block:: python def parse_vcf_line(line): if line.startswith("#"): return None fields = line.strip().split("\t") chrom, pos, var_id, ref, alt, qual, filt, info, fmt = fields[:9] samples = fields[9:] info_dict = {} for item in info.split(";"): if "=" in item: key, value = item.split("=", 1) info_dict[key] = value else: info_dict[item] = True return { "chrom": chrom, "pos": int(pos), "id": var_id, "ref": ref, "alt": alt, "qual": float(qual) if qual != "." else None, "filter": filt, "info": info_dict, "format": fmt, "samples": samples } vcf_line = "chr1\t12345\trs123\tA\tG\t999\tPASS\tAC=2;AF=0.333;DP=30\tGT:AD:DP:GQ" result = parse_vcf_line(vcf_line) print(f"位置: {result['chrom']}:{result['pos']}") print(f"变异: {result['ref']}>{result['alt']}") print(f"等位基因频率: {result['info'].get('AF', 'N/A')}") 种系变异检测 ~~~~~~~~~~~~ 种系变异（germline variant）是从父母遗传获得的变异，存在于个体的所有细胞中。短读长种系SNV/Indel检测的主流工具包括GATK HaplotypeCaller、FreeBayes、BCFtools和Samtools mpileup。GATK HaplotypeCaller采用局部组装策略，在候选变异区域进行de Bruijn图组装，然后通过配对隐马尔可夫模型确定单倍型和变异，是GATK最佳实践流程的核心工具，也是目前最广泛使用的种系变异检测工具。FreeBayes基于贝叶斯统计模型进行变异检测，支持多等位基因位点和池化样本。BCFtools call基于mpileup的碱基计数进行变异检测，速度较快。Samtools mpileup生成碱基堆叠信息，配合bcftools call进行变异检测。 长读长种系变异检测利用长读长的优势来检测短读长难以发现的变异。DeepVariant使用深度神经网络进行变异检测，支持多种测序平台的数据。Clair3是另一个基于深度学习的长读长变异检测工具，速度比DeepVariant快。NanoCaller专门针对Nanopore数据优化。Sniffles2是长读长结构变异检测的标准工具，它可以检测缺失、插入、倒位、重复和易位等多种结构变异类型。单样本与多样本联合检测是群体变异检测的两种策略，单样本检测对每个样本独立进行变异检测，然后合并结果；多样本联合检测同时分析多个样本，利用样本间的信息提高检测灵敏度和基因型准确性。GATK GenomicsDBImport将多个样本的gVCF导入GenomicsDB数据库，然后使用GenotypeGVCFs进行联合基因分型。BCFtools的mpileup模式支持多样本同时检测。群体VCF生成是群体遗传学分析的基础，通常需要将多个样本的变异信息合并为一个群体VCF文件。当你进行大规模群体测序项目的变异检测时，多样本联合检测是推荐的做法，因为它可以提高低频变异的检测灵敏度。 .. code-block:: bash # GATK HaplotypeCaller种系变异检测 gatk HaplotypeCaller \ -R reference.fasta \ -I aligned_dedup.bam \ -O sample.g.vcf.gz \ -ERC GVCF # 多样本联合基因分型 gatk GenomicsDBImport \ --genomicsdb-workspace-path my_database \ --batch-size 50 \ -V sample1.g.vcf.gz -V sample2.g.vcf.gz -V sample3.g.vcf.gz gatk GenotypeGVCFs \ -R reference.fasta \ -V gendb://my_database \ -O cohort.vcf.gz # BCFtools快速变异检测 bcftools mpileup -f reference.fasta -Ou aligned_dedup.bam | \ bcftools call -mv -Oz -o variants.vcf.gz 体细胞变异检测 ~~~~~~~~~~~~~~ 体细胞变异（somatic variant）是在个体生命过程中获得的变异，主要存在于肿瘤细胞中。体细胞变异检测通常需要肿瘤-正常配对样本设计，即同时测序同一个体的肿瘤组织和正常组织，通过比较两者来识别肿瘤特异的体细胞变异。这种配对设计可以有效去除种系变异的干扰，确保检测到的变异是肿瘤特异的。 体细胞SNV/Indel检测的主流工具包括Mutect2、Strelka2、VarScan2和VarDictJava。Mutect2是GATK中的体细胞变异检测工具，它采用与HaplotypeClient类似的局部组装策略，同时使用配对正常样本和群体等位基因频率来过滤种系变异。Strelka2采用差分概率模型比较肿瘤和正常样本的等位基因频率，检测灵敏度较高。VarScan2使用启发式方法比较肿瘤和正常样本的碱基计数，支持多种变异类型。VarDictJava是一个灵活的体细胞变异检测工具，特别适合临床应用。 体细胞拷贝数变异检测使用CNVkit、FACETS和Sequenza等工具。CNVkit通过比较肿瘤和正常样本的读段深度来检测拷贝数变异，同时利用off-target区域的信息提高分辨率。FACETS使用等位基因频率和读段深度联合分析来推断拷贝数和纯度。Sequenza从测序数据估计肿瘤纯度和倍性，然后检测拷贝数变异。体细胞结构变异检测使用Manta、Delly和GRIDSS等工具，这些工具可以检测肿瘤中的缺失、重复、倒位和易位等结构变异。亚克隆变异检测与克隆结构推断使用PyClone和SciClone等工具，它们通过聚类变异的等位基因频率来推断肿瘤的克隆结构，揭示肿瘤的异质性。肿瘤突变负荷（Tumor Mutational Burden，TMB）计算是免疫治疗疗效预测的重要指标，定义为每Mb基因组中体细胞突变的数量。突变特征分析使用SigProfiler和Mutalisk等工具，通过非负矩阵分解从突变谱中提取突变过程特征，不同的突变特征对应不同的致癌因素。当你进行肿瘤基因组分析时，体细胞变异检测是核心步骤，需要根据样本类型和分析目标选择合适的检测工具和参数。 .. code-block:: bash # Mutect2体细胞变异检测 gatk Mutect2 \ -R reference.fasta \ -I tumor.bam -tumor tumor_sample \ -I normal.bam -normal normal_sample \ --germline-resource af-only-gnomad.vcf.gz \ --panel-of-normals pon.vcf.gz \ -O somatic_unfiltered.vcf.gz gatk FilterMutectCalls \ -V somatic_unfiltered.vcf.gz \ -R reference.fasta \ -O somatic_filtered.vcf.gz # Strelka2体细胞变异检测 configureStrelkaSomaticWorkflow.py \ --normalBam normal.bam --tumorBam tumor.bam \ --referenceFasta reference.fasta \ --runDir strelka_run strelka_run/runWorkflow.py -m local -j 16 # CNVkit体细胞拷贝数变异 cnvkit.py batch tumor.bam -n normal.bam \ -f reference.fasta --output-reference my_reference.cnn \ --output-dir cnvkit_output --diagram --scatter 结构变异检测 ~~~~~~~~~~~~ 结构变异（SV）是基因组中大于50bp的序列改变，包括缺失、重复、倒位、易位和插入等类型。结构变异对基因功能的影响通常比SNV更大，可以导致基因融合、基因破坏或基因剂量改变，与许多遗传疾病和肿瘤密切相关。结构变异的检测策略主要分为四类。 基于读长对（read-pair）的SV检测方法利用双端读段之间的距离和方向异常来推断结构变异。当双端读段的比对距离显著大于或小于预期插入片段大小时，提示可能存在缺失或插入；当读段对的方向异常（如反向比对）时，提示可能存在倒位。Lumpy、Manta和Delly是常用的基于读长对的SV检测工具。基于局部组装的SV检测方法将异常读段及其邻近读段提取出来进行局部组装，从组装结果中直接识别结构变异的断点序列。Manta、SvABA和DELLY2支持局部组装策略，这种方法可以提供更精确的断点位置和插入序列信息。基于读长深度（read-depth）的SV检测方法通过比较不同基因组区域的读段覆盖深度来检测拷贝数变异。CNVnator使用均值-方差分割算法检测读段深度的变化，Control-FREEC通过滑动窗口方法检测拷贝数变异。基于双端读长的片段大小分布方法通过分析读段对的插入片段大小分布来识别异常，读段对方向异常可以识别倒位和易位。 基于长读长的SV检测利用长读长可以跨越结构变异断点的优势，直接从单条读段中检测结构变异。Sniffles2是长读长SV检测的标准工具，CuteSV使用聚类策略检测SV，SVision可以检测复杂的嵌套SV。SV断点精确解析与序列内容推断是SV检测的高级分析，通过局部组装或长读长比对可以精确确定断点位置和插入序列的来源。插入来源分析调查未知插入序列的来源，判断它是来自转座子插入、病毒整合还是其他机制。当你需要全面检测基因组中的结构变异时，结合多种检测策略可以提高灵敏度，特别是对于复杂的嵌套结构变异。 .. code-block:: bash # Manta结构变异检测 configManta.py --bam aligned.bam --referenceFasta reference.fasta \ --runDir manta_run manta_run/runWorkflow.py -m local -j 16 # Sniffles2长读长SV检测 sniffles --input aligned_ont.bam \ --vcf sv_calls.vcf \ --reference reference.fasta \ --threads 16 # CNVnator读段深度CNV检测 cnvnator -root sample.root -genome reference.fasta -tree aligned.bam cnvnator -root sample.root -his 1000 -genome reference.fasta cnvnator -root sample.root -stat 1000 cnvnator -root sample.root -partition 1000 cnvnator -root sample.root -call 1000 > cnv_calls.txt 变异过滤与质控 ~~~~~~~~~~~~~~ 变异过滤是去除假阳性变异、保留真实变异的关键步骤。基于质量分数的过滤使用VCF文件中的QUAL、GQ、DP和AD等字段来过滤低质量变异。QUAL是变异位点的Phred质量分数，GQ是基因型质量分数，DP是测序深度，AD是各等位基因的深度。基于硬过滤（hard filtering）使用预设的阈值对变异进行过滤，例如QUAL<30或DP<10的变异被过滤掉。GATK最佳实践推荐使用VQSR（Variant Quality Score Recalibration）代替硬过滤，VQSR利用变异的多个质量指标构建机器学习模型来区分真阳性和假阳性变异。然而，VQSR需要大量已知变异作为训练集，对于非人类物种可能不适用。GATK最佳实践过滤指标包括QD（质量深度比，过滤QD<2）、FS（Fisher精确检验，过滤FS>60）、SOR（链偏倚比，过滤SOR>3）、MQRankSum（比对质量秩和检验，过滤MQRankSum<-12.5）和ReadPosRankSum（读段位置秩和检验，过滤ReadPosRankSum<-8）。 群体分层去除常见变异是将变异与公共数据库中的已知变异进行比较，去除在健康人群中频率较高的变异，这些变异不太可能是致病变异。常用的参考数据库包括gnomAD（基因组聚合数据库，包含超过12万个人的外显子组和基因组数据）和1000 Genomes Project（千人基因组计划）。变异位置重叠与注释过滤根据变异的功能性影响预测来过滤，例如只保留影响蛋白质功能的变异。变异可视化核查使用IGV查看变异位点的比对情况，人工确认变异的真实性。当你进行罕见病基因诊断时，严格的变异过滤是缩小候选变异范围的关键步骤，通常需要结合多种过滤策略逐步缩小候选列表。 .. code-block:: bash # GATK硬过滤 gatk VariantFiltration \ -V variants.vcf.gz \ --filter-expression "QD < 2.0" --filter-name "QD2" \ --filter-expression "FS > 60.0" --filter-name "FS60" \ --filter-expression "MQ < 40.0" --filter-name "MQ40" \ --filter-expression "SOR > 3.0" --filter-name "SOR3" \ --filter-expression "MQRankSum < -12.5" --filter-name "MQRankSum-12.5" \ --filter-expression "ReadPosRankSum < -8.0" --filter-name "ReadPosRankSum-8" \ -O filtered_variants.vcf.gz # GATK VQSR变异质量分数校准 gatk VariantRecalibrator \ -R reference.fasta -V variants.vcf.gz \ --resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap.vcf.gz \ --resource:omni,known=false,training=true,truth=false,prior=12.0 omni.vcf.gz \ --resource:1000G,known=false,training=true,truth=false,prior=10.0 1000G.vcf.gz \ --resource:dbsnp,known=true,training=false,truth=false,prior=2.0 dbsnp.vcf.gz \ -mode SNP -O output.recal --tranches-file output.tranches # bcftools过滤 bcftools filter -e 'QUAL<30 || DP<10 || MQ<40' -Oz -o filtered.vcf.gz variants.vcf.gz 变异注释与优先级排序 ~~~~~~~~~~~~~~~~~~~~ 变异注释是为检测到的变异添加功能信息和临床信息的过程，变异优先级排序是根据注释信息将最可能致病的变异排在前面。这是从数百万变异中找到致病变异的关键步骤。 功能性影响预测使用SnpEff、ANNOVAR和VEP（Variant Effect Predictor）等工具来预测变异对基因功能的影响。SnpEff根据基因注释预测变异的效应类型（如同义突变、错义突变、无义突变、剪接位点突变等）。ANNOVAR是一个灵活的变异注释框架，支持多种数据库的注释。VEP是Ensembl提供的变异注释工具，功能全面且持续更新。保守性评分评估变异位点在进化中的保守程度，高度保守的位点发生变异更可能有害。PhyloP和PhastCons是基于多物种比对的保守性评分，GERP++通过检测进化约束来评估保守性。剪接影响预测使用SpliceAI和dbscSNV来预测变异对mRNA剪接的影响，SpliceAI使用深度学习模型预测剪接位点的创建或破坏，这对于远离经典剪接位点但影响剪接的变异尤为重要。蛋白功能影响预测使用SIFT（预测是否影响蛋白质功能）、PolyPhen-2（预测是否可能有害）、CADD（综合多种特征的评分）和REVEL（集成多个预测工具的结果）来评估错义变异的致病性。 致病性数据库提供已知致病变异的信息。ClinVar是公共的致病变异数据库，记录了变异与疾病的关系和临床意义分类。HGMD（Human Gene Mutation Database）是商业数据库，收录了大量已发表的致病变异。OMIM（Online Mendelian Inheritance in Man）提供基因与遗传病的关系信息。DisGeNET整合了基因与疾病的关联数据。人群频率查询使用gnomAD、1000 Genomes、ExAC和ALFA等数据库来获取变异在人群中的频率信息，频率极低的变异更可能是致病变异。调控区域注释使用ENCODE和Roadmap Epigenomics的数据来注释位于非编码区的变异，JASPAR数据库提供转录因子结合位点的信息。药物基因组学注释使用PharmGKB数据库来注释影响药物反应的变异。结果输出与过滤脚本结合ANNOVAR和tabix等工具，可以根据注释信息对变异进行灵活的过滤和筛选。当你进行临床基因诊断时，变异注释和优先级排序是从海量变异中找到致病原因的核心步骤，需要综合多种证据来源进行判断。 .. code-block:: bash # SnpEff变异功能注释 snpEff -v GRCh38 variants.vcf.gz > annotated_snpeff.vcf # VEP变异注释 vep -i variants.vcf.gz -o annotated_vep.vcf \ --cache --assembly GRCh38 --offline \ --everything --fork 4 # ANNOVAR变异注释 table_annovar.pl variants.vcf humandb/ \ -buildver hg38 -out annotated \ -remove -protocol refGene,clinvar_20240101,gnomad211_exome,dbnsfp42a \ -operation g,f,f,f -nastring . -vcfinput # ClinVar致病性查询 .. code-block:: python import subprocess def query_clinvar(variant_list): results = [] for var in variant_list: chrom, pos, ref, alt = var cmd = f"tabix clinvar.vcf.gz {chrom}:{pos}-{pos}" output = subprocess.run(cmd, shell=True, capture_output=True, text=True) for line in output.stdout.strip().split("\n"): if line and not line.startswith("#"): fields = line.split("\t") if fields[3] == ref and fields[4] == alt: info = dict(item.split("=") for item in fields[7].split(";") if "=" in item) results.append({ "variant": f"{chrom}:{pos}{ref}>{alt}", "clinical_significance": info.get("CLNSIG", "unknown"), "review_status": info.get("CLNREVSTAT", "unknown") }) return results 比较基因组学 ------------ 比较基因组学是通过比较不同物种或不同个体之间的基因组序列来揭示进化关系、功能约束和基因组结构变化的学科。比较基因组学的研究范围从两个物种的全基因组比对到数百个物种的泛基因组分析，提供了理解基因组进化和功能的重要视角。 基因组比对 ~~~~~~~~~~ 全基因组比对是比较基因组学的基础操作，它将两个或多个基因组序列进行全局比对，识别同源区域和差异区域。MUMmer是最常用的全基因组比对工具集，它使用后缀树算法快速寻找最大精确匹配（Maximal Exact Match，MUM），然后通过LIS（最长递增子序列）算法将匹配排序，最后用Smith-Waterman算法填充间隙。MUMmer适用于近缘物种之间的基因组比对，速度非常快。LASTZ是另一个全基因组比对工具，它使用种子-扩展策略，适合远缘物种之间的比对。minimap2在组装模式下可以将组装的contig比对到参考基因组，用于组装评估和基因组比较。Cactus和progressiveCactus是基于图的全基因组比对工具，它们可以同时比对多个基因组，生成多基因组比对结果，这种方法比两两比对更能保留多物种共享的同源关系。 比对后过滤与处理是全基因组比对的重要步骤。MUMmer的delta-filter可以根据比对长度和一致性过滤低质量比对。Chain/Net格式是UCSC基因组浏览器使用的比对结果格式，Chain记录局部比对的片段，Net将Chain组织成层次结构，区分一级同源（orthology）和二级同源（paralogy）。SyMAP是一个基因组共线性可视化和分析工具，可以交互式地查看基因组之间的共线性关系。微共线性分析关注局部区域的同源关系，MCScanX是最常用的微共线性分析工具，它可以检测共线性区块、计算共线性基因对和进行共线性分析。SynMap是CoGe平台上的在线共线性分析工具，DAGchainer通过动态规划算法检测共线性基因对。当你比较两个近缘物种的基因组时，MUMmer加MCScanX是一个经典的分析组合。 .. code-block:: bash # MUMmer全基因组比对 nucmer --maxmatch -p genome_compare reference.fasta query.fasta delta-filter -m -i 90 -l 1000 genome_compare.delta > filtered.delta show-coords -THrd filtered.delta > coords.txt mummerplot --png -p dotplot filtered.delta # MCScanX微共线性分析 # 首先准备BLASTP结果和GFF文件 MCScanX mcscanx_input # Cactus多基因组比对 cactus jobStore seqfile.txt output.maf --binariesMode local 核心与泛基因组分析 ~~~~~~~~~~~~~~~~~~ 泛基因组（pan-genome）是一个物种中所有基因的集合，包括核心基因组（core genome，所有个体共有的基因）、辅助基因组（accessory genome，部分个体拥有的基因）和特异基因簇（unique genes，仅个别个体拥有的基因）。泛基因组构建使用Panaroo、Roary、OrthoFinder和PGAP等工具。Panaroo在去噪和纠错方面表现优秀，适合高质量基因组的泛基因组构建。Roary是最广泛使用的大规模细菌泛基因组分析工具，它基于CD-HIT的快速聚类算法。OrthoFinder不仅可以构建泛基因组，还可以推断基因树的系统发育关系。PGAP是NCBI的参考泛基因组分析流程。基因存在-缺失矩阵是泛基因组分析的核心输出，它记录了每个基因在每个样本中的存在或缺失状态，这个矩阵可以用于下游的群体遗传学分析和关联研究。泛基因组曲线拟合通过Heaps定律拟合泛基因组大小与样本数量的关系，判断泛基因组是开放型（随样本增加持续增长）还是闭合型（随样本增加趋于稳定）。物种内多样性评估通过计算核心基因的SNP距离和pangenome的基因内容距离来评估物种内的遗传多样性。当你研究一个物种的泛基因组时，样本的多样性和数量是关键因素，样本越多样化，发现的辅助基因越多。 .. code-block:: bash # Roary泛基因组构建 roary -e -n -p 16 -i 95 -cd 99 *.gff # Panaroo泛基因组构建（更精确） panaroo -i *.gff -o panaroo_output --clean-mode strict -t 16 # OrthoFinder直系同源分析 orthofinder -f proteins_dir/ -t 16 -a 8 # 泛基因组曲线拟合 .. code-block:: python import numpy as np from scipy.optimize import curve_fit def heaps_law(n, kappa, alpha): return kappa * n ** alpha sample_sizes = np.array([5, 10, 20, 30, 50, 80, 100]) pan_genes = np.array([4500, 5200, 6000, 6500, 7200, 7800, 8100]) popt, pcov = curve_fit(heaps_law, sample_sizes, pan_genes) kappa, alpha = popt print(f"Heaps定律参数: kappa={kappa:.1f}, alpha={alpha:.3f}") if alpha < 1: print("泛基因组类型: 闭合型 (closed)") else: print("泛基因组类型: 开放型 (open)") 基因组重排与结构进化 ~~~~~~~~~~~~~~~~~~~~ 基因组重排是指基因组中大尺度结构的变化，包括倒位、易位、重复和缺失等。同源区块（synteny blocks）检测是识别基因组重排的基础，同源区块是指两个基因组中基因顺序和方向保守的区域。共线性点图（dotplot）是最直观的基因组比较可视化方法，MUMmer的mummerplot工具可以生成点图，横轴和纵轴分别表示两个基因组的位置，点表示同源区域，对角线表示共线性，偏离对角线的模式表示重排事件。倒位在点图上表现为对角线的镜像，易位表现为远离对角线的同源片段，重复表现为同一位置的多个同源片段。基因组进化速率估计通过计算同源区块之间的替换率和重排率来估计基因组的进化速度。基因组大小演化探讨基因组大小变化的机制，包括全基因组复制、转座子扩增和大规模缺失等。当你比较不同物种的基因组结构时，共线性分析可以揭示基因组在进化过程中发生的重排事件，帮助理解物种分化的基因组基础。 直系同源与旁系同源 ~~~~~~~~~~~~~~~~~~ 直系同源（ortholog）和旁系同源（paralog）是理解基因进化关系的两个核心概念。直系同源基因是由物种分化事件产生的同源基因，通常保持相同的功能；旁系同源基因是由基因复制事件产生的同源基因，功能可能发生分化。直系同源基因群集使用OrthoFinder、OrthoMCL和eggNOG等工具。OrthoFinder使用基因树与物种树协调的方法来推断直系同源关系，是目前最准确的直系同源推断工具之一。OrthoMCL使用MCL聚类算法对BLAST结果进行聚类，速度较快但精度略低。eggNOG是一个预计算的直系同源群数据库，覆盖了大量物种。 基因家族扩张与收缩分析使用CAFE和BADGER等工具，CAFE通过随机生灭模型来推断基因家族的扩张和收缩事件，判断哪些基因家族在特定进化分支上发生了显著的大小变化。基因复制（Whole Genome Duplication，WGD）事件检测通过Ks（同义替换率）分布分析来识别，WGD事件会在Ks分布中产生一个峰值，因为同时复制的基因对具有相似的进化时间。基因组内共线性分析也可以检测WGD事件，大规模的区块重复是WGD的证据。物种特异性基因识别是找出仅在某个物种中存在的基因，这些基因可能与该物种的特有性状相关。系统发育树与基因树一致性检验使用ASTRAL等工具来检测基因树与物种树之间的不一致，不一致可能由不完全谱系分选（Incomplete Lineage Sorting，ILS）、基因流或基因复制/丢失等原因造成。当你研究基因的进化关系时，区分直系同源和旁系同源是非常重要的，因为功能推断通常基于直系同源关系。 .. code-block:: bash # OrthoFinder直系同源群集 orthofinder -f proteins_dir/ -t 16 -a 8 -S diamond # CAFE基因家族扩张与收缩分析 cafe5 -i gene_families.txt -t species_tree.nwk -p # Ks分布分析检测WGD事件 .. code-block:: python import numpy as np from scipy.stats import gaussian_kde ks_values = np.loadtxt("ks_pairs.tsv", usecols=2) ks_filtered = ks_values[(ks_values > 0.01) & (ks_values < 5)] kde = gaussian_kde(ks_filtered) x_grid = np.linspace(0, 5, 500) density = kde(x_grid) peaks_x = x_grid[np.argsort(density)[-3:]] print(f"Ks分布峰值位置: {sorted(peaks_x)}") print("如果存在明显的第二个峰，可能表明发生过全基因组复制事件") 比较基因组学可视化 ~~~~~~~~~~~~~~~~~~ 比较基因组学的结果需要通过有效的可视化来展示和解读。共线性图是最常用的比较基因组学可视化方式，Circos可以生成环形基因组图，在图上展示共线性关系、变异分布和注释信息等多层数据，是发表高质量论文的标准工具。Synteny Viewer和Rideogram是R语言中的共线性可视化工具，适合在R环境中进行交互式分析。基因组浏览器比较轨道使用UCSC、Ensembl和JBrowse等浏览器来并排展示多个物种的基因组注释，通过比较轨道可以直观地看到同源区域的注释差异。插入/缺失区域标注在比较基因组中标记物种间存在或缺失的区域，这些区域可能包含物种特异的适应性基因。系统发育基因组学（Phylogenomics）利用基因组规模的数据构建物种树，使用超矩阵（concatenation）或超树（coalescence）方法整合大量基因的系统发育信号。当你需要展示比较基因组学的结果时，Circos图和共线性点图是最常用的可视化方式，它们可以直观地展示基因组之间的结构关系。 .. code-block:: bash # Circos共线性图配置 # 需要准备karyotype, links, highlights等配置文件 circos -conf circos.conf # 使用R绘制共线性图 .. code-block:: r library(tidyr) library(ggplot2) synteny_data <- data.frame( species1_chr = c("chr1", "chr1", "chr2"), species1_start = c(1000, 50000, 2000), species1_end = c(40000, 90000, 30000), species2_chr = c("chrA", "chrA", "chrB"), species2_start = c(2000, 60000, 1000), species2_end = c(42000, 100000, 29000) ) ggplot(synteny_data) + geom_segment(aes(x = species1_start, xend = species2_start, y = 1, yend = 2, color = species1_chr), size = 1.5, alpha = 0.7) + scale_color_brewer(palette = "Set1") + theme_minimal() + labs(title = "基因组共线性图", x = "基因组位置", y = "") 泛基因组学 ---------- 泛基因组学是研究一个物种内所有基因组变异的学科，它超越了单一参考基因组的局限，全面描述物种的基因组多样性。传统的线性参考基因组只代表一个或少数个体的基因组序列，无法涵盖物种内的全部遗传变异，这导致了参考基因组偏倚问题。泛基因组学通过构建包含物种内所有已知序列的泛基因组来解决这个问题，为更准确的比对和变异检测提供基础。 泛基因组概念 ~~~~~~~~~~~~ 泛基因组（pan-genome）的概念最早由Tettelin等人在2007年研究链球菌基因组时提出，它被定义为一个物种中所有基因序列的总和。泛基因组可以分为开放型和闭合型两种。开放型泛基因组随着新个体的加入持续增长，意味着每次测序新个体都可能发现新的基因，这通常发生在高度多样化的物种中，如许多细菌物种。闭合型泛基因组在加入一定数量的个体后趋于稳定，新个体不再贡献新的基因，这通常发生在遗传多样性较低的物种中。判断泛基因组是开放型还是闭合型，可以通过Heaps定律拟合泛基因组曲线来实现，当Heaps指数小于1时为开放型，等于1时为闭合型。 线性参考基因组偏倚是传统基因组分析的根本性问题。当你将测序数据比对到线性参考基因组时，样本中存在但参考基因组中不存在的序列无法被正确比对，这些序列被称为非参考等位基因或存在-缺失变异（Presence/Absence Variation，PAV）。线性参考基因组偏倚会导致变异检测的假阴性和等位基因频率的估计偏差，特别是对于非参考等位基因频率较高的群体。图泛基因组（graph pan-genome）是解决参考偏倚的新方法，它将所有已知的基因组变异整合到一个图结构中，图中的节点表示序列片段，边表示序列之间的连接关系。 图泛基因组构建使用minigraph、vg和PanGenome Graph Builder等工具。minigraph是一个轻量级的图泛基因组构建工具，它可以快速将多个基因组组装整合为一个图结构，适合大规模泛基因组的初步构建。vg（Variation Graph）是一个功能全面的图基因组工具包，支持图构建、图比对和图上的变异检测。GraphAligner将长读长比对到图泛基因组上，vg giraffe是vg中专门为短读长设计的图比对工具。泛基因组变异调用在图上进行变异检测，可以同时检测SNV、Indel和结构变异，减少参考偏倚。泛基因组注释将基因和功能元件的注释映射到图结构上，处理基因的存在-缺失变异。存在-缺失变异（PAV）分析统计每个基因在每个个体中的存在或缺失状态，分析PAV与表型的关联。泛基因组规模的GWAS将PAV和SNV与传统GWAS方法结合，在泛基因组尺度上寻找与性状关联的遗传变异。当你研究一个遗传多样性丰富的物种时，泛基因组分析可以揭示线性参考基因组无法捕获的遗传变异，为育种和群体遗传学研究提供更全面的信息。 .. code-block:: bash # minigraph图泛基因组构建 minigraph -c -xggs reference.fasta assemblies.txt > graph.gfa # vg图泛基因组构建 vg construct -r reference.fasta -v variants.vcf.gz > graph.vg vg index -x graph.xg -g graph.gcsa graph.vg # vg giraffe短读长图比对 vg giraffe -x graph.xg -g graph.gcsa -H graph.min \ -d graph.dist -f reads_R1.fq.gz -f reads_R2.fq.gz \ > aligned.gam # GraphAligner长读长图比对 GraphAligner -g graph.gfa -f long_reads.fasta \ -a aligned.gaf -t 16 重复序列分析与转座子生物学 -------------------------- 重复序列是基因组中最丰富的功能元件，转座子是重复序列的主要组成部分。转座子不仅是基因组的寄生序列，它们还通过插入突变、基因表达调控和基因组结构变异等方式对基因组进化和基因功能产生深远影响。重复序列分析是基因组注释的重要组成部分，转座子生物学则深入探讨转座子的活性、进化和功能。 转座子活性评估 ~~~~~~~~~~~~~~ 转座子活性评估是检测转座子在个体间插入位置多态性的过程。转座子插入多态性（Transposable Element Insertion Polymorphism，TE-IP）是指转座子在同一个体的不同个体中的插入位置不同，这种多态性可以像SNP一样作为遗传标记使用。TE插入多态性的检测方法主要有三种：基于短读长的检测通过识别读段对中一个读段比对到基因组而另一个读段比对到转座子序列的嵌合读段来发现新的TE插入；基于长读长的检测可以直接跨越TE插入位点和侧翼序列，更准确地定位插入位置；基于组装的检测通过比较不同个体的基因组组装来发现TE插入差异。当你研究转座子对表型变异的贡献时，TE插入多态性分析可以揭示转座子插入与基因表达变化或表型变异之间的关联。 .. code-block:: bash # 短读长TE插入多态性检测（TELR） telr call -b aligned.bam -c reference.fasta \ -t te_consensus.fasta -o telr_output # 长读长TE插入检测（TEbreak） tebreak -b aligned_ont.bam -r reference.fasta \ -t te_families.fasta -o tebreak_output # 基于组装的TE插入比较 # 比较不同个体组装中的TE注释差异 转座子家族进化 ~~~~~~~~~~~~~~ 转座子家族进化通过分子系统发育分析来研究转座子家族的进化历史和分类关系。转座子的系统发育树可以揭示转座子家族的扩张时间、分化关系和水平转移事件。构建转座子系统发育树通常使用转座子编码的转座酶或反转录酶序列，因为这些蛋白质序列包含了转座子进化的信息。不同转座子家族的进化速率不同，年轻的转座子家族通常具有更高的序列一致性和活性，而古老的转座子家族积累了较多突变，可能已经失去转座活性。当你分析一个新基因组中的转座子组成时，系统发育分析可以帮助确定转座子家族的分类归属和进化来源。 TE插入时间估计 ~~~~~~~~~~~~~~ TE插入时间估计是推断转座子插入基因组的时间，这对于理解转座子的进化动态非常重要。对于LTR反转录转座子，插入时间可以通过5'LTR和3'LTR之间的序列分歧度来估计。LTR反转录转座子在插入基因组时，5'LTR和3'LTR的序列是完全相同的，随着时间推移，两个LTR独立积累突变，序列分歧度逐渐增大。假设中性进化，插入时间T可以用公式T=K/(2r)来估计，其中K是两个LTR之间的替换距离，r是每个位点每年的替换速率。当你研究植物基因组中LTR反转录转座子的爆发历史时，LTR分歧度分析可以揭示转座子在不同进化时期的活性变化，通常会发现某些时期转座子活性显著增强，形成转座子爆发。 .. code-block:: python def estimate_ltr_age(k_distance, substitution_rate=1.3e-8): years = k_distance / (2 * substitution_rate) return years ltr_pairs = [ ("RLG1_chr1_10000", 0.02), ("RLG1_chr2_50000", 0.15), ("RLC2_chr3_80000", 0.001), ("RLG3_chr5_200000", 0.45), ] print(f"{'LTR名称':<25} {'K距离':<10} {'插入时间(百万年)':<15}") print("-" * 50) for name, k in ltr_pairs: age = estimate_ltr_age(k) print(f"{name:<25} {k:<10.3f} {age/1e6:<15.2f}") TE对基因表达的调控影响 ~~~~~~~~~~~~~~~~~~~~~~ 转座子对基因表达的调控影响是近年来备受关注的研究领域。转座子序列中包含大量的调控元件，如启动子、增强子、沉默子和转录因子结合位点，当转座子插入到基因附近时，这些调控元件可以影响邻近基因的表达。转座子对基因表达的影响可以是正面的（提供新的启动子或增强子）或负面的（干扰原有调控元件或导致表观遗传沉默）。在某些情况下，转座子来源的调控元件被宿主基因组驯化，成为基因调控网络的重要组成部分。例如，人类基因组中约25%的启动子区域含有转座子来源的序列。当你研究基因表达的调控机制时，考虑转座子的贡献可以揭示传统分析可能忽略的调控来源。 TE衍生的小非编码RNA ~~~~~~~~~~~~~~~~~~~~ 转座子可以产生多种小非编码RNA，包括siRNA（小干扰RNA）、piRNA（Piwi相互作用RNA）和miRNA（微小RNA）。这些小RNA在转座子沉默和基因组防御中发挥关键作用。piRNA是动物生殖细胞中特异性表达的小RNA，主要功能是沉默转座子以维持基因组稳定性。piRNA通过指导Piwi蛋白切割转座子mRNA和介导转座子位点的表观遗传修饰来实现转座子沉默。siRNA通过RNA干扰途径沉默双链RNA来源的转座子转录本。某些miRNA来源于转座子序列，它们可能同时调控转座子和宿主基因的表达。当你研究生殖细胞或早期胚胎中的转座子调控时，小RNA分析是理解转座子沉默机制的重要手段。 全基因组转座子图谱构建 ~~~~~~~~~~~~~~~~~~~~~~ 全基因组转座子图谱构建是系统描述基因组中所有转座子位置、类型和状态的过程。转座子图谱包括转座子的基因组坐标、分类归属、序列完整性（完整或片段化）、插入方向和侧翼序列等信息。构建转座子图谱通常使用RepeatMasker的注释结果作为基础，然后通过自定义脚本整合和格式化信息。对于非参考基因组的转座子注释，需要先运行RepeatModeler构建物种特异的重复序列库，然后使用RepeatMasker进行注释。转座子图谱可以用于多种下游分析，如转座子密度分布、转座子与基因的关联分析、转座子的年龄分布等。当你需要全面了解一个基因组中的转座子组成时，全基因组转座子图谱提供了系统性的视图。 转座子相关的结构变异 ~~~~~~~~~~~~~~~~~~~~ 转座子活动是基因组结构变异的重要来源。转座子的插入可以导致基因破坏或基因调控改变；不等交换（unequal crossing over）发生在同源转座子之间，可以导致缺失和重复；转座子介导的重组可以产生倒位和易位；非等位同源重组（Non-Allelic Homologous Recombination，NAHR）发生在同源但非等位的转座子拷贝之间，是基因组结构变异的主要机制之一。当你检测到基因组中的结构变异时，检查断点附近是否存在转座子序列可以帮助判断结构变异的产生机制，许多结构变异的断点位于转座子内部或附近，提示转座子介导的重组是其产生的原因。 .. code-block:: python def analyze_te_sv_association(sv_breakpoints, te_annotations, window=1000): associated = [] for sv_id, sv_chr, sv_pos in sv_breakpoints: for te_chr, te_start, te_end, te_type in te_annotations: if sv_chr == te_chr: if abs(sv_pos - te_start) <= window or abs(sv_pos - te_end) <= window: associated.append({ "sv_id": sv_id, "sv_pos": sv_pos, "te_type": te_type, "te_start": te_start, "te_end": te_end, "distance": min(abs(sv_pos - te_start), abs(sv_pos - te_end)) }) break print(f"SV总数: {len(sv_breakpoints)}") print(f"与TE关联的SV: {len(associated)} ({len(associated)/len(sv_breakpoints)*100:.1f}%)") te_types = {} for item in associated: t = item["te_type"] te_types[t] = te_types.get(t, 0) + 1 for t, count in sorted(te_types.items(), key=lambda x: -x[1]): print(f" {t}: {count}") return associated 单倍型分析 ---------- 单倍型（haplotype）是指同一条染色体上紧密连锁的一组等位基因的组合，它反映了亲本染色体的遗传信息。在二倍体生物中，每个个体拥有两套同源染色体，分别来自父本和母本。传统的变异检测通常不考虑变异之间的相位关系，即不知道哪些变异位于同一条染色体上。单倍型分型（phasing）就是确定变异之间的相位关系，将杂合变异分配到各自的亲本染色体上，重建单倍型序列的过程。单倍型信息对于理解基因组功能、疾病遗传模式和群体历史具有重要意义。 单倍型分型概念 ~~~~~~~~~~~~~~ 单倍型分型的核心目标是将个体的杂合变异分配到两条同源染色体上，确定哪些变异位于同一条染色体上。单倍型分型的方法可以分为两大类：基于实验数据的长读长分型和基于统计推断的分型。长读长分型利用PacBio HiFi或Nanopore等长读长测序技术，由于一条读段可以跨越多个变异位点，读段上的等位基因组合直接反映了单倍型信息。当一条长读段覆盖了多个杂合位点时，这些位点上的等位基因必然来自同一条染色体，因此可以直接确定它们的相位关系。长读长分型的准确性取决于读段长度和覆盖度，读段越长，能够跨越的变异位点越多，分型结果越完整。 统计推断分型不依赖长读长数据，而是利用群体中的连锁不平衡（Linkage Disequilibrium，LD）信息来推断变异之间的相位关系。统计分型的基本原理是：在群体中频繁共遗传的等位基因更可能位于同一条染色体上。统计分型需要一个参考面板（reference panel），即已经完成分型的大量样本的单倍型数据，通过将目标样本的基因型与参考面板中的单倍型进行比较，找到最可能的相位配置。统计分型的准确性取决于参考面板的大小和多样性、标记密度以及目标样本与参考面板的遗传距离。当你只有短读长测序数据时，统计推断分型是获得单倍型信息的主要方法，但它的分型长度受限于连锁不平衡的衰减距离，通常只能获得几十kb到几百kb的相位块。 .. code-block:: bash # WhatsHap长读长分型 whatshap phase --reference reference.fasta \ --output phased.vcf.gz \ variants.vcf.gz aligned_longread.bam # SHAPEIT统计分型 shapeit5 --input variants.vcf.gz \ --reference reference.fasta \ --map genetic_map.txt \ --region chr1:1-250000000 \ --output phased_shapeit.vcf.gz \ --thread 16 # Eagle统计分型 eagle --vcf=variants.vcf.gz \ --geneticMapFile=genetic_map.txt \ --outPrefix=eagle_phased \ --numThreads 16 分型工具 ~~~~~~~~ 单倍型分型工具有多种选择，适用于不同的数据类型和分析需求。WhatsHap是一个专门利用长读长数据进行分型的工具，它采用分阶段动态规划算法，可以高效地利用覆盖多个变异位点的长读段来确定相位关系。WhatsHap特别适合PacBio HiFi和Nanopore数据的分型，它还可以整合家族信息（trio phasing）来提高分型准确性。HapCut2是另一个长读长分型工具，它使用最大切割图算法来寻找最优的相位配置，支持多种测序平台的数据。SHAPEIT是最广泛使用的统计分型工具，它通过隐马尔可夫模型将目标样本的基因型与参考面板中的单倍型进行匹配，推断最可能的相位配置。SHAPEIT5是最新版本，支持大规模群体数据的快速分型。Eagle是另一个统计分型工具，它使用基于数据压缩的快速算法，在保持高准确性的同时大幅提高了分型速度。当你需要对新测序样本进行分型时，如果有长读长数据，优先使用WhatsHap或HapCut2；如果只有短读长数据，SHAPEIT或Eagle配合适当的参考面板是标准选择。 单体型的变异调用与评估 ~~~~~~~~~~~~~~~~~~~~~~ 单体型水平的变异调用是在分型后进行变异检测，考虑变异的相位信息来提高检测准确性。传统的变异检测将每个位点独立处理，忽略了位点之间的关联，而单体型水平的检测利用了相邻位点的相位信息，可以更准确地估计等位基因频率和基因型质量。分型质量的评估通常使用switch error rate和hamming error rate两个指标。Switch error rate是指在分型结果中，相邻两个位点之间相位翻转的比例，它衡量分型的局部准确性。Hamming error rate是指分型结果中错误分配的等位基因占总等位基因数的比例，它衡量分型的全局准确性。N50分型长度是另一个重要指标，表示将所有相位块按长度排序后，累计长度达到总长度50%时的相位块长度，它衡量分型的连续性。当你评估分型结果时，switch error rate应该低于1%，N50分型长度越长越好，理想情况下应该达到染色体级别。 .. code-block:: python def evaluate_phasing(phased_vcf, truth_vcf): switch_errors = 0 total_pairs = 0 for chrom in phased_vcf: prev_phase = None for pos, (ref_allele, alt_allele, phase) in phased_vcf[chrom]: if prev_phase is not None: total_pairs += 1 if phase != prev_phase: switch_errors += 1 prev_phase = phase switch_error_rate = switch_errors / total_pairs if total_pairs > 0 else 0 print(f"Switch error rate: {switch_error_rate:.4f} ({switch_error_rate*100:.2f}%)") print(f"总相位对数: {total_pairs}") print(f"Switch错误数: {switch_errors}") return switch_error_rate 单倍型解析组装 ~~~~~~~~~~~~~~ 单倍型解析组装（haplotype-resolved assembly）是直接组装出两套单倍型序列，而不是组装一个共识序列。传统的基因组组装通常产生一个共识序列，将两条同源染色体的信息混合在一起，丢失了单倍型信息。单倍型解析组装利用Hi-C和长读长数据的组合来实现分相组装。Hi-C数据提供了染色体级别的相位信息，因为同一单倍型上的序列在三维空间中更倾向于相互靠近。长读长数据提供了局部的相位信息，因为一条读段上的变异必然来自同一条染色体。DipAsm是一个利用Hi-C和长读长数据进行单倍型解析组装的工具，它首先使用长读长进行contig组装，然后利用Hi-C数据将contig分配到两个单倍型，最后进行scaffold构建。HiCanu是Canu的HiFi模式，它可以在组装过程中区分两个单倍型，产生分相的组装结果。当你组装一个高杂合度的基因组时，单倍型解析组装可以避免等位基因被错误地组装为独立的contig，产生更准确的基因组表示。 .. code-block:: bash # DipAsm单倍型解析组装 dipasm -t 32 \ --hifi pacbio_hifi.fasta.gz \ --hic hic_R1.fq.gz,hic_R2.fq.gz \ -o dipasm_output # HiCanu分相组装 canu -p asm -d hicaru_output \ genomeSize=3g \ -pacbio-hifi pacbio_hifi.fasta.gz \ maxInputCoverage=200 单倍型特异性表达 ~~~~~~~~~~~~~~~~ 单倍型特异性表达（allele-specific expression，ASE）是指同源染色体上等位基因表达水平的差异。在二倍体生物中，大多数基因的两个等位基因表达水平相似，但某些基因存在表达偏向，一个等位基因的表达量显著高于另一个。ASE的检测需要同时知道基因型和表达量信息，首先确定个体的杂合位点，然后统计RNA-seq数据中每个等位基因对应的读段数量，使用统计检验（如二项检验或beta-binomial检验）判断是否存在显著的等位基因偏向。ASE的生物学意义包括基因组印记（imprinting，由亲本来源决定的等位基因表达差异）、顺式调控变异（影响邻近基因表达的变异）和X染色体失活等。当你研究基因表达的调控机制时，ASE分析可以揭示传统差异表达分析无法检测的等位基因水平的表达差异，特别是对于杂合致病变异的功能影响评估。 .. code-block:: python from scipy.stats import binom_test def test_ase(ref_counts, alt_counts, min_count=10, p_threshold=0.05): total = ref_counts + alt_counts if total < min_count: return None p_value = binom_test(ref_counts, total, 0.5) ref_ratio = ref_counts / total is_ase = p_value < p_threshold and abs(ref_ratio - 0.5) > 0.2 return { "ref_count": ref_counts, "alt_count": alt_counts, "ref_ratio": ref_ratio, "p_value": p_value, "is_ase": is_ase } sites = [(85, 15), (55, 45), (120, 30), (50, 50)] for ref, alt in sites: result = test_ase(ref, alt) if result: status = "ASE" if result["is_ase"] else "非ASE" print(f"Ref:{result['ref_count']} Alt:{result['alt_count']} " f"Ref比例:{result['ref_ratio']:.2f} p={result['p_value']:.4f} {status}") 单倍型网络分析 ~~~~~~~~~~~~~~ 单倍型网络分析是研究群体中单倍型之间进化关系的方法，与系统发育树不同，单倍型网络允许节点之间存在多条路径，更适合表示群体水平的进化关系，特别是当存在重组、平行突变或祖先单倍型仍然存在时。单倍型网络的构建方法包括中位数连接网络（median-joining network）、TCS网络和最小生成树网络等。中位数连接网络是最常用的方法，它通过引入假设的中间单倍型来连接观察到的单倍型，网络中的节点表示单倍型，边表示突变步骤。单倍型网络在群体遗传学和系统地理学中广泛应用，可以揭示群体的扩张历史、迁移路径和种群结构。当你研究一个基因或基因组区域在群体中的变异模式时，单倍型网络可以直观地展示不同单倍型之间的进化关系和频率分布。 .. code-block:: python def build_median_joining_network(haplotypes): from collections import Counter hap_counts = Counter(haplotypes) unique_haps = list(hap_counts.keys()) n = len(unique_haps) distances = {} for i in range(n): for j in range(i+1, n): dist = sum(a != b for a, b in zip(unique_haps[i], unique_haps[j])) distances[(i, j)] = dist print(f"单倍型数量: {n}") for (i, j), dist in sorted(distances.items(), key=lambda x: x[1]): print(f" H{i+1}({unique_haps[i][:20]}...) <-> H{j+1}({unique_haps[j][:20]}...) " f"距离={dist} 频率=({hap_counts[unique_haps[i]]},{hap_counts[unique_haps[j]]})") return distances haplotypes = [ "ATCGATCGATCGATCGATCG", "ATCGATCGATCGATCGATCG", "ATCGATCAATCGATCGATCG", "ATCGATCAATCGATCGATCG", "ATCGATCAATCGATCAATCG", "GTCGATCGATCGATCGATCG", ] build_median_joining_network(haplotypes) 基因组学分析流程与工作流 ------------------------ 基因组学分析通常涉及多个步骤和工具的组合，从原始测序数据到最终的生物学结论，需要经过数据质控、比对、变异检测、注释等多个环节。手动执行这些步骤不仅效率低下，而且容易出错，难以保证分析的可重复性。工作流管理系统通过将分析步骤定义为可执行的流程，实现了基因组学分析的自动化、标准化和可重复。 典型WGS分析流程 ~~~~~~~~~~~~~~~ 全基因组测序（WGS）的标准分析流程包括以下主要步骤：首先对原始FASTQ数据进行质量评估和质控修剪，去除低质量碱基和接头序列；然后将质控后的读段比对到参考基因组，生成BAM文件；对比对结果进行排序和索引，标记PCR重复；进行碱基质量分数校准（Base Quality Score Recalibration，BQSR）；使用校准后的BAM文件进行变异检测，生成VCF文件；最后对变异进行过滤和功能注释。这个流程的每一步都有明确的输入输出和质控标准，GATK最佳实践流程是目前最广泛使用的WGS分析标准。当你处理一批新的WGS数据时，遵循标准化的分析流程可以确保结果的可比性和可重复性，同时减少因分析差异导致的假阳性或假阴性结果。 .. code-block:: bash # 典型WGS分析流程（GATK最佳实践） # 步骤1: 质控 fastp -i sample_R1.fq.gz -I sample_R2.fq.gz \ -o clean_R1.fq.gz -O clean_R2.fq.gz # 步骤2: 比对 bwa mem -t 16 -M reference.fasta clean_R1.fq.gz clean_R2.fq.gz | \ samtools view -bS - > aligned.bam # 步骤3: 排序标记重复 samtools sort -@ 8 -o aligned_sorted.bam aligned.bam picard MarkDuplicates I=aligned_sorted.bam O=dedup.bam M=metrics.txt samtools index dedup.bam # 步骤4: BQSR gatk BaseRecalibrator -I dedup.bam -R reference.fasta \ --known-sites dbsnp.vcf.gz -O recal_data.table gatk ApplyBQSR -I dedup.bam -R reference.fasta \ --bqsr-recal-file recal_data.table -O recalibrated.bam # 步骤5: 变异检测 gatk HaplotypeCaller -R reference.fasta -I recalibrated.bam \ -O sample.g.vcf.gz -ERC GVCF GATK最佳实践 ~~~~~~~~~~~~ GATK最佳实践（Best Practices Pipeline）是Broad研究所制定的一套基因组学分析标准流程，它为种系变异检测和体细胞变异检测分别提供了详细的流程说明和参数推荐。种系变异检测流程包括：数据预处理（质控、比对、排序、标记重复、BQSR）、变异检测（HaplotypeCaller生成gVCF、联合基因分型）、变异过滤（VQSR或硬过滤）和功能注释。体细胞变异检测流程包括：数据预处理、Mutect2变异检测、FilterMutectCalls过滤、CalculateContamination污染估计和功能注释。GATK最佳实践不仅提供了工具和参数的选择，还详细说明了每个步骤的原理和质控标准。当你进行人类基因组的变异检测时，遵循GATK最佳实践可以确保分析结果的可靠性和与其他研究的一致性，这是发表高质量论文的基本要求。 Sentieon商业加速流程 ~~~~~~~~~~~~~~~~~~~~ Sentieon是一个商业化的基因组学分析加速软件，它重新实现了GATK最佳实践流程中的核心算法，在保持结果一致性的同时大幅提高了计算速度。Sentieon的DNAseq模块对应GATK的种系变异检测流程，TNseq模块对应体细胞变异检测流程。Sentieon使用高度优化的C++实现和并行计算技术，速度通常比GATK快5-10倍，内存占用也更低。当你需要处理大规模基因组数据（如数千个样本的群体测序项目）时，Sentieon可以显著缩短分析时间和降低计算成本，但需要注意它是商业软件，需要购买许可证。 工作流语言实现 ~~~~~~~~~~~~~~ 工作流语言是专门用于描述计算流程的领域特定语言，它们将分析步骤、输入输出和依赖关系以代码形式定义，使得流程可以自动化执行和版本控制。Snakemake是Python风格的工作流管理系统，它使用规则（rule）定义每个分析步骤，通过文件名模式匹配自动确定步骤之间的依赖关系，支持本地、集群和云平台的执行。Nextflow是另一个流行的工作流系统，它使用Groovy DSL定义流程，支持多种执行平台（本地、SLURM、AWS、Google Cloud等），具有强大的容器集成能力。CWL（Common Workflow Language）是一个开放标准，它使用YAML/JSON格式描述工作流，强调可移植性和互操作性。WDL（Workflow Definition Language）是Broad研究所开发的工作流语言，语法简洁，与Terra平台紧密集成。当你需要构建一个可重复的基因组学分析流程时，Snakemake和Nextflow是两个最流行的选择，Snakemake更适合Python用户，Nextflow更适合需要跨平台执行的场景。 .. code-block:: python # Snakemake工作流示例 (Snakefile) rule all: input: "results/variants/sample.vcf.gz" rule fastqc: input: r1="data/{sample}_R1.fq.gz", r2="data/{sample}_R2.fq.gz" output: html="qc/{sample}_fastqc.html" shell: "fastp -i {input.r1} -I {input.r2} " "-o clean/{wildcards.sample}_R1.fq.gz " "-O clean/{wildcards.sample}_R2.fq.gz " "--html {output.html}" rule align: input: r1="clean/{sample}_R1.fq.gz", r2="clean/{sample}_R2.fq.gz", ref="reference.fasta" output: bam="results/aligned/{sample}.bam" shell: "bwa mem -t 16 {input.ref} {input.r1} {input.r2} | " "samtools sort -@ 8 -o {output.bam}" rule variant_call: input: bam="results/aligned/{sample}.bam", ref="reference.fasta" output: vcf="results/variants/{sample}.vcf.gz" shell: "gatk HaplotypeCaller -R {input.ref} -I {input.bam} " "-O {output.vcf}" 容器化流程 ~~~~~~~~~~ 容器化技术是解决生物信息学软件环境依赖问题的有效方案。Docker是最广泛使用的容器化平台，它将软件及其所有依赖打包到一个容器镜像中，确保软件在任何环境中都能一致运行。在基因组学分析中，几乎所有常用工具都有Docker镜像可用，这意味着你不需要在本地安装和配置复杂的软件环境，只需拉取对应的Docker镜像即可运行分析。Singularity（现称Apptainer）是专为高性能计算环境设计的容器化平台，它不需要root权限，更适合在共享集群上使用。当你需要在不同的计算环境中部署基因组学分析流程时，使用容器化可以避免环境配置问题，确保分析结果的可重复性。工作流系统如Nextflow和Snakemake都支持直接调用Docker或Singularity容器，使得整个分析流程可以在容器化环境中运行。 .. code-block:: bash # Docker容器化运行BWA docker run --rm -v /data:/data biocontainers/bwa:v0.7.17 \ bwa mem -t 16 /data/reference.fasta /data/R1.fq.gz /data/R2.fq.gz # Singularity容器化运行GATK singularity exec --bind /data:/data docker://broadinstitute/gatk:4.4.0.0 \ gatk HaplotypeCaller -R /data/reference.fasta -I /data/aligned.bam \ -O /data/variants.vcf.gz # Nextflow中使用容器 # nextflow.config # docker.enabled = true # process.container = 'biocontainers/bwa:v0.7.17' 云平台执行 ~~~~~~~~~~ 云计算平台为基因组学分析提供了弹性可扩展的计算资源，特别适合需要大规模并行计算的项目。DNAnexus是一个专门面向基因组学的云平台，它提供了预配置的分析工具和数据管理功能，支持GATK最佳实践流程的一键执行。Seven Bridges是另一个基因组学云平台，它支持CWL工作流，提供了丰富的公共数据集和分析工具。Terra.bio是Broad研究所开发的云平台，原生支持WDL工作流，与GATK和Google Cloud紧密集成，提供了协作分析和数据共享功能。当你需要处理超大规模的基因组数据或需要快速扩展计算能力时，云平台可以按需提供计算资源，避免了本地集群的建设和维护成本。 流程性能调优 ~~~~~~~~~~~~ 基因组学分析流程的性能调优是提高分析效率的重要环节。内存优化需要根据每个工具的内存需求合理分配资源，例如BWA-MEM的内存主要由参考基因组的索引大小决定，GATK的内存需求与输入BAM文件的大小相关。CPU优化需要考虑工具的并行化能力，有些工具（如BWA-MEM、samtools sort）可以高效利用多线程，而有些工具（如某些R脚本）是单线程的。IO优化在处理大量文件时尤为重要，使用管道（pipe）将步骤串联可以避免中间文件的读写，使用压缩格式（如CRAM替代BAM）可以减少存储空间和IO时间。当你优化一个基因组学分析流程时，首先应该识别性能瓶颈，然后针对性地进行优化，通常比对和变异检测是最耗时的步骤，优化这两个步骤可以获得最大的性能提升。 .. code-block:: python import time import subprocess def benchmark_step(step_name, command, iterations=3): times = [] for i in range(iterations): start = time.time() subprocess.run(command, shell=True, capture_output=True) elapsed = time.time() - start times.append(elapsed) print(f" 迭代{i+1}: {elapsed:.1f}秒") avg = sum(times) / len(times) print(f"{step_name} 平均耗时: {avg:.1f}秒") return avg print("=== WGS流程性能基准测试 ===") t1 = benchmark_step("BWA比对", "bwa mem -t 16 ref.fasta R1.fq R2.fq > /dev/null") t2 = benchmark_step("排序", "samtools sort -@ 8 -o /dev/null aligned.bam") t3 = benchmark_step("变异检测", "gatk HaplotypeCaller -R ref.fasta -I aligned.bam -O /dev/null") total = t1 + t2 + t3 print(f"\n总耗时: {total:.1f}秒") print(f"BWA占比: {t1/total*100:.1f}% 排序占比: {t2/total*100:.1f}% 变异检测占比: {t3/total*100:.1f}%") 基因组学分析的伦理与数据安全 ---------------------------- 基因组数据不仅包含个体的遗传信息，还可能揭示其家族成员和群体的遗传特征。基因组学研究的伦理问题涉及隐私保护、知情同意、数据共享和公平使用等多个方面，理解这些伦理规范对于负责任地开展基因组学研究至关重要。 人与保护物种基因组数据的访问控制 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 人类基因组数据包含高度敏感的个人信息，需要严格的访问控制来保护参与者的隐私。受控访问（controlled access）是保护人类基因组数据的主要机制，研究者在获取数据前需要提交研究计划、通过伦理审查并签署数据使用协议。保护物种的基因组数据同样需要访问控制，因为某些濒危物种的位置信息可能被用于非法捕猎，基因组数据可能被用于生物盗猎。当你申请使用受控访问的基因组数据时，需要明确说明研究目的、数据使用范围和数据保护措施，确保数据仅用于经批准的研究目的。 .. code-block:: python def check_data_access_permission(user_id, dataset_id, permissions_db): required_level = permissions_db.get(dataset_id, {}).get("access_level", "open") user_level = permissions_db.get(user_id, {}).get("clearance", "none") approved_studies = permissions_db.get(user_id, {}).get("approved_studies", []) if required_level == "open": return True if required_level == "controlled": return (user_level in ["researcher", "admin"] and dataset_id in approved_studies) if required_level == "restricted": return user_level == "admin" return False permissions = { "user001": {"clearance": "researcher", "approved_studies": ["phs000001", "phs000002"]}, "user002": {"clearance": "admin", "approved_studies": ["all"]}, "phs000001": {"access_level": "controlled"}, "phs000002": {"access_level": "controlled"}, "phs000003": {"access_level": "restricted"}, } print(f"user001访问phs000001: {check_data_access_permission('user001', 'phs000001', permissions)}") print(f"user001访问phs000003: {check_data_access_permission('user001', 'phs000003', permissions)}") 去标识化方法 ~~~~~~~~~~~~ 去标识化（de-identification）是从基因组数据中移除或遮蔽可以直接识别个体身份的信息的过程。直接标识符包括姓名、身份证号、联系方式等，这些信息在数据共享前必须移除。然而，基因组数据本身具有高度可识别性，即使移除了所有直接标识符，基因组序列仍然可以唯一地识别一个人。准标识符如年龄、性别、地理来源等，虽然单独使用时不能识别个体，但组合使用时可能通过交叉比对重新识别个体。k-匿名化和l-多样性是常用的去标识化技术，k-匿名化确保每个准标识符组合至少对应k个个体，l-多样性确保每个等价类中敏感属性至少有l个不同的值。当你需要共享基因组数据时，去标识化是保护参与者隐私的第一步，但需要认识到基因组数据本身的可识别性使得完全匿名化几乎不可能。 .. code-block:: python def k_anonymize(data, quasi_identifiers, k=5): groups = {} for record in data: key = tuple(record[qi] for qi in quasi_identifiers) groups.setdefault(key, []).append(record) anonymized = [] suppressed = 0 for key, group in groups.items(): if len(group) >= k: anonymized.extend(group) else: suppressed += len(group) print(f"原始记录数: {len(data)}") print(f"k-匿名化后保留: {len(anonymized)}") print(f"抑制记录数: {suppressed} ({suppressed/len(data)*100:.1f}%)") return anonymized sample_data = [ {"id": "P001", "age": "30-40", "sex": "M", "region": "East", "diagnosis": "A"}, {"id": "P002", "age": "30-40", "sex": "M", "region": "East", "diagnosis": "B"}, {"id": "P003", "age": "30-40", "sex": "F", "region": "West", "diagnosis": "A"}, ] k_anonymize(sample_data, ["age", "sex", "region"], k=2) 二次使用同意规范 ~~~~~~~~~~~~~~~~ 二次使用同意（consent for secondary use）是指研究参与者同意其数据在原始研究目的之外被使用。基因组学研究的二次使用涉及复杂的伦理问题，因为基因组数据的价值往往在原始研究设计时无法完全预见。广泛同意（broad consent）允许数据在相关领域内被灵活使用，但可能超出参与者的预期。动态同意（dynamic consent）允许参与者在研究过程中随时更新其同意偏好，但增加了管理复杂度。分层同意（tiered consent）让参与者可以选择不同级别的数据共享范围。当你设计基因组学研究的知情同意流程时，需要平衡数据共享的科学价值与参与者自主权的尊重，确保参与者充分理解其数据可能的使用方式。 基因组数据库的受控访问 ~~~~~~~~~~~~~~~~~~~~~~ 基因组数据库的受控访问机制是保护敏感基因组数据的重要屏障。dbGaP（database of Genotypes and Phenotypes）是NIH建立的受控访问数据库，存储了大量人类基因组学和表型数据，研究者需要通过数据访问委员会（Data Access Committee，DAC）的审批才能获取数据。EGA（European Genome-phenome Archive）是欧洲的受控访问基因组数据库，与dbGaP类似，需要通过数据访问申请才能使用数据。这些数据库通常要求研究者签署数据使用协议，承诺仅将数据用于经批准的研究目的，不尝试重新识别参与者，并采取适当的数据安全措施。当你需要使用dbGaP或EGA中的数据时，申请过程可能需要数周时间，建议提前规划。 告知同意书与ICMJE要求 ~~~~~~~~~~~~~~~~~~~~~~ 告知同意书（Informed Consent Form）是基因组学研究中保障参与者权益的核心文件，它必须以参与者能够理解的语言清楚地说明研究目的、数据使用方式、潜在风险和参与者权利。ICMJE（International Committee of Medical Journal Editors）要求临床试验在发表前进行注册，并要求研究者披露利益冲突。对于基因组学研究，告知同意书还需要特别说明基因组数据的共享计划、数据保护措施和参与者撤回同意的权利。当你撰写基因组学研究的告知同意书时，建议参考所在机构伦理委员会的模板，并确保涵盖基因组数据特有的伦理考量。 族群特有的变异与隐私 ~~~~~~~~~~~~~~~~~~~~ 族群特有的变异是指在某些族群中频率较高而在其他族群中罕见的遗传变异。这些变异可以用于推断个体的族群归属，但同时也带来了隐私风险。即使个体数据已经去标识化，通过比对族群特有变异的频率分布，仍然可能推断出个体的族群来源。对于少数族群，这种推断可能特别敏感，因为少数族群的基因组数据较少，更容易被重新识别。此外，族群特有的致病变异可能导致对该族群的歧视或污名化。当你研究特定族群的基因组变异时，需要特别注意保护该族群的集体隐私，避免研究结果的误用。 基因组数据的匿名分享风险 ~~~~~~~~~~~~~~~~~~~~~~~~ 基因组数据的匿名分享是一个持续的伦理挑战。研究表明，即使不包含任何表型信息，仅凭基因组数据本身就可以通过多种方式进行重新识别。一种方法是利用公开的基因族谱数据库（如GEDmatch），通过搜索远亲的基因组数据来缩小候选范围。另一种方法是利用基因组数据中的表型预测信息（如眼睛颜色、头发颜色等外观特征）来辅助识别。随着基因组数据库的不断扩大和表型预测算法的改进，基因组数据的重新识别风险将持续增加。当你考虑分享基因组数据时，需要评估重新识别的风险，采取适当的保护措施，如限制数据的分辨率、使用合成数据或差分隐私技术。 .. code-block:: python import random def differential_privacy_count(true_count, epsilon=1.0): sensitivity = 1 scale = sensitivity / epsilon noise = random.gauss(0, scale) return max(0, round(true_count + noise)) true_counts = {"variant_A": 150, "variant_B": 3, "variant_C": 89} print("差分隐私保护示例 (epsilon=1.0):") print(f"{'变异':<15} {'真实计数':<10} {'隐私计数':<10}") print("-" * 35) for variant, count in true_counts.items(): private_count = differential_privacy_count(count, epsilon=1.0) print(f"{variant:<15} {count:<10} {private_count:<10}") print("\n注意: 稀有变异(variant_B)的计数变化比例更大，隐私保护更强")