adding simple spatial capability

Author: scyrusm

panopticon/analysis.py

@@ -2095,7 +2095,8 @@ def get_cluster_differential_expression(loom,
                                         ident2_downsample_size=None,
                                         min_cluster_size=0,
                                         gene_alternate_name=None,
-                                        gene_subset_mask=None):
+                                        gene_subset_mask=None,
+                                        alternate_test=None):
     """
 
     Parameters
@@ -2242,6 +2243,8 @@ def get_cluster_differential_expression(loom,
         meanexpexpr2 = np.mean(2**data2, axis=1)
         fracexpr1 = np.mean(data1 > 0, axis=1)
         fracexpr2 = np.mean(data2 > 0, axis=1)
+        if alternate_test is not None:
+            alternate_test_pvalues = alternate_test(data1, data2, axis=1).pvalue
 
 
 ###
@@ -2255,6 +2258,8 @@ def get_cluster_differential_expression(loom,
         meanexpexpr2 = []
         fracexpr1 = []
         fracexpr2 = []
+        if alternate_test is not None:
+            alternate_test_pvalues = []
         for igene, gene in enumerate(
                 tqdm(genelist, desc='Computing Mann-Whitney p-values')):
             genes.append(gene)
@@ -2273,6 +2278,8 @@ def get_cluster_differential_expression(loom,
             meanexpexpr2.append(np.mean(2**data2[igene, :]))
             fracexpr1.append((data1[igene, :] > 0).mean())
             fracexpr2.append((data2[igene, :] > 0).mean())
+            if alternate_test is not None:
+                alternate_test_pvalues.append(alternate_test(data1[igene,:], data2[igene,:]).pvalue)
     output = pd.DataFrame(genelist)
     output.columns = ['gene']
     output['pvalue'] = pvalues
@@ -2286,6 +2293,9 @@ def get_cluster_differential_expression(loom,
         meanexpexpr2) / np.log(2)
     output['FracExpr1'] = fracexpr1
     output['FracExpr2'] = fracexpr2
+    if alternate_test is not None:
+        output['AlternateTestPValue'] = [1 if np.isnan(x) else x for x in alternate_test_pvalues]
+        output['AlternateTestQvalue'] = fdrcorrection(output['AlternateTestPValue'],is_sorted=False)[1]
     if gene_alternate_name is not None:
         gene2altname = {
             gene: altname

panopticon/utilities.py

@@ -461,6 +461,7 @@ def cohensd(g1, g2):
 
     return (np.mean(g1) - np.mean(g2)) / s
 
+
 def glassdelta(experimental=None, control=None):
     """Returns Glass' delta for the effect size of group 1 values (g1) over group 2 values (g2). g2 is assumed to be the control group
 
@@ -477,7 +478,8 @@ def glassdelta(experimental=None, control=None):
     
     """
     if (experimental is None) or (control is None):
-        raise Exception("experimental and control should be lists or numpy vectors")
+        raise Exception(
+            "experimental and control should be lists or numpy vectors")
     scontrol = np.std(control, ddof=1)
 
     return (np.mean(experimental) - np.mean(control)) / scontrol
@@ -557,7 +559,9 @@ def convert_10x_h5(path_10x_h5,
                    gene_whitelist=None,
                    output_type='loom',
                    write_chunked=False,
-                   chunk_size=512):
+                   chunk_size=512,
+                   exclude_feature_type=None,
+                   verbose=True):
     """
 
     Parameters
