TCGA数据库更新啦！跟着之前TCGA系列教程进行官方数据下载和合并实操时遇到问题和报错，反馈主要有两点：

1. 原本GDC-client下载的转录组数据文件格式后缀为.gz,更新后后缀为.tsv；

2. 除了格式不同，文件内存放的数据及组织形式也发生了改变。

确实，TCGA悄悄在四月初更新了，但全新的组织方式对于我们来说是更佳便捷的！

一、以TCGA-KIRC为例，进行新版TCGA-KIRC数据下载与整理

第一步，下载临床和表达矩阵

1.下载gdc-client软件。进入https://gdc.cancer.gov/access-data/gdc-data-transfer-tool，选择自己电脑匹配的软件，下载到工作目录解压。

2. 下载表达数据。进入https://portal.gdc.cancer.gov/repository

然后进入Cart

点击Dowload下拉菜单中的Manifest，保存为gdc_manifest_expdata.txt文件到工作目录。这个文件用于下载每个样本表达矩阵。

3. 下载表达数据的Metadata文件，保存为metadata.cart.json文件到工作目录。这个文件用于TCGA_ID和文件名ID转化。

4. 下载临床数据。进入https://portal.gdc.cancer.gov/repository

点击Manifest，保存为gdc_manifest_clinical.txt文件到工作目录。这个文件用于下载临床信息。

接下来就是在Rstudio中下载文件了。

options(stringsAsFactors = F)

library(stringr)

project="TCGA-KIRC"

if(!dir.exists("clinical"))dir.create("clinical")

if(!dir.exists("expdata"))dir.create("expdata")

dir()

# [1] "clinical" "expdata" "gdc-client.exe" "gdc_manifest_clinical.txt"

# [5] "gdc_manifest_expdata.txt" "metadata.cart.json" "step01_prepare_data.Rmd" "TCGA_KIRC.Rproj"

command1 <- "./gdc-client download -m gdc_manifest_clinical.txt -d clinical"

command2 <- "./gdc-client download -m gdc_manifest_expdata.txt -d expdata"

system(command = command1) # download clinical data

system(command = command2) # download expression data

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

#[1] 537

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

#[1] 613

至此，537个临床文件，613个测序文件都下载完了。

第二步，准备临床信息文件

在这一步中，将会合并所用样本的临床信息。

library(XML)

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

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]

# additional_studies tumor_tissue_site histological_type

# 1 Kidney Kidney Clear Cell Renal Carcinoma

# 2 Kidney Kidney Clear Cell Renal Carcinoma

# 3 Kidney Kidney Clear Cell Renal Carcinoma

第三步，准备表达矩阵

在这一步中，将会把所用样本的表达矩阵汇总，随机去除重复基因名称，并以gene_name作为合并后表达矩阵的行名。

count_files = dir("expdata/",pattern = "*.tsv$",recursive = T)

exp = list()

for(i in 1:length(count_files)){

exp[[i]] = read.table(paste0("expdata/",count_files[[i]]),header=T,sep="\t")

exp[[i]] = exp[[i]][-(1:4),] # The first 4 rows contain unneeded information

exp[[i]] = exp[[i]]$unstranded # the fourth column (unstranded) was the count value

}

exp = as.data.frame(do.call(cbind,exp))

dim(exp)

#[1] 60660 613

exp[1:4,1:4]

# V1 V2 V3 V4

# 1 2157 4455 10949 2375

# 2 26 13 33 13

# 3 978 1610 1621 1264

# 4 604 831 462 764

#TCGA ID and file name match

meta = jsonlite::fromJSON("metadata.cart.json")

ID = sapply(meta$associated_entities,

function(x){x$entity_submitter_id})

file2id = data.frame(file_name = meta$file_name,

ID = ID)

count_files2 = stringr::str_split(count_files,"/",simplify = T)[,2]

table(count_files2 %in% file2id$file_name)

# TRUE

# 613

file2id = file2id[match(count_files2,file2id$file_name),]

identical(file2id$file_name,count_files2)

#[1] TRUE

colnames(exp) = file2id$ID

gene_name = data.table::fread(paste0("expdata/",count_files[1]))$gene_name

