转录组数据分析完整教程

一、转录组分析概述

1.1 什么是转录组测序

转录组测序（RNA-Seq）是一种利用高通量测序技术对细胞或组织中所有转录本进行定量和定性分析的方法。通过RNA-Seq，我们可以：

• 定量基因表达水平
• 发现新的转录本和可变剪接
• 识别差异表达基因
• 分析基因融合和SNP

1.2 分析流程概览

标准的RNA-Seq分析流程包括以下主要步骤：

二、Linux环境准备

2.2 Conda环境配置

Conda是生物信息学软件管理的最佳工具，推荐使用Miniconda。

安装Miniconda

# 下载Miniconda安装脚本
wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh

# 运行安装脚本
bash Miniconda3-latest-Linux-x86_64.sh

# 按照提示完成安装，然后激活conda
source ~/.bashrc

# 验证安装
conda --version

配置Bioconda通道

# 添加必要的conda通道（按此顺序）
conda config --add channels defaults
conda config --add channels bioconda
conda config --add channels conda-forge

# 设置通道优先级
conda config --set channel_priority strict

# 创建RNA-Seq专用环境
conda create -n rnaseq_env -c conda-forge python=3.10 r-base=4.3.2 -y
conda activate rnaseq_env

2.3 安装分析工具

先安装 mamba 以加速包管理，然后基于 mamba 安装RNA-Seq分析所需的所有工具：

# 激活环境
conda activate rnaseq_env

# 安装 mamba
conda install -c conda-forge mamba -y

# 基于 mamba 安装质控工具
mamba install -c bioconda fastqc multiqc -y

# 基于 mamba 安装比对工具
mamba install -c bioconda hisat2 star bwa -y

# 基于 mamba 安装SAM/BAM处理工具
mamba install -c bioconda samtools picard -y

# 基于 mamba 安装定量工具
mamba install -c bioconda subread htseq -y

# 基于 mamba 安装其他实用工具
mamba install -c bioconda trimmomatic -y

# 安装 sratoolkit（数据下载工具）
mamba install -c bioconda -c conda-forge -c defaults sratoolkit=3.1.1 -y

# 安装 salmon（转录本定量工具）
mamba install -c bioconda -c conda-forge salmon=1.10.1 -y

# 验证安装
fastqc --version
hisat2 --version
samtools --version
featureCounts --version
salmon --version
prefetch --version

三、数据获取与预处理

3.1 从NCBI下载数据

本教程使用NCBI SRA数据库中的公开数据进行演示。以人类RNA-Seq数据为例：

创建项目目录结构

# 创建项目根目录
mkdir -p ~/rnaseq_project
cd ~/rnaseq_project

# 创建子目录
mkdir -p {raw_data,ref_genome,aligned,counts,results,scripts,logs}

# 查看目录结构
tree

使用SRA Toolkit下载数据

以下是两种下载 PRJNA960807 数据集的方案：

方案一：直接使用 PRJNA960807 项目名称下载

# 进入原始数据目录
cd ~/rnaseq_project/raw_data

# 使用prefetch直接下载整个项目的数据
prefetch --option-file <(echo PRJNA960807)

# 解压下载的数据
for sra_file in *.sra; do
    echo "Processing ${sra_file}..."
    fasterq-dump --split-files ${sra_file}
done

方案二：用户输入项目名称下载

# 进入原始数据目录
cd ~/rnaseq_project/raw_data

# 提示用户输入项目名称
echo "请输入项目名称（例如：PRJNA960807）："
read project_name

# 使用prefetch下载指定项目的数据
prefetch --option-file <(echo ${project_name})

# 解压下载的数据
for sra_file in *.sra; do
    echo "Processing ${sra_file}..."
    fasterq-dump --split-files ${sra_file}
done

方案三：批量下载特定样本（选择部分样本作为示例）

# 进入原始数据目录
cd ~/rnaseq_project/raw_data

# 批量下载PRJNA960807项目的部分样本
for accession in SRX20066670 SRX20066671 SRX20066677 SRX20066676 SRX20066680 SRX20066684; do
    echo "Downloading ${accession}..."
    prefetch ${accession}
    fasterq-dump --split-files ${accession}
done

# 或者使用wget直接下载已发布的FASTQ文件
# 注意：实际下载时需要根据具体的SRA编号构建正确的URL

准备样本信息表

