长读测序全长转录本分析¶

# ===== PacBio Iso-Seq 全长转录本分析流程 =====

# 1. 生成 CCS 读（如果从 subreads 开始）
ccs input.subreads.bam ccs.bam \
    --min-rq 0.9 --min-passes 1 --num-threads 16

# 2. 引物去除和全长分类（isoseq3）
# 识别 5' 和 3' 引物，分类为全长（FL）或非全长
isoseq3 tag ccs.bam tagged.bam --design T-12U-16B
isoseq3 refine tagged.bam primers.fasta flnc.bam \
    --require-polya  # 要求 polyA 尾巴

# flnc = Full-Length Non-Chimeric reads（全长非嵌合读）
# 这些是高质量的全长 cDNA 序列

# 3. 聚类：将相似的全长读聚成转录本
isoseq3 cluster flnc.bam clustered.bam \
    --verbose --use-qvs

# 4. 比对到参考基因组
# 使用 minimap2 的 splice 模式（专门处理 RNA 比对）
minimap2 -a -x splice:hq \
    -uf \             # 强制正义链比对
    --secondary=no \  # 不输出次优比对
    -t 16 \
    hg38.fa clustered.hq.fasta.gz | \
    samtools sort -@ 8 -o aligned.sorted.bam
samtools index aligned.sorted.bam

# 5. 转录本塌陷/去冗余（cDNA_Cupcake 或 IsoQuant）
# 方法A：使用 IsoQuant（推荐，更现代）
isoquant.py \
    --reference hg38.fa \
    --genedb gencode.v44.annotation.gtf \
    --bam aligned.sorted.bam \
    --data_type pacbio_ccs \
    --output isoquant_output \
    --threads 16

# 方法B：cDNA_Cupcake（经典方法）
python collapse_isoforms_by_sam.py \
    --input clustered.hq.fasta.gz \
    --sam aligned.sorted.sam \
    --prefix collapsed \
    --min-coverage 0.85 \  # 至少 85% 的读比对上
    --min-identity 0.95    # 至少 95% 序列一致性

# 6. SQANTI3 质量评估和分类
python sqanti3_qc.py \
    collapsed.gff \
    gencode.v44.annotation.gtf \
    hg38.fa \
    -o sqanti_output \
    --report pdf \
    -t 16 \
    --SR_bam short_read.sorted.bam  # 可选：用短读验证剪接位点

# 7. SQANTI3 过滤（去除低质量转录本）
python sqanti3_filter.py rules \
    sqanti_output_classification.txt \
    sqanti_output_corrected.fasta \
    sqanti_output_corrected.gtf \
    -o filtered

# ===== R: 全长转录本下游分析 =====
library(ggplot2)
library(dplyr)

# 读取 SQANTI3 分类结果
sqanti <- read.delim("sqanti_output_classification.txt")

# 统计各类转录本的比例
category_counts <- sqanti %>%
    count(structural_category) %>%
    mutate(pct = n / sum(n) * 100) %>%
    arrange(desc(n))

# 可视化转录本分类
ggplot(category_counts,
       aes(x = reorder(structural_category, -n), y = n, fill = structural_category)) +
    geom_bar(stat = "identity") +
    geom_text(aes(label = paste0(round(pct, 1), "%")), vjust = -0.5) +
    labs(x = "转录本类别", y = "数量",
         title = "SQANTI3 转录本分类") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))

# 统计新发现的可变剪接事件
novel <- sqanti %>%
    filter(structural_category %in% c("NIC", "NNC"))
cat("新剪接组合（NIC）:", sum(novel$structural_category == "NIC"), "\n")
cat("全新剪接位点（NNC）:", sum(novel$structural_category == "NNC"), "\n")

# 每个基因的转录本亚型数量
isoform_per_gene <- sqanti %>%
    group_by(associated_gene) %>%
    summarise(n_isoforms = n()) %>%
    arrange(desc(n_isoforms))

cat("平均每个基因的转录本数:", mean(isoform_per_gene$n_isoforms), "\n")
cat("转录本最多的基因:", isoform_per_gene$associated_gene[1],
    "(", isoform_per_gene$n_isoforms[1], "个)\n")