gene_name = gene_name[-seq(1,4)] # The first 4 rows contain unneeded information

exp = cbind(gene_name=gene_name,exp)

dim(exp)

#[1] 60660 614

exp = exp[!duplicated(exp$gene_name),]

rownames(exp) = exp$gene_name

exp = exp[,-1]

dim(exp)

#[1] 59427 613

第四步，过滤基因

这里只保留在50%样本中都表达的基因。

#gene filter

exp = exp[apply(exp, 1, function(x) sum(x > 0) > 0.5*ncol(exp)), ] # only keep genes express in half of samples

dim(exp)

#[1] 32809 613

第五步，确定分组信息

TCGA-KIRC数据中包含了541个Tumor和72个Normal样本。

#group information

library(stringr)

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)

# normal tumor

# 72 541

第六步，保存数据。

if(!dir.exists("data"))dir.create("data")

save(exp,clinical,Group,project,file = paste0("data/",project,"_gdc.Rdata")

二、如何使用TCGAbiolinks包下载最新的TCGA数据

第一步：重新安装R包

TCGAbiolinks及TCGAbiolinksGUI.data 包最近都有更新，但在Bioconductor上还没有更新，所以要从git.hub上安装。此外，一些依赖包也需要安装。安装前记得删除之前的TCGAbiolinks及TCGAbiolinksGUI.data 包。

if (!require("devtools"))

install.packages("devtools")

devtools::install_github("BioinformaticsFMRP/TCGAbiolinksGUI.data", host = "https://api.github.com")

BiocManager::install("AnnotationHub")

BiocManager::install("ExperimentHub")

BiocManager::install("SummarizedExperiment")

BiocManager::install("GenomeInfoDbData")

BiocManager::install("GenomicDataCommons")

devtools::install_github("BioinformaticsFMRP/TCGAbiolinks",

host = "https://api.github.com")

第二步：读取基因注释文件

Ginfo<-read.table("overlapTable27.txt", header=T,sep="\t",check.names=F,row.names = 1) #可点击下方下载

第三步：开始下载

此处以TCGA-PAAD为例。更新后的STAR-Counts数据，由于包含了不同类别的数据。所以特别大。

library(TCGAbiolinks)

projects <- c("TCGA-PAAD")#选择PAAD数据集

# RNA Transcriptome

query <- GDCquery(

project = projects,

data.category = "Transcriptome Profiling",

data.type = "Gene Expression Quantification",

workflow.type = "STAR - Counts",#选择STAR-Count数据格式

#barcode = ids

)

GDCdownload(query,

method = "api",

directory = "./data",

files.per.chunk = 10)

第四步：数据转换

rna <- GDCprepare(query, save = FALSE,

summarizedExperiment = T, directory = "./data",

remove.files.prepared = F)

expMatrix<-rna@assays@data@listData$tpm_unstrand

#上面选择了tpm_unstrand数据格式，可以更换为其他格式

rownames(expMatrix)<-rna@rowRanges$gene_id

colnames(expMatrix)<-rna@colData@listData$barcode

colnames(expMatrix) <- substr(colnames(expMatrix),

start = 1,stop = 15)

第五步：基因注释

rownames(Ginfo)<-Ginfo$gene_id

comgene <- intersect(rownames(expMatrix),rownames(Ginfo))

expr <- as.data.frame(expMatrix)[comgene,]

Ginfo <- Ginfo[comgene,]

expr$gene <- Ginfo$genename

expr <- expr[!duplicated(expr$gene),]

Ginfo <- Ginfo[rownames(expr),]

rownames(expr) <- expr$gene

expr <- expr[,-ncol(expr)]

第六步：获得数据的表达谱

第7步：获取临床数据

clinicaldata<-as.data.frame(rna@colData)

第8步：如何选择合适的数据类型来分析

在TCGA的数据中，每个基因的unstranded表达数值=stranded_first+stranded_second。以及官方转换好的FPKM/TMP数据均基于unstranded的数据，因此笔者认为，我们日常常规分析的时候，就仍然使用unstranded的数据就好，stranded的数据可能适用于更精细的研究。

转自：“科研Z库”微信公众号

如有侵权，请联系本站删除！