TCGA数据挖掘（3）：miRNA数据的下载

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
knitr::opts_chunk$set(fig.width = 6,fig.height = 6,collapse = TRUE)
knitr::opts_chunk$set(message = FALSE)

本文的内容是用GDC下载并整理表达矩阵和临床信息数据。

1.从网页选择数据，下载manifest文件

数据存放网站：https://portal.gdc.cancer.gov/

在Repository勾选自己需要的case和file类型。以CHOL为例：
case-Project选择TCGA-CHOL。

file-选择如图：

左右分别是expdata 和clinical的样本选择截图。选好后，点击右侧manifest键下载对应的清单文件。

2.使用gdc-client工具下载

注意：

将gdc-client(mac)或gdc-client.exe(windows)放在工作目录下；

将manifest文件放在工作目录下。

options(stringsAsFactors = F)
library(stringr)
cancer_type="TCGA-KIRC"
if(!dir.exists("clinical"))dir.create("clinical")
if(!dir.exists("expdata"))dir.create("expdata")
dir()
#下面两行命令在terminal完成
#./gdc-client.exe download -m gdc_manifest.2020-08-18_KIRC_clinical.txt -d clinical
#./gdc-client.exe download -m gdc_manifest.2020-08-18_KIRC_expdata.txt -d expdata

length(dir("./clinical/"))
length(dir("./expdata/"))

可以看到，下载的文件是按样本存放的，我们需要得到的是表格，需要将他们批量读入R语言并整理。

3.整理临床信息

library(XML)
result <- xmlParse("./clinical/000b43b4-1fe2-4aaf-9552-7d519eae17a6/nationwidechildrens.org_clinical.TCGA-CJ-4916.xml")
rootnode <- xmlRoot(result)
rootsize <- xmlSize(rootnode)
print(rootnode[1])
#print(rootnode[2])
xmldataframe <- xmlToDataFrame(rootnode[2])
head(t(xmlToDataFrame(rootnode[2])))

xmls = dir("clinical/",pattern = "*.xml$",recursive = T)

td = function(x){
  result <- xmlParse(file.path("clinical/",x))
  rootnode <- xmlRoot(result)
  xmldataframe <- xmlToDataFrame(rootnode[2])
  return(t(xmldataframe))
}

cl = lapply(xmls,td)
cl_df <- t(do.call(cbind,cl))
cl_df[1:3,1:3]
clinical = data.frame(cl_df)
clinical[1:4,1:4]

4.整理表达矩阵

探索数据：先任选两个counts文件读取，并观察geneid的顺序是否一致。

由此可知，他们的geneid顺序是一致的，可以直接cbind，不会导致顺序错乱。

批量读取所有的counts.gz文件。

count_files = dir("expdata/",pattern = "*mirnas.quantification.txt$",recursive = T)

ex = function(x){
  result <- read.table(file.path("expdata/",x),row.names = 1,sep = "\t",header = T)[1]
  return(result)
}
head(ex("000c6bb3-c47d-49a1-ad0f-57093521f922/1260cbd9-170e-4032-b45e-99620b9304da.mirbase21.mirnas.quantification.txt"))

exp = lapply(count_files,ex)
exp <- do.call(cbind,exp)
dim(exp)
exp[1:4,1:4]

发现问题：这样产生出来的表达矩阵没有列名。

解决办法：找到一个文件名与样本ID一一对应的文件。cart-json文件。

meta <- jsonlite::fromJSON("metadata.cart.2020-03-23.json")
colnames(meta)
ids <- meta$associated_entities;class(ids)
ids[[1]]
class(ids[[1]][,1])

可以看到，meta$associated_entities是个列表，这个列表里包含数据框，数据框的第一列内容就是tcga样本id。

注意，换了数据需要自己探索存放在哪一列。不一定是完全一样的，需要确认清楚。

ID = sapply(ids,function(x){x[,1]})
file2id = data.frame(file_name = meta$file_name,
                     ID = ID)

文件名与TCGA样本ID的对应关系已经得到，接下来是将其添加到表达矩阵中，成为行名。需要找到读取文件的顺序，一一对应修改。

head(file2id$file_name)
head(count_files)
count_files2 = stringr::str_split(count_files,"/",simplify = T)[,2]
count_files2[1] %in% file2id$file_name

count_files2的顺序就是列名的顺序，根据它来调整file2id的顺序。此处需要再次理解一下match函数。

file2id = file2id[match(count_files2,file2id$file_name),]
colnames(exp) = file2id$ID
exp[1:4,1:4]

表达矩阵整理完成,需要过滤一下那些在很多样本里表达量都为0的基因。过滤标准不唯一。

dim(exp)
exp = exp[apply(exp, 1, function(x) sum(x > 1) > 9), ]
dim(exp)
exp[1:4,1:4]

分组信息

根据样本ID的第14-15位，给样本分组（tumor和normal）

table(str_sub(colnames(exp),14,15))
group_list = ifelse(as.numeric(str_sub(colnames(exp),14,15)) < 10,'tumor','normal')
group_list = factor(group_list,levels = c("normal","tumor"))
table(group_list)
save(exp,clinical,group_list,cancer_type,file = paste0(cancer_type,"gdc.Rdata"))