CN115083521B - 一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法及*** - Google Patents
一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法及*** Download PDFInfo
- Publication number
- CN115083521B CN115083521B CN202210865067.6A CN202210865067A CN115083521B CN 115083521 B CN115083521 B CN 115083521B CN 202210865067 A CN202210865067 A CN 202210865067A CN 115083521 B CN115083521 B CN 115083521B
- Authority
- CN
- China
- Prior art keywords
- mutation
- cell
- data
- site
- tumor
- 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
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
- G16B20/50—Mutagenesis
-
- 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
- G16B25/00—ICT specially adapted for hybridisation; ICT specially adapted for gene or protein expression
- G16B25/20—Polymerase chain reaction [PCR]; Primer or probe design; Probe optimisation
-
- 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
- G16B30/10—Sequence alignment; Homology search
Landscapes
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Health & Medical Sciences (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Engineering & Computer Science (AREA)
- Chemical & Material Sciences (AREA)
- General Health & Medical Sciences (AREA)
- Biophysics (AREA)
- Theoretical Computer Science (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Medical Informatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Biotechnology (AREA)
- Evolutionary Biology (AREA)
- Analytical Chemistry (AREA)
- Proteomics, Peptides & Aminoacids (AREA)
- Molecular Biology (AREA)
- Genetics & Genomics (AREA)
- Chemical Kinetics & Catalysis (AREA)
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
- Investigating Or Analysing Biological Materials (AREA)
Abstract
本发明公开一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法,包括:获取待测样品的单细胞转录组测序数据并基于单细胞转录组测序数据的分析获得第一数据;获取待测样品的突变位点信息;基于第一数据及突变位点信息进行突变位点的突变分析及肿瘤细胞类群的鉴定;基于肿瘤细胞类群的鉴定及突变位点的突变分析获得待测样品的肿瘤细胞类群统计信息。还公开了对应***、电子设备及计算机可读存储介质,快速利用突变位点鉴定单细胞转录组测序数据中肿瘤细胞类群,包括基于已知肿瘤的单细胞转录组测序数据和突变位点信息,在单细胞水平进行突变分析,分析所有细胞类群的位点突变频率,实现鉴定肿瘤细胞类群,还可以从单细胞水平分析肿瘤细胞的异质性。
Description
技术领域
本发明涉及医疗和生物技术领域,尤其涉及一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法及***。
背景技术
单细胞转录组测序(single cell RNA-Sequence,scRNA-Seq)是近几年出现的在单个细胞水平对转录组进行高通量测序的技术,其一次可以测序几千到上万个细胞转录组表达。随着单细胞转录组测序技术的出现与不断改进,在单细胞分辨率下研究肿瘤的基因组及表达组特征成为可能。单细胞转录组测序在肿瘤研究中可以探索肿瘤异质性、肿瘤耐药机制、免疫治疗等方面,已经广泛应用在各种肿瘤研究中。
通过scRNA-Seq可以获取肿瘤组织中所检测每一个细胞的表达值,通过每个细胞的基因表达值,利用无监督聚类将所检测的细胞分成不同的类别(cluster),通过细胞标志物(marker),获得每个类别(cluster)所属的细胞类型,所属的细胞类型包括免疫细胞(B细胞或T细胞等)、基质细胞、间充质细胞、干细胞或表皮细胞等。根据每个cluster的基因表达情况,获得肿瘤细胞亚群的情况。无监督聚类属于非监督性技术,一般分为两个步骤:首先,基于单细胞转录组测序数据,在特定长度的区域内,为每个细胞估计出拷贝数变异的方向和程度;然后,基于拷贝数变异的相关信息,采用无监督聚类的方法,将所有细胞聚为两类,取拷贝数变异程度较大的一类作为恶性细胞。虽然通过细胞标志物(marker),可以区分免疫细胞和非免疫细胞,但是由于取样的时候不能全部取到肿瘤组织,在肿瘤组织细胞中有一部分是正常细胞,所以非免疫细胞类别(cluster)里面会同时存在肿瘤细胞和正常细胞,通过基因表达很难区分哪一些类别(cluster)细胞是正常细胞,哪一些类别(cluster)是肿瘤细胞。
发明内容
为了解决现有技术中存在的问题,本发明提供了如下技术方案,一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法及***,快速利用突变位点鉴定单细胞转录组测序(scRNA-Seq)数据中肿瘤细胞类群,包括基于已知肿瘤的单细胞转录组测序和突变位点信息,在单细胞水平进行突变分析,分析所有细胞类群(cluster)的位点的突变频率,从而实现鉴定肿瘤细胞类群,分析肿瘤异质性。
本发明一方面提供了一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法,包括:
S1,获取待测样品的单细胞转录组测序数据并基于所述单细胞转录组测序数据的分析获得第一数据;
S2,获取待测样品的突变位点信息;
S3,基于所述第一数据以及突变位点信息进行突变位点的突变分析以及所述肿瘤细胞类群的鉴定;
S4,基于所述肿瘤细胞类群的鉴定以及突变位点的突变分析获得所述待测样品的所述肿瘤细胞类群统计信息。
优选的,所述获取待测样品的单细胞转录组测序数据包括:从因美纳公司Bio-Rad单细胞测序方法(illumina® Bio-Rad® Single-Cell Sequencing Solution)、BD公司Rhapsody单细胞分析***(BD Rhapsody™ Single-Cell Analysis System)、10xgenomics公司铬单细胞测序方法(10x Chromium Single Cell Gene ExpressionSolution)、ICELL8单细胞制备***(ICELL8 Single-Cell System)和/或C1单细胞制备***(C1™)获取待测样品的单细胞转录组测序数据。
优选的,所述第一数据包括:
基因组比对结果文件;
细胞条形码文件(barcode);以及
细胞聚类结果。
优选的,所述基因组比对结果文件为bam文件。
优选的,所述S2,获取待测样品的突变位点信息包括:
获取待测样品的肿瘤位点突变的基因组位置信息、待测样品的脱氧核糖核酸(DNA)体细胞突变数据和热点突变数据,其中所述基因组位置信息对应的基因组与所述基因组比对结果文件中的基因组完全一致。
优选的,所述待测样品的突变位点信息由基因突变检测数据或先验知识获得,所述先验知识包括:
(1)肿瘤外显子测序(WES)或者特定基因组套(panel)测序;
(2)公共数据库中记载的热点突变,所述公共数据库包括癌症基因组图谱(TCGA,The Cancer Genome Atlas)或肿瘤体细胞突变数据库(COSMIC);
(3)文章或者数据库中已经公开的相关肿瘤突变数据。
优选的,所述S3,基于所述第一数据以及突变位点信息进行突变位点的突变分析以及所述肿瘤细胞类群的鉴定包括:
S31,基于所述基因组比对结果文件进行碱基矫正,以矫正测序产生的噪音,从而精准分析突变位点的突变情况,包括:分析所述基因组比对结果文件中每一个细胞在突变位点位置的比对情况,将单细胞转录组测序数据中具有相同唯一分子标签(UMI)的测序读长(reads)片段聚成同一唯一分子标签(UMI)簇;对聚成同一唯一分子标签(UMI)簇的多个测序读长(reads)片段,在突变位点位置的所述比对情况进行判断,所述判断按照碱基矫正程序进行:
(1)如果唯一分子标签(UMI)簇中多个所述测序读长(reads)片段全部为相同的碱基,则确定所述唯一分子标签(UMI)簇为突变碱基;
(2)如果唯一分子标签(UMI)簇中多个所述测序读长(reads)片段包括不同的碱基,且其中比例最大的碱基占比超过80%,则所述唯一分子标签(UMI)簇为最大比例碱基;
(3)如果唯一分子标签(UMI)簇中多个所述测序读长(reads)片段包括不同的碱基,且其中比例最大的碱基占比低于80%,则放弃使用所述唯一分子标签(UMI)簇的信息;
针对所有唯一分子标签(UMI)簇按照所述碱基矫正程序(1)-(3)顺次进行判断,从而基于每个细胞的所有唯一分子标签(UMI)簇进行碱基矫正后获得矫正结果;
S32,基于矫正结果进行突变位点的突变分析,包括:确定参考基因,如果矫正结果中存在唯一分子标签(UMI)簇与参考碱基不一致,则确定细胞在所述突变位点存在突变;或确定多个突变基因,分别统计每个突变位点上所述多个突变碱基的唯一分子标签(UMI)数量,如果任意一个突变碱基的唯一分子标签(UMI)数量大于0,则确定所述突变位点存在细胞突变;
S33,基于所述突变位点的突变分析进行所述肿瘤细胞类群的鉴定。
优选的,所述S4,基于所述肿瘤细胞类群的鉴定以及突变位点的突变分析获得所述待测样品的所述肿瘤细胞类群统计信息包括:
S41,统计待测样品中每一个细胞类群中,携带突变位点的细胞数量;
S42,基于所述携带突变位点的细胞数量确定每一个细胞类群中所携带的肿瘤细胞数量以及类群特异的肿瘤细胞突变谱。
优选的,所述S3与S4之间还包括:对所述突变位点进行多次统计检验,控制所述突变位点的假阴性和假阳性结果的产生,包括:
对于每个突变位点,在所述突变位点对应的基因中,随机选择N个位点数据;
对于所述N个位点按照顺序执行所述S1-S3获得所述N个位点的突变情况,作为所述N个位点的背景值;
基于所述背景值构建背景噪声统计模型;
应用卡方检验(Chi-square test)或费舍尔精确检验(Fisher exact test),基于所述背景噪声统计模型以及所述N个位点的突变情况,排除多个干扰和错误,包括:
计算所述N个位点的突变情况统计显著性,排除测序错误产生的干扰;
基于所述N个位点的非免疫细胞突变情况,排除核糖核酸(RNA)编辑产生的干扰;
统计各个细胞类群所有的突变位点的突变频率,基于费舍尔精确检验,排除建库过程中聚合酶链式反应(PCR)产生的错误;
将免疫细胞对应的细胞类群合并,通过费舍尔精确检验比较各个突变位点非免疫细胞类群和免疫细胞类群的突变细胞比例P,将P值小于第一阈值的确定为肿瘤细胞类群备选集合;统计每个所述肿瘤细胞类群备选集合中所述突变细胞比例P值小于第一阈值的位点数目,基于费舍尔精确检验,确定非免疫细胞比例最高的为最终的肿瘤细胞类群。
优选的,所述第一阈值为0.05。
本发明的第二方面,提供一种单细胞转录组测序数据中肿瘤细胞类群的鉴定***,包括:
测序数据获取模块,获取待测样品的单细胞转录组测序数据并基于所述单细胞转录组测序数据的分析获得第一数据;
突变位点获取模块,用于获取待测样品的突变位点信息;
突变分析及肿瘤细胞类群鉴定模块,用于基于所述第一数据以及突变位点信息进行突变位点的突变分析以及所述肿瘤细胞类群的鉴定;
统计模块,用于基于所述肿瘤细胞类群的鉴定以及突变位点的突变分析获得所述待测样品的所述肿瘤细胞类群统计信息。
本发明的第三方面提供一种电子设备,包括处理器和存储器,所述存储器存储有多条指令,所述处理器用于读取所述指令并执行如第一方面所述的方法。
本发明的第四方面提供一种计算机可读存储介质,所述计算机可读存储介质存储有多条指令,所述多条指令可被处理器读取并执行如第一方面所述的方法。
本发明提供的单细胞转录组测序数据中肿瘤细胞类群的鉴定方法、***和电子设备,具有如下有益效果:
(1)快速利用突变位点鉴定单细胞转录组测序(scRNA-Seq)数据中肿瘤细胞类群,基于已知肿瘤的单细胞转录组测序和突变位点信息,在单细胞水平进行突变分析,分析所有细胞类群(cluster)的位点的突变频率,从而实现鉴定肿瘤细胞类群,分析肿瘤异质性;
(2)可用于任何突变和任何肿瘤,并且能够通过鉴定的肿瘤类群,可以探索肿瘤的异质性;
(3)利用单细胞核糖核酸(RNA)测序数据中唯一分子标签(UMI)信息,构建唯一分子标签(UMI)成簇,矫正测序产生的噪音, 从而精准分析位点突变情况;
(4) 多次统计检验控制假阴性和假阳性结果的产生。
附图说明
图1为本发明所述的单细胞转录组测序数据中肿瘤细胞类群的鉴定方法流程示意图。
图2为本发明提供的单细胞转录组测序数据中肿瘤细胞类群的鉴定***原理结构图。
图3为本发明提供的bam文件格式下比对结果的示意图。
图4所示为本发明提供的待测样品bam文件数据截图。
图5所示是cluster0、cluster4、cluster6突变细胞在所有位点的突变情况。
图6为本发明提供的电子设备一种实施例的结构示意图。
具体实施方式
为了更好的理解上述技术方案,下面将结合说明书附图以及具体的实施方式对上述技术方案做详细的说明。
本发明提供的方法可以在如下的终端环境中实施,该终端可以包括一个或多个如下部件:处理器、存储器和显示屏。其中,存储器中存储有至少一条指令,所述指令由处理器加载并执行以实现下述实施例所述的方法。
处理器可以包括一个或者多个处理核心。处理器利用各种接口和线路连接整个终端内的各个部分,通过运行或执行存储在存储器内的指令、程序、代码集或指令集,以及调用存储在存储器内的数据,执行终端的各种功能和处理数据。
存储器可以包括随机存储器(Random Access Memory,RAM),也可以包括只读存储器(Read-Only Memory,ROM)。存储器可用于存储指令、程序、代码、代码集或指令。
显示屏用于显示各个应用程序的用户界面。
除此之外,本领域技术人员可以理解,上述终端的结构并不构成对终端的限定,终端可以包括更多或更少的部件,或者组合某些部件,或者不同的部件布置。比如,终端中还包括射频电路、输入单元、传感器、音频电路、电源等部件,在此不再赘述。
实施例一
如图1所示,本实施例提供了一种单细胞转录组测序数据中肿瘤细胞类群(cluster)的鉴定方法,包括:
S1,获取待测样品的单细胞转录组测序数据并基于所述单细胞转录组测序数据的分析获得第一数据;
S2,获取待测样品的突变位点信息;
S3,基于所述第一数据以及突变位点信息进行突变位点的突变分析以及所述肿瘤细胞类群的鉴定;
S4,基于所述肿瘤细胞类群的鉴定以及突变位点的突变分析获得所述待测样品的所述肿瘤细胞类群统计信息。
作为优选的实施方式,所述获取待测样品的单细胞转录组测序数据包括:从因美纳公司Bio-Rad单细胞测序方法(illumina® Bio-Rad® Single-Cell SequencingSolution)、BD公司Rhapsody单细胞分析***(BD Rhapsody™ Single-Cell AnalysisSystem)、10x genomics公司铬单细胞测序方法(10x Chromium Single Cell GeneExpression Solution)、ICELL8单细胞制备***(ICELL8 Single-Cell System)和/或C1单细胞制备***(C1™)获取待测样品的单细胞转录组测序数据。
作为优选的实施方式,所述第一数据包括:
基因组比对结果文件;
细胞条形码文件(barcode);以及
细胞聚类结果。
作为优选的实施方式,所述基因组比对结果文件为bam文件。
随着生物信息数据的爆发式增长,存储生物信息的文件格式也多样化起来,不同的文件格式往往有不同的目的:为了软件之间的兼容性以及人类可读性而采用的用于数据加工、解析和处理的格式,比如.tsv,.csv等;为了提高计算机计算效率的数据格式,一般可读性不强的二进制文件,比如本实施例采用的bam文件。bam文件是sam文件的二进制格式,顺序对齐/映像格式文件(sam,Sequence Alignment/Map Format)产生于比对之后的数据输出,记录了比对的具体情况。文件中以tab键分割,包括上下两部分:
头文件部(Header section)和对齐部(Alignments section)
1、头文件部(Header section)
该部分全部以“@”开头,提供基本的软件版本,参考序列信息,排序信息等
@HD行:这一行中有各种不同的标识
标识“VN”用以说明格式版本
标识“SO”用以说明比对排序的情况,有未知(unknown (default))、未分类(unsorted)、队列名(queryname)和坐标(coordinate)四个选项,对于坐标(coordinate)选项,排序的主键是对齐部(Alignments section)的第三列“RNAME”,其顺序由@SQ行的“SN”标识的顺序定义,次要排序键是对齐部(Alignments section)的第四列“S”字段。对于RNAME和POS相等的比对,排列顺序则是任意的;
@SQ行的“SN”标签是参考序列说明,它的值主要是用于对齐部(Alignmentssection)的第三列“RNAME”和第七列“MRNM”比对的记录;
@PG行是使用的程序说明;该行“ID”为程序记录标识符,“PN”为程序名字,“CL”为命令行;
@CO行是任意的说明信息。
2、对齐部(Alignments section)
该部分包含了11列必需字段,无效或者没有的字段一般用“0”或者“*”表示;将比对情况记录成bam文件格式的形式如图3所示。
图3中的对齐部(Alignments section)有6行12列信息详细介绍了6条read的比对情况,其中前11列为必需字段,每列的含义简单汇总如下。
第1列:读长(Read)的名字 (Qname)
第2列:每一个测序读长(read)的比对情况(FLAG),可以用十进制数字(或者十六进制数字)表示,如果比对情况有多个,将多个比对情况所代表的十进制数字加和就是这一行的比对情况(FLAG)。如,上图3中r001的比对情况(FLAG)是99(1+2+32+64),则表示了“该测序读长(read)是对读长(pair read)中的一个”,“对读长(pair read)中每个都能够正确比对上”,“该测序读长(read)的匹配读长(mate read)的反向互补可以比对上”,“该测序读长(read)是对读长(pair read)中的测序读长1(read1)”;r001的另一个比对情况(FLAG)是147(1+2+16+128),则表示“该测序读长(read)是对读长(pair read)中的一个”,“对读长(pair read)中每个都能够正确比对上”,“该测序读长(read)是原测序读长(read)的反向互补”,“该测序读长(read)是对读长(pair read)中的测序读长2(read2)”(也就是说,该测序读长(read)是测序读长2(read2)的反向互补序列)。值得注意的是,r001是对读长(pairread),而且都能比对上,所以r001出现了两次,如果r001的测序读长1(read1)比对到参考序列的2个地方,r001的名字则会出现三次;如果测序读长1(read1)比对上一次,测序读长2(read2)没有比对上,r001仍会出现2次,不过,其中一个r001的第三列为“*”;所以对端(pair-end)测序,测序读长1(read1)文件和测序读长2(read2)文件同时映射(mapping),相同测序读长(reads)的id最少出现2次。
第3列:比对上的参考序列的名字(RNAME),该名字出现在头文件部(Headersection)的@SQ行的SN标识中,如果该测序读长(read)没有比对上,也就是说该测序读长(read)在参考序列上没有坐标,那么这一列则用“”表示,那么这一行的POS和CIGAR列也会是“”。
第4列:比对到的参考序列“RNAME”最左侧的位置坐标(POS read),也是CIGAR中第一个比对标识“M”对应的最左侧碱基在参考序列的位置,未比对上的read在参考序列中没有坐标,此列标识为“0”。
第5列:比对的质量值(MAPQ),计算方法为比对错误率的-10*log10的值,一般是取四舍五入的整数值,如果是255,说明该比对值无效。
第6列:表示测序读长(read)中每个碱基的比对情况标识符CIGAR(CIGAR)。
第7列:该测序读长(read)的匹配读长(mate read)比对上的参考序列的名字(MRNM),该名字出现在头文件部(Header section)的@SQ行的SN标识中,
如果和该测序读长(read)所在行的第三列“RNAME”一样,则用“=”表示,说明这对测序读长(read)比对到了同一条参考序列上;
如果匹配读长(mate read)没有比对上,第七列则用“*”表示;
如果这对测序读长(read)没有比对到同一条参考序列,那么这一列则是匹配读长(mate read)所在行第三列的“RNAME”。
第8列:该测序读长(read)的匹配读长(mate read)比对到的参考序列“RNAME”最左侧的位置坐标(MPOS),也是匹配读长(mate read)的标识(CIGAR)中第一个比对标识“M”对应的最左侧碱基在参考序列的位置,未比对上的测序读长(read)在参考序列中没有坐标,此列标识为“0”。
第9列:表示pair read完全匹配到同一条参考序列时,两个read之间的长度(ISIZE),可理解为测序文库的长度。
第10列:存储的序列(SEQ),没有存储,此列则用“*”标识。该序列的长度一定等于CIGAR标识中“M”,“I”,“S”,“=”,“X”标识的碱基长度之和。
第11列:序列的每个碱基对应一个碱基质量字符(QUAL),每个碱基质量字符对应的ASCII码值减去33(Sanger Phred-33 质量值体系),即为该碱基的测序质量得分(PhredQuality Score)。不同测序质量得分代表不同的碱基测序错误率,如测序质量得分值为20和30分别表示碱基测序错误率为1%和0.1%。
作为优选的实施方式,所述S2,获取待测样品的突变位点信息包括:
获取待测样品的肿瘤位点突变的基因组位置信息、待测样品的DNA体细胞突变数据(该数据通常来自肿瘤外显子测序或者特定基因组套测序)和热点突变,本实施例中待测样品的DNA体细胞突变数据和热点突变作为典型的肿瘤突变数据,采用任意给定的方式获取,均在本发明的保护范围内。其中所述基因组位置信息对应的基因组与所述基因组比对结果文件中的基因组完全一致。
作为优选的实施方式,所述待测样品的突变位点信息由基因突变检测数据或先验知识获得,所述先验知识包括:
(1)肿瘤外显子组测序(WES),最常见的获取肿瘤的体细胞突变信息的技术。根据WES数据计算得到的位点突变信息可以用于鉴定scRNA-Seq肿瘤细胞类群;或者特定基因组套(panel)测序,特定基因组套(panel)测序是高通量基因检测和基因测序发展起来后用的一个词语,它是指在检测中不只是检测一个位点、一个基因。而是同时检测多个基因、多个位点。这些位点和基因需要按照一个标准进行选择和组合,从而构成一个检测组套(panel),因此基因检测组套(panel)可以理解为基因组合、基因集合;特定基因组套(panel)测序是一个基因组合,在基因检测中使用基因组套(panel)所检测的基因比单一的位点要多,比PCR技术检测的序列要长,相对来说,获得的基因信息量要多一些且更为全面;
(2)在没有基因突变检测数据时,利用公共数据库中一些“热点突变”也可以在一定程度帮助鉴定单细胞转录组测序数据中的肿瘤细胞类群。目前癌症基因组图谱(TCGA,The Cancer Genome Atlas)和肿瘤体细胞突变数据库(COSMIC)等数据库(包括但不限于TCGA和COSMIC数据库,本领域技术人员可以根据需要选择收录癌症样本体细胞突变信息的其他数据库,均在本发明的保护范围内)收录了很多癌症样本体细胞突变(somaticmutation)的信息,通过这些数据的信息,获知有些肿瘤存在“热点突变”,热点突变就是在很多癌症样本该位点都发生了突变,比如胰腺癌的KRAS G12位点突变,有报告称该位点高达90%患者存在突变。因此利用配套的位点突变信息或者一些热点突变信息,可以应用以上方法鉴定单细胞转录组测序数据中肿瘤细胞类群,从而揭示哪些类群细胞是肿瘤细胞,哪些类群细胞是正常细胞,以及每个类群内的肿瘤细胞突变谱。
(3)其他已经公开的肿瘤突变数据,例如在文章或者数据库中已经公开的。
作为优选的实施方式,所述S3,基于所述第一数据以及突变位点信息进行突变位点的突变分析以及所述肿瘤细胞类群的鉴定包括:
S31,基于所述基因组比对结果文件进行碱基矫正,以矫正测序产生的噪音,从而精准分析突变位点的突变情况,包括:分析所述基因组比对结果文件中每一个细胞在突变位点位置的比对情况,将单细胞转录组测序数据中具有相同唯一分子标签(UMI)的测序读长(reads)片段聚成同一唯一分子标签(UMI)簇;对聚成同一唯一分子标签(UMI)簇的多个测序读长(reads)片段,在突变位点位置的所述比对情况进行判断,所述判断按照碱基矫正程序进行:
(1)如果唯一分子标签(UMI)簇中多个所述测序读长(reads)片段全部为相同的碱基,则确定该唯一分子标签(UMI)簇为突变碱基;
(2)如果唯一分子标签(UMI)簇中多个所述测序读长(reads)片段包括不同的碱基,且其中比例最大的碱基占比超过80%,则该唯一分子标签(UMI)簇为最大比例碱基;
(3)如果唯一分子标签(UMI)簇中多个所述测序读长(reads)片段包括不同的碱基,且其中比例最大的碱基占比低于80%,则放弃使用该唯一分子标签(UMI)簇的信息;
针对所有唯一分子标签(UMI)簇按照所述碱基矫正程序(1)-(3)顺次进行判断,从而基于每个细胞的所有唯一分子标签(UMI)簇进行碱基矫正后获得矫正结果;
S32,基于矫正结果进行突变位点的突变分析,包括:确定参考基因,如果矫正结果中存在唯一分子标签(UMI)簇与参考碱基不一致,则确定该细胞在该突变位点存在突变;或确定多个突变碱基,分别统计每个突变位点上所述多个突变碱基的唯一分子标签(UMI)数量,如果任意一个突变碱基的唯一分子标签(UMI)数量大于0,则确定在所述突变位点存在细胞突变;
S33,基于所述突变位点的突变分析进行所述肿瘤细胞类群的鉴定,确定单细胞分析细胞的所有类群中包括哪些类群和细胞类型,所有细胞类型中哪些为免疫细胞或非免疫细胞,哪些非免疫细胞中的哪些细胞突变成了肿瘤细胞,进而确定哪些细胞类群为肿瘤细胞类群。
作为优选的实施方式,所述S4,基于所述肿瘤细胞类群的鉴定以及突变位点的突变分析获得所述待测样品的所述肿瘤细胞类群统计信息包括:
S41,统计待测样品中每一个细胞类群中,携带突变位点的细胞数量;
S42,基于所述携带突变位点的细胞数量确定每一个细胞类群中所携带的肿瘤细胞数量以及类群特异的肿瘤细胞突变谱。由于肿瘤的异质性,同一位点的非免疫细胞类别(cluster)里面存在不同突变的情形。探索肿瘤细胞的异质性对于肿瘤的靶向治疗、免疫治疗和疾病认识十分重要。
作为优选的实施方式,由于单细胞转录组测序技术限制,有的位点可能未覆盖到,或者覆盖率很低,因此所述S3与S4之间还包括:对所述突变位点进行多次统计检验,控制所述突变位点的假阴性和假阳性结果的产生,包括:
对于每个突变位点,在所述突变位点对应的基因中,随机选择N个位点数据;
对于所述N个位点顺序执行所述S1-S3获得所述N个位点的突变情况,作为所述N个位点的背景值;
基于所述背景值构建背景噪声统计模型;
应用卡方检验(Chi-square test)或费舍尔精确检验(Fisher exact test),基于所述背景噪声统计模型以及所述N个位点的突变情况,排除多个干扰和错误,包括:
计算所述N个位点的突变情况统计显著性,排除测序错误产生的干扰;
基于所述N个位点的非免疫细胞突变情况,排除核糖核酸编辑(RNA edit)产生的干扰;
统计各个细胞类群所有的突变位点的突变频率,基于费舍尔精确检验,排除一些其他的因素产生的干扰,如建库过程中聚合酶链式反应(PCR)产生的错误。脱氧核糖核酸(DNA)建库方法是分子生物学的实验原理,其本质就是在待测片段两端加上接头的过程。目前脱氧核糖核酸(DNA)建库方法按照接头连接方式不同可分为五类:TA克隆连接接头建库、Swift法建库、转座酶法建库、聚合酶链式反应(PCR)扩增子建库、平末端连接接头建库,聚合酶链式反应(PCR)扩增子建库是捕获建库的一种,适用于临床背景下靶向基因的研究;
为了进一步确定S33中肿瘤细胞类群,将免疫细胞对应的细胞类群合并,通过费舍尔精确检验比较各个位点非免疫细胞类群和免疫细胞类群的突变细胞比例P,以P值小于第一阈值的为肿瘤细胞类群备选集合;第一阈值本实施例设定为0.05,当然本领域技术人员可以根据需要设定适当的阈值范围。统计每个所述肿瘤细胞类群备选集合中所述突变细胞比例P值小于第一阈值的位点数目,基于费舍尔精确检验,确定非免疫细胞比例高的为最终的肿瘤细胞类群。
本实施例采用胰腺癌组织单细胞RNA测序数据以及若干突变位点的数据进行进一步的说明。本优选实施例的方法鉴定胰腺癌单细胞转录组测序中的肿瘤细胞。
数据获取,获取待测样品WES 24个体细胞突变信息,scRNA-Seq 10x cellranger的bam文件,barcode文件和cluster信息文件,通过bam文件中CB tag和UB tag可以知道该测序读长(read)来自哪个细胞(cell)的唯一分子标签(UMI)簇测序读长(read),bam文件数据部分截图4所示;
碱基矫正,过滤掉比对质量低的测序读长(read)(过滤bam文件第五列小或等于0read),过滤掉FLAG(bam文件第二列)为256(非初始对齐,not primary alignment),FLAG为2048(补充对齐,supplementary alignment),NH tag大于1,这些考虑为多比对read,过滤掉FLAG为512(读长失效平台/样品质量核查,read fails platform/vendor qualitychecks),这些测序读长(read)考虑为低质量测序读长(read)。剩余的测序读长(read)分析用来分析每一个细胞在位点突变位置的比对情况,相同的唯一分子标签(UMI)的测序读长(read)聚成一簇,对该簇的测序读长(read)该位置比对情况进行分析,如果全部是相同的碱基那么该唯一分子标签(UMI)簇为该碱基,如果该位置存在不同的碱基,且最大比例的碱基占比超过80%,那么该唯一分子标签(UMI)簇为该最大比例碱基,否则弃掉该唯一分子标签(UMI)簇。所有的唯一分子标签(UMI)簇都以这个规则分析。下表是其中的一个示例,CAAGGCCCATGAACCT-1(细胞barcode)中的ACTTTGTCCT(分子唯一标签(UMI))族矫正为A碱基。
突变分析,对每一个细胞所有的唯一分子标签(UMI)簇的碱基进行分析,统计突变碱基的唯一分子标签(UMI)数量,如果有一个突变碱基唯一分子标签(UMI)数量大于0,那么认为该细胞该位点突变存在突变,下表是其中的一个示例,是待测样本单细胞转录组测序所有细胞在KRAS G12(hg19,chr12: 25398284)位点的分析情况,第一列是单个细胞(barcode),第2列是参考碱基(reference)检测情况,‘C.1’,表示参考碱基是C并且有1个唯一分子标签(UMI)簇检测到,第3列是突变碱基(alle)检测情况,‘A.1’,表示突变碱基是A并且有1个唯一分子标签(UMI)簇检测到。第4列是细胞类型,mut表示该细胞是突变细胞,wild表示是野生型细胞。
肿瘤类群鉴定,下表是单细胞分析细胞类群和细胞注释的结果,一共14个细胞类群(cluster),6种细胞类型(Epithelial_cells、T_cells、Macrophage、 Tissue_stem_cells、B_cell、Endothelial_cells),其中T_cells、B_cell和Macrophage为免疫细胞,其他为非免疫细胞,在非免疫细胞中不知道哪些细胞突变了变成了肿瘤细胞。因此分析体细胞突变位点的信息有助于判断哪些细胞类群(cluster)是肿瘤类群。
细胞类群(Cluster) | 细胞类型(CellType) |
0 | Epithelial_cells |
1 | Epithelial_cells |
2 | T_cells |
3 | Macrophage |
4 | Epithelial_cells |
5 | Tissue_stem_cells |
6 | Epithelial_cells |
7 | B_cell |
8 | Epithelial_cells |
9 | Epithelial_cells |
10 | Endothelial_cells |
11 | T_cells |
12 | Epithelial_cells |
13 | T_cells |
将单细胞转录组测序数据当作常规转录组测序数据,分析并得到了24个突变位点。分析这24个突变位点各个细胞类群的突变情况,下表是14个细胞类群(14 cluster)的检测到的突变细胞数目与总细胞数目信息,比cluster 0(C_0)在1.24020353(1号染色体24020353)位置突变信息为23|535(25为检测到的突变细胞,535为该位点检测的总细胞数量)。
位置(position) | C_0 | C_1 | C_2 | C_3 | C_4 | C_5 | C_6 | C_7 | C_8 | C_9 | C_10 | C_11 | C_12 | C_13 |
1.24020353 | 23|535 | 2|278 | 1|258 | 1|195 | 4|151 | 1|113 | 6|123 | 1|75 | 1|43 | 0|35 | 0|34 | 0|39 | 1|37 | 1|35 |
11.64888468 | 16|607 | 5|383 | 1|316 | 1|244 | 5|181 | 1|140 | 7|137 | 1|78 | 1|53 | 0|53 | 0|48 | 0|50 | 2|45 | 1|38 |
11.8705604 | 27|546 | 1|277 | 1|279 | 1|194 | 9|164 | 3|125 | 10|130 | 1|76 | 2|51 | 0|40 | 0|36 | 1|44 | 3|36 | 1|31 |
12.12063523 | 12|537 | 8|266 | 1|214 | 1|170 | 3|160 | 3|115 | 7|124 | 1|65 | 0|41 | 1|37 | 0|30 | 1|36 | 0|37 | 0|27 |
12.5264245 | 11|561 | 3|298 | 1|63 | 1|73 | 7|189 | 0|36 | 6|138 | 1|16 | 1|38 | 1|17 | 2|15 | 1|11 | 0|21 | 0|14 |
13.53254183 | 19|164 | 3|41 | 1|25 | 1|25 | 5|45 | 1|18 | 6|39 | 1|5 | 0|1 | 0|6 | 0|4 | 1|4 | 0|6 | 0|3 |
15.72491605 | 16|533 | 3|251 | 1|115 | 1|178 | 4|151 | 2|87 | 6|115 | 0|21 | 0|33 | 0|21 | 0|22 | 1|28 | 2|25 | 1|25 |
17.39775844 | 11|472 | 2|198 | 1|90 | 1|118 | 8|191 | 0|53 | 9|135 | 1|31 | 0|47 | 0|25 | 0|19 | 1|14 | 0|28 | 0|13 |
19.1438874 | 12|496 | 1|205 | 1|247 | 1|182 | 4|146 | 2|97 | 7|121 | 0|73 | 0|41 | 0|35 | 0|36 | 1|40 | 0|33 | 1|33 |
19.50002789 | 10|567 | 6|278 | 1|300 | 1|234 | 1|158 | 2|124 | 4|125 | 1|81 | 1|56 | 1|51 | 0|45 | 0|44 | 1|41 | 1|36 |
19.55899369 | 13|622 | 2|464 | 1|353 | 1|276 | 2|195 | 1|149 | 6|143 | 1|86 | 1|68 | 1|57 | 0|52 | 0|50 | 0|47 | 1|39 |
19.58904497 | 7|411 | 3|193 | 1|126 | 1|85 | 3|153 | 0|75 | 6|124 | 1|56 | 0|23 | 0|21 | 1|23 | 1|19 | 1|22 | 0|21 |
2.232576129 | 13|593 | 1|349 | 1|281 | 1|236 | 2|173 | 1|124 | 8|138 | 1|76 | 0|50 | 0|46 | 0|52 | 1|48 | 0|40 | 1|39 |
20.57607355 | 24|567 | 6|332 | 1|272 | 1|217 | 6|158 | 3|115 | 8|124 | 1|57 | 1|44 | 0|35 | 0|43 | 1|44 | 0|31 | 0|36 |
20.62153146 | 49|539 | 4|254 | 1|129 | 1|121 | 10|136 | 2|66 | 14|115 | 0|20 | 1|38 | 0|17 | 1|21 | 1|20 | 1|33 | 1|28 |
4.159631991 | 18|40 | 3|7 | 1|6 | 1|6 | 6|10 | 1|2 | 6|8 | 1|1 | 0|1 | - | 1|1 | 1|2 | - | 1|2 |
6.133138193 | 14|589 | 2|355 | 1|325 | 1|249 | 5|178 | 1|145 | 8|141 | 0|84 | 0|61 | 0|55 | 0|49 | 1|48 | 0|45 | 1|37 |
6.33240475 | 26|508 | 4|219 | 1|214 | 1|161 | 6|145 | 1|97 | 6|116 | 1|70 | 1|38 | 1|30 | 1|32 | 0|31 | 2|33 | 1|26 |
8.146015174 | 19|559 | 2|318 | 1|269 | 1|211 | 4|170 | 2|118 | 9|138 | 1|73 | 0|52 | 0|46 | 1|46 | 0|41 | 1|41 | 1|32 |
8.98725964 | 33|115 | 1|9 | 1|7 | 1|14 | 7|25 | 0|8 | 7|22 | 1|3 | 1|6 | 2|4 | 1|3 | 0|3 | 0|1 | 1|4 |
8.99057271 | 22|505 | 6|253 | 1|271 | 1|164 | 2|149 | 2|105 | 6|125 | 1|75 | 3|40 | 0|41 | 1|39 | 0|38 | 0|30 | 1|31 |
9.130914528 | 25|617 | 6|440 | 1|139 | 1|185 | 1|119 | 0|66 | 2|119 | 1|33 | 0|57 | 1|43 | 0|28 | 0|29 | 2|49 | 1|32 |
9.19378401 | 6|507 | 0|277 | 1|315 | 1|217 | 0|166 | 1|134 | 6|137 | 1|81 | 1|54 | 0|47 | 1|49 | 0|48 | 0|41 | 1|36 |
12.25398284 | 7|12 | 1|1 | 0|2 | 0|1 | 1|1 | 0|5 | 3|6 | 0|1 | 0|1 | 0|1 | 0|1 | 0|1 | - | 0|1 |
由于单细胞转录组测技术限制,有的位点可能未覆盖到,或者覆盖率很低。为了控制假阴性和假阳性结果,对每个突变位点,该位点对应基因中,随机选取100个位点。对每个选取出的未突变位点,进行上述各个步骤分析。分析该基因随机选取的100个位点的突变情况,作为该位点的背景值。下表是所有位点对应的基因随机取100个位点得到的背景值。“细胞(cells)”列为100个随机位点总共检测到的细胞数,“突变细胞(mutation cells)”列为100个随机位点总共检测到的突变的细胞数,“突变比率(mutation percent)”列为突变细胞数目比例,即突变的背景值。不同位点背景值不一样,因此矫正效应也各不相同。
位置(position) | 细胞(cells) | 突变细胞(mutation cells) | 突变比率(mutation percent) |
20.57607355 | 236037 | 3232 | 0.013692769 |
11.64888468 | 101 | 0 | 0 |
17.39775844 | 1111 | 0 | 0 |
12.25398284 | 34744 | 101 | 0.002906977 |
12.5264245 | 202 | 0 | 0 |
9.130914528 | 30502 | 0 | 0 |
8.98725964 | 22523 | 101 | 0.004484305 |
15.72491605 | 156954 | 1212 | 0.007722008 |
20.62153146 | 183315 | 2323 | 0.012672176 |
4.159631991 | 5454 | 0 | 0 |
2.232576129 | 202 | 0 | 0 |
1.24020353 | 101 | 0 | 0 |
11.8705604 | 101 | 0 | 0 |
19.55899369 | 606 | 0 | 0 |
8.99057271 | 101 | 0 | 0 |
8.146015174 | 15049 | 202 | 0.013422819 |
12.12063523 | 202 | 0 | 0 |
19.50002789 | 2020 | 0 | 0 |
6.133138193 | 4141 | 0 | 0 |
19.1438874 | 202 | 0 | 0 |
6.33240475 | 202 | 0 | 0 |
19.58904497 | 101 | 0 | 0 |
9.19378401 | 86860 | 303 | 0.003488372 |
13.53254183 | 13332 | 0 | 0 |
采用费舍尔精确检验(Fisher exact test)判断该cluster检测到突变细胞与背景位点检测突变细胞比例的差别,下表是费舍尔精确检验矫正后各个cluster突变的细胞信息,很多低频的cluster都被矫正为无突变,例如1.24020353该位点的cluster12 1|37矫正为0|37,cluster13 1|35矫正为0|35。
位置(position) | C_0 | C_1 | C_2 | C_3 | C_4 | C_5 | C_6 | C_7 | C_8 | C_9 | C_10 | C_11 | C_12 | C_13 |
1.24020353 | 23|535 | 0|278 | 0|258 | 0|195 | 0|151 | 0|113 | 6|123 | 0|75 | 0|43 | 0|35 | 0|34 | 0|39 | 0|37 | 0|35 |
11.64888468 | 0|607 | 0|383 | 0|316 | 0|244 | 0|181 | 0|140 | 7|137 | 0|78 | 0|53 | 0|53 | 0|48 | 0|50 | 0|45 | 0|38 |
11.8705604 | 27|546 | 0|277 | 0|279 | 0|194 | 9|164 | 0|125 | 10|130 | 0|76 | 0|51 | 0|40 | 0|36 | 0|44 | 3|36 | 0|31 |
12.12063523 | 12|537 | 8|266 | 0|214 | 0|170 | 0|160 | 3|115 | 7|124 | 0|65 | 0|41 | 0|37 | 0|30 | 0|36 | 0|37 | 0|27 |
12.5264245 | 11|561 | 0|298 | 0|63 | 0|73 | 7|189 | 0|36 | 6|138 | 0|16 | 0|38 | 0|17 | 2|15 | 0|11 | 0|21 | 0|14 |
13.53254183 | 19|164 | 3|41 | 1|25 | 1|25 | 5|45 | 1|18 | 6|39 | 1|5 | 0|1 | 0|6 | 0|4 | 1|4 | 0|6 | 0|3 |
15.72491605 | 16|533 | 0|251 | 0|115 | 0|178 | 4|151 | 0|87 | 6|115 | 0|21 | 0|33 | 0|21 | 0|22 | 0|28 | 2|25 | 0|25 |
17.39775844 | 11|472 | 2|198 | 0|90 | 0|118 | 8|191 | 0|53 | 9|135 | 1|31 | 0|47 | 0|25 | 0|19 | 1|14 | 0|28 | 0|13 |
19.1438874 | 12|496 | 0|205 | 0|247 | 0|182 | 4|146 | 0|97 | 7|121 | 0|73 | 0|41 | 0|35 | 0|36 | 0|40 | 0|33 | 0|33 |
19.50002789 | 10|567 | 6|278 | 0|300 | 0|234 | 0|158 | 2|124 | 4|125 | 1|81 | 1|56 | 1|51 | 0|45 | 0|44 | 1|41 | 1|36 |
19.55899369 | 13|622 | 0|464 | 0|353 | 0|276 | 0|195 | 0|149 | 6|143 | 0|86 | 0|68 | 0|57 | 0|52 | 0|50 | 0|47 | 0|39 |
19.58904497 | 0|411 | 0|193 | 0|126 | 0|85 | 0|153 | 0|75 | 6|124 | 0|56 | 0|23 | 0|21 | 0|23 | 0|19 | 0|22 | 0|21 |
2.232576129 | 13|593 | 0|349 | 0|281 | 0|236 | 0|173 | 0|124 | 8|138 | 0|76 | 0|50 | 0|46 | 0|52 | 0|48 | 0|40 | 0|39 |
20.57607355 | 24|567 | 0|332 | 0|272 | 0|217 | 6|158 | 0|115 | 8|124 | 0|57 | 0|44 | 0|35 | 0|43 | 0|44 | 0|31 | 0|36 |
20.62153146 | 49|539 | 0|254 | 0|129 | 0|121 | 10|136 | 0|66 | 14|115 | 0|20 | 0|38 | 0|17 | 0|21 | 0|20 | 0|33 | 0|28 |
4.159631991 | 18|40 | 3|7 | 1|6 | 1|6 | 6|10 | 1|2 | 6|8 | 1|1 | 0|1 | - | 1|1 | 1|2 | - | 1|2 |
6.133138193 | 14|589 | 2|355 | 0|325 | 0|249 | 5|178 | 1|145 | 8|141 | 0|84 | 0|61 | 0|55 | 0|49 | 1|48 | 0|45 | 1|37 |
6.33240475 | 26|508 | 0|219 | 0|214 | 0|161 | 6|145 | 0|97 | 6|116 | 0|70 | 0|38 | 0|30 | 0|32 | 0|31 | 2|33 | 0|26 |
8.146015174 | 19|559 | 0|318 | 0|269 | 0|211 | 0|170 | 0|118 | 9|138 | 0|73 | 0|52 | 0|46 | 0|46 | 0|41 | 0|41 | 0|32 |
8.98725964 | 33|115 | 1|9 | 1|7 | 0|14 | 7|25 | 0|8 | 7|22 | 1|3 | 1|6 | 2|4 | 1|3 | 0|3 | 0|1 | 1|4 |
8.99057271 | 22|505 | 0|253 | 0|271 | 0|164 | 0|149 | 0|105 | 6|125 | 0|75 | 3|40 | 0|41 | 0|39 | 0|38 | 0|30 | 0|31 |
9.130914528 | 25|617 | 6|440 | 1|139 | 1|185 | 1|119 | 0|66 | 2|119 | 1|33 | 0|57 | 1|43 | 0|28 | 0|29 | 2|49 | 1|32 |
9.19378401 | 6|507 | 0|277 | 0|315 | 0|217 | 0|166 | 0|134 | 6|137 | 0|81 | 0|54 | 0|47 | 0|49 | 0|48 | 0|41 | 0|36 |
12.25398284 | 7|12 | 1|1 | 0|2 | 0|1 | 1|1 | 0|5 | 3|6 | 0|1 | 0|1 | 0|1 | 0|1 | 0|1 | - | 0|1 |
为了进一步确定肿瘤细胞类群(cluster),将免疫细胞的细胞类群(cluster)合并,通过费舍尔精确检验比较各个位点非免疫细胞类群(cluster)和免疫细胞类群(cluster)突变细胞比例,以P值小于0.05为肿瘤细胞类群(cluster)。下表是各个位点非免疫细胞与免疫细胞相比费舍尔精确检验平均P值。
位置(position) | C_0 | C_1 | C_4 | C_5 | C_6 | C_8 | C_9 | C_10 | C_12 |
1.24020353 | 0 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 1 |
11.64888468 | 1 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 1 |
11.8705604 | 0 | 1 | 0 | 1 | 0 | 1 | 1 | 1 | 0.0001 |
12.12063523 | 0.0003 | 0.0002 | 1 | 0.006 | 0 | 1 | 1 | 1 | 1 |
12.5264245 | 0.0478 | 1 | 0.0093 | 1 | 0.0066 | 1 | 1 | 0.0057 | 1 |
13.53254183 | 0.188 | 0.5802 | 0.3039 | 0.7308 | 0.1318 | 1 | 1 | 1 | 1 |
15.72491605 | 0.0002 | 1 | 0.007 | 1 | 0.0002 | 1 | 1 | 1 | 0.0039 |
17.39775844 | 0.0971 | 0.5709 | 0.0157 | 1 | 0.0013 | 1 | 1 | 1 | 1 |
19.1438874 | 0.0001 | 1 | 0.0016 | 1 | 0 | 1 | 1 | 1 | 1 |
19.50002789 | 0.0075 | 0.0085 | 1 | 0.1109 | 0.006 | 0.2077 | 0.1916 | 1 | 0.1582 |
19.55899369 | 0 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 1 |
19.58904497 | 1 | 1 | 1 | 1 | 0.0005 | 1 | 1 | 1 | 1 |
2.232576129 | 0 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 1 |
20.57607355 | 0 | 1 | 0.0001 | 1 | 0 | 1 | 1 | 1 | 1 |
20.62153146 | 0 | 1 | 0 | 1 | 0 | 1 | 1 | 1 | 1 |
4.159631991 | 0.2124 | 0.4285 | 0.124 | 0.5439 | 0.0433 | 1 | - | 0.3333 | - |
6.133138193 | 0.0004 | 0.3896 | 0.0038 | 0.4146 | 0 | 1 | 1 | 1 | 1 |
6.33240475 | 0 | 1 | 0.0001 | 1 | 0 | 1 | 1 | 1 | 0.0037 |
8.146015174 | 0 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 1 |
8.98725964 | 0.0207 | 0.6557 | 0.0767 | 1 | 0.0478 | 0.5236 | 0.0889 | 0.3215 | 1 |
8.99057271 | 0 | 1 | 1 | 1 | 0 | 0.0003 | 1 | 1 | 1 |
9.130914528 | 0.0017 | 0.4083 | 0.7158 | 1 | 0.3978 | 1 | 0.3885 | 1 | 0.123 |
9.19378401 | 0.0055 | 1 | 1 | 1 | 0 | 1 | 1 | 1 | 1 |
12.25398284 | 0.0249 | 0.1429 | 0.1429 | 1 | 0.0909 | 1 | 1 | 1 | - |
统计每个细胞类群(cluster) P值<0.05的位点数目,并用费舍尔精确检验跟非免疫细胞比较,最终cluster0、cluster4、cluster6为肿瘤的细胞类群(cluster)。
细胞类群(cluster) | P值<0.05的位点数目(total_P value_less_0.05_site) | 百分比(percent) | P值(P value) |
0.Epithelial_cells | 19 | 0.791666667 | 7.37E-09 |
1.Epithelial_cells | 2 | 0.083333333 | 0.4894 |
4.Epithelial_cells | 9 | 0.375 | 0.001559 |
5.Tissue_stem_cells | 1 | 0.041666667 | 1 |
6.Epithelial_cells | 21 | 0.875 | 1.81E-10 |
8.Epithelial_cells | 1 | 0.041666667 | 1 |
9.Epithelial_cells | 0 | 0 | 1 |
10.Endothelial_cells | 1 | 0.041666667 | 1 |
12.Epithelial_cells | 3 | 0.125 | 0.234 |
肿瘤类群统计,统计该样品中每一个细胞类群(cluster)中,携带突变位点的细胞数量,从而判断该类群携带多少肿瘤细胞。分别统计cluster0、cluster4、cluster6突变的细胞比例,用这些cluster中24个突变检测到的突变细胞除以24个位点检测到的所有细胞。注:由于测序数据深度通常不够,总是存在假阴性现象,这个比例通常小于真实肿瘤细胞比例。因此,该步骤计算得到比例可以作为真实肿瘤细胞比例的下限值。
细胞类群(Cluster) | 肿瘤细胞数量(Tumor cell) | 总细胞数量(Total cell) | 肿瘤细胞占比(Tumor cell percent) |
0.Epithelial_cells | 265 | 622 | 0.426 |
4.Epithelial_cells | 57 | 197 | 0.289 |
6.Epithelial_cells | 68 | 143 | 0.475 |
细胞类群(cluster)特异的肿瘤细胞突变谱分析,由于肿瘤的异质性,同一位点的细胞类群(cluster)里面存在不同突变的情形。探索肿瘤细胞的异质性对于肿瘤的靶向治疗、免疫治疗和疾病认识十分重要。图5所示是cluster0、cluster4、cluster6突变细胞在所有位点的突变情况,每一列表示为一个细胞,每一行为一个位点,黑色表示该细胞在该位点存在突变,白色表示未突变,灰色表示为未检到该细胞该位点的信息。
实施例二
如图2所示,本实施例提供一种单细胞转录组测序数据中肿瘤细胞类群的鉴定***,包括:
测序数据获取模块101,获取待测样品的单细胞转录组测序数据并基于所述单细胞转录组测序数据的分析获得第一数据;
突变位点获取模块102,用于获取待测样品的突变位点信息;
突变分析及肿瘤细胞类群鉴定模块103,用于基于所述第一数据以及突变位点信息进行突变位点的突变分析以及所述肿瘤细胞类群的鉴定;
统计模块104,用于基于所述肿瘤细胞类群的鉴定以及突变位点的突变分析获得所述待测样品的所述肿瘤细胞类群统计信息。
该***可实现上述实施例一提供的鉴定方法,具体的鉴定方法可参见实施例一中的描述,在此不再赘述。
本发明还提供了一种存储器,存储有多条指令,指令用于实现如实施例一的方法。
如图6所示,本发明还提供了一种电子设备,包括处理器301和与处理器301连接的存储器302,存储器302存储有多条指令,指令可被处理器加载并执行,以使处理器能够执行如实施例一的方法。
尽管已描述了本发明的优选实施例,但本领域内的技术人员一旦得知了基本创造性概念,则可对这些实施例做出另外的变更和修改。所以,所附权利要求意欲解释为包括优选实施例以及落入本发明范围的所有变更和修改。显然,本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样,倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内,则本发明也意图包含这些改动和变型在内。
Claims (10)
1.一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法,其特征在于,包括:
S1,获取待测样品的单细胞转录组测序数据并基于所述单细胞转录组测序数据的分析获得第一数据;
S2,获取待测样品的突变位点信息;
S3,基于所述第一数据以及突变位点信息进行突变位点的突变分析以及所述肿瘤细胞类群的鉴定;
S4,基于所述肿瘤细胞类群的鉴定以及突变位点的突变分析获得所述待测样品的所述肿瘤细胞类群统计信息;
所述第一数据包括:
基因组比对结果文件;
细胞条形码文件;以及
细胞聚类结果;
所述S2,获取待测样品的突变位点信息包括:
获取待测样品的肿瘤位点突变的基因组位置信息、待测样品的脱氧核糖核酸体细胞突变数据和热点突变数据,其中所述基因组位置信息对应的基因组与所述基因组比对结果文件中的基因组完全一致;
所述S3,基于所述第一数据以及突变位点信息进行突变位点的突变分析以及所述肿瘤细胞类群的鉴定包括:
S31,基于所述基因组比对结果文件进行碱基矫正,以矫正测序产生的噪音,从而精准分析突变位点的突变情况,包括:分析所述基因组比对结果文件中每一个细胞在突变位点位置的比对情况,将单细胞转录组测序数据中具有相同唯一分子标签的测序读长片段聚成同一唯一分子标签簇;对聚成同一唯一分子标签簇的多个测序读长片段,在突变位点位置的所述比对情况进行判断,所述判断按照碱基矫正程序进行:
(1)如果唯一分子标签簇中多个所述测序读长片段全部为相同的碱基,则确定所述唯一分子标签簇为突变碱基;
(2)如果唯一分子标签簇中多个所述测序读长片段包括不同的碱基,且其中比例最大的碱基占比超过80%,则所述唯一分子标签簇为最大比例碱基;
(3)如果唯一分子标签簇中多个所述测序读长片段包括不同的碱基,且其中比例最大的碱基占比低于80%,则放弃使用所述唯一分子标签簇的信息;
针对所有唯一分子标签簇按照所述碱基矫正程序(1)-(3)顺次进行判断,从而基于每个细胞的所有唯一分子标签簇进行碱基矫正后获得矫正结果;
S32,基于矫正结果进行突变位点的突变分析,包括:确定参考基因,如果矫正结果中存在唯一分子标签簇与参考碱基不一致,则确定细胞在所述突变位点存在突变;或确定多个突变碱基,分别统计每个突变位点上所述多个突变碱基的唯一分子标签数量,如果任意一个突变碱基的唯一分子标签数量大于0,则确定在所述突变位点存在细胞突变;
S33,基于所述突变位点的突变分析进行所述肿瘤细胞类群的鉴定。
2.根据权利要求1所述的一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法,其特征在于,所述获取待测样品的单细胞转录组测序数据包括:从因美纳公司Bio-Rad单细胞测序方法、BD公司Rhapsody单细胞分析***、10x genomics公司铬单细胞测序方法、ICELL8单细胞制备***和/或C1单细胞制备***获取待测样品的单细胞转录组测序数据。
3.根据权利要求1所述的一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法,其特征在于,所述基因组比对结果文件为bam文件。
4.根据权利要求1所述的一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法,其特征在于,所述待测样品的突变位点信息由基因突变检测数据或先验知识获得,所述先验知识包括:
(1)肿瘤外显子测序或者特定基因组套测序;
(2)公共数据库中记载的热点突变,所述公共数据库包括癌症基因组图谱或肿瘤体细胞突变数据库;
(3)文章或者数据库中已经公开的肿瘤突变数据。
5.根据权利要求1所述的一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法,其特征在于,所述S4,基于所述肿瘤细胞类群的鉴定以及突变位点的突变分析获得所述待测样品的所述肿瘤细胞类群统计信息包括:
S41,统计待测样品中每一个细胞类群中,携带突变位点的细胞数量;
S42,基于所述携带突变位点的细胞数量确定每一个细胞类群中所携带的肿瘤细胞数量以及类群特异的肿瘤细胞突变谱。
6.根据权利要求5所述的一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法,其特征在于,所述S3与S4之间还包括:对所述突变位点进行多次统计检验,控制所述突变位点的假阴性和假阳性结果的产生,包括:
对于每个突变位点,在所述突变位点对应的基因中,随机选择N个位点数据;
对于所述N个位点按照顺序执行所述S1-S3获得所述N个位点的突变情况,作为所述N个位点的背景值;
基于所述背景值构建背景噪声统计模型;
应用卡方检验或费舍尔精确检验,基于所述背景噪声统计模型以及所述N个位点的突变情况,排除多个干扰和错误,包括:
计算所述N个位点的突变情况统计显著性,排除测序错误产生的干扰;
基于所述N个位点的非免疫细胞突变情况,排除核糖核酸编辑产生的干扰;
统计各个细胞类群所有的突变位点的突变频率,基于费舍尔精确检验,排除建库过程中聚合酶链式反应产生的错误;
将免疫细胞对应的细胞类群合并,通过费舍尔精确检验比较各个突变位点非免疫细胞类群和免疫细胞类群的突变细胞比例P,将P值小于第一阈值的确定为肿瘤细胞类群备选集合;统计每个所述肿瘤细胞类群备选集合中所述突变细胞比例P值小于第一阈值的位点数目,基于费舍尔精确检验,确定非免疫细胞比例最高的为最终的肿瘤细胞类群。
7.根据权利要求6所述的一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法,其特征在于,所述第一阈值为0.05。
8.一种单细胞转录组测序数据中肿瘤细胞类群的鉴定***,用于实施根据权利要求1-7任一所述的单细胞转录组测序数据中肿瘤细胞类群的鉴定方法,其特征在于,包括:
测序数据获取模块(101),获取待测样品的单细胞转录组测序数据并基于所述单细胞转录组测序数据的分析获得第一数据;
突变位点获取模块(102),用于获取待测样品的突变位点信息;
突变分析及肿瘤细胞类群鉴定模块(103),用于基于所述第一数据以及突变位点信息进行突变位点的突变分析以及所述肿瘤细胞类群的鉴定;
统计模块(104),用于基于所述肿瘤细胞类群的鉴定以及突变位点的突变分析获得所述待测样品的所述肿瘤细胞类群统计信息。
9.一种电子设备,其特征在于,包括处理器和存储器,所述存储器存储有多条指令,所述处理器用于读取所述指令并执行如权利要求1-7任一所述的鉴定方法。
10.一种计算机可读存储介质,其特征在于,所述计算机可读存储介质存储有多条指令,所述多条指令可被处理器读取并执行如权利要求1-7任一所述的鉴定方法。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210865067.6A CN115083521B (zh) | 2022-07-22 | 2022-07-22 | 一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法及*** |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210865067.6A CN115083521B (zh) | 2022-07-22 | 2022-07-22 | 一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法及*** |
Publications (2)
Publication Number | Publication Date |
---|---|
CN115083521A CN115083521A (zh) | 2022-09-20 |
CN115083521B true CN115083521B (zh) | 2022-11-11 |
Family
ID=83243002
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202210865067.6A Active CN115083521B (zh) | 2022-07-22 | 2022-07-22 | 一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法及*** |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN115083521B (zh) |
Families Citing this family (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN116486913B (zh) * | 2023-05-23 | 2023-10-03 | 浙江大学 | 基于单细胞测序从头预测调控突变的***、设备和介质 |
CN116758994B (zh) * | 2023-07-03 | 2024-02-27 | 杭州联川生物技术股份有限公司 | 区分肿瘤细胞和非肿瘤细胞的基因集、方法、介质及设备 |
Family Cites Families (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US6033674A (en) * | 1995-12-28 | 2000-03-07 | Johns Hopkins University School Of Medicine | Method of treating cancer with a tumor cell line having modified cytokine expression |
WO2019129239A1 (en) * | 2017-12-29 | 2019-07-04 | Act Genomics Co., Ltd. | Method and system for sequence alignment and variant calling |
CN109022553B (zh) * | 2018-06-29 | 2019-10-25 | 裕策医疗器械江苏有限公司 | 用于肿瘤突变负荷检测的基因芯片及其制备方法和装置 |
CN112111565A (zh) * | 2019-06-20 | 2020-12-22 | 上海其明信息技术有限公司 | 一种细胞游离dna测序数据的突变分析方法和装置 |
CN110577983A (zh) * | 2019-09-29 | 2019-12-17 | 中国科学院苏州生物医学工程技术研究所 | 高通量单细胞转录组与基因突变整合分析方法 |
CN111321209A (zh) * | 2020-03-26 | 2020-06-23 | 杭州和壹基因科技有限公司 | 一种用于循环肿瘤dna测序数据双端矫正的方法 |
CN115198003B (zh) * | 2020-11-18 | 2023-07-25 | 鲲羽生物科技(江门)有限公司 | 适用于barcode测序的转录组空间位置信息的检测方法及其应用 |
CN113160887B (zh) * | 2021-04-23 | 2022-06-14 | 哈尔滨工业大学 | 一种融合了单细胞tcr测序数据的肿瘤新生抗原筛选方法 |
-
2022
- 2022-07-22 CN CN202210865067.6A patent/CN115083521B/zh active Active
Also Published As
Publication number | Publication date |
---|---|
CN115083521A (zh) | 2022-09-20 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Sheng et al. | Multi-perspective quality control of Illumina RNA sequencing data analysis | |
Krawitz et al. | Microindel detection in short-read sequence data | |
CN115083521B (zh) | 一种单细胞转录组测序数据中肿瘤细胞类群的鉴定方法及*** | |
US10127351B2 (en) | Accurate and fast mapping of reads to genome | |
CN106068330B (zh) | 将已知等位基因用于读数映射中的***和方法 | |
US20190338349A1 (en) | Methods and systems for high fidelity sequencing | |
Anderson et al. | ReCombine: a suite of programs for detection and analysis of meiotic recombination in whole-genome datasets | |
CN113035273B (zh) | 一种快速、超高灵敏度的dna融合基因检测方法 | |
CN114420212B (zh) | 一种大肠杆菌菌株鉴定方法和*** | |
CN113160882A (zh) | 一种基于三代测序的病原微生物宏基因组检测方法 | |
CN112349346A (zh) | 检测基因组区域中的结构变异的方法 | |
US20140288844A1 (en) | Characterization of biological material in a sample or isolate using unassembled sequence information, probabilistic methods and trait-specific database catalogs | |
Han et al. | Novel algorithms for efficient subsequence searching and mapping in nanopore raw signals towards targeted sequencing | |
CN111321209A (zh) | 一种用于循环肿瘤dna测序数据双端矫正的方法 | |
CN109461473B (zh) | 胎儿游离dna浓度获取方法和装置 | |
CN116064755A (zh) | 一种基于连锁基因突变检测mrd标志物的装置 | |
CN109949866B (zh) | 病原体操作组的检测方法、装置、计算机设备和存储介质 | |
Zhang et al. | MaLAdapt reveals novel targets of adaptive introgression from Neanderthals and Denisovans in worldwide human populations | |
JPWO2019132010A1 (ja) | 塩基配列における塩基種を推定する方法、装置及びプログラム | |
CN112885407B (zh) | 一种基于二代测序的微单倍型检测分型***和方法 | |
CN114990202A (zh) | Snp位点在评估基因组异常的应用及评估基因组异常的方法 | |
CN114974432A (zh) | 一种生物标志物的筛选方法及其相关应用 | |
CN114566214A (zh) | 检测基因组缺失***变异的方法及检测装置和计算机可读存储介质与应用 | |
US20170226588A1 (en) | Systems and methods for dna amplification with post-sequencing data filtering and cell isolation | |
US20220399079A1 (en) | Method and system for combined dna-rna sequencing analysis to enhance variant-calling performance and characterize variant expression status |
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 |