6dcfa7d1962d9e9dbfa5cbe651b6e301101287ef hiram Fri Sep 25 12:51:58 2026 -0700 beginning to get a miniMap2 process in place refs #34360 diff --git src/hg/utils/automation/doMiniMap2.pl src/hg/utils/automation/doMiniMap2.pl new file mode 100755 index 00000000000..779ad99ca2f --- /dev/null +++ src/hg/utils/automation/doMiniMap2.pl @@ -0,0 +1,785 @@ +#!/usr/bin/env perl + +# DO NOT EDIT the /cluster/bin/scripts copy of this file -- +# edit ~/kent/src/hg/utils/automation/doMiniMap2.pl instead. + +# doMiniMap2.pl -- same-species / same-haplotype liftOver pipeline, +# modeled directly on doSameSpeciesLiftOver.pl, but using minimap2 +# instead of blat as the alignment engine. +# +# Intended use: two very similar same-species assemblies where blat's +# -fastMap query-chunking dance isn't needed and minimap2's splice-free +# asm5/asm10/asm20 presets do a better job with the larger indels and +# structural differences you see between the two haplotypes of a +# diploid (trio-binned or hifiasm/verkko dual-assembly) genome, or +# between two closely related strain assemblies. +# +# Requires the kent command pafToPsl (PAF+cigar -> PSL) so the existing +# axtChain/chainNet/netChainSubset toolchain can be reused unchanged. + +use Getopt::Long; +use warnings; +use strict; +use FindBin qw($Bin); +use lib "$Bin"; +use HgAutomate; +use HgRemoteScript; +use HgStepManager; +use AssemblyDivergence qw(mashDistance choosePipeline + $mashAsm5Max $mashAsm10Max $mashLastzMin $mashWarnMax); + +# Option variable names, both common and peculiar to this script: +use vars @HgAutomate::commonOptionVars; +use vars @HgStepManager::optionVars; +use vars qw/ + $opt_buildDir + $opt_target2Bit + $opt_targetSizes + $opt_query2Bit + $opt_querySizes + $opt_minimapPreset + $opt_minimapCpu + $opt_chainRam + $opt_chainCpu + $opt_regenerateMash + /; + +# Specify the steps supported with -continue / -stop: +my $stepper = new HgStepManager( + [ { name => 'align', func => \&doAlign }, + { name => 'chain', func => \&doChain }, + { name => 'net', func => \&doNet }, + { name => 'load', func => \&doLoad }, + { name => 'cleanup', func => \&doCleanup }, + ] + ); + +# Option defaults: +my $dbHost = 'hgwdev'; +my $ramG = '32g'; # minimap2 index + 8 threads on a whole genome needs headroom +my $cpu = 1; +# minimapPreset is normally left undef and picked automatically by +# estimateDivergence() (mash distance -> asm5/asm10/asm20); -minimapPreset +# overrides that and skips the mash run entirely. +my $minimapPreset; +my $minimapCpu = 8; # -t N threads given to each minimap2 job +my $chainRam = '16g'; # -chainRam=Ng argument +my $chainCpu = 1; # -chainCpu=N argument +# mash distance thresholds ($mashAsm5Max/$mashAsm10Max/$mashLastzMin/ +# $mashWarnMax) used to pick a preset (or reject the pair outright) come +# from AssemblyDivergence.pm, shared with mashDistance.pl and anything +# else that needs to triage a pair between this pipeline and +# doBlastzChainNet.pl. + +my $base = $0; +$base =~ s/^(.*\/)?//; + +sub usage { + # Usage / help / self-documentation: + my ($status, $detailed) = @_; + # Basic help (for incorrect usage): + print STDERR " +usage: $base fromDb toDb +options: +"; + print STDERR $stepper->getOptionHelp(); + print STDERR <<_EOF_ + -buildDir dir Use dir instead of default + $HgAutomate::clusterData/\$fromDb/$HgAutomate::trackBuild/mm2.\$toDb.\$date + (necessary when continuing at a later date). + -target2Bit /path/target.2bit Full path to target sequence (fromDb) + -query2Bit /path/query.2bit Full path to query sequence (toDb) + -targetSizes /path/target.chrom.sizes Full path to target chrom.sizes (fromDb) + -querySizes /path/query.chrom.sizes Full path to query chrom.sizes (toDb) + -minimapPreset asm5|asm10|asm20 minimap2 -x preset to use. By default + this is chosen automatically via AssemblyDivergence.pm + (same logic as the standalone mashDistance.pl): mash + sketches \$fromDb and \$toDb and picks a preset from + the mash distance between them (< $mashAsm5Max -> asm5, + < $mashAsm10Max -> asm10, < $mashLastzMin -> asm20). + At or above $mashLastzMin this script refuses to run + -- that pair looks too diverged for minimap2's asm* + presets; use doBlastzChainNet.pl instead, or give + -minimapPreset to skip the mash run and force a + specific preset anyway. + -minimapCpu N Threads given to each minimap2 cluster job (-t N), + default: $minimapCpu + -chainRam Ng Cluster ram size for chain step, default: -chainRam=$chainRam + -chainCpu N Cluster CPUs number for chain step, default: -chainCpu=$chainCpu + -regenerateMash Force the mash divergence check to re-sketch and + overwrite even if a cached .msh (GenArk mashSketch/ + cache or \$buildDir/mashDistance.txt) already + exists. Ignored if -minimapPreset is also given, + since that skips the mash run entirely. +_EOF_ + ; + print STDERR &HgAutomate::getCommonOptionHelp('dbHost' => $dbHost, + 'workhorse' => '', + 'fileServer' => '', + 'ram' => $ramG, + 'cpu' => $cpu, + 'bigClusterHub' => ''); + print STDERR " +Automates a same-species/same-haplotype liftOver (minimap2/chain/net) +pipeline, patterned after doSameSpeciesLiftOver.pl but using minimap2 +in place of blat -fastMap: + align: Aligns the assemblies using minimap2 -cx \$minimapPreset on a + big cluster, one job per target sequence against the full + query genome, then converts PAF+cigar to PSL with pafToPsl. + chain: Chains the alignments on a big cluster. + net: Nets the alignments, uses netChainSubset to extract liftOver chains. + load: Installs liftOver chain files, calls hgAddLiftOverChain on $dbHost. + cleanup: Removes or compresses intermediate files. +All operations are performed in the build directory which is +$HgAutomate::clusterData/\$fromDb/$HgAutomate::trackBuild/mm2.\$toDb.\$date unless -buildDir is given. +"; + # Detailed help (-help): + print STDERR " +Assumptions: +1. /scratch/data/\$db/\$db.2bit contains sequence for database/assembly \$db. +2. $HgAutomate::clusterData/\$db/chrom.sizes contains all sequence names and sizes from + \$db.2bit. +3. The \$db.2bit files have already been distributed to cluster-scratch + (/scratch/data/<db>/). +4. pafToPsl is on \$PATH on the machine that runs the align cluster jobs. +5. mash is on \$PATH on the machine that runs this script (used once, up + front, to auto-select -minimapPreset unless it is given explicitly). +6. fromDb and toDb are two haplotypes/assemblies of the same species/ + individual -- this is not a general any-vs-any pipeline. For that, + use doBlastzChainNet.pl. +" if ($detailed); + print "\n"; + exit $status; +} + + +# Globals: +# Command line args: tDb=fromDb, qDb=toDb +my ($tDb, $qDb); + +# Other: +my ($buildDir); +my ($tSeq, $tSizes, $qSeq, $qSizes, $QDb, $fileServer); +my ($liftOverChainDir, $liftOverChainFile, $liftOverChainPath, $dbExists); + +sub checkOptions { + # Make sure command line options are valid/supported. + my $ok = GetOptions(@HgStepManager::optionSpec, + 'buildDir=s', + 'target2Bit=s', + 'targetSizes=s', + 'query2Bit=s', + 'querySizes=s', + 'minimapPreset=s', + 'minimapCpu=i', + 'chainRam=s', + 'chainCpu=i', + 'regenerateMash', + @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); + if ($opt_minimapPreset) { + $minimapPreset = $opt_minimapPreset; + if ($minimapPreset !~ /^asm(5|10|20)$/) { + die "-minimapPreset must be one of asm5, asm10, asm20 (got '$minimapPreset')\n"; + } + } + # else: leave $minimapPreset undef -- estimateDivergence() will set it + # from a mash distance once $tSeq/$qSeq/$buildDir are known. +} + + +sub getClusterSeqs { + # Choose cluster and look for already-installed 2bit files on appropriate + # cluster-scratch storage. Exit with an error message if we can't find them. + my $paraHub = $opt_bigClusterHub ? $opt_bigClusterHub : + &HgAutomate::chooseClusterByBandwidth(); + my ($tSeqScratch, $qSeqScratch); + if ($opt_target2Bit) { + $tSeqScratch = $opt_target2Bit + } else { + my @okFilesystems = + &HgAutomate::chooseFilesystemsForCluster($paraHub, 'in'); + if ( -e "/scratch/data/$tDb/$tDb.2bit" ) { + $tSeqScratch = "/scratch/data/$tDb/$tDb.2bit"; + } else { + foreach my $fs (@okFilesystems) { + &HgAutomate::verbose(1, "checking $fs/$tDb/$tDb.2bit\n"); + if (&HgAutomate::machineHasFile($paraHub, "$fs/$tDb/$tDb.2bit")) { + $tSeqScratch = "$fs/$tDb/$tDb.2bit"; + last; + } + } + } + if (! defined $tSeqScratch) { + die "align: can't find $tDb/$tDb.2bit in " . + join("/, ", @okFilesystems) . "/ -- please distribute.\n"; + } + } + + if ($opt_query2Bit) { + $qSeqScratch = $opt_query2Bit; + } else { + my @okFilesystems = + &HgAutomate::chooseFilesystemsForCluster($paraHub, 'in'); + if ( -e "/scratch/data/$qDb/$qDb.2bit" ) { + $qSeqScratch = "/scratch/data/$qDb/$qDb.2bit"; + } else { + foreach my $fs (@okFilesystems) { + if (&HgAutomate::machineHasFile($paraHub, "$fs/$qDb/$qDb.2bit")) { + $qSeqScratch = "$fs/$qDb/$qDb.2bit"; + last; + } + } + } + if (! defined $qSeqScratch) { + die "align: can't find $qDb/$qDb.2bit in " . + join("/, ", @okFilesystems) . "/ -- please distribute.\n"; + } + } + &HgAutomate::verbose(1, "Using $paraHub, $tSeqScratch and $qSeqScratch\n"); + return ($paraHub, $tSeqScratch, $qSeqScratch); +} # getClusterSeqs + + + +######################################################################### +# * step: align [bigClusterHub] +sub doAlign { + my $runDir = "$buildDir/run.mm2"; + &HgAutomate::mustMkdir($runDir); + + my $pafDir = "$runDir/paf"; + &HgAutomate::checkCleanSlate('align', 'chain', $pafDir, 'run.time'); + + my ($paraHub, $tSeqScratch, $qSeqScratch) = &getClusterSeqs(); + + # Unlike blat -fastMap, minimap2 does not need the query pre-split into + # small chunks or an .ooc repeat mask -- it indexes and aligns whole + # sequences directly. We only partition on the target side, one job + # per target sequence (same trick doSameSpeciesLiftOver.pl uses: pick a + # chunkSize at least as large as the biggest target sequence so nothing + # actually gets split), and align each target job against the *whole* + # query genome in one shot. + + # script for a single job: convert target 2bit spec(s) to fasta, align + # against the pre-built whole-query fasta, convert PAF -> PSL. + # NOTE: partitionSequence.pl may bundle several small target sequences + # into one $runDir/tParts/partNNN.lst file instead of handing us a bare + # 2bit spec directly (it does this whenever a genome has lots of tiny + # scaffolds -- see its own usage message) -- path1 in that case is a + # *relative* path like "tParts/part010.lst", so it has to be resolved + # against $runDir before we cd elsewhere, and expanded into one + # twoBitToFa call per line rather than treated as a single 2bit spec. + my $fh = &HgAutomate::mustOpen(">$runDir/job.sh"); + print $fh <<_EOF_ +#!/bin/bash +set -beEu -o pipefail + +targetSpec=\$1 +outPsl=\$2 + +if [[ "\$targetSpec" == *.lst ]]; then + targetSpec="$runDir/\$targetSpec" +fi + +unset TMPDIR +if [ -d "/data/tmp" ]; then + export TMPDIR="/data/tmp" +elif [ -d "/scratch/tmp" ]; then + export TMPDIR="/scratch/tmp" +else + tmpSz=`df --output=avail -k /tmp | tail -1` + shmSz=`df --output=avail -k /dev/shm | tail -1` + if [ "\$shmSz" -gt "\$tmpSz" ]; then + mkdir -p /dev/shm/tmp + chmod 777 /dev/shm/tmp + export TMPDIR="/dev/shm/tmp" + else + export TMPDIR="/tmp" + fi +fi + +# Use local disk for output, and move the final result to \$outPsl +# when done, to minimize I/O. +tmpDir=`mktemp -d -p \$TMPDIR doMm2.XXXXXX` +pushd \$tmpDir > /dev/null + +: > target.fa +if [[ "\$targetSpec" == *.lst ]]; then + while read -r spec; do + twoBitToFa "\$spec" stdout + done < "\$targetSpec" > target.fa +else + twoBitToFa \$targetSpec target.fa +fi + +minimap2 -cx $minimapPreset --secondary=no -t $minimapCpu \\ + target.fa $runDir/query.fa > tmpOut.paf + +pafToPsl -tSizes=$tSizes -qSizes=$qSizes tmpOut.paf tmpOut.psl + +mv tmpOut.psl \$outPsl + +popd > /dev/null +rm -rf \$tmpDir +_EOF_ + ; + close($fh); + &HgAutomate::run("chmod a+x $runDir/job.sh"); + + &HgAutomate::makeGsub($runDir, + 'job.sh $(path1) {check out line ' . + $pafDir . '/$(file1).psl}'); + + my $paraRun = &HgAutomate::paraRun($ramG, $minimapCpu); + my $whatItDoes = "It performs a cluster run of minimap2 -cx $minimapPreset."; + my $bossScript = newBash HgRemoteScript("$runDir/doAlign.bash", $paraHub, + $runDir, $whatItDoes); + + # Don't allow target sequences to be split -- we align (and chain) whole + # target sequences against the whole query, we never lift target coords + # back up. Use the max target seq size as the chunkSize for + # partitionSequence.pl on the target. + my $tpSize = `awk '{print \$2;}' $tSizes | sort -nr | head -1`; + chomp $tpSize; + my $minTpSize = 10000000; + $tpSize = $minTpSize if ($tpSize < $minTpSize); + + my $gensub2 = &HgAutomate::gensub2(); + $bossScript->add(<<_EOF_ +# Convert the whole query 2bit to fasta once; every target job aligns +# against this same file. +twoBitToFa $qSeqScratch query.fa + +# Compute partition (coordinate ranges) for cluster job. This does +# not need to be run on the build fileserver because it does not actually +# split any sequences -- it merely computes ranges based on the chrom.sizes. +rm -rf tParts +$Bin/partitionSequence.pl $tpSize 0 $tSeqScratch \\ + $tSizes 2000 \\ + -lstDir=tParts > t.lst + +mkdir $pafDir + +$gensub2 t.lst single gsub jobList + +$paraRun +_EOF_ + ); + $bossScript->execute(); +} # doAlign + + +######################################################################### +# * step: chain [smallClusterHub] + +sub doChain { + my $runDir = "$buildDir/run.chain"; + &HgAutomate::mustMkdir($runDir); + + my $pafDir = "$buildDir/run.mm2/paf"; + my $mm2DoneFile = "$buildDir/run.mm2/run.time"; + &HgAutomate::checkCleanSlate('chain', 'net', 'chainRaw'); + &HgAutomate::checkExistsUnlessDebug('align', 'chain', + $pafDir, $mm2DoneFile); + + my ($paraHub, $tSeqScratch, $qSeqScratch) = &getClusterSeqs(); + + # One PSL per target sequence already (align step did not split further), + # so the chain job list is simply the contents of $pafDir/*.psl. + &HgAutomate::run("ls $pafDir/*.psl | xargs -n 1 basename > $runDir/pslParts.lst"); + + # script for a single job: chain one target-sequence's PSL. + my $fh = &HgAutomate::mustOpen(">$runDir/job.sh"); + print $fh <<_EOF_ +#!/bin/bash +set -beEu -o pipefail + +inPsl=\$1 +outChain=\$2 + +unset TMPDIR +if [ -d "/data/tmp" ]; then + export TMPDIR="/data/tmp" +elif [ -d "/scratch/tmp" ]; then + export TMPDIR="/scratch/tmp" +else + tmpSz=`df --output=avail -k /tmp | tail -1` + shmSz=`df --output=avail -k /dev/shm | tail -1` + if [ "\$shmSz" -gt "\$tmpSz" ]; then + mkdir -p /dev/shm/tmp + chmod 777 /dev/shm/tmp + export TMPDIR="/dev/shm/tmp" + else + export TMPDIR="/tmp" + fi +fi + +tmpOut=`mktemp -p \$TMPDIR doMm2.chain.XXXXXX` + +axtChain -verbose=0 -linearGap=medium -psl $pafDir/\$inPsl \\ + $tSeqScratch $qSeqScratch stdout \\ +| chainBridge -linearGap=medium stdin $tSeqScratch $qSeqScratch \\ + \$tmpOut +mv \$tmpOut \$outChain +chmod 664 \$outChain +_EOF_ + ; + close($fh); + &HgAutomate::run("chmod a+x $runDir/job.sh"); + + &HgAutomate::makeGsub($runDir, + 'job.sh $(path1) ' . + '{check out line+ chainRaw/$(path1).chain}'); + my $whatItDoes = "It does a cluster run to chain the minimap2 alignments."; + my $bossScript = newBash HgRemoteScript("$runDir/doChain.bash", $paraHub, + $runDir, $whatItDoes); + my $paraRun = &HgAutomate::paraRun($chainRam, $chainCpu); + my $gensub2 = &HgAutomate::gensub2(); + $bossScript->add(<<_EOF_ +mkdir chainRaw + +$gensub2 pslParts.lst single gsub jobList +$paraRun +_EOF_ + ); + $bossScript->execute(); +} # doChain + + +######################################################################### +# * step: net [workhorse] +sub doNet { + my $runDir = "$buildDir/run.chain"; + my @outs = ("$runDir/$tDb.$qDb.all.chain.gz", + "$runDir/$tDb.$qDb.noClass.net.gz"); + &HgAutomate::checkCleanSlate('net', 'load', @outs); + &HgAutomate::checkExistsUnlessDebug('chain', 'net', "$runDir/chainRaw/"); + + my $whatItDoes = +"It nets the chained minimap2 alignments and runs netChainSubset to produce +liftOver chains."; + my $mach = &HgAutomate::chooseWorkhorse(); + my $bossScript = newBash HgRemoteScript("$runDir/doNet.bash", $mach, + $runDir, $whatItDoes); + my $chromBased = (`wc -l < $tSizes` <= $HgAutomate::splitThreshold); + my $lump = $chromBased ? "" : "-lump=100"; + $bossScript->add(<<_EOF_ +unset TMPDIR +if [ -d "/data/tmp" ]; then + export TMPDIR="/data/tmp" +elif [ -d "/scratch/tmp" ]; then + export TMPDIR="/scratch/tmp" +else + tmpSz=`df --output=avail -k /tmp | tail -1` + shmSz=`df --output=avail -k /dev/shm | tail -1` + if [ "\$shmSz" -gt "\$tmpSz" ]; then + mkdir -p /dev/shm/tmp + chmod 777 /dev/shm/tmp + export TMPDIR="/dev/shm/tmp" + else + export TMPDIR="/tmp" + fi +fi +# Use local scratch disk... this can be quite I/O intensive: +tmpDir=`mktemp -d -p \$TMPDIR doMm2.net.XXXXXX` + +# Merge up the hierarchy and assign unique chain IDs: +chainMergeSort chainRaw/*.chain \\ +| chainSplit $lump \$tmpDir/chainSplit stdin +endsInLf \$tmpDir/chainSplit/*.chain + +mkdir \$tmpDir/netSplit \$tmpDir/overSplit +for f in \$tmpDir/chainSplit/*.chain; do + split=\$(basename "\$f" .chain) + chainNet \$f \\ + $tSizes $qSizes \\ + \$tmpDir/netSplit/\$split.net /dev/null + netChainSubset \$tmpDir/netSplit/\$split.net \$f stdout \\ + | chainStitchId stdin \$tmpDir/overSplit/\$split.chain +done +endsInLf \$tmpDir/netSplit/*.net +endsInLf \$tmpDir/overSplit/*.chain + +cat \$tmpDir/chainSplit/*.chain | gzip -c > $tDb.$qDb.all.chain.gz +cat \$tmpDir/netSplit/*.net | gzip -c > $tDb.$qDb.noClass.net.gz + +cat \$tmpDir/overSplit/*.chain | gzip -c > $buildDir/$liftOverChainFile +# make quickLift chain: +chainSwap $buildDir/$liftOverChainFile stdout \\ + | chainToBigChain stdin $buildDir/$tDb.$qDb.quick.chain.txt \\ + $buildDir/$tDb.$qDb.quick.link.txt + +rm -rf \$tmpDir/ +_EOF_ + ); + $bossScript->execute(); +} # doNet + + +######################################################################### +# * step: load [dbHost] +sub doLoad { + my $runDir = "$buildDir"; + &HgAutomate::checkExistsUnlessDebug('net', 'load', + "$buildDir/$liftOverChainFile"); + + my $whatItDoes = +"It makes links from $HgAutomate::gbdb/ and goldenPath/ (download area) to the liftOver +chains file, and calls hgAddLiftOverChain to register the $HgAutomate::gbdb location."; + my $bossScript = newBash HgRemoteScript("$runDir/doLoad.bash", $dbHost, + $runDir, $whatItDoes); + + if ($dbExists) { + $bossScript->add(<<_EOF_ +# Link to standardized location of liftOver files: +mkdir -p $liftOverChainDir +rm -f $liftOverChainPath +ln -s $buildDir/$liftOverChainFile $liftOverChainPath +tmpFile=`mktemp -t -p /tmp tmpMd5.XXXXXX` +grep -v $liftOverChainFile $liftOverChainDir/md5sum.txt > \$tmpFile || true +md5sum $buildDir/$liftOverChainFile | sed -e "s#$buildDir/##;" >> \$tmpFile +sort \$tmpFile > $liftOverChainDir/md5sum.txt +rm -f \$tmpFile + +# Link from download area: +mkdir -p $HgAutomate::goldenPath/$tDb/liftOver +rm -f $HgAutomate::goldenPath/$tDb/liftOver/$liftOverChainFile +ln -s $liftOverChainPath $HgAutomate::goldenPath/$tDb/liftOver/ + +# Link from genome browser fileserver: +mkdir -p $HgAutomate::gbdb/$tDb/liftOver +rm -f $HgAutomate::gbdb/$tDb/liftOver/$liftOverChainFile +ln -s $liftOverChainPath $HgAutomate::gbdb/$tDb/liftOver/ + +# Add an entry to liftOverChain table in central database (specified in +# ~/.hg.conf) so that hgLiftOver will know that this is available: +hgAddLiftOverChain $tDb $qDb +_EOF_ + ); + } else { + $bossScript->add(<<_EOF_ +hgLoadChain -test -noBin -tIndex $tDb chain$QDb $buildDir/$liftOverChainFile +wget --no-check-certificate -O bigChain.as 'https://raw.githubusercontent.com/ucscGenomeBrowser/kent/refs/heads/master/src/hg/lib/bigChain.as' +wget --no-check-certificate -O bigLink.as 'https://raw.githubusercontent.com/ucscGenomeBrowser/kent/refs/heads/master/src/hg/lib/bigLink.as' +sed 's/.000000//' chain.tab | awk 'BEGIN {OFS="\\t"} {print \$2, \$4, \$5, \$11, 1000, \$8, \$3, \$6, \$7, \$9, \$10, \$1}' > chain${QDb}.tab +bedToBigBed -type=bed6+6 -as=bigChain.as -tab chain${QDb}.tab $tSizes chain${QDb}.bb +awk 'BEGIN {OFS="\\t"} {print \$1, \$2, \$3, \$5, \$4}' link.tab | sort -k1,1 -k2,2n > chain${QDb}Link.tab +bedToBigBed -type=bed4+1 -as=bigLink.as -tab chain${QDb}Link.tab $tSizes chain${QDb}Link.bb + +bedToBigBed -type=bed6+6 -as=bigChain.as -tab $tDb.$qDb.quick.chain.txt $qSizes $tDb.$qDb.quick.bb +bedToBigBed -type=bed4+1 -as=bigLink.as -tab $tDb.$qDb.quick.link.txt $qSizes $tDb.$qDb.quickLink.bb + +totalBases=`ave -col=2 $tSizes | grep "^total" | awk '{printf "%d", \$2}'` +basesCovered=`bigBedInfo chain${QDb}Link.bb | grep "basesCovered" | cut -d' ' -f2 | tr -d ','` +percentCovered=`echo \$basesCovered \$totalBases | awk '{printf "%.3f", 100.0*\$1/\$2}'` +printf "%d bases of %d (%s%%) in intersection\\n" "\$basesCovered" "\$totalBases" "\$percentCovered" > fb.$tDb.chain.${QDb}Link.txt + +qBases=`ave -col=2 $qSizes | grep "^total" | awk '{printf "%d", \$2}'` +qCovered=`bigBedInfo $tDb.$qDb.quickLink.bb | grep "basesCovered" | cut -d' ' -f2 | tr -d ','` +qPerCent=`echo \$qCovered \$qBases | awk '{printf "%.3f", 100.0*\$1/\$2}'` +printf "%d bases of %d (%s%%) in intersection\\n" "\$qCovered" "\$qBases" "\$qPerCent" > fb.$tDb.quick${QDb}Link.txt +rm -f link.tab chain.tab bigChain.as bigLink.as chain${QDb}.tab chain${QDb}Link.tab $tDb.$qDb.quick.chain.txt $tDb.$qDb.quick.link.txt + +_EOF_ + ); + } + $bossScript->execute(); +} # doLoad + + +######################################################################### +# * step: cleanup [fileServer] +sub doCleanup { + my $runDir = "$buildDir"; + my $whatItDoes = "It cleans up or compresses intermediate files."; + $fileServer = &HgAutomate::chooseFileServer($runDir); + my $bossScript = newBash HgRemoteScript("$runDir/doCleanup.bash", $fileServer, + $runDir, $whatItDoes); + $bossScript->add(<<_EOF_ +rm -f run.mm2/query.fa +rm -rf run.mm2/paf/ +rm -rf run.chain/chainRaw/ +rm -f mashSketch.a.msh mashSketch.b.msh +_EOF_ + ); + $bossScript->execute(); +} # doCleanup + + +sub getSeqAndSizes { + if ($opt_target2Bit) { + $tSeq = $opt_target2Bit + } else { + # Test assumptions about 2bit and chrom.sizes files. + $tSeq = "/scratch/data/$tDb/$tDb.2bit"; + if (! -e $tSeq) { + # allow it to exist here too: + my $fs = "$HgAutomate::clusterData"; + &HgAutomate::verbose(1, "checking $fs/$tDb/$tDb.2bit\n"); + if (-e "$fs/$tDb/$tDb.2bit") { + $tSeq = "$fs/$tDb/$tDb.2bit"; + } + } + } + + if ($opt_targetSizes) { + $tSizes = $opt_targetSizes; + } else { + $tSizes = "$HgAutomate::clusterData/$tDb/chrom.sizes"; + } + + if ($opt_query2Bit) { + $qSeq = $opt_query2Bit; + } else { + $qSeq = "/scratch/data/$qDb/$qDb.2bit"; + if (! -e $qSeq) { + # allow it to exist here too: + my $fs = "$HgAutomate::clusterData"; + &HgAutomate::verbose(1, "checking $fs/$qDb/$qDb.2bit\n"); + if (-e "$fs/$qDb/$qDb.2bit") { + $qSeq = "$fs/$qDb/$qDb.2bit"; + } + } + } + + if ($opt_querySizes) { + $qSizes = $opt_querySizes; + } else { + $qSizes = "$HgAutomate::clusterData/$qDb/chrom.sizes"; + } + + my $problem = 0; + foreach my $file ($tSeq, $tSizes, $qSeq, $qSizes) { + if (! -e $file) { + warn "Error: cannot find required file \"$file\"\n"; + $problem = 1; + } + } + if ($problem && !$opt_debug) { + warn "Run $base -help for a description of expected files.\n"; + exit 1; + } +} + + +sub estimateDivergence { + # Pick -minimapPreset automatically via AssemblyDivergence.pm's mash + # distance between tDb and qDb, unless the user already gave + # -minimapPreset explicitly. Cheap (mash sketch+dist on a whole genome + # is seconds to low minutes), so we just run it locally rather than + # turning it into its own cluster step. + if ($minimapPreset) { + &HgAutomate::verbose(1, + "Using explicit -minimapPreset=$minimapPreset (mash not run).\n"); + return; + } + my $mashFile = "$buildDir/mashDistance.txt"; + if (-e $mashFile && ! $opt_regenerateMash) { + # Reuse a previous estimate so -continue steps stay consistent with + # whatever the align step (if already run) actually used. + open(my $fh, "<", $mashFile) || die "Can't read $mashFile: $!\n"; + my $line = <$fh>; + close($fh); + if ($line =~ /^mashDistance=(\S+)\s+minimapPreset=(\S+)/) { + $minimapPreset = $2; + &HgAutomate::verbose(1, "Reusing cached mash distance $1 -> " . + "-minimapPreset=$minimapPreset from $mashFile\n"); + return; + } + } + if ($opt_debug) { + $minimapPreset = 'asm5'; + &HgAutomate::verbose(1, "-debug: skipping mash, using -minimapPreset=$minimapPreset\n"); + return; + } + &HgAutomate::verbose(1, + "Estimating $tDb/$qDb divergence with mash to pick -minimapPreset ...\n"); + # $dbHost matters here: hgcentraltest (used to recognize a plain UCSC + # db's own .2bit) is only reachable from hgwdev, never from whatever + # host might eventually run this -- pass our own $dbHost through + # rather than let AssemblyDivergence.pm assume one. + my $dist = eval { &AssemblyDivergence::mashDistance($tSeq, $qSeq, $buildDir, $opt_regenerateMash, $dbHost); }; + if ($@) { + $minimapPreset = 'asm5'; + warn "estimateDivergence: $@" . + "falling back to -minimapPreset=$minimapPreset. Pass " . + "-minimapPreset explicitly to pick it yourself and silence this.\n"; + return; + } + my ($pipeline, $preset, $warning) = &AssemblyDivergence::choosePipeline($dist); + warn "$warning\n" if ($warning); + if ($pipeline ne 'minimap2') { + die "estimateDivergence: mash distance $dist between $tDb and $qDb " . + "looks too diverged for doMiniMap2.pl -- use doBlastzChainNet.pl " . + "instead, or re-run with -minimapPreset to force minimap2 anyway.\n"; + } + $minimapPreset = $preset; + open(my $fh, ">", $mashFile) || die "Can't write $mashFile: $!\n"; + print $fh "mashDistance=$dist minimapPreset=$minimapPreset\n"; + close($fh); + &HgAutomate::verbose(1, + "mash distance $dist between $tDb and $qDb -> -minimapPreset=$minimapPreset\n"); +} # estimateDivergence + +######################################################################### +# main + +# Prevent "Suspended (tty input)" hanging: +&HgAutomate::closeStdin(); + +&checkOptions(); + +&usage(1) if (scalar(@ARGV) != 2); +($tDb, $qDb) = @ARGV; + +# may be working on a 2bit file that does not have a database browser +$dbExists = 0; +$dbExists = 1 if (&HgAutomate::databaseExists($dbHost, $tDb)); + +&getSeqAndSizes(); +$QDb = ucfirst($qDb); +$liftOverChainDir = "$HgAutomate::clusterData/$tDb/$HgAutomate::trackBuild/liftOver"; +$liftOverChainFile = "${tDb}To${QDb}.over.chain.gz"; +$liftOverChainPath = "$liftOverChainDir/$liftOverChainFile"; +$chainRam = $opt_chainRam ? $opt_chainRam : $chainRam; +$chainCpu = $opt_chainCpu ? $opt_chainCpu : $chainCpu; +$minimapCpu = $opt_minimapCpu ? $opt_minimapCpu : $minimapCpu; +$ramG = $opt_ram ? $opt_ram : $ramG; + +my $date = `date +%Y-%m-%d`; +chomp $date; +$buildDir = $opt_buildDir ? $opt_buildDir : + "$HgAutomate::clusterData/$tDb/$HgAutomate::trackBuild/mm2.$qDb.$date"; + +if (! -d $buildDir) { + if ($stepper->stepPrecedes('align', $stepper->getStartStep())) { + die "$buildDir does not exist; try running again with -buildDir.\n"; + } + &HgAutomate::mustMkdir($buildDir); +} + +&estimateDivergence(); + +$stepper->execute(); + +my $stopStep = $stepper->getStopStep(); +my $upThrough = ($stopStep eq 'cleanup') ? "" : + " (through the '$stopStep' step)"; + +&HgAutomate::verbose(1, + "\n *** All done!$upThrough\n"); +&HgAutomate::verbose(1, + " *** Steps were performed in $buildDir\n"); +if ($stepper->stepPrecedes('net', $stopStep)) { + &HgAutomate::verbose(1, + " *** Test installation ($HgAutomate::gbdb, goldenPath, hgLiftover " . + "operation) on $dbHost.\n"); +} +&HgAutomate::verbose(1, "\n");