763. 单分子RNA直接测序¶

传统RNA-seq（cDNA方法）：
  RNA → 反转录为cDNA → 扩增 → 测序
  丢失的信息：
  ① RNA修饰全部丢失（反转录不保留修饰信息）
  ② poly(A)尾长度不准（反转录可能不完整）
  ③ PCR扩增引入偏差
  ④ 链方向可能丢失

直接RNA测序（dRNA-seq）：
  RNA → 直接穿过纳米孔 → 读取电流信号
  保留的信息：
  ① RNA修饰（m6A/m5C等改变电流信号）
  ② 真实的poly(A)尾长度
  ③ 无扩增偏差
  ④ 天然链方向

缺点：
  ① 通量低（比cDNA低10-100倍）
  ② 需要大量起始RNA（~500ng poly(A)+ RNA）
  ③ 错误率相对较高（~5-10%）
  ④ 成本高

# ===== ONT直接RNA测序数据处理 =====

# Step 1: Basecalling（原始信号→序列）
# 使用Dorado（ONT最新basecaller，2024+推荐）
dorado basecaller \
  rna004_130bps_sup@v5.1.0 \          # 模型名称(RNA004化学)
  pod5_directory/ \                     # 原始POD5文件目录
  --emit-fastq \                       # 输出FASTQ
  --modified-bases m6A \               # 同时检测m6A修饰
  > sample_drna.fastq                  # 输出序列

# 或输出BAM格式（包含修饰信息）
dorado basecaller \
  rna004_130bps_sup@v5.1.0 \
  pod5_directory/ \
  --reference reference.fa \           # 参考基因组
  --modified-bases m6A \               # 检测m6A
  > sample_drna.bam                    # 输出BAM（含修饰标记）

# Step 2: 比对
minimap2 \
  -ax splice \                         # 转录组比对模式（考虑剪接）
  -uf \                                # 正链方向（dRNA是正链）
  -k14 \                               # k-mer大小
  --secondary=no \                     # 不输出次级比对
  reference.fa \                       # 参考基因组
  sample_drna.fastq \                  # 输入FASTQ
  | samtools sort -@ 4 -o sample_drna_sorted.bam  # 排序

samtools index sample_drna_sorted.bam  # 建索引

# Step 3: 基本QC
NanoPlot \
  --bam sample_drna_sorted.bam \       # 输入BAM
  --outdir nanoplot_output \           # 输出目录
  -t 8 \                              # 线程
  --title "Direct RNA-seq QC"          # 标题

# ===== m6A修饰检测方法 =====

# 方法一：Dorado内置修饰检测（推荐，2025最新）
# Dorado basecaller同时输出修饰信息
dorado basecaller \
  rna004_130bps_sup@v5.1.0 \
  pod5_directory/ \
  --reference reference.fa \
  --modified-bases m6A \               # 检测m6A
  > sample_with_mods.bam

# 从BAM中提取修饰信息
dorado summary sample_with_mods.bam > summary.tsv  # 总结

# 用modkit提取修饰位点
modkit pileup \
  sample_with_mods.bam \               # 含修饰的BAM
  m6a_sites.bed \                      # 输出BED格式的修饰位点
  --ref reference.fa \                 # 参考基因组
  --filter-threshold 0.5               # 修饰概率阈值

# ===== 方法二：m6Anet检测m6A =====
# pip install m6anet

# Step 1: 准备输入数据
# m6Anet需要nanopolish eventalign的输出
# 先运行nanopolish
import subprocess

# nanopolish index
subprocess.run([
    "nanopolish", "index",
    "-d", "fast5_directory/",          # FAST5文件目录
    "sample_drna.fastq"                # FASTQ文件
])

# nanopolish eventalign
subprocess.run([
    "nanopolish", "eventalign",
    "--reads", "sample_drna.fastq",     # 读段
    "--bam", "sample_drna_sorted.bam",  # BAM文件
    "--genome", "reference.fa",         # 参考基因组
    "--signal-index",                   # 输出信号索引
    "--scale-events",                   # 标准化事件
    "--summary", "summary.txt",         # 摘要文件
    "-t", "8"                           # 线程
], stdout=open("eventalign.txt", "w"))

