m6A修饰分析MeRIP¶

一句话概述：m6A（N6-甲基腺苷）是mRNA上最常见的化学修饰，MeRIP-seq通过抗体免疫沉淀富集m6A片段后测序来定位修饰位点，exomePeak2/MACS2是常用的peak calling工具，2024年m6ACali用机器学习校准减少假阳性。

核心知识点速览¶

概念	白话解释
m6A	RNA上腺苷（A）的第6位氮原子加了个甲基（-CH3）
MeRIP-seq	用m6A抗体"钓"出带修饰的RNA片段，再测序
表观转录组	类比表观基因组，研究RNA上的化学修饰
Writer	写入m6A的酶（METTL3/METTL14/WTAP复合物）
Eraser	去除m6A的酶（FTO、ALKBH5）
Reader	识别m6A的蛋白（YTHDF1/2/3、YTHDC1/2、IGF2BP1/2/3）
DRACH基序	m6A最常出现的序列上下文：D=A/G/U, R=A/G, A=修饰位点, C, H=A/C/U
exomePeak2	专门为MeRIP-seq设计的R/Bioconductor peak calling包
MACS2	原为ChIP-seq设计的peak calling工具，可适配MeRIP-seq
m6ACali	2024年机器学习校准框架，减少抗体非特异性结合的假阳性

一、m6A修饰基础¶

1.1 m6A的生物学意义¶

m6A是mRNA上最丰富的内部修饰（占所有甲基化修饰的~80%）

m6A的"写入-擦除-读取"系统：
  Writers（甲基转移酶）: METTL3 + METTL14 + WTAP → 催化m6A安装
  Erasers（去甲基化酶）: FTO、ALKBH5 → 去除m6A
  Readers（识别蛋白）:   YTHDF1/2/3 → 影响RNA命运

m6A的功能：
  1. 影响mRNA稳定性（YTHDF2促进降解）
  2. 调控翻译效率（YTHDF1促进翻译）
  3. 影响可变剪接（YTHDC1）
  4. 调控RNA核质转运
  5. 影响miRNA加工
  6. 参与DNA损伤修复

临床意义：
  - m6A异常与多种癌症相关
  - FTO是肥胖和2型糖尿病的风险基因
  - m6A调节免疫应答和肿瘤微环境

1.2 MeRIP-seq原理¶

MeRIP-seq（Methylated RNA Immunoprecipitation Sequencing）流程：

1. 提取总RNA → 打断成100-200nt片段
2. 取一部分作为Input（背景对照）
3. 剩余部分用m6A抗体做免疫沉淀（IP）
4. 洗涤去除非特异结合
5. 提取IP和Input的RNA
6. 建库测序（通常是RNA-seq建库流程）

分析思路：
  IP样本中富集的区域 vs Input背景 → m6A修饰位点（peak）
  类似ChIP-seq的peak calling思路

局限性：
  - 分辨率100-200nt（一个peak内可能有多个m6A位点）
  - 抗体有非特异性结合（假阳性）
  - 需要大量起始RNA（通常>5μg总RNA）

二、MeRIP-seq数据分析流程¶

2.1 数据预处理¶

# 1. 质控
fastqc IP_R1.fq.gz IP_R2.fq.gz Input_R1.fq.gz Input_R2.fq.gz  # 质控报告

# 2. 接头去除
trim_galore --paired \
  --fastqc \                  # 去接头后再质控
  IP_R1.fq.gz IP_R2.fq.gz    # IP样本
trim_galore --paired \
  --fastqc \
  Input_R1.fq.gz Input_R2.fq.gz  # Input样本

# 3. 比对到基因组（用STAR或HISAT2）
STAR --runThreadN 8 \
  --genomeDir star_index/ \          # STAR索引
  --readFilesIn IP_R1_trimmed.fq.gz IP_R2_trimmed.fq.gz \  # IP reads
  --readFilesCommand zcat \          # gzip文件
  --outSAMtype BAM SortedByCoordinate \  # 输出排序BAM
  --outFileNamePrefix IP_             # 前缀

