RNA-seq 差异表达分析流程：从 fastq 到 DEG

1
2
3
4
5

# FastQC 质控
fastqc -t 8 -o qc/ *.fastq.gz

# MultiQC 汇总
multiqc qc/ -o multiqc_report/

1
2
3
4
5
6
7

trim_galore --paired \
  --quality 20 \
  --length 35 \
  --max_n 0 \
  --trim-n \
  -o trimmed/ \
  sample_R1.fastq.gz sample_R2.fastq.gz

1
2
3
4
5
6

hisat2 -p 8 \
  -x genome_index/genome \
  -1 trimmed/sample_R1_val_1.fq.gz \
  -2 trimmed/sample_R2_val_2.fq.gz \
  -S align/sample.sam \
  2> align/sample_hisat2.log

1
2
3

samtools view -bS align/sample.sam | \
  samtools sort -@ 8 -o align/sample.sorted.bam
samtools index align/sample.sorted.bam

1
2
3
4
5

featureCounts -T 8 \
  -p -B -C \
  -a annotation.gtf \
  -o counts.txt \
  *.sorted.bam

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22

library(DESeq2)
library(clusterProfiler)

# 读取计数矩阵
countData <- read.table("counts.txt", header = TRUE, row.names = 1)
colData <- data.frame(
  condition = factor(c(rep("Control", 3), rep("Treatment", 3)))
)

# 创建 DESeq2 对象
dds <- DESeqDataSetFromMatrix(
  countData = countData,
  colData = colData,
  design = ~ condition
)

# 运行差异分析
dds <- DESeq(dds)
res <- results(dds, alpha = 0.05)

# 筛选显著差异基因
sig_genes <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

# 基因 ID 转换
gene_list <- sig_genes$log2FoldChange
names(gene_list) <- rownames(sig_genes)

# GO 富集
go_enrich <- enrichGO(
  gene = names(gene_list),
  OrgDb = org.Hs.eg.db,
  keyType = "ENSEMBL",
  ont = "BP",
  pAdjustMethod = "BH",
  pvalueCutoff = 0.05
)