# get default parameters the papermill way.
input_adata = "../results/05_prepare_de_analysis/adata.h5ad"
output_adata = "../results/50_prepare_analysis_cd39/adata.h5ad"
results_dir = "../results/50_prepare_analysis_cd39/"
cpus = 32

# Parameters
input_adata = "input_adata.h5ad"
output_adata = "adata.h5ad"
cpus = "8"
results_dir = "."

import scanpy as sc
import matplotlib.pyplot as plt
import matplotlib
import numpy as np
import sys

sys.path.extend(("lib", "../lib"))
from jupytertools import *

fix_logging(sc.settings)
matplotlib.rcParams.update({"figure.autolayout": True, "figure.max_open_warning": 0})

# setup R integration
import rpy2.rinterface_lib.callbacks
import anndata2ri
from rpy2.robjects import pandas2ri
from rpy2.robjects.packages import importr

pandas2ri.activate()
anndata2ri.activate()
%load_ext rpy2.ipython

%%R
library(conflicted)
conflict_prefer("Position", "base")
library(dplyr)
library(ggplot2)
library(ggpubr)
library(ggbeeswarm)
library(edgeR)
library(magrittr)
options(max.print=100)
options(repr.matrix.max.cols=50, repr.matrix.max.rows=6)

R[write to console]: [conflicted] Will prefer base::Position over any other package

R[write to console]: Loading required package: magrittr

R[write to console]: Loading required package: limma

adata = sc.read_h5ad(input_adata)

sc.pl.umap(
    adata,
    color=["cell_type", "cluster", "ENTPD1", "FOXP3", "CD4"],
    legend_loc="on data",
    size=10,
    cmap="magma",
    ncols=2,
)

/home/sturm/scratch/projects/2020/kortekaas2020_paper/work/conda/run_notebook-b99176dfb05e0522988294512373169f/lib/python3.7/site-packages/IPython/core/pylabtools.py:128: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  fig.canvas.print_figure(bytes_io, **kw)

Define PD1+ CD8+ cells

We perform unsupervised clustering at high resolution and manually pick the clusters expressing NKG2a.

sc.pp.neighbors(adata, n_neighbors=20)

computing neighbors
    using 'X_pca' with n_pcs = 50

    finished

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

running Leiden clustering

    finished

sc.pl.umap(
    adata,
    color=["leiden", "ENTPD1"],
    legend_loc="on data",
    size=10,
    ncols=2,
    cmap="magma",
)

/home/sturm/scratch/projects/2020/kortekaas2020_paper/work/conda/run_notebook-b99176dfb05e0522988294512373169f/lib/python3.7/site-packages/IPython/core/pylabtools.py:128: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  fig.canvas.print_figure(bytes_io, **kw)

cluster_map = {
    "2" : "Treg",
    "19": "Treg",
    "21": "CD4",
    "9": "CD8"
}
adata.obs["cd39_status"] = [cluster_map[x] if x in cluster_map else "na" for x in adata.obs["leiden"]]

sc.pl.umap(
    adata,
    color=["cd39_status", "ENTPD1"],
    size=10,
    cmap="magma",
    legend_loc="on data",
    ncols=2,
)

... storing 'cd39_status' as categorical

adata.write_h5ad(output_adata, compression="lzf")

Prepare objects for edgeR

adata_edger = adata.copy()
# edgeR expects raw counts, normalization does not have any effect. We can therefore simply undo log1p
adata_edger.X = np.expm1(adata_edger.X)

de_dir = results_dir
!mkdir -p {de_dir}

%%R -i adata_edger -i de_dir
adata0 = adata_edger
dim(adata0)

[1]

%%R
adata = adata0[, colData(adata0)$cd39_status != "na"]
colData(adata)$cd39_status %<>% droplevels()
design = model.matrix(~0 + cd39_status + patient + n_genes + mt_frac, data=colData(adata))
contrasts = makeContrasts(
    cd4_cd8 = cd39_statusCD4 - cd39_statusCD8,
    cd4_treg = cd39_statusCD4 - cd39_statusTreg,
    cd8_treg = cd39_statusCD8 - cd39_statusTreg,
    levels = colnames(design)
)
save(adata, design, contrasts, file=paste0(de_dir, '/cd39_status.rda'), compress=FALSE)