|
| 1 | +#!/usr/bin/env python |
| 2 | +# https://github.com/broadinstitute/gatk-sv/blob/main/scripts/inputs/preprocess_gtf.py |
| 3 | + |
| 4 | +""" |
| 5 | +Preprocess GENCODE basic GTF to extract canonical protein-coding transcripts for functional consequence annotation. |
| 6 | +""" |
| 7 | + |
| 8 | +import argparse |
| 9 | +import gzip |
| 10 | + |
| 11 | + |
| 12 | +CHROM_FIELD = 0 |
| 13 | +ELEMENT_FIELD = 2 |
| 14 | +ATTRIBUTES_FIELD = 8 |
| 15 | +TRANSCRIPT_TYPES = {"protein_coding", "nonsense_mediated_decay"} |
| 16 | +CANONICAL = {"MANE_Plus_Clinical", "MANE_Select", "Ensembl_canonical"} |
| 17 | + |
| 18 | + |
| 19 | +# Flexibly open .gz or uncompressed file to read |
| 20 | +def _open(filename): |
| 21 | + if filename.endswith(".gz"): |
| 22 | + return gzip.open(filename, 'rt') |
| 23 | + else: |
| 24 | + return open(filename, 'r') |
| 25 | + |
| 26 | + |
| 27 | +# Extract transcript type and canonical status |
| 28 | +def parse_attributes(field): |
| 29 | + # format: key1 "value1"; key2 "value2"; |
| 30 | + # keys may be repeated so cannot convert directly to dictionary |
| 31 | + attributes_list = [tuple(x.replace('"', '').split(' ')) for x in field.rstrip(";").split("; ")] |
| 32 | + protein = False |
| 33 | + canonical = False |
| 34 | + for key, val in attributes_list: |
| 35 | + if key == "tag" and val in CANONICAL: |
| 36 | + canonical = True |
| 37 | + elif key == "transcript_type" and val in TRANSCRIPT_TYPES: |
| 38 | + protein = True |
| 39 | + return protein, canonical |
| 40 | + |
| 41 | + |
| 42 | +def process(gtf, outfile): |
| 43 | + with _open(gtf) as inp, open(outfile, 'w') as out: |
| 44 | + gene_line = "" |
| 45 | + for line in inp: |
| 46 | + if line.startswith("#"): |
| 47 | + continue |
| 48 | + fields = line.rstrip('\n').split('\t') |
| 49 | + |
| 50 | + # Drop mitochondria |
| 51 | + if fields[CHROM_FIELD] == 'chrM': |
| 52 | + continue |
| 53 | + |
| 54 | + # Store gene line to print if transcript is eligible |
| 55 | + if fields[ELEMENT_FIELD] == "gene": |
| 56 | + gene_line = line |
| 57 | + continue |
| 58 | + |
| 59 | + # Select protein-coding and canonical transcripts only |
| 60 | + protein, canonical = parse_attributes(fields[ATTRIBUTES_FIELD]) |
| 61 | + if protein and canonical: |
| 62 | + out.write(gene_line + line) |
| 63 | + gene_line = "" # only print gene line before first transcript line |
| 64 | + |
| 65 | + |
| 66 | +def main(): |
| 67 | + parser = argparse.ArgumentParser() |
| 68 | + parser.add_argument('gtf', help="Input GTF from GENCODE") |
| 69 | + parser.add_argument('outfile', help="Output filename") |
| 70 | + args = parser.parse_args() |
| 71 | + |
| 72 | + process(args.gtf, args.outfile) |
| 73 | + |
| 74 | + |
| 75 | +if __name__ == '__main__': |
| 76 | + main() |
0 commit comments