當前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

Meerkat软件

發布時間：2023/12/20 编程问答 63 豆豆

生活随笔收集整理的這篇文章主要介紹了 Meerkat软件小編覺得挺不錯的,現在分享給大家,幫大家做個參考.

一、準備工作

meerkat 0.189版本和以前的版本相比，支持bwa mem 輸出的bam文件，還支持全外顯子數據count SV。
meerkat原理

1.1 需要準備的軟件

unix/Linux系統（自帶）
CMake（自帶）
PERL 5.8.1及以上（自帶）
BioPERL 1.5.0及以上（自行安裝）
R 2.3.1及以上（自帶）
samtools 0.1.5到0.1.19(不支持新版本samtools)
BWA 0.6.2.（已經可以支持新版本的BWA，但是 split read alignment 的時候必須用0.6.2版本）
NCBI blast 2.2.24及以上（自行安裝）
Primer32.2.0及以上（自行安裝）

1.2 需要準備的文件

參考基因組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;

bwa index 對基因組文件建立的index
samtools faidx 對基因組文件建立的index

samtools faidx hg19ref_order.fa

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

UCSC下載Repeat annotation 基因注釋文件也可以在這里輸出

1.3 照著manual 安裝。

下載meerkat壓縮包，解壓。進入meerkat文件夾。

build mybamtools, 生成lib文件夾，文件夾包含著需要鏈接的動態庫

cd ./src/ tar xjvf mybamtools.tbz cd mybamtools mkdir build cd build cmake .. make

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

結果報錯,作如下調試

