预后因子鉴定

点击查看代码

#!~/soft/miniconda3/envs/py3/bin/Rscript

#加载R包
suppressPackageStartupMessages(library(data.table))
suppressPackageStartupMessages(library(rtracklayer))
suppressPackageStartupMessages(library(rtracklayer))
suppressPackageStartupMessages(library(limma))
suppressPackageStartupMessages(library(edgeR))
suppressPackageStartupMessages(library(DESeq2))
suppressPackageStartupMessages(library(survival))
suppressPackageStartupMessages(library(survminer))
suppressPackageStartupMessages(library(glmnet))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(GGally))
suppressPackageStartupMessages(library(rms))
suppressPackageStartupMessages(library(survivalROC))
suppressPackageStartupMessages(library(plotROC))
suppressPackageStartupMessages(library(tidyr))
suppressPackageStartupMessages(library(optparse))


option_list <- list(
  make_option(c("-c", "--count"), help="the htseq-count from ucsc-xena which normlized with log(count + 1 ) " ),
  make_option(c("-m", "--metadata"), help="the clinical info from ucsc-xena "),
  make_option(c("-g","--gene list"),help = "gene list"),
  make_option(c("-r","--reference"),help = "gtf ")
)
opt_parser = OptionParser(option_list=option_list);
opt <- parse_args(opt_parser);



#args <- commandArgs()
# logcount <- fread(opt$c,sep = '\t',header = T)
# clinical =  fread(opt$m,sep = '\t',header = T)
# logcount <- as.data.frame(logcount)
# rownames(logcount) <- logcount[,1]
# logcount <- logcount[,-1]
# genes <- row.names(logcount)
# count <- 2^logcount-1
# save(count,file = 'count.rda')

# ## 做差异基因时筛选配对的肿瘤-癌旁样本
# exp = count
# 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.1,]
# group_list=ifelse(as.numeric(substr(colnames(data),14,15)) < 10,'tumor','normal')
# tcga_normal <- data[,group_list == 'normal']
# tcga_tumor <- data[, group_list == 'tumor']
# tcga_tumor_need <- tcga_tumor[,((substr(colnames(tcga_tumor),1,12))%in%(substr(colnames(tcga_normal),1,12)))]
# same <- intersect(row.names(tcga_normal),row.names(tcga_tumor_need))
# mydata <- cbind(tcga_tumor_need[same,],tcga_normal[same,])
# integerData = apply(mydata,2,as.integer);rownames(integerData) = rownames(mydata)
# integerData = as.data.frame(integerData)
# integerData1 = dplyr::add_rownames(integerData,'geneid')
# expr_df_nopoint <- integerData1 %>% 
#   tidyr::separate(geneid,into = c("gene_id"),sep="\\.") 
  
# load("gtf_df.Rda") #Homo_sapiens.GRCh38.107.chr.gtf
# g2s <- gtf_df %>% 
#   dplyr::filter(type=="gene")%>% 
#    dplyr::select(c(gene_name,gene_id))
# mRNA_exprSet <- gtf_df %>% 
#   dplyr::filter(type=="gene") %>% #筛选gene,和编码指标
#   dplyr::select(c(gene_name,gene_id)) %>% 
#   dplyr::inner_join(expr_df_nopoint,by ="gene_id")
# mRNA_exprSet = mRNA_exprSet[,-2]
# ## 重复的基因名取均值
# expr_mean=aggregate(.~gene_name,mean,data=mRNA_exprSet)
# expr =  tibble::column_to_rownames(expr_mean,'gene_name')
# group_list <- c(rep('tumor',ncol(tcga_tumor_need)),rep('normal',ncol(tcga_normal)))
# expr2 = apply(expr,2,as.integer);rownames(expr2) = rownames(expr)
# save(expr2, group_list, file = 'exprSet_by_group_list.Rdata')
# data = expr2
# group_list = factor(group_list)
# design <- model.matrix(~0+group_list)
# rownames(design) = colnames(data)
# colnames(design) <- levels(group_list)
# DGElist <- DGEList( counts = data, group = group_list)
# keep.exprs<- edgeR::filterByExpr(data,group=group_list)
# DGElist <- DGElist[ keep.exprs, , keep.lib.sizes = FALSE ]

# DGElist <- calcNormFactors( DGElist )
# DGElist <- estimateGLMCommonDisp(DGElist, design)
# DGElist <- estimateGLMTrendedDisp(DGElist, design)
# DGElist <- estimateGLMTagwiseDisp(DGElist, design)

