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

时间:2024-01-03 18:00:50

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

在这个目录下

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

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...