酶切实验一共进行了7组实验，并通过凝胶电泳检测结果，同时，最终得出最佳酶切组合为：MluCI + AvaII，用stacks软件提供的process_radtags程序对测序公司返回的数据进行拆包，发现process_radtags并不支持avaII限制性内切酶，查看process_radtags的使用手册，该软件支持的内切酶列表如下：
'aciI', 'ageI', 'aluI', 'apaLI', 'apeKI', 'apoI', 'aseI', 'bamHI', 'bbvCI', 'bfaI', 'bfuCI', 'bgIII', 'bsaHI', 'bspDI', 'bstYI', 'btgI', 'cac8I', 'claI', 'csp6I', 'ddeI', 'dpnII', 'eaeI', 'ecoRI', 'ecoRV', 'ecoT22I', 'haeIII', 'hinP1I', 'hindIII', 'hpaII', 'hpyCH4IV', 'kpnI', 'mluCI', 'mseI', 'mslI', 'mspI', 'ncoI', 'ndeI', 'nheI', 'nlaIII', 'notI', 'nsiI', 'nspI', 'pacI', 'pspXI', 'pstI', 'rsaI', 'sacI', 'sau3AI', 'sbfI', 'sexAI', 'sgrAI', 'speI', 'sphI', 'taqI', 'xbaI', or 'xhoI'.
由于stacks是开源的软件，我们可以对其程序进行修改，将avaII的酶切位点信息加入到程序中。
第一步下载stacks
stacks2.54下载地址：
https://catchenlab.life.illinois.edu/stacks/
# tar -xzf stacks-2.54.tar.gz
# cd stacks-2.54/
第二步，查找avaII限制性内切酶的酶切位点
结果为：
5'···G^GWCC···3'
3'···CCWG^G···5'
其中W代表A或者T
第三步，编辑renc.cc文件将酶切位点的信息加入其中
vi ./stacks-2.54/src/renz.cc
一共要修改四处地方：
第一处添加：const char *avaII[] = {"GACC","GTCC", "GGTC","GGAC"}; // 第1/4处添加，添加酶切位点信息，输入酶切序列和反向互补序列
第二处添加：renz["avaII"] = avaII; // 第2/4处添加，G^GWCC, avaII
第三处添加：renz_cnt["avaII"] = 2; // 第3/4处添加,输入不同的酶切序列数量，注意：反向互补序列不算
第四处添加：renz_len["avaII"] = 4; // 第4/4处添加,输入酶切序列长度
具体添加位置见如下代码：
本人修改过的renz.h文件已经存入网盘
#include "renz.h"
using namespace std;
//
// First line of static array contains each enzyme's cut sites. Second
// line is the reverse complement of each cut site.
//
const char *avaII[] = {"GACC","GTCC",
"GGTC","GGAC"}; // 第1/4处添加，添加酶切位点信息，输入酶切序列和反向互补序列
const char *aciI[] = {"CGC", "CGG", // C/CGC, AciI
"GCG", "CCG"};
const char *ageI[] = {"CCGGT", // A/CCGGT, AgeI
"ACCGG"};
······此处省略代码
void
initialize_renz(map<string, const char **> &renz, map<string, int> &renz_cnt, map<string, int> &renz_len) {
renz["avaII"] = avaII; // 第2/4处添加，G^GWCC, avaII
renz["sbfI"] = sbfI; // CCTGCA/GG, SbfI
renz["pstI"] = pstI; // CTGCA/G, PstI
renz["notI"] = notI; // GC/GGCCGC, NotI
······此处省略代码
renz["pspXI"] = pspXI; // VC/TCGAGB, PspXI
renz["hpyCH4IV"] = hpyCH4IV; // A/CGT, HpyCH4IV
renz_cnt["avaII"] = 2; // 第3/4处添加,输入不同的酶切序列数量，注意：反向互补序列不算
renz_cnt["sbfI"] = 1;
renz_cnt["pstI"] = 1;
renz_cnt["notI"] = 1;
······此处省略代码
renz_cnt["pspXI"] = 3;
renz_cnt["hpyCH4IV"] = 1;
renz_len["avaII"] = 4; // 第4/4处添加,输入酶切序列长度
renz_len["sbfI"] = 6;
renz_len["pstI"] = 5;
renz_len["notI"] = 6;
wq保存并退出
# cd ..
# ./configure
出现如下错误：
checking for functional regular expressions... configure: error: Regular expressions are not functional, you need g++ 4.9.0 or greater.
原来是gcc版本低了，输入 gcc -v 查看当前gcc版本
Thread model: posix
gcc version 4.8.5 20150623 (Red Hat 4.8.5-44) (GCC)
当前我的gcc版本是4.8.5，要求的gcc版本大于4.9.0
接下来是升级gcc版本
# yum install centos-release-scl
# yum install devtoolset-7 // 升级到gcc7
# scl enable devtoolset-7 bash // 激活gcc版本
# gcc -v
Using built-in specs.
COLLECT_GCC=gcc
COLLECT_LTO_WRAPPER=/opt/rh/devtoolset-7/root/usr/libexec/gcc/x86_64-redhat-linux/7/lto-wrapper
Target: x86_64-redhat-linux
Configured with: ../configure --enable-bootstrap --enable-languages=c,c++,fortran,lto --prefix=/opt/rh/devtoolset-7/root/usr --mandir=/opt/rh/devtoolset-7/root/usr/share/man --infodir=/opt/rh/devtoolset-7/root/usr/share/info --with-bugurl=http://bugzilla.redhat.com/bugzilla --enable-shared --enable-threads=posix --enable-checking=release --enable-multilib --with-system-zlib --enable-__cxa_atexit --disable-libunwind-exceptions --enable-gnu-unique-object --enable-linker-build-id --with-gcc-major-version-only --enable-plugin --with-linker-hash-style=gnu --enable-initfini-array --with-default-libstdcxx-abi=gcc4-compatible --with-isl=/builddir/build/BUILD/gcc-7.3.1-20180303/obj-x86_64-redhat-linux/isl-install --enable-libmpx --enable-gnu-indirect-function --with-tune=generic --with-arch_32=i686 --build=x86_64-redhat-linux
Thread model: posix
gcc version 7.3.1 20180303 (Red Hat 7.3.1-5) (GCC)
成功的将gcc升级到7.3.1
继续之前出错的步骤
# ./configure
# make
出现如下结果，表示编译成功
g++ -g -O2 -fopenmp -o ustacks src/ustacks.o htslib/libhts.a libcore.a -lz -lgomp
g++ -g -O2 -fopenmp -o cstacks src/cstacks.o htslib/libhts.a libcore.a -lz -lgomp
g++ -g -O2 -fopenmp -o process_radtags src/process_radtags.o htslib/libhts.a libcore.a libclean.a -lz -lgomp
g++ -g -O2 -fopenmp -o sstacks src/sstacks.o htslib/libhts.a libcore.a -lz -lgomp
g++ -g -O2 -fopenmp -o process_shortreads src/process_shortreads.o htslib/libhts.a libcore.a libclean.a -lz -lgomp
g++ -g -O2 -fopenmp -o kmer_filter src/kmer_filter.o htslib/libhts.a libcore.a libclean.a -lz -lgomp
g++ -g -O2 -fopenmp -o clone_filter src/clone_filter.o htslib/libhts.a libcore.a libclean.a -lz -lgomp
g++ -g -O2 -fopenmp -o populations src/populations.o htslib/libhts.a libcore.a libpop.a -lz -lgomp
g++ -g -O2 -fopenmp -o phasedstacks src/phasedstacks.o htslib/libhts.a libcore.a -lz -lgomp
g++ -g -O2 -fopenmp -o tsv2bam src/tsv2bam.o htslib/libhts.a libcore.a -lz -lgomp
g++ -g -O2 -fopenmp -o gstacks src/gstacks.o src/debruijn.o src/SuffixTree.o htslib/libhts.a libcore.a -lz -lgomp
make[2]: Leaving directory '/home/packages/stacks-2.54'
make[1]: Leaving directory '/home/packages/stacks-2.54'
接下来我们使用修改过的process_radtags拆包RAD-seq原始文件
# /home/packages/stacks-2.54/process_radtags -P -p ./raw/ -b ./pool-1.1.txt -o ./samples/ -c -q -r --inline_index --renz_1 mluCI --renz_2 avaII
Processing paired-end data.
Using Phred+33 encoding for quality scores.
Found 13 paired input file(s).
Searching for single-end, inlined and paired-end, indexed barcodes.
Loaded 48 barcodes (5bp / 6bp).
Will attempt to recover barcodes with at most 1 / 1 mismatches.
Processing file 1 of 13 [51_R1_001.fastq.gz]
Reading data from:
./raw/51_R1_001.fastq.gz and
./raw/51_R2_001.fastq.gz
Processing RAD-Tags...1M...2M...3M...
可以看到添加的酶切位点已经能够成功被识别
总结：
经历上述添加酶切位点信息的过程可以发现，将所有的限制性内切酶都添加进去，对stacks的作者来说，是一件非常简单的事情。但stacks的作者并没有选择这么做的原因是：他不能保证加入其他酶切信息能够得到准确的结果，毕竟他能够测试的数据集也是有限的。如果添加的酶切位点造成意想不到的结果，建议向作者反馈异常信息，提升软件的适用范围，帮助更多的人。
创作不易，转载请注明出处！
