Nanopore直接RNA测序（RNA修饰检测/polyA长度/全长转录本）¶

# ===== Dorado basecalling（推荐，已完全取代Guppy） =====
# Dorado是ONT唯一维护的basecaller（最新v1.4.0，2026年2月）
# 注意：Dorado v1.0+已移除Fast5和RNA002支持，仅支持POD5格式和RNA004

# 下载Dorado（从ONT官网或GitHub releases）
# https://github.com/nanoporetech/dorado/releases

# 下载RNA模型（最新v5.3.0）
dorado download --model rna004_130bps_sup@v5.3.0

# Basecalling（SUP模型，最高准确率）
dorado basecaller rna004_130bps_sup@v5.3.0 pod5_dir/ \
  --emit-fastq > basecalled.fastq

# 或输出BAM格式（推荐，保留信号信息）
dorado basecaller rna004_130bps_sup@v5.3.0 pod5_dir/ \
  --reference genome.fa > aligned.bam

# 可用RNA模型速度/准确率等级：
# rna004_130bps_fast@v5.3.0  - 最快，准确率较低
# rna004_130bps_hac@v5.3.0   - 速度与准确率均衡
# rna004_130bps_sup@v5.3.0   - 最高准确率（推荐）

# SUP模型支持的RNA修饰检测（无需额外工具）：
# m5C_2OmeC, m6A_DRACH, inosine_m6A_2OmeA, pseU_2OmeU, 2OmeG
# 使用方法：dorado basecaller rna004_130bps_sup@v5.3.0 pod5_dir/ \
#   --modified-bases m6A_DRACH

# ===== Guppy basecalling（已停止维护，仅旧数据使用） =====
# Guppy已被Dorado取代，新项目不建议使用
# 如果有旧Fast5数据，需先用pod5 convert将Fast5转为POD5格式
# pip install pod5
# pod5 convert fast5 fast5_dir/ --output converted.pod5
guppy_basecaller \
  --input_path fast5_dir/ \
  --save_path output_dir/ \
  --config rna_r9.4.1_70bps_hac.cfg \
  --device cuda:0 \
  --recursive \
  --compress_fastq

# ===== 质控 =====
# 使用NanoPlot
pip install NanoPlot

NanoPlot --fastq basecalled.fastq \
  --outdir nanoplot_qc/ \
  --threads 8 \
  --plots dot kde

# 查看关键指标：
# - Mean read quality (推荐 >Q8 for dRNA)
# - Mean read length (dRNA一般500-2000nt)
# - Number of reads
# - N50 read length

# 使用pycoQC（支持直接读取sequencing_summary.txt）
pip install pycoQC
pycoQC -f sequencing_summary.txt -o pycoQC_report.html

# ===== 质量过滤 =====
# 使用NanoFilt
pip install nanofilt

NanoFilt --quality 7 --length 100 --maxlength 100000 \
  < basecalled.fastq > filtered.fastq

# ===== Minimap2比对（推荐用于长读长） =====
# 比对到参考基因组
minimap2 -ax splice \
  --secondary=no \
  -uf \
  -k14 \
  --MD \
  -t 16 \
  reference_genome.fa \
  filtered.fastq | \
  samtools sort -@ 8 -o aligned.sorted.bam

samtools index aligned.sorted.bam

# 比对到转录组（如果只需定量）
minimap2 -ax map-ont \
  -t 16 \
  transcriptome.fa \
  filtered.fastq | \
  samtools sort -@ 8 -o aligned_transcriptome.sorted.bam

# ===== 比对统计 =====
samtools flagstat aligned.sorted.bam
# 关键指标：mapping rate(dRNA通常>85%)

# ===== 全长转录本鉴定（FLAIR） =====
# FLAIR: Full-Length Alternative Isoform analysis of RNA
pip install flair-brookslab

# 步骤1：Correct - 校正splice junctions
flair correct \
  -q aligned.sorted.bed \
  -g reference_genome.fa \
  -f annotation.gtf \
  -o corrected \
  --threads 16

# 步骤2：Collapse - 合并reads为isoforms
flair collapse \
  -r filtered.fastq \
  -q corrected_all_corrected.bed \
  -g reference_genome.fa \
  -f annotation.gtf \
  -o collapsed \
  --threads 16 \
  --generate_map

