In [1]:
%matplotlib inline
%load_ext autoreload
%autoreload 2
%load_ext line_profiler
In [2]:
import scanpy as sc
import random
from unicoord import scu
from unicoord.visualization import draw_loss_curves
import torch
from line_profiler import LineProfiler
In [3]:
sc.settings.verbosity = 3             # verbosity: errors (0), warnings (1), info (2), hints (3)
sc.logging.print_header()
# sc.settings.set_figure_params(dpi=80, facecolor='white')
sc.settings.set_figure_params(vector_friendly=False)

scanpy==1.9.1 anndata==0.7.5 umap==0.3.10 numpy==1.17.3 scipy==1.7.3 pandas==1.3.5 scikit-learn==0.22 statsmodels==0.10.2 python-igraph==0.7.1 louvain==0.6.1 leidenalg==0.7.0

load hECA lung data

In [21]:
adata = sc.read_h5ad('h5ad/Lung.Adult.pp.h5ad')
In [22]:
sc.pl.embedding(adata, 'X_umap',legend_loc='on data', legend_fontsize=10,
                color=['leiden','donor_id','seq_tech','cell_type'], ncols=2)
In [5]:
adata = adata.raw.to_adata()
sc.pp.normalize_total(adata, target_sum=1e4 ,exclude_highly_expressed= True)
sc.pp.log1p(adata)
adata

normalizing counts per cell The following highly-expressed genes are not considered during normalization factor computation:
['CCL21', 'FTL', 'SCGB1A1', 'SCGB3A1', 'SCGB3A2', 'SFTPA2', 'SFTPC']
    finished (0:00:00)
Out[5]:
AnnData object with n_obs × n_vars = 54615 × 20770
    obs: 'user_id', 'study_id', 'cell_id', 'organ', 'region', 'subregion', 'seq_tech', 'sample_status', 'donor_id', 'donor_gender', 'donor_age', 'original_name', 'cl_name', 'hcad_name', 'tissue_type', 'cell_type', 'marker_gene', 'cid', 'RNA_snn_res.0.4', 'seurat_clusters', 'nCount_RNA', 'nFeature_RNA', 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'leiden'
    var: 'n_cells', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts'
    uns: 'cell_type_colors', 'hvg', 'leiden', 'leiden_colors', 'neighbors', 'pca', 'study_id_colors', 'tissue_type_colors', 'umap', 'log1p'
    obsm: 'X_pca', 'X_umap'
In [6]:
adata.obs.donor_id.value_counts()
Out[6]:
AdultLung3        8930
AdultLung1_21Y    8415
390C              8162
356C              7663
367C              7320
AdultLung2_21Y    6006
343B              4385
368C              3381
Braga2019_2        167
Braga2019_1         92
Braga2019_3         52
Braga2019_4         42
Name: donor_id, dtype: int64
In [7]:
adata.obs.seq_tech.value_counts()
Out[7]:
10X              31264
Microwell-seq    23351
Name: seq_tech, dtype: int64

supervise batch

In [8]:
scu.model_unicoord_in_adata(adata, n_cont=30, n_diff=0, n_clus = [],
                            obs_fitting=['seq_tech'])
In [9]:
scu.train_unicoord_in_adata(adata, epochs=100, chunk_size=20000, slot = "cur")

/home/gaohx/anaconda3/envs/scvi/lib/python3.7/site-packages/pandas/core/arrays/categorical.py:2631: FutureWarning: The `inplace` parameter in pandas.Categorical.remove_unused_categories is deprecated and will be removed in a future version.
  res = method(*args, **kwargs)

training chunk 1 / 3 of the data

training chunk 2 / 3 of the data

training chunk 3 / 3 of the data



















In [10]:



    


fig = draw_loss_curves(adata.uns['unc_stuffs']['loss'])
# if save_figs:
#     fig.savefig(os.path.join(savePath, 'img', 'fig1_lossCurves.png'))
fig.show()














































In [11]:



    


scu.embed_unicoord_in_adata(adata, chunk_size=5000)

sc.pp.neighbors(adata, use_rep='unicoord')

sc.tl.leiden(adata, resolution=0.5)

sc.tl.umap(adata)


























computing neighbors

















/home/gaohx/anaconda3/envs/scvi/lib/python3.7/site-packages/numba/typed_passes.py:271: NumbaPerformanceWarning: 
The keyword argument 'parallel=True' was specified but no transformation for parallel execution was possible.

To find out why, try turning on parallel diagnostics, see http://numba.pydata.org/numba-doc/latest/user/parallel.html#diagnostics for help.

File "../../../../../home/gaohx/anaconda3/envs/scvi/lib/python3.7/site-packages/umap/rp_tree.py", line 135:
@numba.njit(fastmath=True, nogil=True, parallel=True)
def euclidean_random_projection_split(data, indices, rng_state):
^

  state.func_ir.loc))
