764. RNA修饰m6A检测分析¶

一句话概述：检测mRNA上的m6A甲基化修饰，揭示基因表达的"表观转录组"调控层——就像发现mRNA上也有"便利贴"标记，告诉细胞"这条信息要加急处理/马上销毁/换个地方翻译"。

核心知识点速查表¶

概念	白话解释	关键工具
m6A	N6-甲基腺嘌呤，mRNA上最丰富的修饰	调控RNA代谢
DRACH	m6A的共识序列(D=A/G/U, R=A/G)	约70%的m6A在此
Writer	m6A甲基转移酶(添加m6A)	METTL3/METTL14
Eraser	m6A去甲基化酶(去除m6A)	FTO/ALKBH5
Reader	m6A识别蛋白(读取m6A)	YTHDF1/2/3
MeRIP-seq	抗体富集+测序检测m6A	最经典方法
GLORI	化学法检测m6A(单碱基分辨率)	2023年Nature

一、原理（白话版）¶

1.1 m6A是什么？¶

m6A = 腺嘌呤(A)的第6位氮上加了一个甲基(-CH3)

m6A的角色（类比"便利贴"）：
  📝 加速降解：YTHDF2读到m6A → 招募降解机器 → mRNA被快速降解
  📝 促进翻译：YTHDF1读到m6A → 招募翻译起始因子 → 翻译效率↑
  📝 影响剪接：YTHDC1读到m6A → 改变剪接模式
  📝 核输出：YTHDC1促进m6A-mRNA从核到胞质的运输

m6A的动态调控：
  Writer(写入器): METTL3/METTL14/WTAP → 添加m6A
  Eraser(擦除器): FTO/ALKBH5 → 去除m6A
  Reader(阅读器): YTHDF1/2/3, YTHDC1/2 → 识别并执行功能

m6A的生物学重要性：
  - 胚胎发育必需
  - 干细胞分化调控
  - 免疫应答调节
  - 肿瘤发生发展（METTL3在多种癌症中异常）

1.2 主要检测方法¶

方法	分辨率	原理	定量
MeRIP-seq	~100bp峰	抗体富集m6A片段	半定量
miCLIP	单碱基	交联+免疫沉淀	半定量
GLORI	单碱基	化学转化(C→T)	定量(化学计量)
m6A-SAC-seq	单碱基	酶法标记	定量
Nanopore dRNA	单碱基	电流信号差异	定量
SELECT	单位点	引物延伸	定量验证

二、MeRIP-seq分析流程¶

2.1 数据处理¶

# ===== MeRIP-seq数据处理 =====
# MeRIP-seq有两组数据：Input(总RNA) 和 IP(m6A抗体富集)

# Step 1: 质控与比对
# 对Input和IP分别处理
for sample in input_rep1 input_rep2 ip_rep1 ip_rep2; do
    # 质控
    fastp \
      -i ${sample}_R1.fastq.gz \        # 输入R1
      -I ${sample}_R2.fastq.gz \        # 输入R2
      -o ${sample}_clean_R1.fq.gz \     # 输出R1
      -O ${sample}_clean_R2.fq.gz \     # 输出R2
      -h ${sample}_fastp.html           # QC报告

    # 比对到基因组
    STAR \
      --runThreadN 8 \                  # 线程数
      --genomeDir star_index/ \         # STAR索引
      --readFilesIn ${sample}_clean_R1.fq.gz ${sample}_clean_R2.fq.gz \
      --readFilesCommand zcat \         # 解压缩
      --outSAMtype BAM SortedByCoordinate \  # 排序BAM
      --outFileNamePrefix ${sample}_ \  # 输出前缀
      --outFilterMultimapNmax 1         # 只保留唯一比对

    # 建索引
    samtools index ${sample}_Aligned.sortedByCoord.out.bam
done

2.2 Peak calling（核心步骤）¶

