还不知道TCGA数据如何下载处理?快跟上小果步伐去学习!!
想要分析肿瘤数据的小伙伴,都少不了TCGA数据的使用,这个数据库呢全名是癌症基因组图谱计划。.作为目前最大的癌症基因信息数据库,包含了许多肿瘤的数据,所以小伙伴都很想用TCGA数据去做分析,但是刚入手不知道从哪里下载,小果这里教大家如何去下载并使用R语言去处理数据。
这里小果使用UCSC xena数据库 这个数据库包含又TCGA数据库肿瘤信息,并且数据已经整理好,最重要的是他是开源的!!来看一下:#xena官网
#https://xenabrowser.net/datapages/
这是他的主页面 这里小果使用宫颈癌示例数据去演示,宫颈癌在TCGA数据库中是CESC我们找到这个 点击进入
这里面就是他的数据了,我们主要使用的呢就是counts数据还以FPKMS数据,这里就展示Counts数据的处理,FPKM数据和counts数据处理形式是一样的,
我们点击进入
这里就会有他的信息,我们下载就会得到
这样一个表达文件。下面我们就去处理这个数据吧!
公众号回复111获取代码,代码编号:231107
#install.packages("tidyverse")
library(tidyverse)#每次重新打开R都要library一下
##文件的读取
#读取tsv文件
counts1 = read.table(file = 'TCGA-CESC.htseq_counts.tsv', sep = 't', header = TRUE)
得到一下的矩阵
我们需要对矩阵进行处理,将第一列放在文本框中
rownames(counts1) <- counts1[,1] #Alt <-
#x <- counts1[,1:3]
counts1 = counts1[,-1]
接下来就是数据处理,我们要得到肿瘤的样本还以正常的样本数据,含有01A的是肿瘤样本,11A的为正常,因此我们使用下面table函数去寻找并计算,
#table函数
table(substr(colnames(counts1),14,16))
#c("01A","11A")
#%in%符号用于判断是否属于
counts1 <- counts1[,substr(colnames(counts1),14,16)%in% c("01A","11A")]
table(substr(colnames(counts1),14,16))
#保留列名前15位
rownames(counts1) <- substr(rownames(counts1),1,15)
ceiling(1.2)
ceiling(3.8)
counts <- ceiling(2^(counts1)-1)
这样一来列的名称中含有小数点的数字没有了,这样方便我们去对照寻找基因名
我们将
write.table(counts,"counts.txt",sep = "t",row.names = T,col.names = NA,quote = F)
那么我们有了数据,肯定要做差异表达,下面小果带大家去学习对上述矩阵进行差异表达分析。
counts <- read.table("counts.txt",sep = "t",row.names = 1,check.names = F,stringsAsFactors = F,header = T)
Ginfo_0 <- read.table("gene_length_Table.txt",sep = "t",check.names = F,stringsAsFactors = F,header = T,row.names = 1)
就是上述注释文件,包含我们想要的基因信息
我们只要RNA信息,并且找到Gene名替换矩阵前面的第一列信息:
Ginfo <- Ginfo_0[which(Ginfo_0$genetype == "protein_coding"),] #只要编码RNA
#美元符号代表提取列
#取行名交集intersect
comgene <- intersect(rownames(counts),rownames(Ginfo))
counts <- counts[comgene,]
#class(counts)#判断数据类型
#class(comgene)
Ginfo <- Ginfo[comgene,]
a <- rownames(counts)
b <- rownames(Ginfo)
identical(a,b)#判断是否相同
counts$Gene <- as.character(Ginfo$genename) #新增Gene Symbol
counts <- counts[!duplicated(counts$Gene),] #去重复
rownames(counts) <- counts$Gene #将行名变为Gene Symbol
counts <- counts[,-ncol(counts)] #去除最后一列
write.table(counts, file = "CESC_counts_mRNA_all.txt",sep = "t",row.names = T,col.names = NA,quote = F)
这里我们#保存癌症患者的counts
tumor <- colnames(counts)[substr(colnames(counts),14,16) == "01A"]
counts_01A <- counts[,tumor]#并且找到只要肿瘤的矩阵,为了我们后续差异分析用
write.table(counts_01A, file = "CESC_counts_mRNA_01A.txt",sep = "t",row.names = T,col.names = NA,quote = F)
下面开始我们的差异分析吧!
#差异分析
library(tidyverse)
#安装BiocManager
#if(!require(DESeq2))BiocManager::install('DESeq2')
library(DESeq2)#这里我们用这个包DESeq2,做差异分析的包
后面就跟着小果代码跑跑跑就行了
counts = counts[apply(counts, 1, function(x) sum(x > 1) > 32), ]
conditions=data.frame(sample=colnames(counts),
group=factor(ifelse(substr(colnames(counts),14,16) == "01A","T","N"),levels = c("N","T"))) %>%
column_to_rownames("sample")
dds <- DESeqDataSetFromMatrix(
countData = counts,
colData = conditions,
design = ~ group)
dds <- DESeq(dds)#这一步有点慢,小伙伴多多等待不要动
resultsNames(dds)
res <- results(dds)
save(res,file = "CESC_DEG.rda")#一定要保存!
res_deseq2 <- as.data.frame(res)%>%
arrange(padj) %>%
dplyr::filter(abs(log2FoldChange) > 4, padj < 0.05)#根据自己需要
这样结果就跑完了,我们得到基因数1194,小伙伴觉得差异基因太少或者太多,自行去修改
##dplyr::filter(abs(log2FoldChange) > 4, padj < 0.05)#根据自己需要这里调整阈值的选择。一定在合理范围之内
以上就是我们对TCGA数据的处理了,小伙伴有没有学会呢,多多理解代码的意义,不要跑错了
往期推荐