# 步骤3：Quantify - 定量isoform表达
flair quantify \
  -r reads_manifest.tsv \
  -i collapsed.isoforms.fa \
  --threads 16 \
  -o quantified

# 步骤4：DiffExp - 差异isoform分析
flair diffExp \
  -q quantified.counts.tsv \
  -o diff_results \
  --threads 16

# ===== FLAMES（另一个全长转录本工具） =====
# FLAMES支持单细胞和bulk
# R包：BiocManager::install("FLAMES")
# Python: pip install FLAMES-py

# ===== StringTie2 长读长模式 =====
stringtie aligned.sorted.bam \
  -L \
  -G annotation.gtf \
  -o stringtie_output.gtf \
  -p 16 \
  -A gene_abundances.tab

# ===== 方法1：m6Anet（深度学习方法，推荐） =====
pip install m6anet

# 步骤1：数据预处理（需要原始信号）
# 确保BAM文件包含move table或使用Nanopolish
m6anet dataprep \
  --input_dir eventalign_output/ \
  --output_dir m6anet_dataprep/ \
  --n_processes 16

# 步骤2：预测m6A位点
m6anet inference \
  --input_dir m6anet_dataprep/ \
  --output_dir m6anet_results/ \
  --n_processes 16 \
  --num_iterations 1000

# 对RNA004数据，使用专门训练的模型：
# m6anet inference \
#   --input_dir m6anet_dataprep/ \
#   --out_dir m6anet_results/ \
#   --pretrained_model HEK293T_RNA004 \
#   --n_processes 16

# 结果文件：data.site_proba.csv
# 包含每个DRACH motif位点的m6A概率

# ===== 方法2：Nanopolish eventalign =====
# 需要原始fast5文件和basecalled序列

# 建立索引
nanopolish index \
  -d fast5_dir/ \
  basecalled.fastq

# 比对
minimap2 -ax splice -uf -k14 reference.fa basecalled.fastq | \
  samtools sort -o aligned.bam
samtools index aligned.bam

# Event alignment（将信号对齐到参考序列）
nanopolish eventalign \
  --reads basecalled.fastq \
  --bam aligned.bam \
  --genome reference.fa \
  --signal-index \
  --scale-events \
  --threads 16 > eventalign.txt

# ===== 方法3：xPore（比较两个条件） =====
pip install xpore

# xPore通过比较两组样本的信号差异检测修饰
# 需要control（如METTL3 KO）和treatment（WT）
xpore dataprep \
  --eventalign eventalign_control.txt \
  --out_dir xpore_control/

xpore dataprep \
  --eventalign eventalign_treatment.txt \
  --out_dir xpore_treatment/

xpore diffmod \
  --config xpore_config.yml \
  --n_processes 16

# ===== 方法4：ELIGOS（检测多种修饰） =====
pip install eligos

# ELIGOS通过比较dRNA信号与预期信号检测修饰
eligos2 rna_mod \
  -i aligned.bam \
  -reg regions.bed \
  -ref reference.fa \
  -o eligos_output/ \
  -p 16 \
  --oddR 1.5 \
  --esb 0.1 \
  --pval 0.05

# ===== 方法5：Tombo（ONT旧工具） =====
pip install ont-tombo

# 重新归正信号
tombo resquiggle fast5_dir/ reference.fa \
  --processes 16 \
  --rna

# 检测修饰（与标准模型比较）
tombo detect_modifications de_novo \
  --fast5-basedirs fast5_dir/ \
  --statistics-file-basename tombo_stats \
  --processes 16

# 输出可视化
tombo plot most_significant \
  --fast5-basedirs fast5_dir/ \
  --statistics-filename tombo_stats.tombo.stats \
  --plot-standard-model \
  --num-bases 10 \
  --num-regions 5

# ===== Nanopolish polya =====
# 需要原始fast5信号文件
nanopolish polya \
  --reads basecalled.fastq \
  --bam aligned.bam \
  --genome transcriptome.fa \
  --threads 16 > polya_results.tsv