# ===== MACS2 Peak Calling =====
# m6A MeRIP-seq的peak calling类似ChIP-seq
macs2 callpeak \
  -t ip_rep1_sorted.bam ip_rep2_sorted.bam \     # IP样本（处理组）
  -c input_rep1_sorted.bam input_rep2_sorted.bam \  # Input（对照）
  -f BAMPE \                            # 配对端BAM
  --nomodel \                           # 不建模（RNA不是组蛋白）
  --keep-dup all \                      # 保留所有重复
  -g hs \                               # 基因组大小（hs=human）
  -q 0.05 \                             # FDR阈值
  --outdir macs2_output \               # 输出目录
  -n m6a_peaks                          # 输出名称

# ===== exomePeak2（专门为MeRIP-seq设计的peak caller）=====
# R语言
# BiocManager::install("exomePeak2")

# ===== R语言：exomePeak2分析 =====
library(exomePeak2)  # 加载exomePeak2

# 运行peak calling
peaks <- exomePeak2(
  bam_ip = c("ip_rep1.bam", "ip_rep2.bam"),       # IP样本BAM
  bam_input = c("input_rep1.bam", "input_rep2.bam"),  # Input BAM
  gff = "annotation.gtf",                           # 基因注释
  genome = "hg38",                                   # 基因组版本
  paired_end = TRUE                                  # 配对端
)

# 结果在 exomePeak2_output/ 目录
# 包含BED格式的peak文件和差异分析结果

2.3 差异m6A分析¶

# ===== 差异m6A分析：处理 vs 对照 =====
library(exomePeak2)

# 差异peak calling
diff_peaks <- exomePeak2(
  bam_ip = c("treat_ip_1.bam", "treat_ip_2.bam"),     # 处理组IP
  bam_input = c("treat_input_1.bam", "treat_input_2.bam"),  # 处理组Input
  bam_ip_treated = c("ctrl_ip_1.bam", "ctrl_ip_2.bam"),  # 对照组IP
  bam_input_treated = c("ctrl_input_1.bam", "ctrl_input_2.bam"),  # 对照组Input
  gff = "annotation.gtf",
  genome = "hg38",
  paired_end = TRUE
)

# 输出：差异甲基化位点
# log2FC > 0: 处理组m6A增加(hyper-methylation)
# log2FC < 0: 处理组m6A减少(hypo-methylation)

2.4 Motif分析与可视化¶

# ===== Python分析m6A peak =====
import pandas as pd  # 导入pandas
import matplotlib.pyplot as plt  # 导入matplotlib
import pybedtools  # 导入pybedtools

# 读取MACS2 peak结果
peaks = pd.read_csv(
    "macs2_output/m6a_peaks_peaks.narrowPeak",
    sep="\t", header=None,
    names=["chr", "start", "end", "name", "score", "strand",
           "signalValue", "pvalue", "qvalue", "peak"]
)
print(f"检测到 {len(peaks)} 个m6A peaks")

# Peak在基因上的分布（metagene分析）
# m6A富集在3'UTR和Stop codon附近是经典特征
# 用Guitar或MetaPlotR做metagene plot

# ===== DRACH motif搜索 =====
# 提取peak序列
peak_bed = pybedtools.BedTool("macs2_output/m6a_peaks_peaks.narrowPeak")

# 用HOMER做motif分析
# findMotifsGenome.pl peaks.bed hg38 motif_output/ -rna -size 100

# 手动统计DRACH motif
import re  # 导入正则

drach_pattern = re.compile(r'[AGT][AG]AC[ACT]')  # DRACH motif
# D=[AGT], R=[AG], A=A, C=C, H=[ACT]

# 统计peak中DRACH出现频率
# 预期：m6A peaks中DRACH motif显著富集

# ===== R语言：Guitar包做metagene图 =====
library(Guitar)  # 加载Guitar包
library(TxDb.Hsapiens.UCSC.hg38.knownGene)

# 加载peak
peaks_gr <- import("m6a_peaks.bed")  # BED转GRanges

