manhattan.py

   1 # Copyright (C) The Lightning Authors. All rights reserved.
   2 #
   3 # SPDX-License-Identifier: AGPL-3.0
   4
   5 import csv
   6 import os
   7 import sys
   8
   9 import matplotlib
  10 import numpy
  11 import pandas
  12 import qmplot
  13
  14 (_,
  15  input_path,
  16  output_path,
  17  ) = sys.argv
  18 columns = numpy.load(os.path.join(input_path, 'onehot-columns.npy'))
  19
  20 # pvalue maps tag# => [pvalue1, pvalue2, ...] (one het p-value and one hom p-value for each tile variant)
  21 pvalue = {}
  22 for i in range(columns.shape[1]):
  23     tag = columns[0,i]
  24     x = pvalue.get(tag, [])
  25     x.append(pow(10, -columns[4,i] / 1000000))
  26     pvalue[tag] = x
  27
  28 # tilepos maps tag# => (chromosome, position)
  29 tilepos = {}
  30 for dirent in os.scandir(input_path):
  31     if dirent.name.endswith('.annotations.csv'):
  32         with open(dirent, 'rt', newline='') as annotations:
  33             for annotation in csv.reader(annotations):
  34                 # 500000,0,2,=,chr1,160793649,,,
  35                 if annotation[3] == "=":
  36                     tilepos[int(annotation[0])] = (annotation[4], int(annotation[5]))
  37
  38 series = {"#CHROM": [], "POS": [], "P": []}
  39 for tag, chrpos in sorted(tilepos.items(), key=lambda item: (item[1][0][-1] > '9', item[1][0].lstrip('chr').zfill(2), item[1][1])):
  40     for p in pvalue.get(tag, []):
  41         series["#CHROM"].append(chrpos[0])
  42         series["POS"].append(chrpos[1])
  43         series["P"].append(p)
  44
  45 qmplot.manhattanplot(data=pandas.DataFrame(series),
  46                      suggestiveline=2e-10,  # Turn off suggestiveline
  47                      genomewideline=2e-11,  # Turn off genomewidel
  48                      sign_line_cols=["#D62728", "#2CA02C"],
  49                      marker=".",
  50                      alpha = 0.6,
  51                      hline_kws={"linestyle": "--", "lw": 1.3},
  52                      title="Tile Variant Manhattan Plot",
  53                      # xtick_label_set=xtick,
  54                      xlabel="Chromosome",
  55                      ylabel=r"$-log_{10}{(P)}$",
  56                      xticklabel_kws={"rotation": "vertical"})
  57 matplotlib.pyplot.savefig(output_path, bbox_inches="tight")