5c6501db6f2650e1f23fb8d196b4c9c0b3fd8d10 hiram Wed Jul 29 15:11:42 2026 -0700 wire in the new procedure for doNcbiGene refs #37958 diff --git src/hg/utils/automation/doNcbiGene.pl src/hg/utils/automation/doNcbiGene.pl index 80f7ced1911..4484f38bdf7 100755 --- src/hg/utils/automation/doNcbiGene.pl +++ src/hg/utils/automation/doNcbiGene.pl @@ -1,384 +1,392 @@ #!/usr/bin/env perl # DO NOT EDIT the /cluster/bin/scripts copy of this file -- # edit ~/kent/src/hg/utils/automation/doNcbiGene.pl instead. use Getopt::Long; use warnings; use strict; use File::stat; use FindBin qw($Bin); use lib "$Bin"; use HgAutomate; use HgRemoteScript; use HgStepManager; use AsmHub; # Option variable names, both common and peculiar to this script: use vars @HgAutomate::commonOptionVars; use vars @HgStepManager::optionVars; use vars qw/ $opt_buildDir $opt_assemblySource $opt_chromSizes $opt_namesFile $opt_liftFile /; # Specify the steps supported with -continue / -stop: my $stepper = new HgStepManager( [ { name => 'ncbiGene', func => \&doNcbiGene }, { name => 'cleanup', func => \&doCleanup }, ] ); # Option defaults: my $dbHost = 'hgwdev'; my $workhorse = 'hgwdev'; my $defaultWorkhorse = 'hgwdev'; my $defaultFileServer = 'hgwdev'; my $fileServer = 'hgwdev'; my $base = $0; $base =~ s/^(.*\/)?//; sub usage { my ($status, $detailed) = @_; print STDERR " -usage: $base [options] asmId +usage: $base [options] asmId db required arguments: - asmId - assembly identifier, e.g. GCF_000001405.32_GRCh38.p6 + asmId - assembly identifier at NCBI, e.g. + GCF_000001405.32_GRCh38.p6 -- used to locate the + source gff3/lift/remove.dups.list files. + db - name used for the track's own output files, e.g. + bigDataUrl bbi/db.ncbiGene.bb -- may equal asmId + for a plain GenArk build, or be a custom -dbName. options: "; print STDERR $stepper->getOptionHelp(); print STDERR <<_EOF_ -buildDir dir Use dir instead of default (current directory). This *is* the runDir -- typically the hub's trackData/ncbiGene -- no further nesting is added. - -assemblySource dir Directory holding \${asmId}_genomic.gff.gz (and - optionally \${asmId}.remove.dups.list). + -assemblySource dir Directory holding \${asmId}_genomic.gff.gz. -chromSizes path Path to the assembly's chrom.sizes file. - -namesFile path Path to the hub's \$asmId.names.tab (built by the + -namesFile path Path to the hub's \$db.names.tab (built by the gatewayPage step), e.g. - \$buildDir/../../html/\$asmId.names.tab -- used + \$buildDir/../../html/\$db.names.tab -- used for the archived-version description page. -liftFile path Optional lift file translating NCBI names to UCSC names, e.g. \$buildDir/../../sequence/\$asmId.ncbiToUcsc.lift _EOF_ ; print STDERR &HgAutomate::getCommonOptionHelp('dbHost' => $dbHost, 'workhorse' => $defaultWorkhorse, 'fileServer' => $defaultFileServer); print STDERR " Automates construction of the 'ncbiGene' track from an assembly's own NCBI GFF3 gene predictions, for assembly hub (GenArk) builds. Steps: - ncbiGene: if a previous \${asmId}.ncbiGene.bb exists and the source gff + ncbiGene: if a previous \$db.ncbiGene.bb exists and the source gff is newer, archive the previous version under archive/<date>/ (keyed by the previous build's own gff-derived mtime), then translate the current gff3 into a bigGenePred track and swap it into place atomically. cleanup: compress intermediate files "; print "\n"; exit $status; } # Globals: -my ($asmId); +my ($asmId, $db); my ($buildDir, $assemblySource, $chromSizes, $namesFile, $liftFile); my ($secondsStart, $secondsEnd); sub checkOptions { my $ok = GetOptions(@HgStepManager::optionSpec, 'buildDir=s', 'assemblySource=s', 'chromSizes=s', 'namesFile=s', 'liftFile=s', @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); $workhorse = $opt_workhorse if ($opt_workhorse); $fileServer = $opt_fileServer if ($opt_fileServer); } # same mtime-compare doAssemblyHub.pl has used all along sub needsUpdate($$) { my ($source, $result) = @_; if (-s $result) { return (stat($source)->mtime > stat($result)->mtime) ? 1 : 0; } return 1; } ######################################################################### # archive whatever the previous build left behind, keyed by the date # that was stamped onto it (the existing 'touch -r $gffFile' convention # means the old .bb's mtime already *is* the previous source gff's date, # so it doubles as a version string for free -- no separate release # metadata needed, unlike ncbiRefSeq's NCBI-supplied $verString). sub archivePriorVersion { - my $priorBb = "$buildDir/$asmId.ncbiGene.bb"; + my $priorBb = "$buildDir/$db.ncbiGene.bb"; return if (! -s "$priorBb"); my $priorMtime = stat($priorBb)->mtime; my ($mday,$mon,$year) = (localtime($priorMtime))[3,4,5]; my $priorVersion = sprintf("%04d-%02d-%02d", $year+1900, $mon+1, $mday); my $archiveDir = "$buildDir/archive/$priorVersion"; if ( -d "$archiveDir" ) { &HgAutomate::verbose(1, "# ncbiGene: archive/$priorVersion already exists, not re-archiving\n"); return; } &HgAutomate::mustMkdir($archiveDir); foreach my $ext (qw( bb gtf.gz ix ixx stats.txt )) { - my $f = "$buildDir/$asmId.ncbiGene.$ext"; - rename($f, "$archiveDir/$asmId.ncbiGene.$ext") if (-e "$f"); + my $f = "$buildDir/$db.ncbiGene.$ext"; + rename($f, "$archiveDir/$db.ncbiGene.$ext") if (-e "$f"); } - my $geneAttrs = "$buildDir/$asmId.geneAttrs.ncbi.txt.gz"; - rename($geneAttrs, "$archiveDir/$asmId.geneAttrs.ncbi.txt.gz") if (-e "$geneAttrs"); - my $fb = "$buildDir/fb.$asmId.ncbiGene.txt"; - rename($fb, "$archiveDir/fb.$asmId.ncbiGene.txt") if (-e "$fb"); + my $geneAttrs = "$buildDir/$db.geneAttrs.ncbi.txt.gz"; + rename($geneAttrs, "$archiveDir/$db.geneAttrs.ncbi.txt.gz") if (-e "$geneAttrs"); + my $fb = "$buildDir/fb.$db.ncbiGene.txt"; + rename($fb, "$archiveDir/fb.$db.ncbiGene.txt") if (-e "$fb"); writeArchiveHub($archiveDir, $priorVersion); &HgAutomate::verbose(1, "# ncbiGene: archived previous version to archive/$priorVersion/\n"); } # archivePriorVersion ######################################################################### # a single hub.txt (useOneFile on) that puts just the archived ncbiGene # track up for viewing. No twoBitPath/organism/defaultPos are needed in -# the genome stanza: $asmId is already a published GenArk genome, so -# 'genome $asmId' with no twoBitPath just attaches this one extra track -# to that already-known assembly (see trackHubGenomeReadRa() in -# hg/lib/trackHub.c -- twoBitPath is only required when *introducing* a -# new genome). +# the genome stanza: the true NCBI accession is already a published +# GenArk genome, so 'genome $accession' with no twoBitPath just attaches +# this one extra track to that already-known assembly (see +# trackHubGenomeReadRa() in hg/lib/trackHub.c -- twoBitPath is only +# required when *introducing* a new genome). sub writeArchiveHub { my ($archiveDir, $priorVersion) = @_; - my $haveIx = ( -s "$archiveDir/$asmId.ncbiGene.ix" ); + my $haveIx = ( -s "$archiveDir/$db.ncbiGene.ix" ); writeArchiveHtml($archiveDir, $priorVersion); # 'genome' wants the bare accession (GCF_937001465.1), not the full # asmId with its _assemblyName suffix (GCF_937001465.1_mOrcOrc1.1) -- # that suffix is what NCBI adds for the gff/download file names, but - # it is not the name the genome is registered under. + # it is not the name the genome is registered under. This must come + # from the true NCBI $asmId, not the output-naming $db, since $db may + # be an arbitrary -dbName. my @parts = split('_', $asmId); my $accession = "$parts[0]_$parts[1]"; my $hubTxt = "$archiveDir/hub.txt"; open(my $fh, ">", $hubTxt) or die "can not write $hubTxt: $!"; print $fh <<_HUB_; hub ${asmId}_ncbiGene_${priorVersion} shortLabel $asmId ncbiGene archive $priorVersion longLabel Archived NCBI GFF3 gene track for $asmId, superseded $priorVersion when NCBI updated the source annotation useOneFile on email genome-www\@soe.ucsc.edu genome $accession track ncbiGene shortLabel NCBI GenBank ($priorVersion) longLabel Gene models submitted to GenBank, ENA, DDBJ -- archived version $priorVersion visibility pack color 0,80,150 altColor 150,80,0 colorByStrand 0,80,150 150,80,0 type bigGenePred -bigDataUrl $asmId.ncbiGene.bb -html $asmId.ncbiGene +bigDataUrl $db.ncbiGene.bb +html $db.ncbiGene _HUB_ if ($haveIx) { - print $fh "searchIndex name\nsearchTrix $asmId.ncbiGene.ix\n"; + print $fh "searchIndex name\nsearchTrix $db.ncbiGene.ix\n"; } close($fh); } # writeArchiveHub ######################################################################### # the description page for the archived hub above, built by handing the # same paths off to AsmHub::ncbiGeneDescription() that asmHubNcbiGene.pl # uses for the live track -- one source of truth for what this page says, # just with an extra banner noting it is a superseded version. sub writeArchiveHtml { my ($archiveDir, $priorVersion) = @_; - my $bbPath = "$archiveDir/$asmId.ncbiGene.bb"; - my $statsPath = "$archiveDir/$asmId.ncbiGene.stats.txt"; + my $bbPath = "$archiveDir/$db.ncbiGene.bb"; + my $statsPath = "$archiveDir/$db.ncbiGene.stats.txt"; my $archiveNote = "This is an archived copy of the ncbiGene track as it stood on $priorVersion, before NCBI's source GFF3 annotation was updated to a newer version."; my $html = AsmHub::ncbiGeneDescription($bbPath, $statsPath, $chromSizes, $namesFile, $asmId, $archiveNote); - my $htmlPath = "$archiveDir/$asmId.ncbiGene.html"; + my $htmlPath = "$archiveDir/$db.ncbiGene.html"; open(my $fh, ">", $htmlPath) or die "can not write $htmlPath: $!"; print $fh $html; close($fh); } # writeArchiveHtml ######################################################################### # * step: ncbiGene [workhorse] sub doNcbiGene { my $gffFile = "$assemblySource/${asmId}_genomic.gff.gz"; if ( ! -s "${gffFile}" ) { &HgAutomate::verbose(1, "# step ncbiGene: no gff file found at:\n# $gffFile\n"); return; } - if ( ! needsUpdate($gffFile, "$buildDir/$asmId.ncbiGene.bb") ) { + if ( ! needsUpdate($gffFile, "$buildDir/$db.ncbiGene.bb") ) { &HgAutomate::verbose(1, "# ncbiGene step previously completed\n"); return; } &HgAutomate::mustMkdir($buildDir); archivePriorVersion(); my $whatItDoes = "translate NCBI GFF3 gene definitions into a track"; my $bossScript = newBash HgRemoteScript("$buildDir/doNcbiGene.bash", $workhorse, $buildDir, $whatItDoes); my $dupList = ""; - if ( -s "${assemblySource}/${asmId}.remove.dups.list" ) { - $dupList = " | (grep -v -f \"${assemblySource}/${asmId}.remove.dups.list\" || true)"; + # this list is curated by hand into the build's own download/ area, + # not something NCBI supplies in $assemblySource + if ( -s "${buildDir}/../../download/${asmId}.remove.dups.list" ) { + $dupList = " | (grep -v -f \"${buildDir}/../../download/${asmId}.remove.dups.list\" || true)"; } $bossScript->add(<<_EOF_ -export asmId=$asmId +export asmId=$db export gffFile=$gffFile export chromSizes=$chromSizes function cleanUp() { rm -f \$asmId.ncbiGene.genePred.gz \$asmId.ncbiGene.genePred rm -f \$asmId.geneAttrs.ncbi.txt } if [ \$gffFile -nt \$asmId.ncbiGene.bb ]; then ln -sf \$gffFile ./ (gff3ToGenePred -warnAndContinue -useName \\ -attrsOut=\$asmId.geneAttrs.ncbi.txt \$gffFile stdout \\ 2>> \$asmId.ncbiGene.log.txt || true) | genePredFilter \\ -chromSizes=\$chromSizes stdin stdout \\ $dupList | gzip -c > \$asmId.ncbiGene.genePred.gz genePredCheck \$asmId.ncbiGene.genePred.gz zcat \$asmId.ncbiGene.genePred.gz > ncbiGene.\$asmId.gp genePredToGtf -utr file ncbiGene.\$asmId.gp stdout | gzip -c > \$asmId.ncbiGene.gtf.gz rm -f ncbiGene.\$asmId.gp export howMany=`genePredCheck \$asmId.ncbiGene.genePred.gz 2>&1 | grep "^checked" | awk '{print \$2}'` if [ "\${howMany}" -eq 0 ]; then printf "# ncbiGene: no gene definitions found in \$gffFile\\n" cleanUp exit 0 fi export ncbiGenePred="\$asmId.ncbiGene.genePred.gz" _EOF_ ); if ( -s "$liftFile" ) { $bossScript->add(<<_EOF_ liftUp -extGenePred -type=.gp stdout \\ $liftFile warn \\ \$asmId.ncbiGene.genePred.gz | gzip -c \\ > \$asmId.ncbiGene.ucsc.genePred.gz ncbiGenePred="\$asmId.ncbiGene.ucsc.genePred.gz" _EOF_ ); } $bossScript->add(<<_EOF_ ~/kent/src/hg/utils/automation/gpToIx.pl \$ncbiGenePred \\ > \$asmId.gpToIx.txt ~/kent/src/hg/utils/automation/gffAttrsToIx.py \$asmId.geneAttrs.ncbi.txt \\ \$ncbiGenePred > \$asmId.attrsToIx.txt sort -u \$asmId.gpToIx.txt \$asmId.attrsToIx.txt > \$asmId.ncbiGene.ix.txt if [ -s \$asmId.ncbiGene.ix.txt ]; then ixIxx \$asmId.ncbiGene.ix.txt \$asmId.ncbiGene.ix.new \$asmId.ncbiGene.ixx.new fi rm -f \$asmId.ncbiGene.ix.txt \$asmId.gpToIx.txt \$asmId.attrsToIx.txt genePredToBigGenePred \$ncbiGenePred stdout \\ | sort -k1,1 -k2,2n > \$asmId.ncbiGene.bed (bedToBigBed -type=bed12+8 -tab -as=\$HOME/kent/src/hg/lib/bigGenePred.as \\ -extraIndex=name \$asmId.ncbiGene.bed \\ \$chromSizes \$asmId.ncbiGene.bb.new || true) if [ ! -s "\$asmId.ncbiGene.bb.new" ]; then printf "# ncbiGene: failing bedToBigBed\\n" 1>&2 exit 255 fi touch -r\$gffFile \$asmId.ncbiGene.bb.new bigBedInfo \$asmId.ncbiGene.bb.new | egrep "^itemCount:|^basesCovered:" \\ | sed -e 's/,//g' > \$asmId.ncbiGene.stats.txt LC_NUMERIC=en_US /usr/bin/printf "# ncbiGene %s %'d %s %'d\\n" `cat \$asmId.ncbiGene.stats.txt` | xargs echo # basesCovered comes straight out of the bigBedInfo call above -- # no need to separately rebuild exons via bedToExons/bedSingleCover.pl export totalBases=`ave -col=2 \$chromSizes | grep total | awk '{printf "%d", \$NF}'` export basesCovered=`grep basesCovered \$asmId.ncbiGene.stats.txt | awk '{printf "%s", \$NF}'` export percentCovered=`echo \$basesCovered \$totalBases | awk '{printf "%.3f", 100.0*\$1/\$2}'` printf "%d bases of %d (%s%%) in intersection\\n" "\$basesCovered" "\$totalBases" "\$percentCovered" > fb.\$asmId.ncbiGene.txt # atomic swap -- a client mid-read of the .bb/.ix/.ixx over http never # sees a half-written file, unlike the old build-in-place behavior mv -f \$asmId.ncbiGene.bb.new \$asmId.ncbiGene.bb if [ -s \$asmId.ncbiGene.ix.new ]; then mv -f \$asmId.ncbiGene.ix.new \$asmId.ncbiGene.ix mv -f \$asmId.ncbiGene.ixx.new \$asmId.ncbiGene.ixx fi else printf "# ncbiGene step previously completed\\n" 1>&2 fi _EOF_ ); $bossScript->execute(); } # doNcbiGene ######################################################################### # * step: cleanup [fileServer] sub doCleanup { my $whatItDoes = "compress intermediate files"; my $bossScript = new HgRemoteScript("$buildDir/doCleanup.csh", $fileServer, $buildDir, $whatItDoes); $bossScript->add(<<_EOF_ -gzip -f $asmId.geneAttrs.ncbi.txt $asmId.ncbiGene.log.txt +gzip -f $db.geneAttrs.ncbi.txt $db.ncbiGene.log.txt _EOF_ ); $bossScript->execute(); } # doCleanup ######################################################################### # main &HgAutomate::closeStdin(); &checkOptions(); -&usage(1) if (scalar(@ARGV) != 1); +&usage(1) if (scalar(@ARGV) != 2); $secondsStart = `date "+%s"`; chomp $secondsStart; -($asmId) = @ARGV; +($asmId, $db) = @ARGV; $assemblySource = $opt_assemblySource or die "ERROR: -assemblySource is required\n"; $chromSizes = $opt_chromSizes or die "ERROR: -chromSizes is required\n"; $namesFile = $opt_namesFile or die "ERROR: -namesFile is required\n"; $liftFile = $opt_liftFile ? $opt_liftFile : ""; $buildDir = $opt_buildDir ? $opt_buildDir : `pwd`; chomp $buildDir; $stepper->execute(); $secondsEnd = `date "+%s"`; chomp $secondsEnd; my $elapsedSeconds = $secondsEnd - $secondsStart; my $elapsedMinutes = int($elapsedSeconds/60); $elapsedSeconds -= $elapsedMinutes * 60; my $stopStep = $stepper->getStopStep(); my $upThrough = ($stopStep eq 'cleanup') ? "" : " (through the '$stopStep' step)"; &HgAutomate::verbose(1, "\n *** All done !$upThrough Elapsed time: ${elapsedMinutes}m${elapsedSeconds}s\n"); &HgAutomate::verbose(1, " *** Steps were performed in $buildDir\n"); &HgAutomate::verbose(1, "\n");