20260910 110123 Rna Seq 差异表达分析全流程实现详解

# 对6个样本的R1/R2共12个FASTQ文件运行FastQC
fastqc -t 8 -o qc/raw/ raw_data/*.fastq.gz

multiqc qc/raw/ -o qc/multiqc/ -f

from Bio import SeqIO
from collections import defaultdict

# 确认编码格式：2011年后的Illumina数据均为Phred+33（'fastq'）
# 千万不要猜编码——从仪器元数据确认

# 计算每个read位置的平均质量，识别3'端衰减
position_quals = defaultdict(list)
for record in SeqIO.parse('sample1_R1.fastq', 'fastq'):
    for i, q in enumerate(record.letter_annotations['phred_quality']):
        position_quals[i].append(q)

for pos in sorted(position_quals)[:20]:
    quals = position_quals[pos]
    print(f'Position {pos}: mean={sum(quals)/len(quals):.1f}')

# NovaSeq RNA-seq 预处理：轻修剪（aligner 会 soft-clip，不要 --dedup）
fastp -i raw_data/sample1_R1.fq.gz -I raw_data/sample1_R2.fq.gz \
      -o clean/sample1_R1.fq.gz -O clean/sample1_R2.fq.gz \
      --detect_adapter_for_pe \
      --trim_poly_g \
      -q 20 -l 50 \
      -w 8 \
      -h clean/sample1.html -j clean/sample1.json

fastqc -t 8 -o qc/trimmed/ clean/*.fastq.gz
multiqc qc/ -o qc/compare/ -f   # 同时对比 raw vs trimmed

# sjdbOverhang = readlength - 1 = 150 - 1 = 149
STAR --runMode genomeGenerate \
     --genomeDir star_index/ \
     --genomeFastaFiles GRCh38.fa \
     --sjdbGTFfile gencode.v44.annotation.gtf \
     --sjdbOverhang 149 \
     --runThreadN 16

for sample in sample1 sample2 sample3 sample4 sample5 sample6; do
  STAR --genomeDir star_index/ \
       --readFilesIn clean/${sample}_R1.fq.gz clean/${sample}_R2.fq.gz \
       --readFilesCommand zcat \
       --outFileNamePrefix aligned/${sample}_ \
       --outSAMtype BAM SortedByCoordinate \
       --quantMode GeneCounts \
       --twopassMode Basic \
       --runThreadN 8
done

featureCounts -a gencode.v44.annotation.gtf \
              -o counts.txt \
              -s 2 \
              -p --countReadPairs \
              -T 8 \
              aligned/sample1_Aligned.sortedByCoord.out.bam \
              aligned/sample2_Aligned.sortedByCoord.out.bam \
              aligned/sample3_Aligned.sortedByCoord.out.bam \
              aligned/sample4_Aligned.sortedByCoord.out.bam \
              aligned/sample5_Aligned.sortedByCoord.out.bam \
              aligned/sample6_Aligned.sortedByCoord.out.bam

cat counts.txt.summary

library(DESeq2)
library(edgeR)

counts <- read.csv('counts.txt', sep='\t', row.names=1, skip=1)
counts <- counts[, 6:11]  # 去掉前5列元数据列，只保留计数列
coldata <- data.frame(condition = factor(c('ctrl','ctrl','ctrl','treat','treat','treat')),
                      row.names = colnames(counts))

dds <- DESeqDataSetFromMatrix(countData = counts, colData = coldata, design = ~condition)

colSums(counts(dds))        # 每样本总reads
colSums(counts(dds) > 0)    # 每样本检出基因数

# 用 edgeR 的 filterByExpr 做有原则的过滤，优于硬编码 CPM>1
keep <- filterByExpr(counts(dds), group = dds$condition)
dds <- dds[keep, ]
cat("保留基因数:", sum(keep), "，过滤基因数:", sum(!keep), "\n")

# VST（方差稳定变换）用于 QC 可视化，不是用于 DE 检验
# n>30 用 vst()，n<30 用 rlog()——本例 n=6 偏小，但 vst 通常也能用
vsd <- vst(dds, blind = TRUE)   # blind=TRUE 用于无监督QC

plotPCA(vsd, intgroup = 'condition')

# 样本距离热图
pheatmap::pheatmap(as.matrix(dist(t(assay(vsd)))), main = 'Sample distance (VST)')

dds <- DESeq(dds)
cooks <- assays(dds)[['cooks']]   # 注意：不是 results(dds)$cooksd
boxplot(log10(cooks), las = 2, main = "Cook's distance per sample")

# 已知批次：放入 design 公式，让引擎在原始计数层面校正
design(dds) <- ~ batch + condition   # condition 在最后 = 感兴趣的对比

# 绝不要对 DE 输入数据使用 removeBatchEffect() 或 ComBat 后再喂给 DESeq2
# 那样会双重校正

library(DESeq2)

# 重新构建（如果上面改了 design）
dds <- DESeqDataSetFromMatrix(countData = counts, colData = coldata, design = ~condition)
dds <- DESeq(dds)

# 提取结果：显式命名对比，alpha=0.05
res <- results(dds, contrast = c('condition', 'treat', 'ctrl'), alpha = 0.05)
summary(res)

# LFC 收缩（apeglm）：降噪但不改变 p 值
res_shrunk <- lfcShrink(dds, coef = 'condition_treat_vs_ctrl', type = 'apeglm')

library(DESeq2)
library(EnhancedVolcano)
library(ggplot2)

plotDispEsts(dds)

res_df <- as.data.frame(res)
ggplot(res_df, aes(x = pvalue)) + geom_histogram(bins = 50) + theme_bw()

# 未收缩 vs 收缩对比
par(mfrow = c(1, 2))
plotMA(res, ylim = c(-5, 5), main = 'Unshrunken')
plotMA(res_shrunk, ylim = c(-5, 5), main = 'apeglm shrunk')

library(EnhancedVolcano)
EnhancedVolcano(res_shrunk,
    lab = rownames(res_shrunk),
    x = 'log2FoldChange',
    y = 'pvalue',          # 注意：用未收缩的 p 值
    pCutoff = 0.05,
    FCcutoff = 1,
    max.overlaps = Inf)

# 取 top 50 差异基因
top_genes <- head(order(res$padj), 50)
mat <- assay(vsd)[top_genes, ]
pheatmap::pheatmap(mat, scale = 'row',     # 行z-score：看pattern
                   annotation_col = coldata,
                   show_rownames = FALSE)

# 查看某个关键基因的表达模式
plotCounts(dds, gene = 'TP53', intgroup = 'condition')

library(clusterProfiler)
library(org.Hs.eg.db)

# 用 DESeq2 的 stat（Wald 统计量）作为排名指标
# 关键：用带符号的、方差校准的统计量，不用原始 p 值
gl <- res$stat
names(gl) <- rownames(res)
gl <- sort(gl, decreasing = TRUE)

# 去掉 NA
gl <- gl[!is.na(gl)]

# GO Biological Process
set.seed(42)   # 必须设种子，保证可重复
gse_go <- gseGO(geneList = gl,
                OrgDb = org.Hs.eg.db,
                ont = 'BP',
                minGSSize = 10,
                maxGSSize = 500,
                seed = TRUE)
gse_go <- setReadable(gse_go, OrgDb = org.Hs.eg.db)

# KEGG 通路
gse_kegg <- gseKEGG(geneList = gl,
                    organism = 'hsa',
                    minGSSize = 10,
                    maxGSSize = 500,
                    seed = TRUE)

head(gse_go@result, 10)

# 气泡图
dotplot(gse_go, showCategory = 20) + ggtitle('GSEA GO BP')

# GSEA 运行得分曲线
library(enrichplot)
gseaplot(gse_go, by = 'count', geneSetID = 1)

原始 FASTQ
  │
  ├─ Step 1: FastQC + MultiQC 质控报告
  │          （RNA-seq 在 per-base content / duplication 报红是正常的）
  │
  ├─ Step 2: fastp 一步预处理
  │          （自动接头修剪 + poly-G去除 + 质量过滤，绝不 --dedup）
  │
  ├─ Step 3: STAR 比对
  │          （sjdbOverhang=149, twopassMode, --quantMode GeneCounts）
  │
  ├─ Step 4: featureCounts 定量
  │          （-s 2 反向链, -p 按fragment计数, 验证Assigned>70%）
  │
  ├─ Step 5: 计数矩阵 QC
  │          （filterByExpr过滤, VST变换, PCA聚类, Cook's离群值, 批次检查）
  │
  ├─ Step 6: DESeq2 差异检验
  │          （Wald检验, alpha=0.05, apeglm LFC收缩, 显式contrast）
  │
  ├─ Step 7: 可视化
  │          （离散度图, p值直方图, MA图, 火山图, 热图, plotCounts）
  │
  └─ Step 8: GSEA 通路富集
             （用 stat 排名, clusterProfiler gseGO/gseKEGG, NES+FDR判定）