f9ce722d011fb87c9f61a4d892b9140ef7fa9047 hiram Mon Sep 12 20:38:56 2022 -0700 now working through step chromAlias and fixup that procedure to work correctly on ucsc names refs #29811 diff --git src/hg/utils/automation/asmHubChromAlias.pl src/hg/utils/automation/asmHubChromAlias.pl index 99c7af2..0f5318d 100755 --- src/hg/utils/automation/asmHubChromAlias.pl +++ src/hg/utils/automation/asmHubChromAlias.pl @@ -1,29 +1,29 @@ #!/usr/bin/env perl use strict; use warnings; my $argc = scalar(@ARGV); if ($argc != 1) { printf STDERR "usage: asmHubChromAlias.pl asmId > asmId.chromAlias.tab\n\n"; printf STDERR "where asmId is something like: GCF_000001735.4_TAIR10.1\n"; printf STDERR "Outputs a tab file suitable for processing with ixIxx to\n"; printf STDERR "create an index file to use in an assembly hub.\n"; printf STDERR "This command assumes it is in a work directory in the\n"; printf STDERR "assembly hub: .../asmId/trackData/chromAlias/\n"; - printf STDERR "and .../asmId/downloads/ and .../asmId/sequence/ have\n"; + printf STDERR "and ../../download/ and ../../sequence/ have\n"; printf STDERR "been created and processed for the hub build.\n"; exit 255; } my %ucscToRefSeq; # key is UCSC sequence name, value is RefSeq name my %ucscToGenbank; # key is UCSC sequence name, value is GenBank name my $ucscNames = 0; # == 1 if sequence is UCSC names, == 0 if NCBI names my $dupCount = 0; my %dupToSequence; # key is duplicate name, value is target sequence name # to manage duplicates coming in from assembly report my $sequenceCount = 0; my %sequenceSizes; # key is sequence name, value is sequence size my $asmId = shift; @@ -160,31 +160,32 @@ chomp $line; my ($ucscName, $seqName) = split('\s+', $line); $ncbiToUcsc{$seqName} = $ucscName; $ucscToNcbi{$ucscName} = $seqName; ++$nameCount; $ucscNames = 1 if (defined($sequenceSizes{$ucscName})); if ($isRefSeq) { $ucscToRefSeq{$ucscName} = $seqName; } else { $ucscToGenbank{$ucscName} = $seqName; } } close (FH); my $dupsNotFound = 0; -my $dupsList = "../../download/$asmId.dups.txt.gz"; +my $dupsList = `ls ../../download/*.dups.txt.gz 2> /dev/null|head -1`; +chomp $dupsList if (defined($dupsList)); if ( -s "$dupsList" ) { open (FH, "zcat $dupsList | awk '{print \$1, \$3}'|") or die "can not read $dupsList"; while (my $line = ) { chomp $line; my ($dupAlias, $dupTarget) = split('\s+', $line); $dupToSequence{$dupAlias} = $dupTarget; if ($ucscNames) { if (!defined($ncbiToUcsc{$dupTarget})) { printf STDERR "# ERROR: can not find dupTarget: $dupTarget in ncbiToUcsc for dupAlias: $dupAlias\n"; $dupsNotFound += 1; } else { printf STDERR "# addAlias(\"ucsc\", $dupAlias, $ncbiToUcsc{$dupTarget});\n"; } } else { printf STDERR "# addAlias($asmSource, $dupAlias, $dupTarget);\n"; @@ -247,79 +248,90 @@ } if (!defined($alias)) { printf STDERR "ERROR: can not find alias name for '%s' sequence\n", $seqName; exit 255; } addAlias($ucscNames ? "ucsc" : $asmSource, $alias, $seqName); } # foreach my $sequence (sort keys %chrNames) { # printf "%s\t%s\n", $chrNames{$sequence}, $sequence; # } ### next set of names are the equivalents declared by NCBI ### if they exist my %chr2acc; # key is sequence name, value is NCBI chromosome name -my $asmStructCount = `ls ../../download/${asmId}_assembly_structure/*/*/chr2acc 2> /dev/null | wc -l`; +my $asmStructCount = `ls ../../download/*_assembly_structure/*/*/chr2acc 2> /dev/null | wc -l`; chomp $asmStructCount; if ( $asmStructCount > 0 ) { # printf STDERR "# second name set processing chr2acc files\n"; - open (FH, "grep -h -v '^#' ../../download/${asmId}_assembly_structure/*/*/chr2acc|") or die "can not grep chr2acc files"; + open (FH, "grep -h -v '^#' ../../download/*_assembly_structure/*/*/chr2acc|") or die "can not grep chr2acc files"; while (my $line = ) { chomp $line; my ($alias, $seqName) = split('\t', $line); $alias =~ s/ /_/g; # some assemblies have spaces in chr names ... $alias =~ s/:/_/g; # one assembly had : in a chr name $chr2acc{$seqName} = $alias; if ($ucscNames) { $seqName = $ncbiToUcsc{$seqName}; } if (defined($dupToSequence{$seqName})) { # avoid duplicates printf STDERR "# skipping duplicate name $seqName\n"; } else { addAlias("ncbi", $alias, $seqName); } } close (FH); } # foreach my $sequence (sort keys %chr2acc) { # printf "%s\t%s\n", $chr2acc{$sequence}, $sequence; # } my $dbgCount = 0; # printf STDERR "# third set processing assembly_report\n"; # column 1 is the 'assembly' name # column 5 is the GenBank-Accn, column 7 is the RefSeq-Accn -open (FH, "grep -v '^#' ../../download/${asmId}_assembly_report.txt | cut -d\$'\t' -f1,5,7|") or die "can not grep assembly_report"; +open (FH, "grep -v '^#' ../../download/*_assembly_report.txt | cut -d\$'\t' -f1,5,7|") or die "can not grep assembly_report"; while (my $line = ) { chomp $line; ++$dbgCount; - my ($asmName, $gbkName, $refSeqName) = split('\s+', $line); + my ($asmName, $gbkName, $refSeqName) = split('\t', $line); + $asmName =~ s/ /_/g; # some assemblies have spaces in chr names ... + $asmName =~ s/:/_/g; # one assembly had : in chr name if (defined($dupToSequence{$asmName})) { # avoid duplicates printf STDERR "# skipping duplicate name $asmName\n"; next; } elsif (defined($dupToSequence{$gbkName})) { # avoid duplicates printf STDERR "# skipping duplicate name $gbkName\n"; next; } elsif (defined($dupToSequence{$refSeqName})) { # avoid duplicates printf STDERR "# skipping duplicate name $refSeqName\n"; next; } - printf STDERR "# '%s'\t'%s'\t'%s'\n", $asmName, $gbkName, $refSeqName if ($dbgCount < 5); + printf STDERR "# asmRpt: '%s'\t'%s'\t'%s'\n", $asmName, $gbkName, $refSeqName if ($dbgCount < 5); # next if ($refSeqName eq "na"); # may not be any RefSeq name # next if ($gbkName eq "na"); # may not be any GenBank name + # fill in ncbiToUcsc for potentially the 'other' NCBI name + if (defined($ncbiToUcsc{$refSeqName}) && !defined($ncbiToUcsc{$gbkName})) { + $ncbiToUcsc{$gbkName} = $ncbiToUcsc{$refSeqName}; + $ucscToNcbi{$ncbiToUcsc{$refSeqName}} = $gbkName; + } + if (defined($ncbiToUcsc{$gbkName}) && !defined($ncbiToUcsc{$refSeqName})) { + $ncbiToUcsc{$refSeqName} = $ncbiToUcsc{$gbkName}; + $ucscToNcbi{$ncbiToUcsc{$gbkName}} = $refSeqName; + } if ($refSeqName ne "na") { my $seqName = $refSeqName; if (! $isRefSeq) { $seqName = $gbkName; } if ($ucscNames) { $seqName = $ncbiToUcsc{$seqName}; } if (!defined($seqName)) { if (defined($aliasOut{"refseq"})) { if (defined($aliasOut{"refseq"}->{$refSeqName})) { $seqName = $aliasOut{"refseq"}->{$refSeqName}; } } }