<jupyter_start><jupyter_text>IntroductionThis notebook performs the data normalisation of the collection of pancreas data. This collection consists of the following data sets:* CelSeq (GSE81076)* CelSeq2 (GSE85241)* Fluidigm C1 (GSE86469)* SMART-Seq2 (E-MTAB-5061)* inDrop (GSE84133)* SMARTer (GSE81608)For all but the SMARTer dataset, we perform *scran* normalisation. Subsequently, we log-scale all data sets. Import packages.<jupyter_code>import sys
# I have to add paths to the following github repos
sys.path.append('~/Python/scanpy/')
sys.path.append('~/Python/anndata/')
sys.path.append('~/Python/anndata2ri/')
import anndata
import scanpy as sc
import scipy as sci
sc.logging.print_versions()
sc.settings.verbosity = 3
import numpy as np
import scipy as sp
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import rcParams
from matplotlib import colors
import seaborn as sb
import datetime
now = datetime.datetime.now()
today = now.strftime("%Y%m%d")
import rpy2.rinterface_lib.callbacks
import logging

from rpy2.robjects import pandas2ri
import anndata2ri
rpy2.rinterface_lib.callbacks.logger.setLevel(logging.ERROR)

# Automatically convert rpy2 outputs to pandas dataframes
pandas2ri.activate()
anndata2ri.activate()
%load_ext rpy2.ipython<jupyter_output><empty_output><jupyter_text>Load data In the following, import the `csv` files (created in R) to convert them into an `anndata` file.<jupyter_code>path_dir = '../data/'
adata = sc.read_csv(filename = path_dir + 'human_pancreas.csv').T
metadata = pd.read_csv(path_dir + 'human_pancreas_meta.csv', index_col=0)
metadata
adata.obs['tech'] = metadata['tech'].values
adata.obs['celltype'] = metadata['celltype'].values<jupyter_output><empty_output><jupyter_text>Save to file.<jupyter_code>sc.write(adata=adata, filename = path_dir + 'human_pancreas.h5ad')
adata = sc.read(path_dir + 'human_pancreas.h5ad')<jupyter_output><empty_output><jupyter_text>Normalize with scran per data set Here, we split the data by dataset (i.e. tech variable). All but the 'smarter' dataset (which is in RPKM) will be normalised with scran.<jupyter_code>%%R

# Load all the R libraries we will be using in the notebook
library(scran)
#library(RColorBrewer)
#library(slingshot)
#library(monocle)
#library(gam)
#library(clusterExperiment)
#library(ggplot2)
#library(plyr)
#library(MAST)<jupyter_output><empty_output><jupyter_text>Perform clustering.<jupyter_code>a,b = np.unique(adata.obs['tech'], return_index=True)
protocol_id = a[np.argsort(b)]
adatas =[]
for protocol in enumerate(a[np.argsort(b)]):
    print(protocol[1])
    adata_pp = adata[adata.obs['tech']==protocol[1]].copy()
    if protocol[1]!='smarter':
        sc.pp.normalize_per_cell(adata_pp, counts_per_cell_after=1e6)
        sc.pp.log1p(adata_pp)
        sc.pp.pca(adata_pp, n_comps=15, svd_solver='arpack')
        sc.pp.neighbors(adata_pp)
        sc.tl.louvain(adata_pp, key_added='groups', resolution=0.5)
        #Preprocess variables for scran normalization
    adatas.append(adata_pp)<jupyter_output>celseq
normalizing by total count per cell
    finished (0:00:00): normalized adata.X and added    'n_counts', counts per cell before normalization (adata.obs)
computing PCA with n_comps = 15
    finished (0:00:01)
computing neighbors
    using 'X_pca' with n_pcs = 15
    finished: added to `.uns['neighbors']`
    'distances', distances for each pair of neighbors
    'connectivities', weighted adjacency matrix (0:00:00)
running Louvain clustering
    using the "louvain" package of Traag (2017)
    finished: found 7 clusters and added
    'groups', the cluster labels (adata.obs, categorical) (0:00:00)
celseq2
normalizing by total count per cell
    finished (0:00:00): normalized adata.X and added    'n_counts', counts per cell before normalization (adata.obs)
computing PCA with n_comps = 15
    finished (0:00:02)
computing neighbors
    using 'X_pca' with n_pcs = 15
    finished: added to `.uns['neighbors']`
    'distances', distances for each pair of neighbors
    'connectivities', wei[...]<jupyter_text>We have to assign the clusters manually in order to call the size factor computation in R.<jupyter_code>#Preprocess variables for scran normalization
