CN116194995A - 在低覆盖度的下一代测序数据中识别染色体空间不稳定性如同源修复缺陷的方法 - Google Patents
在低覆盖度的下一代测序数据中识别染色体空间不稳定性如同源修复缺陷的方法 Download PDFInfo
- Publication number
- CN116194995A CN116194995A CN202180060052.6A CN202180060052A CN116194995A CN 116194995 A CN116194995 A CN 116194995A CN 202180060052 A CN202180060052 A CN 202180060052A CN 116194995 A CN116194995 A CN 116194995A
- Authority
- CN
- China
- Prior art keywords
- coverage
- data
- chromosome
- hrd
- sample
- 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.)
- Pending
Links
Images
Classifications
-
- 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
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N20/00—Machine learning
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/02—Neural networks
- G06N3/04—Architecture, e.g. interconnection topology
- G06N3/0464—Convolutional networks [CNN, ConvNet]
-
- 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
- G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
-
- 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
- G16B35/00—ICT specially adapted for in silico combinatorial libraries of nucleic acids, proteins or peptides
-
- 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
- G16B40/00—ICT specially adapted for biostatistics; ICT specially adapted for bioinformatics-related machine learning or data mining, e.g. knowledge discovery or pattern finding
- G16B40/20—Supervised data analysis
-
- 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
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Health & Medical Sciences (AREA)
- Theoretical Computer Science (AREA)
- Medical Informatics (AREA)
- Bioinformatics & Cheminformatics (AREA)
- General Health & Medical Sciences (AREA)
- Biophysics (AREA)
- Biotechnology (AREA)
- Bioinformatics & Computational Biology (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Evolutionary Biology (AREA)
- Software Systems (AREA)
- Data Mining & Analysis (AREA)
- Chemical & Material Sciences (AREA)
- Molecular Biology (AREA)
- Evolutionary Computation (AREA)
- Artificial Intelligence (AREA)
- Proteomics, Peptides & Aminoacids (AREA)
- Analytical Chemistry (AREA)
- Bioethics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Databases & Information Systems (AREA)
- Mathematical Physics (AREA)
- General Physics & Mathematics (AREA)
- Computing Systems (AREA)
- General Engineering & Computer Science (AREA)
- Genetics & Genomics (AREA)
- Epidemiology (AREA)
- Public Health (AREA)
- Library & Information Science (AREA)
- Biochemistry (AREA)
- Biomedical Technology (AREA)
- Computational Linguistics (AREA)
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
- Investigating Or Analysing Biological Materials (AREA)
- Apparatus For Radiation Diagnosis (AREA)
- Image Processing (AREA)
Abstract
基因组数据分析仪可被配置为用机器学习模型,如经训练的卷积神经网络来检测并表征肿瘤样本中基因组不稳定性的存在。基因组数据分析仪可以使用全基因组测序读数作为输入数据,即使是在各种临床肿瘤学设置中经常采用的高通量测序工作流程中的低测序覆盖度下。基因组数据分析仪可以排列来自染色体臂或全染色体的比对读数数据覆盖度,以形成可能是图像形式的覆盖度数据信号阵列。经训练的机器学习模型可以处理覆盖度数据信号阵列,以确定肿瘤样本中是否存在染色体空间不稳定性(CSI),例如由同源修复或重组缺陷(HRD)引起的基因组不稳定性。后者的指示可以指导对肿瘤的首选抗癌治疗的选择。
Description
技术领域
本文描述的方法一般涉及基因组分析,更具体而言,涉及一种基因组信息在检测和治疗癌症中的应用。
背景技术
作为抗癌治疗应答的预测器的肿瘤基因组学
除了固有的种系突变外,癌细胞还经常携带体细胞染色体大尺度的畸变,如某些等位基因或基因组区域的拷贝数变异、重复或缺失。其中一些变体可能导致丧失一些基因组功能,例如肿瘤抑制机制,特别是同源重组修复(HRR或HR)功能,从而使癌症更具侵袭性。识别这些基因组突变在个性化癌症药物治疗的最近进展中具有关键意义,因为它们已被证明可以预测患有细胞过度增殖症的受试者对某些抗癌疗法的应答。受试者可以是人或动物。抗癌治疗的实例包括这样的烷基化剂,例如,但不限于,铂基的化疗剂、卡铂、顺铂、异丙铂、奈达铂、奥沙利铂、吡铂、氮芥、苯丁酸氮芥、美法仑、环磷酰胺、异环磷酰胺、雌莫司汀、卡莫司汀、洛莫司汀、福莫司汀、链脲菌素、白消安、哌泊溴烷、甲基苄肼、达卡巴嗪、塞替派、替莫唑胺和/或其他抗肿瘤的铂配合化合物;DNA损伤剂或辐射;蒽环类药物,例如,但不限于,表柔比星或多柔比星;拓扑异构酶I抑制剂,例如,但不限于,喜树碱、拓扑替康、伊立替康;和/或PARP(聚ADP-核糖聚合酶)抑制剂。利用肿瘤DNA同源修复缺陷(HRD)选择性地破坏癌细胞的PARP抑制剂的实例是:奥拉帕尼、卢卡帕尼、尼拉帕尼(MK 4827)和他拉唑帕尼(BMN-673),它们在美国和欧洲已经被批准用于某些癌症类型;PARP抑制剂的其他实例包括:伊尼帕利(iniparib)、CEP 9722(-14-甲氧基-9-[(4-甲基哌嗪-1-基)甲基]-9,19-二氮杂五环[10.7.0.0^{2,6}.0^{7,11}.0^{13,18}]十九-1(12),2(6),7(11),13(18),14,16-己烯-8,10-二酮)、3-氨基苯甲酰胺、维拉匹里(velapirib)、帕米帕利(Pamiparib)或E7016(10-((4-羟基哌啶-1-基)甲基)色烯并[4,3,2-de]酞嗪-3(2H)-酮)。PARP抑制剂间接诱导大量的DNA双链断裂(DSB),这防止HRD肿瘤细胞的发展,而正常细胞一般能够通过HR修复这些断裂(Keung等人,2019,Journal of Clinical Medicine 8(4),p.435)。
HRD基因组分析测试
因此,一种对DNA样本(如来自肿瘤的样本)进行分类的方法可能有助于对可能的癌症类型进行诊断,或者由于基因组测序和分析而可以根据患者肿瘤样本中DNA拷贝数严重畸变的特征对其进行最适合的癌症治疗。特别是,确定癌症是否是同源修复缺陷(HRD)可能对治疗计划有相当大的帮助。在过去的十年中,为此目的已经确定了不同的基因组突变特征。欧洲专利EP2609216B1披露了使用全球染色体畸变分数(GCAS)来预测抗癌治疗如PARP抑制剂、放疗或包括铂类化疗药物的化疗的结果。欧洲专利EP2817630B1提出检测端粒等位基因不平衡(TAI)事件的数量,如果该数量高于已知耐铂的类似癌症的参考值,则选择含铂的疗法。居里研究所和INSERM的斯特恩、马尼和波波娃的欧洲专利EP2859118B1披露了一种通过计数每个基因组中至少跨越3兆碱基的片段的数量来预测HRD的方法,该片段的数量相当于染色体拷贝数的大规模转换(LST)。来自Myriad Genetics的Abkevich等人的欧洲专利EP2794907B1披露了对至少一对人类染色体中长于11兆碱基但短于整个染色体的杂合性缺失(LOH)区域的总数进行计数,并将该数字与参考数值相比较,以预测患者对各种可能的癌症治疗的应答;EP2981624B1公开了使用LOH、TAI和LST指标,而EP3180447B1公开了将它们的值合并成为HRD分数,将其与参考数值进行比较,以预测患者对各种可能的癌症治疗的应答。后一种检测方法目前被Myriad Genetics myChoice CDx测定用作测试卵巢癌患者的HRD,其在中位覆盖度至少为500X的许多基因上采用基于定制杂交的靶向富集。这种HRD分数最近也被证明是关于高级别浆液性卵巢癌(HGSOC)进展的可能的预后预测因素(Takaya等人,“Homologous recombination deficiency status-based classificationof high-grade serous ovarian carcinoma”,Nature Research Scientific Reports(2020)10:2757)。
与Myriad myChoice测定相似,Foundation Medicine的商业测定FoundationFocus CDxBRCALOH专门分析BRCA1和BRCA2基因,以确定卵巢癌患者是否会对卢卡帕尼(一种PARP抑制剂)疗法有应答。后者测定也使用了基于杂交的定制捕获许多基因的所有编码外显子,中位覆盖度为500X,但只使用了从全基因组拷贝数特征和SNP的小等位基因比例组合中评估的LOH分数。
WGS数据的HRD分析方法
更一般而言,WO2017191074提出了根据从分析不同的基因组改变,包括碱基替换、重排和印迹特征(HRDetect分数),计算出的概率分数来表征肿瘤DNA样本的HRD状态。如于2017年3月13日在线发表的“HRDetect is apredictor of BRCA1 and BRCA2 deficiencybased on mutational signatures”,H.Davies等人,Nature Medicine中所述,通过后者的HRDetect分数,来自全基因组观察的突变特征,全基因组测序(WGS)被证明是检测HRD的可能方法,灵敏度为98.7%。与WGS相比,当只应用全外显子组测序(WES)时,HRDetect的灵敏度明显下降,介于46.8%和73%之间,这取决于生物信息学算法的具体调整。
除了乳腺癌和卵巢癌,正如“Genomic aberration based molecular signaturesefficiently characterize homologous recombination deficiency in prostatecancer”,Sztupinszki等人最近也开始研究使用LOH、TAI和LST作为WGS和WES数据中的HRD签名n个前列腺肿瘤的指标(scarHRD)。
由于上述方法都依赖于SNV和/或INDEL的调用,它们在下一代测序(NGS)工作流程中需要高覆盖深度(通常至少30X)。WGS和大样本中的高覆盖度要求大大增加了临床实践中的分析成本,包括湿法实验室实验和干法实验室数据处理的开销。因此,使用低通量WGS(LP-WGS-1x至5x)或甚至超低通量WGS(ULP-WGS-低至0.1x)的替代方法可能对临床肿瘤学更有利。在“ShallowHRD:Detection of Homologous Recombination Deficiency fromshallow Whole Genome Sequencing“,Bioinformatics,Apr 21 2020,Eeckhoutte等人中,描述了ShallowHRD,这是一种表征LST状态的软件方法,类似于Manie,Stern和Pova的方法,但专门针对浅层WGS数据,覆盖度低至约1X。他们的方法只是使用染色体臂内拷贝数改变的计数。这个计数可以估计为在覆盖度数据信号中除去小于3Mb的片段后,至少10Mb的相邻片段之间的大规模转换的数量。据作者所说,少于15个转换能够识别HRD阴性肿瘤,而超过19个转换能够识别HRD阳性肿瘤,其灵敏度和特异性与现有技术的scarHRD方法相似。虽然这种方法为在临床肿瘤学中使用更具成本效益的NGS工作流程提供了希望,但我们注意到它的灵敏度仍然低于HRDetect和SNP阵列上LST指标的测量所取得的结果。这可能是由于它不能整合现有技术方法中的LOH和TAI特征。此外,它只能够对那些染色体改变大且频繁,足以造成十几个大规模转变的肿瘤进行分类,而不考虑这些转变一方面在每个染色体臂内如何分布,另一方面在人类基因组的多条染色体之间如何分布。
最近在肿瘤学方面的一些研究工作表明,功能失调的DNA的高突变性可能会在染色体的不同区域引起异质性改变。在“Regulation of mitotic recombination betweenDNA repeats in centromeres”,Nucleic Acids Research,2017,Vol.45,No.19,Zafaret等人中,在“The dark side of centromeres,types,causes and consequences ofstructural abnormalities implicating centromeric DNA”,Barra et al.,NatureCommunications(2018)9:4340,Barra等人报道了在衍生自结直肠癌和腺癌的癌细胞系中,观察到染色体重排和断裂的比例明显,这可能是由于肿瘤中着丝粒区的固有脆弱性,这在其他癌症中也很常见。Barra等人强调缺乏模型和技术来明确分析这些区域的高度重复序列。
因此,需要改进的基因组分析方法,这些方法可以很容易地、经济有效地部署在常规肿瘤学临床实践的自动化NGS工作流程中,同时适应最近的癌症基因组学发现,特别是适应于肿瘤基因组中某些染色体畸变的空间特征的有希望的发现。特别需要改进的基因组分析方法,整合最近的癌症基因组学发现,并允许识别与诊断、预后、治疗选择和患者管理有关的特征,包括HRD。更好的是,这些方法和技术也可以在较低的覆盖度下适用,以减少常规临床设置中的整体分析成本,同时与现有技术相比在检测HR事件方面保持合理的准确性。此外,改进的基因组分析方法还可以通过采用机器学习技术,而不是依赖现有技术中的LOH、LST或TAI指标等标量HRD指标的明确测量和阈值来促进HR缺陷肿瘤的表征,从而可以更容易地利用越来越多的临床数据对人类和动物的各种不同癌症的肿瘤进行分类。
发明内容
提出了一种确定受试者DNA样本的同源重组缺陷(HRD)状态的方法,所述方法包括:获得一组待分析的受试者DNA样本的全基因组测序读数,将这一组受试者DNA样本的测序读数与参考基因组进行比对,其中,所述参考基因组被分为多个区块(bin),每个区块属于待分析的全基因组染色体中染色体臂的同一基因组区域,对沿着每个染色体臂的每个区块中的比对读数的数值进行计算并归一化以获得所述染色体臂的覆盖度信号,将染色体臂的覆盖度信号排列成受试者DNA样本的覆盖度数据信号阵列,将覆盖度数据信号阵列输入到经训练的机器学习模型中,其中,所述模型已通过使用一组已知同源重组缺陷状态的样本进行训练,以区分来自同源重组缺陷状态为阳性的样本的覆盖度数据信号阵列和来自同源重组缺陷状态为阴性的样本的覆盖度数据信号阵列,从而确定受试者DNA样本的同源重组缺陷分数(HRD分数),以及根据所述经训练的机器学习模型的HRD分数来确定受试者DNA样本的阴性、阳性或不确定的同源重组缺陷(HRD)状态。在可能的实施方式中,这一组测序读数可以是从全基因组测序中获得的,其中,读数深度覆盖度至多为30x;或者这一组测序读数可以是从低通量全基因组测序中获得的,其中,读数深度覆盖度至少为0.1X且至多为5X。在可能的实施方式中,对沿着每个染色体臂的每个区块中的比对读数的数值进行计算并归一化以获得所述染色体臂的覆盖度信号可以包括对每个样本的覆盖度信号进行归一化,和/或通过GC含量进行归一化以应用GC-偏移校正。在可能的实施方式中,所述染色体臂的覆盖度信号可以被排列成一维覆盖度数据信号矢量或二维覆盖度数据信号图像。在可能的实施方式中,所述染色体臂的覆盖度信号可以通过按行比对每条染色体的覆盖度数据信号和每条染色体臂的着丝粒区块,即,与染色体臂的着丝粒区域相邻的最近区块,而被排列成二维覆盖度数据信号图像。在可能的实施方式中,所述机器学习模型可以是使用一组具有已知同源重组缺陷状态的肿瘤数据样本作为训练标签进行预先训练的。所述训练数据集可以用通过结合数据样本的染色体的数据与已知同源重组缺陷状态标签而产生的人工样本数据进行增强。生成数据增强样本可以是为了表示在真实样本数据集中观察到的纯度-倍性比分布。在可能的实施方式中,参考基因组可以被划分为第一组至多100kbp的区块,并还包括在排列每个染色体臂上的覆盖度信号之前,将100kbp区块坍塌成第二组至少500kbp的较大区块的步骤。第一组区块的区块可以具有至多100kbp的均匀尺寸,而第二组区块的区块可以具有2.5至3.5Mbp的尺寸并且通过池化第一组区块中25至35个100kbp的区块得到。
提出了一种确定患者DNA样本的同源重组缺陷(HRD)状态的体外方法,所述方法包括:
-提供来自患者样本的DNA片段;
-构建包含与一组染色体重叠的所述片段的文库;
-对所述文库进行测序,以达到至多30X的全基因组测序覆盖度,优选达到至少0.1X且至多5X的基因组测序覆盖度;
并且基于根据本文公开的方法获得的经训练的机器学习模型的分析来确定患者样本的HRD状态。所述患DNA样本可以是肿瘤细胞外游离DNA(cfDNA)、新鲜冷冻组织(FFT)或福尔马林固定石蜡包埋(FFPE)样本。所述患样本的HRD分数或HRD状态可以是对癌症治疗方案的肿瘤应答的预测因子。所述的癌症治疗方案可以选自由烷基化剂、铂类化疗剂、卡铂、顺铂、异丙铂、奈达铂、奥沙利铂、吡铂、氮芥、苯丁酸氮芥、美法仑、环磷酰胺、异环磷酰胺、雌莫司汀、卡莫司汀、洛莫司汀、福莫司汀、链脲菌素、白消安、哌泊溴烷、甲基苄肼、达卡巴嗪、塞替派、替莫唑胺和/或其他抗肿瘤的铂配合化合物、DNA损伤剂、放射治疗剂、蒽环类药物、表柔比星、多柔比星;拓扑异构酶I抑制剂、喜树碱、拓扑替康、伊立替康、PARP(聚ADP-核糖聚合酶)抑制剂、奥拉帕尼、卢卡帕尼、尼拉帕尼、他拉唑帕尼、伊尼帕利、CEP 9722、MK4827、BMN-673、3-氨基苯甲酰胺、维拉匹里、帕米帕利或E7016所组成的组。
提出了一种选择癌症患者用铂类化疗剂、DNA损伤剂、蒽环类药物、拓扑异构酶I抑制剂、PARP抑制剂进行治疗的方法,包括根据本文公开的方法检测肿瘤患者样本为HRD阳性的步骤。在可能的实施方式中,所述患者可以患有选自高级别浆液性卵巢癌、前列腺癌、乳腺癌或胰腺癌的癌症。
提出了一种训练机器学习算法以确定受试者DNA样本的同源重组缺陷(HRD)状态的方法,所述方法包括向机器学习监督训练算法输入来自已知同源重组缺陷状态为阳性的样本的覆盖度数据信号阵列和来自已知同源重组缺陷状态为阴性的样本的覆盖度数据信号阵列。
经训练的机器学习模型可以是随机森林模型、神经网络模型、深度学习分类器或卷积神经网络模型。神经网络模型经训练的机器学习模型可以是卷积神经网络,所述卷积神经网络被训练成在其输出端产生阳性或阴性的HRD状态的单标签二元分类,或阳性、阴性或不确定的HRD状态的单标签多类别分类,或代表HRD状态的标量HRD分数。所述机器学习模型可以使用通过从共享相同HRD状态及相同归一化纯度和倍性比的一组真实样本的染色体上取样数据产生的数据增强集在半监督模式下训练。
提出了一种用于表征DNA样本的方法,所述方法包括从患者样本中获得一组与所述参考基因组比对的测序读数;将参考基因组的至少两个染色体的碱基对位置(bp)划分为一组区块,每个区块对应于至多20兆字节对(20Mbp)的基因组区域,每个区块仅包含来自单个染色体臂的覆盖度数据;从比对的读数中估计已测序患者样本的每个区块的覆盖度数据;将覆盖度数据排列成多维阵列,该阵列包括沿一个维度的染色体或染色体臂,以及沿另一个维度的所述染色体或染色体臂的一组区块,其特征在于,每个染色体或染色体臂的着丝粒区块或端粒区块被排列成多维阵列的空间排列;将多维阵列输入到经训练的机器学习模型;以及在经训练的机器学习模型的输出端处生成染色体空间不稳定性(CSI)指标。在可能的实施方式中,CSI指标可以是所述患者样本具有同源重组(HR)缺陷的高(HRD+)或低(HRD-)可能性的指标。患者样本可以是肿瘤样本,而患者样本的CSI指标可以是对癌症治疗方案的肿瘤应答的预测因子,癌症治疗方案包括铂类化疗剂、DNA损伤剂、蒽环类药物、拓扑异构酶I抑制剂或PARP抑制剂。
在一个实施方式中,机器学习模型是通过使用来自具有已知基因组不稳定状态,如HRD状态(阳性或阴性)的样本的标记的覆盖度数据在监督或半监督模式下进行训练的。
在一个实施方式中,基因组不稳定性的染色体空间不稳定性(CSI)指标是基于在一个DNA样本中至少一个染色体臂的至少一个区域具有基因组不稳定性的所述DNA样本与在另一个DNA样本中至少一个染色体臂的至少一个区域没有基因组不稳定性的所述另一个DNA样本之间的差异,这可以由机器学习模型学习。
在可能的实施方式中,第一组区块可以在排列多维阵列之前通过相对于每个染色体臂的长度调整较大区块的尺寸来进一步坍塌成较大区块,使得区块尺寸沿染色体臂保持恒定。第一组区块可以有至多100kbp的均一尺寸,而塌陷的区块可以有至少500kbp的尺寸。覆盖度数据可以被归一化和/或可以被分割,以便在排列多维阵列之前推断出离散的拷贝数值。在可能的实施方式中,多维阵列的空元素可以用最接近的区块的值或预定的值来填充。在可能的实施方式中,经训练的机器学习模型可以是随机森林模型或神经网络模型,例如卷积神经网络,其被训练成在其输出端产生阳性或阴性的CSI状态的单标签二元分类,或阳性、阴性或不确定的CSI状态的单标签多类别分类,或代表CSI状态的标量HRD分数。患者样本可以是肿瘤样本,机器学习模型可以在监督或半监督模式下使用一组根据目标应用用CSI状态标记的真实样本来训练。在可能的实施方式中,CSI状态标签可以是根据HRDetect方法和/或使用通过对一组共享相同HRD状态并具有类似肿瘤内容的真实样本的染色体进行采样而产生的人工样本的HRD状态。人工样本的选择可以是为了重新创造与真实样本数据集相同的纯度-倍性分布。
在一个特定的实施方式中,患者样本可以是肿瘤样本,并且机器学习模型可以通过使用具有HRD状态的标记数据在监督或半监督模式下进行训练。训练数据可以由从具有HRD状态的临床样本中获得的并排列成多维阵列的覆盖度数据组成。训练数据标签可以从HRD检测方法中获得,例如HRDetect方法。
在一个特定的实施方式中,考虑到由样本特定属性(如纯度和倍性)引入的偏差的数据增强策略,能够用于增加用于训练机器学习算法的多维阵列的多样性和数量。
附图说明
图1表示根据本公开的某些实施方式的下一代测序系统。
图2表示根据本公开的某些实施方式的基因组分析工作流程。
图3示意性示出了低覆盖度数据准备预处理工作流程的可能的实施方式,包括将小覆盖度区块坍塌和平滑成大区块,以更好地证明沿不同染色体臂可观察到的更具体的空间事件。
图4描绘了沿着人类DNA样本中的22条非性染色体的归一化读数覆盖度数据信号的实例,其可以是用低通全基因组测序实验的覆盖度数据信号中的第一组小区块位初步测量的。
图5绘制了归一化读数覆盖度数据信号的实例,其可以用第二组较大区块从覆盖度数据信号图4中计算出来的。
图6示出了排列成对齐它们的着丝粒区域的整套人类基因组染色体(详见x轴),以及它们各自的p臂和q臂染色体臂长(y轴)。
图7示意了根据所提出的覆盖度数据制备方法的某些实施方式,在空间上被重新排列为染色体按行排列而且覆盖区块沿列排列的二维矩阵的归一化覆盖度数据,使得在染色体空间不稳定性分析之前着丝粒区块在单列中竖直排列,其中,a)HRD阴性肿瘤DNA样本,和b)HRD阳性肿瘤DNA样本。
图8示意性地示出了可能的染色体空间不稳定性分析仪(例如,用于HRD状态分析)的内部数据工作流程架构。
图9示出了根据所提出的CSI分析仪数据处理模块的某些实施方式的卷积神经网络的可能架构。
图10a)和图10b)示出了560乳腺癌基因组数据库的2个样本的归一化覆盖度数据。
图11示出了来自560乳腺癌基因组数据库的202个BAM文件样本的关于22条非性染色体的归一化覆盖度数据。
图12将现有技术的HRD分数方法与提议的CSI指标分类器的HRDetect分数进行比较。
图13绘制了一系列69个具有不同BRCA缺失(可能与HR缺失有关)测试样本的HRD分数、提议的CSI指标分数和HRDetect分数的结果。
图14示出了来自至少两个分析的示例性染色体(染色体A和染色体B)的每个染色体臂(p臂和q臂)相对于至少一个与着丝粒区域相邻的区块或一个与端粒区域相邻的区块的覆盖度信号图像的示例性比对情况。图片A是至少两个示范性分析染色体(染色体A和染色体B)的每个染色体臂(p臂和q臂)的覆盖度信号图像,有点的矩形为与端粒区相邻的区块,竖直条纹矩形为与着丝粒区相邻的区块;图片B是覆盖度信号图像关于与着丝粒区域相邻的区块的比对情况;图片C是覆盖度信号图像关于与端粒区相邻的区块的比对情况;图片D是覆盖度信号图像关于与着丝粒区相邻的区块的比对情况;图片E是覆盖度信号图像关于与端粒区相邻的区块的比对情况。
图15示出了如实施例2所述的从结合的lpWGS和定向捕获所获得的示例性输出数据。图15A示出了从结合低通WGS(lpWGS)和靶向测序的工作流程中产生的NGS数据所获得的原始覆盖度,以100kB的间隔进行计算。出于说明目的,在靶向实验中富集的基因组区域的原始覆盖度用空标记表示,而lpWGS区域用实心标记表示。图15B示出了在非目标(仅lpWGS)区域(左,图15B的B1)和目标(捕获富集)区域(右,图15B的B2)的映射读数的排列。图15C示出了仅限于lpWGS区域的归一化覆盖度(左图,图15C的C1),以及在捕获富集区域识别的变体的等位基因部分(右图,图15C的C2)。
图16示出了如实施例3所述的人工样本的制备。图片A示出了跨越随机选择以生成人工样本的3个原始肿瘤样本的基因组(X轴,归一化覆盖度图)的归一化覆盖度(Y轴)。每个样本的肿瘤含量(或纯度,T)、倍性(P)和纯度倍性比(PPR)在样本的相应的归一化覆盖度图上方标出。图片B示出了对样本进行于电脑(in-silico)中稀释后的归一化覆盖度图,以匹配在该组的原始样本中观察到的最低PPR,图片C示出了通过结合稀释后的样本覆盖度图的数据所组装的人工样本的归一化覆盖度图。原始样本、稀释样本和人工样本的归一化覆盖点的颜色被保留下来。
具体实施方式
本公开至少部分基于以下发现:本公开的机器学习训练的分析仪,旨在处理与基础染色体排列有关的肿瘤样本测序数据覆盖度,可以提取一般的染色体空间不稳定性(CSI)的指标,特别是肿瘤样本的同源修复缺陷(HRD)状态的指标。
现在将通过参考更详细的实施方式来描述所提出的方法和系统。然而,所提出的方法和系统可以以不同的形式体现出来,不应该被理解为局限于本文所述的实施方式。相反,提供这些实施方式是为了使本公开彻底和完整,并向本领域的技术人员充分表达其范围。
定义
“DNA样本”是指来自生物体的核酸样本,例如可以从实体肿瘤或液体中提取。该生物体可以是人类、动物、植物、真菌或微生物。核酸可以在固体样本中发现,如福尔马林固定的石蜡包埋(FFPE)样本。另外,核酸也可以在有限的数量或低浓度,如血液或血浆中的循环肿瘤DNA中找到。
“DNA片段”是指由高分子量DNA的片段化产生的短的DNA片段。片段化可能是在样本生物体内自然发生的,也可能是通过应用于DNA样本的DNA片段化方法人为产生的,例如通过机械剪切、超声处理、酶解片段化和其他方法。在片段化后,DNA片段可以进行末端修复,以确保每个分子拥有钝端。为了提高连接效率,可以在片段化的DNA的每个3'钝端添加腺嘌呤,使DNA片段能够连接到具有互补的dT-突出端的接头。
“DNA产品”是指通过操作、延伸、连接、加倍、扩增、复制、编辑和/或切割DNA片段以使其适应下一代测序工作流程而产生的工程化的DNA片段。
“DNA接头产品”是指将DNA片段与DNA接头连接以使其与下一代测序工作流程兼容而产生的DNA产品。
“DNA文库”是指DNA产品或DNA接头产品的集合,以使DNA片段与下一代测序工作流程兼容。
“池”是指来自相同或不同生物体的多个DNA样本(例如48个样本、96个样本或更多),这些样本可被多路复用到单个高通量测序分析中。每个样本都可以通过独特的样本条形码在池中识别。
“核苷酸序列”或“多核苷酸序列”是指核苷酸,如胞嘧啶(由序列串中的C字母表示)、胸腺嘧啶(由序列串中的T字母表示)、腺嘌呤(由序列串中的A字母表示)、鸟嘌呤(由序列串中的G字母表示)以及尿嘧啶(由序列串中的U字母表示)的任何聚合物或低聚物。它可以是DNA或RNA,或其组合。它可以永久地或暂时地以单链或双链的形式存在。除非另有说明,核酸序列是按5'到3'的方向从左到右书写的。
“扩增”是指多核苷酸扩增反应,以产生从一个或多个母体序列复制的多个多核苷酸序列。扩增可以是通过各种方法,例如聚合酶链式反应(PCR)、线性聚合酶链反应、基于核酸序列的扩增、滚动循环扩增和其他方法产生的。
“测序”是指从DNA文库中读出核苷酸的序列,以产生能够由生物信息学计算机在生物信息学工作流程中处理的一组测序读数。高通量测序(HTS)或下一代测序(NGS)是指多个序列的实时并行测序,通常每个序列在50到几千个碱基对之间。典型的NGS技术包括Illumina、Ion Torrent Systems、Oxford Nanopore Technologies、Complete Genomics、Pacific Biosciences、BGI和其他公司的技术。根据实际的技术,NGS测序可能需要用测序接头或引物进行样本制备以促进进一步的测序步骤,以及扩增步骤,以便对单一母体分子的多个实例进行测序,例如在通过合成测序的情况下,在输送到流动池之前用PCR扩增进行测序。
“测序深度”或“测序覆盖度”或“测序的深度”是指基因组被测序的次数。在目标富集工作流程中,整个基因组中只有一小部分感兴趣的区域被测序,因此增加测序深度而不面临太大数据存储和数据处理的开销可能是合理的。在一些不需要沿基因组高分辨率的基因组分析应用中,例如在检测拷贝数改变中,低通(LP)覆盖度(1x-10x)或甚至超低通(ULP)覆盖度(<1x-并非所有位置都被测序)在信息技术基础设施成本方面可能更有效,但这些工作流程需要更复杂的生物信息学方法和技术来处理从测序仪和比对仪输入的不太可靠的数据。此外,除了与数据存储和处理有关的较高成本外,还需要通过平衡覆盖深度和常规临床工作流程中可能平行检测的样本数量来优化实验性NGS运行的操作成本,即用测序用样本加载测序仪。事实上,下一代测序仪在一次实验中(即在给定的运行中)能产生的读数总数仍然有限。覆盖度越低,用于基因组分析的每个样本的读数就越少,而在下一代测序运行中可复用的样本数量就越多。
“比对(aligning)”或“比对(alignment)”或“比对仪(aligner)”是指根据应用情况在生物信息学工作流程中,对测序读数与参考基因组序列进行逐个碱基的映射和比对。例如,在靶向富集应用中,根据实验扩增过程中使用的杂交捕获探针,测序读数预期映射到特定的靶向基因组区域,比对可以相对于相应的序列进行专门搜索,该序列由基因组坐标定义,如参考基因组中的染色体编号、起始位置和终止位置。如生物信息学实践中已知的,在一些实施方式中,本文采用的“比对”方法还可能包括某些预处理步骤,以促进测序读数的映射和/或从读数中去除不相关的数据,例如通过去除非配对读数,和/或通过修剪作为读数末端的接头序列,和/或其他读数预处理过滤手段。具有不同坐标系统(绝对或相对位置索引,基于0或基于1等)的示例性生物信息学数据表示方法包括BED格式、GTF格式、GFF格式、SAM格式、BAM格式、VCF格式、BCF格式、Wiggle格式、GenomicRanges格式、BLAST格式、GenBank/EMBL特征表格式,以及其他。
“覆盖度”或“序列读数覆盖度”或“读数覆盖度”是指与基因组位置或一组基因组位置对齐的测序读数的数量。一般而言,具有较高覆盖度的基因组区域与下游基因组表征中的较高可靠性有关,特别是在调用变体时。
“区块(Bin)”、“基因组片段”、“分区”、“基因组部分”或“染色体部分”指的是基因组中感兴趣的连续区域。由于这样的区域可能包括变体,区块一般是指基因组的位置或区域,而不是指固定的DNA序列。在生物信息学方法和过程中,区块可以通过其沿参考基因组的起点和终点基因组坐标来识别,并且区块的长度可以用从起点到终点坐标的碱基数(b,kb,Mb)或碱基对(bp,kbp,Mbp)来衡量。一般而言,区块可以对应于整个染色体、染色体的片段、参考基因组的片段、多个染色体的部分、多个染色体、多个染色体的部分和/或其组合。优选的是,区块是通过将参考基因组分割成基因组部分(例如,按尺寸、片段、连续区域、任意定义尺寸的连续区域等进行分割)而获得的染色体的部分。可以使用生物信息学中已知的任何合适的标准对基因组部分进行选择、分类、过滤和/或从考虑中删除。沿着基因组,区块可以有相同的、均一长度或不同的、可变的长度。
“覆盖度区块计数”或“覆盖度计数”或“计数”指的是测序读数或对端读数(代表DNA片段)的数量,其被映射到区块或部分地与区块重叠。计数可以从部分或全部的原始序列读数和/或预处理的序列读数映射或比对到(即有关的)基因组部分而得出。根据本领域已知的各种生物信息学方法,一些组的读数可以在计数前以读数为单位或以对端读数为单位进行加权、移除、过滤、归一化、调整、坍塌、合并、添加和/或减去,或通过其组合进行处理。在一些实施方式中,计数可与不确定度或质量值相关。在一些实施方式中,读数或对端读数可能不完全包含在单个区块内,而是在两个相邻的区块内重叠;在这种情况下,读数可能被计算在具有最大重叠的区块内。
“覆盖度数据信号”或“覆盖度信号”或“覆盖度数据信号阵列”是指与各自计数相关的区块的集合,它可以排列成一维矢量(从而获得覆盖度信号矢量或覆盖度信号阵列或一维热图),二维矩阵(从而获得覆盖度信号矩阵或覆盖度信号图像或二维热图),或任何合适的拓扑结构。覆盖度数据信号可以选择性地被归一化,以消除技术偏差。覆盖度数据信号阵列可以排列成二维表示,在本文被称为覆盖度数据信号图像或覆盖度数据信号多维阵列。
“染色体臂”是指染色体的两个部分(臂)中的任何一个,它们由染色体的着丝粒相互结合。p臂是指最短的臂,而q臂是指最长的臂。每条臂都以端粒为终点。在中间着丝粒染色体中,p臂和q臂的尺寸相似。在亚中间着丝粒染色体中,p臂比q臂短。在近端着丝粒染色体中,p臂非常短。在端着丝粒的染色体中,p臂不再存在或短到在检查染色体时不再可见。人类的正常细胞不携带端着丝粒染色体;然而,它们可能在某些肿瘤细胞中被发现。
染色体的“着丝粒的区域”或“着丝粒区”是指紧挨着染色体的着丝粒的区域。
染色体覆盖度数据信号的“着丝粒区块”是指在与染色体臂相关的覆盖区块集合中与着丝粒区域最接近的区块的覆盖区块。至多可以有两个着丝粒区块与任何给定的人类染色体的基因组覆盖度数据相关:沿着测序数据基因组坐标系统,一个在着丝粒的左边,一个在着丝粒的右边。
染色体的“端粒区域”是指染色体的端粒旁边的区域。
染色体覆盖度数据信号的“端粒区块”是指在与该染色体相关的覆盖区块集合中分别与该染色体起始和末端的两个端粒区域中的任何一个最接近的覆盖区块。对于任何给定的人类染色体,至多可以有两个独立的、遥远的端粒区块与基因组覆盖度数据相关联,一个是p臂,另一个是q臂。
“染色体不稳定性”或“CIN”是指一类基因组不稳定,其中染色体是不稳定的,使得整个染色体或染色体的部分被复制或删除。更具体而言,CIN是指整个染色体或部分染色体的增加或丢失的速度增加。这种基因组改变可能特别发生在肿瘤细胞中,涉及整个染色体的增加或丢失或结构畸变,如染色体的严重重排。
“染色体空间不稳定性”或“CSI”是指染色体不稳定性,其中增加或丢失的速度增加可根据事件在基因组上的空间分布来进行表征。可以理解的是,CSI信息可以在经训练的机器学习模型的输出端获得,或者直接作为机器学习模型为给定表型计算的CSI分数,例如标量数字获得;或者作为CSI状态获得,它可以被标记为正值“CSI+”,负值“CSI-”,或者可能不确定或未确定的值“未确定的CSI”。也可以根据CSI分数推断出CSI的阳性或阴性状态。
“同源重组缺陷状态”或“HRD状态”是指同源重组途径的分类,并且与任何导致同源重组途径缺陷的细胞状态/事件有关。HRD状态可分为阳性(HRD+),其中同源重组途径存在缺陷;或可分为阴性(HRD-),其中同源重组途径不存在缺陷;或可分为未确定的其他状态(HRD不确定,HRD未知)。
“机器学习模型”是指数据模型或数据分类器,其为使用数据科学领域中已知的监督、半监督或无监督学习技术训练出来的,而不是明确的统计模型。数据输入可以表示为一维信号(矢量)、二维信号(矩阵),或更普遍的多维阵列信号(例如张量,或表示为其红、绿、蓝颜色分解平面的3*2D信号的RGB彩色图像-3个矩阵),和/或其组合。多维阵列在数学上的定义是沿至少两个维度排列的数据结构,每个维度记录了大于1的值。
在深度学习分类器的情况下,数据输入通过一系列的数据处理层进一步处理,以隐含地捕获隐藏的数据结构、数据签名和基本模式。由于使用了多个数据处理层,深度学习促进了自动数据处理的普遍化,使之成为复杂模式检测和数据分析任务的多样性。机器学习模型可以在有监督、半监督或无监督的学习框架内进行训练。在监督学习框架内,模型根据输入和匹配输出的实例对学习函数,来映射来自输入数据集的输出结果。用于监督学习的机器学习模型的实例包括:支持矢量机(SVM)、回归分析、线性回归、逻辑回归、朴素贝叶斯(Bayes)、线性判别分析、决策树、K近邻算法、随机森林、人工神经网络(ANN),如卷积神经网络(CNN)、循环神经网络(RNN)、全连接神经网络、长短期记忆(LSTM)模型以及其他;和/或其组合。在无监督学习框架内训练的模型会推断出确定数据集隐藏结构的函数,而不需要关于数据的预先知识。本领域已知的无监督机器学习模型的实例包括:聚类,如k均值聚类、混合模型聚类、层次聚类;异常检测方法;主成分分析(PCA)、独立成分分析(ICA)、T型分布的随机邻接嵌入(t-SNE);生成模型;和/或无监督的神经网络;自动编码器;和/或其组合。半监督学习(SSL)是机器学习框架,在这个框架内,人们可以使用标记的和未标记的数据来训练模型。数据增强方法可以选择用于从稀缺组的真实数据样本中产生人工数据样本,并增加用于模型训练的数据数量和多样性。无标签的数据当与少量的标签数据一起使用时与其他框架相比可以产生相当大的学习准确性的提高。当只有部分可用数据被标记时,这种方法特别有趣。
“卷积神经网络”或“CNN”是指一种机器学习模型,它使用多个数据处理层,称为卷积层,以最适合解决分类或回归任务的方式来表示输入数据。在训练过程中,使用本领域已知的优化算法,如反向传播算法来执行随机梯度下降,对每个CNN层优化权重参数。在运行时,所得到的训练后的CNN可以非常有效地处理输入数据,例如,在学习分类任务的情况下,将其分类为正确的数据输出标签,并尽可能减少假阳性和假阴性。卷积神经网络也可以与循环神经网络相结合,以产生深度学习分类器。
基因组分析系统
现在将参照图1进一步详细描述示例性的基因组分析系统和工作流程。对于DNA分析技术领域的技术人员来说显而易见的是,基因组分析工作流程包括:在实验室(也称为“湿式实验室”)进行的初步实验步骤,以产生DNA分析数据,如下一代测序工作流程中的原始测序读数;以及使用生物信息系统(也称为“干式实验室”)对DNA分析数据进行的后续数据处理步骤,以进一步确定终端用户感兴趣的信息,如详细识别DNA变体和相关注释。根据实际应用、实验室设置和生物信息学平台,DNA分析工作流程的各种实施方式是可能的。图1描述了NGS系统的实例,包括湿式实验室系统,其中DNA样本首先用DNA文库制备方案100进行实验制备,该方案可以产生、适应测序和扩增DNA片段,以方便NGS测序仪110处理。在下一代测序工作流程中,产生的DNA分析数据可以作为FASTQ格式的原始测序读数的数据文件产生。然后,该工作流程还可以包括干式实验室基因组数据分析系统120,其将根据提议的方法制备的DNA样本池的原始测序读数作为输入,并应用一系列数据处理步骤来表征输入样本的某些基因组特征。示例性的基因组数据分析系统120是索菲亚数据驱动医学平台(Sophia DDM),其在2020年已经被全球1000多家医院用于自动识别和表征基因组变体并报告给终端用户,但也可以使用其他系统。例如,在国际PCT专利申请WO2017/220508中描述了可由基因组数据分析系统120应用于基因组变体分析的数据处理步骤的不同详细可能的实施方式,但其他实施方式也是可能的。
如图1所示,基因组数据分析仪120可以包括序列比对模块121,其将原始NGS测序数据与参考基因组进行比较,例如医学应用中的人类基因组,或兽医应用中的动物基因组。在传统的基因组数据分析系统中,所得的比对数据可由变体调用模块(未示出)进一步过滤和分析,以检索变体信息,如SNP和INDEL多态性。变体调用模块可以被配置为执行不同的变体调用算法。然后,所检测到的变体信息可由基因组数据分析模块120输出为基因组变体报告,供终端用户进一步处理,例如用可视化工具,和/或由进一步的变体注释处理模块(未示出)进一步处理。在可能的实施方式中,基因组数据分析系统120可以包括自动数据处理模块,如覆盖度数据准备模块122以准备空间排列的覆盖度数据信号,和染色体空间不稳定性(CSI)分析模块123以分析空间排列的覆盖度数据信号并获得CSI信息,然后可以报告给终端用户,例如用可视化工具报告给终端用户,或报告给另一个下游过程(未示出)。在可能的实施方式中,染色体空间不稳定性(CSI)分析模块123可以适于分析空间排列的覆盖度数据信号,并从CSI信息中得出CSI分数,如HRD分数,或CSI状态,如HRD状态。然后,HRD分数、HRD状态和/或CSI信息可以报告给终端用户,例如用可视化工具报告给终端用户,或报告给另一个下游过程(未示出)。
数据处理工作流程
基因组数据分析仪120可以通过采用和结合不同的数据处理方法来处理测序数据以产生基因组数据分析报告。
序列比对模块121可以被配置为执行不同的比对算法。可以使用标准的原始数据比对算法,如为快速处理大量基因组数据测序读数而优化的Bowtie2或BWA,但也可以使用其他实施方式。比对结果可以表示为一个或几个BAM或SAM格式的文件,如生物信息学领域的技术人员所知,但也可以使用其他格式,例如压缩格式或为保序加密优化的格式,这取决于基因组数据分析仪120对存储优化和/或基因组数据隐私执行的要求。
图2示出了自动数据处理模块的可能工作流程,以从生物信息学渠道中的比对文件(如BAM文件)的分析中产生CSI报告。在可能的实施方式中,覆盖度数据准备模块122可以被调整为从比对的读数中产生200覆盖度数据信号,作为与每个基因组区块对齐的读数或对端读数的原始数量,并将所得数据210排列为空间排列的覆盖度数据信号。空间排列的覆盖度数据可表示为矩阵,即一组矩阵,或表示为覆盖度多维阵列数据(张量),以方便CSI分析仪123对二维图像,即多平面二维图像,或多维阵列(张量)输入信号的处理。然后,染色体空间不稳定性(CSI)分析仪模块123可以进一步分析、分类和/或归类220空间排列的覆盖度数据信息,以向终端用户报告230该数据是否表现出某些特性,这些特性是肿瘤样本的可能的染色体空间不稳定性(CSI)的特征。在可能的实施方式中,染色体空间不稳定性(CSI)分析模块123可以产生CSI分数。在另一个可能的实施方式中,染色体空间不稳定性(CSI)分析模块123可以产生CSI状态,如阴性状态、阳性状态或未确定状态。在可能的实施方式中,CSI状态可以根据CSI分数来推断,但其他实施方式也是可能的。在优选的实施方式中,CSI可以将肿瘤样本的同源修复缺陷(HRD)表征为HRD分数或HRD状态,但也可以表征与染色体空间不稳定性有关的其他基因组特征,例如由BRCA1或BRCA2缺陷引起的HRD。然后,这份报告可能有助于对癌症亚型和/或其预后进行肿瘤基因组学诊断。该报告还有助于选择专门针对所分析的肿瘤的医疗方法,例如使用某些癌症治疗方法,利用HRR途径的缺陷来减少或停止肿瘤细胞的增殖。
基因组数据分析仪120可以是计算机系统或计算机系统的一部分,包括中央处理单元(本文中的CPU,“处理器”或“计算机处理器”),存储器,如RAM和存储单元如硬盘,以及通过通信网络,例如互联网或本地网络与其他计算机系统进行通信的通信接口。基因组数据分析仪计算系统、环境和/或配置的实例包括但不限于:个人计算机系统、服务器计算机系统、瘦客户机、厚客户机、手持或笔记本电脑设备、多处理器系统、基于微处理器的系统、机顶盒、可编程消费电子产品、网络PC、微型计算机系统、大型计算机系统、图形处理单元(GPU)等。在一些实施方式中,计算机系统可包括一个或多个计算机服务器,其与众多其他通用或特殊用途的计算系统一起运行,并且可实现分布式计算,如在基因组数据农场中的云计算。在一些实施方式中,基因组数据分析仪120可以被整合到大规模的并行系统中。在一些实施方式中,基因组数据分析仪120可以直接集成到下一代测序系统中。
基因组数据分析仪120计算机系统可以在计算机系统可执行指令(如由计算机系统执行的程序模块)的一般情况下进行调整。一般而言,程序模块可以包括执行特定任务或实现特定抽象数据类型的例程、程序、对象、组件、逻辑、数据结构等。正如那些熟练掌握计算机编程技术的人所熟知的那样,程序模块可以使用本地操作系统和/或文件系统功能、独立的应用程序;浏览器或应用程序插件、小程序等;商业或开源库和/或库工具,因为可以用Python、Biopython、C/C++或其他编程语言进行编程;自定义脚本,如Perl或Bioperl脚本。
指令可在分布式云计算环境中执行,其中任务由通过通信网络连接的远程处理设备执行。在分布式云计算环境中,程序模块可以位于本地和远程计算机系统存储介质中,包括存储器存储设备。
因此,可以理解,这里描述的方法是计算机实现的方法。
覆盖度数据准备-区块化
在可能的实施方式中,首先将参考基因组坐标划分为一组P区块,将BAM或SAM文件(测序读数比对后)中的比对读数作为输入,使用本领域已知的生物信息学方法,计数映射到每个区块的DNA片段的数量,以产生200个沿P区块的覆盖度数据信号。对于生物信息学技术领域的技术人员来说显而易见的是,在一些实施方式中,一些区块可以被忽略、过滤或与其他区块合并以提高输入数据的质量。
在可能的简单实施方式中,沿参考基因组的区块可以使用均一尺寸,但也可以选择异质的尺寸。区块的尺寸可被选择成使其为每个染色体的覆盖度数据信号的分析保留足够的空间细节,同时便于用本文所述的机器学习分析方法进行覆盖度数据信号分析。优选的是,区块的尺寸可以为300bp至20Mpb,但其他实施方式也是可能的。
在参考人类基因组的情况下,染色体从最大的染色体(chr1)到最小的染色体(chr22)依次排序,可能排除了性染色体。人类染色体按长度递减进行编号,从249Mbp的chr1到51Mbp的chr22。例如,通过将参考基因组按3Mbp的区块划分,可以得到一维阵列(矢量),其中最长的1号染色体有84个区块,最短的22号染色体有4个区块。表1示出了理论上着丝粒区各自的起点和终点位置,染色体总长度及其类型,以及人类基因组中每个染色体的p臂(短臂)和q臂(长臂)各自的长度(Mbp),参照在用于读数比对的参考基因组的坐标系统中的GRCh37参考基因组(hs37d5版本)。请注意,在目前的高通量测序技术下,并不是所有的人类基因组区块都可以包含BAM文件中的比对读数,这是因为很难将短的NGS读数准确地映射到人类基因组的某些重复序列或有问题的区域,例如近端着丝粒染色体13、14、15、21、22或Y的p臂。
表1.在GRCh37参考基因组坐标系统(hs37d5版本)中表示的人类染色体的着丝粒区位置、染色体长度(bp)和染色体臂长度(Mbp)的开始和终止位置
对生物信息学领域的技术人员来说显而易见的是,在可能的实施方式中,沿全基因组测序数据的区块尺寸和位置也可以根据BAM文件的实际覆盖度进行调整。一般而言,我们可以认为P个区块的集合中的每个区块可以由对应于参考坐标系中的区块{区块_开始,区块_终止}1<=区块_指数<=P的起始和结束位置的2个值来定义。因此,区块的尺寸可以是可变的,即区块_尺寸=区块_终止-区块_开始。在可能的实施方式中,区块的尺寸可以为0.5Mbp至1.5Mbp,或从1.5Mbp至2.5Mbp,或从2.5Mbp至3.5Mbp,或从3.5Mbp至4.5Mbp,或从4.5Mbp至5.5Mbp,或从5.5Mbp至6.5Mbp,或从6.5Mbp至7.5Mbp,或从7.5Mbp至8.5Mbp,或从8.5Mbp至9.5Mbp。在另一个可能的实施方式中,可以选择区块的尺寸,以便为每个染色体或每个染色体臂获得固定数量的区块。其他实施方式也是可能的。
然后可以通过计数固定的区块集合的每个区块中的比对读数的数量来获得每个染色体或每个染色体臂的覆盖度信号。在可能的实施方式中,得到的每个染色体或每个染色体臂的高分辨率覆盖度数据信号阵列可作为下游CSI分析的高维度数据输入。在可能的替代实施方式中,可使用100kbp的第一尺寸来产生第一组高分辨率覆盖度区块。产生的高分辨率覆盖度数据信号可以进一步坍塌成1到20Mbp的第二组较大的区块,例如通过结合约30个100kbp的大初始区块来产生较大尺寸(2.5Mbp至3.5Mbp)的区块,从而降低下游CSI分析的数据的维度。为此可采用生物信息学领域的技术人员已知的各种方法,如平滑、使用中位数或平均值来降低维度、池化、抽样和其他方法。在优选的实施方式中,坍塌将每个区块唯一地分配给特定的染色体臂,而不将来自不同染色体臂的覆盖度数据混入区块。然后,覆盖度数据从一个染色体臂清晰地分离到下一个染色体臂,这可能有助于在染色体臂水平上分析染色体空间不稳定性事件。
在可能的实施方式中,坍塌可以包括根据每个染色体臂的长度来调整坍塌区块的尺寸,以减少染色体臂边界处的边界效应。例如,在染色体臂的长度不能被所考虑的初始区块尺寸所分割的情况下,最后得到的坍塌区块通常会比其他的短,造成边界效应(例如,染色体臂任何一端的着丝粒区块或端粒区块可能不包括整个3Mbp的区域,并可能导致这些区域的扭曲的覆盖度,而可能是分析CSI事件的关键)。减少这种影响的策略是调整相对于染色体臂长度的区块化尺寸,使所考虑的染色体臂长度可以被区块化尺寸分割,具有可忽略的剩余部分。这就保证了在所考虑的染色体臂上的所有坍塌区块都有相同的实际尺寸。
在可能的进一步实施方式中,坍塌可以包括沿染色体臂单独调整坍塌区块的开始或结束位置,从而调整实际尺寸,使坍塌区块的可变尺寸之和仍然等于染色体臂的长度,同时沿染色体臂或多或少地均匀分布。
图3示意性地示出了染色体300,包括对应于着丝粒的中心位置320、对应于p臂端粒的起始位置310和对应于q臂端粒的末端位置330。该染色体包括从起始位置310开始的第一端粒区315,围绕着丝粒的着丝粒区325,以及结束于末端位置330的第二端粒区335。对于生物信息学领域的技术人员来说显而易见的是,着丝粒区域325的一部分351在覆盖这个区域325的现有参考基因组中没有相关的参考序列,因此不能被映射。在可能的实施方式中,对应于有限覆盖度区域或低质量区域的高分辨率区块也可以根据相应区域的基因组知识被移除(例如图3中的区域370)。在优选的实施方式中,可以从BAM文件中为每个染色体臂单独构建一组不同的第一高分辨率覆盖度区块350、352。然后,每组高分辨率的区块可以进一步坍塌成第二组更大范围的区块360、362,用于每个臂。区块也可以根据相应区域的基因组知识跨越更短或更大的基因组区域(例如,2.5Mbp的区块371,然后是3.5Mbp的区块372,平均到3Mbp的区块),和/或优化两个染色体臂各自末端的区块(端粒区块372,375,这是靠近端粒310、330的最近的区块;或着丝粒区块373,374,这是靠近着丝粒320的最近的区块)。
例如,来自每个染色体臂的初始100kbp的区块可以被独立平滑,以获得3Mbp的区块(即30个100kbp的区块的聚集)。由于每条染色体臂的区块总数不可能被30整除,染色体臂的最后平滑区块通常包含少于30个原始区块。如果保留这些区域,这可能会给这些区域带来太多重要性,反之,如果删除这些区域,则会损失太多信息。为了减少边界效应,可以为每个染色体臂独立地坍塌较大的可变尺寸的区块。在可能的实施方式中,可以选择导致最小的非完整区块余量和最接近目标尺寸的区块尺寸。在可能的实施方式中,最后一个区块如果不完整,可以拒绝(因为现在可以忽略不计),但其他实施方式也是可能的。表2第3列示出了每个染色体臂的提议的坍塌区块尺寸的示例性分布,最初的一组100kbp的区块沿2*22染色体臂坍塌成目标3Mbp尺寸的较大区块。第4列示出了所得的被拒绝的边界区块的尺寸,与第2列中默认的剩余高分辨率区块的数量相比较,当所有臂的固定尺寸为3Mpb时,将形成臂边界的最后一个坍塌区块。
表2.示例性的区块坍塌分布,从100kbp的区块到3MBp的目标坍塌区块尺寸。
覆盖度数据准备-归一化
对于生物信息学领域的技术人员来说显而易见的是,在一些实施方式中,原始覆盖度数据信号可以通过本领域已知的方法进一步归一化、过滤或平滑化。在可能的实施方式中,归一化可以包括每个样本的归一化,通过将覆盖度数据信号除以整个样本的平均覆盖度信号,以摆脱与样本测序实验有关的可能导致每个样本覆盖度差异的偏差。归一化也可以包括通过GC含量的归一化来应用GC-偏移校正。其他的实施方式也是可能的,如逐个区块归一化,线性和非线性最小二乘回归,GC LOESS(GC归一化),LOWESS,PERUN(每个样本的归一化)、RM、GCRM、cQn和/或其组合。
在可能的进一步的实施方式中,归一化的数据可以进一步分割成区块内的子区域,以确定覆盖度数据信号中的同质区域,可能对应于不同的拷贝数(通常每个染色体2个)。覆盖度分段法可以从给定样本的覆盖度曲线推断出相对于其他基因组区域的离散片段的拷贝数。为此,覆盖度曲线被分解为段(即基因组的一部分),其中的覆盖度被认为是恒定的。然后,在上一步中定义的每个片段都可以与以覆盖度信号振幅定义的离散水平相联系。在这个阶段,拷贝数变异(CNV)事件已经可以被检测到,但与每个覆盖度水平相关的绝对拷贝数仍然是未知的。
在样本肿瘤内容和倍性已知或可以从数据中推断的情况下,通过分割确定的离散水平可以映射到反映肿瘤中存在的拷贝数的整数(例如,CN=1,CN=2,……)。对生物信息学技术领域的技术人员来说显而易见的是,无论是否有绝对的拷贝数参考,这种分割方法都可以减少噪音,并便于识别样本的拷贝数概况。然而,在存在噪声的情况下,分割可能会产生错误的结果,因为它抑制了由拷贝数变化导致的覆盖度波动,特别是在低覆盖深度测序时。因此,在优选的实施方式中,本文所述的基因组分析方法可以在没有分割步骤的情况下运行。在另一个实施方式中,在应用本文公开的机器学习方法之前,可以应用分割步骤对数据进行预处理。
在可能的归一化和/或分割步骤之后,覆盖度数据通常由一组一维矢量组成,其中每个矢量对应染色体。在可能的实施方式中,覆盖度数据也可以由一组一维矢量组成,每个矢量对应于染色体臂。在可能的实施方式中,矢量的尺寸可以是可变的,反映了不同的染色体有不同的尺寸。在另一个可能的实施方式中,区块的尺寸可变,以便每个染色体矢量(或每个染色体臂矢量)具有相同的尺寸。
图4描绘了全基因组NGS测序实验的示例性覆盖度数据1D信号,该实验以在人类基因组前22条染色体(不包括性染色体)上的1X覆盖度,100kbp的区块,按样本和GC含量归一化后进行。图5示出了相同的覆盖度数据信号在每条染色体上的排列,坍塌后的区块为3Mbp。
覆盖度数据的空间排列
在优选的实施方式中,覆盖度数据因此可以进一步在空间上排列210,以形成多维数据阵列,其中覆盖度数据被组织成沿着阵列结构中的任何它们的端粒或着丝粒的区块对齐每个染色体或染色体臂。多维阵列结构可能特别适合受益于最近为非基因组应用(如图像分类器)开发的高效基因组分析机器学习训练模型,而不是像现有技术中的HRD检测方法那样依赖简单的人类工程化的特征。
换句话说,并与所提供的端粒或着丝粒区块的定义相一致,覆盖度数据信号矢量可进一步在空间上排列210,以形成覆盖度数据信号图像,其中覆盖度数据被组织成沿至少一个与着丝粒区域相邻的着丝粒区块或至少一个与二维阵列结构中的端粒区域相邻的端粒区块对齐每个染色体或染色体臂的区块。
在可能的实施方式中,覆盖度数据信号矢量可进一步在空间上排列210,以形成一维覆盖度数据阵列,其中覆盖度数据被组织为沿基因组的每个染色体或染色体臂的区块的尾部到头部的连接。
在另一个优选的实施方式中,覆盖度数据信号矢量可进一步在空间上排列210以形成多维覆盖度数据信号阵列,其中覆盖度数据被组织成沿阵列结构中与着丝粒区域相邻的至少两个区块对齐每个染色体或染色体臂的区块,其中一个区块在着丝粒的左边,另一个区块在着丝粒的右边,并且该多维数据阵列是二维阵列,具体是二维图像。
在另一个优选的实施方式中,覆盖度数据可以进一步在空间上排列210以形成一个多维数据阵列,其中覆盖度数据被组织成将每个染色体或染色体臂的区块沿着阵列结构中邻近端粒区域的至少两个区块对齐,其中,一个区块位于与该染色体相关的覆盖区块集内的染色体起点,另一个区块位于与该染色体相关的覆盖区块集内的染色体终点,并且该多维数据阵列是二维阵列,特别是二维图像。
在另一个优选的实施方式中,覆盖度数据可进一步在空间上排列210以形成多维数据阵列,其中覆盖度数据被组织成沿着阵列结构中与着丝粒区域相邻的至少一个或至少两个区块或与端粒区域相邻的至少一个或至少两个区块比对每个染色体或染色体臂的区块,其中,这个多维数据阵列是二维阵列,具体是二维图像。
在这种情况下,与着丝粒或端粒区域相邻的区块意味着它是最接近这些区域的区块。
在进一步的优选实施方式中,每个染色体臂的区块是独立计算的。
在一个实施方式中,方法可以包括但不限于这样步骤,其中对于每个染色体臂的覆盖度数据可以通过在2.5Mbp至3.5Mbp的窗口上取中位归一化覆盖度进一步归一化。对于每个染色体臂来说,确切的区块分辨率可以在这两个值之间独立选择,从而在染色体臂的末端得到最小的可能的非完整区块。最后不完整的区块被丢弃在机器学习模型的输入中,如CNN。因此,每个p臂的最后一个区块包含的覆盖度信号在其位置上最接近其相应的染色体着丝粒。同样,每个q臂的第一个区块包含的覆盖度信号是最接近其相应染色体着丝粒的位置。同时,每个p臂的第一个区块包含覆盖度信号就其位置而言是最接近其相应的p臂染色体端粒区域。同样,每个q臂的最后一个区块包含的覆盖度信号是最接近其相应q臂染色体端粒区域的位置。
因此,在图3中可以理解和描述,着丝粒区域可能不包括在这些区块中。因此,p臂的最后一个区块和q臂的第一个区块反而标志着一个染色体的着丝粒区域。
组织覆盖度数据以比对本文所述的区块的效果允许检测检测到的事件之间的空间关系,如染色体臂内的基因组不稳定性。机器学习模型,如CNN,可以使用检测到的事件的空间分布信息,因此当准确的状态检测得益于这种类型的空间排列时,有可能提高其输出预测。
图6示出了根据着丝粒位置排列的染色体的图形表示。从表1和图6可以看出,最长的染色体chr1的长度是249Mbp,最长的p臂的长度对于chr1是124.7Mbp,而最长的q臂在chr2中实际上更长,为147.9Mbp。在可能的实施方式中,在给定的恒定的区块_尺寸中,覆盖度数据可以被适应于(上限(chr_p-臂_max_len)/区块_尺寸+上限(chr_q-臂_max_len)/区块_尺寸)*N的二维数组,其中N是表示染色体的数量,chr_p-臂_max_len是最长的p臂长度,chr_q-臂_max_len是要联合分析的N个染色体集合中最长的q臂长度。作为说明性的实例,对于22条非性染色体的臂上的全基因组测序覆盖度数据按其着丝粒区块排列,以区块_尺寸=3Mbp的恒定尺寸的区块,尺寸为92*22(或22*92)的阵列可以由基因组数据分析仪120在空间上排列210作为CSI分析仪123的数据输入。最长的p臂124.7Mbp(chr1)的覆盖度数据可能确实适合于42个3Mbp的区块,而最长的q臂147.9Mbp(chr2)可能适合于50个3Mbp的区块。在可能的实施方式中,多维阵列可以被排列为代表2*N个染色体臂在待分析的N个染色体组中最长的染色体臂的区块。例如,覆盖度数据可以在空间上排列成尺寸为44*50的阵列,用于分析22条染色体,其恒定的区块尺寸为3Mbp。每个染色体臂的覆盖度数据可以在阵列中排列成行(或列),以便染色体臂的着丝粒区块在阵列中保持一致。在可能的实施方式中,可以应用图3上示意性的自适应区块来去除不相关的高分辨率区块,并调整沿染色体基因组的可变区块尺寸和位置。例如,尺寸为84*22(或22*84)的阵列可以由基因组数据分析仪120在空间上排列210,作为CSI分析仪123的数据输入。其他实施方式也是可能的。
由于区块的数量随着染色体和/或染色体臂的长度而变化,阵列的一些元素需要通过填充覆盖度数据中没有的虚拟数据来填补。对于数据科学领域的技术人员来说,可以有不同的选择来填充多维阵列的空条目(或元素)。在可能的实施方式中,空元素可以由常量值来填充。在另一个可能的实施方式中,空元素可以从预定义的常量值的掩码阵列中填充。在另一个可能的实施方式中,空元素可以根据行或列中被真实数据填充的部分的覆盖度数据来填充。在可能的实施方式中,空的元素可以通过重复覆盖度值,从沿着染色体或染色体臂行(或列)的区块的列(或行)的最后一个可用区块向左或向右(分别为高于或低于)进行填充。
在另一个可能的实施方式中,可以选择区块的单个尺寸,以促进阵列的密集填充。在具有特殊基因组意义的区域可以选择较小的区块尺寸(对应于更好的分辨率),以增加某些染色体臂的区块的总数。
图7示出了归一化的覆盖度数据在空间上重新排列为适合用机器学习模型进行CSI分析的二维阵列图像的实例,分别来自a)HRD阴性肿瘤DNA样本和b)HRD阳性肿瘤DNA样本。染色体按行排列,覆盖区块沿列排列,使得在染色体空间不稳定性分析之前,根据所提出的覆盖度数据制备方法的某些实施方式,着丝粒区块在单列中竖直排列。像素越亮,绘制的区块的归一化覆盖度计数值C区块越高,反之,像素越暗,该数值越低。
染色体空间不稳定性(CSI)分析仪
基因组数据分析系统250可进一步将染色体空间不稳定性分析仪123应用于空间排列的覆盖度数据多维阵列,以自动分析和报告表征每个DNA样本中两个或更多个染色体或染色体臂的某些空间不稳定性模式的一种或多种基因组特性。在肿瘤样本分析的情况下,CSI分析仪模块123可以报告染色体畸变的负担,如沿每条染色体臂的特定基因组区域的等位基因的大量缺失或复制。
在可能的实施方式中,CSI分析仪模块123是计算机实施的算法,它被训练来区分具有基因组不稳定性(被标记为HRD阳性状态)的DNA样本和没有基因组不稳定性(被标记为HRD阴性状态)的另一个DNA样本。
在可能的实施方式中,在肿瘤样本的情况下,CSI分析仪模块123可以识别特征(signature),该特征可以识别具有类似特征的样本,从而可能对给定的治疗有良好的应答。
CSI分析仪模块123可相应地识别、分类和报告指标、标量分数或特征指标的组合,以描述肿瘤样本的一个或多个基因组特性。在优选的实施方式中,CSI分析仪模块123可以相应地识别、分类和报告肿瘤样本的状态,如HRD-阴性或HRD-阳性,或者可能是未确定(不确定)的。在可能的实施方式中,CSI分析模块123可以报告标量分数作为肿瘤样本的HRD可能性的指标。对于本领域技术人员来说显而易见的是,这种生物信息学方法随后大大促进了对肿瘤样本中癌细胞基因组改变的详细了解,并使个性化医学治疗适应患者的推断出的癌细胞生物学的具体情况。
在一个实施方式中,基因组特性报告可以在图形用户界面上显示给终端用户。在另一个可能的实施方式中,基因组特性报告可以被制作成文本文件,以便进一步自动处理。其他实施方式也是可能的。
在优选的实施方式中,CSI分析仪模块123可以被调整为将机器学习模型应用于空间排列的覆盖度数据多维阵列。在可能的实施方式中,CSI分析仪模块123可以包括经训练的神经网络或经训练的神经网络的组合作为数据处理模块来分析220作为经训练的神经网络输入的空间排列的覆盖度数据多维阵列,并产生作为经训练的神经网络输出的CSI结果。
现在将参考图8进一步详细描述示例性的CSI分析仪123及其数据处理工作流程。对于机器学习领域的技术人员来说,显然可以使用本领域已知的不同方法和架构来处理准备的覆盖度数据多维阵列输入,并得出肿瘤样本的CSI状态或CSI分数。
示例性的CNN机器学习模型
为了对多维阵列(如二维图像)进行分类,其中染色体臂已沿其各自的着丝粒位置排列,以突出表征该染色体组的CSI状态的视觉模式,可通过采用图像模式识别领域中已知的方法和技术,如TensorFlow、Caffe、Caffe2、Pytorch、Theano或其他,来构建和训练卷积神经网络(CNN)。
图8示意性地示出了卷积神经网络(CNN)的实例,它可以处理患者肿瘤DNA样本的空间排列覆盖度数据矩阵图像,以将其分类为HRD阳性(HRD+)或HRD阴性(HRD-)。在图8的实例中,经训练的CNN将空间排列的覆盖度数据矩阵图像800作为输入,该图像以22条非性染色体的行排列,从顶部的chr1到底部的chr22,并且以区块的列排列(在图8的示例性图示中,约3Mbp尺寸,沿染色体的着丝粒区域排列,基本填充为0覆盖度数,但其他实施方式也是可能的)。输入的图像800被输入到有几个处理层的第一卷积神经网络810,以量化一组中间特征,这些特征本身被输入到第二神经网络分类器820,负责提取两个分数(一个为HRD-,一个为HRD+)作为其两个输出。然后,CSI分析仪可以通过比较和阈值化HRD-和HRD+的输出值来得出HRD状态(HRD-,HRD+,或不确定)。
更一般而言,对于本领域技术人员而言显而易见的是,CSI分析仪123可以被调整为采用各种CNN架构。卷积网络可以包括第一系列串行连接的卷积数据处理层810,它可以进一步馈送第二系列全连接层720以输出各种预测结果。例如,CNN结构可以包括一个或多个卷积层,布置为用于处理一维、二维或多维数据;可选地,一个或多个最大、中值或平均集合层,布置为用于处理一维、二维或多维数据;一个或多个多重滤波层,以便在训练期间进行正则化;可选地,一个或多个批量正则化层;扁平化层,接着是一个或多个全连接层;或其组合。可以使用本领域已知的不同激活方法,例如ReLU、SELU,或更普遍的防止梯度消失的任何激活函数。不同的输出激活也可以使用本领域已知的方法,例如,sigmoid函数方法,以产生与单标签二元分类相关的概率(例如,正面状态的概率,负面状态的概率=1-正面状态的概率);或用softmax方法产生与单标签多类分类输出相关的概率(例如,阳性状态、阴性状态或不确定状态的互斥概率),或用回归方法,如线性方法或ReLU分类方法产生标量值,该标量值不打算被解释为概率,例如HRD分数。
在优选的实施方式中,CSI分析仪123通过对CNN的输出进行阈值处理,产生最终的输出决定,如每个样本的HRD+或HRD-状态分析。阈值可以根据应用进行调整,以便根据终端用户的需求优化CSI分析的灵敏度和/或特异性,例如用于诊断、预后或指导选择最有效的癌症治疗。在可能的实施方式中,CSI分析仪123还可以根据终端用户的需要进一步报告230CNN的输出值。
图9示出了一个可能的详细CNN架构的实例,其可由CSI分析仪123采用,用于分析220空间排列的覆盖度数据。对于网络内部的每个卷积层,过滤器的数量以及过滤器的尺寸被表示为[N x H x W],其中N是过滤器的数量,H是过滤器的高度,W是过滤器的宽度。每个卷积层都通过激活函数(正则化线性单元,ReLU),然后再进行批量归一化(BN)。中间卷积层可以与最大池化层交错,然后按照通常的CNN结构,由平均池化层和扁平化层进行交错。在这个实例中,CNN的输出层被扁平化,通过批量归一化层,Dropout层,并使用sigmoid激活函数给单节点密集层(即全连接层)作为输入。
染色体空间不稳定性(CSI)分析仪训练
对于有监督和半监督的CSI分析仪训练,可以使用一组标记的输入-输出对作为训练集。例如,为了训练CSI分析仪以对HRD状态进行分类,可以使用公共领域HRD标记的数据样本子集,使网络能够通过优化过程(如反向传播)沿其不同层调整其参数(如权重),使其输出误差最小化,该优化过程在训练阶段被特别应用。对于机器学习领域的技术人员来说,显然可以使用不同的损失函数来测量误差,例如,对于具有单标签二元分类输出的模型,可以使用二进制交叉熵测量,对于具有单标签多类别分类输出的模型,可以使用分类交叉熵测量,或基于均方根误差、绝对误差或其他误差测量的回归方法。优化可以采用现有的基于随机下降梯度的优化器,如ADAM或RMPSprop。规则化可以采用放弃法、早期停止法和/或L-1或L-2参数规则化。其他实施方式也是可能的。
随着更多与染色体不稳定性有关的临床病症被当前和未来的肿瘤学研究分类,额外的数据集可用于训练机器学习模型,旨在检测由HRD或其他现象引起的基因组疤痕,如对特定治疗的应答。在可能的实施方式中,HRD+和HRD-的二元状态,以及可选择的“不确定”状态可作为多类基础真实标签使用。在另一个可能的实施方式中,来自现有技术方法(如HRDetect)的HRD分数可作为标量输出标签。根据本文所公开的方法的任何可能的实施方式,只要在训练时和运行时对多维阵列格式实施方式的选择保持不变,输入的覆盖度信号数据阵列也可以在空间上排列为多维数据。
在可能的实施方式中,用于训练机器学习模型的额外数据集可以包括来自患者DNA样本的覆盖度数据信号阵列,其中,这些患者是已知的同源重组缺陷状态,和/或其中这些患者接受了癌症治疗方案,并且治疗结果是已知的(应答与无应答)。
在可能的实施方式中,用于训练机器学习模型的额外数据集可以包括来自原代或永生化细胞系实验或肿瘤衍生的类器官实验产生的DNA样本的覆盖度数据信号阵列,特别是确定预测有HRD的类器官是否被癌症治疗方案消除或其生长减少的实验。
在可能的实施方式中,可以采用半监督学习框架来训练机器学习模型。事实上,获取标记的数据可能具有挑战性,而且成本很高,特别是在HRD状态预测的背景下,这使得半监督学习成为这一应用的有希望的框架。由于模型的过拟合程度由其复杂性和接受的训练量决定,为模型提供更多的训练实例可以帮助减少过拟合。数据增强包括从现有数据中创建新的人工训练样本。在可能的实施方式中,人工样本可以通过从共享相同HRD状态并具有类似肿瘤内容(或仿生预处理以模仿类似肿瘤内容)和倍性的现有真实样本中抽取染色体来生成。
在可能的实施方式中,可以通过从共享相同HRD状态和相同归一化纯度和倍性比的可用真实样本中取样染色体来生成人工样本,其中,纯度是样本中肿瘤细胞的百分比,倍性是整个基因组的平均染色体拷贝数,从而表征肿瘤细胞平均有多少个完整的染色体组。
在一个实施方式中,机器学习模型是在来自标记的真实样本的覆盖度数据信号阵列数据上训练的。
在另一个实施方式中,机器学习模型是在覆盖度数据信号阵列数据上训练的,这些数据来自于标记的真实样本和标记的数据增强样本,这些样本是通过从具有相同HRD状态的可用真实样本中取样染色体并确保所取样的染色体的纯度/倍性比率相同而产生的。
因此,任何机器学习模型都可以学习将HRD状态为阳性的DNA样本和HRD状态为阴性的另一个DNA样本之间的差异量化为CSI分数,然后根据CSI分数结果对测试的DNA样本进行分类。
在可能的实施方式中,机器学习模型在监督或半监督模式下被训练,至少使用一组根据HRDetect方法标记有HRD状态的真实样本。优选的是,学习模型在监督模式下使用至少一组根据HRDetect方法标记有HRD状态的真实样本进行训练。
在可能的实施方式中,机器学习模型是在监督或半监督模式下训练的,使用的是通过对一组具有相同HRD状态的真实样本的染色体进行抽样并考虑纯度倍性比的差异而产生的人工样本。在可能的实施方式中,机器学习模型是在监督模式下训练的,使用的是通过对一组共享相同HRD状态的真实样本的染色体取样产生的人工样本,并考虑到纯度倍性比的差异。
在一个实施方式中,经训练的机器学习模型计算出受试者DNA样本的CSI分数,并使用该CSI分数来分类受试者DNA样本属于具有某种CSI状态的样本组的概率。
在一个实施方式中,经训练的机器学习模型计算出受试者DNA样本的CSI分数,并使用该CSI分数来确定受试者DNA样本属于一组样本的HRD状态。
在一个实施方式中,提供了一种确定DNA样本的CSI状态的基于计算机的方法,其中,CSI状态是HRD状态。
在一个实施方式中,提供了一种确定DNA样本的CSI状态的基于计算机的方法,其中,DNA样本是来自癌症样本的DNA。
在一个实施方式中,提供了确定DNA样本CSI状态的基于计算机的方法,其中,低通全基因组测序覆盖度至少为0.1X,且最高为30X,如1X至10X,如0.1X至5X,或如0.1X至1X。
在一个实施方式中,提供的是确定DNA样本的CSI状态的基于计算机的方法,其中,低通量基因组测序至少为0.1X且至多30X,例如1X至10X,例如0.1X至5X或例如0.1X至1X,获得的是一组包含来自至少2个和至多22个染色体的读数的受试者DNA样本的测序读数。
在一个实施方式中,提供了确定DNA样本的CSI状态的基于计算机的方法,其中,执行的是DNA样本的测序读数与人类参考基因组的比对。
在一个实施方式中,提供了确定DNA样本的CSI状态的基于计算机的方法,其中,执行的是通过GC含量进行归一化以应用GC-偏移校正。
在一个实施方式中,提供了确定DNA样本的CSI状态的基于计算机的方法,所述CSI状态为HRD状态。
在一个实施方式中,提供了确定DNA样本的CSI状态的基于计算机的方法,其中,机器学习模型通过使用已知基因组不稳定状态如HRD+/HRD-状态的样本的训练数据集进行预先训练,因此该机器学习模型被训练为区分具有某种CSI状态如HRD状态的样本的覆盖度谱图特性的样本。
诊断/预后方法
在一个实施方式中,提供了一种表征从患者获得的DNA样本的体外方法,该方法包括。
-从患者样本中分离出DNA片段;
-构建包含与一组染色体重叠的所述片段的测序文库;
-对该文库进行测序,使之达到至多30X的测序覆盖度;
-将得到的测序读数与人类参考基因组进行比对;
-并用本发明的方法生成受试者样本的CSI分数。
在一个实施方式中,提供了一种表征从患者获得的DNA样本的体外方法,该方法包括。
-从患者样本中分离出DNA片段;
-构建包含与一组染色体重叠的所述片段的测序文库;
-对该文库进行测序,使之达到至多30X的测序覆盖度;
-将得到的测序读数与人类参考基因组进行比对;
-并根据用本发明的方法得到的CSI分数,确定患者样本的HRD状态。
一种选择癌症患者用PARP抑制剂治疗的方法,包括用本发明所述的任何方法检测肿瘤患者样本的HRD阳性的步骤。
在一个实施方式中,患者患有癌症,可通过使用本文所述的任何方法将其分类为HRD阳性(HRD+)或HRD阴性(HRD-)。
在一个实施方式中,提供了一种表征从患者获得的DNA样本的体外方法,其中,患者样本是肿瘤样本,患者样本的染色体空间不稳定性指标是肿瘤对癌症治疗方案的应答的预测因子,该治疗方案包括铂类化疗剂、DNA损伤剂、蒽环类药物、拓扑异构酶I抑制剂,或PARP抑制剂。
在一个实施方式中,提供了一种表征从患者获得的DNA样本的体外方法,包括基于染色体空间不稳定性指标识别肿瘤患者样本是否同源重组(HR)缺陷的步骤,其中,是同源重组(HR)缺陷的可能性大,表明PARP抑制剂可用于治疗癌症的方法。PARP抑制剂可单独使用或与其他治疗方法联合使用。
在一个实施方式中,提供了一种选择癌症患者接受铂类化疗剂、DNA损伤剂、蒽环类药物、拓扑异构酶I抑制剂或PARP抑制剂治疗的方法,包括根据本文所述的确定DNA癌症样本CSI状态的基于计算机的方法,检测患者DNA癌症样本为HRD阳性的步骤。在一个实施方式中,用于选择癌症受试者进行治疗的方法是体外方法。
在一个实施方式中,提供的是确定本文所述的DNA癌症样本的CSI状态的基于计算机的方法,其中,该方法进一步包括:
当确定该DNA癌症样本是同源重组缺陷时,通过向测试受试者施用聚ADP核糖聚合酶(PARP)抑制剂治疗癌症。
根据一个方面,样本是患者的样本,其形式为组织、新鲜冷冻组织(FFT)、血液或任何体液,或细胞学标本/制备物(FFPE、涂片)等。根据一个方面,样本是患者的肿瘤样本,包括FFPE样本。
在另一个特定的实施方式中,患者患有癌症,尤其是高级别浆液性卵巢癌、前列腺癌、乳腺癌、胰腺癌等。
根据一个实施方式,提供了一种根据本文所述方法诊断癌症类型为HDR+或HDR-的方法,其中,从患有或怀疑患有癌症的受试者那里获得样本。
湿式实验室工作流程中的整合
-CSI分析仪的数据生成-WGS工作流程
在可能的实施方式中,提供了一种确定受试者DNA样本的CSI状态的方法,该方法包括:
a)提供包括核酸的样本材料;
b)准备用于全基因组测序的第一核酸测序文库;
c)对步骤b)中获得的第一核酸测序文库进行测序;
d)分析步骤c)中获得的核酸序列,其中,
根据本文所述的确定受试者DNA样本的CSI状态的基于计算机的方法,分析(在步骤c中获得的)核酸序列的序列。
-CSI分析仪的数据生成-WGS工作流程和变体调用程序目标富集工作流程
在可能的实施方式中,提供了一种确定受试者DNA样本的CSI状态的方法,该方法包括:
a)提供包括核酸的样本材料;
b)准备用于全基因组测序的第一核酸测序文库;
c)对步骤b)中获得的核酸测序文库进行测序;
d)分析步骤c)中获得的核酸序列,其中,
根据本文所述的确定受试者DNA样本的CSI状态的基于计算机的方法,分析(在步骤c中获得的)核酸序列的序列。
其中,该方法还包括:
e)从同一样本材料中制备第二核酸测序文库;
f)对第二核酸测序文库进行靶向富集;
g)对步骤f)中获得的靶向富集文库进行测序;
h)分析步骤g)中获得的核酸序列,其中,
根据任何已知的变体调用方法分析来自富集的核酸文库(在步骤f中获得)的序列。
对于基因组分析领域的技术人员来说显而易见的是,第一和第二文库可以在同一时间,或在不同时间,以任何顺序从同一样本中独立制备。在可能的实施方式中,第一和第二文库可以通过将一个文库分成两个子集获得,一个子集(第一文库-用于CSI分析的WGS文库)经过步骤c)和d),而另一个子集(第二文库-用于变体调用分析的靶向富集库)经过子集e)至h)。
在另一个的实施方式中,提供了一种确定受试者DNA样本的CSI状态的方法,该方法包括:
a)提供包括核酸的样本材料;
b)准备第一核酸测序文库;
c)对核酸测序文库进行靶向富集,获得第二富集的核酸测序文库;
d)对步骤b)和步骤c)中获得的核酸测序文库进行测序;
e)分析步骤d)中获得的核酸序列,其中,
在排除与c)中富集的区域比对的序列后,根据本文所述的确定受试者DNA样本CSI状态的基于计算机的方法,分析来自第一个非富集核酸序列文库(在步骤b中获得)的非靶向基因组区域的序列
和
根据任何已知的变体调用方法分析来自第二富集的核酸文库(在步骤c中获得)的序列。
对基因组学领域的技术人员来说显而易见的是,上述工作流程可以使用不同的基因组分析步骤。
提供包括核酸的样本材料的步骤可以根据任何已知的DNA提取方法进行。如果多个样本在一个实验中一起测序,可以采用基因组条形码化方法来识别下游分析中的样本材料。
制备核酸测序文库的步骤可以按照任何已知的制备测序文库的方法进行。可选地,制备核酸测序文库的步骤b)可包括进一步的核酸序列扩增步骤b.0.)。扩增核酸可以按照任何已知的方法进行,如聚合酶链式反应(PCR)。
进行靶向富集的步骤可以根据任何已知的靶向富集DNA-测序或RNA-测序的方法进行,如靶向杂交捕获(即基于杂交捕获的靶向测序)或基于扩增子的方法(扩增子测序)。
根据一个实施方式,靶向富集是基于捕获的靶向富集。根据该实施方式,进行靶向富集的步骤包括将至少一个探针或探针组与来自已知可能携带感兴趣变体的基因组区域(例如但不限于BRCA1和BRCA2区域)的靶向核酸杂交(“靶向区域”或“富集区域”),洗掉非靶向核酸并富集核酸。可以理解的是,术语“探针”或可互换的“诱饵”或“(探针)核酸分子”或“捕获探针”或“(DNA/RNA)寡核苷酸(捕获)探针”是指可以与靶核酸分子杂交的核酸分子。任何已知的探针设计都可以使用。术语“靶核酸”是指可被所用探针捕获的基因或转录物中的核酸区域。优选的是靶核酸是与HRR途径相关的基因,如选自但不限于BRCA1和BRCA2。
可以理解的是,在步骤b)中获得的核酸测序文库可以分为至少两部分,其中,一部分将不进行靶向富集,因此仍为非富集核酸文库;而另一部分将进行靶向,并形成富集核酸文库。
可以理解的是,非富集的核酸文库和富集的核酸文库可以一起或分别负载到测序仪。因此,根据所选择的工作流程,可以单独或一起获得核酸序列。
当非富集的核酸文库和富集的核酸文库一起被负载到测序仪中时,这些文库以预先确定的浓度负载。这些相对浓度是靶向(富集)和非靶向区域所需覆盖度的函数,使得非富集文库的核酸覆盖度为0.1X至10X,0.5X至10X,至少0.2X,或优选至少0.1X且至多30X;而富集文库的核酸覆盖度至少为30X,或30X至100X,或100X至500X,或500X至1000X,或1000X至5000X,优选至少4000X。
核酸文库的测序步骤可以根据任何已知的方法和使用已知的测序仪来进行。
在可能的实施方式中,来自非富集核酸文库的序列分析和来自富集核酸样本的序列分析是分别进行的。
当两个文库一起测序时,根据本文所述的基于计算机的方法对来自非富集核酸文库的序列进行分析以获得CSI状态,可以包括进一步过滤出靶向富集文库中覆盖的区域的步骤。这些区域可由CSI分析仪根据任何已知的屏蔽策略进行屏蔽或过滤,例如:
m1)仅通过将测序读数与参考基因组在基因组区域进行比对,而这些基因组区域未被靶向富集文库覆盖;或
m2)通过在比对文件(例如BAM或SAM文件)中只选择未被靶向富集文库覆盖的比对基因组区域;或
m3)通过从CSI分数测定中排除高分辨率区块,即与靶向基因组区域重叠的区块(例如,通过屏蔽与靶向富集文库覆盖的基因组区域重叠的区块)。
在非富集文库和富集文库的联合测序中,屏蔽从富集文库中获得的基因组的靶向富集区域的步骤确保靶向富集引入的覆盖度差异不会影响CSI分析。
在可能的实施方式中,对来自富集核酸文库的序列的分析包括根据已知的方法进行变体调用。从该分析中获得的有关变体等位基因分数(VAF)的数据是额外信息,可用于进一步表征患者DNA样本。在可能的实施方式中,样本的HRD状态可以从CSI状态、变体调用结果(例如,如果与同源修复缺陷有关的基因中的一个或多个变体被鉴定出来,例如,但不限于BRCA1或BRCA2基因中的某些变体),或从其任何组合(例如,如果CSI状态被鉴定为阳性或如果与同源修复缺陷有关的基因中的一个或多个变体被鉴定)得出。在可能的实施方式中,可以先进行变体调用分析,只有在变体调用结果没有鉴定出样本的HRD状态时,才根据本文所述的基于计算机的方法对样本进行CSI分析,以便从沿样本的全基因组的基因组不稳定性中鉴定可能的HRD状态。
实验-实施例1
在第一个实验中,图9的CNN已经使用在560名乳腺癌患者中包含的133个新鲜冷冻肿瘤样本子集中测得的低通WGS覆盖度数据进行训练,其中,BRCA缺陷状态和HRD状态(由HRDetect分数定义)是公众可得的,参考Nik-Zainal等人,“Landscape of somaticmutations in 560breast cancer whole-genome sequences”,Nature 534,47–54(2016)”-来自可在European Genome-phenome Archive EGA(https://www.ebi.ac.uk/ega/studies/EGAS00001001178)获得的Wellcome Trust Sanger Institute和International Cancer Genome Consortium ICGC的数据集,并且参考“HRDetect is apredictor of BRCA1 and BRCA2 deficiency based on mutational signatures”,H.Davies等人,Nature Medicine,2017年3月13日在线发布。这个由133个公众可得的真实样本组成的原始训练数据集被增强为3083个人工训练样本。进行数据增强是为了保持在原始数据集中观察到的相同的纯度-倍性分布。通过数据扩增获得的人工样本只用于训练和验证。原始数据集的很大一部分(来自560个乳腺癌研究的202个样本中的69个)被保留用于测试,没有参与数据增强过程。
每个肿瘤样本的原始BAM文件首先被下采样为1000万对端读数,以模拟低通WGS,覆盖度信号首先使用第一组100kbp的高分辨率覆盖区块进行预计算,并通过执行样本和GC正常化进行预处理。由此产生的归一化覆盖度信号(y轴)分别显示在图10a)中的HRD-样本和图10b)中的HRD+样本的全基因组(不包括性染色体)(x轴)。X轴给出了每个染色体边界的竖直虚线的基因组坐标。
使用图3的自适应区块化策略,将归一化的覆盖度数据进一步坍塌为2.5Mpb至3.5Mbp的平均3Mbp目标区块尺寸。图11示出了在应用根据本公开的一些实施方式的空间排列方法之前,对133个原始训练样本按其CSI状态排序的所得覆盖度数据为一维信号,每个样本一行,HR缺陷的样本在Y轴的顶部。X轴示出了每个染色体边界的竖直虚线的基因组坐标。
训练数据集中样本的归一化覆盖度数据被进一步在空间上排列,形成由84个区块*22条染色体组成的二维阵列,用于训练CNN模型以预测HRD状态。22条染色体的覆盖区块被绘制成22行,从chr1到chr22,相对于它们的着丝粒区块排列。对于较短的染色体臂,通过复制同一行中最近的端粒区块的值来填补空区块。图7a)示出了所产生的HRD样本的阵列,图7b)示出了所产生的HRD+样本的阵列。
图12示出了将经训练的CNN应用于202个样本的完整数据集所得到的结果(左边),其包括未用于训练的69个测试样本(表示为在CNN输出端测量的“BRCAness”状态指示,然后通过softmax转换为概率),并且示出了现有技术HRD分数(通过平均LOH、LST和TAI分数得到)的结果(右边),以HRDetect分数为基准的分类器。从图12所示的直方图中可以看出,无论对整个202个真实样本集进行测试还是仅对未用于训练机器学习模型的69个样本集进行测试,所提出的方法在预测HRDetect阳性样本方面都优于现有技术的HRD分数方法。
图13示出了测试集中69个样本中的每个样本的三个HRDetect分数(顶部)、HRD分数(底部)和在经训练的CNN(中间)输出处获得的SOPHiA CSI(此处表示为BRCAness)分数。在每个单独的图片中,样本按分数结果的增加从左到右排序。没有BRCA缺陷的样本显示为浅灰色。BRCA缺陷的样本显示为深灰色(不同的灰色水平用来表示导致BRCA缺陷的突变类型)。如图13所证实的,所有BRCA缺陷样本的HRDetect分数和SOPHiA BRCAness CSI分数都非常高。这一结果表明,CSI指标在识别BRCA缺陷样本方面优于HRD分数。
实验-实施例2-在单一湿式实验室工作流程中,对高等级卵巢浆液性腺癌Ovkate
细胞进行lpWGS与基于捕获的靶标板相结合
对OVKATE(RRID:CVCL_3110)细胞(https://web.expasy.org/cellosaurus/CVCL_3110)进行了低通量全基因组测序和基于捕获的靶标板。根据制造商的说明,使用DNeasy血液和组织试剂盒(Qiagen)提取DNA。然后使用SOPHiA基因文库制备试剂盒制备全基因组测序文库。全基因组文库的一部分被用来进行靶标富集,使用涵盖HRR途径相关基因的探针小组,并使用SOPHiA基因捕获方案进行。然后,富集的文库和全基因组文库一起在NextseqMid流动池(Illumina)上测序。两种类型的文库在流式细胞上被负载和平衡,以便在基因组上达到大约1-2X的覆盖度,在捕获板的目标区域达到超过1000X的覆盖度。
图15示出了分析lpWGS和定向捕获后获得的数据的示例性工作流程。
从培养的Ovkate细胞中提取的DNA按照结合低通WGS和靶向富集的NGS工作流程进行处理。对NGS数据进行生物信息学处理以产生BAM文件,其含有与参考基因组比对的读数信息。映射到靶向区域的读数可以进行变体调用,并测量富集区域,图15C(右,图15C中的C2),的等位基因比例。在掩盖了通过捕获富集方法靶向的区域后,来自lp-WGS区域的读数可以计算出全基因组的覆盖水平概况。后者被归一化以消除偏差,图15C(左,图15C中的C1),并用CNN处理以计算出HRD分数。Ovkate细胞被分类为HRD阳性(分数15.8)。
因此,结合lpWGS和基于捕获的靶标测序可以平行检测样本CSI状态和感兴趣基因的突变HRD状态。
实验-实施例3-获得增强的训练数据
在一个实验中,首先从Nik-Zainal等人,“Landscape of somatic mutationsin560breast cancer whole-genome sequences”,Nature 534,47–54(2016)的公众可得数据中选择了169个具有已知HRD状态的新鲜冷冻组织的真实样本。对于每个人工样本,可将具有相同HRD状态的随机数量的真实样本中的染色体随机组合,以创建数据增强(DA)训练样本。创建人工样本的样本数量可以是从指数分布中抽取的随机数,其使用:N=K*exp(-K*x),x是一个随机数,K=1/3。使用这样的方法可以确保只有有限数量的样本数据被组合成人工样本。
接下来,确定池中纯度/倍性比最低的样本,降低所有其他样本的纯度,使池中所有样本的纯度/倍性比等于这个比。这样可以防止引入从具有不同肿瘤纯度和倍性的样本中随机选择的染色体之间的覆盖幅度的差异(图16A)。为了达到这个目的,通过添加来自正常样本的测序读数,于电脑中(in-silico)降低纯度(图16B)。这一策略确保了在所有DA样本中,对于数据增强的样本,在所有染色体上观察到给定倍性的覆盖度差异幅度是恒定的(图16C),但在DA训练数据的纯度/倍性比中引入了偏差。特别是这些样本的纯度和倍性的分布往往低于对原始样本集的观察。应用Metropolis-Hastings和Gibbs抽样方法来考虑这个潜在的干扰因素,并确保由此产生的4403个保留DA样本的纯度/倍性分布与原始样本相匹配。
这样一组DA训练数据可用于训练机器学习模型。对于机器学习领域的技术人员来说显而易见的是,当可能需要支持某些类别的方法如HRD的性能的大量数据不可用时,这样的数据增强策略可以促进CSI分类器机器学习模型的训练。数据增强可以在电脑中准备,以增加训练数据集的大小和多样性。成功的数据增强策略可以保留增强数据集的关键属性,从而导致人工数据和真实数据无法区分。
所提出的方法的进一步优势和好处
所提出的机器学习方法进一步带来了额外优势,因为它们本身就能适应越来越多的个性化医学数据的可用性,特别是在肿瘤学实践中。它们不需要明确的生物学模型(例如优选是在某些染色体臂区,如端粒或着丝粒周围发生的疑似事件)来建立适合指导诊断、预后和/或治疗的预测特征。当采用半监督训练框架时,数据增强有利于开发更稳健的数据模型,对噪声数据或DNA样本中的低肿瘤含量不太敏感。随着更多数据的可用,模型可以被重新训练,而不需要重新架构运行中的基因组分析系统和工作流程,因为只有模型参数发生变化。最初为特定应用(如乳腺癌)开发的训练模型可以转移到其他应用(如卵巢癌、前列腺癌或胰腺癌)和/或其他样本类型(FFPE、FFT、cfDNA、ctDNA)。它们也可用于预测不同状态,包括对不同疗法和治疗的应答。
本文提出的方法适用于lpWGS测序数据,因此与使用高覆盖度(>30X)WGS或SNP阵列的现有技术方法相比,更便宜、更容易在标准临床实践中实施和部署。
其他实施方式和应用
虽然上面已经描述了各种实施方式,但应该理解,这些方案是以举例的方式提出的,而不是限制性。对相关领域的技术人员来说,显然可以在不偏离精神和范围的情况下对其进行形式和细节的各种改变。事实上,在阅读了上述描述之后,对于相关技术技术人员来说,如何实现替代性的实施方式是显而易见的。
虽然所提出的方法的示范性实施方式和应用已经与靶向下一代测序基因组分析有关,但对那些生物信息学领域的技术人员来说显而易见的是,它们也可以适用于从其他肿瘤样本基因组分析工作流程中检测和分类HR缺陷,例如使用SNP阵列和阵列CGH湿式实验室工作流程中的基因组数据。SNP阵列、阵列CGH或下一代测序(NGS)技术可用于产生随拷贝数变化的计数。此外,低通量WGS(0.1X至5X)或大型靶向富集测序文库,例如可从WES(全外显子组测序)或CES(毛细管电泳)板中获得,可单独使用或结合使用以产生提出的机器学习方法的输入数据。另外,低通量WGS可以与小型靶向板(基于扩增子的或基于捕获的)在单一或单独的湿室工作流程中结合。可能的话,从靶向测序方法发出的脱靶读数也可以作为主要的或补充的输入数据进入多维阵列,同时仍然很适合经训练的神经网络作为CSI分析仪123架构中的输入数据处理组件。
虽然所提出的方法的示例性实例和应用已被描述为与覆盖度数据信号阵列有关,该信号阵列被排列为以行排列的染色体覆盖区块的图像,并按其着丝粒区块竖直排列,但对于生物信息学领域的技术人员来说,显然,各种其他实施例也是可能的。图14示出了来自至少两个分析的染色体(染色体A和染色体B)的每个染色体臂(p臂和q臂)相对于至少一个与着丝粒区域或端粒区块相邻的区块的覆盖度信号图像的实例性比对。
在可能的实施方式中,染色体臂可表示为阵列中的行,所有臂的着丝粒区块可沿阵列的第一列对齐。在另一个可能的实施方式中,染色体臂可以在阵列中表示为行,所有臂的端粒区块可以沿阵列的第一列对齐。在另一个可能的实施方式中,染色体臂可表示为阵列中的行,所有臂的着丝粒区块可沿阵列的最后一列对齐。在另一个可能的实施方式中,染色体臂可表示为阵列中的行,所有臂的端粒区块可沿阵列的最后一列对齐。
在可能的实施方式中,染色体臂可表示为阵列中的列,所有染色体臂的着丝粒区块可沿阵列的第一行对齐。在另一个可能的实施方式中,染色体臂可表示为阵列中的列,所有臂的端粒区块可沿阵列的第一行对齐。在另一个可能的实施方式中,染色体臂可表示为阵列中的列,所有臂的着丝粒区块可沿阵列的最后一行对齐。在另一个可能的实施方式中,染色体臂可以表示为阵列中的列,所有臂的端粒区块可以沿阵列的最后一行对齐。
在可能的实施方式中,整个染色体可表示为阵列中的列,所有染色体的着丝粒区块可沿阵列的行对齐。着丝粒区块的对齐行可以在阵列的中心,也可以根据染色体组的p臂和q臂各自的区块长,从中心向上或向下移动。与染色体的p臂相关的区块和与q臂相关的区块可以分别在着丝粒区块行的上方和下方,或者反过来在着丝粒区块行的下方和上方。
在可能的实施方式中,整个染色体可表示为阵列中的行,所有染色体的着丝粒区块可沿阵列的列对齐。着丝粒区块的对齐列可以在阵列的中间,也可以根据染色体组的p臂和q臂各自的区块长从中间向右或向左偏移。与染色体的p臂相关的区块和与q臂相关的区块可以分别从着丝粒区块列向右和向左移动,或者反过来从着丝粒区块列向左和向右移动。
在可能的实施方式中,整个染色体可以在阵列中表示为行,所有染色体的两个端粒区块之一可以沿阵列的列对齐。在可能的实施方式中,端粒区块的对齐列可以是阵列的第一列。在另一个可能的实施方式中,端粒区块的对齐列可以是阵列的最后一列。在可能的实施方式中,p臂端粒区块可用于对齐。在另一个可能的实施方式中,q臂端粒区块可用于对齐。在另一个可能的实施方式中,可以为每条染色体单独选择p臂或q臂端粒区块的对齐。
在可能的实施方式中,整个染色体可以在阵列中表示为列,所有染色体的两个端粒区块中的一个可以沿着阵列的行进行对齐。在可能的实施方式中,端粒区块的对齐行可以是阵列的第一行。在另一个可能的实施方式中,端粒区块的对齐行可以是阵列的最后一行。在可能的实施方式中,p臂端粒区块可用于对齐。在另一个可能的实施方式中,q臂端粒区块可用于对齐。在另一个可能的实施方式中,可以为每条染色体单独选择p臂或q臂端粒区块的对齐。
基因组数据分析仪120计算机系统(本文中也为“系统”)120可被编程或以其他方式配置,以实现除本文所述CSI和HRD分析仪系统和方法之外的不同基因组数据分析方法,例如接收和/或结合测序数据、调用拷贝数改变和/或注释变体以进一步表征肿瘤样本。机器学习模型可以采用扩展的多维阵列输入,除了归一化覆盖度数据二维阵列输入外,还包括诸如GC含量、区块尺寸、映射性、映射质量、变体等位基因分数(VAF)等信息。机器学习模型还可以采用额外的标量输入,提供肿瘤含量信息(纯度)、样本倍性信息。基因组数据分析系统也可用于评估样本的质量(即测量FFPE降解的程度)。基因组数据分析系统还可用于根据癌症的CSI状态水平对其进行分类。基因组数据分析仪系统还可用于根据癌症的CSI状态和与免疫逃逸事件的关联对其进行分层(Bakhoum等人,2018,Cell 174(6),p.1347–1360)。
虽然所提出的方法的示例性实施方式和应用已经用CNN作为训练的机器学习模型进行了描述,但在可能的实施方式中(未图示),可以交替使用随机森林机器学习模型。随机森林或随机决策森林是一种用于分类、回归和其他任务的集合学习方法,它通过在训练时构建大量的决策树并输出类的模式(分类)或单个树的平均预测(回归)来操作。随机决策森林纠正了决策树对其训练集过度拟合的习惯,这使得它成为高维输入数据的良好候选方案。
除非另有定义,本文使用的所有技术和科学术语与本发明所属领域的普通技术人员通常理解的含义相同。本文描述中使用的术语仅用于描述特定的实施方式,并不意味着是限制性的。在说明书和所附权利要求书中,除非上下文明确指出,单数形式的:“一个/一种(a)”、“一个/一种(an)”和“该”也包括复数形式。
除非有相反的说明,以下说明书和所附权利要求书中列出的数字参数是近似值,可根据寻求获得的所需特性而变化,因此可由术语“约”来修改。至少,而且不是为了限制权利要求书范围内等价理论的应用,每个数字参数应根据有效数字的数量和普通的四舍五入方法来解释。
尽管提出广泛范围的数字范围和参数是近似值,但具体实例中提出的数值是尽可能精确的报告。然而,任何数值都必然包含某些误差,这些误差是由各自测试测量中发现的标准偏差造成的。本说明书中给出的每一个数值范围都将包括落在这种较宽的数值范围内的每一个较窄的数值范围,就像这种较窄的数值范围都明确地写在这里一样。
对于数字数据通信领域的技术人员来说是显而易见的,这里描述的方法可以无差别地应用于各种数据结构,如数据文件或数据流。因此,术语“数据”、“数据集”、“数据结构”、“数据字段”、“文件”或“流”可以在本说明书中无差别地使用。
尽管上面的详细描述包含了许多具体的细节,但这些不应该被理解为限制了本发明的实施方式范围,而只是提供了几个实施方式中的一些说明。
此外,应该理解的是,任何突出功能和优点的图都只是出于举例的目的。所公开的方法具有足够的灵活性和可配置性,因此它们可以以所示以外的方式加以利用。
Claims (21)
1.一种确定受试者DNA样本的同源重组缺陷(HRD)状态的基于计算机的方法,所述方法包括:
-获得一组待分析的受试者DNA样本的全基因组测序读数,
-将这一组受试者DNA样本的测序读数与参考基因组进行比对,其中,所述参考基因组被分为多个区块,每个区块属于待分析的全基因组染色体中染色体臂的同一基因组区域,
-对沿着每个染色体臂的每个区块中的比对读数的数值进行计数并归一化以获得所述染色体臂的覆盖度信号,
-将染色体臂的覆盖度信号排列成受试者DNA样本的覆盖度数据信号阵列,
-将覆盖度数据信号阵列输入到经训练的机器学习模型中,其中,所述模型已使用一组已知同源重组缺陷状态的样本进行训练,以区分来自同源重组缺陷状态为阳性的样本的覆盖度数据信号阵列和来自同源重组缺陷状态为阴性的样本的覆盖度数据信号阵列,
从而确定受试者DNA样本的同源重组缺陷分数(HRD分数),以及
根据所述经训练的机器学习模型的HRD分数来确定受试者DNA样本的阴性、阳性或不确定的同源重组缺陷(HRD)状态。
2.根据权利要求1所述的方法,其中,一组测序读数是从全基因组测序中获得的,其中,读数深度覆盖度至多为30X。
3.根据权利要求2所述的方法,其中,一组测序读数是从低通量全基因组测序中获得的,其中,读数深度覆盖度至少为0.1X且至多为5X。
4.根据权利要求1至3所述的方法,其中,对沿着每个染色体臂的每个区块中的比对读数的数值进行计数并归一化以获得所述染色体臂的覆盖度信号包括对每个样本的覆盖度信号进行归一化,和/或通过GC含量进行归一化以应用GC-偏移校正。
5.根据权利要求1至4所述的方法,其中,所述染色体臂的覆盖度信号被排列成一维覆盖度数据信号矢量或二维覆盖度数据信号图像。
6.根据权利要求5所述的方法,其中,所述染色体臂的覆盖度信号通过按行比对每条染色体的覆盖度数据信号和每条染色体臂的着丝粒区块,即,与染色体臂的着丝粒区域相邻的最近区块,而被排列成二维覆盖度数据信号图像。
7.根据权利要求1至6所述的方法,其中,所述机器学习模型是使用一组具有已知同源重组缺陷状态的肿瘤数据样本作为训练标签进行预先训练的。
8.根据权利要求7所述的方法,其中,所述训练数据集用通过结合数据样本的染色体的数据与已知同源重组缺陷状态标签而产生的人工样本数据进行增强。
9.根据权利要求8所述的方法,其中,生成数据增强的样本以表示如在真实样本数据集中观察到的纯度-倍性比分布。
10.根据权利要求1所述的方法,其中,所述参考基因组被划分为第一组至多100kbp的区块,并还包括在排列每个染色体臂上的覆盖度信号之前,将100kbp区块坍塌成第二组至少500kbp的较大区块的步骤。
11.根据权利要求10所述的方法,其中,第一组区块的区块具有至多100kbp的均匀尺寸,而第二组区块的区块具有2.5至3.5Mbp的尺寸并且通过池化第一组区块中25至35个100kbp的区块得到。
12.一种确定患者DNA样本的同源重组缺陷(HRD)状态的体外方法,所述方法包括:
-提供来自患者样本的DNA片段;
-构建包含与一组染色体重叠的所述片段的文库;
-对所述文库进行测序,以达到至多30X的全基因组测序覆盖度,优选达到至少0.1X且至多5X的基因组测序覆盖度;
并且基于根据权利要求1获得的经训练的机器学习模型的分析来确定患者样本的HRD状态。
13.根据前述权利要求中任何一项所述的方法,其中,所述患者DNA样本是肿瘤细胞外游离DNA(cfDNA)、新鲜冷冻组织(FFT)或福尔马林固定石蜡包埋(FFPE)样本。
14.根据前述权利要求中任何一项所述的方法,其中,所述患者样本的HRD分数或HRD状态是对癌症治疗方案的肿瘤应答的预测因子。
15.根据权利要求15所述的方法,其中,所述癌症治疗方案选自由烷基化剂、铂类化疗剂、卡铂、顺铂、异丙铂、奈达铂、奥沙利铂、吡铂、氮芥、苯丁酸氮芥、美法仑、环磷酰胺、异环磷酰胺、雌莫司汀、卡莫司汀、洛莫司汀、福莫司汀、链脲菌素、白消安、哌泊溴烷、甲基苄肼、达卡巴嗪、塞替派、替莫唑胺和/或其他抗肿瘤的铂配合化合物、DNA损伤剂、放射治疗剂、蒽环类药物、表柔比星、多柔比星;拓扑异构酶I抑制剂、喜树碱、拓扑替康、伊立替康、PARP(聚ADP-核糖聚合酶)抑制剂、奥拉帕尼、卢卡帕尼、尼拉帕尼、他拉唑帕尼、伊尼帕利、CEP9722、MK 4827、BMN-673、3-氨基苯甲酰胺、维拉匹里、帕米帕利或E7016所组成的组。
16.一种选择癌症患者用铂类化疗剂、DNA损伤剂、蒽环类药物、拓扑异构酶I抑制剂、PARP抑制剂进行治疗的方法,包括根据权利要求1所述的方法检测肿瘤患者样本为HRD阳性的步骤。
17.根据权利要求16所述的方法,其中,所述患者患有选自高级别浆液性卵巢癌、前列腺癌、乳腺癌或胰腺癌的癌症。
18.一种训练机器学习算法以确定受试者DNA样本的同源重组缺陷(HRD)状态的方法,所述方法包括:
向机器学习监督训练算法输入来自已知同源重组缺陷状态为阳性的样本的覆盖度数据信号阵列和来自已知同源重组缺陷状态为阴性的样本的覆盖度数据信号阵列。
19.根据前述权利要求中任何一项所述的方法,其中,经训练的机器学习模型是随机森林模型、神经网络模型、深度学习分类器或卷积神经网络模型。
20.根据权利要求19所述的方法,其中,神经网络模型经训练的机器学习模型是卷积神经网络,所述卷积神经网络被训练成在其输出端产生阳性或阴性的HRD状态的单标签二元分类,或阳性、阴性或不确定的HRD状态的单标签多类别分类,或代表HRD状态的标量HRD分数。
21.根据前述权利要求中任何一项所述的方法,其中,所述机器学习模型已使用数据增强集在半监督模式下训练,所述数据增强集通过从共享相同HRD状态及相同归一化纯度和倍性比的一组真实样本的染色体取样数据而产生。
Applications Claiming Priority (3)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
EP20187813.9 | 2020-07-27 | ||
EP20187813.9A EP3945525A1 (en) | 2020-07-27 | 2020-07-27 | Methods for identifying chromosomal spatial instability such as homologous repair deficiency in low coverage next-generation sequencing data |
PCT/EP2021/071073 WO2022023381A1 (en) | 2020-07-27 | 2021-07-27 | Methods for identifying chromosomal spatial instability such as homologous repair deficiency in low coverage next-generation sequencing data |
Publications (1)
Publication Number | Publication Date |
---|---|
CN116194995A true CN116194995A (zh) | 2023-05-30 |
Family
ID=71833187
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202180060052.6A Pending CN116194995A (zh) | 2020-07-27 | 2021-07-27 | 在低覆盖度的下一代测序数据中识别染色体空间不稳定性如同源修复缺陷的方法 |
Country Status (9)
Country | Link |
---|---|
US (3) | US20220028481A1 (zh) |
EP (2) | EP3945525A1 (zh) |
JP (1) | JP2023535962A (zh) |
KR (1) | KR20230045009A (zh) |
CN (1) | CN116194995A (zh) |
AU (1) | AU2021314892A1 (zh) |
BR (1) | BR112023000014A2 (zh) |
CA (1) | CA3185856A1 (zh) |
WO (1) | WO2022023381A1 (zh) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN117373678A (zh) * | 2023-12-08 | 2024-01-09 | 北京望石智慧科技有限公司 | 基于突变签名的疾病风险预测模型构建方法及分析方法 |
Families Citing this family (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20230265525A1 (en) * | 2022-02-22 | 2023-08-24 | Genegeniedx Corp | Methods for detecting homologous recombination deficiency in cancer patients |
EP4297037A1 (en) * | 2022-06-24 | 2023-12-27 | Seqone | Device for determining an indicator of presence of hrd in a genome of a subject |
WO2024050366A1 (en) * | 2022-08-30 | 2024-03-07 | Foundation Medicine, Inc. | Systems and methods for classifying and treating homologous repair deficiency cancers |
CN115330603B (zh) * | 2022-10-17 | 2023-01-20 | 湖南自兴智慧医疗科技有限公司 | 基于深度学习卷积神经网络的人类染色体图像摆正方法 |
CN116129123B (zh) * | 2023-02-27 | 2024-01-05 | 中国矿业大学 | 基于不确定度校准和区域分解的端到端染色体分割方法 |
CN116646010B (zh) * | 2023-07-27 | 2024-03-29 | 深圳赛陆医疗科技有限公司 | 人源性病毒检测方法及装置、设备、存储介质 |
Family Cites Families (12)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
ES2588979T3 (es) | 2010-08-24 | 2016-11-08 | Dana-Farber Cancer Institute, Inc. | Procedimientos para la predicción de una respuesta contra el cáncer |
EP2794907B2 (en) | 2011-12-21 | 2022-11-23 | Myriad Genetics, Inc. | Methods and materials for assessing loss of heterozygosity |
CA2864481C (en) | 2012-02-23 | 2020-07-14 | The Children's Hospital Corporation | Methods for predicting anti-cancer response |
EP2859118B1 (en) | 2012-06-07 | 2017-11-22 | Institut Curie | Methods for detecting inactivation of the homologous recombination pathway (brca1/2) in human tumors |
ES2777228T3 (es) | 2013-04-05 | 2020-08-04 | Myriad Genetics Inc | Métodos para evaluar la deficiencia de recombinación homóloga y predecir la respuesta al tratamiento del cáncer |
ES2800673T3 (es) | 2014-08-15 | 2021-01-04 | Myriad Genetics Inc | Métodos y materiales para evaluar una deficiencia de recombinación homóloga |
WO2017191076A1 (en) * | 2016-05-01 | 2017-11-09 | Genome Research Limited | Method of characterising a dna sample |
JP7224185B2 (ja) * | 2016-05-01 | 2023-02-17 | ゲノム・リサーチ・リミテッド | Dnaサンプルを特徴付ける方法 |
US11923049B2 (en) | 2016-06-22 | 2024-03-05 | Sophia Genetics S.A. | Methods for processing next-generation sequencing genomic data |
US11482303B2 (en) * | 2018-06-01 | 2022-10-25 | Grail, Llc | Convolutional neural network systems and methods for data classification |
US11164655B2 (en) * | 2019-12-10 | 2021-11-02 | Tempus Labs, Inc. | Systems and methods for predicting homologous recombination deficiency status of a specimen |
JP2023526252A (ja) * | 2020-05-14 | 2023-06-21 | ガーダント ヘルス, インコーポレイテッド | 相同組換え修復欠損の検出 |
-
2020
- 2020-07-27 EP EP20187813.9A patent/EP3945525A1/en not_active Withdrawn
-
2021
- 2021-07-27 WO PCT/EP2021/071073 patent/WO2022023381A1/en active Application Filing
- 2021-07-27 AU AU2021314892A patent/AU2021314892A1/en active Pending
- 2021-07-27 JP JP2023505760A patent/JP2023535962A/ja active Pending
- 2021-07-27 US US17/386,255 patent/US20220028481A1/en active Pending
- 2021-07-27 EP EP21754754.6A patent/EP4189685A1/en active Pending
- 2021-07-27 CA CA3185856A patent/CA3185856A1/en active Pending
- 2021-07-27 BR BR112023000014A patent/BR112023000014A2/pt unknown
- 2021-07-27 CN CN202180060052.6A patent/CN116194995A/zh active Pending
- 2021-07-27 KR KR1020237002728A patent/KR20230045009A/ko unknown
- 2021-11-23 US US17/534,368 patent/US20220084626A1/en active Pending
-
2022
- 2022-03-07 US US17/688,791 patent/US20220310199A1/en active Pending
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN117373678A (zh) * | 2023-12-08 | 2024-01-09 | 北京望石智慧科技有限公司 | 基于突变签名的疾病风险预测模型构建方法及分析方法 |
CN117373678B (zh) * | 2023-12-08 | 2024-03-05 | 北京望石智慧科技有限公司 | 基于突变签名的疾病风险预测模型构建方法及分析方法 |
Also Published As
Publication number | Publication date |
---|---|
EP4189685A1 (en) | 2023-06-07 |
BR112023000014A2 (pt) | 2023-02-07 |
US20220310199A1 (en) | 2022-09-29 |
KR20230045009A (ko) | 2023-04-04 |
US20220084626A1 (en) | 2022-03-17 |
JP2023535962A (ja) | 2023-08-22 |
US20220028481A1 (en) | 2022-01-27 |
EP3945525A1 (en) | 2022-02-02 |
CA3185856A1 (en) | 2022-02-03 |
WO2022023381A1 (en) | 2022-02-03 |
AU2021314892A1 (en) | 2023-03-09 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN116194995A (zh) | 在低覆盖度的下一代测序数据中识别染色体空间不稳定性如同源修复缺陷的方法 | |
JP7247253B2 (ja) | 経験的バリアントスコア(evs)ベースの深層学習バリアントコーラ | |
KR102662206B1 (ko) | 심층 학습 기반 비정상 스플라이싱 검출 | |
KR102562419B1 (ko) | 심층 신경망에 기반한 변이체 분류자 | |
Amaratunga et al. | Exploration and analysis of DNA microarray and protein array data | |
JP2022521492A (ja) | 相同組換え欠損を推定するための統合された機械学習フレームワーク | |
CN112048559A (zh) | 基于m6A相关的IncRNA网络胃癌预后的模型构建及临床应用 | |
KR102628141B1 (ko) | 서열-특정 오류(sse)를 유발시키는 서열 패턴을 식별하기 위한 심층 학습-기반 프레임워크 | |
KR20210110241A (ko) | 인간백혈구항원 하플로타입 기반 다중 분류 인공지능 모델을 이용한 면역항암제 적응증 및 반응 예측 시스템 및 방법 | |
Lock et al. | Bayesian genome-and epigenome-wide association studies with gene level dependence | |
Gupta et al. | A new deep learning technique reveals the exclusive functional contributions of individual cancer mutations | |
US20230316054A1 (en) | Machine learning modeling of probe intensity | |
Patel | Cancer risk prediction with next generation sequencing data using machine learning | |
Narain et al. | Biometry in Genomics | |
Bakewell et al. | Toward'smart'DNA microarrays: algorithms for improving data quality and statistical inference |
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 |