CN113862351B - 体液样本中鉴定胞外rna生物标志物的试剂盒及方法 - Google Patents

体液样本中鉴定胞外rna生物标志物的试剂盒及方法 Download PDF

Info

Publication number
CN113862351B
CN113862351B CN202010618721.4A CN202010618721A CN113862351B CN 113862351 B CN113862351 B CN 113862351B CN 202010618721 A CN202010618721 A CN 202010618721A CN 113862351 B CN113862351 B CN 113862351B
Authority
CN
China
Prior art keywords
exrna
data
gene
mapping
matrix
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010618721.4A
Other languages
English (en)
Other versions
CN113862351A (zh
Inventor
鲁志
史斌斌
李洋
陈旭鹏
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Tsinghua University
Original Assignee
Tsinghua University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Tsinghua University filed Critical Tsinghua University
Priority to CN202010618721.4A priority Critical patent/CN113862351B/zh
Publication of CN113862351A publication Critical patent/CN113862351A/zh
Application granted granted Critical
Publication of CN113862351B publication Critical patent/CN113862351B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • C—CHEMISTRY; METALLURGY
    • C12—BIOCHEMISTRY; BEER; SPIRITS; WINE; VINEGAR; MICROBIOLOGY; ENZYMOLOGY; MUTATION OR GENETIC ENGINEERING
    • C12Q—MEASURING OR TESTING PROCESSES INVOLVING ENZYMES, NUCLEIC ACIDS OR MICROORGANISMS; COMPOSITIONS OR TEST PAPERS THEREFOR; PROCESSES OF PREPARING SUCH COMPOSITIONS; CONDITION-RESPONSIVE CONTROL IN MICROBIOLOGICAL OR ENZYMOLOGICAL PROCESSES
    • C12Q1/00—Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions
    • C12Q1/68—Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions involving nucleic acids
    • C12Q1/6876—Nucleic acid products used in the analysis of nucleic acids, e.g. primers or probes
    • C12Q1/6883—Nucleic acid products used in the analysis of nucleic acids, e.g. primers or probes for diseases caused by alterations of genetic material
    • C12Q1/6886—Nucleic acid products used in the analysis of nucleic acids, e.g. primers or probes for diseases caused by alterations of genetic material for cancer
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/30—Detection of binding sites or motifs
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B50/00—ICT programming tools or database systems specially adapted for bioinformatics
    • C—CHEMISTRY; METALLURGY
    • C12—BIOCHEMISTRY; BEER; SPIRITS; WINE; VINEGAR; MICROBIOLOGY; ENZYMOLOGY; MUTATION OR GENETIC ENGINEERING
    • C12Q—MEASURING OR TESTING PROCESSES INVOLVING ENZYMES, NUCLEIC ACIDS OR MICROORGANISMS; COMPOSITIONS OR TEST PAPERS THEREFOR; PROCESSES OF PREPARING SUCH COMPOSITIONS; CONDITION-RESPONSIVE CONTROL IN MICROBIOLOGICAL OR ENZYMOLOGICAL PROCESSES
    • C12Q2600/00—Oligonucleotides characterized by their use
    • C12Q2600/158—Expression markers
    • C—CHEMISTRY; METALLURGY
    • C12—BIOCHEMISTRY; BEER; SPIRITS; WINE; VINEGAR; MICROBIOLOGY; ENZYMOLOGY; MUTATION OR GENETIC ENGINEERING
    • C12Q—MEASURING OR TESTING PROCESSES INVOLVING ENZYMES, NUCLEIC ACIDS OR MICROORGANISMS; COMPOSITIONS OR TEST PAPERS THEREFOR; PROCESSES OF PREPARING SUCH COMPOSITIONS; CONDITION-RESPONSIVE CONTROL IN MICROBIOLOGICAL OR ENZYMOLOGICAL PROCESSES
    • C12Q2600/00—Oligonucleotides characterized by their use
    • C12Q2600/178—Oligonucleotides characterized by their use miRNA, siRNA or ncRNA

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Chemical & Material Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • General Health & Medical Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Organic Chemistry (AREA)
  • Biophysics (AREA)
  • Genetics & Genomics (AREA)
  • Biotechnology (AREA)
  • Pathology (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Medical Informatics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Evolutionary Biology (AREA)
  • Zoology (AREA)
  • Wood Science & Technology (AREA)
  • Immunology (AREA)
  • Molecular Biology (AREA)
  • Bioethics (AREA)
  • Microbiology (AREA)
  • Oncology (AREA)
  • Hospice & Palliative Care (AREA)
  • Databases & Information Systems (AREA)
  • Biochemistry (AREA)
  • General Engineering & Computer Science (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本发明公开了体液样本中鉴定胞外RNA生物标志物的试剂盒及方法。本发明提供了获得癌症患者体液exRNA生物标志物的方法，本发明研究出针对体液exRNA特点的预处理、降噪、归一化方法，同时利用机器学习为主的最新人工智能技术，对鉴定可靠的exRNA生物标志物起到至关重要的作用。本发明的方法包括了表达值定量、归一化、特征选择、构建机器学习模型，可以更方便地让用户从测序数据出发得到生物标志物。

Description

体液样本中鉴定胞外RNA生物标志物的试剂盒及方法
技术领域
本发明属于生命科学领域，尤其涉及体液样本中鉴定胞外RNA生物标志物的试剂盒及方法。
背景技术
早期诊断癌症患者的五年生存率要比晚期诊断癌症患者高5-10倍，因此尽早的发现癌症可以提供最大的治愈机会(Aravanis et al.,2017)。常见的可以用作无创检测的生物标志物大多是蛋白质分子或者特定序列的DNA(如ctDNA)，这些标志物因为自身特性在灵敏度、误诊率和组织溯源上有着难以克服的困难，检测的经济成本也比较大，难以普及。例如2018年最新发表在Science上的CancerSEEK，通过整合61个DNA扩增片段和8种蛋白质，能够从血液中准确地检测出8种常见类型的癌症(Cohen et al.,2018)，但是该方法的检测成本需要500美元，在确定癌症类型上的表现也比较差(最低的只有-40％的溯源准确度)。而RNA由于其在中心法则中所处的特殊地位，越来越多的研究发现其在疾病发生发展中可以作为一种更有优势的标志物：RNA标志物与DNA和蛋白标志物相比，具有更好的敏感性和组织特异性(Lopez et al.,2015)(表1)。
表1.三种类型的生物标志物的特征比较
Figure BDA0002562158410000011
在多种如血清、唾液以及尿液的体液中等，可以检测到一类非侵入性细胞外RNA(extracellular RNA,exRNA)(Redzic et al.,2014)。exRNA的类型主要包括mRNA片段和多种非编码RNA：如miRNA、piRNA、snRNA、lncRNA、circular RNA等(图1)，他们具有不同的剪切形式，修饰种类和细胞内外分布等。因为这些分子由于具备自身稳定结构(Chen,2016；Memczak et al.,2013)、类细胞膜结构和蛋白质的保护，加上某些RNA具有特定的结构，exRNA在多种体液(血清、唾液、尿液等)中可以抵抗体液中RNase的降解，从而稳定存在(Paula Maria Godoy,2018)。这些exRNA可以成为一类有效的生物标志物，服务于人体健康状况检测和疾病的诊断，如癌症的早期诊断、肿瘤生长状况监测、以及预后辅助诊断(Xi etal.,2017)。
国际上已有一些规模较大的研究团队和商业组织，开始将exRNA(extracellularRNA)作为生物标志物进行探索性的研究。最近，美国NIH下属的转化科学国家中心(NCATS)启动了exRNA研究项目ERCC(Extracellular RNA Communication Consortium)(Ainszteinet al.,2015)，下设多个研究方向，在包括1)exRNA治疗方案；2)分子标志物；3)作用机制等多方面同时资助多个课题组开展研究。在之前的体液exRNA研究中，miRNA受到广泛关注。miRNA可在多种细胞中进行内源性表达，并分泌至多种体液中(血液、唾液和尿液)稳定存在。基于这些特征，miRNA可以作为非侵入性生物标志物，成为包括癌症在内的人类疾病的理想候选生物标志物之一(Nakamura et al.,2016)。各种新exRNA的发现及研究进展迅速，可以作为生物标志物的exRNA远不止miRNA这一类。同时，exRNA通常会被包裹进入外泌体(exosomes)、微囊泡(microvesicles,MVs)(Valadi et al.,2007)等囊泡结构中、以及非囊泡结构的核糖核酸蛋白复合体(RNPs)(Vickers et al.,2011)等不同组分中，有着不同的粒径大小。针对exRNA不同种类、不同组分分布等特点的鉴定和细致研究却相对稀少。此外，现有鉴定潜在exRNA生物标志物的生物信息学方法也局限在比较表达值的差异，而忽略了其他突变、调控特征。
RNA生物信息学的发展和挑战：为了解析体液exRNA测序得到的高通量、大噪音、片段化的数据，发展针对exRNA的生物信息分析技术是关键。体液中exRNA的测序数据具有非均一化、动态性强、易降解、呈现数据碎片化、杂音大、呈现批次效应等特点(图2A、图2B)，尚缺乏专业和标准化的生物信息学分析方法(Byron et al.,2016)。例如，不同批次提取的体液样本之间存在很大差异(图2C)(Yuan et al.,2016)。与此同时，不同个体间基因表达的差异也不容忽视，不同临床患者的基因表达存在很大的个体差异(图2C)(Wang et al.,2015)，如果不对数据进行前处理与归一化，将为差异表达的分析带来了很大的困难，对后续分析的结果带来很大的误差，虽然可以通过排秩(rank-based)的方法去除了个体差异的影响(Li et al.,2017)，但是仍有很大提升空间。除此以外，RNA分子除了能够反映基因组变异的信息，同时后转录调控过程使得RNA分子具有广泛的多态性；血液样本中的RNA分子不同于普通细胞或组织样本的存在形式，往往受到降解作用的影响，以碎片的形式存在；对多重标志物信息的整合模型缺乏研究。
现有的方法存在如下缺点1、大部分已有的exRNA数据分析流程没有考虑exRNA样本间的异质性和批次效应等问题，导致分析结果受到干扰。2、exRNA测序数据存在不均匀的片段化现象，按照一般RNA测序的表达值计算方法，容易忽略一条RNA中不同片段之间表达值的差异。3、现有的exRNA数据分析流程大多专注于解决RNA表达值定量问题。
由于针对解决上述问题生物信息学工具的缺乏，研究者往往使用非针对性的工具进行数据分析，这会相应的引入大量系统误差、导致鉴定不到可靠的exRNA生物标志物，为后续的验证与功能机制的研究带来困难。
发明内容
本发明一个目的是提供A和B所示的物质的用途。
本发明提供了A和B所示的物质在制备获得癌症患者体液exRNA生物标志物产品中的应用：
A)用于获得癌症病人和正常人体液不同样本的exRNA测序数据的物质；
B)记载对所述不同样本的exRNA测序数据进行数据分析的可读载体；
所述对不同样本的exRNA测序数据进行数据分析的方法如下：
1)对所述不同样本的exRNA测序数据进行预处理，得到基因计数矩阵；
所述预处理为如下1)-A或1)-B；
所述1)-A：若exRNA小于或等于50bp，则所述预处理在所述顺序映射中，依次包括去接头、顺序映射、exRNA结构域提取和基因计数矩阵的构建；
所述顺序映射包括如下：
1)-A-1)，将去接头后数据映射到人类rRNA数据库；
1)-A-2)，将1)-A-1)处理后未映射的数据继续映射到UniVec数据库；
1)-A-3)，将1)-A-2)处理后未映射的数据依次按照如下顺序映射到关注的RNA类别上：lncRNA、miRNA、mRNA、piRNA、snoRNA、snRNA、srpRNA、tRNA、TUCP RNA和Y RNA；
1)-A-4)，将1)-A-3)处理后未映射的数据映射到人类完整的基因组序列上；
1)-A-5)，将1)-A-4)处理后未映射的数据映射到circBase数据库中定义的环状RNA上；
所述exRNA结构域提取包括如下：将所述1)-A-4)映射后数据的每个转录本划分为20个nt的bin，然后计算每个bin的平均读数覆盖率，并过滤掉平均读数覆盖率低于5的bin；从第一个有效bin的起始位置开始进行双向搜索以获得局部最大的读数覆盖率；一旦找到局部最大值，则从该局部最大值的位置向5’端和3’端两个方向进行峰拓展，直到该位置的读数覆盖率小于5或小于局部最大值的1/2；删除短于10nt的峰，并合并相邻的峰；再针对每个样品独立执行了上述峰提取程序，并计算每个峰在所有样品中重复出现的次数，峰重复率大于10％的区域被保留并合并为可信峰，这些可信峰又称为结构域特征，即为exRNA结构域；
所述基因计数矩阵的构建为将所述exRNA结构域、1)-A-3)映射到miRNA的数据和1)-A-5映射后数据生成癌症病人基因计数矩阵和正常人基因计数矩阵；
所述1)-B：若exRNA大于50bp，则所述预处理依次包括去接头、顺序映射和基因计数矩阵的构建；
所述顺序映射包括如下：先将所述去接头数据映射到到人类rRNA数据库；再将前一步未映射的数据映射到人类hg38基因组；再将前一步未映射的数据映射到circBase数据库中定义的环状RNA；
所述基因计数矩阵的构建为将所述1)-B中的顺序映射中的映射到人类hg38基因组和所述映射到circBase数据库的数据生成癌症病人基因计数矩阵和正常人基因计数矩阵；
2)将1)得到的各个样本的癌症病人基因计数矩阵和正常人基因计数矩阵分别使用如下5种样本库大小归一化方法：CPM、TMM、UQ、RLE和CPM-top进行处理，得到各个样本的5种归一化处理后矩阵；再将所述各个样本的5种归一化处理后矩阵分别用limma、ComBat、和RUV方法进行批次效应校正，得到各个样本15种批次校正后矩阵；检测所述各个样本15种批次校正后矩阵，再从中选择mkNN和UCA两个指标均最接近于1的归一化处理和批次效应校正组合方法作为目标归一化方法和目标批次效应校正方法；最后采用该目标归一化方法和目标批次效应校正方法对1)得到的各个样本的癌症病人基因计数矩阵和正常人基因计数矩阵进行处理，得到各个样本的归一化和批次效应去除后的exRNA基因表达矩阵；
3)将所述各个样本的归一化和批次效应去除后的exRNA基因表达矩阵进行对数转换以让分布更加集中而接近正态分布；再对每个特征或每一行进行标准化，得到各个标准化处理得到的exRNA基因表达矩阵；再用如下7种特征选择方法对所述各个标准化处理得到的exRNA基因表达矩阵进行特征选择，得到各个癌症病人和正常人的7种exRNA基因表达矩阵特征选择结果；再将所述各个癌症病人和正常人的7种exRNA基因表达矩阵特征选择结果分别用如下逻辑回归、随机森林和决策树规则进行再次特征选择，得到癌症病人21种exRNA基因表达矩阵特征选择结果和正常人的21种exRNA基因表达矩阵特征选择结果；最后从所述癌症病人21种exRNA基因表达矩阵特征选择结果和正常人的21种exRNA基因表达矩阵特征选择结果中选出特征稳定性评估KI值以及四智能预测模型结果评估AUROC值均接近1的作为癌症病人目标exRNA候选生物标志物；
所述7种特征选择方法为(1)差异表达，(2)基于逻辑斯蒂回归的wrapper方法，(3)基于随机森林的wrapper方法，(4)ReliefF，(5)SURF，(6)MultiSURF和(7)SIS。
上述mkNN指标的公式如下：
Figure BDA0002562158410000051
其中B是批次的种类数，b代表批次的编号，Nb是批次编号为b的样本数，k代表选取的最近邻的样本数，N代表所有样本数；而
Figure BDA0002562158410000052
是每个样本的k个最近邻样本中与该样本来自同一批次的样本数的平均值。
上述中若exRNA测序数据获得中的测序文库构建引入内参spike-in RNA序列，则先映射到spike-in上。
上述A和B所示的物质在制备鉴定癌症患者体液样本中胞外RNA生物标志物也是本发明保护的范围。
上述A和B所示的物质在制备鉴定或辅助鉴定癌症患者产品中的应用也是本发明保护的范围。
上述A和B所示的物质在制备区分或辅助区分癌症患者和非癌症患者产品中的应用也是本发明保护的范围。
上述A和B所示的物质在制备预测待测样品是否来源于癌症患者产品中的应用也是本发明保护的范围。
本发明另一个目的是提供一种获得癌症患者体液exRNA生物标志物的系统。
本发明提供的系统，包括第一个目的中A和B所示的物质。
本发明还有一个目的是提供一种获得癌症患者体液exRNA生物标志物的方法。
本发明提供的方法，包括如下步骤：
1)获得癌症患者和非癌症患者体液样本的exRNA测序数据；
在本发明的实施例中，癌症病人和正常人(非癌症病人)血液中的exRNA测序数据(体液样本exRNA测序数据)可以从GEO(Gene Expression Omnibus)数据库中直接获得；也可以通过如下方法获得：1、分别提取癌症病人和正常人血液中的exRNA；2、测序：使用
Figure BDA0002562158410000053
Multiplex Small RNA Library Prep Set试剂盒(E7330S)上述1得到的exRNA构建测序文库；再使用Illumina HiSeq X10测序平台对exRNA构建测序文库进行测序，得到癌症病人体液样本exRNA测序数据和正常人体液样本exRNA测序数据。
2)对所述不同样本的exRNA测序数据进行预处理，得到基因计数矩阵；
所述预处理为如下2)-A或2)-B；
所述2)-A：若exRNA小于等于50bp，则所述预处理在所述顺序映射中，依次包括去接头、顺序映射、exRNA结构域提取和基因计数矩阵的构建；
所述顺序映射包括如下：
2)-A-1)，将去接头后数据映射到人类rRNA数据库；
2)-A-2)，将2)-A-1)处理后未映射的数据继续映射到UniVec数据库；
2)-A-3)，将2)-A-2)处理后未映射的数据依次按照如下顺序映射到关注的RNA类别上：lncRNA、miRNA、mRNA、piRNA、snoRNA、snRNA、srpRNA、tRNA、TUCP RNA和Y RNA；
2)-A-4)，将2)-A-3)处理后未映射的数据映射到人类完整的基因组序列上；
2)-A-5)，将2)-A-4)处理后未映射的数据映射到circBase数据库中定义的环状RNA上；
所述exRNA结构域提取(具体代码见实施例)包括如下：将所述2)-A-4)映射后数据的每个转录本划分为20个nt的bin，然后计算每个bin的平均读数覆盖率，并过滤掉平均读数覆盖率低于5的bin；从第一个有效bin的起始位置开始进行双向搜索以获得局部最大的读数覆盖率；一旦找到局部最大值，则从该局部最大值的位置向5’端和3’端两个方向进行峰拓展，直到该位置的读数覆盖率小于5或小于局部最大值的1/2；删除短于10nt的峰，并合并相邻的峰；5、针对每个样品独立执行了上述峰提取程序，并计算每个峰在所有样品中重复出现的次数，峰重复率大于10％的区域被保留并合并为可信峰，这些可信峰又称为结构域特征，即为exRNA结构域；
所述基因计数矩阵的构建为将所述exRNA结构域、2)-A-3)映射到miRNA的数据和所述2)-A-5映射后数据生成癌症病人基因计数矩阵和正常人基因计数矩阵；
所述2)-B：若exRNA大于50bp，则所述预处理依次包括去接头、顺序映射和基因计数矩阵的构建；
所述顺序映射包括如下：先将去接头数据映射到到人类rRNA数据库；再将未映射的数据映射到人类hg38基因组；再将未映射的数据映射到circBase数据库中定义的环状RNA；
所述基因计数矩阵的构建为将所述2)-B中的顺序映射中的映射到人类hg38基因组和所述映射到circBase数据库的数据生成癌症病人基因计数矩阵和正常人基因计数矩阵；
3)将2)得到的各个样本的基因计数矩阵分别使用如下5种样本库大小归一化方法：CPM、TMM、UQ、RLE和CPM-top进行处理，得到各个样本的5种归一化处理后矩阵；再将所述各个样本的5种归一化处理后矩阵分别用limma、ComBat、和RUV方法进行批次效应校正，得到各个样本15种批次校正后矩阵；检测所述各个样本15种批次校正后矩阵，再从中选择mkNN和UCA两个指标均最接近于1的归一化处理和批次效应校正组合方法作为目标归一化方法和目标批次效应校正方法；最后采用该目标归一化方法和目标批次效应校正方法对2)得到的各个样本的基因计数矩阵进行处理，得到各个样本的归一化和批次效应去除后的exRNA基因表达矩阵；
4)将所述各个样本的归一化和批次效应去除后的exRNA基因表达矩阵进行对数转换以让分布更加集中而接近正态分布；再对每个特征或每一行进行标准化；再用如下7种特征选择方法对上述各个标准化处理得到的exRNA基因表达矩阵进行特征选择，得到各个癌症病人和正常人的7种exRNA基因表达矩阵特征选择结果(图2D)；再将所述各个癌症病人和正常人的7种exRNA基因表达矩阵特征选择结果分别用如下逻辑回归、随机森林和决策树规则进行再次特征选择，得到癌症病人21种exRNA基因表达矩阵特征选择结果和正常人的21种exRNA基因表达矩阵特征选择结果；最后从所述癌症病人21种exRNA基因表达矩阵特征选择结果和正常人的21种exRNA基因表达矩阵特征选择结果中选出特征稳定性评估(KI，Kuncheva index)以及四智能预测模型结果评估(AUROC)均接近1的作为癌症病人目标exRNA候选生物标志物(图2E)；
所述7种特征选择方法对包括(1)差异表达，(2)基于逻辑斯蒂回归的wrapper方法，(3)基于随机森林的wrapper方法，(4)ReliefF，(5)SURF，(6)MultiSURF，(7)SIS。
由上述方法制备的鉴定癌症患者体液样本中胞外RNA生物标志物也是本发明保护的范围。上述标志物在是制备鉴定或辅助鉴定癌症患者产品中的应用也是本发明保护的范围。上述标志物在是制备区分或辅助区分癌症患者和非癌症患者产品中的应用也是本发明保护的范围。上述标志物在是制备预测待测样品是否来源于癌症患者产品中的应用也是本发明保护的范围。
本发明还提供一种预测待测患者是否为癌症样本的试剂盒。
本发明提供的试剂盒，上述的胞外RNA生物标志物。
在本发明中，存在如下名词解释：
1、Biomarker：生物标志物
2、ctDNA(Circulating tumor DNA)：循环肿瘤DNA
3、exRNA(extracellular RNA)：胞外RNA
4、miRNA：微小RNA
5、piRNA：piwi-interacting RNA
6、snRNA：核小RNA
7、lncRNA：长非编码RNA
8、circular RNA：环化RNA
9、rRNA：核糖体RNA
10、RNase：核糖核酸酶
11、Domain Calling:结构域的检测
12、Rank-based：排秩
13、Feature selection：特征选择
14、ROC：接收者操作特征曲线(receiver operating characteristic curve)，或者叫ROC曲线
15、AUC(area under the ROC curve)：曲线下面积
本发明研究出针对体液exRNA特点的预处理、降噪、归一化方法，同时利用机器学习为主的最新人工智能技术，对鉴定可靠的exRNA生物标志物起到至关重要的作用。
在所述顺序映射中，若exRNA测序数据获得中的测序文库构建引入内参spike-inRNA序列，则先映射到spike-in上。
针对体液exRNA测序数据的非均一化、动态性强、易降解、呈现数据碎片化、杂音大、呈现批次效应等特点，本发明的发明人开发了exSEEK软件，搭建了一套统一的数据处理流程，以期在最大程度上去解决exRNA测序数据中的问题，并通过特征选择(featureselection)的方法挑选潜在与癌症诊断和预后相关的exRNA生物标志物。利用鉴定到这些潜在生物标志物，通过搭建人工智能预测模型来预测癌症的患病风险与预后效果。该数据处理流程(图3)，大致可以分为三个主要部分：1、数据的预处理，体液样本的测序数据将通过质检步骤，挑选出质量较好的样本，并过滤掉rRNA、接头序列等测序片段。随后，测得的片段，将被比对和注释到不同的RNA类型，特别发展了结构域的检测的方法来解决exRNA数据的碎片化问题，对于小RNA测序数据，使用长RNA的结构域与miRNA，piRNA的全长特征组合的矩阵作为表达矩阵。2、数据噪音和干扰因素的去除，通过过滤、多种归一化和去除批次效应的方法，去除掉数据中存在的噪音和干扰因素。3、生物标志物的鉴定和预测模型的建立，通过多种特征选择的方法，挑选与疾病预前、预后相关的RNA生物标志物，并建立基于机器学习的预测模型。本发明的方法包括了表达值定量、归一化、特征选择、构建机器学习模型，可以更方便地让用户从测序数据出发得到生物标志物。exSEEK软件还提供了大量可定制的参数，便于用户探索不同的参数组合，找到适合的参数设置。
附图说明
图1为extracellular RNA(exRNA)的生成过程与分类。
图2为exRNA测序数据存在批次化严重和异质性等问题；A exRNA来源；B exRNA数据存在的四个问题：碎片化，稀疏化，异质性和批次效应；C问题及其对应解决方案；D特征选择和机器学习框架；E挑选出的生物标志物和模型效果评估。
图3为exRNA测序数据的处理流程；exSEEK主要包括三个部分：原始测序reads经过质量控制，序列比对和基因计数，处理为基因计数矩阵；基因计数矩阵经过样本库大小归一化和批次效应去除，得到表达值矩阵；从表达值矩阵进行机器学习和特征选择，寻找生物标志物和构建预测模型。
图4为结构域的检测方法以及结果统计，图4a为以AC006548.3(lncRNA)为例展示结构域的检测的结果，细胞外游离RNA有清晰的峰，而组织RNA则无此特征；图4b为exRNA数据相比于组织RNA数据有明显的峰；图4c为结构域的长度集中在40个碱基左右；图4d为通过icSHAPE和图5a为RBP结合位点数据分析。
图5为结构域的类型统计和差异分析；图5a为RBP结合位点数据分析，可以看到找到的结构域具有统计学上显著的结构稳定性和与RBP结合的趋势(p-value<＝0.01)；图5b为差异表达RNA的比例；图5c为各种类型RNA数量和表达量的统计。
图6为数据的归一化、去除批次效应；(a)展示使用的两种度量指标UCA(非监督聚类正确性)和mKNN(m-K最近邻)所选择的归一化与去批次效应的方法组合，右上角的方法是希望挑选的；(b)展示了使用最好的归一化方法RLE后数据的归一化情况；(c)展示了进一步去除批次效应后主成分分析的结果；(d)展示了去除批次效应模型显著地消除了批次效应；(e)展示了数据的异质性得到了消除，并且RNA表达值的大小关系得以保留。
图7为特征选择方法评估挑选；展示使用的特征选择方法，通过测试多种不同的特征选择方法，使用AUC和KI两个指标进行分类效果和稳定性的评估，筛选最好的特征选择方法选择的特征作为生物标志物。
图8为exRNA生物标志物预测模型的ROC曲线；已知的其他生物标志物相比，exSEEK模型挑选出的生物标志物可以取得很好的预测结果。
图9为挑选出的RNA的聚类热图和差异表达结果；挑选出的RNA作为生物标志物可以明显地区分正常人和早期癌症患者；点的位置代表该RNA在模型中的权重，大小代表表达值的对数，条形图的长度代表fold change。
具体实施方式
下述实施例中所使用的实验方法如无特殊说明，均为常规方法。
下述实施例中所用的材料、试剂等，如无特殊说明，均可从商业途径得到。
实施例1、获得癌症样本的exRNA生物标志物的方法
I、癌症病人血液中的exRNA测序数据的获得
癌症病人和正常人(非癌症病人)血液中的exRNA测序数据(体液样本exRNA测序数据)可以从GEO(Gene Expression Omnibus)数据库中直接获得；
也可以通过如下方法获得：
1、分别提取癌症病人和正常人血液中的exRNA；
2、测序；
使用
Figure BDA0002562158410000101
Multiplex Small RNA Library Prep Set试剂盒(E7330S)上述1得到的exRNA构建测序文库；再使用Illumina HiSeq X10测序平台对exRNA构建测序文库进行测序，得到癌症病人体液样本exRNA测序数据和正常人体液样本exRNA测序数据。
上述建库时可以引入内参spike-in RNA序列。
II、通过测序数据获得癌症样本的exRNA生物标志物
一、exRNA测序数据的预处理
1、RNA的测序数据去接头和映射处理
体液样本exRNA测序数据将通过质检步骤，挑选出质量较好的样本，并过滤掉rRNA、接头序列等测序片段，随后，测得的片段，将被比对和注释到不同的RNA类型，得到exRNA测序数据reads映射的RNA类型的文件。
根据不同的建库方法获得的exRNA测序数据可以为小RNA(小于或者等于50bp)，也可以为长RNA(大于50bp)；根据RNA长度的不同，选择如下A或B的方法分别对上述I获得的癌症病人体液样本exRNA测序数据和正常人体液样本exRNA测序数据进行去接头和映射处理：
A、针对小RNA的测序数据处理方法
1.1、去接头
使用cutadapt软件对上述I得到的癌症病人体液样本exRNA测序数据和正常人体液样本exRNA测序数据读数(reads)进行去3’端的接头(adapter)，得到去除接头reads。
cutadapt软件设置如下参数：平均质量得分低于30或exRNA测序数据读数(reads)长度短于16nt。
1.2、reads顺序映射(mapping)
将上述1.1得到的去除接头reads用Bowtie2软件按照如下顺序进行映射，得到的reads映射到多种类型RNA的注释文件，目的是使测序的RNAreads比对和注释到不同的RNA类型。
Bowtie2软件的设置参数应选择“--sensitive”选项，并且不建议使用“--local”参数进行局部匹配，因为这会增加reads错误mapping的比例。
上述按照如下顺序进行映射依次包括如下：
1)将上述1.1得到的去除接头reads映射到spike-in上；
若建库时没有引入内参spike-in RNA序列，则直接进行如下2)；
2)将1)处理后未映射的reads继续mapping(映射)到人类rRNA数据库(来自NCBIRefSeq数据库)；
3)将2)处理后未映射的reads继续mapping到UniVec数据库；
采用2)和3)中2个数据库的原因是由于某些游离RNA-seq数据可能受到载体(vector)的污染。
4)将3)处理后未映射的reads依次按照如下顺序mapping(映射)到关注的RNA类别上：lncRNA、miRNA、mRNA、piRNA、snoRNA、snRNA、srpRNA、tRNA、TUCP RNA和Y RNA。
由于本发明对体液中lncRNA的表达谱及其潜在的功能非常感兴趣，所以将lncRNA(来自GENCODE v27和MiTranscriptome数据库)设定为最高优先级以此来增加lncRNA峰检测的灵敏度；miRNA的序列信息来自于miRBase数据库；piRNA序列信息来自于piRNABank数据库；其余上述RNA类别均来自GENCODE v27数据库。
根据GENCODE V27注释文件的每个记录中的“gene_type”属性，从中提取其他各类RNA的基因组位置，并将每个转录本的所有外显子的序列连接起来，这一做法的好处是可以对跨外显子的reads进行mapping。对于具有多个转录本的RNA，仅提取其序列最长的转录本的位置信息以避免mapping过程中reads归属的模糊性(ambiguity)问题。
5)将4)处理后未映射的reads mapping到人类完整的基因组序列上(UCSC hg38序列信息)，以获取来自启动子(promoter)、增强子(enhancer)和重复区域(repeats)的reads信息。
启动子和增强子的具体位置由ChromHMM软件从9个细胞系的数据中提取。ChromHMM软件处理后的结果可以从UCSC Genome Browser上下载并通过CrossMap软件转换为human hg38坐标。重复序列的具体位置从UCSC Genome Browser的“rmsk”表来获取。
6)将5)处理后未映射的reads映射到circBase数据库中定义的环状RNA(circRNA)；
无法比对到circRNA的reads被定义为无法被比对的读数(unmapped reads)，舍去，保留上述所有映射数据。
B、体液长RNA测序数据的处理
2.1去接头
使用cutadapt软件对上述I得到的exRNA测序数据读数(reads)去除3’端的adapter接头，得到去除接头reads。
cutadapt软件设置如下参数：平均质量得分低于30或exRNA测序数据读数(reads)长度短于16nt。
2.2reads顺序映射(mapping)
将上述2.1得到的去除接头reads用Bowtie2软件按照如下顺序进行映射,得到的reads映射到多种类型RNA的注释文件，目的是使测序的RNAreads比对和注释到不同的RNA类型。
Bowtie2软件的设置参数应选择“--sensitive”选项，并且不建议使用“--local”参数进行局部匹配，因为这会增加reads错误mapping的比例。
上述按照如下顺序进行映射依次包括如下：
1)将上述2.1得到的去除接头reads映射到spike-in上；
若建库时没有引入内参spike-in RNA序列，则直接进行如下2)；
2)将1)处理后未映射的reads继续mapping到人类rRNA数据库(来自NCBI RefSeq数据库)；
3)将2)处理后未映射的reads继续mapping到人类hg38基因组；
使用STAR软件进行mapping,关键参数为--outSAMtypeBAM Unsorted；--outReadsUnmappedFastx；--outSAMmultNmax 1；--seedPerWindowNmax50。
4)将3)处理后未映射的reads映射到circBase数据库中定义的环状RNA(circRNA)；
无法比对到circRNA的reads被定义为无法被比对的读数(unmapped reads)，舍去，保留上述所有映射数据。
circRNA的mapping与小RNA使用的策略相似，但需要额外定义配对的两条reads在基因组上的位置顺序，即要求从read-1的5’-末端读取到read2的3'-末端所得的序列与circRNA的反向拼接(back-splicing junction)序列相同。
在上述两步mapping步骤中，只有在配对的两条reads均得到mapping的情况下，这对reads才能被用于后续的分析。由于文库构建的过程中可能存在PCR过度扩增的问题，因此测序数据中可能会包含大量重复的reads。为此，使用picardMarkDuplicates去除重复的配对reads。
2、小RNA(小于50bp)的exRNA结构域提取(峰提取)
针对于小RNA的测序数据进行此步骤，若是长RNA的测序数据直接进入下面步骤3；
鉴于体液中的exRNA存在严重的碎片化及reads分布不均匀的问题，想通过一些生信方法来找到信噪比较高的可信的片段化区域，并将其定义为exRNA结构域(domain)。域检测所使用的方法类似于CLIP-seq或ChIP-seq数据集中的峰值检测(peak calling)方法。然而，exRNA数据的peak calling方法与CLIP-seq或ChIP-seq数据的方法有以下几个方面的不同:(1)ChIP-seq实验通常会利用对照组(control)DNA文库来标准化实验组DNA文库中reads的非均匀覆盖，然而，大多数exRNA数据集中都不包含这样的对照组；(2)CLIP-seq常用的peak calling工具，如PARalyzer和CIMS，都是利用交联引起的特征事件来确定单核苷酸分辨率下的RBP结合位点，然而，无法从exRNA数据集中获取RBP的结合信息，因此如何确定exRNA数据的domain边界是一个关键点的问题；(3)传统的peaking calling工具，例如Piranha，使用每个read的起始位置来计算峰值的覆盖度(coverage)，这可能导致该类软件找出的峰值位置偏离真正的峰值位置。
基于以上问题，发明人开发了一种专门针对exRNA数据的域检测算法，其总体设计思想与Piranha软件类似，但是本发明的方法可以做到更高灵的敏度以及找到更加准确的域位置。
exRNA结构域提取的具体方法过程如下：
为了从小的exRNA-seq数据中识别出长RNA的重复短片段，提出了一种基于局部最大(local maximum)读数覆盖率(reads coverage)的结构域提取(峰提取)的方法，并将该方法应用到1中得到的每个转录本。
1)首先，将上述1中1.2中的5)映射到人类完整的基因组序列后数据的每个转录本划分为20个nt的bin，然后计算每个bin的平均读数覆盖率，并过滤掉平均读数覆盖率低于5的bin；
2)、从第一个有效bin(覆盖率大于5)的起始位置开始进行双向搜索以获得局部最大的读数覆盖率；
3、一旦找到局部最大值，则从该局部最大值的位置向5’端和3’端两个方向进行峰拓展，直到该位置的读数覆盖率小于5或小于局部最大值的1/2；
4、删除短于10nt的峰，并合并相邻的峰；
5、针对每个样品独立执行了上述峰提取程序，并计算每个峰在所有样品中重复出现的次数，峰重复率大于10％的区域被保留并合并为可信峰，这些可信峰又称为结构域特征，即为exRNA结构域。
exRNA结构域提取采用的具体代码如下：
Figure BDA0002562158410000141
Figure BDA0002562158410000151
Figure BDA0002562158410000161
Figure BDA0002562158410000171
Figure BDA0002562158410000181
Figure BDA0002562158410000191
Figure BDA0002562158410000201
3、基因计数矩阵的构建
对于短RNA来说，将上述1中1.2的4)得到的miRNA映射数据、上述1中1.2的6)得到的映射数据和上述2得到的exRNA结构域，使用带有“-t exon-g gene_id-M-s”参数的featureCounts软件生成癌症病人基因计数矩阵和正常人基因计数矩阵。
对于长RNA来说，将上述1中2.2的3)和4)以及使用带有“-t exon-g gene_id-M-s”参数的featureCounts软件生成癌症病人基因计数矩阵和正常人基因计数矩阵。
基因计数矩阵中，行为转录本或结构域(特征)，每一列为样本。
二、数据的归一化和去除批次效应
不同于传统的RNA-seq测序数据，使用微量建库得到的exRNA测序数据往往具有大量噪音与碎片化的特点。与此同时，不同的病人样本、取样和保存时间、建库方法等因素都会对数据引入潜在的噪音。通过归一化和去除批次效应的模型可以部分消除这一影响。
1、样本库大小进行归一化
不同测序样本中每种基因的原始reads计数由多方面因素决定，除了RNA真正的浓度之外，还有RNA提取效率、RT-PCR对RNA片段的选择、PCA扩增倍数、样本稀释倍数、测序深度等影响。这些因素对每个基因reads数的影响可认为是等倍数的，但在不同样本之间存在明显的差异。这种影响所有基因的倍数又可被称为样本库大小(library size)。在对不同样本进行整合分析前，我们需要对样本库大小归一化。样本库大小归一化的目标是对所有基因的原始reads数乘以某个倍数后，让它们归一化后的样本库大小相同。
将上述一得到的各个癌症病人基因计数矩阵和正常人基因计数矩阵，分别使用了4种常见的样本库大小归一化方法：CPM,TMM,UQ和RLE，以及1种修改过的CPM(CPM-top)进行处理，得到各个样本的5种归一化处理后矩阵。
CPM的修改版本(CPM-top)：对于每个样本，去除最高表达的k个基因，然后用剩余基因的reads数的和作为样本库大小。最后，对所有基因都除以该样本库大小值，计算出归一化的表达值。
以上5个归一化方法中，UQ、TMM、RLE均用R语言的edgeR函数库(Robinson et al.,2010)中相应归一化函数实现，CPM和CPM-top用R语言实现。
2、批次效应的校正
尝试用两种RNA-seq和一种单细胞RNA-seq中常见的批次算法对上述1的5个归一化方法处理后矩阵进行批次效应校正。
将上述1得到的各个样本的5种归一化处理后矩阵分别用如下3种方法：limma),ComBat)和RUV)进行批次效应校正，得到各个样本15种批次校正后矩阵。
limma和ComBat两个方法均需要以已知的批次信息作为输入，通过线性回归的方法对每个基因的表达值进行校正，消除在不同批次之间的差异。limma的removeBatchEffect函数是一个较为简单的批次校正算法。limma的基本假设是，每个基因的表达值是固定因素和随机因素的加和效应。不同样本之间批次的差异和其他的干扰因素属于一类固定效应，也是需要消除的效应。另外，每个基因的表达值由于生物学扰动，实验技术中的噪声的影响也存在差异，可假设服从一定的概率分布(例如正态分布)。对RNA-seq测序reads数而言，一般先转换为对数，让分布变得相对集中，接近正态分布。limma以每个基因的表达值作为因变量，以批次作为协变量进行回归，然后用表达值减去固定效应值(也可认为是残差)，就得到校正过的基因表达值。
ComBat与limma的思想较为类似，也是基于线性回归，假设批次效应对表达值的影响是加和效应。不同的是，ComBat认为对每个基因批次校正的值与该基因的表达值之间存在关系。由于样本数较少，如果数据中出现极端值，那么线性回归系数将受到较大影响，产生过度校正，产生新的批次效应。ComBat为了解决过度校正问题，采取了经验贝叶斯方法(empirical Bayesian method)，假设批次校正值服从某种先验分布，而先验分布的参数本身又从数据中估计得来。经验贝叶斯方法的应用在小样本数据中更加抵抗数据中的噪声，因此可能是比limma更准确的方法。
前两种批次校正方法都需要提供批次信息作为方法的输入。如果批次信息未知或者批次信息与生物学类别存在较大重合，那么前面两种依赖批次信息的方法将不再适合。批次信息与生物学类别重合时，若进行批次校正，那么生物学类别造成的差异也随之消除，导致校正后的表达值丢失大部分有用信息。不依赖批次信息的代表方法是RUV(removeunwanted variation)。RUV采用的方法是，寻找对照基因，即在不同条件之间没有差异的基因，通过因子分析(factor analysis)的方法，估计未知的批次效应。RUV细分为三种方法：RUVg,RUVs和RUVr。其中，RUVg和RUVs都是在提供负对照样本(不受需要研究的生物学因素影响的样本)和负对照基因的情况下使用的。而RUVr则不需要已知负对照基因或者负对照样本，而是先通过差异表达回归去除生物学因素后，用残差进行因子分析批次效应。使用RUV时，需要预先指定位置批次效应因子的个数k。k的选择与批次效应的复杂度有关，可认为反映的是描述批次效应所需要的独立变量数。
实验中，对有批次信息的样本，分别用上述三种方法进行批次效应去除，而对于没有批次信息或者批次与生物学因素重合的样本，只采用RUV进行批次效应去除。
3、表达值矩阵前处理效果评估
对于归一化和批次效应校正(统称表达值矩阵前处理)，希望达到两个方面的结果：表达值在批次之间没有差异；表达值在生物学因素之间有差异(Lazar et al.,2013)。从这两个方面，分别采用下面描述的mkNN和UCA两个指标来评估。
将上述2得到各个样本15种批次校正后矩阵用如下mkNN和UCA两个指标来评估，选择mkNN和UCA两个指标均最接近于1的组合方法作为目标归一化方法和目标批次效应校正方法，采用该目标归一化方法和目标批次效应校正方法对上述一得到的各个癌症病人基因计数矩阵和正常人基因计数矩阵进行处理，得到各个癌症病人和正常人归一化和批次效应去除后的exRNA基因表达矩阵。
3.1、mkNN指标
基于一篇单细胞测序数据整合方法(Butler et al.,2018)中采用的alignmentscore提出mkNN指标。原始的alignment score假设单细胞测序样本来自多个批次，同时属于不同的细胞类型。那么选一个细胞，根据基因表达谱寻找与其最接近的k个细胞，然后统计其中与该细胞来自相同批次的比例。如果存在明显的批次效应，那么最邻近的k个细胞与该细胞来此同一个批次的比例将较大(接近1)。相反，如果没有明显的批次效应，k个最邻近细胞的批次的比例应该与所有细胞中来自该批次的比例相近，即相当于从总体中随机挑选k个细胞。最邻近的k个细胞所占比例可通过假设无批次效应情况下的期望进行归一化，得到一个大部分范围在0-1之间的值，称为alignment score：
Figure BDA0002562158410000231
其中N表示全部细胞的个数，而x则是每个细胞的k个最邻近的细胞中与该细胞来自同一批次的比例的平均值。Alignment score取值小于1，而可能小于0。完全没有批次情况下，alignment score接近0。
Alignment score的很大局限性是，只能应用于两个批次的情况。
将alignment score扩展到多个批次，提出了mKNN指标：
Figure BDA0002562158410000232
其中B是批次的种类数，b代表批次的编号，Nb是批次编号为b的样本数，k代表选取的最近邻样本的数量，N代表所有样本数；而
Figure BDA0002562158410000233
是每个样本的k个最近邻样本中与该样本来自同一批次的样本数的平均值。实际使用中，为了方便，用1–mKNN代替mKNN，使得批次效应越明显的数据集该指标越接近0。
3.2、UCA指标
采用了一种用于评估非监督学习的指标UCA(unsupervised clusteringaccuracy,参考文献Lopez,R.,Regier,J.,Cole,M.B.,Jordan,M.I.,and Yosef,N.(2018).Deep generative modeling for single-cell transcriptomics.Nature Methods 15,1053–1058.)(使用python中的sklearn.cluster.KMeans软件包进行计算)。
本质上，UCA反映的是聚类效果或者不同类别样本之间的区分度。采用UCA而不是之间用分类器的好处是，UCA需要的计算量更小，而且采用非监督学习降低了过拟合的风险。原始矩阵首先通过数据降维，然后用K-means算法进行聚类。聚类的结果相当于预测的类别标签，而期望评估非监督学习与真实聚类的类别之间的一致性。真实类别和聚类的类别之间的匹配可看成是二分图的最大匹配为题，可采用匈牙利算法，每个样本所属类别对应上，然后建立混淆矩阵，计算出准确度。准确度计算与机器学习的分类问题中的准确度计算方法相同。
三、特征选择
1)、首先对上述二得到的各个癌症病人和正常人归一化和批次效应去除后的exRNA基因表达矩阵进行对数转换以让分布更加集中而接近正态分布；
2)、然后对每个特征(每一行)进行标准化(即减去特征在样本之间的平均值，然后除以标准差)；
3)再使用如下7种特征选择方法对上述各个标准化处理得到的exRNA基因表达矩阵进行特征选择，得到各个癌症病人和正常人的7种exRNA基因表达矩阵特征选择结果；
7种特征选择方法为(1)差异表达(使用DESeq2软件默认参数得到)，(2)基于逻辑斯蒂回归(logistic regression)的wrapper方法(scikit-learn软件包)，(3)基于随机森林(RandomForest)的wrapper方法(scikit-learn软件包)，(4)ReliefF(scikit-rebate软件包)，(5)SURF(scikit-rebate软件包)，(6)MultiSURF(scikit-rebate软件包)和(7)SIS(R软件包)。
4)将上述各个癌症病人和正常人的7种exRNA基因表达矩阵特征选择结果分别用如下逻辑回归(logistic regression)、随机森林(random forest)和决策树规则进行再次特征选择，得到21种exRNA基因表达矩阵特征选择结果。
在特征选择的过程中，由于样本数量少和类别不均衡，没有采用K倍交叉验证，而是每次完全随机把所有样本划分为训练集和测试集，并且保证训练集和测试集中正负样本的比例与总体近似相同(使用scikit-learn软件包中的StratifiedShuffleSplit)。整体的交叉验证采用的是嵌套的策略：模型的评估是在所有样本中用50次交叉验证得到的(外部循环)，而模型的超参数优化只在训练集上进行5折交叉验证来优化(内部循环)。在内部循环中，有10％的样本被选作测试集而剩余90％为训练集。50次外部循环中测试集上的AUROC的平均值作为最终汇报的模型表现。
在外部循环中，只有训练集的样本被用作内部循环。在采用的稳健特征选择中，每次随机抽取90％的的样本作为训练集，用于拟合模型和特征选择。特征选择可采用逻辑回归(logistic regression)或随机森林(random forest)。模型拟合完成后，每个特征的权重或重要性可计算出来。例如在逻辑回归模型中，回归系数绝对值越大的特征越重要。在随机森林模型中，被越多决策树规则使用的特征越重要。特征按照重要性从大到小进行排序后，选择所需数量(k)的特征。重复多次样本抽取后，将每个特征被抽取的次数从大到小进行排序，抽取最大的k个特征作为该轮外部循环选择的特征。
模型超参数的优化在特征选择之前进行。超参数采用的策略是网格搜索的方法。例如对逻辑回归而言，有一个参数C需要进行优化，在10e-5到10e5之间搜索。对随机森林而言，有两个可优化的超参数：最大树深度(max_depth)和决策树个数(n_estimators)。在25到75之间选择决策树个数，而在3到5之间选择最大树深度。本部分内容使用python的sklearn库实现。
5)从21种exRNA基因表达矩阵特征选择结果中选出特征稳定性评估(KI，Kunchevaindex，参考文献Abeel,T.,Helleputte,T.,Peer,Y.V.de,Dupont,P.,and Saeys,Y.(2010).Robust biomarker identification for cancer diagnosis with ensemblefeature selection methods.Bioinformatics 26,392–398.)以及四智能预测模型结果评估(AUROC)均接近1的作为目标exRNA候选生物标志物。
用KI(Kuncheva index)评估特征选择在不同样本选取之间的稳定性。对于一个包含n个样本和N个特征的表达值矩阵，对其进行不放回采用，得到80％的样本(0.8n)。采用某种特征选择方法在该子样本集上，选择d个特征，将其表示为fi。不放回抽样进行50次，每次采取不同的随机初始值。那么在任意两次抽取的特征集fi,fj之间，可计算出KI值：
Figure BDA0002562158410000251
其中d2项反映了完全随机抽取相同数量特征的稳定性。如果进行M次不放回采样，特征选择的平均KI值(AKI)为：
Figure BDA0002562158410000252
如果M足够大，那么对于完全随机的特征选择，AKI接近于0，意味着特征选择非常不稳定。如果每次不放回采样后，总是选择固定的特征，那么AKI接近1，意味着特征选择非常稳定。
用AUROC衡量其模型的准确度和精确度(AUROC越接近于1说明模型作为exRNA癌症诊断标志物的性能越好)。
实施例2、肝癌样本中exRNA测序数据处理获得exRNA生物标志物的验证
I、癌症病人血液中的exRNA测序数据的获得
收集了GEO(Gene Expression Omnibus)中四项研究的数据：GSE123972，GSE113994，GSE53080和GSE94582(表2)；GSE123972数据集(短RNA)包含了来自肝癌(hepatocellular carcinoma,HCC)患者的30个血浆样本，来自健康供体(HD)的10个样本和来自慢性乙型肝炎感染者(chronic hepatitis B infection,CHB)的3个样本。
收集了GEO中研究的数据：GSE113994，GSE53080，GSE113994，GSE53080和GSE94582均来自健康供体，其中GSE113994包含了216个来自健康供体(HD)的样品。GSE113994包含的216个健康供体(HD)的样品中存在许多低质量的测序结果，因此过滤掉了测序深度少于2,000,000个读数的样本，或映射到NCBI UniVec数据库的读数大于总读数20％的样本，亦或映射到非人类基因组的读数大于总读数的20％的样本。
表2测试数据集的样本统计
Figure BDA0002562158410000261
为了更好地解释机器学习模型找到的体液中片段化RNA的生物学功能，以及揭示这些片段RNA的潜在产生机制，对这些RNA进行了二级结构和RBP相互作用的分析。为此，下载并处理了GEO数据库中的RNA二级结构高通量测序(icSHAPE)数据(GSE74353)，该icSHPAE数据来源于人HEK293T细胞系。并从POSATR2(http://lulab.life.tsinghua.edu.cn/postar2/)数据库中下载并处理了各类RBP结合位点的数据。
II、通过测序数据获得癌症样本的exRNA生物标志物
一、exRNA测序数据的预处理
1、RNA的测序数据去接头和映射处理
按照实施例1的II、一的1的A方法对进行处理；
2、小RNA(小于50bp)的exRNA结构域提取(峰提取):按照实施例1的II、一的2的方法处理；
使用结构域检测算法得到的结构域示例如图4a所示，可以发现总体水平上exRNA的结构域特征明显强于组织中RNA，如图4b所示。图4c显示结构域的长度集中在40个碱基左右；图4d通过icSHAPE和图5a为RBP结合位点数据分析，可以看到找到的结构域具有统计学上显著的结构稳定性和与RBP结合的趋势(p-value<＝0.01)，证明找到的结构域作为特征具有稳定存在的可能；图5c各种类型RNA数量和表达量的统计，可以发现除了miRNA，其他种类的RNA也有较多的种类和一定的比例；图5b差异表达RNA的比例，可以见到多种RNA均有差异表达，有作为生物标志物的潜力。
3、基因计数矩阵的构建:按照实施例1的II、一的3的方法处理；
二、数据的归一化和去除批次效应
1、样本库大小的进行归一化
2、批次效应的校正
3、表达值矩阵前处理效果评估
将上述2的5个归一化方法和3个批次效应校正方法处理的结果用如下mkNN和UCA两个指标来评估，选择mkNN和UCA两个指标均最接近于1的组合方法(以RLE作为归一化方法，以COMBAT作为批次效应校正方法)以作为目标归一化方法和目标批次效应校正方法，采用该目标归一化方法和目标批次效应校正方法对上述一步骤3得到的exRNA基因表达矩阵进行处理，得到归一化和批次效应去除后的exRNA基因表达矩阵。
数据的总体处理效果如图6c所示。使用上述实施例1二中1-3部分的指标UCA和mKNN作为评价聚类效果和去除批次效应效果的度量方式，来综合挑选一个归一化方法和去除批次效应方法的组合，如图6a所示。通过综合评估可以选择最佳的矩阵处理方法，在本实施例中为RLE和Limma的组合(mkNN和UCA两个指标分别为0.794，0.686)，如图6b所示，首先使用RLE可以对样本进行样本库大小的归一化，用RLE图分析了归一化之前和之后的结果。RLE分析中，在大部分基因在样本间没有差异表达的前提假设下，所有样本的RLE分布的中位数应该在0附近。归一化之前，不同样本RLE分布按照数据集来源呈现异质性，例如GSE123972相对于其他基因差异明显偏高，很可能是因为测序深度高于其他数据集的缘故。用RLE进行归一化后，大部分样本的RLE中位数都位于0附近，RLE的中位数距离0的偏差在样本之间的平均值从4.25降低到了0.45，可认为很大程度消除了数据异质性。然后如图6c所示可以看到使用Limma进行批次效应去除的效果，从PCA图可以观察到，去除批次效应前，正常样本按照数据集来源明显区分到不同的类别，而正常与癌症样本没有明显区分开。应用limma进行批次效应校正后，不同来源的正常样本较为均匀地混合到一起(mKNN从0.055上升到0.792)，同时正常样本与癌症样本的区分更加明显了。如图6d所示观察方差分解分析的结果，发现，批次效应去除之前，批次占有比类别更大的贡献，而去除批次效应后大大降低了。图6e展示了数据的异质性得到了消除，并且RNA表达值的大小关系得以保留。
三、特征选择
采用实施例1的II、三的方法，对上述二得到的归一化和批次效应去除后的exRNA基因表达矩阵进行特征选择，选出特征稳定性评估(KI，Kuncheva index)以及四智能预测模型结果评估(AUROC)均接近1的作为目标exRNA候选生物标志物。
使用多种机器学习模型组合来进行特征选择和分类，如图7所示，使用AUC和KI两个指标进行分类准确性以及特征选择稳定性的衡量指标。挑选出AUC和KI两个指标综合起来最好的模型作为特征选择和分类模型。在本例中使用DiffExp_TTest作为特征选择模型，用RandomForest作为分类模型。
使用ROC曲线评估挑选出的生物标志物和已知生物标志物的分类效果，并且获得灵敏性，特异性和精度等指标，可以得出比经典的生物标志物更好的分类效果，如图8所示。
将挑选出的生物标志物绘制热力图也可用于展示分类效果，同时可以获得挑选出的生物标志物对智能预测模型的贡献，表达值等信息，如图9所示。从图9可以看出，标志物miR-122-5p和miR-122-3p为目标exRNA候选生物标志物，这些标志物为现有文献报道的肝癌标志物，表明，本发明的方法正确。

Claims (7)

1.A和B所示的物质在制备获得癌症患者体液exRNA生物标志物产品中的应用：
A)用于获得癌症病人和正常人体液不同样本的exRNA测序数据的物质；
B)记载对所述不同样本的exRNA测序数据进行数据分析的可读载体；
所述对不同样本的exRNA测序数据进行数据分析的方法如下：
1)对所述不同样本的exRNA测序数据进行预处理，得到基因计数矩阵；
所述预处理为如下1)-A或1)-B；
所述1)-A：若exRNA小于或等于50bp，则所述预处理在顺序映射中，依次包括去接头、顺序映射、exRNA结构域提取和基因计数矩阵的构建；
所述顺序映射包括如下：
1)-A-1)，将去接头后数据映射到人类rRNA数据库；
1)-A-2)，将1)-A-1)处理后未映射的数据继续映射到UniVec数据库；
1)-A-3)，将1)-A-2)处理后未映射的数据依次按照如下顺序映射到关注的RNA类别上：lncRNA、miRNA、mRNA、piRNA、snoRNA、snRNA、srpRNA、tRNA、TUCP RNA和Y RNA；
1)-A-4)，将1)-A-3)处理后未映射的数据映射到人类完整的基因组序列上；
1)-A-5)，将1)-A-4)处理后未映射的数据映射到circBase数据库中定义的环状RNA上；
所述exRNA结构域提取包括如下：将所述1)-A-4)映射后数据的每个转录本划分为20个nt的bin，然后计算每个bin的平均读数覆盖率，并过滤掉平均读数覆盖率低于5的bin；从第一个有效bin的起始位置开始进行双向搜索以获得局部最大的读数覆盖率；一旦找到局部最大值，则从该局部最大值的位置向5’端和3’端两个方向进行峰拓展，直到该位置的读数覆盖率小于5或小于局部最大值的1/2；删除短于10nt的峰，并合并相邻的峰；再针对每个样品独立执行了上述峰提取程序，并计算每个峰在所有样品中重复出现的次数，峰重复率大于10％的区域被保留并合并为可信峰，这些可信峰又称为结构域特征，即为exRNA结构域；
所述基因计数矩阵的构建为将所述exRNA结构域、1)-A-3)映射到miRNA的数据和所述1)-A-5映射后数据生成癌症病人基因计数矩阵和正常人基因计数矩阵；
所述1)-B：若exRNA大于50bp，则所述预处理依次包括去接头、顺序映射和基因计数矩阵的构建；
所述顺序映射包括如下：先将所述去接头数据映射到到人类rRNA数据库；再将前一步未映射的数据映射到人类hg38基因组；再将前一步未映射的数据映射到circBase数据库中定义的环状RNA；
所述基因计数矩阵的构建为将所述1)-B中的顺序映射中的映射到人类hg38基因组和所述映射到circBase数据库的数据生成癌症病人基因计数矩阵和正常人基因计数矩阵；
2)将1)得到的各个样本的癌症病人基因计数矩阵和正常人基因计数矩阵分别使用如下5种样本库大小归一化方法：CPM、TMM、UQ、RLE和CPM-top进行处理，得到各个样本的5种归一化处理后矩阵；再将所述各个样本的5种归一化处理后矩阵分别用limma、ComBat、和RUV方法进行批次效应校正，得到各个样本15种批次校正后矩阵；检测所述各个样本15种批次校正后矩阵，再从中选择mkNN和UCA两个指标均最接近于1的归一化处理和批次效应校正组合方法作为目标归一化方法和目标批次效应校正方法；最后采用该目标归一化方法和目标批次效应校正方法对1)得到的各个样本的癌症病人基因计数矩阵和正常人基因计数矩阵进行处理，得到各个样本的归一化和批次效应去除后的exRNA基因表达矩阵；
3)将所述各个样本的归一化和批次效应去除后的exRNA基因表达矩阵进行对数转换以让分布更加集中而接近正态分布；再对每个特征或每一行进行标准化，得到各个标准化处理得到的exRNA基因表达矩阵；再用如下7种特征选择方法对所述各个标准化处理得到的exRNA基因表达矩阵进行特征选择，得到各个癌症病人和正常人的7种exRNA基因表达矩阵特征选择结果；再将所述各个癌症病人和正常人的7种exRNA基因表达矩阵特征选择结果分别用逻辑回归、随机森林和决策树规则进行再次特征选择，得到癌症病人21种exRNA基因表达矩阵特征选择结果和正常人的21种exRNA基因表达矩阵特征选择结果；最后从所述癌症病人21种exRNA基因表达矩阵特征选择结果和正常人的21种exRNA基因表达矩阵特征选择结果中选出特征稳定性评估KI值以及智能预测模型结果评估AUROC值均接近1的作为癌症病人目标exRNA候选生物标志物；
所述7种特征选择方法为(1)差异表达，(2)基于逻辑斯蒂回归的wrapper方法，(3)基于随机森林的wrapper方法，(4)ReliefF，(5)SURF，(6)MultiSURF和(7)SIS。
2.权利要求1中A和B所示的物质在制备鉴定癌症患者体液样本中胞外RNA生物标志物产品中的应用。
3.权利要求1中A和B所示的物质在制备鉴定或辅助鉴定癌症患者产品中的应用。
4.权利要求1中A和B所示的物质在制备区分或辅助区分癌症患者和非癌症患者产品中的应用。
5.权利要求1中A和B所示的物质在制备预测待测样品是否来源于癌症患者产品中的应用。
6.一种获得癌症患者体液exRNA生物标志物的系统，包括权利要求1中A和B所示的物质。
7.一种获得癌症患者体液exRNA生物标志物的方法，包括如下步骤：
1)获得癌症患者和非癌症患者体液样本的exRNA测序数据；
2)对不同样本的exRNA测序数据进行预处理，得到基因计数矩阵；
所述预处理为如下2)-A或2)-B；
所述2)-A：若exRNA小于等于50bp，则所述预处理在顺序映射中，依次包括去接头、顺序映射、exRNA结构域提取和基因计数矩阵的构建；
所述顺序映射包括如下：
2)-A-1)，将去接头后数据映射到人类rRNA数据库；
2)-A-2)，将2)-A-1)处理后未映射的数据继续映射到UniVec数据库；
2)-A-3)，将2)-A-2)处理后未映射的数据依次按照如下顺序映射到关注的RNA类别上：lncRNA、miRNA、mRNA、piRNA、snoRNA、snRNA、srpRNA、tRNA、TUCP RNA和Y RNA；
2)-A-4)，将2)-A-3)处理后未映射的数据映射到人类完整的基因组序列上；
2)-A-5)，将2)-A-4)处理后未映射的数据映射到circBase数据库中定义的环状RNA上；
所述exRNA结构域提取包括如下：将所述2)-A-4)映射后数据的每个转录本划分为20个nt的bin，然后计算每个bin的平均读数覆盖率，并过滤掉平均读数覆盖率低于5的bin；从第一个有效bin的起始位置开始进行双向搜索以获得局部最大的读数覆盖率；一旦找到局部最大值，则从该局部最大值的位置向5’端和3’端两个方向进行峰拓展，直到该位置的读数覆盖率小于5或小于局部最大值的1/2；删除短于10nt的峰，并合并相邻的峰；针对每个样品独立执行了上述峰提取程序，并计算每个峰在所有样品中重复出现的次数，峰重复率大于10％的区域被保留并合并为可信峰，这些可信峰又称为结构域特征，即为exRNA结构域；
所述基因计数矩阵的构建为将所述exRNA结构域、2)-A-3)映射到miRNA的数据和所述2)-A-5映射后数据生成癌症病人基因计数矩阵和正常人基因计数矩阵；
所述2)-B：若exRNA大于50bp，则所述预处理依次包括去接头、顺序映射和基因计数矩阵的构建；
所述顺序映射包括如下：先将去接头数据映射到到人类rRNA数据库；再将未映射的数据映射到人类hg38基因组；再将未映射的数据映射到circBase数据库中定义的环状RNA；
所述基因计数矩阵的构建为将所述2)-B中的顺序映射中的映射到人类hg38基因组和所述映射到circBase数据库的数据生成癌症病人基因计数矩阵和正常人基因计数矩阵；
3)将2)得到的各个样本的基因计数矩阵分别使用如下5种样本库大小归一化方法：CPM、TMM、UQ、RLE和CPM-top进行处理，得到各个样本的5种归一化处理后矩阵；再将所述各个样本的5种归一化处理后矩阵分别用limma、ComBat、和RUV方法进行批次效应校正，得到各个样本15种批次校正后矩阵；检测所述各个样本15种批次校正后矩阵，再从中选择mkNN和UCA两个指标均最接近于1的归一化处理和批次效应校正组合方法作为目标归一化方法和目标批次效应校正方法；最后采用该目标归一化方法和目标批次效应校正方法对2)得到的各个样本的基因计数矩阵进行处理，得到各个样本的归一化和批次效应去除后的exRNA基因表达矩阵；
4)将所述各个样本的归一化和批次效应去除后的exRNA基因表达矩阵进行对数转换以让分布更加集中而接近正态分布；再对每个特征或每一行进行标准化；再用如下7种特征选择方法对上述各个标准化处理得到的exRNA基因表达矩阵进行特征选择，得到各个癌症病人和正常人的7种exRNA基因表达矩阵特征选择结果；再将所述各个癌症病人和正常人的7种exRNA基因表达矩阵特征选择结果分别用逻辑回归、随机森林和决策树规则进行再次特征选择，得到癌症病人21种exRNA基因表达矩阵特征选择结果和正常人的21种exRNA基因表达矩阵特征选择结果；最后从所述癌症病人21种exRNA基因表达矩阵特征选择结果和正常人的21种exRNA基因表达矩阵特征选择结果中选出特征稳定性评估以及智能预测模型结果评估均接近1的作为癌症病人目标exRNA候选生物标志物；
所述7种特征选择方法为(1)差异表达，(2)基于逻辑斯蒂回归的wrapper方法，(3)基于随机森林的wrapper方法，(4)ReliefF，(5)SURF，(6)MultiSURF和(7)SIS。
CN202010618721.4A 2020-06-30 2020-06-30 体液样本中鉴定胞外rna生物标志物的试剂盒及方法 Active CN113862351B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010618721.4A CN113862351B (zh) 2020-06-30 2020-06-30 体液样本中鉴定胞外rna生物标志物的试剂盒及方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010618721.4A CN113862351B (zh) 2020-06-30 2020-06-30 体液样本中鉴定胞外rna生物标志物的试剂盒及方法

