差异分析生存分析热图绘制基于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()