a) RepeatMaskerをubuntu18.04にインストール
  b) RepeatMaskerの出力結果をRで解析する

1 RepeatMaskerインストール

1.0.2 Text::Soundex

  • perlモジュール, RepeatMaskerを実施するときにエラーでたのでインストール

1.0.5 Dfam

  • ファイルサイズ大(今回のバージョンでは解凍して83G)

1.0.8 RepeatMasker実行

  • -pa: パラレル
  • -species: human,mouse,rattus,"ciona savignyi",arabidopsis, etc.
  • -s: Slow search[0-5%高感度、デフォルトより2-3倍低速]
  • -q: Quick search[5-10%低感度、デフォルトより2-5倍高速(ほとんどの状況ではクイック検索で問題無し)]
  • -qq:ラッシュジョブ[0-5%高感度、デフォルトより4->10倍高速]
  • -nolow: Does not mask low_complexity DNA or simple repeats
  • -dir: 出力ディレクトリ. 指定しない場合currentディレクトリにできる。既に出力ファイルがある場合、古いフィアルは新しいディレクトリに移される。

2 RepeatMasker出力結果

2.1 ショートサマリーファイル(.tbl)

  • すべてのサブクラスがリストされていないため、クラスの合計はサブクラスの合計よりも多い。
  • 指定したspeciesによって表の形が異なる。サブクラスの集計はoutファイルで行う。
  • tblファイルに示されたsuperfamilyはWicker et al.2007で定義されている。
  • 10bpより短い配列はアノテーションされていないがマスクはされている。
##  [1] "SINEs:             8133      1155881 bp    9.06 %"
##  [2] "      ALUs         3028       355297 bp    2.79 %"
##  [3] "      MIRs          720        88733 bp    0.70 %"
##  [4] "LINEs:             1900       808260 bp    6.34 %"
##  [5] "      LINE1        1412       727607 bp    5.71 %"
##  [6] "      LINE2         384        67080 bp    0.53 %"
##  [7] "      L3/CR1         79        10322 bp    0.08 %"
##  [8] "LTR elements:      5510      1311103 bp   10.28 %"
##  [9] "      ERVL          492       131805 bp    1.03 %"
## [10] "      ERVL-MaLRs   2327       654365 bp    5.13 %"
## [11] "      ERV_classI    197        52951 bp    0.42 %"
## [12] "      ERV_classII  2459       449807 bp    3.53 %"
## [13] "DNA elements:      1002       202308 bp    1.59 %"
## [14] "     hAT-Charlie    633       124619 bp    0.98 %"
## [15] "     TcMar-Tigger   161        35931 bp    0.28 %"
## [16] "Unclassified:        44        15941 bp    0.13 %"
## [17] "Total interspersed repeats:  3493493 bp   27.40 %"
## [18] "Small RNA:          143        12069 bp    0.09 %"
## [19] "Satellites:         598        65894 bp    0.52 %"
## [20] "Simple repeats:       0            0 bp    0.00 %"
## [21] "Low complexity:       0            0 bp    0.00 %"

2.2 アノテーションファイル(.out)

2.2.1 アノテーションファイルをタブ区切りテキストファイルに変換する(R)。

  • 一見tab区切りテキストに見えるが違うのでreadLinesで読み込む
  • 14カラムあるのが正しいようだが、便宜上16列にする。
  • 最後列のアスタリスク（*）は、ドメインの一部（<80％）に含まれるスコアの高い一致。
  • IDが重複しているものはタンデムに並んでいるもの?
  • 1スタートみたいなのでbedに変換するとき注意
# RepeatMaskerのアノテーションファイルをタブ区切りデータに変換する
rmout <- function(in_f){
  # outファイル読み込み 
  rmlines <- readLines(in_f)
  
  # 編集 
  library(dplyr)
  rmlist <- rmlines[-1:-3] %>% # ヘーッダー部分を除外
    sapply(., function(x) strsplit(x, " {1,}")) %>%  # 空白文字で分割
    setNames(., NULL) %>% 
    lapply(., function(x){ 
      if (!x[length(x)] == "*") {　# 最後列の*の有無で分岐
        x[length(x) + 1] <- ""
      } 
      if (x[1] == "") { # 先頭の空文字除去
        x[-1]
      } else {
        x
      }
    }) 
    
  
  # データフレームに変換, class/familyを分割する, カラム名は適当 
    tag <- c("SW_score", "perc_div", "perc_del", "perc_ins", "query", "query_begin", "query_end",
           "(query_left)","strand", "matching_repeat", "class_family","repeat_begin","repeat_end",
           "(repeat_left)", "ID", "high_score")
    rmdat <- data.frame(do.call(rbind, rmlist), stringsAsFactors = F) %>%
    setNames(., tag)

  return(rmdat)

}

# RepeatMaskerのアノテーションファイルをタブ区切りデータに変換
in_rmout <- "~/db/repeat/rmout_210408/rmtest.fa.out"
rmdat <- rmout(in_rmout)

# 一部表示
knitr::kable(head(rmdat), align = 'c', caption = "repeatmasker output")

