zoukankan      html  css  js  c++  java
  • 生物结构变异分析软件meerkat 0.189使用笔记(一)

      一、准备工作

        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

          结果报错如下,

           

        作如下调试   

    make clean (清除刚才的安装)
    #修改makefile
    #from: ... -lbamtools -lbamtools-utils     
    #to: ... -lbamtools -lbamtools-utils -lz
    make

         编译成功,但是运行./bamreader 继续报错

       

       解决方法如下

    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

  • 相关阅读:
    AcWing 1027. 方格取数 dp
    AcWing 1014. 登山 dp
    acwing 482. 合唱队形 dp
    LeetCode 1463. 摘樱桃II dp
    LeetCode 100. 相同的树 树的遍历
    LeetCode 336. 回文对 哈希
    LeetCode 815. 公交路线 最短路 哈希
    算法问题实战策略 DARPA大挑战 二分
    算法问题实战策略 LUNCHBOX 贪心
    AcWing 1100. 抓住那头牛 BFS
  • 原文地址:https://www.cnblogs.com/ywliao/p/6438732.html
Copyright © 2011-2022 走看看