3a103914b717293d5a2723c2306db20a889a31c7 hiram Fri Feb 14 14:47:59 2020 -0800 now with new path names for download access refs #23891 diff --git src/hg/makeDb/doc/asmHubs/mkAsmStats.pl src/hg/makeDb/doc/asmHubs/mkAsmStats.pl index 06f6e58..c8613fc 100755 --- src/hg/makeDb/doc/asmHubs/mkAsmStats.pl +++ src/hg/makeDb/doc/asmHubs/mkAsmStats.pl @@ -8,30 +8,32 @@ if ($argc != 2) { printf STDERR "mkAsmStats Name asmName\n"; printf STDERR "e.g.: mkAsmStats Mammals mammals\n"; exit 255; } my $Name = shift; my $asmHubName = shift; my $home = $ENV{'HOME'}; my $toolsDir = "$home/kent/src/hg/makeDb/doc/asmHubs"; my $commonNameList = "$asmHubName.asmId.commonName.tsv"; my $commonNameOrder = "$asmHubName.commonName.asmId.orderList.tsv"; my @orderList; # asmId of the assemblies in order from the *.list files # the order to read the different .list files: +my %betterName; # key is asmId, value is better common name than found in + # assembly_report my $assemblyTotal = 0; # complete list of assemblies in this group 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 @@ -91,63 +93,73 @@ ### 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; my $doneMsg = ""; if ($asmCount < $assemblyTotal) { $doneMsg = sprintf(" (%d build completed, %.2f %% finished)", $asmCount, $percentDone); } +if ($assemblyTotal > 1) { print <<"END" </tbody> <tfoot><tr><th>TOTALS:</th><td align=center colspan=3>total assembly count ${assemblyTotal}${doneMsg}</td> <td align=right>$commaSeqCount</td> <td align=right>$commaNuc</td> <td align=right>$commaGapCount</td> <td align=right>$commaGapSize</td> <td colspan=1> </td> </tr></tfoot> </table> END +} 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/asmStatsPrimates.html'>Primates</a></th>\n" if ($asmHubName ne "primates"); printf "<th><a href='../mammals/asmStatsMammals.html'>Mammals</a></th>\n" if ($asmHubName ne "mammals"); printf "<th><a href='../birds/asmStatsBirds.html'>Birds</a></th>\n" if ($asmHubName ne "birds"); printf "<th><a href='../fish/asmStatsFish.html'>Fish</a></th>\n" if ($asmHubName ne "fish"); printf "<th><a href='../vertebrate/asmStatsVertebrate.html'>other vertebrates</a></th>\n" if ($asmHubName ne "vertebrate"); 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) = @_; my ($sequenceCount, $totalSize) = split('\s+', `ave -col=2 $chromSizes | egrep "^count|^total" | awk '{printf "%d\\n", \$NF}' | xargs echo`); return ($sequenceCount, $totalSize); } @@ -174,39 +186,39 @@ 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() { foreach my $asmId (reverse(@orderList)) { + 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); - $accessionDir .= "/" . $asmId; - my $buildDir = "/hive/data/genomes/asmHubs/refseqBuild/$accessionDir"; + my $buildDir = "/hive/data/genomes/asmHubs/refseqBuild/$accessionDir/$asmId"; my $asmReport="$buildDir/download/${asmId}_assembly_report.txt"; next if (! -s "$asmReport"); - my ($gcPrefix, $asmAcc, $asmName) = split('_', $asmId, 3); my $chromSizes = "${buildDir}/${asmId}.chrom.sizes"; my $twoBit = "${buildDir}/trackData/addMask/${asmId}.masked.2bit"; next if (! -s "$twoBit"); my $faSizeTxt = "${buildDir}/${asmId}.faSize.txt"; if ( ! -s "$faSizeTxt" ) { printf STDERR "twoBitToFa $twoBit stdout | faSize stdin > $faSizeTxt\n"; print `twoBitToFa $twoBit stdout | faSize stdin > $faSizeTxt`; } my ($gapSize, $maskPerCent) = maskStats($faSizeTxt); $overallGapSize += $gapSize; my ($seqCount, $totalSize) = asmCounts($chromSizes); $overallSeqCount += $seqCount; # my $totalSize=`ave -col=2 $chromSizes | grep "^total" | awk '{printf "%d", \$NF}'`; $overallNucleotides += $totalSize; my $gapCount = gapStats($buildDir, $asmId); @@ -237,59 +249,61 @@ $bioSample =~ s/.*:\s+//; } } elsif ($line =~ m/BioProject:/) { if ($bioProject =~ m/notFound/) { ++$itemsFound; $bioProject = $line; $bioProject =~ s/.*:\s+//; } } elsif ($line =~ m/Organism name:/) { 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+\(.*//; } } elsif ($line =~ m/Taxid:/) { if ($taxId =~ m/notFound/) { ++$itemsFound; $taxId = $line; $taxId =~ s/.*:\s+//; } } } close (FH); - my $hubUrl = "https://hgdownload.soe.ucsc.edu/hubs/$accessionDir"; + 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.ucsc.edu/cgi-bin/hgGateway?hubUrl=%s/%s.hub.txt&genome=%s&position=lastDbPos' target=_blank>%s</a></td>\n", $hubUrl, $asmId, $asmId, $commonName; - printf " <td align=center><a href='https://hgdownload.soe.ucsc.edu/hubs/%s/genomes/%s/' target=_blank>%s</a></td>\n", $asmHubName, $asmId, $sciName; - printf " <td align=left><a href='https://www.ncbi.nlm.nih.gov/assembly/%s_%s/' target=_blank>%s</a></td>\n", $gcPrefix, $asmAcc, $asmId; + printf "<td align=center><a href='https://genome.ucsc.edu/cgi-bin/hgGateway?hubUrl=%s/hub.txt&genome=%s&position=lastDbPos' target=_blank>%s</a></td>\n", $hubUrl, $asmId, $commonName; + printf " <td align=center><a href='%s/' target=_blank>%s</a></td>\n", $hubUrl, $sciName; + printf " <td align=left><a href='https://www.ncbi.nlm.nih.gov/assembly/%s/' target=_blank>%s</a></td>\n", $accessionId, $asmId; printf " <td align=right>%s</td>\n", commify($seqCount); printf " <td align=right>%s</td>\n", commify($totalSize); printf " <td align=right>%s</td>\n", commify($gapCount); printf " <td align=right>%s</td>\n", commify($gapSize); printf " <td align=right>%.2f</td>\n", $maskPerCent; printf "</tr>\n"; } } ############################################################################## ### main() ############################################################################## open (FH, "<$toolsDir/${commonNameOrder}") or die "can not read ${commonNameOrder}"; while (my $line = <FH>) { chomp $line; my ($commonName, $asmId) = split('\t', $line); push @orderList, $asmId; + $betterName{$asmId} = $commonName; ++$assemblyTotal; } close (FH); startHtml(); startTable(); tableContents(); endTable(); endHtml();