Publications (2)

Publication Number Publication Date
CN113862351A CN113862351A (zh) 2021-12-31
CN113862351B true CN113862351B (zh) 2023-04-07

Family

ID=78981680

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010618721.4A Active CN113862351B (zh) 2020-06-30 2020-06-30 体液样本中鉴定胞外rna生物标志物的试剂盒及方法

Country Status (1)

Country Link
CN (1) CN113862351B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114822722B (zh) * 2022-04-27 2022-12-06 吴嘉瑞 一种治疗肝炎的中药药效物质筛选方法
CN117594133A (zh) * 2024-01-19 2024-02-23 普瑞基准科技（北京）有限公司 用于判别子宫病变类别的生物标志物的筛选方法及其应用

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108707663A (zh) * 2018-04-19 2018-10-26 深圳华大基因股份有限公司 用于癌症样本miRNA测序定量结果评价的试剂、制备方法和应用
CN109517881A (zh) * 2018-06-13 2019-03-26 清华大学 一种体液微量游离rna的高通量测序文库构建方法
CN110706749A (zh) * 2019-09-10 2020-01-17 至本医疗科技（上海）有限公司 一种基于组织器官分化层次关系的癌症类型预测系统和方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6232073B1 (en) * 1999-02-05 2001-05-15 Virginia Commonwealth University Nucleic acid marker for cancer

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108707663A (zh) * 2018-04-19 2018-10-26 深圳华大基因股份有限公司 用于癌症样本miRNA测序定量结果评价的试剂、制备方法和应用
CN109517881A (zh) * 2018-06-13 2019-03-26 清华大学 一种体液微量游离rna的高通量测序文库构建方法
CN110706749A (zh) * 2019-09-10 2020-01-17 至本医疗科技（上海）有限公司 一种基于组织器官分化层次关系的癌症类型预测系统和方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
RNA Biomarkers: Frontier of Precision Medicine for Cancer;Xi X等;《Noncoding RNA》;20170220;第3卷(第1期);第1-17页 *
基于主成分分析和神经网络的癌症驱动基因预测模型;周莉;《中国优秀硕士学位论文全文数据库 医药卫生科技辑》;20180115(第01期);第1-51页 *
通过RNA-蛋白相互作用的联合聚类鉴定RNA调控元件;李洋;《中国博士学位论文全文数据库 基础科学辑》;20200615(第06期);第1-112页 *

