学习笔记 — TCGA 差异表达分析及可视化

一、TCGA数据下载 （LIHC为例）

数据下载的方式和之前学习的临床数据的下载类似，先进入官网 https://portal.gdc.cancer.gov/

新版TCGA数据库下载流程：

Cohort Builder → Program (TCGA)、 Project (LIHC) → 点击 Repository → 侧边栏筛选：Experimental Strategy (RNA-Seq)  → Data Category (transcriptome profiling) → Data Type (Gene Expression Quantification)

可以看到筛选出来的数据包括了424个文件，一共371例患者数据， 全部加入Cart，再下载 Sample Sheet 和 Cart 两个文件。

二、TCGA数据数据整理

在RNA-Seq的基础理论篇中学习了RNA-Seq的基本原理以及数据清洗的几种常用格式，如：CPM、TPM、RPKM/FPKM。接下来就利用TCGA数据库的LIHC数据进行数据整理，获得我们进行差异分析需要的counts、TPM等数据。

### 解压数据，创建存储文件夹
setwd("TCGA-LIHC") # 设置工作路径
dir.create('RawMatrix/') # 新建文件夹存储下载的原始数据
tar_file <- "./gdc_download_20240606_135942.082516.tar.gz" 
extract_dir <- "./RawMatrix"
untar(tar_file, exdir = extract_dir) # 导入tar.gz，并解压文件
dir.create('RawData/') # 新建文件夹存储count/TPM/差异表达矩阵等txt格式
dir.create('RawData/csv/') # 新建文件夹存储csv格式的矩阵

### 数据整理
library(data.table)
library(dplyr)
sample_sheet <- fread("./gdc_sample_sheet.2024-06-06.tsv") # 读取样本信息
sample_sheet$Barcode <- substr(sample_sheet$`Sample ID`,1,15) # 取ID前15字符作为barcode
sample_sheet1 <- sample_sheet %>% filter(!duplicated(sample_sheet$Barcode)) # 去重
sample_sheet2 <- sample_sheet1 %>% filter(grepl("01$|11$|06$",sample_sheet1$Barcode)) # Barcode的最后两位：01表示肿瘤样本，11表示正常样本，06表示转移样本

TCGA_LIHC_Exp <- fread("./RawMatrix/00fb4b52-e6a4-4ad9-bbed-584e25851aca/ba056a7d-5370-4fe9-af1e-2e3de42e205f.rna_seq.augmented_star_gene_counts.tsv") # 任意读取一个文件
TCGA_LIHC_Exp <- TCGA_LIHC_Exp[!1:4,c("gene_id","gene_name","gene_type")] # 创建包含"gene_id","gene_name","gene_type"的数据框，用于合并表达数据

### 将所有样本合并成一个数据框
for (i in 1:nrow(sample_sheet2)) {
  
  folder_name <- sample_sheet2$`File ID`[i]
  file_name <- sample_sheet2$`File Name`[i]
  sample_name <- sample_sheet2$Barcode[i]
  
  data1 <- fread(paste0("./RawMatrix/",folder_name,"/",file_name))
  #unstranded代表count值；如果要保存TPM，则改为tpm_unstranded
  data2 <- data1[!1:4,c("gene_id","gene_name","gene_type","unstranded")] 
  colnames(data2)[4] <- sample_name
  
  TCGA_LIHC_Exp <- inner_join(TCGA_LIHC_Exp,data2)
  
}

### 根据需要的表达比例筛选满足条件的基因
zero_percentage <- rowMeans(TCGA_LIHC_Exp[, 4:ncol(TCGA_LIHC_Exp)] == 0) 
TCGA_LIHC_Exp1 <- TCGA_LIHC_Exp[zero_percentage < 0.6, ] # 筛选出表达超过60%的基因

install.packages("BiocManager")
library(BiocManager)
BiocManager::install("edgeR")
library(edgeR)
TCGA_LIHC_Exp1 = avereps(TCGA_LIHC_Exp[,-c(1:3)],ID = TCGA_LIHC_Exp$gene_name) # 对重复基因名取平均表达量，并将基因名作为行名
TCGA_LIHC_Exp1 <- TCGA_LIHC_Exp1[rowMeans(TCGA_LIHC_Exp1)>100,] # 根据需要去除低表达基因，这里设置的平均表达量100为阈值

### 创建样本分组
library(stringr)
tumor <- colnames(TCGA_LIHC_Exp1)[substr(colnames(TCGA_LIHC_Exp1),14,15) == "01"]
normal <- colnames(TCGA_LIHC_Exp1)[substr(colnames(TCGA_LIHC_Exp1),14,15) == "11"]
tumor_sample <- TCGA_LIHC_Exp1[,tumor]
normal_sample <- TCGA_LIHC_Exp1[,normal]
exprSet_by_group <- cbind(tumor_sample,normal_sample)
gene_name <- rownames(exprSet_by_group)
exprSet <- cbind(gene_name, exprSet_by_group)  # 将gene_name列设置为数据框的行名，合并后又添加一列基因名

### 存储counts和TPM数据
fwrite(exprSet,"./RawData/TCGA_LIHC_Count.txt") # txt格式
write.csv(exprSet, "./RawData/csv/TCGA_LIHC_Count.csv", row.names = FALSE) # csv格式
fwrite(exprSet,"./RawData/TCGA_LIHC_Tpm.txt") # txt格式
write.csv(exprSet, "./RawData/csv/TCGA_LIHC_Tpm.csv", row.names = FALSE) # csv格式

三、差异表达分析

这里分享用edgeR包进行差异分析的操作：

