Stacks分析RAD-Seq数据的内部原理
-
一，背景知识
1.Restriction enzymes
限制酶是DNA切割酶。 每种酶识别一个或几个靶序列，并在这些序列处或附近切割DNA。 DNA连接酶（DNA ligase）是一种DNA连接酶。 如果两个DNA的末端匹配，则连接酶可以将它们连接起来，形成一个单一的，不间断的DNA分子。这种酶最早在细菌中发现，细菌使用这种机制来切割外部病毒，达到防御目的，因此叫限制酶。
酶切点根据不同的酶的特性，会产生粘性端和钝端等情况。
2.RAD-Seq
RAD-Seq就是首先使用酶切，将整个基因组酶切，然后两端加接头构建文库，上机进行单端或者双端测序。
二，Stacks分析原理
Stacks是分析RAD-Seq数据的常用流程，流程采用模块化设计，从RAD原始序列预处理到Call SNP以及最后的下游分析，涉及多个子模块，整体过程如下：
本文从源码层面梳理一下几个核心模块的内部实现：
2.1 数据准备
我们使用 ddrage软件来生产ddrad-seq的模拟数据，使用默认参数产生三个个体的数据
2.2 原理分析
1.process_radtags
对原始序列做质控.从Rad-seq图中可以看出，建库后的序列包括adaptor+barcode(index)+酶切点，adaptor用于引导测序，barcode用于标识单个样本，也叫index。通常情况下illumina会自动把adaptor和barcode切除，所以大部分情况radseq的序列头都是酶切后的序列。举一个实际例子，来看看RAD-Seq拿到的序列是什么样子：
假设某实验用的两个酶是SacI 和MseI ，
SacI 酶：
Fastq中的R1 reads头是：
MseI 酶：
Fastq中的R2 reads头是：
1.1 检查barcode、酶切点的完整性
1.2 根据barcode拆分样本
1.3 计算序列平均质量，丢弃低于90%平均质量的序列处理完原始数据后，根据是否有参考序列或者是否希望用参考序列，会走两个完全不同的流程。
2.无参考基因组的de novo过程：
Stacks 1.0 [3] 中对de novo 过程做了详细描述
2.1 ustacks
上图中蓝色框中的A~F步骤详细说明了ustacks的步骤，具体如下：
nohup ./stack2_ustacks.sh &
A 将单个样本中的序列分类为stack stack中的序列被称为primary序列
B 将stack拆分为kmer字典 并通过查询该字典找出相似的stack
C 相似的stack成为匹配 用图中的连线连接 匹配的节点合并为loci 不在stack中的序列被称为二级序列(secondary reads)
D 将二级序列也和loci进行比较 用以增加比较深度
E 构建一致性序列(consensus sequence)来记录分型和snp信息
2.2 cstacks
2.3 sstacks
maximum likelihood statistical model
3.有参考基因组的过程：讨论一：是否使用参考基因组？
参考:
1.N. Rochette, A. Rivera‐Colón, and J. Catchen. Stacks 2: Analytical methods for paired‐end sequencing improve RADseq‐based population genomics. Molecular Ecology, 28(21):4737-4754. 2019.
2.https://www.floragenex.com/sbg-ddrad-seq
3.https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3276136/
Stacks: Building and Genotyping Loci De Novo From Short-Read Sequences -
编译调试版本
./configure CFLAGS='-g -O0' CXXFLAGS='-g -O0'