# 整形
rmmdat <- rmdat %>% mutate(len = as.integer(.$query_end) - as.integer(.$query_begin) + 1) %>% 
  select(class_family, query_begin, query_end, len) %>%
  group_by(class_family) %>% 
  summarise(total_len = sum(len), n = n()) %>% 
  mutate(class = sapply(strsplit(class_family, "/"), "[", 1), 
         family = sapply(strsplit(class_family, "/"), "[", 2)) %>%
  mutate(family = ifelse(is.na(family), "unclassified", family)) %>%
  select(class, family, total_len, n) %>% 
  arrange(class) 

# classの並びを見やすい様に変える
cl <- c("SINE", "LINE", "LTR", "DNA", "RC", "Retroposon", 
  "rRNA", "scRNA", "snRNA", "srpRNA", "tRNA", "Satellite", "Unknown")
rmmdat <- rmmdat %>% 
  mutate(class = factor(class, levels = cl)) %>%
  arrange(class)

# 一部表示
rmmdat %>% group_by(class) %>% top_n(3) %>% 
  knitr::kable(., align = 'c', caption = "class/family summarise")
## Selecting by n
repeatmasker output
SW_score perc_div perc_del perc_ins query query_begin query_end (query_left) strand matching_repeat class_family repeat_begin repeat_end (repeat_left) ID high_score
732 16.6 5.8 6.4 Contig1 70 243 (12751134) + Lx7 LINE/L1 7101 7273 (423) 1
250 15.5 5.2 0.0 Contig1 260 317 (12751060) + Lx7 LINE/L1 7617 7677 (19) 1
748 17.0 1.9 4.7 Contig1 855 1008 (12750369) + B1_Mur4 SINE/Alu 1 150 (0) 2
272 14.7 2.9 0.0 Contig1 1693 1760 (12749617) C RLTR17B_Mm LTR/ERVK (479) 323 254 3
408 13.2 2.5 5.1 Contig1 2307 2386 (12748991) C B1_Mur3 SINE/Alu (8) 141 64 4
616 14.7 1.0 1.0 Contig1 2441 2543 (12748834) C B1_Mur3 SINE/Alu (0) 149 47 5
class/family summarise
class family total_len n
SINE Alu 367380 3047
SINE B2 271798 1584
SINE B4 423806 2543
LINE CR1 10322 82
LINE L1 728349 2025
LINE L2 67117 447
LTR ERVK 463312 2877
LTR ERVL 135899 604
LTR ERVL-MaLR 668096 2849
DNA hAT-Charlie 125272 691
DNA hAT-Tip100 22356 104
DNA TcMar-Tigger 36033 177
RC Helitron 151 1
Retroposon L1-dep 13091 31
rRNA unclassified 940 11
scRNA unclassified 6352 79
snRNA unclassified 1446 13
srpRNA unclassified 643 6
tRNA unclassified 1436 24
Satellite unclassified 65380 617
Satellite Y-chromosome 1744 26
Unknown unclassified 3040 18

2.2.2 アノテーションファイルをgff3に変換(bash)

  • RepeatMaskerのユーティリティプログラムrmOutToGFF3.plを使ってgff3に変換出来る。

2.3 配列の切り出し1

  • RepeatMaskerのユーティリティプログラムrmOut2Fasta.plを使えば.outファイルからfasta配列をきりだせる。
  • NCBI RefSeqのゲノムデータの付置データとして*rm.out.gzのみが置いてある場合があるので、そのファイルを使って配列を切り出すことが出来る。
>Contig1:70-243 ( Lx7:7101-7273 orient=+ )
GCCTGCACAGATCTATCTAAGCCAGATGGGGCACCAGTGCTGAGAAGGGC
AATGGACAGAAGCCCAACCCAGAATTTAGCTCCAGTTGATAACCACTTGA
GAAYGAAATATTATTTTTCTACAACAAGAGTCTTACTGAGTCTTCTAACA
AACCACACCTAAGGGCAGGTCCCA
>Contig1:260-317 ( Lx7:7617-7677 orient=+ )
TGGGTGGGAGGAAGGAGTTGTGGGAAGGGAAACTGTGATCAGAATACATT
GTAGGAAA

2.4 配列の切り出し2

  • リピートマスカーを-gffオプションをつけて実行し、gffファイルをもとに配列を切り出す。
  • gffreadでは切り出せなかった。書き換えればできるかも
  • gff2bedでgffをbedに変換 -> bedtools getfastaで切り出した

3 環境.

## R version 3.5.2 (2018-12-20)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] ja_JP.UTF-8/ja_JP.UTF-8/ja_JP.UTF-8/C/ja_JP.UTF-8/ja_JP.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] dplyr_0.8.5
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.7        knitr_1.33        magrittr_2.0.1    tidyselect_1.1.1 
##  [5] R6_2.5.0          rlang_0.4.11      fansi_0.5.0       highr_0.9        
##  [9] stringr_1.4.0     tools_3.5.2       xfun_0.24         utf8_1.2.2       
## [13] jquerylib_0.1.4   htmltools_0.5.1.1 ellipsis_0.3.2    yaml_2.2.1       
## [17] digest_0.6.27     assertthat_0.2.1  tibble_3.1.3      lifecycle_1.0.0  
## [21] crayon_1.4.1      purrr_0.3.4       vctrs_0.3.8       sass_0.4.0       
## [25] glue_1.4.2        evaluate_0.14     rmarkdown_2.9     stringi_1.7.3    
## [29] compiler_3.5.2    bslib_0.2.5.1     pillar_1.6.1      jsonlite_1.7.2   
## [33] pkgconfig_2.0.3