学徒和学员已经陆续出师,是时候把生信技能树的舞台交给后辈了! 呃提笔来写这个总结,貌似是上周发生的事了,但是不太想回忆,学习新东西会抱着满满的热情,但是翻旧账就会让人心里不是很舒畅...特别是没能完成的旧事重提。唔但是也算是慢慢习惯了完成一件事会记下点个人感悟,nice~这是个好习惯(最开始当然也只是Jimmy大大的任务模式啦)。 这里感谢老师费心帮助我完成这个,还找了位前辈小蚂蚁山神(就是昨天完成任务的那个:学徒数据挖掘之IgG4-RD患者的WGCNA分析)来帮助我,非常感谢大家的帮助:happy: 写完这个,又发生了一件特别让我感动的事情,不管多小的人物类似于我,Jimmy大大都会认真的回答你的小小疑问,遇到一个好老师真的很重要,超开心能认识这么棒的老师和一些志同道合的小伙伴。虽然学习生信的初衷可能只是为了发一篇简单的小文章好找工作一点,但是不知不觉中每天学习一点生信技能树的公众号,论坛和各种这个团队优秀伙伴的简书笔记成为了一个习惯,习惯着等任务,没有任务时自己找事情做,B站不再是天天刷吃播视频哈哈哈,B站真的能学习,老师的那么多视频得刷到啥时候啊,还不断更新~收藏的东西也不会吃灰:six::six::six:。害,就是我太懒了,不愿意拓展知识。 首先接到完成文章中一个生存分析的图的任务: 数据文章中说明使用了cbioportal上的TCGA的数据:
原文中的图为这样:
基本术语:
用在线xena下载数据,直接下载临床信息,全部都是整理好的,分14个数据集的和19个数据集的,19的那个。 xena数据库中会有两个数据:
有了这么多数据文件,接下来就是该秀出我在生信技能树学习班获得的R语言知识啦: rm(list=ls()) options(stringsAsFactors = F) ##读取数据 mRNA_HiSeqV2 = read.table('HiSeqV2.gz',header = T,sep = '\t',check.names = F,row.names = 1) dim(mRNA_HiSeqV2) mRNA_HiSeqV2[1:4,1:4]
#查看NA的数据 dim(mRNA_HiSeqV2) exp = mRNA_HiSeqV2 exp = exp[apply(exp, 1, function(x) sum(x > 1) > 9), ] dim(exp) exp[1:4,1:4]
##临床信息(在这里没有用到) mRNA_clinical = read.table('LIHC_clinicalMatrix',header = T,sep = '\t',row.names = 1) dim(mRNA_clinical) mRNA_clinical[1:4,1:4]
##生存相关信息 mRNA_survival = read.table('LIHC_survival.txt.gz',header = T,sep = '\t',row.names = 1) dim(mRNA_survival) mRNA_survival <- mRNA_survival[, -1] mRNA_survival[1:4,1:4]
save(mRNA_survival,mRNA_clinical,exp,group_list,file='LIHC.Rdata')
#01–09是癌症,10–19是正常,20–29是癌旁 根据样本ID的第14-15位,给样本分组(tumor和normal) library(stringr) table(str_sub(colnames(exp),14,15)) group_list = ifelse(as.numeric(str_sub(colnames(exp),14,15)) < 10,'tumor','normal') group_list = factor(group_list,levels = c("normal","tumor")) table(group_list)
##xena下载的数据是经过log的count值,所以要还原才可以进行差异分析 exprSet = exp exprSet <- 2^exprSet-1 dim(exprSet) exprSet[1:4,1:4]
#并且还要取整数,现在就变成count值了 exprSet <- floor(exprSet) exprSet[1:4,1:4]
save(mRNA_survival,mRNA_clinical,exp,group_list,exprSet,file='LIHC1.Rdata')
rm(list=ls()) options(stringsAsFactors = F) load('LIHC1.Rdata') library("stringr")
先看以count数据的结果 OS tumor_sample <- colnames(exprSet)[substr( colnames(exprSet),14,15) < 10] mRNA_survival = mRNA_survival[,-1] pheno = rownames(mRNA_survival)[substr(rownames(mRNA_survival),14,15) < 10] fin_tumor = intersect(tumor_sample,pheno)
# 提取表达矩阵和临床信息 exprSet = exprSet[,fin_tumor] pheno_fina = mRNA_survival[fin_tumor,]
### 这里又归一化了? exprSet = log2(exprSet+1)
代码超级长,绝大部分都是生信技能树jimmy老师写的,我只是一个代码搬运工以及修修补补。
# 开始生存分析 library(ggrisk) library(survival) library(survminer)
# 先做无病生存期,首先把DFI不明(NA)的样本给去掉 DFI = pheno_fina[,3:6] dfi = DFI[!is.na(DFI$DFI),] dfi_count = exprSet[,rownames(dfi)] dd = scale(as.numeric( dfi_count["LHPP",]))
# 绘图 mySurv<-Surv(dfi$DFI.time, dfi$DFI) ### 以LHPP的表达量z-score后,<= -1为界,哦z-score是这样用的啊 dfi_group<-ifelse(dd <= -1 ,"low","high") fi <-survfit(mySurv~dfi_group,data = dfi) survminer::ggsurvplot(fi, data = dfi, risk.table = TRUE, pval = T, )
# 接下来绘制os的 os = pheno_fina[,c(1,2,7,8)] os = os[!is.na(os$OS),] os_count = exprSet[,rownames(os)]
dd = scale(as.numeric( os_count["LHPP",]))
mySurv<-Surv(os$OS.time, os$OS) dfi_group<-ifelse(dd <= -1 ,"low","high") fi <-survfit(mySurv~dfi_group,data = dfi) survminer::ggsurvplot(fi, data = dfi, risk.table = TRUE, pval = T, )
这大概算是完成了
原文没有指明,他用的表达数据是Count,标准化后的Count,FPKM还是TPM,也 没有说是否对数据进行Log2(dat + 1)处理,不过表达量本身仅仅是用来把病人分组,何种归一化并不会影响基因表达量高低这个排序。 写到最后如果你时间比较充裕,而且确实想把生存分析学会,学好,我在生信技能树多次分享过生存分析的细节也值得你认真读完;
“ID转换和生存分析”群的钉钉群号: 35371384,如果你下载钉钉软件,搜索进群这个能力都没有,我还是建议你放弃学生信哈! 如果你也想开启自己的生物信息学数据处理生涯,加入生信技能树小圈子,还等什么呢,赶快行动起来吧! 发邮件(jmzeng1314@163.com)给生信技能树创始人jimmy就有惊喜哦!当然了,不能是辣鸡或者骚扰邮件啦,带上自己的简历和想学习交流的诚心吧! |
|