CN110714078A

CN110714078A - 一种用于ii期结直肠癌复发预测的标记基因及应用

Info

Publication number: CN110714078A
Application number: CN201910932402.8A
Authority: CN
Inventors: 丁克峰; 陆玮; 肖乾; 李军
Original assignee: Zhejiang University ZJU
Current assignee: Zhejiang University ZJU
Priority date: 2019-09-29
Filing date: 2019-09-29
Publication date: 2020-01-21
Anticipated expiration: 2039-09-29
Also published as: CN110714078B

Abstract

本发明公开了一种用于II期结直肠癌复发预测的标记基因及应用，本发明采用基因芯片荟萃分析技术鉴定与II期结直肠癌患者复发相关的公共差异表达基因，该模型预测II期结直肠癌患者5年复发风险的AUC值为0.806，在测试集中能将高复发风险和低复发风险的患者显著分开；本发明使用了Lasso Cox回归模型，既起到了建立模型的作用又起到了变量重要性筛选的作用，使模型的变量维度大幅度降低，将有助于降低基因表达检测的成本，有利于该模型在临床应用中的推广。

Description

一种用于II期结直肠癌复发预测的标记基因及应用

(一)技术领域

本发明涉及生物信息学领域，具体涉及采用基因芯片荟萃分析技术鉴定II期结直肠癌复发相关的公共差异表达基因技术，构建一种基于公共差异表达基因与Lasso Cox模型建立的II期结直肠癌患者肿瘤复发预测模型。

(二)背景技术

结直肠癌是全世界最常见的恶性肿瘤之一，在各类恶性肿瘤中其发病率排名第三，死亡率排名第四。近几年来，随着人民经济水平的提高、生活方式的改变，我国的结直肠癌发病率正逐渐升高并且呈年轻化趋势，每年新发的结直肠癌患者数目增幅大约在4％左右。

结直肠癌的治疗采用以外科手术治疗为主，同时辅以化学治疗、放射治疗、靶向治疗、免疫治疗等综合的治疗方式。对于早期结直肠癌患者，高质量的结直肠癌根治性手术能为患者带来显著的获益。但是部分早期结直肠癌患者根治性手术后出现局部复发、转移性复发，这些复发患者的预后往往较差，并且有研究发现结直肠癌患者根治手术后的RFS(无复发生存期)越短，其总生存期也越短。因此，早期结直肠癌根治手术后复发风险的预测具有重要意义，将有助于指导结直肠癌患者手术后辅助治疗方案的制定。

II期结直肠癌由于患者的预后差异较大，部分IIB、IIC期患者的预后甚至差于IIIA期患者，一直是近年来的关注热点。NCCN(美国国立综合癌症网络)指南提出II期结直肠癌的8个高危因素：①T4期肿瘤②肿瘤伴穿孔③肿瘤伴梗阻④淋巴管、血管侵犯⑤神经侵犯⑥检出淋巴结<12个⑦肿瘤为低分化或未分化⑧切缘阳性。同时近年来研究发现MSI-H(微卫星高度不稳定)或dMMR(错配修复蛋白缺失)是II期结直肠癌的低危因素。具有任一高危因素的II期结直肠癌应被视为高危患者，而具有低危因素的II期结直肠癌患者应被视为低危患者，若既不属于高危患者也不属于低危患者应被视为普危患者。通常认为，低危II期患者和普危II期患者预后较好，所以辅助化疗的获益很小，而高危II期患者预后较差，术后复发风险较高，可以从辅助化疗中获益。然而目前有许多临床研究发现，现在对高危II期的定义并不准确，许多高危II期患者并无复发，而一些普危II期患者却发生了复发转移，这可能与传统的高危因素仅考虑了患者的临床病理特征，没有考虑肿瘤本身的生物学特征有关，而目前的基因芯片技术、高通量测序技术可以让研究者们更好的挖掘肿瘤的基因表达信息，进而反映肿瘤的生物学特征。

(三)发明内容

本发明的目的是根据II期结直肠癌患者肿瘤组织的基因表达信息，采用基因芯片荟萃分析技术鉴定II期结直肠癌复发相关的公共差异表达基因，并采用Lasso Cox模型建立了II期结直肠癌患者复发预测模型。

本发明采用的技术方案为：

