一种CRISPR目标区间变异检测的分析方法

技术领域

本发明涉及基因检测分析技术领域，尤其涉及一种CRISPR目标区间变异检测的分析方法。

背景技术

Clustered regularly interspaced short palindromic repeats(CRISPR)——规律成簇的间隔短回文重复和CRISPR-associated protein 9(Cas9)，共同组成了一套系统，是细菌用来抵御外源DNA入侵的天然防御系统。目前利用这个原理构建了很强的RNA引导的DNA靶向平台——主要用来基因组编辑、转录扰乱、表观遗传调控等。

CRISPR/Cas9从2013年以来，凭借其特有的敲除靶序列效果，迅速风靡科研界。与传统的siRNA文库和shRNA文库比较，有更加严苛的筛选标准，更有效去除假阳性，而且有着成本低廉、构建简单、操作方便的诸多优点。

而现有的CRISPR分析主要是集中在针对目标区间单位点碱基突变的分析，而忽略了完整序列在成功转导后的比率，从而造成假阳性的情况。

发明内容

本发明的提供一种CRISPR目标区间变异检测的分析方法。

本发明的方案是：

一种CRISPR目标区间变异检测的分析方法，包括下列步骤：

1)将高质量测序数据使用flash软件将两条reads按照overlap进行合并，再使用cd_hit软件进行聚类，获得聚类结果与每一类代表性序列；

2)使用mafft软件将聚类的代表性序列与目标区间进行序列的多重比对，扫描目标区间，根据扫描结果查找出突变；

3)根据步骤2)中找出的突变进行分类，再结合步骤1)中的聚类结果中每一个cluster的reads数除以总聚类的reads数，得到每个cluster条目对应分类的突变率；

4)将步骤1)中的聚类结果的reads数与步骤2)中的扫描结果统计每种cluster的丰度，然后结合步骤3)中突变率，筛选，获得分析结果。

作为优选的技术方案，步骤1)中的所述聚类为当长序列包含短序列时把短序列归为长序列。

作为优选的技术方案，所述步骤1)中高质量测序数据获得方法：

将原始二代下机数据按照barcode拆分成单个样本的fastq文件，将得到的所述fastq文件使用fastp软件过滤去接头处理得到高质量测序数据。

作为优选的技术方案，所述步骤1)中高质量测序数据使用flash软件将两条reads根据overlap合并成一条长reads；将得到的所述长reads的fastq文件转化为fasta格式的文件，使用cd-hit软件对转化后的所述fasta格式的文件内序列进行聚类，所述聚类结果为各个cluster的reads数以及各个cluster的代表性序列；

作为优选的技术方案，使用mafft软件将聚类结果中的各个cluster的代表性序列与目标区间进行序列进行两两比对。

作为优选的技术方案，所述步骤4)中筛选为筛选高频cluster与目标区间高频的突变位点。

由于采用了上述技术方案一种CRISPR目标区间变异检测的分析方法，包括下列步骤，1)将高质量测序数据使用flash软件将两条reads按照overlap进行合并，再使用cd_hit软件进行聚类，获得聚类结果与每一类代表性序列；2)使用mafft软件将聚类的代表性序列与目标区间进行序列的多重比对，扫描目标区间，根据扫描结果查找出突变；3)根据步骤2)中找出的突变进行分类，再结合步骤1)中的聚类结果中每一个cluster的reads数除以总聚类的reads数，得到每个cluster条目对应分类的突变率；4)将步骤1)中的聚类结果的reads数与步骤2)中的扫描结果统计每种cluster的丰度，然后结合步骤3)中突变率，筛选，获得分析结果。

本发明的优点：

本发明避免了因忽略完整序列在成功转导后的比率造成假阳性的情况出现，成本低廉、构建简单、操作方便，通过本发明的分析得到的结果更接近真实情况，提高了准确率。

附图说明

图1为本发明的框架流程图；

图2为本发明实施例2的统计比对率结果图；

图3为本发明实施例2的比对结果扫描突变位点图。

具体实施方式

为了弥补以上不足，本发明提供了一种CRISPR目标区间变异检测的分析方法以解决上述背景技术中的问题。

一种CRISPR目标区间变异检测的分析方法，包括下列步骤：

1)将高质量测序数据使用flash软件将两条reads按照overlap进行合并，再使用cd_hit软件进行聚类，获得聚类结果与每一类代表性序列；

2)使用mafft软件将聚类的代表性序列与目标区间进行序列的多重比对，扫描目标区间，根据扫描结果查找出突变；

3)根据步骤2)中找出的突变进行分类，再结合步骤1)中的聚类结果中每一个cluster的reads数除以总聚类的reads数，得到每个cluster条目对应分类的突变率；

4)将步骤1)中的聚类结果的reads数与步骤2)中的扫描结果统计每种cluster的丰度，然后结合步骤3)中突变率，筛选，获得分析结果。

步骤1)中的所述聚类为当长序列包含短序列时把短序列归为长序列。

所述步骤1)中高质量测序数据获得方法：

