6a6d7a78f1de5456925fc8f4f675795cbb61ee35
max
  Tue Nov 29 03:16:02 2022 -0800
adding tests for bigGuessDb, refs #30316

diff --git src/utils/bigGuessDb/bigGuessDb src/utils/bigGuessDb/bigGuessDb
index 419b62a..6e9241a 100755
--- src/utils/bigGuessDb/bigGuessDb
+++ src/utils/bigGuessDb/bigGuessDb
@@ -6,82 +6,88 @@
 import logging, optparse, sys
 from collections import defaultdict
 from os.path import join, isfile, expanduser
 import os, gzip, subprocess
 
 # ==== functions =====
 def parseArgs():
     " setup logging, parse command line arguments and options. -h shows auto-generated help page "
     parser = optparse.OptionParser("usage: %prog [options] filename - given a bigBed or "\
             "bigWig file, " \
         "guess the assembly based on the chrom sizes")
 
     parser.add_option("-d", "--debug", dest="debug", action="store_true", \
         help="show debug messages")
     parser.add_option("", "--index", dest="index", action="store_true", \
-        help="go /hive/data/genomes and build a file ~/.guessAssembly.txt with all chromSizes")
+        help="go /hive/data/genomes and build an index file with all chromSizes")
     parser.add_option("-b", "--best", dest="best", action="store_true", \
-        help="only print a single string, the best matching assembly")
+        help="only print a single string, the best matching assembly, or 'emptyFile' or 'notFound'. ")
+    parser.add_option("", "--indexFile", dest="indexFname", action="store", \
+        help="Use specified index file, default is %default", default="~/.guessAssembly.txt.gz")
     (options, args) = parser.parse_args()
 
     if args==[] and not options.index:
         parser.print_help()
         exit(1)
 
     if options.debug:
         logging.basicConfig(level=logging.DEBUG)
         logging.getLogger().setLevel(logging.DEBUG)
     else:
         logging.basicConfig(level=logging.INFO)
         logging.getLogger().setLevel(logging.INFO)
 
     return args, options
 
 def parseSizes(inFname, doSubset):
     " given a chrom.sizes file, return the 10 longest and 10 shortest chrom names "
     logging.info("Reading %s",inFname)
     sizes = list()
     for line in open(inFname):
         if line.startswith("#"):
             continue
         chrom, size = line.rstrip("\n").split("\t")[:2]
         sizes.append( (int(size), chrom) )
+
+    if len(sizes)==0:
+        logging.error("bigBed/bigWig file is empty. Cannot guess assembly.")
+        return None
+
     sizes.sort()
 
     if not doSubset:
         return sizes
 
-    someSizes = sizes[:10]
-    someSizes.extend(sizes[-10:])
+    someSizes = sizes[-20:] # small chroms carry less information and have fewer features
     return someSizes
 
 def writeSizes(allSizes, outFname):
     " write all sizes to the index file "
     ofh = gzip.open(outFname, "wt") # "write" "text"
     for db, dbSizes in allSizes.items():
         sizeParts = ["%s=%d" % (chrom, size) for size,chrom in dbSizes]
         sizeStr = ",".join(sizeParts)
         ofh.write("%s\t%s\n" % (db, sizeStr))
     ofh.close()
     logging.info("Wrote %s", outFname)
 
 def buildIndex(inDir, outFname):
     """ go over all direct subdirectories of inDir and find a chrom.sizes file,
     compact it to format db -> list of (chrom,size) and write to outFname """
     allSizes = dict()
     for db in os.listdir(inDir):
-        if "Patch" in db or db == "sonMus0" or db.startswith("braNey"):
+        if "Patch" in db or db == "sonMus0" or db.startswith("braNey") or db.endswith(".lowec"):
             continue
 
         subDir = join(inDir, db)
         chromFname = join(subDir, "chrom.sizes")
         if not isfile(chromFname):
             chromFname = join(subDir, db+".sizes")
 
         if not isfile(chromFname):
             print("not found "+chromFname)
             continue
 
         doSubset = True
         if db.startswith("hg") or db.startswith("mm"):
             doSubset = False
 
@@ -155,30 +161,39 @@
     dbMatches = list(dbChromMatch.items()) # dbMatches is now a list of db -> list of chromosomes
     dbMatches.sort(key=sortBySecondLen)
     return dbMatches
 
 def printAllMatches(dbChromMatch):
     " print all matching dbs as a tsv "
     print("#db\tmatchCount\tmatchList")
     for db, chromList in dbChromMatch:
         print("\t".join([db, str(len(chromList)), ",".join(chromList)]))
 
 # ----------- main --------------
 def main():
     " entry point to script "
     args, options = parseArgs()
 
-    indexFname = expanduser("~/.guessAssembly.txt.gz")
+    indexFname = expanduser(options.indexFname)
     if options.index:
         buildIndex("/hive/data/genomes", indexFname)
     else:
         inFname = args[0]
         fileSizes = bigSizes(inFname)
 
+        if len(fileSizes) == 0:
+            logging.debug("%s is empty. Cannot determine assembly." % inFname)
+            hits = ( ("emptyFile", 0),  )
+        else:
             sizeIndex = readSizeIndex(indexFname)
             hits = findBestDb(sizeIndex, fileSizes)
+
         if options.best:
+            if len(hits)==0:
+                print("notFound")
+            else:
+                #if (len(hits[0][1]) >= 2): # need more than a single match, as chrM often matches
                 print(hits[0][0])
         else:
             printAllMatches(hits)
 
 main()