2022年4月份TCGA进行了更新，测序数据变成了STAR-count文件。以往count文件是htseq.counts.gz格式，现在变成star_gene_counts.tsv格式了。
对比下count文件内容：
image.png
可以看得出来，新版的count文件第一列为gene_id，第二列为gene_name，第四列unstranded也就是count value。另外前面1-4行是注释信息，在后续处理时需要去除。
image.png
接下来以TCGA-KIRC为例，进行新版TCGA-KIRC数据下载与整理。
第一步，下载临床和表达矩阵
-
下载gdc-client软件。进入https://gdc.cancer.gov/access-data/gdc-data-transfer-tool，选择自己电脑匹配的软件，下载到工作目录解压。
image.png
-
下载表达数据。进入https://portal.gdc.cancer.gov/repository
image.png
image.png
然后进入Cart
image.png
点击Dowload下拉菜单中的Manifest，保存为gdc_manifest_expdata.txt文件到工作目录。这个文件用于下载每个样本表达矩阵。 -
下载表达数据的Metadata文件，保存为metadata.cart.json文件到工作目录。这个文件用于TCGA_ID和文件名ID转化。
image.png
-
下载临床数据。进入https://portal.gdc.cancer.gov/repository
image.png
image.png
点击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
第五步，确定分组信息
#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
TCGA-KIRC数据中包含了541个Tumor和72个Normal样本。
第六步，保存数据
if(!dir.exists("data"))dir.create("data")
save(exp,clinical,Group,project,file = paste0("data/",project,"_gdc.Rdata"))
致谢：
本文代码参考了“生信技能树”、“果子学生信”公众号中的代码。
王雄 2022-7-10
网友评论