RepeatMasker

Example procedure to operate RepeatMasker on a typical genome assembly sequence.

Partition

The source tree script simplePartition.pl can be used to partition the genome. The output of the script is merely a listing of twoBit file references representing the partitions. A specific bit of sequence is obtained from the twoBit file at the time it is needed for the RepeatMasker operation. UCSC uses partitions of 500,000 bases when using the cross_match alignment engine with RepeatMasker. That size of partition often results in a single job run-time of about 1.5 hours. When using the rmblast alignment engine with RepeatMasker, the typical single job run-time is on the order of 10 minutes. The resulting masking is very similar. This example constructs partition lists in the specified RMPart directory hierarchy:

mkdir -p /data/genomes/ricCom1/bed/repeatMasker/
cd /data/genomes/ricCom1/bed/repeatMasker
simplePartition.pl ../../ricCom1.unmasked.2bit 500000 /data/genomes/ricCom1/RMPart

Cluster run

The template file used here:

#LOOP
./RMRun.csh /data/genomes/ricCom1/RMPart/$(path1).out
#ENDLOOP

with the gensub2 command constructs the command list jobList to run on the cluster.

mkdir -p /data/genomes/ricCom1/bed/repeatMasker/run.cluster
cd /data/genomes/ricCom1/bed/repeatMasker/run.cluster
ln -s /data/genomes/ricCom1/RMPart ./RMPart
gensub2/data/genomes/ricCom1/RMPart/partitions.lst single template jobList

This example constructs 880 commands to run in the jobList file, typically:

./RMRun.csh /data/genomes/ricCom1/RMPart/000/000.lst.out
./RMRun.csh /data/genomes/ricCom1/RMPart/001/001.lst.out
./RMRun.csh /data/genomes/ricCom1/RMPart/002/002.lst.out
... etc.

The RMRun.csh script performs the RepeatMasker operation on each partition:

#!/bin/csh -ef

set finalOut = $1

set inLst = $finalOut:r
set inLft = $inLst:r.lft
set alignOut = $finalOut:r.align
set catOut = $finalOut:r.cat

# Use local disk for output, and move the final result to $outPsl
# when done, to minimize I/O.
set tmpDir = `mktemp -d -p /scratch/tmp doRepeatMasker.cluster.XXXXXX`
pushd $tmpDir

# Initialize local library
RepeatMasker -engine rmblast -pa 1 -species 'Ricinus communis' /dev/null

foreach spec (`cat $inLst`)
  # Remove path and .2bit filename to get just the seq:start-end spec:
  set base = `echo $spec | sed -r -e 's/^[^:]+://'`
  # If $spec is the whole sequence, twoBitToFa removes the :start-end part,
  # which causes liftUp to barf later.  So tweak the header back to
  # seq:start-end for liftUp's sake:
  twoBitToFa $spec stdout \
  | sed -e "s/^>.*/>$base/" > $base.fa
  RepeatMasker -engine rmblast -pa 1 -align -species 'Ricinus communis' $base.fa
  if (-e $base.fa.cat) then
    mv $base.fa.cat $catOut
  endif
end

# Lift up (leave the RepeatMasker header in place because we'll liftUp
# again later):
liftUp -type=.out stdout $inLft error *.fa.out > tmpOut__out