# 输出包含每条read的polyA长度估计
# 列：readname, contig, position, leader_start, adapter_start, 
#     polya_start, transcript_start, read_rate, polya_length, qc_tag

# 过滤高质量结果
awk '$NF == "PASS"' polya_results.tsv > polya_pass.tsv

# ===== tailfindr（R包） =====
# BiocManager::install("tailfindr")
# 支持DNA和RNA的polyA测量

# ===== Dorado内置polyA估计 =====
dorado basecaller rna004_130bps_sup@v5.3.0 pod5_dir/ \
  --estimate-poly-a > basecalled_with_polya.bam

# polyA长度信息存储在BAM的pt tag中
samtools view basecalled_with_polya.bam | \
  awk '{for(i=12;i<=NF;i++) if($i~/^pt:/) print $1, $i}'

# ===== R中分析polyA长度 =====

# R分析polyA结果
library(ggplot2)
library(dplyr)

# 读取nanopolish polya结果
polya <- read.table("polya_pass.tsv", header = TRUE, sep = "\t")

# 基本统计
summary(polya$polya_length)

# 按基因聚合
polya_by_gene <- polya %>%
  group_by(contig) %>%
  summarise(
    median_polya = median(polya_length),
    mean_polya = mean(polya_length),
    n_reads = n(),
    sd_polya = sd(polya_length)
  ) %>%
  filter(n_reads >= 10)  # 至少10条reads

# 全局polyA长度分布
ggplot(polya, aes(x = polya_length)) +
  geom_histogram(bins = 100, fill = "steelblue") +
  xlim(0, 300) +
  labs(x = "Poly(A) Tail Length (nt)", y = "Count",
       title = "Distribution of Poly(A) Tail Length") +
  theme_classic()

# 两组比较（如不同条件下的polyA长度变化）
polya$group <- ifelse(grepl("control", polya$readname), "Control", "Treatment")
ggplot(polya, aes(x = polya_length, fill = group)) +
  geom_density(alpha = 0.5) +
  xlim(0, 300) +
  labs(x = "Poly(A) Tail Length", y = "Density") +
  theme_classic()

# 特定基因的polyA长度分布
gene_of_interest <- "ACTB"
polya_gene <- polya[polya$contig == gene_of_interest, ]
ggplot(polya_gene, aes(x = polya_length)) +
  geom_histogram(bins = 30) +
  labs(title = paste("Poly(A) length of", gene_of_interest))

# ===== IsoQuant（isoform定量和发现） =====
pip install isoquant

# 已知isoform定量
isoquant.py \
  --reference reference_genome.fa \
  --genedb annotation.gtf \
  --fastq filtered.fastq \
  --data_type nanopore \
  --output isoquant_output/ \
  --threads 16 \
  --complete_genedb

# 新isoform发现
isoquant.py \
  --reference reference_genome.fa \
  --genedb annotation.gtf \
  --fastq filtered.fastq \
  --data_type nanopore \
  --output isoquant_novel/ \
  --threads 16 \
  --transcript_quantification unique_only

# ===== SQANTI3（isoform质量评估） =====
# 对发现的isoform进行分类和质量评估
python sqanti3_qc.py \
  collapsed.isoforms.gtf \
  annotation.gtf \
  reference_genome.fa \
  --output sqanti_output \
  --dir sqanti_dir/ \
  -t 16

# SQANTI3将isoform分类为：
# FSM: Full Splice Match（完全匹配已知）
# ISM: Incomplete Splice Match（部分匹配）
# NIC: Novel In Catalog（新组合，已知junction）
# NNC: Novel Not in Catalog（含新junction）
# Fusion, Genic, Antisense, Intergenic等

# ===== 差异isoform使用（DTU分析） =====
# 使用DRIMSeq + stageR

# R中进行DTU（Differential Transcript Usage）分析
library(DRIMSeq)
library(stageR)

# 读取isoform定量数据
counts <- read.table("isoquant_output/transcript_counts.tsv",
                     header = TRUE, row.names = 1)

# 准备DRIMSeq数据
d <- dmDSdata(counts = counts, samples = sample_info)

# 过滤低表达isoform
d <- dmFilter(d, min_samps_gene_expr = 3, min_samps_feature_expr = 3,
              min_gene_expr = 10, min_feature_expr = 5)

