|
| 1 | +#!/usr/bin/python3 |
| 2 | + |
| 3 | +import math |
| 4 | +import sys |
| 5 | +import argparse |
| 6 | + |
| 7 | +from random import seed |
| 8 | +from random import random |
| 9 | + |
| 10 | +from reportlab.lib import colors |
| 11 | +from reportlab.lib.units import cm |
| 12 | +from Bio.Graphics import GenomeDiagram |
| 13 | +from Bio.SeqFeature import SeqFeature, FeatureLocation |
| 14 | + |
| 15 | +colors_list = [colors.red, colors.blue, colors.yellow, colors.brown, colors.green, colors.darkmagenta, colors.greenyellow, colors.lavender, colors.purple, colors.navy, colors.orange, colors.paleturquoise, colors.powderblue, colors.rosybrown, colors.turquoise, colors.black] |
| 16 | + |
| 17 | +parser = argparse.ArgumentParser(description='Create an image with the distribution of a list of genes. Lines in the input TSV file must have four columns: (1) the group to wich the gene belongs to (each group is drawn in a different color), (2) the name of the gene, (3) the start coordinate, and (4) the end coordinate.') |
| 18 | + |
| 19 | +parser.add_argument('INPUT_FILE', help='input TSV file with gene names and coordinates') |
| 20 | +parser.add_argument('-o', '--output', default='image', help='output image name (without file extension)') |
| 21 | +parser.add_argument('-f', '--format', default='pdf', help='output image format: PS, PDF, SVG, JPG, BMP, GIF, PNG, TIFF or TIFF (default is PDF)') |
| 22 | +parser.add_argument('-s', '--start', help='start coordinate to draw (if not provided, it is automatically determined from the input data)') |
| 23 | +parser.add_argument('-e', '--end', help='end coordinate to draw (if not provided, it is automatically determined from the input data)') |
| 24 | +parser.add_argument('-p','--preserve-strand', help='use this flag to preserve gene strands', action='store_true') |
| 25 | +parser.add_argument('-b','--breaks', default='20', help='number of axis breaks to represent the scale') |
| 26 | +parser.add_argument('-r','--random-seed', default='1', help='random seed for color generation') |
| 27 | + |
| 28 | +arg = parser.parse_args() |
| 29 | + |
| 30 | +fh = open(arg.INPUT_FILE) |
| 31 | +lines = fh.readlines() |
| 32 | + |
| 33 | +series_features = {} |
| 34 | +series_indexes = {} |
| 35 | + |
| 36 | +series_index = 0 |
| 37 | +for line in lines: |
| 38 | + if not ('#' in line): |
| 39 | + temp_line = line.upper().replace('\n', '') |
| 40 | + sp_line = temp_line.split(sep='\t') |
| 41 | + series = sp_line[0] |
| 42 | + gene_name = sp_line[1] |
| 43 | + |
| 44 | + if int(sp_line[2]) > int(sp_line[3]): |
| 45 | + start = int(sp_line[3]) |
| 46 | + end = sp_line [2] |
| 47 | + if arg.preserve_strand: |
| 48 | + strand = -1 |
| 49 | + else: |
| 50 | + strand = 1 |
| 51 | + else: |
| 52 | + start = int(sp_line[2]) |
| 53 | + end = int(sp_line[3]) |
| 54 | + strand = 1 |
| 55 | + |
| 56 | + if series in series_features: |
| 57 | + series_features[series].append((gene_name, strand, start, end)) |
| 58 | + else: |
| 59 | + series_features[series] = [(gene_name, strand, start, end)] |
| 60 | + series_indexes[series] = series_index |
| 61 | + series_index = series_index + 1 |
| 62 | + |
| 63 | +start = sys.maxsize |
| 64 | +end = -1 |
| 65 | + |
| 66 | +gdd = GenomeDiagram.Diagram("diagram", tracklines = False, y = 0.4) |
| 67 | +gd_track_for_features = gdd.new_track(1, scale = True, height = 1, scale_smallticks = 0) |
| 68 | +gds_features = gd_track_for_features.new_set() |
| 69 | + |
| 70 | +seed(int(arg.random_seed)) |
| 71 | + |
| 72 | +for series in series_features.keys(): |
| 73 | + if series_indexes[series] < len(colors_list): |
| 74 | + current_color = colors_list[series_indexes[series]] |
| 75 | + else: |
| 76 | + current_color = colors.Color(random(), random(), random()) |
| 77 | + |
| 78 | + for i in range(0, len(series_features[series])): |
| 79 | + current_feature = series_features[series][i] |
| 80 | + feature = SeqFeature(FeatureLocation(int(current_feature[2]), int(current_feature[3])), strand=current_feature[1]) |
| 81 | + gds_features.add_feature(feature, name="{}".format(current_feature[0]), label=True, color=current_color) |
| 82 | + |
| 83 | + if int(current_feature[2]) < start: |
| 84 | + start = int(current_feature[2]) |
| 85 | + if int(current_feature[3]) > end: |
| 86 | + end = int(current_feature[3]) |
| 87 | + |
| 88 | + |
| 89 | +if not arg.start == None: |
| 90 | + start = int(arg.start) |
| 91 | +else: |
| 92 | + start = start - 1000 |
| 93 | + |
| 94 | +if not arg.end == None: |
| 95 | + end = int(arg.end) |
| 96 | +else: |
| 97 | + end = end + 1000 |
| 98 | + |
| 99 | +if start > end: |
| 100 | + sys.exit("start must be less than end") |
| 101 | + |
| 102 | +gd_track_for_features.start = start |
| 103 | +gd_track_for_features.end = end |
| 104 | +gd_track_for_features.scale_largetick_interval = (end-start) / int(arg.breaks) |
| 105 | + |
| 106 | +gdd.draw(format="linear", fragments=1, start=start, end=end) |
| 107 | +gdd.write(arg.output+'.'+arg.format, arg.format) |
0 commit comments