肿瘤转录组多模式信息分析平台PipeOne及其构建方法

技术领域

本发明涉及肿瘤转录组多模式信息分析领域，特别涉及肿瘤转录组多模式信息分析平台PipeOne及其构建方法。

背景技术

转录组是特定细胞或者组织中产生的所有RNA转录本的集合，该术语有时也仅指所有的mRNA。转录组学是应用高通量的方法如微阵列或高通量RNA测序技术(RNA-seq)，来研究转录组的学科，了解转录组是揭示细胞和组织在发育或疾病条件下分子组成的关键

研究发现超过90％的人类基因会经历可变剪切,生成不同的 mRNA异构体，这大大增加了人类转录组和蛋白质组的复杂性。可变剪切控制失调可能导致疾病

基因融合是通常由染色体重排产生的嵌合基因，一些融合基因可以驱动癌症，是作为治疗靶标和诊断的生物标志物。单核苷酸多态性 (single-nucleotidepolymorphism，SNP)是基因组中特定位置的一个核苷酸替代。RNA-seq数据也可以用SNP和基因融合的检测

除了mRNA外，从RNA-seq中还可以检测到大量的非编码RNA，如长链非编码RNA(lncRNA)。相当一部分的lncRNA已被证明在不同的细胞活动中具有功能，并在各种癌症中失调。例如，TUG1可以与 Polycomb阻遏物复合物相互作用并调节神经胶质瘤细胞的自我更新

转座子是基因组中可移动的DNA，占人类基因组的近一半

细胞中的不同生物过程可以相互影响。例如，RNA编辑可影响选择性剪切和环状RNA的生物发生

RNA-seq简介

Illumina短读长(read)cDNA测序是应用最广泛的RNA测序技术(short-read RNA-seq)，在SRA数据库中公开的RNA-seq数据中的95％

相比于第二代的short-read RNA-seq，第三代RNA测序技术如 long-read RNA-seq和Direct RNA-seq可以测得更长(1-50kb)，可以直接完整的将转录本序列测出来，不过三代测序技术测序通量和准确率较低，在基因定量上不如第二代测序

如没有特殊声明，本研究所指的RNA-seq都是短读长RNA-seq。 RNA-seq的开发和应用早在十多年前就已经开始了，在转录组学领域得到广泛应用，并成为生命科学研究的一种常用手段，其中最普遍的应用是差异表达分析，即鉴定在两类不同样品间有表达水平差异的基因

流程管理系统Nextflow和Snakemake简介

生物信息学的分析经常需要组合使用多款独立的软件，例如，基因组注释套件Sanger Companion pipeline

流程管理系统Nextflow

Snakemake

转录组分析平台简介

从RNA-seq原始数据获取某种信息，往往需要配合多款软件进行，比如计算基因表达水平，就需要从质量控制，序列比对，再到计算表达量，以及最后生成完整的数据表格，这些需要考虑每个软件的运行环境，文件的输入输出，中间文件处理等问题，当样本量很大时，还要考虑并行计算，批处理的方法来加快任务。为此，研究者发展了开发了很多RNA-seq分析流程，比如最近三年新发布的流程就包括 RNACocktail

此外还有专门鉴定特定RNA的流程，比如鉴定lncRNA的 lncPipe

虽然已经出现不少处理RNA-seq的流程，但是他们往往只关注转录组的一两种模式信息，即计算基因表达水平以进行差异表达分析，或者是可变剪切分析，比如UTAP，BioJupies,RASflow。有的可以分析多种模式的信息，比如RNACocktail，VIPER，hppRNA，aRNApipe，但一般只进行原始数据处理，或者下游分析只是对每种信息模式分别进行分析，没有将多种模式信息进行整合分析。还有的在安装和运行上需要有较高的编程水平才能驾驭，比如RNACocktail和aRNApipe。因此亟需一款安装方便，运行简单，能够整合多种模式信息的RNA-seq 分析流程。

