一、准备工作
meerkat 0.189版本和以前的版本相比，支持bwa mem 输出的bam文件，还支持全外显子数据count SV。
meerkat原理：参见http://compbio.med.harvard.edu/Meerkat/
1.1 需要准备的软件
1. unix/Linux系统（自带）
2. CMake（自带）
3. PERL 5.8.1及以上（自带）
4. BioPERL 1.5.0及以上（自行安装）
5. R 2.3.1及以上（自带）
6. samtools 0.1.5到0.1.19(不支持新版本samtools)
7. BWA 0.6.2.（已经可以支持新版本的BWA，但是 split read alignment 的时候必须用0.6.2版本）
8. NCBI blast 2.2.24及以上（自行安装）
9. Primer32.2.0及以上（自行安装）
1.2 需要准备的文件
1.参考基因组fasta文件（单独放在文件夹），运行perl脚本，用BioPerl的Bio:DB::Fasta进行处理
#!/bin/perl
use Bio::DB::Fasta;
# Create database from a directory of Fasta files
my $db = Bio::DB::Fasta->new('/home/ywliao/utilities/UCSC/hg19_FA');
my @ids = $db->get_all_primary_ids;
2.bwa index 对基因组文件建立的index（实验室已有）
3. samtools faidx 对基因组文件建立的index
samtools faidx hg19ref_order.fa
4.UCSC下载的参考基因注释文件，knowGene.txt 用sort refGene.txt -k 3,3 -k 5,5n > refGene_sorted.txt命令进行sort
sort knownGene.txt -k 3,3 -k 5,5n > hg19_knowGene_sorted.txt
5.UCSC下载Repeat annotation。（基因注释文件也可以在这里输出)
1.3 照着manual 安装。
下载meerkat压缩包，解压。进入meerkat文件夹。
1.build mybamtools, 生成lib文件夹，文件夹包含着需要链接的动态库
cd ./src/ tar xjvf mybamtools.tbz cd mybamtools mkdir build cd build cmake .. make
2.build bamreader
tar xjvf bamreader.tbz cd bamreader # Edit Makefile and set BTROOT to the path to which mybamtools was extracted. #vi Makefile #BTROOT = /home/ywliao/bin/Meerkat/src/mybamtoolsmake mv ./bamreader ../../bin
结果报错如下，
作如下调试
|
1
2
3
4
|
make clean (清除刚才的安装）#修改makefile#from: ... -lbamtools -lbamtools-utils #to: ... -lbamtools -lbamtools-utils -lz<br>make |
编译成功，但是运行./bamreader 继续报错
解决方法如下
| 1 | export LD_LIBRARY_PATH=/home/ywliao/bin/Meerkat/src/mybamtools/lib |
将mybamtools/lib路径加入LD_LIBRARY_PATH变量即可。
3.build dre
tar xjvf dre.tbz cd dre #Edit Makefile and set BTROOT to the path to which mybamtools was extracted. #vi Makefile #BTROOT = /home/ywliao/bin/Meerkat/src/mybamtools/
make mv ./dre ../../bin/
4.build sclus
tar xjvf sclus.tbz cd sclus make mv ./sclus ../../bin/
二、预处理
manual明确指出不建议用默认参数
perl ./scripts/pre_process.pl [options]
-b FILE 已经sort和index的bam文件
-k INT 过滤掉的最小的覆盖度（过滤覆盖过多reads的位置，默认500;过滤mapped到着丝粒的reads，通过它显示出的覆盖次数，在肿瘤样品中应该观察拷贝数，应设置一个更高的数值，比如1500，以至于不忽略这些事件
-r INT 被用于计算分布的插入长度的幅度（默认1000),会生成一个pdf的分布图，显示插入片段长度的分布，0关掉这个函数
-n INT 每个read group被用于计算插入片段大小分布的reads数，0 使用全部reads，默认1000
-l INT 提取配对的softclip reads，或者其他配对的，但是有某一个mapped不上或者都mapped不上的reads，默认1。对于插入片段很小的，在sv断点上就会有reads覆盖，这样得到的reads就会部分比对或者比对不上。运行的时候，对于一个末端mapped上，另一个read末端部分比对上的reads对，会把部分比对read的unmapped部分提取出来和mapped的read组成人为的read对；对于一个末端比对上，一个末端unmapped的reads对，那么unmapped read 的起始和末端的序列分别提取和mapped的read组成两对人为的read对；-c 参数就是控制提取的部分的大小，这样人为的reads对重新mapping 到参考基因组。如果插入片段小但是read的长度长，那么就会很大的增加敏感性。对于短长度的read，应设置为0。对于bwa mem 出来的基因组，不需要重新mapping，所以可以关掉这一参数，在meerkat.pl中也一样。
-q INT 削减reads数，等同于bwa 的-q参数，默认15
-c INT 设置开始和末端剪下softclip 或者unmapped 的read的bp数，这些剪下的reads 用来比对参考基因组，寻找更小事件。应轻微小于1/2的read长度，默认参数适合长读长的人类基因组。
-s INT 设置开始和末端剪下softclip 或者unmapped 的read的bp数，这些剪下的reads 用来split reads mapping，必须和下一步meerkat的-s参数设置一样。在不牺牲mapping能力的情况下，值可以设的小一点。应该设置1/5到1/3的read长度。
-u INT 处理uu reads 对(双unmapped reads对，分成4对。默认0。如果测序质量够好，并且基因组没有什么重复的话，对于识别小事件非常有用，人类基因组建议关闭函数。
-f INT clip 比对时允许输出到XA标签的备择比对数量，默认100
-N INT clip和split reads必须Ns阈值，默认是5
-t INT bwa align用到的线程数
-R STR 包含黑名单reads的文件，一个group id 一行，如果对于一个group的单一比对reads少于30%，推荐不出这个group，如果group的
-I STR bwa_index路径,bwa index 生成的参考基因index路径，不是文件，用于bwa align，如果l（L发音）参数设为1的话应设置
-A STR 参考基因的fasta.fai文件，用于bwa align(查看代码发现就是上文提到的samtools建立的参考基因的fai文件）
-S STR samtools路径，如果不存在于环境变量的话
-W STR bwa途径，如果不存在于环境变量的话
-P STR 指定运行步骤，[ all | is | cl ]，默认全部
is:提取unmapped,softclip reads，计算插入片段分布
cl: map soft clip 配对reads 到参考基因组，如果使用多线程的话，应分步，cl1运行多线程，cl2运行单线程
-h help
manual中给出的例子。
1. 50bp的reads,<10x TCGA 基因组
建议使用-s 18 -l 0 -q 0
估计50bp片段过小，所以-s 选项 以1/3 read 长度，短长度reads，-l设置为0，估计测序深度不深，所以 并不trimming reads,所以-q 设置为0
2. 101bp reads, 20-30x and 60-80x TCGA 基因组
建议使用-s 20 -k 1500 -q 15
如果是bwa mem出来的文件，建议使用-s 20 -k 1500 -q 15 -l 0
75-101bp的reads，-s 选项应该设置为1/5 read 长度，20,因为人类癌症基因，所以-k选项设为1500，测序深度够深，所以可以设置过滤的basequality为15。bwa mem mapping的数数据-l设置为0
3. TCGA WES 数据
建议使用-s 20 -k 10000 -q 5
-k 10000表示10000的copy number的reads也会留下，-q 5,就是过滤的basequality为5
这次我们实验室分析的数据，150bp 读长，测序深度8x,bwa mem 肿瘤数据，我选择参数为-s 30 -k 1500 -q 0 -l 0
perl /home/ywliao/bin/Meerkat/scripts/pre_process.pl -b $inputfile -k 1500 -t 20 -l 0 -q 0 -P is -A $hg19_fai -W $bwa_dir -s 30 -S $samtools_dir perl /home/ywliao/bin/Meerkat/scripts/pre_process.pl -b $inputfile -k 1500 -t 20 -l 0 -q 0 -P cl1 -A $hg19_fai -W $bwa_dir -s 30 -S $samtools_dir perl /home/ywliao/bin/Meerkat/scripts/pre_process.pl -b $inputfile -k 1500 -t 20 -l 0 -q 0 -P cl2 -A $hg19_fai -W $bwa_dir -s 30 -S $samtools_dir
参考资料
meerkat manual :http://gensoft.pasteur.fr/docs/Meerkat/0.189/Manual_0.189.pdf