# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: python3.7 (sc-tutorial)
#     language: python
#     name: sc-tutorial
# ---

# # Introduction
#
# This 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.

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
# -

# # Load data

# In the following, import the `csv` files (created in R) to convert them into an `anndata` file.

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

# Save to file.

sc.write(adata=adata, filename = path_dir + 'human_pancreas.h5ad')

adata = sc.read(path_dir + 'human_pancreas.h5ad')

# # 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. 

# + language="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)
# -

# Perform clustering.

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)


# We have to assign the clusters manually in order to call the size factor computation in R.

#Preprocess variables for scran normalization
input_groups = adatas[0].obs['groups']
data_mat = adata[adata.obs['tech']==protocol_id[0]].X.T

# + magic_args="-i data_mat -i input_groups -o size_factors" language="R"
#
# 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

# + magic_args="-i data_mat -i input_groups -o size_factors" language="R"
#
# 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

# + magic_args="-i data_mat -i input_groups -o size_factors" language="R"
#
# 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

# + magic_args="-i data_mat -i input_groups -o size_factors" language="R"
#
# 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

# + magic_args="-i data_mat -i input_groups -o size_factors" language="R"
#
# 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

# Now, let us concatenate all temporary anndata objects to one file in order to concatenate the size factors to one vector.

adata_tmp = adatas[0].concatenate(adatas[1:])

# Assign size factors to original data set collection.

adata.obs['size_factors'] = adata_tmp.obs['size_factors'].values

# Check size factor distribution.

sb.distplot(adata.obs['size_factors'], bins=50, kde=False)
plt.show()

# Free memory.

del adata_tmp
del adata_pp
del adatas

# Finally, normalise data.

#Normalize adata 
adata.X /= adata.obs['size_factors'].values[:,None]
sc.pp.log1p(adata)

# Save normalised data to file.

sc.write(adata=adata, filename=path_dir + 'human_pancreas_norm.h5ad')

StarcoderdataJupyterScriptsDedupFiltered