问题来源:
探寻:NCR2 表达对GBM肿瘤芯片数据中cell-cycle 和 cytokines 通路基因的相关关系
- 通过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']]
- 下载这些基因在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
- 画出散点图
#定义一个画散点图的函数
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)