转录组可变polyA分析¶

# APA示意图
# Gene X:
#   ====[CDS]=====[3'UTR_short]pA1............[3'UTR_long]pA2
#                              ↑                           ↑
#                         proximal pA                  distal pA
#
# Short isoform: CDS + short 3'UTR (使用proximal pA)
# Long isoform:  CDS + long 3'UTR  (使用distal pA)
#
# PDUI = reads_using_distal / (reads_using_proximal + reads_using_distal)
# PDUI高 → 倾向使用distal pA → 长3'UTR
# PDUI低 → 倾向使用proximal pA → 短3'UTR

# === DaPars2 分析流程 ===

# 1. 准备BedGraph文件（从BAM生成）
# 对每个样本
genomeCoverageBed -ibam sample1.sorted.bam -bg -split > sample1.bedgraph
genomeCoverageBed -ibam sample2.sorted.bam -bg -split > sample2.bedgraph

# 2. 准备基因注释（从GTF提取3'UTR）
# DaPars需要的格式：gene_symbol, chr, strand, 3UTR_start, 3UTR_end
python dapars_extract_anno.py -b hg38_refseq.bed -o extracted_3UTR.bed

# 3. 准备配置文件
cat > dapars2_config.txt << 'EOF'
# Annotated 3'UTR file
Annotated_3UTR=extracted_3UTR.bed

# BedGraph files (Group1: Control, Group2: Treatment)
Group1_Tophat_aligned_Wig=control1.bedgraph,control2.bedgraph,control3.bedgraph
Group2_Tophat_aligned_Wig=treat1.bedgraph,treat2.bedgraph,treat3.bedgraph

Output_directory=dapars2_output/
Output_result_file=APA_results

# Parameters
Num_least_in_group1=1
Num_least_in_group2=1
Coverage_cutoff=30
FDR_cutoff=0.05
PDUI_cutoff=0.2
Fold_change_cutoff=0.59
EOF

# 4. 运行DaPars2
python DaPars2_main.py dapars2_config.txt

# DaPars2结果解析
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# 读取结果
results = pd.read_csv("dapars2_output/APA_results_All_Prediction_Results.txt", sep="\t")

# 关键列说明：
# PDUI_Group1, PDUI_Group2: 各组的平均PDUI
# PDUI_Group_diff: Group2 - Group1
# P_val, adjusted.P_val: 统计检验结果
# Pass_Filter: 是否通过多重检验

# 筛选显著APA事件
sig_apa = results[results['Pass_Filter'] == 'Y'].copy()
print(f"Significant APA events: {len(sig_apa)}")

# 3'UTR缩短事件（PDUI下降 → proximal pA使用增加）
shortened = sig_apa[sig_apa['PDUI_Group_diff'] < -0.2]
# 3'UTR延长事件
lengthened = sig_apa[sig_apa['PDUI_Group_diff'] > 0.2]

print(f"3'UTR shortening: {len(shortened)} genes")
print(f"3'UTR lengthening: {len(lengthened)} genes")

# 可视化PDUI变化
plt.figure(figsize=(8, 6))
plt.scatter(results['PDUI_Group1'], results['PDUI_Group2'],
            alpha=0.3, s=5, c='grey')
plt.scatter(shortened['PDUI_Group1'], shortened['PDUI_Group2'],
            alpha=0.7, s=10, c='blue', label='Shortened')
plt.scatter(lengthened['PDUI_Group1'], lengthened['PDUI_Group2'],
            alpha=0.7, s=10, c='red', label='Lengthened')
plt.plot([0, 1], [0, 1], 'k--', alpha=0.5)
plt.xlabel("PDUI (Control)")
plt.ylabel("PDUI (Treatment)")
plt.legend()
plt.title("Alternative Polyadenylation Changes")
plt.savefig("APA_scatter.pdf")

# === QAPA 分析流程 ===

# 1. 构建QAPA参考库
# 下载PolyASite数据库
wget https://polyasite.unibas.ch/download/atlas/2.0/GRCh38.96/atlas.clusters.2.0.GRCh38.96.bed.gz

# 生成QAPA注释文件
qapa build --db ensembl_identifiers.txt \
    --gencode_fasta gencode.v38.transcripts.fa \
    -o qapa_library

# 2. 使用Salmon定量
salmon index -t qapa_library.fa -i qapa_index --threads 8

# 对每个样本定量
for sample in sample1 sample2 sample3; do
    salmon quant -i qapa_index \
        -l A \
        -1 ${sample}_R1.fq.gz -2 ${sample}_R2.fq.gz \
        -o salmon_${sample} \
        --threads 8
done

# 3. QAPA定量APA使用率
qapa quant \
    --db qapa_library.metadata.txt \
    salmon_sample1/quant.sf \
    salmon_sample2/quant.sf \
    salmon_sample3/quant.sf \
    > qapa_pau_results.txt

# === QAPA结果分析（R）===
library(dplyr)
library(ggplot2)

# 读取QAPA PAU (Poly(A) site Usage) 结果
pau <- read.table("qapa_pau_results.txt", header = TRUE, sep = "\t")

# PAU值说明：每行是一个pA位点，值为该pA位点在对应基因所有pA位点中的使用比例
# Proximal PAU + Distal PAU = 1（对于有两个pA位点的基因）

# 筛选有两个及以上pA位点的基因
multi_pa_genes <- pau %>%
  group_by(Gene) %>%
  filter(n() >= 2) %>%
  ungroup()

# 计算各样本的"长UTR使用率"（取最远端pA位点的PAU）
# 按 APA_ID 排序，每个基因的最后一个是distal
distal_pau <- multi_pa_genes %>%
  group_by(Gene) %>%
  slice_tail(n = 1) %>%
  ungroup()

