afee0fcf2a7b7de2eb37e1420dc126607373453e hiram Wed May 27 17:31:45 2020 -0700 rework scripts to allow subset hubs to be created from a list of assemblies refs #23734 diff --git src/hg/makeDb/doc/asmHubs/mkGenomes.pl src/hg/makeDb/doc/asmHubs/mkGenomes.pl index 0bf16e8..a00425b 100755 --- src/hg/makeDb/doc/asmHubs/mkGenomes.pl +++ src/hg/makeDb/doc/asmHubs/mkGenomes.pl @@ -1,93 +1,94 @@ #!/usr/bin/env perl use strict; use warnings; use File::Basename; my $argc = scalar(@ARGV); -if ($argc != 2) { - printf STDERR "mkGenomes.pl Name asmName\n"; - printf STDERR "e.g.: mkGenomes.pl Mammals mammals\n"; +if ($argc != 1) { + printf STDERR "mkGenomes.pl [two column name list] > .../hub/genomes.txt\n"; + printf STDERR "e.g.: mkGenomes.pl vgp.primary.assemblies.tsv > .../vgp/genomes.txt\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 "column 2: common name for species, columns separated by tab\n"; + printf STDERR "result will write a local asmId.genomes.txt file for each hub\n"; + printf STDERR "and a local asmId.hub.txt file for each hub\n"; + printf STDERR "and a local asmId.groups.txt file for each hub\n"; + printf STDERR "and the output to stdout will be the overall genomes.txt\n"; + printf STDERR "index file for all genomes in the given list\n"; exit 255; } -my $Name = shift; -my $asmHubName = shift; - -my %betterName; # key is asmId, value is common name -my $srcDocDir = "${asmHubName}AsmHub"; -my $destDir = "/hive/data/genomes/asmHubs/$asmHubName"; +my $commonNameOrder = 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"; - -open (FH, "<$toolsDir/${commonNameList}") or die "can not read $toolsDir/${commonNameList}"; -while (my $line = <FH>) { - next if ($line =~ m/^#/); - chomp $line; - my ($asmId, $name) = split('\t', $line); - $betterName{$asmId} = $name; -} -close (FH); +my %commonName; # key is asmId, value is common name my @orderList; # asmId of the assemblies in order from the *.list files # the order to read the different .list files: my $assemblyCount = 0; 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); + my ($asmId, $commonName) = split('\t', $line); + if (defined($commonName{$asmId})) { + printf STDERR "ERROR: duplicate asmId: '%s'\n", $asmId; + printf STDERR "previous name: '%s'\n", $commonName{$asmId}; + printf STDERR "duplicate name: '%s'\n", $commonName; + exit 255; + } + $commonName{$asmId} = $commonName; push @orderList, $asmId; + printf STDERR "orderList[$assemblyCount] = $asmId\n"; ++$assemblyCount; } close (FH); my $buildDone = 0; my $orderKey = 0; -foreach my $asmId (reverse(@orderList)) { +foreach my $asmId (@orderList) { ++$orderKey; my ($gcPrefix, $accession, undef) = split('_', $asmId); my $accessionId = sprintf("%s_%s", $gcPrefix, $accession); 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"; if ($gcPrefix eq "GCA") { $buildDir = "/hive/data/genomes/asmHubs/genbankBuild/$accessionDir/$asmId"; } my $asmReport="$buildDir/download/${asmId}_assembly_report.txt"; my $trackDb = "$buildDir/$asmId.trackDb.txt"; if ( ! -s "${trackDb}" ) { printf STDERR "# %03d not built yet: %s\n", $orderKey, $asmId; + printf STDERR "# '%s'\n", $trackDb; next; } if ( ! -s "${asmReport}" ) { printf STDERR "# %03d missing assembly_report: %s\n", $orderKey, $asmId; next; } ++$buildDone; printf STDERR "# %03d genomes.txt %s/%s\n", $buildDone, $accessionDir, $accessionId; my $descr=`grep -i "organism name:" $asmReport | head -1 | sed -e 's#.*organism name: *##i; s# (.*\$##;'`; chomp $descr; my $orgName=`grep -i "organism name:" $asmReport | head -1 | sed -e 's#.* name: .* (##; s#).*##;'`; chomp $orgName; - $orgName = $betterName{$asmId} if (exists($betterName{$asmId})); printf "genome %s\n", $accessionId; printf "trackDb ../%s/%s/trackDb.txt\n", $accessionDir, $accessionId; printf "groups groups.txt\n"; printf "description %s\n", $orgName; printf "twoBitPath ../%s/%s/%s.2bit\n", $accessionDir, $accessionId, $accessionId; printf "organism %s\n", $descr; my $chrName=`head -1 $buildDir/$asmId.chrom.sizes | awk '{print \$1}'`; chomp $chrName; my $bigChrom=`head -1 $buildDir/$asmId.chrom.sizes | awk '{print \$NF}'`; chomp $bigChrom; my $oneThird = int($bigChrom/3); my $tenK = $oneThird + 10000; $tenK = $bigChrom if ($tenK > $bigChrom); my $defPos="${chrName}:${oneThird}-${tenK}"; @@ -116,56 +117,56 @@ close (GF); my $localHubTxt = "$buildDir/${asmId}.hub.txt"; open (HT, ">$localHubTxt") or die "can not write to $localHubTxt"; printf HT "hub %s genome assembly\n", $accessionId; printf HT "shortLabel %s\n", $orgName; printf HT "longLabel %s/%s/%s genome assembly\n", $orgName, $descr, $asmId; printf HT "genomesFile genomes.txt\n"; printf HT "email hclawson\@ucsc.edu\n"; printf HT "descriptionUrl html/%s.description.html\n", $asmId; close (HT); my $localGroups = "$buildDir/${asmId}.groups.txt"; open (GR, ">$localGroups") or die "can not write to $localGroups"; print GR <<_EOF_ name user -label Custom +label Custom Tracks priority 1 defaultIsClosed 1 name map -label Mapping +label Mapping and Sequencing priority 2 defaultIsClosed 0 name genes -label Genes +label Genes and Gene Predictions priority 3 defaultIsClosed 0 name rna -label mRNA +label mRNA and EST priority 4 defaultIsClosed 0 name regulation label Regulation priority 5 defaultIsClosed 0 name compGeno -label Comparative +label Comparative Genomics priority 6 defaultIsClosed 0 name varRep label Variation priority 7 defaultIsClosed 0 name x label Experimental priority 10 defaultIsClosed 1 _EOF_ ; close (GR);