Also Published As

Publication number Publication date
CN113862351A (zh) 2021-12-31

Similar Documents

Publication Publication Date Title
US20220351805A1 (en) Systems and methods for detecting cellular pathway dysregulation in cancer specimens
CN109072309B (zh) 癌症进化检测和诊断
JP2022521492A (ja) 相同組換え欠損を推定するための統合された機械学習フレームワーク
WO2019232435A1 (en) Convolutional neural network systems and methods for data classification
CN113257350B (zh) 基于液体活检的ctDNA突变程度分析方法和装置、ctDNA性能分析装置
CN112086129B (zh) 预测肿瘤组织cfDNA的方法及系统
CN106156543B (zh) 一种肿瘤ctDNA信息统计方法
US20030017481A1 (en) Methods for classifying samples and ascertaining previously unknown classes
Larsson et al. Comparative microarray analysis
CN112951327B (zh) 药物敏感预测方法、电子设备及计算机可读存储介质
CN112218957A (zh) 用于确定在无细胞核酸中的肿瘤分数的系统及方法
US20200126637A1 (en) Methods for identifying agents with desired biological activity
CN112289376B (zh) 一种检测体细胞突变的方法及装置
EP3973080A1 (en) Systems and methods for determining whether a subject has a cancer condition using transfer learning
CN113862351B (zh) 体液样本中鉴定胞外rna生物标志物的试剂盒及方法
EP4035161A1 (en) Systems and methods for diagnosing a disease condition using on-target and off-target sequencing data
AU2020215312A1 (en) Method of predicting survival rates for cancer patients
CN111733251A (zh) 一种特征miRNA表达谱组合及肾透明细胞癌早期预测方法
CN116312800A (zh) 一种基于血浆中循环rna全转录组测序的肺癌特征识别方法、装置和存储介质
Hobbs et al. Biostatistics and bioinformatics in clinical trials
CN110462056A (zh) 基于dna测序数据的样本来源检测方法、装置和存储介质
Lauria Rank-based miRNA signatures for early cancer detection
Olorunshola Classifying Different Cancer Types Based on Transcriptomics Data Using Machine Learning Algorithms
WO2022140616A1 (en) Taxonomy-independent cancer diagnostics and classification using microbial nucleic acids and somatic mutations
Yim Leveraging of single molecule sequencing methods for less invasive cancer detection

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant