80484198a43b2152023cadd0b6abd26785ab69ea hiram Wed Jul 29 14:53:53 2026 -0700 now archiving versions of NCBI GenBank gene track when updated refs #37958 diff --git src/hg/utils/automation/doNcbiGene.pl src/hg/utils/automation/doNcbiGene.pl new file mode 100755 index 00000000000..80f7ced1911 --- /dev/null +++ src/hg/utils/automation/doNcbiGene.pl @@ -0,0 +1,384 @@ +#!/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 +required arguments: + asmId - assembly identifier, e.g. GCF_000001405.32_GRCh38.p6 + +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). + -chromSizes path Path to the assembly's chrom.sizes file. + -namesFile path Path to the hub's \$asmId.names.tab (built by the + gatewayPage step), e.g. + \$buildDir/../../html/\$asmId.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 + 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 ($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"; + 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 $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"); + + 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). +sub writeArchiveHub { + my ($archiveDir, $priorVersion) = @_; + my $haveIx = ( -s "$archiveDir/$asmId.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. + 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 +_HUB_ + if ($haveIx) { + print $fh "searchIndex name\nsearchTrix $asmId.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 $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"; + 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") ) { + &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)"; + } + + $bossScript->add(<<_EOF_ +export asmId=$asmId +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 +_EOF_ + ); + $bossScript->execute(); +} # doCleanup + +######################################################################### +# main + +&HgAutomate::closeStdin(); + +&checkOptions(); +&usage(1) if (scalar(@ARGV) != 1); + +$secondsStart = `date "+%s"`; +chomp $secondsStart; + +($asmId) = @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");