海洋生物构造变异分析软件meerkat 0.189运用笔记(二)

一、 运行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再度好

微生物 1

 

 

 微生物 2

 微生物 3

 

 

         -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
:http://gensoft.pasteur.fr/docs/Meerkat/0.189/Manual\_0.189.pdf**
**

相关文章