/home/gaohx/anaconda3/envs/scvi/lib/python3.7/site-packages/umap/nndescent.py:92: NumbaPerformanceWarning: 
The keyword argument 'parallel=True' was specified but no transformation for parallel execution was possible.

To find out why, try turning on parallel diagnostics, see http://numba.pydata.org/numba-doc/latest/user/parallel.html#diagnostics for help.

File "../../../../../home/gaohx/anaconda3/envs/scvi/lib/python3.7/site-packages/umap/utils.py", line 409:
@numba.njit(parallel=True)
def build_candidates(current_graph, n_vertices, n_neighbors, max_candidates, rng_state):
^

  current_graph, n_vertices, n_neighbors, max_candidates, rng_state
/home/gaohx/anaconda3/envs/scvi/lib/python3.7/site-packages/numba/typed_passes.py:271: NumbaPerformanceWarning: 
The keyword argument 'parallel=True' was specified but no transformation for parallel execution was possible.

To find out why, try turning on parallel diagnostics, see http://numba.pydata.org/numba-doc/latest/user/parallel.html#diagnostics for help.

File "../../../../../home/gaohx/anaconda3/envs/scvi/lib/python3.7/site-packages/umap/nndescent.py", line 47:
    @numba.njit(parallel=True)
    def nn_descent(
    ^

  state.func_ir.loc))

















    finished: added to `.uns['neighbors']`
    `.obsp['distances']`, distances for each pair of neighbors
    `.obsp['connectivities']`, weighted adjacency matrix (0:00:10)
running Leiden clustering
    finished: found 52 clusters and added
    'leiden', the cluster labels (adata.obs, categorical) (0:00:14)
computing UMAP
    finished: added
    'X_umap', UMAP coordinates (adata.obsm) (0:00:55)



















In [12]:



    


sc.pl.embedding(adata, 'X_umap',legend_loc='on data', legend_fontsize=10,
                color=['leiden','donor_id','seq_tech','cell_type'], ncols=2)


















































not supervise cell type












In [13]:



    


scu.model_unicoord_in_adata(adata, n_cont=100, n_diff=0, n_clus = [],
                            obs_fitting=[])

















In [14]:



    


scu.train_unicoord_in_adata(adata, epochs=100, chunk_size=20000, slot = "cur")


























/home/gaohx/anaconda3/envs/scvi/lib/python3.7/site-packages/pandas/core/arrays/categorical.py:2631: FutureWarning: The `inplace` parameter in pandas.Categorical.remove_unused_categories is deprecated and will be removed in a future version.
  res = method(*args, **kwargs)

















training chunk 1 / 3 of the data

















 
 






















training chunk 2 / 3 of the data

















 
 






















training chunk 3 / 3 of the data

















 
 








































In [15]:



    


fig = draw_loss_curves(adata.uns['unc_stuffs']['loss'])
# if save_figs:
#     fig.savefig(os.path.join(savePath, 'img', 'fig1_lossCurves.png'))
fig.show()














































In [16]:



    


scu.embed_unicoord_in_adata(adata, chunk_size=5000)

















In [17]:



    


sc.pp.neighbors(adata, use_rep='unicoord')


























computing neighbors

















/home/gaohx/anaconda3/envs/scvi/lib/python3.7/site-packages/numba/typed_passes.py:271: NumbaPerformanceWarning: 
The keyword argument 'parallel=True' was specified but no transformation for parallel execution was possible.

To find out why, try turning on parallel diagnostics, see http://numba.pydata.org/numba-doc/latest/user/parallel.html#diagnostics for help.

File "../../../../../home/gaohx/anaconda3/envs/scvi/lib/python3.7/site-packages/umap/nndescent.py", line 47:
    @numba.njit(parallel=True)
    def nn_descent(
    ^

  state.func_ir.loc))

















    finished: added to `.uns['neighbors']`
    `.obsp['distances']`, distances for each pair of neighbors
    `.obsp['connectivities']`, weighted adjacency matrix (0:00:08)



















In [18]:



    


sc.tl.leiden(adata, resolution=0.5)


























running Leiden clustering
    finished: found 44 clusters and added
    'leiden', the cluster labels (adata.obs, categorical) (0:00:18)



















In [19]:



    


sc.tl.umap(adata)


























computing UMAP
    finished: added
    'X_umap', UMAP coordinates (adata.obsm) (0:00:58)



















In [20]:



    


sc.pl.embedding(adata, 'X_umap',legend_loc='on data', legend_fontsize=10,
                color=['leiden','donor_id','seq_tech','cell_type'], ncols=2)