[发明专利]一种新的脑胶质瘤分子分型方法

权利要求书

1.一种新的脑胶质瘤分子分型方法，其特征在于：包括如下步骤：

步骤一：收集数据

从数据库中获取基因表达谱和临床信息数据，五个胶质瘤队列纳入本研究：TCGA-GBMLGG队列(n＝892)，三个CGGA队列(mRNA-array(n＝301)，mRNAseq_325(n＝325)和mRNAseq_693(n＝693))和Rembrandt队列(n＝475)，从公开数据库收集了四个独立的免疫治疗队列，包括：①Roh队列：抗CTLA-4、抗PD-1治疗队列；②GSE100797：过继性T细胞治疗队列；③GSE78220：抗PD-1治疗队列；④GSE35640：抗MAGE-A3治疗队列，根据recistv1.1标准，完全缓解(CR)或部分缓解(PR)患者被视为免疫治疗有反应者，疾病稳定(SD)或疾病进展(PD)患者被视为免疫治疗无反应者，不可评估(NE)患者则被剔除；

步骤二：突变特征谱

突变特征库(第二版)聚焦于碱基置换突变，突变点的碱基置换包含六种类型：CA，CG，CT，TA，TC和TG，突变点两侧(5'和3'端)各可搭配四种碱基(A、T、C、G)，最终可产生96种可能的突变类型(6种突变位点的碱基替换类型×4种，5'碱基×4种，3'碱基)，在体细胞中，各种机制引起DNA损伤，进而发生体细胞突变，使得细胞基因组不断变化，各种类型的突变不断积累，最终，在96种突变类型上具有不同的积累，形成一个个独特的突变积累组合，检测到的每种组合即为一种“突变特征”；

步骤三：亚型识别

(一)、数据获取：从COSMIC网站获取获取每种突变特征的特征数据信息，从TCGA-LGG和TCGA-GBM获得的体细胞突变数据在去除沉默突变后用于构建每个个体的突变特征谱，参考基因组为h38

(二)、R包：DeconstructSigs和NMFpackage

(三)、方法

①移除突变数据中的沉默突变

②利用计算机技术将移除沉默突变的突变数据转化为突变环境矩阵

③使用DeconstructSigs包，分析每个样本中30个signature的组成比例，参考的signature是COSMIC，cutoff值设置为0.06，标准化方式使用“exome2genome”，最终我们得到了一个矩阵(行为30个signature，列为每个样本，cell的值为signature在每个样本的比例，所有signature加起来刚好为1)

④使用NMF包进行提取和聚类分析，设置潜在的ranks＝2:5，运算执行次数设置为50，method设置为‘lee’，最终通过从cophenetic系数和轮廓系数决定最佳rank＝4，如图1，即根据TCGA-GBMLGG队列中每个患者的突变特征谱，将胶质瘤分为四种分子亚型最优

⑤非负矩阵分解的一个特性是倾向于产生观察数据的稀疏表示，导致双聚类的自然应用，通过少量特征表征样本组，在NMF模型中，根据对每个样本贡献最大的基础组分(即在系数矩阵的每列中具有最大系数的基础组分)对样本进行分组，然后，通过根据基础矩阵计算的基础特异性评分选择的一组特征对每组样本进行表征，上述过程由NMF包实现，根据所有患者的突变特征谱，构建NMF模型，并通过extractFeatures函数(方法设置为“max”)对最基本的特定特征进行提取，最终，患者分为4个基础组，并提取出11种可表征每组样本的关键突变特征(mutational signature 1，3，5，8，12，13，15，16，21，26和30)，结果如图2所示，每种亚型都有特异的突变特征变量，而后，根据提取出的这11种最基础的突变特征，进行NMF聚类分析，将TCGA-GBMLGG队列的所有患者分为四种亚型，命名为C1，C2，C3和C4，如图3所示；

步骤四：胶质瘤风险指数(glioma risk index，GRI)的构建

①数据分析，TCGA-GBMLGG胶质瘤队列作为训练集进行建模；三个CGGA队列(mRNA-array(n＝301)，mRNAseq_325(n＝325)和mRNAseq_693(n＝693))

