获取表达矩阵以及绘制生存曲线

继续加油！

Sys.setenv(LANGUAGE = 'en')
options(stringsAsFactors = FALSE)
rm(list=ls())
library(tidyverse)
library(GEOquery)
library(dplyr)
##在GEO中下载GSE数据，以GSE12417为例
#普遍下载
gset <- getGEO('GSE12417',destdir = '.',getGPL = F)#有时候网速不好，存在数据下载不完整的情况
#换种方式下载
remotes::install_github("jmzeng1314/GEOmirror")
eSet=geoChina('GSE12417')
#还是不行的话 再换吧 哈哈
library(remotes)
url='https://gitee.com/jmzeng/annoprobe.git'
install_git(url)
gset=geoChina('GSE12417')  
class(gset) #"list" 8.4M

gset.png

gset <- gset[[2]]
class(gset)  #ExpressionSet
#提取表达矩阵和临床信息都是在gset为ExpressionSet基础上进行的
expr <- exprs(gset) #提取表达矩阵
exp <- as.data.frame(expr) #将表达矩阵转换为数据框，方便操作
##探针和id的转换
library(data.table)
anno <- fread('GPL96-57554.txt',sep = '\t',header = T,data.table = F)
gene_symbol <- anno$`Gene Symbol`
a <- strsplit(gene_symbol,split=' /// ',fixed = T) #list
a1 <- sapply(a, function(x){x[1]})
##将id和symbol放在一个数据框内
z <- data.frame(anno$ID,a1)
##合并z和表达矩阵
exp1 <- merge(z,exp,by.x = 1,by.y=0)
##将a1作为exp1的行名
rownames(exp1) <- exp1$a1 ##存在重复值
exp2 <- distinct(exp1,a1,.keep_all = T)#去重复
rownames(exp2) <- exp2$a1 #存在缺失值
exp3 <- na.omit(exp2)
rownames(exp3) <- exp3$a1
exp4 <- exp3[,-c(1,2)]  ##最终的表达矩阵
##绘制生存曲线
#提取临床信息
pd <- pData(gset)
p <- pd$characteristics_ch1  #提取含有生存时间和生存状态的列
head(p)
p1 <- strsplit(p,split='; ',fixed = T)
os.time1 <- sapply(p1, function(x){x[3]})
os.time2 <- strsplit(os.time1,split = ' ',fixed = T)
os.time3 <- sapply(os.time2, function(x){x[3]})          
os.time <- as.numeric(os.time3)
status1 <- sapply(p1, function(x){x[4]})
head(status1)
status2 <- strsplit(status1,split=': ',fixed = T)
status3 <- sapply(status2, function(x){x[2]})
status <- as.numeric(status3)
#将os.time和status合并为数据框
x <- data.frame(os.time,status)
rownames(x) <- rownames(pd)
#选取某一个基因进行绘制
exp4.e <- exp4['EZH2',]
exp4.e <- t(exp4.e) #转置
exp4.e <- as.data.frame(exp4.e)
class(exp4.e) #data.frame
##将exp4.e和x组合（含os.time和status）
exp4.e.x <- merge(exp4.e,x,by=0)
exp4.e.x$EZH2<- ifelse(exp4.e.x$EZH2>median(exp4.e.x$EZH2),'high','low')
library(survminer)
library(survival)
fit <- survfit(Surv(os.time,status)~EZH2,data = exp4.e.x)
ggsurvplot(fit,conf.int = F,pval = T,risk.table = F,
           linetype = 'strata')

ggsurvplot(fit,conf.int = F,pval = T,risk.table = T,
           risk.table.col='strata',
           linetype = c("solid", "dashed"),
           surv.median.line = 'hv',
           palette = c('red','blue'))
?ggsurvplot

生存曲线.png