STAR --runThreadN 8 \
  --genomeDir star_index/ \
  --readFilesIn Input_R1_trimmed.fq.gz Input_R2_trimmed.fq.gz \
  --readFilesCommand zcat \
  --outSAMtype BAM SortedByCoordinate \
  --outFileNamePrefix Input_

# 4. 索引BAM
samtools index IP_Aligned.sortedByCoord.out.bam     # IP索引
samtools index Input_Aligned.sortedByCoord.out.bam   # Input索引

# 5. 去重复（推荐）
picard MarkDuplicates \
  INPUT=IP_Aligned.sortedByCoord.out.bam \
  OUTPUT=IP_dedup.bam \
  METRICS_FILE=IP_dup_metrics.txt \
  REMOVE_DUPLICATES=true        # 去除PCR重复
samtools index IP_dedup.bam     # 重新索引

2.2 MACS2 Peak Calling¶

# MACS2原本为ChIP-seq设计，需要调整参数用于MeRIP-seq

macs2 callpeak \
  -t IP_dedup.bam \            # IP样本（treatment）
  -c Input_dedup.bam \         # Input样本（control）
  -f BAMPE \                   # 双端测序格式
  --nomodel \                  # 不建立模型（RNA不同于DNA）
  --extsize 150 \              # 片段延伸长度（~RNA片段大小）
  --gsize 1e8 \                # 有效基因组大小（转录组~100Mb）
  -q 0.05 \                    # FDR阈值
  --keep-dup all \             # 保留所有reads
  -n m6A_peaks \               # 输出前缀
  --outdir macs2_output/       # 输出目录

# MACS2输出文件：
# m6A_peaks_peaks.narrowPeak — peak位置和统计信息
# m6A_peaks_summits.bed — peak顶点位置
# m6A_peaks_treat_pileup.bdg — IP信号
# m6A_peaks_control_lambda.bdg — Input信号

# 重要：MeRIP-seq用MACS2的注意事项
# 1. 必须加--nomodel（RNA片段分布不同于ChIP-seq）
# 2. --gsize用转录组大小（~1e8），不用基因组大小
# 3. --extsize根据实际RNA片段大小调整

2.3 exomePeak2 Peak Calling（推荐）¶

# exomePeak2是专为MeRIP-seq设计的R包
# BiocManager::install("exomePeak2")

library(exomePeak2)

# 运行peak calling
result <- exomePeak2(
  bam_ip = c("IP_rep1.bam", "IP_rep2.bam"),        # IP样本BAM文件
  bam_input = c("Input_rep1.bam", "Input_rep2.bam"), # Input样本BAM
  txdb = TxDb.Hsapiens.UCSC.hg38.knownGene,         # 转录组数据库
  genome = "hg38",                                    # 基因组版本
  paired_end = TRUE                                   # 双端测序
)

# 获取peak结果
peaks <- results(result)  # 提取peak GRanges对象
length(peaks)             # peak总数

# 导出为BED文件
export(peaks, "m6A_peaks.bed", format = "BED")  # 导出

# exomePeak2的优势：
# 1. 专门为RNA修饰设计（考虑转录组结构）
# 2. 支持差异修饰分析（两组样本比较）
# 3. GC含量校正
# 4. 基于负二项分布的统计模型

2.4 差异m6A分析¶

# 比较两个条件下m6A修饰的差异
library(exomePeak2)

# 差异peak calling
diff_result <- exomePeak2(
  bam_ip = c("KO_IP_1.bam", "KO_IP_2.bam"),          # 处理组IP
  bam_input = c("KO_Input_1.bam", "KO_Input_2.bam"),  # 处理组Input
  bam_ip_treated = c("WT_IP_1.bam", "WT_IP_2.bam"),   # 对照组IP
  bam_input_treated = c("WT_Input_1.bam", "WT_Input_2.bam"),  # 对照组Input
  txdb = TxDb.Hsapiens.UCSC.hg38.knownGene,
  genome = "hg38",
  paired_end = TRUE
)

# 提取差异修饰位点
diff_peaks <- results(diff_result)
# log2FC > 0: 处理组m6A增加（hyper-methylation）
# log2FC < 0: 处理组m6A减少（hypo-methylation）

三、m6ACali机器学习校准（2024）¶

# m6ACali解决MeRIP-seq的抗体非特异性结合问题
# 不需要IVT对照，从数据特征直接判断真假peak

