- 加载包
library(TCGAbiolinks)
library(dplyr)
library(DT)
library(SummarizedExperiment)
- 下载临床数据
clinical_data <- GDCquery_clinic(project = "TCGA-LIHC", type = "clinical")
- 下载表达谱数据
Expr_df <- GDCquery(project = "TCGA-LIHC",
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
workflow.type = "HTSeq - FPKM")
GDCdownload(Expr_df, method = "api", files.per.chunk = 100)
expdat <- GDCprepare(query = Expr_df)
Expr_matrix=assay(expdat)
但此时得到的是ensemblID，需要转成我们熟悉的symbol ID
- 基因ID的转换
library(clusterProfiler)
gene<-bitr(rownames(Expr_matrix),"ENSEMBL","SYMBOL","org.Hs.eg.db")
Expr_matrix<-cbind(rownames(Expr_matrix),Expr_matrix)
colnames(Expr_matrix)[1]<-"ENSEMBL"
df<-merge(gene,Expr_matrix,by="ENSEMBL")
- 整理分组信息
group <- strsplit(colnames(df_50)[-1],"[-]")
class<-sapply(group,function(I){I[4]})
control<-grepl("11",class)
control<-which(control==TRUE)
class[control]<-"normal"
class[-control]<-"cancer"
- 整理TCGAbiolinks临床数据
整体思路如下：
- clinical data是submitter_id，只有三个字段，即“TCGA-3Z-A93Z”，而表达谱数据有7个"TCGA-3Z-A93Z-01A-01R-0864-07"，需要把表达谱数据ID拆分成3个字段与submitter_id对上
- vital_status为生存分析中的OS，将Alive和Dead改成0和1
- 关于OS.time: 当患者dead时，days_to_death为OS.time；当患者alive时，days_to_last_follow_up为OS.time
- 此外还有性别，年龄，种族等临床信息，大家自行选择
#拆分TCGA表达谱7个字段的ID，并组成3个
newid<-lapply(strsplit(rownames(df_50),'-'),function(i){paste0(i[1:3],collapse = '-')})
newid<-sapply(1:611,function(i){newid[[i]]})
df_50<-cbind(df_50,newid)
colnames(clinical_data1)[1]<-'newid'
df_OS<-merge(clinical_data1,df_50,by="newid")
df_OS$OS.time<-df_OS$OS.time/365
rownames(df_OS)<-df_OS[,1]
#去重，ID存在一对多的情况，因为懒我就robust的直接删掉了，你可以自己进行高标准选择
df_OS_dropdu<-df_OS[-which(duplicated(df_OS[,1])),]
realdata<-df_OS_dropdu
rownames(realdata)<-realdata[,1]
realdata<-realdata[,-1]
