2022-09-29 新版TCGA数据下载

2022-09-30 本文已影响0人学习生信的小兔子

参考：公众号果子学生信 https://mp.weixin.qq.com/s/rdFnq6jCMIjmrWI8A8fS5g

GDC网站

https://portal.gdc.cancer.gov/

点击repository

点击case

确定数据下载的组织

点击file选择下载数据的类型

接下来把数据选择的数据加入到购物车，购物车里面数量会变成文件数目，当前是1226

下载数据

此时点击购物车，就会进入下载页面，因为当前数据挺大的，不建议直接下载而是采用他推荐的GDC Data Transfer Tool来进行而使用这个工具，需要两个文件，一个是Manifest，一个是metadata

这两个文件，会根据下载时间的不同，生成对应日期的名称

下载GDC Data Transfer Tool工具

解压，放在同一个文件夹中

R下载数据

### 1.创建数据下载文件夹
if(!file.exists("rawdata")) dir.create("rawdata")

### 2.生成下载代码
### 这个就是之前下载的manifest文件，每个人的不一样，注意修改
manifest <- "gdc_manifest_20220929_023209.txt"
### 这个是下载的文件夹
rawdata <- "rawdata"
command <- sprintf("./gdc-client download -m %s -d %s",
                   manifest,
                   rawdata)
### 3.执行下载操作
system(command = command)

这里我下载了3个小时。。。

整理数据

这时候所有数据都在rawdata这个文件夹中，总共有1226个文件夹，每一个里面点开后有一个文件

看看这个文件（44beea06-28a7-490a-b6bb-4e023230e51c.rna_seq.augmented_star_gene_counts.tsv）的样子

其中第一列 gene_id我们是需要的，第四列unstranded 就是传统的counts 数据，以前的TCGA数据只提供两列。
其中第二列是转换过后的基因名称，第三列是基因类型，包括编码和非编码信息，这个可以保留一份用来做基因注释
我们只需要第1列和第4列

1.把所有的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")  
}

合并以后的文件都放在同一个文件夹

2.提取文件名称和TCGAID 的对应关系【metadta文件中】

#提取文件名称和TCGAID 的对应关系
metadata <- jsonlite::fromJSON("metadata.cart.2022-09-29.json")
library(dplyr)
metadata_id <- metadata %>% 
  dplyr::select(c(file_name,associated_entities))

metadata

metadataid

## 提取对应的文件
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了

3.批量读取数据

##第三步，批量读取数据
myfread <- function(files){
  data = data.table::fread(paste0("data_in_one/",files))
  data = data[-seq(1,4),]
  data = data$unstranded
}

files <- dir("data_in_one")
system.time(test <- myfread(files[1]))测试读取一个文件的所需时间

system.time(f <- lapply(files,myfread))  #lapply批量读取
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
test=expr_df[1:20,1:20]

最后添加基因名称，因为这些基因名称和类型，在每一个样本中都有，所以我们随便读取一个数据，然后添加即可。

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)

顺利完成，该数据格式，第一列是基因，后面都是样本，转录组下游分析，无论是公共数据，还是自测数据，都应该获取这样的数据
总代码·

### 1.从GDC下载counts数据
## https://portal.gdc.cancer.gov/repository
### 1.创建数据下载文件夹
if(!file.exists("rawdata")) dir.create("rawdata")

### 2.生成下载代码
### 这个就是之前下载的manifest文件，每个人的不一样，注意修改
manifest <- "gdc_manifest_20220929_023209.txt"
### 这个是下载的文件夹
rawdata <- "rawdata"
command <- sprintf("./gdc-client download -m %s -d %s",
                   manifest,
                   rawdata)
### 3.执行下载操作
system(command = command)



#合并文件
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")  
}


#提取文件名称和TCGAID 的对应关系
metadata <- jsonlite::fromJSON("metadata.cart.2022-09-29.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")

##第三步，批量读取数据
myfread <- function(files){
  data = data.table::fread(paste0("data_in_one/",files))
  data = data[-seq(1,4),]
  data = data$unstranded
}

files <- dir("data_in_one")
system.time(test <- myfread(files[1]))

system.time(f <- lapply(files,myfread))

expr_df <- as.data.frame(do.call(cbind,f))
test=expr_df[1:20,1:20]


rownames(naid_df) <- naid_df[,1]
naid_df <- naid_df[files,]
colnames(expr_df) <- naid_df$TCGA_id
test=expr_df[1:20,1:20]


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)

save(expr_df,naid_df,file="exp.Rdata")