把GDC下载的多个TCGA文件批量读入R

多个单独文件合并到单一文件夹

但是下载是一个个压缩文件，更要命的是，这些压缩文件还藏在一个个单独的文件夹中。

我在以前的一个帖子里讲述了，如何把这些压缩文件放在同一个文件夹中。
TCGA提取非编码RNA并完成下游分析

直觉告诉我们，一个个复制粘贴是可行的。对！有时候我们就是要凭着自己的直觉做事，这一次我们用更加简单的代码来实现这个功能(文末有免费操作视频)。

首先所有的原始数据存在rawdata这个文件夹中，我们现在创建一个新的文件夹叫data_in_one,用来存放所有的压缩文件。

dir.create('data_in_one')


如果要用for循环来做这个事情，诀窍只有一个：

你要清晰地定义，如何做一件事，你能做一件事，就能做多件事。

实现一个，很简单

• 1. 打开rawdata下面的第一个文件夹

• 2. 找到该文件夹下面的压缩文件

• 3. 复制粘贴到新的文件夹data_in_one中

好了，根据这个思路，我们顺利地写出了这个循环

for (dirname in dir('rawdata/')){  
  ## 使用list.files函数找到rawdata里面单个文件夹下面的压缩文件
  file <- list.files(paste0(getwd(),'/rawdata/',dirname),pattern = '*.counts')  #找到对应文件夹中的内容，pattern可以是正则表达式
  ## 使用file.copy函数复制粘贴压缩文件到data_in_one
  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')

好了，打完收工，下次我们再来讲基于这个数据的差异分析以及如何标准化。