TCGA 单样品合并

书接上回：R语言下载TCGA数据 非R包

我们已经获取到TCGA单样品数据，接下来就是合并和处理合并后重复基因名问题。

01.合并前准备

library(rjson)
library(dplyr) 

#从原始数据中把*.tsv文件复制到data_in_one文件夹中，便于后续操作

if(!file.exists("data_in_one")) dir.create("data_in_one")
for (dirname in dir("rawdata/")){  
  mydir <- paste0(getwd(),"/rawdata/",dirname)#注意这里和上书提到的位置一致
  file <- list.files(mydir,pattern = "_counts")
  myfile <- paste0(mydir,"/",file)
  file.copy(myfile,"data_in_one")  
}

02.获取TCGA样品编号

#读取json文件，提取样品TCGA编号名称

metadata <- jsonlite::fromJSON("metadata.cart.2023-11-22.json")
metadata_id <- metadata %>% 
  dplyr::select(c(file_name,associated_entities)) 

naid_df <- data.frame()
for (i in 1:nrow(metadata)){
  naid_df[i,1] <- metadata_id$file_name[i]
  naid_df[i,2] <- metadata_id$associated_entities[i][[1]]$entity_submitter_id
}

colnames(naid_df) <- c("filename","TCGA_id")

03.获取不同矩阵类型

#提取你想要的类型 counts或者fpkm值

myfread <- function(files){
  data = data.table::fread(paste0("data_in_one/",files))
  data = data[-seq(1,4),]
  data = data$fpkm_unstranded#博主这边是fpkm值，想获取counts矩阵$后边换成unstranded
}

files <- dir("data_in_one")
system.time(f <- lapply(files,myfread))
expr_df <- as.data.frame(do.call(cbind,f))

#查看expr_df，发现我们的矩阵已经获取成功了

04.添加样品TCG编号和基因名

#但是发现样品名称和基因symbol没有，需要进一步处理

#先获取样品名称

rownames(naid_df) <- naid_df[,1]
naid_df <- naid_df[files,]
colnames(expr_df) <- naid_df$TCGA_id
 

#在获取基因名称

gene_id <- data.table::fread(paste0("data_in_one/",files[1]))$gene_name#这里是基因symbol如果想看到ENS开头的名称，$gene_id就可
gene_id <- gene_id[-seq(1,4)]
expr_df <- cbind(gene_id=gene_id,expr_df)

#结果是这样的基因symbol结果

#基因ENS结果是这样的

05.去除重复基因名

#如果做的是ENSG下边就不用看了，输出结果就行，如果是symbol需要处理一下重复基因名，可参考博主 R去除重复基因名博客

matrix0 <- order(rowMeans(expr_df[,-1]),decreasing = T)
expr_ordered=expr_df[matrix0,]
keep=!duplicated(expr_ordered$gene_id)
expr_max=expr_ordered[keep,]

06.输出

#最终结果输出

expr_max=as.matrix(expr_max)
rownames(expr_max) <- expr_max[,1]
exp=expr_max[,2:ncol(expr_max)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
lea=rbind(id=colnames(data),data)
write.table(lea,file = "TCGA.txt",sep = "\t",quote = F,col.names = F)
 

到这里矩阵中无重复基因名，用于接下来分析，今天的分享就到这里啦，博主抽空再把获取临床数据，SNP数据和CNV数据获取更新一下，至于时间嘛就待定了，欲知后事如何，请待下回分享~