火山图
功能分析
library("org.Hs.eg.db")
rt=read.csv("diff.csv",sep=",",check.names=F,header=T)
rt<-rt[,1]
rt<-as.data.frame(rt)
genes=as.vector(rt[,1])
entrezIDs <- mget(genes, org.Hs.egSYMBOL2EG, ifnotfound=NA)
entrezIDs <- as.character(entrezIDs)
out=cbind(rt,entrezID=entrezIDs)
write.table(out,file="id.txt",sep="\t",quote=F,row.names=F)
###############GO分析
library("clusterProfiler")
library("org.Hs.eg.db")
library("enrichplot")
library("ggplot2")
rt=read.table("id.txt",sep="\t",header=T,check.names=F)
rt=rt[is.na(rt[,"entrezID"])==F,]
cor=rt$cor
gene=rt$entrezID
names(cor)=gene
#GO富集分析
kk <- enrichGO(gene = gene,
OrgDb = org.Hs.eg.db,
pvalueCutoff =0.05,
qvalueCutoff = 0.05,
ont="all",
readable =T)
write.table(kk,file="GO.txt",sep="\t",quote=F,row.names = F)
#柱状图
tiff(file="GO.tiff",width = 26,height = 20,units ="cm",compression="lzw",bg="white",res=600)
barplot(kk, drop = TRUE, showCategory =10,split="ONTOLOGY") + facet_grid(ONTOLOGY~., scale='free')
dev.off()
#####################KEGG分析
library("clusterProfiler")
library("org.Hs.eg.db")
library("enrichplot")
library("ggplot2")
setwd("C:\\Users\\lexb4\\Desktop\\CCLE\\09.KEGG")
rt=read.table("id.txt",sep="\t",header=T,check.names=F)
rt=rt[is.na(rt[,"entrezID"])==F,]
cor=rt$cor
gene=rt$entrezID
names(cor)=gene
#kegg富集分析
kk <- enrichKEGG(gene = gene, organism = "hsa", pvalueCutoff =0.05, qvalueCutoff =0.05)
write.table(kk,file="KEGG.txt",sep="\t",quote=F,row.names = F)
#柱状图
tiff(file="barplot.tiff",width = 20,height = 12,units ="cm",compression="lzw",bg="white",res=600)
barplot(kk, drop = TRUE, showCategory = 20)
dev.off()
6.3 批量KM曲线
输入数据KM.TXT，数据是FPKM.
KM.txt
library(survival)
rt=read.table("KM.txt",header=T,sep="\t",check.names=F)
outTab=data.frame()
for(gene in colnames(rt[,4:ncol(rt)])){
a=rt[,gene]<=median(rt[,gene])
diff=survdiff(Surv(futime, fustat) ~a,data = rt)
pValue=1-pchisq(diff$chisq,df=1)
outTab=rbind(outTab,cbind(gene=gene,pvalue=pValue))
fit <- survfit(Surv(futime, fustat) ~ a, data = rt)
summary(fit)
#只将p<0.05的基因画出来。
if(pValue<0.05){
if(pValue<0.001){
pValue="<0.001"
}else{
pValue=round(pValue,3)
pValue=paste0("=",pValue)
}
pdf(file=paste(gene,".survival.pdf",sep=""),
width=6,
height=6)
plot(fit,
lwd=2,
col=c("red","blue"),
xlab="Time (year)",
mark.time=T, #显示删失数据
ylab="Survival rate",
main=paste(gene,"(p", pValue ,")",sep=""))
legend("topright",
c("High","Low"),
lwd=2,
col=c("red","blue"))
dev.off()
}
}
write.table(outTab,file="survival.xls",sep="\t",row.names=F,quote=F)
6.4 差异基因进行COX模型构建
用FPKM数据做COX模型的构建，用OS相关的基因做分析。
输入数据exptime.txt
exptime.txt
单因素cox
pFilter=0.01 #定义单因素显著性
library(survival) #引用包
rt=read.table("expTime.txt",header=T,sep="\t",check.names=F,row.names=1) #读取输入文件
sigGenes=c("futime","fustat")
outTab=data.frame()
for(i in colnames(rt[,3:ncol(rt)])){
cox <- coxph(Surv(futime, fustat) ~ rt[,i], data = rt)
coxSummary = summary(cox)
coxP=coxSummary$coefficients[,"Pr(>|z|)"]
outTab=rbind(outTab,
cbind(id=i,
HR=coxSummary$conf.int[,"exp(coef)"],
HR.95L=coxSummary$conf.int[,"lower .95"],
HR.95H=coxSummary$conf.int[,"upper .95"],
pvalue=coxSummary$coefficients[,"Pr(>|z|)"])
)
if(coxP<pFilter){
sigGenes=c(sigGenes,i)
}
}
write.table(outTab,file="uniCox.xls",sep="\t",row.names=F,quote=F)
uniSigExp=rt[,sigGenes]
uniSigExp=cbind(id=row.names(uniSigExp),uniSigExp)
write.table(uniSigExp,file="uniSigExp.txt",sep="\t",row.names=F,quote=F)
LASSO回归
library("glmnet")
library("survival")
rt=read.table("uniSigExp.txt",header=T,sep="\t",row.names=1,check.names=F) #读取文件
rt$futime[rt$futime<=0]=1
x=as.matrix(rt[,c(3:ncol(rt))])
y=data.matrix(Surv(rt$futime,rt$fustat))
fit <- glmnet(x, y, family = "cox", maxit = 1000)
pdf("lambda.pdf")
plot(fit, xvar = "lambda", label = TRUE)
dev.off()
cvfit <- cv.glmnet(x, y, family="cox", maxit = 1000)
pdf("cvfit.pdf")
plot(cvfit)
abline(v=log(c(cvfit$lambda.min,cvfit$lambda.1se)),lty="dashed")
dev.off()
coef <- coef(fit, s = cvfit$lambda.min)
index <- which(coef != 0)
actCoef <- coef[index]
lassoGene=row.names(coef)[index]
lassoGene=c("futime","fustat",lassoGene)
lassoSigExp=rt[,lassoGene]
lassoSigExp=cbind(id=row.names(lassoSigExp),lassoSigExp)
write.table(lassoSigExp,file="lassoSigExp.txt",sep="\t",row.names=F,quote=F)
多因素COX模型的构建
library(survival)
library(survminer)
rt=read.table("lassoSigExp.txt",header=T,sep="\t",check.names=F,row.names=1) #读取train输入文件
rt[,"futime"]=rt[,"futime"]/365
#使用train组构建COX模型
multiCox=coxph(Surv(futime, fustat) ~ ., data = rt)
multiCox=step(multiCox,direction = "both")
multiCoxSum=summary(multiCox)
#输出模型相关信息
outTab=data.frame()
outTab=cbind(
coef=multiCoxSum$coefficients[,"coef"],
HR=multiCoxSum$conf.int[,"exp(coef)"],
HR.95L=multiCoxSum$conf.int[,"lower .95"],
HR.95H=multiCoxSum$conf.int[,"upper .95"],
pvalue=multiCoxSum$coefficients[,"Pr(>|z|)"])
outTab=cbind(id=row.names(outTab),outTab)
write.table(outTab,file="multiCox.xls",sep="\t",row.names=F,quote=F)
#绘制森林图
pdf(file="forest.pdf",
width = 8, #图片的宽度
height = 5, #图片的高度
)
ggforest(multiCox,
main = "Hazard ratio",
cpositions = c(0.02,0.22, 0.4),
fontsize = 0.7,
refLabel = "reference",
noDigits = 2)
dev.off()
#输出train组风险文件
riskScore=predict(multiCox,type="risk",newdata=rt) #利用train得到模型预测train样品风险
coxGene=rownames(multiCoxSum$coefficients)
coxGene=gsub("`","",coxGene)
outCol=c("futime","fustat",coxGene)
medianTrainRisk=median(riskScore)
risk=as.vector(ifelse(riskScore>medianTrainRisk,"high","low"))
write.table(cbind(id=rownames(cbind(rt[,outCol],riskScore,risk)),cbind(rt[,outCol],riskScore,risk)),
file="riskTrain.txt",
sep="\t",
quote=F,
row.names=F)
####################下面是验证集的风险文件输出#################
#输出test组风险文件
rtTest=read.table("test.txt",header=T,sep="\t",check.names=F,row.names=1) #读取test输入文件
rtTest[,"futime"]=rtTest[,"futime"]/365
riskScoreTest=predict(multiCox,type="risk",newdata=rtTest) #利用train得到模型预测test样品风险
riskTest=as.vector(ifelse(riskScoreTest>medianTrainRisk,"high","low"))
write.table(cbind(id=rownames(cbind(rtTest[,outCol],riskScoreTest,riskTest)),cbind(rtTest[,outCol],riskScore=riskScoreTest,risk=riskTest)),
file="riskTest.txt",
sep="\t",
quote=F,
row.names=F)
6.5 风险图形绘制
风险生存曲线
library(survival)
#绘制train组生存曲线
rt=read.table("riskTrain.txt",header=T,sep="\t")
diff=survdiff(Surv(futime, fustat) ~risk,data = rt)
pValue=1-pchisq(diff$chisq,df=1)
pValue=signif(pValue,4)
pValue=format(pValue, scientific = TRUE)
fit <- survfit(Surv(futime, fustat) ~ risk, data = rt)
summary(fit) #查看五年生存率
pdf(file="survivalTrain.pdf",width=5.5,height=5)
plot(fit,
lwd=2,
col=c("red","blue"),
xlab="Time (year)",
ylab="Survival rate",
main=paste("Survival curve (p=", pValue ,")",sep=""),
mark.time=T)
legend("topright",
c("high risk", "low risk"),
lwd=2,
col=c("red","blue"))
dev.off()
#绘制test组生存曲线
rt=read.table("riskTest.txt",header=T,sep="\t")
diff=survdiff(Surv(futime, fustat) ~risk,data = rt)
pValue=1-pchisq(diff$chisq,df=1)
pValue=signif(pValue,4)
pValue=format(pValue, scientific = TRUE)
fit <- survfit(Surv(futime, fustat) ~ risk, data = rt)
summary(fit) #查看五年生存率
pdf(file="survivalTest.pdf",width=5.5,height=5)
plot(fit,
lwd=2,
col=c("red","blue"),
xlab="Time (year)",
ylab="Survival rate",
main=paste("Survival curve (p=", pValue ,")",sep=""),
mark.time=T)
legend("topright",
c("high risk", "low risk"),
lwd=2,
col=c("red","blue"))
dev.off()
KM曲线
风险ROC曲线
library(survivalROC)
rt=read.table("riskTrain.txt",header=T,sep="\t",check.names=F,row.names=1)
pdf(file="rocTrain.pdf",width=6,height=6)
par(oma=c(0.5,1,0,1),font.lab=1.5,font.axis=1.5)
roc=survivalROC(Stime=rt$futime, status=rt$fustat, marker = rt$riskScore,
predict.time =1, method="KM")
plot(roc$FP, roc$TP, type="l", xlim=c(0,1), ylim=c(0,1),col='red',
xlab="False positive rate", ylab="True positive rate",
main=paste("ROC curve (", "AUC = ",round(roc$AUC,3),")"),
lwd = 2, cex.main=1.3, cex.lab=1.2, cex.axis=1.2, font=1.2)
abline(0,1)
dev.off()
rt=read.table("riskTest.txt",header=T,sep="\t",check.names=F,row.names=1)
pdf(file="rocTest.pdf",width=6,height=6)
par(oma=c(0.5,1,0,1),font.lab=1.5,font.axis=1.5)
roc=survivalROC(Stime=rt$futime, status=rt$fustat, marker = rt$riskScore,
predict.time =1, method="KM")
plot(roc$FP, roc$TP, type="l", xlim=c(0,1), ylim=c(0,1),col='red',
xlab="False positive rate", ylab="True positive rate",
main=paste("ROC curve (", "AUC = ",round(roc$AUC,3),")"),
lwd = 2, cex.main=1.3, cex.lab=1.2, cex.axis=1.2, font=1.2)
abline(0,1)
dev.off()
ROC曲线
