163_全长转录本分析IsoSeq¶

# === isoseq3流水线 ===
# 输入: PacBio SMRT Cell数据 (subreads.bam)

# Step 1: CCS生成 (circular consensus sequencing)
# pbccs已集成到isoseq3中
ccs subreads.bam ccs.bam \
  --min-rq 0.99 \
  --min-passes 3 \
  --num-threads 32

# Step 2: Primer去除和全长reads识别
# 寻找5'和3'引物 → 判定全长
isoseq3 tag ccs.bam flt.bam --design T-12U-16B
# 或经典模式:
lima ccs.bam primers.fasta fl.bam --isoseq

# Step 3: 去除poly-A和concatemer
isoseq3 refine fl.bam primers.fasta flnc.bam \
  --require-polya  # 要求有poly-A(确保全长)

# FLNC = Full-Length Non-Chimeric reads

# Step 4: 聚类(可选，用于降低冗余)
isoseq3 cluster flnc.bam polished.bam \
  --verbose --use-qvs

# Step 5: 比对到参考基因组
# 使用minimap2 (splice-aware模式)
minimap2 -ax splice:hq \
  -uf \
  --secondary=no \
  reference.fa \
  polished.hq.fasta.gz \
  | samtools sort -@ 8 -o aligned.bam

samtools index aligned.bam

# 或直接对FLNC比对(跳过cluster步骤，保留更多信息)
minimap2 -ax splice:hq -uf --secondary=no \
  reference.fa flnc.bam \
  | samtools sort -@ 8 -o flnc_aligned.bam

# === ONT cDNA测序处理 ===
# 使用pychopper识别全长reads(有5'和3'引物)
pychopper -r report.pdf \
  -S stats.tsv \
  -t 16 \
  input.fastq.gz \
  full_length.fastq.gz

# 比对到基因组 (splice-aware)
minimap2 -ax splice \
  -uf \
  -k14 \
  --secondary=no \
  --junc-bed known_junctions.bed \
  reference.fa \
  full_length.fastq.gz \
  | samtools sort -@ 8 -o ont_aligned.bam

samtools index ont_aligned.bam

# === ONT direct RNA处理 ===
# Direct RNA不需要pychopper(无PCR引物)
# 但需要注意reads是3'→5'方向
minimap2 -ax splice \
  -uf \
  -k14 \
  --secondary=no \
  reference.fa \
  drna_reads.fastq.gz \
  | samtools sort -@ 8 -o drna_aligned.bam

# === 安装IsoQuant ===
pip install isoquant
# 或 conda install -c bioconda isoquant

# === 运行IsoQuant ===
isoquant.py \
  --reference reference.fa \
  --genedb annotation.gtf \
  --bam aligned.bam \
  --data_type pacbio_ccs \
  --output isoquant_output/ \
  --threads 16 \
  --complete_genedb

# 参数说明:
# --data_type: pacbio_ccs / nanopore / assembly
# --complete_genedb: 输出完整基因数据库(含新发现)
# --genedb: 已知注释(GTF/GFF3)

# === ONT数据 ===
isoquant.py \
  --reference reference.fa \
  --genedb annotation.gtf \
  --fastq full_length.fastq.gz \
  --data_type nanopore \
  --output isoquant_ont/ \
  --threads 16

# === 输出文件 ===
# isoquant_output/
#   *.transcript_models.gtf     : 重构的转录本模型
#   *.read_assignments.tsv      : 每条read的转录本归属
#   *.transcript_counts.tsv     : 转录本计数
#   *.gene_counts.tsv           : 基因计数
#   *.transcript_tpm.tsv        : TPM量化

# === 安装SQANTI3 ===
conda create -n sqanti3 python=3.9
conda activate sqanti3
pip install sqanti3

