核糖体图谱Ribo-seq¶

Ribo-seq的核心思想：
  翻译中的核糖体像"火车"沿着mRNA"铁轨"行驶
  核糖体"覆盖"的mRNA片段受保护，不被核酸酶降解

实验步骤：
  1. 用翻译抑制剂（如CHX）冻住核糖体在mRNA上的位置
  2. 用核酸酶（RNase I）消化未被保护的RNA
  3. 纯化核糖体-RNA复合物
  4. 提取RPF（约28-32nt的小片段）
  5. 建库测序
  6. 同时做RNA-seq获取mRNA丰度

关键质控指标：
  - RPF长度分布应集中在28-32nt
  - 三核苷酸周期性（triplet periodicity）：reads按3nt周期分布
  - CDS区富集：>80%的reads应该在编码区
  - 起始密码子和终止密码子处有明显的reads堆积

1. 翻译组学（Translatome）
   哪些mRNA正在被翻译？翻译到什么程度？

2. 翻译效率（Translation Efficiency, TE）
   TE = RPF_RPKM / mRNA_RPKM
   TE高: mRNA翻译活跃（可能有IRES或优势密码子）
   TE低: mRNA存在翻译抑制（可能有uORF或RNA二级结构）

3. 新ORF发现
   lncRNA上的小ORF → 小蛋白/微肽
   5'UTR中的uORF → 调控主ORF翻译
   内含子保留区的ORF → 新蛋白变体

4. 密码子使用偏好
   核糖体在某些密码子上停留更久 → 翻译暂停
   稀有密码子 → 翻译速度减慢
   与tRNA丰度相关

# 1. 质控
fastqc RPF_sample.fq.gz     # Ribo-seq数据质控
fastqc mRNA_sample.fq.gz    # 配对RNA-seq质控

# 2. 去接头（Ribo-seq reads通常很短，约28-32nt）
cutadapt \
  -a AGATCGGAAGAGCACACGTCT \  # 3'接头序列（根据建库方法调整）
  -m 25 \                     # 最小保留长度25nt
  -M 35 \                     # 最大保留长度35nt（RPF应该28-32nt）
  -o RPF_trimmed.fq.gz \      # 输出文件
  RPF_sample.fq.gz            # 输入文件

# 3. 去除rRNA reads（Ribo-seq中rRNA污染通常很高）
# 建立rRNA索引
bowtie2-build rRNA_sequences.fa rRNA_index  # 建rRNA索引

# 比对到rRNA，保留未比对的reads
bowtie2 \
  -x rRNA_index \              # rRNA索引
  -U RPF_trimmed.fq.gz \       # 输入
  --un-gz RPF_norRNA.fq.gz \   # 未比对到rRNA的reads（保留）
  -S /dev/null \               # 比对结果丢弃
  -p 8                         # 线程数

echo "去除rRNA前后reads数量比较："
echo "去除前: $(zcat RPF_trimmed.fq.gz | wc -l) / 4"
echo "去除后: $(zcat RPF_norRNA.fq.gz | wc -l) / 4"

# 4. 比对到基因组
STAR --runThreadN 8 \
  --genomeDir star_index/ \              # STAR索引
  --readFilesIn RPF_norRNA.fq.gz \       # 去rRNA后的reads
  --readFilesCommand zcat \              # gzip文件
  --outSAMtype BAM SortedByCoordinate \  # 输出排序BAM
  --outFilterMismatchNmax 2 \            # 最多允许2个错配
  --outFilterMultimapNmax 1 \            # 只保留唯一比对
  --alignIntronMax 20000 \               # 最大内含子长度
  --outFileNamePrefix RPF_               # 前缀

samtools index RPF_Aligned.sortedByCoord.out.bam  # 索引BAM

# Ribo-seq特有的质控指标
library(riboWaltz)  # Ribo-seq专用R包