### 差异表达分析- edgeR包

library(edgeR)
group_list <- c(rep('tumor',ncol(tumor_sample)),rep('normal',ncol(normal_sample)))
group_list <- factor(group_list)
design <- model.matrix(~0+group_list)
rownames(design) <- colnames(exprSet_by_group)
colnames(design) <- levels(group_list) # 创建分组

# 创建一个DGEList对象，用于存储差异表达分析所需的数据
DGElist <- DGEList(counts = exprSet_by_group, group = group_list)

# 使用cpm值对低表达量的基因进行过滤。
keep_gene <- rowSums(cpm(DGElist) > 1 ) >= 2
DGElist <- DGElist[keep_gene, keep.lib.sizes = FALSE] # 保留符合条件的基因

# 校正测序深度
DGElist <- calcNormFactors( DGElist ) 

# 估算离散度，该步骤会对DGElist进行添加或更新，一般看CommonDisp就行
DGElist <- estimateGLMCommonDisp(DGElist, design) # 共同离散度
DGElist <- estimateGLMTrendedDisp(DGElist, design) # 趋势离散度
DGElist <- estimateGLMTagwiseDisp(DGElist, design) # 基因特异的离散度

# 拟合广义线性模型
fit <- glmFit(DGElist, design) # 偏离分布模型的基因即差异表达基因DEGs
results <- glmLRT(fit, contrast = c(-1, 1)) # 似然比检验，contrast = c(-1, 1)即对分组的两个条件检验，这里是tumor和normal

# 提取差异表达的top基因
LIHC_nrDEG_edgeR <- topTags(results, n = nrow(DGElist)) # n = nrow(DGElist)即全部保存
LIHC_nrDEG_edgeR <- as.data.frame(LIHC_nrDEG_edgeR)

# 保存原始差异基因矩阵文件
fwrite(LIHC_nrDEG_edgeR,"./RawData/LIHC_nrDEG_edgeR.txt", row.names = TRUE) # csv格式
write.csv(LIHC_nrDEG_edgeR, "./RawData/csv/LIHC_nrDEG_edgeR.csv", row.names = T) # csv格式

# 筛选差异基因
library(data.table)
library(dplyr)
LIHC_Match_DEG <- LIHC_nrDEG_edgeR # 或者从输出的文件里读取
LIHC_Match_DEG$log10FDR <- -log10(LIHC_Match_DEG$FDR)
colnames(LIHC_Match_DEG)[1] <- "gene_name"
LIHC_Match_DEG <- LIHC_Match_DEG %>% 
  mutate(DEG = case_when(logFC > 2 & FDR < 0.05 ~ "Up",
                         abs(logFC) < 2 | FDR > 0.05 ~ "None",
                         logFC < -2 & FDR < 0.05 ~ "Down")) # 打标签：logFC > 2 & FDR < 0.05：上调基因，logFC < -2 & FDR < 0.05：下调基因，其它认为无显著差异

# 保存添加标签后的基因
fwrite(LIHC_Match_DEG,"./RawData/LIHC_Match_DEG.txt") # txt格式
write.csv(LIHC_Match_DEG, "./RawData/csv/LIHC_Match_DEG.csv", row.names = F) # csv格式

四、可视化之火山图 (Volcano Plot) 

LIHC_Match_DEG <- fread("./RawData/LIHC_Match_DEG.txt") # 读取打完标签的基因列表

down_gene <- LIHC_Match_DEG[LIHC_Match_DEG$DEG == "Down", ]
up_gene <- LIHC_Match_DEG[LIHC_Match_DEG$DEG == "Up", ]

uptop <- rownames(up_gene)[1:10]  # 上调的前10基因
downtop <- rownames(down_gene)[1:10] # 下调的前10基因

LIHC_Match_DEG$label <- ifelse(LIHC_Match_DEG$gene_name %in% c(uptop,dwtop), LIHC_Match_DEG$gene_name, "") # 后面画图时用来突出显著表达的前10个基因

# 加载需要用到的程序包
library(data.table)
library(ggplot2)
library(ggprism)
library(ggrepel)

# 画图 volcano plot
ggplot(LIHC_Match_DEG, aes(x = logFC, y = log10FDR, colour = DEG)) +
  geom_point(alpha = 0.85, size = 1.5) + # 设置点的透明度和大小
  scale_color_manual(values = c('steelblue', 'gray', 'brown')) + # 调整点的颜色
  xlim(c(-11, 11)) +  # 调整x轴的范围
  geom_vline(xintercept = c(-2, 2), lty = 4, col = "black", lwd = 0.8) + # x轴辅助线
  geom_hline(yintercept = -log10(0.05), lty = 4, col = "black", lwd = 0.8) + # y轴辅助线
  labs(x = "logFC", y = "-log10FDR") + # x、y轴标签
  ggtitle("TCGA LIHC DEG") +  # 图表标题
  theme(plot.title = element_text(hjust = 0.5), legend.position = "right", legend.title = element_blank()) +  # 设置图表标题和图例位置
  geom_label_repel(data = LIHC_Match_DEG, aes(label = label),  # 添加标签
                   size = 3, box.padding = unit(0.5, "lines"),
                   point.padding = unit(0.8, "lines"),
                   segment.color = "black",
                   show.legend = FALSE, max.overlaps = 10000) +  # 标签设置
  theme_prism(border = TRUE)

输出结果如下：

后续将继续学习富集分析，感兴趣的小伙伴麻烦点赞、关注、收藏哦~

欢迎大家在评论区留言讨论，如有不对敬请指出~

下次见~