# Step 2: m6Anet数据准备
# m6anet dataprep --eventalign eventalign.txt --out_dir m6anet_input/ --n_processes 8

# Step 3: 运行m6Anet推断
# m6anet inference --input_dir m6anet_input/ --out_dir m6anet_output/ --n_processes 8

# Step 4: 分析结果
import pandas as pd  # 导入pandas

results = pd.read_csv("m6anet_output/data.site_proba.csv.gz")  # 读取结果
print(f"检测到 {len(results)} 个候选m6A位点")

# 筛选高置信度m6A位点
# probability_modified > 0.9 → 高置信度
high_conf = results[results["probability_modified"] > 0.9]
print(f"高置信度m6A位点: {len(high_conf)}")

# DRACH motif分析（m6A的经典序列模体）
# D=A/G/U, R=A/G, A=中心A, C=C, H=A/C/U
drach_sites = high_conf[high_conf["kmer"].str.match(r"[AGT][AG]AC[ACT]")]
print(f"DRACH motif中的m6A: {len(drach_sites)} ({len(drach_sites)/len(high_conf)*100:.1f}%)")

# ===== poly(A)尾长度分析 =====
# 使用nanopolish polya或Dorado
import pandas as pd  # 导入pandas
import matplotlib.pyplot as plt  # 导入matplotlib

# 方法一：nanopolish polya
# nanopolish polya --reads sample.fastq --bam sample.bam --genome ref.fa > polya.tsv

# 读取poly(A)结果
polya = pd.read_csv("polya.tsv", sep="\t")  # 读取制表符分隔

# 过滤高质量估计
polya_pass = polya[polya["qc_tag"] == "PASS"]  # 只要PASS的
print(f"有效poly(A)估计: {len(polya_pass)} 条读段")
print(f"平均poly(A)长度: {polya_pass['polya_length'].mean():.1f} nt")

# 可视化poly(A)分布
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 左图：全局分布
axes[0].hist(polya_pass["polya_length"], bins=100, edgecolor="black", alpha=0.7)
axes[0].set_xlabel("Poly(A) Tail Length (nt)")
axes[0].set_ylabel("Read Count")
axes[0].set_title("Global Poly(A) Length Distribution")
axes[0].axvline(x=polya_pass["polya_length"].median(), color="red",
                linestyle="--", label=f"Median: {polya_pass['polya_length'].median():.0f} nt")
axes[0].legend()

# 右图：按基因的poly(A)长度
# 合并基因注释信息
top_genes = polya_pass.groupby("contig")["polya_length"].mean().nlargest(10)
axes[1].barh(range(len(top_genes)), top_genes.values)
axes[1].set_yticks(range(len(top_genes)))
axes[1].set_yticklabels(top_genes.index)
axes[1].set_xlabel("Mean Poly(A) Length (nt)")
axes[1].set_title("Top 10 Genes by Poly(A) Length")

plt.tight_layout()
plt.savefig("polya_analysis.png", dpi=300)
plt.show()

# ===== 直接RNA测序速查 =====

# Basecalling + m6A检测
dorado basecaller rna004_model pod5_dir/ \
    --reference ref.fa --modified-bases m6A > out.bam

# 比对
minimap2 -ax splice -uf -k14 ref.fa reads.fq | samtools sort > sorted.bam

# m6Anet
nanopolish eventalign --reads reads.fq --bam sorted.bam --genome ref.fa
m6anet dataprep --eventalign eventalign.txt --out_dir prep/
m6anet inference --input_dir prep/ --out_dir results/

# poly(A)分析
nanopolish polya --reads reads.fq --bam sorted.bam --genome ref.fa

# 关键点：
# dRNA优势 = 修饰检测 + poly(A) + 无PCR偏差
# dRNA劣势 = 低通量 + 高成本 + 高错误率
# m6A DRACH motif: [AGU][AG]AC[ACU]
# 高置信m6A: probability > 0.9