#單因素cox回歸用于篩選與預后相關的免疫基因

#文件要求熄云，帶有臨床生存與基因表達值的文件clinical+exp.txt

#生信課視頻上的代碼：不知道為什么就是運行不出來席函，于是我放棄了

library(survival)

pFilter=0.05

rt=read.table("clinical+exp.txt",header = T,sep="\t",check.names = F,row.names = 1)

outTab=data.frame()

sigGenes=c("funtime","funstat")

for (i in colnames(rt[,3:ncol(rt)])){

? ? cox<-coxph(Surv(funtime,funstat)~rt[,i],data = rt)

? ? coxSummary=summary(cox)

? ? coxP=coxSummary$coefficients[,"Pr(>|z|)"]

? ? ? ? ? if(coxP<pFilter){

? ? ? ? ? ? ? ? ? ? ? ? ? sigGenes=c(sigGenes,i)

? ? ? ? ? ? ? ? ? ? ? ? ? outTab=rbind(outTab,cbind(id=i,

? ? ? ? ? ? ? ? ? ? ? ? ? HR=coxSummary$conf.int[,"exp(coef)"],

? ? ? ? ? ? ? ? ? ? ? ? ? HRL=coxSummary$conf.int[,"lower .95"],

? ? ? ? ? ? ? ? ? ? ? ? ? HRH=coxSummary$conf.int[,"upper .95"],

? ? ? ? ? ? ? ? ? ? ? ? ? pValue=coxSummary$coefficients[,"Pr(>|z|)"]

? ? }

? ? ? }

#網友提供的方法

#需要多一個變量名的文件，是指差異表達的基因名仪芒，后來想想要不是有些搗亂的基因名炫贤，是不是用于做單因素分析的基因就是所需的covariates

rt=read.table("clinical_exp_gaibian.txt",header = T,sep="\t",check.names = F,row.names = 1)

covariates<-read.csv("names.csv")

covariates<-as.vector(covariates)

covariates

covariates=t(covariates)

time<-rt$funtime

time<-as.vector(time)

status<-rt$funstat

#有個小問題，就是變量命名里不能有-例如HLA-A有缆，會出現找不到HLA的錯誤，這里把這些帶-的都刪掉了

univ_formulas <- sapply(covariates,

? ? ? ? ? ? ? ? ? ? ? ? function(x) as.formula(paste('Surv(time, status)~', x)))

univ_models <- lapply( univ_formulas, function(x){coxph(x, data = rt)})

univ_results <- lapply(univ_models,

? ? ? ? ? ? ? ? ? ? ? ? function(x){

? ? ? ? ? ? ? ? ? ? ? ? x <- summary(x)

? ? ? ? ? ? ? ? ? ? ? ? p.value<-signif(x$wald["pvalue"])

? ? ? ? ? ? ? ? ? ? ? ? wald.test<-signif(x$wald["test"])

? ? ? ? ? ? ? ? ? ? ? ? beta<-signif(x$coef[1]);#coeficient beta

? ? ? ? ? ? ? ? ? ? ? ? HR <-signif(x$coef[2]);#exp(beta)#温亲，digits=2代表保留兩位有效數字

? ? ? ? ? ? ? ? ? ? ? ? #HR <-signif(x$coef[2],digits = 2);

? ? ? ? ? ? ? ? ? ? ? ? HR.confint.lower <- signif(x$conf.int[,"lower .95"])#,3,是3個有效數字

? ? ? ? ? ? ? ? ? ? ? ? HR.confint.upper <- signif(x$conf.int[,"upper .95"])

? ? ? ? ? ? ? ? ? ? ? ? HR <- paste0(HR, " (",

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? HR.confint.lower, "-", HR.confint.upper, ")")

? ? ? ? ? ? ? ? ? ? ? ? res<-c(beta, HR, wald.test, p.value)

? ? ? ? ? ? ? ? ? ? ? ? names(res)<-c("beta", "HR (95% CI for HR)", "wald.test",

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? "p.value")

? ? ? ? ? ? ? ? ? ? ? ? return(res)

? ? ? ? ? ? ? ? ? ? ? ? #return(exp(cbind(coef(x),confint(x))))

? ? ? ? ? ? ? ? ? ? ? })

res <- t(as.data.frame(univ_results, check.names = FALSE))

as.data.frame(res)

write.csv(res,"cox_singel.csv")

#LSAAO回歸

#文件要求：與單因素回歸的文件格式相似

LASSO是為了避免過度擬合和刪除高度相關的基因

library(glmnet)

library(survival)

setwd("C:/Users/仙女/Documents/liver cancer/LASSO")

rt=read.table("Lasso.txt",header = T,sep="\t",check.names = F,row.names = 1)

rt$funtime[rt$funtime<=0]=0.003

#LASSO回歸中生存時間不能是0棚壁，否則會報錯

x=as.matrix(rt[,c(3:ncol(rt))])

y=data.matrix(Surv(rt$funtime,rt$funstat))

fit<-glmnet(x,y,family = "cox",maxit = 1000)

pdf(file = "lambda.pdf",onefile = FALSE)

plot(fit,xvar="lambda",label=TRUE)

dev.off

cvfit<-cv.glmnet(x,y,family="cox",maxit=1000)

#做交叉驗證并輸出圖形

pdf("cvfit.pdf",onefile = FALSE)

#onefile=FALSE很重要，因為如果沒有這個輸出圖應該在pdf的第二頁栈虚，前面有個空白的pdf頁

plot(cvfit)

abline(v=log(c(cvfit$lambda.min,cvfit$lambda.1se)),lty="dashed")

dev.off()

#maxit=1000是隨機模擬1000次袖外，這代表這次做完之后把腳本關閉重新再做得出的基因是不相同的

coef<-coef(fit,s=cvfit$lambda.min)

index<-which(coef !=0)

actCoef<-coef[index]

lassoGene=rownames(coef)[index]

lassoGene=c("funtime","funstat",lassoGene)

lassoSigExp=rt[,lassoGene]

lassoSigExp=cbind(id=row.names(lassoSigExp),lassoSigExp)

write.table(lassoSigExp,file = "lassoSigExp2.txt",sep="\t",row.names = F,quote=F)

#雖然我的lambda圖沒有做成

#多變量cox回歸及建模人群預測模型的構建

#文件要求：帶有生存數據的基因表達

setwd("C:/Users/仙女/Documents/liver cancer/multi_cox")

rt=read.table("lassoSigExp.txt",header = T,sep="\t",check.names = F,row.names = 1)

rt[,"funtime"]=rt[,"funtime"]/365

#這里生存時間要轉換成年

#使用建模數據構建模型

multiCox=coxph(Surv(funtime,funstat)~.,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)

#其實每運行一次產生的結果都不一樣，即使是相同的lassoSigExp文件

#結果顯示有的基因HR的p值并不小于0.05魂务，這也是可以的曼验，它可能跟別的基因互補一起預測生存

#輸出病人的風險值

riskScore=predict(multiCox,type = "risk",newdata = rt)

coxGene=row.names(multiCoxSum$coefficients)

coxGene=gsub("`","",coxGene)

outCol=c("funtime","funstat",coxGene)

risk=as.vector(ifelse(riskScore>median(riskScore),"high","low"))

write.table(cbind(id=row.names(cbind(rt[,outCol],riskScore,risk)),cbind(rt[,outCol],riskScore,risk)),

? ? ? ? ? ? file = "risk.txt",

? ? ? ? ? ? sep="\t",

? ? ? ? ? ? quote=F,

? ? ? ? ? ? row.names = F)

#可是沒有預測的具體公式，好像是summary(multiCox)的coef值

#模型驗證

#輸入驗證人群的風險文件

rtTest=read.table("test.txt",header = T,sep="\t",check.names = F,row.names = 1)

rtTest[,"funtime"]=rtTest[,"funtime"]/365

riskScoreTest=predict(multiCox,type = "risk",newdata = rtTest)

riskTest=as.vector(ifelse(riskScoreTest>medianTrainRisk,"high","low"))

write.table(cbind(id=rownames(cbind(rtTest[,outCol],riskScoreTest,riskTest)),cbind(rtTest[,outCol],riskScoreTest,riskTest)),

? ? ? ? ? ? ? ? ? file="riskTest.txt",

? ? ? ? ? ? ? ? ? sep="\t",

? ? ? ? ? ? ? ? ? quote=F,

? ? ? ? ? ? ? ? ? row.names=F)