765. 翻译组学Ribo-seq分析¶

实验步骤：
  1. 用RNase消化未被核糖体保护的mRNA → 只保留核糖体"脚下"的~28-30nt片段
  2. 用蔗糖梯度离心纯化核糖体+片段复合物
  3. 提取RNA片段 → 建库 → 测序

生物学意义：
  mRNA转录 ≠ 蛋白翻译
  - 有些mRNA高表达但翻译很少（翻译抑制）
  - 有些mRNA低表达但翻译很多（翻译增强）
  - Ribo-seq直接看哪些mRNA正在被翻译

翻译效率(TE) = Ribo-seq信号 / RNA-seq信号
  TE > 1: 翻译效率高于平均
  TE < 1: 翻译效率低于平均
  TE变化: 翻译调控（不改变mRNA水平的基因表达调控）

三碱基周期性(3-nt periodicity)：
  核糖体每次移动3个碱基(一个密码子)
  → Ribo-seq读段在CDS中显示3nt的周期性模式
  → 这是正在翻译的"指纹信号"
  → 用于验证数据质量和识别翻译ORF

P-site偏移(P-site offset)：
  RPF的5'端与核糖体P-site之间有固定偏移
  通常：28nt RPF → P-site在5'端第13个碱基
  精确的P-site定位 → 单密码子分辨率的翻译图谱

# ===== Ribo-seq基础数据处理 =====

# Step 1: 质控与去接头
fastp \
  -i sample_ribo.fastq.gz \            # 输入Ribo-seq reads
  -o sample_ribo_clean.fq.gz \         # 输出清洁reads
  --length_required 25 \               # 最短25nt
  --length_limit 35 \                  # 最长35nt（RPF在28-30nt）
  -h sample_ribo_qc.html               # QC报告

# Step 2: 去除rRNA（Ribo-seq中rRNA污染严重）
bowtie2 \
  -x rRNA_index \                      # rRNA参考索引
  -U sample_ribo_clean.fq.gz \         # 输入reads
  --un-gz sample_norRNA.fq.gz \        # 未比对的reads（非rRNA）
  -S /dev/null \                       # 丢弃比对结果
  -p 8 \                              # 8线程
  --very-sensitive                     # 高敏感度

echo "去除rRNA后reads数:"
zcat sample_norRNA.fq.gz | wc -l | awk '{print $1/4}'

# Step 3: 比对到基因组
STAR \
  --runThreadN 8 \
  --genomeDir star_index/ \            # STAR索引
  --readFilesIn sample_norRNA.fq.gz \  # 非rRNA reads
  --readFilesCommand zcat \
  --outSAMtype BAM SortedByCoordinate \
  --outFilterMultimapNmax 1 \          # 只保留唯一比对
  --outFilterMismatchNmax 2 \          # 最多2个错配
  --alignIntronMax 50000 \             # 最大内含子长度
  --outFileNamePrefix sample_ribo_

samtools index sample_ribo_Aligned.sortedByCoord.out.bam

# ===== 检查Ribo-seq数据质量 =====
import pysam  # 导入pysam
import matplotlib.pyplot as plt  # 导入matplotlib
import numpy as np  # 导入numpy
from collections import Counter  # 导入Counter

# 读取BAM
bam = pysam.AlignmentFile("sample_ribo_sorted.bam", "rb")  # 打开BAM

# 统计RPF长度分布
lengths = []  # 存储读段长度
for read in bam:
    if not read.is_unmapped and read.mapping_quality >= 20:  # 过滤
        lengths.append(read.query_length)  # 记录长度

length_dist = Counter(lengths)  # 长度分布

# 可视化RPF长度
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 左图：长度分布
x = sorted(length_dist.keys())
y = [length_dist[k] for k in x]
axes[0].bar(x, y, color="steelblue", edgecolor="black")
axes[0].set_xlabel("RPF Length (nt)")
axes[0].set_ylabel("Read Count")
axes[0].set_title("RPF Length Distribution")
axes[0].axvline(x=28, color="red", linestyle="--", label="28nt")
axes[0].axvline(x=30, color="red", linestyle="--", label="30nt")
axes[0].legend()

