二、TCGA数据下载

TCGA癌症类型

1.数据下载类型

(1) gdc-client（TCGA官方渠道）
(2) Xena（国外一学校提供的全但不好下，多尝试）
(3) gdcRNAtools（gdc-client的衍生工具可能更新跟不上）

1.没有count数据的情况

2.gdc-client数据下载步骤

2.1

(1)下载manifest文件

1）bing → 搜索GDC tcga → 首页点击Respository

（1）

2）勾选自己需要的case和file类型。以CHOL为例： case——Project选择TCGA-CHOL，其他按需求选择

file-选择如图：

（2）

3）选好后，点击右侧manifest键下载对应的清单文件

（3）

(2)下载表达数据（组学数据）

1）case选择不变，file选择如下：

（1）

2）选好后，点击右侧manifest键下载对应的清单文件

（2）

(3)下载数据

1）将manifest文件放在工作目录下

options(stringsAsFactors = F)
library(stringr)
cancer_type="TCGA-CHOL"#前缀
if(!dir.exists("clinical"))dir.create("clinical")#判断文件是否在工作目录下，不存在就新建
if(!dir.exists("expdata"))dir.create("expdata")
dir()#展示工作目录下都有哪些文件

2）在terminal下输入口令

./gdc-client --help
#（每次运行要把./gdc-client粘贴到工作目录下, 
# window电脑有的需要加exe有的不需要，MAC不需要）

./gdc-client download -m gdc_manifest_cl.2020-03-23.txt -d clinical 
./gdc-client download -m gdc_manifest_expdata.2020-03-23.txt -d expdata

回到R运行命令，查看下载文件个数

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

3）整理临床信息

（3）

library(XML)
result <- xmlParse("./clinical/142aea0e-7a7b-4ac4-9dbb-0f62e2379599/nationwidechildrens.org_clinical.TCGA-W5-AA2O.xml")#读取xml下的一个文件，/后的文件名是tab补全的,两下tab

rootnode <- xmlRoot(result)#提取节点
xmlSize(rootnode)#输出几个节点

print(rootnode[1])#输出第一个节点
print(rootnode[2])#真正有用的在这个节点
xmldataframe <- xmlToDataFrame(rootnode[2])#把xml转成数据框
head(t(xmlToDataFrame(rootnode[2])))

xmls = dir("clinical/",pattern = "*.xml$",recursive = T)#读取全部xml文件，recursive层级的匹配即下一级文件夹也算
cl = list()#用循环读取
for(i in 1:length(xmls)){
  result <- xmlParse(paste0("clinical/",xmls[[i]]))
  rootnode <- xmlRoot(result)
  cl[[i]] = xmlToDataFrame(rootnode[2])
}#读文件、取节点、转数据框与上面单个文件步骤相同
clinical <- do.call(rbind,cl)#简化，按行拼接在一起，变成了整个数据框
clinical[1:3,1:3]#clinical是为生存分析准备的，差异分析太详细了，但如果分组按临床信息也可以用

4）整理表达矩阵

#探索数据：先任选两个counts文件读取，并观察geneid的顺序是否一致。
#options(stringsAsFactors = F)#4.0以前版本，凡是要读取和生成数据框都要运行一下
x = read.table("expdata/03aee74e-4e37-4a58-a720-c90e807d2f40/be5bc6a0-9720-47ac-953e-fa8d0c32cd82.htseq.counts.gz",
               row.names = 1,sep = "\t")#如果数据没有列名和多列后面可写成sep = "\t",header = T)[1]
x2 = read.table("expdata/10d08172-48d2-49e7-b760-721163492cc1/c1071bcd-5a0c-4e09-a578-fc4b6dbe26ad.htseq.counts.gz",
                row.names = 1,sep = "\t")
identical(rownames(x),rownames(x2))#识别是否一致


#批量读取所有的counts.gz文件
count_files = dir("expdata/",pattern = "*.htseq.counts.gz$",recursive = T)
exp = list()

for(i in 1:length(count_files)){
  exp[[i]] <- read.table(paste0("expdata/",count_files[[i]]),row.names = 1,sep = "\t")
}#如果数据没有列名和多列后面可写成sep = "\t",header = T)[1]

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


#发现表达矩阵
没有列名。
#解决办法：找到一个文件名与样本ID一一对应的文件，cart-json文件。
在GDC选TCGA、CHOL、RNA-Seq、Counts  → 点Add All File to Cart  
→  cart-json文件 → 点  Metadata下载
meta <- jsonlite::fromJSON("metadata.cart.2020-03-23.json")
colnames(meta)

