1. sci666首页
  2. 实用技巧
  3. 生物信息学

使用TCGAbiolinks批量下载TCGA的表达量数据

TCGA的数据下载,需要去GDC,使用官方工具,
下载的数据每个样本都在单个文件夹中,需要一系列的处理,才能得到表达矩阵。
而数据分析和挖掘就是从表达矩阵开始的。

以往我写贴讲课都是用的手工的方法,目的是为了教会大家批量操作,而批量是我学习编程到现在收货最大的部分。
还要一个工具,集TCGA数据下载处理可视化为一身,就是TCGAbiolinks。
但是因为网络原因,我重来没用过。直到有一天,唐医生帮我彻底解决了网络问题,我才能像一个正常生信工作人员一样生活。
我现在无论是SRA,还是github,都没有网络限制了。

然后我就用上了TCGAbiolinks,最终我的评价就是,

TCGAbiolinks用来下载数据十分方便,用来分析和可视化达不到要求。

现在以下载是RNAseq的数据为例,写一个批量的呈现,下载所有33个癌症的数据变成Rdata格式,并分享。

对包和项目的探索

有哪些项目的数据可以下载呢?

library(TCGAbiolinks)
projects <- getGDCprojects()
使用TCGAbiolinks批量下载TCGA的表达量数据

这个返回的表格列举的可以下载的项目,其中每个项目都还有很多组学可以下载。

提取TCGA的项目

library(dplyr)
projects <- projects %>% 
as.data.frame() %>% 
select(project_id,tumor) %>% 
  filter(grepl(pattern="TCGA",project_id))

返回的是33行。

使用TCGAbiolinks批量下载TCGA的表达量数据

测试单个项目的数据下载

我们以卵巢癌为例,下载表达量数据。
卵巢癌的代号是TCGA-OV
表达矩阵获取四部曲
第一,使用GDCquery函数来查询下载信息

query.exp <- GDCquery(project = "TCGA-OV", 
                  data.category = "Transcriptome Profiling",
                  data.type = "Gene Expression Quantification",
                  workflow.type = "HTSeq - Counts")

第二,使用GDCdownload函数来下载数据

GDCdownload(query.exp)

下载的数据,会默认在当前工作目录创建GDCdata文件夹,并存放在里面。

使用TCGAbiolinks批量下载TCGA的表达量数据

第三,使用GDCprepare批量读取数据

pre.exp <- GDCprepare(query = query.exp, save = TRUE, save.filename = "TCGA_OV_RNAseq.Rdata")

此时如果想把读取的数据保存为Rdata,下次直接加载,就把
save参数设置为TRUE
save.filename给一个名称。
这时候,这个pre.exp是个对象,需要从里面提取表达量的信息

第四,使用assay函数提取表达量信息

countsdata <- SummarizedExperiment::assay(pre.exp)

现在就获得了行是基因,列是样本的表达量矩阵啦

使用TCGAbiolinks批量下载TCGA的表达量数据

批量下载TCGA的数据

能做一个就能做多个,能做多个的前提是能做一个,而且要清晰定义一次操作的所有步骤。

上面四步已经做完了,现在可以批量做了。

projects这个表达已经获取到了,循环他就可以。

使用TCGAbiolinks批量下载TCGA的表达量数据
dir.create("TCGA_RNA_data")
for (i in 1:nrow(projects)) {
## 0.运行信息
  print(paste0("Downloading number ",i,",project name: ",projects$project_id[i]))
## 1.查询信息
  query.exp = GDCquery(project = projects$project_id[i], 
                        data.category = "Transcriptome Profiling",
                        data.type = "Gene Expression Quantification",
                        workflow.type = "HTSeq - Counts")
## 2.正式下载
  GDCdownload(query.exp)
## 3.多个数据合并
  pre.exp = GDCprepare(query = query.exp)
## 4.提取表达量数据
  countsdata = SummarizedExperiment::assay(pre.exp)
## 5.保存数据
  save(countsdata,file = paste0("TCGA_RNA_data/",projects$project_id[i],".Rdata"))
}

经过一段时间的等待,下载完毕,

使用TCGAbiolinks批量下载TCGA的表达量数据

总共大概4个G,自己可以试试这个技能。

这些可能会帮助到你: 问答社区 | 共享百度SVIP | 留言建议

欢迎入群交流:生信分析群: 732179952 · Meta分析群: 797345521 · 医学科研交流群: 797345521

发表评论

登录后才能评论