将原始二代下机数据按照barcode拆分成单个样本的fastq文件，将得到的所述fastq文件使用fastp软件过滤去接头处理得到高质量测序数据。

所述步骤1)中高质量测序数据使用flash软件将两条reads根据overlap合并成一条长reads；将得到的所述长reads的fastq文件转化为fasta格式的文件，使用cd-hit软件对转化后的所述fasta格式的文件内序列进行聚类，所述聚类结果为各个cluster的reads数以及各个cluster的代表性序列；

使用mafft软件将聚类结果中的各个cluster的代表性序列与目标区间进行序列进行两两比对。

所述步骤4)中筛选为筛选高频cluster与目标区间高频的突变位点。

为了使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体实施例，进一步阐述本发明。

实施例1：

1)将高质量测序数据使用flash软件将两条reads按照overlap进行合并，再使用cd_hit软件进行聚类，获得聚类结果；

2)使用mafft软件将聚类结果与目标区间进行序列的多重比对，扫描目标区间，根据扫描结果查找出突变；

3)根据步骤2)中找出的突变进行分类，再结合步骤1)中的聚类结果中每一个cluster的reads数除以总聚类的reads数，得到每个cluster条目对应分类的突变率；

4)将步骤1)中的聚类结果的reads数与步骤2)中的扫描结果统计每种cluster的丰度，然后结合步骤3)中突变率，筛选，获得分析结果。

作为优选的技术方案，步骤1)中的所述聚类为当长序列包含短序列时把短序列归为长序列。

所述步骤1)中高质量测序数据获得方法：

将原始二代下机数据按照barcode拆分成单个样本的fastq文件，将得到的所述fastq文件使用fastp软件过滤去接头处理得到高质量测序数据。

所述步骤1)中高质量测序数据使用flash软件将两条reads根据overlap合并成一条长reads；将得到的所述长reads的fastq文件转化为fasta格式的文件，使用cd-hit软件对转化后的所述fasta格式的文件内序列进行聚类，所述聚类结果为各个cluster的reads数以及各个cluster的代表性序列；

使用mafft软件将聚类结果中的各个cluster的代表性序列与目标区间进行序列进行两两比对。

所述步骤4)中筛选为筛选高频cluster与目标区间高频的突变位点。

实施例2

1、以小鼠特定目标基因扩增为例，使用传统的bwa比对然后统计目标区域的突变会发现在前130bp序列中有8个高频突变位点，步骤如下：

一)将test样本高质量的数据使用flash进行搭载：

flash test_HQ_R1.fq.gz test_HQ_R2.fq.gz-m 10-M 300-c-q-x 0.1-t

5|pigz-9-p 2>test.merged.fastq.gz

二)将搭载后的reads使用bwa软件与参考目标序列进行比对：

bwa mem-t 4-R

'@RG\tID:'${name}'\tSM:'test'\tPL:illumina\tLB:PE400\tPU:unit'

ref.fasta test.merged.fastq.gz>test.sam

三)使用picard将sam文件转化为bam并进行排序生成索引：

/picard.jar SortSam INPUT＝test.sam OUTPUT＝test.sorted.bam

SORT_ORDER＝coordinate CREATE_INDEX＝true

TMP_DIR＝./temp/picardSortSam/test

VALIDATION_STRINGENCY＝LENIENT

四)使用iDES软件中的bam2freq.pl脚本将比对bam文件转化文本文件，统计比对率。结果如图1所述：

Perl bam2freq.pl-a-q 0test.sorted.bam ref.fasta test.bed

2、使用聚类后的reads与参考序列进行多重比对，从最多的前8类聚类reads结果可以看出，在前130bp序列中有一些测序结果中有效reads发生突变位点在32，其他位置的突变是多种reads共同作用结果包括可能被放大的随机错误，步骤如下：

一)将test样本高质量的数据使用flash进行搭载：

flash test_HQ_R1.fq.gz test_HQ_R2.fq.gz-m 10-M 300-c-q

-x 0.1-t 5|pigz-9-p 2>test.merged.fastq.gz

二)将fastq文件转化为fasta文件：

Zcat test.merged.fastq.gz|awk'{if(NR％4＝＝1||

NR％4＝＝2)print$0}'test.fastq|sed's/@/>/'>test.fasta

三)使用cd-hit对test.fasta文件进行聚类(按照100％相似)：

cd-hit-est-i tes.fasta-c 1-o PSN19CV922-1-aS 1

四)将test样本的各个cluster代表序列与ref.fasta进行两辆比对(cluster代表序列用sequence_n.fa表示n代表类别)：

for i in sequence_*.fa；do mafft--globalpair--maxiterate

1000--clustalout$i>$i.out；done

五)依据比对结果扫描突变位点，如图2所示：

六)结合聚类reads数与扫描结果统计每种cluster的丰度；

七)从结果中可以看出在前130bp序列中有一些测序结果中有效reads发生突变位点在32，其他位置的突变是多种reads共同作用结果包括可能被放大的随机错误；

八)将cluster的代表性reads在32好位点发生突变的丰度汇总，统计整体发生突变的碱基比率。

以上显示和描述了本发明的基本原理、主要特征及本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。