ids <- meta$associated_entities;class(ids)
#entities列中都是数据框，数据框里有TCGA的样本ID(例：TCGA-CJ-4923-01A-01R-1425-13)
ids[[1]]
ids[[1]][,1]#取出来，不一定第几行，换数据要改

ID = sa pply(ids,function(x){x[,1]})#批量提取样本ID,与上面ids[[1]][,1]同步改动
file2id = data.frame(file_name = meta$file_name,
                     ID = ID)#取出file_name和样本ID的对应关系

#将两个文件的file_name变成名字一样只有顺序不一样
head(file2id$file_name,2)
head(count_files,2)
count_files2 = stringr::str_split(count_files,"/",simplify = T)[,2]#str_split按/做切割
table(count_files2 %in% file2id$file_name)#验证是不是数据一样

#调整顺序添加列名
file2id = file2id[match(count_files2,file2id$file_name),]#调整file2id的顺序，match放在中括号逗号前面代表调整行的顺序
identical(file2id$file_name,count_files2)#检查顺序
colnames(exp) = file2id$ID
exp[1:4,1:4]

#过滤一下那些在很多样本里表达量都为0的基因。过滤标准不唯一
dim(exp)
#exp = exp[rowSums(exp)>0,]
exp = exp[apply(exp, 1, function(x) sum(x > 1) > 9), ]
#保留一行里大于一的样本数量是否大于九，九为样本数量比较小的一组的数量
#另一种标准：将九换成ncol(exp)*0.75,可选乘多少
#最宽松标准：exp=exp[rouSums(exp)>0],所有样本表达量和不为零
dim(exp)
exp[1:4,1:4]
exp = as.matrix(exp)

#根据样本ID的第14-15位，给样本分组（tumor和normal）
table(str_sub(colnames(exp),14,15))#查看每个分组各有多少个
Group = ifelse(as.numeric(str_sub(colnames(exp),14,15)) < 10,'tumor','normal')#分组命名
Group = factor(Group,levels = c("normal","tumor"))#把分组变成因子，前面是对照组
table(Group)#看和前面数量是否一致

save(exp,clinical,Group,cancer_type,file =paste0(cancer_type,"_gdc.Rdata"))
#保存，exp表达矩阵、clinical临床信息、group根据表达矩阵列名实现的分组、cancer_type癌症名

3.Xena数据下载

搜索Xena UCSC → launch Xena → DATA SETS → Ctrl+F搜索需要的癌症简称 → 选GDC数据（xena上的数据log过要使用得复原回去）

# 1.xena
if(F){
  download.file(url = "https://gdc.xenahubs.net/download/TCGA-CHOL.htseq_counts.tsv.gz",destfile = "counts.tsv.gz")
  download.file(url = "https://gdc.xenahubs.net/download/TCGA-CHOL.GDC_phenotype.tsv.gz",destfile = "phenotype.tsv.gz")
  download.file(url = "https://gdc.xenahubs.net/download/TCGA-CHOL.survival.tsv.gz",destfile = "survival.tsv.gz")
}
dat = read.table("counts.tsv.gz",check.names = F,row.names = 1,header = T)
#逆转log
dat = as.matrix(2^dat - 1)
dat[1:4,1:4]
as.character(dat[1:100,1:10]) #有一些小数

# 用apply转换为整数矩阵
exp = apply(dat, 2, as.integer)
exp[1:4,1:4] #行名消失
rownames(exp) = rownames(dat)#找回行名

clinical = read.table("phenotype.tsv.gz",fill = T,header = T,sep = "\t")
surv = read.table("survival.tsv.gz",header = T)
clinical[1:4,1:4]
surv[1:4,1:4]

# 2.GDCRNATools
针对GDC的R包可以自己学习
# http://bioconductor.org/packages/devel/bioc/vignettes/GDCRNATools/inst/doc/GDCRNATools.html

4.GEO数据下载

library(GEOquery)
eSet = getGEO("GSE162550",destdir = F,getGPL = F)#按芯片方式下载
rm(list = ls())
dat=
read.table("GSE162550_gene_sample_count_with_symbol.xls",
                 fill = T,sep = "\t",header = T)
table(!duplicated(dat$Symbol))

o = order(rowSums(dat[,4:9]),decreasing = T)
dat = dat[o,] 
dat = dat[!duplicated(dat$Symbol),]
dat = dat[dat$Symbol!="---",]

exp = dat[,4:9]
rownames(exp) = dat$Symbol

• 是学习生信技能树课程的笔记哦