# === 运行SQANTI3 QC ===
python sqanti3_qc.py \
  isoquant_output/transcript_models.gtf \
  annotation.gtf \
  reference.fa \
  --output sqanti3_output/ \
  --cage_peak cage_peaks.bed \
  --polyA_motif_list polyA_list.txt \
  --short_reads illumina.bam \
  -t 16

# 输入:
# 1. 新发现的转录本GTF (来自IsoQuant/isoseq3/TAMA)
# 2. 参考注释GTF
# 3. 参考基因组FASTA
# 可选: CAGE峰(验证5'端)、poly-A信号、短读长比对(验证junction)

# === SQANTI3转录本分类 ===
# FSM (Full Splice Match): 完全匹配已知转录本
# ISM (Incomplete Splice Match): 匹配已知转录本但5'/3'端缺失
# NIC (Novel In Catalog): 新剪接组合但所有junction已知
# NNC (Novel Not in Catalog): 含至少一个新junction
# Antisense: 反义转录本
# Genic Intron: 完全位于内含子内
# Genic Genomic: 部分覆盖外显子但无剪接
# Intergenic: 基因间区域
# Fusion: 跨越多个基因的融合转录本

# === SQANTI3结果分析 (R) ===
library(ggplot2)
library(dplyr)

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

# 转录本类别分布
category_counts <- sqanti %>%
  group_by(structural_category) %>%
  summarise(n = n()) %>%
  mutate(pct = n / sum(n) * 100)

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) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "SQANTI3 Transcript Classification", x = "", y = "Count")

# 过滤高质量转录本
hq_transcripts <- sqanti %>%
  filter(structural_category %in% c("FSM", "ISM", "NIC", "NNC")) %>%
  filter(all_canonical == "canonical") %>%  # 经典剪接位点
  filter(!is.na(FL_count) & FL_count >= 3)  # 至少3条全长reads支持

cat(sprintf("高质量转录本: %d / %d (%.1f%%)\n", 
            nrow(hq_transcripts), nrow(sqanti), 
            nrow(hq_transcripts)/nrow(sqanti)*100))

library(IsoformSwitchAnalyzeR)

# === 分析可变剪接事件 ===
# IsoformSwitchAnalyzeR整合全长转录本数据
switchList <- importIsoformExpression(
  isoformCountMatrix = "transcript_counts.tsv",
  designMatrix = data.frame(
    sampleID = c("sample1", "sample2", "sample3", "control1", "control2", "control3"),
    condition = c(rep("Tumor", 3), rep("Normal", 3))
  )
)

# 添加注释
switchList <- importRdata(
  isoformCountMatrix = count_matrix,
  designMatrix = design,
  isoformExonAnnoation = "transcript_models.gtf",
  isoformNtFasta = "transcript_sequences.fa"
)

# 差异isoform使用分析
switchList <- isoformSwitchTestDEXSeq(switchList, dIFcutoff = 0.1)

# 分析开放阅读框
switchList <- analyzeORF(switchList)

# 分析isoform switch的功能后果
switchList <- analyzeSwitchConsequences(switchList)

# 可视化特定基因的isoform switch
switchPlot(switchList, gene = "BRCA1")

# === 完整Iso-Seq分析流程 ===
# PacBio Iso-Seq
isoseq3 tag ccs.bam flt.bam --design T-12U-16B
isoseq3 refine flt.bam primers.fasta flnc.bam --require-polya
minimap2 -ax splice:hq -uf reference.fa flnc.bam | samtools sort -o aligned.bam
isoquant.py --reference reference.fa --genedb anno.gtf --bam aligned.bam \
  --data_type pacbio_ccs --output results/ --threads 16
python sqanti3_qc.py results/transcript_models.gtf anno.gtf reference.fa

# ONT全长转录本
pychopper input.fq.gz full_length.fq.gz
minimap2 -ax splice -uf reference.fa full_length.fq.gz | samtools sort -o ont.bam
isoquant.py --reference reference.fa --genedb anno.gtf --bam ont.bam \
  --data_type nanopore --output ont_results/ --threads 16