# 比较条件间差异
# 假设前3列为Control，后3列为Treatment
distal_pau$mean_ctrl <- rowMeans(distal_pau[, c("sample1", "sample2", "sample3")])
distal_pau$mean_treat <- rowMeans(distal_pau[, c("sample4", "sample5", "sample6")])
distal_pau$delta_PAU <- distal_pau$mean_treat - distal_pau$mean_ctrl

# t检验
distal_pau$pvalue <- apply(distal_pau, 1, function(row) {
  ctrl_vals <- as.numeric(row[c("sample1", "sample2", "sample3")])
  treat_vals <- as.numeric(row[c("sample4", "sample5", "sample6")])
  t.test(ctrl_vals, treat_vals)$p.value
})
distal_pau$padj <- p.adjust(distal_pau$pvalue, method = "BH")

# 显著APA基因
sig_apa <- distal_pau %>% filter(padj < 0.05, abs(delta_PAU) > 0.1)

# === APA功能影响分析 ===
library(clusterProfiler)
library(org.Hs.eg.db)

# 1. 对3'UTR缩短基因做功能富集
shortened_genes <- sig_apa$Gene[sig_apa$delta_PAU < -0.1]

ego_short <- enrichGO(gene = shortened_genes,
                       OrgDb = org.Hs.eg.db,
                       keyType = "SYMBOL",
                       ont = "BP", pvalueCutoff = 0.05)
dotplot(ego_short, title = "3'UTR Shortened Genes - GO Enrichment")

# 2. 失去的miRNA结合位点分析
# 下载TargetScan预测结果
# 对每个缩短基因，找在proximal-distal之间区域的miRNA位点
library(GenomicRanges)

# 缩短区域（从proximal pA到distal pA）
lost_regions <- GRanges(
  seqnames = shortened_info$chr,
  ranges = IRanges(shortened_info$proximal_pa, shortened_info$distal_pa),
  strand = shortened_info$strand
)

# miRNA结合位点
mirna_sites <- import("targetscan_mirna_sites.bed")

# 找到在缩短区域中的miRNA位点
lost_mirna <- subsetByOverlaps(mirna_sites, lost_regions)
lost_mirna_summary <- table(lost_mirna$miRNA)
print(sort(lost_mirna_summary, decreasing = TRUE)[1:20])

# 3. APA与基因表达变化的关联
# 假设：3'UTR缩短 → 逃避miRNA降解 → 蛋白表达上调
apa_de_merge <- merge(sig_apa, de_results, by.x = "Gene", by.y = "gene_symbol")
cor.test(apa_de_merge$delta_PAU, apa_de_merge$log2FoldChange)
# 预期负相关：PDUI降低(缩短) 与 表达上调 相关

# === QuantSeq 3' mRNA-Seq 分析 ===

# QuantSeq只产生3'端的单端reads

# 1. 比对（使用STAR，注意单端模式）
STAR --runThreadN 16 \
    --genomeDir star_index \
    --readFilesIn sample_R1.fq.gz \
    --readFilesCommand zcat \
    --outSAMtype BAM SortedByCoordinate \
    --outFileNamePrefix quantseq_sample_

# 2. 使用polyAtools进行pA位点识别
# 识别reads的polyA痕迹（internal priming过滤）
python polyAtools.py \
    --bam quantseq_sample_.bam \
    --genome hg38.fa \
    --output polya_sites.bed \
    --internal-priming-filter

# 3. 使用APAtrap分析QuantSeq数据
# APAtrap专为3'端测序设计
Rscript APAtrap.R \
    --bam quantseq_sample_.bam \
    --annotation hg38.gtf \
    --output apatrap_results.txt

# 4. 使用MAAPER (Model-based Analysis of APA using PAS-seq and RNA-seq)
Rscript MAAPER.R \
    --pas_bam pas_seq.bam \
    --rna_bam rnaseq.bam \
    --gtf annotation.gtf \
    --output maaper_results.txt

# === APA调控因子分析 ===

# 关键APA调控因子
apa_regulators <- c(
  # CFIm复合物（促进distal pA使用 → 长UTR）
  "NUDT21",  # CFIm25
  "CPSF6",   # CFIm68
  # CstF复合物（促进proximal pA使用 → 短UTR）
  "CSTF2",   # CstF64
  "CSTF2T",  # CstF64tau
  "CSTF3",   # CstF77
  # CPSF复合物
  "CPSF1", "CPSF2", "CPSF3", "CPSF4",
  # 其他
  "PABPN1",  # nuclear poly(A) binding protein
  "FIP1L1",  # Fip1
  "PCF11"    # CF II
)

# 检查调控因子的差异表达
regulator_de <- de_results[de_results$gene_symbol %in% apa_regulators, ]
print(regulator_de[, c("gene_symbol", "log2FoldChange", "padj")])

# CFIm25下调 → 全局3'UTR缩短（已在多种肿瘤中验证）
# CstF64上调 → 全局3'UTR缩短

# 调控因子表达与PDUI变化的关联
cor_with_regulators <- sapply(apa_regulators, function(reg) {
  if (reg %in% rownames(expression_matrix)) {
    cor(expression_matrix[reg, ], mean_pdui_per_sample, use = "complete.obs")
  } else NA
})

# DaPars快速流程
genomeCoverageBed -ibam sample.bam -bg -split > sample.bedgraph
python DaPars2_main.py config.txt

# QAPA快速流程
qapa build --db identifiers.txt --gencode_fasta transcripts.fa -o library
salmon index -t library.fa -i index
salmon quant -i index -l A -1 R1.fq.gz -2 R2.fq.gz -o quant_out
qapa quant --db library.metadata.txt quant_out/quant.sf > pau_results.txt