a56d88e05670b759ff3b32829542537ccc790c57 lrnassar Tue Apr 28 19:18:20 2026 -0700 Address CR feedback on insight + tp53 hub scripts. refs #37418 Drop duplicated bash() wrappers in favor of subprocess.run / check_output with list args, eliminating shell=True, embedded-quote concerns, and stderr-into-stdout merging. Centralize common operations as run_sort_bed/run_liftOver in tp53FuncLib alongside existing run_bedToBigBed. Switch HTML escaping to stdlib html.escape() consistently. insightHCIPriors mouseover (previously unescaped) now escapes HGVS fields, addressing the specific c.123A>G case Jonathan flagged. Replace invalid </br> tags with <br> across all five affected mouseover sites. diff --git src/hg/makeDb/scripts/insight/insightPVS1.py src/hg/makeDb/scripts/insight/insightPVS1.py index 6cd94fcceee..b42efa10ec9 100644 --- src/hg/makeDb/scripts/insight/insightPVS1.py +++ src/hg/makeDb/scripts/insight/insightPVS1.py @@ -10,42 +10,33 @@ classifications (PVS1, PVS1_Moderate, or PVS1_n.a.) based on their position. Coordinates are computed dynamically from codon positions by querying NCBI RefSeq transcript coordinates from hgsql for both hg38 and hg19. Transcripts used: MLH1: NM_000249.4 (chr3, + strand) MSH2: NM_000251.3 (chr2, + strand) MSH6: NM_000179.3 (chr2, + strand) PMS2: NM_000535.7 (chr7, - strand) Author: Generated for InSiGHT VCEP Date: 2025 """ -import subprocess +import html import os - -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) - bashStdout = rawBashOutput.stdout - except subprocess.CalledProcessError as e: - raise RuntimeError("command '{}' return with error (code {}): {}".format(e.cmd, e.returncode, e.output)) - return(bashStdout) +import subprocess # ============================================================================ # Configuration # ============================================================================ OUTPUT_DIR = "/hive/users/lrnassar/insightHub/pvs1" # Transcripts to query TRANSCRIPTS = { 'MLH1': 'NM_000249.4', 'MSH2': 'NM_000251.3', 'MSH6': 'NM_000179.3', 'PMS2': 'NM_000535.7', } # PVS1 regions defined by codon ranges @@ -90,31 +81,31 @@ uint thickStart; "Same as chromStart" uint thickEnd; "Same as chromEnd" uint reserved; "RGB value (use R,G,B string in input file)" string rule; "Codon position rule" string acmgCode; "ACMG classification code" string _mouseOver; "Field only used as mouseOver" )""" # ============================================================================ # Functions for coordinate computation # ============================================================================ def get_transcript_info(db, accession): """Query hgsql to get transcript information from ncbiRefSeq""" query = f"SELECT name, chrom, strand, txStart, txEnd, cdsStart, cdsEnd, exonStarts, exonEnds FROM ncbiRefSeq WHERE name='{accession}'" - result = bash(f'hgsql {db} -Ne "{query}"') + result = subprocess.check_output(["hgsql", db, "-Ne", query], text=True) if not result.strip(): raise ValueError(f"Transcript {accession} not found in {db}.ncbiRefSeq") fields = result.strip().split('\t') # Parse exon starts and ends (comma-separated, trailing comma) exon_starts = [int(x) for x in fields[7].rstrip(',').split(',')] exon_ends = [int(x) for x in fields[8].rstrip(',').split(',')] return { 'name': fields[0], 'chrom': fields[1], 'strand': fields[2], 'txStart': int(fields[3]), @@ -232,33 +223,36 @@ # Clamp aa_end to protein length if aa_end > protein_length: aa_end = protein_length if strand == '+': segments = aa_to_genomic_plus(aa_start, aa_end, cds_regions) else: segments = aa_to_genomic_minus(aa_start, aa_end, cds_regions) for seg_start, seg_end, exon_num in segments: if seg_start >= seg_end: print(f" WARNING: Skipping invalid segment for {gene} {name}: {seg_start}-{seg_end}") continue - # HTML-encode special characters for mouseOver (UCSC browser can't handle ≤ ≥) - rule_html = rule.replace('≤', '≤').replace('≥', '≥').replace('>', '>').replace('<', '<') - mouse_over = f"<b>Name: </b>{name}</br><b>Gene: </b>{gene}</br><b>Rule: </b>{rule_html}</br><b>ACMG Code: </b>{acmg_code}" + # HTML-escape; UCSC mouseover doesn't render raw ≤/≥ so map to entities. + rule_html = html.escape(rule).replace('≤', '≤').replace('≥', '≥') + mouse_over = (f"<b>Name: </b>{html.escape(name)}<br>" + f"<b>Gene: </b>{html.escape(gene)}<br>" + f"<b>Rule: </b>{rule_html}<br>" + f"<b>ACMG Code: </b>{html.escape(acmg_code)}") bed_line = f"{chrom}\t{seg_start}\t{seg_end}\t{name}\t0\t.\t{seg_start}\t{seg_end}\t{color}\t{rule}\t{acmg_code}\t{mouse_over}" bed_lines.append(bed_line) return bed_lines def create_track(db, output_dir): """Create BED and bigBed files for a given genome assembly""" print(f"\n{'='*70}") print(f"Processing {db}") print(f"{'='*70}") # Query transcript info from hgsql print(f"\nQuerying transcript coordinates from {db}.ncbiRefSeq...") transcripts_info = {} for gene, accession in TRANSCRIPTS.items(): @@ -269,40 +263,43 @@ transcripts_info[gene] = tx_info # Generate BED entries print("\nGenerating BED entries...") bed_lines = generate_bed_entries(db, transcripts_info) print(f" Generated {len(bed_lines)} region segments") # Write BED file bed_file = os.path.join(output_dir, f"InSiGHTPVS1_{db}.bed") print(f"\nWriting BED file: {bed_file}") with open(bed_file, 'w') as f: f.write('\n'.join(bed_lines) + '\n') # Sort BED file print("Sorting BED file...") - bash(f"sort -k1,1 -k2,2n {bed_file} -o {bed_file}") + subprocess.run(["sort", "-k1,1", "-k2,2n", bed_file, "-o", bed_file], check=True) # Create bigBed as_file = os.path.join(output_dir, "InSiGHTPVS1.as") bb_file = os.path.join(output_dir, f"InSiGHTPVS1{db.capitalize()}.bb") chrom_sizes = f"/cluster/data/{db}/chrom.sizes" print(f"\nCreating bigBed file: {bb_file}") try: - bash(f"bedToBigBed -as={as_file} -type=bed9+3 -tab {bed_file} {chrom_sizes} {bb_file}") + subprocess.run( + ["bedToBigBed", "-as=" + as_file, "-type=bed9+3", "-tab", + bed_file, chrom_sizes, bb_file], + check=True) print(f" Successfully created: {bb_file}") except Exception as e: print(f" ERROR creating bigBed: {e}") return bed_file, bb_file # ============================================================================ # Main execution # ============================================================================ if __name__ == "__main__": print("=" * 70) print("InSiGHT VCEP PVS1 Decision Track Generator") print("=" * 70) # Create output directory if needed