【手把手教学】进阶版TCGA数据整理与分析
之前小云和小伙伴分享了如何使用TCGA数据库，同时以胰腺癌数据为实例进行了数据下载，那么接下来就要对下载的数据根据需要进行整理
TCGA数据整理
# 数据读取 -------------------------------------------------------------------
首先读取下载数据的记录，这样免得大家二次下载，省去等待时间
load("TCGA_PAAD.Rdata")
# 加载所需要的包 -------------------------------------------------------------
library(dplyr) #数据整理
#只取protein coding，编码蛋白的mRNA进行分析
table(expr_count$gene_type)
mRNA <- expr_count[expr_count$gene_type=="protein_coding",]
# count数据清洗 ---------------------------------------------------------------
mRNA <- mRNA[,-1]
dim(mRNA)#查看一共多少行
mRNA$mean <- rowMeans(mRNA[,2:184])#计算每一行表达量的平均值
mRNA <- arrange(mRNA,desc(mean))#按mean值大小降序排列
mRNA <- mRNA %>% distinct(gene_name, .keep_all = T)#删除重复Gene的行
mRNA <- select(mRNA,-c(mean))#删除之前生成的mean一行
#将肿瘤和非肿瘤样本区分开排序
rownames(mRNA) <- mRNA$gene_name
mRNA <- mRNA[,-1]
table(substr(colnames(mRNA),14,16))
那么-01A和-11A，这里的字母A又是什么含义呢？
从这个示意图上可以看到，有时候我们可以从一个病例身上取多个样本，不论是肿瘤样本，还是癌旁正常对照，然后存放在不同的管子里面。这里的A，B，C就表示样本的顺序。官方文档的解释如下。
我们可以看到后缀中数字与样本类型的对应关系，如下图
解释完楼，来我们继续~
Tumor <- grep('01A',colnames(mRNA))
Tumor #肿瘤样本所处位置
Tumor_mRNA <- mRNA[,Tumor]
Normal <- grep('11A',colnames(mRNA))
Normal #正常样本所处位置
Normal_mRNA <- mRNA[,Normal]
mRNA_count <- cbind(Normal_mRNA,Tumor_mRNA)
write.csv(mRNA_count,file = "mRNA_count.csv")
# TPM数据清洗 -----------------------------------------------------------------
mRNA_TPM <- expr_TPM[expr_TPM$gene_type=="protein_coding",]
mRNA_TPM <- mRNA_TPM[,-1]
dim(mRNA_TPM)#查看一共多少行
mRNA_TPM <- mRNA_TPM %>%
mutate(mean=rowMeans(.[,-1])) %>%
arrange(desc(mean)) %>%
distinct(gene_name,.keep_all = T) %>%
select(-mean)
#将肿瘤样本和非肿瘤样本分开排序
rownames(mRNA_TPM) <- mRNA_TPM$gene_name
mRNA_TPM <- mRNA_TPM[,-1]
table(substr(colnames(mRNA_TPM),14,16))
Tumor <- grep('01A',colnames(mRNA_TPM))
Tumor #肿瘤样本所处位置
Tumor_mRNA_TPM <- mRNA_TPM[,Tumor]
Normal <- grep('11A',colnames(mRNA_TPM))
Normal #正常样本所处位置
Normal_mRNA_TPM <- mRNA_TPM[,Normal]
mRNA_TPM_new <- cbind(Normal_mRNA_TPM,Tumor_mRNA_TPM)
# 保存文件用于下次的分析 -------------------------------------------------------------
save(mRNA_count,mRNA_TPM_new,file = "count_and_TPM.Rdata")
DESeq2差异分析
数据整理完成后，一般我们会获得一个原始 read count表达矩阵，其中行是基因，列是样品。常用的差异分析工具包括limma，edgeR，DESeq2。小图就以R包DESeq2的差异表达基因分析方法为例做个简单演示，因为小图今天下载的是count数据。
简单地说，DESeq2将对原始reads count进行建模，使用标准化因子（scale factor）来解释库深度的差异。然后DESeq2估计基因的离散度，并缩小这些估计值以生成更准确的离散度估计，从而对reads count进行建模。最后，DESeq2拟合负二项分布的模型，并使用Wald检验或似然比检验进行假设检验。
来，马上开始实操~
setwd("D:/ DESEq2差异表达")
# 清空变量 --------------------------------------------------------------------
rm(list=ls(all=TRUE))
options(stringsAsFactors = F)
library(DESeq2)
expr <- expr[rowMeans(expr) > 1,]
rownames(expr)=mRNA_count[,1]
# 创建分组信息 ------------------------------------------------------------------
table(substr(colnames(expr),14,16))#查看肿瘤样本和正常样本数量
Tumor <- grep('01A',colnames(expr))#查看肿瘤样本所处位置
Tumor
group <- factor(c(rep("Normal",times=4),rep("Tumor",times=178)))#创建分组因子变量
Data <- data.frame(row.names = colnames(expr), #创建分组数据框
group = group)
View(Data)
#基因表达矩阵(expr)、样本分组数据框(Data)都已经准备完毕
# 开始进行差异表达分析 --------------------------------------------------------------
#第一步：构建DEseq2对象(dds)
dds <- DESeqDataSetFromMatrix(countData = expr,
colData = Data,
design = ~ group)
#第二步：开始差异分析
dds2 <- DESeq(dds)
res <- results(dds2, contrast=c("group", "Tumor", "Normal"))#肿瘤在前，对照在后
##或者res= results(dds)
res <- res[order(res$pvalue),]
summary(res)
my_result <- as.data.frame(res)#转成容易查看的数据框
my_result <- na.omit(my_result)#删除倍数为0的值
#第三步：保存差异分析的结果
library(dplyr)
my_result$Gene_symbol<-rownames(my_result)
my_result <- my_result %>% dplyr::select('Gene_symbol',colnames(my_result)[1:dim(my_result)[2]-1],everything())
rownames(my_result) <- NULL
write.csv(my_result,file="my_result_DESeq2.csv")
# DEG的筛选 ------------------------------------------------------------------
my_result$regulate <- ifelse(my_result$pvalue > 0.05, "unchanged",
ifelse(my_result$log2FoldChange > 0, "up-regulated",
ifelse(my_result$log2FoldChange < 0, "down-regulated", "unchanged")))
table(my_result$regulate)
#可以把上调基因和下调基因取出放在一块
DEG_deseq2 <-subset(my_result, pvalue < 0.05 & abs(log2FoldChange) > 0)
upgene <- DEG_deseq2[DEG_deseq2$regulate=='up-regulated',]
downgene <- DEG_deseq2[DEG_deseq2$regulate=='down-regulated',]
write.csv(DEG_deseq2,file= "DEG_deseq2_0.05.csv")
好了今天的进阶版TCGA数据整理分享就到这啦，小伙伴们有遇到问题欢迎随时来找小云分享讨论呀！
