plotting annotations

genda/eQTL/plotting.py

@@ -1,19 +1,20 @@
 """ Plotting functions for eQTLs
 """
+from collections import defaultdict
+import re
+
 import numpy as np
 import pandas as pd
 import matplotlib
 import statsmodels.api as sm
-
 import matplotlib.pyplot as plt
 from statsmodels.graphics import utils as smutils
 from statsmodels.graphics.regressionplots import abline_plot
 
-from genda.plotting import should_not_plot, add_gene_bounderies
+from genda.plotting import (should_not_plot, add_gene_bounderies,
+        snp_arrow)
 from genda import calculate_minor_allele_frequency, calculate_ld
 
-from IPython import embed
-
 
 def var_boxplot(ax, x, colors=None):
     """
@@ -185,8 +186,6 @@ def plot_eQTL(meQTL, gene_name, annotation, dosage, ax=None,
     snp_pv = adj_pv.iloc[iix[0]]
     color1 = calculate_ld(dosage_sub,
             snp)[dosage_sub.index].values
-    print(snp)
-    print(color1[0:5])
     if ax is None:
         ax_orig = False 
         fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(16, 6), 
@@ -204,13 +203,7 @@ def plot_eQTL(meQTL, gene_name, annotation, dosage, ax=None,
     ### Actual scatter #############################
     im = ax.scatter(pos, adj_pv, s=dosage_maf, c = color1)
     #:TODO make the arrow into a funciton
-    arrow_start = ylim - max(adj_pv/6.0)/2
-    arrow_length = max(adj_pv/6.0/2) - max(adj_pv/6.0/2)/8
-    ax.arrow(snpx, arrow_start, 0, - arrow_length, 
-            head_width=0.01, head_length=0.1, fc='k',
-            ec='k')
-    ax.text(snpx-0.05 , arrow_start + max(adj_pv/6.0)/5.0, 
-            snp, style='italic')
+    ax = snp_arrow(snpx, snp_pv, snp, ax)
     ax.set_ylabel(r'$-log_{10}$ eQTL p-value')
     ax.set_xlabel(r'Position (Mb)')
     if should_not_plot(gene_annot):

genda/plotting/__init__.py

@@ -1,65 +1,9 @@
-from collections import defaultdict
-import matplotlib.pyplot as plt
-import matplotlib.patches as patches
-import re
 #from scipy.stats import linregress, pearsonr
 #from matplotlib import figure
-import pandas as pd
-import numpy as np
-import pysam
 
 from matplotlib.path import Path
-from genda.formats.dosages import grab_gene_location
+from .annotation import (snp_arrow)
 from .plotting_utils import *
 
 
-def make_rectangle(start, end, y1, y2):
-    verts = [(start, y1),
-            (start, y2),
-            (end, y2),
-            (end, y1),
-            (start, y1),
-            ]
-    codes = [
-            Path.MOVETO,
-            Path.LINETO,
-            Path.LINETO,
-            Path.LINETO,
-            Path.CLOSEPOLY,
-            ]
-    return (Path(verts, codes))
 
-
-def create_path(gtf_iterator, gene, height=0.5, x_scale=1):
-    """ Create a path for a given transcript
-
-    :TODO change name
-    """
-    transcripts = defaultdict(list)
-    t_rxe = re.compile('transcript_id "(ENST[0-9]+)"')
-    t_add = 0
-    for i in gtf_iterator:
-        mverts = []
-        i = i.split("\t")
-        if (i[3] == gene or gene in i[9]) and i[7] == 'exon':
-            verts = [(int(i[1])/x_scale, t_add),
-                     (int(i[1])/x_scale, height + t_add),
-                     (int(i[2])/x_scale, height + t_add),
-                     (int(i[2])/x_scale, t_add),
-                     (int(i[1])/x_scale, t_add),
-                    ] 
-
-	    codes = [
-		    Path.MOVETO,
-		    Path.LINETO,
-		    Path.LINETO,
-		    Path.LINETO,
-		    Path.CLOSEPOLY,
-		    ]
-	    transcript = t_rxe.search(i[9]).group(1)
-	    if transcript not in transcripts.keys():
-	        t_add += 1
-	    else: pass
-	    transcripts[transcript].append(Path(verts, codes))
-        else: pass
-    return(transcripts)

genda/plotting/annotation.py

@@ -0,0 +1,36 @@
+
+
+def snp_arrow(pos, height, snpid, ax):
+    """ Draws an arrow pointing at a SNP in 
+    a Manhattan plot
+
+    pos - an object with pos and height
+    height - height
+    ax - matplotlib axes
+
+    """
+
+    rel_width = 25
+
+    p_width = ax.get_xlim()
+    p_height = ax.get_ylim()
+    '''
+    spacer = (p_height[1] - p_height[0])/rel_width
+    a_head_width = (p_width[1] - p_width[0])/rel_width
+    a_head_length = p_height[1] - p_height[0]/rel_width
+    arrow_length = p_height - height - spacer
+    ax.arrow(pos, p_height[1], 0, 
+            arrow_length, 
+            head_width = a_head_width,
+            head_length = a_head_length, 
+            fc = 'k')
+    '''
+
+    ax.annotate(snpid, xy=(pos, height),
+            xytext = (pos, p_height[1]),
+            arrowprops=dict(facecolor='black', shrink=0.05))
+
+
+    return ax
+
+

genda/plotting/plotting_utils.py

@@ -1,8 +1,8 @@
 """ Functions for adding to an axis.  Most functions here require passing in an
 axis object from matplotlib
 """
-
-#from collections import defaultdict
+from collections import defaultdict
+import re
 import matplotlib.patches as patches
 from matplotlib.path import Path
 import numpy as np
@@ -11,6 +11,41 @@
 from genda.formats.dosages import grab_gene_location
 
 
+def create_path(gtf_iterator, gene, height=0.5, x_scale=1):
+    """ Create a path for a given transcript
+
+    :TODO change name
+    """
+    transcripts = defaultdict(list)
+    t_rxe = re.compile('transcript_id "(ENST[0-9]+)"')
+    t_add = 0
+    for i in gtf_iterator:
+        mverts = []
+        i = i.split("\t")
+        if (i[3] == gene or gene in i[9]) and i[7] == 'exon':
+            verts = [(int(i[1])/x_scale, t_add),
+                     (int(i[1])/x_scale, height + t_add),
+                     (int(i[2])/x_scale, height + t_add),
+                     (int(i[2])/x_scale, t_add),
+                     (int(i[1])/x_scale, t_add),
+                    ] 
+
+	    codes = [
+		    Path.MOVETO,
+		    Path.LINETO,
+		    Path.LINETO,
+		    Path.LINETO,
+		    Path.CLOSEPOLY,
+		    ]
+	    transcript = t_rxe.search(i[9]).group(1)
+	    if transcript not in transcripts.keys():
+	        t_add += 1
+	    else: pass
+	    transcripts[transcript].append(Path(verts, codes))
+        else: pass
+    return(transcripts)
+
+
 def should_not_plot(x):
     """
     """