# 设计矩阵
design <- model.matrix(~ condition, data = DRIMSeq::samples(d))

# 差异分析
d <- dmPrecision(d, design = design)
d <- dmFit(d, design = design)
d <- dmTest(d, coef = "conditionTreatment")

# 结果
res_gene <- DRIMSeq::results(d)  # 基因水平
res_txp <- DRIMSeq::results(d, level = "feature")  # 转录本水平

# stageR两阶段检验
pScreen <- res_gene$pvalue
names(pScreen) <- res_gene$gene_id
pConfirmation <- matrix(res_txp$pvalue, ncol = 1)
rownames(pConfirmation) <- res_txp$feature_id

stageObj <- stageRTx(pScreen = pScreen, pConfirmation = pConfirmation,
                     pScreenAdjusted = FALSE, tx2gene = tx2gene_df)
stageObj <- stageWiseAdjustment(stageObj, method = "dtu", alpha = 0.05)
dtu_results <- getAdjustedPValues(stageObj, order = TRUE)

#!/bin/bash
# ============================================
# Nanopore Direct RNA-seq Analysis Pipeline
# ============================================

# === 0. 环境设置 ===
THREADS=16
GENOME="reference_genome.fa"
ANNOTATION="annotation.gtf"
TRANSCRIPTOME="transcriptome.fa"
RAW_DIR="pod5_pass/"
OUTPUT="drna_analysis/"
mkdir -p $OUTPUT/{basecall,qc,align,isoform,modification,polya}

# === 1. Basecalling ===
dorado basecaller rna004_130bps_sup@v5.3.0 $RAW_DIR \
  --estimate-poly-a \
  --emit-fastq > $OUTPUT/basecall/basecalled.fastq

# === 2. Quality Control ===
NanoPlot --fastq $OUTPUT/basecall/basecalled.fastq \
  --outdir $OUTPUT/qc/ \
  --threads $THREADS

# === 3. Filtering ===
NanoFilt --quality 7 --length 100 \
  < $OUTPUT/basecall/basecalled.fastq \
  > $OUTPUT/basecall/filtered.fastq

# === 4. Alignment ===
minimap2 -ax splice -uf -k14 --MD -t $THREADS \
  $GENOME $OUTPUT/basecall/filtered.fastq | \
  samtools sort -@ 8 -o $OUTPUT/align/aligned.sorted.bam
samtools index $OUTPUT/align/aligned.sorted.bam

# === 5. Alignment stats ===
samtools flagstat $OUTPUT/align/aligned.sorted.bam > $OUTPUT/align/flagstat.txt

# === 6. Isoform identification (IsoQuant) ===
isoquant.py \
  --reference $GENOME \
  --genedb $ANNOTATION \
  --bam $OUTPUT/align/aligned.sorted.bam \
  --data_type nanopore \
  --output $OUTPUT/isoform/ \
  --threads $THREADS

# === 7. m6A detection (m6Anet) ===
# 需要先做eventalign
nanopolish index -d $RAW_DIR $OUTPUT/basecall/filtered.fastq

nanopolish eventalign \
  --reads $OUTPUT/basecall/filtered.fastq \
  --bam $OUTPUT/align/aligned.sorted.bam \
  --genome $TRANSCRIPTOME \
  --signal-index --scale-events \
  --threads $THREADS > $OUTPUT/modification/eventalign.txt

m6anet dataprep \
  --eventalign $OUTPUT/modification/eventalign.txt \
  --out_dir $OUTPUT/modification/m6anet_prep/ \
  --n_processes $THREADS

m6anet inference \
  --input_dir $OUTPUT/modification/m6anet_prep/ \
  --out_dir $OUTPUT/modification/m6anet_results/ \
  --n_processes $THREADS

# === 8. Poly(A) length ===
nanopolish polya \
  --reads $OUTPUT/basecall/filtered.fastq \
  --bam $OUTPUT/align/aligned.sorted.bam \
  --genome $TRANSCRIPTOME \
  --threads $THREADS > $OUTPUT/polya/polya_results.tsv

echo "Pipeline completed!"