Add colored Manhattan plots and fdr threshold

Users can now specify to generate colored Manhattan plots and also specify which fdr threshold should be plotted.

Changelog

@@ -1,3 +1,7 @@
+1.7.2
+--
+* Enhancement: Users can specifiy if the chromosome should be colored in different colors and which fdr threshold to include in the Manhattan plots
+
 1.7.1
 --
 * Bugfix: loading a gwas result from csv will properly check if it contains pvalues or scores

pygwas/__init__.py

@@ -1,3 +1,3 @@
-__version__ = '1.7.1'
-__updated__ = "13.06.2017"
+__version__ = '1.7.2'
+__updated__ = "20.09.2017"
 __date__ = "20.8.2014"

pygwas/core/plotting.py

@@ -14,7 +14,9 @@
 
 SUPPORTED_FORMAT=('png','pdf')
 
-def plot_gwas_result(gwas_result,output_file,chrs=None,mac=15, marker_size=10):
+def plot_gwas_result(gwas_result,output_file,chrs=None,mac=15, marker_size=10, fdr = 'all', color_map=None):
+    if not color_map:
+        color_map = ['#4F94CD', '#36648B']
     chr_map = None
     if chrs is not None:
         chr_map = set(chrs)
@@ -34,7 +36,6 @@ def plot_gwas_result(gwas_result,output_file,chrs=None,mac=15, marker_size=10):
     data = _get_data(gwas_result)
     offset = 0
     markersize = marker_size
-    color_map = ['#4F94CD', '#36648B']
     plt.figure(figsize=(11, 3.8))
     plt.axes([0.045, 0.15, 0.99, 0.61])
 
@@ -71,12 +72,18 @@ def plot_gwas_result(gwas_result,output_file,chrs=None,mac=15, marker_size=10):
 
     score_range = max_score - min_score
     padding = 0.05*(score_range)
-    bonf_handle, = plt.plot([0, x_range], [bonferroni_threshold, bonferroni_threshold], color='r', linestyle="--", linewidth=1, alpha=.5)
-    if bh_thres is not None:
-        bh_handle, = plt.plot([0, x_range], [bh_thres, bh_thres], color='b', linestyle='--', linewidth=1, alpha=.5)
-        plt.figlegend((bonf_handle, bh_handle), ('Bonferroni', 'Benjamini Hochberg'), 'upper right')
-    else:
-        plt.figlegend((bonf_handle,), ('Bonferroni',), 'upper right')
+    fdr_labels = set()
+    fdr_handles = ()#
+    if fdr in ('all','bonferroni'):
+        handle, = plt.plot([0, x_range], [bonferroni_threshold, bonferroni_threshold], color='r', linestyle="--", linewidth=1, alpha=.5)
+        fdr_handles = fdr_handles + (handle,)
+        fdr_labels.add("Bonferroni")
+    if bh_thres is not None and fdr in ('all', 'benjamini_hochberg'):
+        handle, = plt.plot([0, x_range], [bh_thres, bh_thres], color='b', linestyle='--', linewidth=1, alpha=.5)
+        fdr_handles = fdr_handles + (handle,)
+        fdr_labels.add("Benjamini Hochberg")
+    if len(fdr_labels) > 0:
+        plt.figlegend(fdr_handles, fdr_labels, 'upper right')
     plt.axis([-x_range * 0.01, x_range * 1.01, min_score - padding, max_score + padding])
     ax = plt.gca()
     ax.xaxis.set_tick_params(direction='out')

pygwas/pygwas.py

@@ -102,6 +102,8 @@ def get_parser(program_license,program_version_message):
     plotter_parser.add_argument('-m','--macs',dest='macs',default=15,type=int,help='Minor Allele Count filter (Default: 15)')
     plotter_parser.add_argument('-s','--size',dest='marker_size',default=10,type=int,help='Size of the markers in the Manhattan plot (Default: 10)')
     plotter_parser.add_argument("-o",'--output',dest='output',required=True,help='The output image file')
+    plotter_parser.add_argument("-f",'--fdr',dest='fdr',default="all", help='The FDR threshold to plot. Default[all]', choices=["all", "bonferroni","benjamini_hochberg"])
+    plotter_parser.add_argument('--colored',dest='colored',action='store_true', help='Flag to use different colors for the chromosomes (Default:False)')
     plotter_parser.add_argument(dest="file", help="GWAS result file (.hdf5 or .csv)", metavar="FILE")
     plotter_parser.set_defaults(func=plot)
 
@@ -159,6 +161,7 @@ def main():
     program_build_date = str(__updated__)
     program_version_message = '%%(prog)s %s (%s)' % (program_version, program_build_date)
     program_shortdesc = "The main module for running Genome Wide Association studies"
+
     program_license = '''%s
 
   Created by Ümit Seren on %s.
@@ -268,7 +271,12 @@ def plot(args):
         gwas_result = result.load_from_hdf5(args['file'])
     else:
         gwas_result = result.load_from_csv(args['file'])
-    plotting.plot_gwas_result(gwas_result,args['output'],chrs,args['macs'],marker_size=marker_size)
+    fdr = args['fdr']
+    colored = args['colored']
+    color_map = None
+    if colored:
+        color_map = ['b', 'g', 'r', 'c', 'm']
+    plotting.plot_gwas_result(gwas_result,args['output'],chrs,args['macs'],marker_size=marker_size, fdr=fdr, color_map=color_map)
 
 
 def qq_plot(args):