dff254260c6b5cb3ab0f2257775fce47bbbfa3ea hiram Wed Mar 25 11:56:34 2020 -0700 now staging on hgdownload.soe.ucsc.edu refs #24748 diff --git src/hg/makeDb/doc/asmHubs/trackData.pl src/hg/makeDb/doc/asmHubs/trackData.pl index 326e842..937d6d0 100755 --- src/hg/makeDb/doc/asmHubs/trackData.pl +++ src/hg/makeDb/doc/asmHubs/trackData.pl @@ -27,102 +27,103 @@ 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 } -# ($itemCount, $percentCover) = oneTrackData($trackFile, $sizeNoGaps, $trackFb); +# ($itemCount, $percentCover) = oneTrackData($asmId, $track, $trackFile, $totalSize, $trackFb, $runDir); # might have a track feature bits file (trackFb), maybe not sub oneTrackData($$$$$$) { my ($asmId, $trackName, $file, $genomeSize, $trackFb, $runDir) = @_; # printf STDERR "# %s\n", $file; my $itemCount = 0; my $percentCover = 0; if (! -s "${file}") { if ($trackName eq "gapOverlap") { if (-s "${runDir}/$asmId.gapOverlap.bed.gz" ) { my $lineCount=`zcat "${runDir}/$asmId.gapOverlap.bed.gz" | head | wc -l`; chomp $lineCount; if (0 == $lineCount) { return("0", "0 %"); } else { return("n/a", "n/a"); } } + } elsif ($trackName eq "gap") { + return("0", "0 %"); } else { return("n/a", "n/a"); } } if ($file =~ m/.bw$/) { my $bigWigInfo = `bigWigInfo "$file" | egrep "basesCovered:|mean:" | awk '{print \$NF}' | xargs echo | sed -e 's/,//g;'`; chomp $bigWigInfo; my ($bases, $mean) = split('\s+', $bigWigInfo); $percentCover = sprintf("%.2f %%", 100.0 * $bases / $genomeSize); $itemCount = sprintf ("%.2f", $mean); # printf STDERR "# bigWigInfo %s %s %s\n", $itemCount, $percentCover, $file; } else { 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); # 56992654 bases of 2616369673 (2.178%) in intersection if ( -s "${trackFb}" ) { -printf STDERR "# $trackFb\n"; my ($itemBases, undef, undef, $noGapSize, undef) = split('\s+', `cat $trackFb`, 5); $percentCover = sprintf("%.2f %%", 100.0 * $itemBases / $noGapSize); } # printf STDERR "# bigBedInfo %s %s %s\n", $itemCount, $percentCover, $file; } return ($itemCount, $percentCover); } ############################################################################## ### start the HTML output ############################################################################## 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"; } print <<"END" <!DOCTYPE HTML 4.01 Transitional> -<!--#set var="TITLE" value="$Name genomes assembly hubs" --> +<!--#set var="TITLE" value="$Name genomes assembly hubs, track statistics" --> <!--#set var="ROOT" value="../.." --> <!--#include virtual="\$ROOT/inc/gbPageStartHardcoded.html" --> -<h1>$Name Genomes assembly hubs</h1> +<h1>$Name Genomes assembly hubs, track statistics</h1> <p> Assemblies from NCBI/Genbank/Refseq sources, $subSetMessage. </p> -<h3>See also: <a href='index.html'>hub access</a></h3><br> +<h3>See also: <a href='index.html'>hub access</a>, <a href='asmStats$Name.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 } ############################################################################## ### start the table output ############################################################################## sub startTable() { print <<"END" @@ -131,31 +132,32 @@ <th>common name<br>link to genome browser</th> <th class="sorttable_numeric">gc5 base</th> <th class="sorttable_numeric">AGP<br>gap</th> <th class="sorttable_numeric">all<br>gaps</th> <th class="sorttable_numeric">assembly<br>sequences</th> <th class="sorttable_numeric">rmsk</th> <th class="sorttable_numeric">TRF<br>simpleRepeat</th> <th class="sorttable_numeric">window<br>Masker</th> <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> <th class="sorttable_numeric">genes<br>ncbi</th> <th class="sorttable_numeric">ncbiRefSeq</th> <th class="sorttable_numeric">xenoRefGene</th> - <th class="sorttable_numeric">augustus</th> + <th class="sorttable_numeric">augustus<br>genes</th> + <th class="sorttable_numeric">Ensembl<br>genes</th> </tr></thead><tbody> END } ############################################################################## ### end the table output ############################################################################## sub endTable() { my $commaNuc = commify($overallNucleotides); my $commaSeqCount = commify($overallSeqCount); my $commaGapSize = commify($overallGapSize); my $commaGapCount = commify($overallGapCount); my $percentDone = 100.0 * $asmCount / $assemblyTotal; @@ -171,42 +173,53 @@ } else { print <<"END" </tbody> </table> END } } # sub endTable() ############################################################################## ### end the HTML output ############################################################################## sub endHtml() { if ($asmHubName ne "viral") { - printf "<p>\nOther assembly hubs available:<br>\n<table border='1'><thead>\n<tr>"; - - printf "<th><a href='../primates/trackData.html'>Primates</a></th>\n" - if ($asmHubName ne "primates"); - printf "<th><a href='../mammals/trackData.html'>Mammals</a></th>\n" - if ($asmHubName ne "mammals"); - printf "<th><a href='../birds/trackData.html'>Birds</a></th>\n" - if ($asmHubName ne "birds"); - printf "<th><a href='../fish/trackData.html'>Fish</a></th>\n" - if ($asmHubName ne "fish"); - printf "<th><a href='../vertebrate/trackData.html'>other vertebrates</a></th>\n" - if ($asmHubName ne "vertebrate"); + printf "<p>\n<table border='1'><thead>\n<tr>"; + printf "<th>Assembly hubs index pages: </th>\n"; + printf "<th><a href='../primates/index.html'>Primates</a></th>\n"; + printf "<th><a href='../mammals/index.html'>Mammals</a></th>\n"; + printf "<th><a href='../birds/index.html'>Birds</a></th>\n"; + printf "<th><a href='../fish/index.html'>Fish</a></th>\n"; + printf "<th><a href='../vertebrate/index.html'>other vertebrates</a></th>\n"; + + printf "</tr><tr>\n"; + printf "<th>Hubs assembly statistics: </th>\n"; + printf "<th><a href='../primates/asmStatsPrimates.html'>Primates</a></th>\n"; + printf "<th><a href='../mammals/asmStatsMammals.html'>Mammals</a></th>\n"; + printf "<th><a href='../birds/asmStatsBirds.html'>Birds</a></th>\n"; + printf "<th><a href='../fish/asmStatsFish.html'>Fish</a></th>\n"; + printf "<th><a href='../vertebrate/asmStatsVertebrate.html'>other vertebrates</a></th>\n"; + + printf "</tr><tr>\n"; + printf "<th>Hubs track statistics: </th>\n"; + printf "<th><a href='../primates/trackData.html'>Primates</a></th>\n"; + printf "<th><a href='../mammals/trackData.html'>Mammals</a></th>\n"; + printf "<th><a href='../birds/trackData.html'>Birds</a></th>\n"; + printf "<th><a href='../fish/trackData.html'>Fish</a></th>\n"; + printf "<th><a href='../vertebrate/trackData.html'>other vertebrates</a></th>\n"; printf "</tr></thead>\n</table>\n</p>\n"; } print <<"END" </div><!-- closing gbsPage from gbPageStartHardcoded.html --> </div><!-- closing container-fluid from gbPageStartHardcoded.html --> <!--#include virtual="\$ROOT/inc/gbFooterHardcoded.html"--> <script type="text/javascript" src="/js/sorttable.js"></script> </body></html> END } sub asmCounts($) { my ($chromSizes) = @_; @@ -236,34 +249,36 @@ my ($buildDir, $asmId) = @_; my $gapBed = "$buildDir/trackData/allGaps/$asmId.allGaps.bed.gz"; my $gapCount = 0; if ( -s "$gapBed" ) { $gapCount = `zcat $gapBed | awk '{print \$3-\$2}' | ave stdin | grep '^count' | awk '{print \$2}'`; } chomp $gapCount; return ($gapCount); } ############################################################################## ### tableContents() ############################################################################## sub tableContents() { - my @trackList = qw(gc5Base gap allGaps assembly rmsk simpleRepeat windowMasker gapOverlap tandemDups cpgIslandExtUnmasked cpgIslandExt ncbiGene ncbiRefSeq xenoRefGene augustus); + my @trackList = qw(gc5Base gap allGaps assembly rmsk simpleRepeat windowMasker gapOverlap tandemDups cpgIslandExtUnmasked cpgIslandExt ncbiGene ncbiRefSeq xenoRefGene augustus ensGene); + my $asmCounted = 0; foreach my $asmId (reverse(@orderList)) { + my $tracksCounted = 0; my ($gcPrefix, $asmAcc, $asmName) = split('_', $asmId, 3); my $accessionId = sprintf("%s_%s", $gcPrefix, $asmAcc); my $accessionDir = substr($asmId, 0 ,3); $accessionDir .= "/" . substr($asmId, 4 ,3); $accessionDir .= "/" . substr($asmId, 7 ,3); $accessionDir .= "/" . substr($asmId, 10 ,3); my $buildDir = "/hive/data/genomes/asmHubs/refseqBuild/$accessionDir/$asmId"; # next if ($asmId ne "GCF_000001405.39_GRCh38.p13"); my $asmReport="$buildDir/download/${asmId}_assembly_report.txt"; if (! -s "$asmReport") { printf STDERR "# no assembly report:\n# %s\n", $asmReport; next; } my $chromSizes = "${buildDir}/${asmId}.chrom.sizes"; my $twoBit = "${buildDir}/trackData/addMask/${asmId}.masked.2bit"; @@ -308,70 +323,83 @@ if ($sciName =~ m/notFound/) { ++$itemsFound; $commonName = $line; $sciName = $line; $commonName =~ s/.*\(//; $commonName =~ s/\)//; $commonName = $betterName{$asmId} if (exists($betterName{$asmId})); $sciName =~ s/.*:\s+//; $sciName =~ s/\s+\(.*//; } } } close (FH); my $hubUrl = "https://hgdownload.soe.ucsc.edu/hubs/$accessionDir/$accessionId"; printf "<tr><td align=right>%d</td>\n", ++$asmCount; - printf "<td align=center><a href='https://genome-test.gi.ucsc.edu/h/%s' target=_blank>%s<br>%s</a></td>\n", $accessionId, $commonName, $accessionId; + printf "<td align=center><a href='https://genome.ucsc.edu/h/%s' target=_blank>%s<br>%s</a></td>\n", $accessionId, $commonName, $accessionId; foreach my $track (@trackList) { my $trackFile = "$buildDir/bbi/$asmId.$track"; my $trackFb = "$buildDir/trackData/$track/fb.$asmId.$track.txt"; my $runDir = "$buildDir/trackData/$track"; my ($itemCount, $percentCover); if ( "$track" eq "gc5Base" ) { $trackFile .= ".bw"; } else { $trackFile .= ".bb"; } + my $customKey = ""; 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 = `zcat "$faOut" | wc -l`; 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`); - $percentCover = sprintf("%.2f %%", $percentCover); chomp $percentCover; + $customKey = sprintf("%.2f", $percentCover); + $percentCover = sprintf("%.2f %%", $percentCover); } } else { ($itemCount, $percentCover) = oneTrackData($asmId, $track, $trackFile, $totalSize, $trackFb, $runDir); + if (($percentCover =~ m/%/) || ($percentCover !~ m#n/a#)) { + $customKey = $percentCover; + $customKey =~ s/[ %]+//; } + } + if (length($customKey)) { + printf " <td align=right sorttable_customkey='%s'>%s<br>(%s)</td>\n", $customKey, $itemCount, $percentCover; + } else { printf " <td align=right>%s<br>(%s)</td>\n", $itemCount, $percentCover; } + $tracksCounted += 1 if ($itemCount ne "n/a"); + } printf "</tr>\n"; + $asmCounted += 1; + printf STDERR "# %03d\t%02d tracks\t%s\n", $asmCounted, $tracksCounted, $asmId; } } ############################################################################## ### main() ############################################################################## open (FH, "<$toolsDir/${commonNameOrder}") or die "can not read ${commonNameOrder}"; while (my $line = <FH>) { next if ($line =~ m/^#/); chomp $line; my ($commonName, $asmId) = split('\t', $line); push @orderList, $asmId; $betterName{$asmId} = $commonName; ++$assemblyTotal;