使用R语言下载TCGA数据库癌症基因表达数据小例子

Posted 2021-04-29 小明的数据分析笔记本

参考资料

• 简书文章 TCGA癌症缩写、癌症中英文对照

• CDSN文章 通过整理TCGA数据，探索某癌症的癌组织和正常组织的差异基因

• TCGA Workflow: Analyze cancer genomics and epigenomics data using Bioconductor packages

• TCGAbiolinks包下载TCGA数据进行表达差异分析-乳腺癌案例

代码

数据下载

BiocManager::install("TCGAbiolinks")
library(TCGAbiolinks)
?GDCquery
x<-getGDCprojects()$project_id
x == "TCGA-COAD"
x == "TCGA-READ"
TCGAbiolinks:::getProjectSummary("TCGA-COAD")
query<-GDCquery(project = "TCGA-COAD",
                legacy=F,
                experimental.strategy ="RNA-Seq",
                data.category = "Transcriptome Profiling",
                data.type="Gene Expression Quantification",
                workflow.type="HTSeq - Counts")
GDCdownload(query)

这里遇到的问题是：所有数据都下载下来了，如何将数据整合成表达量矩阵的形式呢？

更新

上面遇到的问题今天找到了解决办法：可以使用SummarizedExperiment包中的函数assay()函数将表达矩阵提取出来；colData()函数好像是获得一些样本信息

查看TCGAbiolinks包的帮助文档
browseVignettes("TCGAbiolinks")
尝试重复帮助文档中的第3项内容 Downloading and preparing files for analysis
这篇帮助文档中用到的例子是TCGA-GBM——Glioblastoma multiforme——多形成性胶质细胞瘤

获取基因表达矩阵

query_3 <- GDCquery(project = "TCGA-GBM",
                    data.category = "Gene expression",
                    data.type = "Gene expression quantification",
                    platform = "Illumina HiSeq", 
                    file.type  = "results",
                    experimental.strategy = "RNA-Seq",
                    barcode = barcode,
                    legacy = TRUE)
GDCdownload(query_3)
dataPrep1<-GDCprepare(query=query_3)
library(SummarizedExperiment)
expdat<-assay(dataPrep1)
dim(expdat)
head(expdat)

备注：file.type 的另外一个参数是 normalized_results

DESeq2基因差异表达分析

library(stringr)
group_list<-ifelse(str_sub(colnames(expdat),14,15) == "01","tumor","normal")
group_list
coldata<-data.frame(row.names = colnames(expdat),
                    condition=group_list)
coldata
library(DESeq2)
expdat<-round(expdat,0)
dds<-DESeqDataSetFromMatrix(countData = expdat,
                            colData = coldata,
                            design = ~ condition)
keep<-rowSums(counts(dds))>=10
dds<-dds[keep,]
dds<-DESeq(dds)
res<-results(dds,contrast = c("condition","tumor","normal"))
resOrdered <- res[order(res$pvalue),]
DEG <- as.data.frame(resOrdered)
DEG <- na.omit(DEG)
logFC_cutoff<-with(DEG,mean(abs(log2FoldChange))+2*sd(abs(log2FoldChange)))
logFC_cutoff
DEG$change<-as.factor(ifelse(DEG$pvalue<0.05&abs(DEG$log2FoldChange)>logFC_cutoff,
                             ifelse(DEG$log2FoldChange>logFC_cutoff,"UP","DOWN"),
                             "NOT"))
library(ggplot2)
this_title <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
                     '\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,
                     '\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',]))
ggplot(data=DEG,aes(x=log2FoldChange,y=-log10(pvalue),color=change))+
  geom_point(alpha=0.4,size=1.75)+
  labs(x="log2 fold change")+ ylab("-log10 pvalue")+
  ggtitle(this_title)+theme_bw(base_size = 10)+
  theme(plot.title = element_text(size=15,hjust=0.5))+
  scale_color_manual(values=c('blue','black','red'))

Rplot.png

使用R语言包 clusterProfiler 差异表达基因的GO富集分析

先看一下这个包的帮助文档

browseVignettes("clusterProfiler")
help(package="clusterProfiler")
library(clusterProfiler)
data(geneList,package="DOSE")
names(geneList)[1:100]

准备自己的输入数据

df<-DEG[which(DEG$change=="UP"),]
dim(df)
rownames(df)
gene_list<-str_split(rownames(df),"\\|")
gene_list
class(gene_list)
gene_list<-as.data.frame(matrix(unlist(a),ncol=2,byrow = T))
head(gene_list)
yy<-enrichGO(gene=gene_list$V2,
             OrgDb='org.Hs.eg.db',
             ont="BP",
             pvalueCutoff = 0.01)
dotplot(yy)

Rplot01.png

有时间试着用今天分析的数据准备一下GOplot包的输入数据。

公众号二维码.jpg

中国加油！武汉加油！