6a48e07c2fbefbf0a0e5e74662dd768ef93fca88 hiram Wed Apr 29 11:27:30 2020 -0700 adding ncbiRmsk option to use NCBI rm.out.gz file instead of cluster run and record versionInfo.txt of RepeatMasker operation refs #23734 diff --git src/hg/utils/automation/doRepeatMasker.pl src/hg/utils/automation/doRepeatMasker.pl index 286e999..8e05e16 100755 --- src/hg/utils/automation/doRepeatMasker.pl +++ src/hg/utils/automation/doRepeatMasker.pl @@ -15,74 +15,85 @@ use HgRemoteScript; use HgStepManager; # Hardcoded command path: my $RepeatMaskerPath = "/hive/data/staging/data/RepeatMasker"; my $RepeatMasker = "$RepeatMaskerPath/RepeatMasker"; my $RepeatMaskerEngine = "-engine crossmatch -s"; # Let parasol pick defaults my $parasolRAM = ""; # Option variable names, both common and peculiar to this script: use vars @HgAutomate::commonOptionVars; use vars @HgStepManager::optionVars; use vars qw/ $opt_buildDir + $opt_ncbiRmsk + $opt_liftSpec $opt_species $opt_unmaskedSeq $opt_customLib $opt_useHMMER $opt_useRMBlastn $opt_splitTables $opt_noSplit $opt_updateTable /; # Specify the steps supported with -continue / -stop: my $stepper = new HgStepManager( [ { name => 'cluster', func => \&doCluster }, { name => 'cat', func => \&doCat }, { name => 'mask', func => \&doMask }, { name => 'install', func => \&doInstall }, { name => 'cleanup', func => \&doCleanup }, ] ); # Option defaults: my $dbHost = 'hgwdev'; my $unmaskedSeq = "\$db.unmasked.2bit"; my $defaultSpecies = 'scientificName from dbDb'; +my $ncbiRmsk = ""; # path to NCBI *_rm.out.gz file to use NCBI result + # for repeat masking result + # in the assembly as calculated by NCBI +my $liftSpec = ""; # lift file from NCBI coordinates to UCSC coordinates + # usually required with ncbiRmsk my $base = $0; $base =~ s/^(.*\/)?//; sub usage { # Usage / help / self-documentation: my ($status, $detailed) = @_; # Basic help (for incorrect usage): print STDERR " usage: $base db options: "; print STDERR $stepper->getOptionHelp(); print STDERR <<_EOF_ -buildDir dir Use dir instead of default $HgAutomate::clusterData/\$db/$HgAutomate::trackBuild/RepeatMasker.\$date (necessary when continuing at a later date). -species sp Use sp (which can be quoted genus and species, or a common name that RepeatMasker recognizes. Default: $defaultSpecies. + -ncbiRmsk path/file_rm.out.gz - Use the repeat masker result supplied in + the assembly as calculated by NCBI + -liftSpec path/file.lift - Use this lift file to lift the NCBI coordinates + to UCSC coordinates, usually used with ncbiRmsk -unmaskedSeq seq.2bit Use seq.2bit as the unmasked input sequence instead of default ($unmaskedSeq). -customLib lib.fa Use custom repeat library instead of RepeatMaskers\'s. -useRMBlastn Use NCBI rmblastn instead of crossmatch -useHMMER Use hmmer instead of crossmatch ( currently for human only ) -updateTable load into table name rmskUpdate (default: rmsk) -splitTables split the _rmsk tables (default is not split) -noSplit default behavior, this option no longer required. _EOF_ ; print STDERR &HgAutomate::getCommonOptionHelp('dbHost' => $dbHost, 'workhorse' => '', 'bigClusterHub' => ''); print STDERR " Automates UCSC's RepeatMasker process for genome database \$db. Steps: @@ -107,51 +118,63 @@ exit $status; } # Globals: # Command line args: db my ($db); # Other: my ($buildDir, $chromBased, $updateTable, $secondsStart, $secondsEnd); sub checkOptions { # Make sure command line options are valid/supported. my $ok = GetOptions(@HgStepManager::optionSpec, 'buildDir=s', 'species=s', + 'ncbiRmsk=s', + 'liftSpec=s', 'unmaskedSeq=s', 'customLib=s', 'useRMBlastn', 'useHMMER', 'splitTables', 'noSplit', 'updateTable', @HgAutomate::commonOptionSpec, ); &usage(1) if (!$ok); &usage(0, 1) if ($opt_help); &HgAutomate::processCommonOptions(); my $err = $stepper->processOptions(); usage(1) if ($err); $dbHost = $opt_dbHost if ($opt_dbHost); } ######################################################################### # * step: cluster [bigClusterHub] sub doCluster { + return if ($ncbiRmsk); my $runDir = "$buildDir/run.cluster"; + if ( -d "$runDir" ) { + if ( -s "$runDir/run.time" ) { + &HgAutomate::verbose(1, + "\ncluster step previously completed\n"); + return; + } else { + die "\nERROR: cluster step may not be completed properly\n"; + } + } &HgAutomate::mustMkdir($runDir); my $paraHub = $opt_bigClusterHub ? $opt_bigClusterHub : &HgAutomate::chooseClusterByBandwidth(); if (! -e $unmaskedSeq) { die "Error: required file $unmaskedSeq does not exist."; } my @okIn = grep !/scratch/, &HgAutomate::chooseFilesystemsForCluster($paraHub, "in"); my @okOut = &HgAutomate::chooseFilesystemsForCluster($paraHub, "out"); if (scalar(@okOut) > 1) { @okOut = grep !/$okIn[0]/, @okOut; } my $inHive = 0; $inHive = 1 if ($okIn[0] =~ m#/hive/data/genomes#); @@ -277,58 +300,73 @@ ; close($fh); &HgAutomate::makeGsub($runDir, "./RMRun.csh {check out line $partDir/\$(path1).out}"); my $whatItDoes = "It computes a logical partition of unmasked 2bit into 500k chunks and runs it on the cluster with the most available bandwidth."; my $bossScript = new HgRemoteScript("$runDir/doCluster.csh", $paraHub, $runDir, $whatItDoes); $bossScript->add(<<_EOF_ chmod a+x dummyRun.csh chmod a+x RMRun.csh -./dummyRun.csh # Record RM version used: +printf "The repeat files provided for this assembly were generated using RepeatMasker.\\ + Smit, AFA, Hubley, R & Green, P.,\\ + RepeatMasker Open-4.0.\\ + 1996-2010 <http://www.repeatmasker.org>.\\ +\\ +VERSION:\\n" > ../versionInfo.txt + +./dummyRun.csh | grep -v "dev/null" >> ../versionInfo.txt + +grep 'version of RepeatMasker\$' $RepeatMasker >> ../versionInfo.txt +grep RELEASE $RepeatMaskerPath/Libraries/RepeatMaskerLib.embl >> ../versionInfo.txt +printf "# RepeatMasker engine: %s\\n" "${RepeatMaskerEngine}" >> ../versionInfo.txt + ls -ld $RepeatMaskerPath $RepeatMasker grep 'version of RepeatMasker\$' $RepeatMasker grep RELEASE $RepeatMaskerPath/Libraries/RepeatMaskerLib.embl echo "# RepeatMasker engine: $RepeatMaskerEngine" _EOF_ ); if ($opt_useRMBlastn) { $bossScript->add(<<_EOF_ +printf "# using rmblastn:\\n" >> ../versionInfo.txt echo "# useRMBlastn: rmblastn:" -grep RMBLAST_DIR $RepeatMaskerPath/RepeatMaskerConfig.pm | awk '{print \$NF}' +grep RMBLAST_DIR $RepeatMaskerPath/RepeatMaskerConfig.pm | grep rmblastn-2 | awk '{print \$NF}' >> ../versionInfo.txt _EOF_ ); } if ($opt_useHMMER) { $bossScript->add(<<_EOF_ +printf "# using Dfam library and HMMER3:\\n" >> ../versionInfo.txt echo "# useHMMER: Dfam library: " ls -ld $RepeatMaskerPath/Libraries/Dfam.hmm -grep Release: $RepeatMaskerPath/Libraries/Dfam.hmm +grep Release: $RepeatMaskerPath/Libraries/Dfam.hmm >> ../versionInfo.txt echo "# useHMMER: HMMER3: " -grep -m 1 ^HMMER3 $RepeatMaskerPath/Libraries/Dfam.hmm +grep -m 1 ^HMMER3 $RepeatMaskerPath/Libraries/Dfam.hmm >> ../versionInfo.txt _EOF_ ); } if (length($repeatLib) > 0) { $bossScript->add(<<_EOF_ +printf "# RepeatMasker library options: %s\\n" "${repeatLib}" >> ../versionInfo.txt echo "# RepeatMasker library options: '$repeatLib'" _EOF_ ); } if ( ! $opt_unmaskedSeq && ! $inHive) { $bossScript->add(<<_EOF_ mkdir -p $clusterSeqDir rsync -av $unmaskedSeq $clusterSeq _EOF_ ); } if ($opt_unmaskedSeq) { $bossScript->add(<<_EOF_ rm -rf $partDir $Bin/simplePartition.pl $clusterSeq 500000 $partDir @@ -339,107 +377,147 @@ rm -rf $partDir $Bin/simplePartition.pl $clusterSeq 500000 $partDir rm -f $buildDir/RMPart ln -s $partDir $buildDir/RMPart _EOF_ ); } my $binPara = "/parasol/bin/para"; if ( ! -s "$binPara" ) { # allow PATH to find the para command $binPara = "para"; } my $gensub2 = &HgAutomate::gensub2(); $bossScript->add(<<_EOF_ +printf "\\nPARAMETERS:\\ +$RepeatMasker $RepeatMaskerEngine -align $repeatLib\\n" >> ../versionInfo.txt + $gensub2 $partDir/partitions.lst single gsub jobList +exit 0 $binPara $parasolRAM make jobList $binPara check $binPara time > run.time cat run.time _EOF_ ); if (! $opt_unmaskedSeq && ! $inHive) { $bossScript->add(<<_EOF_ rm -f $clusterSeq _EOF_ ); } $bossScript->execute(); } # doCluster ######################################################################### # * step: cat [fileServer] sub doCat { my $runDir = "$buildDir"; + if (! -s "${ncbiRmsk}" ) { &HgAutomate::checkExistsUnlessDebug('cluster', 'cat', "$buildDir/run.cluster/run.time"); - + } + if ( -s "$runDir/doCat.bash" ) { + if ( ! (-s "$runDir/$db.sorted.fa.out.gz" || -s "$runDir/$db.sorted.fa.out" ) ) { + die "\nERROR: cat step may not be completed properly\n"; + } else { + &HgAutomate::verbose(1, + "\ncat step previously completed\n"); + return; + } + } my $whatItDoes = "It concatenates .out files from cluster run into a single $db.sorted.fa.out.\n" . "liftUp (with no lift specs) is used to concatenate .out files because it\n" . "uniquifies (per input file) the .out IDs which can then be used to join\n" . "fragmented repeats in the Nested Repeats track."; my $fileServer = &HgAutomate::chooseFileServer($runDir); my $bossScript = newBash HgRemoteScript("$runDir/doCat.bash", $fileServer, $runDir, $whatItDoes); - # Use symbolic link created in cluster step: - my $partDir = "$buildDir/RMPart"; - my $indent = ""; - my $path = ""; - my $levels = $opt_debug ? 1 : `cat $partDir/levels`; - chomp $levels; $bossScript->add("# Debug mode -- don't know the real number of levels\n") if ($opt_debug); + if ($ncbiRmsk) { + my $zippedSource = 0; # asssume not a gz zipped file + $zippedSource = 1 if ($ncbiRmsk =~ m/.gz$/); + my $liftOpts = ""; + if (-s "${liftSpec}" ) { + $liftOpts = "| liftUp -type=.out stdout $liftSpec carry stdin"; + } + $bossScript->add(<<_EOF_ +export db="${db}" +export ncbiRmOutFile="${ncbiRmsk}" +cat /dev/null > "\${db}.sorted.fa.out" +_EOF_ + ); + if ($zippedSource) { + $bossScript->add(<<_EOF_ +zcat "\${ncbiRmOutFile}" | tail -n +4 | sort -k5,5 -k6,6n $liftOpts >> "\${db}.sorted.fa.out" +_EOF_ + ); + } else { + $bossScript->add(<<_EOF_ +tail -n +4 "\${ncbiRmOutFile}" sort -k5,5 -k6,6n $liftOpts >> "\${db}.sorted.fa.out" +_EOF_ + ); + } + } else { # not using the NCBI supplied file, use local cluster run result + my $partDir = "$buildDir/RMPart"; + my $levels = $opt_debug ? 1 : `cat $partDir/levels`; + chomp $levels; + # Use symbolic link created in cluster step: + my $indent = ""; + my $path = ""; # Make the appropriate number of nested foreach's to match the number # of levels from the partitioning. If $levels is 1, no foreach required. for (my $l = $levels - 2; $l >= 0; $l--) { my $dir = ($l == ($levels - 2)) ? $partDir : "\$d" . ($l + 1); $bossScript->add(<<_EOF_ ${indent}for d$l in $dir/??? do _EOF_ ); $indent .= ' '; $path .= "\$d$l/"; } for (my $l = 0; $l <= $levels - 2; $l++) { $bossScript->add(<<_EOF_ ${indent} liftUp ${path}cat.out /dev/null carry ${path}???/*.out ${indent} liftUp ${path}cat.align /dev/null carry ${path}???/*.align ${indent}done _EOF_ ); $indent =~ s/ //; $path =~ s/\$d\d+\/$//; } $bossScript->add(<<_EOF_ liftUp $db.fa.out /dev/null carry $partDir/???/*.out liftUp $db.fa.align /dev/null carry $partDir/???/*.align head -3 $db.fa.out > $db.sorted.fa.out tail -n +4 $db.fa.out | sort -k5,5 -k6,6n >> $db.sorted.fa.out _EOF_ ); + } $bossScript->add(<<_EOF_ # Use the ID column to join up fragments of interrupted repeats for the # Nested Repeats track. Try to avoid the Undefined id errors. -($Bin/extractNestedRepeats.pl $db.fa.out 2> nestRep.err || true) | sort -k1,1 -k2,2n > $db.nestedRepeats.bed +($Bin/extractNestedRepeats.pl $db.sorted.fa.out 2> nestRep.err || true) | sort -k1,1 -k2,2n > $db.nestedRepeats.bed if [ -s "nestRep.err" ]; then export lineCount=`(grep "Undefined id, line" nestRep.err || true) | cut -d' ' -f6 | wc -l` if [ "\${lineCount}" -gt 0 ]; then if [ "\${lineCount}" -gt 10 ]; then printf "ERROR: too many Undefined id lines (> 10) reported by extractNestedRepeats.pl" 1>&2 exit 255 fi export sedExpr=`grep "Undefined id, line" nestRep.err | cut -d' ' -f6 | sed -e 's/\$/d;/;' | xargs echo` export sedCmd="sed -i.broken -e '\$sedExpr'" eval \$sedCmd $db.fa.out ($Bin/extractNestedRepeats.pl $db.fa.out 2> nestRep2.err || true) | sort -k1,1 -k2,2n > $db.nestedRepeats.bed if [ -s "nestRep2.err" ]; then printf "ERROR: the attempt of cleaning nestedRepeats did not work" 1>&2 exit 255 fi @@ -584,70 +662,72 @@ | splitFileByColumn -col=5 stdin /cluster/data/$db -chromDirs \\ -ending=.fa.out -head=/tmp/rmskHead.txt _EOF_ ); $bossScript->execute(); } } # doInstall ######################################################################### # * step: cleanup [fileServer] sub doCleanup { my $runDir = "$buildDir"; my $whatItDoes = "It cleans up or compresses intermediate files."; my $fileServer = &HgAutomate::chooseFileServer($runDir); - my $bossScript = new HgRemoteScript("$runDir/doCleanup.csh", $fileServer, + my $bossScript = newBash HgRemoteScript("$runDir/doCleanup.bash", $fileServer, $runDir, $whatItDoes); $bossScript->add(<<_EOF_ -if (-e $db.fa.align) then +if [ -e "$db.fa.align" ]; then gzip $db.fa.align -endif +fi rm -fr RMPart/* rm -fr RMPart -if ( -d /hive/data/genomes/$db/RMPart ) then +if [ -d "/hive/data/genomes/$db/RMPart" ]; then rmdir /hive/data/genomes/$db/RMPart -endif +fi _EOF_ ); $bossScript->execute(); } # doCleanup ######################################################################### # main # Prevent "Suspended (tty input)" hanging: &HgAutomate::closeStdin(); # Make sure we have valid options and exactly 1 argument: &checkOptions(); &usage(1) if (scalar(@ARGV) != 1); $secondsStart = `date "+%s"`; chomp $secondsStart; ($db) = @ARGV; # Now that we know the $db, figure out our paths: my $date = `date +%Y-%m-%d`; chomp $date; $buildDir = $opt_buildDir ? $opt_buildDir : "$HgAutomate::clusterData/$db/$HgAutomate::trackBuild/RepeatMasker.$date"; $unmaskedSeq = $opt_unmaskedSeq ? $opt_unmaskedSeq : "$HgAutomate::clusterData/$db/$db.unmasked.2bit"; my $seqCount = `twoBitInfo $unmaskedSeq stdout | wc -l`; $chromBased = ($seqCount <= $HgAutomate::splitThreshold) && $opt_splitTables; $updateTable = $opt_updateTable ? "Update" : ""; +$ncbiRmsk = $opt_ncbiRmsk ? $opt_ncbiRmsk : ""; +$liftSpec = $opt_liftSpec ? $opt_liftSpec : ""; # Do everything. $stepper->execute(); # Tell the user anything they should know. my $stopStep = $stepper->getStopStep(); my $upThrough = ($stopStep eq 'cleanup') ? "" : " (through the '$stopStep' step)"; $secondsEnd = `date "+%s"`; chomp $secondsEnd; my $elapsedSeconds = $secondsEnd - $secondsStart; my $elapsedMinutes = int($elapsedSeconds/60); $elapsedSeconds -= $elapsedMinutes * 60;