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");