From b0e57c72714822ad0bdf30218b703ad68e19f921 Mon Sep 17 00:00:00 2001 From: lucas_miranda Date: Mon, 22 Mar 2021 17:57:09 +0100 Subject: [PATCH 1/6] Added color and ax options to manhattan_plot in util.py --- fastlmm/util/util.py | 104 +++++++++++++++++++++++++++++++++---------- 1 file changed, 80 insertions(+), 24 deletions(-) diff --git a/fastlmm/util/util.py b/fastlmm/util/util.py index a492bd0f..fac32d7c 100644 --- a/fastlmm/util/util.py +++ b/fastlmm/util/util.py @@ -434,11 +434,22 @@ def _color_list(chr_list,rle): result = [index_to_color[chr_to_index[chr]%len(index_to_color)] for chr in chr_list] return result -def manhattan_plot(chr_pos_pvalue_array,pvalue_line=None,plot_threshold=1.0,vline_significant=False,marker="o", chromosome_starts=None, xaxis_unit_bp=True, alpha=0.5): +def manhattan_plot( + chr_pos_pvalue_array, + pvalue_line=None, + plot_threshold=1.0, + vline_significant=False, + marker="o", + chromosome_starts=None, + xaxis_unit_bp=True, + alpha=0.5, + color=None, + ax=None, +): """ Function to create a Manhattan plot. See http://en.wikipedia.org/wiki/Manhattan_plot. - :param chr_pos_pvalue_array: an n x 3 numpy array. The three columns are the chrom number + :param chr_pos_pvalue_array: an n x 3 numpy array. The three columns are the chrom number (as a number), the position, and pvalue. :type chr_pos_pvalue_array: numpy array :param pvalue_line: (Default: None). If given, draws a line at that PValue. @@ -453,7 +464,7 @@ def manhattan_plot(chr_pos_pvalue_array,pvalue_line=None,plot_threshold=1.0,vlin :param chromosome_starts: chromosome, cumulative start position, cumulative stop position cumulative chromosome starts, for plotting. If None (default), this is estimated from data :type chromosome_starts: [Nchrom x 3] ndarray - :param xaxis_unit_bp: If true, plot cumulative position in basepair units on x axis. If False, only + :param xaxis_unit_bp: If true, plot cumulative position in basepair units on x axis. If False, only use rank of SNP positions. (default: True) :type xaxis_unit_bp: Boolean :param alpha: alpha (opaqueness) for P-value markers in scatterplot (default 0.5) @@ -461,6 +472,12 @@ def manhattan_plot(chr_pos_pvalue_array,pvalue_line=None,plot_threshold=1.0,vlin :rtype: chromosome_starts [Nchrom x 3] ndarray: chromosome, cumulative start position, cumulative stop position cumulative chromosome starts used in plotting. + :param alpha: alpha (opaqueness) for P-value markers in scatterplot (default 0.5) + :type alpha: number + :param color: color to use for the final scatterplot. If None (default) chromosomes are alternately colored + :type color: string + :param ax: ax where to plot the figure. If None (default) a fresh canvas is used + :type ax: matplotlib.ax :Example: @@ -481,42 +498,81 @@ def manhattan_plot(chr_pos_pvalue_array,pvalue_line=None,plot_threshold=1.0,vlin # create a copy of the data and sort it by chrom and then position array = np.array(chr_pos_pvalue_array) if plot_threshold: - array = array[array[:,2]<=plot_threshold] + array = array[array[:, 2] <= plot_threshold] else: plot_threshold = 1.0 - array=array[np.argsort(array[:,1]),:] #sort by ChrPos - array=array[np.argsort(array[:,0],kind='mergesort'),:] #Finally, sort by Chr (but keep ChrPos in case of ties) - rle = list(_run_length_encode(array[:,0])) - - if xaxis_unit_bp: #compute and use cumulative basepair positions for x-axis + array = array[np.argsort(array[:, 1]), :] # sort by ChrPos + array = array[ + np.argsort(array[:, 0], kind="mergesort"), : + ] # Finally, sort by Chr (but keep ChrPos in case of ties) + rle = list(_run_length_encode(array[:, 0])) + + if xaxis_unit_bp: # compute and use cumulative basepair positions for x-axis if chromosome_starts is None: chromosome_starts = _compute_x_positions_chrom(array) chr_pos_list = _compute_x_positions_snps(array, chromosome_starts) - plt.xlim([0,chromosome_starts[-1,2]+1]) - plt.xticks(chromosome_starts[:,1:3].mean(1),chromosome_starts[:,0]) - else: #use rank indices for x-axis + plt.xlim([0, chromosome_starts[-1, 2] + 1]) + plt.xticks(chromosome_starts[:, 1:3].mean(1), chromosome_starts[:, 0]) + else: # use rank indices for x-axis chr_pos_list = np.arange(array.shape[0]) - xTickMarks = [str(int(item)) for item,count in rle] - plt.xlim([0,array.shape[0]]) + xTickMarks = [str(int(item)) for item, count in rle] + plt.xlim([0, array.shape[0]]) plt.xticks(list(_rel_to_midpoint(rle)), xTickMarks) - y = -np.log10(array[:,2]) + y = -np.log10(array[:, 2]) max_y = y.max() - if pvalue_line and vline_significant: #mark significant associations (ones that pass the pvalue_line) by a red vertical line: - idx_significant = array[:,2] Date: Mon, 22 Mar 2021 19:26:41 +0100 Subject: [PATCH 2/6] Improved ax options to manhattan_plot in util.py to avoid redundancy --- fastlmm/util/util.py | 61 +++++++++++++++++--------------------------- 1 file changed, 23 insertions(+), 38 deletions(-) diff --git a/fastlmm/util/util.py b/fastlmm/util/util.py index fac32d7c..84438fd8 100644 --- a/fastlmm/util/util.py +++ b/fastlmm/util/util.py @@ -495,6 +495,11 @@ def manhattan_plot( """ import matplotlib.pyplot as plt + # Create an empty canvas if none is provided + # this way the function can only work on 'ax' + if ax is None: + _, ax = plt.subplots() + # create a copy of the data and sort it by chrom and then position array = np.array(chr_pos_pvalue_array) if plot_threshold: @@ -511,12 +516,12 @@ def manhattan_plot( if chromosome_starts is None: chromosome_starts = _compute_x_positions_chrom(array) chr_pos_list = _compute_x_positions_snps(array, chromosome_starts) - plt.xlim([0, chromosome_starts[-1, 2] + 1]) + ax.set_xlim([0, chromosome_starts[-1, 2] + 1]) plt.xticks(chromosome_starts[:, 1:3].mean(1), chromosome_starts[:, 0]) else: # use rank indices for x-axis chr_pos_list = np.arange(array.shape[0]) xTickMarks = [str(int(item)) for item, count in rle] - plt.xlim([0, array.shape[0]]) + ax.set_xlim([0, array.shape[0]]) plt.xticks(list(_rel_to_midpoint(rle)), xTickMarks) y = -np.log10(array[:, 2]) max_y = y.max() @@ -529,7 +534,7 @@ def manhattan_plot( y_significant = y[idx_significant] chr_pos_list_significant = chr_pos_list[idx_significant] for i in range(len(chr_pos_list_significant)): - plt.axvline( + ax.axvline( x=chr_pos_list_significant[i], ymin=0.0, ymax=y_significant[i], @@ -537,41 +542,21 @@ def manhattan_plot( alpha=0.8, ) - if ax is None: - - plt.scatter( - chr_pos_list, - y, - marker=marker, - c=(_color_list(array[:, 0], rle) if color is None else color), - edgecolor="none", - s=y / max_y * 20 + 0.5, - alpha=alpha, - ) - plt.xlabel("chromosome") - plt.ylabel("-log10(P value)") - - if pvalue_line: - plt.axhline(-np.log10(pvalue_line), linestyle="--", color="gray") - plt.ylim([-np.log10(plot_threshold), None]) - - else: - - ax.scatter( - chr_pos_list, - y, - marker=marker, - c=(_color_list(array[:, 0], rle) if color is None else color), - edgecolor="none", - s=y / max_y * 20 + 0.5, - alpha=alpha, - ) - ax.set_xlabel("chromosome") - ax.set_ylabel("-log10(P value)") - - if pvalue_line: - ax.axhline(-np.log10(pvalue_line), linestyle="--", color="gray") - ax.set_ylim([-np.log10(plot_threshold), None]) + ax.scatter( + chr_pos_list, + y, + marker=marker, + c=(_color_list(array[:, 0], rle) if color is None else color), + edgecolor="none", + s=y / max_y * 20 + 0.5, + alpha=alpha, + ) + ax.set_xlabel("chromosome") + ax.set_ylabel("-log10(P value)") + + if pvalue_line: + ax.axhline(-np.log10(pvalue_line), linestyle="--", color="gray") + ax.set_ylim([-np.log10(plot_threshold), None]) return chromosome_starts From 021a4c0456cd78ead51a0577d8bf6afdea138723 Mon Sep 17 00:00:00 2001 From: lucas_miranda Date: Mon, 22 Mar 2021 22:03:16 +0100 Subject: [PATCH 3/6] Improved ax options to manhattan_plot in util.py to avoid redundancy --- fastlmm/util/util.py | 24 +++++++++++++++++++++--- 1 file changed, 21 insertions(+), 3 deletions(-) diff --git a/fastlmm/util/util.py b/fastlmm/util/util.py index 84438fd8..27c67d92 100644 --- a/fastlmm/util/util.py +++ b/fastlmm/util/util.py @@ -498,7 +498,7 @@ def manhattan_plot( # Create an empty canvas if none is provided # this way the function can only work on 'ax' if ax is None: - _, ax = plt.subplots() + ax = plt.axes(label="manhattan") # create a copy of the data and sort it by chrom and then position array = np.array(chr_pos_pvalue_array) @@ -517,12 +517,30 @@ def manhattan_plot( chromosome_starts = _compute_x_positions_chrom(array) chr_pos_list = _compute_x_positions_snps(array, chromosome_starts) ax.set_xlim([0, chromosome_starts[-1, 2] + 1]) - plt.xticks(chromosome_starts[:, 1:3].mean(1), chromosome_starts[:, 0]) + + # If canvas already exists, add new ticks to it and keep only unique values + if ax.get_xticks()[0] == 0.0: + ax.set_xticks([]) + + xticks = np.concatenate([ax.get_xticks(), chromosome_starts[:, 1:3].mean(1)]) + old_ticklabels = [float(i.get_text()) for i in ax.get_xticklabels()] + + unique_chr, unique_idx = np.unique( + np.concatenate([old_ticklabels, chromosome_starts[:, 0]]), return_index=True + ) + + xticks = xticks[unique_idx] + xticklabels = unique_chr + + ax.set_xticks(np.sort(xticks)) + ax.set_xticklabels(np.sort(xticklabels)) + else: # use rank indices for x-axis chr_pos_list = np.arange(array.shape[0]) xTickMarks = [str(int(item)) for item, count in rle] ax.set_xlim([0, array.shape[0]]) - plt.xticks(list(_rel_to_midpoint(rle)), xTickMarks) + ax.set_xticks(list(_rel_to_midpoint(rle))) + ax.set_xticklabels(xTickMarks) y = -np.log10(array[:, 2]) max_y = y.max() From 3fd43e13fc5cd390365999177355be57598f802f Mon Sep 17 00:00:00 2001 From: lucas_miranda Date: Mon, 22 Mar 2021 22:16:52 +0100 Subject: [PATCH 4/6] Improved ax options to manhattan_plot in util.py to avoid redundancy --- fastlmm/util/util.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/fastlmm/util/util.py b/fastlmm/util/util.py index 27c67d92..7b863b81 100644 --- a/fastlmm/util/util.py +++ b/fastlmm/util/util.py @@ -532,6 +532,9 @@ def manhattan_plot( xticks = xticks[unique_idx] xticklabels = unique_chr + print(unique_chr) + print(unique_idx) + ax.set_xticks(np.sort(xticks)) ax.set_xticklabels(np.sort(xticklabels)) @@ -541,6 +544,7 @@ def manhattan_plot( ax.set_xlim([0, array.shape[0]]) ax.set_xticks(list(_rel_to_midpoint(rle))) ax.set_xticklabels(xTickMarks) + y = -np.log10(array[:, 2]) max_y = y.max() From 44f07381c714a5afe91fc3e9680fbea9d56f13f6 Mon Sep 17 00:00:00 2001 From: lucas_miranda Date: Mon, 22 Mar 2021 22:39:10 +0100 Subject: [PATCH 5/6] Improved ax options to manhattan_plot in util.py to avoid redundancy --- fastlmm/util/util.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/fastlmm/util/util.py b/fastlmm/util/util.py index 7b863b81..d4a6dcea 100644 --- a/fastlmm/util/util.py +++ b/fastlmm/util/util.py @@ -469,15 +469,14 @@ def manhattan_plot( :type xaxis_unit_bp: Boolean :param alpha: alpha (opaqueness) for P-value markers in scatterplot (default 0.5) :type alpha: number - - :rtype: chromosome_starts [Nchrom x 3] ndarray: chromosome, cumulative start position, cumulative stop position - cumulative chromosome starts used in plotting. :param alpha: alpha (opaqueness) for P-value markers in scatterplot (default 0.5) :type alpha: number :param color: color to use for the final scatterplot. If None (default) chromosomes are alternately colored :type color: string :param ax: ax where to plot the figure. If None (default) a fresh canvas is used :type ax: matplotlib.ax + :rtype: chromosome_starts [Nchrom x 3] ndarray: chromosome, cumulative start position, cumulative stop position + cumulative chromosome starts used in plotting. :Example: @@ -532,9 +531,6 @@ def manhattan_plot( xticks = xticks[unique_idx] xticklabels = unique_chr - print(unique_chr) - print(unique_idx) - ax.set_xticks(np.sort(xticks)) ax.set_xticklabels(np.sort(xticklabels)) From e8b6c12c964cd960fddc40dc7538e89f0a5d472f Mon Sep 17 00:00:00 2001 From: lucas_miranda Date: Tue, 23 Mar 2021 13:40:22 +0100 Subject: [PATCH 6/6] Simplyfied code --- fastlmm/util/util.py | 19 ++----------------- 1 file changed, 2 insertions(+), 17 deletions(-) diff --git a/fastlmm/util/util.py b/fastlmm/util/util.py index d4a6dcea..a07bf847 100644 --- a/fastlmm/util/util.py +++ b/fastlmm/util/util.py @@ -516,23 +516,8 @@ def manhattan_plot( chromosome_starts = _compute_x_positions_chrom(array) chr_pos_list = _compute_x_positions_snps(array, chromosome_starts) ax.set_xlim([0, chromosome_starts[-1, 2] + 1]) - - # If canvas already exists, add new ticks to it and keep only unique values - if ax.get_xticks()[0] == 0.0: - ax.set_xticks([]) - - xticks = np.concatenate([ax.get_xticks(), chromosome_starts[:, 1:3].mean(1)]) - old_ticklabels = [float(i.get_text()) for i in ax.get_xticklabels()] - - unique_chr, unique_idx = np.unique( - np.concatenate([old_ticklabels, chromosome_starts[:, 0]]), return_index=True - ) - - xticks = xticks[unique_idx] - xticklabels = unique_chr - - ax.set_xticks(np.sort(xticks)) - ax.set_xticklabels(np.sort(xticklabels)) + ax.set_xticks(chromosome_starts[:, 1:3].mean(1)) + ax.set_xticklabels(chromosome_starts[:, 0]) else: # use rank indices for x-axis chr_pos_list = np.arange(array.shape[0])