2018全球癌症统计报告显示，全球每年癌症新增病例1800万例，新增死亡960万人，我国成为癌症发病和死亡人数最高的国家

RNA测序(RNA-seq)已广泛用于功能基因组学的研究，可以通过不同的生物信息学分析手段从RNA-seq数据中获得各种信息，包括基因表达水平，可变剪接(AS)，选择性聚腺苷酸化(APA)，基因融合，RNA编辑和单核苷酸多态性(SNP)。超过90％的人类基因经历了可变剪切

除了mRNA外，从RNA-seq中还可以检测到大量的非编码RNA，如 circRNA

细胞中的不同生物过程可以相互影响。例如，RNA编辑可影响可变剪切和环状RNA的生物发生

发明内容

为解决上述现有技术存在的问题，本发明的目的在于提供一种肿瘤转录组多模式信息分析平台PipeOne及其构建方法。

为达到上述目的，本发明的技术方案为：

一种肿瘤转录组多模式信息分析平台PipeOne，该平台包含三个模块；

其中，模块一是RNA-seq的原始数据处理，应用不同的生物信息学软件到FASTQ文件，最终得到多个格式一致的RNA-seq多模式的信息表格；

模块二使用机器学习方法随机森林应用于RNA-seq多模式信息数据来区分癌和癌旁组织，通过有监督的学习，找出有对区分癌和癌旁有贡献的特征及特征权重；

模块三将RNA-seq多模式的信息进行非负矩阵分解，将分解得到的结果用于癌症亚型分析，分好的亚型再进行有监督的随机森林分析，以对特征进行权重排序。

进一步的，所述模块一具体构建流程为：

RNA-seq的FASTQ格式文件作为输入数据，流程的第一步先确定数据是否为处理过的干净数据，如果不是则用FASTP进行原始数据清洗和质控；FASTP清洗后的FASTQ文件经过8个不同的Nextflow脚本进行处理，包括：

1.转录组重构及新lncRNA的查找，并得到mRNA和lncRNA的表达水平。用HISAT(版本2.1.0)

[1]外显子数大于1,或者是单个外显子的转录本长度大于2000 nt；

[2]长度大于或等于200nt；

[3]至少2个样品中TPM>0.1；

[4]CPAT、CPPred和PLEK都认为该转录本没有编码潜力。

2.将所有lncRNAs(GENCODE lncRNA,LNCipedia lncRNA和新 lncRNA)，根据它们与蛋白质编码基因在基因组上的距离可以分为5 类，即反义(antisense),内含子区(intronic),顺势调控 (cis-regulatory),双向(divergent),基因间(intergenic)。Antisense lncRNAs在其反义链上会与基因的外显子至少有1nt重叠。cis-regulatorylncRNAs和编码基因处在同一DNA链，基因组上的距离范围在2000nt以内。divergentlncRNA与蛋白质编码基因的基因组距离范围和cis-regulatory lncRNA一样，但是在反义链上。其余的与任何编码基因不重叠lncRNA认为是intergenic lncRNA。

3.我们用三个环状RNA检测软件来检测circRNA。(如果提供的 RNA-seq是totalRNA-seq)。主要有四个步骤，(1)比对：RNA-seq 数据经过BWA

4.逆转座子组(retrotranscriptome)的表达水平计算。用 bowtie2

