• windowmasker 标记基因组中的重复序列和低复杂度序列


    下载地址:ftp://ftp.ncbi.nlm.nih.gov/pub/agarwala/windowmasker/

    在这个目录下

    其中windowmasker 为linux 平台的可执行文件

    windowmasker 只需要根据基因组序列本身,就可以识别并标记高重复序列和低复杂度序列,

    其有两种工作模式, 第一种为WinMasker模式, 用于识别重复序列;第二种为DUST模式,用于识别低复杂度序列;

    windowsmasker 的处理过程分为两步:

    第一步先生成一个count文件,通过指定 -mk_counts 开启,在这一步中会进行4次处理,其中前3步是可选的,

      pass1 检查输入文件中的duplicate reads, 通过指定-checkdup true 开启;

      pass2计算基因组总的碱基数;

      pass3 计算阈值;

      pass4 生成 count 文件;

    count文件的格式有4种可选,分别是ascii,   binary, oascii, obinary ,  默认为ascii, 通过-sformat 参数指定, ascii 格式是文本格式,是人类可读的;binary 格式是二进制格式,在第二步中load 的更快; 而oascii 和 obinary 是通过哈希优化过的结构, 能够提升2.5-4倍的运行速度,同时消耗的内存更多,他们二者之间的区别在于,oascii 人类可读, obinarty 不可读;

    备注: windowmasker 提供了一个count 文件的转换功能, 通过指定 -convert 开启, -in 指定原始的 count 文件, -out  指定生成的新的count 文件, -sformat 指定新生成文件的格式, -smen 指定生成文件时所用的最大内存, 只有输出文件格式为oascii 和 obinary 时才起作用, 其值为正整数, 默认为512, 单位数M, 即512M的内存, 当这个值不能满足要求时, 会报错;

    第二步根据count文件标记输入序列中的重复序列, 通过指定-ustat 参数开启,该参数的值是第一步生成的count 文件;在这一步中也可以加入-dust 使用DUST 模式同时标记低复杂度序列,只有同时为重复序列和低复杂度序列的 区域被标记。 这一步运行完后生成的文件有两种格式,分别是interval和fasta, 默认为interval, 通过-outfmt 格式指定,interval 的格式输出内容为标记的序列的区间, 而fasta 格式用小写字母标记序列;

     准备输入文件: 以人类hg38 的Y染色体为例:

    wget http://ftp.ncbi.nih.gov/genomes/H_sapiens/CHR_Y/hs_ref_GRCh38.p2_chrY.fa.gz
    gunzip hs_ref_GRCh38.p2_chrY.fa.gz
    sed 's/>.*/>Y/' hs_ref_GRCh38.p2_chrY.fa  > hg38.chrY.fa

    查看fasat 文件内容:

    head hg38.chrY.fa
    >Y
    CTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTCTGAAAGTGGACCTATCAGCAGGATGTGGGTG
    GGAGCAGATTAGAGAATAAAAGCAGACTGCCTGAGCCAGCAGTGGCAACCCAATGGGGTCCCTTTCCATA
    CTGTGGAAGCTTCGTTCTTTCACTCTTTGCAATAAATCTTGCTATTGCTCACTCTTTGGGTCCACACTGC
    CTTTATGAGCTGTGACACTCACCGCAAAGGTCTGCAGCTTCACTCCTGAGCCAGTGAGACCACAACCCCA
    CCAGAAAGAAGAAACTCAGAACACATCTGAACATCAGAAGAAACAAACTCCGGACGCGCCACCTTTAAGA
    ACTGTAACACTCACCGCGAGGTTCCGCGTCTTCATTCTTGAAGTCAGTGAGACCAAGAACCCACCAATTC
    CAGACACACTAGGACCCTGAGACAACCCCTAGAAGAGCACCTGGTTGATAACCCAGTTCCCATCTGGGAT
    TTAGGGGACCTGGACAGCCCGGAAAATGAGCTCCTCATCTCTAACCCAGTTCCCCTGTGGGGATTTAGGG
    GACCAGGGACAGCCCGTTGCATGAGCCCCTGGACTCTAACCCAGTTCCCTTCTGGAATTTAGGGGCCCTG

    运行第一步,得到count文件,输出格式采用obinary

    windowmasker -mk_counts -in hg38.chrY.fa -infmt fasta -out chrY.count -sformat obinary
    computing the genome length
    pass 1
    pass 2
    optimizing the data structure
    

    运行第二步,标记重复序列

    1) 生成interval 格式的输出

    windowmasker -ustat chrY.count -in hg38.chrY.fa -out hg38.chrY.masked -outfmt interval
    head hg38.chrY.masked 
    >Y
    0 - 39
    52 - 71
    167 - 202
    232 - 261
    278 - 297
    344 - 368
    377 - 416
    823 - 840
    1431 - 1452

    2) 生成fasta 格式的输出

    windowmasker -ustat chrY.count -in hg38.chrY.fa -out hg38.chrY.fasta.masked -outfmt fasta
    head hg38.chrY.fasta.masked 
    >Y
    ctaaccctaaccctaaccctaaccctaaccctaaccctctGAAAGTGGACCTatcagcag
    gatgtgggtgggAGCAGATTAGAGAATAAAAGCAGACTGCCTGAGCCAGCAGTGGCAACC
    CAATGGGGTCCCTTTCCATACTGTGGAAGCTTCGTTCTTTCACTCTTtgcaataaatctt
    gctattgctcactctttgggtccACACTGCCTTTATGAGCTGTGACACTCACcgcaaagg
    tctgcagcttcactcctgagccAGTGAGACCACAACCCcaccagaaagaagaaactcaGA
    ACACATCTGAACATCAGAAGAAACAAACTCCGGACGCGCCACCTttaagaactgtaacac
    tcaccgcgaGGTTCCGCgtcttcattcttgaagtcagtgagaccaagaacccaccaaTTC
    CAGACACACTAGGACCCTGAGACAACCCCTAGAAGAGCACCTGGTTGATAACCCAGTTCC
    CATCTGGGATTTAGGGGACCTGGACAGCCCGGAAAATGAGCTCCTCATCTCTAACCCAGT
    

    利用convert进行格式转换

    windowmasker -convert -in chrY.count -out chrY.count.ascii -sformat ascii 
    reading counts...
    converting counts...
    converting parameters...
    final processing...

      

  • 相关阅读:
    [LC] 71. Simplify Path
    [LC] 225. Implement Stack using Queues
    [Coding Made Simple / LeetCode 1235] Maximum Profit in Job Scheduling
    [Coding Made Simple] Cutting Rod for max profit
    [Coding Made Simple] Longest Common Substring
    [GeeksForGeeks] Convert an array to reduced form
    [GeeksForGeeks] Find if there is a pair with a given sum in a sorted and rotated array.
    [Coding Made Simple] Optimal Binary Search Tree
    [GeeksForGeeks] Write a program to delete a tree
    [GeeksForGeeks] Check if two trees are Isomorphic
  • 原文地址:https://www.cnblogs.com/xudongliang/p/5072961.html
Copyright © 2020-2023  润新知