一种单细胞ATAC-seq数据分析方法
专利摘要:本发明提供一种单细胞ATAC-seq数据分析方法,包括以下步骤:步骤S1,对测序原始数据进行数据分析与质控;步骤S2,比对分析;步骤S3,插入片段分析;步骤S4,富集区域Peak分析;步骤S5,单细胞亚群分类;步骤S6,对Peak相关基因进行注释和富集;步骤S7,TF-motif分析;步骤S8,亚群可及性差异分析;步骤S9,差异可及性位点相关基因分析,对鉴定出的差异TF-motif所在peak区域最邻近的转录起始位点所对应的基因注释等步骤。本发明构建了一个全面、分析内容丰富的单细胞ATAC-seq数据分析流程,分析结果揭示了大量的生物信息,方便人们深入挖掘蕴藏在单细胞水平内的生物学现象和特征,分析流程及结果以html的形式进行可视化展示,分析内容层次明了,结果展现形式多样,增加了报告的可读性。
专利说明:一种单细胞ATAC-seq数据分析方法 技术领域 本发明属于生物信息学技术领域,具体涉及一种单细胞ATAC-Seq数据的生物分析技术领域。 背景技术 众所周知,大多数基因组中的染色质是紧密盘绕在细胞核内的,仅有一小部分区域比较松散,这部分无核小体的裸露DNA区域称为开放染色质(open chromatin),DNA复制和基因转录往往在这些区域发生。ATAC-seq(Assay for Transposase-AccessibleChromatin with high-throughput sequencing,染色质易开放区测序)是利用Tn5转座酶切割染色质开放区,并加上测序引物进行高通量测序,通过生物信息分析鉴定转录因子结合位点和核小体区域位置,从而为研究基因调控、DNA印记等提供有效的方法。 目前,单细胞技术迅速发展,该技术能够很好地解决群体细胞异质性这一难题,因此在各组学中都得到了充分地应用。2013年,Tang F等人首次在单细胞内应用RRBS-seq技术,开启了单细胞表观遗传组学研究的大门。2015年,Buenrostro等人发表了单细胞ATAC-seq技术,即一种对转座酶易感染色质进行高通量测序从而绘制出个体细胞的易感基因组的方法。从成千上万个细胞的总体上紧密类似的图谱中获得来自数百个单细胞ATAC-seq图谱,能够清楚地了解细胞间的变异,发现变异不仅与特定的反式因子和顺式元件系统密切相关,而且也是与细胞变异诱导或者抑制相关的反式因子的组合。单细胞ATAC-seq技术可以在解决细胞异质性的同时,研究细胞转录调控的机制,成为了单细胞表观遗传学的一大突破。目前主要应用在肿瘤异质性研究、基因调控网络分析、细胞谱系示踪和生物标记物发现等方面。 然而,目前对于单细胞表观遗传组测序数据的分析流程还尚未有统一的标准,公开号为CN 107368701A的中国专利公开了一种大批量单细胞ATAC-seq数据质量控制和分析方法,包括对测序片段水平和多细胞水平的质量控制、单细胞层面的质量控制、细胞聚类和细胞特异峰的探测,最终给用户提供一份质量控制的报告文档。目前,还缺乏比较全面、个性化的能够揭示更详细信息的单细胞ATAC-seq数据分析流程,并且分析报告通过html的形式进行展示,这样更加便捷,条目清晰,同时,报告中设置有超链接,便于用户能够对报告进行透彻理解,增加了报告的可读性。 发明内容 本发明的目的在于克服上述现有技术的不足之处而提供一种单细胞ATAC-seq数据分析方法。 为实现上述目的,本发明采取的技术方案为: 一种单细胞ATAC-seq数据分析方法,包括以下步骤: S1.对测序原始数据进行数据分析与质控,得到高质量的数据进行后续分析; S2.比对分析,将修剪后的reads比对到参考基因组,分析确定唯一比对且MAPQ>30的reads数量,并去除比对到线粒体上的reads,将过滤后的reads用于后续分析; S3.插入片段分析,利用两端reads的比对信息可以计算出各样品的文库插入片段长度,从而绘制出插入片段分布图; S4.Peak分析,主要包括样本有效细胞信息结果统计、Peak区域片段分布统计及最终peak-barcode矩阵分析; S5.单细胞亚群分类,对数据进行LSA降维处理,然后进行k-medoids聚类,最后将转换后的矩阵利用t-SNE进行可视化。 S6.对Peak相关基因注释和富集; S7.TF-motif分析,主要包括TF-motif鉴定和TF-motif显著程度分析。 S8:亚群可及性差异分析通过差异TF-motif统计,差异TF-motif聚类以及差异TF-motif在各亚群之间的分布,获得各个差异TF motif在各细胞亚群中的富集分布。 S9:差异可及性位点相关基因分析,对鉴定出的差异TF-motif所在peak区域最邻近的转录起始位点(peak位于该TSS上游1000bp或下游100bp区域)所对应的基因注释。 进一步地,所述步骤S1的数据分析及质控是通过以下方法实现:对测序原始数据中出现错误的barcodes进行过滤和校正,将每个barcodes序列与数据库中已知的barcodes序列进行比对,查找与已知barcodes碱基错配≤2bp的barcodes,并根据读barcodes的丰度和错配碱基的质量值进行评分,分值大于90%的barcodes被认为是正确的barcodes。 进一步地,所述步骤S2比对分析包括比对信息、测序饱和度分析以及TTS周围信号分布统计。 进一步地,所述比对信息包括以下步骤:统计结果以表格的形式进行展示,结果列明了样本名称、比对到目标区域(转录起始位点、DNase超敏感区域、增强子或启动子区域)片段占所有片段的百分比、比对到转录起始位点片段占所有片段的百分比等内容。 进一步地,所述测序饱和度分析包括以下步骤:分析结果以饱和曲线图进行表示,横坐标为每个细胞平均reads数,纵坐标为每个细胞唯一比对的平均片段数。 进一步地,所述TTS周围信号分布统计包括以下步骤:统计TSS上下游各1000bp窗口内的信号强度,每个位点碱基存在于切割位点的片段数目定义为每个位点的信号值,标准化后绘制检测样本TTS周围信号分布图。 进一步地,所述步骤S6Peak相关基因注释,利用bedtools将peak最邻近的转录起始位点所对应的基因作为该peak相关基因。 进一步地,所述步骤S6Peak相关基因富集包括GO富集分析和KO富集分析。 进一步地,所述步骤S9包括差异可及性位点相关基因注释、差异可及性位点相关基因GO功能富集分析和差异可及性位点相关基因KO功能富集分析。 本发明的有益效果:本发明围绕单细胞ATAC-seq数据进行研究,构建了一个全面、分析内容丰富的单细胞ATAC-seq数据分析流程。分析结果揭示了大量的生物信息,进一步方便人们深入挖掘蕴藏在单细胞水平内的生物学现象和特征。最终分析流程及结果以html的形式进行可视化展示,分析内容层次明了,结果展现形式多样,增加了报告的可读性。 附图说明 图1为本发明单细胞ATAC-seq数据分析流程。 具体实施方式 为了更加简洁明了的展示本发明的技术方案、目的和优点,下面结合具体实施例及其附图对本发明做进一步的详细描述。 实施例 本实施例单细胞ATAC-seq数据分析流程,包括如下步骤: 步骤S1:对测序原始数据进行数据分析与质控,得到高质量的数据进行后续分析。主要是利用Cell Ranger软件对测序原始数据中出现错误的barcodes进行过滤和校正。将每个barcodes序列与数据库中已知的barcodes序列进行比对,查找与已知barcodes碱基错配≤2bp的barcodes,并根据读barcodes的丰度和错配碱基的质量值进行评分,分值大于90%的barcodes被认为是正确的barcodes。 步骤S2:比对分析,使用cutadapt识别reads末端引物反向互补序列,并从reads序列中去除,然后,使用BWA-MEM将修剪后的reads比对到参考基因组,随后利用Duplication分析确定唯一比对且MAPQ>30的reads数量,并去除比对到线粒体上的reads,将过滤后的reads用于后续分析。该分析主要包括了以下3部分。 步骤S2.1:比对信息统计,统计结果以表格的形式进行展示,结果列明了样本名称、比对到目标区域(转录起始位点、DNase超敏感区域、增强子或启动子区域)片段占所有片段的百分比、比对到转录起始位点片段占所有片段的百分比等内容。 步骤S2.2:测序饱和度分析,分析结果以饱和曲线图进行表示,横坐标为每个细胞平均reads数,纵坐标为每个细胞唯一比对的平均片段数。 步骤S2.3:TTS周围信号分布统计,统计TSS上下游各1000bp窗口内的信号强度,每个位点碱基存在于切割位点的片段数目定义为每个位点的信号值,标准化后绘制检测样本TTS周围信号分布图。横坐标是相对于TSS的位置信息,纵坐标是相对信号强度。这个图有助于评估文库的信噪比,因为与基因组的基因间区和内含子区相比,TSSs及其周围的启动子区具有较高的染色质易接近性。 步骤S3:插入片段分析,由于Tn5转座酶优先攻击染色质开放区,一般来说,大多数DNA都是不包含或只包含一个核小体的短片段,同时也有一些包含多个核小体的长片段,在含量分布上会呈现明显的片段分布特征。利用两端reads的比对信息可以计算出各样品的文库插入片段长度,从而绘制出插入片段分布图。 步骤S4:富集区域Peak分析,主要包括样本有效细胞信息结果统计、Peak区域片段分布统计及最终peak-barcode矩阵分析,得出差异可及性位点,具体步骤如下: 步骤S4.1:样本有效细胞信息结果统计,首先对peak区域进行鉴定,然后利用peak区域片段数进行有效细胞鉴定,将最终筛选得到的有效细胞进行基本信息统计。其中,peak区域鉴定,利用片段末端位点计算基因组中每个碱基对的转座事件的发生数目,设置401bp的滑动窗口在基因组上进行滑窗,并拟合类ZINBA混合模型,设置1/5比值比作为阈值,大于该阈值的区域被认为是peak信号(富含开放染色质),小于阈值的则被认为是背景噪声。因此,并非所有的酶切位点都位于peak区域内。最后将相邻500bp内的peak合并,得到最终peak区域。接着,利用鉴定的peak区域中的片段count值,拟合两个负二项分布的混合模型来捕捉细胞信号和背景噪声,从而达到区别有效细胞barcode与非细胞barcode,过滤掉非细胞barcode数据后用于后续分析。根据样本有效细胞信息结果绘制出有效细胞回收鉴定图和有效细胞片段分布图。 步骤S4.2:Peak区域片段分布统计,根据peak分析绘制得到每个细胞peak区域所含片段数分布散点图。 步骤S4.3:peak-barcode矩阵分析,根据peak和barcode信息,最终生成一个peak-barcode矩阵。 步骤S5:单细胞亚群分类,首先对数据进行LSA降维处理,然后进行k-medoids聚类,最后将转换后的矩阵利用t-SNE进行可视化。包括单细胞亚群分类信息统计和单细胞亚群分类可视化展示。其中,单细胞亚群分类信息统计结果用统计表和绘制各亚群细胞数量柱状图进行展示。 步骤S6:Peak相关基因注释和富集: 步骤S6.1:Peak相关基因注释主要包括:Peak相关基因注释,利用bedtools,将peak最邻近的转录起始位点(peak位于该TSS上游1000bp或下游100bp区域)所对应的基因作为该peak相关基因。对该样本的peak相关基因进行统计,结果用peak相关基因注释表进行展示。 步骤S6.2:GO富集分析,对上述的peak相关基因进行GO富集分析,富集分析结果一般用柱形图、富集气泡图和条形图进行展示。 步骤S6.3:KO富集分析,对上述的peak相关基因进行KO富集分析,KEGG是有关Pathway的主要公共数据库。Pathway显著性富集分析以KEGG Pathway为单位,应用超几何检验,找出与整个基因组背景相比,在基因中显著性富集的Pathway。通过Pathway显著性富集能确定基因参与的最主要生化代谢途径和信号转导途径。分析结果通常是利用P值(或Q值)最小的前20个pathway来作图,绘制KO富集气泡图和KO富集条形图。 步骤S7:TF-motif分析,主要包括TF-motif鉴定和TF-motif显著程度分析。其中,TF-motif鉴定首先从JASPAR数据库获得TF motif的位置权重矩阵(PWM),再利用MOODS扫描每个peak,查找每个peak与TF匹配的motif(默认p值≤1E-5)。TF-motif显著程度分析主要是针对peak区域中结合TF的motif,分析其富集的显著性程度,主要是通过计算每个细胞内特定TF motif的开放区域reads数目与总开放区域reads数目的比值,并进行深度均一化与Z-Score标准化。 步骤S8:亚群可及性差异分析,主要包括以下3个步骤: 步骤S8.1:差异TF-motif统计,分别计算特定亚群及其剩余亚群中每个TF motif所在peak区域内的酶切位点数目的平均值,利用edgR进行差异分析,并进行Benjamini-Hochberg多重检验校正,默认取亚群中FDR值Top100的TF motif作为该亚群的差异TFmotif。 步骤S8.2:差异TF-motif聚类,结果以样本中各亚群显著TF-motif聚类热图进行呈现。 步骤S8.3:差异TF-motif在各亚群之间的分布,根据TF motif的富集分数,利用表达分布热图展示各个差异TF motif在各细胞亚群中的富集分布。 步骤S9:差异可及性位点相关基因分析,差异可及性位点相关基因是指利用bedtools,对鉴定出的差异TF-motif所在peak区域最邻近的转录起始位点(peak位于该TSS上游1000bp或下游100bp区域)所对应的基因注释。主要包括差异可及性位点相关基因注释、差异可及性位点相关基因GO功能富集分析和差异可及性位点相关基因KO功能富集分析3部分。 以上所述实施例仅表达了本发明的几种实施方式,其描述较为具体和详细,但并不能因此而理解为对发明专利范围的限制。应当指出的是,对于本领域的普通技术人员来说,在不脱离本发明构思的前提下,还可以做出若干变形和改进,这些都属于本发明的保护范围。因此,本发明专利的保护范围应以所附权利要求为准。
权利要求:
1.一种单细胞ATAC-seq数据分析方法,其特征在于,包括以下步骤:
S1.对测序原始数据进行数据分析与质控,得到高质量的数据进行后续分析;
S2.比对分析,将修剪后的reads比对到参考基因组,分析确定唯一比对且MAPQ>30的reads数量,并去除比对到线粒体上的reads,将过滤后的reads用于后续分析;
S3.插入片段分析,利用两端reads的比对信息可以计算出各样品的文库插入片段长度,从而绘制出插入片段分布图;
S4.Peak分析,主要包括样本有效细胞信息结果统计、Peak区域片段分布统计及最终peak-barcode矩阵分析;
S5.单细胞亚群分类,对数据进行LSA降维处理,然后进行k-medoids聚类,最后将转换后的矩阵利用t-SNE进行可视化;
S6.对Peak相关基因进行注释和富集;
S7.TF-motif分析主要包括TF-motif鉴定和TF-motif显著程度分析;
S8:亚群可及性差异分析通过差异TF-motif统计,差异TF-motif聚类以及差异TF-motif在各亚群之间的分布,获得各个差异TF motif在各细胞亚群中的富集分布;
S9:差异可及性位点相关基因分析,对鉴定出的差异TF-motif所在peak区域最邻近的转录起始位点所对应的基因注释。
2.如权利要求1所述的分析方法,其特征在于,所述步骤S1的数据分析及质控是通过以下方法实现:对测序原始数据中出现错误的barcodes进行过滤和校正,将每个barcodes序列与数据库中已知的barcodes序列进行比对,查找与已知barcodes碱基错配≤2bp的barcodes,并根据读barcodes的丰度和错配碱基的质量值进行评分,分值大于90%的barcodes被认为是正确的barcodes。
3.如权利要求1所述的分析方法,其特征在于,所述步骤S2比对分析包括比对信息、测序饱和度分析以及TTS周围信号分布统计。
4.如权利要求2所述的分析方法,其特征在于,所述比对信息包括以下步骤:统计结果以表格的形式进行展示,结果列明了样本名称、比对到目标区域片段占所有片段的百分比、比对到转录起始位点片段占所有片段的百分比等内容。
5.如权利要求2所述的分析方法,其特征在于,所述测序饱和度分析包括以下步骤:分析结果以饱和曲线图进行表示,横坐标为每个细胞平均reads数,纵坐标为每个细胞唯一比对的平均片段数。
6.如权利要求2所述的分析方法,其特征在于,所述TTS周围信号分布统计包括以下步骤:统计TSS上下游各1000bp窗口内的信号强度,每个位点碱基存在于切割位点的片段数目定义为每个位点的信号值,标准化后绘制检测样本TTS周围信号分布图。
7.如权利要求1所述的分析方法,其特征在于,所述步骤S6 Peak相关基因注释,利用bedtools将peak最邻近的转录起始位点所对应的基因作为该peak相关基因。
8.如权利要求1所述的分析方法,其特征在于,所述步骤S6Peak相关基因富集包括GO富集分析和KO富集分析。
9.如权利要求1所述的分析方法,其特征在于,所述步骤S9包括差异可及性位点相关基因注释、差异可及性位点相关基因GO功能富集分析和差异可及性位点相关基因KO功能富集分析。
公开号:CN110544509
申请号:CN201910768671.5A
发明人:夏昊强 高川 周煌凯 张羽 陶勇 罗玥 陈飞钦 曾川川
拥有者:广州基迪奥生物科技有限公司
申请日:2019-08-20
公开日:2019-12-06
全文下载