input_groups = adatas[0].obs['groups']
data_mat = adata[adata.obs['tech']==protocol_id[0]].X.T
%%R -i data_mat -i input_groups -o size_factors

size_factors = computeSumFactors(data_mat, clusters=input_groups, min.mean=0.1)
# Visualize the estimated size factors
adatas[0].obs['size_factors'] = size_factors

sc.pl.scatter(adatas[0], 'size_factors', 'n_counts')

sb.distplot(size_factors, bins=50, kde=False)
plt.show()
#Preprocess variables for scran normalization
input_groups = adatas[1].obs['groups']
data_mat = adata[adata.obs['tech']==protocol_id[1]].X.T
%%R -i data_mat -i input_groups -o size_factors

size_factors = computeSumFactors(data_mat, clusters=input_groups, min.mean=0.1)
# Visualize the estimated size factors
adatas[1].obs['size_factors'] = size_factors

sc.pl.scatter(adatas[1], 'size_factors', 'n_counts')


sb.distplot(size_factors, bins=50, kde=False)
plt.show()
#Preprocess variables for scran normalization
input_groups = adatas[2].obs['groups']
data_mat = adata[adata.obs['tech']==protocol_id[2]].X.T
%%R -i data_mat -i input_groups -o size_factors

size_factors = computeSumFactors(data_mat, clusters=input_groups, min.mean=0.1)
# Visualize the estimated size factors
adatas[2].obs['size_factors'] = size_factors

sc.pl.scatter(adatas[2], 'size_factors', 'n_counts')


sb.distplot(size_factors, bins=50, kde=False)
plt.show()
#Preprocess variables for scran normalization
input_groups = adatas[3].obs['groups']
data_mat = adata[adata.obs['tech']==protocol_id[3]].X.T
%%R -i data_mat -i input_groups -o size_factors

size_factors = computeSumFactors(data_mat, clusters=input_groups, min.mean=0.1)
# Visualize the estimated size factors
adatas[3].obs['size_factors'] = size_factors

sc.pl.scatter(adatas[3], 'size_factors', 'n_counts')


sb.distplot(size_factors, bins=50, kde=False)
plt.show()
#Preprocess variables for scran normalization
input_groups = adatas[4].obs['groups']
data_mat = adata[adata.obs['tech']==protocol_id[4]].X.T
%%R -i data_mat -i input_groups -o size_factors

size_factors = computeSumFactors(data_mat, clusters=input_groups, min.mean=0.1)
# Visualize the estimated size factors
adatas[4].obs['size_factors'] = size_factors

sc.pl.scatter(adatas[4], 'size_factors', 'n_counts')


sb.distplot(size_factors, bins=50, kde=False)
plt.show()
# Set size factor to 1 for RPKM values
adatas[5].obs['size_factors'] = 1<jupyter_output><empty_output><jupyter_text>Now, let us concatenate all temporary anndata objects to one file in order to concatenate the size factors to one vector.<jupyter_code>adata_tmp = adatas[0].concatenate(adatas[1:])<jupyter_output><empty_output><jupyter_text>Assign size factors to original data set collection.<jupyter_code>adata.obs['size_factors'] = adata_tmp.obs['size_factors'].values<jupyter_output><empty_output><jupyter_text>Check size factor distribution.<jupyter_code>sb.distplot(adata.obs['size_factors'], bins=50, kde=False)
plt.show()<jupyter_output><empty_output><jupyter_text>Free memory.<jupyter_code>del adata_tmp
del adata_pp
del adatas<jupyter_output><empty_output><jupyter_text>Finally, normalise data.<jupyter_code>#Normalize adata 
adata.X /= adata.obs['size_factors'].values[:,None]
sc.pp.log1p(adata)<jupyter_output><empty_output><jupyter_text>Save normalised data to file.<jupyter_code>sc.write(adata=adata, filename=path_dir + 'human_pancreas_norm.h5ad')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12198338388500593, 'fraction_numerical': 0.019781089278649612, 'mean_word_length': 6.449901768172888, 'pattern_counts': {'":': 0, '<': 48, '<?xml version=': 0, '>': 48, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '291982', 'license': 'MIT', 'n_tokens_mistral': 2821, 'n_tokens_neox': 2628, 'n_words': 621, 'path': 'notebooks/data_preprocessing/pancreas/02_normalize_human_pancreas.ipynb', 'repo_name': 'Hrovatin/scib', 'size': '7583'}