735f7be7f5d6854830d5ef7cc5387fc29027778d lrnassar Thu Jun 19 18:21:58 2025 -0700 Making improvements to the panelApp otto job, primarily in support of the new panelApp Australia track. 1. Adding our standard check that will cause an exit if the item count is more than 10% different. 2. Tweak Max's code to fix an issue when specifying the specific panels to extract, this is used for the Australia genes track. 3. Remove the check that only made track items if the versions was >.99, we now make track items out of everything. 4. Max had accidentally specified the two 'onlyPanels' for Genomics England, but it should be Australia. Refs #35758 diff --git src/hg/utils/otto/panelApp/doPanelApp.py src/hg/utils/otto/panelApp/doPanelApp.py index 9d634e6af7b..07ff72860f3 100755 --- src/hg/utils/otto/panelApp/doPanelApp.py +++ src/hg/utils/otto/panelApp/doPanelApp.py @@ -1,27 +1,30 @@ #!/hive/data/outside/otto/panelApp/venv/bin/python3 + from datetime import date -import os -import shutil -import requests -import time -import json import pandas as pd -import sys -import re -import gzip -import logging +import gzip, logging, re, sys, json, time, requests, shutil, os, subprocess + +def bash(cmd): + """Run the cmd in bash subprocess""" + try: + rawBashOutput = subprocess.run(cmd, check=True, shell=True,\ + stdout=subprocess.PIPE, universal_newlines=True, stderr=subprocess.STDOUT) + bashStdoutt = rawBashOutput.stdout + except subprocess.CalledProcessError as e: + raise RuntimeError("command '{}' return with error (code {}): {}".format(e.cmd, e.returncode, e.output)) + return(bashStdoutt) def getArchDir(db): " return hgwdev archive directory given db " dateStr = date.today().strftime("%Y-%m-%d") archDir = "/usr/local/apache/htdocs-hgdownload/goldenPath/archive/%s/panelApp/%s" % (db, dateStr) if not os.path.isdir(archDir): os.makedirs(archDir) return archDir def writeBb(hg19Table, hg38Table, subTrack): " sort the pandas tables, write to BED and convert " for db in ["hg19", "hg38"]: archDir = getArchDir(db) bedFname = "current/%s/%s.bed" % (db, subTrack) @@ -52,44 +55,62 @@ def updateGbdbSymlinks(country): " update the symlinks in /gbdb. Not really necessary but kept this code just in case. " if country == "Australia": subtracks = ["genesAus", "tandRepAus", "cnvAus"] else: subtracks = ["genes", "tandRep", "cnv"] for db in ["hg19", "hg38"]: archDir = getArchDir(db) for subTrack in subtracks: if db=="hg19" and "cnv" in subTrack: continue # no cnv on hg19 cmd = "ln -sf `pwd`/current/%s/%s.bb /gbdb/%s/panelApp/%s.bb" % (db, subTrack, db, subTrack) assert(os.system(cmd)==0) +def checkIfFilesTooDifferent(oldFname,newFname): + # Exit if the difference is more than 10% + + oldItemCount = bash('bigBedInfo '+oldFname+' | grep "itemCount"') + oldItemCount = int(oldItemCount.rstrip().split("itemCount: ")[1].replace(",","")) + + newItemCount = bash('bigBedInfo '+newFname+' | grep "itemCount"') + newItemCount = int(newItemCount.rstrip().split("itemCount: ")[1].replace(",","")) + + if abs(newItemCount - oldItemCount) > 0.1 * max(newItemCount, oldItemCount): + sys.exit(f"Difference between itemCounts greater than 10%: {newItemCount}, {oldItemCount}") + else: + print(oldFname+" vs. new count: "+str(oldItemCount)+" - "+str(newItemCount)) + def flipFiles(country): " rename the .tmp files to the final filenames " if country == "Australia": subtracks = ["genesAus", "tandRepAus", "cnvAus"] else: subtracks = ["genes", "tandRep", "cnv"] for db in ["hg19", "hg38"]: archDir = getArchDir(db) for subTrack in subtracks: if db=="hg19" and "cnv" in subTrack: # no cnvs for hg19 yet continue oldFname = "current/%s/%s.bb.tmp" % (db, subTrack) newFname = "current/%s/%s.bb" % (db, subTrack) + + #Check if files are more than 10% different + checkIfFilesTooDifferent(oldFname,newFname) + os.replace(oldFname, newFname) def getAllPages(url, results=[]): " recursively download all pages. Stack should be big enough " try: myResponse = requests.get(url) if (myResponse.ok): jData = json.loads(myResponse.content.decode()) # If jData is empty create, else append if "error" in jData.keys() or not "results" in jData.keys(): raise Exception("Error in keys when downloading %s" % url) if "count" in jData and not "page" in url: print("API says that there are %d results for url %s" % (jData["count"], url)) results.extend(jData["results"]) @@ -559,30 +580,35 @@ logging.error("HTTP error on %s, retrying after 1 minute (trial %d)" % (entityUrl, count)) time.sleep(60) # Wait 1 minute before trying again count += 1 # Count the number of tries before failing if count > 10: # Quit afer 10 failed attempts assert False, "Cannot get URL after 10 attempts" jsonData = myResponse.content #jData = myResponse.json() jData = json.loads(jsonData.decode()) jsonFh.write(jsonData) jsonFh.write("\n".encode()) res = jData['results'] + + #filter by onlyPanels early, if specified + if onlyPanels is not None: + res = [entry for entry in res if entry.get('panel', {}).get('id') in onlyPanels] + num_gene_variant = len(res) count = 0 while count != num_gene_variant: temp_attribute_dictionary = dict() string_dict_key = 'gene_{}'.format(continuous_count) gene_range_37 = None gene_range_38 = None try: ensembl_genes_GRch37_82_location = res[count]['gene_data']['ensembl_genes']['GRch37']['82']['location'] location_37 = ensembl_genes_GRch37_82_location.split(':') chromo_37 = 'chr'+location_37[0] gene_range_37 = location_37[1].split('-') # on hg19, we have added a chrMT sequence later. @@ -760,37 +786,34 @@ temp_attribute_dictionary[attribute] = 'Gain-of-function' elif mode[0] == 'O' or mode[0] == 'o': temp_attribute_dictionary[attribute] = 'Other' else: temp_attribute_dictionary[attribute] = mode except: temp_attribute_dictionary[attribute] = '' #----------------------------------------------------------------------------------------------------------- panel_list = ['id','name', 'disease_group', 'disease_sub_group', 'status', 'version_created'] for attribute in panel_list: try: x = res[count]['panel'][attribute] - if attribute=="id" and onlyPanels is not None: - if int(x) not in onlyPanels: - continue if not x: temp_attribute_dictionary[attribute] = '' else: - temp_attribute_dictionary[attribute] = res[count]['panel'][attribute] + temp_attribute_dictionary[attribute] = x except: temp_attribute_dictionary[attribute] = '' #----------------------------------------------------------------------------------------------------------- version_num = 0.0 try: version_num = float(res[count]['panel']['version']) temp_attribute_dictionary['version'] = version_num except: temp_attribute_dictionary['version'] = '' #----------------------------------------------------------------------------------------------------------- try: @@ -846,61 +869,59 @@ ''' Replace all tab in value with spaces and removes new lines ''' for key, item in temp_attribute_dictionary.items(): try: if isinstance(item, int): pass elif isinstance(item, float): pass else: temp_attribute_dictionary[key] = item.replace('\t', ' ').strip().strip("\n").strip("\r") except: pass - # Version Threshold = 0.99 - max_num = float(0.99) - - if version_num > max_num: if temp_attribute_dictionary['label'] not in repeat19 and gene_range_37 is not None: # Removes Repeats repeat19.append(temp_attribute_dictionary['label']) blockSizes = int(gene_range_37[1]) - int(gene_range_37[0]) hg19_dict[string_dict_key] = [chromo_37, int(gene_range_37[0]), gene_range_37[1], temp_attribute_dictionary['label'], score, strand, gene_range_37[0], gene_range_37[1], rgb, blockCount, blockSizes, blockStarts, temp_attribute_dictionary['gene_symbol'], temp_attribute_dictionary['biotype'], temp_attribute_dictionary['hgnc_id'], temp_attribute_dictionary['gene_name'], temp_attribute_dictionary['omim_gene'], ensembl_genes_GRch38_90_ensembl_id, temp_attribute_dictionary['entity_type'], temp_attribute_dictionary['entity_name'], temp_attribute_dictionary['confidence_level'], temp_attribute_dictionary['penetrance'], temp_attribute_dictionary['mode_of_pathogenicity'], temp_attribute_dictionary['publications'], temp_attribute_dictionary['evidence'], temp_attribute_dictionary['phenotypes'], temp_attribute_dictionary['mode_of_inheritance'], temp_attribute_dictionary['tags'], temp_attribute_dictionary['id'], temp_attribute_dictionary['name'], temp_attribute_dictionary['disease_group'], temp_attribute_dictionary['disease_sub_group'], temp_attribute_dictionary['status'], temp_attribute_dictionary['version'], temp_attribute_dictionary['version_created'], temp_attribute_dictionary['relevant_disorders'], temp_attribute_dictionary['mouseOverField']] if temp_attribute_dictionary['label'] not in repeat38 and gene_range_38 is not None: # Remove Repeats repeat38.append(temp_attribute_dictionary['label']) blockSizes = int(gene_range_38[1]) - int(gene_range_38[0]) hg38_dict[string_dict_key] = [chromo_38, int(gene_range_38[0]), gene_range_38[1], temp_attribute_dictionary['label'], score, strand, gene_range_38[0], gene_range_38[1], rgb, blockCount, blockSizes, blockStarts, temp_attribute_dictionary['gene_symbol'], temp_attribute_dictionary['biotype'], temp_attribute_dictionary['hgnc_id'], temp_attribute_dictionary['gene_name'], temp_attribute_dictionary['omim_gene'], ensembl_genes_GRch38_90_ensembl_id, temp_attribute_dictionary['entity_type'], temp_attribute_dictionary['entity_name'], temp_attribute_dictionary['confidence_level'], temp_attribute_dictionary['penetrance'], temp_attribute_dictionary['mode_of_pathogenicity'], temp_attribute_dictionary['publications'], temp_attribute_dictionary['evidence'], temp_attribute_dictionary['phenotypes'], temp_attribute_dictionary['mode_of_inheritance'], temp_attribute_dictionary['tags'], temp_attribute_dictionary['id'], temp_attribute_dictionary['name'], temp_attribute_dictionary['disease_group'], temp_attribute_dictionary['disease_sub_group'], temp_attribute_dictionary['status'], temp_attribute_dictionary['version'], temp_attribute_dictionary['version_created'], temp_attribute_dictionary['relevant_disorders'], temp_attribute_dictionary['mouseOverField']] + + count = count + 1 continuous_count = continuous_count + 1 print('Genes with missing coordinates in one assembly (written to missing_genes.txt):') print(genes_missing_info) print('Genes with missing coordinates in both assemblies (written to missing_genes.txt):') missSyms = [] for miss in genes_no_location: missSyms.append(miss["gene_symbol"]) print(",".join(missSyms)) missOfh = open("missing_genes.txt", "w") missOfh.write("* Not found in one assembly:\n") missOfh.write("\n".join(genes_missing_info)) @@ -937,46 +958,46 @@ "Penetranace", "Mode of Pathogenicity", "Publications", "Evidence", "Phenotypes", "Mode of Inheritance", "Tags", "Panel ID", "Panel Name", "Disease Group", "Disease Subgroup", "Status", "Panel Version", "Version Created", "Relevant Disorders", "MouseOverField"] return pd_19_table, pd_38_table def main(): " create the 2 x three BED files and convert each to bigBed and update the archive " # the script uses relative pathnames, so make sure we're always in the right directory os.chdir("/hive/data/outside/otto/panelApp") # First update panelApp England # Gene panels track - hg19Bed, hg38Bed = downloadGenes("https://panelapp.genomicsengland.co.uk/api/v1/panels/?format=json", onlyPanels=[137, 126]) + hg19Bed, hg38Bed = downloadGenes("https://panelapp.genomicsengland.co.uk/api/v1/panels/?format=json") writeBb(hg19Bed, hg38Bed, "genes") # STRs track hg19Bed, hg38Bed = downloadTandReps("https://panelapp.genomicsengland.co.uk/api/v1/strs/?format=json") writeBb(hg19Bed, hg38Bed, "tandRep") # CNV track hg38Bed = downloadCnvs('https://panelapp.genomicsengland.co.uk/api/v1/regions/?format=json') # no hg19 CNV data yet from PanelApp - still true as of 5/20/2025 writeBb(None, hg38Bed, "cnv") flipFiles('England') updateGbdbSymlinks('England') # Now update panelApp Australia # Genes track - hg19Bed, hg38Bed = downloadGenes("https://panelapp-aus.org/api/v1/panels/?format=json") + hg19Bed, hg38Bed = downloadGenes("https://panelapp-aus.org/api/v1/panels/?format=json", onlyPanels=[137, 126]) writeBb(hg19Bed, hg38Bed, "genesAus") # STRs track hg19Bed, hg38Bed = downloadTandReps("https://panelapp-aus.org/api/v1/strs/?format=json") writeBb(hg19Bed, hg38Bed, "tandRepAus") # CNVs track hg38Bed = downloadCnvs('https://panelapp-aus.org/api/v1/regions/?format=json') # no hg19 CNV data yet from PanelApp - still true as of 5/20/2025 writeBb(None, hg38Bed, "cnvAus") flipFiles('Australia') updateGbdbSymlinks('Australia') print("PanelApp otto update: OK") main()