nucmer全基因组序列比对

作者: DumplingLucky | 来源:发表于2021-05-03 07:23 被阅读0次

    1. MUMmer简介

    MUMmer(Maximal Unique Match mer)软件中使用率最高的命令nucmer能用于对两个基因组assemblies进行比较。现在最新版本的MUMmer4相比于MUMmer3,能用于更大的基因组序列比较,运行速度更快。

    2. MUMmer4软件安装

    先安装高版本autoconf、automake和yaggo

    wget http://ftp.gnu.org/gnu/autoconf/autoconf-2.69.tar.gz -P ~/software_packages/
    tar zxf ~/software_packages/autoconf-2.69.tar.gz
    cd autoconf-2.69
    ./configure --prefix=/opt/sysoft/autoconf-2.69
    make -j 24 && make install
    echo 'PATH=/opt/sysoft/autoconf-2.69/bin/:$PATH' >> ~/.bashrc
    source ~/.bashrc
    cd .. && rm -rf autoconf-2.69
    
    wget http://ftp.gnu.org/gnu/automake/automake-1.16.tar.gz -P ~/software_packages/
    tar zxf ~/software_packages/automake-1.16.tar.gz
    cd automake-1.16
    ./configure --prefix=/opt/sysoft/automake-1.16
    make -j 24 && make install
    echo 'PATH=/opt/sysoft/automake-1.16/bin/:$PATH' >> ~/.bashrc
    source ~/.bashrc
    cd .. && rm -rf automake-1.16
    
    wget https://github.com/gmarcais/yaggo/releases/download/v1.5.10/yaggo-1.5.10.tgz -P ~/software_packages/
    tar zxf ~/software_packages/yaggo-1.5.10.tgz -C /opt/sysoft/
    make
    echo 'PATH=/opt/sysoft/yaggo-1.5.10:$PATH' >> ~/.bashrc
    source ~/.bashrc
    

    再安装MUMmer4

    $ wget https://github.com/mummer4/mummer/archive/v4.0.0beta2.tar.gz -O -P ~/software_packages/mummer-4.0.0beta2.tar.gz
    tar zxf ~/software_packages/mummer-4.0.0beta2.tar.gz
    cd mummer-4.0.0beta2
    autoreconf -fi
    ./configure --prefix=/opt/biosoft/mummer-4.0.0beta2
    source ~/.bashrc.gcc    要使用高版本GCC进行编译
    make -j 24 && make install
    cp -r MANUAL.md docs /opt/biosoft/mummer-4.0.0beta2/
    cd .. && rm -rf mummer-4.0.0beta2
    cd /opt/biosoft/mummer-4.0.0beta2/docs/
    latex maxmat3man.tex
    dvipdf maxmat3man.dvi
    echo 'PATH=$PATH:/opt/biosoft/mummer-4.0.0beta2/bin/' >> ~/.bashrc
    source ~/.bashrc
    

    若需要使用MUMmer调用gunplot绘图,则需要安装高版本的gunplot

    wget https://sourceforge.net/projects/gnuplot/files/gnuplot/5.2.2/gnuplot-5.2.2.tar.gz -P ~/software_packages/
    tar zxf ~/software_packages/gnuplot-5.2.2.tar.gz 
    cd gnuplot-5.2.2
    ./configure --prefix=/opt/sysoft/gnuplot-5.2.2 && make -j 24 && make install
    cd .. && rm -rf gnuplot-5.2.2
    echo 'PATH=/opt/sysoft/gnuplot-5.2.2/bin/:$PATH' >> ~/.bashrc
    source ~/.bashrc
    

    3. mumer命令

    mummer命令是MUMmer软件的核心程序。其它程序,如nucmer就是perl编写的调用此核心程序进行数据分析的脚本。由于该程序是核心程序,MUMmer软件给出了该命令的详细说明文档。

    mummer命令用于计算query assembly和reference assembly中>=20bp的匹配序列(maximal matches),该段maximal matches序列在query和reference基因组序列中完全一致,且延长到最长。

    mumer命令的使用:

    mumer用法
    $ mummer [options] <ref.fasta> <query1.fasta> [query2.fasta] ...
        mummer只能输入一个reference assembly文件,最多支持输入32个query 
    assembly文件。程序将结果输出到标准输出。
    常用示例:
    $ mummer -mumreference -l 20 -b -n -c -F ref.fa query.fa > out.txt
    
    常用参数:
    -mum
    -mumreference
    -maxmatch
        mummer程序用于计算获得maximal matches。程序通过suffix array算法将
    reference assembly构建索引数据库,然后再将query assembly和数据库比对,
    搜索并得到maximal matches。
        有很多maximal matches可能属于重复序列,而重复序列的比对可能没太大意义,
    且重复序列的比对极大消耗计算量。因此,需要考虑是否需要得到重复序列的比对结果。
        以上3个参数则对应重复序列的处理方法:默认是-mumreference参数,表示构建
    数据库时,去掉了reference assemly中的重复序列,得到的结果中,maximal 
    matches在reference assembly中是unique的,这种方法计算速度最快;-mum参数
    则是在前者基础上,对query assembly先去除重复序列,再进行比对,得到的结果中,
    maximal matches在reference assembly和query assembly中都是unique的,
    这种方法计算速度较慢;-maxmatch参数则完全不考虑重复序列,该方法计算速度最慢,
    同时生成的结果文件非常大。
    
    -l <int>    default: 20
        设置maximal match的最小长度。
    
    -b
    -r
        默认设置下,mummer进行maximal matches搜索时,要求query和reference
    上的maximal matches序列完全一致;若加上-r参数,则要求query和reference
    上的maximal matches序列反向互补;若加上-b参数,则以上两种情况都考虑。
    
    -n
        mummer对query和reference序列进行比较时,能识别任意字符,且对字符大小写
    不敏感,因此,也适合于蛋白序列的比较。若加入该参数,则仅识别ATCG四种碱基字符,
    所有其他字符都认为是不能匹配的。适合于assembly中含有W/Y/N等简并碱基的情况。
    
    -c
        若加入参数-b或-r,则会考虑maximal matches序列反向互补匹配的情况。这时,
    mummer将query序列转换成反向互补序列,在输出文件中,>一行后面有reverse字样,
    其结果中第三列的值也是该反向互补序列的比对位点,表示query反向互补序列从该位点
    开始往后是maximal matches序列;若加入-c参数,则将第三列的值转换成query序列
    的位点(query序列长度 - 第三列值 + 1),表示从该位点往前是maximal matches
    序列。
    
    -F
        程序默认设置下,一般会生成4列数据,第一列表示reference的序列ID,若
    reference序列仅有一条,则不会给出该列数据。加入该参数,则强制性给出该列
    数据。
    

    mummer命令的结果解读:

    1. 结果中>开头行表示某条query序列,之后的行表示该条query序列的详细比对结果;
    若该行有reverse,则表示该query序列的反向互补序列的结果。
    2. 第一列是reference序列ID;第二列表示reference起始位点;第三列表示query
    序列或其反向互补序列的起始位点;第四列表示maximal matches序列长度。
    
    对一段3000bp的序列作为reference序列,对其841-1800区段进行反向互补后,作为
    query序列
    
    1. 不加 -b 参数,则不会对query序列反向互补序列进行分析
    $ mummer -n -F ref.fa query.fa
    > query_seq
      ref_seq         1         1       840
      ref_seq      1801      1801      1140
    
    2. 加 -b 不加 -c 参数,额外给出反向互补序列的结果
    $ mummer -b -n -F ref.fa query.fa
    > query_seq
      ref_seq         1         1       840
      ref_seq      1801      1801      1140
    > query_seq Reverse
      ref_seq       841      1141       960
    注意 Reverse 表示反向互补结果,query反向互补序列从其1141位点往后共960bp属于
    maximal match。
    
    3. 加 -b 和 -c 参数
    $ mummer -b -c -n -F ref.fa query.fa
    > query_seq
      ref_seq         1         1       840
      ref_seq      1801      1801      1140
    > query_seq Reverse
      ref_seq       841      1800       960
    
    注意 Reverse 表示反向互补结果,query序列从其1800位点往前960bp属于maximal 
    match。
    

    4. nucmer命令

    nucmer命令用于对nucleotide序列进行比对。相比于mummer命令,nucmer命令能将相邻的maximal matches连起来作为cluster,然后对cluster两端进行延伸,形成大的匹配区域;并且计算SNP数量和INDEL间的距离。

    nucmer命令是MUMmer软件使用最多的命令,常用于相似性很高的核酸序列比较,特别是通一
    个物种间的基因组序列比较,或不同De novo组装软件得到的assemblies之间的比较。

    nucmer命令的使用:

    用法和常用示例:
    $ nucmer [options] ref:path query:path
    
    $ nucmer -c 200 -g 200 ref.fa query.fa
    常用参数:
    -p <string>    default: out
        设置输出文件前缀。
    --num
    --mumreference
    --maxmatch
        mummer程序用于计算获得maximal matches。程序通过suffix array算法将
    reference assembly构建索引数据库,然后再将query assembly和数据库比对,
    搜索并得到maximal matches。
        有很多maximal matches可能属于重复序列,而重复序列的比对可能没太大意义,
    且重复序列的比对极大消耗计算量。因此,需要考虑是否需要得到重复序列的比对结果。
        以上3个参数则对应重复序列的处理方法:默认是-mumreference参数,表示构建
    数据库时,去掉了reference assemly中的重复序列,得到的结果中,maximal 
    matches在reference assembly中是unique的,这种方法计算速度最快;-mum参数
    则是在前者基础上,对query assembly先去除重复序列,再进行比对,得到的结果中,
    maximal matches在reference assembly和query assembly中都是unique的,
    这种方法计算速度较慢;-maxmatch参数则完全不考虑重复序列,该方法计算速度最慢,
    同时生成的结果文件非常大。
        nucmer命令调用mummer程序进行计算,沿用了以上3个参数,只是参数的使用由
    一个中划线变成了两个中划线。
        若比较两个基因组,希望得到更大的比对区域,则考虑使用--maxmatch参数;若
    想整合两个assemblies,则最好不对重复序列进行比对,考虑使用--mumreference
    参数。
    
    -f | --forward
    -r | --reverse
        默认设置下,程序对正负链都进行比对;加入-f参数则仅对正链进行比对;加入-r
    参数则仅对负链进行比对。
    
    -l <int>    default: 20
        设置maximal match的最小长度。
    
    -g <int>    default: 90
        将两个相邻的maximal matches连起来,作为一个cluster,要求这两个matches
    之间的距离不超过该参数指定的碱基数。
    
    -c <int>    default: 65
        一个cluster中所有maximal matches的长度总和要不小于该参数值。
    
    -d <float>    default: .12
    -D <int>    default: 5
        若cluster中相邻的两个maximal matches和中间Gap区段中,变异的碱基比率
    要低于-d参数值;变异的个数要少于-D参数值。
    --extend
    --noextend
        设置是否对clustering两端进行延伸。默认设置是--extend。
    
    -b <int>    default: 200
        对cluster两端进行延伸(cluster侧翼存在一定长度的匹配区域)时,不能跨越
    长度超过200bp的低匹配得分区域。
    
    -t | --threads <int>    
        设置程序运行所使用的线程数。默认值是使用计算机所有CPU线程数。nucmer多
    线程效率不高:例如,80线程服务器运行一个nucmer命令,平均负载不到40,真实的
    CPU线程资源消耗不到20。
    
    --save <string>
        加入该参数,则能生成reference序列的索引数据库,该参数值是生成的索引数
    据库文件前缀。
    
    --load <string>
        输入reference索引数据库前缀。使用前一个参数和此参数,在多次利用nucmer
    命令进行分析时,能节约计算时间。
    

    nucmer程序输出文件为out.delta。该文件示例内容:

    /home/username/reference.fasta /home/username/query.fasta
    NUCMER
    >ref_seq1 query_seq1 500 2000000
    88 198 1641558 1641668 0 0 0
    0
    167 4877 1 4714 15 15 0
    2456
    1
    -11
    769
    950
    1
    1
    -142
    -1
    0
    >ref_seq2 query_seq4 50000 30000
    5198 22885 5389 23089 18 18 0
    -6
    -32
    -1
    -1
    -1
    7
    1130
    0
    

    out.delta文件格式解析:

    1. 文件第一行是两个需要比较的基因组序列fasta文件的绝对路径。
    2. 第二行是比对类型,可以是NUCMER或PROMER。
    3. 第三行往后是比对结果。每个比对结果以行首为 > 号开头,整行为 0 结尾。
    4. 比对结果第一行以 > 号开头,后面包含空格分割的4列,分别是:reference序列
       ID、query序列ID、该reference序列长度、该query序列长度。
    5. 比对结果第二行包含空格分割的7列数值,分别是:reference序列比对起始位点、
       reference序列比对结束位点、query序列比对起始位点、query序列比对结束位
       点、错配位点数、不相似位点数(在蛋白序列比对中有意义,核酸序列比对中和前
       一个数值一致)、非字母字符个数(在蛋白序列比对中表示终止密码子个数,在核
       酸序列比对中一般是0)
    6. 第三行往后是Delta值,表示离下一个IDNEL之间的距离,正数表示在reference
       中是Inserion,负数表示在reference中是Deletion。最后一行是0。
       例如: reference  acg.tagctgag$
                 query  .cggtag.tgag$
       Delta = (1, -3, 4, 0)。
    

    参考:
    http://www.chenlianfu.com/?p=2559

    相关文章

      网友评论

        本文标题:nucmer全基因组序列比对

        本文链接:https://www.haomeiwen.com/subject/lrjhrltx.html