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"))
参考资料
- TCGA的28篇教程- 指定癌症查看感兴趣基因的表达量
- 【r<-统计|绘图】使用R进行生存分析——一文打尽https://www.jianshu.com/p/4ad9ba730719
请关注“恒诺新知”微信公众号,感谢“R语言“,”数据那些事儿“,”老俊俊的生信笔记“,”冷🈚️思“,“珞珈R”,“生信星球”的支持!