# fit <- glmFit(DGElist, design)
# results <- glmLRT(fit, contrast = c(-1, 1)) 
# nrDEG_edgeR <- topTags(results, n = nrow(DGElist))
# nrDEG_edgeR <- as.data.frame(nrDEG_edgeR)
# head(nrDEG_edgeR)
# padj = 0.05 # 自定义
# #print(padj)
# foldChange= 1 # 自定义
# #print(foldChange)
# nrDEG_edgeR_signif  = nrDEG_edgeR[(nrDEG_edgeR$FDR < padj & 
#                        (nrDEG_edgeR$logFC > foldChange | nrDEG_edgeR$logFC < (-foldChange))),]
# nrDEG_edgeR_signif = nrDEG_edgeR_signif[order(nrDEG_edgeR_signif$logFC),]
# save(nrDEG_edgeR_signif,file = 'nrDEG_edgeR_signif.Rdata')


load(opt$r) #Homo_sapiens.GRCh38.107.chr.gtf
logcount <- fread(opt$c,sep = '\t',header = T)
clinical =  fread(opt$m,sep = '\t',header = T)
gene.list = fread(opt$g,header =T) #格式为一列，列名为 'gene'


logcount <- as.data.frame(logcount)
rownames(logcount) <- logcount[,1]
logcount <- logcount[,-1]
genes <- row.names(logcount)
count <- 2^logcount-1
save(count,file = 'count.rda')

exp = count
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)>10,]
group_list=ifelse(as.numeric(substr(colnames(data),14,15)) < 10,'tumor','normal')
tcga_normal <- data[,group_list == 'normal']
tcga_tumor <- data[, group_list == 'tumor']
tcga_tumor = as.data.frame(tcga_tumor)
tcga_tumor = dplyr::add_rownames(tcga_tumor,'geneid')

expr_df_nopoint <- tcga_tumor %>% 
  tidyr::separate(geneid,into = c("gene_id"),sep="\\.") 
g2s <- gtf_df %>% 
  dplyr::filter(type=="gene")%>% 
   dplyr::select(c(gene_name,gene_id))
mRNA_exprSet <- gtf_df %>% 
  dplyr::filter(type=="gene") %>% #筛选gene,和编码指标
  dplyr::select(c(gene_name,gene_id)) %>% 
  dplyr::inner_join(expr_df_nopoint,by ="gene_id")
mRNA_exprSet = mRNA_exprSet[,-2]
dim(mRNA_exprSet)
expr_mean=aggregate(.~gene_name,mean,data=mRNA_exprSet)
expr =  tibble::column_to_rownames(expr_mean,'gene_name')

clinical = clinical[,c(1,2,4)] #不同来源的clinical信息可能不同 ，这是基因ucsc xena 三列为sampleid OS OS.time
meta = clinical[clinical$sample %in% colnames(expr),]
tumor=expr[rownames(expr) %in% gene.list$gene,]
tumor=tumor[,colnames(tumor) %in% meta$sample]
mySurv=with(meta,Surv(OS.time, OS))
log_rank_p <- apply(tumor , 1 , function(gene){
  # gene=exprSet[1,]
  meta$group=ifelse(gene>median(gene),'high','low')  
  data.survdiff=survdiff(mySurv~group,data=meta)
  p.val = 1 - pchisq(data.survdiff$chisq, length(data.survdiff$n) - 1)
  return(p.val)
})
save(tumor,meta,file = 'coxinput.rda')
load('coxinput.rda')
tumor_cpm = log2(tumor + 1) #对原始count取对数
meta$OS.time=meta$OS.time/365
t_tumor = t(tumor_cpm)
t_tumor = as.data.frame(t_tumor)
t_tumor = dplyr::add_rownames(t_tumor,'sample')
cox.input = merge(meta,t_tumor,by = 'sample')
cox.input = tibble::column_to_rownames(cox.input,var = 'sample')
cox.input = as.data.frame(cox.input)

pFilter=0.05     
#单因素过滤条件
outTab=data.frame()
sigGenes=c("OS.time","OS")
for(gene in colnames(cox.input[,3:ncol(cox.input)])){
    cox=coxph(Surv(OS.time, OS) ~ cox.input[,gene], data = cox.input)
    coxSummary = summary(cox)
    coxP=coxSummary$coefficients[,"Pr(>|z|)"]
        
    if(coxP<pFilter){
        sigGenes=c(sigGenes,gene)
        outTab=rbind(outTab,
                   cbind(gene=gene,
                   HR=coxSummary$conf.int[,"exp(coef)"],
                   HR.95L=coxSummary$conf.int[,"lower .95"],
                   HR.95H=coxSummary$conf.int[,"upper .95"],
                   pvalue=coxP) )
    }
}
write.table(outTab,file="uniCox.txt",sep="\t",row.names=F,quote=F)
surSigExp=cox.input[,sigGenes]
write.table(surSigExp,file="uniSigExp.txt",sep="\t",row.names=F,quote=F)

