表观遗传学分析 ================== 表观遗传学是研究在不改变DNA序列的情况下，基因表达发生可遗传变化的一门学科。与基因组学研究DNA序列本身不同，表观遗传学关注的是DNA甲基化、组蛋白修饰、染色质可及性和三维结构等调控层次，它们共同决定了哪些基因在特定细胞类型和特定发育阶段被激活或沉默。表观遗传修饰具有动态性和可逆性，受环境因素影响显著，这使得它们在发育调控、疾病发生和衰老过程中扮演着关键角色。当你在研究中发现基因组层面的变异无法完全解释表型差异时，表观遗传学分析往往能够提供重要的补充信息。本章将系统介绍DNA甲基化、染色质可及性、组蛋白修饰、染色质互作、RNA修饰等核心表观遗传学分析内容，涵盖从实验技术到数据分析的完整流程。 6.1 DNA甲基化分析 ----------------- DNA甲基化是最早被发现也是研究最为深入的表观遗传修饰之一。在哺乳动物中，DNA甲基化主要发生在CpG二核苷酸的胞嘧啶5号碳原子上，形成5-甲基胞嘧啶（5mC）。CpG位点在基因组中的分布并不均匀，大约70-80%的CpG分散在基因组中且通常处于甲基化状态，而剩余的CpG聚集形成所谓的CpG岛，主要位于基因启动子区域，通常保持非甲基化状态。当启动子区域的CpG岛发生异常甲基化时，会导致基因转录沉默，这一机制在肿瘤抑制基因的失活中尤为常见。DNA甲基化分析的目标就是定量测量基因组中各个CpG位点的甲基化水平，并比较不同样本之间的差异，从而揭示甲基化在基因调控和疾病中的作用。 6.1.1 甲基化检测技术 ^^^^^^^^^^^^^^^^^^^^^ 重亚硫酸盐测序（BS-seq）是DNA甲基化检测的金标准方法。其核心原理是利用重亚硫酸盐（bisulfite）处理DNA，将未甲基化的胞嘧啶（C）转化为尿嘧啶（U），而甲基化的胞嘧啶（5mC）则保持不变。经过PCR扩增后，U被读作T，而5mC仍然读作C。通过比较处理前后的序列变化，可以在单碱基分辨率下确定每个CpG位点的甲基化状态。当你需要获得全基因组范围内精确的甲基化图谱时，BS-seq是最可靠的选择，但其成本较高，因为重亚硫酸盐处理会导致DNA降解，需要较高的测序深度来补偿。 简化代表性重亚硫酸盐测序（RRBS）通过限制性内切酶（如MspI）切割富含CpG的片段，然后对这些片段进行重亚硫酸盐测序，从而以较低的成本覆盖基因组中大部分的CpG岛和启动子区域。RRBS特别适合于大样本量的甲基化研究，当你关注的主要是启动子区域的甲基化变化而非全基因组范围时，RRBS在成本和覆盖度之间提供了良好的平衡。 甲基化DNA免疫沉淀测序（MeDIP-seq）利用抗5mC抗体特异性地富集甲基化的DNA片段，然后对富集的片段进行测序。与BS-seq不同，MeDIP-seq不能提供单碱基分辨率的甲基化信息，而是反映基因组区域的整体甲基化水平。MeDIP-seq的优势在于成本较低且不需要重亚硫酸盐处理，适合于大尺度的甲基化模式筛查。当你需要进行初步的甲基化组学探索或处理大量样本时，MeDIP-seq是一个实用的选择。 甲基化敏感限制性内切酶测序（MRE-seq）利用对甲基化敏感的限制性内切酶（如HpaII）只能切割非甲基化CpG位点的特性，通过酶切和测序来识别未甲基化的区域。MRE-seq通常与MeDIP-seq联合使用，前者检测非甲基化区域，后者检测甲基化区域，两者互补可以提供更全面的甲基化图谱。 全基因组重亚硫酸盐测序（WGBS）是BS-seq在全基因组层面的实现，它对整个基因组进行重亚硫酸盐处理和深度测序，理论上可以覆盖基因组中所有的CpG位点。WGBS是甲基化检测最全面的方法，但成本也最高，通常每个样本需要30倍以上的基因组覆盖度才能获得可靠的甲基化调用。当你需要绘制完整的甲基化图谱或发现新的甲基化区域时，WGBS是首选方法。 靶向甲基化测序结合了目标区域捕获技术和重亚硫酸盐测序，通过探针捕获感兴趣的基因组区域（如特定基因的启动子、增强子或已知的差异甲基化区域），然后对捕获的区域进行BS-seq。这种方法在保持单碱基分辨率的同时大幅降低了测序成本，特别适合于验证WGBS发现的候选区域或在临床研究中聚焦于特定的甲基化标志物。 氧化重亚硫酸盐测序（oxBS-seq）是区分5mC和5hmC（5-羟甲基胞嘧啶）的关键技术。5hmC是TET蛋白氧化5mC的产物，在脑组织中含量丰富，具有重要的调控功能。oxBS-seq先用氧化剂将5hmC氧化为5fC（5-甲酰胞嘧啶），再经重亚硫酸盐处理使5fC转化为U，从而单独检测5mC的水平。将常规BS-seq的结果减去oxBS-seq的结果，即可得到5hmC的分布。当你研究脑组织或胚胎干细胞等5hmC含量较高的样本时，oxBS-seq能够揭示常规BS-seq无法区分的羟甲基化信息。 Tet辅助重亚硫酸盐测序（TAB-seq）是另一种区分5mC和5hmC的技术，与oxBS-seq相反，TAB-seq通过TET蛋白将5mC氧化为5caC，同时保护5hmC不被转化，从而单独检测5hmC的水平。TAB-seq在单碱基分辨率下检测5hmC，适合于需要精确定位5hmC位点的研究。 纳米孔直接RNA/DNA甲基化检测（Nanopore）利用牛津纳米孔测序技术的电信号变化来直接检测DNA或RNA上的甲基化修饰，无需重亚硫酸盐处理。纳米孔测序通过测量DNA链穿过纳米孔时引起的电流变化来读取序列信息，甲基化的碱基会产生与未甲基化碱基不同的电信号特征。这种方法的优势在于可以直接检测多种修饰类型（5mC、5hmC、6mA等），且能够获得长读长信息，有助于在重复区域和结构变异区域进行甲基化分析。当你需要研究等位基因特异性甲基化或复杂基因组区域的甲基化模式时，纳米孔技术提供了独特的优势。 PacBio HiFi甲基化检测利用PacBio单分子实时测序（SMRT）的动力学信息来检测DNA甲基化。在SMRT测序中，聚合酶合成DNA链时，遇到甲基化碱基会导致聚合速率发生变化，这种动力学特征可以被用来推断甲基化的存在。HiFi读长兼具高准确性和长读长的优势，使得PacBio在复杂区域的甲基化检测中表现出色。 单细胞重亚硫酸盐测序（scBS-seq，scWGBS）将BS-seq技术应用于单个细胞，能够揭示细胞群体中甲基化的异质性。由于单个细胞的DNA量极少，scBS-seq需要经过全基因组扩增，这会引入覆盖度不均匀和扩增偏倚等问题。尽管如此，scBS-seq在研究早期胚胎发育、肿瘤异质性和罕见细胞类型的甲基化特征方面具有不可替代的价值。当你需要了解细胞间甲基化差异而非群体平均水平时，单细胞甲基化测序是必要的工具。 甲基化阵列（Illumina 450K和EPIC 850K）是基于探针杂交的甲基化检测平台，通过设计针对特定CpG位点的探针来定量甲基化水平。450K阵列覆盖约45万个CpG位点，而EPIC 850K阵列覆盖约85万个CpG位点，主要增加了增强子区域的覆盖。甲基化阵列的优势在于成本低、通量高、数据质量稳定，是大规模EWAS研究中最常用的平台。当你需要处理数百甚至数千个样本的甲基化分析时，阵列平台在成本效益和数据一致性方面具有明显优势。 .. code-block:: bash # 不同甲基化检测技术的典型测序深度要求 # WGBS: ~30x 基因组覆盖度 # RRBS: ~10-30x 有效覆盖度（仅覆盖CpG富集区域） # MeDIP-seq: ~20-40M reads # EPIC 850K阵列: 无需测序，芯片扫描 # RRBS文库制备示例（使用MspI酶切） # MspI识别CCGG位点，切割后产生含CpG的片段 # 典型片段大小: 40-220bp 6.1.2 甲基化数据分析流程 ^^^^^^^^^^^^^^^^^^^^^^^^^ 重亚硫酸盐比对是甲基化数据分析的第一步，其核心挑战在于重亚硫酸盐处理后基因组序列发生了C到T的转换，使得常规的比对工具无法直接使用。Bismark是最流行的重亚硫酸盐比对工具，它通过将参考基因组和测序读段分别进行C-to-T和G-to-A的转换，然后使用Bowtie2或HISAT2进行三方向比对（原始读段、C-to-T转换读段、G-to-A转换读段），最终确定读段的最佳比对位置和原始链方向。BSMAP是另一个常用的比对工具，它使用容错匹配策略来处理重亚硫酸盐转换。BS-Seeker2同样基于三方向转换策略，但支持更多的比对器后端。BWA-meth是BWA-MEM的重亚硫酸盐适配版本，速度较快且内存占用较低。GemBS是一个集成了比对、甲基化调用和下游分析的综合流程。当你处理大规模的WGBS数据时，Bismark和BWA-meth是效率较高的选择，而GemBS适合需要标准化流程的项目。 .. code-block:: bash # Bismark比对流程 # 1. 构建重亚硫酸盐基因组索引 bismark_genome_preparation /path/to/reference_genome/ # 2. 比对（双端测序） bismark --genome /path/to/reference_genome/ \ -1 sample_R1.fq.gz -2 sample_R2.fq.gz \ --parallel 8 --output_dir bismark_output/ # 3. 去重 deduplicate_bismark --paired bismark_output/sample_R1_bismark_bt2_pe.bam # 4. 甲基化提取 bismark_methylation_extractor \ --paired-end --gzip --comprehensive \ --bedGraph --CX_context \ bismark_output/sample_R1_bismark_bt2_pe.deduplicated.bam 甲基化调用是从比对结果中提取每个CpG位点甲基化信息的过程。Bismark methylation extractor是Bismark流程自带的甲基化调用工具，它从BAM文件中统计每个CpG位点上支持甲基化（C）和非甲基化（T）的读段数量。MethylDackel是另一个快速的甲基化调用工具，它可以处理Bismark或其他比对器产生的BAM文件，提取每个位点的甲基化和非甲基化计数。甲基化调用的结果通常以每个CpG位点的甲基化读段数和总读段数的形式呈现。 甲基化水平的计算通常使用beta值，其公式为 $\beta = mC / (mC + umC)$ ，其中mC是支持甲基化的读段数，umC是支持非甲基化的读段数。beta值的取值范围是0到1，0表示完全非甲基化，1表示完全甲基化。另一种常用的度量是M值，其公式为 $M = log2(\beta / (1 - \beta))$ ，M值在统计检验中更为稳定，适合用于差异甲基化分析。当你进行下游的统计检验时，M值通常比beta值更适合，因为M值近似服从正态分布；而在结果展示和生物学解释时，beta值更为直观。 .. code-block:: r # 甲基化水平计算 mC <- 8 umC <- 12 beta <- mC / (mC + umC) cat(sprintf("Beta值: %.4f\n", beta)) M <- log2(beta / (1 - beta)) cat(sprintf("M值: %.4f\n", M)) # 批量计算beta值和M值 meth_counts <- data.frame( mC = c(8, 15, 2, 20, 5), umC = c(12, 5, 18, 0, 15) ) meth_counts$beta <- with(meth_counts, mC / (mC + umC)) meth_counts$M <- with(meth_counts, log2(beta / (1 - beta))) print(meth_counts) 过滤低覆盖度CpG位点是为了确保甲基化水平估计的可靠性。覆盖度过低的位点其甲基化水平的估计方差较大，容易引入假阳性或假阴性的差异甲基化结果。通常建议每个CpG位点至少有5-10倍的覆盖度才纳入分析，对于WGBS数据可能需要更高的阈值。此外，还需要过滤掉位于已知SNP位点的CpG，因为SNP导致的C-to-T变化会被误认为甲基化状态的变化。当你处理阵列数据时，还需要过滤检测p值大于0.01的探针和交叉反应性探针。 CpG岛注释是将甲基化位点映射到基因组功能区域的过程。CpG岛通常定义为长度至少200bp、GC含量大于50%、观察值与期望值之比大于0.6的基因组区域。除了CpG岛本身，还需要注释CpG岛的上下游岸区（shore，距离CpG岛2kb以内）和架子区（shelf，距离CpG岛2-4kb），因为这些区域的甲基化变化同样具有重要的生物学意义。研究表明，组织特异性和疾病相关的差异甲基化更多地发生在CpG岛的岸区而非岛内。 甲基化数据的存储格式有多种。BEDGraph格式记录每个基因组区域的甲基化水平，适合在基因组浏览器中可视化。bigWig是BEDGraph的二进制索引版本，加载速度更快，适合大规模数据的可视化。methylation BED格式记录每个CpG位点的甲基化计数和百分比。CX report是Bismark输出的详细甲基化报告，包含所有胞嘧啶上下文（CpG、CHG、CHH）的甲基化信息。当你需要在不同工具之间交换甲基化数据时，了解这些格式的特点和转换方法非常重要。 甲基化可视化是理解甲基化模式的重要手段。IGV（Integrative Genomics Viewer）可以直接加载BAM文件和BEDGraph文件，逐个位点查看甲基化状态。UCSC Genome Browser提供了自定义track功能，可以将甲基化数据作为自定义轨道添加到基因组浏览器中。MethBank是一个专门用于DNA甲基化数据的数据库和可视化平台，提供了多种物种的参考甲基化组数据。当你需要在候选基因区域检查甲基化模式时，IGV是最方便的本地可视化工具。 差异甲基化分析是比较不同组别之间甲基化水平的统计方法。DSS（Dispersion Shrinkage for Sequencing data）使用贝叶斯方法估计甲基化水平的离散度，然后进行Wald检验或似然比检验来识别差异甲基化位点，特别适合基于测序的甲基化数据。limma包虽然最初设计用于基因表达分析，但也可以应用于甲基化阵列数据的差异分析，通过经验贝叶斯方法提高小样本情况下的统计功效。methylKit提供了从比对到差异分析的完整流程，支持多种归一化和统计检验方法。RnBeads是一个综合的甲基化分析框架，支持阵列和测序数据，集成了质量控制、预处理、差异分析和注释等功能。当你处理WGBS或RRBS数据时，DSS是推荐的选择；对于阵列数据，limma和RnBeads更为常用。 .. code-block:: r # 使用DSS进行差异甲基化分析 library(DSS) library(bsseq) # 构建BSseq对象 positions <- c(100, 200, 300, 400, 500) chr <- rep("chr1", 5) # 对照组样本 meth_ctrl1 <- c(8, 15, 2, 20, 5) total_ctrl1 <- c(20, 20, 20, 20, 20) meth_ctrl2 <- c(7, 14, 3, 18, 6) total_ctrl2 <- c(20, 20, 20, 20, 20) # 处理组样本 meth_treat1 <- c(15, 8, 12, 19, 16) total_treat1 <- c(20, 20, 20, 20, 20) meth_treat2 <- c(16, 7, 13, 20, 15) total_treat2 <- c(20, 20, 20, 20, 20) BSobj <- BSseq( chr = chr, pos = positions, M = cbind(meth_ctrl1, meth_ctrl2, meth_treat1, meth_treat2), Cov = cbind(total_ctrl1, total_ctrl2, total_treat1, total_treat2), sampleNames = c("ctrl1", "ctrl2", "treat1", "treat2") ) # 差异甲基化分析 dmlTest <- DMLtest(BSobj, group1 = c("ctrl1", "ctrl2"), group2 = c("treat1", "treat2")) head(dmlTest) # 调用差异甲基化位点（DML） dmls <- callDML(dmlTest, p.threshold = 0.001) head(dmls) # 调用差异甲基化区域（DMR） dmrs <- callDMR(dmlTest, p.threshold = 0.01) head(dmrs) DMR（差异甲基化区域）检测与合并是将相邻的差异甲基化位点合并为连续区域的过程。单个CpG位点的差异甲基化可能由技术噪音引起，而连续多个CpG位点的一致性差异更可能具有生物学意义。DMR检测算法通常基于滑动窗口或区域分割策略，将相邻的差异甲基化位点合并，并评估整个区域的统计显著性。合并后的DMR需要满足一定的条件，如最小长度（通常至少包含3个CpG位点）、最小甲基化差异（如delta beta > 0.2）和统计显著性（如FDR < 0.05）。当你报告差异甲基化结果时，DMR通常比单个DMP更具生物学可解释性。 差异甲基化位点（DMP）统计关注的是单个CpG位点水平的差异分析。DMP分析输出每个CpG位点的统计检验结果，包括甲基化差异（delta beta或diffMeth）、p值和校正后的FDR值。DMP的筛选标准通常包括FDR < 0.05和绝对甲基化差异 > 0.2（对于beta值）。当你需要精确定位甲基化变化的位置时，DMP分析提供了最细粒度的信息，但需要注意多重检验校正的问题，因为人类基因组中约有2800万个CpG位点。 基因组区域甲基化特征分析关注不同功能区域的甲基化模式。启动子区域的甲基化通常与基因沉默相关，当启动子CpG岛发生高甲基化时，转录因子无法结合，基因表达被抑制。基因体（gene body）的甲基化则呈现不同的模式，高甲基化的基因体往往与活跃转录的基因相关，这种正相关关系的机制尚不完全清楚，可能与抑制基因内转录起始和转座子活性有关。增强子区域的甲基化变化对基因表达的影响更为复杂，低甲基化的增强子通常处于活跃状态。重复元件（如LINE、SINE、LTR）的甲基化维持了基因组的稳定性，其去甲基化可能导致转座子激活和基因组不稳定。 6.1.3 甲基化功能解析 ^^^^^^^^^^^^^^^^^^^^^ 启动子甲基化与基因表达负相关分析是甲基化功能研究中最基本的分析之一。在大多数情况下，启动子区域CpG岛的高甲基化与基因表达下调相关，这是因为甲基化的CpG会阻碍转录因子的结合，同时招募甲基化结合蛋白（如MeCP2），进而招募组蛋白去乙酰化酶等抑制性复合物，形成紧密的染色质结构。当你同时拥有甲基化数据和RNA-seq数据时，可以通过计算启动子甲基化水平与基因表达量之间的相关系数来验证这一负相关关系。需要注意的是，这种负相关并非绝对的，基因体甲基化与基因表达可能呈正相关，而增强子甲基化的关系则更为复杂。 甲基化定量性状位点（mQTL）分析是将DNA甲基化水平作为数量性状，寻找与甲基化变异相关的遗传变异位点。mQTL分析的基本思路与eQTL分析类似，通过线性回归模型检验SNP基因型与CpG位点甲基化水平之间的关联。mQTL可以分为cis-mQTL（SNP与CpG位点在基因组上距离较近，通常在1Mb以内）和trans-mQTL（SNP与CpG位点距离较远或位于不同染色体）。cis-mQTL的效应通常较强且容易检测，而trans-mQTL的效应较弱但可能揭示跨染色体的调控网络。当你需要理解遗传变异如何通过影响甲基化进而影响表型时，mQTL分析提供了关键的中间环节。 甲基化与转录因子结合位点关系分析探讨甲基化如何影响转录因子的DNA结合活性。许多转录因子含有CpG在内的识别序列，甲基化可以直接阻碍这些转录因子的结合。例如，CTCF是一个对甲基化高度敏感的绝缘子蛋白，其结合位点的甲基化会导致CTCF结合丧失，进而影响染色质环结构和基因表达。另一方面，一些蛋白（如MeCP2、Kaiso等）特异性地结合甲基化的CpG位点，招募抑制性复合物。当你研究特定转录因子的调控机制时，检查其结合位点的甲基化状态可以提供重要的功能线索。 等位基因特异性甲基化（ASM）是指基因组上同一CpG位点在两条同源染色体上具有不同甲基化状态的现象。ASM包括基因组印记（imprinting）相关的ASM和遗传变异导致的ASM。基因组印记是由亲本来源决定的ASM，印记基因只表达来自特定亲本的等位基因。遗传变异导致的ASM通常与附近的SNP相关，即cis-mQTL效应。当你使用长读长测序或等位基因特异性比对时，可以检测到ASM现象，这对于理解等位基因特异性表达和印记基因的调控至关重要。 甲基化年龄推断（DNAm年龄，时钟）是利用特定CpG位点的甲基化水平来预测个体生物学年龄的技术。DNA甲基化随年龄发生规律性变化，某些CpG位点的甲基化水平与年龄高度相关。Horvath时钟使用353个CpG位点，适用于多种组织类型；Hannum时钟使用71个CpG位点，主要基于血液样本开发。DNAm年龄与实际年龄的差值被称为年龄加速（AgeAccel），正值的年龄加速意味着个体的表观遗传年龄大于实际年龄，可能与加速衰老或疾病风险增加相关。当你研究衰老相关疾病或环境暴露对衰老的影响时，表观遗传时钟提供了一个客观的生物学年龄度量。 甲基化与疾病关联（EWAS，表观组关联研究）是将甲基化作为暴露因素或中间表型，研究其与疾病或性状之间关联的方法。EWAS的设计类似于GWAS，但以CpG位点为单位检验甲基化水平与表型的关联。EWAS面临的主要挑战包括细胞类型异质性（不同细胞类型的甲基化模式差异很大）、多重检验负担（数百万个CpG位点）和因果方向的确定（甲基化变化是疾病的原因还是结果）。当你设计EWAS研究时，需要特别注意样本量估算、细胞组成校正和结果验证策略。 组织特异性甲基化模式是指不同组织类型具有独特的甲基化特征。每种组织都有其特征性的甲基化图谱，这些特征主要反映了组织特异性基因的调控状态。组织特异性差异甲基化区域（tDMR）可以用来推断未知样本的组织来源，这在法医学和液体活检中有重要应用。当你分析来自不同组织的甲基化数据时，需要首先确认样本的组织来源，因为组织混杂是EWAS中最常见的混杂因素之一。 .. code-block:: r # 启动子甲基化与基因表达相关性分析 set.seed(42) n_genes <- 100 promoter_meth <- runif(n_genes, 0, 1) gene_expr <- -0.6 * promoter_meth + rnorm(n_genes, 0.5, 0.2) gene_expr <- pmax(gene_expr, 0) cor_test <- cor.test(promoter_meth, gene_expr, method = "spearman") cat(sprintf("Spearman相关系数: %.4f, p值: %.4e\n", cor_test$estimate, cor_test$p.value)) # 简单的mQTL分析模拟 genotypes <- sample(c(0, 1, 2), n_genes, replace = TRUE, prob = c(0.5, 0.4, 0.1)) meth_levels <- 0.3 + 0.15 * genotypes + rnorm(n_genes, 0, 0.1) mQTL_result <- lm(meth_levels ~ genotypes) summary(mQTL_result) 6.2 染色质可及性分析 --------------------- 染色质可及性是指染色质DNA被蛋白质和酶（如转录因子、RNA聚合酶等）接触和结合的能力。在真核细胞中，DNA缠绕在组蛋白八聚体上形成核小体，核小体进一步折叠形成高级染色质结构。当染色质处于开放状态时，DNA暴露出来，转录因子可以结合并调控基因表达；当染色质处于紧密状态时，DNA被包裹在核小体内部，转录因子无法接触。染色质可及性分析的目标就是绘制基因组中开放染色质区域的图谱，识别活跃的调控元件，并推断转录因子的结合活性。 6.2.1 检测技术 ^^^^^^^^^^^^^^^ ATAC-seq（Assay for Transposase-Accessible Chromatin）是目前最主流的染色质可及性检测技术。它利用Tn5转座酶优先切割和标记开放染色质区域的DNA，同时加入测序接头，一步完成开放区域的切割和文库构建。ATAC-seq的优势在于实验流程简单快速（仅需约3小时）、所需起始细胞量少（低至500个细胞），且能够同时提供核小体定位信息。当你需要快速获得样本的染色质可及性图谱时，ATAC-seq是首选方法。 DNase-seq（DNase I超敏感位点测序）利用DNase I酶优先切割开放染色质区域的特性来识别DNase I超敏感位点（DHS）。DNase-seq是最早用于全基因组染色质可及性检测的技术之一，ENCODE项目产生了大量的人类和小鼠DNase-seq数据作为参考。与ATAC-seq相比，DNase-seq需要更多的起始细胞和更复杂的实验流程，但它在某些特定应用中仍有价值，特别是当你需要与ENCODE的历史数据进行比较时。 MNase-seq（微球菌核酸酶测序）利用微球菌核酸酶（MNase）优先消化核小体之间的连接DNA，保留核小体包裹的DNA片段。与ATAC-seq和DNase-seq检测开放区域不同，MNase-seq主要用来确定核小体的精确位置。通过分析MNase消化后保留的约147bp DNA片段，可以在全基因组范围内绘制核小体占据图谱。当你需要研究核小体定位与基因调控的关系时，MNase-seq提供了独特的核小体分辨率信息。 NOMe-seq（核小体占据与甲基化测序）结合了GpC甲基转移酶和重亚硫酸盐测序，同时检测核小体占据和内源性DNA甲基化状态。GpC甲基转移酶在开放染色质区域对GpC二核苷酸进行甲基化，而核小体包裹的GpC则被保护。通过重亚硫酸盐测序同时读取GpC甲基化（反映染色质可及性）和CpG甲基化（反映内源性甲基化），NOMe-seq可以在同一分子上同时获得两种表观遗传信息。当你需要研究染色质可及性与DNA甲基化的关系时，NOMe-seq提供了单分子分辨率的信息。 scATAC-seq（单细胞ATAC-seq）将ATAC-seq技术应用于单个细胞，能够揭示细胞群体中染色质可及性的异质性。scATAC-seq的数据具有极度稀疏的特点，每个细胞通常只捕获基因组中少数开放区域的信号。常用的scATAC-seq平台包括10X Genomics Chromium和Sci-ATAC-seq等。当你需要识别细胞类型特异的调控元件或构建单细胞染色质可及性图谱时，scATAC-seq是必要的工具。 基于转座酶的索引技术（sci-ATAC）通过组合索引策略大幅提高了单细胞ATAC-seq的通量。sci-ATAC使用多轮分池和标记，以组合方式为每个细胞赋予独特的条形码，从而在无需物理分离每个细胞的情况下实现单细胞分辨率的染色质可及性分析。当你需要处理大量细胞但成本有限时，sci-ATAC提供了高性价比的解决方案。 空间ATAC-seq是空间组学技术在染色质可及性领域的延伸，它能够在保留组织空间位置信息的同时检测染色质可及性。空间ATAC-seq使得研究者可以在组织切片的原位观察不同区域的染色质开放状态，对于理解组织微环境中细胞状态的异质性具有重要意义。当你需要将染色质可及性信息与组织空间结构关联时，空间ATAC-seq提供了独特的视角。 6.2.2 ATAC-seq 数据分析 ^^^^^^^^^^^^^^^^^^^^^^^^ ATAC-seq数据的比对与常规ChIP-seq或RNA-seq有所不同。由于Tn5转座酶切割的DNA片段包括核小体游离片段和核小体结合片段，比对时需要使用能够处理可变长度片段的比对器。Bowtie2是最常用的ATAC-seq比对工具，它支持局部比对和末端到末端比对模式，可以有效地将ATAC-seq读段比对到参考基因组。minimap2也可以用于ATAC-seq比对，特别是在处理长读长ATAC-seq数据时。当你设置比对参数时，建议使用--very-sensitive模式并允许最多2个错配，以平衡灵敏度和特异性。 .. code-block:: bash # ATAC-seq比对流程 # 1. 比对到参考基因组 bowtie2 -x /path/to/genome_index \ -1 sample_R1.fq.gz -2 sample_R2.fq.gz \ --very-sensitive -X 2000 \ -p 8 -S sample.sam # 2. 转换为BAM并排序 samtools view -bS sample.sam | samtools sort -o sample.sorted.bam samtools index sample.sorted.bam # 3. 去除线粒体reads samtools idxstats sample.sorted.bam | cut -f 1 | grep -v "chrM" | \ xargs samtools view -b sample.sorted.bam > sample.noMT.bam samtools index sample.noMT.bam 去除线粒体reads是ATAC-seq数据处理中特别重要的一步。线粒体DNA没有核小体结构，完全处于开放状态，因此Tn5转座酶会大量切割线粒体DNA。在典型的ATAC-seq数据中，线粒体reads可能占总reads的20-60%，如果不加以去除会严重影响有效数据的比例。去除线粒体reads的方法包括使用samtools过滤chrM染色体的比对结果，或在比对前将线粒体基因组从参考中排除。当你发现线粒体reads比例异常高（超过50%）时，可能需要检查实验质量或考虑使用线粒体DNA含量较低的细胞类型。 片段长度分布分析是ATAC-seq数据质量控制的核心指标。ATAC-seq产生的DNA片段长度分布呈现特征性的周期性模式：核小体游离片段（nucleosome-free fragment，NFR）长度小于100bp，代表转录因子结合的开放区域；单核小体片段长度约180-247bp，代表被单个核小体包裹的DNA；双核小体片段长度约315-473bp；三核小体片段长度约558-615bp。一个高质量的ATAC-seq实验应该显示出清晰的核小体游离峰和单核小体峰。当你观察到片段长度分布缺乏周期性模式或核小体游离峰不明显时，可能提示实验质量不佳。 .. code-block:: r # ATAC-seq片段长度分布分析 library(ggplot2) # 模拟片段长度分布数据 fragment_lengths <- c( rnorm(5000, mean = 50, sd = 20), rnorm(3000, mean = 200, sd = 30), rnorm(1500, mean = 400, sd = 40), rnorm(500, mean = 600, sd = 50) ) fragment_lengths <- abs(fragment_lengths) fragment_lengths <- fragment_lengths[fragment_lengths > 0 & fragment_lengths < 800] df <- data.frame(length = fragment_lengths) ggplot(df, aes(x = length)) + geom_histogram(binwidth = 5, fill = "steelblue", alpha = 0.7) + geom_vline(xintercept = c(100, 200, 400, 600), linetype = "dashed", color = "red") + annotate("text", x = c(50, 200, 400, 600), y = rep(Inf, 4), vjust = 2, label = c("NFR", "单核小体", "双核小体", "三核小体")) + labs(x = "片段长度(bp)", y = "计数", title = "ATAC-seq片段长度分布") + theme_minimal() 插入位点偏移校正是ATAC-seq数据分析中一个微妙但重要的步骤。Tn5转座酶以二聚体形式结合DNA，两个单体在插入位点各偏移9bp，导致实际切割位置比Tn5结合中心偏移+4bp和-5bp。为了获得转录因子结合位点的精确位置，需要对比对结果进行偏移校正：正链reads的5'端向3'方向移动+4bp，负链reads的5'端向5'方向移动-5bp。这种校正对于后续的足迹分析和基序富集分析尤为重要。当你进行转录因子足迹分析时，Tn5偏移校正是必须的预处理步骤。 峰值调用（peak calling）是识别基因组中染色质开放区域的关键步骤。MACS2是最常用的峰值调用工具，对于ATAC-seq数据建议使用--nomodel --shift -100 --extsize 200参数，这些参数适应ATAC-seq片段的特征。HOMER的findPeaks命令也可以用于ATAC-seq峰值调用，它提供了丰富的注释功能。Genrich是一个专门为ATAC-seq设计的峰值调用器，它可以直接处理BAM文件并考虑Tn5偏移。当你调用ATAC-seq峰时，需要确保使用适当的对照组或使用输入DNA作为背景，虽然ATAC-seq通常不需要专门的input对照。 峰的质量评估使用多个指标来衡量ATAC-seq实验和峰值调用的质量。FRiP（Fraction of Reads in Peaks）是指落在峰值区域内的reads占总reads的比例，是衡量信噪比的关键指标，高质量的ATAC-seq数据FRiP应大于0.2。峰值数（Number of Peaks）反映了检测到的开放区域数量，人类细胞通常应检测到30000-100000个峰。NRF（Non-redundant Fraction）是衡量文库复杂度的指标，计算方式为唯一比对位置数除以总比对reads数。PBC（PCR Bottleneck Coefficient）是衡量PCR扩增偏倚的指标，计算方式为最大非冗余位置数除以非冗余位置数。当你发现FRiP过低或峰值数异常少时，可能需要检查实验质量或调整峰值调用参数。 峰注释是将检测到的开放区域映射到基因组功能元件的过程。注释通常关注峰与最近转录起始位点（TSS）的距离，以及峰落在何种基因组元件中（启动子、外显子、内含子、基因间区、增强子等）。ATAC-seq峰在启动子区域的富集反映了活跃转录基因的开放染色质状态，而远距离的峰可能对应于增强子等远端调控元件。当你需要将染色质可及性与基因调控联系起来时，峰注释是连接开放区域与靶基因的桥梁。 差异可及性分析比较不同条件之间染色质开放程度的变化。最常用的方法是将峰值区域作为特征，统计每个峰区域内的reads计数，然后使用DESeq2或edgeR等差异分析工具进行统计检验。DiffBind是专门为染色质数据设计的差异分析R包，它封装了DESeq2和edgeR的分析流程，并提供了峰值合并和注释功能。当你进行差异可及性分析时，需要注意峰值区域的定义在不同样本之间可能不同，通常需要先合并所有样本的峰值，然后在统一的峰值集上进行计数和比较。 .. code-block:: r # 使用DiffBind进行差异可及性分析 library(DiffBind) # 创建样本信息表 samples <- dba.sampleSheet( SampleID = c("ctrl1", "ctrl2", "ctrl3", "treat1", "treat2", "treat3"), Condition = c("Control", "Control", "Control", "Treatment", "Treatment", "Treatment"), bamReads = c("ctrl1.bam", "ctrl2.bam", "ctrl3.bam", "treat1.bam", "treat2.bam", "treat3.bam"), Peaks = c("ctrl1_peaks.narrowPeak", "ctrl2_peaks.narrowPeak", "ctrl3_peaks.narrowPeak", "treat1_peaks.narrowPeak", "treat2_peaks.narrowPeak", "treat3_peaks.narrowPeak") ) # 创建DBA对象并计算一致性峰值集 dba_obj <- dba(sampleSheet = samples) dba_obj <- dba.count(dba_obj) # 差异分析 dba_obj <- dba.contrast(dba_obj, categories = DBA_CONDITION, minMembers = 2) dba_obj <- dba.analyze(dba_obj) # 提取差异可及性区域 diff_regions <- dba.report(dba_obj, th = 0.05) head(diff_regions) 基序富集分析是在差异可及性区域中寻找转录因子结合基序的富集模式，从而推断哪些转录因子可能驱动了观察到的染色质可及性变化。HOMER的findMotifsGenome命令可以快速在基因组区域中搜索已知的和de novo的转录因子结合基序。MEME套件提供了更全面的基序分析工具，包括MEME（de novo基序发现）、AME（已知基序富集检验）和CentriMo（中心富集分析）。当你发现一组差异可及性区域时，基序富集分析可以帮助你识别潜在的调控转录因子，缩小后续实验验证的范围。 足迹分析（footprinting）是在ATAC-seq数据中检测转录因子结合位点上reads覆盖度下降的模式。当转录因子结合到DNA上时，它会保护该区域免受Tn5切割，在开放区域的信号中形成一个局部的"凹陷"，即足迹。HINT-ATAC是一个专门为ATAC-seq数据设计的足迹分析工具，它考虑了Tn5偏移和核小体信号的影响。TOBIAS是另一个足迹分析工具，它通过比较不同条件下的足迹深度变化来推断转录因子活性的变化。当你需要从ATAC-seq数据中推断转录因子的实际结合状态而不仅仅是基序的存在时，足迹分析提供了更直接的信息。 核小体定位分析利用ATAC-seq数据中的核小体相关片段来推断核小体的精确位置。NucTools是一个用于核小体定位分析的工具，它可以分析核小体占据图谱和核小体定位的周期性模式。ATAC-seq数据中单核小体片段（约180-247bp）的末端对应核小体的边界，通过分析这些片段的分布可以推断核小体的位置和占据频率。当你研究核小体定位与基因表达调控的关系时，ATAC-seq数据提供了比MNase-seq更简便的核小体信息来源。 scATAC-seq分析流程需要处理数据的极度稀疏性。Signac是Seurat生态系统的单细胞ATAC-seq分析R包，它提供了从预处理到降维聚类的完整流程，并且可以与Seurat的scRNA-seq分析无缝整合。ArchR是一个高性能的单细胞染色质可及性分析框架，支持大规模数据集的处理，提供了峰值调用、基序富集、轨迹推断等功能。Cicero专注于分析单细胞染色质可及性数据中的共可及性（co-accessibility）模式，通过识别远距离开放区域之间的共可及性关系来推断增强子-启动子的连接。当你需要将scATAC-seq与scRNA-seq数据整合分析时，Signac+Seurat或ArchR是推荐的选择。 .. code-block:: r # 使用Signac进行scATAC-seq分析 library(Signac) library(Seurat) # 读取10X Genomics scATAC-seq数据 counts <- Read10X_h5("filtered_peak_bc_matrix.h5") fragment_file <- "fragments.tsv.gz" # 创建ChromatinAssay对象 chrom_assay <- CreateChromatinAssay( counts = counts, sep = c(":", "-"), fragments = fragment_file ) # 创建Seurat对象 atac_obj <- CreateSeuratObject( counts = chrom_assay, assay = "peaks" ) # 计算核小体信号和TSS富集分数 atac_obj <- NucleosomeSignal(atac_obj) atac_obj <- TSSEnrichment(atac_obj) # 质量过滤 atac_obj <- subset( atac_obj, nCount_peaks > 1000 & nCount_peaks < 50000 & nucleosome_signal < 2 & TSS.enrichment > 2 ) # 降维和聚类 atac_obj <- RunTFIDF(atac_obj) atac_obj <- FindTopFeatures(atac_obj, min.cutoff = "q0") atac_obj <- RunSVD(atac_obj) atac_obj <- RunUMAP(atac_obj, reduction = "lsi", dims = 2:30) atac_obj <- FindNeighbors(atac_obj, reduction = "lsi", dims = 2:30) atac_obj <- FindClusters(atac_obj, resolution = 0.5) 6.3 组蛋白修饰分析 ------------------ 组蛋白修饰是表观遗传调控的另一个核心层次。组蛋白的N端尾巴可以发生多种共价修饰，包括甲基化、乙酰化、磷酸化、泛素化等，这些修饰改变了染色质的结构和蛋白质相互作用，从而影响基因的表达状态。不同类型的组蛋白修饰标记了不同的染色质状态：某些修饰（如H3K4me3、H3K27ac）标记活跃的基因调控区域，而另一些修饰（如H3K27me3、H3K9me3）则标记沉默的染色质区域。组蛋白修饰分析的主要技术是染色质免疫沉淀测序（ChIP-seq），它通过特异性抗体富集目标修饰的基因组区域，然后通过测序确定这些区域的位置。 6.3.1 ChIP-seq 基础 ^^^^^^^^^^^^^^^^^^^^^ 染色质免疫沉淀测序（ChIP-seq）的原理是利用特异性抗体识别并结合目标蛋白或修饰，然后通过免疫沉淀将蛋白质-DNA复合物从细胞裂解物中分离出来，最后对富集的DNA片段进行测序。ChIP-seq实验的成功与否在很大程度上取决于抗体的质量和实验条件的优化。当你计划一个ChIP-seq实验时，需要仔细考虑抗体选择、交联条件、对照设置和测序深度等因素。 抗体选择是ChIP-seq实验中最关键的决策之一。抗体的特异性直接决定了实验结果的可靠性，非特异性抗体会导致背景信号升高和假阳性峰的出现。在选择抗体时，应该优先考虑经过ChIP-seq验证的抗体，并查阅ENCODE等数据库中该抗体的历史使用数据。抗体的批次间差异也是一个需要关注的问题，不同批次的同一抗体可能表现不同。当你选择抗体时，建议参考ENCODE推荐的抗体列表，并在可能的情况下使用已发表研究中验证过的抗体。 交联条件的选择取决于目标蛋白的类型。XChIP（交联ChIP）使用甲醛等交联剂将蛋白质与DNA固定在一起，适用于大多数转录因子和组蛋白修饰的研究。Native ChIP（非交联ChIP）不使用交联剂，直接利用微球菌核酸酶消化染色质，适用于组蛋白修饰的研究，因为组蛋白与DNA的结合比转录因子更为稳定。当你研究组蛋白修饰时，Native ChIP通常能提供更低的背景和更精确的信号；而研究转录因子时，XChIP是必需的，因为转录因子与DNA的结合较弱，不交联会在操作过程中丢失。 输入对照（input DNA）是ChIP-seq实验中必不可少的阴性对照。input DNA是未经免疫沉淀的DNA样本，它反映了染色质片段化和测序过程中的偏倚，如开放染色质区域更容易被片段化和测序。在峰值调用时，input对照用于估计背景信号水平，从而区分真实的富集信号和背景噪音。当你进行ChIP-seq实验时，每个样本都应配备相应的input对照，否则峰值调用的可靠性会大打折扣。 ChIP-seq实验设计还需要考虑峰数量期望和测序深度。不同类型的组蛋白修饰产生不同特征的峰：转录因子和活性启动子标记（如H3K4me3）产生尖锐的窄峰，通常需要20-40M reads即可获得良好的覆盖；而广泛修饰（如H3K27me3、H3K36me3）覆盖较大的基因组区域，产生宽峰，需要更高的测序深度（40-60M reads）。当你设计实验时，需要根据目标修饰的峰类型确定合适的测序深度。 6.3.2 ChIP-seq 数据分析 ^^^^^^^^^^^^^^^^^^^^^^^^^ ChIP-seq数据的比对通常使用Bowtie2或BWA进行。比对参数的设置与ATAC-seq类似，使用--very-sensitive模式以获得最佳的比对效果。比对后需要检查比对率、唯一比对率等质量指标，比对率过低可能提示样本质量问题或参考基因组选择不当。 重复去除是ChIP-seq数据处理中的重要步骤。PCR扩增产生的重复reads不代表真实的生物学信号，需要被标记或去除。Picard的MarkDuplicates工具是最常用的去重工具，它可以标记PCR重复而不删除它们，保留信息用于后续的质量评估。当你处理ChIP-seq数据时，去重后的唯一比对reads数是评估文库复杂度的重要指标。 .. code-block:: bash # ChIP-seq数据处理流程 # 1. 比对 bowtie2 -x /path/to/genome_index \ -U chip_sample.fq.gz \ --very-sensitive -p 8 \ -S chip_sample.sam # 2. 转换并排序 samtools view -bS chip_sample.sam | samtools sort -o chip_sample.sorted.bam samtools index chip_sample.sorted.bam # 3. 去除PCR重复 picard MarkDuplicates \ INPUT=chip_sample.sorted.bam \ OUTPUT=chip_sample.dedup.bam \ METRICS_FILE=chip_sample.metrics.txt \ REMOVE_DUPLICATES=true samtools index chip_sample.dedup.bam 峰值调用是ChIP-seq数据分析的核心步骤，其目标是从比对数据中识别基因组上被目标蛋白或修饰显著富集的区域。MACS2是最广泛使用的峰值调用工具，对于转录因子和H3K4me3等产生窄峰的修饰，使用默认的narrow peak模式；对于H3K27me3和H3K36me3等产生宽峰的修饰，使用--broad参数调用broad peak。SICER专门设计用于检测广泛组蛋白修饰的富集区域，它利用空间聚类算法识别跨越数kb到数十kb的富集区域。HOMER的findPeaks命令提供了多种峰值调用算法和丰富的注释功能。MUSIC是另一个支持多尺度峰值检测的工具，适合于同时存在窄峰和宽峰的数据。 .. code-block:: bash # MACS2峰值调用 # 窄峰（转录因子，H3K4me3等） macs2 callpeak -t chip_sample.dedup.bam \ -c input_sample.dedup.bam \ -f BAM -g hs \ -n sample_narrow \ --outdir peaks/ \ -q 0.01 # 宽峰（H3K27me3，H3K36me3等） macs2 callpeak -t chip_sample.dedup.bam \ -c input_sample.dedup.bam \ -f BAM -g hs \ -n sample_broad \ --outdir peaks/ \ --broad -q 0.1 峰值质量控制使用ENCODE项目制定的标准来评估ChIP-seq数据的质量。FRiP（Fraction of Reads in Peaks）应大于0.01，即至少1%的reads落在峰值区域内，这个阈值对于组蛋白修饰可能更高。NSC（Normalized Strand Cross-correlation）和RSC（Relative Strand Cross-correlation）是基于链间交叉相关分析的质量指标，NSC应大于1.05，RSC应大于0.8。当你的数据未达到这些标准时，可能需要重新评估实验质量或调整分析参数。 峰可视化帮助研究者直观地理解组蛋白修饰的分布模式。IGV可以逐个基因地查看ChIP-seq信号的分布，适合于检查个别候选基因区域的修饰状态。deepTools提供了批量生成热图和metaplot的功能，可以在TSS等参考点周围绘制所有基因的信号分布，揭示组蛋白修饰的整体模式。当你需要在大量基因上总结修饰分布特征时，deepTools的热图和metaplot比逐基因查看更有效率。 .. code-block:: bash # deepTools可视化 # 生成bigWig文件 bamCoverage -b chip_sample.dedup.bam \ -o chip_sample.bw \ --binSize 10 \ --normalizeUsing RPGC \ --effectiveGenomeSize 2913022398 \ --extendReads # 绘制TSS周围的metaplot computeMatrix reference-point \ -S chip_sample.bw \ -R genes.bed \ --referencePoint TSS \ -b 3000 -a 3000 \ -o matrix_TSS.gz plotProfile -m matrix_TSS.gz \ -o profile_TSS.png \ --perGroup # 绘制热图 plotHeatmap -m matrix_TSS.gz \ -o heatmap_TSS.png 差异结合分析比较不同条件之间组蛋白修饰或转录因子结合的差异。DiffBind是专门为ChIP-seq差异分析设计的R包，它使用DESeq2或edgeR作为统计后端，对峰值区域的reads计数进行差异检验。MAnorm通过比较两个样本之间的信号强度来识别差异结合区域，不需要生物学重复。CSAR是另一个差异分析工具，基于泊松分布模型检验富集信号的差异。当你有生物学重复时，DiffBind是推荐的选择；当只有单个样本时，MAnorm可以作为替代方案。 重塑特定修饰的特征是理解组蛋白修饰功能的重要步骤。H3K4me3在活跃转录基因的启动子区域形成尖锐的峰，其信号强度与基因表达水平正相关。H3K27ac在活跃增强子和启动子区域富集，是区分活跃增强子与poised增强子的关键标记。H3K36me3沿着活跃转录的基因体延伸，其分布范围与转录延伸距离相关。当你分析特定修饰的ChIP-seq数据时，理解其典型的分布特征有助于判断数据质量和生物学意义。 宽峰与尖峰的区别是ChIP-seq数据分析中需要特别注意的问题。窄峰（narrow peak）通常由转录因子和某些组蛋白修饰（如H3K4me3、H3K27ac）产生，峰的宽度通常在数百bp范围内，信号集中且尖锐。宽峰（broad peak）由广泛分布的组蛋白修饰（如H3K27me3、H3K36me3、H3K9me3）产生，覆盖数kb到数十kb的基因组区域，信号相对平缓。这两种峰类型需要不同的峰值调用策略和参数设置，MACS2的narrow peak模式和broad peak模式分别适用于这两种情况。当你同时分析多种组蛋白修饰时，需要根据每种修饰的峰类型选择相应的分析策略。 超级增强子鉴定是ChIP-seq分析中的一个重要应用。超级增强子是一类覆盖范围大、信号强度高的增强子簇，它们在细胞身份决定和疾病发生中发挥着关键作用。ROSE（Rank Ordering of Super-Enhancers）算法利用H3K27ac的ChIP-seq信号来鉴定超级增强子，其基本思路是将相邻的H3K27ac峰缝合为增强子区域，然后按照信号强度排序，通过斜率变化点区分超级增强子和普通增强子。当你研究细胞类型特异性调控或肿瘤中的增强子重编程时，超级增强子鉴定可以揭示关键的调控元件。 共结合分析通过比较多个ChIP-seq数据集的峰值重叠来研究不同蛋白或修饰之间的协同关系。例如，同时分析转录因子和组蛋白修饰的ChIP-seq数据，可以发现转录因子结合位点与特定组蛋白修饰的共定位模式。共结合分析通常使用bedtools intersect等工具计算峰值重叠，并使用统计检验评估重叠的显著性。当你需要理解多个调控因子之间的协同调控机制时，共结合分析提供了重要的线索。 6.3.3 组蛋白修饰的功能解读 ^^^^^^^^^^^^^^^^^^^^^^^^^^^ 活性启动子标记（H3K4me3和H3K9ac）是识别活跃转录基因的关键表观遗传特征。H3K4me3在活跃基因的转录起始位点（TSS）附近形成特征性的峰，通常跨越TSS上下游约1-2kb的区域。H3K9ac与H3K4me3经常共定位，两者的存在标志着启动子处于开放和活跃的状态。当你需要判断一个基因是否在特定细胞类型中活跃表达时，H3K4me3和H3K9ac的ChIP-seq信号是比DNA甲基化更直接的活性标记。 活性增强子标记（H3K27ac和H3K4me1）是识别活跃增强子的核心特征。H3K4me1在增强子区域富集，但单独的H3K4me1不能区分活跃增强子和poised增强子。H3K27ac是区分活跃增强子（H3K27ac阳性）和poised增强子（H3K27ac阴性、H3K4me1阳性）的关键标记。当你需要识别细胞类型特异性的活跃增强子时，H3K27ac的ChIP-seq数据是最重要的依据。 转录延伸标记（H3K36me3）沿着活跃转录的基因体分布，从TSS下游延伸到转录终止位点。H3K36me3的分布范围与基因的长度和转录活性相关，长基因和高表达基因具有更强的H3K36me3信号。H3K36me3还与剪接调控和DNA甲基化维持相关，Setd2介导的H3K36me3可以招募DNMT3B进行基因体甲基化。当你需要评估基因的转录延伸状态时，H3K36me3是最合适的标记。 多梳抑制标记（H3K27me3和H3K9me2/me3）标记了沉默的染色质区域。H3K27me3由多梳抑制复合物2（PRC2）催化，覆盖大范围的基因组区域，在发育调控基因的启动子区域尤为常见。H3K9me2和H3K9me3主要与组成型异染色质相关，标记着丝粒、端粒和重复元件等区域。H3K27me3和H3K4me3同时存在的启动子被称为二价启动子（bivalent promoter），在胚胎干细胞中常见，代表基因处于"准备"状态，可以快速激活或永久沉默。当你研究干细胞分化或肿瘤中的基因沉默机制时，多梳标记的分析尤为重要。 基因体甲基化与H3K36me3的关联是一个有趣的表观遗传交叉调控现象。基因体内的CpG甲基化与H3K36me3呈正相关，活跃转录的基因同时具有高水平的基因体甲基化和H3K36me3。这种关联的分子机制是H3K36me3可以招募DNMT3B到基因体区域进行甲基化，而基因体甲基化又可以抑制异常的基因内转录起始。当你同时分析甲基化和组蛋白修饰数据时，基因体区域H3K36me3与甲基化的正相关是一个可预期的特征。 染色质状态分割与注释是将多种组蛋白修饰的ChIP-seq数据整合，将基因组划分为具有不同功能状态的区域。ChromHMM是一种基于隐马尔可夫模型（HMM）的方法，它学习多种组蛋白修饰的组合模式，将基因组分割为不同数目的状态（如15状态或25状态模型），每种状态对应特定的染色质功能（如活跃启动子、强增强子、弱增强子、转录延伸、多梳抑制等）。Segway使用动态贝叶斯网络进行类似的分割，但建模方式不同。HMMstate是另一种染色质状态分割工具。当你拥有多种组蛋白修饰的ChIP-seq数据时，染色质状态分割可以提供基因组功能的全面注释。 染色质状态图是ENCODE项目定义的标准化的基因组功能注释。ENCODE 25状态模型和15状态模型使用不同组合的组蛋白修饰来定义染色质状态，每种状态有明确的生物学含义和命名。例如，活跃启动子状态的特征是H3K4me3+H3K27ac+H3K9ac，强增强子状态的特征是H3K4me1+H3K27ac，多梳抑制状态的特征是H3K27me3。当你需要将基因组功能注释与自己的数据关联时，ENCODE的染色质状态图提供了标准化的参考。 .. code-block:: r # 使用ChromHMM进行染色质状态分割（命令行示例） # ChromHMM需要Java运行环境 # 1. 准备输入文件：将每个样本的BAM转换为二值化的信号文件 # BinarizeBam命令将BAM文件转换为ChromHMM格式 # 2. 学习染色质状态模型 # LearnModel命令执行HMM模型训练 # system("java -jar ChromHMM.jar LearnModel -p 8 binarized_dir output_dir 15 hg38") # 3. 结果解读 # emissions_15.txt: 各状态的组蛋白修饰发射概率 # transitions_15.txt: 状态之间的转移概率 # sample_15_segments.bed: 基因组分割结果 6.4 染色质互作分析 ------------------ 染色质在细胞核中并非线性排列，而是通过折叠和环化形成复杂的三维结构。染色质互作分析研究的是基因组上远距离位点之间的物理接触，这些互作对于基因调控至关重要——增强子可以跨越数十万甚至数百万碱基对与启动子形成环状结构来激活基因表达。染色质构象捕获技术及其衍生方法使得全基因组范围内的染色质互作检测成为可能，为理解基因调控的空间维度提供了强有力的工具。 6.4.1 染色质构象捕获技术 ^^^^^^^^^^^^^^^^^^^^^^^^^ 3C（染色质构象捕获）是最早的染色质互作检测技术，它通过甲醛交联固定细胞内蛋白质-DNA复合物，然后用限制性内切酶消化，在近距离连接后检测特定的一对基因组位点之间的互作频率。3C是"一对一"的检测模式，只能检测预先选定的两个位点之间的互作，适合于验证已知的增强子-启动子互作。 4C（环状染色质构象捕获）采用"一对全"的检测模式，以一个特定的位点（viewpoint）为锚点，检测其与基因组上所有其他位点的互作。4C通过反向PCR扩增viewpoint周围的连接片段，然后进行高通量测序。当你需要研究一个特定基因或增强子与全基因组的互作网络时，4C提供了全面的信息。 5C（3C-碳拷贝）采用"多对多"的检测模式，通过设计覆盖目标区域的引物池来同时检测多个位点之间的互作。5C适合于研究特定基因组区域（如一个基因簇或一个拓扑关联结构域）内部的所有互作关系，但引物设计的复杂性限制了其应用范围。 Hi-C是"全对全"的染色质互作检测技术，它使用生物素标记的核苷酸富集连接片段，然后进行全基因组测序，理论上可以检测基因组上所有位点对之间的互作频率。Hi-C数据通常以接触矩阵（contact matrix）的形式呈现，矩阵中每个元素表示对应两个基因组区域之间的互作频率。Hi-C是研究全基因组染色质组织结构的最强大工具，当你需要绘制完整的染色质三维结构图谱时，Hi-C是首选方法。 Capture Hi-C结合了Hi-C和靶向捕获技术，通过设计探针捕获包含特定区域（如所有基因启动子）的连接片段，从而以更高的分辨率和深度检测目标区域的互作。Capture Hi-C特别适合于研究启动子与远端调控元件之间的互作，在疾病遗传学研究中被广泛用于将GWAS信号与靶基因连接起来。 ChIA-PET（染色质免疫沉淀-配对末端标签测序）结合了ChIP和Hi-C的原理，先通过免疫沉淀富集特定蛋白（如RNA聚合酶II、CTCF等）介导的染色质互作，然后检测这些蛋白相关的互作网络。ChIA-PET提供了蛋白特异性的染色质互作信息，当你需要研究特定蛋白在染色质三维结构中的作用时，ChIA-PET是理想的选择。 HiChIP和PLAC-seq是ChIA-PET的改进版本，它们在Hi-C文库制备后进行免疫沉淀，而非在连接前进行，这使得实验流程更简单且信噪比更高。HiChIP和PLAC-seq产生的数据质量通常优于ChIA-PET，且需要的细胞量更少。当你需要研究特定蛋白相关的染色质互作但样本量有限时，HiChIP和PLAC-seq比ChIA-PET更实用。 Micro-C使用微球菌核酸酶代替限制性内切酶消化染色质，实现了核小体分辨率的染色质互作检测。与Hi-C相比，Micro-C能够检测更精细的染色质结构，包括相邻核小体之间的互作。当你需要研究核小体尺度的染色质组织时，Micro-C提供了Hi-C无法达到的分辨率。 SPRITE（Split-Pool Recognition of Interactions by Tag Extension）是一种不依赖连接反应的染色质互作检测技术，它通过多轮分池和标记来识别同一复合物中的多个DNA分子。SPRITE可以检测包含多个位点的复合互作，而传统的Hi-C只能检测两两互作。当你需要研究涉及多个基因组位点的复杂染色质组装时，SPRITE提供了独特的多路互作信息。 单细胞Hi-C将染色质互作检测应用于单个细胞，能够揭示细胞群体中染色质三维结构的异质性。单细胞Hi-C的数据极度稀疏，每个细胞只能捕获少量的互作信息，但通过聚合大量细胞的数据可以重建染色质结构的细胞类型特异性特征。 6.4.2 Hi-C 数据分析 ^^^^^^^^^^^^^^^^^^^^^ Hi-C数据的比对需要将双端reads分别比对到参考基因组，然后根据比对位置确定每对reads代表的基因组互作。HiC-Pro是一个完整的Hi-C数据处理流程，包括比对、过滤、接触矩阵构建和归一化。HiCAlign专注于Hi-C数据的快速比对。Juicer是Aiden实验室开发的Hi-C分析流程，提供了从原始数据到可视化的端到端解决方案。distiller是Juicer的下一代版本，使用更模块化的设计。当你需要处理大规模Hi-C数据时，HiC-Pro和Juicer是两个最成熟的选择。 .. code-block:: bash # 使用Juicer处理Hi-C数据 # Juicer自动化流程 juicer.sh -g hg38 -s MboI -z /path/to/genome.fa \ -p /path/to/chrom.sizes \ -y /path/to/restriction_sites.txt \ -d /path/to/fastq_dir/ \ -D /path/to/juicer_dir/ \ -t 16 # 使用HiC-Pro HiC-Pro -i /path/to/fastq_dir/ \ -o /path/to/output/ \ -c config_hicpro.txt 过滤噪音是Hi-C数据处理中必不可少的步骤。自连接（self-ligation）是指同一条限制性片段的两端被连接在一起的情况，这些读对不包含互作信息。未连接（dangling end）是指限制性酶切后未被连接的末端，这些读对也不代表真实的互作。PCR重复需要被去除以避免扩增偏倚。当你处理Hi-C数据时，这些噪音来源的过滤可以显著提高数据质量。 接触矩阵构建是将过滤后的互作读对按照基因组位置分箱（binning），构建表示互作频率的矩阵。bin大小的选择影响数据的分辨率和信噪比：较小的bin（如5kb）提供更高的分辨率但噪音较大，较大的bin（如1Mb）信号更平滑但丢失细节。常用的bin大小包括5kb、10kb、25kb、50kb、100kb、500kb和1Mb，具体选择取决于测序深度和分析目的。当你需要在分辨率和信噪比之间取得平衡时，25kb-50kb的bin大小通常是TAD分析的良好起点。 ICE归一化（迭代纠正经验偏差消除）是Hi-C数据归一化的标准方法。ICE假设基因组上每个bin的互作总数应该相等，任何偏差都来源于技术噪音。ICE通过迭代地除以行和列的边际来消除系统偏差，直到矩阵收敛。归一化后的接触矩阵更适合进行下游的TAD检测和环识别。当你比较不同样本的Hi-C数据时，归一化是确保可比性的关键步骤。 质量控制评估Hi-C数据的多个方面。顺式（cis）互作比率是指同一条染色体内的互作占总互作的比例，高质量的Hi-C数据顺式比率通常较高。反式（trans）互作比率是不同染色体之间的互作比例。有效读对比例是经过过滤后保留的读对占总读对的比例。受限制片段比对效率反映了比对到预期限制性酶切片段的reads比例。当你评估Hi-C数据质量时，这些指标可以帮助你判断数据是否适合进行下游分析。 可视化是理解Hi-C数据的重要手段。Juicebox是一个交互式的Hi-C数据可视化软件，可以加载.hic格式的接触矩阵，以热图形式展示染色质互作，并支持TAD和环的标注。HiGlass是一个基于Web的Hi-C数据可视化平台，支持多数据集的同步浏览和比较。Circos可以将Hi-C互作以弦图的形式展示在染色体圈图上，适合展示染色体间的反式互作。当你需要在全基因组尺度上浏览Hi-C数据时，Juicebox是最直观的工具。 拓扑关联结构域（TAD）识别是Hi-C数据分析的核心内容之一。TAD是染色质上自我互作频率高于与外部区域互作频率的连续区域，大小通常在数百kb到数Mb之间。TAD在细胞类型之间相对保守，被认为是基因调控的基本单位，增强子通常在TAD内部与其靶基因互作，而TAD边界阻止了跨TAD的异常激活。insulation score方法通过计算每个bin的对角线方向上互作频率的局部最小值来识别TAD边界。Directionality Index方法计算每个bin向上游和下游互作的不对称性，TAD边界处方向性指数发生显著变化。HiTC是R语言中的Hi-C数据分析包，提供了TAD检测和可视化功能。TADbit是一个综合的Hi-C分析框架，支持多种TAD检测算法。当你需要识别基因组上的TAD结构时，这些工具提供了互补的检测策略。 隔板（boundaries）检测关注的是TAD边界的精确位置和特征。TAD边界通常富集CTCF结合位点、活跃启动子和tRNA基因等特征。CTCF蛋白通过其方向性的结合位点形成环挤压（loop extrusion）的结构基础，cohesin复合物沿着染色质滑动直到遇到方向正确的CTCF位点停止，从而形成TAD边界。当你研究TAD边界的生物学特征时，将边界位置与CTCF ChIP-seq和基因注释数据进行整合分析可以揭示边界形成的分子机制。 染色质环检测识别基因组上两个远距离位点之间具有显著高互作频率的互作对。FitHiC2使用样条拟合来估计期望的互作频率分布，然后检验观察到的互作频率是否显著高于期望值。HICCUPS是Juicer软件包中的环检测算法，它使用局部背景模型和多重检验校正来识别显著的染色质环。Mustache使用类似于图像处理中的尺度空间方法来检测接触矩阵中的显著互作峰。当你需要识别增强子-启动子环等功能性染色质互作时，这些工具提供了统计严谨的检测方法。 高级互作分析包括A/B区室识别。A/B区室是染色质在更大尺度上的组织形式，A区室对应开放、基因密集、早期复制的常染色质，B区室对应紧密、基因稀疏、晚期复制的异染色质。A/B区室可以通过对归一化后的接触矩阵进行PCA分析来识别，第一主成分的正值通常对应A区室，负值对应B区室。当你需要研究染色质的大尺度组织变化时，A/B区室分析提供了比TAD更宏观的视角。 互作网络分析将染色质互作数据构建为网络模型，其中基因组区域是节点，互作频率是边的权重。网络分析可以识别高度连接的枢纽区域（hub），这些区域可能在染色质组织中发挥核心作用。当你需要从系统层面理解染色质三维结构时，网络分析提供了超越两两互作的全局视角。 .. code-block:: r # 使用HiTC进行TAD检测 library(HiTC) # 构建Hi-C数据对象 # 假设已有归一化的接触矩阵 # xgi和ygi为基因组区间对象 # 计算方向性指数 # di <- directionalityIndex(hic_obj) # TAD检测 # tads <- detectTAD(hic_obj, method = "directionality") # 可视化 # mapC(hic_obj) 6.4.3 环与互作注释 ^^^^^^^^^^^^^^^^^^^ 启动子-增强子循环是染色质互作中最具生物学意义的一类互作。Capture Hi-C通过靶向捕获启动子区域的连接片段，可以高分辨率地检测启动子与远端调控元件的互作。CHiCAGO（Capture Hi-C Analysis of Genomic Organisation）是专门为Capture Hi-C数据设计的统计方法，它考虑了距离衰减效应和技术噪音来识别显著的启动子-增强子互作。当你需要将增强子与其靶基因配对时，Capture Hi-C结合CHiCAGO分析提供了最直接的实验证据。 互作与基因表达关联分析将Hi-C数据与RNA-seq数据整合，研究染色质互作变化与基因表达变化之间的关系。当增强子-启动子互作增强时，基因表达通常上调；当互作减弱或丧失时，基因表达可能下调。这种整合分析需要同时考虑互作频率的变化和表达水平的变化，使用相关性分析或回归模型来评估两者之间的关联。当你拥有匹配的Hi-C和RNA-seq数据时，这种整合分析可以揭示三维基因组重组对基因调控的影响。 超级增强子与基因的互作映射将超级增强子的多个组成元件与其靶基因连接起来。超级增强子通常跨越数十kb的区域，包含多个增强子模块，这些模块可能通过不同的互作与同一个或多个靶基因形成环状结构。HiChIP和PLAC-seq等技术可以高效地检测超级增强子相关的染色质互作，帮助你理解超级增强子如何通过三维基因组结构调控基因表达。 疾病相关的非编码变异与靶基因通过互作连接是Hi-C数据在遗传学中的重要应用。GWAS研究发现了大量位于非编码区的疾病相关变异，但这些变异的靶基因往往难以确定，因为它们可能通过远距离的染色质互作调控远端基因。通过将GWAS信号与Hi-C或Capture Hi-C数据整合，可以将非编码变异映射到其物理互作的启动子区域，从而推断候选靶基因。当你需要解读GWAS发现的非编码信号时，染色质互作数据提供了从变异到基因的功能桥梁。 进化保守的互作网络是染色质三维结构在物种间的保守性研究。比较不同物种的Hi-C数据可以发现，TAD结构在哺乳动物之间具有较高程度的保守性，而具体的环互作则可能存在物种差异。进化保守的互作往往与重要的基因调控功能相关，而物种特异的互作可能反映了谱系特异的调控创新。当你进行跨物种的表观遗传学比较时，染色质三维结构的保守性和差异性分析可以揭示调控元件的进化动态。 6.5 表观转录组学 ----------------- 表观转录组学研究的是RNA分子上的化学修饰，这些修饰不影响RNA的碱基序列，但可以影响RNA的剪接、稳定性、翻译效率和定位等性质。RNA修饰与DNA甲基化和组蛋白修饰不同，它们直接作用于转录产物，在转录后调控层面发挥着重要作用。近年来，随着检测技术的进步，特别是m6A（N6-甲基腺苷）检测方法的发展，表观转录组学已经成为表观遗传学研究中一个蓬勃发展的方向。 RNA甲基化修饰类型多样，其中m6A是最常见和研究最深入的RNA内部修饰，在哺乳动物mRNA中平均每个转录本含有3-5个m6A位点。m5C（5-甲基胞嘧啶）在RNA中也有分布，tRNA和rRNA中含量较高，mRNA中的m5C功能正在被揭示。m1A（N1-甲基腺苷）主要分布在tRNA和rRNA中，mRNA5'UTR区域的m1A可能与翻译起始调控相关。假尿嘧啶（Ψ）是最丰富的RNA修饰之一，主要存在于tRNA和rRNA中，影响RNA的结构和功能。5hmC在RNA中也有报道，但研究相对较少。当你研究RNA修饰时，需要了解不同修饰类型的分布特征和功能差异。 m6A检测技术中最常用的是MeRIP-seq（又称m6A-seq），它利用抗m6A抗体免疫沉淀富集含m6A的RNA片段，然后进行高通量测序。MeRIP-seq可以在转录组范围内检测m6A的分布区域，但分辨率约为100-200nt，无法定位到单个碱基。当你需要进行大规模的m6A分布筛查时，MeRIP-seq是最成熟和常用的方法。 m6A单碱基分辨率方法包括miCLIP、m6A-REF-seq和MAZTER-seq等。miCLIP利用m6A抗体与RNA交联后诱导的突变或截断来定位m6A位点，可以达到单核苷酸分辨率。m6A-REF-seq利用m6A修饰对逆转录酶的影响来推断修饰位点。MAZTER-seq利用MazF核酸酶特异性切割m6A之前的ACA序列，通过分析切割模式来定位m6A。当你需要精确定位m6A位点时，这些单碱基分辨率方法比MeRIP-seq提供了更精确的信息。 m6A数据分析使用专门的工具来识别m6A富集区域。exomePeak是最早开发的m6A峰检测工具，它使用滑动窗口策略在MeRIP-seq数据中识别m6A富集峰。m6Aviewer提供了更全面的可视化和分析功能，支持差异m6A分析。当你处理MeRIP-seq数据时，这些工具可以帮助你从IP和input样本的比较中识别m6A修饰区域。 m6A富集与转录后调控的关系是表观转录组学研究的核心问题。m6A富集在终止密码子附近和3'UTR区域，这种特征性的分布暗示了其在mRNA代谢中的调控作用。m6A可以影响可变剪接，m6A修饰的外显子倾向于被跳过。m6A还影响RNA稳定性，通过招募YTHDF2等阅读器蛋白促进mRNA降解。在翻译效率方面，m6A可以通过招募eIF3等翻译起始因子促进帽非依赖性翻译。当你发现基因的m6A水平发生变化时，需要综合考虑其对剪接、稳定性和翻译的多重影响。 m6A差异修饰分析比较不同条件之间m6A修饰水平的变化。差异m6A分析的方法类似于ChIP-seq的差异分析，通过比较IP信号的富集程度来识别差异修饰区域。常用的工具包括exomePeak2、MACS2（应用于MeRIP-seq数据）和DRME等。当你研究m6A在疾病或发育中的动态变化时，差异修饰分析是关键的步骤。 RNA修饰与表型关联分析将m6A等修饰的变化与生物学表型联系起来。例如，肿瘤中m6A写入酶（如METTL3）和擦除酶（如FTO、ALKBH5）的表达异常与肿瘤进展相关，特定基因m6A修饰的变化可能影响其表达并驱动肿瘤发生。当你需要验证m6A修饰的功能重要性时，将修饰变化与基因表达和表型变化进行关联分析是必要的。 f5C和ac4C等其他修饰技术也有各自专用的分析工具。f5C（5-甲酰胞嘧啶）在RNA中的检测可以使用特定的化学转化方法结合测序。ac4C（N4-乙酰胞嘧啶）由NAT10催化，在mRNA的CDS区域富集，可以促进翻译。这些修饰的分析工具相对m6A而言还不够成熟，但随着检测技术的进步，相应的分析方法也在不断发展。 .. code-block:: r # m6A差异修饰分析示例 # 使用exomePeak2进行MeRIP-seq数据分析 # library(exomePeak2) # 输入文件: IP的BAM文件, Input的BAM文件, 基因注释GTF文件 # result <- exomePeak2( # IP_bam = "IP_sample.bam", # Input_bam = "Input_sample.bam", # GENE_ANNO_GTF = "genes.gtf" # ) # 提取m6A峰 # peaks <- result$peak # head(peaks) 6.6 表观遗传重编程与发育 ------------------------- 表观遗传重编程是指在发育过程中表观遗传标记的大规模擦除和重建，这对于确保正确的基因表达程序至关重要。从受精到胚胎发育，从原始生殖细胞到体细胞重编程，表观遗传修饰经历了动态的变化，这些变化决定了细胞的命运和分化方向。理解表观遗传重编程的机制对于发育生物学和再生医学具有重要意义。 受精后DNA甲基化重编程是一个大规模的表观遗传擦除和重建过程。在受精后，父本基因组经历主动去甲基化，通过TET蛋白介导的氧化反应迅速去除大部分5mC；母本基因组则经历被动去甲基化，在卵裂过程中由于DNMT1的排斥而逐渐稀释甲基化标记。这种不对称的去甲基化过程确保了合子基因组的表观遗传重置，为全能性的建立创造了条件。在随后的植入前发育中，甲基化水平逐渐重建，到植入后阶段恢复到体细胞水平。当你研究早期胚胎发育时，理解这种甲基化重编程的时序性是解读甲基化数据的关键。 原始生殖细胞中甲基化擦除是另一种重要的表观遗传重编程事件。在原始生殖细胞（PGC）的发育过程中，基因组经历近乎完全的去甲基化，以擦除亲本的表观遗传记忆，为下一代重新建立印记和性别特异的甲基化模式做准备。PGC的去甲基化涉及被动和主动两种机制，其中TET蛋白介导的羟甲基化在主动去甲基化中发挥核心作用。印记控制区域（ICR）的甲基化在PGC中也被擦除，随后在配子发生过程中根据性别重新建立。当你研究印记基因或生殖细胞发育时，PGC甲基化重编程的知识是理解这些过程的基础。 体细胞重编程（iPSC）中甲基化和组蛋白重塑特征反映了从分化状态到多能性状态的转变。在诱导多能干细胞（iPSC）的重编程过程中，分化相关的基因启动子经历去甲基化，而多能性基因的启动子从沉默状态被激活。组蛋白修饰也发生大规模重塑，分化标记（如H3K9me3）被擦除，多能性标记（如H3K4me3）被建立。然而，iPSC的重编程并不完全，某些区域可能保留供体细胞的表观遗传记忆，这被称为表观遗传残留。当你使用iPSC模型时，需要意识到这种表观遗传残留可能影响iPSC的分化倾向。 DMR与印记控制区（ICR）的关系是表观遗传重编程研究中的核心问题。印记控制区是调控印记基因表达的DNA区域，其甲基化状态由亲本来源决定。ICR的甲基化在配子发生过程中建立，在受精后和体细胞分裂中被维持。当ICR的甲基化发生异常时，会导致印记基因的表达紊乱，引起如Prader-Willi综合征、Angelman综合征和Beckwith-Wiedemann综合征等印记疾病。当你研究印记基因的调控时，ICR的甲基化状态是最重要的分析对象。 等位基因特异性甲基化（印记基因）是基因组印记的分子基础。印记基因的两个等位基因具有不同的甲基化状态，只有来自特定亲本的等位基因被表达。这种等位基因特异性的表达模式依赖于ICR的差异性甲基化，ICR通过调控染色质结构和CTCF结合来控制印记基因簇的表达。当你分析印记基因时，需要使用能够区分等位基因的分析方法，如等位基因特异性甲基化分析。 X染色体失活是雌性哺乳动物中一条X染色体被沉默的过程，以确保X连锁基因的剂量补偿。Xist长非编码RNA是X染色体失活的启动因子，它从失活的X染色体上转录并顺式包裹整条染色体，招募PRC2等抑制性复合物建立H3K27me3修饰和DNA甲基化。X染色体失活是一个多层表观遗传调控的过程，涉及组蛋白修饰、DNA甲基化和染色质重塑的协同作用。当你研究雌性样本的表观遗传数据时，需要注意X染色体失活对甲基化和组蛋白修饰信号的影响。 6.7 表观基因组关联研究（EWAS） ------------------------------- 表观基因组关联研究（EWAS）是将表观遗传变异与疾病或性状关联起来的研究方法。与GWAS关注遗传变异不同，EWAS关注的是表观遗传变异（主要是DNA甲基化）与表型的关联。EWAS的优势在于表观遗传变异比遗传变异更动态，可以反映环境暴露和生活方式的影响，因此能够捕获GWAS无法检测的疾病风险因素。然而，EWAS也面临细胞异质性、因果方向不确定和多重检验负担等挑战。 EWAS设计通常采用病例对照的研究框架，使用甲基化芯片（如EPIC 850K）或WGBS进行全基因组甲基化检测。样本量的确定需要考虑效应大小、甲基化变异和期望的统计功效，通常需要数百到数千个样本才能获得可靠的关联结果。当你设计EWAS研究时，需要特别注意匹配病例和对照的年龄、性别、组织类型和细胞组成，因为这些因素都是甲基化的重要混杂因素。 细胞异质性校正是EWAS分析中最关键的步骤之一。不同细胞类型具有截然不同的甲基化模式，如果病例和对照之间的细胞组成存在差异，这种差异可能被误认为与疾病相关的甲基化变化。Houseman方法是最常用的细胞组成估计方法，它利用参考面板中各细胞类型的特征性甲基化模式来估计混合样本中各细胞类型的比例。估计的细胞比例可以作为协变量纳入关联分析模型中进行校正。当你处理血液样本的EWAS数据时，白细胞比例的校正是必须的步骤。 差异甲基化位点（DMP）关联分析是EWAS的核心分析内容。对于每个CpG位点，使用线性回归或混合模型检验甲基化水平与疾病状态的关联，同时校正年龄、性别、细胞比例等协变量。对于连续型性状（如血压、BMI等），使用甲基化水平与表型值的线性关联模型。对于二分类性状（如疾病/健康），使用逻辑回归模型。当你进行大规模的EWAS分析时，需要使用高效的计算工具来处理数百万个CpG位点的检验。 表观变异解释遗传力是指DNA甲基化变异可以作为遗传变异和表型之间的中介，帮助解释GWAS无法完全解释的遗传力。mQTL分析表明，相当比例的甲基化变异受遗传变异调控，这些受遗传调控的甲基化位点可能通过影响基因表达来贡献表型变异。当你需要理解遗传变异如何通过表观遗传机制影响疾病时，将mQTL与EWAS结果整合可以提供更完整的因果链。 EWAS多重检验校正是控制假阳性的必要步骤。由于EWAS同时检验数百万个CpG位点，需要使用FWER（家族错误率）或FDR（错误发现率）方法进行多重检验校正。Bonferroni校正最为严格，但可能过于保守；FDR校正（如Benjamini-Hochberg方法）在控制假阳性的同时保留了更多的统计功效。当你报告EWAS结果时，应该同时报告校正前和校正后的p值。 复制与meta分析是验证EWAS结果可靠性的重要步骤。由于EWAS容易受到批次效应和混杂因素的影响，独立队列的复制验证是确认关联真实性的关键。Meta分析可以整合多个EWAS研究的结果，提高统计功效并减少单个研究的偏倚。当你发现EWAS关联信号时，在独立队列中进行复制验证是必要的后续步骤。 EWAS结果与遗传多态性（mQTL）整合可以帮助理解甲基化变化的遗传基础。如果一个EWAS关联的CpG位点同时也是mQTL位点，那么遗传变异可能通过影响甲基化进而影响疾病风险。这种整合分析可以使用共定位（colocalization）方法来评估EWAS信号和mQTL信号是否共享同一个因果变异。当你需要区分甲基化变化是疾病的原因还是结果时，mQTL整合分析提供了遗传学层面的证据。 孟德尔随机化用于甲基化-疾病因果推断是一种利用遗传变异作为工具变量来推断甲基化与疾病之间因果关系的统计方法。如果遗传变异与甲基化水平相关（mQTL），且遗传变异仅通过甲基化影响疾病风险（排除多效性），那么可以推断甲基化与疾病之间存在因果关系。当你需要确定EWAS发现的甲基化变化是否具有因果作用时，孟德尔随机化分析提供了比观察性关联更强的因果证据。 6.8 表观遗传时钟与衰老 ----------------------- 表观遗传时钟是基于DNA甲基化水平预测生物学年龄的统计模型。随着年龄增长，基因组中特定CpG位点的甲基化水平发生规律性变化，这些变化可以被用来构建年龄预测模型。表观遗传时钟不仅是衰老研究的重要工具，还为研究环境暴露、疾病状态和生活方式对生物学衰老的影响提供了量化指标。 基于CpG甲基化水平的年龄预测模型中，Horvath时钟和Hannum时钟是最著名的两个模型。Horvath时钟使用353个CpG位点，基于多种组织和细胞类型的甲基化数据训练，具有泛组织适用性，可以准确预测大多数组织类型的年龄。Hannum时钟使用71个CpG位点，主要基于血液样本训练，在血液样本中预测精度较高。这两个模型都使用弹性网络回归（elastic net regression）进行特征选择和模型训练。当你需要评估样本的表观遗传年龄时，Horvath时钟是最通用的选择。 泛组织时钟与组织特异性时钟各有其适用场景。泛组织时钟（如Horvath时钟）的优势在于可以在不同组织类型之间进行比较，但其在特定组织中的精度可能不如组织特异性时钟。组织特异性时钟针对特定组织（如血液、大脑、皮肤等）训练，在该组织中具有更高的预测精度。PhenoAge时钟和GrimAge时钟是第二代表观遗传时钟，它们不仅预测日历年龄，还整合了与衰老相关的临床指标，对死亡率和发病率的预测能力更强。当你研究特定组织的衰老时，组织特异性时钟可能提供更精确的估计。 表观遗传年龄加速（AgeAccel）是指表观遗传年龄与日历年龄之间的差值。正的AgeAccel表示个体的表观遗传年龄大于日历年龄，暗示加速衰老；负的AgeAccel表示表观遗传年龄小于日历年龄，暗示延缓衰老。AgeAccel的计算方式有多种，包括残差法（表观遗传年龄对日历年龄回归的残差）和差值法（表观遗传年龄减去日历年龄）。当你研究影响衰老的因素时，AgeAccel是核心的因变量。 衰老相关DMR（年龄相关甲基化漂移）描述了随年龄变化的甲基化模式。某些CpG位点随年龄发生渐进性高甲基化，而另一些位点则发生渐进性低甲基化，这种现象被称为"甲基化漂移"。年龄相关的高甲基化主要发生在CpG岛，而低甲基化主要发生在CpG岛外的区域。这种漂移可能与细胞分裂过程中DNMT1保真度的下降和被动去甲基化有关。当你分析衰老相关的甲基化数据时，区分年龄相关的高甲基化和低甲基化区域有助于理解衰老的表观遗传机制。 时钟模型建立使用岭回归或弹性网络等正则化回归方法。弹性网络回归结合了L1（Lasso）和L2（Ridge）正则化，可以在高维数据（数十万个CpG位点）中进行特征选择和模型拟合。模型训练需要大样本量的甲基化数据，使用交叉验证来评估模型的预测性能。当你需要构建特定人群或特定组织的表观遗传时钟时，需要收集足够大的训练数据集并进行严格的交叉验证。 表观遗传时钟在疾病和环境暴露中的应用日益广泛。在疾病研究中，多种疾病（如癌症、心血管疾病、神经退行性疾病等）与AgeAccel相关，疾病状态可能导致表观遗传年龄加速。在环境暴露研究中，吸烟、空气污染、早期生活压力等因素与AgeAccel正相关，而健康生活方式（如运动、饮食）可能与AgeAccel负相关。当你需要量化环境暴露或疾病对生物学衰老的影响时，表观遗传时钟提供了一个客观的度量工具。 .. code-block:: r # 使用R计算表观遗传年龄 # 使用DNAmAge包或手动实现Horvath时钟 # library(DNAmAge) # 假设已有甲基化beta值矩阵 # beta_matrix: 行为CpG位点，列为样本 # 计算Horvath时钟年龄 # horvath_age <- agep(beta_matrix, method = "horvath") # 计算AgeAccel # chronological_age <- c(45, 52, 38, 61, 55) # age_accel <- horvath_age - chronological_age # cat("表观遗传年龄加速:\n") # print(age_accel) 6.9 染色质动力学与蛋白-DNA相互作用 ------------------------------------ 染色质动力学研究的是染色质结合蛋白与DNA的动态相互作用，这些相互作用决定了染色质的状态和基因的调控模式。与传统的ChIP-seq相比，新一代的染色质分析技术（如CUT&RUN和CUT&Tag）提供了更高的信噪比和更低的细胞需求，使得在有限样本中研究蛋白-DNA互作成为可能。 CUT&RUN（靶向切割与释放）利用蛋白A/G融合的微球菌核酸酶（pA-MNase）在抗体靶向的位置切割DNA，释放与目标蛋白结合的DNA片段。与ChIP-seq相比，CUT&RUN在完整细胞中进行，避免了交联和超声片段化带来的偏倚，信噪比显著提高，且需要的细胞数更少（可低至100个细胞）。当你需要在有限样本中检测转录因子或组蛋白修饰的结合位点时，CUT&RUN是比ChIP-seq更高效的选择。 CUT&Tag（转座子融合蛋白）利用蛋白A融合的Tn5转座酶（pA-Tn5）在抗体靶向的位置切割并标记DNA，同时加入测序接头。CUT&Tag的实验流程比CUT&RUN更简单，可以直接在细胞或细胞核中进行，且产生的文库可以直接用于测序。CUT&Tag特别适合于低输入量和单细胞应用，当你需要进行高通量的蛋白-DNA互作分析时，CUT&Tag提供了最简便的方案。 DamID（DNA腺嘌呤甲基转移酶鉴定）利用大肠杆菌DNA腺嘌呤甲基转移酶（Dam）与目标蛋白融合表达，在目标蛋白结合的DNA区域引入GATC位点的腺嘌呤甲基化（m6A），然后通过DpnI酶切和测序检测这些甲基化位点。DamID不需要抗体和交联，可以在活细胞中实时检测蛋白-DNA互作。当你需要研究蛋白-DNA互作的动态变化或避免抗体相关的偏倚时，DamID提供了独特的优势。 染色质蛋白结合位点的检测可以通过目标蛋白的ChIP-seq来实现，这些数据与组蛋白修饰的ChIP-seq分析类似，但针对的是特定的染色质结合蛋白（如CTCF、cohesin、RNA聚合酶II等）。转录因子的ChIP-seq产生窄峰，分析流程与组蛋白修饰的窄峰分析类似。当你研究特定蛋白在基因组上的结合模式时，ChIP-seq仍然是最经典的方法。 共占位分析通过比较两种蛋白的Peaks重叠来研究它们在基因组上的协同结合。如果两种蛋白的Peaks显著重叠，说明它们可能在相同的基因组区域协同调控。共占位分析需要使用统计检验（如超几何检验）来评估重叠的显著性，因为随机重叠也可能发生。当你需要研究两个蛋白是否在相同的调控元件上协同作用时，共占位分析提供了基因组层面的证据。 .. code-block:: bash # CUT&Tag数据处理流程 # 1. 比对（与ATAC-seq类似，使用Bowtie2） bowtie2 -x /path/to/genome_index \ -1 sample_R1.fq.gz -2 sample_R2.fq.gz \ --very-sensitive -X 2000 \ -p 8 -S sample.sam # 2. 转换并排序 samtools view -bS sample.sam | samtools sort -o sample.sorted.bam samtools index sample.sorted.bam # 3. 去除线粒体reads和重复 samtools idxstats sample.sorted.bam | cut -f 1 | grep -v "chrM" | \ xargs samtools view -b sample.sorted.bam > sample.noMT.bam picard MarkDuplicates INPUT=sample.noMT.bam OUTPUT=sample.dedup.bam \ METRICS_FILE=sample.metrics.txt REMOVE_DUPLICATES=true # 4. 峰值调用（使用MACS2或SEACR） # SEACR特别适合CUT&RUN/CUT&Tag数据 # bash SEACR_1.3.sh sample.bedgraph 0.01 norm stringent sample_peaks 6.10 表观编辑与干扰 -------------------- 表观编辑是指利用基因工程工具定向修改基因组特定位点的表观遗传状态，而不改变DNA序列本身。CRISPR-dCas9系统的出现使得靶向表观编辑成为可能，研究者可以精确地在特定基因位点引入或去除甲基化标记、组蛋白修饰等表观遗传变化，从而直接验证表观遗传修饰的功能。表观编辑技术不仅在基础研究中用于功能验证，在疾病治疗中也展现了潜在的应用前景。 dCas9融合表观修饰酶是表观编辑的基本策略。将催化失活的Cas9（dCas9）与表观修饰酶融合，通过sgRNA引导到目标位点，可以在该位点引入特定的表观遗传修饰。dCas9-DNMT3A融合蛋白可以在目标位点引入DNA甲基化，沉默基因表达。dCas9-TET1融合蛋白可以在目标位点去除DNA甲基化，激活基因表达。dCas9-p300融合蛋白具有组蛋白乙酰转移酶活性，可以在目标位点引入H3K27ac修饰，激活增强子和启动子。当你需要验证特定CpG位点的甲基化对基因表达的影响时，dCas9-DNMT3A或dCas9-TET1提供了直接的实验工具。 CRISPRoff/CRISPRon是专门用于DNA甲基化编辑的表观编辑系统。CRISPRoff将dCas9与DNMT3A和DNMT3L融合，可以在目标基因的启动子区域引入稳定的DNA甲基化，实现长期的基因沉默。这种沉默效应即使在移除编辑工具后仍然可以维持，并且可以通过细胞分裂传递给子代细胞。CRISPRon将dCas9与TET1融合，可以去除CRISPRoff引入的甲基化，恢复基因表达。当你需要实现可遗传的基因沉默时，CRISPRoff提供了比传统的CRISPRi更持久的解决方案。 组蛋白修饰编辑使用dCas9与组蛋白修饰酶的融合蛋白来改变目标位点的组蛋白修饰状态。dCas9-HDAC融合蛋白可以去除组蛋白乙酰化标记，抑制基因表达。dCas9-p300融合蛋白可以添加H3K27ac标记，激活基因表达。dCas9-LSD1融合蛋白可以去除H3K4me2标记，抑制增强子活性。当你需要研究特定组蛋白修饰对基因调控的影响时，组蛋白修饰编辑提供了直接的因果关系证据。 瞬时表观调控技术追求的是可逆和可控的表观遗传修饰改变。与基因组编辑不同，表观编辑的目标是实现可逆的修饰变化，使研究者能够观察修饰引入和去除对基因表达的影响。一些瞬时调控技术使用小分子诱导系统（如脱落酸诱导的dCas9-表观修饰酶系统），可以在特定时间点启动或终止表观编辑。当你需要研究表观遗传修饰的时间依赖性效应时，瞬时调控技术比组成型表达系统更合适。 编辑后评估是验证表观编辑效果的关键步骤。靶向甲基化测序（如靶向BS-seq或靶向RRBS）可以精确测量编辑位点和周围区域的甲基化水平变化。ChIP-seq可以评估组蛋白修饰编辑的效果。RNA-seq可以评估编辑对基因表达的影响。当你完成表观编辑实验后，需要使用这些分析方法来确认编辑的特异性和效率，并评估是否存在脱靶效应。 .. code-block:: r # 表观编辑效果评估示例 # 比较编辑前后目标位点的甲基化水平 before_edit <- c(0.1, 0.15, 0.12, 0.08, 0.11) after_edit <- c(0.75, 0.82, 0.78, 0.71, 0.80) # t检验评估甲基化变化 t_result <- t.test(before_edit, after_edit, paired = TRUE) cat(sprintf("甲基化变化: %.3f -> %.3f\n", mean(before_edit), mean(after_edit))) cat(sprintf("p值: %.4e\n", t_result$p.value)) # 评估基因表达变化 expr_before <- c(5.2, 4.8, 5.5, 4.9, 5.1) expr_after <- c(1.2, 0.8, 1.5, 0.9, 1.1) expr_result <- t.test(expr_before, expr_after, paired = TRUE) cat(sprintf("基因表达变化: %.2f -> %.2f\n", mean(expr_before), mean(expr_after))) cat(sprintf("p值: %.4e\n", expr_result$p.value)) 6.11 表观遗传数据分析的最佳实践 --------------------------------- 表观遗传数据分析涉及多种数据类型和复杂的分析流程，遵循最佳实践可以确保分析结果的可靠性和可重复性。从实验设计到数据处理、从统计分析到结果解读，每个环节都有需要注意的关键问题。本节总结表观遗传数据分析中的最佳实践建议，帮助研究者避免常见的陷阱并提高分析质量。 处理不同批次的样本时加入批内对照是控制批次效应的基本策略。表观遗传数据对技术变异特别敏感，不同批次之间的系统差异可能被误认为生物学差异。在每个批次中加入共同的对照样本（如参考DNA或细胞系），可以用来评估和校正批次效应。当你必须分批处理样本时，确保每个批次包含各处理组的样本，并加入批内对照以便后续的批次校正。 表观数据与基因表达数据联合分析是多组学整合的重要方向。将甲基化或染色质可及性数据与RNA-seq数据整合，可以揭示表观遗传修饰对基因表达的调控关系。常用的整合方法包括MOFA（Multi-Omics Factor Analysis）、Seurat的多模态整合、LIGER和mixOmics等。这些方法可以在不同组学数据之间识别共享的变异模式，揭示跨层次的调控关系。当你拥有匹配的表观遗传和转录组数据时，多组学整合分析可以提供比单一组学更全面的生物学洞见。 高通量表观数据的可重复性是确保结果可信的关键。ENCODE项目制定了严格的数据质量标准和分析流程，为表观遗传数据的可重复性提供了参考。使用标准化的分析流程（如ENCODE管道）可以减少分析变异，提高结果的可比性。当你报告表观遗传分析结果时，应该详细描述分析流程和参数设置，以便其他研究者能够重复你的分析。 使用Anndata存储甲基化和染色质可及性矩阵是单细胞表观数据管理的推荐做法。Anndata是Python中Scanpy生态系统的数据结构，可以高效地存储高维稀疏矩阵和丰富的注释信息。对于单细胞ATAC-seq和单细胞甲基化数据，Anndata格式提供了灵活的数据存储和访问方式。当你处理大规模的单细胞表观数据时，使用Anndata可以简化数据管理并提高分析效率。 工作流管理使用标准化的分析流程可以确保分析的可重复性和可移植性。nf-core项目提供了一系列基于Nextflow的高质量分析流程，包括nf-core/methylseq（甲基化测序分析）、nf-core/atacseq（ATAC-seq分析）和nf-core/chipseq（ChIP-seq分析）。这些流程经过社区审查和测试，集成了最佳实践的分析步骤，支持容器化部署，可以在不同的计算环境中运行。当你需要处理大规模的表观遗传数据时，使用nf-core流程可以节省流程开发的时间并确保分析质量。 .. code-block:: bash # 使用nf-core流程进行表观遗传数据分析 # nf-core/methylseq: 甲基化测序分析 nextflow run nf-core/methylseq \ -profile docker \ --input samplesheet.csv \ --genome GRCh38 \ --aligner bismark \ --outdir results/methylseq # nf-core/atacseq: ATAC-seq分析 nextflow run nf-core/atacseq \ -profile docker \ --input samplesheet.csv \ --genome GRCh38 \ --outdir results/atacseq # nf-core/chipseq: ChIP-seq分析 nextflow run nf-core/chipseq \ -profile docker \ --input samplesheet.csv \ --genome GRCh38 \ --outdir results/chipseq