gencodeGTFtoTSS.py

#! /usr/bin/env python3
#
# Original author: Peter Orchard
# Modified: 2021-03-02 (by Vivek)

import sys
import gzip
import argparse
import pathlib

parser = argparse.ArgumentParser(
    description="Create TSS Bed file from GTF annotations."
)

parser.add_argument("gtf", type=pathlib.Path, help="Path to the GTF file")
parser.add_argument(
    "-r",
    "--right",
    dest="right",
    type=int,
    default=0,
    help="Extend intervals by distance (bp) towards right (downstream). Default: 0",
)
parser.add_argument(
    "-l",
    "--left",
    dest="left",
    type=int,
    default=0,
    help="Extend intervals by distance (bp) towards left (upstream). Default: 0",
)

args = parser.parse_args()

gtf = args.gtf

open_fun, mode = open, "r"

if gtf.suffix == ".gz":
    open_fun = gzip.open
    mode = "rt"

print(
    f"Extending TSS by {args.left} bp upstream, {args.right} bp downstream.",
    file=sys.stderr,
)

with open_fun(gtf, mode) as f:
    for line in f:
        if line.startswith("#"):
            continue

        (
            chrom,
            source,
            feature_type,
            start,
            end,
            score,
            strand,
            phase,
            info,
        ) = line.rstrip().split("\t")

        if feature_type != "transcript":
            continue

        info = [i.split(" ") for i in info.rstrip(";").split("; ")]
        info = {i[0]: i[1].replace('"', "") for i in info}
        tss_start = (
            (int(start) - args.left) if strand == "+" else (int(end) - 1 - args.left)
        )
        tss_end = tss_start + (1 + args.right + args.left)
        gene_name = info["gene_name"]
        gene_status = info["gene_status"]
        gene_type = info["gene_type"]

        if gene_status == "NOVEL":
            continue

        if gene_type not in [
            "lncRNA",
            "rRNA",
            "protein_coding",
            "retained_intron",
            "processed_transcript",
            "non_coding",
            "ambiguous_orf",
            "lincRNA",
            "macro_lncRNA",
            "bidirectional_promoter_lncRNA",
        ]:
            continue

        print(
            "{chrom}\t{tss_start}\t{tss_end}\t{gene_name}\t.\t{strand}".format(
                **locals()
            ),
            file=sys.stdout,
        )