前言：之前做过一些数据分析与数据挖掘相关的工作，最近抽空将之前做的内容简单整理一下，方便查看，主要使用R语言和PERL脚本语言，使用TCGA和ICGC数据库中的临床数据，做类似的分析可以参考一下，如果想查看详细内容与数据可以通过本人的Gitee及Github仓库下载，链接于篇尾附上。
一、标题：Effect of HSP90AB1 on the local immune response of hepatocellular carcinoma and it realtionship to prognosis（HSP90β对肝癌局部免疫的影响及对肝癌患者预后的影响）
二、部分代码及结果展示：
1、整理TCGA数据库肝细胞癌临床数据的部分PERL脚本
use strict; #use warnings; use XML::Simple; opendir(RD, ".") or die $!; my @dirs=readdir(RD); closedir(RD); open(WF,">clinical.xls") or die $!; print WF "IdtfutimetfustattAgetGendertGradetStagetTtMtNn"; foreach my $dir(@dirs){ #print $dir . "n"; next if($dir eq '.'); next if($dir eq '..'); #print $dir . "n"; if(-d $dir){ opendir(RD,"$dir") or die $!; while(my $xmlfile=readdir(RD)){ if($xmlfile=~/.xml$/){ #print "$dir\$xmlfilen"; my $userxs = XML::Simple->new(KeyAttr => "name"); my $userxml=""; if(-f "$dir/$xmlfile"){ $userxml = $userxs->XMLin("$dir/$xmlfile"); }else{ $userxml = $userxs->XMLin("$dir$xmlfile"); } # print output #open(WF,">dumper.txt") or die $!; #print WF Dumper($userxml); #close(WF); my $disease_code=$userxml->{'admin:admin'}{'admin:disease_code'}{'content'}; #get disease code my $disease_code_lc=lc($disease_code); my $patient_key=$disease_code_lc . ':patient'; #ucec:patient my $follow_key=$disease_code_lc . ':follow_ups'; my $patient_barcode=$userxml->{$patient_key}{'shared:bcr_patient_barcode'}{'content'}; #TCGA-AX-A1CJ my $gender=$userxml->{$patient_key}{'shared:gender'}{'content'}; #male/female my $age=$userxml->{$patient_key}{'clin_shared:age_at_initial_pathologic_diagnosis'}{'content'}; my $race=$userxml->{$patient_key}{'clin_shared:race_list'}{'clin_shared:race'}{'content'}; #white/black my $grade=$userxml->{$patient_key}{'shared:neoplasm_histologic_grade'}{'content'}; #G1/G2/G3 my $clinical_stage=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:clinical_stage'}{'content'}; #stage I my $clinical_T=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:tnm_categories'}{'shared_stage:clinical_categories'}{'shared_stage:clinical_T'}{'content'}; my $clinical_M=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:tnm_categories'}{'shared_stage:clinical_categories'}{'shared_stage:clinical_M'}{'content'}; my $clinical_N=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:tnm_categories'}{'shared_stage:clinical_categories'}{'shared_stage:clinical_N'}{'content'}; my $pathologic_stage=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:pathologic_stage'}{'content'}; #stage I my $pathologic_T=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:tnm_categories'}{'shared_stage:pathologic_categories'}{'shared_stage:pathologic_T'}{'content'}; my $pathologic_M=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:tnm_categories'}{'shared_stage:pathologic_categories'}{'shared_stage:pathologic_M'}{'content'}; my $pathologic_N=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:tnm_categories'}{'shared_stage:pathologic_categories'}{'shared_stage:pathologic_N'}{'content'}; $gender=(defined $gender)?$gender:"unknow"; $age=(defined $age)?$age:"unknow"; $race=(defined $race)?$race:"unknow"; $grade=(defined $grade)?$grade:"unknow"; $clinical_stage=(defined $clinical_stage)?$clinical_stage:"unknow"; $clinical_T=(defined $clinical_T)?$clinical_T:"unknow"; $clinical_M=(defined $clinical_M)?$clinical_M:"unknow"; $clinical_N=(defined $clinical_N)?$clinical_N:"unknow"; $pathologic_stage=(defined $pathologic_stage)?$pathologic_stage:"unknow"; $pathologic_T=(defined $pathologic_T)?$pathologic_T:"unknow"; $pathologic_M=(defined $pathologic_M)?$pathologic_M:"unknow"; $pathologic_N=(defined $pathologic_N)?$pathologic_N:"unknow"; my $survivalTime=""; my $vital_status=$userxml->{$patient_key}{'clin_shared:vital_status'}{'content'}; my $followup=$userxml->{$patient_key}{'clin_shared:days_to_last_followup'}{'content'}; my $death=$userxml->{$patient_key}{'clin_shared:days_to_death'}{'content'}; if($vital_status eq 'Alive'){ $survivalTime="$followupt0"; } else{ $survivalTime="$deatht1"; } for my $i(keys %{$userxml->{$patient_key}{$follow_key}}){ eval{ $followup=$userxml->{$patient_key}{$follow_key}{$i}{'clin_shared:days_to_last_followup'}{'content'}; $vital_status=$userxml->{$patient_key}{$follow_key}{$i}{'clin_shared:vital_status'}{'content'}; $death=$userxml->{$patient_key}{$follow_key}{$i}{'clin_shared:days_to_death'}{'content'}; }; if($@){ for my $j(0..5){ #假设最多有6次随访 my $followup_for=$userxml->{$patient_key}{$follow_key}{$i}[$j]{'clin_shared:days_to_last_followup'}{'content'}; my $vital_status_for=$userxml->{$patient_key}{$follow_key}{$i}[$j]{'clin_shared:vital_status'}{'content'}; my $death_for=$userxml->{$patient_key}{$follow_key}{$i}[$j]{'clin_shared:days_to_death'}{'content'}; if( ($followup_for =~ /d+/) || ($death_for =~ /d+/) ){ $followup=$followup_for; $vital_status=$vital_status_for; $death=$death_for; my @survivalArr=split(/t/,$survivalTime); if($vital_status eq 'Alive'){ if($followup>$survivalArr[0]){ $survivalTime="$followupt0"; } } else{ if($death>$survivalArr[0]){ $survivalTime="$deatht1"; } } } } } my @survivalArr=split(/t/,$survivalTime); if($vital_status eq 'Alive'){ if($followup>$survivalArr[0]){ $survivalTime="$followupt0"; } } else{ if($death>$survivalArr[0]){ $survivalTime="$deatht1"; } } } print WF "$patient_barcodet$survivalTimet$aget$gendert$gradet$pathologic_staget$pathologic_Tt$pathologic_Mt$pathologic_Nn"; } } close(RD); } } close(WF); 2、使用R语言分析正常组与肿瘤组中HSP90AB1的表达情况
#if (!requireNamespace("BiocManager", quietly = TRUE)) # install.packages("BiocManager") #BiocManager::install("limma") #install.packages("ggplot2") #install.packages("ggpubr") #引用包 library(limma) library(ggplot2) library(ggpubr) expFile="symbol.txt" #表达输入文件 gene="VCAN" #基因的名称 setwd("C:\Users\lexb4\Desktop\geneImmune\07.diff") #设置工作目录 #读取基因表达文件,并对数据进行处理 rt=read.table(expFile, header=T, sep="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=t(data[gene,,drop=F]) #正常和肿瘤数目 group=sapply(strsplit(rownames(data),"\-"), "[", 4) group=sapply(strsplit(group,""), "[", 1) group=gsub("2", "1", group) conNum=length(group[group==1]) #正常组样品数目 treatNum=length(group[group==0]) #肿瘤组样品数目 Type=c(rep(1,conNum), rep(2,treatNum)) #差异分析 exp=cbind(data, Type) exp=as.data.frame(exp) colnames(exp)=c("gene", "Type") exp$Type=ifelse(exp$Type==1, "Normal", "Tumor") exp$gene=log2(exp$gene+1) #设置比较组 group=levels(factor(exp$Type)) exp$Type=factor(exp$Type, levels=group) comp=combn(group,2) my_comparisons=list() for(i in 1:ncol(comp)){my_comparisons[[i]]<-comp[,i]} #绘制boxplot boxplot=ggboxplot(exp, x="Type", y="gene", color="Type", xlab="", ylab=paste0(gene, " expression"), legend.title="Type", palette = c("blue","red"), add = "jitter")+ stat_compare_means(comparisons=my_comparisons,symnum.args=list(cutpoints = c(0, 0.001, 0.01, 0.05, 1), symbols = c("***", "**", "*", "ns")),label = "p.signif") #输出图片 pdf(file=paste0(gene,".diff.pdf"), width=5, height=4.5) print(boxplot) dev.off() 
3、使用R语言分析不同类型免疫细胞在肝细胞癌中的表达水平及相关关系
#install.packages("corrplot") library(corrplot) #引用包 immFile="CIBERSORT-Results.txt" #免疫细胞浸润的结果文件 pFilter=0.05 #免疫细胞浸润结果过滤条件 setwd("C:\Users\Administrator\Desktop\geneimmune\10immunePlot") #设置工作目录 #读取免疫细胞浸润的结果文件，并对数据进行整理 immune=read.table(immFile, header=T, sep="t", check.names=F, row.names=1) immune=immune[immune[,"P-value"]<pFilter,] immune=as.matrix(immune[,1:(ncol(immune)-3)]) data=t(immune) #绘制柱状图 col=rainbow(nrow(data), s=0.7, v=0.7) pdf(file="barplot.pdf", width=22, height=10) par(las=1,mar=c(8,5,4,16),mgp=c(3,0.1,0),cex.axis=1.5) a1=barplot(data,col=col,yaxt="n",ylab="Relative Percent",xaxt="n",cex.lab=1.8) a2=axis(2,tick=F,labels=F) axis(2,a2,paste0(a2*100,"%")) axis(1,a1,labels=F) par(srt=60,xpd=T);text(a1,-0.02,colnames(data),adj=1,cex=0.6);par(srt=0) ytick2=cumsum(data[,ncol(data)]);ytick1=c(0,ytick2[-length(ytick2)]) legend(par('usr')[2]*0.98,par('usr')[4],legend=rownames(data),col=col,pch=15,bty="n",cex=1.3) dev.off() #删除正常样品 group=sapply(strsplit(colnames(data),"\-"), "[", 4) group=sapply(strsplit(group,""), "[", 1) group=gsub("2", "1", group) data=data[,group==0,drop=F] #绘制免疫细胞相关性的图形 pdf(file="corrplot.pdf", width=13, height=13) par(oma=c(0.5,1,1,1.2)) immune=immune[,colMeans(immune)>0] M=cor(immune) corrplot(M, method = "color", order = "hclust", tl.col="black", addCoef.col = "black", number.cex = 0.8, col=colorRampPalette(c("blue", "white", "red"))(50) ) dev.off()
4、使用R语言分析正常组及肝癌组中不同免疫细胞浸润水平
#install.packages("pheatmap") #install.packages("vioplot") #引用包 library(vioplot) library(pheatmap) input="CIBERSORT-Results.txt" #免疫细胞浸润文件 pFilter=0.05 #免疫细胞浸润结果过滤条件 setwd("C:\Users\Administrator\Desktop\生信文章\geneimmune\11heatmap\vioplot-high") #设置工作目录 #读取免疫结果文件，并对数据进行整理 immune=read.table("CIBERSORT-Results.txt", header=T, sep="t", check.names=F, row.names=1) immune=immune[immune[,"P-value"]<pFilter,,drop=F] immune=as.matrix(immune[,1:(ncol(immune)-3)]) data=t(immune) #正常和肿瘤数目 group=sapply(strsplit(colnames(data),"\-"), "[", 4) group=sapply(strsplit(group,""), "[", 1) group=gsub("2", "1", group) conNum=length(group[group==1]) #正常组样品数目 treatNum=length(group[group==0]) #肿瘤组样品数目 #定义热图的注释文件 Type=c(rep("Normal",conNum), rep("Tumor",treatNum)) names(Type)=colnames(data) Type=as.data.frame(Type) #绘制热图 pdf(file="heatmap.pdf", width=12, height=6) pheatmap(data, annotation=Type, color = colorRampPalette(c(rep("green",1), rep("black",1), rep("red",3)))(100), cluster_cols =F, show_colnames=F, fontsize = 8, fontsize_row=7, fontsize_col=5) dev.off() #绘制小提琴图 data=t(data) outTab=data.frame() pdf(file="vioplot.pdf", width=13, height=8) par(las=1, mar=c(10,6,3,3)) x=c(1:ncol(data)) y=c(1:ncol(data)) xMax=ncol(data)*3-2 plot(x,y, xlim=c(0,xMax),ylim=c(min(data),max(data)+0.02), main="", xlab="", ylab="Fraction", pch=21, col="white", xaxt="n") #对每个免疫细胞循环，绘制小提琴图，正常样品用绿色表示，肿瘤样品用红色表示 for(i in 1:ncol(data)){ if(sd(data[1:conNum,i])==0){ data[1,i]=0.00001 } if(sd(data[(conNum+1):(conNum+treatNum),i])==0){ data[(conNum+1),i]=0.00001 } conData=data[1:conNum,i] treatData=data[(conNum+1):(conNum+treatNum),i] vioplot(conData,at=3*(i-1),lty=1,add = T,col = 'green') vioplot(treatData,at=3*(i-1)+1,lty=1,add = T,col = 'red') wilcoxTest=wilcox.test(conData, treatData) p=wilcoxTest$p.value if(p<pFilter){ cellPvalue=cbind(Cell=colnames(data)[i], pvalue=p) outTab=rbind(outTab, cellPvalue) } mx=max(c(conData,treatData)) lines(c(x=3*(i-1)+0.2,x=3*(i-1)+0.8),c(mx,mx)) text(x=3*(i-1)+0.5, y=mx+0.02, labels=ifelse(p<0.001, paste0("p<0.001"), paste0("p=",sprintf("%.03f",p))), cex = 0.8) } legend("topright", c("Normal", "Tumor"), lwd=5,bty="n",cex=1.2, col=c("green","red")) text(seq(1,xMax,3),-0.05,xpd = NA,labels=colnames(data),cex = 1,srt = 45,pos=2) dev.off() #输出免疫细胞和p值表格文件 write.table(outTab,file="diff.result.txt",sep="t",row.names=F,quote=F)
5、不同拷贝子数目的HSP90β对中性粒细胞和CD8阳性T细胞在肝癌局部浸润的影响
6、HSP90β基因的表达水平、拷贝子水平及甲基化水平与不同淋巴细胞数量之间的关系
7、HSP90β基因与免疫调节基因之间的关系
8、与HSP90β相关的免疫调节基因的蛋白互作网络
9、GO及相关通路分析
10、整理患者基因表达水平与临床生存信息
#if (!requireNamespace("BiocManager", quietly = TRUE)) # install.packages("BiocManager") #BiocManager::install("limma") library(limma) #引用包 expFile="symbol.txt" #表达数据文件 cliFile="time.txt" #临床数据文件 geneFile="gene.txt" #基因列表文件 setwd("C:\Users\Administrator\Desktop\geneimmune\24mergeTime") #工作目录（需修改） #读取表达文件，并对输入文件整理 rt=read.table(expFile, header=T, sep="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,] #读取免疫基因的表达量 gene=read.table(geneFile, header=F, sep="t", check.names=F) sameGene=intersect(as.vector(gene[,1]), row.names(data)) data=data[sameGene,] #删掉正常样品 group=sapply(strsplit(colnames(data),"\-"), "[", 4) group=sapply(strsplit(group,""), "[", 1) group=gsub("2", "1", group) data=data[,group==0] colnames(data)=gsub("(.*?)\-(.*?)\-(.*?)\-(.*?)\-.*", "\1\-\2\-\3", colnames(data)) data=t(data) data=avereps(data) #读取生存数据 cli=read.table(cliFile, header=T, sep="t", check.names=F, row.names=1) #读取临床文件 #数据合并并输出结果 sameSample=intersect(row.names(data), row.names(cli)) data=data[sameSample,] cli=cli[sameSample,] out=cbind(cli,data) out=cbind(id=row.names(out),out) write.table(out,file="expTime.txt",sep="t",row.names=F,quote=F) 11、筛选HSP90β蛋白互作网络中预后相关的免疫调节基因，绘制森林图
#install.packages('survival') library(survival) #引用包 coxPfilter=0.05 #显著性过滤标准 inputFile="expTime.txt" #输入文件 setwd("C:\Users\Administrator\Desktop\geneimmune\25uniCox") #设置工作目录 rt=read.table(inputFile, header=T, sep="t", check.names=F, row.names=1) #读取输入文件 rt$futime=rt$futime/365 rt[,3:ncol(rt)]=log2(rt[,3:ncol(rt)]+1) #对基因进行循环，找出预后相关的基因 outTab=data.frame() sigGenes=c("futime","fustat") for(i in colnames(rt[,3:ncol(rt)])){ #cox分析 cox <- coxph(Surv(futime, fustat) ~ rt[,i], data = rt) coxSummary = summary(cox) coxP=coxSummary$coefficients[,"Pr(>|z|)"] #保留预后相关的基因 if(coxP<coxPfilter){ sigGenes=c(sigGenes,i) 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|)"]) ) } } #输出单因素的结果 write.table(outTab,file="uniCox.txt",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) ############绘制森林图函数############ bioForest=function(coxFile=null, forestFile=null, forestCol=null){ #读取输入文件 rt <- read.table(coxFile, header=T, sep="t", check.names=F, row.names=1) 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=6.5, height=5) 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(1.5-0.5*0.2,n:1,pVal,adj=1,cex=text.cex);text(1.5-0.5*0.2,n+1,'pvalue',cex=text.cex,font=2,adj=1) text(3,n:1,Hazard.ratio,adj=1,cex=text.cex);text(3,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="uniCox.txt", forestFile="forest.pdf", forestCol=c("red","green")) 12、使用筛选出的基因构建预后模型
#install.packages("glmnet") #install.packages("survival") #install.packages('survminer') #引用包 library(glmnet) library(survival) library(survminer) inputFile="uniSigExp.txt" #单因素显著基因的表达输入文件 setwd("C:\Users\lexb4\Desktop\geneImmune\26.model") #设置工作目录 rt=read.table(inputFile, header=T, sep="t", row.names=1, check.names=F) #读取输入文件 #COX模型构建 multiCox=coxph(Surv(futime, fustat) ~ ., data = rt) multiCox=step(multiCox, direction="both") multiCoxSum=summary(multiCox) #输出模型相关信息 outMultiTab=data.frame() outMultiTab=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|)"]) outMultiTab=cbind(id=row.names(outMultiTab),outMultiTab) write.table(outMultiTab, file="multiCox.txt", sep="t", row.names=F, quote=F) #输出风险文件 score=predict(multiCox, type="risk", newdata=rt) coxGene=rownames(multiCoxSum$coefficients) coxGene=gsub("`", "", coxGene) outCol=c("futime", "fustat", coxGene) risk=as.vector(ifelse(score>median(score), "high", "low")) outTab=cbind(rt[,outCol], riskScore=as.vector(score), risk) write.table(cbind(id=rownames(outTab),outTab), file="risk.txt", sep="t", quote=F, row.names=F) #绘制森林图 pdf(file="multi.forest.pdf", width=10, height=6, onefile=FALSE) ggforest(multiCox, data=rt, main = "Hazard ratio", cpositions = c(0.02,0.22, 0.4), fontsize = 0.7, refLabel = "reference", noDigits = 2) dev.off()
13、绘制该预后模型高低风险组的生存曲线
#install.packages("survival") #install.packages("survminer") #引用包 library(survival) library(survminer) setwd("C:\Users\lexb4\Desktop\geneImmune\27.survival") #设置工作目录 #定义生存曲线的函数 bioSurvival=function(inputFile=null, outFile=null){ #读取输入文件 rt=read.table(inputFile, header=T, sep="t", check.names=F) #比较高低风险组生存差异，得到显著性p值 diff=survdiff(Surv(futime, fustat) ~ risk, data=rt) pValue=1-pchisq(diff$chisq, df=1) if(pValue<0.001){ pValue="p<0.001" }else{ pValue=paste0("p=",sprintf("%.03f",pValue)) } fit <- survfit(Surv(futime, fustat) ~ risk, data = rt) #print(surv_median(fit)) #绘制生存曲线 surPlot=ggsurvplot(fit, data=rt, conf.int=T, pval=pValue, pval.size=6, surv.median.line = "hv", legend.title="Risk", legend.labs=c("High risk", "Low risk"), xlab="Time(years)", break.time.by = 1, palette=c("red", "blue"), risk.table=TRUE, risk.table.title="", risk.table.col = "strata", risk.table.height=.25) pdf(file=outFile, onefile=FALSE, width=6.5, height=5.5) print(surPlot) dev.off() } #调用函数，绘制生存曲线 bioSurvival(inputFile="risk.txt", outFile="survival.pdf")
14、绘制不同的风险曲线
#install.packages("pheatmap") library(pheatmap) #引用包 setwd("C:\Users\Administrator\Desktop\geneimmune\28riskPlot") #设置工作目录 #定义风险曲线的函数 bioRiskPlot=function(inputFile=null, riskScoreFile=null, survStatFile=null, heatmapFile=null){ rt=read.table(inputFile, header=T, sep="t", check.names=F, row.names=1) #读取输入文件 rt=rt[order(rt$riskScore),] #按照风险打分对样品排序 #绘制风险曲线 riskClass=rt[,"risk"] lowLength=length(riskClass[riskClass=="low"]) highLength=length(riskClass[riskClass=="high"]) lowMax=max(rt$riskScore[riskClass=="low"]) line=rt[,"riskScore"] line[line>10]=10 pdf(file=riskScoreFile, width=7, height=4) plot(line, type="p", pch=20, xlab="Patients (increasing risk socre)", ylab="Risk score", col=c(rep("green",lowLength),rep("red",highLength)) ) abline(h=lowMax,v=lowLength,lty=2) legend("topleft", c("High risk", "Low Risk"),bty="n",pch=19,col=c("red","green"),cex=1.2) dev.off() #绘制生存状态图 color=as.vector(rt$fustat) color[color==1]="red" color[color==0]="green" pdf(file=survStatFile, width=7, height=4) plot(rt$futime, pch=19, xlab="Patients (increasing risk socre)", ylab="Survival time (years)", col=color) legend("topleft", c("Dead", "Alive"),bty="n",pch=19,col=c("red","green"),cex=1.2) abline(v=lowLength,lty=2) dev.off() #绘制风险热图 rt1=rt[c(3:(ncol(rt)-2))] rt1=t(rt1) annotation=data.frame(type=rt[,ncol(rt)]) rownames(annotation)=rownames(rt) pdf(file=heatmapFile, width=7, height=4) pheatmap(rt1, annotation=annotation, cluster_cols = FALSE, cluster_rows = FALSE, show_colnames = F, scale="row", color = colorRampPalette(c(rep("green",3), "white", rep("red",3)))(50), fontsize_col=3, fontsize=7, fontsize_row=8) dev.off() } #调用函数，绘制风险曲线 bioRiskPlot(inputFile="risk.txt", riskScoreFile="riskScore.pdf", survStatFile="survStat.pdf", heatmapFile="heatmap.pdf")
15、绘制不同风险因素森林图比较，并进行预后分析
#install.packages('survival') library(survival) #引用包 setwd("C:\Users\Administrator\Desktop\geneimmune\29indep") #设置工作目录 ############绘制森林图函数############ bioForest=function(coxFile=null, forestFile=null, forestCol=null){ #读取输入文件 rt <- read.table(coxFile, header=T, sep="t", check.names=F, row.names=1) 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=6.5, height=4.5) 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(1.5-0.5*0.2,n:1,pVal,adj=1,cex=text.cex);text(1.5-0.5*0.2,n+1,'pvalue',cex=text.cex,font=2,adj=1) text(3.1,n:1,Hazard.ratio,adj=1,cex=text.cex);text(3.1,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=3) abline(v=1, col="black", lty=2, lwd=2) boxcolor = ifelse(as.numeric(hr) > 1, forestCol, forestCol) points(as.numeric(hr), n:1, pch = 15, col = boxcolor, cex=2) axis(1) dev.off() } ############绘制森林图函数############ #定义独立预后分析函数 indep=function(riskFile=null,cliFile=null,uniOutFile=null,multiOutFile=null,uniForest=null,multiForest=null){ risk=read.table(riskFile, header=T, sep="t", check.names=F, row.names=1) #读取风险文件 cli=read.table(cliFile, header=T, sep="t", check.names=F, row.names=1) #读取临床文件 #数据合并 sameSample=intersect(row.names(cli),row.names(risk)) risk=risk[sameSample,] cli=cli[sameSample,] rt=cbind(futime=risk[,1], fustat=risk[,2], cli, riskScore=risk[,(ncol(risk)-1)]) #单因素独立预后分析 uniTab=data.frame() for(i in colnames(rt[,3:ncol(rt)])){ cox <- coxph(Surv(futime, fustat) ~ rt[,i], data = rt) coxSummary = summary(cox) uniTab=rbind(uniTab, 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|)"]) ) } write.table(uniTab,file=uniOutFile,sep="t",row.names=F,quote=F) bioForest(coxFile=uniOutFile, forestFile=uniForest, forestCol="green") #多因素独立预后分析 uniTab=uniTab[as.numeric(uniTab[,"pvalue"])<1,] rt1=rt[,c("futime", "fustat", as.vector(uniTab[,"id"]))] multiCox=coxph(Surv(futime, fustat) ~ ., data = rt1) multiCoxSum=summary(multiCox) multiTab=data.frame() multiTab=cbind( 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|)"]) multiTab=cbind(id=row.names(multiTab),multiTab) write.table(multiTab,file=multiOutFile,sep="t",row.names=F,quote=F) bioForest(coxFile=multiOutFile, forestFile=multiForest, forestCol="red") } #独立预后分析 indep(riskFile="risk.txt", cliFile="clinical.txt", uniOutFile="uniCox.txt", multiOutFile="multiCox.txt", uniForest="uniForest.pdf", multiForest="multiForest.pdf")
16、绘制ROC曲线
#install.packages("survival") #install.packages("survminer") #install.packages("timeROC") #引用包 library(survival) library(survminer) library(timeROC) riskFile="risk.txt" #风险输入文件 cliFile="clinical.txt" #临床数据文件 setwd("C:\Users\Administrator\Desktop\geneimmune\30ROC") #修改工作目录 #读取风险输入文件 risk=read.table(riskFile, header=T, sep="t", check.names=F, row.names=1) risk=risk[,c("futime", "fustat", "riskScore")] #读取临床数据文件 cli=read.table(cliFile, header=T, sep="t", check.names=F, row.names=1) #合并数据 samSample=intersect(row.names(risk), row.names(cli)) risk1=risk[samSample,,drop=F] cli=cli[samSample,,drop=F] rt=cbind(risk1, cli) #定义颜色 bioCol=rainbow(ncol(rt)-1, s=0.9, v=0.9) #绘制ROC曲线 predictTime=3 #定义预测年限 aucText=c() pdf(file="ROC.pdf", width=6, height=6) #绘制风险得分的ROC曲线 i=3 ROC_rt=timeROC(T=rt$futime, delta=rt$fustat, marker=rt[,i], cause=1, weighting='aalen', times=c(predictTime),ROC=TRUE) plot(ROC_rt, time=predictTime, col=bioCol[i-2], title=FALSE, lwd=2) aucText=c(paste0("Risk", ", AUC=", sprintf("%.3f",ROC_rt$AUC[2]))) abline(0,1) #对临床数据进行循环，绘制临床数据的ROC曲线 for(i in 4:ncol(rt)){ ROC_rt=timeROC(T=rt$futime, delta=rt$fustat, marker=rt[,i], cause=1, weighting='aalen', times=c(predictTime),ROC=TRUE) plot(ROC_rt, time=predictTime, col=bioCol[i-2], title=FALSE, lwd=2, add=TRUE) aucText=c(aucText, paste0(colnames(rt)[i],", AUC=",sprintf("%.3f",ROC_rt$AUC[2]))) } #绘制联合的ROC曲线 multiCox=coxph(Surv(futime, fustat) ~ ., data = rt) score=predict(multiCox, type="risk", newdata=rt) ROC_rt=timeROC(T=rt$futime, delta=rt$fustat, marker=score,cause=1, weighting='aalen', times=c(predictTime),ROC=TRUE) plot(ROC_rt, time=predictTime, col=bioCol[ncol(rt)-1], title=FALSE, lwd=2, add=TRUE) aucText=c(aucText, paste0("Risk+Clinical", ", AUC=", sprintf("%.3f",ROC_rt$AUC[2]))) #绘制图例，得到ROC曲线下的面积 legend("bottomright", aucText,lwd=2,bty="n",col=bioCol[1:(ncol(rt)-1)]) dev.off()
17、绘制列线图与校准曲线
#install.packages("rms") library(rms) #引用包 riskFile="risk.txt" #风险输入文件 cliFile="clinical.txt" #临床数据文件 setwd("C:\Users\Administrator\Desktop\生信文章\geneimmune\31Nomo") #修改工作目录 #读取风险输入文件 risk=read.table(riskFile, header=T, sep="t", check.names=F, row.names=1) risk=risk[,c("futime", "fustat", "riskScore")] #读取临床数据文件 cli=read.table(cliFile, header=T, sep="t", check.names=F, row.names=1) #合并数据 samSample=intersect(row.names(risk), row.names(cli)) risk1=risk[samSample,,drop=F] cli=cli[samSample,,drop=F] rt=cbind(risk1, cli) paste(colnames(rt)[3:ncol(rt)],collapse="+") #数据打包 dd <- datadist(rt) options(datadist="dd") #生成函数 f <- cph(Surv(futime, fustat) ~ riskScore+Age+Gender+Grade+Stage+T+M+N, x=T, y=T, surv=T, data=rt, time.inc=1) surv <- Survival(f) #建立nomogram nom <- nomogram(f, fun=list(function(x) surv(1, x), function(x) surv(2, x), function(x) surv(3, x)), lp=F, funlabel=c("1-year survival", "2-year survival", "3-year survival"), maxscale=100, fun.at=c(0.99, 0.9, 0.8, 0.7, 0.5, 0.3,0.1,0.01)) #nomogram可视化 pdf(file="Nomogram.pdf",height=8.5,width=9.5) plot(nom) dev.off() #calibration curve time=3 #预测年限 f <- cph(Surv(futime, fustat) ~ riskScore+Age+Gender+Grade+Stage+T+M+N, x=T, y=T, surv=T, data=rt, time.inc=time) cal <- calibrate(f, cmethod="KM", method="boot", u=time, m=75, B=1000) pdf(file="calibration.pdf", width=9.5, height=8.5) plot(cal, xlab=paste0("Nomogram-Predicted Probability of ", time, "-Year OS"), ylab=paste0("Actual ", time, "-Year OS(proportion)"), col="red", sub=T) dev.off()
18、最后在ICGC肿瘤数据库中再次验证该模型的准确性，代码与以上类似
三、总结：
这里只展示了研究的部分内容，有部分研究结果是使用从在线数据库分析获取的，比如中性粒细胞和CD8阳性T细胞的变化情况。研究整体说明的问题是抑制HSP90β可能会通过调节筛选出的免疫基因来改善肝细胞癌患者的预后，后续可以增加部分验证实验来证明研究结果，通过抑制HSP90β来研究不同免疫基因的改变情况。更详细的研究内容可通过以下访问链接获取：
Gitee码云：
Github：
![洛谷 P11345 [KTSC 2023 R2] 基地简化 题解](http://www.itfaba.com/wp-content/themes/kemi/timthumb.php?src=http://www.itfaba.com/wp-content/themes/kemi/img/random/1.jpg&w=218&h=124&zc=1)
