github-actions[bot]
commented
2 years ago 单基因生信分析--差异分析&生存分析 by 医学僧的科研日记
写在前边
单基因功能试验仍然是小白入门的主流，但很多人往往都会遇到无生存信息，随访又难的问题，老老实实随访会搭上大部分时间精力；鳞次栉比的网页工具也能实现预后曲线，但其形式单一，不同网站差异很大。通过公众号分享，代码能轻松解决以上问题，但对菜鸟来讲不太容易领悟。这里，小王子通过自己的理解，总结出以下缩减版的单基因差异表达以及生存分析模块，希望对你有所帮助！
01
文件准备
如上图，我们需要从UCSC Xena（https://xenabrowser.net/datapages/）网站下载所需的表达、注释以及生存文件（以下以LIHC文件为例），然后用以下代码进行预处理：
设置工作目录library(data.table)library(magrittr)expr <- fread('TCGA-LIHC.htseq_fpkm.tsv.gz',h=T,check.names = F)ann <- read.table('gencode.v22.annotation.gene.probeMap',h=T)expr2 <- merge(ann[,1:2],expr,by=1) %>% .[,-1]library(limma)expr2 <- avereps(expr2[,-1],ID=expr2$gene) %>% as.data.frame() ##如果有重复的行，取均值##这里得到的expr2中既有mRNA,也有lncRNA和miRNA，所以我们要将其区分开来library(biomaRt)ensembl <- useMart(biomart = "ENSEMBL_MART_ENSEMBL",dataset = "hsapiens_gene_ensembl",host = "www.ensembl.org")listAttributes(ensembl)[63,]feature_info <- getBM(attributes = c("gene_biotype","hgnc_symbol"),filters = "hgnc_symbol",values = rownames(expr2), mart = ensembl)table(feature_info$gene_biotype)write.csv(feature_info,'annotation.csv',row.names = F)mRNA <- expr2[feature_info$hgnc_symbol[feature_info$gene_biotype=='protein_coding'],]lncRNA <- expr2[feature_info$hgnc_symbol[feature_info$gene_biotype=='lncRNA'],]miRNA <- expr2[feature_info$hgnc_symbol[feature_info$gene_biotype=='miRNA'],]save(mRNA,lncRNA,miRNA,file='expression.Rdata')
生存文件预处理：
rm(list = ls())sur <- read.table('TCGA-LIHC.survival.tsv.gz',h=T,check.names = F)##这里的OS.time是以天为单位的，转化为以年为单位sur$OS.time <- sur$OS.time/365sur <- sur[,-3]save(sur,file = 'survival.Rdata')
02
差异分析
接下来，我们对肿瘤和正常样本中单基因的表达量进行差异分析（详细内容可以参考我们往期内容“为图形添加p值”）
rm(list = ls())load('expression.Rdata')library(magrittr)gene <- mRNA["AASS",]%>%t()%>%as.data.frame() ##AASS为我们的目的基因gene$type <- sapply(gene,function(x){ifelse(substr(rownames(gene),14,15)!=11, "T","N")})%>%as.factor()table(gene$type)##癌和正常组织已区分好，接下来按照往期推送“为图形添加p值”开始画图library(ggplot2)library(ggpubr)ggplot(aes(type,AASS,color=type),data=gene)+ ###横轴跟纵轴的标签geom_dotplot(binaxis="y", stackdir='centerwhole', method="dotdensity",binpositions="all",stackgroups = F,stackratio=.5,dotsize = .4)+stat_compare_means(method= "t.test",label.x = 1.4, label.y= 7)+ ###p值的位置，根据需要可调stat_boxplot(geom = "errorbar",width=.2)+ ###误差线宽度，根据需要可调geom_boxplot(alpha=0, width=0.6)+ ###中间箱图宽度 ，根据需要可调theme_classic()
03
预后分析
接下来，我们根据前期文章“R-生存曲线三步走”以及“R-生存分析一行代码直达”两篇推送，对LIHC进行生存分析（也可以利用GEPIA、Kaplan-Meier Plotter等网页工具进行生存分析，不过囿于不同网站结果差异很大，以及无法根据需要进行修饰，推荐利用R进行生存分析）
rm(list = ls())library(survival)library(survminer)load('expression.Rdata')load('survival.Rdata')gene <- 'AASS' ###目的基因data <- mRNA[gene,]surdata <- merge(sur,t(data),by.x = 1,by.y = 0)cut <- surv_cutpoint(surdata,'OS.time','OS','AASS') ###需要更换目的基因print(paste0('The optimal cutpoint is ',cut$cutpoint$cutpoint,'.'))surdata$level<- ifelse(surdata$AASS>cut$cutpoint$cutpoint,'High','Low') ###需要更换目的基因fit <- survfit(Surv(OS.time,OS)~level,data = surdata)surv_pvalue(fit)$pvalggsurvplot(fit,pval = T,pval.method = T,risk.table = T)
2.方法二：根据中位值分组
rm(list = ls())library(survival)library(survminer)load('expression.Rdata')load('survival.Rdata')gene <- 'AASS' ###目的基因data <- mRNA[gene,]surdata <- merge(sur,t(data),by.x = 1,by.y = 0)surdata$level <- ifelse(surdata[,4]>median(surdata[,4]),'High','Low')fit <- survfit(Surv(OS.time,OS)~level,data = surdata)surv_pvalue(fit)$pvalggsurvplot(fit,pval = T,pval.method = T,risk.table = T)
rm(list = ls())library(survival)library(survminer)load('expression.Rdata')load('survival.Rdata')gene <- 'AASS' ###目的基因data <- mRNA[gene,]surdata <- merge(sur,t(data),by.x = 1,by.y = 0)surdata$level <- ifelse(surdata[,4]>mean(surdata[,4]),'High','Low')fit <- survfit(Surv(OS.time,OS)~level,data = surdata)surv_pvalue(fit)$pvalggsurvplot(fit,pval = T,pval.method = T,risk.table = T)
写在后边
总的来说，生存分析的三种方式不同的分组方式会得到不同的效果，最直接的反应就是p值的不同，就LIHC中的AASS基因来讲，中位值和平均值差异均不显著，但最佳截断点却可以！这就给我们的小论文带来了曙光。另外，等R水平从菜鸟进阶成功后还可以根据不同比例分组进行预后分析！总之，只要基因有意义，p值总会有意义的，具体的方法可以参考前期推送“R-生存曲线三步走”以及“R-生存分析一行代码直达”。
https://mp.weixin.qq.com/s/7EdwrhtsfrgKYoTAL8_WrA