写在前面：本文为微信公众号:生信星球的数据挖掘线上班的随堂笔记，感谢小洁老师的付出！
TCGA
TCGA能做的
0.准备工作
0-1 Terminal
Terminal操作系统分为两个部分，一部分称作内核，另一部分成为用户交互界面。终端就是连接内核与交互界面的这座桥，它允许用户在交互界面上打开一个叫做「Terminal 终端」的应用程序，在其中输入命令，系统会直接给出反馈。
- Git软件：方便管理命令行。使用时右击文件夹[get bush here]
- Rstudio里有terminal，可连接git。操作方式：tools → terminal → git bash
0-2 Rmd
- Rmd文件介绍：
---
【头文件】
output: html_document
---
···
【代码块】
···
【正文部分】（无需#）
- 插入代码块 ctrl+alt+i
-
运行代码
- 基本语法同markdown
- 渲染导出文件
- 代码块选项
- echo = FALSE 该代码块运行但不显示
-
eval = FALSE 该代码块不运行但显示
1.TCGA数据下载方式
-
TCGA癌症类型
- 根据mRNA，芯片等等得到预后的结果。
1. TCGA gdc
- 官方工具
1. 网页选择：manifest（清单文件）——file+case
- 数据：counts...+json——RNA-seq+miRNA-seq+exon...
- counts文件需整理为表达矩阵。
- json储存样本文件的详细信息。
- 临床信息：xml文件——病人ID，生死，分期，死亡时间等等
-
病人临床信息多于RNA-seq/组学，因为有的人没全测。
样本ID第14-15位为tumor或normal情况 - 下载表达矩阵（.htseq.counts.gz）和clinical数据（.xml）
terminal操作
- ./gdc-client download --help #工作目录下打开gdc-client并执行下载功能-的help文件
-
-d 为文件夹名,-m 文件信息
#下面两行命令在terminal完成，需要更改-m后面的文件名
./gdc-client.exe download -m gdc_manifest_cl.2020-03-23.txt -d clinical
./gdc-client.exe download -m gdc_manifest_mrna.2020-03-23.txt -d mrna
- gdc和manifast需要在同一个目录下。
- 使用'length(dir("./clinical/"))'在控制台查看下载进度。
- 1.从网页选择数据，下载manifest文件,数据存放网站,然后在Repository勾选自己需要的case和file类型。
- 2.使用gdc-client工具下载
注意：将gdc-client(mac)或gdc-client.exe(windows)放在工作目录下；将manifest文件放在工作目录下。
options(stringsAsFactors = F)
cancer_type="TCGA-CHOL"
if(!dir.exists("clinical"))dir.create("clinical")
if(!dir.exists("mrna"))dir.create("mrna")
dir()
#下面两行命令在terminal完成
#./gdc-client download -m gdc_manifest_cl.2020-03-23.txt -d clinical
#./gdc-client download -m gdc_manifest_mrna.2020-03-23.txt -d mrna
length(dir("./clinical/"))
length(dir("./mrna/"))
3.整理临床信息
library(XML)
result <- xmlParse("./clinical/142aea0e-7a7b-4ac4-9dbb-0f62e2379599/nationwidechildrens.org_clinical.TCGA-W5-AA2O.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)
#列出了文件夹所有的xml文件
#以下为以上的数据处理的函数td
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.整理表达矩阵
count_files = dir("mrna/",pattern = "*.htseq.counts.gz$",recursive = T)
ex = function(x){
result <- read.table(file.path("mrna/",x),row.names = 1,sep = "\t")
return(result)
}
head(ex("03aee74e-4e37-4a58-a720-c90e807d2f40/be5bc6a0-9720-47ac-953e-fa8d0c32cd82.htseq.counts.gz"))
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])
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的基因。过滤标准不唯一。这里取出至少有9个样本表达量大于1的基因们。
dim(exp)
exp = exp[apply(exp, 1, function(x) sum(x > 1) > 9), ]
dim(exp)
exp[1:4,1:4]
5.分组信息
- 本样本根据样本ID的第14-15位，给样本分组（tumor和normal），随数据更改——如本次为"TCGA-W5-AA36-01A-11R-A41I-07"。
table(stringr::str_sub(colnames(exp),14,15))
group_list = ifelse(as.numeric(substr(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"))
