Guideseq分析流程

Guideseq是一种检测体内基因编辑脱靶的测序方法，基本原理是CRISPR系统剪切基因组后，在进行基因组修复时可能会将一段设计好的双链DNA，即tag，插入到基因组中，通过tag设计引物可将基因组断裂位点两侧的序列扩增出来建库，之后进行二代测序，通过Guideseq的分析流程即可鉴定出脱靶位点。
具体细节可以看他们的论文：
GUIDE-seq enables genome-wide profiling of off-target cleavage by CRISPR-Cas nucleases | Nature Biotechnology
下图是Guideseq的原理图：

Guideseq

我就不介绍Guideseq的具体分析流程怎么跑了，可以参考这里aryeelab/guideseq: Analysis pipeline for the GUIDE-seq assay. For news, see: (github.com)，不过目前有2个版本，老版本只支持python2，新版本支持了python3，具体代码有略微差异。

Guideseq的分析共分为7步，分别是Demultiplex、Umitag、Consolidate、Align、IdentifyOfftargetSites、FilterBackgroundSites以及Visualize。
符合Guideseq分析的标准数据包含read1、read2、index1和index2四个文件，其中read1和read2就是测序reads，index1为样本index，长度为8nt，index2为8nt的样本index+8nt的umi。
这里记录一下Guideseq的分析思路，即每一步具体做了什么，内容根据python2版本的代码整理。

1. Demultiplex

首先是Demultiplex，这一步的作用是根据index将测序数据按样本进行拆分，因此，如果你的测序数据是多个样本混合一起上机测序的，可以直接跑Guideseq流程。如果你的数据已经按照样本拆分好了，那么就可以跳过这一步。
这一步的具体操作是：
1） 将index1和index2的[1:8]的序列(这里是0base，因此各取7nt的序列)组合起来作为样本index；
2）根据样本index将测序文件按样本拆分，拆分后的单样本文件reads数必须>=min_reads (默认10000)，否则存放到undetermined文件；
3）这一步未对原始数据进行更改，仅拆分。

2. Umitag

这一步是对拆分好的数据进行umi标记，用于下一步合并PCR重复的reads。具体操作是:
1）在read1和read2的read id最后加上molecular_barcode，其组成是index2[8:16]_read1[0:6]_read2[0:6]；
2）将reads根据molecular_barcode进行排序。

3. Consolidate

此步的目的是合并PCR重复的reads，因为最终要对鉴定的脱靶位点reads数进行统计，所以要消除PCR的效应。具体操作是：
1）首先将molecular_barcode相同的reads划分进同一个bin；
2）对同一个bin中的reads逐个碱基进行合并：
i) 对每个位置的碱基进行统计，只有质量值大于min_qual (默认15)的碱基才被统计，同时记录每个位置A、T、C、G、N的最高质量值；
ii) 每个位置取出现次数最多的那个碱基，若此碱基出现的频率大于min_freq (默认0.9)，则将该位点合并成此碱基，质量值取该碱基在该位置的最高质量值，否则该位点合并为N，质量值为0；
3）readname 变成 "molecular_barcode_n 原始readname中的第二个元素"，n为被合并的reads数。

4.Align

使用bwa mem进行比对。

5.IdentifyOfftargetSites

这一步根据比对后的sam文件，鉴定脱靶位点。
1）取mapq>=50且flag包含128 (PE测序read2)且flag不包含2048 (嵌合比对)的比对结果，获取read的比对起始位点（根据sam文件中TLEN插入片段长度的正负号不同计算方法也不同）；
2）获得所有比对起始位点及其对应的read数，将同一染色体上距离<=10bp的相邻位点归入一个window，对所有window的read count进行统计，并筛选出满足以下条件之一的window：
i）比对到正链的read count和比对到负链的read count均大于0；
ii）以tag_primer1开头的read count和以tag_primer2开头的read count 均大于0；
3）对于通过筛选的window，确定window中总read count最大的位点，以此为中心，前后各延伸25nt，截取基因组序列；
4）将截取的基因组序列与target sequence进行比对，若最优比对的mismatch<=8，则此基因组位点为脱靶位点，以该window中总read count进行量化。

6. FilterBackgroundSites

一般来说，做实验时会有一个阴性对照，这一步的目的是过滤掉阴性对照中检测到的脱靶位点，但实际的代码是将实验组与对照组的交集保存到单独的文件中，在最后的可视化中并没有将这些位点扣除，这里有些奇怪，不知道我是否理解有误。

7. Visualize

最后就是将鉴定到的脱靶位点进行可视化。
放一张Guideseq论文中的可视化图

image.png