TCGA之获取指定癌症的感兴趣基因的表达量

本文使用cdgsr包获取感兴趣基因的mRNA表达量，并绘制相关关系图和生存曲线。相关资料来源于一篇Cell文献，原文章是Natural Killer Cells Control Tumor Growth by Sensing a Growth Factor（PMID:29275861）。

原文相关图示如下，在GBM（胶质细胞瘤）中，NCR2的表达同细胞因子相关基因、细胞循环基因呈现相关关系，NCR2相关的两个生存曲线。

首先是载入需要的R包

# 载入必要的包
library("cgdsr")
library("KEGG.db")
library("tidyverse")
library("survival")
library("survminer")

细胞循环通路和细胞因子通路的所有基因ID可以通过KEGG.db包获取（查询资料得知KEGG.db包已经很久没人维护了，推荐使用KEGGREST，这里仍然使用KEGG.db）。

#ls("package:KEGG.db")
cellcycle_genes=KEGGPATHID2EXTID[['hsa04110']] #细胞循环通路基因
cytokine_genes=KEGGPATHID2EXTID[['hsa04060']] #细胞因子通路基因

使用cgdsr包获取GBM中相应的样本中的mRNA表达量数据：

###获取GBM（神经胶质瘤）的表达量数据###
# 创建CGDS对象
mycgds <- CGDS("http://www.cbioportal.org/")
cancerstudies <- getCancerStudies(mycgds)%>% filter(grepl("gbm",.[,1]))%>%select(1)
cancerstudy <- cancerstudies[5,] # 数据集是"gbm_tcga"

# 获取样本列表gbm_tcga_mrna_U133
caseL <- getCaseLists(mycgds,cancerstudy)%>%select(1:4) # 选择mRNA表达数据中样本量最大的一组样本
caselist <- caseL[9,1] #gbm_tcga_mrna_U133
# 要获取的数据类型：mRNA数据
geneticP <- getGeneticProfiles(mycgds,cancerstudy) #U133
geneticprofilter <- geneticP[4,1]

# 根据样本类型和数据类型，使用getProfileData函数获取相应基因的表达量数据
geneticDatacellcycle <- getProfileData(mycgds,cellcycle_genes,geneticprofilter,caselist)
geneticDatacytokine <- getProfileData(mycgds,cytokine_genes,geneticprofilter,caselist)
NCR2_exp <- getProfileData(mycgds,"NCR2",geneticprofilter,caselist)

绘制相关性点图

#绘图函数,基本上全是Jimmy大神写的函数，多加了一个拟合曲线方程
cov_plot <- function(x,y){
  plot(x,y,xlab="NCR2",ylab="Gene")
  model<-lm(y~x)
  summary(model)
  intercept <- format(model$coefficients[1], digits = 4)
  slop <- format(model$coefficients[2], digits = 4)
  abline(intercept,slop, lty=1, lwd=2, col="red")
  label <- paste0("y=",as.numeric(slop),"x+",as.numeric(intercept))
  text(0.95*max(x), 0.95*max(y), label, col = "red")
  r= format(cor(x,y),digits = 4)
  p= format(cor.test(x,y)$p.value,digits = 4)
  title(main = paste0('p value=',p), 
        sub  = paste0('correlation=',r))
}
#绘制NCR2和细胞循环通路第一个基因ANAPC1之间的相关关系
cov_plot(NCR2_exp$NCR2,geneticDatacellcycle$ANAPC1)

结果如下，添加了拟合曲线方程，其他基本上是一样的。

其他基因也是类似的，调整相应的基因名字即可，如果要输出所有的图，可以写一个循环。

#"ANAPC1"  "ANAPC10" "ANAPC11" "ANAPC13"等基因
cov_plot(NCR2_exp$NCR2,geneticDatacellcycle$ANAPC1)
cov_plot(NCR2_exp$NCR2,geneticDatacellcycle$ANAPC10)
cov_plot(NCR2_exp$NCR2,geneticDatacellcycle$ANAPC11)
cov_plot(NCR2_exp$NCR2,geneticDatacellcycle$ANAPC13)

原文还有两个个图（I、J）是生存曲线，但是我并没有看懂这两个图是什么意思，下面的生存曲线是NCR2高表达（Q4）与低表达（Q1）对GBM病人总生存率的关系。

#library("survival") 生存分析需要survival包
#首先获取临床数据，并取出OS_MONTHS、OS_STATUS两列
#再根据NCR2的表达量，将其分为"High" "Medium" "Low"
#最后合并两个表格，只取NCR2的值是"High"和"Low"
NCR2_clinical <- getClinicalData(mycgds, caselist)%>%select(starts_with("OS")) %>%
  transmute(time=OS_MONTHS*30,status=ifelse(OS_STATUS == "LiVING",0,1)) 
NCR2_HL <- NCR2_exp %>% 
  mutate(NCR2=ifelse(NCR2>quantile(NCR2)[2],ifelse(NCR2>quantile(NCR2)[4],"High","Medium"),"Low"))
NCR2_clinical_data <- cbind(NCR2_clinical,NCR2_HL) %>%filter(NCR2_HL!="Medium")


#KM生存分析，并绘图即可
with(NCR2_clinical_data,{
  NCR2_km <<-survfit(Surv(time,status)~NCR2)
})

#library("survminer")
ggsurvplot(NCR2_km,data = NCR2_clinical_data, pval = T,
           title="Kaplan-Meier Curve for NCR2 expression", 
           legend.title="NCR2", legend.labs=c("High", "Low"))