书接上文：小白必会之：2024最新TCGA数据下载与处理手把手教学
获取表达矩阵后接上DESeq2差异分析流程
Step1. 数据预处理
废话不多说直接上代码：
首先读取metadata,获取样本信息
setwd("D:/work/bininfo/SXFX955/TCGA/TCGA-COAD/")#设置工作路径,即下载的数据的路径
rm(list = ls())
library(dplyr)
# 1. 读取metadata,获取样本信息
metadat <- read.table(file = "gdc_sample_sheet.2024-08-11.tsv",sep = '\t',header = T)
count <- read.csv('TCGA_COAD_count.csv',row.names = 1)
由于R语言的特性,'-'短横杠在读取的时候会变成‘.’，为了和metadat中SampleID对应,我们把'-'改成下划线’_‘,并重新保存,否则每次读取都要重新命名
colnames(count) <- gsub('\\.','_',colnames(count))
write.csv(count,file = 'TCGA_COAD_count.csv')
count <- read.csv('TCGA_COAD_count.csv',row.names = 1)
# metadat也改一下并重新生成一个meta文件
meta <- metadat[,c("Case.ID","Sample.ID","Sample.Type")]
meta$Case.ID <- gsub('-','_',meta$Case.ID)
meta$Sample.ID <- gsub('-','_',meta$Sample.ID)
查看一下分组信息
table(meta$Sample.Type)
这里我们把转移和复发的去掉(就两例)
rmid <- meta[(meta$Sample.Type=='Metastatic' | meta$Sample.Type=='Recurrent Tumor'),'Sample.ID']
count <- count[,-match(rmid,colnames(count))]
meta <- meta[(meta$Sample.Type!='Metastatic' & meta$Sample.Type!='Recurrent Tumor'),]
table(meta$Sample.Type)
名字太长了我们改一下
meta$group <- plyr::mapvalues(meta$Sample.Type,
from = c('Primary Tumor', 'Solid Tissue Normal'),
to = c('Tumor','Normal'))
检查一下count矩阵的列名是否都在meta的sample
identical(sort(meta$Sample.ID),sort(colnames(count))) #顺序不一样我们sort一下
小插曲
这里发现居然meta样本信息中sampleid和count列名对不上，需要仔细检查一下哪里不对
首先match函数看下有多少个id对不上
match(colnames(count),meta$Sample.ID) %>% is.na() %>% table()
查看是哪10个文件
colnames(count)[which(is.na(match(colnames(count),meta$Sample.ID)))]
拿出一个进行匹配
grep('TCGA_A6_5656_01A',meta$Sample.ID,value = T)
发现居然有重复的id,继续在count矩阵中查看是怎么回事
head(count[,grep('TCGA_A6_2674_01A',colnames(count))])
回到原始metadat中查看,注意把下划线'_'改回短横杠'-'
metadat[grep('TCGA-A6-2674-01A',metadat$Sample.ID),]
可以发现确实是不同的文件,但是Case.ID, Sample.ID, Sample.Type是一样的。
这里就需要提及到TCGA样本名重复的问题了，其实也可以理解,这是因为同一个病人会有多个样本,关于这个问题网上有博主写过推文：
想了解更多的可以参照https://www.jianshu.com/p/f07721d3bd57和https://www.jianshu.com/p/74c36463a97e
这里就不提及更多了
我们简单的去除这些重复的文件（更为妥善的解决办法可以参考上述链接，本篇主题是TCGA数据差异分析），代码如下
# 简单去除重复
count <- count[,grep('_1$',colnames(count),invert = T)]
meta <- meta[!duplicated(meta$Sample.ID),]
rownames(meta) <- meta$Sample.ID
小插曲结束
通过查看count矩阵得知基因是Ensembl ID,且后面有小数点(代表的是该Ensembl ID的版本)
因此我们需要将小数点去除
首先将ENSEMBL ID后的version去掉并赋值给新的一列命名为ENSEMBL列
# 首先将ENSEMBL ID后的version去掉并赋值给新的一列
count$ENSEMBL <- sapply(rownames(count),function(x){strsplit(x,'\\.')[[1]][1]})
然后判断去掉小数点后是否有重复ENSEMBL ID
count$ENSEMBL[which(table(count$ENSEMBL)>1)]
count[which(count$ENSEMBL=='ENSG00000124333'),c(1:3)]
rowSums(count[which(count$ENSEMBL=='ENSG00000124333'),c(1:512)])
还挺多，但我们发现重复的这个ID，所有样本的count值为0，我们可以简单去除这些ID
# 去除所有为0的ID
count <- count[rowSums(count[,c(1:512)])>0,]
# 再次检测是否有重复ID
count$ENSEMBL[which(table(count$ENSEMBL)>1)]
最后将count的列名设置为ENSEMBL列的ID
rownames(count) <- count$ENSEMBL
# 去除ENSEMBL列
count <- count[,-513]
去除后数据预处理就完成了，接下来进行第二步：DESeq2差异分析流程
Step2. DESeq2差异分析流程
代码如下
library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = count, colData = meta[colnames(count),], design = ~ group)
dds
# 过滤低表达基因
keep <- rowSums(counts(dds)>10) >= 10 #至少10个样本表达值大于10
table(keep)
dds <- dds[keep,]
# 注意Tumor在前，Normal在后,即为 Tumor/Normal
res <- results(dds, contrast=c("group","Tumor","Normal"))
resDE_res <- as.data.frame(res)
DE_res$label <- ifelse(DE_res$padj<0.05 & abs(DE_res$log2FoldChange)>= 1,
ifelse(DE_res$log2FoldChange >= 0,'up','down'),'no')
table(DE_res$label)
# 保存结果到DEresult文件夹
write.csv(DE_res,file = './DEresult/DE_res.csv')
画火山图进行可视化
# 火山图
library(ggplot2)
DE_res$label <- factor(DE_res$label,levels = c("up","down","no"))
up_label <- paste0("up:",table(DE_res$label)["up"])
down_label <- paste0("down:",table(DE_res$label)["down"])
no_label <- paste0("no_sig:",table(DE_res$label)["no"])#padj pvalue
ggplot(DE_res)+geom_point(aes(log2FoldChange,-log10(padj),color=label))+ scale_color_manual(values = c("up"="#d62728","down"="#1f77b4","no"="grey"), labels=c(up_label,down_label,no_label)) + theme_bw(base_size = 14)
ggsave('./DEresult/volcanoplot.png',width = 5,height = 4)
画PCA图并可视化
# 首先使用DESeq2的VST转换用于标准化
vsd <- vst(dds, blind=FALSE)
# DESeq2的plotPCA函数直接出PCA图
pcaplot <- DESeq2::plotPCA(vsd, intgroup = "group", ntop = 500, returnData = FALSE)
# 更改一下点的大小
pcaplot[["layers"]][[1]][["aes_params"]][["size"]] = 1
pcaplot +
# 自定义点的颜色
scale_color_manual(values = c("Tumor"="#d62728","Normal"="#1f77b4")) +
# 添加置信椭圆
stat_ellipse(aes(fill=group),type="norm",geom="polygon",alpha=0.2,color=NA) +
# 自定义椭圆颜色
scale_fill_manual(values = c("Tumor"="#d62728","Normal"="#1f77b4")) +
theme_bw(base_size = 14)
ggsave('./DEresult/PCAplot.png',width = 5,height = 4)
Step3 基因symbol转换
由于count矩阵行名为Ensembl ID，我们需要将Ensembl ID转换为大家熟知的Gene Symbol，两种方法
方法一：使用clusterProfiler包的bitr函数
# 将Ensembl ID转换为Gene Symbol ####
# 方法一
library(clusterProfiler)
library(org.Hs.eg.db)
DE_res$ENSEMBL <- rownames(DE_res)
IDdf <- bitr(DE_res$ENSEMBL,fromType = "ENSEMBL",
toType = "SYMBOL",OrgDb = org.Hs.eg.db, drop = F)
head(IDdf)
# 合并IDdf和DE_res
DE_res_symbol <- merge(DE_res,IDdf,by='ENSEMBL')
up_genes <- na.omit(subset(DE_res_symbol,label=='up')$SYMBOL)
dw_genes <- na.omit(subset(DE_res_symbol,label=='down')$SYMBOL)
any(duplicated(up_genes))
方法二：对于ensembl id推荐使用官方的biomaRt包
# 方法二
library("biomaRt")
# 用useMart函数链接到人类的数据库
ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
# 除此以外还有很多：使用 listDatasets(ensembl) 查看
attributes <- listAttributes(ensembl)
IDdf2 <- getBM(attributes = c("ensembl_gene_id","hgnc_symbol"),
filters = "ensembl_gene_id",
values = count$ENSEMBL,
mart = ensembl)
head(IDdf2)
# 合并IDdf和DE_res
DE_res_symbol2 <- merge(DE_res,IDdf2,by='ENSEMBL')
可能有小伙伴会问为何不在做差异分析之前进行转换，主要是因为基因ID的转换会出现一对多和多对一的情况。
如果需要转换后再做差异分析，我们可以采取不同的策略，比如一个Ensembl ID对应多个Gene Symbol时，可以简单粗暴随机取一个，而多个Ensembl ID对应一个Gene Symbol时，可以‘加和’或者‘取平均’。
至于在做差异分析之前是否要进行转换，仁者见仁智者见智。不过国外专注于生物信息学领域的问答平台：Biostars上早就有人讨论过这个问题：https://www.biostars.org/p/352492/
高赞回答是这样说的：
简而言之就是不同的Ensembl ID代表了不同的genomic loci，理论上说它们的功能可能是有差异的（否则取多个id名字干嘛对吧），高赞答主认为既然有差异不能简单认为它们是一样的，对数据进行加和取平均。当然话又说回来，既然是同一个基因，那么上面说的差异可能也没有那么大，对于你的结果来说影响其实也很小。所以我的看法是，两个字：‘随你’
over ~
