查看感兴趣基因的表达量及其相关关系

问题来源：
探寻：NCR2 表达对GBM肿瘤芯片数据中cell-cycle 和 cytokines 通路基因的相关关系

1. 通过KEGG.db 来下载KEGG相关通路的基因

if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("KEGG.db", version = "3.8")

library(KEGG.db)

去官网查看细胞周期相关的通路：
http://www.genome.jp/kegg-bin/show_pathway?hsa04110
细胞因子相关的通路：
http://www.genome.jp/kegg-bin/show_pathway?hsa04060
用KEGG来下载相应的基因列表：

cellcycle_genes=KEGGPATHID2EXTID[['hsa04110']]
cytokine_genes=KEGGPATHID2EXTID[['hsa04060']]

2. 下载这些基因在GBM芯片数据中的表达量

library(cgdsr)
#library(DT)
mycgds = CGDS("http://www.cbioportal.org/")   
all_TCGA_studies <- getCancerStudies(mycgds)  #cbioportal 中相关的研究
all_TCGA_studies[grepl('gbm',all_TCGA_studies$cancer_study_id),]  #找出与gbm研究相关的研究
all_tables <- getCaseLists(mycgds, "gbm_tcga")
all_dataset <- getGeneticProfiles(mycgds, "gbm_tcga")

#最后确定要
my_dataset <- 'gbm_tcga_mrna_U133'
my_table <- "gbm_tcga_mrna_U133" 
#下载那两个基因list表达图谱
cellcycle_expr <- getProfileData(mycgds, cellcycle_genes, my_dataset, my_table)
cytokine_expr <- getProfileData(mycgds, cytokine_genes, my_dataset, my_table)
#下载单个基因的表达向量
NCR2_expr <- getProfileData(mycgds, 'NCR2', my_dataset, my_table)

关于getProfileData 函数的说明：
mycgds : A CGDS objects
cellcycle_genes: 基因名字的一个向量
mydataset： genetic profile ID 的向量 或者单个 genetic profile 字符串
前三个参数都是必须参数
my_table : A case list ID

3. 画出散点图

#定义一个画散点图的函数
cor_plot <- function(x,y){
  plot(x,y, xlab = 'NCR2', ylab = 'gene')
  model = lm(y ~ x)
  summary(model) 
  int =  model$coefficient["(Intercept)"]
  slope =model$coefficient["x"]
  abline(int, slope,
         lty=1, lwd=2, 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))
}

#执行函数画出散点图
cor_plot(x=NCR2_expr$NCR2,y=cellcycle_expr$CCNB1)
cor_plot(x=NCR2_expr$NCR2,y=cellcycle_expr$BUB1B)
cor_plot(x=NCR2_expr$NCR2,y=cytokine_expr$CCL4)