set nonomatch
set alignFiles = ( *.align )
if ( ${#alignFiles} && -e $alignFiles[1] ) then
  liftUp -type=.align stdout $inLft error *.align \
    > tmpOut__align
else
  touch tmpOut__align
endif

# Move final result into place:
mv tmpOut__out $finalOut
mv tmpOut__align $alignOut

popd
rm -rf $tmpDir

If you have a parasol cluster, this operation is performed.

para -ram=2g make jobList
para check
para time > run.time
cat run.time
Completed: 880 of 880 jobs
CPU time in finished jobs:     202665s    3377.75m    56.30h    2.35d  0.006 y
IO & Wait Time:                119279s    1987.98m    33.13h    1.38d  0.004 y
Average job time:                 366s       6.10m     0.10h    0.00d
Longest finished job:            6789s     113.15m     1.89h    0.08d
Submission to last job:          6938s     115.63m     1.93h    0.08d

Collect cluster run results

Gather the partition results into a single file:

cd /data/genomes/ricCom1/bed/repeatMasker
liftUp ricCom1.fa.out /dev/null carry /data/genomes/ricCom1/bed/repeatMasker/RMPart/???/*.out
liftUp ricCom1.fa.align /dev/null carry /data/genomes/ricCom1/bed/repeatMasker/RMPart/???/*.align
head -3 ricCom1.fa.out > ricCom1.sorted.fa.out
tail -n +4 ricCom1.fa.out | sort -k5,5 -k6,6n >> ricCom1.sorted.fa.out

In addition, using the UCSC script: extractNestedRepeats.pl for the nestedRepeats track:

# Use the ID column to join up fragments of interrupted repeats for the
# Nested Repeats track.
extractNestedRepeats.pl ricCom1.fa.out | sort -k1,1 -k2,2n > ricCom1.nestedRepeats.bed

Masking sequence

The unmasked.2bit file is masked with the RepeatMasker results:

cd /data/genomes/ricCom1/bed/repeatMasker
twoBitMask /data/genomes/ricCom1/ricCom1.unmasked.2bit ricCom1.sorted.fa.out ricCom1.rmsk.2bit

Measure the masking:

twoBitToFa ricCom1.rmsk.2bit stdout | faSize stdin > faSize.rmsk.txt

Results in the faSize.rmsk.txt file:

350621860 bases (13662715 N's 336959145 real 326566002 upper 10393143 lower) in 25763 sequences in 1 files
Total size: mean 13609.5 sd 106411.1 min 202 (EQ999533) max 4693355 (EQ973772) median 1094
%2.96 masked total, %3.08 masked real

This example indicates a very low masking result from RepeatMasker. The Window Masker result masks the sequence much better and will be used for the eventual final masking.

Convert RM .out to bigBed

The standard RepeatMasker track in the genome browser is a special unique track type that displays several different aspects of the RepeatMasker data in side by side displays. Such a track type is not available in a track hub, but an approximation of the data display can be displayed with several different bigBed tracks in the track hub.

First, convert the RepeatMasker .out file format to the UCSC track table format:

 hgLoadOut -tabFile=ricCom1.rmsk.tab -nosplit test ricCom1.sorted.fa.out

Split the resulting table rows by their RepeatMasker class attribute:

mkdir rmskClass
sort -k12,12 ricCom1.rmsk.tab | splitFileByColumn -ending=tab -col=12 -tab stdin rmskClass

This produces a number of files in the directory ./rmskClass/ - for example

ls -ogrt rmskClass
total 12096
-rw-rw-r-- 1  694267 Apr 18 10:31 DNA.tab
-rw-rw-r-- 1    8842 Apr 18 10:31 DNA?.tab
-rw-rw-r-- 1  129866 Apr 18 10:31 LINE.tab
-rw-rw-r-- 1    6208 Apr 18 10:31 LINE?.tab
-rw-rw-r-- 1  690106 Apr 18 10:31 LTR.tab
-rw-rw-r-- 1  857390 Apr 18 10:31 Low_complexity.tab
-rw-rw-r-- 1   14911 Apr 18 10:31 Other.tab
... etc.

Convert those files in the same groupings as the standard UCSC RepeatMasker class displays to bigBed file types:

for T in SINE LINE LTR DNA Simple Low_complexity Satellite
do
 ./toBed6+10.pl rmskClass/${T}*.tab | sort -k1,1 -k2,2n > ${DB}.rmsk.${T}.bed
done
./toBed6+10.pl rmskClass/rRNA.tab | sort -k1,1 -k2,2n > ${DB}.rmsk.RNA.bed
./toBed6+10.pl rmskClass/Other.tab rmskClass/RC*.tab | sort -k1,1 -k2,2n > ricCom1.rmsk.Other.bed
for T in DNA LINE LTR Low_complexity Other RNA SINE Satellite Simple
do
  bedToBigBed -tab -type=bed6+10 -as=rmsk16.as ricCom1.rmsk.${T}.bed \
   ../../chrom.sizes ricCom1.rmsk.${T}.bb
done

Where the script toBed6+10.pl is:

#!/usr/bin/env perl

use strict;
use warnings;

my $sixteenSamples = 0;
my $fifteenSamples = 0;

my $argc = scalar(@ARGV);
if ($argc < 1) {
  printf STDERR "usage: toBed6+10.pl rmskClass/file.tab > ricCom1.class.tab\n";
  printf STDERR "e.g.: ./toBed6+10.pl rmskClass/rRNA.tab > ricCom1.rmsk.RNA.tab\n";
  exit 255;
}

while (my $file = shift) {

open (FH, "<$file") or die "can not read $file";
while (my $line = <FH>) {
  next if ($line =~ m/^\s+SW\s+perc.*/);
  next if ($line =~ m/^score\s+div.*/);
  next if ($line =~ m/^$/);
  chomp $line;
  my @a = split('\t', $line);
  shift @a;
  die "ERROR: not finding 16 fields in a line" if (scalar(@a) != 16);
  printf "%s\t%d\t%d\t%s\t0\t%s\t%d\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n",
    $a[4], $a[5], $a[6], $a[9], $a[8], $a[0],
    $a[1], $a[2], $a[3], $a[7], $a[10], $a[11], $a[12], $a[13], $a[14];
}
close (FH);
}

And the rmsk16.as bigBed structure definition is:

table repeatMasker
"RepeatMasker .out record"
    (
    string chrom;    "Genomic sequence name"
    uint chromStart;     "Start in genomic sequence"
    uint chromEnd;       "End in genomic sequence"
    string name;     "Name of repeat"
    uint score;       "always 0 place holder"
    char[1] strand;     "Relative orientation + or -"
    uint swScore;       "Smith Waterman alignment score"
    uint milliDiv;      "Base mismatches in parts per thousand"
    uint milliDel;      "Bases deleted in parts per thousand"
    uint milliIns;      "Bases inserted in parts per thousand"
    int  genoLeft;      "-#bases after match in genomic sequence"
    string repClass;    "Class of repeat"
    string repFamily;   "Family of repeat"
    int repStart;       "Start (if strand is +) or -#bases after match (if strand is -) in repeat sequence"
    uint repEnd;        "End in repeat sequence"
    int repLeft;        "-#bases after match (if strand is +) or start (if strand is -) in repeat sequence"
    )

Composite RMSK trackDb.txt

Please note the composite track arrangement to display the resulting class files in a fashion resembling the standard UCSC RepeatMasker track display:

http://genome-test.soe.ucsc.edu/~hiram/hubs/Plants/araTha1/trackDb.txt

In the section starting:

track repeatMasker_
compositeTrack on
shortLabel RepeatMasker
longLabel Repeating Elements by RepeatMasker
group varRep
priority 149.1
visibility dense
type bed 3 .
noInherit on
html ../trackDescriptions/repeatMasker

See also

Please note specific details for running each of the masking operations:

• Window Masker
• TRF Simple Repeats