转录组学分析 ================== 转录组学是研究细胞或组织中所有RNA分子的学科，它反映了基因在特定时间、特定条件下的表达状态。与基因组学研究DNA序列不同，转录组学关注的是基因表达的动态变化，这使得我们能够理解细胞如何响应环境刺激、发育过程和疾病状态。RNA测序（RNA-seq）技术是当前转录组学研究的主流方法，它通过高通量测序获得样本中所有RNA的序列信息，进而进行表达量定量、差异表达分析和功能注释等后续分析。本章将从实验设计开始，系统介绍RNA-seq数据的处理流程和分析方法。 5.1 RNA-seq 实验设计 -------------------- RNA-seq实验设计的合理性直接决定了后续数据分析的可靠性和结论的可信度。一个良好的实验设计需要在统计功效、成本控制和数据质量之间取得平衡。生物学重复和技术重复是两个容易混淆但本质不同的概念。生物学重复是指来自不同个体或独立培养的样本，它们之间的变异反映的是真实的生物变异。技术重复是指同一样本经过独立的文库制备和测序过程产生的数据，它们之间的变异主要来源于实验操作和仪器误差。在差异表达分析中，生物学重复是必需的，因为它提供了估计组内变异的基础；而技术重复通常不是必需的，除非你需要评估实验流程的稳定性。当你计划一项RNA-seq实验时，至少需要3个生物学重复才能进行有意义的统计分析，更多的重复（如5-6个）可以提高检测低表达基因差异表达的灵敏度。 重复数的确定应该基于统计功效分析（power analysis）。统计功效是指在真实存在差异的情况下能够检测到该差异的概率，它与效应大小、组内变异、样本量和显著性水平相关。在RNA-seq实验中，常用的功效分析工具包括RNASeqPower、powsimR和ssize等R包。这些工具可以根据预期的差异倍数（如2倍变化）、目标检测到的最小表达水平、期望的功效（通常为0.8）和假阳性率（通常为0.05），计算出所需的生物学重复数。当你的研究对象是异质性较大的群体（如临床患者样本）时，可能需要更多的重复来克服较大的组内变异；而当使用遗传背景一致的模型动物或细胞系时，较少的重复可能就足够了。 测序策略的选择涉及单端测序和双端测序两种模式。单端测序只从片段的一端读取序列，成本较低且操作简单，适用于基因表达定量等不需要精确配对信息的分析。双端测序从片段的两端分别读取序列，产生成对的读段，虽然成本更高但提供了更丰富的信息。双端测序的优势在于：可以提高比对准确性，特别是在跨越剪接点或重复区域时；可以准确估计插入片段长度分布；有助于区分多映射读段的真实来源；对于可变剪接和新转录本发现等应用更为重要。当你需要进行复杂的转录组分析（如融合基因检测、新转录本鉴定）时，双端测序是更好的选择；如果仅用于常规的差异表达分析，单端测序在控制成本的前提下也能满足需求。 测序深度的选择取决于你的分析目的。对于基因表达谱分析和差异表达检测，每个样本20-30 million reads（单端）或15-25 million paired-end reads通常是足够的，这个深度可以保证大多数中等和高表达基因被可靠地检测到。如果你需要鉴定稀有转录本或进行全面的转录本组装，则需要更高的深度（50-100 million reads甚至更多）。测序深度还与基因组复杂度相关：人类和小鼠等哺乳动物的转录组比细菌或酵母更复杂，因此需要更高的测序深度来覆盖所有的转录本。此外，polyA富集策略捕获的是mRNA，其复杂度低于rRNA去除策略捕获的总RNA，因此在相同深度下polyA富集的有效覆盖更高。 链特异性文库（strand-specific library）保留了原始RNA分子的链方向信息，这对于正确识别反义转录本、重叠基因和双向启动子区域的转录活动至关重要。非链特异性文库在反转录过程中丢失了方向信息，无法区分一条read来自正义链还是反义链的转录本。dUTP法是最常用的链特异性建库方法之一，它在第二链cDNA合成时掺入dUTP，使得第二链在PCR扩增时不能有效扩增，从而保留第一链的方向信息。Illumina链特异性试剂盒采用类似的原理，通过化学修饰实现链特异性的保留。RNA降解预处理法利用了RNA降解的方向性特征来实现链特异性。当你研究的基因位点存在重叠的正负链转录本时（如某些lncRNA与mRNA重叠的情况），链特异性文库可以避免错误地将reads分配给错误的转录本，从而提高定量的准确性。 混池测序（multiplexing）是在同一个测序通道中同时运行多个样本的技术，通过在每个样本的文库上添加独特的索引序列（index/barcode）来区分不同来源的数据。合理的index设置需要注意避免索引跳跃（index hopping）问题，即Illumina平台上的某些模式中，一个样本的索引可能被错误地分配给另一个样本的读段。为了解决这个问题，推荐使用unique dual index（UDI）方案，即每对读段都有两个唯一的索引。当你需要在有限的预算内处理大量样本时，混池测序可以显著降低单个样本的成本，但需要注意不要过度混池导致每个样本的实际测序深度不足。 批次效应是影响RNA-seq数据质量的重要因素，它指的是由于实验在不同时间、由不同操作人员、使用不同试剂批次等因素引起的系统性偏差。批次效应可以在数据中引入与生物学因素无关的变异来源，如果不加以控制可能导致错误的结论。随机化是减轻批次效应的基本策略，即确保每个处理组的样本均匀地分布在各个批次中，而不是将同一处理组的所有样本集中在一个批次中完成。当你无法完全消除批次效应时，应该在实验设计中包含批次信息，以便在后续数据分析中使用ComBat等方法进行校正。 Spike-in对照是外源添加的已知浓度的RNA分子，用于评估和校正实验中的技术变异。ERCC（External RNA Controls Consortium） spike-in是一套广泛使用的合成RNA混合物，包含了92个不同浓度梯度的外源RNA。Spike-in对照的主要用途包括：监测反转录效率的一致性、评估文库制备的偏差、作为归一化的外部参考、检测样本间的技术变异。需要注意的是，spike-in的使用并不总是必要的，对于大多数标准的差异表达分析，基于内参基因的归一化方法已经足够。当你需要进行跨物种比较或研究极端条件下的全局表达变化时，spike-in对照可能更有价值。 样本收集与RNA提取质控是整个RNA-seq流程的第一步，也是最容易出问题的环节之一。RNA的质量直接影响后续的分析结果，特别是对于需要全长转录本信息的应用（如全长RNA-seq）。RNA完整性通常用RIN值（RNA Integrity Number）来衡量，范围从1（完全降解）到10（完整无损）。一般来说，RIN值大于7的RNA适合大多数RNA-seq应用，而RIN值大于8则更适合需要高质量RNA的应用。28S/18S rRNA比例是另一种评估RNA质量的指标，完整RNA的比例应接近2:1。当你的样本难以获取新鲜组织时（如临床FFPE样本），即使RIN值较低也可以进行RNA-seq分析，但需要对结果进行谨慎解读。 文库制备方法的选择取决于你的研究目标和样本类型。PolyA富集方法通过oligo(dT)引物选择性富集带polyA尾的mRNA，这种方法获得的文库主要由编码蛋白的mRNA组成，背景噪音较低，适合基因表达定量分析。rRNA去除方法使用探针杂交或酶切方式去除rRNA，保留了包括mRNA、lncRNA、circRNA和其他非编码RNA在内的总RNA成分，适合于全面研究转录组的方法。全长文库制备方法（如SMARTer、Oxford Nanopore全长RNA-seq）可以获得完整的转录本序列，对于新转录本发现和可变剪接分析特别有价值。当你主要关注蛋白编码基因的表达变化时，polyA富集是经济高效的选择；如果你想研究非编码RNA或发现新的转录本，rRNA去除或全长方法是更好的选择。 单细胞RNA-seq（scRNA-seq）文库平台的发展使我们能够在单个细胞的分辨率下研究基因表达。10X Genomics是目前最流行的商业化scRNA-seq平台，它基于微液滴技术，每次运行可以同时处理数千到数万个细胞，具有高通量和相对标准化的优势。Drop-seq是一个开源的低成本替代方案，原理类似但需要更多的手工优化。Smart-seq2是一种基于孔板的全长scRNA-seq方法，虽然通量较低（通常数百个细胞）但可以获得完整的转录本信息，更适合于可变剪接和等位基因特异性表达的研究。当你需要研究细胞群体的异质性或发现稀有的细胞亚群时，10X Genomics是首选平台；如果你需要高分辨率的转录本信息，Smart-seq2可能更合适。 低起始量RNA-seq方案解决了只有极少量RNA可用时的测序需求。Smart-seq系列方法利用模板转换机制，可以从皮克级的总RNA开始构建文库，非常适合单细胞或少量的显微切割组织。SMARTer系列商业试剂盒进一步简化了这一过程，提供了针对超低输入量的优化方案。这些方法的共同特点是使用了特殊的引物设计和扩增策略，以最大限度地减少扩增偏倚并保持表达的相对定量关系。当你只能获取极少量的样本材料（如激光显微切割的组织切片、流式分选的少量细胞）时，这些低起始量方案是实现RNA-seq的唯一途径。 .. code-block:: r # RNA-seq统计功效分析 library(RNASeqPower) # 参数设置: 深度=20M, 离散度=0.4, alpha=0.05, power=0.8, 目标检测2倍变化 result <- RNASeqPower(nreps = 3:10, depth = 20, disp = 0.4, alpha = 0.05, power = 0.8, effect = log2(2)) print("各重复数对应的可检测最小效应:") print(result) # 推荐每组至少3-6个生物学重复 # 测序深度估算 estimate_depth <- function(target_genes = 20000, avg_reads_per_gene = 30) { total_reads <- target_genes * avg_reads_per_gene depth_millions <- total_reads / 1e6 cat(sprintf("推荐测序深度: ~%dM reads\n", ceiling(depth_millions))) cat(sprintf("双端模式: ~%dM paired-end reads\n", ceiling(depth_millions / 2))) return(depth_millions) } estimate_depth() 5.2 RNA-seq 数据预处理 ----------------------- RNA-seq数据预处理是将原始测序数据转换为可用于下游分析的干净比对数据的过程，这一步骤的质量直接影响最终结果的可靠性。原始数据的质量评估是预处理的第一步，FastQC是最常用的质量评估工具，它可以生成关于碱基质量分布、GC含量、接头污染、重复序列等多方面的报告。MultiQC可以将多个FastQC报告汇总为一个综合报告，便于批量比较多个样本的质量状况。当你拿到一批新的RNA-seq数据时，首先运行FastQC检查是否存在明显的质量问题（如末端碱基质量下降、接头未去除、异常GC含量等），这些问题需要在后续步骤中加以处理。 接头修剪是去除测序读段末端的接头序列的过程。接头序列是在文库制备过程中加入的，它们的存在会干扰比对结果并降低有效读段的长度。cutadapt是一个轻量级高效的接头修剪工具，支持精确匹配和容错匹配模式，可以处理各种类型的接头序列。Trimmomatic是另一个广泛使用的工具，除了接头修剪外还可以进行滑动窗口质量过滤和多步连续处理。fastp集成了质量控制、接头修剪和质量过滤于一体，运行速度快且自动生成HTML格式的质量报告。当你处理Illumina平台的RNA-seq数据时，需要根据所使用的建库试剂盒确定正确的接头序列，常见的 TruSeq 接头序列是 AGATCGGAAGAGC。 低质量碱基修剪是为了去除测序质量较差的碱基，这些碱基可能引入错误的比对或变异调用。Phred质量分数是衡量每个碱基测序可靠性的指标，Q30表示99.9%的正确率，Q20表示99%的正确率。滑动窗口（sliding window）是一种常用的质量过滤策略，它在读段上滑动一个固定大小的窗口，计算窗口内的平均质量分数，如果平均值低于阈值则截断该位置之后的所有碱基。这种策略比单纯按照固定位置截断更加灵活，因为质量下降的位置在不同读段之间可能有差异。当你设置质量过滤参数时，需要权衡保留尽可能多的有用数据和去除不可靠碱基之间的关系，过于严格的过滤会损失过多数据，过于宽松则会引入噪声。 去除rRNA reads是RNA-seq预处理的可选步骤，因为即使在polyA富集后仍可能有少量rRNA残留，或者在rRNA去除策略中未能完全去除。SortMeRNA是一个专门用于快速分类和去除rRNA的工具，它使用k-mer方法和Burrows-Wheeler变换来加速搜索过程。bbsplit属于BBMap套件的一部分，可以将reads分配到多个参考序列上，从而分离出比对到rRNA参考的reads。Bowtie2比对到rRNA参考是另一种方法，先将reads比对到rRNA数据库，然后提取未比对的reads用于后续分析。当你的数据显示rRNA比例异常偏高（超过总reads的5-10%）时，执行这一步可以显著减少无效数据并节省计算资源。 去除globin reads主要应用于血液样本的RNA-seq分析。Globin mRNA（血红蛋白mRNA）在全血样本中占主导地位，可以达到总mRNA的50-70%，这会严重稀释其他基因的表达信号并浪费大量的测序深度。针对小鼠和人类的globin去除方案可以使用特定的探针或抗体进行免疫耗竭，也可以在计算阶段通过比对到globin基因参考序列来过滤掉相关的reads。当你研究全血或外周血单核细胞（PBMC）的转录组时，globin去除几乎是必须的步骤，否则大部分测序资源将被globin mRNA占据。 重复序列过滤在某些情况下是有用的步骤。PCR扩增过程中的过度扩增会产生大量完全相同的duplicate reads，这些reads不代表原始分子数量的增加，而是PCR偏倚的结果。然而，在RNA-seq分析中，duplicate removal需要谨慎对待，因为高表达基因天然会产生大量相同的reads，盲目去除会导致表达量的低估。Picard的MarkDuplicates工具可以标记潜在的PCR duplicate，但在RNA-seq中通常建议仅标记而不删除，或者仅在非常高的duplicate rate时才考虑去除。当你使用UMI（Unique Molecular Identifier）技术的文库时，可以通过UMI信息准确地识别和去重真正的PCR duplicates。 比对前归一化是指在进行正式比对之前，先进行一轮快速的粗略比对以检测可能的污染源。例如，如果你的样本可能受到细菌或病毒污染，可以先快速比对到相应的参考基因组，评估污染的程度。这一步可以帮助你决定是否需要调整后续分析策略，或者是否需要重新考虑样本的处理方式。当你处理未知来源或高风险污染可能的样本（如环境样本、临床感染样本）时，比对前归一化可以作为质量控制的有益补充。 .. code-block:: bash # RNA-seq数据预处理流程 # 1. 质量评估 fastqc -t 8 -o qc/ raw_data/*.fastq.gz # 2. 汇总QC报告 multiqc -o qc/ qc/ # 3. 接头修剪和质量过滤 (fastp) fastp -i sample_R1.fastq.gz -I sample_R2.fastq.gz \ -o clean_R1.fq.gz -O clean_R2.fq.gz \ --detect_adapter_for_pe \ --cut_front --cut_tail \ --cut_window_size 4 --cut_mean_quality 20 \ --length_required 36 \ --json fastp.json --html fastp.html # 4. 去除rRNA reads (SortMeRNA) sortmerna --ref /path/to/silva-bac-16s-id.fasta,silva-bac-16s-id-index \ --reads clean_R1.fq.gz \ --reads clean_R2.fq.gz \ --other non_rRNA_1.fq.gz \ --other non_rRNA_2.fq.gz \ --paired_in rRNA_pairs.fq.gz \ --fastx --log -a 8 5.3 比对到参考基因组与转录组 ------------------------------ RNA-seq数据的比对与基因组DNA测序的比对有一个根本性的区别：RNA-seq读段可能跨越剪接点，这意味着一个读段的两个部分可能对应基因组上不相邻的外显子区域。传统的短读长比对器（如BWA-MEM）假设读段与参考序列基本连续，无法处理这种跨越较大间隙的比对情况。剪接感知比对器（splice-aware aligner）专门设计用于处理RNA-seq数据，它们能够在比对过程中识别和处理剪接连接（splice junction），将读段正确地定位到基因组上。 STAR（Spliced Transcripts Alignment to a Reference）是目前最广泛使用的RNA-seq比对器之一，它基于两阶段比对策略：第一阶段使用最大可扩展种子（maximally mappable seed）快速找到候选位置，第二阶段进行无gap和有gap的精确比对。STAR的速度非常快，可以在数小时内完成人类基因组的全转录组比对，同时保持了较高的灵敏度和准确性。HISAT2是TopHat2/HISAT的继承者，它采用分层索引FM-index（Hierarchical Graph FM Index），将基因组索引组织成层次结构，使得比对速度更快且内存占用更低。Subread套件中的align功能也支持剪接感知比对，特别适合内存受限的环境。RSEM（RNA-Seq by Expectation Maximization）结合了比对和定量功能，可以同时输出比对结果和基因表达量估计。Salmon在轻比对模式下可以直接从FASTQ文件进行定量，无需生成完整的BAM文件，大大提高了处理速度。 比对索引构建是使用任何比对器的第一步。对于STAR来说，索引构建需要提供参考基因组FASTA文件和可选的基因注释GTF文件，包含GTF信息可以使STAR在比对时利用已知的剪接位点，提高比对效率和准确性。STAR的索引构建参数中，--sjdbOverhang参数应该设置为读段长度减1，这是因为STAR使用剪接位点的锚定序列来识别剪接连接。HISAT2的索引同样基于参考基因组和基因注释，但其层次化的索引结构使得索引文件较小且加载迅速。当你更换参考基因组版本或更新基因注释时，必须重新构建索引以确保比对结果的一致性。 双端比对一致性处理是双端RNA-seq数据特有的问题。理想情况下，一对双端reads应该比对到基因组的相近位置且方向相反（--fr表示forward-reverse配对模式），插入片段长度应在预期范围内。然而，由于剪接的存在，双端reads可能跨越一个或多个外显子-内含子边界，导致插入片段长度看起来异常大。大多数剪接感知比对器会自动处理这种情况，但仍需注意检查配对信息是否正确保留。当你观察到大量异常配对的reads时，可能是比对参数设置不当或存在结构变异等问题。 比对后过滤是根据比对质量筛选可靠的比对结果。唯一比对（uniquely mapped）是指read只有一个最佳的比对位置，这类reads的定位最为可靠。多比对（multi-mapped）是指read有多个同等优秀的比对位置，通常发生在重复序列或基因家族成员之间。对于基因表达定量而言，如何处理多比对reads是一个重要的问题：简单地丢弃它们会导致某些基因家族的表达量被低估，而将它们随机分配到一个位置又可能引入噪声。常用的解决方案包括使用加权分配算法（如RSEM、eXpress）根据多比对位置的表达证据按比例分配reads，或者使用EM算法迭代优化分配方案。当你关注基因家族或重复区域的表达时，需要特别注意多比对reads的处理策略。 比对质量评估是验证比对结果可靠性的关键步骤。Qualimap RNA-seq QC模块专门针对RNA-seq数据提供详细的统计信息，包括比对率、外显子/内含子/基因间区的分布、5'/3'偏向性、覆盖度均一性等。RSeQC是一套全面的RNA-seq质量控制工具集合，它可以评估基因体覆盖均匀性、剪接位点质量、链特异性准确性、重复率等多个维度。Picard的CollectRnaSeqMetrics工具可以计算核糖体RNA比率、编码区覆盖率等重要指标。当你完成一批RNA-seq样本的比对后，运行这些QC工具并汇总结果可以帮助你识别异常样本并及时采取补救措施。 比对输出格式方面，SAM（Sequence Alignment/Map）是文本格式的比对结果，便于人工查看但不适合大规模存储。BAM（Binary Alignment Map）是SAM的二进制压缩格式，是实际存储和传输比对结果的标准格式。CRAM（Compressed Reference-oriented Alignment Map）是更新的压缩格式，它存储相对于参考序列的差异信息而非完整的序列，可以实现更高的压缩比（通常比BAM小40-60%）。当你需要长期存储大量RNA-seq比对数据时，CRAM是更节省空间的选择，但需要注意使用CRAM时需要随时访问参考基因组文件。 长读长RNA-seq比对需要专门的工具来处理长读长的特点。minimap2是目前最流行的长读长比对工具，它支持PacBio和Nanopore数据，具有极高的比对速度和不错的准确性。uLTRA是专门为长读长RNA-seq设计的剪接感知比对器，它可以更好地处理长读长跨越多个剪接点的情况。MAJIQ-LR结合了MAJIQ的剪接量化能力和长读长的优势，用于检测和定量可变剪接事件。当你使用Oxford Nanopore或PacBio Iso-Seq进行全长转录本研究时，选择合适的长读长比对工具至关重要。 无参考转录组比对适用于没有高质量参考基因组的物种或希望发现全新转录本的情况。de novo组装（从头组装）是指不依赖参考基因组，仅根据reads之间的overlap关系将它们组装成contig或转录本。常用的de novo组装工具包括Trinity、SOAPdenovo-Trans和Oases等。组装完成后，可以将组装得到的转录本作为参考进行后续的定量和分析。当你研究非模式物种或希望发现参考注释中不存在的转录本时，无参考转录组组装是必要的选择，但需要注意de novo组装的计算成本较高且结果质量受测序深度和数据质量的影响较大。 .. code-block:: bash # STAR比对流程 # 1. 构建基因组索引 STAR --runThreadN 16 \ --runMode genomeGenerate \ --genomeDir /path/to/STAR_index \ --genomeFastaFiles GRCh38.fa \ --sjdbGTFfile gencode.v44.annotation.gtf \ --sjdbOverhang 149 # 2. 剪接感知比对 STAR --runThreadN 16 \ --genomeDir /path/to/STAR_index \ --readFilesIn clean_R1.fq.gz clean_R2.fq.gz \ --readFilesCommand zcat \ --outFileNamePrefix sample1_ \ --outSAMtype BAM SortedByCoordinate \ --quantMode GeneCounts \ --twopassMode Basic # 3. HISAT2比对（替代方案） hisat2-build GRCh38.fa hisat2_index/grch38 hisat2 -p 16 -x hisat2_index/grch38 \ -1 clean_R1.fq.gz -2 clean_R2.fq.gz \ -S sample1.sam # 4. Salmon伪比对定量（无需完整比对） salmon index -t transcripts.fa -i salmon_index --type quasi -k 31 salmon quant -i salmon_index -l A \ -1 clean_R1.fq.gz -2 clean_R2.fq.gz \ -p 8 --validateMappings -o salmon_out 5.4 转录本定量 -------------- 转录本定量是RNA-seq分析的核心环节，它的目标是从比对后的数据中估计每个基因或转录本的表达丰度。定量的准确性直接影响后续差异表达分析和功能解释的可靠性。根据是否依赖完整的基因组比对，定量方法可以分为基于比对的方法和基于轻比对（伪比对）的方法两大类。 5.4.1 基于比对的定量 ~~~~~~~~~~~~~~~~~~~~ 基于比对的定量方法首先将reads比对到参考基因组或转录组，然后统计每个基因或转录本对应的reads数量。HTSeq-count是Python编写的计数工具，它提供了三种计数模式来处理reads与基因特征的重叠关系。Union模式是最常用的模式，如果一个read与多个基因特征重叠，只要其中任何一个特征被完全包含在read中，该read就被计入该特征。Intersection_strict模式要求read与基因特征的交集必须是非空的且交集的碱基数等于read的长度。Intersection_non_empty模式要求read与基因特征的交集非空即可。featureCounts是Subread套件的一部分，它是目前最广泛使用的计数工具之一，运行速度快且支持多种参数配置。 基因水平定量和转录本水平定量是两种不同粒度的定量策略。基因水平定量将所有属于同一基因的转录本的reads合并计数，得到每个基因的总表达量，这是差异表达分析中最常用的形式。转录本水平定量分别估计每个转录本的表达丰度，可以揭示同一基因内部不同转录本的相对比例变化，对于可变剪接分析尤为重要。然而，转录本水平的定量面临更大的挑战，因为不同转录本之间存在大量共享的外显子和reads，使得解析各转录本的具体贡献变得更加困难。 多比对reads的处理是基于比对定量的核心难题之一。当一个read可以等同地比对到多个位置时（例如基因家族的 paralogs 或重复元件），简单的丢弃或随机分配都会带来偏差。RSEM（RNA-Seq by Expectation Maximization）使用期望最大化算法来迭代地估计每个转录本的表达水平，并根据当前估计值为多比对reads分配概率权重。eXpress采用了类似在线EM的策略，可以在单次扫描数据的过程中完成定量，速度更快。这些方法的核心思想是：如果一个read可以比对到转录本A和B，那么它应该根据A和B的相对表达量按比例贡献给两者，而不是全部给其中一个。 片段（fragment）计数与read计数的区别在双端测序中尤为重要。一个片段（fragment）在文库制备过程中产生两个配对的reads（paired-end reads），这两个reads代表的是同一个原始cDNA分子的两端信息。正确的做法是统计片段的数量而不是read的数量，因为一个片段的两个reads不应该被当作两个独立的观察值来计数。HTSeq-count和featureCounts都提供了-p参数来指定使用片段计数模式，此时配对的reads会被作为一个整体来处理。当你使用双端数据进行定量时，务必确认使用片段计数模式，否则表达量会被大约高估一倍。 .. code-block:: r # 基于比对的定量 library(Rsubread) # featureCounts计数（推荐） fc <- featureCounts(files = list.files("bam_dir", "\\.bam$"), annot.inbuilt = "hg38", isGTFAnnotationFile = TRUE, GTF.featureType = "exon", GTF.attrType = "gene_id", useMetaFeatures = TRUE, nthreads = 8) counts <- fc$counts print(head(counts)) # HTSeq-count替代方案 # htseq-count -f bam -r pos -s reverse -t exon -i gene_id sample.bam annotation.gtf 5.4.2 基于比对的轻量定量（伪比对） ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 伪比对（pseudo-alignment或quasi-mapping）是一种革命性的转录本定量方法，它跳过了传统比对中耗时的完整比对步骤，而是通过快速近似的方式确定每个read最可能来自哪些转录本。Salmon和Kallisto是这一领域的代表性工具，它们的运行速度比传统流程快数十倍甚至上百倍，同时保持着相当的准确性。 伪比对的基本原理依赖于k-mer索引。工具预先为参考转录本构建k-mer索引，对于每个输入read，提取其中的k-mers并在索引中查找匹配的转录本。通过分析k-mer的共现模式和位置一致性，可以快速推断read最可能来自哪个转录本以及在该转录本上的大致位置。这个过程不需要进行完整的动态规划比对，因此速度极快。Salmon使用选择性映射算法（selective alignment）在初始快速定位后对候选位置进行局部验证，兼顾了速度和准确性。 转录本量化校正是伪比对工具的重要特性。原始的计数或映射速率会受到多种系统偏倚的影响，包括GC含量偏倚（GC含量极端的转录本可能被系统性低估或高估）、片段长度分布偏倚（较短的转录本产生较少的片段）、序列特异性偏倚（随机六聚体引发或片段化的偏好性）等。Salmon内置了丰富的偏倚校正模型，可以在定量过程中同时估计和校正这些偏倚因素的影响。Kallisto也提供了基本的偏倚校正选项。当你追求最高精度的定量结果时，启用Salmon的完整偏倚校正模型（--gcBias --seqBias）是推荐的设置。 准比对（quasi-mapping）是Salmon使用的核心技术术语，它描述了一种介于完全不比对和完全比对之间的中间策略。准比对不仅记录read与转录本的对应关系，还记录了大概的映射位置信息，这些信息足以用于后续的偏倚校正和有效性评估，但又避免了完整比对的高昂计算代价。 Salmon与Kallisto都是优秀的伪比对工具，但它们在一些细节上有所不同。Salmon提供了更丰富的偏倚校正选项和更细致的映射质量评估，支持选择性映射以提高准确性，并且可以处理GC偏倚和序列偏倚等多种偏倚类型。Kallisto的设计理念更为简洁，强调极致的速度和最小的内存占用，其bootstrapping功能可以快速生成表达量估计的不确定性区间。两者在大多数情况下给出高度一致的定量结果，选择哪一个很大程度上取决于个人偏好和具体需求。当你需要更精细的控制和更全面的偏倚校正时，Salmon可能是更好的选择；当你优先考虑简洁性和极致速度时，Kallisto值得考虑。 5.4.3 归一化方法 ~~~~~~~~~~~~~~~~~~ 原始计数（raw count）是定量工具的直接输出，表示每个基因或转录本被多少个reads或fragments覆盖。然而，原始计数不能直接用于样本间比较，因为它受到测序深度、基因长度和RNA组成等多种因素的干扰。归一化的目的是消除这些技术性因素带来的变异，使样本间的表达量具有可比性。 RPKM（Reads Per Kilobase per Million reads）和FPKM（Fragments Per Kilobase per Million reads）是最早广泛使用的归一化方法之一。RPKM = (10^6 × C) / (N × L / 1000)，其中C是某基因的原始计数，N是样本总reads数，L是基因长度（单位kb）。FPKM是RPKM的双端版本，用片段数代替reads数。这两种方法同时校正了测序深度（除以百万reads）和基因长度（除以kb），使得不同长度基因的表达量可以相互比较。然而，RPKM/TPKM的主要局限在于它们不适合跨样本比较，因为分母中的总reads数会受到少数极高表达基因的强烈影响，当一个样本中有几个超高表达基因时，其他所有基因的RPKM/FPKM都会被人为压低。 TPM（Transcripts Per Million）是对RPKM/FPKM的重要改进。TPM的计算分为两步：首先将每个基因的计数除以其长度得到每kb的reads数（RPK），然后将所有基因的RPK求和，最后将每个基因的RPK除以RPK总和再乘以10^6。关键区别在于TPM的归一化因子是所有基因RPK的和，这意味着所有样本的TPM之和恒定为10^6，使得TPM值在样本之间具有内在可比性。当你需要在多个样本之间比较基因的相对表达水平时，TPM是比FPKM更合适的选择。 TMM（Trimmed Mean of M-values）是edgeR包推荐的归一化方法。TMM的基本思想是：在大多数基因没有差异表达的前提下，找到一个缩放因子使得两个样本之间的表达量分布尽可能一致。具体做法是计算每个基因的对数比值（M值）和对数均值（A值），然后剔除极端的M值和A值（trimming），用剩余基因的加权平均M值作为归一化因子。TMM对于组成型差异（compositional differences，即少数基因的大幅度变化影响整体分布）具有良好的鲁棒性。当你使用edgeR进行差异表达分析时，TMM归一化是默认且推荐的选择。 DESeq2的中位数归一化（median-of-ratios）是该包的核心归一化方法。对于每个基因，计算其计数与所有样本中位数的比值，然后取所有基因比值的中位数作为该样本的归一化因子。这种方法假设大多数基因没有差异表达，因此中位数比值应该接近1。DESeq2的归一化因子被称为size factor，后续的所有分析都在此基础上进行。DESeq2的中位数归一化对于离群值具有较强的抵抗力，因为它使用中位数而非均值。当你使用DESeq2进行差异表达分析时，通常会自动应用中位数归一化，无需手动干预。 分位数归一化（Quantile normalization）强制使所有样本的表达量分布完全相同，即将每个样本的表达量排序后替换为所有样本在该排名处的平均值。这种方法最初开发用于微阵列数据分析，在RNA-seq中的应用较少见，因为它可能过度修正真实的生物学差异。相对对数表达（RLE）是通过计算每个基因与其在各样本中位数的比值来评估归一化效果的方法，常用于检测异常样本。上四分位数归一化（upper quartile normalization）使用上四分位数而非总数或中位数作为归一化基准。无参考的归一化方法使用一组已知稳定表达的管家基因（housekeeping genes）作为参照来进行归一化，这在缺乏完整参考基因组的情况下可能有用，但需要事先验证所选管家基因在你的实验条件下确实稳定表达。 .. code-block:: r # 归一化方法比较 library(edgeR) # 模拟计数数据 set.seed(123) counts <- matrix(rnbinom(6000, mu = rpois(100, 500), size = 0.1), nrow = 100, ncol = 6) rownames(counts) <- paste0("Gene", 1:100) colnames(counts) <- c("Ctrl1", "Ctrl2", "Ctrl3", "Treat1", "Treat2", "Treat3") gene_lengths <- sample(500:15000, 100) # TMM归一化 (edgeR) y <- DGEList(counts = counts) y <- calcNormFactors(y, method = "TMM") cat("TMM归一化因子:\n"); print(y$samples$norm.factors) # 计算RPKM/FPKM rpkm <- rpkm(y, gene.length = gene_lengths) fpkm <- rpkm / 2 # 双端数据除以2 cat("\nRPKM示例 (前5个基因):\n"); print(head(rpkm[, 1:2])) # 计算TPM (TPM更适合跨样本比较) tpm_function <- function(counts_matrix, lengths) { rpk <- t(t(counts_matrix) / (lengths / 1000)) tpm <- t(t(rpk) / colSums(rpk)) * 1e6 return(tpm) } tpm_vals <- tpm_function(counts, gene_lengths) cat("\nTPM每列总和检验 (应接近1e6):\n"); print(colSums(tpm_vals)) 5.5 差异表达分析 ---------------- 差异表达分析的目标是识别在不同实验条件（如处理vs对照、疾病vs健康、不同时间点）之间表达水平发生显著变化的基因。这是RNA-seq分析中最核心也是最常用的分析步骤，其结果直接指导后续的功能注释和机制探索。 5.5.1 经典工具与方法 ~~~~~~~~~~~~~~~~~~~~~~~ 基于负二项分布的方法是RNA-seq差异表达分析的主流框架。RNA-seq计数数据的特点是离散的、非负整数值，且均值与方差相关（即过离散现象），负二项分布恰好能很好地描述这种数据特性。DESeq2是最广泛使用的基于负二项分布的工具，它通过估计每个基因的离散度参数来建模均值-方差关系，然后拟合广义线性模型并进行统计检验。edgeR同样是基于负二项分布的经典工具，与DESeq2相比，edgeR在某些场景下（如样本量很小时）可能表现更好。DSS（Dispersion Shrinkage for Sequencing data）和baySeq也是基于负二项分布或类似分布的方法，各有特色。 基于线性模型的方法提供了更灵活的实验设计能力。limma-voom将limma（原本为微阵列数据开发的线性模型框架）与voom变换相结合，voom变换将计数数据转换为适合线性模型的log-CPM值的同时估计每个基因的均值-方差关系权重。sleuth是为伪比对工具（如Kallisto）的bootstrap输出设计的差异分析方法，它使用线性模型框架并考虑了定量不确定性。当你需要处理复杂的实验设计（如多因素交叉设计、配对样本、时间序列）时，基于线性模型的方法通常比单纯的负二项分布方法更具灵活性。 基于非参数方法不依赖于数据服从特定分布的假设，因此在某些情况下可能更稳健。NOIseq使用噪声分布来模拟技术噪声，通过比较观测到的差异与噪声分布来判断显著性。SAMseq是SAM（Significance Analysis of Microarrays）的RNA-seq版本，使用基于排列检验的非参数方法。当你对数据分布的假设不太确定，或者样本量太小不足以可靠地估计分布参数时，非参数方法可以作为备选方案。 基于贝叶斯方法将先验信息纳入差异表达分析中。baySeq使用经验贝叶斯方法估计不同表达模式的 posterior probability。EBSeq（Empirical Bayes Analysis of differential gene expression for RNA-seq data）专门用于识别差异表达和差异外显子/isoform使用。贝叶斯方法在小样本情况下往往表现更好，因为先验信息可以弥补数据不足的问题。当你处理的实验组样本量很少（如每组仅有2-3个重复）时，基于贝叶斯的方法可能给出更稳定的估计。 5.5.2 DESeq2 详解（关键步骤） ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ DESeq2是目前最流行的RNA-seq差异表达分析工具，理解其工作流程对于正确使用和解释结果至关重要。离散度估计是DESeq2的第一步，也是最具特色的一步。离散度（dispersion）参数描述了方差相对于均值的偏离程度。DESeq2首先为每个基因单独估计离散度（gene-wise dispersion），然后拟合一条离散度与均值关系的趋势曲线（fitted dispersion trend），最后将基因特异的离散度向趋势曲线收缩（shrinkage），得到最终的离散度估计。收缩估计的强度取决于基因的信息量：高表达基因的离散度估计较为可靠，收缩较少；低表达基因的离散度估计不确定，收缩较多。这种经验贝叶斯收缩策略有效地解决了小样本情况下离散度估计不稳定的问题。 拟合负二项广义线性模型是DESeq2的核心统计框架。DESeq2使用广义线性模型（GLM）来描述计数数据与实验变量之间的关系，链接函数为log link，误差分布为负二项分布。最简单的模型公式为 ~ condition，表示表达量仅受实验条件的影响。更复杂的模型可以包含多个协变量，如 ~ batch + condition 表示同时考虑批次效应和实验条件，~ time + treatment + time:treatment 表示包含时间、处理及其交互作用。模型矩阵的设计决定了DESeq2可以检验哪些对比（contrasts）。 Wald检验和似然比检验（LRT）是DESeq2提供的两种统计检验方法。Wald检验用于检验单个系数是否为零（如某个处理效应是否存在），它计算系数估计值与其标准误的比值，在零假设下渐近服从标准正态分布。似然比检验用于比较两个嵌套模型（如完整模型 vs 缺少某个因子的简化模型），它计算两个模型对数似然之差的两倍，在零假设下渐近服从卡方分布。LRT特别适用于筛选性分析（screening analysis），即寻找在任何水平上有变化的基因而不指定具体的对比方向。当你有明确的两组比较需求时，Wald检验更常用；当你想探索多个水平之间的整体差异时，LRT更合适。 独立过滤（independent filtering）是DESeq2提高统计功效的重要策略。多重检验校正（如Benjamini-Hochberg FDR）需要同时对所有基因进行校正，基因数量越多，每个基因需要达到的阈值就越严格。独立过滤的思想是：对于那些在任何样本中表达量都很低的基因，它们几乎不可能达到显著性阈值，提前将这些基因从检验中排除可以减少多重检验的负担，从而放宽剩余基因的FDR阈值。DESeq2默认使用每个基因的平均归一化计数作为过滤标准，过滤掉平均计数过低（如小于10）的基因。独立过滤不会增加假阳性率，因为在过滤标准与检验统计量独立的条件下，它只是减少了不必要的检验次数。 批次效应包含在模型设计中是处理批次影响的推荐方法。如果在实验设计阶段未能完全随机化，或者事后发现有显著的批次效应，可以将batch作为一个协变量加入到DESeq2的设计公式中。例如，design = ~ batch + condition 告诉DESeq2在检验condition效应之前先校正batch的影响。这种方法比先进行数学上的批次校正（如ComBat）再进行差异分析更为严谨，因为它将批次效应明确地纳入统计模型中，可以正确地估计自由度和不确定性。当你知道样本的批次信息时，始终建议将其包含在设计公式中。 多因素设计允许你同时考察多个实验变量的主效应和交互效应。例如，在一个包含时间和处理两因素的时间序列实验中，设计公式 ~ time + treatment + time:treatment 可以检验处理效应是否随时间变化。交互项 time:treatment 的系数告诉你处理组与对照组之间的差异在不同时间点上是否不同。DESeq2的结果函数允许你通过name参数指定要提取哪个系数的检验结果，或通过contrast参数自定义对比。当你研究的问题涉及多个因素时，仔细思考设计公式的写法和需要检验的对比非常重要，不恰当的设计可能导致无法回答你想问的科学问题。 收缩估计（shrinkage）用于log2 fold change（LFC）是DESeq2的一个强大功能。普通的极大似然估计给出的LFC对于低表达基因或低计数基因可能极其不稳定（如一个基因在对照组中计数为0而在处理组中计数为1，LFC将为无穷大）。DESeq2提供了apeglm和ashr两种收缩方法，它们将不稳定的LFC估计值向零或其他先验值收缩，使得最终的LFC排序更加可靠。收缩后的LFC特别适合用于基因排序和后续的功能富集分析，因为它们避免了极端值主导排序结果。当你准备基于LFC进行基因排名或绘制火山图时，务必使用收缩后的LFC（如results()函数中的coef="apeglm"选项）。 .. code-block:: r # DESeq2完整差异表达分析流程 library(DESeq2) # 1. 创建DESeqDataSet sample_info <- data.frame( row.names = colnames(counts), condition = factor(c(rep("Control", 3), rep("Treatment", 3))), batch = factor(c("B1", "B2", "B1", "B2", "B1", "B2")) ) dds <- DESeqDataSetFromMatrix( countData = counts, colData = sample_info, design = ~ batch + condition ) # 2. 预过滤低表达基因 keep <- rowSums(counts(dds)) >= 10 dds <- dds[keep, ] # 3. 运行DESeq2标准分析（包含离散度估计、模型拟合、Wald检验） dds <- DESeq(dds) # 4. 提取结果（使用apeglm收缩LFC） res <- results(dds, name = "condition_Treatment_vs_Control", alpha = 0.05, coef="apeglm") resOrdered <- res[order(res$padj), ] # 5. 结果摘要 summary(resOrdered) cat("\n显著差异基因数 (FDR<0.05):", sum(resOrdered$padj < 0.05, na.rm=TRUE), "\n") cat("上调基因数:", sum(resOrdered$log2FoldChange > 1 & resOrdered$padj < 0.05, na.rm=TRUE), "\n") cat("下调基因数:", sum(resOrdered$log2FoldChange < -1 & resOrdered$padj < 0.05, na.rm=TRUE), "\n") 5.5.3 edgeR 详解（关键步骤） ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ edgeR是与DESeq2齐名的经典RNA-seq差异表达分析工具，它同样基于负二项分布框架，但在一些实现细节上有所不同。估计共同离散度（common dispersion）是edgeR离散度估计的第一步，它假设所有基因共享同一个离散度值，通过条件极大似然估计得到这个共同值。共同离散度给出了一个粗糙的全局离散度估计，作为后续精细化估计的起点。在实际数据中，不同基因的离散度显然不可能完全相同，因此共同离散度只是一个初步估计。 估计趋势离散度（trended dispersion）是edgeR的第二步。edgeR观察到离散度通常与基因的表达水平相关：低表达基因倾向于有更高的离散度（更大的相对变异），高表达基因的离散度趋于稳定。趋势离散度为每个表达水平预测了一个期望的离散度值，形成了一条离散度-均值趋势线。这条趋势线的形状和位置取决于数据本身的特性，edgeR通过加权回归来拟合这条曲线。趋势离散度的引入使得edgeR可以利用基因表达量信息来改善离散度估计。 基因特异离散度（tagwise dispersion）是edgeR离散度估计的最后一步。在趋势离散度的基础上，edgeR为每个基因估计一个特异的离散度值，并通过经验贝叶斯方法将其向趋势线收缩。收缩的程度取决于每个基因的信息量：拥有更多样本和更高计数的基因，其离散度估计更可信，收缩较少；反之则收缩较多。这与DESeq2的收缩策略类似，但edgeR的三步估计过程（common → trended → tagwise）使其离散度建模更加透明和可控。 准似然（quasi-likelihood）F检验是edgeR推荐的统计检验方法。与传统的似然比检验不同，准似然方法通过估计额外的dispersion参数来考虑模型拟合的不确定性，这使得它在小样本情况下更加保守和稳健。QL F检验类似于t检验相对于z检验的关系：当样本量大时两者趋近一致，但当样本量小时QL方法更能防止假阳性。当你使用edgeR进行分析时，推荐使用qlFTest而非传统的exactTest或LRT，除非你有特殊理由偏好后者。 使用TMM归一化是edgeR的标准配置。edgeR的calcNormFactors函数实现了TMM归一化，它会计算每个样本的归一化因子并将其存储在DGEList对象中。后续的所有分析（离散度估计、模型拟合、检验）都会自动使用这些归一化因子。TMM归一化与edgeR的负二项分布框架配合良好，是edgeR工作流的标配组合。当你创建edgeR的DGEList对象时，记得调用calcNormFactors进行TMM归一化。 .. code-block:: r # edgeR完整差异表达分析流程 library(edgeR) # 1. 创建DGEList对象 group <- factor(c(rep("Control", 3), rep("Treatment", 3))) y <- DGEList(counts = counts, group = group) # 2. TMM归一化 y <- calcNormFactors(y, method = "TMM") cat("TMM归一化因子:\n"); print(y$samples$norm.factors) # 3. 离散度估计 (common -> trended -> tagwise) y <- estimateDisp(y, design = model.matrix(~group), robust = TRUE) plotBCV(y, main = "离散度-均值关系") # 可视化 # 4. 准似然F检验 (推荐) design <- model.matrix(~group) fit <- glmQLFit(y, design) qlf <- glmQLFTest(fit, coef = 2) # 检验Treatment vs Control # 5. 提取结果 topTags(qlf, n = 10) res_edgeR <- topTags(qlf, n = Inf)$table sig_genes <- subset(res_edgeR, FDR < 0.05 & abs(logFC) > 1) cat("\nedgeR显著差异基因数:", nrow(sig_genes), "\n") 5.5.4 limma-voom ~~~~~~~~~~~~~~~~ limma-voom是将微阵列分析领域成熟的limma框架适配到RNA-seq数据的桥梁。voom变换的核心思想是将计数数据转换为适合线性模型分析的形式，同时保留计数数据的均值-方差关系信息。具体而言，voom对每个基因计算log2(counts-per-million)，然后拟合log-CPM与均值之间的关系，得到每个观测值的精确权重。这些权重随后传递给limma的线性模型过程，使得高表达基因（方差估计更可靠）获得更高的权重，低表达基因（方差估计不确定）获得较低的权重。 从计数转换为log-CPM是voom的第一步。首先使用TMM或其他方法对原始计数进行归一化，然后加上一个小的偏移量（offset，通常为0.5或1）以避免取对数时出现log(0)，最后乘以10^6并取log2。转换后的log-CPM值近似服从正态分布，满足了线性模型的前提假设。当你查看voom的输出图时，会看到mean-variance trend图显示方差随平均表达量的变化趋势，这是voom的关键诊断信息。 均值-方差关系建模是voom的核心创新。voom使用局部加权回归（lowess/loess）拟合log-CPM均值与方差之间的关系曲线，然后根据每个基因的平均表达量从曲线上读取预测的方差，取倒数作为该基因所有观测值的权重。这种建模方式巧妙地捕捉了RNA-seq数据中方差随表达量变化的特性，使得limma的经验贝叶斯平滑方法可以有效运作。当你运行voom后，务必检查mean-variance plot确保趋势合理，如果出现异常模式可能提示数据存在问题。 使用limma的线性模型和经验贝叶斯是voom之后的自然延续。一旦获得了加权后的log-CPM值，就可以像分析微阵列数据一样使用lmFit拟合线性模型、eBayes进行经验贝叶斯平滑、topTable提取差异表达结果。limma框架的最大优势在于其灵活的实验设计能力： contrasts矩阵可以定义任意复杂的对比组合，duplicateCorrelation可以处理相关性样本（如重复测量），removeBatchEffect可以方便地去除批次效应。当你需要处理复杂的实验设计时，limma-voom往往比DESeq2/edgeR提供更直观和灵活的操作界面。 适用于复杂实验设计是limma-voom的一大优势。无论是配对样本、重复测量、多因素析因设计还是纵向数据，limma都可以通过设计矩阵和对比矩阵的组合来优雅地处理。例如，配对样本只需在设计中添加subject作为block factor，时间序列分析可以使用splines或多项式来建模时间趋势，多批次的实验可以通过ComBat或设计矩阵中的batch term来校正。当你面对一个实验设计复杂的项目时，limma-voom可能是最省心的选择。 .. code-block:: r # limma-voom完整流程 library(limma) library(edgeR) y <- DGEList(counts = counts, group = group) y <- calcNormFactors(y) v <- voom(y, design = model.matrix(~group), plot = TRUE) design <- model.matrix(~ 0 + group) colnames(design) <- levels(group) fit <- lmFit(v, design) contrast.matrix <- makeContrasts(Treatment - Control, levels = design) fit2 <- contrasts.fit(fit, contrast.matrix) fit2 <- eBayes(fit2) topTable(fit2, coef = "Treatment - Control", number = 10) 5.5.5 差异表达后处理 ~~~~~~~~~~~~~~~~~~~~~~~ 多重检验校正是差异表达分析中不可或缺的步骤。当你同时检验数千甚至上万个基因时，即使在没有真实差异的情况下，也会有相当数量的基因仅仅因为随机波动而显示出统计学显著性。Benjamini-Hochberg（BH）方法控制的False Discovery Rate（FDR）是RNA-seq分析中最常用的多重检验校正方法，它控制的是所有被声明为显著的结果中错误发现的比例的期望值。Bonferroni校正更为严格，它控制Family-Wise Error Rate（FWER），即犯至少一次错误的概率，但由于过于保守在RNA-seq中较少使用。q-value是Storey提出的一种改进方法，它估计了真实原假设的比例（π0），在保持FDR控制的同时可能比BH方法更有力。当你报告差异表达结果时，FDR < 0.05 是最常见的显著性阈值。 log2倍数变化（LFC）阈值用于筛选生物学意义显著的差异表达基因。统计学显著性（低FDR）不一定意味着生物学重要性：一个基因可能有统计显著的变化但实际变化幅度很小（如1.05倍），这种变化可能在生物学上是无关紧要的。常见的LFC阈值包括 ``|LFC|`` > 1（即2倍变化）和 ``|LFC|`` > 2（即4倍变化），具体选择取决于研究领域和实验性质。有些研究者主张不预设LFC阈值，而是将所有基因按LFC排序供下游分析使用。当你设定LFC阈值时，需要在灵敏度和特异性之间做出权衡：过于宽松的阈值会引入大量弱效应基因，过于严格的阈值可能遗漏重要的调控基因。 增强列表是差异表达分析的最终输出产品，它应该包含足够的信息供后续解读使用。一个完善的差异表达结果表通常包括：基因ID、基因符号（symbol）、基因描述（description）、对照组平均表达量、处理组平均表达量、log2 fold change、p值、校正后p值（padj/FDR）、baseMean（平均表达量）等信息。你可以使用biomaRt或AnnotationDbi等R包从数据库中获取基因符号和描述信息。当你准备提交论文或与合作者分享结果时，一个信息完整、格式清晰的增强列表是必不可少的。 结果可视化是展示和传达差异表达分析发现的重要手段。MA图（MA plot）以log2 fold change（M值）为纵轴、平均表达量（A值）为横轴，展示了差异表达的整体格局，显著差异的基因通常以颜色突出显示。火山图（volcano plot）以-log10(p-value)为纵轴、log2 fold change为横轴，同时展示了统计显著性和效应大小，是最常用的差异表达可视化方式。热图（heatmap）和聚类热图（clustered heatmap）展示了差异基因在各样本中的表达模式，可以直观地看出样本分组和基因聚类的一致性。EnhancedVolcano、pheatmap、ComplexHeatmap等R包提供了高质量的绑图功能。当你撰写论文或做学术报告时，精心设计的可视化图表可以极大地提升结果的可读性和说服力。 LFC收缩对排序的影响是一个容易被忽视但很重要的问题。未经收缩的LFC估计值对于低表达基因可能极度不稳定，如果直接使用这些LFC进行基因排序，排序结果可能被极端值扭曲。收缩后的LFC（如DESeq2+apeglm或ashr得到的LFC）更加平滑和可靠，基于收缩LFC的排序更能反映真实的生物学差异。当你进行基因集富集分析（GSEA）时，输入基因列表的排序方式会直接影响富集结果，因此使用收缩后的LFC进行排序是推荐的做法。 差异表达基因集输出是差异表达分析的最终交付物。除了完整的增强列表外，你可能还需要导出特定阈值的基因子集（如FDR < 0.05且|LFC| > 1的基因列表），用于后续的功能富集分析、通路分析或实验验证。导出格式可以是CSV、TSV或Excel，确保包含足够的元数据信息以便追踪。当你完成差异表达分析后，保存好完整的分析脚本和中间对象（如dds、dge等R对象），这对于结果的可重现性和后续的审稿修改都非常重要。 .. code-block:: r # 差异表达结果可视化 library(EnhancedVolcano) library(pheatmap) res_df <- as.data.frame(resOrdered) res_df$gene <- rownames(res_df) # 火山图 EnhancedVolcano(res_df, lab = res_df$gene, x = "log2FoldChange", y = "padj", pCutoff = 0.05, FCcutoff = 1, title = "差异表达火山图", subtitle = "DESeq2 + apeglm", pointSize = 2.0, labSize = 3.0) # 热图 (top50差异基因) top50_genes <- head(rownames(resOrdered)[order(resOrdered$padj)], 50) norm_counts <- assay(vst(dds))[top50_genes, ] pheatmap(norm_cols, annotation_col = sample_info["condition", drop=FALSE], show_rownames = FALSE, cluster_rows = TRUE, main = "Top50 差异基因热图") # 导出增强列表 sig_res <- subset(res_df, padj < 0.05 & abs(log2FoldChange) > 1) write.csv(sig_res, "DEG_results_FDR0.05_LFC1.csv", row.names = FALSE) cat("导出显著差异基因数:", nrow(sig_res), "\n") 5.6 功能富集分析 ------------------- 差异表达分析给出了哪些基因发生了变化，而功能富集分析试图回答这些变化了的基因在功能上有什么共同点。通过将差异表达基因与已知的基因功能注释数据库进行比较，富集分析可以发现显著富集的生物学通路、分子功能和细胞组分，从而帮助理解实验处理或疾病状态的潜在生物学机制。 5.6.1 过表示分析（ORA） ~~~~~~~~~~~~~~~~~~~~~~~~ 过表示分析（Over-Representation Analysis，ORA）是最经典和最直观的富集分析方法。ORA的基本逻辑是：在一组差异表达基因中，某个功能类别（如GO term或KEGG pathway）中的基因数量是否显著多于随机期望？如果某种功能的基因在差异表达列表中出现得过于频繁，说明该功能可能与所研究的生物学过程相关。ORA的统计基础是超几何检验或Fisher精确检验，它们都是在给定总体中各类别基因数量的前提下，计算观察到当前富集程度或更极端情况的概率。 超几何检验和Fisher精确检验是ORA的两种统计方法，它们在大多数情况下给出几乎相同的结果。超几何检验从有限总体不放回抽样的角度建立模型，而Fisher精确检验从2x2列联表的角度出发。两者的区别主要在于理论出发点不同，实际应用中选择哪一种影响不大。当你使用clusterProfiler等R包进行ORA时，通常不需要手动选择检验方法，包内部已经做了合适的默认设置。 背景基因集选择是ORA中一个关键但经常被忽视的问题。背景基因集（background/universe set）代表了你在理论上可能检测到差异表达的所有基因。常见的选择包括：所有在样本中被检测到表达的基因（如count > 10的基因）、基因组中的所有蛋白质编码基因、或者芯片/RNA-seq平台覆盖的所有基因。不同的背景基因集选择可能导致截然不同的富集结果。一般推荐使用所有表达基因作为背景，因为这反映了你的实验实际能够检测的范围。当你进行ORA时，务必明确说明你所使用的背景基因集是什么，并在不同分析中保持一致。 多重检验校正在ORA中同样必要，因为你可能同时检验数千个GO terms或数百条通路。Benjamini-Hochberg FDR校正仍然是首选方法。考虑到GO terms之间的层级结构和高度重叠性，一些研究者还开发了考虑GO结构的校正方法（如topGO中的elim和weight算法），它们可以减少冗余的显著结果。当你报告ORA结果时，FDR < 0.05 或 FDR < 0.01 是常见的显著性阈值，同时建议报告富集倍数（enrichment fold）和涉及的基因数量。 代表性工具方面，clusterProfiler是Bioconductor中最流行和功能最全面的富集分析R包，它支持GO、KEGG、Reactome、DO（疾病本体论）等多种数据库的富集分析，并提供丰富的可视化功能（dotplot、cnetplot、emapplot等）。GOstats和topGO是较早的GO富集分析工具，topGO特别擅长处理GO的层级结构以减少冗余。当你需要一个一站式解决富集分析和可视化的方案时，clusterProfiler通常是最佳选择。 .. code-block:: r # GO/KEGG过表示分析 (ORA) library(clusterProfiler) library(org.Hs.eg.db) gene_list <- rownames(sig_res) go_bp <- enrichGO(gene = gene_list, OrgDb = org.Hs.eg.db, keyType = "SYMBOL", ont = "BP", pAdjustMethod = "BH", qvalueCutoff = 0.05, readable = TRUE) kegg_res <- enrichKEGG(gene = gene_list, organism = "hsa", pAdjustMethod = "BH", qvalueCutoff = 0.05) cat("GO BP显著富集项数:", nrow(go_bp@result), "\n") 基因集数据库是富集分析的信息来源。Gene Ontology（GO）包括三个分支：生物学过程（Biological Process，BP）描述基因参与的生物学活动，细胞组分（Cellular Component，CC）描述基因产物所在的细胞位置，分子功能（Molecular Function，MF）描述基因产物的分子活性。KEGG（Kyoto Encyclopedia of Genes and Genomes）是代谢和信号转导通路数据库，在生物医学研究中应用极为广泛。Reactome是另一个高质量的通路数据库，涵盖了更广泛的物种和更详细的反应层级。WikiPathways是社区驱动的通路数据库，内容更新较快。Hallmark基因集是MSigDB中的一个精选集合，包含50个经过精心整理的代表核心生物学过程的基因集，特别适合作为GSEA的输入。当你选择基因集数据库时，需要考虑研究物种的覆盖情况和数据库的维护更新频率。 5.6.2 基因集富集分析（GSEA） ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 基因集富集分析（Gene Set Enrichment Analysis，GSEA）与ORA的根本区别在于：GSEA使用所有基因的排序信息而非仅使用差异表达基因列表。GSEA的基本原理是：如果某个基因集中的基因确实在某种生物学过程中协同发挥作用，那么它们在整个基因排序列表中应该倾向于聚集在顶部（上调）或底部（下调），而不是随机散布。GSEA通过计算一个游走统计量（running enrichment score）来量化这种聚集程度，并通过排列检验（permutation test）评估其统计显著性。 预排序GSEA（pre-ranked GSEA）是最常用的GSEA模式。首先需要将所有基因按照某种指标排序，通常使用log2 fold change或t统计量作为排序依据。排序后的基因列表作为GSEA的输入，GSEA从列表底部开始向上遍历，遇到基因集中的基因时增加富集分数，遇到非基因集基因时减少分数。富集分数的最大偏离值（正或负）即为该基因集的富集得分（Enrichment Score，ES），其绝对值大小反映了富集程度。标准化后的ES（NES）考虑了基因集大小的因素，使得不同大小的基因集之间可以比较。当你已经有了差异表达分析的结果时，预排序GSEA是最直接的GSEA应用方式。 经典GSEA软件是由Broad Institute开发的Java应用程序，提供了图形用户界面和命令行两种使用方式。它内置了大量预定义的基因集集合（MSigDB），支持多种排列检验策略（基因排列 vs 样本排列），并生成详细的报告和可视化图表。经典GSEA软件的输出包括富集图（enrichment plot，显示ES沿排序列表的走势）、前导基因（leading-edge subset，对富集贡献最大的基因子集）和详细统计表格。当你首次学习GSEA或需要生成符合发表标准的图表时，经典GSEA软件是一个很好的起点。 fgsea（fast Gene Set Enrichment Analysis）是GSEA的高效R语言实现。fgsea的核心优势在于速度：它使用优化的算法使得即使是大规模的基因集数据库（如MSigDB的全部集合）也可以在数秒到数分钟内完成分析。fgsea接受预排序的基因列表和基因集数据库作为输入，输出NES、p-value、FDR等统计量以及leading-edge基因信息。fgsea与tidyverse生态系统兼容良好，便于整合到现代R分析流程中。当你需要在R中进行大规模GSEA分析时，fgsea是首选工具。 岭分析（ridgeplot）是GSEA的一种可视化变体，它展示了基因集中各基因在排序分布中的密度曲线。与经典的enrichment plot相比，ridgeplot可以同时展示多个基因集的分布模式，便于比较不同基因集的富集特征。clusterProfiler的ridgeplot函数和fgsea的plotEnrichment函数都可以生成岭图。当你需要在一个图中比较多个基因集的富集模式时，ridgeplot是非常有效的可视化手段。 .. code-block:: r # GSEA (fgsea) library(fgsea); library(msigdbr) gene_ranks <- setNames(resOrdered$log2FoldChange, rownames(resOrdered)) gene_ranks <- sort(gene_ranks, decreasing = TRUE) hallmark <- msigdbr(species = "Homo sapiens", category = "H") pathways <- split(hallmark$gene_symbol, hallmark$gs_name) fgsea_res <- fgsea(pathways, stats = gene_ranks, nperm = 1000, minSize = 15, maxSize = 500) cat("GSEA显著通路数:", nrow(subset(fgsea_res, padj < 0.05)), "\n") 5.6.3 其它富集方法 ~~~~~~~~~~~~~~~~~~ 功能归类是将语义相似的GO terms聚合成更高级别的功能描述，以减少冗余和提高可解释性。REVIGO是一个常用的GO术语精简工具，它基于GO术语之间的语义相似性（semantic similarity）对显著富集的GO terms进行聚类和代表性术语选择。通过REVIGO处理后的GO富集结果更加简洁明了，避免了报告中充斥着大量高度重叠的GO术语的问题。当你得到数百个显著的GO terms需要解读时，REVIGO可以帮助你提炼出核心的功能主题。 通路拓扑学分析超越了简单的基因计数，考虑了基因在通路中的位置和相互作用关系。SPIA（Signaling Pathway Impact Analysis）结合了ORA的证据和通路拓扑信息（如基因在通路中的位置、相互作用类型等），给出综合考虑的通路显著性评分。PathNet将通路视为网络节点，通过分析差异基因在网络中的传播效应来评估通路的重要性。当你希望更深入地理解通路层面的调控机制而非仅仅是基因富集时，通路拓扑学分析方法值得尝试。 网络富集分析将富集结果以网络的形式展现，便于理解基因集之间的关联和重叠。Enrichment Map是Cytoscape的一个插件，它自动构建富集结果网络：节点代表基因集，边的粗细代表基因集之间的重叠程度，节点颜色表示富集方向。这种可视化方式可以直观地发现功能模块和聚类模式。当你需要展示大量基因集之间的关系并识别功能簇时，Enrichment Map是非常强大的可视化工具。 基因集变异分析（GSVA）和单样本GSEA（ssGSEA）将基因集富集的概念扩展到了单个样本层面。GSVA为每个样本的每个基因集计算一个富集分数，从而将基因表达矩阵转换为基因集活性矩阵。ssGSEA是GSVA的一种特定实现方式。这使得我们可以对基因集级别的活动进行差异表达分析、聚类和分类等操作。GSVA在免疫学研究中特别有用，例如计算不同免疫细胞类型的signature activity来推测肿瘤微环境的免疫浸润情况。当你需要从基因集角度描述单个样本的特征时，GSVA/ssGSEA是合适的方法。 富集结果可视化是展示富集分析发现的最后一步。dotplot（点图）以每个基因集为一个点，横轴为富集倍数或NES，纵轴为基因集名称，点的颜色和大小可以编码额外信息（如基因数量、FDR值），是展示ORA或GSEA结果最紧凑和常用的方式。cnetplot（概念网络图）展示了基因集与其成员基因之间的连接关系，可以直观地看到哪些基因驱动了特定基因集的富集。emapplot（富集图谱）展示了基因集之间的重叠关系网络，类似Enrichment Map的简化版。ridgeplot（岭图）展示了基因集中基因的排序分布密度。clusterProfiler统一提供了这些绑图函数，使得富集结果的可视化变得十分便捷。当你准备富集分析结果的图表时，建议组合使用多种可视化方式以全面展示分析发现。 .. code-block:: r # 富集结果可视化 library(enrichplot) # dotplot: GO BP富集结果点图 dotplot(go_bp, showCategory = 20, title = "GO BP富集分析 (Dotplot)") + theme(axis.text.y = element_text(size = 8)) # cnetplot: 基因-通路网络图 cnetplot(go_bp, showCategory = 5, foldChange = gene_ranks, circular = FALSE) # emapplot: 富集项关系网络 emapplot(go_bp, showCategory = 15) # ridgeplot: GSEA基因集分布岭图 ridgeplot(fgsea_res[1:6, ]) 5.7 可变剪接分析 ---------------- 可变剪接（alternative splicing）是真核基因表达调控的重要机制，通过不同的剪接方式从一个基因产生多种mRNA异构体（isoform）。人类基因组中超过95%的多外显子基因经历可变剪接，可变剪接的异常与众多疾病密切相关。RNA-seq数据为可变剪接分析提供了前所未有的机会，使我们能够在全基因组范围内系统地检测和定量剪接事件。 5.7.1 剪接事件类型 ~~~~~~~~~~~~~~~~~~ 外显子跳跃（Exon Skipping，SE）是最常见的可变剪接类型，指某个外显子在成熟mRNA中被跳过不包括在内。在外显子跳跃事件中，上游外显子的供体位点（donor site）直接与下游外显子的受体位点（acceptor site）相连，中间的外显子被完全跳过。外显子跳跃可以改变蛋白质的结构域组成，甚至引入移码突变导致提前终止密码子的产生。当你分析可变剪接数据时，SE事件通常占总可变剪接事件的很大比例，是需要重点关注的类型。 内含子保留（Intron Retention，RI）是指某个内含子未被切除而保留在成熟mRNA中。在哺乳动物中，RI曾经被认为主要是剪接缺陷的副产品，但现在认识到许多RI事件具有重要的调控功能，如通过引入提前终止密码子触发无义介导的mRNA降解（NMD）。RI事件的检测需要特别注意，因为未完全剪接的前体mRNA和基因组DNA污染都可能表现为假的RI信号。当你观察到RI事件时，需要谨慎验证其真实性，排除技术性假阳性的可能。 可变5'剪接位点（Alternative 5' Splice Site，A5SS）和可变3'剪接位点（Alternative 3' Splice Site，A3SS）是指使用不同的供体或受体位点导致外显子边界的变化。A5SS使用不同的上游剪接供体位点，改变了外显子的5'端边界；A3SS使用不同的下游剪接受体位点，改变了外显子的3'端边界。这两种事件通常导致外显子长度的缩短或延长，可能影响蛋白质的功能域或调控元件。A5SS和A3SS的检测精度依赖于剪接位点附近有足够的reads覆盖。 互斥外显子（Mutually Exclusive Exons，MXE）是指两个或多个外显子中只能有一个被包含在成熟的mRNA中，它们互相排斥。MXE事件在蛋白质多样性产生中扮演重要角色，例如在神经元受体和离子通道基因中频繁出现，通过切换不同的外显子来改变蛋白质的性质。MXE事件的检测需要同时考虑多个外显子之间的协调变化模式，比单一类型的剪接事件检测更具挑战性。 可变起始外显子和可变终止外显子是指转录起始于不同的第一个外显子或终止于不同的最后一个外显子。这些事件通常与可变启动子（alternative promoter）使用和多聚腺苷酸化位点（alternative polyadenylation site）选择相关联，它们影响了mRNA的5'端和3'端UTR区域，进而可能影响mRNA的稳定性、翻译效率和亚细胞定位。当你研究基因表达的调控层面时，这些事件提供了超越编码区变化的重要视角。 5.7.2 基于比对的剪接检测 ~~~~~~~~~~~~~~~~~~~~~~~~ 从比对中提取剪接点是检测可变剪接的基础。剪接点（splice junction）是外显子-内含子边界处形成的连接，在RNA-seq比对中以N操作符表示（如在CIGAR字符串中出现的75N500表示跨越了一个500bp的内含子）。STAR、HISAT2等剪接感知比对器会在比对过程中自动识别和记录剪接点信息，输出SJ.out.tab文件包含所有检测到的剪接点及其支撑读段数。从BAM文件中提取剪接信息也可以使用regtools、junction_extract等专用工具。当你开始可变剪接分析时，第一步就是从比对结果中全面地收集剪接点信息。 rMATs（Replicate Multivariate Analysis of Transcript Splicing）是目前最广泛使用的差异剪接分析工具之一。rMATs可以检测五种主要的可变剪接事件类型（SE、RI、A5SS、A3SS、MXE），支持配对和不配对样本设计，使用似然比检验判断剪接事件是否在条件间存在显著差异。rMATs的输出包括每个剪接事件的PSI（Percent Spliced In）值、P-value和FDR。rMATs需要输入BAM文件和基因注释GTF文件，运行时需要较大的内存（通常需要30GB以上的人类基因组分析）。当你需要进行全面的差异可变剪接分析时，rMATs是首选工具之一。 .. code-block:: bash # rMATs差异剪接分析 python rmats.py --b1 control_bams.txt --b2 treat_bams.txt \ --gtf gencode.v44.annotation.gtf -o rmats_output \ --readLength 150 --nthread 16 --libType fr-firststrand # SUPPA2剪接分析 suppa.py generateEvents -i gencode.v44.annotation.gtf -o events -f ioe # leafcutter剪接聚类 leafcutter_cluster.jar -i bam_list.txt -o leafcutter_out -m 50 MAJIQ（Modeling Alternative Junction Inclusion Quantification）采用了一种不同于传统事件分类的方法。MAJIQ首先构建局部剪接图（local splice graph），将一个基因区域内所有检测到的剪接路径建模为图结构，然后量化每条路径的使用变化（表示为ΔPsi，Delta PSI）。MAJIQ的优势在于它可以发现新颖的剪接事件而不局限于预先定义的事件类型，并且提供了可视化工具Voila用于交互式浏览剪接变化。当你希望发现新的剪接事件或研究复杂剪接模式时，MAJIQ提供了更灵活的分析框架。 MISO（Mixture of Isoforms）是基于贝叶斯推断的剪接分析工具。MISO为每个剪接事件建立一个统计模型，使用马尔可夫链蒙特卡洛（MCMC）方法估计剪接比例的后验分布。MISO的优势在于它提供了完整的不确定性量化（置信区间），而不仅仅是点估计，这对于判断剪接变化的统计可靠性很有帮助。MISO主要用于分析预先定义的剪接事件类型，其事件库需要从参考注释中生成。当你关注剪接变化的不确定性估计时，MISO的贝叶斯框架提供了独特价值。 SUPPA2（Speedy Processing of Annotations for Alternative Splicing Analysis）采用了一种轻量级的策略：它不直接操作BAM文件，而是基于转录本定量结果（如Salmon/Kallisto的输出）来计算PSI值和检测差异剪接事件。SUPPA2首先从GTF注释中推导所有可能的剪接事件及其对应的转录本隶属关系，然后根据转录本表达量计算每个事件的PSI值。SUPPA2的速度极快（因为不需要重新处理BAM文件），且与伪比对流程无缝集成。当你已经在使用Salmon/Kallisto进行定量并希望快速添加可变剪接分析时，SUPPA2是最便捷的选择。 leafcutter采用了与众不同的思路：它不依赖于基因注释中预定义的剪接事件，而是直接从比对数据中聚类内含子剪接模式。leafcutter首先提取所有检测到的剪接 junction，然后基于它们的共表达模式进行聚类，每个聚类代表一个潜在的差异剪接区域。leafcutter的优势在于它可以发现注释中不存在的剪接变化，且不受限于五种经典事件类型的框架。leafcutter特别适合于发现新的剪接变异和研究群体水平的剪接多态性。当你想要一种注释无关的（annotation-free）剪接分析方法时，leafcutter提供了独特而有价值的视角。 Whippet是基于图索引的快速可变剪接分析工具。Whippet构建了一个剪接图索引，将基因组上的所有可能剪接路径编码为图结构，然后通过快速图匹配来定量和比较剪接事件。Whippet兼具速度和灵活性，支持差异剪接检测、新事件发现和PSI值计算等功能。Whippet的索引构建是一次性的前期投入，后续的查询和分析都非常快速。当你需要平衡分析速度和检测全面性时，Whippet是一个值得考虑的现代工具。 5.7.3 基于组装的方法 ~~~~~~~~~~~~~~~~~~~~ StringTie和Cufflinks是基于组装的转录本重构工具，它们从RNA-seq比对结果中从头组装转录本结构。StringTie使用网络流算法在剪接图上寻找最优路径来重构转录本，其输出是重构的GTF文件，包含每个样本中检测到的转录本结构及表达量估计。StringTie可以合并多个样本的组装结果（StringTie --merge）以获得更全面的转录本注释。Cufflinks是较早的组装工具（现已不再积极维护），其配套的Cuffdiff工具曾广泛用于差异转录本表达分析。当你需要发现新的转录本或完善现有的基因注释时，StringTie组装是推荐的方法。 比较转录本表达量是组装后分析的核心任务。Cuffdiff是Cufflinks套件的一部分，专门用于基于组装转录本的差异表达分析。Ballgown是Bioconductor中专门处理转录本级别表达数据的R包，它可以导入StringTie或Cufflinks的输出，进行转录本级别的差异表达检验和可视化。与基因级别的差异表达分析相比，转录本级别的分析面临更大的挑战：转录本之间的共享reads使得定量更加不确定，需要更大的样本量才能达到相当的统计功效。当你对特定基因的不同isoform感兴趣时，基于组装的方法提供了直接的定量途径。 使用长读长RNA-seq改良剪接检测是近年来的重要发展方向。长读长（如PacBio Iso-Seq和Oxford Nanopore direct RNA-seq）可以跨越完整的转录本，直接观测到剪接连接的模式，无需像短读长那样进行推断。TALON（Transcript Annotation from Long read ONtologies）是一个专为长读长RNA-seq设计的转录本鉴定和定量工具，它可以处理Iso-Seq和Nanopore数据，并与短读长的注释进行整合。FLAIR（Full-Length Alternative Isoform analysis of RNA）专注于从长读长数据中发现和定量可变剪接异构体。当你有机会使用长读长RNA-seq数据时，这些工具可以大幅提高剪接检测的准确性和全面性。 5.7.4 剪接定量指标 ~~~~~~~~~~~~~~~~~~ 百分比剪接指数（PSI，Ψ，读作psi）是描述剪接事件中使用某一剪接方式的转录本比例的核心指标。对于一个典型的外显子跳跃事件，PSI = inclusion reads / (inclusion reads + skipping reads)，表示包含该外显子的转录本占总相关转录本的比例。PSI值的范围为0到1，PSI=1表示该外显子总是被包含（完全inclusion），PSI=0表示该外显子总是被跳跃（完全skipping），PSI=0.5表示两种剪接方式各占一半。PSI值在不同条件之间的变化（ΔPSI）是衡量差异剪接程度的直观指标。当你报告可变剪接结果时，PSI值和ΔPSI是最重要的定量指标之一。 差异剪接统计检验用于判断观察到的PSI变化是否具有统计学意义。不同的工具采用不同的统计框架：rMATs使用似然比检验，MISO使用贝叶斯后验概率，MAJIQ使用后验置信区间，leafcutter使用似然 ratio 检验。无论哪种方法，都需要考虑PSI估计的不确定性，特别是对于低表达的剪接事件，PSI的估计方差较大，统计检验的功效较低。当你解释差异剪接结果时，不仅要看p值或FDR，还要关注ΔPSI的大小和PSI估计的置信区间宽度。 差异剪接可视化是传达剪接分析发现的重要手段。Sashimi plot（也称junction plot）是RNA-seq数据可视化的标志性图像，它展示了基因组区域上reads的覆盖度和剪接连接的弧线，弧线的粗细表示跨越该剪接点的reads数量。Sashimi plot可以直观地展示不同样本或条件之间的剪接模式差异。IGV（Integrative Genomics Viewer）提供了交互式的剪接点视图，可以逐个检查剪接点的支撑reads和序列上下文。ggsashimi和sgviz等R包可以编程生成高质量的Sashimi plot。当你需要展示具体的剪接事件时，Sashimi plot是最具说服力的可视化方式。 剪接事件功能注释帮助理解可变剪接变化的生物学后果。并非所有的剪接变化都具有功能意义：有些变化发生在UTR区域可能不影响蛋白质序列，有些变化可能只是噪音。最重要的功能影响之一是剪接变化是否引入了提前终止密码子（premature termination codon，PTC），PTC可能触发无义介导的mRNA降解（NMD），从而导致该转录本的有效表达量下降。其他功能性影响包括：是否改变了蛋白质的功能域、是否影响了磷酸化位点或蛋白互作界面、是否移除了关键的调控元件等。当你对差异剪接结果进行优先级排序时，功能影响程度应该是重要的考量因素。 利用可变剪接数据库可以为你的分析结果提供注释和参考。VastDB是一个专门的可变剪接数据库，收录了人和小鼠的大量组织和发育阶段的可变剪接事件及其PSI值。ASevent数据库收集了来自多个研究的可变剪接事件及其与疾病的关联。SpliceDB提供了剪接位点和剪接因子的综合信息。当你检测到一个新的差异剪接事件时，查询这些数据库可以了解该事件是否已被报道、在哪些组织中活跃、以及是否有已知的疾病关联。 .. code-block:: r # PSI计算与差异剪接可视化 library(ggplot2) # PSI计算函数 calc_psi <- function(inclusion_reads, skipping_reads) { total <- inclusion_reads + skipping_reads psi <- inclusion_reads / total return(list(PSI = psi, CI_lower = binom.test(inclusion_reads, total)$conf.int[1], CI_upper = binom.test(inclusion_reads, total)$conf.int[2])) } # 模拟外显子跳跃事件数据 events <- data.frame( Event = paste0("SE", 1:6), Ctrl_incl = c(85, 52, 30, 45, 70, 15), Ctrl_skip = c(15, 48, 70, 55, 30, 85), Treat_incl = c(20, 80, 35, 48, 25, 75), Treat_skip = c(80, 20, 65, 52, 75, 25) ) events$Ctrl_PSI <- sapply(1:nrow(events), function(i) calc_psi(events$Ctrl_incl[i], events$Ctrl_skip[i])$PSI) events$Treat_PSI <- sapply(1:nrow(events), function(i) calc_psi(events$Treat_incl[i], events$Treat_skip[i])$PSI) events$dPSI <- events$Treat_PSI - events$Ctrl_PSI print("外显子跳跃事件PSI分析:") print(events[, c("Event", "Ctrl_PSI", "Treat_PSI", "dPSI")]) sig_as <- subset(events, abs(dPSI) > 0.2) cat("\n显著差异剪接事件数 (|dPSI|>0.2):", nrow(sig_as), "\n") 5.8 非编码RNA分析 ----------------- 非编码RNA（non-coding RNA，ncRNA）是指不编码蛋白质但在细胞中具有重要调控功能的RNA分子。随着高通量测序技术的发展，人们发现基因组中只有不到2%的序列编码蛋白质，而大部分转录产物是非编码RNA。非编码RNA包括长链非编码RNA（lncRNA）、微小RNA（miRNA）、环状RNA（circRNA）、piRNA、snoRNA等多种类型，它们在基因表达调控、染色质重塑、RNA加工等过程中发挥关键作用。 5.8.1 长链非编码RNA（lncRNA） ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 长链非编码RNA（long non-coding RNA，lncRNA）是指长度超过200个核苷酸且不具备蛋白质编码能力的RNA分子。lncRNA的鉴别需要满足几个标准：首先，转录本长度必须大于200nt，这是区分lncRNA和小非编码RNA的基本界限；其次，开放阅读框（ORF）较短，通常小于100个氨基酸，且缺乏完整的起始密码子和终止密码子；第三，编码潜力低，即不具备翻译成功能蛋白质的能力；第四，外显子数量通常较少，多数lncRNA只有1-2个外显子。当你从RNA-seq数据中鉴定新的lncRNA时，需要综合运用这些标准进行筛选。 编码潜力预测是lncRNA鉴定的核心步骤。CPC2（Coding Potential Calculator 2）是基于序列特征和机器学习的编码潜力预测工具，它使用支持向量机模型根据ORF长度、等电点、分子量等特征判断转录本是否具有编码能力。CPAT（Coding-Potential Assessment Tool）使用逻辑回归模型，基于ORF长度、Fickett得分和六聚体得分进行预测，速度较快。PLEK（Predictor of Long non-coding RNAs and messenger RNAs based on k-mer scheme）基于k-mer频率特征使用支持向量机进行分类。LGC（Long non-coding RNA Gene Calculator）通过计算GC含量和ORF长度等简单特征进行快速筛选。当你需要从大量新转录本中筛选lncRNA候选时，建议组合使用多个工具取交集以提高预测准确性。 lncRNA的定量方法与mRNA类似，但需要特异性的注释文件。由于lncRNA的注释在标准基因注释（如GENCODE）中可能不完整，你可能需要整合来自NONCODE、LNCipedia等专门数据库的lncRNA注释。定量工具如featureCounts、HTSeq-count、Salmon等同样适用于lncRNA，只需提供包含lncRNA的GTF注释文件。当你进行lncRNA定量时，需要注意lncRNA的表达量通常低于mRNA，因此可能需要更深的测序深度或更宽松的过滤阈值。 lncRNA的差异表达分析与mRNA相同，可以使用DESeq2、edgeR或limma-voom等工具。由于lncRNA表达量普遍较低，在离散度估计和统计检验时需要特别注意。建议使用独立过滤时适当放宽阈值，避免将大量低表达的lncRNA过早排除。当你发现lncRNA差异表达后，后续的功能注释和验证是更大的挑战，因为大多数lncRNA的功能尚不明确。 lncRNA与其他RNA的共表达网络分析是推断lncRNA功能的重要策略。基于"guilt-by-association"原则，如果一个lncRNA与已知功能的mRNA高度共表达，那么该lncRNA可能参与相似的生物学过程。WGCNA（Weighted Gene Co-expression Network Analysis）是构建共表达网络的常用工具，它将基因按表达模式聚类成模块，然后分析模块与表型的关联。当你研究lncRNA的功能时，共表达网络分析可以提供重要的线索和假设。 lncRNA的亚细胞定位对其功能模式有重要影响。核内lncRNA通常参与染色质调控、转录控制和RNA加工，而细胞质lncRNA可能参与mRNA稳定性调控、翻译控制或作为miRNA海绵。lncLocator和iLoc-lncRNA是基于序列特征预测lncRNA亚细胞定位的工具。当你研究一个新的lncRNA时，预测其亚细胞定位可以帮助推断其可能的作用机制。 lncRNA数据库为研究提供了重要的参考资源。NONCODE是目前最全面的lncRNA数据库之一，收录了多个物种的lncRNA注释和表达信息。LNCipedia提供人类lncRNA的序列和注释信息，并支持批量查询和下载。lncRNAdb收集了文献中功能已验证的lncRNA及其生物学功能描述。LNCat是一个整合多数据源的lncRNA注释平台，提供lncRNA的序列、结构、表达和功能信息。当你需要注释新发现的lncRNA或查找特定lncRNA的信息时，这些数据库是重要的起点。 .. code-block:: bash # lncRNA编码潜力预测 (CPC2) CPC2.py -i transcripts.fasta -o cpc2_output.txt # CPAT编码潜力预测 cpat.py -g transcripts.fasta -d human_Hexamer.tsv -x human_logitModel.RData \ -o cpat_output.txt --min-orf 50 --top-orf 5 # lncRNA定量 (使用Salmon) salmon quant -i lncrna_index -l A -1 R1.fq.gz -2 R2.fq.gz \ -p 8 --validateMappings -o lncrna_quant # WGCNA共表达网络分析 # R代码见下文 .. code-block:: r # WGCNA共表达网络分析 library(WGCNA) datExpr <- t(normalized_counts) powers <- c(1:20) sft <- pickSoftThreshold(datExpr, powerVector = powers) softPower <- sft$powerEstimate adjacency <- adjacency(datExpr, power = softPower) TOM <- TOMsimilarity(adjacency) dissTOM <- 1 - TOM geneTree <- hclust(as.dist(dissTOM), method = "average") dynamicMods <- cutreeDynamic(dendro = geneTree, minClusterSize = 30) moduleColors <- labels2colors(dynamicMods) cat("识别到的模块数:", length(unique(moduleColors)), "\n") 5.8.2 微小RNA（miRNA） ~~~~~~~~~~~~~~~~~~~~~~~ 微小RNA（microRNA，miRNA）是一类长度约22个核苷酸的小非编码RNA，通过碱基配对识别靶mRNA的3'UTR区域，导致mRNA降解或翻译抑制。miRNA在基因表达调控中扮演着核心角色，据估计人类基因组中编码了超过2000个miRNA，调控着超过60%的蛋白编码基因。 miRNA测序文库制备需要特殊的实验流程。由于miRNA分子短小，需要使用特定的small RNA建库试剂盒，通过连接接头而非随机引物反转录来捕获小RNA分子。文库制备过程中需要使用特定大小的PAGE胶或磁珠筛选18-30nt的片段，去除接头二聚体和其他大小的RNA。当你进行miRNA测序时，需要注意起始RNA的质量和数量要求，以及文库制备过程中可能产生的接头二聚体污染。 接头修剪是miRNA-seq数据处理的关键步骤。由于miRNA序列很短，测序读长通常覆盖整个miRNA加上部分接头序列，因此必须准确去除接头。cutadapt是最常用的接头修剪工具，它可以精确识别并去除3'端的接头序列。fastx_clipper是FASTX-Toolkit的一部分，也可以用于接头修剪。当你处理miRNA-seq数据时，需要确保接头序列正确，并设置适当的最小保留长度（通常18nt）以避免过度修剪。 比对到miRBase是miRNA定量和注释的标准流程。miRBase是最权威的miRNA数据库，收录了经过实验验证的miRNA序列及其注释。比对工具可以使用Bowtie或STAR，设置允许错配数通常为0或1，因为miRNA序列很短且高度保守。当你进行miRNA分析时，建议使用最新版本的miRBase以获得最完整的miRNA注释。 miRNA定量与mRNA类似，可以使用featureCounts从比对后的BAM文件中计数。需要注意的是，由于miRNA序列短且可能存在多个miRNA共享相似序列的情况，多映射reads的处理需要特别关注。常用的策略包括：只保留唯一映射的reads、按比例分配多映射reads、或使用专门的miRNA定量工具如miRDeep2。当你进行miRNA定量时，建议同时统计每个miRNA的reads数和归一化后的表达量（如RPM，reads per million）。 miRNA的差异表达分析与mRNA类似，可以使用DESeq2、edgeR等工具。由于miRNA数量相对较少（通常几百到几千个），多重检验校正的负担较小，统计功效相对较高。当你发现差异表达的miRNA后，下一步通常是预测其靶基因并分析其可能的功能影响。 miRNA靶基因预测是理解miRNA功能的关键步骤。TargetScan是最广泛使用的靶基因预测工具之一，它基于种子序列匹配和进化保守性进行预测，并提供预测的上下文得分。miRanda考虑了miRNA-mRNA双链的热力学稳定性和序列互补性。miRDB使用机器学习方法基于序列特征进行预测。RNAhybrid基于RNA二级结构预测最小自由能来评估结合可能性。当你进行靶基因预测时，建议使用多个工具取交集以提高预测可靠性，并结合表达数据进行验证。 miRNA-mRNA整合分析是揭示miRNA调控机制的重要策略。如果miRNA上调而其预测靶基因下调，这支持miRNA对靶基因的负向调控关系。更严谨的分析可以使用相关性分析或因果推断方法来验证miRNA-mRNA的调控关系。当你研究特定生物学过程中的miRNA调控时，整合miRNA和mRNA表达数据可以构建更完整的调控网络。 miRNA功能富集分析帮助理解差异miRNA的生物学意义。miRPath和DIANA-mirPath是专门用于miRNA富集分析的工具，它们基于miRNA的靶基因进行通路富集分析。分析结果可以揭示差异miRNA可能影响的生物学通路和功能类别。当你完成miRNA差异表达分析后，功能富集分析是解读其生物学意义的重要步骤。 .. code-block:: bash # miRNA-seq数据处理流程 # 1. 接头修剪 cutadapt -a TGGAATTCTCGGGTGCCAAGG -m 18 -o trimmed.fq.gz raw.fq.gz # 2. 比对到miRBase bowtie -v 1 -m 1 --best --strata mirbase_index trimmed.fq.gz \ -S mirna_aligned.sam # 3. 定量 (featureCounts) featureCounts -a mirbase.gtf -t miRNA -o mirna_counts.txt mirna_aligned.bam # 4. miRDeep2新miRNA发现 miRDeep2.pl trimmed.fq.gz genome.fa mirna_aligned.arf \ mature_miRNAs.fa hairpin.fa > mirdeep2_report.txt .. code-block:: r # miRNA靶基因预测与整合分析 library(multiMiR) miRNA_list <- c("hsa-miR-21-5p", "hsa-miR-155-5p", "hsa-miR-122-5p") targets <- get_multimir(mirna = miRNA_list, target = TRUE, summary = TRUE, table = "validated") cat("验证的miRNA-靶基因对数:", nrow(targets@validated), "\n") # miRNA-mRNA相关性分析 cor_matrix <- cor(t(miRNA_expr), t(mRNA_expr), method = "spearman") significant_pairs <- which(abs(cor_matrix) > 0.5, arr.ind = TRUE) cat("显著相关的miRNA-mRNA对数:", nrow(significant_pairs), "\n") 5.8.3 其它非编码RNA ~~~~~~~~~~~~~~~~~~~ 环状RNA（circular RNA，circRNA）是一类通过反向剪接形成的共价闭合环状RNA分子。circRNA具有高度稳定性，不易被RNA外切酶降解，在基因表达调控中可能作为miRNA海绵发挥作用。CIRI2是基于成对嵌合读段检测circRNA的工具，它扫描BAM文件中的嵌合比对信号来识别反向剪接位点。circExplorer通过分析剪接点信息来识别circRNA。find_circ是最早开发的circRNA检测工具之一，通过搜索未比对reads中的反向剪接信号来发现circRNA。当你研究circRNA时，需要注意使用RNase R处理样本可以富集circRNA并去除线性RNA。 piRNA（PIWI-interacting RNA）是一类长度约24-32nt的小RNA，主要在生殖细胞中表达，与PIWI蛋白结合参与转座子沉默和基因组稳定性维持。piRNA的注释和定量可以使用small RNA-seq数据，通过比对到piRBase等专门的piRNA数据库进行识别。当你研究生殖系统或肿瘤中的piRNA时，这些数据库提供了重要的参考资源。 snoRNA（small nucleolar RNA）是一类指导rRNA和snRNA化学修饰的小RNA，主要位于核仁中。snoRNA可以从small RNA-seq或常规RNA-seq数据中检测，使用RNAcentral等综合数据库进行注释。snoRNA的异常表达可能与某些疾病相关，特别是与核糖体功能相关的疾病。 tRNA衍生片段（tRNA-derived fragments，tRFs）和rRNA衍生片段（rRNA-derived fragments，rRFs）是近年来发现的新型小RNA分子。它们来源于tRNA或rRNA的特异性切割，可能具有基因表达调控功能。这些片段的分析需要专门的工具和数据库，如tRFdb和MINTbase。当你研究应激反应或癌症等条件下的小RNA组学时，这些新型小RNA值得关注。 .. code-block:: bash # circRNA检测 (CIRI2) CIRI2.pl -I aligned.bam -O circRNA_output.txt -F reference.fa # circExplorer检测 circExplorer.py parse -t STAR aligned.bam > circExplorer_output.txt # find_circ检测 find_circ.py -p ref_name -s sample_name unmapped.bam find_circ_output/ 5.9 单细胞转录组学 ------------------ 单细胞RNA测序（single-cell RNA-seq，scRNA-seq）技术使我们能够在单个细胞的分辨率下研究基因表达，揭示细胞群体的异质性和稀有细胞类型。与bulk RNA-seq测量细胞群体的平均表达不同，scRNA-seq捕获每个细胞的转录组信息，为理解发育、疾病和免疫等复杂生物学过程提供了前所未有的视角。 5.9.1 单细胞RNA-seq预处理 ~~~~~~~~~~~~~~~~~~~~~~~~~ 测序数据拆分是scRNA-seq数据处理的第一步。10X Genomics等平台使用细胞条形码（cell barcode）标记每个细胞的reads，需要根据barcode序列将混合的测序数据拆分到各个细胞。Cell Ranger是10X Genomics官方提供的处理流程，它自动完成barcode识别、比对和定量。对于其他平台，可能需要使用自定义脚本或UMI-tools等工具进行barcode处理。当你处理scRNA-seq数据时，首先需要了解所使用平台的数据格式和barcode设计。 比对到参考基因组可以使用STAR或Cell Ranger封装的STAR流程。由于scRNA-seq数据通常使用UMI（Unique Molecular Identifier）标记原始分子，比对时需要保留UMI信息用于后续去重。Cell Ranger的比对流程针对scRNA-seq数据进行了优化，包括对剪接位点的特殊处理和对多映射reads的策略。当你使用非10X平台的数据时，可能需要手动配置STAR参数以适应scRNA-seq的特点。 基因表达定量在scRNA-seq中涉及UMI计数去重。由于PCR扩增会产生重复的reads，UMI标记可以识别来自同一原始分子的reads并去除PCR重复。featureCounts、HTSeq-count等工具可以用于scRNA-seq的基因计数，但需要设置适当的参数处理UMI。Cell Ranger流程自动完成UMI去重和计数。当你进行定量时，需要注意区分UMI计数和read计数，UMI计数更能反映真实的分子数量。 细胞-基因计数矩阵构建是scRNA-seq分析的起点。矩阵的行代表基因，列代表细胞，每个元素是该基因在该细胞中的UMI计数。这个矩阵通常非常稀疏，因为每个细胞只表达基因组中的一部分基因，且scRNA-seq的检测深度有限。当你构建计数矩阵后，需要进行质量控制过滤掉低质量细胞和可能的空液滴。 .. code-block:: bash # Cell Ranger处理流程 (10X Genomics) cellranger mkfastq --run=/path/to/bcl --output-dir=fastq_output cellranger count --id=sample1 --transcriptome=/path/to/refdata \ --fastqs=/path/to/fastq_output --sample=sample1 --expect-cells=5000 5.9.2 质量控制 ~~~~~~~~~~~~~~ 每个细胞的UMI总数（nCount_RNA）是衡量细胞质量的基本指标。高质量细胞应该有适中的UMI总数，过低可能表示细胞裂解不完全或RNA丢失，过高可能是双细胞或多细胞。典型的哺乳动物单细胞UMI总数在1000到50000之间，具体取决于测序深度和细胞类型。当你设置过滤阈值时，需要根据数据分布和实验设计确定合理的范围。 每个细胞检测到的基因数量（nFeature_RNA）是另一个重要指标。基因数量与UMI数量相关，但可以独立反映细胞的转录活性。低质量细胞通常检测到的基因数量很少，而高活性细胞可能检测到数千个基因。当你进行质量控制时，可以绘制nFeature_RNA与nCount_RNA的散点图来评估数据质量。 线粒体基因比例是识别低质量细胞的关键指标。受损或凋亡细胞会释放细胞质RNA，而线粒体RNA由于位于线粒体内而相对保留，导致线粒体基因比例升高。通常，线粒体基因比例超过10-20%的细胞被认为是低质量细胞需要过滤。当你计算线粒体比例时，需要识别所有线粒体编码的基因（通常以MT-或mt-开头）。 核糖体基因比例是另一个可选的质量指标。核糖体基因在某些细胞类型中高表达，但异常高的比例可能表示细胞状态异常或技术问题。当你研究特定细胞类型时，需要了解其正常的核糖体基因表达水平。 细胞周期评分可能混淆细胞类型鉴定，因为处于不同细胞周期阶段的细胞可能聚类到不同的簇。Seurat提供了细胞周期评分功能，基于已知的细胞周期相关基因表达计算S期和G2M期得分。当你发现聚类结果可能受细胞周期影响时，可以考虑在归一化时回归掉细胞周期效应。 双重或多重细胞的检测是scRNA-seq质量控制的重要环节。在微液滴平台中，有时两个或多个细胞会被封装到同一个液滴中，表现为一个"细胞"具有异常高的UMI数和基因数。DoubletFinder和Scrublet是检测双细胞的常用工具，它们基于模拟双细胞的特征来识别可能的doublets。当你处理高细胞密度的样本时，双细胞检测尤为重要。 过滤死细胞与空液滴可以使用EmptyDrops方法。EmptyDrops通过比较每个液滴的表达谱与背景空液滴的表达谱，识别真正含有细胞的液滴。这对于区分低表达细胞和空液滴特别有用，避免错误地过滤掉真实的低表达细胞。当你使用10X Genomics数据时，Cell Ranger已经进行了初步的细胞识别，但你可能需要进一步调整阈值。 .. code-block:: r # Seurat单细胞质量控制 library(Seurat) pbmc <- Read10X(data.dir = "filtered_feature_bc_matrix/") seurat_obj <- CreateSeuratObject(counts = pbmc, project = "PBMC") seurat_obj[["percent.mt"]] <- PercentageFeatureSet(seurat_obj, pattern = "^MT-") VlnPlot(seurat_obj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3) seurat_obj <- subset(seurat_obj, subset = nFeature_RNA > 200 & nFeature_RNA < 7500 & percent.mt < 20) cat("过滤后细胞数:", ncol(seurat_obj), "\n") # 双细胞检测 (DoubletFinder) library(DoubletFinder) seurat_obj <- SCTransform(seurat_obj) seurat_obj <- RunPCA(seurat_obj) seurat_obj <- FindNeighbors(seurat_obj, dims = 1:30) sweep.res.list <- paramSweep_v3(seurat_obj, PCs = 1:30, sct = TRUE) sweep.stats <- summarizeSweep(sweep.res.list, GT = FALSE) bcmvn <- find.pK(sweep.stats) optimal_pK <- bcmvn$pK[which.max(bcmvn$BCmetric)] 5.9.3 归一化与批次校正 ~~~~~~~~~~~~~~~~~~~~~~ 平滑分位数归一化（sctransform）是Seurat推荐的归一化方法。sctransform使用正则化负二项回归对UMI计数进行建模，去除测序深度的影响，同时保留生物学变异。与简单的对数归一化相比，sctransform可以更好地处理技术噪音，特别是对于低表达基因。当你使用Seurat进行scRNA-seq分析时，SCTransform函数是推荐的归一化起点。 基于对数归一化（LogNormalize）是更简单的方法，将每个细胞的计数除以该细胞的总计数，乘以缩放因子（通常为10000），然后取自然对数。这种方法计算简单但可能引入偏倚，特别是对于低计数的基因。当你需要快速处理数据或与历史数据比较时，LogNormalize是一个可行的选择。 基于正则化负二项回归的SCnorm是另一种归一化方法，它考虑了基因特异的测序深度依赖性，对不同表达水平的基因使用不同的归一化策略。当你发现数据存在复杂的测序深度效应时，SCnorm可能提供更准确的归一化。 批次效应去除是整合多个样本或批次数据的关键步骤。Harmony是一种快速的批次校正方法，它在PCA空间中迭代地调整细胞坐标以去除批次效应。Seurat的Integration方法使用典型相关分析（CCA）识别共享的细胞状态，然后通过锚点（anchor）对齐不同批次的细胞。fastMNN基于相互最近邻（mutual nearest neighbors）进行批次校正。BBKNN（Batch Balanced k-Nearest Neighbors）在构建邻居图时考虑批次信息。当你整合多个批次的数据时，需要选择适合你数据特点的批次校正方法。 数据缩放（ScaleData）是归一化后的标准步骤，它将每个基因的表达值中心化为均值为0，并可选地标准化为方差为1。缩放确保所有基因在下游分析（如PCA）中具有相同的权重。当你使用Seurat时，ScaleData函数还可以回归掉不需要的变异来源（如线粒体比例、细胞周期得分）。 .. code-block:: r # Seurat归一化与批次校正 seurat_obj <- SCTransform(seurat_obj, vars.to.regress = "percent.mt") # Harmony批次校正 library(harmony) seurat_obj <- RunPCA(seurat_obj) seurat_obj <- RunHarmony(seurat_obj, group.by.vars = "batch") # Seurat Integration批次整合 pbmc.list <- SplitObject(seurat_obj, split.by = "batch") pbmc.list <- lapply(pbmc.list, SCTransform) features <- SelectIntegrationFeatures(pbmc.list, nfeatures = 3000) pbmc.anchors <- FindIntegrationAnchors(pbmc.list, anchor.features = features) pbmc.integrated <- IntegrateData(pbmc.anchors) 5.9.4 降维与可视化 ~~~~~~~~~~~~~~~~~~ 主成分分析（PCA）是scRNA-seq降维的标准方法。PCA将高维基因表达矩阵投影到低维主成分空间，保留最大的方差。在scRNA-seq中，通常使用高度可变基因（highly variable genes）进行PCA，以减少计算量并突出生物学变异。当你进行PCA时，需要选择适当数量的主成分用于后续分析。 t-SNE（t-distributed Stochastic Neighbor Embedding）是一种非线性降维方法，特别适合可视化高维数据的局部结构。t-SNE保持局部邻域关系，使得相似的细胞在二维空间中聚集在一起。然而，t-SNE的全局结构不可靠，不同簇之间的距离没有意义。当你使用t-SNE可视化时，需要记住它主要用于展示而非定量分析。 UMAP（Uniform Manifold Approximation and Projection）是另一种非线性降维方法，近年来在scRNA-seq中广泛使用。UMAP比t-SNE更好地保留了全局结构，运行速度更快，且参数更直观。UMAP已经成为scRNA-seq可视化的首选方法。当你比较不同样本或条件时，UMAP图可以直观地展示细胞群体的分布。 扩散图（Diffusion Map）是一种保留数据流形结构的降维方法，特别适合分析连续的生物学过程（如细胞分化）。扩散图基于扩散距离，可以更好地捕捉细胞的连续状态变化。当你研究发育轨迹或细胞分化时，扩散图是有价值的可视化工具。 选择主成分数量是降维的关键决策。肘部图（ElbowPlot）展示每个主成分解释的方差，通常在某个点后方差下降趋于平缓，该点可以作为主成分数量的选择依据。JackStraw是一种统计检验方法，评估每个主成分是否显著高于随机期望。当你不确定应该使用多少主成分时，建议尝试不同的数量并评估结果的稳定性。 .. code-block:: r # Seurat降维与可视化 seurat_obj <- RunPCA(seurat_obj, npcs = 50) ElbowPlot(seurat_obj, ndims = 50) seurat_obj <- RunUMAP(seurat_obj, dims = 1:30) seurat_obj <- RunTSNE(seurat_obj, dims = 1:30) DimPlot(seurat_obj, reduction = "umap", label = TRUE) 5.9.5 聚类与细胞类型注释 ~~~~~~~~~~~~~~~~~~~~~~~~ 图聚类是scRNA-seq细胞聚类的主流方法。Louvain算法是一种社区检测算法，通过优化模块度来识别细胞群体。Leiden算法是Louvain的改进版本，保证识别的社区是连通的，且运行速度更快。Seurat的FindClusters函数默认使用Louvain或Leiden算法。当你进行聚类时，需要选择适当的分辨率参数，分辨率越高识别的簇越多。 聚类分辨率选择影响最终识别的细胞群体数量。较低的分辨率（如0.4-0.6）识别较少的大簇，适合粗粒度的细胞类型划分；较高的分辨率（如0.8-1.2）识别更多的小簇，可以区分更细的细胞亚型。当你不确定合适的分辨率时，建议尝试多个分辨率并比较结果的生物学合理性。 标记基因识别是定义细胞类型的基础。Seurat的FindAllMarkers函数可以识别每个簇相对于其他簇的差异表达基因，使用Wilcoxon秩和检验、MAST或其他统计方法。标记基因应该是该簇特异性高表达的基因，用于定义细胞身份。当你识别标记基因后，需要结合生物学知识解释每个簇的细胞类型。 细胞类型自动注释工具可以加速细胞类型鉴定过程。SingleR基于参考数据集的基因表达模式自动注释细胞类型。CellAssign使用概率模型将细胞分配到预定义的细胞类型。ScType整合了多个细胞标记数据库进行注释。当你有合适的参考数据集时，自动注释可以提供快速且相对准确的细胞类型预测。 人工注释基于已知标记基因是细胞类型鉴定的金标准。PanglaoDB和CellMarker数据库收录了大量细胞类型的标记基因。你可以使用小提琴图、特征图（FeaturePlot）或点图（DotPlot）可视化标记基因在不同簇中的表达，从而判断每个簇的细胞类型。当你进行人工注释时，需要综合多个标记基因的信息，避免仅依赖单一基因。 差异表达基因列表可视化帮助展示细胞类型的特征。热图（DoHeatmap）可以展示多个标记基因在多个簇中的表达模式。小提琴图（VlnPlot）展示单个基因在单个或多个簇中的表达分布。点图（DotPlot）同时展示基因表达水平和表达细胞比例。特征图（FeaturePlot）在UMAP或t-SNE图上展示基因表达的空间分布。当你准备发表scRNA-seq结果时，这些可视化图表是展示细胞类型特征的重要手段。 .. code-block:: r # Seurat聚类与标记基因识别 seurat_obj <- FindNeighbors(seurat_obj, dims = 1:30) seurat_obj <- FindClusters(seurat_obj, resolution = 0.8) markers <- FindAllMarkers(seurat_obj, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25) top5_markers <- markers %>% group_by(cluster) %>% top_n(n = 5, wt = avg_log2FC) # 细胞类型自动注释 (SingleR) library(SingleR) library(celldex) ref <- HumanPrimaryCellAtlasData() predictions <- SingleR(test = seurat_obj@assays$RNA@data, ref = ref, labels = ref$label.main) seurat_obj$SingleR_label <- predictions$labels # 可视化 VlnPlot(seurat_obj, features = c("CD3D", "CD14", "MS4A1"), ncol = 3) FeaturePlot(seurat_obj, features = c("CD3D", "CD14"), cols = c("grey", "blue")) DotPlot(seurat_obj, features = top5_markers$gene) + RotatedAxis() 5.9.6 轨迹推断与伪时间分析 ~~~~~~~~~~~~~~~~~~~~~~~~~~ 单细胞轨迹推断旨在重建细胞的连续状态变化，如发育分化、应激响应或疾病进展。Monocle3是最广泛使用的轨迹分析工具之一，它使用UMAP或DDRTree降维后构建最小生成树来推断细胞轨迹。Slingshot是一种更灵活的轨迹推断方法，它可以在任何降维结果上拟合平滑曲线。scVelo实现了RNA velocity分析，利用未剪接和剪接mRNA的比例推断细胞的未来状态。当你研究细胞分化或发育过程时，轨迹分析可以揭示细胞状态的连续变化。 分支点与细胞命运决策是轨迹分析的核心发现。在细胞分化过程中，细胞可能到达决策点后向不同方向分化，形成轨迹的分支。识别这些分支点可以帮助理解细胞命运决定的机制。当你分析轨迹时，需要识别分支点并分析不同分支的细胞特征。 差异动力学分析使用tradeSeq等工具识别沿轨迹变化的基因。这些基因可能在特定分支点改变表达模式，或在伪时间上呈现特定的变化趋势。当你发现轨迹上的差异表达基因后，可以进行功能富集分析理解其生物学意义。 推断分化路径需要定义轨迹的起点和终点。起点通常是已知的祖细胞或干细胞，终点是分化的细胞类型。你可以基于生物学知识或标记基因表达来定义起点和终点。当你定义好轨迹后，可以计算每个细胞的伪时间值，表示其在分化路径上的进度。 基因表达随伪时间的变化可以揭示分化过程中的调控机制。你可以可视化单个基因沿伪时间的表达变化，或聚类具有相似变化模式的基因。当你发现沿伪时间变化的基因后，可以进一步分析其调控网络和功能。 .. code-block:: r # Monocle3轨迹分析 library(monocle3) cds <- as.cell_data_set(seurat_obj) cds <- learn_graph(cds) cds <- order_cells(cds, root_cells = root_cluster_cells) plot_cells(cds, color_cells_by = "pseudotime") # 沿伪时间的差异表达基因 graph_test_res <- graph_test(cds, neighbor_graph = "principal_graph") sig_genes <- subset(graph_test_res, q_value < 0.05) cat("沿伪时间变化的基因数:", nrow(sig_genes), "\n") 5.9.7 细胞-细胞相互作用 ~~~~~~~~~~~~~~~~~~~~~~~ 配体-受体对分析是推断细胞间通讯的主要方法。CellPhoneDB是一个常用的工具，它基于已知配体-受体对的数据库，通过统计检验识别细胞类型之间显著富集的相互作用。CellChat提供了更全面的分析框架，包括信号流可视化和通路分析。NicheNet专注于预测配体对靶基因的调控效应。当你研究肿瘤微环境或免疫细胞相互作用时，细胞间通讯分析可以揭示重要的调控机制。 空间临近性推断对于非空间转录组数据是必要的。由于缺乏真实的空间信息，这些工具基于细胞类型的共现模式和相互作用潜力推断可能的空间关系。当你使用非空间数据时，需要谨慎解读推断的空间关系。 通讯网络可视化帮助理解复杂的细胞间相互作用。CellChat提供了多种可视化方式，包括信号流图、网络图和和弦图。这些可视化可以直观地展示哪些细胞类型是主要的信号发送者或接收者。当你报告细胞间通讯结果时，网络可视化是有效的展示方式。 .. code-block:: r # CellChat细胞间通讯分析 library(CellChat) cellchat <- createCellChat(object = seurat_obj) CellChatDB <- CellChatDB.human cellchat <- subsetData(cellchat) cellchat <- identifyOverExpressedGenes(cellchat) cellchat <- identifyOverExpressedInteractions(cellchat) cellchat <- computeCommunProb(cellchat) cellchat <- filterCommunication(cellchat, min.cells = 10) cellchat <- computeCommunProbPathway(cellchat) cellchat <- aggregateNet(cellchat) netVisual_circle(cellchat@net$count, vertex.weight = groupSize) 5.9.8 单细胞多组学整合 ~~~~~~~~~~~~~~~~~~~~~~ CITE-seq技术结合了单细胞RNA测序和蛋白质标记检测，使用DNA条形码标记的抗体同时测量细胞表面蛋白和转录组。这使得可以更准确地定义细胞类型，特别是当转录组信息不足以区分细胞亚型时。当你需要高分辨率的细胞类型鉴定时，CITE-seq提供了额外的蛋白质层面信息。 scATAC-seq与scRNA-seq的整合可以同时揭示基因调控和基因表达信息。通过整合染色质可及性和基因表达数据，可以推断转录因子活性和基因调控网络。当你研究基因调控机制时，多组学整合提供了更全面的视角。 scRNA-seq与空间转录组的整合可以将单细胞分辨率的信息映射到组织空间位置。Seurat的锚点方法可以用于整合scRNA-seq和空间转录组数据，推断空间位置上每个点的细胞类型组成。当你需要理解细胞在组织中的空间分布时，这种整合分析非常有价值。 .. code-block:: r # CITE-seq数据处理 library(Seurat) cbmc <- Read10X_h5("filtered_feature_bc_matrix.h5") cbmc_rna <- cbmc$`Gene Expression` cbmc_adt <- cbmc$`Antibody Capture` cbmc_seurat <- CreateSeuratObject(counts = cbmc_rna) cbmc_seurat[["ADT"]] <- CreateAssayObject(counts = cbmc_adt) cbmc_seurat <- NormalizeData(cbmc_seurat, assay = "ADT", method = "CLR") 5.10 转录本从头组装（de novo 转录组组装） ----------------------------------------- 转录本从头组装是指在缺乏高质量参考基因组的情况下，仅利用RNA-seq读段之间的重叠关系重构转录本序列。这对于非模式生物、新发现的物种或参考基因组不完整的物种尤为重要。 无参考基因组物种的转录组分析需要从头组装作为起点。当你的研究物种没有可用的参考基因组，或者参考基因组质量较差时，从头组装是获得转录本序列的唯一途径。组装得到的转录本可以作为后续定量、差异表达分析和功能注释的基础。 短读长组装器利用de Bruijn图方法从短读长数据构建转录本。Trinity是目前最广泛使用的短读长转录组组装工具，它使用Inchworm、Chrysalis和Butterfly三个模块逐步构建转录本。SOAPdenovo-Trans是SOAPdenovo的转录组版本，针对转录本的特性进行了优化。Trans-ABySS是ABySS的转录组组装模式。当你使用短读长数据进行从头组装时，Trinity是首选工具。 长读长组装器可以利用PacBio Iso-Seq或Nanopore全长RNA-seq数据直接获得完整转录本。Iso-Seq是PacBio提供的全长转录本测序和分析流程，可以获得从转录起始到终止的完整序列。TALON是专门用于长读长转录本分析的工具，可以整合多个样本的长读长数据。IDP（Isoform Detection and Prediction）结合长读长和短读长数据进行转录本鉴定。当你有长读长数据时，可以获得更准确的转录本结构。 组装后冗余去除是提高组装质量的重要步骤。由于组装可能产生大量重叠或相似的转录本，需要去除冗余以获得非冗余的转录本集合。CD-HIT-EST基于序列相似性聚类去除冗余转录本。EvidentialGene使用更复杂的策略评估转录本质量并选择代表性序列。当你完成组装后，冗余去除可以简化后续分析。 组装完整性和冗余评估使用TransRate和BUSCO等工具。TransRate基于读段比对质量评估组装的完整性和准确性。BUSCO（Benchmarking Universal Single-Copy Orthologs）通过检查保守的单拷贝直系同源基因是否存在来评估组装完整性。当你评估组装质量时，这些工具提供了客观的量化指标。 转录本注释包括预测ORF、同源搜索和功能注释。TransDecoder可以预测转录本的开放阅读框，识别可能的编码区域。BLAST或DIAMOND可以将转录本与已知数据库比对进行同源注释。eggNOG-mapper、InterProScan等工具可以进行功能注释。当你获得组装的转录本后，注释是理解其功能的关键步骤。 组装后定量可以使用Salmon或Kallisto直接对组装的转录本进行定量。Salmon的索引模式可以为组装的转录本建立索引，然后进行伪比对定量。这使得在无参考基因组的情况下也能进行表达量分析。当你需要比较不同条件下的基因表达时，定量是必要的步骤。 差异表达在无参考下进行可以使用RSEM结合Trinity的差异分析流程。Trinity提供了完整的差异表达分析流程，包括使用RSEM进行定量、使用edgeR或DESeq2进行差异分析。当你完成从头组装后，可以使用这个流程进行差异表达分析。 .. code-block:: bash # Trinity从头转录组组装 Trinity --seqType fq --max_memory 100G --CPU 16 \ --left left.fq.gz --right right.fq.gz \ --output trinity_out_dir # 组装质量评估 TransRate --assembly trinity_out_dir/Trinity.fasta \ --left left.fq.gz --right right.fq.gz \ --output transrate_output # BUSCO评估 busco -i Trinity.fasta -l metazoa_odb10 -m transcriptome -o busco_output # 冗余去除 cd-hit-est -i Trinity.fasta -o Trinity_nr.fasta -c 0.95 -n 10 # 组装后定量 salmon index -t Trinity_nr.fasta -i trinity_index salmon quant -i trinity_index -l A -1 R1.fq.gz -2 R2.fq.gz \ -p 8 -o salmon_quant 5.11 等位基因特异性表达 ----------------------- 等位基因特异性表达（Allele-Specific Expression，ASE）是指同源染色体上两个等位基因表达水平的不平衡现象。在二倍体生物中，大多数基因的两个等位基因表达水平相近，但某些情况下会出现一个等位基因显著高于另一个的表达偏向。 等位基因识别需要变体信息（SNP）。只有当一个基因座存在杂合SNP时，才能区分两个等位基因并检测ASE。因此，ASE分析需要先进行基因型分析或使用已知的SNP信息。当你进行ASE分析时，需要确保有足够的杂合位点覆盖。 读取比对时保留等位基因来源是ASE分析的关键。常规比对可能将reads随机分配到参考等位基因或替代等位基因，引入偏倚。WASP工具通过在比对前将reads中的SNP位点替换为替代等位基因，比较替换前后比对结果来识别和去除比对偏倚。GATK的ASEReadCounter可以统计每个等位基因的reads数量。当你进行ASE分析时，需要使用专门的工具处理比对偏倚问题。 等位基因不平衡检验使用Fisher精确检验或二项式检验。对于每个杂合位点，统计参考等位基因和替代等位基因的reads数量，然后检验是否偏离1:1的期望比例。当你进行ASE检验时，需要考虑多重检验校正和效应大小阈值。 亲本来源特异表达在杂交物种RNA-seq中可以检测。当使用来自不同物种或品系的亲本杂交产生的F1代时，可以根据SNP区分来自父本和母本的等位基因，检测印记基因或亲本偏向表达。当你研究基因组印记或亲本效应时，杂交系统提供了清晰的遗传标记。 等位基因表达的环境依赖性是ASE研究的重要方向。某些等位基因的偏向表达可能在特定条件下（如应激、发育阶段、组织类型）发生变化，揭示顺式调控变异的功能影响。当你研究基因型-表型关系时，ASE分析可以揭示遗传变异如何影响基因表达。 .. code-block:: bash # GATK ASE分析 gatk ASEReadCounter -R reference.fa -I sample.bam \ -V variants.vcf.gz -O ase_counts.table # WASP去除比对偏倚 python run_wasp.py --is_paired_end --output_dir wasp_output \ samples.bam snps_file.vcf .. code-block:: r # ASE统计检验 ase_data <- read.table("ase_counts.table", header = TRUE) ase_data$p_value <- apply(ase_data, 1, function(row) { ref_count <- as.numeric(row["refCount"]) alt_count <- as.numeric(row["altCount"]) if (ref_count + alt_count < 10) return(NA) binom.test(c(ref_count, alt_count), p = 0.5)$p.value }) ase_data$FDR <- p.adjust(ase_data$p_value, method = "BH") sig_ase <- subset(ase_data, FDR < 0.05 & abs(refCount - altCount) > 5) cat("显著ASE位点数:", nrow(sig_ase), "\n") 5.12 RNA编辑分析 ---------------- RNA编辑是指转录后RNA序列的改变，使得RNA序列与对应的DNA模板序列不同。最常见的RNA编辑类型是腺苷到肌苷（A-to-I）编辑，由ADAR酶催化，肌苷在翻译和剪接中被识别为鸟苷。 RNA编辑类型以A-to-I最为常见，在人类转录组中广泛存在，特别是在Alu重复序列区域。C-to-U编辑由APOBEC酶催化，在载脂蛋白B等特定基因中发生。其他类型的RNA编辑相对罕见。当你研究RNA编辑时，A-to-I编辑是主要关注对象。 检测工具包括REDItools、JACUSA和RNAEditor。REDItools是一套用于检测和定量RNA编辑的工具，它比较RNA-seq和DNA-seq数据识别编辑位点。JACUSA使用更复杂的统计模型检测RNA编辑，支持多种实验设计。RNAEditor提供了完整的分析流程，包括编辑检测、注释和可视化。当你进行RNA编辑分析时，需要同时有RNA-seq和匹配的DNA-seq数据。 区分编辑位点与基因组变异是RNA编辑分析的核心挑战。单核苷酸多态性（SNP）也会表现为RNA与DNA序列的差异，需要通过比较多个样本或使用数据库过滤已知的SNP位点。当你检测到潜在的编辑位点时，需要排除SNP的可能性。 编辑频率计算表示编辑发生的比例。对于每个编辑位点，编辑频率等于编辑reads数除以总reads数。编辑频率可以从0（无编辑）到1（完全编辑），反映了编辑的效率和调控。当你报告RNA编辑结果时，编辑频率是重要的定量指标。 功能影响包括改变氨基酸序列、影响剪接和改变microRNA结合。编码区的A-to-I编辑可能改变密码子，导致氨基酸替换（如谷氨酰胺变为精氨酸）。剪接位点附近的编辑可能影响剪接位点识别。3'UTR区域的编辑可能改变miRNA结合位点，影响基因表达调控。当你发现RNA编辑位点后，评估其功能影响是重要的后续步骤。 .. code-block:: bash # REDItools RNA编辑检测 python REDItoolDnaVcf.py -i rna.bam -j dna.vcf -f reference.fa \ -o reditools_output # JACUSA RNA编辑检测 java -jar jacusa.jar call-2 -r reference.fa -p 8 \ sample1.rna.bam sample2.rna.bam -o jacusa_output.txt # RNAEditor完整流程 RNAEditor.py -r reference.fa -g annotation.gtf \ -b rna.bam -d dna.bam -o rnaeditor_output 5.13 融合基因检测 ----------------- 融合基因是由两个不同基因的片段连接形成的嵌合基因，通常由染色体易位、缺失或倒位等结构变异产生。融合基因在肿瘤中常见，某些融合基因（如BCR-ABL）是重要的诊断标志物和治疗靶点。 肿瘤融合基因检测工具包括STAR-Fusion、FusionCatcher、Arriba和FusionMap。STAR-Fusion基于STAR比对器检测的嵌合读段识别融合基因。FusionCatcher专门设计用于肿瘤样本，整合了多种检测策略。Arriba是快速的融合基因检测工具，针对RNA-seq数据优化。FusionMap可以同时处理RNA-seq和DNA-seq数据。当你研究肿瘤转录组时，融合基因检测是重要的分析内容。 嵌合比对识别是融合基因检测的基础。跨越融合断点的reads会比对到两个不同的基因位置，形成嵌合比对。STAR等比对器可以识别并报告这些嵌合比对。当你进行融合基因检测时，高质量的比对是前提条件。 融合断点边界验证需要检查断点附近的reads覆盖和序列质量。真实的融合基因应该在断点两侧都有足够的reads支持，且断点位置符合剪接位点规律（通常位于外显子边界）。当你发现候选融合基因时，需要仔细检查断点质量。 融合基因翻译产物预测可以推断融合蛋白的结构和功能。融合可能保留、破坏或创造新的功能域，影响蛋白质的活性或定位。当你发现融合基因后，预测其蛋白质产物有助于理解其致癌机制。 公共数据库如ChiTars和FusionDB收录了已知的融合基因及其相关信息。ChiTars专注于肿瘤相关的融合基因，FusionDB提供更全面的融合基因注释。当你发现新的融合基因时，可以查询这些数据库了解是否已有报道。 .. code-block:: bash # STAR-Fusion融合基因检测 STAR-Fusion --genome_lib_dir /path/to/CTAT_resource_lib \ --left_fq R1.fq.gz --right_fq R2.fq.gz \ --output_dir star_fusion_output --CPU 8 # Arriba快速融合检测 arriba -x aligned.bam -o fusions.tsv -O fusions.discarded.tsv \ -a reference.fa -g annotation.gtf # FusionCatcher肿瘤融合检测 fusioncatcher.py -i R1.fq.gz,R2.fq.gz -o fusioncatcher_output \ -d /path/to/bowtie_index 5.14 转录组进化分析 ------------------- 转录组进化分析研究基因表达在进化过程中的变化规律，揭示基因表达调控的演化机制和适应性意义。 同源转录本识别是跨物种转录组比较的基础。直系同源基因（orthologs）是不同物种中由共同祖先基因进化而来的基因，通常保留相似的功能。旁系同源基因（paralogs）是基因复制产生的同源基因，可能分化出不同的功能。当你进行跨物种转录组分析时，准确识别直系同源基因对是关键步骤。 跨物种转录组比较可以研究表达量进化和Phylostratigraphy。表达量进化研究基因表达水平在进化过程中的变化速率和模式。Phylostratigraphy通过将基因按其进化起源分组，研究基因组的进化历史。当你研究物种分化或适应性进化时，转录组比较提供了功能层面的视角。 基因表达分化可以使用多种度量方法。对于编码序列，dN/dS比值可以衡量选择压力；对于表达量，可以使用表达差异距离、表达分歧指数等方法量化表达分化程度。当你比较不同物种或群体的转录组时，表达分化分析可以识别快速进化的基因。 系统发育转录组学将转录组数据置于系统发育框架中分析。通过构建基因表达的系统发育树，可以研究表达模式的进化历史。当你研究多个相关物种时，系统发育方法可以区分系统发育信号和适应性变化。 新基因起源鉴定关注孤儿基因（orphan genes）。孤儿基因是指在特定类群中缺乏明显同源物的基因，可能代表新起源的基因。研究孤儿基因的表达模式和功能有助于理解新基因的起源和进化。当你研究基因组创新时，孤儿基因是重要的研究对象。 .. code-block:: r # 跨物种表达分化分析 library(orthologr) orthologs <- orthologs(query_file = "species1_proteins.fa", subject_file = "species2_proteins.fa") # 表达分化计算 expr_divergence <- function(expr1, expr2, ortholog_pairs) { divergence <- sapply(1:nrow(ortholog_pairs), function(i) { gene1 <- ortholog_pairs$gene1[i] gene2 <- ortholog_pairs$gene2[i] if (gene1 %in% rownames(expr1) && gene2 %in% rownames(expr2)) { mean1 <- mean(expr1[gene1, ]) mean2 <- mean(expr2[gene2, ]) return(abs(log2(mean1 + 1) - log2(mean2 + 1))) } return(NA) }) return(divergence) } 5.15 空间转录组学 ----------------- 空间转录组学技术在保留组织空间位置信息的同时测量基因表达，使得可以在组织原位研究细胞类型分布和细胞间相互作用。 空间转录组技术包括10X Visium、Slide-seq、MERFISH和seqFISH等。10X Visium使用带有空间条形码的捕获阵列，每个spot捕获该位置的转录本，分辨率约为55微米。Slide-seq使用更密集的珠阵列，分辨率可达10微米。MERFISH和seqFISH是基于成像的技术，使用荧光标记的探针逐个检测转录本，可以达到单分子分辨率。当你选择空间转录组技术时，需要权衡分辨率、通量和成本。 数据格式包括SPOT或细胞坐标加表达矩阵。对于基于阵列的技术，数据以spot为单位，每个spot有空间坐标和基因表达量。对于基于成像的技术，数据以单细胞为单位，每个细胞有空间坐标和基因表达矩阵。当你处理空间转录组数据时，需要了解数据的空间坐标系统。 空间表达模式可视化使用Squidpy、Giotto和Seurat等工具。这些工具可以在组织切片图像上叠加基因表达信息，生成空间特征图。当你展示空间转录组结果时，可视化是直观展示空间模式的重要手段。 空间变异基因检测使用SpatialDE和Trendsceek等工具。这些方法识别在空间上具有非随机表达模式的基因，可能标记特定的组织区域或功能状态。当你研究组织结构时，空间变异基因可以揭示区域特异性功能。 空间区域分割与细胞类型mapping将空间位置划分为功能区域并注释细胞类型。可以使用聚类方法识别空间上连续的表达区域，然后基于标记基因注释细胞类型。对于低分辨率数据，可以使用解卷积方法推断每个spot的细胞类型组成。当你分析空间转录组数据时，区域分割和细胞类型注释是核心分析步骤。 空间细胞间相互作用扩展了CellChat等工具到空间数据。通过考虑细胞的空间临近性，可以更准确地推断细胞间通讯。当你研究肿瘤微环境或发育组织时，空间细胞间相互作用分析可以揭示重要的调控网络。 .. code-block:: r # Seurat空间转录组分析 library(Seurat) st_data <- Load10X_Spatial(data.dir = "spatial_output/") st_data <- SCTransform(st_data, assay = "Spatial") st_data <- RunPCA(st_data) st_data <- FindNeighbors(st_data, dims = 1:30) st_data <- FindClusters(st_data, resolution = 0.8) SpatialDimPlot(st_data, label = TRUE) # 空间变异基因检测 library(SpatialDE) spatial_results <- SpatialDE::run(counts, coordinates) sig_genes <- subset(spatial_results, qvalue < 0.05) cat("空间变异基因数:", nrow(sig_genes), "\n") 5.16 时序RNA-seq分析 --------------------- 时序RNA-seq分析研究基因表达随时间的变化规律，揭示动态生物学过程（如发育、应激响应、药物处理）中的转录调控机制。 时间序列设计需要合理的时间间隔和采样点数量。时间间隔应该足够密集以捕捉快速的变化，同时覆盖足够长的时间跨度以观察完整的过程。采样点数量通常需要5个以上才能进行有意义的时间趋势分析。当你设计时序RNA-seq实验时，需要根据生物学过程的时间尺度确定采样策略。 差异表达分析适用于时间序列的方法包括limma时间趋势、DESeq2 LRT和maSigPro。limma可以使用自然样条或多项式建模时间效应。DESeq2的似然比检验可以比较包含时间效应的模型与不包含时间效应的模型。maSigPro专门为时间序列和多重条件设计，可以识别时间依赖的差异表达。当你分析时序数据时，需要选择适合实验设计的统计方法。 聚类表达谱使用Mfuzz和STEM等工具。Mfuzz使用模糊c均值聚类将基因按表达模式分组。STEM（Short Time-series Expression Miner）专门针对短时间序列设计，提供预定义的表达模式模板。当你发现大量时间依赖的差异基因后，聚类可以帮助识别具有相似变化模式的基因群。 发现阶段性表达基因是时序分析的重要目标。某些基因在特定时间窗口高表达，标记特定的发育阶段或响应阶段。识别这些阶段性表达基因可以帮助理解过程的时序控制。当你研究发育或分化过程时，阶段性表达基因是重要的标记。 转录调控网络动态利用时间延迟推断调控关系。转录因子的表达变化通常先于其靶基因的变化，通过分析表达变化的时间延迟可以推断调控方向和强度。当你研究转录调控时，时间序列数据提供了推断因果关系的独特机会。 .. code-block:: r # maSigPro时间序列差异表达 library(maSigPro) design <- make.design.matrix(sample_info, degree = 3) fit <- p.vector(counts, design, Q = 0.05) tstep <- T.fit(fit, step.method = "backward") sig_genes <- get.siggenes(tstep, rsq = 0.6, vars = "groups") cat("时间依赖差异基因数:", length(sig_genes$genes), "\n") # Mfuzz聚类 library(Mfuzz) expr_set <- new("ExpressionSet", exprs = as.matrix(normalized_counts)) expr_set <- standardise(expr_set) c <- mfuzz(expr_set, c = 6, m = 1.25) mfuzz.plot(expr_set, cl = c) 5.17 RNA-seq分析工作流与最佳实践 --------------------------------- RNA-seq分析工作流将各个分析步骤整合为自动化流程，确保分析的可重复性和标准化。工作流管理系统使得复杂的分析流程可以在不同计算环境中一致运行。 nf-core/rnaseq是基于Nextflow的标准化RNA-seq分析流程，由nf-core社区维护。它整合了从质控、比对、定量到差异表达的完整分析链，支持多种比对器和定量工具的选择。流程自动生成MultiQC报告和R Markdown分析报告。当你需要一个开箱即用的RNA-seq分析流程时，nf-core/rnaseq是推荐的选择。 Snakemake RNA-seq流程使用Python风格的规则定义，适合Python用户构建自定义流程。Snakemake的优势在于其直观的规则定义方式和强大的并行执行能力。当你需要定制化分析流程时，Snakemake提供了灵活的框架。 使用Docker/Singularity容器保证可重复性。容器将软件及其所有依赖打包，确保分析在不同计算环境中产生一致的结果。nf-core/rnaseq和大多数现代工作流都支持容器化执行。当你分享分析流程或在不同平台运行时，容器化是确保可重复性的关键。 报告生成使用MultiQC汇总质量控制结果，使用R Markdown生成分析报告。MultiQC整合多个样本的QC报告为单一文档，便于批量评估数据质量。R Markdown可以将代码、结果和解释整合为可重复的报告文档。当你完成RNA-seq分析后，完整的报告文档对于结果分享和审稿回复都很重要。 快速原型示例数据集可以从GEO和ArrayExpress获取。这些公共数据库收录了大量已发表的RNA-seq数据集，可以用于测试分析流程或方法开发。当你开发新方法或学习RNA-seq分析时，公共数据集是宝贵的资源。 跨实验比较与meta分析整合多个研究的差异表达结果。由于不同研究可能使用不同的实验设计、分析流程和报告标准，meta分析需要仔细处理异质性和批次效应。当你需要综合多个研究的证据时，meta分析可以提供更稳健的结论。 De Bruijn图组装器用于基因组组装，在转录组从头组装中也有应用。SPAdes、ABySS、SOAPdenovo2和Velvet是常用的de Bruijn图组装器。SPAdes针对细菌和低复杂度基因组优化，ABySS适合大规模基因组组装，SOAPdenovo2是SOAPdenovo的升级版本，Velvet是最早的de Bruijn图组装器之一。 最优k-mer选择使用KmerGenie等工具。k-mer大小影响组装的连续性和准确性，较小的k-mer可以跨越重复区域但可能引入错误连接，较大的k-mer更准确但可能断裂。KmerGenie通过分析k-mer频率分布推荐最优的k值。当你进行de novo组装时，选择合适的k-mer大小是重要的参数调优步骤。 读长纠错使用BayesHammer和Musket等工具。BayesHammer是SPAdes内置的纠错模块，使用贝叶斯模型纠正测序错误。Musket使用k-mer频率信息进行纠错。纠错可以提高组装质量，特别是对于高错误率的数据。当你使用错误率较高的测序数据时，纠错是提高组装质量的有效手段。 脱泡去除与图简化是de Bruijn图组装的后续处理步骤。气泡（bubble）是由于测序错误或变异产生的图结构异常，需要去除以获得线性contig。图简化通过移除低覆盖度的路径和合并线性区域来简化图结构。当你理解de Bruijn图组装原理时，这些步骤是理解组装过程的关键。 Scaffolding方法使用SSPACE、Opera和SOAPdenovo Scaffolding等工具。Scaffolding利用配对端信息或长距离连接信息将contig连接成更长的scaffold，中间的gap用N表示。当你需要更长的组装序列时，scaffolding是必要的步骤。 Gap闭合使用GapFiller、GapCloser和Sealer等工具。这些工具利用配对端reads或长读长数据填充scaffold中的gap区域，提高组装完整性。当你需要完整的基因组序列时，gap闭合是最后的组装优化步骤。 .. code-block:: bash # nf-core/rnaseq工作流 nextflow run nf-core/rnaseq -profile docker \ --input samplesheet.csv --outdir results \ --genome GRCh38 --aligner star_salmon # Snakemake RNA-seq流程示例 # Snakefile内容见下文 .. code-block:: python # Snakemake RNA-seq流程 (Snakefile) rule all: input: "results/differential_expression/DEG_results.csv" rule fastqc: input: "data/{sample}.fastq.gz" output: "qc/{sample}_fastqc.html" shell: "fastqc {input} -o qc/" rule align: input: r1 = "data/{sample}_R1.fastq.gz", r2 = "data/{sample}_R2.fastq.gz" output: "results/aligned/{sample}.bam" shell: "STAR --runThreadN 8 --genomeDir index " "--readFilesIn {input.r1} {input.r2} --outFileNamePrefix results/{wildcards.sample} " "--outSAMtype BAM SortedByCoordinate" rule deseq2: input: expand("results/counts/{sample}.txt", sample=SAMPLES) output: "results/differential_expression/DEG_results.csv" script: "scripts/run_deseq2.R" .. code-block:: bash # de Bruijn图组装器 # SPAdes组装 spades.py -1 R1.fq.gz -2 R2.fq.gz -o spades_output -t 16 -m 64 # ABySS组装 abyss-pe k=64 name=abyss in='R1.fq.gz R2.fq.gz' np=16 # KmerGenie最优k-mer选择 kmergenie R1.fq.gz,R2.fq.gz -o kmergenie_output # Gap闭合 gapfiller.pl -l library.txt -s contigs.fa -b gapfiller_output