make clean (清除剛才的安裝） #修改makefile #from: ... -lbamtools -lbamtools-utils #to: ... -lbamtools -lbamtools-utils -lz<br>make

編譯成功，但是運行./bamreader 繼續報錯

解決方法如下

export LD_LIBRARY_PATH=/home/ywliao/bin/Meerkat/src/mybamtools/lib

將mybamtools/lib路徑加入LD_LIBRARY_PATH變量即可。

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/

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_dirperl /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_dirperl /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

前面已經依序安裝了meerkat 的環境和meerkat，運行了預處理一步，在相對應的bam文件目錄下生成了大批文件，因此，當要用meerkat處理某個bam文件時，應先將該bam文件移動到專有的一個文件夾，manual中也建議這樣用。
預處理生成的文件包括：
黑名單文件.gz
isinfo文件：包括插入大小信息
pdf文件：插入大小的分布圖，unmapped reads長度的分布圖，softclip reads長度分布圖
pre.log文件：日志文件，包括輸入的參數，輸入輸出信息，reads數，unallign reads數等
softclips.fq.gz： softlip reads文件，完整的read
softclips.rdist：softclip 的reads長度分布信息記錄
unmapped.fq.gz：unmapped的reads 文件，完整的read
unmapped.rdist： unmapped的reads長度的分布信息
sr1.fq.gz : 從softclip read 或者 unmaped read 切下來的指定bp 的reads
sr2.fq.gz : 與sr1.fq配對的reads
一個文件夾包括1_1fq.gz，1_2fq.gz ，里面有不同長度的reads。每個read group 人工的reads 對

1.1 meerkat.pl perl ./scripts/meerkat.pl [options] -b FIle 已經sort過的bam文件 -k int [0/1]是否使用預處理產生的黑名單文件，默認1 -d FLT call discordant mate pairs 的標準差閾值，默認3.這個選項控制怎樣尋找discordant read對，等同于定義最大的concordant fragment(配對的reads定位到的片段），計算方法是：中位插入大小+d x sd,如果插入大小分布圖狹窄并對稱，就使用默認參數，例如下面一二圖。如果分布圖偏寬，或者帶著長尾，三四圖，參數應設為5，對于高深度(>30x)，盡管分布圖狹窄，但是使用5會輕微好一點。如果峰窄，但是仍然帶著一個尾部(圖五)，好像大部分TCGA數據那樣，在meekat.pl這一步使用5比3更好 <img src="e31c18eb-9872-4c9a-ae8d-9a862aca9f97_files/4a514de1-be66-47e7-abd8-bd7a250cbd77.jpg" style="vertical-align: bottom; max-width: 100%;"><img src="e31c18eb-9872-4c9a-ae8d-9a862aca9f97_files/7811316c-a85a-4ced-9e23-8475c7156f7e.jpg" style="vertical-align: bottom; max-width: 100%;"><img src="e31c18eb-9872-4c9a-ae8d-9a862aca9f97_files/198fa504-d0b6-4df8-8f66-d79418680015.jpg" style="vertical-align: bottom; max-width: 100%;">-c FLT 集簇discordant mate pair標準差閾值，默認和-d相同。控制怎樣集簇，構建斷點的置信區間。等同于定義覆蓋斷點的最大片段（中位插入大小+c x sd)。如果-d 選項是5或者小于5，使用用-c相同的參數，如果-d 很大，比如10，那么使用更小的-c，比如5。如果數據有很高的測序深度，或者有廣泛的拷貝的變化，那么使用5而不是3來避免不能構建斷點的置信區間。 -p FLT call一個事件支持的配對數閾值，默認2 -o INT 需要支持全長reads對的數目，默認是0，設定這個選項會降低小復雜事件的敏感度。如果使用-a 1 -l 1推薦使用-o 1 來避免重復序列導致的過多人工產物 -q INT call一個事件支持的split reads數，默認是1 -z INT 事件數的最大值，默認1,000,000,000 -s INT 從unmapped reads 開始和末端切下來的bp數，必須和與處理的-s 參數設置一樣-m INT [0/1]，如果設定為1，使用meerkat 去去除duplicates，如果設為0，使用Picard 的 flag 'd‘ marked ，或者其他工具去除duplicates.bwa mem 的數據必須用picard mark duplicate ，所以要設為0. -a INT [0/1],處理非單一比對，默認1。如果測序質量不好，或者測序深度不夠，將參數設為0，這樣全部成對的reads都被作為單一比對使用。這依賴bwa mapping 時生成的XT標簽。如果沒有XT標簽，使用選項Q. -u INT [0/1],使用bam文件的全部比對，如果BAM文件不是由BWA產生的，或者由bwa產生但是沒有XT標簽，那么開這個選項，開了這個選項會強制關閉-a選項。默認是0。開了這個選項之后應使用-Q 參數過濾低mapping reads,推薦-Q設置10，對于bwa比對過的bam，也可以繼續過濾。 -Q INT 被使用的reads的最小mapping quality，默認是0 -g INT 在主要的bam文件使用的備擇mapping數，默認使用全部。備擇mapping 數之前通過bwa -N 參數輸出到XA標簽中。bwa mem 默認輸出備擇mapping。 -f INT 在clipped 比對中考慮的備擇mapping數，默認使用全部。 -l INT [0/1],是否考慮clipped 比對，默認1.和預處理一樣，對于bwa mem 出來的基因組，不需要重新mapping. -t INT 在bwa比對中用到的核數，默認1 -R STR 包括黑名單read group的文件，每一行一個read group ID -F STR fasta文件夾路徑 -S STR samtools路徑，如果samtools不在環境變量中的話 -W STR bwa路徑，如果bwa不在環境變量中的話 -B STR blastall 和 formatdb 的路徑，如果不在環境變量的話 -P STR 指定運行順序，dc|cl|mpd|alg|srd|rf|all，默認all，每一步運行都需要上一步的結果dc: extract discordant read pairs，提取discordant read對cl: construct clusters of discordant read pairs，將discordant read對建簇mpd: call events based on read pairs，基于read 對call 事件alg: align split reads to candidate break point regions，比對split read 到候選的斷點區域。如果你有多核心的話，可以將alg步驟切為兩步,alg1和alg2，alg1運行多線程，alg2運行單線程。srd: confirm events based on split reads and filter results，基于split reads和過濾的結果對事件進行證明rf: refine break points by local alignments，通過區域比對定義斷點all: run all above steps ,運行全部步驟 -h 幫助

manual 中的舉例：
50bp reads，<10x 的TCGA 基因組使用-s 18 -d 5 -a 0 -l 0 -q 1，猜測：reads 長度較小，所以取1/3 read 長度，-s 取18， TCGA 基因組,插入分布狹窄帶尾，所以-d 設為5，測序深度較低，reads 長度較短，所以盡量保留數據，-a 設為0, -l 設為0，-q 設的較低，1。
75-101bp reads, 30-40x TCGA 基因組使用-s 20 -d 5 -p 3 -o 1 -a 0 -u 1 -Q 10，猜測：reads長度較長，取1/5read長度，-s 取20，TCGA 基因組,插入分布狹窄帶尾，所以-d 設為5，長度較長，深度較深，因此降低敏感度，增加特異度，所以-p 設為3，-o 設為1，由于沒有XT tag和XA tag ，因此-a 1 選項無法運行，因此設為-u 1 -a 0 -Q 10 。如果這是bwa mem的數據的話，參數應設為-s 20 -d 5 -p 3 -o 1 -m 0 -l 0，bwa mem 數據不需要重新比對softclips -l 0，必須用picard去除duplicate，-m設為0,估計這個是早版本的bwa，因此比對時可以生成XT標簽，-a 默認為1。
101bp reads, 60-80x TCGA 基因組使用-s 20 -d 5 -p 5 -o 3，75-101bp 使用-s 20，TCGA基因組使用-d 5，測序深度深，-o 設置更高3。
如果腫瘤基因組60x，相對應的正常基因組30x，那么使用60x的參數，常常用配對的正常組織中的black.list.gz 重命名并放到腫瘤bam文件處理的文件夾，替換cancer的blacklist.gz文件。

1.2 mechanism.pl perl ./scripts/mechanism.pl [options] -b FILE sort過的bam文件 -o INT [0/1]在TE包括rmsk類型 \"Other\"，默認1 -t INT TE的最大值，默認100000 -z INT 被處理的SV的數量限制，默認1000000000 -R STR repeat注釋路徑，能夠從UCSC下載 -h help

二、參照

manual中給出的數據，HapMap個體NA18507（42x深度，100bp讀長，500bp插入大小），用10核bwa比對花費1.5天和10GB的內存。30x深度的腫瘤數據，要多于兩天并且超過30G的內存，如果測序質量不太好，比如很多的嵌合比對和許多非單一mapped的reads，就會花費更多的時間和內存。、

三、結果

運行完預處理和meerkat.pl后，能夠得到兩個文件.intra.refined.typ.sorted和inter.refined.typ.sorted，運行完mechanism.pl后，會得到.variants文件，否則，就該回去看一下設置是否出現問題。運行的腫瘤基因組文件的時候，一定要把正常組織的blacklist文件替換腫瘤基因組的blacklist文件，blacklist文件可以在預處理中生成。如果在預處理中出現錯誤信息“differing read lengths“，沒有關系，僅僅是告訴你在一些read group中長度不一致。

四、包含的其他程序

4.1 轉換variant 文件到vcf格式 perl ./scripts/meerkat2vcf.pl -i variantfile -o vcffile -H headerfile -F /db/hg18/hg18_fasta/ -F選項可以產生一個頭文件 4.2 融合基因注釋 perl ./scripts/fusions.pl -i variantfile -G /db/hg18/refGene_hg18_sorted.txt 4.3 為變異位點設計引物使用primer.pl -i FILE 輸入文件，fusion.pl產生的融合文件-o FILE 輸出文件 -p STR 引物固定序列-c INT 列數補償，默認0，如果第一列存在如何事件的樣品名稱，那么設為1 -f INT 側翼區域，默認500，設計引物所在的區域-s STR 被“，”分隔的引物大小，默認20,23,25,27 -m STR 引物最小，最優，最大Tm值，","分隔，默認50，60，65 -n INT 對每一個引物片段，設計的音物個數，默認5 -r INT 掩蓋repeat，默認0 -q INT 輸出設計引物的側翼序列 -F STR fasta文件文件夾路徑 -P STR primer_core程序文件夾路徑 -L STR bla文件夾t路徑 -V STR blat服務器，例如fServer start 10.11.240.76 17777/reference/hg18/hg18.2bit -stepSize=5，服務器路徑應該為10.11.240.76 -T STR blat端口，上面例子中的17777 -h help

全部音物都是由Primer3生成，對于每一個事件，挑出.1和.2，不同的取向認為是不同的事件，所以取引物的時候直接拷貝出來，不需要額外的反向互補，如果序列是小寫字母，那么說明引物在重復序列。理論上，你應該用 1 blat hit 挑出小寫的引物。如果blat hit 是0，意味著由太多hit了，所以不要挑選這種引物。有時候，即使引物是在重復序列上（小寫字母），但是在基因組上仍然是單一比對的，（1 blat hit)，因為是重復元件的變異，挑選這種引物是可以的。如果由一些引物PCR沒有結果，你可以挑選2個正向引物，兩個反向引物來同時進行4 個PCR 反應。引物設計的普遍規則仍然要使用，比如，你應該挑選TM 值相差不大并且GC含量不太極端的。 4.4 計算等位基因頻率使用discon.pl,這個腳本會告訴你全部斷點的discordant 和 concordant的read對. -i FILE 輸入variants 文件，必須 -o FILE 輸入文件，必須 -D INT 從bam文件中計算discordant pair的數目，默認0，基因分型的時候打開這個選項 -B FILE bam文件，必須 -C FILE 由Meerkat 產生的cluster文件，必須-I FILE Meerkat 運行時產生的isinfo文件 -K FILE 包含要忽略的read group的文件名，一個read ID 一行，和 meerkat.pl的R參數一樣 -S FILE samtools文件夾路徑，如果samtools不再環境變量中 -d FLT call discordant read 對的標準差閾值，默認3 -Q INT 使用的reads 最小的mapping quality，默認0 -h help

perl ./scripts/discon.pl -d 5 -Q 10 -i $somaticg -o $somaticg_rp -B $cancer_bam -C $cancer_cluster -I $cancer_isinfo -K $cancer_blacklistrg -S /home/ly55/opt/samtools/ 結果文件里面每一個事件有一個 RP tag ，A_B_C_D格式A：全長的discordant read 對數B：從部分比對的reads中discordant的reads數C：第一斷點的concordant reads 數D：第二斷點的concordant reads數A+B應該等與Meerkat得到的總的reads數

參考資料

meerkat manual

轉載于:https://www.cnblogs.com/raisok/p/10918187.html

總結

以上是生活随笔為你收集整理的Meerkat软件的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。

上一篇：网校系统是怎样搭建的？
下一篇：牛客 NC24858 [USACO 20