@@ -618,6 +622,30 @@ def convert_10x_h5(path_10x_h5,
         ca[labelkey] = [label] * len(barcodes)
 
     m = filtered_feature_bc_matrix.m
+
+    feature_types = np.array(filtered_feature_bc_matrix.feature_ref.get_feature_types_excluding_deprecated_probes() )
+    if exclude_feature_type is not None:
+        if type(exclude_feature_type)==str:
+            exclude_feature_type = [exclude_feature_type]
+        elif type(exclude_feature_type)==tuple:
+            exclude_feature_type = list(exclude_feature_type)
+
+        try:
+            iterator = iter(exclude_feature_type)
+        except TypeError:
+            raise Exception("exclude_feature_type must be iterable or str, is of type {}".format(type(exclude_feature_type)))
+        feature_type_mask = ~np.isin(feature_types, exclude_feature_type)
+        feature_types = np.array(feature_types)[np.array(feature_type_mask)]
+        m = m[feature_type_mask, :]
+        features = list(np.array(features)[feature_type_mask])
+        features_common_names = list(np.array(features_common_names)[feature_type_mask])
+
+        if verbose:
+            print("Removing the following feature types: {}".format(
+                exclude_feature_type))
+    if verbose:
+        print(
+            "Including the following feature types: {}".format(np.unique(feature_types)))
     if gene_whitelist is not None:
         if len(gene_whitelist) > 0:
             mask = np.isin(features, gene_whitelist)

panopticon/visualization.py

@@ -1030,14 +1030,14 @@ def position_to_xytext_habt_va(position, effect_size, negpval, maxx, maxy,
                 "\'{}\':  invalid position character selection".format(
                     position))
         return xytext, habt, va
+
     negpvals = []
     for gene in genemarklist:
         genedf = diffex[diffex[gene_column] == gene]
         negpval = -np.log(genedf.iloc[0][pval_column]) / np.log(10)
         negpvals.append(negpval)
     genemarklist = list(np.array(genemarklist)[np.argsort(negpvals)][::-1])
 
-
     if positions != 'side':
         if type(positions) == dict:
             positions = [positions[key] for key in genemarklist]
@@ -1052,7 +1052,8 @@ def position_to_xytext_habt_va(position, effect_size, negpval, maxx, maxy,
             effect_size = genedf.iloc[0][effect_size_col]
             ax.scatter(effect_size, negpval, marker='.', color='k')
             xytext, habt, va = position_to_xytext_habt_va(
-                position, effect_size, negpval, maxx, maxy, gene_label_offset_scale)
+                position, effect_size, negpval, maxx, maxy,
+                gene_label_offset_scale)
             anno = ax.annotate(
                 gene, (effect_size, negpval),
                 xytext,
@@ -1073,37 +1074,46 @@ def position_to_xytext_habt_va(position, effect_size, negpval, maxx, maxy,
             if gene in gene_position_dict_for_side_annotations.keys():
                 position = gene_position_dict_for_side_annotations[gene]
                 xytext, habt, va = position_to_xytext_habt_va(
-                    position, effect_size, negpval, maxx, maxy, gene_label_offset_scale)
-                anno = ax.annotate(
-                    gene, (effect_size, negpval),
-                    xytext,
-                    va=va,
-                    ha=habt,
-                    path_effects=[pe.withStroke(linewidth=2, foreground="white")])
+                    position, effect_size, negpval, maxx, maxy,
+                    gene_label_offset_scale)
+                anno = ax.annotate(gene, (effect_size, negpval),
+                                   xytext,
+                                   va=va,
+                                   ha=habt,
+                                   path_effects=[
+                                       pe.withStroke(linewidth=2,
+                                                     foreground="white")
+                                   ])
             else:
-    
+
                 if effect_size < no_effect_line:
-                    xytext = (left_edge + .03 * maxx * side_annotation_gene_label_offset_scale ,
-                              top_edge - lcounter)
+                    xytext = (
+                        left_edge +
+                        .03 * maxx * side_annotation_gene_label_offset_scale,
+                        top_edge - lcounter)
                     habt = 'left'
-                    va='center'
+                    va = 'center'
                     lcounter += counterscale
                 else:
-                    xytext = (right_edge - .03 * maxx * side_annotation_gene_label_offset_scale,
-                              top_edge - rcounter)
+                    xytext = (
+                        right_edge -
+                        .03 * maxx * side_annotation_gene_label_offset_scale,
+                        top_edge - rcounter)
                     habt = 'right'
-                    va='center'
+                    va = 'center'
                     rcounter += counterscale
-                anno = ax.annotate(
-                    gene, (effect_size, negpval),
-                    xytext=xytext,
-                    va=va,
-                    ha=habt,
-                    arrowprops=dict(facecolor='black',
-                                    width=0.1,
-                                    headwidth=0,
-                                    alpha=0.25),
-                    path_effects=[pe.withStroke(linewidth=2, foreground="white")])
+                anno = ax.annotate(gene, (effect_size, negpval),
+                                   xytext=xytext,
+                                   va=va,
+                                   ha=habt,
+                                   arrowprops=dict(facecolor='black',
+                                                   width=0.1,
+                                                   headwidth=0,
+                                                   alpha=0.25),
+                                   path_effects=[
+                                       pe.withStroke(linewidth=2,
+                                                     foreground="white")
+                                   ])
             if draggable_annotations:
                 anno.draggable()
             ax.scatter(effect_size, negpval, marker='.', color='k')
@@ -2492,11 +2502,15 @@ def plot_color_coded_embedding(loom,
                                y_ca,
                                category_ca=None,
                                category_as_continuum=False,
+                               use_gex_as_ca=False,
+                               gene_ra='gene_common_name',
+                               gex_layer='log2(TP10k+1)',
                                fig=None,
                                ax=None,
                                color_palette='colorblind',
                                legend=True,
-                               on_figure_annotation=False):
+                               on_figure_annotation=False,
+                               s=2):
     import numpy as np
     import seaborn as sns
     if fig is not None:
@@ -2507,10 +2521,15 @@ def plot_color_coded_embedding(loom,
             raise Exception("Both or neither of fig, ax may be None")
         fig, ax = plt.subplots(figsize=(4, 4))
     if category_as_continuum:
+        if use_gex_as_ca:
+            igene = np.where(loom.ra[gene_ra] == category_ca)[0][0]
+            c = loom[gex_layer][igene, :]
+        else:
+            c = loom.ca[category_ca]
         g = ax.scatter(loom.ca[x_ca],
                        loom.ca[y_ca],
-                       s=2,
-                       c=loom.ca[category_ca],
+                       s=s,
+                       c=c,
                        cmap=color_palette)
         plt.colorbar(g, label=category_ca)
     else:
@@ -2531,7 +2550,7 @@ def plot_color_coded_embedding(loom,
         ax.scatter(
             loom.ca[x_ca][shuffle],
             loom.ca[y_ca][shuffle],
-            s=2,
+            s=s,
             c=[category2color[x] for x in loom.ca[category_ca][shuffle]])
         from matplotlib.lines import Line2D
         legend_elements = [
@@ -2583,3 +2602,70 @@ def plot_color_coded_embedding(loom,
     ax.set_ylabel(y_ca, loc='bottom', fontsize=14)
 
     return fig, ax
+
+
+def gsea_plot(ranking,
+              pathway2genelist_dict,
+              left_label='Enriched for genes at beginning of ranking',
+              right_label='Enriched for genes at end of ranking',
+              figsize=(9/2,7/2)
+              ):
+    import matplotlib.pyplot as plt
+    import seaborn as sns
+    from panopticon.analysis import get_enrichment_score
+    import matplotlib
+
+    palette = sns.color_palette('colorblind', 12)
+    fig, axes = plt.subplots(len(pathway2genelist_dict.keys()) + 1,
+                             1,
+                             figsize=figsize,
+                             height_ratios=[20] +
+                             [1] * len(pathway2genelist_dict.keys()),
+                             sharex=True)
+    mins = []
+    maxs = []
+    for ikey, key in enumerate([y for y in pathway2genelist_dict.keys()]):
+        es = get_enrichment_score(ranking,
+                                  pathway2genelist_dict[key],
+                                  presorted=True,
+                                  return_es_curve=True,
+                                  return_pvalue=True,
+                                  use_fgsea=True)
+        axes[0].plot(es.enrichment_score_curve,
+                     label=key,
+                     lw=3,
+                     color=palette[ikey])
+        maxs.append(np.max(es.enrichment_score_curve))
+        mins.append(np.min(es.enrichment_score_curve))
+
+        for x in np.where(np.isin(ranking, pathway2genelist_dict[key]))[0]:
+            axes[ikey + 1].axvline(x)
+        axes[ikey + 1].set_ylabel(key + '\n' * 5,
+                                  rotation=0,
+                                  ha='left',
+                                  va='bottom')
+        axes[ikey + 1].set_yticks([])
+        axes[ikey + 1].yaxis.set_label_position("right")
+        axes[ikey + 1].yaxis.tick_right()
+        for side in ['top', 'bottom', 'right', 'left']:
+            axes[ikey + 1].spines[side].set_visible(False)
+        axes[ikey + 1].set_xticks([])
+        if es.p_value == 0:
+            key_with_pval = key + ', p<{0:.5g}'.format(1 / 10000)
+        else:
+            key_with_pval = key + ', p={0:.5g}'.format(es.p_value)
+        axes[ikey + 1].set_ylabel(key_with_pval,
+                                  rotation=0,
+                                  ha='left',
+                                  va='center')
+
+    axes[0].legend(bbox_to_anchor=(1, 1))
+    axes[0].set_ylabel('running enrichment score')
+    axes[-1].set_xlabel('rank in gene list\n' + r'$\leftarrow$ ' + left_label +
+                        ' ' * 20 + right_label + r' $\rightarrow$')
+    axes[0].set_ylim([np.min(mins), np.max(maxs)])
+    axes[0].spines['top'].set_position(('data', 0))
+    axes[0].spines['bottom'].set_position(('data', 0))
+    axes[0].spines['right'].set_visible(False)
+    plt.tight_layout()
+    return fig, axes