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/insightClinDomains.py src/hg/makeDb/scripts/insight/insightClinDomains.py index 3cd4301be7c..6e6fc894944 100644 --- src/hg/makeDb/scripts/insight/insightClinDomains.py +++ src/hg/makeDb/scripts/insight/insightClinDomains.py @@ -1,333 +1,326 @@ #!/usr/bin/env python3 """ InSiGHT VCEP Clinically Relevant Protein Domains Track Generator This script generates UCSC Genome Browser tracks (bigBed format) for the InSiGHT Variant Curation Expert Panel (VCEP) clinically relevant protein domains for Lynch syndrome genes: MLH1, MSH2, MSH6, and PMS2. Coordinates are computed dynamically from amino acid 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 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) - # ============================================================================ # Configuration # ============================================================================ OUTPUT_DIR = "/hive/users/lrnassar/insightHub/clinDomains" # Transcripts to query TRANSCRIPTS = { 'MLH1': 'NM_000249.4', 'MSH2': 'NM_000251.3', 'MSH6': 'NM_000179.3', 'PMS2': 'NM_000535.7', } # Domain data from proteinDomains.txt (gene, domain_name, aa_start, aa_end) DOMAINS = [ ('MLH1', 'ATPase', 20, 147), ('MLH1', 'MutSa interaction', 1, 329), ('MLH1', 'EXO1 interaction', 410, 650), ('MLH1', 'NLS', 469, 478), ('MLH1', 'MLH3/PMS1/PMS2 interaction', 506, 743), ('MLH1', 'NES', 585, 595), ('MLH1', 'NES', 650, 663), ('PMS2', 'ATPase', 1, 365), ('PMS2', 'NLS', 625, 632), ('PMS2', 'MLH1 interaction', 675, 855), ('MSH2', 'DNA binding', 1, 125), ('MSH2', 'Connector', 125, 300), ('MSH2', 'Lever', 300, 457), ('MSH2', 'Clamp', 457, 554), ('MSH2', 'Lever', 554, 620), ('MSH2', 'ATPase', 620, 856), ('MSH2', 'Helix-turn-helix', 856, 934), ('MSH2', 'MutLa interaction', 204, 225), ('MSH2', 'MSH3/MSH6 interaction', 378, 625), ('MSH2', 'MutLa interaction', 673, 686), ('MSH2', 'MSH3/MSH6 interaction', 875, 934), ('MSH2', 'EXO1 stabilisation & interaction', 266, 671), ('MSH6', 'PCNA interaction', 4, 11), ('MSH6', 'PWWP', 89, 194), ('MSH6', 'NLS', 246, 249), ('MSH6', 'NLS', 298, 313), ('MSH6', 'DNA binding', 360, 405), ('MSH6', 'Connector', 405, 701), ('MSH6', 'Lever', 701, 968), ('MSH6', 'Clamp', 968, 981), ('MSH6', 'Lever', 981, 1132), ('MSH6', 'ATPase', 1132, 1360), ('MSH6', 'MSH2 interaction', 326, 575), ('MSH6', 'MSH2 interaction', 1302, 1360), ] # AutoSQL definition for the BED9+4 format AUTOSQL = """table InSiGHTclinDomains "InSiGHT VCEP clinically relevant protein domains for Lynch syndrome genes (MLH1, MSH2, MSH6, PMS2)" ( string chrom; "Reference sequence chromosome or scaffold" uint chromStart; "Start position in chromosome" uint chromEnd; "End position in chromosome" string name; "Domain name" uint score; "Not used, all 0" char[1] strand; "Not used, all ." uint thickStart; "Same as chromStart" uint thickEnd; "Same as chromEnd" uint reserved; "RGB value (use R,G,B string in input file)" string geneSymbol; "Gene symbol" string NMaccession; "NCBI NM isoform accession" string AAlocation; "Amino acid location of domain" 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]), 'txEnd': int(fields[4]), 'cdsStart': int(fields[5]), 'cdsEnd': int(fields[6]), 'exonStarts': exon_starts, 'exonEnds': exon_ends, } def build_cds_regions(tx): """Build list of CDS regions from transcript info""" cds_regions = [] for i in range(len(tx['exonStarts'])): ex_start = tx['exonStarts'][i] ex_end = tx['exonEnds'][i] # Skip exons outside CDS if ex_end <= tx['cdsStart'] or ex_start >= tx['cdsEnd']: continue # Clip exon to CDS boundaries region_start = max(ex_start, tx['cdsStart']) region_end = min(ex_end, tx['cdsEnd']) cds_regions.append((region_start, region_end, i+1)) # For minus strand, reverse order (CDS starts from high genomic coords) if tx['strand'] == '-': cds_regions = cds_regions[::-1] return cds_regions def aa_to_genomic_plus(aa_start, aa_end, cds_regions): """Convert AA range to genomic segments for + strand genes""" nt_start = (aa_start - 1) * 3 + 1 # 1-based nucleotide position nt_end = aa_end * 3 segments = [] cumulative = 0 for start, end, exon_num in cds_regions: region_len = end - start region_nt_start = cumulative + 1 region_nt_end = cumulative + region_len if region_nt_end >= nt_start and region_nt_start <= nt_end: overlap_nt_start = max(nt_start, region_nt_start) overlap_nt_end = min(nt_end, region_nt_end) genomic_start = start + (overlap_nt_start - region_nt_start) genomic_end = start + (overlap_nt_end - region_nt_start) + 1 segments.append((genomic_start, genomic_end, exon_num)) cumulative += region_len return segments def aa_to_genomic_minus(aa_start, aa_end, cds_regions): """Convert AA range to genomic segments for - strand genes""" nt_start = (aa_start - 1) * 3 + 1 nt_end = aa_end * 3 segments = [] cumulative = 0 for start, end, exon_num in cds_regions: region_len = end - start region_nt_start = cumulative + 1 region_nt_end = cumulative + region_len if region_nt_end >= nt_start and region_nt_start <= nt_end: overlap_nt_start = max(nt_start, region_nt_start) overlap_nt_end = min(nt_end, region_nt_end) offset_from_end_start = overlap_nt_start - region_nt_start offset_from_end_end = overlap_nt_end - region_nt_start genomic_end = end - offset_from_end_start genomic_start = end - offset_from_end_end - 1 if genomic_start > genomic_end: genomic_start, genomic_end = genomic_end, genomic_start segments.append((genomic_start, genomic_end, exon_num)) cumulative += region_len return segments def generate_bed_entries(db, transcripts_info): """Generate BED entries for all domains""" bed_lines = [] for gene, domain, aa_start, aa_end in DOMAINS: tx = transcripts_info[gene] chrom = tx['chrom'] strand = tx['strand'] accession = tx['name'] cds_regions = build_cds_regions(tx) 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} {domain}: {seg_start}-{seg_end}") continue aa_loc = f"{aa_start}-{aa_end}" - mouse_over = f"<b>Domain: </b>{domain}</br><b>Gene: </b>{gene}</br><b>Transcript: </b>{accession}</br><b>Amino acid loc:</b> {aa_loc}" + mouse_over = f"<b>Domain: </b>{domain}<br><b>Gene: </b>{gene}<br><b>Transcript: </b>{accession}<br><b>Amino acid loc:</b> {aa_loc}" bed_line = f"{chrom}\t{seg_start}\t{seg_end}\t{domain}\t0\t.\t{seg_start}\t{seg_end}\t230,3,131\t{gene}\t{accession}\t{aa_loc}\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(): print(f" {gene}: {accession}") transcripts_info[gene] = get_transcript_info(db, accession) # Generate BED entries print("\nGenerating BED entries...") bed_lines = generate_bed_entries(db, transcripts_info) print(f" Generated {len(bed_lines)} domain segments") # Write BED file bed_file = os.path.join(output_dir, f"InSiGHTclinDomains_{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, "InSiGHTclinDomains.as") bb_file = os.path.join(output_dir, f"InSiGHTclinDomains{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+4 -tab {bed_file} {chrom_sizes} {bb_file}") + subprocess.run( + ["bedToBigBed", "-as=" + as_file, "-type=bed9+4", "-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 Clinically Relevant Protein Domains Track Generator") print("=" * 70) # Create output directory if needed os.makedirs(OUTPUT_DIR, exist_ok=True) # Write AutoSql file as_file = os.path.join(OUTPUT_DIR, "InSiGHTclinDomains.as") print(f"\nWriting AutoSql file: {as_file}") with open(as_file, 'w') as f: f.write(AUTOSQL) # Process both assemblies output_files = {} for db in ['hg38', 'hg19']: try: bed_file, bb_file = create_track(db, OUTPUT_DIR) output_files[db] = {'bed': bed_file, 'bb': bb_file} except Exception as e: print(f"\nERROR processing {db}: {e}") # Print summary print("\n" + "=" * 70) print("Output Files") print("=" * 70) print(f" AutoSql file: {as_file}") for db, files in output_files.items(): print(f"\n {db}:") print(f" BED file: {files['bed']}") if os.path.exists(files['bb']): print(f" BigBed file: {files['bb']}") print("\n" + "=" * 70) print("Custom Track Line (for bigBed)") print("=" * 70) print(""" track type=bigBed name="InSiGHT Domains" \\ description="InSiGHT VCEP Clinically Relevant Protein Domains" \\ visibility=pack itemRgb=on \\ bigDataUrl=<URL_TO_YOUR_BB_FILE> """) print("=" * 70) print("Done!") print("=" * 70)