diff --git a/pygenometracks/readGwas.py b/pygenometracks/readGwas.py new file mode 100644 index 00000000..ea8c77e7 --- /dev/null +++ b/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) diff --git a/pygenometracks/tests/generateAllOutput.sh b/pygenometracks/tests/generateAllOutput.sh index 6e5055a3..1eee9419 100644 --- a/pygenometracks/tests/generateAllOutput.sh +++ b/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 diff --git a/pygenometracks/tests/test_data/gwas.ini b/pygenometracks/tests/test_data/gwas.ini new file mode 100644 index 00000000..f7390f6b --- /dev/null +++ b/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] diff --git a/pygenometracks/tests/test_data/gwas_1.gwas b/pygenometracks/tests/test_data/gwas_1.gwas new file mode 100644 index 00000000..dae8e6b2 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/pygenometracks/tests/test_data/gwas_2.gwas b/pygenometracks/tests/test_data/gwas_2.gwas new file mode 100644 index 00000000..ce7eb4c6 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/pygenometracks/tests/test_data/master_gwas.png b/pygenometracks/tests/test_data/master_gwas.png new file mode 100644 index 00000000..6b229356 Binary files /dev/null and b/pygenometracks/tests/test_data/master_gwas.png differ diff --git a/pygenometracks/tests/test_data/master_gwas_chrY.png b/pygenometracks/tests/test_data/master_gwas_chrY.png new file mode 100644 index 00000000..96642596 Binary files /dev/null and b/pygenometracks/tests/test_data/master_gwas_chrY.png differ diff --git a/pygenometracks/tests/test_gwasTrack.py b/pygenometracks/tests/test_gwasTrack.py new file mode 100644 index 00000000..376c48fa --- /dev/null +++ b/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) diff --git a/pygenometracks/tracks/GwasTrack.py b/pygenometracks/tracks/GwasTrack.py new file mode 100644 index 00000000..07b29b5b --- /dev/null +++ b/pygenometracks/tracks/GwasTrack.py @@ -0,0 +1,153 @@ +import numpy as np +from intervaltree import Interval, IntervalTree +from tqdm import tqdm + +from ..readGwas import ReadGwas +from ..utilities import change_chrom_names, count_lines, opener +from .GenomeTrack import GenomeTrack + +DEFAULT_GWAS_COLOR = '#ff7f00' + + +class GwasTrack(GenomeTrack): + SUPPORTED_ENDINGS = ['.gwas', '.linear', '.logistic', '.assoc', '.qassoc'] # this is used by make_tracks_file to guess the type of track based on file name + TRACK_TYPE = 'gwas' + OPTIONS_TXT = GenomeTrack.OPTIONS_TXT + f""" +# File containing the data. We expect an IGV .gwas format file with the columns: CHR, BP, SNP and P. +# Optionally, extra annotation columns can be added. +file = +# Indicate if your file has a header: +file_has_header = false +# Each SNP will be plotted as a 'o' and you can control color/size etc... +# Inside color +#color = red +# Border color +#border_color = black +# Line width +#line_width = 0.5 +# Size +#marker_size = 45 +# set show_data_range to false to hide the text on the upper-left showing the data range +show_data_range = true +# the default for min_value and max_value is 'auto' which means that the scale will go +# roughly from the minimum value found in the region plotted to the maximum value found. +min_value = 0 +#max_value = auto +# Optional. If not given is guessed from the file ending. +file_type = {TRACK_TYPE} + """ + + DEFAULTS_PROPERTIES = {'max_value': None, + 'min_value': None, + 'show_data_range': True, + 'orientation': None, + 'color': DEFAULT_GWAS_COLOR, + 'border_color': 'black', + 'line_width': 0.5, + 'marker_size': 45, + 'file_has_header': False} + + NECESSARY_PROPERTIES = ['file'] + SYNONYMOUS_PROPERTIES = {'max_value': {'auto': None}, + 'min_value': {'auto': None}} + POSSIBLE_PROPERTIES = {} + BOOLEAN_PROPERTIES = ['file_has_header', 'show_data_range'] + STRING_PROPERTIES = ['title', 'file_type', 'file', 'color', 'border_color'] + FLOAT_PROPERTIES = {'max_value': [- np.inf, np.inf], + 'min_value': [- np.inf, np.inf], + 'height': [0, np.inf], + 'marker_size': [0, np.inf], + 'line_width': [0, np.inf]} + INTEGER_PROPERTIES = {} + + def __init__(self, *args, **kwarg): + super(GwasTrack, self).__init__(*args, **kwarg) + self.interval_tree = self.process_gwas(self.properties['region']) + + def process_gwas(self, plot_regions=None): + """Read the gwas file and store values in a IntervalTree + + :param list plot_regions: list of plotted regions (like [(chrom1, start1, end1), (chrom2, start2, end2)]), defaults to None + :return None + """ + + total_length = count_lines(opener(self.properties['file']), + asBed=True) + if self.properties['file_has_header']: + total_length -= 1 + gwas_file_h = ReadGwas(opener(self.properties['file']), + has_header=self.properties['file_has_header']) + + valid_intervals = 0 + interval_tree = {} + + if plot_regions is not None: + chroms_to_plot = set([v[0] for v in plot_regions]) + else: + chroms_to_plot = None + + for record in tqdm(gwas_file_h, total=total_length): + + if plot_regions is not None and record.chromosome not in chroms_to_plot: + continue + + if record.chromosome not in interval_tree: + interval_tree[record.chromosome] = IntervalTree() + + interval_tree[record.chromosome].add(Interval(record.position, + record.position + 1, record)) + valid_intervals += 1 + + try: + gwas_file_h.file_handle.close() + except AttributeError: + pass + + if valid_intervals == 0: + self.log.warning("No valid intervals were found in file " + f"{self.properties['file']} for regions" + f"{plot_regions}.\n") + + return interval_tree + + def plot(self, ax, chrom_region, start_region, end_region): + """ + Plot a scatter plot for the GWAS data. + The p-values are transformed as -log10(pvalue), so the y-axis will show the exponents of the p-values. + + :param ax: matplotlib axis + :param chrom_region: chromosome name + :param start_region: start position of the region + :param end_region: end position of the region + :return: None + """ + if chrom_region not in self.interval_tree.keys(): + chrom_region_before = chrom_region + chrom_region = change_chrom_names(chrom_region) + if chrom_region not in self.interval_tree.keys(): + self.log.warning("*Warning*\nNo interval was found when " + "overlapping with both " + f"{chrom_region_before}:{start_region}-{end_region}" + f" and {chrom_region}:{start_region}-{end_region}" + " inside the gwas file. " + "This will generate an empty track!!\n") + self.adjust_ylim(ax) + return + + gwas_overlap = \ + self.interval_tree[chrom_region][start_region:end_region] + + # Fill in the position and pvalues lists with data from the GWAS file + position = [region.begin for region in gwas_overlap] + # Notice the -log10 transformation + y_values = [-np.log10(region.data.pvalue) if region.data.pvalue > 0 else 0 + for region in gwas_overlap] + + # Plot the scatterplot + ax.scatter(position, y_values, + s=self.properties['marker_size'], + color=self.properties['color'], marker='o', + edgecolors=self.properties['border_color'], + linewidths=self.properties['line_width']) + + self.adjust_ylim(ax)