生物结构变异分析软件meerkat 0.189使用笔记（一）
一、準備工作
meerkat 0.189版本和以前的版本相比，支持bwa mem 輸出的bam文件，還支持全外顯子數(shù)據(jù)count SV。
meerkat原理：參見http://compbio.med.harvard.edu/Meerkat/
1.1 需要準備的軟件
1. unix/Linux系統(tǒng)（自帶）
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.（已經(jīng)可以支持新版本的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文件夾，文件夾包含著需要鏈接的動態(tài)庫
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
結(jié)果報錯如下，
作如下調(diào)試
|
1 2 3 4 |
makeclean (清除剛才的安裝）#修改makefile#from: ... -lbamtools -lbamtools-utils#to: ... -lbamtools -lbamtools-utils -lz<br>make
|
編譯成功，但是運行./bamreader 繼續(xù)報錯
解決方法如下
| 1 | exportLD_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/
二、預(yù)處理
manual明確指出不建議用默認參數(shù)
perl ./scripts/pre_process.pl [options]
-b FILE 已經(jīng)sort和index的bam文件
-k INT 過濾掉的最小的覆蓋度（過濾覆蓋過多reads的位置，默認500;過濾mapped到著絲粒的reads，通過它顯示出的覆蓋次數(shù)，在腫瘤樣品中應(yīng)該觀察拷貝數(shù)，應(yīng)設(shè)置一個更高的數(shù)值，比如1500，以至于不忽略這些事件
-r INT 被用于計算分布的插入長度的幅度（默認1000),會生成一個pdf的分布圖，顯示插入片段長度的分布，0關(guān)掉這個函數(shù)
-n INT 每個read group被用于計算插入片段大小分布的reads數(shù)，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 參數(shù)就是控制提取的部分的大小，這樣人為的reads對重新mapping 到參考基因組。如果插入片段小但是read的長度長，那么就會很大的增加敏感性。對于短長度的read，應(yīng)設(shè)置為0。對于bwa mem 出來的基因組，不需要重新mapping，所以可以關(guān)掉這一參數(shù)，在meerkat.pl中也一樣。
-q INT 削減reads數(shù)，等同于bwa 的-q參數(shù)，默認15
-c INT 設(shè)置開始和末端剪下softclip 或者unmapped 的read的bp數(shù)，這些剪下的reads 用來比對參考基因組，尋找更小事件。應(yīng)輕微小于1/2的read長度，默認參數(shù)適合長讀長的人類基因組。
-s INT 設(shè)置開始和末端剪下softclip 或者unmapped 的read的bp數(shù)，這些剪下的reads 用來split reads mapping，必須和下一步meerkat的-s參數(shù)設(shè)置一樣。在不犧牲mapping能力的情況下，值可以設(shè)的小一點。應(yīng)該設(shè)置1/5到1/3的read長度。
-u INT 處理uu reads 對(雙unmapped reads對，分成4對。默認0。如果測序質(zhì)量夠好，并且基因組沒有什么重復(fù)的話，對于識別小事件非常有用，人類基因組建議關(guān)閉函數(shù)。
-f INT clip 比對時允許輸出到XA標簽的備擇比對數(shù)量，默認100
-N INT clip和split reads必須Ns閾值，默認是5
-t INT bwa align用到的線程數(shù)
-R STR 包含黑名單reads的文件，一個group id 一行，如果對于一個group的單一比對reads少于30%，推薦不出這個group，如果group的
-I STR bwa_index路徑,bwa index 生成的參考基因index路徑，不是文件，用于bwa align，如果l（L發(fā)音）參數(shù)設(shè)為1的話應(yīng)設(shè)置
-A STR參考基因的fasta.fai文件，用于bwa align(查看代碼發(fā)現(xiàn)就是上文提到的samtools建立的參考基因的fai文件）
-S STR samtools路徑，如果不存在于環(huán)境變量的話
-W STR bwa途徑，如果不存在于環(huán)境變量的話
-P STR 指定運行步驟，[ all | is | cl ]，默認全部
is:提取unmapped,softclip reads，計算插入片段分布
cl: map soft clip 配對reads 到參考基因組，如果使用多線程的話，應(yīng)分步，cl1運行多線程，cl2運行單線程
-h help
manual中給出的例子。
1. 50bp的reads,<10x TCGA 基因組
建議使用-s 18 -l 0 -q 0
估計50bp片段過小，所以-s 選項 以1/3 read 長度，短長度reads，-l設(shè)置為0，估計測序深度不深，所以 并不trimming reads,所以-q 設(shè)置為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 選項應(yīng)該設(shè)置為1/5 read 長度，20,因為人類癌癥基因，所以-k選項設(shè)為1500，測序深度夠深，所以可以設(shè)置過濾的basequality為15。bwa mem mapping的數(shù)數(shù)據(jù)-l設(shè)置為0
3. TCGA WES 數(shù)據(jù)
建議使用-s 20 -k 10000 -q 5
-k 10000表示10000的copy number的reads也會留下，-q 5,就是過濾的basequality為5
這次我們實驗室分析的數(shù)據(jù)，150bp 讀長，測序深度8x,bwa mem 腫瘤數(shù)據(jù)，我選擇參數(shù)為-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
總結(jié)
以上是生活随笔為你收集整理的生物结构变异分析软件meerkat 0.189使用笔记（一）的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 更新npm至最新版本
- 下一篇： 怎么创建具有真实纹理的CG场景岩石？