get_gene_structures.cgi

#!/usr/bin/python3
import base64
import cgi
import json
import urllib.request
from PIL import Image, ImageDraw

print('Access-Control-Allow-Origin: *')
print('Content-Type: application/json\n')  # HTML is following

# ----- CONSTANTS -----
EXON_IMG_WIDTH = 450
EXON_IMG_HEIGHT = 8

# ----- VARIABLES -----
exon_and_CDS = {"start": [], "end": []}
mRNA = {"start": [], "end": []}

expression_score = []
variants = []
variant_count = 0

# ----- GET THE LOCUS OF INTEREST -----
gene_id = cgi.FieldStorage().getvalue('locus')

# ----- GETS MAPPING INFO FOR THE GENE ID -----
map_info = json.loads(urllib.request.urlopen(
    "https://bar.utoronto.ca/webservices/bar_araport/gene_structure_by_locus.php?locus=" + gene_id).read())

printout = ""
printout = printout + "{"
printout = printout + "\"locus\" : \"" + gene_id + "\", "
printout = printout + "\"locus_start\" : \"" + str(map_info['features'][0]['start']) + "\", "
printout = printout + "\"locus_end\" : \"" + str(map_info['features'][0]['end']) + "\", "
printout = printout + "\"splice_variants\" : ["

# Get start/end of the LOCUS
start = int(map_info['features'][0]['start'])
end = int(map_info['features'][0]['end'])
strand = int(map_info['features'][0]['strand'])

i = 0
for subfeature in map_info['features'][0]['subfeatures']:
    # Need a comma if it is not the first element
    if i == 0:
        printout = printout + "{"
    else:
        printout = printout + ", {"

    # To return a count of the number of variants returned by Araport
    variant_count = variant_count + 1

    # Extract variant
    variants.append(subfeature['subfeatures'])

    # Return all exons' coordinates
    printout = printout + "\"exon_coordinates\" : ["

    # Define the colours
    white = (255, 255, 255)
    black = (0, 0, 0)

    red = (220, 20, 60)
    orange = (255, 140, 0)
    blue = (0, 0, 255)
    # TO DO: Fix the green and dark green shades...
    green = (166, 220, 166)
    darkgreen = (0, 125, 0)

    # Generate gene structure image based on the information in map_info.
    # Create an image for gene structure
    exon_graph_image = Image.new("RGB", (EXON_IMG_WIDTH, EXON_IMG_HEIGHT), white)
    exongraph = ImageDraw.Draw(exon_graph_image)

    count = 0  # Need a comma if it is not the first element...
    for region in variants[i]:
        # We want to return regions marked as type = exons and type = CDS
        if region['type'] == 'exon' or region['type'] == 'CDS':
            exon_and_CDS["start"].append(int(region['start']))
            exon_and_CDS["end"].append(int(region['end']))
            if count == 0:
                printout = printout + "{" + "\"exon_start\" : " + str(int(region['start'])) + ", \"exon_end\" : " + str(
                    int(region['end'])) + "}"
            else:
                printout = printout + ", {" + "\"exon_start\" : " + str(
                    int(region['start'])) + ", \"exon_end\" : " + str(int(region['end'])) + "}"
            count = count + 1  # To add a comma only...
        # We want to graph all types of features in the gene structure image
        if region['type'] == 'exon':
            exongraph.rectangle((
                (int(float(region['start'] - start) / (end - start) * EXON_IMG_WIDTH), 0),
                (int(float(region['end'] - start) / (end - start) * EXON_IMG_WIDTH), EXON_IMG_HEIGHT)), darkgreen)
        elif region['type'] == 'CDS':
            exongraph.rectangle((
                (int(float(region['start'] - start) / (end - start) * EXON_IMG_WIDTH), 0),
                (int(float(region['end'] - start) / (end - start) * EXON_IMG_WIDTH), EXON_IMG_HEIGHT)), darkgreen)
        elif region['type'] == 'five_prime_UTR':
            exongraph.rectangle((
                (int(float(region['start'] - start) / (end - start) * EXON_IMG_WIDTH), 0),
                (int(float(region['end'] - start) / (end - start) * EXON_IMG_WIDTH), EXON_IMG_HEIGHT)), green)
        elif region['type'] == 'three_prime_UTR':
            exongraph.rectangle((
                (int(float(region['start'] - start) / (end - start) * EXON_IMG_WIDTH), 0),
                (int(float(region['end'] - start) / (end - start) * EXON_IMG_WIDTH), EXON_IMG_HEIGHT)), green)

    # Nucleotide padding
    nucleotidePadding = 100
    exongraph.rectangle(((0, 0), ((EXON_IMG_WIDTH / nucleotidePadding), EXON_IMG_HEIGHT)), white)
    exongraph.rectangle((((EXON_IMG_WIDTH - (EXON_IMG_WIDTH / nucleotidePadding)), 0),
                         (EXON_IMG_WIDTH, EXON_IMG_HEIGHT)), white)

    # Line in the middle
    exongraph.rectangle(((0, EXON_IMG_HEIGHT / 2), (EXON_IMG_WIDTH, EXON_IMG_HEIGHT / 2)), black)

    # Insert appropriate arrows to indicate direction of the gene
    if strand == -1:
        exongraph.polygon(
            ((0, EXON_IMG_HEIGHT / 2), (3, EXON_IMG_HEIGHT - (EXON_IMG_HEIGHT - 1)), (3, EXON_IMG_HEIGHT - 1)), black)
    elif strand == +1:
        exongraph.polygon(((EXON_IMG_WIDTH, EXON_IMG_HEIGHT / 2),
                                 (EXON_IMG_WIDTH - 3, EXON_IMG_HEIGHT - (EXON_IMG_HEIGHT - 1)),
                                 (EXON_IMG_WIDTH - 3, EXON_IMG_HEIGHT - 1)), black)

    f = open("get_exon_base64_exongraph.png", "wb")
    exon_graph_image.save(f)
    f.close()

    printout = printout + "], " + "\"start\" : " + str(start) + ", " + "\"end\" : " + str(
        end) + ", " + "\"gene_structure\" : "
    printout = printout + "\""

    with open("get_exon_base64_exongraph.png", "rb") as fl:
        printout = printout + base64.b64encode(fl.read()).decode("utf-8")
    printout = printout + "\""
    fl.close()

    i = i + 1
    printout = printout + "}"

printout = printout + "]" + ", \"variant_count\" : \"" + str(variant_count) + "\"" + "}"
print(printout.replace('\n', ' '))