# 绘制metagene图
GuitarPlot(
  gfeinger = peaks_gr,                 # peak坐标
  txdb = TxDb.Hsapiens.UCSC.hg38.knownGene,  # 转录本数据库
  saveToPDFprefix = "m6a_metagene"     # 保存PDF
)
# 预期结果：m6A peak在3'UTR和CDS末端显著富集

三、GLORI化学法（单碱基分辨率）¶

# ===== GLORI分析流程概述 =====
# GLORI(glyoxal and nitrite-mediated deamination of unmethylated adenosines)
# 原理：未甲基化的A被化学转化为I(读为G)，m6A不被转化(仍为A)
# 类比：类似亚硫酸盐测序检测DNA甲基化

# Step 1: 比对（需要特殊处理A→G转化）
# 类似BS-seq，需要用支持A→G转化的比对器

# Step 2: 计算甲基化水平
# 每个A位点：m6A level = A读段数 / (A读段数 + G读段数)
# m6A level > 10% 且有统计学显著性 → 认为是m6A位点

# GLORI的优势：
# ① 单碱基分辨率（MeRIP-seq只能到~100bp峰）
# ② 定量（化学计量比，不是抗体富集的半定量）
# ③ 不依赖抗体（减少批次差异）

四、常见报错与解决¶

报错信息	原因	解决方案
`MACS2: too few peaks`	IP效率低	检查富集效率(>2倍)
`exomePeak2: no peaks`	数据量不够	增加测序深度(≥30M reads)
`No DRACH enrichment`	抗体非特异性	检查抗体质量，考虑换抗体
`Input/IP不匹配`	样本搞混	检查样本对应关系
`Peak太宽`	片段化不够	确保RNA片段在100-150nt
`低重复性`	技术噪声大	至少2个生物学重复

五、面试高频问题¶

Q1: m6A的Writer/Eraser/Reader各是什么？¶

A: Writer(甲基转移酶)：METTL3-METTL14复合物是核心催化组件，WTAP/VIRMA等辅助蛋白决定位置特异性。Eraser(去甲基化酶)：FTO和ALKBH5能去除m6A，使修饰可逆。Reader(识别蛋白)：YTHDF1促翻译、YTHDF2促降解、YTHDF3协同1和2、YTHDC1影响剪接和核输出。

Q2: MeRIP-seq和GLORI的核心区别？¶

A: MeRIP-seq用抗体富集m6A片段，分辨率约100bp（峰水平），半定量。GLORI用化学转化区分m6A和非修饰A，分辨率为单碱基，定量（化学计量比）。MeRIP-seq更成熟、文章多；GLORI更精确但实验难度大。

Q3: m6A在肿瘤中的作用？¶

A: m6A是把"双刃剑"：①促癌：METTL3在AML/肝癌等中过表达，m6A促进癌基因mRNA稳定性或翻译；②抑癌：在某些癌症中m6A促进抑癌基因表达。关键取决于m6A修饰的靶基因和识别它的Reader蛋白。这使得m6A成为潜在的治疗靶点。

六、速查表¶

# ===== m6A分析速查 =====

# MeRIP-seq流程
STAR比对 → MACS2/exomePeak2 peak calling → motif分析 → metagene图

# Peak calling
macs2 callpeak -t ip.bam -c input.bam -f BAMPE \
    --nomodel --keep-dup all -q 0.05

# R: exomePeak2
exomePeak2(bam_ip=ip_bams, bam_input=input_bams,
           gff="anno.gtf", genome="hg38")

# m6A经典特征：
# DRACH motif: [AGU][AG]AC[ACU]
# 富集位置: 3'UTR + Stop codon附近
# 每个mRNA平均3-5个m6A位点

# Writer/Eraser/Reader:
# Writer: METTL3/14 + WTAP
# Eraser: FTO, ALKBH5
# Reader: YTHDF1(翻译) YTHDF2(降解) YTHDC1(剪接)