拿去！TCGA改版后转录组数据的批量下载及分享

之前我们讲了TCGA改版后转录组数据下载和整理的方法，
TCGA改版后转录组数据的下载以及整理
但是你实际操作的时候会发现，下载的数据大而且慢。
以TCGA 的BRCA数据为例，总共1226个样本，下载的数据量大概是5个G，按照1Mb/s 的数据，大概是85分钟。
而整理后的数据只有81Mb，所以我们下载了很多冗余数据。

考虑到有些区域的网速受限，我们不仅要教大家如何打鱼，今天还要把池塘里的鱼全部捞上来送给大家。
大概就是授人以渔同时授人以鱼，这也是我一贯的理念，想学就毫无保留的教，不想学就挥金如土的送。

本次批量的方法用的是TCGAbiolinks，上一个版本的数据可以参考
使用TCGAbiolinks批量下载TCGA的表达量数据(老版)

既然用别人的包，那就耐心等待，看作者的心情和精力做事，我去github上看了一些，作者已经更新了新的代码。

批量处理的诀窍在于，你要能够做好单次能做的事情，测试好单次数据的流程。
需要注意的是，数据改版来的突然，所以更新的代码需要使用开发版本的R包

## 1.安装
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install("BioinformaticsFMRP/TCGAbiolinksGUI.data")
BiocManager::install("BioinformaticsFMRP/TCGAbiolinks")
## 2.加载
library(TCGAbiolinks)
## 3.查看版本
packageVersion("TCGAbiolinks")

我当前使用的版本是2.25.0，看看你的是不是落后了很多代。

测试流程

下面测试TCGA GBM数据的两个样本

## 1.查询
query <- GDCquery(
  project = "TCGA-GBM",
  data.category = "Transcriptome Profiling",
  data.type = "Gene Expression Quantification", 
  workflow.type = "STAR - Counts",
  barcode = c("TCGA-14-0736-02A-01R-2005-01", "TCGA-06-0211-02A-02R-2005-01")
)
## 2.下载
GDCdownload(query = query)
## 3.整理
data <- GDCprepare(query = query)

你看这个方法有多简单，不需要鼠标点点点，直接下载，不仅下载，而且他还能够自动把数据整理成表格。
也就是之前我们的那些手艺活，一行代码就搞定了。当前下载的原始数据会自动在工作目录建立新的文件夹GDCdata

值得强调的是，这个GDCprepare 会自动储存单个文件中的六个数据，这在运行代码的过程中，可以看得出来。

而且，这六个数据可以通过代码方便的取出来。都是整理好的。
通过提取这些表格信息可以看到

names(data@assays)

确实有六种

[1] "unstranded"       "stranded_first"   "stranded_second"  "tpm_unstrand"    
[5] "fpkm_unstrand"    "fpkm_uq_unstrand"

可以通过以下代码，按照顺序分别提取不同的数据, assay 函数默认提取的是第一个

library(SummarizedExperiment)
test1 <- assay(data,1)
test2 <- assay(data,2)
test3 <- assay(data,3)
test4 <- assay(data,4)
test5 <- assay(data,5)
test6 <- assay(data,6)

提取后的counts数据如下，是不是超级方便，如果需要的保存一下，就是以后要用到的表达矩阵了。

正式运行

先提取TCGA的项目编号

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

总共33个项目

建立一个新的文件夹用于存储整理好的数据

  if(!file.exists("TCGA_RNA_data")) dir.create("TCGA_RNA_data")

批量执行

  ## 批量运行
  for (i in 1:nrow(projects)) {
    ## 运行信息
    print(paste0("Downloading number ",i,",project name: ",projects$project_id[i]))
    ## 下载信息
    query.exp = GDCquery(project = projects$project_id[i], 
                         data.category = "Transcriptome Profiling",
                         data.type = "Gene Expression Quantification", 
                         workflow.type = "STAR - Counts",
                         )
    ## 正式下载
    GDCdownload(query.exp)
    ## 多个数据合并
    pre.exp = GDCprepare(query = query.exp)
    ## 提取表达量数据
    countsdata = SummarizedExperiment::assay(pre.exp)
    ## 保存数据
    saveRDS(countsdata,file = paste0("TCGA_RNA_data/",projects$project_id[i],".rds"))
  }

速度时好时坏，最终下载的原始数据有43个G，整理好的数据大概700多Mb

如果想要使用某个数据，可以直接加载，因为当时保存的数据是矩阵，所以需要额外转换一下。

countsdata <- readRDS(file = "TCGA_RNA_data_STAR/TCGA-BRCA.rds")
countsdata <- as.data.frame(countsdata)
test <- countsdata[1:100,1:100]

现在就可以畅游数据的海洋了

我把这些数据弄好了链接，在果子学生信公众号回复"果子新版批量TCGA"自助获取吧。
之前老版本的数据就回复"果子批量TCGA"

那如果我需要TPM的数据怎么办，假如你像我一样原始数据(就是一个个文件)下载好了，我们只要改一下流程中assay部分就行，其余的以此类推。

## 提取表达量数据
countsdata = SummarizedExperiment::assay(pre.exp,4)

但是实际上没有必要，因为ounts数据转换为TPM格式十分方便。
又因为TCGA改版后，一套数据同时提供了counts，TPM，FPKM三个格式，我们有了阳性对照，就可以测试自己的转换流程是否正确了。

而获取数据分析的阳性对照，是自我迭代的神器。