知识的价值不在于占有,而在于使用。

生信自学网-速科生物-生物信息学数据库挖掘视频教程

当前位置: 主页 > ICGC >

差异分析生存分析热图绘制基于ICGC数据库

时间:2019-07-30 11:33来源:生信自学网 作者:乐伟 点击:
在这里,做差异用到的limma R包,大家对这个包应该非常熟悉了,只是还没用来做ICGC数据的差异。

得到ICGC表达矩阵之后,我们就可以做后面的差异分析了
在这里,做差异用到的limma R包,大家对这个包应该非常熟悉了,只是还没用来做ICGC数据的差异。
1、下面我们来看下做差异的代码:

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("limma", version = "3.8")

library("limma")

setwd("C:\\Users\\lexb4\\Desktop\\ICGCexp\\04.diff")              #设置工作目录
inputFile="sampleExp.txt"                                         #输入文件
fdrFilter=0.05                                                    #fdr临界值
logFCfilter=2                                                     #logFC临界值
conNum=202                                                        #normal组样品数目
treatNum=240                                                      #tumor组样品数目

#读取输入文件
outTab=data.frame()
grade=c(rep(1,conNum),rep(2,treatNum))
rt=read.table(inputFile,sep="\t",header=T,check.names=F)
rt=as.matrix(rt)
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data=avereps(data)
data=data[rowMeans(data)>0.5,]
更多原创代码请联系生信自学网微信:18520221056
关注微信公众号“biowolf_cn"

然后我们就得到了差异结果,在下面的结果表格中,第一列就是基因名称,conMean基因在正常组表达量的均值,treatMean基因在癌症表达了的均值;logFC大家应该非常熟悉,就是log2的肿瘤组比正常组,如果logFC大于零,说明该基因在肿瘤组是上调的;pValue统计检验的P值,fdr叫做后的P值。

2、接下来可以绘制一下热图,热图用到的输入文件就是我们做差异得到的diffGeneExp.txt,用到的包是pheatmap:

install.packages("pheatmap")

setwd("C:\\Users\\lexb4\\Desktop\\ICGCexp\\05.pheatmap")      #设置工作目录
rt=read.table("diffGeneExp.txt",sep="\t",header=T,row.names=1,check.names=F)
rt=log2(rt+0.001)
library(pheatmap)
Type=c(rep("con",202),rep("treat",243))    #修改对照和处理组样品数目
names(Type)=colnames(rt)
Type=as.data.frame(Type)
获取生信自学网原创代码,可以直接购买课程


3、生存分析
前面我们提到ICGC数据库部分癌症提供临床数据,那么我们可以用下载的临床数据,结合做差异得到的表达数据,整合成做生存的文件:

有了这样一个含有基因表达、生存时间、生存状态的表格,我们就可以做生存分析,生存分析用到的是survival R包,而且是批量做生存曲线,并且得到了每个基因的生存曲线P值,方便大家分析:

install.packages("survival")

setwd("C:\\Users\\lexb4\\Desktop\\ICGCexp\\07.survival")    #工作目录(需修改)
library(survival)
pFilter=0.05                                                #由于图形太多,只对p小于0.05的基因绘图
rt=read.table("survivalInput.txt",header=T,sep="\t",check.names=F)    #读取输入文件
rt$futime=rt$futime/365                                     #如果以月为单位,除以30;以年为单位,除以365
outTab=data.frame()



购买《ICGC数据库挖掘》课程



加生信自学网群
责任编辑:伏泽
作者申明:本文版权属于生信自学网(微信号:18520221056)未经授权,一律禁止转载!
BioWolf二维码生成器
------分隔线----------------------------
GEO芯片数据库挖掘生信视频教程
推荐内容
TCGA数据库挖掘文章套路生信视频教程
中药复方网络药理学文章套路生信视频教程