New gwas track

pygenometracks/readGwas.py

@@ -0,0 +1,80 @@
+# -*- coding: utf-8 -*-
+import collections
+
+from .utilities import InputError, to_string
+
+
+class ReadGwas(object):
+    """
+    Reads a GWAS file. The expected fields are:
+    chromosome, position, name, and pvalue.
+
+    Example:
+    gwas = ReadGwas(open("file.gwas", 'r'))
+    for record in gwas:
+        print(record.chromosome, record.position, record.pvalue)
+    """
+
+    def __init__(self, file_handle, has_header=False):
+        """
+        :param file_handle: file handle
+        """
+        self.file_handle = file_handle
+        self.line_number = 0
+
+        # Define the fields for GWAS
+        self.fields = ['chromosome', 'position', 'name', 'pvalue']
+        self.GwasRecord = collections.namedtuple('GwasRecord', self.fields)
+
+        # Skip the header line if present
+        if has_header:
+            next(self.file_handle)
+            self.line_number += 1
+
+    def __iter__(self):
+        return self
+
+    def get_no_comment_line(self):
+        """
+        Skips comment lines starting with '#' or empty lines.
+        :return: a valid line
+        """
+        line = next(self.file_handle)
+        line = to_string(line)
+        if line.startswith("#") or line.strip() == '':
+            line = self.get_no_comment_line()
+
+        self.line_number += 1
+        return line
+
+    def __next__(self):
+        """
+        :return: GwasRecord object
+        """
+        line = self.get_no_comment_line()
+        return self.get_gwas_record(line)
+
+    def get_gwas_record(self, gwas_line):
+        """
+        Processes each line from a GWAS file and returns a namedtuple object.
+
+        :param gwas_line: a single line from the GWAS file
+        :return: GwasRecord object
+        """
+        line_data = gwas_line.strip()
+        line_data = to_string(line_data)
+        line_data = line_data.split("\t")
+
+        if len(line_data) < 4:
+            raise InputError(f"Line {self.line_number} does not have 4 fields: {gwas_line}."
+                             f"We expect at least 4 field, corresponding to: chromosome, position, name, pvalue.")
+
+        try:
+            chromosome = line_data[0]
+            position = int(line_data[1])
+            name = line_data[2]
+            pvalue = float(line_data[3])
+        except ValueError as e:
+            raise InputError(f"Error parsing line {self.line_number}: {gwas_line}\n{e}")
+
+        return self.GwasRecord(chromosome, position, name, pvalue)

pygenometracks/tests/generateAllOutput.sh

@@ -181,6 +181,10 @@ pgt --tracks ./pygenometracks/tests/test_data/browser_tracks_hic_small_test_squa
 pgt --tracks ./pygenometracks/tests/test_data/browser_tracks_hic_inbetween.ini --BED ./pygenometracks/tests/test_data/chrY_regions.bed --trackLabelFraction 0.23 --width 38 --dpi 130 -o ./pygenometracks/tests/test_data/master_hic_inbetween.png
 pgt --tracks ./pygenometracks/tests/test_data/mcool_hic_matrix_square.ini --region X:2500000-3500000 --trackLabelFraction 0.23 --width 38 --dpi 130 -o ./pygenometracks/tests/test_data/master_mcool_hic_matrix_square.png
 
+# gwas
+pgt --tracks ./pygenometracks/tests/test_data/gwas.ini --region X:3000000-3200000 --trackLabelFraction 0.2 --dpi 130 -o ./pygenometracks/tests/test_data/master_gwas.png
+pgt --tracks ./pygenometracks/tests/test_data/gwas.ini --region chrY:3000000-3200000 --trackLabelFraction 0.2 --dpi 130 -o ./pygenometracks/tests/test_data/master_gwas_chrY.png
+
 # demo
 pgt --tracks ./pygenometracks/tests/test_data/demo2.ini --region chrX:3320000-3370000 -o ./pygenometracks/tests/test_data/demo2.png
 pgt --tracks ./pygenometracks/tests/test_data/demo.ini --region chrX:3000000-3500000 -o ./pygenometracks/tests/test_data/demo.png

pygenometracks/tests/test_data/gwas.ini

@@ -0,0 +1,29 @@
+
+[gwas]
+file = gwas_1.gwas
+height = 4
+title = test_1 default values
+
+[spacer]
+
+[gwas_2]
+file = gwas_2.gwas
+file_has_header = True
+height = 2
+title = test_2 file_has_header = true color = #50E3C2 border_color = red line_width = 2 marker_size = 90 show_data_range = false
+color = #50E3C2
+border_color = red
+line_width = 2
+marker_size = 90
+show_data_range = false
+
+[spacer]
+
+[gwas_2]
+file = gwas_1.gwas
+height = 4
+title = test_1 default values min_value = 0 max_value = 15
+min_value = 0
+max_value = 15
+
+[x-axis]

pygenometracks/tests/test_data/gwas_1.gwas

@@ -0,0 +1,12 @@
+X	3002145	rs7823451	4.7e-08
+X	3005892	rs6610293	1.3e-06
+X	3008734	rs9034128	9.2e-05
+X	3031022	rs4420981	6.1e-10
+X	3042567	rs2093847	2.8e-07
+X	3056789	rs7734012	5.5e-09
+X	3072341	rs1182736	3.9e-06
+X	3079880	rs5502918	7.4e-05
+X	3112450	rs3847291	1.8e-08
+X	3145670	rs6029184	9.9e-07
+X	3180230	rs2938471	3.2e-10
+X	3197650	rs8840123	6.3e-08