set.seed(1234)
rt = surSigExp
x=as.matrix(rt[,c(3:ncol(rt))])
y=data.matrix(Surv(rt$OS.time,rt$OS))
fit <- glmnet(x, y, family = "cox", maxit = 1000)
pdf("lasso.lambda.pdf")
plot(fit, xvar = "lambda", label = TRUE)
dev.off()
cvfit <- cv.glmnet(x, y, family="cox", maxit = 1000)
pdf("lasso.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("OS.time","OS",lassoGene)
lassoSigExp=rt[,lassoGene]
lassoSigExp=cbind(id=row.names(lassoSigExp),lassoSigExp)
write.table(lassoSigExp,file="lasso.SigExp.txt",sep="\t",row.names=F,quote=F)
#COX模型构建
rt=read.table("lasso.SigExp.txt",header=T,sep="\t",check.names=F,row.names=1)
multiCox=coxph(Surv(OS.time, OS) ~ ., 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)
outTab=gsub("`","",outTab)
write.table(outTab,file="multi.Cox.txt",sep="\t",row.names=F,quote=F)
#输出病人风险值
riskScore=predict(multiCox, type="risk", newdata=rt)
coxGene=rownames(multiCoxSum$coefficients)
coxGene=gsub("`", "", coxGene)
outCol=c("OS.time", "OS", coxGene)
riskOut=cbind(rt[,outCol], riskScore)
riskOut=cbind(id=rownames(riskOut), riskOut)
write.table(riskOut, file="riskScore.txt", sep="\t", quote=F, row.names=F)
bioForest=function(coxFile=null, forestFile=null, forestCol=null){
    #读取输入文件
    rt <- read.table(coxFile,header=T,sep="\t",row.names=1,check.names=F)
    gene <- rownames(rt)
    hr <- sprintf("%.3f",rt$"HR")
    hrLow  <- sprintf("%.3f",rt$"HR.95L")
    hrHigh <- sprintf("%.3f",rt$"HR.95H")
    Hazard.ratio <- paste0(hr,"(",hrLow,"-",hrHigh,")")
    pVal <- ifelse(rt$pvalue<0.001, "<0.001", sprintf("%.3f", rt$pvalue))
        
    #输出图形
    pdf(file=forestFile, width=9, height=7)
    n <- nrow(rt)
    nRow <- n+1
    ylim <- c(1,nRow)
    layout(matrix(c(1,2),nc=2),width=c(3,2.5))
        
    #绘制森林图左边的临床信息
    xlim = c(0,3)
    par(mar=c(4,2.5,2,1))
    plot(1,xlim=xlim,ylim=ylim,type="n",axes=F,xlab="",ylab="")
    text.cex=0.8
    text(0,n:1,gene,adj=0,cex=text.cex)
    text(2.08-0.5*0.2,n:1,pVal,adj=1,cex=text.cex);text(2.08-0.5*0.2,n+1,'pvalue',cex=text.cex,font=2,adj=1)
    text(3.12,n:1,Hazard.ratio,adj=1,cex=text.cex);text(3.12,n+1,'Hazard ratio',cex=text.cex,font=2,adj=1,)
        
    #绘制森林图
    par(mar=c(4,1,2,1),mgp=c(2,0.5,0))
    xlim = c(0,max(as.numeric(hrLow),as.numeric(hrHigh)))
    plot(1,xlim=xlim,ylim=ylim,type="n",axes=F,ylab="",xaxs="i",xlab="Hazard ratio")
    arrows(as.numeric(hrLow),n:1,as.numeric(hrHigh),n:1,angle=90,code=3,length=0.05,col="darkblue",lwd=2.5)
    abline(v=1,col="black",lty=2,lwd=2)
    boxcolor = ifelse(as.numeric(hr) > 1, forestCol[1], forestCol[2])
    points(as.numeric(hr), n:1, pch = 15, col = boxcolor, cex=1.6)
    axis(1)
    dev.off()
}
#绘制模型森林图
bioForest(coxFile="multi.Cox.txt", forestFile="model.multiForest.pdf", forestCol=c("red","green"))

#绘制模型中基因单因素森林图
uniRT=read.table("uniCox.txt",header=T,sep="\t",row.names=1,check.names=F)
uniRT=uniRT[coxGene,]
uniRT=cbind(id=row.names(uniRT), uniRT)
write.table(uniRT, file="unicox.forest.txt", sep="\t", row.names=F, quote=F)
bioForest(coxFile="unicox.forest.txt", forestFile="model.uniForest.pdf", forestCol=c("red","green"))