70b908a5c6475d6c2770a3305f9319d88d0cb62e hiram Mon Sep 28 22:31:27 2026 -0700 readjusting output results and feature bit file locations refs #34360 diff --git src/hg/utils/automation/doMiniMap2.pl src/hg/utils/automation/doMiniMap2.pl index b552c8624bf..48408288837 100755 --- src/hg/utils/automation/doMiniMap2.pl +++ src/hg/utils/automation/doMiniMap2.pl @@ -1,781 +1,860 @@ #!/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_minimapSecondaryN + $opt_minimapSecondaryRatio $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; +my $ramG = '24g'; # minimap2 index + 8 threads on a whole genome needs headroom # 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 +# Secondary alignments are kept (not '--secondary=no') so axtChain/chainNet +# see the same kind of overlapping candidate alignments a lastz/blastz run +# would give them, and can use genome-scale synteny context to arbitrate +# ambiguous/duplicated regions -- rather than trusting minimap2's own local, +# per-query-region primary/secondary call alone. Defaults match minimap2's +# own -N/-p defaults; only override via -minimapSecondaryN/-Ratio if you +# have a specific reason to. +my $minimapSecondaryN = 5; # -N INT to minimap2 +my $minimapSecondaryRatio = 0.8; # -p FLOAT to minimap2 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 + $HgAutomate::clusterData/\$fromDb/$HgAutomate::trackBuild/miniMap2\$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 + -ram Ng Cluster ram size for minimap2, default: -ram=$ramG + -minimapSecondaryN N minimap2 -N: max secondary alignments retained per + primary hit, default: $minimapSecondaryN. Secondary + alignments are kept (not -secondary=no) so + axtChain/chainNet get the same kind of overlapping + candidate alignments a lastz run would give them, + and can resolve ambiguous/duplicated regions using + genome-scale synteny context instead of relying on + minimap2's own local primary/secondary call alone. + -minimapSecondaryRatio F minimap2 -p: minimum secondary-to-primary score + ratio to keep a secondary at all, default: + $minimapSecondaryRatio -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, + 'minimapCpu' => $minimapCpu, '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', + 'minimapSecondaryN=i', + 'minimapSecondaryRatio=f', '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 \\ +minimap2 -cx $minimapPreset -N $minimapSecondaryN -p $minimapSecondaryRatio -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_ +# Record how the preset was chosen (mash distance, if it ran -- see +# estimateDivergence()) and exactly which minimap2 and what arguments were +# used, for reproducibility -- the per-job command line is otherwise +# buried in job.sh on whatever cluster node happened to run it. +if [ -s $buildDir/mashDistance.txt ]; then + cat $buildDir/mashDistance.txt > $buildDir/version.txt +else + printf "minimapPreset=$minimapPreset (explicit -minimapPreset, mash not run)\\n" > $buildDir/version.txt +fi +printf "minimap2 " >> $buildDir/version.txt +minimap2 --version >> $buildDir/version.txt +printf "minimap2 -cx $minimapPreset -N $minimapSecondaryN -p $minimapSecondaryRatio -t $minimapCpu target.fa query.fa > out.paf\\n" >> $buildDir/version.txt + # 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"; + my $runDir = "$buildDir/axtChain/run"; &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 $runDir = "$buildDir/axtChain"; my @outs = ("$runDir/$tDb.$qDb.all.chain.gz", - "$runDir/$tDb.$qDb.noClass.net.gz"); + "$runDir/noClass.net"); &HgAutomate::checkCleanSlate('net', 'load', @outs); - &HgAutomate::checkExistsUnlessDebug('chain', 'net', "$runDir/chainRaw/"); + &HgAutomate::checkExistsUnlessDebug('chain', 'net', "$runDir/run/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 \\ +chainMergeSort run/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 +# noClass.net is scaffolding for netChainSubset above, not a kept product +# (matches doBlastzChainNet.pl's own treatment -- see cleanup step): +cat \$tmpDir/netSplit/*.net > noClass.net -cat \$tmpDir/overSplit/*.chain | gzip -c > $buildDir/$liftOverChainFile +cat \$tmpDir/overSplit/*.chain | gzip -c > $runDir/$liftOverChainFile # make quickLift chain: -chainSwap $buildDir/$liftOverChainFile stdout \\ - | chainToBigChain stdin $buildDir/$tDb.$qDb.quick.chain.txt \\ - $buildDir/$tDb.$qDb.quick.link.txt +chainSwap $runDir/$liftOverChainFile stdout \\ + | chainToBigChain stdin $runDir/$tDb.$qDb.quick.chain.txt \\ + $runDir/$tDb.$qDb.quick.link.txt rm -rf \$tmpDir/ _EOF_ ); $bossScript->execute(); } # doNet ######################################################################### # * step: load [dbHost] sub doLoad { - my $runDir = "$buildDir"; + my $runDir = "$buildDir/axtChain"; &HgAutomate::checkExistsUnlessDebug('net', 'load', - "$buildDir/$liftOverChainFile"); + "$runDir/$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."; +chains file, calls hgAddLiftOverChain to register the $HgAutomate::gbdb location, and +builds the quickLift chain/link bigBed pair (for any target/query pair, regardless of +whether \$tDb is a real UCSC database)."; my $bossScript = newBash HgRemoteScript("$runDir/doLoad.bash", $dbHost, $runDir, $whatItDoes); + $bossScript->add(<<_EOF_ +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' +_EOF_ + ); + if ($dbExists) { $bossScript->add(<<_EOF_ # Link to standardized location of liftOver files: mkdir -p $liftOverChainDir rm -f $liftOverChainPath -ln -s $buildDir/$liftOverChainFile $liftOverChainPath +ln -s $runDir/$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 +md5sum $runDir/$liftOverChainFile | sed -e "s#$runDir/##;" >> \$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' +# Full-chain bigChain pair, from the complete alignment (all.chain.gz) -- +# the standard genome-browser Chain/Net-style track pair: +hgLoadChain -test -noBin -tIndex $tDb chain$QDb $runDir/$tDb.$qDb.all.chain.gz 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 +printf "%d bases of %d (%s%%) in intersection\\n" "\$basesCovered" "\$totalBases" "\$percentCovered" > $buildDir/fb.$tDb.chain${QDb}Link.txt + +rm -f link.tab chain.tab chain${QDb}.tab chain${QDb}Link.tab + +# LiftOver-subset bigChain pair, from the netChainSubset-extracted chain +# (over.chain.gz) -- the strict, single-best-path-per-region subset: +hgLoadChain -test -noBin -tIndex $tDb chainLiftOver$QDb $runDir/$liftOverChainFile +sed 's/.000000//' chain.tab | awk 'BEGIN {OFS="\\t"} {print \$2, \$4, \$5, \$11, 1000, \$8, \$3, \$6, \$7, \$9, \$10, \$1}' > chainLiftOver${QDb}.tab +bedToBigBed -type=bed6+6 -as=bigChain.as -tab chainLiftOver${QDb}.tab $tSizes chainLiftOver${QDb}.bb +awk 'BEGIN {OFS="\\t"} {print \$1, \$2, \$3, \$5, \$4}' link.tab | sort -k1,1 -k2,2n > chainLiftOver${QDb}Link.tab +bedToBigBed -type=bed4+1 -as=bigLink.as -tab chainLiftOver${QDb}Link.tab $tSizes chainLiftOver${QDb}Link.bb + +liftOverBasesCovered=`bigBedInfo chainLiftOver${QDb}Link.bb | grep "basesCovered" | cut -d' ' -f2 | tr -d ','` +liftOverPercentCovered=`echo \$liftOverBasesCovered \$totalBases | awk '{printf "%.3f", 100.0*\$1/\$2}'` +printf "%d bases of %d (%s%%) in intersection\\n" "\$liftOverBasesCovered" "\$totalBases" "\$liftOverPercentCovered" > $buildDir/fb.$tDb.chainLiftOver${QDb}.txt + +rm -f link.tab chain.tab chainLiftOver${QDb}.tab chainLiftOver${QDb}Link.tab +_EOF_ + ); + } + + # quickLift bigBed pair: built for any target/query pair regardless of + # $dbExists -- it only needs the quick.chain.txt/quick.link.txt doNet + # already produced, plus $qSizes, no dependency on \$tDb being a real + # database. + $bossScript->add(<<_EOF_ +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 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 +printf "%d bases of %d (%s%%) in intersection\\n" "\$qCovered" "\$qBases" "\$qPerCent" > $buildDir/fb.$tDb.quick${QDb}Link.txt +rm -f bigChain.as bigLink.as $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 -rf axtChain/run/chainRaw/ +rm -f axtChain/noClass.net +# mashSketch.{a,b}.msh only ever exist here if estimateDivergence() fell +# back to a one-off sketch (tSeq/qSeq didn't resolve to a real GenArk +# accession or UCSC db -- see AssemblyDivergence.pm's sketch()); for a +# normal assembly the sketch lives in its own permanent, shared cache +# under /hive/data/genomes/, never here, so this is a no-op in that case: 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"); my $dist = eval { &AssemblyDivergence::mashDistance($tSeq, $qSeq, $buildDir, $opt_regenerateMash); }; 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; +$minimapSecondaryN = $opt_minimapSecondaryN ? $opt_minimapSecondaryN : $minimapSecondaryN; +$minimapSecondaryRatio = $opt_minimapSecondaryRatio ? $opt_minimapSecondaryRatio : $minimapSecondaryRatio; $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"; + "$HgAutomate::clusterData/$tDb/$HgAutomate::trackBuild/miniMap2${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");