(一)
技术领域：
：本发明涉及一种ii期结直肠癌复发预测模型的构建方法及预后标记基因。(二)
背景技术：
：：结直肠癌的治疗采用以外科手术治疗为主，同时辅以化学治疗、放射治疗、靶向治疗、免疫治疗等综合的治疗方式。对于早期结直肠癌患者，高质量的结直肠癌根治性手术能为患者带来显著的获益。但是部分早期结直肠癌患者根治性手术后出现局部复发、转移性复发，这些复发患者的预后往往较差，并且有研究发现结直肠癌患者根治手术后的rfs(无复发生存期)越短，其总生存期也越短。因此，早期结直肠癌根治手术后复发风险的预测具有重要意义，将有助于指导结直肠癌患者手术后辅助治疗方案的制定。ii期结直肠癌由于患者的预后差异较大，部分iib、iic期患者的预后甚至差于iiia期患者，一直是近年来的关注热点。nccn(美国国立综合癌症网络)指南提出ii期结直肠癌的8个高危因素：①t4期肿瘤②肿瘤伴穿孔③肿瘤伴梗阻④淋巴管、血管侵犯⑤神经侵犯⑥检出淋巴结<12个⑦肿瘤为低分化或未分化⑧切缘阳性。同时近年来研究发现msi-h(微卫星高度不稳定)或dmmr(错配修复蛋白缺失)是ii期结直肠癌的低危因素。具有任一高危因素的ii期结直肠癌应被视为高危患者，而具有低危因素的ii期结直肠癌患者应被视为低危患者，若既不属于高危患者也不属于低危患者应被视为普危患者。通常认为，低危ii期患者和普危ii期患者预后较好，所以辅助化疗的获益很小，而高危ii期患者预后较差，术后复发风险较高，可以从辅助化疗中获益。然而目前有许多临床研究发现，现在对高危ii期的定义并不准确，许多高危ii期患者并无复发，而一些普危ii期患者却发生了复发转移，这可能与传统的高危因素仅考虑了患者的临床病理特征，没有考虑肿瘤本身的生物学特征有关，而目前的基因芯片技术、高通量测序技术可以让研究者们更好的挖掘肿瘤的基因表达信息，进而反映肿瘤的生物学特征。(三)技术实现要素：本发明目的是根据ii期结直肠癌患者肿瘤组织的基因表达信息，鉴定ii期结直肠癌复发相关基因，并采用机器学习方法构建了ii期结直肠癌复发预测模型。本发明采用的技术方案是：本发明提供一种用于ii期结直肠癌复发预测的预后标记基因，所述预后标记基因包括nvl、acer、jup、plat、junb、igdcc3、ankrd27、nup50、git1、prkcdbp、tdo2、loc100287896、mcmbp、sned1、esm1、marco、fam46c、lmo2、haus8、ttc13、adnp2、rsbn1l、raver2、snora21、pnrc1和pja1。本发明还提供一种所述预后标记在构建ii期结直肠癌复发预测模型中的应用，所述模型在r语言程序中的表达式：v.obj_2ndscore＝pred(v.obj_2nd,data)，其中v.obj_2nd为根据ii期结直肠癌患者复发相关的预后标记基因的表达值，采用随机生存森林模型方法构建ii期结直肠癌患者肿瘤复发预测模型，v.obj_2ndscore为肿瘤复发风险评分。本发明所述ii期结直肠癌复发预测模型的构建方法，所述方法包括以下步骤：(1)获取基因表达数据集：获取ii期结直肠癌肿瘤样本，获取各个样本的mrna表达数据，其检测技术包括但不限于基因芯片技术、高通量转录组测序技术、实时荧光定量qpcr技术等，并通过聚类分析检测离群值并排除离群值；(2)鉴定与ii期结直肠癌患者复发相关的公共差异表达基因：采用基因芯片荟萃分析技术，以logrank检验计算每一个数据集中每一个基因的p值；然后，采用minp方法合并每一个数据集中每一个基因的p值，最终得到每一个基因的pminp值；采用benjamini-hochberg方法对每一个基因的pminp值进行校正，得到每一个基因的fdr值(错误发现率)，根据fdr<0.1的标准筛选与ii期结直肠癌患者复发相关的公共差异表达基因(即预后标记基因)；(3)差异表达基因的信号通路富集分析：采用metascape数据库(/gp/index.html#/main/step1)对步骤(2)鉴定得到的与ii期结直肠癌患者复发相关的公共差异表达基因进行信号通路富集分析，并以p＝0.01作为阈值筛选差异表达基因显著富集的信号通路，即p值小于0.01的信号通路；(4)模型的构建：根据步骤(2)ii期结直肠癌患者复发相关的公共差异表达基因的表达值，采用随机生存森林模型方法构建ii期结直肠癌患者肿瘤复发预测模型v.obj_2nd，所述模型在r语言程序中的表达式：v.obj_2ndscore＝pred(v.obj_2nd,data)其中，v.obj_2ndscore为肿瘤复发风险评分，本发明所述模型中v.obj_2ndscore的大小即代表着肿瘤复发风险概率的大小，通过r语言程序的stats包的predict函数，输入模型(v.obj_2nd)和患者的基因表达数据(data)，即可获得患者肿瘤复发风险评分。根据ii期结直肠癌患者的肿瘤复发风险评分的中位数(优选中位数为6.10)，将患者分为高风险组(复发风险评分高于中位数)与低风险组(复发风险评分低于中位数)。与现有技术相比，本发明有益效果主要体现在：本发明提供了一种根据ii期结直肠癌患者肿瘤的基因表达信息预测患者复发风险的模型，该模型预测ii期结直肠癌患者5年复发风险的auc值为0.993，在测试集中能将高复发风险和低复发风险的患者显著分开(hr＝1.901，95％ci＝1.124-3.215)；关于模型的建立与选择，本发明使用了随机生存森林模型，并且根据变量所在的最大子树的最小深度值进行变量筛选，挑选重要的变量重新建立模型，使模型的变量维度大幅度降低；在测试集中使用随机生存森林模型将患者分成高复发风险组与低复发风险组后，两组间患者的无复发生存期具有显著性差异。(四)附图说明图1：数据集筛选流程。图2：gse14333数据集的层次聚类。图3：差异表达基因显著富集的信号通路。柱状图颜色反映p值大小，颜色越深，p值越小；柱状图右侧为信号通路名称。图4：随机生存森林模型v.obj_2nd的调优参数组合。横坐标为终末节点的数据点个数，纵坐标为每个节点的候选变量数目，颜色深度反映袋外错误率大小，颜色越深，袋外错误率越大。图5：随机生存森林模型v.obj_2nd中基因表达值与患者无复发生存期之间的非线性关系。横坐标为每个基因表达值，纵坐标为患者的无复发生存率。图6：训练集的时间依赖的roc曲线。分别为1年、3年、5年时，随机生存森林模型的时间依赖的roc曲线图。图7：测试集中根据随机生存森林模型v.obj_2nd预测的高、低复发风险组的生存曲线。线条b为根据lassocox模型预测的高复发风险组；线条a为根据lassocox模型预测的低复发风险组。(五)具体实施方式下面结合具体实施例对本发明进行进一步描述，但本发明的保护范围并不仅限于此：实施例1基于ii期结直肠癌患者肿瘤组织基因表达信息的肿瘤复发预测模型(1)获取基因表达数据集基因表达数据指多个个体样本的mrna表达数据表达数据，其检测技术包括但不限于基因芯片技术、高通量转录组测序技术、实时荧光定量qpcr技术等。通过检索ncbi(美国国立生物技术信息中心)的高通量基因表达数据库(https://www.ncbi.nlm.nih.gov/gds/)获取基因表达数据集，检索式为：("colorectalcancer"[allfields]or"coloncancer"[allfields]or"rectalcancer"[allfields])and"expressionprofilingbyarray"[filter]，共获取981个基因表达数据集，并根据纳入标准与排除标准对数据集进行筛选。数据集的纳入标准为：(1)样本类型为ii期结直肠癌患者手术后肿瘤组织样本；(2)患者的rfs数据可获取；(3)检测技术为基因表达谱芯片。以上3项标准全部满足的数据集将被纳入后续分析。数据集的排除标准为：(1)样本类型非ii期结直肠癌患者手术后肿瘤组织样本；(2)仅提供患者复发与否信息，而未提供患者的rfs信息；(3)数据集样本量小于40。最终纳入分析的样本为6个数据集(gse14333、gse17538、gse33113、gse39582、gse24551)中的651例ii期结直肠癌患者手术后肿瘤组织样本，数据集的筛选流程见图1。以上3项标准中任意1项不满足的数据集将被排除。下载上述651例样本的基因芯片原始数据，采用r语言程序的oligo包读取数据并采用rma(robustmultiarrayaverage)算法进行标准化预处理(目的是消除由于实验技术等其他因素所导致的基因表达量的变化，并且使各个样本的基因表达数据处于可比的水平)，获取各个样本的基因表达矩阵，基因表达矩阵横坐标为每个样本的编号，纵坐标为探针集的编号，其r语言程序代码如下(以数据集gse14333为例，其余数据集只需要将gse14333替换为其他数据集的名称即可)。为了得到每个样本的基因表达值，采用david数据库(http://david.ncifcrf.gov/)将探针集的编号转换为基因名，得到每个数据集的基因表达矩阵。library(oligo)library(pd.hg.u133.plus.2)genecels＝list.celfiles(listgzipped＝t,full.name＝t)affygenefs<-read.celfiles(genecels)genecore<-rma(affygenefs)write.csv([email protected]$exprs,file＝"gse14333.csv")由于实验误差等因素存在，基因芯片数据可能存在离群值，离群值的出现将为后续建立模型、测试模型带来困难。根据基因表达矩阵计算样本间的pearson相关系数，并采用(1-pearson相关系数)/2作为样本间的距离进行层次聚类(采用r语言程序的nbclust包的hclust函数进行聚类)，以聚类的高度0.2作为阈值检测离群值，其中一个数据集gse14333中样本的层次聚类图见图2。本实施例中6个数据集的651例样本的基因表达矩阵都通过了离群值检测，并被纳入下一步分析，离群值检测的r语言程序代码如下(以数据集gse14333为例，其余数据集只需要将gse14333替换为其他数据集的名称即可)。library(nbclust)gse14333＝read.table("gse14333.txt",sep＝"\t",header＝true)gse14333＝gse14333[-c(1,2),]row.names(gse14333)＝gse14333[,1]gse14333＝gse14333[,-c(1,2)]dd<-as.dist((1-cor(gse14333))/2)cluster＝hclust(dd)tiff(filename＝"gse14333.gif",width＝4500,height＝3500,res＝300)plot(cluster,main＝"clusterdendrogramofgse14333",xlab＝na)dev.off()cutree(cluster,h＝0.2)(2)鉴定与ii期结直肠癌患者复发相关的公共差异表达基因传统的肿瘤患者复发相关基因鉴定、复发预测模型建立往往基于单个数据集，而这种方法具有样本量缺乏、差异表达基因的可重复性较差等缺点，而本发明的创新点之一在于采用基因芯片荟萃分析技术，鉴定与ii期结直肠癌患者复发相关的公共差异表达基因。采用r语言程序的metade包读取各个样本的基因表达矩阵，将患者的rfs时间、复发状态与基因表达矩阵进行匹配，并进行后续的基因芯片荟萃分析。由于基因表达矩阵存在多个探针集编号对应一个基因名的情况，采用r语言程序metade包的metade.match函数根据探针集编号的iqr值(四分位数间距值)筛选探针集，并选取同一个基因名iqr值最大的探针集代表该基因的表达值。采用合并p值的方法进行基因芯片荟萃分析，该方法的实现依赖于r语言程序metade包的metade.rawdata函数。首先，采用logrank检验计算每一个数据集中每一个基因的p值；然后，采用minp方法合并每一个数据集中每一个基因的p值，最终得到每一个基因的pminp值。由于对每一个基因都采用了一次假设检验，为了减少假设检验过程中的α错误，采用benjamini-hochberg方法对每一个基因的pminp值进行校正，得到每一个基因的fdr值(错误发现率)。根据fdr<0.1的标准筛选与ii期结直肠癌患者复发相关的公共差异表达基因，共鉴定得到479个与ii期结直肠癌患者复发相关的公共差异表达基因。其r语言程序代码如下：rm(list＝ls())library(metade)study.names<-c("gse14333","gse17538","gse24551","gse33113","gse39582","gse92921")crc.raw<-metade.read(study.names,skip＝rep(2,6),via＝"txt",matched＝false,log＝false)crc.raw_match＝metade.match(crc.raw,pool.replicate＝"iqr")dim(crc.raw[[1]][[1]])crc.merged<-metade.merge(crc.raw_match)dim(crc.merged[[1]][[1]])crc.filtered<-metade.filter(crc.merged,c(0.2,0))dim(crc.filtered[[1]][[1]])crc.filtered$gse14333$censoring.status＝crc.raw$gse14333$censoring.statuscrc.filtered$gse17538$censoring.status＝crc.raw$gse17538$censoring.statuscrc.filtered$gse24551$censoring.status＝crc.raw$gse24551$censoring.statuscrc.filtered$gse33113$censoring.status＝crc.raw$gse33113$censoring.statuscrc.filtered$gse39582$censoring.status＝crc.raw$gse39582$censoring.statuscrc.filtered$gse92921$censoring.status＝crc.raw$gse92921$censoring.statusstart＝sys.time()set.seed()metade.res.minp<-metade.rawdata(crc.filtered,ind.method＝rep("logrank",6),meta.method＝c("minp"),paired＝rep(f,7),asymptotic＝t)b<-sys.time()-startprint(b)view(metade.res.minp$meta.analysis$fdr)count.denumber(metade.res.minp,p.cut＝c(0.01,0.05),q.cut＝c(0.1))draw.denumber(metade.res.minp,1,fdr＝t)write.table(metade.res.minp$meta.analysis$fdr,file＝"deg.txt",sep＝"\t")(3)差异表达基因的信号通路富集分析基因富集分析常用于对基因芯片、转录组测序等高通量技术获得的差异表达基因进行功能注释、通路注释，进而有助于探索生物学表型、临床疾病内在的分子机制。采用metascape数据库(/gp/index.html#/main/step1)对鉴定得到的479个与ii期结直肠癌患者复发相关的公共差异表达基因进行信号通路富集分析，并以p＝0.01作为阈值筛选差异表达基因显著富集的信号通路，见图3。(4)在训练集中构建ii期结直肠癌患者肿瘤复发预测模型为了构建ii期结直肠癌患者肿瘤复发预测模型，将6个数据集的651例样本按照60％、40％的比例随机分为训练集、测试集，并在训练集中构建肿瘤复发预测模型。由于基因表达数据具有高数据维度、高相关性等特征，传统的cox回归模型不适用于高维数据，而随机生存森林模型可以很好的解决这一问题。除此之外，基因表达值对患者复发风险的影响可能是非线性的，传统的cox回归模型本质上是线性模型，可能会降低模型的预测精度，而随机生存森林可以发现基因表达值对患者复发风险的非线性效应。关于随机生存森林模型的选择也是本发明的创新点之一。随机生存森林模型的重要参数为每个节点的候选变量数目与终末节点的数据点个数，通过计算各种参数组合条件下训练集中的袋外错误率，确定使袋外错误率最低的参数组合，见图4。随机生存森林模型可以根据变量所在的最大子树的最小深度值进行变量筛选，其筛选阈值为所有变量所在的最大子树的最小深度值的平均值。经过变量筛选后，还可以利用筛选得到的变量重复上述模型构建步骤，直至所有变量所在的最大子树的最小深度值在阈值以下。经过2次随机生存森林模型的变量筛选，最终得到26个基因表达值构成的预测模型v.obj_2nd，所使用的软件为r语言程序的randomforestsrc包，具体r语言程序代码如下：library("randomforestsrc")library("survival")library("timeroc")data＝rbind(gse14333_bn,gse17538_bn,gse24551_bn,gse33113_bn,gse39582_bn,gse92921_bn)data＝data.frame(data)colnames(data)[480]＝"time"set.seed()random_sample＝sample(651)train_index＝random_sample[1:390]test_index＝random_sample[391:651]train_data＝data[train_index,]test_data＝data[test_index,]tune_data＝tune(surv(time,status)～.,data＝train_data,ntreetry＝1000,trace＝t)print(tune_data$rf)v.obj<-rfsrc(surv(time,status)～.,data＝train_data,nodesize＝2,mtry＝34,ntree＝1000,importance＝true)print(v.obj)plot(v.obj)pred＝predict(v.obj,train_data)pred#1stsimplify#v.max$topvarstrain_data_1st＝train_data[,c(v.max$topvars,"time","status")]test_data_1st＝test_data[,c(v.max$topvars,"time","status")]tune_data_1st＝tune(surv(time,status)～.,data＝train_data_1st,ntreetry＝1000,trace＝t)print(tune_data_1st$rf)v.obj_1st<-rfsrc(surv(time,status)～.,data＝train_data_1st,nodesize＝1,mtry＝11,ntree＝1000,importance＝true)print(v.obj_1st)plot(v.obj_1st)plot.survival(v.obj_1st,cens.model＝"rfsrc")pred_1st＝predict(v.obj_1st,train_data_1st)pred_1st#2ndsimplify#v.max_1st$topvarstrain_data_2nd＝train_data[,c(v.max_1st$topvars,"time","status")]test_data_2nd＝test_data[,c(v.max_1st$topvars,"time","status")]tune_data_2nd＝tune(surv(time,status)～.,data＝train_data_2nd,ntreetry＝1000,trace＝t)print(tune_data_2nd$rf)v.obj_2nd<-rfsrc(surv(time,status)～.,data＝train_data_2nd,nodesize＝2,mtry＝10,ntree＝1000,importance＝true)print(v.obj_2nd)plot(v.obj_2nd)plot.variable(v.obj_2nd,surv.type＝"surv",time＝60,plots.per.page＝6)pred_2nd＝predict(v.obj_2nd,train_data_2nd)pred_2nd#minimal_depthv.max_2nd<-max.subtree(v.obj_2nd)print(round(v.max_2nd$order,3))#theminimaldepthisthefirstorderdepthprint(round(v.max_2nd$order[,1],3))minimal_depth_2nd＝round(v.max_2nd$order[,1],3)print(v.max_2nd$threshold)v.max_2nd$topvars#test#pred_test_2nd＝predict(v.obj_2nd,test_data_2nd)#tdroc#roc_train_2nd＝timeroc(t＝train_data_2nd$time,delta＝train_data_2nd$status,marker＝pred_2nd$predicted,cause＝1,iid＝true,times＝c(12,36,60))plot(roc_train_2nd,time＝12,add＝f)plot(roc_train_2nd,time＝36,add＝f)plot(roc_train_2nd,time＝60,add＝f)所述模型在r语言程序中的表达式：v.obj_2ndscore＝pred(v.obj_2nd,data)其中，v.obj_2ndscore为肿瘤复发风险评分，本发明所述模型中v.obj_2ndscore的大小即代表着肿瘤复发风险概率的大小，通过r语言程序的stats包的predict函数，输入模型(v.obj_2nd)和患者的基因表达数据(data)，即可获得患者肿瘤复发风险评分。v.obj_2nd模型中26个基因的名称及其所在的最大子树的最小深度值见表1。表1.随机生存森林模型v.obj_2nd中26个基因的最大子树的最小深度值通过变量的边缘效应图展示了基因表达值与患者无复发生存期之间的非线性关系，见图5。在训练集中采用随机生存森林模型v.obj_2nd构建肿瘤复发预测模型后，采用时间依赖的roc曲线(受试者工作曲线)的auc(曲线下面积)评估模型的预测效果。auc的取值范围为0到1之间，auc越大说明模型的预测效果越好。随机生存森林模型v.obj_2nd的auc1年＝0.995、auc3年＝0.999、auc5年＝0.993，见图6。(5)在测试集中评估模型预测效果获得了测试集中所有样本的基因表达矩阵，并提取随机生存森林模型v.obj_2nd所纳入的26个基因表达值(nvl、acer、jup、plat、junb、igdcc3、ankrd27、nup50、git1、prkcdbp、tdo2、loc100287896、mcmbp、sned1、esm1、marco、fam46c、lmo2、haus8、ttc13、adnp2、rsbn1l、raver2、snora21、pnrc1、pja1)。将构建的随机生存森林预测模型v.obj_2nd导入r语言程序，调用r语言程序的predict函数并将模型、测试集中样本的上述26基因表达矩阵作为输入参数，即可获得模型在测试集中的肿瘤复发风险评分。根据ii期结直肠癌患者的肿瘤复发风险评分的中位数(在本实施例中中位数为6.10)，将患者分为高风险组(复发风险评分高于6.10)与低风险组(复发风险评分低于6.10)。分别绘制两组患者的kaplan-meier生存曲线，采用logrank检验两组间的无复发生存期是否具有显著差异，并采用mantel-haenszel方法计算hr值(风险比)及其95％ci(置信区间)。测试集按照上述方法分成高复发风险组与低复发风险组后，两组间logrank检验的p值＝0.0165，高复发风险组相比于低复发风险组的hr＝1.901，95％ci＝1.124-3.215，见图7。以上所述了本发明的一个实施例，本领域的普通技术人员可以理解，应当指出，在不脱离本发明的原理和宗旨的情况下可以对这些实施例、方法学、选取的模型、编程软件进行多种变化、修改、替换和补充，这些变化、修改、替换和补充也应该视为本发明的保护范围。当前第1页1 2 3 当前第1页1 2 3
技术特征：
1.一种用于ii期结直肠癌复发预测的预后标记基因，其特征在于所述标记基因包括nvl、acer、jup、plat、junb、igdcc3、ankrd27、nup50、git1、prkcdbp、tdo2、loc100287896、mcmbp、sned1、esm1、marco、fam46c、lmo2、haus8、ttc13、adnp2、rsbn1l、raver2、snora21、pnrc1和pja1。
2.一种权利要求1所述预后标记在构建ii期结直肠癌复发预测模型中的应用。
3.如权利要求2所述的应用，其特征在于所述模型在r语言程序中的表达式：v.obj_2ndscore＝pred(v.obj_2nd,data)，其中v.obj_2nd为根据ii期结直肠癌患者复发相关的预后标记基因的表达值，采用随机生存森林模型方法构建ii期结直肠癌患者肿瘤复发预测模型，v.obj_2ndscore为肿瘤复发风险评分。
4.如权利要求2所述的应用，其特征在于模型的构建方法包括以下步骤：(1)获取基因表达数据集：获取ii期结直肠癌肿瘤样本，获取各个样本的mrna表达数据，并通过聚类分析检测离群值并排除离群值；(2)鉴定与ii期结直肠癌患者复发相关的公共差异表达基因：采用基因芯片荟萃分析技术，以logrank检验计算每一个数据集中每一个基因的p值；然后，采用minp方法合并每一个数据集中每一个基因的p值，最终得到每一个基因的pminp值；采用benjamini-hochberg方法对每一个基因的pminp值进行校正，得到每一个基因的fdr值，根据fdr<0.1的标准筛选与ii期结直肠癌患者复发相关的公共差异表达基因；(3)差异表达基因的信号通路富集分析：采用metascape数据库对步骤(2)鉴定得到的与ii期结直肠癌患者复发相关的公共差异表达基因进行信号通路富集分析，并以p＝0.01作为阈值筛选差异表达基因显著富集的信号通路，即p值小于0.01的信号通路；(4)模型的构建：根据步骤(2)ii期结直肠癌患者复发相关的公共差异表达基因的表达值，采用随机生存森林模型方法构建ii期结直肠癌患者肿瘤复发预测模型v.obj_2nd，
所述模型在r语言程序中的表达式：v.obj_2ndscore＝pred(v.obj_2nd,data)。
技术总结
本发明公开了一种用于II期结直肠癌复发预测的预后标记基因及应用，本发明提供了一种根据II期结直肠癌患者肿瘤的基因表达信息预测患者复发风险的模型，该模型预测II期结直肠癌患者5年复发风险的AUC值为0.993，在测试集中能将高复发风险和低复发风险的患者显著分开；关于模型的建立与选择，本发明使用了随机生存森林模型，并且根据变量所在的最大子树的最小深度值进行变量筛选，挑选重要的变量重新建立模型，使模型的变量维度大幅度降低；在测试集中使用随机生存森林模型将患者分成高复发风险组与低复发风险组后，两组间患者的无复发生存期具有显著性差异。
技术研发人员：丁克峰;陆玮;肖乾;李军
受保护的技术使用者：浙江大学
技术研发日：.09.29
技术公布日：.02.14
