25ea393c310876076793bbad0e093f83161e9d78 hiram Tue Aug 18 14:39:58 2026 -0700 improve efficiency of the trackData.pl an now using faSize to get sizes out of the 2bit files no redmine diff --git src/hg/makeDb/doc/asmHubs/trackData.pl src/hg/makeDb/doc/asmHubs/trackData.pl index c8aba832782..50f58aad6ab 100755 --- src/hg/makeDb/doc/asmHubs/trackData.pl +++ src/hg/makeDb/doc/asmHubs/trackData.pl @@ -1,75 +1,67 @@ #!/usr/bin/env perl use strict; use warnings; use File::Basename; use FindBin qw($Bin); use lib "$Bin"; use commonHtml; use File::stat; my $argc = scalar(@ARGV); -if ($argc < 3) { - printf STDERR "usage: trackData.pl Name asmHubName [two column name list] > trackData.html\n"; - printf STDERR "e.g.: trackData.pl Mammals mammals mammals.asmId.commonName.tsv > trackData.html\n"; +if ($argc != 5) { + printf STDERR "usage: trackData.pl Name asmHubName orderList prodOutFile testOutFile\n"; + printf STDERR "e.g.: trackData.pl Mammals mammals mammals.asmId.commonName.tsv trackData.html testTrackData.html\n"; printf STDERR "the name list is found in \$HOME/kent/src/hg/makeDb/doc/asmHubs/\n"; - printf STDERR "\nthe two columns are 1: asmId (accessionId_assemblyName)\n"; + printf STDERR "\nthe two columns in the name list are 1: asmId (accessionId_assemblyName)\n"; printf STDERR "column 2: common name for species, columns separated by tab\n"; + printf STDERR "\nWrites both the production and the '-test' variant of the track\n"; + printf STDERR "statistics page in a single pass over the assembly list -- each\n"; + printf STDERR "assembly's per-track stats (bigWigInfo/bigBedInfo/faSize etc.) are\n"; + printf STDERR "measured once and reused for both pages, instead of running the\n"; + printf STDERR "whole scan twice.\n"; exit 255; } my $home = $ENV{'HOME'}; my $toolsDir = "$home/kent/src/hg/makeDb/doc/asmHubs"; my $sciNameOverrideFile = "$toolsDir/sciNameOverride.txt"; my %sciNameOverride; # key is accession, value is corrected scientific name my %taxIdOverride; # key is accession, value is corrected taxId # keys for both of those can also be the asmId if ( -s "${sciNameOverrideFile}" ) { open (my $sn, "<", "${sciNameOverrideFile}") or die "can not read ${sciNameOverrideFile}"; while (my $line = <$sn>) { next if ($line =~ m/^#/); next if (length($line) < 2); chomp $line; my ($accO, $asmIdO, $sciNameO, $taxIdO) = split('\t', $line); $sciNameOverride{$accO} = $sciNameO; $sciNameOverride{$asmIdO} = $sciNameO; $taxIdOverride{$accO} = $taxIdO; $taxIdOverride{$asmIdO} = $taxIdO; } close ($sn); } - -my $testOutput = 0; -my $spliceOut = -1; - -if ($argc > 2) { - for (my $i = 0; $i < $argc; ++$i) { - if ($ARGV[$i] =~ /-test/) { - $testOutput = 1; - $spliceOut = $i; - } - } -} -if ($spliceOut != -1) { - splice @ARGV, $spliceOut, 1; -} my $Name = shift; my $asmHubName = shift; my $inputList = shift; +my $prodOutFile = shift; +my $testOutFile = shift; my $orderList = $inputList; if ( ! -s "$orderList" ) { $orderList = $toolsDir/$inputList; } my @orderList; # asmId of the assemblies in order from the orderList file my %commonName; # key is asmId, value is a common name, perhaps more appropriate # than found in assembly_report file # assembly_report my $vgpIndex = 0; $vgpIndex = 1 if ($Name =~ m/vgp/i); my $hprcIndex = 0; $hprcIndex = 1 if ($Name =~ m/hprc/i); my $brcIndex = 0; $brcIndex = 1 if ($Name =~ m/brc/i); @@ -78,30 +70,48 @@ my $asmCount = 0; # count of assemblies completed and in the table my $overallNucleotides = 0; my $overallSeqCount = 0; my $overallGapSize = 0; my $overallGapCount = 0; ############################################################################## # from Perl Cookbook Recipe 2.17, print out large numbers with comma delimiters: ############################################################################## sub commify($) { my $text = reverse $_[0]; $text =~ s/(\d\d\d)(?=\d)(?!\d*\.)/$1,/g; return scalar reverse $text } +############################################################################## +# capture(&): run a block of code that 'print's/'printf's, and return +# everything it printed as a string, instead of letting it go to STDOUT. +# This lets startHtml()/startTable()/endTable()/endHtml() keep their +# original print-based bodies unchanged, while the caller decides which +# of the two output files (or both) a given fragment belongs in. +############################################################################## +sub capture(&) { + my ($code) = @_; + my $buf = ''; + open(my $fh, '>', \$buf) or die "capture: $!"; + my $old = select($fh); + $code->(); + select($old); + close($fh); + return $buf; +} + # ($itemCount, $percentCover) = bigWigMeasure($trackFile, $genomeSize); sub bigWigMeasure($$) { my ($file, $genomeSize) = @_; my $bigWigInfo = `bigWigInfo "$file" | egrep "basesCovered:|mean:" | awk '{print \$NF}' | xargs echo | sed -e 's/,//g;'`; chomp $bigWigInfo; my ($bases, $mean) = split('\s+', $bigWigInfo); my $itemCount = sprintf ("%.2f", $mean); my $percentCover = sprintf("%.2f %%", 100.0 * $bases / $genomeSize); return ($itemCount, $percentCover); } # $percentCover = pcFbFile($trackFb); sub pcFbFile($) { my ($trackFb) = @_; my ($itemBases, undef, undef, $noGapSize, undef) = split('\s+', `cat $trackFb`, 5); @@ -139,31 +149,31 @@ my $bigBedInfo = `bigBedInfo "$file" | egrep "itemCount:|basesCovered:" | awk '{print \$NF}' | xargs echo | sed -e 's/,//g;'`; chomp $bigBedInfo; my ($items, $bases) = split('\s', $bigBedInfo); $itemCount = commify($items); $percentCover = sprintf("%.2f %%", 100.0 * $bases / $genomeSize); if ( -s "${trackFb}" ) { $percentCover = pcFbFile($trackFb); } # printf STDERR "# bigBedInfo %s %s %s\n", $itemCount, $percentCover, $file; } } return ($itemCount, $percentCover); } # sub oneTrackData($$$$$$) ############################################################################## -### start the HTML output +### start the HTML output -- identical for the prod and -test pages ############################################################################## sub startHtml() { my $timeStamp = `date "+%F"`; chomp $timeStamp; my $subSetMessage = "subset of $asmHubName only"; if ($asmHubName eq "vertebrate") { $subSetMessage = "subset of other ${asmHubName}s only"; } if ($vgpIndex) { my $vgpSubset = "(set of primary assemblies)"; if ($orderList =~ m/vgp.alternate/) { $vgpSubset = "(set of alternate/haplotype assemblies)"; @@ -249,42 +259,76 @@ } my $indexUrl = "index"; my $asmStats = "asmStats"; print <<"END" <h3>See also: <a href='$indexUrl.html'>hub access</a>, <a href='$asmStats.html'>assembly statistics</a></h3><br> <h3>Data resource links</h3> NOTE: <em>Click on the column headers to sort the table by that column</em><br> The <em>link to genome browser</em> will attach only that single assembly to the genome browser.<br> The numbers are: item count (percent coverage)<br> Except for the gc5Base column which is: overall GC % average (percent coverage) END } # sub startHtml() -# order of columns in the table +############################################################################## +# buildTrackList($testOutput, $asmHubName): the order of columns in the +# table, for either the production ($testOutput = 0) or -test ($testOutput +# = 1) variant. This used to be a single global array that tableContents() +# mutated in place with these same splices, gated on a global $testOutput -- +# calling this twice (once per variant) instead reproduces both variants +# exactly, without needing two separate runs of the script. +############################################################################## +sub buildTrackList($$) { + my ($testOutput, $asmHubName) = @_; # eliminated the ncbiGene track -my @trackList = qw(ncbiRefSeq xenoRefGene augustus ensGene gc5Base allGaps assembly rmsk simpleRepeat windowMasker cpgIslandExtUnmasked); -### XXX beware, this trackList is going to be edited below to add or -### remove elements depending upon the situation + my @list = qw(ncbiRefSeq xenoRefGene augustus ensGene gc5Base allGaps assembly rmsk simpleRepeat windowMasker cpgIslandExtUnmasked); + if ($testOutput) { # add extra columns during 'test' output +# 0 1 2 3 4 5 6 +# 7 8 9 10 11 12 13 +# 14 +# my @trackList = qw(ncbiRefSeq xenoRefGene augustus ensGene gc5Base gap allGaps assembly rmsk simpleRepeat windowMasker gapOverlap tandemDups cpgIslandExtUnmasked cpgIslandExt); +# 0 1 2 3 4 5 6 +# 7 8 9 10 +# my @trackList = qw(ncbiRefSeq xenoRefGene augustus ensGene gc5Base allGaps assembly rmsk simpleRepeat windowMasker cpgIslandExtUnmasked); + splice @list, 11, 0, "cpgIslandExt"; + splice @list, 10, 0, "tandemDups"; + splice @list, 10, 0, "gapOverlap"; + splice @list, 5, 0, "gap"; + } + if ("viral" eq $asmHubName) { + splice @list, 3, 1; + splice @list, 2, 1; + splice @list, 1, 1; + } + if ($testOutput || ("viral" eq $asmHubName)) { # add extra columns during 'test' output + splice @list, 1, 0, "ncbiGene"; + } + if ("viral" eq $asmHubName) { + splice @list, 0, 1; + } + return @list; +} ############################################################################## ### start the table output ############################################################################## -sub startTable() { +sub startTable($) { + my ($testOutput) = @_; -# coordinate the order of these column headings with the @trackList listed above +# coordinate the order of these column headings with buildTrackList() above print '<table class="sortable" border="1"> <thead style="position:sticky; top:0;"><tr><th>count</th> <th>common name<br>link to genome browser</th> '; print '<th class="sorttable_numeric">ncbiRefSeq</th> ' if ("viral" ne $asmHubName); print " <th class=\"sorttable_numeric\">ncbiGene</th>\n" if ($testOutput || ("viral" eq $asmHubName)); print ' <th class="sorttable_numeric">xenoRefGene</th> <th class="sorttable_numeric">augustus<br>genes</th> <th class="sorttable_numeric">Ensembl<br>genes</th> ' if ("viral" ne $asmHubName); @@ -304,69 +348,64 @@ <th class="sorttable_numeric">window<br>Masker</th> '; if ($testOutput) { print ' <th class="sorttable_numeric">gap<br>Overlap</th> <th class="sorttable_numeric">tandem<br>Dups</th> <th class="sorttable_numeric">cpg<br>unmasked</th> <th class="sorttable_numeric">cpg<br>island</th> '; } else { print " <th class=\"sorttable_numeric\">cpg<br>islands</th>\n"; } print "</tr></thead><tbody>\n"; -} # sub startTable() +} # sub startTable($) ############################################################################## ### end the table output ############################################################################## -sub endTable() { - -my $commaNuc = commify($overallNucleotides); -my $commaSeqCount = commify($overallSeqCount); -my $commaGapSize = commify($overallGapSize); -my $commaGapCount = commify($overallGapCount); +sub endTable($$$) { + my ($assemblyTotal, $asmCount, $columnCount) = @_; my $percentDone = 100.0 * $asmCount / $assemblyTotal; my $doneMsg = ""; if ($asmCount < $assemblyTotal) { $doneMsg = sprintf(" (%d build completed, %.2f %% finished)", $asmCount, $percentDone); } -my $columnCount = scalar(@trackList); my $colSpanFill = $columnCount - 1; if ($assemblyTotal > 1) { print <<"END" </tbody> <tfoot><tr><th>TOTALS:</th><td style='text-align: center;' colspan=$colSpanFill>total assembly count ${assemblyTotal}${doneMsg}</td> </tr></tfoot> </table> END } else { print <<"END" </tbody> </table> END } -} # sub endTable() +} # sub endTable($$$) ############################################################################## -### end the HTML output +### end the HTML output -- identical for the prod and -test pages ############################################################################## sub endHtml() { &commonHtml::otherHubLinks($vgpIndex, $asmHubName); &commonHtml::htmlFooter($vgpIndex, $asmHubName); } # sub endHtml() sub asmCounts($) { my ($chromSizes) = @_; my ($sequenceCount, $totalSize) = split('\s+', `/cluster/bin/x86_64/ave -col=2 $chromSizes | egrep "^count|^total" | awk '{printf "%d\\n", \$NF}' | xargs echo`); return ($sequenceCount, $totalSize); } sub maskStats($) { @@ -391,58 +430,169 @@ my $gapBed = "$buildDir/trackData/allGaps/$asmId.allGaps.bed.gz"; my $gapCount = 0; if ($asmId !~ m/^GC/) { $gapBed = "/hive/data/genomes/$asmId/$asmId.N.bed"; if ( -s "$gapBed" ) { $gapCount = `awk '{print \$3-\$2}' $gapBed | /cluster/bin/x86_64/ave stdin | grep '^count' | awk '{print \$2}'`; } } elsif ( -s "$gapBed" ) { $gapCount = `zcat $gapBed | awk '{print \$3-\$2}' | /cluster/bin/x86_64/ave stdin | grep '^count' | awk '{print \$2}'`; } chomp $gapCount; return ($gapCount); } ############################################################################## -### tableContents() +# computeTrackCell($asmId, $track, $buildDir, $totalSize) +# returns (itemCount, percentCover, customKey) for one track of one +# assembly. This is the expensive part (bigWigInfo/bigBedInfo/hgsql/etc.) +# and is testOutput-independent, so it only needs to run once per assembly +# per track no matter how many page variants reference that track. +# +# NOTE: the original script additionally retried a "still n/a" ensGene or +# ncbiRefSeq track as ebiGene/ncbiGene, but on production output only. That +# retry re-checked the *same* file path that had just been found missing +# (only $runDir and the diagnostic track-name argument differed, and +# oneTrackData() only consults $runDir for the unrelated 'gapOverlap' case) +# -- so it was a guaranteed no-op and is not reproduced here. ############################################################################## -sub tableContents() { - my $asmCounted = 0; - if ($testOutput) { # add extra columns during 'test' output -# 0 1 2 3 4 5 6 -# 7 8 9 10 11 12 13 -# 14 -# my @trackList = qw(ncbiRefSeq xenoRefGene augustus ensGene gc5Base gap allGaps assembly rmsk simpleRepeat windowMasker gapOverlap tandemDups cpgIslandExtUnmasked cpgIslandExt); -# 0 1 2 3 4 5 6 -# 7 8 9 10 -# my @trackList = qw(ncbiRefSeq xenoRefGene augustus ensGene gc5Base allGaps assembly rmsk simpleRepeat windowMasker cpgIslandExtUnmasked); - splice @trackList, 11, 0, "cpgIslandExt"; - splice @trackList, 10, 0, "tandemDups"; - splice @trackList, 10, 0, "gapOverlap"; - splice @trackList, 5, 0, "gap"; +sub computeTrackCell($$$$) { + my ($asmId, $track, $buildDir, $totalSize) = @_; + my $trackFile = "$buildDir/bbi/$asmId.$track"; + my $trackFb = "$buildDir/trackData/$track/fb.$asmId.$track.txt"; + # no ensGene file ? Then look for ebiGene file + if ($track eq "ensGene" && ! -s $trackFb) { + if ( -d "$buildDir/trackData/ebiGene" ) { + $trackFb = "$buildDir/trackData/ebiGene/fb.ebiGene.txt" if ( -d "$buildDir/trackData/ebiGene/fb.ebiGene.txt"); + $trackFile = "$buildDir/bbi/$asmId.ebiGene"; } - if ("viral" eq $asmHubName) { - splice @trackList, 3, 1; - splice @trackList, 2, 1; - splice @trackList, 1, 1; } - if ($testOutput || ("viral" eq $asmHubName)) { # add extra columns during 'test' output - splice @trackList, 1, 0, "ncbiGene"; + my $runDir = "$buildDir/trackData/$track"; + my ($itemCount, $percentCover); + my $customKey = ""; + if ($asmId !~ m/^GC/) { + $itemCount = "n/a"; + $percentCover = "n/a"; + if ($track eq "ncbiRefSeq") { + my $refSeqDir=`ls -d /hive/data/genomes/$asmId/bed/ncbiRefSeq.20* | tail -1`; + chomp $refSeqDir; + if ( -d "${refSeqDir}" ) { + my $trackFb = "${refSeqDir}/fb.ncbiRefSeq.$asmId.txt"; + if ( -s "${trackFb}" ) { + $itemCount = `hgsql -N -e 'select count(*) from $track;' $asmId 2> /dev/null`; + chomp $itemCount; + $percentCover = pcFbFile($trackFb); } - if ("viral" eq $asmHubName) { - splice @trackList, 0, 1; } + } elsif ($track eq "gc5Base") { + my $bwFile = "/gbdb/$asmId/bbi/gc5Base.bw"; + $bwFile = "/gbdb/$asmId/bbi/gc5BaseBw/gc5Base.bw" if (! -s "${bwFile}"); + ($itemCount, $percentCover) = bigWigMeasure($bwFile, $totalSize); + } elsif ($track eq "rmsk") { + my $rmskStats = "/hive/data/genomes/$asmId/bed/repeatMasker/$asmId.rmsk.stats"; + if (! -s "${rmskStats}") { + my $faOut = "/hive/data/genomes/$asmId/bed/repeatMasker/$asmId.sorted.fa.out.gz"; + if ( -s "$faOut") { + my $items = `zgrep -c ^ "$faOut"`; + chomp $items; + $itemCount = commify($items); + my $masked = `grep masked "/hive/data/genomes/$asmId/bed/repeatMasker/faSize.rmsk.txt" | awk '{print \$4}' | sed -e 's/%//;'`; + chomp $masked; + $percentCover = sprintf("%.2f %%", $masked); + open (RS, ">$rmskStats") or die "can now write to $rmskStats"; + printf RS "%s\t%s\n", $itemCount, $percentCover; + close (RS); + } else { + $itemCount = "n/a"; + $percentCover = "n/a"; + } + } else { + ($itemCount, $percentCover) = split('\s+', `cat $rmskStats`); + chomp $percentCover; + $customKey = sprintf("%.2f", $percentCover); + $percentCover = sprintf("%.2f %%", $percentCover); + } + } # elsif ($track eq "rmsk") + } else { # working on an assembly hub + if ( "$track" eq "gc5Base" ) { + $trackFile .= ".bw"; + } else { + $trackFile .= ".bb"; + } + if ( "$track" eq "rmsk") { + my $rmskStats = "$buildDir/trackData/repeatMasker/$asmId.rmsk.stats"; + if (! -s "${rmskStats}") { + my $faOut = "$buildDir/trackData/repeatMasker/$asmId.sorted.fa.out.gz"; + if ( -s "$faOut") { + my $items = `zgrep -c ^ "$faOut"`; + chomp $items; + $itemCount = commify($items); + my $masked = `grep masked "$buildDir/trackData/repeatMasker/faSize.rmsk.txt" | awk '{print \$4}' | sed -e 's/%//;'`; + chomp $masked; + $percentCover = sprintf("%.2f %%", $masked); + open (RS, ">$rmskStats") or die "can now write to $rmskStats"; + printf RS "%s\t%s\n", $itemCount, $percentCover; + close (RS); + } else { + $itemCount = "n/a"; + $percentCover = "n/a"; + } + } else { + ($itemCount, $percentCover) = split('\s+', `cat $rmskStats`); + chomp $percentCover; + $customKey = sprintf("%.2f", $percentCover); + $percentCover = sprintf("%.2f %%", $percentCover); + } + } else { # not the rmsk track + ($itemCount, $percentCover) = oneTrackData($asmId, $track, $trackFile, $totalSize, $trackFb, $runDir); + } # else not the rmsk track + } # else if ($asmId !~ m/^GC/) + if (($percentCover =~ m/%/) || ($percentCover !~ m#n/a#)) { + $customKey = $percentCover; + $customKey =~ s/[ %]+//; + } + return ($itemCount, $percentCover, $customKey); +} # sub computeTrackCell($$$$) + +# render one <td> cell from a (itemCount, percentCover, customKey) triple +sub renderCell($$$) { + my ($itemCount, $percentCover, $customKey) = @_; + if (length($customKey)) { + return sprintf(" <td style='text-align: right;' sorttable_customkey='%s'>%s<br>(%s)</td>\n", $customKey, $itemCount, $percentCover); + } elsif ($itemCount eq "n/a") { + return " <td style='text-align: right;'>n/a</td>\n"; + } else { + return sprintf(" <td style='text-align: right;'>%s<br>(%s)</td>\n", $itemCount, $percentCover); + } +} + +############################################################################## +### tableContentsBoth() +### walks @orderList exactly once, measuring each assembly's tracks exactly +### once, and returns the table body HTML for both the production and +### -test pages (plus each page's column count, for endTable()'s colspan). +############################################################################## +sub tableContentsBoth() { + my @prodTrackList = buildTrackList(0, $asmHubName); + my @testTrackList = buildTrackList(1, $asmHubName); + my %inUnion; + my @unionTracks = grep { !$inUnion{$_}++ } (@prodTrackList, @testTrackList); + + my $prodBody = ""; + my $testBody = ""; + my $asmCounted = 0; + foreach my $asmId (@orderList) { my $gcPrefix = "GCx"; my $asmAcc = "asmAcc"; my $asmName = "asmName"; my $accessionId = "GCx_098765432.1"; my $accessionDir = ""; my $configRa = "n/a"; my $tracksCounted = 0; my $buildDir = "/hive/data/genomes/asmHubs/refseqBuild/$accessionDir/$asmId"; my $asmReport="$buildDir/download/${asmId}_assembly_report.txt"; my $chromSizes = "${buildDir}/${asmId}.chrom.sizes"; my $twoBit = "${buildDir}/trackData/addMask/${asmId}.masked.2bit"; my $faSizeTxt = "${buildDir}/${asmId}.faSize.txt"; if ($asmId !~ m/^GC/) { $configRa = "/hive/data/genomes/$asmId/$asmId.config.ra"; @@ -464,47 +614,47 @@ ($gcPrefix, $asmAcc, $asmName) = split('_', $asmId, 3); $accessionId = sprintf("%s_%s", $gcPrefix, $asmAcc); $accessionDir = substr($asmId, 0 ,3); $accessionDir .= "/" . substr($asmId, 4 ,3); $accessionDir .= "/" . substr($asmId, 7 ,3); $accessionDir .= "/" . substr($asmId, 10 ,3); $buildDir = "/hive/data/genomes/asmHubs/refseqBuild/$accessionDir/$asmId"; if ($gcPrefix eq "GCA") { $buildDir = "/hive/data/genomes/asmHubs/genbankBuild/$accessionDir/$asmId"; } $asmReport="$buildDir/download/${asmId}_assembly_report.txt"; $chromSizes = "${buildDir}/${asmId}.chrom.sizes"; $twoBit = "${buildDir}/trackData/addMask/${asmId}.masked.2bit"; $faSizeTxt = "${buildDir}/${asmId}.faSize.txt"; } -# my $trackDb="$buildDir/${asmId}.trackDb.txt"; -# next if (! -s "$trackDb"); # assembly build not complete if (! -s "$asmReport") { printf STDERR "# no assembly report:\n# %s\n", $asmReport; next; } if (! -s "$twoBit") { printf STDERR "# no 2bit file:\n# %s\n", $twoBit; - printf "<tr><td style='text-align: right;'>%d</td>\n", ++$asmCount; - printf "<td style='text-align: center;'>%s</td>\n", $accessionId; - printf "<th colspan=15 style='text-align: center;'>missing masked 2bit file</th>\n"; - printf "</tr>\n"; + my $missingRow = sprintf("<tr><td style='text-align: right;'>%d</td>\n", ++$asmCount); + $missingRow .= sprintf("<td style='text-align: center;'>%s</td>\n", $accessionId); + $missingRow .= "<th colspan=15 style='text-align: center;'>missing masked 2bit file</th>\n"; + $missingRow .= "</tr>\n"; + $prodBody .= $missingRow; + $testBody .= $missingRow; next; } if ( ! -s "$faSizeTxt" ) { - printf STDERR "twoBitToFa $twoBit stdout | faSize stdin > $faSizeTxt\n"; - print `twoBitToFa $twoBit stdout | faSize stdin > $faSizeTxt`; + printf STDERR "faSize $twoBit > $faSizeTxt\n"; + print `faSize $twoBit > $faSizeTxt`; } my ($gapSize, $maskPerCent, $sizeNoGaps) = maskStats($faSizeTxt); $overallGapSize += $gapSize; my ($seqCount, $totalSize) = asmCounts($chromSizes); $overallSeqCount += $seqCount; $overallNucleotides += $totalSize; my $gapCount = gapStats($buildDir, $asmId); $overallGapCount += $gapCount; my $sciName = "notFound"; $sciName = $sciNameOverride{$accessionId} if (defined($sciNameOverride{$accessionId})); my $commonName = "notFound"; my $asmDate = "notFound"; my $itemsFound = 0; open (FH, "<$asmReport") or die "can not read $asmReport"; while (my $line = <FH>) { @@ -520,186 +670,77 @@ } } elsif ($line =~ m/Organism name:/) { if ($sciName =~ m/notFound/) { ++$itemsFound; $commonName = $line; $sciName = $line; $commonName =~ s/.*\(//; $commonName =~ s/\)//; $commonName = $commonName{$asmId} if (exists($commonName{$asmId})); $sciName =~ s/.*:\s+//; $sciName =~ s/\s+\(.*//; } } } close (FH); - my $hubUrl = "https://hgdownload.soe.ucsc.edu/hubs/$accessionDir/$accessionId"; + + # the browser/hub links are the only thing that differ between the + # two pages besides the columns themselves my $browserName = $commonName; - my $browserUrl = "https://genome.ucsc.edu/h/$accessionId"; + my $prodBrowserUrl = "https://genome.ucsc.edu/h/$accessionId"; + my $testBrowserUrl = "https://genome-test.gi.ucsc.edu/h/$accessionId"; if ($asmId !~ m/^GC/) { - $hubUrl = "https://hgdownload.soe.ucsc.edu/goldenPath/$asmId/bigZips"; - $browserUrl = "https://genome.ucsc.edu/cgi-bin/hgTracks?db=$asmId"; + $prodBrowserUrl = "https://genome.ucsc.edu/cgi-bin/hgTracks?db=$asmId"; + $testBrowserUrl = "https://genome-test.gi.ucsc.edu/cgi-bin/hgTracks?db=$asmId"; $browserName = "$commonName ($asmId)"; - if ($testOutput) { - $browserUrl = "https://genome-test.gi.ucsc.edu/cgi-bin/hgTracks?db=$asmId"; - $hubUrl = "https://hgdownload-test.gi.ucsc.edu/goldenPath/$asmId/bigZips"; - } - } elsif ($testOutput) { - $browserUrl = "https://genome-test.gi.ucsc.edu/h/$accessionId"; - } - printf "<tr><td style='text-align: right;'>%d</td>\n", ++$asmCount; - printf "<td style='text-align: center;'><a href='%s' target=_blank>%s<br>%s</a></td>\n", $browserUrl, $browserName, $accessionId; - foreach my $track (@trackList) { - my $trackFile = "$buildDir/bbi/$asmId.$track"; - my $trackFb = "$buildDir/trackData/$track/fb.$asmId.$track.txt"; - # no ensGene file ? Then look for ebiGene file - if ($track eq "ensGene" && ! -s $trackFb) { - if ( -d "$buildDir/trackData/ebiGene" ) { - $trackFb = "$buildDir/trackData/ebiGene/fb.ebiGene.txt" if ( -d "$buildDir/trackData/ebiGene/fb.ebiGene.txt"); - $trackFile = "$buildDir/bbi/$asmId.ebiGene"; - } - } - my $runDir = "$buildDir/trackData/$track"; - my ($itemCount, $percentCover); - my $customKey = ""; - if ($asmId !~ m/^GC/) { - $itemCount = "n/a"; - $percentCover = "n/a"; - if ($track eq "ncbiRefSeq") { - my $refSeqDir=`ls -d /hive/data/genomes/$asmId/bed/ncbiRefSeq.20* | tail -1`; - chomp $refSeqDir; - if ( -d "${refSeqDir}" ) { - my $trackFb = "${refSeqDir}/fb.ncbiRefSeq.$asmId.txt"; - if ( -s "${trackFb}" ) { - $itemCount = `hgsql -N -e 'select count(*) from $track;' $asmId 2> /dev/null`; - chomp $itemCount; - $percentCover = pcFbFile($trackFb); - } } - } elsif ($track eq "gc5Base") { - my $bwFile = "/gbdb/$asmId/bbi/gc5Base.bw"; - $bwFile = "/gbdb/$asmId/bbi/gc5BaseBw/gc5Base.bw" if (! -s "${bwFile}"); - ($itemCount, $percentCover) = bigWigMeasure($bwFile, $totalSize); - } elsif ($track eq "rmsk") { - my $rmskStats = "/hive/data/genomes/$asmId/bed/repeatMasker/$asmId.rmsk.stats"; - if (! -s "${rmskStats}") { - my $faOut = "/hive/data/genomes/$asmId/bed/repeatMasker/$asmId.sorted.fa.out.gz"; - if ( -s "$faOut") { - my $items = `zgrep -c ^ "$faOut"`; - chomp $items; - $itemCount = commify($items); - my $masked = `grep masked "/hive/data/genomes/$asmId/bed/repeatMasker/faSize.rmsk.txt" | awk '{print \$4}' | sed -e 's/%//;'`; - chomp $masked; - $percentCover = sprintf("%.2f %%", $masked); - open (RS, ">$rmskStats") or die "can now write to $rmskStats"; - printf RS "%s\t%s\n", $itemCount, $percentCover; - close (RS); - } else { - $itemCount = "n/a"; - $percentCover = "n/a"; - } - } else { - ($itemCount, $percentCover) = split('\s+', `cat $rmskStats`); - chomp $percentCover; - $customKey = sprintf("%.2f", $percentCover); - $percentCover = sprintf("%.2f %%", $percentCover); - } - } # elsif ($track eq "rmsk") -########### need to figure which tables can be measured here -# x else x -# $itemCount = `hgsql -N -e 'select count(*) from $track;' $asmId 2> /dev/null`; -# chomp $itemCount; -# if (length($itemCount) < 1) { -# $itemCount = "n/a"; -# $percentCover = "n/a"; -# } else { -# $percentCover = `featureBits $asmId $track 2>&1 | cut -d' ' -f5 | tr -d ')('`; -# chomp $percentCover; -# $customKey = $percentCover; -# $customKey =~ s/[ %]+//; -# } - } else { # working on an assembly hub - if ( "$track" eq "gc5Base" ) { - $trackFile .= ".bw"; - } else { - $trackFile .= ".bb"; - } - if ( "$track" eq "rmsk") { - my $rmskStats = "$buildDir/trackData/repeatMasker/$asmId.rmsk.stats"; - if (! -s "${rmskStats}") { - my $faOut = "$buildDir/trackData/repeatMasker/$asmId.sorted.fa.out.gz"; - if ( -s "$faOut") { - my $items = `zgrep -c ^ "$faOut"`; - chomp $items; - $itemCount = commify($items); - my $masked = `grep masked "$buildDir/trackData/repeatMasker/faSize.rmsk.txt" | awk '{print \$4}' | sed -e 's/%//;'`; - chomp $masked; - $percentCover = sprintf("%.2f %%", $masked); - open (RS, ">$rmskStats") or die "can now write to $rmskStats"; - printf RS "%s\t%s\n", $itemCount, $percentCover; - close (RS); - } else { - $itemCount = "n/a"; - $percentCover = "n/a"; - } - } else { - ($itemCount, $percentCover) = split('\s+', `cat $rmskStats`); - chomp $percentCover; - $customKey = sprintf("%.2f", $percentCover); - $percentCover = sprintf("%.2f %%", $percentCover); - } - } else { # not the rmsk track - ($itemCount, $percentCover) = oneTrackData($asmId, $track, $trackFile, $totalSize, $trackFb, $runDir); - if (0 == $testOutput) { # only on the production stats page - # if track ensGene does not exist, try the ebiGene track - if ($track eq "ensGene" && $itemCount eq "n/a") { - $runDir = "$buildDir/trackData/ebiGene"; - $trackFile = "$buildDir/bbi/$asmId.$track.bb"; - ($itemCount, $percentCover) = oneTrackData($asmId, "ebiGene", $trackFile, $totalSize, $trackFb, $runDir); - } elsif ($track eq "ncbiRefSeq" && $itemCount eq "n/a") { - # if track ncbiRefSeq does not exist, try the ncbiGene track - $runDir = "$buildDir/trackData/ncbiGene"; - $trackFile = "$buildDir/bbi/$asmId.$track.bb"; - ($itemCount, $percentCover) = oneTrackData($asmId, "ncbiGene", $trackFile, $totalSize, $trackFb, $runDir); - } - } - } # else not the rmsk track - } # else if ($asmId !~ m/^GC/) - if (($percentCover =~ m/%/) || ($percentCover !~ m#n/a#)) { - $customKey = $percentCover; - $customKey =~ s/[ %]+//; + + ++$asmCount; + my $prodRow = sprintf("<tr><td style='text-align: right;'>%d</td>\n", $asmCount); + $prodRow .= sprintf("<td style='text-align: center;'><a href='%s' target=_blank>%s<br>%s</a></td>\n", $prodBrowserUrl, $browserName, $accessionId); + my $testRow = sprintf("<tr><td style='text-align: right;'>%d</td>\n", $asmCount); + $testRow .= sprintf("<td style='text-align: center;'><a href='%s' target=_blank>%s<br>%s</a></td>\n", $testBrowserUrl, $browserName, $accessionId); + + # measure every track needed by either page exactly once + my %cell; # key is track name, value is [itemCount, percentCover, customKey] + foreach my $track (@unionTracks) { + $cell{$track} = [ computeTrackCell($asmId, $track, $buildDir, $totalSize) ]; } - if (length($customKey)) { - printf " <td style='text-align: right;' sorttable_customkey='%s'>%s<br>(%s)</td>\n", $customKey, $itemCount, $percentCover; - } else { - if ($itemCount eq "n/a") { - printf " <td style='text-align: right;'>n/a</td>\n"; - } else { - printf " <td style='text-align: right;'>%s<br>(%s)</td>\n", $itemCount, $percentCover; + + foreach my $track (@prodTrackList) { + my ($itemCount, $percentCover, $customKey) = @{$cell{$track}}; + $tracksCounted += 1 if ($itemCount ne "n/a"); + $prodRow .= renderCell($itemCount, $percentCover, $customKey); } + foreach my $track (@testTrackList) { + my ($itemCount, $percentCover, $customKey) = @{$cell{$track}}; + $testRow .= renderCell($itemCount, $percentCover, $customKey); } - $tracksCounted += 1 if ($itemCount ne "n/a"); - } # foreach my $track (@trackList) - printf "</tr>\n"; + $prodRow .= "</tr>\n"; + $testRow .= "</tr>\n"; + $prodBody .= $prodRow; + $testBody .= $testRow; + $asmCounted += 1; if ($asmId =~ m/^GC/) { printf STDERR "# %03d\t%02d tracks\t%s\n", $asmCounted, $tracksCounted, $asmId; } else { printf STDERR "# %03d\t%02d tracks\t%s_%s (%s)\n", $asmCounted, $tracksCounted, $accessionId, $asmName, $asmId; } } -} + return ($prodBody, $testBody, scalar(@prodTrackList), scalar(@testTrackList)); +} # sub tableContentsBoth() ############################################################################## ### main() ############################################################################## # if there is a 'promoted' list, it has been taken out of the 'orderList' # so will need to stuff it back in at the correct ordered location my %promotedList; # key is asmId, value is common name my $promotedList = dirname(${orderList}) . "/promoted.list"; my @promotedList; # contents are asmIds, in order by lc(common name) my $promotedIndex = -1; # to walk through @promotedList; if ( -s "${promotedList}" ) { open (FH, "<${promotedList}" ) or die "can not read ${promotedList}"; while (my $line = <FH>) { @@ -727,20 +768,30 @@ if (lc($checkInsertName) lt lc($commonName)) { push @orderList, $checkInsertAsmId; $commonName{$checkInsertAsmId} = $checkInsertName; ++$assemblyTotal; printf STDERR "# inserting '%s' before '%s' at # %03d\n", $checkInsertName, $commonName, $assemblyTotal; ++$promotedIndex; # only doing one at this time # TBD: will need to improve this for more inserts } } push @orderList, $asmId; $commonName{$asmId} = $commonName; ++$assemblyTotal; } close (FH); -startHtml(); -startTable(); -tableContents(); -endTable(); -endHtml(); +my $header = capture { startHtml() }; +my ($prodBody, $testBody, $prodCols, $testCols) = tableContentsBoth(); +my $prodHead = capture { startTable(0) }; +my $testHead = capture { startTable(1) }; +my $prodFoot = capture { endTable($assemblyTotal, $asmCount, $prodCols) }; +my $testFoot = capture { endTable($assemblyTotal, $asmCount, $testCols) }; +my $footer = capture { endHtml() }; + +open(my $pfh, '>', $prodOutFile) or die "can not write $prodOutFile"; +print $pfh $header, $prodHead, $prodBody, $prodFoot, $footer; +close($pfh); + +open(my $tfh, '>', $testOutFile) or die "can not write $testOutFile"; +print $tfh $header, $testHead, $testBody, $testFoot, $footer; +close($tfh);