update 2

example/cosmx_analysis.py

@@ -0,0 +1,169 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Nov 12 13:18:24 2021
+@author: Ajit Johnson Nirmal
+COSMX analysis
+"""
+
+import scimap as sm
+import scanpy as sc
+import numpy as np
+import pandas as pd
+
+
+# %% Impot data
+# import cosmx data
+data_path = {
+"D:/cosmx/Lung5_Rep1/Lung5_Rep1-Flat_files_and_images/Lung5_Rep1_exprMat_file.csv": "D:/cosmx/Lung5_Rep1/Lung5_Rep1-Flat_files_and_images/Lung5_Rep1_metadata_file.csv",
+"D:/cosmx/Lung5_Rep2/Lung5_Rep2-Flat_files_and_images/Lung5_Rep2_exprMat_file.csv": "D:/cosmx/Lung5_Rep2/Lung5_Rep2-Flat_files_and_images/Lung5_Rep2_metadata_file.csv",
+"D:/cosmx/Lung5_Rep3/Lung5_Rep3-Flat_files_and_images/Lung5_Rep3_exprMat_file.csv": "D:/cosmx/Lung5_Rep3/Lung5_Rep3-Flat_files_and_images/Lung5_Rep3_metadata_file.csv",
+"D:/cosmx/Lung6/Lung6-Flat_files_and_images/Lung6_exprMat_file.csv": "D:/cosmx/Lung6/Lung6-Flat_files_and_images/Lung6_metadata_file.csv",
+"D:/cosmx/Lung9_Rep1/Lung9_Rep1-Flat_files_and_images/Lung9_Rep1_exprMat_file.csv": "D:/cosmx/Lung9_Rep1/Lung9_Rep1-Flat_files_and_images/Lung9_Rep1_metadata_file.csv",
+"D:/cosmx/Lung9_Rep2/Lung9_Rep2-Flat_files_and_images/Lung9_Rep2_exprMat_file.csv": "D:/cosmx/Lung9_Rep2/Lung9_Rep2-Flat_files_and_images/Lung9_Rep2_metadata_file.csv",
+"D:/cosmx/Lung12/Lung12-Flat_files_and_images/Lung12_exprMat_file.csv": "D:/cosmx/Lung12/Lung12-Flat_files_and_images/Lung12_metadata_file.csv",
+"D:/cosmx/Lung13/Lung13-Flat_files_and_images/Lung13_exprMat_file.csv": "D:/cosmx/Lung13/Lung13-Flat_files_and_images/Lung13_metadata_file.csv",
+            }
+
+data_path = {
+"D:/cosmx/Lung9_Rep1/Lung9_Rep1-Flat_files_and_images/Lung9_Rep1_exprMat_file.csv": "D:/cosmx/Lung9_Rep1/Lung9_Rep1-Flat_files_and_images/Lung9_Rep1_metadata_file.csv",
+  }
+
+adata = cosmx_to_scimap (data_path,
+                     unique_CellId=True,
+                     remove_NegPrb=True,
+                     random_sample=50000,
+                     log=False)
+
+adata.write('D:/cosmx/combined.h5ad')
+adata.write('D:/cosmx/combined_subset.h5ad')
+adata.write('D:/cosmx/Lung9_Rep2.h5ad')
+
+# read
+adata = sc.read('D:/cosmx/combined_subset.h5ad')
+adata = sc.read('D:/cosmx/Lung9_Rep2.h5ad')
+
+#%% Process
+
+adata.obs['fov'] = adata.obs['fov'].astype('category')
+
+# basic filtering
+sc.pp.filter_cells(adata, min_genes=30)
+sc.pp.filter_genes(adata, min_cells=10)
+
+
+
+# Norm
+sc.pp.normalize_total(adata, target_sum=10000, exclude_highly_expressed=True)
+
+# take square root of data
+#adata.X = np.sqrt(adata.X)
+sc.pp.log1p(adata)
+
+# save raw
+adata.raw = adata
+
+# Find MVGs
+sc.pp.highly_variable_genes(adata, n_top_genes=150)
+#sc.pl.highly_variable_genes(adata)
+adata = adata[:, adata.var.highly_variable]
+
+
+# custom normalization
+#adata.raw = adata
+#adata.X = adata.X / adata.X.sum(axis=1, keepdims=True) # Divide each element by total
+
+
+# PCA
+sc.tl.pca(adata, svd_solver='arpack')
+#sc.pl.pca(adata, color='CD74')
+#sc.pl.pca_variance_ratio(adata, log=True)
+
+# Network
+sc.pp.neighbors(adata, n_neighbors=50, n_pcs=10)
+adata.obs['log_CD3D'] = np.log1p(adata.obs['Mean.CD3'])
+
+# find genes
+#adata.var.index[adata.var.index.str.contains('CEACAM8', case=False)]
+
+#UMAP
+#sc.tl.tsne(adata)
+#sc.pl.tsne(adata, use_raw=False, color=['log_CD3D', 'ITGAX', 'CD163', 'CXCL8', 'CD3D', 'CD4','KRT7'], cmap='vlag')
+
+
+sc.tl.umap(adata,min_dist=0.4, spread=2.5)
+sc.pl.umap(adata, use_raw=True, color=['log_CD3D', 'ITGAX', 'CD163', 'CXCL8', 'CD3D', 'FCGR3A','KRT7','leiden'], cmap='vlag')
+
+#cluster
+sc.tl.leiden(adata, resolution=0.3)
+sc.pl.umap(adata, color=['leiden', 'imageid'])
+
+#marker genes
+sc.tl.rank_genes_groups(adata, 'leiden', method='t-test')
+sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)
+
+
+# real space
+plotly (adata,phenotype='leiden',image_id=None,x='CenterX_global_px',y='CenterY_global_px',size=3)
+
+# spatial lag analysisi
+adata = sm.tl.spatial_expression (adata, x_coordinate='CenterX_global_px',
+                                  y_coordinate='CenterY_global_px',
+                                  method='knn', knn=50, imageid='imageid', 
+                                  use_raw=True,subset=None,
+                                  label='spatial_expression_knn')
+
+
+# cluster
+adata = sm.tl.spatial_cluster (adata, k= 5, method = 'kmeans', df_name='spatial_expression_knn') # results will be saved under adata.obs['spatial_kmeans']
+
+#marker genes
+sc.tl.rank_genes_groups(adata, 'spatial_kmeans', method='t-test')
+sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)
+
+
+# view
+plotly (adata,phenotype='spatial_kmeans',image_id=None,x='CenterX_global_px',y='CenterY_global_px',size=3)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

scimap/plotting/lasso_selector.py

@@ -0,0 +1,96 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Nov  5 14:17:09 2019
+@author: Ajit Johnson Nirmal
+Program: Lasso Selector for selecting regions of interest
+"""
+
+# Load packages
+from __future__ import print_function
+#from six.moves import input
+import numpy as np
+from matplotlib.widgets import LassoSelector
+from matplotlib.path import Path
+from IPython import get_ipython
+get_ipython().run_line_magic('matplotlib', 'auto')
+import matplotlib.pyplot as plt
+
+class SelectFromCollection(object):
+    """Select indices from a matplotlib collection using `LassoSelector`.
+
+    Selected indices are saved in the `ind` attribute. This tool highlights
+    selected points by fading them out (i.e., reducing their alpha values).
+    If your collection has alpha < 1, this tool will permanently alter them.
+
+    Note that this tool selects collection objects based on their *origins*
+    (i.e., `offsets`).
+
+    Parameters
+    ----------
+    ax : :class:`~matplotlib.axes.Axes`
+        Axes to interact with.
+
+    collection : :class:`matplotlib.collections.Collection` subclass
+        Collection you want to select from.
+
+    alpha_other : 0 <= float <= 1
+        To highlight a selection, this tool sets all selected points to an
+        alpha value of 1 and non-selected points to `alpha_other`.
+    
+    Example:        
+        marker = 'CD45'
+        x = adata.obs['X_position']
+        y = adata.obs['Y_position']
+        
+        m_idx = adata.var.index.tolist().index(marker) # Get the index of marker of interest
+        tmp_dataframe = pd.DataFrame(adata.X)
+        hue = np.array(tmp_dataframe[m_idx])
+        
+        # Plotting    
+        fig, ax = plt.subplots()
+        pts = ax.scatter(x, y, s=1,c=hue,cmap='viridis')
+        ax.invert_yaxis()
+        
+        # Function call to do the Lasso selection
+        selector = SelectFromCollection(ax, pts)
+        # Return indeces of the selected points
+        tumor_idx= selector.ind
+        len(tumor_idx)
+        
+        # Update adata
+        adata.obs.loc[adata.obs.index[tumor_idx], 'roi-1'] = "roi-1"
+        adata.obs['roi-1'].value_counts() # Checking
+        
+        """
+
+    def __init__(self, ax, collection,alpha_other=0.3):
+        self.canvas = ax.figure.canvas
+        self.collection = collection
+        self.alpha_other = alpha_other
+
+        self.xys = collection.get_offsets()
+        self.Npts = len(self.xys)
+
+        # Ensure that we have separate colors for each object
+        self.fc = collection.get_facecolors()
+        if len(self.fc) == 0:
+            raise ValueError('Collection must have a facecolor')
+        elif len(self.fc) == 1:
+            self.fc = np.tile(self.fc, self.Npts).reshape(self.Npts, -1)
+
+        self.lasso = LassoSelector(ax, onselect=self.onselect)
+        self.ind = []
+
+    def onselect(self, verts):
+        path = Path(verts)
+        self.ind = np.nonzero([path.contains_point(xy) for xy in self.xys])[0]
+        self.fc[:, -1] = self.alpha_other
+        self.fc[self.ind, -1] = 1
+        self.collection.set_facecolors(self.fc)
+        self.canvas.draw_idle()
+
+    def disconnect(self):
+        self.lasso.disconnect_events()
+        self.fc[:, -1] = 1
+        self.collection.set_facecolors(self.fc)
+        self.canvas.draw_idle()