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

大家好，欢迎来到IT知识分享网。

一、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__.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)

输出结果如下：

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

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

下次见~