8e2ce67692e814db90b2d2938915712ff9b0ee51 jeltje.van.baren Wed Apr 16 18:07:03 2025 -0700 Adding recount3 #34886 diff --git src/hg/makeDb/outside/recount3/junctionsToBed.py src/hg/makeDb/outside/recount3/junctionsToBed.py new file mode 100755 index 00000000000..0600aa0cbb9 --- /dev/null +++ src/hg/makeDb/outside/recount3/junctionsToBed.py @@ -0,0 +1,125 @@ +#!/usr/bin/env python3 +import sys +import csv +import os +import re +import csv +import argparse +import textwrap +from collections import namedtuple +from math import log +sys.path.append('/hive/groups/browser/pycbio/lib') +from pycbio.hgdata.bed import Bed, BedBlock, BedReader, intArraySplit + + + +def main(): + parser = argparse.ArgumentParser( + formatter_class=argparse.RawDescriptionHelpFormatter, + description=textwrap.dedent('''\ +Converts gff to bed format, adding expression values and CDS from separate files. +Also outputs a decorator bed file, see https://genome.ucsc.edu/goldenPath/help/decorator.html + + ''')) + #group = parser.add_argument_group('required arguments') + parser.add_argument('--junctions', required=True, type=str, help='annotated gff') + parser.add_argument('--bed', required=True, type=str, help='output bed') + parser.add_argument('--decorator', required=True, type=str, help='output decorator bed') + parser.add_argument('--compilation', required=True, type=str, + choices=['ccle', 'gtexv2', 'srav3h', 'tcgav2', 'srav1m'], + help='one of ccle, gtexv2, srav3h, tcgav2, srav1m') + args = parser.parse_args() + + junction_to_bed(args.bed, args.junctions, args.decorator, args.compilation) + + +def makeSpliceDecorator(bedObj, donor, acc, outFH): + '''Create decorators for splice consensus''' + #chr1 1999 2000 green_circle 0 + 1999 2000 0,255,0,255 1 1 0 chr1:1000-2000:feature glyph 0,255,0,255 Circle + #colorblind color scheme Paul Tol's Muted from + # https://www.nceas.ucsb.edu/sites/default/files/2022-06/Colorblind%20Safe%20Color%20Schemes.pdf + # dark blue = GT donor, AG acceptor (046,037,133), teal= GC donor (093,168,153), faded red = AG/AC (194,106,119) + if bedObj.strand == '.': + # do not decorate + return + # assign colors + splice_map_donor = {'GT': '046,037,133','GC': '093,168,153','AT': '194,106,119'} + splice_map_acc = {'AG': '046,037,133', 'AC': '194,106,119'} + leftColor = splice_map_donor.get(donor, 'default') + rightColor = splice_map_acc.get(acc, 'default') + # correct for negative strand + if bedObj.strand == '-': + leftColor, rightColor = rightColor, leftColor + # decorate the first and last two bases of each block + leftChromStart = bedObj.chromStart + leftChromEnd = leftChromStart + 2 + leftBlocks = [BedBlock(bedObj.chromStart, bedObj.chromStart+2)] + rightChromStart = bedObj.chromEnd - 2 + rightChromEnd = bedObj.chromEnd + rightBlocks = [BedBlock(bedObj.chromEnd-2, bedObj.chromEnd)] + # format the target location (required for format) + target='{}:{}-{}:{}'.format(bedObj.chrom, bedObj.chromStart, bedObj.chromEnd, bedObj.name) + extraCols = [target, 'block', leftColor, 'Circle'] + # create and write bed objects, one for each side + leftObj = Bed(bedObj.chrom, leftChromStart, leftChromEnd, name='decl'+bedObj.name, + strand=bedObj.strand, blocks=leftBlocks, numStdCols=12, itemRgb=leftColor, + extraCols=extraCols) + leftObj.write(outFH) + extraCols = [target, 'block', rightColor, 'Circle'] + rightObj = Bed(bedObj.chrom, rightChromStart, rightChromEnd, name='decr'+bedObj.name, + strand=bedObj.strand, blocks=rightBlocks, numStdCols=12, itemRgb=rightColor, + extraCols=extraCols) + rightObj.write(outFH) + + +def junction_to_bed(outputfile, junctions, decoratorfile, compilation): + '''Parses intron lines, extracts ID info from last field, assigns score''' + urlbase = '<a href="https://snaptron.cs.jhu.edu/snaptron-studies/jxn2studies?compilation='+compilation + urlfill = '&jid={jid}&coords={chrom}:{start}-{end}" target="_blank">link to snaptron</a>' +#'https://snaptron.cs.jhu.edu/snaptron-studies/jxn2studies?compilation=<COMPILATION_ID>&jid=<JXN_ID>&coords=<CHROMOSOME:START-END>' +# url = '<a href="https://snaptron.cs.jhu.edu/gtexv2/snaptron?regions={chrom}:{start}-{end}" target="_blank">{text}</a>' + decFile = open(decoratorfile, 'wt') + with open(outputfile, 'wt') as outfile: + writer = csv.writer(outfile, delimiter='\t', lineterminator=os.linesep) + # extract all exons from the gff, keep exons grouped by transcript + maxScore=1 + for line in open(junctions, 'rt'): + if line.startswith('#'): + continue + [snaptron_id, chrom, start, end, length, strand, annotated, + left_motif, right_motif, left_annotated, right_annotated, samples, + samples_count, coverage_sum, coverage_avg, coverage_median, + source_dataset_id] = line.rstrip().split('\t') + # skip non-genome chromosomes + if any(chrom.startswith(prefix) for prefix in ('ERCC', 'SIRV', 'GL', 'JH', 'chrEBV')): + continue + # turn left and right motifs into donor and acceptor + donor = left_motif.upper() + acc = right_motif.upper() + if strand == '?': + strand = '.' + elif strand == '-': + complement = {'A': 'T', 'T': 'A', 'C': 'G', 'G': 'C'} + donor = ''.join(complement[base] for base in reversed(acc)) + acc = ''.join(complement[base] for base in reversed(left_motif.upper())) + # link out to sample information for each junction + link = urlbase + urlfill.format(jid=snaptron_id, chrom=chrom, start=start, end=end) + # correct coordinates + start = int(start)-1 + end = int(end) + # score scaling + score = int(10*log(int(coverage_sum))) # score field only accepts integers + maxScore = max(score, maxScore) + extraCols = [coverage_sum, samples_count, donor, acc, link] + # create bed object + bedObj = Bed(chrom, start, end, name=snaptron_id, strand=strand, score=score, + thickStart=start, thickEnd=end, numStdCols=9, extraCols=extraCols) + # create splice site decorator file + makeSpliceDecorator(bedObj, donor, acc, decFile) + bedObj.write(outfile) + decFile.close() + print('maxScore', maxScore) + + +if __name__ == "__main__": + main()