# 读取BAM文件
reads_list <- bamtolist(
  bamfolder = "bam_folder/",  # BAM文件目录
  annotation = annotation_df   # 基因注释数据框
)

# 1. RPF长度分布（应该集中在28-32nt）
rlength_distr(reads_list,
              sample = "RPF_sample",
              cl = 99)  # 99%置信区间

# 典型的好数据：28nt有一个明显的主峰
# 如果长度分布太散 → RNA降解或建库问题

# 2. 三核苷酸周期性（最重要的质控）
# reads应该按3nt周期分布在CDS区
metaprofile_psite(reads_list,
                  annotation = annotation_df,
                  sample = "RPF_sample",
                  length_filter_mode = "periodicity")

# 三核苷酸周期性来源：
# 核糖体每翻译一个密码子（3nt）前进一步
# 所以RPF reads的5'端位置应该呈3nt周期分布

# 3. P-site偏移计算
# P-site是核糖体实际翻译的位置
# 需要从RPF reads的5'端推算P-site位置
psite_offset <- psite(reads_list,
                      flanking = 6,      # 侧翼区域
                      extremity = "5end") # 基于5'端计算

# 4. 基因组区域分布
# CDS区应该>80%，5'UTR和3'UTR比例较低
region_psite(reads_list,
             annotation = annotation_df,
             sample = "RPF_sample")

# 翻译效率 = RPF / mRNA
# 需要配对的Ribo-seq和RNA-seq数据

library(DESeq2)

# 准备计数矩阵
# RPF计数（来自Ribo-seq）
rpf_counts <- featureCounts(
  files = "RPF_sample.bam",
  annot.ext = "annotation.gtf",
  isGTFAnnotationFile = TRUE,
  isPairedEnd = FALSE           # Ribo-seq通常是单端
)$counts

# mRNA计数（来自RNA-seq）
mrna_counts <- featureCounts(
  files = "mRNA_sample.bam",
  annot.ext = "annotation.gtf",
  isGTFAnnotationFile = TRUE,
  isPairedEnd = TRUE            # RNA-seq通常是双端
)$counts

# 方法1：简单比值法
rpf_rpkm <- rpkm(rpf_counts, gene_lengths)    # 计算RPF的RPKM
mrna_rpkm <- rpkm(mrna_counts, gene_lengths)   # 计算mRNA的RPKM
TE <- rpf_rpkm / mrna_rpkm                     # 翻译效率

# 方法2：DESeq2差异翻译分析（推荐，更严格）
# 合并RPF和RNA-seq计数矩阵
count_matrix <- cbind(rpf_counts, mrna_counts)  # 合并

# 创建条件信息
col_data <- data.frame(
  condition = rep(c("treated", "control"), each = 4),  # 实验条件
  assay = rep(c("RPF", "mRNA"), times = 4),            # 数据类型
  row.names = colnames(count_matrix)
)

# DESeq2分析
dds <- DESeqDataSetFromMatrix(
  countData = count_matrix,
  colData = col_data,
  design = ~ condition + assay + condition:assay  # 交互项=翻译效率变化
)
dds <- DESeq(dds)

# 交互项结果就是翻译效率的差异
te_results <- results(dds, name = "conditiontreated.assayRPF")

# 密码子占用率反映核糖体在每个密码子上停留的时间
library(riboWaltz)

# 计算每个密码子的P-site密度
codon_usage <- codon_usage_psite(
  reads_list,
  annotation = annotation_df,
  sample = "RPF_sample",
  fastapath = "hg38.fa"       # 参考基因组（提取密码子序列）
)