# 创建样本信息文件
cat > ~/rnaseq_project/samples.txt << 'EOF'
sample_id	group	fq1	fq2
normal_1	normal	SRX20066670_1.fastq.gz	SRX20066670_2.fastq.gz
normal_2	normal	SRX20066671_1.fastq.gz	SRX20066671_2.fastq.gz
normal_3	normal	SRX20066677_1.fastq.gz	SRX20066677_2.fastq.gz
stress_24h_1	stress	SRX20066676_1.fastq.gz	SRX20066676_2.fastq.gz
stress_4h_1	stress	SRX20066680_1.fastq.gz	SRX20066680_2.fastq.gz
stress_12h_1	stress	SRX20066684_1.fastq.gz	SRX20066684_2.fastq.gz
EOF

3.2 数据解压与质控

解压FASTQ文件

# 进入数据目录
cd ~/rnaseq_project/raw_data

# 解压gz文件（如果需要）
# gunzip *.gz

# 查看FASTQ文件格式
zcat SRX20066670_1.fastq.gz | head -12

# FASTQ格式说明：
# 第1行：序列标识符（以@开头）
# 第2行：碱基序列
# 第3行：+（分隔符）
# 第4行：质量分数（Phred+33编码）

使用Trimmomatic进行质控修剪

# 创建质控后数据目录
mkdir -p ~/rnaseq_project/clean_data

# 运行Trimmomatic进行质控
cd ~/rnaseq_project/raw_data

for sample in SRX20066670 SRX20066671 SRX20066677 SRX20066676 SRX20066680 SRX20066684; do
    echo "Processing ${sample}..."
    trimmomatic PE -threads 4 \
        ${sample}_1.fastq.gz ${sample}_2.fastq.gz \
        ../clean_data/${sample}_1_paired.fq.gz ../clean_data/${sample}_1_unpaired.fq.gz \
        ../clean_data/${sample}_2_paired.fq.gz ../clean_data/${sample}_2_unpaired.fq.gz \
        ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \
        LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36
done

Trimmomatic参数说明：

• PE：双端测序模式
• -threads 4：使用4个线程
• ILLUMINACLIP：去除接头序列
• LEADING:3：去除开头质量值低于3的碱基
• TRAILING:3：去除末尾质量值低于3的碱基
• SLIDINGWINDOW:4:15：4碱基窗口平均质量低于15则切除
• MINLEN:36：丢弃长度小于36bp的reads

3.3 FastQC质量检查

FastQC是用于高通量测序数据质量控制的工具，可以生成详细的质控报告。

# 创建质控报告目录
mkdir -p ~/rnaseq_project/qc_reports

# 对原始数据进行质控
cd ~/rnaseq_project/raw_data
fastqc *.fastq.gz -t 4 -o ../qc_reports/

# 对质控后的数据进行质控
cd ~/rnaseq_project/clean_data
fastqc *_paired.fq.gz -t 4 -o ../qc_reports/

# 使用MultiQC汇总所有质控报告
cd ~/rnaseq_project
multiqc qc_reports/ -o qc_reports/multiqc_report

FastQC报告解读

四、参考基因组准备

4.1 下载参考基因组

本教程以人类基因组GRCh38为例，可以从Ensembl或NCBI下载。

# 进入参考基因组目录
cd ~/rnaseq_project/ref_genome

# 从Ensembl下载人类基因组（GRCh38）
wget -c ftp://ftp.ensembl.org/pub/release-110/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz

# 解压基因组文件
gunzip Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz

# 查看基因组文件
head -20 Homo_sapiens.GRCh38.dna.primary_assembly.fa

# 统计染色体信息
grep "^>" Homo_sapiens.GRCh38.dna.primary_assembly.fa | head -10

4.2 构建HISAT2索引

HISAT2是一种快速且灵敏的RNA-Seq比对工具，需要预先构建基因组索引。

# 创建索引目录
mkdir -p ~/rnaseq_project/ref_genome/hisat2_index
cd ~/rnaseq_project/ref_genome

# 构建HISAT2索引（约需30-60分钟，需要8GB以上内存）
# 使用8个线程加速
hisat2-build -p 8 \
    Homo_sapiens.GRCh38.dna.primary_assembly.fa \
    hisat2_index/GRCh38

# 验证索引文件
ls -lh hisat2_index/

索引构建完成后，会生成以下文件：

• GRCh38.1.ht2 到 GRCh38.8.ht2：HISAT2索引文件

下载HISAT2预构建索引（可选）

如果自行构建索引耗时太长，可以下载预构建的索引：

# 下载预构建的人类基因组索引
cd ~/rnaseq_project/ref_genome
wget https://genome-idx.s3.amazonaws.com/hisat/grch38_genome.tar.gz
tar -xzf grch38_genome.tar.gz

# 重命名索引文件
mv grch38 hisat2_index

4.3 下载基因注释文件

基因注释文件（GTF/GFF）包含基因的位置信息，是表达量定量的必需文件。

# 从Ensembl下载GTF注释文件
cd ~/rnaseq_project/ref_genome