5.融合事件检测用Arriba (https://github.com/suhrig/arriba)。参与融合事件的基因标记1，反之为0，生成融合事件表格。

6.RNA编辑水平计算。编辑位点和编辑率用SPRINT

7.可变剪切事件的检测。Reads经过STAR

8.单核苷酸多态性(SNP)检测。在剪切事件得到的STAR 2-pass 比对bam文件用picard标记重复，在用GATK

最后用Python脚本汇总这8个脚本生成的数据，得到16(15，如果RNA-seq不是total RNA)个csv(逗号分隔符)格式的表格，用以下一步的计算。PipeOne模块一的八个脚本可以一次运行，也可以单独运行。

进一步的，所述模块二用于寻找区分癌和癌旁的特征组合，其构建流程为：

首先根据数据提取方差最高的前K个(top Kvariance，默认1000) 的特征，然后将样本数据分为训练集(training set)和测试集 (testing set)，用户可以自己确定训练集和测试集的数据；训练集用Python机器学习包Scikit-learn进行随机森林训练，留一法交叉验证(leave one out cross validation)进行模型评估；模型超参数选择随机森林参数，超参数的默认范围为(用户可以使用参数自行调整)：

n_estimators(决策树的数量):[3,5,7,10,20,30,50, 100]，

max_depth(决策树最大深度):[2,3,4,7,10],

min_samples_split(拆分内部节点所需的最少样本数):[2,3, 4,5,7],

class_weight(类权重):[None,"balanced", "balanced_subsample"]，根据留一法交叉验证的结果，选择最精确的模型，然后计算特征的权重(feature importance,orFeature weight)。

根据特征权重值，选择前K(10,20,50,100,200,all)个特征再进行重复训练，并用测试集数据进行验证。

和模块二一样，每个信息矩阵V取方差最大前K个(top K variance，默认1000)的特征，每个矩阵V在不同参数下进行非负矩阵分解(Non-negative Matrix Factoriztion,NMF)，分解成Wx H矩阵，H矩阵是不同V矩阵的共同分解矩阵。再用各H矩阵进行 Kmeans聚类,选择k＝3-8(默认，用户可以自行选择)。聚类的结果用silhouette值进行评估。用R包survival和survminer对分类结果进行预后分析，选择以生存差异最显著的分类结果为最优分类结果。根据分类结果，再进行有监督的随机森林训练，评估模型的方法采用留一法交叉验证。计算特征对分类模型的贡献值(Feature Importance,or Feature weight)。

相对于现有技术，本发明的有益效果为：

本发明肿瘤转录组多模式信息分析平台，是一站式且易于使用的 RNA-seq多模式信息分析平台，可从大量RNA-seq样本中提取和整合多模式信息，从而系统地识别出疾病的关键分子特征并对癌症进行亚型分析。PipeOne是使用Nextflow(工作流管理系统)和Docker(跨平台容器)构建的，并用已发布工具构建八个脚本和两个自制的分析组件打包在一起，从而使用户免于繁琐地安装各种工具。

附图说明

图1为PipeOne概览图；

图2为PipeOne模块一的技术路线图；

图3为lncRNA鉴定及lncRNA和mRNA的表达水平计算的流程图；

图4为环状RNA鉴定流程图；

图5为模块二的路线图；

图6为模块三的路线图；

图7为PipeOne对肾癌相关因子进行排序图；

图8为特征基因富集的前10种疾病图；

图9为基因UMOD和AQP2在肿瘤(TCGA,红色)和正常组织(GTEx, 黑色)中的表达水平图。图片由数据库GEPIA2生成；

图10为PipeOne确定了三种与临床相关的肾癌亚型图；

图11为PipeOne构建的自定义参考注释图。

具体实施方式

下面结合附图和具体实施方式对本发明技术方案做进一步详细描述：

如图1所示，PipeOne包含三个模块，模块一(Module 1)是 RNA-seq的原始数据处理，应用不同的生物信息学软件到FASTQ文件，最终得到多个格式一致的RNA-seq多模式(multi-modal)的信息表格。模块二(Module 2)使用机器学习方法随机森林应用于RNA-seq多模式(multi-modal)信息数据来区分癌和癌旁组织，通过有监督的学习，找出有对区分癌和癌旁有贡献的特征及特征权重。模块三 (Module 3)将RNA-seq多模式的信息进行非负矩阵分解 (Non-negative matrix factorization)，将分解得到的结果用于癌症亚型分析，分好的亚型再进行有监督的随机森林分析，以对特征进行权重排序。如图1所示。

PipeOne模块一的技术路线

RNA-seq的FASTQ格式文件作为输入数据，流程的第一步先确定数据是否为处理过的干净数据(用户提供参数--cleaned true)，如果不是则用FASTP

1.转录组重构及新lncRNA的查找，并得到mRNA和lncRNA的表达水平。用HISAT(版本2.1.0)

[1]外显子数大于1,或者是单个外显子的转录本长度大于2000 nt；

[2]长度大于或等于200nt；

[3]至少2个样品中TPM>0.1；

[4]CPAT、CPPred和PLEK都认为该转录本没有编码潜力。

本研究将所有lncRNAs(GENCODE lncRNA,LNCipedia lncRNA和新lncRNA)，根据它们与蛋白质编码基因在基因组上的距离可以分为 5类，即反义(antisense),内含子区(intronic),顺势调控 (cis-regulatory),双向(divergent),基因间(intergenic)。Antisense lncRNAs在其反义链上会与基因的外显子至少有1nt重叠。cis-regulatorylncRNAs和编码基因处在同一DNA链，基因组上的距离范围在2000nt以内。divergentlncRNA与蛋白质编码基因的基因组距离范围和cis-regulatory lncRNA一样，但是在反义链上。其余的与任何编码基因不重叠lncRNA认为是intergenic lncRNA。

2.我们用三个环状RNA检测软件来检测circRNA。(如果提供的 RNA-seq是totalRNA-seq)。主要有四个步骤，(1)比对：RNA-seq 数据经过BWA

3.选择性腺苷酸化事件(alternative polyadenylation,APA)。用QAPA

4.逆转座子组(retrotranscriptome)的表达水平计算。用 bowtie2

5.融合事件检测用Arriba (https://github.com/suhrig/arriba)。参与融合事件的基因标记1，反之为0，生成融合事件表格。

6.RNA编辑水平计算。编辑位点和编辑率用SPRINT

7.可变剪切事件的检测。Reads经过STAR

8.单核苷酸多态性(SNP)检测。在剪切事件得到的STAR 2-pass 比对bam文件用picard标记重复，在用GATK

最后用Python脚本汇总这8个脚本生成的数据，得到16(15，如果RNA-seq不是total RNA)个csv(逗号分隔符)格式的表格，用以下一步的计算。PipeOne模块一的八个脚本可以一次运行，也可以单独运行。

PipeOne模块二的技术路线

模块二用于寻找区分癌和癌旁的特征组合。首先根据数据提取方差最高的前K个(top K variance，默认1000)的特征，然后将样本数据分为训练集(training set)和测试集(testing set)，用户可以自己确定训练集和测试集的数据。训练集用Python机器学习包Scikit-learn

n_estimators(决策树的数量):[3,5,7,10,20,30,50, 100]，

max_depth(决策树最大深度):[2,3,4,7,10],

min_samples_split(拆分内部节点所需的最少样本数):[2,3, 4,5,7],

class_weight(类权重):[None,"balanced", "balanced_subsample"]，根据留一法交叉验证的结果，选择最精确的模型，然后计算特征的权重(feature importance,orFeature weight)。

根据特征权重值，选择前K(10,20,50,100,200,all)个特征再进行重复训练，并用测试集数据进行验证。

和模块二一样，每个信息矩阵V取方差最大前K个(top K variance，默认1000)的特征，每个矩阵V在不同参数下进行非负矩阵分解(Non-negative Matrix Factoriztion,NMF)，分解成Wx H矩阵，H矩阵是不同V矩阵的共同分解矩阵。再用各H矩阵进行 Kmeans聚类,选择k＝3-8(默认，用户可以自行选择)。聚类的结果用silhouette值进行评估。用R包survival和survminer对分类结果进行预后分析，选择以生存差异最显著的分类结果为最优分类结果。根据分类结果，再进行有监督的随机森林训练，评估模型的方法采用留一法交叉验证。计算特征对分类模型的贡献值(Feature Importance,or Feature weight)。

PipeOne使用的软件和数据库

PipeOne使用生物信息学流程管理框架Nextflow

表1 PipeOne使用的软件

表2 PipeOne使用的参考基因组和数据库

测试数据

本研究用TCGA(The Cancer Genome Atlas)数据库的肾乳头状细胞癌(Kidneyrenal papillary cell carcinoma，KIRP)的RNA-seq (polyA+)数据，总共有320个样本，其中32个是癌旁样本，288 个是癌症样本。

PipeOne概览

相比现有的可以处理多模式数据的RNA-seq流程，PipeOne在原始数据处理上功能更多，而在下游的分析中PipeOne可以将多种模式信息整合在一起用于癌症相关特征排序(模块二)和亚型分析(模块三)，这个两个功能在目前的流程里尚未发现。PipeOne可谓开创了利用RNA-seq产生的多模式信息来整合研究癌症的方法。PipeOne的模块一在原始数据处理上几乎涵盖了目前RNA-seq数据处理的所有模式如表2所示。虽然暂不具备如同RNACocktail

在安装上RNACocktail将运行环境整合在Docker容器内，或者手动安装各种软件。VIPER支持Conda

在运行上，RNACocktail要求用户有良好的编程能力；PipeOne 则使用优秀的流程管理系统Nextflow

表2 PipeOne和其他可以处理RNA-seq多模式数据的流程比较

PipeOne识别肾癌的关键基因

为了证明PipeOne在寻找疾病相关特征的有效性，本研究将 PipeOne应用到来自TCGA的肾乳头状细胞癌(Kidney renal papillary cell carcinoma，KIRP)的RNA-seq数据，来探究PipeOne 是否可以找到和肾病相关的基因。我们先通过PipeOne模块一对肾癌的RNA-seq原始数据进行处理，得到13个不同模式的信息表格，即 mRNA，lncRNA，逆转录组，多聚腺苷酸化，Fusion，RNA编辑，SNP 各一个信息表格和可变剪切事件的6表格(可变3’剪切位点，可变 5’剪切位点，外显子跳跃,内含子保留,多内含子跳跃,互斥外显子)。取每个模式信息数据方差最大的前1000个特征(互斥外显子事件只检测有708个特征)，总共12708个特征做后续分析。

将PipeOne模块二应用到KIRP数据，随机预留36个样本(20个肿瘤样本，12个癌旁样本)做testset，剩下做tranningset(252 肿瘤样本，20个癌旁样本)。训练和测试结果如表3所示，特异性 (specificity)和灵敏度(sensitivity)都很高。

如图7所示，排名前十的(K＝10)特征包括四个mRNA(CYFIP2， DCN，UMOD和SERPINA5)，三个lncRNA(AP000775.1，lnc-GAPVD1-3 和lnc-SFRP1-2)，一个AS事件(BTBD11)和一个逆转录转座子 (HERVFRD_6p24.2)，这支持了本研究的假设，即多种类型的特征导致了该疾病的发生，因此整合多模式信息有助于更好地解释癌症患者的致病因素。

(A)与KIRP相关的特征，包括mRNA和lncRNA表达水平，逆转录转座子表达水平以及可变剪接。列表示在PipeOne模块二中使用的不同 K值。颜色代表功能的重要性，当K＝all时，按功能重要性对行进行排序。(B)对于每个K值，PipeOne获得了不同类型和数量的特征。AS:可变剪切,APA:可变多聚腺苷酸化，Retro:逆转座子。Weight: 相对权重，将每个模型的特征权重值归一化处理以方便作图显示。

表3随机森林模型在KIRP数据上的表现

使用ToppGene

通过分析发现PipeOne模块二的确可以找到和肾癌相关的已知基因，此外可以推断未知的lncRNA和逆转座子在肾癌中很可能扮演重要角色，可以通过实验验证，进一步研究他们的功能作用。

PipeOne揭示了与临床相关的肾癌亚型

本研究进一步将PipeOne模块三应用于KIRP进行亚型分析。识别出三个重要的亚型(亚型-0/ST0，亚型-1/ST1和亚型-2/ST2； P＝6.3E-4，log-rank test)(图2-8，补充表S1)，其中ST0生存率最差，ST2生存率最佳。对亚型分型贡献最高的前30个特征包括 13个mRNA，5个lncRNA，7个AS事件和5个反转录转座子元件(图 9)。这些因素(按功能重要性排序)包括SEPT2，TUG1(lncRNA)， CTNNA1，ITGB1，MT-CO1，MT-CO2，ERVLB4(逆转座子)，NCAPD2(AS 事件)，STAT1，TNFRSF11B和DDX5。SEPT2编码与膜丝相关的septin 家族成员，并可能调节细胞极性和组织发育。Septins与多种肾脏疾病有关，但其具体功能尚待阐明

如图10所示，(A)三种具有明显不同生存结果的亚型。ST0(亚型_0) 的预后较差，ST2(亚型_2)的存活率最高。(B)各种类型的因素有助于(A)中的亚型分析结果，包括mRNASEPT2，lncRNA TUG1和逆转录转座子ERVLB4(来自两个不同的基因座)。(C)将PipeOne亚型(ST0，ST1，ST2)与已报道的分子亚型(C1，C2a，C2b，C2c-CIMP，没有亚型注释的127个样品表示为NA)进行比较

分析传统RNA-seq数据集可以获得各种信息，如本研究所示，至少有七种模式(如果是total RNA-seq数据则为八种)特征可用。但是，以往的研究大多只关注这些特征中的一两个模式，因此可能会错失对疾病潜在机制的理解。本研究提出了一个一站式肿瘤转录组分析平台，该平台整合了多种已发布的RNA-seq分析工具，可以对每个测序样品进行全面分析。此外，本研究还开发了两个新模块，分别用于寻找与疾病相关的特征和癌症亚型分型。本研究将PipeOne应用在肾脏癌中的分析中，证明了其在利用多模式信息分析癌症样本中的有效性，并将与肾脏癌的发展与进展相关的潜在基因进行了优先排序。

PipeOne模块一重新构建转录组注释以提供全面的lncRNA注释，即GENCODElncRNA和LNCipedia lncRNA和检测的新lncRNA(和 GENCODE lncRNA、LNCipedia lncRNA没有重复)。PipeOne模块一的四个组件可以利用这个自定义的注释，即lncRNA定量，circRNA注释，融合基因检测和剪切事件。对于KIRP样品，PipeOne模块一鉴定了4,320个新lncRNA基因(图2-10C)的7,337个新的lncRNA转录本(主要是在基因间的，反义的和双向的lncRNA，图2-10A-B)，其表达水平与GENCODE和LNCipedia报道的lncRNA相似。

PipeOne的不足之处在于它目前仅关注传统的RNA-seq数据。但是，在此基础上，可以很容易地扩展和完善它，纳入处理其他高通量数据的功能，例如，DNA测序(DNA-seq)，通过微阵列或测序进行的 DNA甲基化谱，以及通过质谱法进行的蛋白质组学。另一个是，本研究仅通过对预后分析结果评估了亚型分析的效果。将来提供更多选择，结合其他方式评估PipeOne预测的亚型，例如，使用患者药物反应结果。使没有预后信息的其他疾病的亚型分析成为可能。此外，可以把PipeOne应用于其他癌症类型的研究，例如肝癌和乳腺癌以及许多复杂的疾病。

如图11所示，(A)与GENCODE v32和LNCipedia v5.2相比，共有7,337个lncRNA转录本是新的。(B)新的lncRNA转录本主要是基因间的，反义的和双向的。(C)共有4320个lncRNA基因是新的。(D) 新的lncRNA基因的表达水平分布与来自GENCODE和LNCipedia的相似。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何不经过创造性劳动想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书所限定的保护范围为准。