# 可视化密码子占用率
# 某些密码子占用率高 = 核糖体在那里停留久 = 翻译暂停
library(ggplot2)
ggplot(codon_usage, aes(x = reorder(codon, -occupancy), y = occupancy)) +
  geom_bar(stat = "identity", fill = "steelblue") +  # 柱状图
  labs(x = "Codon", y = "Occupancy",
       title = "Codon Occupancy (Ribosome Pausing)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, size = 6))  # 旋转标签

# 高占用率密码子通常是：
# - 稀有密码子（对应tRNA少）
# - 脯氨酸密码子（CCC, CCU等，翻译慢）
# - 起始密码子（AUG）和终止密码子附近

# Ribo-seq最独特的应用：发现新的翻译事件

# 1. 使用ORFquant检测活跃翻译的ORF
# BiocManager::install("ORFquant")
library(ORFquant)

# 从Ribo-seq数据中识别翻译的ORF
orf_results <- find_orfs(
  bam_file = "RPF_sample.bam",      # Ribo-seq BAM
  annotation = annotation_gtf,       # 基因注释
  genome = BSgenome.Hsapiens.UCSC.hg38,  # 基因组
  min_reads = 10,                    # 最少10条reads
  min_length = 10                    # 最短10个密码子
)

# 分类新发现的ORF
# uORF: 5'UTR中的上游ORF（调控主ORF翻译）
# dORF: 3'UTR中的下游ORF
# sORF: lncRNA上的小ORF（编码微肽）
# overlapping ORF: 与已知CDS重叠的ORF

# 2. 使用RibORF（另一个常用工具）
# RibORF用三核苷酸周期性和reads均匀性评分ORF的翻译活性

RiboToolkit是免费的Web工具（http://rnabioinfor.tch.harvard.edu/RiboToolkit）

适合不会编程的初学者，提供：
1. 数据预处理（去接头、去rRNA、比对）
2. RPF质控（长度分布、三核苷酸周期性）
3. 基于RPF的表达定量
4. 翻译效率分析
5. 密码子占用率分析
6. 差异翻译分析
7. 功能注释
8. 活跃翻译ORF鉴定

使用方法：
  上传fastq文件 → 选择物种和参数 → 自动运行 → 下载结果

RiboSeq.Org（https://riboseq.org）是2025年NAR发布的综合平台

包含多个互联工具：
1. GWIPS-viz — 全基因组核糖体图谱可视化
2. Trips-Viz — 转录组层面的Ribo-seq可视化
3. RiboGalaxy — 数据处理（Galaxy界面）
4. RDP（RiboSeq Data Portal）— 数据门户

RDP目前收录：
  - 14,840个测序文库
  - 969个研究项目
  - 96个物种
  所有数据预处理完毕，可直接使用

RiboParser/RiboShiny是2025-2026年新发布的分析框架

特点：
1. 优化的P-site检测方法（基于核糖体结构和起始/终止位点）
2. 全面的质控、基因级分析、密码子级分析
3. 支持Ribo-seq变体分析
4. 适用于非模式生物
5. 比RiboToolkit更高的准确性
6. RiboShiny提供交互式可视化

# Ribo-seq分析流程
去接头(cutadapt) → 去rRNA(bowtie2) → 比对(STAR)
→ 质控(riboWaltz/RiboseQC) → P-site计算
→ 翻译效率(DESeq2) → 密码子占用率 → 新ORF发现

# 质控标准
RPF长度: 集中在28-32nt
三核苷酸周期性: CDS区有明显3nt周期
CDS比例: >80%的reads在编码区
rRNA比例: <20%（越低越好）
起始/终止密码子: 有明显的reads峰

# 核心公式
翻译效率(TE) = RPF_RPKM / mRNA_RPKM
TE > 1: 翻译活跃
TE < 1: 翻译抑制
TE ≈ 1: 翻译与转录一致

# 工具选择
初学者/无编程: RiboToolkit（Web）
全面分析: RiboSeq.Org（Web+数据库）
自定义分析: riboWaltz（R包）
差异翻译: DESeq2（交互项模型）
新ORF发现: ORFquant/RibORF

# 密码子占用率含义
高占用 = 核糖体停留久 = 翻译暂停
原因: 稀有密码子、tRNA不足、RNA结构、蛋白折叠