TCGA改版后转录组数据的下载以及整理
自己常用的数据库改版，用户能做的事情就两件，
第一，等到天荒地老，
第二，没有枪没有炮我们自己造
拿TCGA转录组的数据举例，他更新后采用了STAR 作为比对工具，在每一个样本的一个文件中提供了基因名称，基因类型，以及counts，fpkm， tpm三种数据格式。
这本来是好事，但是因为很多朋友在一开始并没有学习如何整理这些数据，或者说不会写代码，这些新的数据就处理不了了。
这时候，只要不着急，等待就是最好的方案，使用perl语言脚本的，可以等等看看作者会不会更新一个新的流程。
使用TCGAbiolinksR包下载的朋友也可以等等，因为只要有人发现下载不了数据，就会给作者提要求，作者应该很快就更新。
另外一种方案就是自己动手，改一下代码就可以直接运行了。迟迟没有动手是因为我之前已经下载好了所有的TCGA数据，并没有更新的需求。
使用TCGAbiolinks批量下载TCGA的表达量数据(老版)
我第一次给别人讲TCGA的时候，因为教学需要以及面子上面过不去(主要是面子上过不去)，硬着头皮写了一个R语言版本的教程。
写第一遍是很难受的，但是积累了一些解决报错的经验，克服了自创流程的心理压力，算是入了R语言的门。
以下，以从GDC数据库下载TCGA的BRCA转录组数据为例，更新一下既往的流程，以一个更快的速度来整理这些数据。
首先通过网址进入GDC网站
https://portal.gdc.cancer.gov/
点击repository 进入仓库，
从cases里面确定数据下载的组织
然后files里面选择下载数据的类型，是转录组数据，选择只是在出现多个选项的时候才选，如果只有一个选项，不勾选也没有关系
比如，当前数据格式tsv只有一种，不勾选没有关系
接下来把数据选择的数据加入到购物车，购物车里面数量会变成文件数目，当前是1226
数据下载
此时点击购物车，就会进入下载页面，因为当前数据挺大的，不建议直接下载而是采用他推荐的GDC Data Transfer Tool来进行
而使用这个工具，需要两个文件，一个是Manifest，一个是metadata
这两个文件，会根据下载时间的不同，生成对应日期的名称
这个GDC Data Transfer Tool工具下载也提供了链接，点击进去，往下看，根据自己的电脑系统下载不同的文件。
我们建议两个平台都下载一下，然后都解压，放在同一个文件夹中，而我们接下来的流程就可以不区分windows还是mac，可以直接运行了。
我在想一般人不会这么做，因为只有一台电脑啊，但是我因为要授课，所以希望大家都能用一样的流程，而不是windows讲一遍，然后mac再讲一遍。
然后创建一个新的r project，把这些文件放进去，然后再创建一些文件夹，这些是习惯，非必须。
比如我喜欢把代码都放在一个文件夹，叫做scripts
进入这个项目中，开始下载数据
### 1.从GDC下载counts数据
## https://portal.gdc.cancer.gov/repository
### 1.创建数据下载文件夹
if(!file.exists("rawdata")) dir.create("rawdata")
### 2.生成下载代码
### 这个就是之前下载的manifest文件，每个人的不一样，注意修改
manifest <- "gdc_manifest_20220509_093925.txt"
### 这个是下载的文件夹
rawdata <- "rawdata"
command <- sprintf("./gdc-client download -m %s -d %s",
manifest,
rawdata)
### 3.执行下载操作
system(command = command)
然后思考片刻，会出现下载进度条，我这里的下载速度大概是1MB/s, 不知道你那里是什么情况。
这个代码是洲更给我改进的，他的方便之处在于，无视操作系统，你不要下载一个终端工具进入终端去操作了，那会带来额外的负担。
我很喜欢这个代码的另外一个理由是，我们在教学时不需要在这里额外解释，什么叫做终端，什么叫做cd，什么叫做pwd，这些都是冗余信息。是当前不需要了解的概念。
当然每个人在这里下载的时候，都可能出现一些状况，如果有问题，可以重复运行一下，下载完毕后会显示
数据整理
这时候所有数据都在rawdata这个文件夹中，总共有1226个文件夹，每一个里面点开后有一个文件
这时候，如果第一次弄，一定要打开文件看看，里面有什么，我们需要什么。
其中第一列 gene_id我们是需要的，第四列unstranded 就是传统的counts 数据，以前的TCGA数据只提供两列。
其中第二列是转换过后的基因名称，第三列是基因类型，包括编码和非编码信息，这个可以保留一份用来做基因注释
现在我们知道了，我们只是需要第1列和第4列，为什么以及剩余的都是些什么，后面再来讲。
搞清楚了情况，就很容易写个代码来提取了。
第一步，把所有的tsv文件，拷贝到新的文件夹data_in_one
if(!file.exists("data_in_one")) dir.create("data_in_one")
### 使用for循环来批量做，回顾for循环的要点
for (dirname in dir("rawdata/")){
## 要查看的单个文件夹的绝对路径
mydir <- paste0(getwd(),"/rawdata/",dirname)
##找到对应文件夹中的文件并提取名称，pattern表示模式，可以是正则表达式
file <- list.files(mydir,pattern = "_counts")
## 当前文件的绝对路径是
myfile <- paste0(mydir,"/",file)
#复制这个文件到目的文件夹
file.copy(myfile,"data_in_one")
}
现在所有的tsv文件都到一个文件夹中了
但是这种文件名称实在太难受了，他有个文件名称和TCGA id的对应关系，藏在了metadata中，这个在一开始下载了的
第二步，提取文件名称和TCGAID 的对应关系
metadata <- jsonlite::fromJSON("metadata.cart.2022-05-09.json")
library(dplyr)
metadata_id <- metadata %>%
dplyr::select(c(file_name,associated_entities))
## 提取对应的文件
naid_df <- data.frame()
for (i in 1:nrow(metadata)){
naid_df[i,1] <- metadata_id$file_name[i]
naid_df[i,2] <- metadata_id$associated_entities[i][[1]]$entity_submitter_id
}
colnames(naid_df) <- c("filename","TCGA_id")
得到的对应关系如下，以后每一个文件都知道TCGA id了
第三步，批量读取数据
先测试读取一个文件的代码
我们当前设想，给一个文件名称，读取数据，选择第四列，删掉前四行就行。
第四列就是我们想要的数据，而不要前四行，是因为一开始打开文件时发现前四行是附带信息么没有用处。
myfread <- function(files){
data = data.table::fread(paste0("data_in_one/",files))
data = data[-seq(1,4),]
data = data$unstranded
}
测试一下函数, 在我的电脑上大概0.03s一个,测试时间是为了知道我们总共会话多少时间，比如，2个小时，意味着你有2个小时不能关电脑
files <- dir("data_in_one")
system.time(test <- myfread(files[1]))
lapply批量读取，最终27s完成
system.time(f <- lapply(files,myfread))
完成后是个列表，需要转换为数据框
expr_df <- as.data.frame(do.call(cbind,f))
最终得到60660行，1226列的数据框，其中60660这个数字很重要，他代表的是当前的基因是60660个，而1226就是样本数
当前的数据是裸的，没有行应该是基因，但是没有，列应该是TCGA ID也没有，接下来添加一下。
先添加列名, 对应关系之前已经提取到了，只是要限定以下行的顺序，需要跟读入的样本顺序一致
rownames(naid_df) <- naid_df[,1]
naid_df <- naid_df[files,]
colnames(expr_df) <- naid_df$TCGA_id
最后添加基因名称，因为这些基因名称和类型，在每一个样本中都有，所以我们随便读取一个数据，然后添加即可。
gene_id <- data.table::fread(paste0("data_in_one/",files[1]))$gene_id
gene_id <- gene_id[-seq(1,4)]
expr_df <- cbind(gene_id=gene_id,expr_df)
顺利完成，记住这个数据格式，第一列是基因，后面都是样本，转录组下游分析，无论是公共数据，还是自测数据，都应该获取这样的数据。
如果感兴趣，可以看看跟以前的流程相比，也没有多出什么内容
把GDC下载的多个TCGA文件批量读入R
而批量处理这样的技巧，我们还有个合集，以及一个配套的课程，感兴趣的可以看看。
视频课程: R语言中的批量操作，裂变你的技能
我们下一次会讲一下，如何使用TCGAbiolink 这个R包来获取最新的转录组数据，并且将把全部TCGA癌症的转录组数据下载下来，提供给大家。
总体而言，对于TCGA这种光谱的数据来说，他的任何格式，或者评分，免疫浸润数据，干性数据，都是不用愁的，因为只要有人在用，就会提供下载，我们只要等待就行。
像刚才讲的这些技能，其中用到的数据增删改查，都是R语言的常用技能，是可以举一反三的。最终，我们学习一个新的技能，是为了缓解焦虑，而不是增加烦恼。