wget -c ftp://ftp.ensembl.org/pub/release-110/gtf/homo_sapiens/Homo_sapiens.GRCh38.110.gtf.gz

# 解压GTF文件
gunzip Homo_sapiens.GRCh38.110.gtf.gz

# 查看GTF文件格式
head -20 Homo_sapiens.GRCh38.110.gtf

# GTF格式说明（9列）：
# 1. 染色体名称
# 2. 来源
# 3. 特征类型（gene, transcript, exon, CDS等）
# 4. 起始位置
# 5. 结束位置
# 6. 得分
# 7. 链方向
# 8. 相位
# 9. 属性（gene_id, transcript_id, gene_name等）

提取基因信息

# 提取所有基因信息
grep -w "gene" Homo_sapiens.GRCh38.110.gtf | \
    awk -F'\t' '{print $9}' | \
    grep -o 'gene_id "[^"]*"' | \
    sed 's/gene_id "//;s/"$//' | \
    sort -u | wc -l

# 创建基因ID到基因名的映射文件
grep -w "gene" Homo_sapiens.GRCh38.110.gtf | \
    awk -F'\t' '{print $9}' | \
    grep -oP 'gene_id "\K[^"]+|gene_name "\K[^"]+' | \
    paste - - > gene_id_to_name.txt

head gene_id_to_name.txt

五、序列比对

5.1 HISAT2比对

HISAT2是基于图的比对工具，特别适合RNA-Seq数据的比对，能够处理剪接位点。

# 创建比对结果目录
cd ~/rnaseq_project

# 比对脚本
for sample in SRX20066670 SRX20066671 SRX20066677 SRX20066676 SRX20066680 SRX20066684; do
    echo "Aligning ${sample}..."
    
    hisat2 -p 8 \
        -x ref_genome/hisat2_index/GRCh38 \
        -1 clean_data/${sample}_1_paired.fq.gz \
        -2 clean_data/${sample}_2_paired.fq.gz \
        --dta \
        --dta-cufflinks \
        --new-summary \
        --summary-file logs/${sample}_align_summary.txt \
        | samtools sort -@ 4 -o aligned/${sample}.sorted.bam -
    
    # 建立BAM索引
    samtools index aligned/${sample}.sorted.bam
done

HISAT2主要参数说明：

• -p 8：使用8个线程
• -x：索引文件前缀
-1/-2：双端测序的FASTQ文件
• --dta：输出适合转录组组装的结果
• --dta-cufflinks：兼容Cufflinks的比对结果
• --new-summary：输出新的汇总格式

查看比对统计

# 查看比对汇总
cat logs/SRX20066670_align_summary.txt

# 使用samtools统计比对率
samtools flagstat aligned/SRX20066670.sorted.bam

5.2 STAR比对（可选）

STAR是另一种常用的RNA-Seq比对工具，速度更快但内存需求更高。

构建STAR索引

# 创建STAR索引目录
mkdir -p ~/rnaseq_project/ref_genome/star_index

# 构建STAR索引（需要约30GB内存）
STAR --runMode genomeGenerate \
    --genomeDir ref_genome/star_index \
    --genomeFastaFiles ref_genome/Homo_sapiens.GRCh38.dna.primary_assembly.fa \
    --sjdbGTFfile ref_genome/Homo_sapiens.GRCh38.110.gtf \
    --runThreadN 8 \
    --genomeSAindexNbases 14

运行STAR比对

# STAR比对
for sample in SRX20066670 SRX20066671 SRX20066677 SRX20066676 SRX20066680 SRX20066684; do
    echo "Aligning ${sample} with STAR..."
    
    STAR --runThreadN 8 \
        --genomeDir ref_genome/star_index \
        --readFilesIn clean_data/${sample}_1_paired.fq.gz clean_data/${sample}_2_paired.fq.gz \
        --readFilesCommand zcat \
        --outFilterType BySJout \
        --outFilterMultimapNmax 20 \
        --alignSJoverhangMin 8 \
        --alignSJDBoverhangMin 1 \
        --outFilterMismatchNmax 999 \
        --outFilterMismatchNoverReadLmax 0.04 \
        --alignIntronMin 20 \
        --alignIntronMax 1000000 \
        --alignMatesGapMax 1000000 \
        --outSAMtype BAM SortedByCoordinate \
        --outFileNamePrefix aligned/${sample}_
    
    # 重命名输出文件
    mv aligned/${sample}_Aligned.sortedByCoord.out.bam aligned/${sample}.sorted.bam
    samtools index aligned/${sample}.sorted.bam
done

5.3 比对结果评估

# 统计每个样本的比对情况
echo -e "Sample\tTotal\tMapped\tUnique\tMulti" > aligned/mapping_stats.txt

for sample in SRX20066670 SRX20066671 SRX20066677 SRX20066676 SRX20066680 SRX20066684; do
    total=$(samtools view -c aligned/${sample}.sorted.bam)
    mapped=$(samtools view -c -F 4 aligned/${sample}.sorted.bam)
    unique=$(samtools view -c -q 30 aligned/${sample}.sorted.bam)
    multi=$((mapped - unique))
    
    echo -e "${sample}\t${total}\t${mapped}\t${unique}\t${multi}" >> aligned/mapping_stats.txt
done

cat aligned/mapping_stats.txt

六、比对后处理

6.1 SAM转BAM格式

SAM格式是纯文本格式，占用空间大。BAM是其二进制压缩格式，更适合存储和处理。

# SAM转BAM并排序（如果在比对时未完成）
# samtools view: 转换格式
# -b: 输出BAM格式
# -S: 输入是SAM格式
# -h: 包含header

samtools view -bS aligned/sample.sam > aligned/sample.bam

# 按坐标排序
samtools sort aligned/sample.bam -o aligned/sample.sorted.bam

# 建立索引（用于快速访问）
samtools index aligned/sample.sorted.bam

6.2 BAM文件排序与索引

# 查看BAM文件头信息
samtools view -H aligned/SRX20066670.sorted.bam

# 查看特定区域的比对
samtools view aligned/SRX20066670.sorted.bam chr1:1000000-1001000 | head

# 统计覆盖度
samtools depth aligned/SRX20066670.sorted.bam | head

# 统计每个染色体的reads数
samtools idxstats aligned/SRX20066670.sorted.bam

6.3 去除重复序列

PCR重复会影响表达量定量的准确性，需要使用Picard或samtools去除。

# 使用Picard去除重复
for sample in SRX20066670 SRX20066671 SRX20066677 SRX20066676 SRX20066680 SRX20066684; do
    echo "Removing duplicates for ${sample}..."
    
    picard MarkDuplicates \
        I=aligned/${sample}.sorted.bam \
        O=aligned/${sample}.dedup.bam \
        M=logs/${sample}_dup_metrics.txt \
        REMOVE_DUPLICATES=true \
        CREATE_INDEX=true
done

# 或者使用samtools标记重复（不删除，只标记）
samtools markdup aligned/SRX20066670.sorted.bam aligned/SRX20066670.markdup.bam

七、表达量定量

7.1 FeatureCounts计数

FeatureCounts是Subread包中的基因计数工具，速度快且内存效率高。

# 创建计数目录
cd ~/rnaseq_project

# 运行FeatureCounts
featureCounts -T 8 \
    -p \
    -t exon \
    -g gene_id \
    -a ref_genome/Homo_sapiens.GRCh38.110.gtf \
    -o counts/gene_counts.txt \
    aligned/*.sorted.bam

# 参数说明：
# -T 8: 使用8个线程
# -p: 双端测序模式
# -t exon: 统计exon层面的计数
# -g gene_id: 按gene_id汇总
# -a: 注释文件路径
# -o: 输出文件

查看计数结果

# 查看计数结果的前几行
head -20 counts/gene_counts.txt

# 提取基因ID和计数列，生成表达矩阵
cut -f1,7- counts/gene_counts.txt | \
    grep -v "^#" | \
    sed 's/aligned\///g;s/.sorted.bam//g' > counts/count_matrix.txt

head counts/count_matrix.txt

7.2 HTSeq计数

HTSeq是另一个常用的计数工具，与DESeq2兼容性良好。

# 使用HTSeq进行计数
for sample in SRX20066670 SRX20066671 SRX20066677 SRX20066676 SRX20066680 SRX20066684; do
    echo "Counting ${sample}..."
    
    htseq-count \
        -f bam \
        -r pos \
        -s no \
        -t exon \
        -i gene_id \
        aligned/${sample}.sorted.bam \
        ref_genome/Homo_sapiens.GRCh38.110.gtf \
        > counts/${sample}_counts.txt
done

# 合并所有样本的计数结果
paste counts/*_counts.txt | \
    awk 'BEGIN{OFS="\t"} {printf "%s", $1; for(i=2;i<=NF;i+=2) printf "\t%s", $i; print ""}' \
    > counts/htseq_count_matrix.txt

7.3 表达矩阵构建

# 创建样本信息表
cat > counts/sample_info.txt << 'EOF'
sample_id	group
SRX20066670	normal
SRX20066671	normal
SRX20066677	normal
SRX20066676	stress
SRX20066680	stress
SRX20066684	stress
EOF

# 查看最终的表达矩阵
echo "表达矩阵维度："
wc -l counts/count_matrix.txt

echo -e "\n表达矩阵前10行："
head counts/count_matrix.txt

八、R语言分析流程

8.1 R环境配置

安装必要的R包

# 安装Bioconductor
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

# 安装核心分析包
BiocManager::install(c(
    "DESeq2",           # 差异表达分析
    "edgeR",            # 另一种差异分析方法
    "limma",            # 线性模型分析
    "tximport",         # 转录本定量导入
    "GenomicFeatures",  # 基因组特征注释
    "org.Hs.eg.db",     # 人类基因注释
    "AnnotationDbi"     # 注释数据库接口
))

# 安装可视化包
install.packages(c(
    "ggplot2",          # 绘图
    "pheatmap",         # 热图
    "RColorBrewer",     # 配色方案
    "EnhancedVolcano",  # 火山图
    "ggrepel",          # 标签防重叠
    "dplyr",            # 数据处理
    "tidyr",            # 数据整理
    "readr"             # 数据读取
))

8.2 数据导入与预处理

# 加载必要的包
library(DESeq2)
library(ggplot2)
library(pheatmap)
library(RColorBrewer)
library(dplyr)
library(readr)

# 设置工作目录
setwd("~/rnaseq_project")

# 读取计数矩阵
countData <- read.table("counts/count_matrix.txt", 
                        header = TRUE, 
                        row.names = 1, 
                        sep = "\t")

# 查看数据结构
dim(countData)
head(countData)

# 读取样本信息
coldata <- read.table("counts/sample_info.txt", 
                      header = TRUE, 
                      row.names = 1, 
                      sep = "\t")

# 确保样本顺序一致
coldata <- coldata[match(colnames(countData), rownames(coldata)), , drop = FALSE]

# 查看样本信息
print(coldata)

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

# 查看对象
print(dds)

# 预过滤：去除低表达的基因
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]

# 设置参考水平（对照组）
dds$group <- relevel(dds$group, ref = "control")

# 运行DESeq2分析
dds <- DESeq(dds)

# 查看结果名称
resultsNames(dds)

8.3 DESeq2差异分析

# 获取差异分析结果
res <- results(dds, contrast = c("group", "treatment", "control"))

# 查看结果摘要
summary(res)

# 查看结果表
head(res)

# 添加基因名称注释（可选）
library(org.Hs.eg.db)

# 将ENSEMBL ID转换为基因名
gene_ids <- rownames(res)
gene_symbols <- mapIds(org.Hs.eg.db, 
                       keys = gene_ids, 
                       column = "SYMBOL", 
                       keytype = "ENSEMBL",
                       multiVals = "first")

res$gene_symbol <- gene_symbols

# 按padj排序
res_ordered <- res[order(res$padj),]
head(res_ordered)

# 筛选显著差异表达基因
sig_genes <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
print(paste("显著差异表达基因数量:", nrow(sig_genes)))

# 保存结果
write.csv(as.data.frame(res_ordered), 
          file = "results/deseq2_results.csv", 
          row.names = TRUE)

write.csv(as.data.frame(sig_genes), 
          file = "results/significant_genes.csv", 
          row.names = TRUE)

标准化表达值

# 获取标准化后的计数（VST或rlog转换）
# VST（方差稳定转换）- 适合大数据集
vsd <- vst(dds, blind = FALSE)

# rlog（正则化对数转换）- 适合小数据集
# rld <- rlog(dds, blind = FALSE)

# 获取标准化表达矩阵
norm_counts <- assay(vsd)
write.csv(norm_counts, file = "results/normalized_counts.csv")

# 计算FPKM/TPM（需要基因长度信息）
# 这里简化处理，使用DESeq2的size factor标准化结果

九、数据可视化

9.1 PCA主成分分析

# PCA分析
pcaData <- plotPCA(vsd, intgroup = "group", returnData = TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))

# 绘制PCA图
pca_plot <- ggplot(pcaData, aes(x = PC1, y = PC2, color = group)) +
    geom_point(size = 3) +
    geom_text(aes(label = name), vjust = -1, size = 3) +
    xlab(paste0("PC1: ", percentVar[1], "% variance")) +
    ylab(paste0("PC2: ", percentVar[2], "% variance")) +
    ggtitle("PCA Plot of RNA-Seq Samples") +
    theme_bw() +
    theme(plot.title = element_text(hjust = 0.5))

print(pca_plot)

# 保存图片
ggsave("results/pca_plot.png", pca_plot, width = 8, height = 6, dpi = 300)

9.2 火山图绘制

# 安装EnhancedVolcano包（如果尚未安装）
if (!require("EnhancedVolcano", quietly = TRUE))
    BiocManager::install("EnhancedVolcano")

library(EnhancedVolcano)

# 绘制火山图
volcano_plot <- EnhancedVolcano(
    res,
    lab = res$gene_symbol,
    x = 'log2FoldChange',
    y = 'padj',
    title = 'Treatment vs Control',
    pCutoff = 0.05,
    FCcutoff = 1,
    pointSize = 2.0,
    labSize = 3.0,
    colAlpha = 0.5,
    legendPosition = 'right',
    legendLabSize = 10,
    legendIconSize = 3.0,
    drawConnectors = TRUE,
    widthConnectors = 0.5,
    colConnectors = 'black'
)

print(volcano_plot)

# 保存图片
ggsave("results/volcano_plot.png", volcano_plot, width = 10, height = 8, dpi = 300)

# 使用ggplot2自定义火山图
res_df <- as.data.frame(res)
res_df$significance <- "Not Significant"
res_df$significance[res_df$padj < 0.05 & res_df$log2FoldChange > 1] <- "Up-regulated"
res_df$significance[res_df$padj < 0.05 & res_df$log2FoldChange < -1] <- "Down-regulated"

volcano_gg <- ggplot(res_df, aes(x = log2FoldChange, y = -log10(padj), color = significance)) +
    geom_point(alpha = 0.6, size = 1.5) +
    scale_color_manual(values = c("Down-regulated" = "blue", 
                                   "Not Significant" = "gray", 
                                   "Up-regulated" = "red")) +
    theme_bw() +
    xlab("Log2 Fold Change") +
    ylab("-Log10 Adjusted P-value") +
    ggtitle("Volcano Plot: Differential Gene Expression") +
    geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "gray") +
    geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "gray")

print(volcano_gg)
ggsave("results/volcano_plot_custom.png", volcano_gg, width = 8, height = 6, dpi = 300)

9.3 热图绘制

# 选择前50个显著差异表达基因绘制热图
top_genes <- head(order(res$padj), 50)
mat <- assay(vsd)[top_genes, ]

# 获取基因名
row_names <- res$gene_symbol[top_genes]
row_names[is.na(row_names)] <- rownames(mat)[is.na(row_names)]
rownames(mat) <- row_names

# 数据标准化（按行）
mat_scaled <- t(scale(t(mat)))

# 定义颜色
annotation_col <- data.frame(
    Group = coldata$group
)
rownames(annotation_col) <- colnames(mat)

ann_colors <- list(
    Group = c(control = "#1f77b4", treatment = "#ff7f0e")
)

# 绘制热图
heatmap_plot <- pheatmap(
    mat_scaled,
    annotation_col = annotation_col,
    annotation_colors = ann_colors,
    cluster_rows = TRUE,
    cluster_cols = TRUE,
    show_rownames = TRUE,
    show_colnames = TRUE,
    fontsize_row = 6,
    fontsize_col = 10,
    main = "Heatmap of Top 50 Differentially Expressed Genes",
    color = colorRampPalette(rev(brewer.pal(9, "RdBu")))(100)
)

# 保存热图
png("results/heatmap_top50.png", width = 10, height = 12, units = "in", res = 300)
print(heatmap_plot)
dev.off()

9.4 MA图与散点图

# MA图
ma_plot <- plotMA(res, main = "MA Plot", ylim = c(-5, 5))

# 使用ggplot2绘制更美观的MA图
res_df$baseMean_log <- log2(res_df$baseMean + 1)

ma_gg <- ggplot(res_df, aes(x = baseMean_log, y = log2FoldChange)) +
    geom_point(aes(color = padj < 0.05), alpha = 0.5, size = 1) +
    scale_color_manual(values = c("TRUE" = "red", "FALSE" = "gray")) +
    theme_bw() +
    xlab("Log2 Mean Expression") +
    ylab("Log2 Fold Change") +
    ggtitle("MA Plot") +
    geom_hline(yintercept = 0, linetype = "dashed", color = "blue")

print(ma_gg)
ggsave("results/ma_plot.png", ma_gg, width = 8, height = 6, dpi = 300)

# 样本间相关性散点图
plot(countData[,1], countData[,2], 
     log = "xy", 
     xlab = colnames(countData)[1], 
     ylab = colnames(countData)[2],
     main = "Sample Correlation Scatter Plot",
     pch = 20, 
     col = rgb(0,0,0,0.3))

# 样本相关性热图
cor_matrix <- cor(countData, method = "spearman")
pheatmap(cor_matrix, 
         main = "Sample Correlation Heatmap",
         display_numbers = TRUE)

十、完整脚本汇总

10.1 Linux分析流程脚本

#!/bin/bash
# RNA-Seq分析完整流程脚本
# 使用方法: bash rnaseq_pipeline.sh

# 设置变量
PROJECT_DIR=~/rnaseq_project
GENOME_URL="ftp://ftp.ensembl.org/pub/release-110/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz"
GTF_URL="ftp://ftp.ensembl.org/pub/release-110/gtf/homo_sapiens/Homo_sapiens.GRCh38.110.gtf.gz"
THREADS=8

# 创建目录结构
echo "=== 创建项目目录 ==="
mkdir -p ${PROJECT_DIR}/{raw_data,clean_data,ref_genome/{hisat2_index,star_index},aligned,counts,results,scripts,logs,qc_reports}

# 下载参考基因组和注释
echo "=== 下载参考基因组 ==="
cd ${PROJECT_DIR}/ref_genome
if [ ! -f "Homo_sapiens.GRCh38.dna.primary_assembly.fa" ]; then
    wget -c ${GENOME_URL}
    gunzip Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
fi

echo "=== 下载基因注释 ==="
if [ ! -f "Homo_sapiens.GRCh38.110.gtf" ]; then
    wget -c ${GTF_URL}
    gunzip Homo_sapiens.GRCh38.110.gtf.gz
fi

# 构建索引
echo "=== 构建HISAT2索引 ==="
if [ ! -f "hisat2_index/GRCh38.1.ht2" ]; then
    hisat2-build -p ${THREADS} \
        Homo_sapiens.GRCh38.dna.primary_assembly.fa \
        hisat2_index/GRCh38
fi

# 数据质控和比对（循环处理每个样本）
process_sample() {
    local sample=$1
    echo "=== 处理样本: ${sample} ==="
    
    # FastQC质控
    fastqc ${PROJECT_DIR}/raw_data/${sample}_*.fastq.gz \
        -t 2 -o ${PROJECT_DIR}/qc_reports/
    
    # Trimmomatic质控
    trimmomatic PE -threads 4 \
        ${PROJECT_DIR}/raw_data/${sample}_1.fastq.gz \
        ${PROJECT_DIR}/raw_data/${sample}_2.fastq.gz \
        ${PROJECT_DIR}/clean_data/${sample}_1_paired.fq.gz \
        ${PROJECT_DIR}/clean_data/${sample}_1_unpaired.fq.gz \
        ${PROJECT_DIR}/clean_data/${sample}_2_paired.fq.gz \
        ${PROJECT_DIR}/clean_data/${sample}_2_unpaired.fq.gz \
        ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \
        LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36
    
    # HISAT2比对
    hisat2 -p ${THREADS} \
        -x ${PROJECT_DIR}/ref_genome/hisat2_index/GRCh38 \
        -1 ${PROJECT_DIR}/clean_data/${sample}_1_paired.fq.gz \
        -2 ${PROJECT_DIR}/clean_data/${sample}_2_paired.fq.gz \
        --dta \
        --summary-file ${PROJECT_DIR}/logs/${sample}_align_summary.txt \
        | samtools sort -@ 4 -o ${PROJECT_DIR}/aligned/${sample}.sorted.bam -
    
    # 建立索引
    samtools index ${PROJECT_DIR}/aligned/${sample}.sorted.bam
    
    # 去除重复
    picard MarkDuplicates \
        I=${PROJECT_DIR}/aligned/${sample}.sorted.bam \
        O=${PROJECT_DIR}/aligned/${sample}.dedup.bam \
        M=${PROJECT_DIR}/logs/${sample}_dup_metrics.txt \
        REMOVE_DUPLICATES=true \
        CREATE_INDEX=true
}

# 导出函数以便并行处理
export -f process_sample
export PROJECT_DIR THREADS

# 样本列表
SAMPLES=("SRX20066670" "SRX20066671" "SRX20066677" "SRX20066676" "SRX20066680" "SRX20066684")

# 顺序处理（或使用parallel并行处理）
for sample in "${SAMPLES[@]}"; do
    process_sample ${sample}
done

# MultiQC汇总
echo "=== 生成MultiQC报告 ==="
multiqc ${PROJECT_DIR}/qc_reports/ -o ${PROJECT_DIR}/qc_reports/multiqc_report

# FeatureCounts计数
echo "=== 基因表达计数 ==="
featureCounts -T ${THREADS} \
    -p -t exon -g gene_id \
    -a ${PROJECT_DIR}/ref_genome/Homo_sapiens.GRCh38.110.gtf \
    -o ${PROJECT_DIR}/counts/gene_counts.txt \
    ${PROJECT_DIR}/aligned/*.dedup.bam

# 提取表达矩阵
cut -f1,7- ${PROJECT_DIR}/counts/gene_counts.txt | \
    grep -v "^#" | \
    sed 's/.dedup.bam//g' > ${PROJECT_DIR}/counts/count_matrix.txt

echo "=== 分析完成 ==="
echo "结果文件位于: ${PROJECT_DIR}/results/"
echo "请运行R脚本进行差异分析和可视化"

10.2 R分析完整代码

#!/usr/bin/env Rscript
# RNA-Seq差异分析完整R脚本
# 使用方法: Rscript deseq2_analysis.R

# 加载包 ----
suppressPackageStartupMessages({
    library(DESeq2)
    library(ggplot2)
    library(pheatmap)
    library(RColorBrewer)
    library(dplyr)
    library(readr)
    library(EnhancedVolcano)
    library(org.Hs.eg.db)
})

# 设置工作目录 ----
setwd("~/rnaseq_project")

# 创建结果目录 ----
dir.create("results/figures", recursive = TRUE, showWarnings = FALSE)

# 读取数据 ----
message("读取计数矩阵...")
countData <- read.table("counts/count_matrix.txt", 
                        header = TRUE, 
                        row.names = 1, 
                        sep = "\t")

coldata <- read.table("counts/sample_info.txt", 
                      header = TRUE, 
                      row.names = 1, 
                      sep = "\t")

# 确保样本顺序一致 ----
coldata <- coldata[match(colnames(countData), rownames(coldata)), , drop = FALSE]
coldata$group <- factor(coldata$group)

# 创建DESeq2对象 ----
message("创建DESeq2对象...")
dds <- DESeqDataSetFromMatrix(
    countData = countData,
    colData = coldata,
    design = ~ group
)

# 预过滤 ----
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]
message(paste("过滤后基因数:", nrow(dds)))

# 设置参考水平 ----
dds$group <- relevel(dds$group, ref = "control")

# 运行DESeq2 ----
message("运行DESeq2分析...")
dds <- DESeq(dds)

# 获取结果 ----
res <- results(dds, contrast = c("group", "treatment", "control"))
res$gene_symbol <- mapIds(org.Hs.eg.db, 
                          keys = rownames(res), 
                          column = "SYMBOL", 
                          keytype = "ENSEMBL",
                          multiVals = "first")

# 保存结果 ----
write.csv(as.data.frame(res[order(res$padj),]), 
          file = "results/deseq2_results.csv")

sig_genes <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
write.csv(as.data.frame(sig_genes), 
          file = "results/significant_genes.csv")

message(paste("显著差异基因数:", nrow(sig_genes)))

# 标准化 ----
vsd <- vst(dds, blind = FALSE)
write.csv(assay(vsd), file = "results/normalized_counts.csv")

# 可视化 ----
message("生成可视化图表...")

# 1. PCA图
pcaData <- plotPCA(vsd, intgroup = "group", returnData = TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))

pca_plot <- ggplot(pcaData, aes(x = PC1, y = PC2, color = group)) +
    geom_point(size = 3) +
    xlab(paste0("PC1: ", percentVar[1], "% variance")) +
    ylab(paste0("PC2: ", percentVar[2], "% variance")) +
    ggtitle("PCA Plot") +
    theme_bw()

ggsave("results/figures/pca_plot.png", pca_plot, width = 8, height = 6, dpi = 300)

# 2. 火山图
volcano_plot <- EnhancedVolcano(
    res,
    lab = res$gene_symbol,
    x = 'log2FoldChange',
    y = 'padj',
    title = 'Treatment vs Control',
    pCutoff = 0.05,
    FCcutoff = 1
)

ggsave("results/figures/volcano_plot.png", volcano_plot, width = 10, height = 8, dpi = 300)

# 3. 热图
top_genes <- head(order(res$padj), 50)
mat <- assay(vsd)[top_genes, ]
row_names <- res$gene_symbol[top_genes]
row_names[is.na(row_names)] <- rownames(mat)[is.na(row_names)]
rownames(mat) <- row_names
mat_scaled <- t(scale(t(mat)))

annotation_col <- data.frame(Group = coldata$group)
rownames(annotation_col) <- colnames(mat)

png("results/figures/heatmap_top50.png", width = 10, height = 12, units = "in", res = 300)
pheatmap(mat_scaled,
         annotation_col = annotation_col,
         cluster_rows = TRUE,
         cluster_cols = TRUE,
         main = "Top 50 DEGs Heatmap",
         color = colorRampPalette(rev(brewer.pal(9, "RdBu")))(100))
dev.off()

# 4. MA图
png("results/figures/ma_plot.png", width = 8, height = 6, units = "in", res = 300)
plotMA(res, main = "MA Plot", ylim = c(-5, 5))
dev.off()

message("分析完成！结果保存在 results/ 目录")