pygenometracks/tests/test_data/gwas_2.gwas

@@ -0,0 +1,13 @@
+CHR	BP	SNP	P	BETA	SE	MAF
+X	3001200	rs100001	2.5e-09	0.62	0.09	0.21
+X	3003450	rs100002	7.8e-08	0.55	0.10	0.23
+X	3007890	rs100003	1.1e-06	0.41	0.08	0.27
+X	3009990	rs100004	4.3e-05	0.33	0.08	0.30
+X	3032100	rs100005	9.2e-10	0.71	0.11	0.18
+X	3045600	rs100006	3.4e-08	0.59	0.10	0.20
+X	3060200	rs100007	6.7e-07	0.48	0.09	0.25
+X	3078500	rs100008	1.9e-05	0.36	0.08	0.29
+X	3100400	rs100009	8.8e-06	-0.40	0.09	0.26
+X	3150000	rs100010	2.2e-07	-0.52	0.10	0.22
+X	3195000	rs100011	1.5e-09	-0.68	0.11	0.19
+X	3199800	rs100012	5.0e-08	-0.60	0.10	0.21

pygenometracks/tests/test_gwasTrack.py

@@ -0,0 +1,119 @@
+import os.path
+from tempfile import NamedTemporaryFile
+
+import matplotlib as mpl
+from get_matplotlib_CI_version import get_CI_mpl_version
+from matplotlib.testing.compare import compare_images
+
+import pygenometracks.plotTracks
+
+mpl.use('agg')
+
+ROOT = os.path.join(os.path.dirname(os.path.abspath(__file__)),
+                    "test_data")
+
+tracks = """
+[gwas]
+file = gwas_1.gwas
+height = 4
+title = test_1 default values
+
+[spacer]
+
+[gwas_2]
+file = gwas_2.gwas
+file_has_header = True
+height = 2
+title = test_2 file_has_header = true color = #50E3C2 border_color = red line_width = 2 marker_size = 90 show_data_range = false
+color = #50E3C2
+border_color = red
+line_width = 2
+marker_size = 90
+show_data_range = false
+
+[spacer]
+
+[gwas_2]
+file = gwas_1.gwas
+height = 4
+title = test_1 default values min_value = 0 max_value = 15
+min_value = 0
+max_value = 15
+
+[x-axis]
+"""
+
+with open(os.path.join(ROOT, "gwas.ini"), 'w') as fh:
+    fh.write(tracks)
+
+tolerance = 13  # default matplotlib pixed difference tolerance
+default_mpl_version = get_CI_mpl_version()
+
+
+def test_gwas_track():
+
+    if mpl.__version__ != default_mpl_version:
+        my_tolerance = 26
+    else:
+        my_tolerance = tolerance
+
+    outfile = NamedTemporaryFile(suffix='.png', prefix='gwas_test_',
+                                 delete=False)
+    ini_file = os.path.join(ROOT, "gwas.ini")
+    region = "X:3000000-3200000"
+    expected_file = os.path.join(ROOT, 'master_gwas.png')
+    args = f"--tracks {ini_file} --region {region} " \
+           "--trackLabelFraction 0.2 --dpi 130 " \
+           f"--outFileName {outfile.name}".split()
+    pygenometracks.plotTracks.main(args)
+    res = compare_images(expected_file,
+                         outfile.name, my_tolerance)
+    assert res is None, res
+
+    os.remove(outfile.name)
+
+
+def test_gwas_track_chrX():
+
+    if mpl.__version__ != default_mpl_version:
+        my_tolerance = 15
+    else:
+        my_tolerance = tolerance
+
+    outfile = NamedTemporaryFile(suffix='.png', prefix='gwas_test_',
+                                 delete=False)
+    ini_file = os.path.join(ROOT, "gwas.ini")
+    region = "chrX:3000000-3200000"
+    expected_file = os.path.join(ROOT, 'master_gwas.png')
+    args = f"--tracks {ini_file} --region {region} " \
+           "--trackLabelFraction 0.2 --dpi 130 " \
+           f"--outFileName {outfile.name}".split()
+    pygenometracks.plotTracks.main(args)
+    res = compare_images(expected_file,
+                         outfile.name, my_tolerance + 14)  # 14 corresponds to the 'chr' on the x axis
+    assert res is None, res
+
+    os.remove(outfile.name)
+
+
+def test_gwas_track_chrY():
+
+    if mpl.__version__ != default_mpl_version:
+        my_tolerance = 26
+    else:
+        my_tolerance = tolerance
+
+    outfile = NamedTemporaryFile(suffix='.png', prefix='gwas_test_',
+                                 delete=False)
+    ini_file = os.path.join(ROOT, "gwas.ini")
+    region = "chrY:3000000-3200000"
+    expected_file = os.path.join(ROOT, 'master_gwas_chrY.png')
+    args = f"--tracks {ini_file} --region {region} " \
+           "--trackLabelFraction 0.2 --dpi 130 " \
+           f"--outFileName {outfile.name}".split()
+    pygenometracks.plotTracks.main(args)
+    res = compare_images(expected_file,
+                         outfile.name, my_tolerance)
+    assert res is None, res
+
+    os.remove(outfile.name)