# m6ACali的核心思路：
# 1. 训练集：用已验证的m6A位点（如miCLIP/m6A-Atlas）
# 2. 特征：提取每个peak的序列特征（外显子长度、基序、GC含量等）
# 3. 分类：用机器学习模型区分真实m6A peak和假阳性
# 4. 校准：对exomePeak2/MACS2的结果进行二次过滤

# m6ACali验证了326个人类MeRIP-seq实验
# 校准后peak中DRACH基序的富集度显著提高
# 发现假阳性常出现在：短外显子、短mRNA、高覆盖区域

四、m6A位点注释和可视化¶

# m6A peak的基因组区域注释
library(Guitar)  # MeRIP-seq专用可视化包

# Guitar图：显示m6A在mRNA上的分布
# m6A主要富集在3'UTR和终止密码子附近

# 读取peak数据
peaks_gr <- import("m6A_peaks.bed")  # 导入peak

# 绘制Guitar图
GuitarPlot(
  peaks_gr,
  txdb = TxDb.Hsapiens.UCSC.hg38.knownGene,
  genome = "hg38",
  saveToPDFprefix = "m6A_distribution"  # 保存PDF
)
# Guitar图显示的是m6A peak在mRNA结构上的分布
# 典型模式：3'UTR和CDS末端富集

# DRACH基序分析
library(Biostrings)
library(BSgenome.Hsapiens.UCSC.hg38)

# 提取peak序列
peak_seqs <- getSeq(BSgenome.Hsapiens.UCSC.hg38, peaks_gr)  # 获取序列

# 搜索DRACH基序（D=A/G/T, R=A/G, A, C, H=A/C/T）
drach_count <- vcountPattern("DRACH", peak_seqs,
                              fixed = FALSE)  # 模糊匹配
drach_pct <- sum(drach_count > 0) / length(drach_count) * 100
cat("含DRACH基序的peak占比:", round(drach_pct, 1), "%\n")
# 正常应该>60%，如果太低说明假阳性多

# 功能富集分析（m6A修饰基因的GO/KEGG）
library(clusterProfiler)
library(org.Hs.eg.db)

# 提取m6A修饰基因
m6a_genes <- unique(peaks_gr$gene_name)  # 获取基因列表

# GO富集
ego <- enrichGO(
  gene = m6a_genes,
  OrgDb = org.Hs.eg.db,
  keyType = "SYMBOL",
  ont = "BP",                  # 生物过程
  pvalueCutoff = 0.05
)
dotplot(ego, showCategory = 15)  # 气泡图

五、新一代m6A检测方法¶

传统MeRIP-seq的局限性和新方法：

| 方法 | 分辨率 | 需要抗体 | 输入量 |
|------|--------|----------|--------|
| MeRIP-seq | 100-200nt | 是 | >5μg |
| miCLIP/eCLIP | 单碱基 | 是 | >10μg |
| MAZTER-seq | 单碱基 | 否 | ~100ng |
| DART-seq | ~50nt | 否 | <1μg |
| m6A-SAC-seq | 单碱基 | 否 | ~50ng |
| Nanopore m6Anet | 单碱基 | 否 | 直接RNA |

Nanopore直接RNA测序：
  - m6Anet：深度学习模型从Nanopore信号直接检测m6A
  - 不需要抗体，不需要片段化
  - 可以检测单分子级别的修饰
  - 但通量和准确率仍在提升中

常见报错与解决¶

报错信息	原因	解决方案
`MACS2: too few peaks`	IP富集效率低	检查IP/Input比值，降低q值阈值
`exomePeak2: txdb mismatch`	注释版本和比对基因组不一致	确保TxDb和BAM用同一基因组版本
`No DRACH motif enrichment`	假阳性太多	用m6ACali校准或增加严格过滤
`Too many peaks (>50000)`	FDR阈值太宽松	提高stringency，用q<0.01
`Guitar: no peak in mRNA`	peak没有注释到转录组	检查peak坐标是否为基因组坐标
`Differential: no sig peaks`	样本差异不大或重复数少	增加生物学重复数（≥3）

速查表¶