# 右图：metagene图（起始密码子附近的读段分布）
# 预期：在AUG起始密码子处有明显的peak
# 且显示3nt周期性
axes[1].set_xlabel("Distance from Start Codon (nt)")
axes[1].set_ylabel("Read Density")
axes[1].set_title("Metagene Profile at Start Codon")

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

# 质量标准：
# ① RPF长度峰值在28-30nt → 说明RNase消化成功
# ② CDS区域占比>60% → reads主要在编码区
# ③ 明显的3nt周期性 → 真实的翻译信号

# ===== 计算翻译效率(TE) =====
import pandas as pd  # 导入pandas

# 读取Ribo-seq和RNA-seq的计数
ribo_counts = pd.read_csv("ribo_gene_counts.tsv", sep="\t", index_col=0)
rna_counts = pd.read_csv("rnaseq_gene_counts.tsv", sep="\t", index_col=0)

# 归一化为RPKM
def rpkm(counts, gene_lengths):
    """计算RPKM"""
    total = counts.sum()  # 总reads数
    rpkm = counts / gene_lengths * 1e9 / total  # RPKM公式
    return rpkm

# 计算TE
# TE = Ribo-seq RPKM / RNA-seq RPKM
te = ribo_counts["rpkm"] / rna_counts["rpkm"]  # 翻译效率
te = te.replace([np.inf, -np.inf], np.nan).dropna()  # 去无穷值

print(f"平均TE: {te.mean():.2f}")
print(f"TE > 2的基因数: {(te > 2).sum()}")  # 翻译效率高的基因
print(f"TE < 0.5的基因数: {(te < 0.5).sum()}")  # 翻译效率低的基因

# 差异翻译效率分析
# 使用Xtail或RiboDiff
# Xtail: 用负二项模型检测条件间的TE变化

# ===== R语言：Xtail差异翻译分析 =====
library(xtail)  # 加载Xtail

# 准备数据
ribo_matrix <- as.matrix(ribo_counts)  # Ribo-seq count矩阵
rna_matrix <- as.matrix(rna_counts)    # RNA-seq count矩阵
condition <- c("ctrl", "ctrl", "treat", "treat")  # 条件

# 运行Xtail
result <- xtail(
  rna_matrix,     # RNA-seq counts
  ribo_matrix,    # Ribo-seq counts
  condition,      # 条件
  bins = 10000    # bin数
)

# 提取差异翻译基因
summary_result <- resultsTable(result)
sig_genes <- summary_result[summary_result$pvalue_final < 0.05, ]
print(paste("差异翻译基因数:", nrow(sig_genes)))

# ===== RiboTISH：发现新的翻译ORF =====
# pip install ribotish

# 预测翻译起始位点和ORF
ribotish predict \
  -b sample_ribo_sorted.bam \          # Ribo-seq BAM
  -g annotation.gtf \                  # 基因注释
  -f reference.fa \                    # 参考基因组
  -o ribotish_orfs.txt \               # 输出ORF
  --framebest \                        # 选最佳阅读框
  -p 8                                 # 线程

# 输出列包括：
# GenomePos: ORF在基因组的位置
# Start/Stop: 起始/终止密码子位置
# Tid: 转录本ID
# Symbol: 基因名
# TisType: 起始类型(annotated/novel)
# RiboPvalue: 核糖体信号的p值

# 筛选高置信新ORF
awk '$NF < 0.01 && $6 == "novel"' ribotish_orfs.txt \
  > novel_orfs_sig.txt  # p<0.01的新ORF

# ===== Ribo-seq速查 =====

# 数据处理
fastp -i ribo.fq.gz --length_required 25 --length_limit 35
bowtie2 -x rRNA_idx -U clean.fq.gz --un-gz norRNA.fq.gz
STAR --readFilesIn norRNA.fq.gz --outFilterMultimapNmax 1

# 质量标准
# RPF长度峰: 28-30nt ✓
# CDS占比: >60% ✓
# 3nt周期性: 明显 ✓
# rRNA比例: <20% ✓

# 翻译效率
# TE = Ribo-seq_RPKM / RNA-seq_RPKM
# 差异翻译: Xtail或RiboDiff

# ORF发现
ribotish predict -b ribo.bam -g anno.gtf -f ref.fa

# 关键概念：
# RPF = 核糖体保护片段(~28-30nt)
# P-site = 精确翻译位置
# TE = 翻译效率
# smORF = 小ORF(<100aa的微蛋白)