- 加載包
library(TCGAbiolinks)
library(dplyr)
library(DT)
library(SummarizedExperiment)
- 下載臨床數(shù)據(jù)
clinical_data <- GDCquery_clinic(project = "TCGA-LIHC", type = "clinical")
- 下載表達(dá)譜數(shù)據(jù)
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)
但此時(shí)得到的是ensemblID扯躺，需要轉(zhuǎn)成我們熟悉的symbol ID
- 基因ID的轉(zhuǎn)換
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臨床數(shù)據(jù)
整體思路如下：
- clinical data是submitter_id脂崔，只有三個(gè)字段巢掺，即“TCGA-3Z-A93Z”，而表達(dá)譜數(shù)據(jù)有7個(gè)"TCGA-3Z-A93Z-01A-01R-0864-07"胖翰，需要把表達(dá)譜數(shù)據(jù)ID拆分成3個(gè)字段與submitter_id對(duì)上
- vital_status為生存分析中的OS，將Alive和Dead改成0和1
- 關(guān)于OS.time: 當(dāng)患者dead時(shí)壕曼，days_to_death為OS.time琳拨；當(dāng)患者alive時(shí)，days_to_last_follow_up為OS.time
- 此外還有性別弛随，年齡瓢喉，種族等臨床信息，大家自行選擇
#拆分TCGA表達(dá)譜7個(gè)字段的ID舀透，并組成3個(gè)
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存在一對(duì)多的情況，因?yàn)閼形揖蛂obust的直接刪掉了盐杂，你可以自己進(jìn)行高標(biāo)準(zhǔn)選擇
df_OS_dropdu<-df_OS[-which(duplicated(df_OS[,1])),]
realdata<-df_OS_dropdu
rownames(realdata)<-realdata[,1]
realdata<-realdata[,-1]