# MeRIP-seq分析流程
质控(FastQC) → 去接头(trim_galore) → 比对(STAR/HISAT2)
→ 去重(Picard) → Peak Calling(exomePeak2/MACS2)
→ 差异分析(exomePeak2) → 注释(Guitar) → 功能富集

# Peak Calling工具选择
推荐首选: exomePeak2（专为MeRIP-seq设计）
替代方案: MACS2（需要调参：--nomodel --gsize 1e8）
质量校准: m6ACali（减少抗体非特异性假阳性）
整合流水线: MeRIPseqPipe（Nextflow自动化）

# MACS2关键参数（MeRIP-seq适配）
--nomodel: 必须（RNA片段分布不同于DNA）
--gsize 1e8: 转录组有效大小
--extsize 150: RNA片段大小
-f BAMPE: 双端测序

# m6A质量检查
DRACH基序富集: >60%的peak含DRACH
分布模式: 3'UTR和终止密码子附近富集（Guitar图）
peak数量: 通常10,000-30,000个
IP/Input比值: 通常2-5倍富集

# m6A调控因子
Writers: METTL3/METTL14/WTAP（安装m6A）
Erasers: FTO/ALKBH5（去除m6A）
Readers: YTHDF1/2/3, YTHDC1/2, IGF2BP1/2/3（识别m6A）

面试高频问题¶

Q1：m6A修饰的生物学功能是什么？¶

答：m6A是mRNA上最丰富的内部修饰，通过"Writer-Eraser-Reader"三元系统动态调控。Writer（METTL3/14/WTAP复合物）安装m6A，Eraser（FTO/ALKBH5）去除m6A，Reader蛋白识别m6A后介导下游效应：YTHDF2促进mRNA降解，YTHDF1促进翻译，YTHDC1影响剪接和核输出，IGF2BP1/2/3稳定mRNA。m6A调控mRNA的整个生命周期——从剪接、核输出、翻译到降解。

Q2：MeRIP-seq和ChIP-seq在分析上有什么区别？¶

答：两者都用抗体富集+测序的策略，但关键区别：①MeRIP-seq的目标是RNA而非DNA，需要比对到转录组而非基因组；②MeRIP-seq的分辨率更低（100-200nt vs ChIP-seq的~200bp），因为RNA需要先打断；③MeRIP-seq用MACS2时需要调参——--nomodel（RNA片段分布不同）、--gsize 1e8（转录组大小而非基因组大小）；④推荐使用exomePeak2等专门为RNA修饰设计的工具。

Q3：如何判断MeRIP-seq数据质量？¶

答：四个维度：①IP效率——IP样本reads应该明显富集在m6A位点，IP/Input比值通常2-5倍；②DRACH基序——真实m6A peak应该>60%含有DRACH共识序列（D=A/G/U, R=A/G, A, C, H=A/C/U）；③分布模式——Guitar图应显示m6A在3'UTR和终止密码子附近明显富集；④重复一致性——至少2个生物学重复的peak重叠率应>50-60%（已知MeRIP-seq重复间overlap仅30-60%，这是该技术的固有局限）。

Q4：exomePeak2和MACS2用于MeRIP-seq有什么区别？¶

答：exomePeak2是专为MeRIP-seq设计的Bioconductor R包，考虑了转录组结构（只在外显子区域call peak），有GC含量校正，使用负二项分布统计模型，支持差异m6A分析。MACS2是通用的ChIP-seq peak caller，用于MeRIP-seq需要大量参数调整（--nomodel, --gsize等），但灵活性更高。2024年m6ACali显示两者都存在抗体非特异性结合导致的假阳性，可以用机器学习校准提高准确性。

Q5：m6A研究的最新趋势是什么？¶

答：①单碱基分辨率方法——m6A-SAC-seq、GLORI等化学方法取代抗体，实现单碱基精度且低输入量；②Nanopore直接RNA测序——m6Anet等深度学习模型从Nanopore电流信号直接检测m6A，不需要任何化学处理；③单细胞m6A——scDART-seq等方法实现单细胞水平的m6A检测；④m6A与疾病——METTL3抑制剂等表观转录组靶向药物进入临床试验。MeRIP-seq虽然仍是主流，但正在被更精确的方法补充和替代。