把GDC下载的多个TCGA文件批量读入R
上一次我们已经从GDC下载了TCGA肿瘤数据库的数据
多个单独文件合并到单一文件夹
但是下载是一个个压缩文件，更要命的是，这些压缩文件还藏在一个个单独的文件夹中。
我在以前的一个帖子里讲述了，如何把这些压缩文件放在同一个文件夹中。
TCGA提取非编码RNA并完成下游分析
直觉告诉我们，一个个复制粘贴是可行的。对！有时候我们就是要凭着自己的直觉做事，这一次我们用更加简单的代码来实现这个功能(文末有免费操作视频)。
首先所有的原始数据存在rawdata这个文件夹中，我们现在创建一个新的文件夹叫data_in_one,用来存放所有的压缩文件。
dir.create("data_in_one")
如果要用for循环来做这个事情，诀窍只有一个：
你要清晰地定义，如何做一件事，你能做一件事，就能做多件事。
实现一个，很简单
1. 打开rawdata下面的第一个文件夹
2. 找到该文件夹下面的压缩文件
3. 复制粘贴到新的文件夹data_in_one中
好了，根据这个思路，我们顺利地写出了这个循环
for (dirname in dir("rawdata/")){
file <- list.files(paste0(getwd(),"/rawdata/",dirname),pattern = "*.counts")
file.copy(paste0(getwd(),"/rawdata/",dirname,"/",file),"data_in_one")
}
很快，所有的文件被复制到了新的文件夹，在我的电脑上用时1秒，很爽！
接下来就好办了，我们批量把这个数据读入R语言即可，但是问题是，TCGA数据是有TCGA barcode的，类似于下面这个：
"TCGA-06-0138-01A-02R-1849-01"
但是我们现在只有这样的文件名
"76c6e110-9a94-4fcf-a534-a21ba4698f86.htseq.counts.gz"
如果就这样把数据读入R，我们会分不清样本间的区别，无法进行下游分组操作。
所以，现在要找到文件名称和TCGA id之间的对应关系。
找出文件名对应的TCGA id
这个对应关系在上次下载的metadata文件中，这个文件是json格式的，很复杂，需要专门的函数读取，我对比了好几种函数，发现jsonlite中的fromJSON函数是最好的。
metadata <- jsonlite::fromJSON("metadata.cart.2019-01-14.json")
我们再用for循环提取对应的两者对应关系
naid_df <- data.frame()
for (i in 1:nrow(metadata)){
naid_df[i,1] <- metadata$file_name[i]
naid_df[i,2] <- metadata$associated_entities[i][[1]]$entity_submitter_id
}
colnames(naid_df) <- c("filename","TCGA_id")
速度很快，效果不错
批量读取数据
现在我们可以放心地读取文件了，在以前，我喜欢把他们解压缩了再读取，但是洲跟提醒我R语言可以直接读取gz压缩文件，那就好办了，先读一个试试效果
test <- data.table::fread(paste0("data_in_one/",naid_df$filename[1]))
数据有两列，一列是ensemble id，一列是对应基因的counts数目。
那么接下来的思路是，我先创建一个大的数据框，然后把每次读取出来的第二列数据变成一列即可。
expr_df <- data.frame(matrix(NA,nrow(test),nrow(naid_df)))
for (i in 1:nrow(naid_df)) {
print(i)
expr_df[,i]= data.table::fread(paste0("data_in_one/",naid_df$filename[i]))[,2]
}
给读入的数据添加列名和基因名称
每一个文件读取时都对应了一个TCGA id，所以用对应的TCGA id 给获得的数据命名即可
colnames(expr_df) <- naid_df$TCGA_id
读取任意单个文件，把他的第一列合并到大的数据框上，作为基因名称
gene_id <- data.table::fread(paste0("data_in_one/",naid_df$filename[1]))$V1
expr_df <- cbind(gene_id=gene_id,expr_df)
去除后5行,保存数据成Rdata格式
这个时坑，只有自己掉进去才知道，最后5行不是我们需要的，可以用tail这个函数查看
tail(expr_df$gene_id,10)
去除最后5行
expr_df <- expr_df[1:(nrow(expr_df)-5),]
看一下现在的数据
保持数据为Rdata格式
save(expr_df,file = "expr_df.Rdata")
以后要用的时候，load一下即可。
load(file = "expr_df.Rdata")
好了，打完收工，下次我们再来讲基于这个数据的差异分析以及如何标准化。
整个用到的文件以及操作讲解的视频我已经上传，回复"TCGA0123"即可免费获取。
为了不让大家内疚，没有开赞赏，点击右下方"好看"即可。