本发明提供一种用于II期结直肠癌复发预测的标记基因，所述标记基因包括：PAOX、SIGLEC7、PHAX、XCR1、TM4SF4、TRIOBP、MCMBP、HCFC1R1、ADNP2、NUP50、GTF2A2、BCCIP、FLJ90680、NVL、ESM1、GABRR2、FAM166A、USP14、JUNB、UBAP2、AP5B1、FAM46C、LDB3和JUP。

本发明还提供一种所述标记基因在构建II期结直肠癌复发预测模型中的应用，所述模型的数学表达式：lasso_coxscore＝∑(基因表达值*回归系数)。

所述基因表达值为标记基因mRNA表达值。

所述回归系数见表1所示：

表1.lasso_cox回归模型中24个基因的回归系数

本发明所述II期结直肠癌复发预测模型的构建方法：(1)获取基因表达数据集：获取II期结直肠癌肿瘤样本的mRNA表达数据，其检测技术包括但不限于基因芯片技术、高通量转录组测序技术、实时荧光定量qPCR技术等，并通过聚类分析检测离群值并排除离群值；(2)鉴定与II期结直肠癌患者复发相关的公共差异表达基因：采用基因芯片荟萃分析技术，以Logrank检验计算步骤(1)每一个数据集中每一个基因的p值；然后，采用minP方法合并每一个数据集中每一个基因的p值，最终得到每一个基因的p_minP值；采用Benjamini-Hochberg方法对每一个基因的p_minP值进行校正，得到每一个基因的FDR值(错误发现率)；根据FDR<0.1的标准筛选与II期结直肠癌患者复发相关的公共差异表达基因；(3)公共差异表达基因的信号通路富集分析：采用Metascape数据库(http://metascape.org/gp/index.html#/main/step1)对步骤(2)鉴定得到的与II期结直肠癌患者复发相关的公共差异表达基因进行信号通路富集分析，并以p＝0.01作为阈值筛选差异表达基因显著富集的信号通路，即p值小于0.01的信号通路；(4)模型构建：根据步骤(2)II期结直肠癌患者复发相关的公共差异表达基因的表达值，采用Lasso Cox方法构建II期结直肠癌患者肿瘤复发预测lasso_cox模型；

所述模型的数学表达式：lasso_coxscore＝∑(基因表达值*回归系数)

所述模型在R语言程序中的表达式：lasso_cox score＝predict(lasso_cox,data)其中，基因表达值为步骤(2)II期结直肠癌患者复发相关的公共差异表达基因的表达值，回归系数为lasso_cox模型中的回归系数，采用R语言程序的glmnet包中的glmnet函数计算得到。lasso_cox score为肿瘤复发风险评分，本发明所述模型中lasso_cox score评分的大小即代表着肿瘤复发风险概率的大小，通过R语言程序的stats包的predict函数，输入模型(lasso_cox)和患者的基因表达数据(data)，即可获得患者肿瘤复发风险评分。根据II期结直肠癌患者的肿瘤复发风险评分的中位数(优选-2.748)，将患者分为高风险组(复发风险评分高于中位数)与低风险组(复发风险评分低于中位数)。

与现有技术相比，本发明有益效果主要体现在：本发明的创新点之一在于采用基因芯片荟萃分析技术，鉴定与II期结直肠癌患者复发相关的公共差异表达基因；与传统的仅根据II期结直肠癌患者的临床病理特征判断复发风险的方法相比，本发明提供了一种根据II期结直肠癌患者肿瘤的基因表达信息预测患者复发风险的模型，该模型预测II期结直肠癌患者5年复发风险的AUC值为0.806，在测试集中能将高复发风险和低复发风险的患者显著分开(HR＝2.052，95％CI＝1.219-3.455)；本发明使用了Lasso Cox回归模型，既起到了建立模型的作用又起到了变量重要性筛选的作用，使模型的变量维度大幅度降低，将有助于降低基因表达检测的成本，有利于该模型在临床应用中的推广。

(四)附图说明

图1：数据集筛选流程。

图2：GSE14333数据集的层次聚类。

图3：差异表达基因显著富集的信号通路。柱状图颜色反映p值大小，颜色越深，p值越小；柱状图右侧为信号通路名称。

图4：lasso_cox回归模型中正则化参数λ与部分似然估计偏差关系。横坐标为正则化参数λ的自然对数值，纵坐标为部分似然估计偏差

图5：训练集的时间依赖的ROC曲线。分别为1年、3年、5年时lasso_cox回归模型的时间依赖的ROC曲线图。

图6：测试集中根据lasso_cox模型预测的高、低复发风险组的生存曲线。线条b为根据lasso_cox模型预测的高复发风险组；线条a为根据lasso_cox模型预测的低复发风险组。

(五)具体实施方式

下面结合具体实施例对本发明进行进一步描述，但本发明的保护范围并不仅限于此：

实施例1、基于II期结直肠癌患者肿瘤组织基因表达信息的肿瘤复发预测模型的构建及验证

(1)、获取基因表达数据集

采用基因表达谱芯片技术获取基因表达数据集，即获取多个个体样本的mRNA表达值，具体为：

通过检索NCBI(美国国立生物技术信息中心)的高通量基因表达数据库(https://www.ncbi.nlm.nih.gov/gds/)获取基因表达数据集，检索式为：("colorectal cancer"[All Fields]OR"colon cancer"[All Fields]OR"rectal cancer"[All Fields])AND"Expression profiling by array"[Filter]，共获取981个基因表达数据集(mRNA)，并根据纳入标准与排除标准对数据集进行筛选。

数据集的纳入标准为：(1)样本类型为II期结直肠癌患者手术后肿瘤组织样本；(2)患者的RFS数据可获取；(3)检测技术为基因表达谱芯片。以上3项标准全部满足的数据集将被纳入后续分析。

数据集的排除标准为：(1)样本类型非II期结直肠癌患者手术后肿瘤组织样本；(2)仅提供患者复发与否信息，而未提供患者的RFS信息；(3)数据集样本量小于40。以上3项标准中任意1项不满足的数据集将被排除。

最终纳入分析的样本为6个数据集(GSE14333、GSE17538、GSE33113、GSE39582、GSE24551)，涉及651例II期结直肠癌患者手术后肿瘤组织样本，数据集的筛选流程见图1。

下载上述651例样本的基因芯片原始数据，采用R语言程序的oligo包读取数据，采用RMA(Robust Multiarray Average)算法对数据进行标准化预处理(其目的是消除由于实验技术等其他因素所导致的基因表达量的变化，并且使各个样本的基因表达数据处于可比的水平)，获得基因表达矩阵，基因表达矩阵横坐标为每个样本的编号，纵坐标为探针集的编号(基因芯片原始数据自带的探针集编号)。采用DAVID数据库(http://david.ncifcrf.gov/)将探针集编号转换成基因名，得到每个数据集的基因表达矩阵。

R语言程序代码如下(以数据集GSE14333为例，其余数据集只需要将GSE14333替换为其他数据集的名称即可)：

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(geneCore@assayData$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.tif",width＝4500,height＝3500,res＝300)

plot(cluster,main＝"Cluster dendrogram of GSE14333",xlab＝NA)

dev.off()

cutree(cluster,h＝0.2)

(2)鉴定与II期结直肠癌患者复发相关的公共差异表达基因

传统的肿瘤患者复发相关基因鉴定、复发预测模型建立往往基于单个数据集，而这种方法具有样本量缺乏、差异表达基因的可重复性较差等缺点，而本发明的创新点之一在于采用基因芯片荟萃分析技术，鉴定与II期结直肠癌患者复发相关的公共差异表达基因，即标记基因。采用R语言程序的MetaDE包读取基因表达矩阵，将患者的RFS时间、复发状态与基因表达矩阵进行匹配，并进行后续的基因芯片荟萃分析。由于基因表达矩阵存在多个探针集编号对应一个基因名的情况，采用R语言程序MetaDE包的MetaDE.match函数，根据探针集编号的IQR值(四分位数间距值，通过MetaDE包的MetaDE.match函数计算得到)筛选探针集，并选取同一个基因名IQR值最大的探针集代表该基因的表达值。

采用合并p值的方法进行基因芯片荟萃分析，该方法的实现依赖于R语言程序MetaDE包的MetaDE.rawdata函数。首先，采用Logrank检验计算每一个数据集中每一个基因的p值；然后，采用minP方法合并每一个数据集中每一个基因的p值，最终得到每一个基因的p_minP值。由于对每一个基因都采用了一次假设检验，为了减少假设检验过程中的α错误，采用Benjamini-Hochberg方法对每一个基因的p_minP值进行校正，得到每一个基因的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.status

crc.filtered$GSE17538$censoring.status＝crc.raw$GSE17538$censoring.status

crc.filtered$GSE24551$censoring.status＝crc.raw$GSE24551$censoring.status

crc.filtered$GSE33113$censoring.status＝crc.raw$GSE33113$censoring.status

crc.filtered$GSE39582$censoring.status＝crc.raw$GSE39582$censoring.status

crc.filtered$GSE92921$censoring.status＝crc.raw$GSE92921$censoring.status

start＝Sys.time()

set.seed(2018)

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()-start

print(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数据库(http://metascape.org/gp/index.html#/main/step1)对鉴定得到的479个与II期结直肠癌患者复发相关的公共差异表达基因进行信号通路富集分析，并以p＝0.01作为阈值筛选差异表达基因显著富集的信号通路，见图3。

(4)在训练集中构建II期结直肠癌患者肿瘤复发预测模型

为了构建II期结直肠癌患者肿瘤复发预测模型，将6个数据集的651例样本按照60％、40％的比例随机分为训练集、测试集，并在训练集中构建肿瘤复发预测模型。由于基因表达数据具有高数据维度、高相关性等特征，传统的Cox回归模型将无法适用，对于高维生存数据预测模型的选择包括但不限于Lasso Cox回归模型，采用Lasso Cox回归方法构建肿瘤复发预测lasso_cox模型，所使用的软件为R语言程序的glmnet包，具体R语言程序代码如下：

library("glmnet")

library("survival")

library("timeROC")

data＝rbind(GSE14333_BN,GSE17538_BN,GSE24551_BN,GSE33113_BN,GSE39582_BN,GSE92921_BN)

set.seed(2018)

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,]

x_train＝train_data[,1:479]

follow_up_train＝train_data[,480:481]

zero_follow_up_train＝(follow_up_train[,1]＝＝0)

follow_up_train[zero_follow_up_train,1]＝0.01

y_train＝Surv(follow_up_train[,1],follow_up_train[,2])

set.seed(2018)

fit1_cv＝cv.glmnet(x_train,y_train,family＝"cox",nfolds＝20)

plot(fit1_cv)

title("Lasso-Cox",line＝2.5)

cox＝coxph(y_train～x_train)

lasso_cox＝glmnet(x_train,y_train,family＝"cox",lambda＝fit1_cv$lambda.min)

lasso_cox$beta

ROC_train＝timeROC(T＝follow_up_train[,1],delta＝follow_up_train[,2],marker＝predict(lasso_cox,x_train),cause＝1,iid＝TRUE,times＝c(12,36,60))

plot(ROC_train,time＝12,add＝F)

plot(ROC_train,time＝36,add＝F)

plot(ROC_train,time＝60,add＝F)

ROC$AUC_train

x_test＝test_data[,1:479]

follow_up_test＝test_data[,480:481]

zero_follow_up_test＝(follow_up_test[,1]＝＝0)

follow_up_test[zero_follow_up_test,1]＝0.01

ROC_test＝timeROC(T＝follow_up_test[,1],delta＝follow_up_test[,2],marker＝predict(lasso_cox,x_test),cause＝1,iid＝TRUE,times＝c(12,36,60))

plot(ROC_test,time＝12,add＝F)

plot(ROC_test,time＝36,add＝F)

plot(ROC_test,time＝60,add＝F)

ROC_test$AUC

write.table(cbind(predict(lasso_cox,x_test),follow_up_test),file＝"test.csv",sep＝",")

Lasso Cox回归模型与传统的Cox回归模型相比，最大的不同在于Lasso Cox回归模型引入了回归系数的L1范数惩罚项的权重λ，又称之为正则化参数λ。通过调整参数λ值，可以使得某些变量的回归系数等于0(使所有基因中除了表1中24个基因之外的基因的回归系数等于0)，这样就达到了变量选择与简化模型的目的。

最优的λ值是根据在训练集中采用20折交叉验证的方法确定的，在该λ取值时模型的部分似然估计偏差达到最小值，见图4，并得出在该λ取值时479个与II期结直肠癌患者复发相关的公共差异表达基因的回归系数，其中455个基因表达值的回归系数等于0，24个基因表达值的回归系数不为0，这24个基因名及其表达值的回归系数见表1。

表1.lasso_cox回归模型中24个基因的回归系数

每个基因表达值的回归系数值表示该基因表达量每变化1个单位，II期结直肠癌患者肿瘤复发风险评分的改变值。若回归系数为正数，则表示该基因表达值升高时肿瘤复发风险增加；类似的，若回归系数为负数，则表示该基因表达值升高时肿瘤复发风险降低。肿瘤复发风险评分(即lasso_cox回归模型)数学计算公式为：肿瘤复发风险评分lasso_coxscore＝∑(基因表达值*回归系数)。

在训练集中采用Lasso Cox回归模型方法构建肿瘤复发预测lasso_cox模型后，采用时间依赖的ROC曲线(受试者工作曲线)的AUC(曲线下面积)评估模型的预测效果。AUC的取值范围为0到1之间，AUC越大说明模型的预测效果越好。lasso_cox回归模型的AUC_1年＝0.825、AUC_3年＝0.821、AUC_5年＝0.806，说明该模型在训练集中预测患者无复发生存期的效果较好，见图5。

(5)在测试集中评估模型预测效果

本实施例获得了测试集中所有样本的基因表达矩阵，并提取Lasso Cox模型所纳入的24个基因表达值(PAOX、SIGLEC7、PHAX、XCR1、TM4SF4、TRIOBP、MCMBP、HCFC1R1、ADNP2、NUP50、GTF2A2、BCCIP、FLJ90680、NVL、ESM1、GABRR2、FAM166A、USP14、JUNB、UBAP2、AP5B1、FAM46C、LDB3、JUP)。根据这24个基因的表达值与lasso_cox模型中24个基因的回归系数，即可获得模型在测试集中的肿瘤复发风险评分。根据II期结直肠癌患者的肿瘤复发风险评分的中位数(在本实施例中中位数为-2.748)，将患者分为高风险组(复发风险评分高于-2.748)与低风险组(复发风险评分低于-2.748)。分别绘制两组患者的Kaplan-Meier生存曲线，采用Logrank检验两组间的无复发生存期是否具有显著差异，并采用Mantel-Haenszel方法计算HR值(风险比)及其95％CI(置信区间)。测试集按照上述方法分成高复发风险组与低复发风险组后，两组间Logrank检验的p值＝0.0068，高复发风险组相比于低复发风险组的HR＝2.052，95％CI＝1.219-3.455，即模型能够显著的区分高复发风险和低复发风险的II期结直肠癌患者，见图6。

以上所述了本发明的一个实施例，本领域的普通技术人员可以理解，应当指出，在不脱离本发明的原理和宗旨的情况下可以对这些实施例、方法学、选取的模型、编程软件进行多种变化、修改、替换和补充，这些变化、修改、替换和补充也应该视为本发明的保护范围。

Claims

1.一种用于II期结直肠癌复发预测的标记基因，其特征在于所述标记基因包括：PAOX、SIGLEC7、PHAX、XCR1、TM4SF4、TRIOBP、MCMBP、HCFC1R1、ADNP2、NUP50、GTF2A2、BCCIP、FLJ90680、NVL、ESM1、GABRR2、FAM166A、USP14、JUNB、UBAP2、AP5B1、FAM46C、LDB3和JUP。

2.一种权利要求1所述标记基因在构建II期结直肠癌复发预测模型中的应用。

3.如权利要求2所述的应用，其特征在于所述模型的数学表达式：lasso_coxscore＝∑(基因表达值*回归系数)，所述基因表达值为标记基因mRNA表达值。

4.如权利要求2所述的应用，其特征在于所述回归系数见表1所示：

表1.lasso_cox回归模型中24个基因的回归系数

5.如权利要求2所述的应用，其特征在于模型的构建方法：(1)获取基因表达数据集：获取II期结直肠癌肿瘤样本的mRNA表达数据，并通过聚类分析检测离群值并排除离群值；(2)鉴定与II期结直肠癌患者复发相关的公共差异表达基因：采用基因芯片荟萃分析技术，以Logrank检验计算步骤(1)每一个数据集中每一个基因的p值；然后，采用minP方法合并每一个数据集中每一个基因的p值，最终得到每一个基因的p_minP值；采用Benjamini-Hochberg方法对每一个基因的p_minP值进行校正，得到每一个基因的FDR值；根据FDR<0.1的标准筛选与II期结直肠癌患者复发相关的公共差异表达基因；(3)公共差异表达基因的信号通路富集分析：采用Metascape数据库对步骤(2)鉴定得到的与II期结直肠癌患者复发相关的公共差异表达基因进行信号通路富集分析，并以p＝0.01作为阈值筛选差异表达基因显著富集的信号通路，即p值小于0.01的信号通路；(4)模型构建：根据步骤(2)II期结直肠癌患者复发相关的公共差异表达基因的表达值，采用Lasso Cox方法构建II期结直肠癌患者肿瘤复发预测lasso_cox模型；

所述模型的数学表达式：lasso_coxscore＝∑(基因表达值*回归系数)。