②筛选4种亚型之间的共同差异表达基因(differentially expressed genes，DEGs)：分别将每一种亚型组与其他三个亚型组配对比较，使用edgeR软件包进行基因表达差异分析，标准为校正p值0.05和|log2 FC|1，结果：鉴定出四组DEGs后取交集，共识别出708个DEGs

③对708个DEGs做单因素cox回归分析，|1-HR|0.5和P-adjust0.05的基因作为预后相关基因纳入下一步分析，(HR：危险比；P-adjust：校正P值)，结果：一共提取到226个基因纳入下一步分析

④将这226个基因两两组合，形成基因对(gene pair)，每对基因包含两个基因，A和B，表示为A|B，在一个样本中，若基因A的表达值高于基因B，那么该A|B基因对的值标记为1，反之为0，这样的赋值设计的优势在于，只需要关注两个基因mRNA表达之间的数学关系，完全忽略了不同平台之间的批次效应，不需要定义截断值(cut-off值)，增加了临床适用性，在TCGA队列中，通过上述赋值方法对每个样本中的所有基因对进行赋值，并剔除在80％以上样本中评分为全为0或全为1的基因对，最终得到一个由样本和基因对构成的二进制0/1矩阵，用于下一步骤的分析

⑤根据上述0/1矩阵，对其中包含的基因对进行Lasso回归，以降维并建模，最优模型由惩罚系数λ决定，当惩罚系数lambda＝0.07094148时，模型最优，最优模型包含了由36个基因组成的44个基因对基于这44个基因对

⑥GRI计算公式设计如下：

GRI＝∑β_i×GPV(i)

其中i为Lasso回归得到的关键基因对，GPV是i的赋值(0/1)，β是i对应的Lasso回归系数，最终，GRI计算公式为：GRI＝0.022×GPV(AGXT|BPIFB4)+0.002×GPV(AGXT|STMND1)+0.040×GPV(C5orf46|CSAG3)+0.031×GPV(CD70|FMO1)+0.070×GPV(DCSTAMP|FMO1)+-0.087×GPV(EDARADD|MAGED4)+0.051×GPV(EMP3|SOCS2)+0.234×GPV(EN1|FAT2)+0.076×GPV(EN1|PXDNL)+0.016×GPV(EN1|TDO2)+-0.020×GPV(ESR2|MAGED4)+0.138×GPV(FAM92B|FCAMR)+-0.285×GPV(FAT2|ITGBL1)+-0.160×GPV(FAT2|TMEM71)+-0.244×GPV(FAT2|WISP1)+-0.058×GPV(FMO1|HIST1H2BH)+-0.021×GPV(FMO1|HOXD11)+0.401×GPV(GPR1|KLRC1)+0.156×GPV(HCRT|PRSS48)+0.115×GPV(HCRT|SLC6A18)+0.055×GPV(HIST1H2AJ|SLC6A18)+0.019×GPV(HIST1H3C|RIPPLY3)+0.161×GPV(HIST1H3F|USP29)+0.270×GPV(HIST1H3G|PRSS48)+0.030×GPV(HIST1H3G|SLC6A18)+0.040×GPV(HIST1H3G|SLCO1B1)+0.098×GPV(HIST1H4B|SLCO1B1)+0.038×GPV(HIST1H4D|SLC6A18)+0.260×GPV(HIST1H4D|SPRR2A)+0.177×GPV(HOXA6|SLC6A18)+0.022×GPV(HOXA6|SPRR2A)+0.026×GPV(HOXD11|POTEF)+0.130×GPV(HOXD11|PRSS48)+0.038×GPV(HOXD11|TCF23)+0.007×GPV(HOXD11|UCN2)+0.207×GPV(IGFBP2|SLC29A1)+0.005×GPV(IL36B|SLCO1B1)+0.057×GPV(MAGED4|MOCOS)+-0.041×GPV(METTL1|PLA2G5)+0.093×GPV(NPIPA7|SLC6A18)+0.243×GPV(PAEP|SLCO1B1)+0.115×GPV(PLEKHN1|TNFSF11)+0.016×GPV(POTEI|POTEJ)+0.004×GPV(RBP1|SOCS2)

⑦用survminer包确定GRI最佳截断点0.8321341，根据此截断值，可将患者分为高、低GRI两组，Kaplan-Meier生存分析表明高GRI患者的预后生存情况比低GRI患者差。