Load data

library(conflicted)
conflict_prefer("Position", "base")
library(dplyr)
conflict_prefer("filter", "dplyr")
conflict_prefer("select", "dplyr")
library(ggplot2)
library(ggpubr)
library(ggbeeswarm)
library(readr)
library(tidyr)
library(reticulate)
library(clusterProfiler)
library(foreach)
library(doMC)
library(org.Hs.eg.db)
library(cowplot)
library(stringr)
library(kableExtra)
library(EnhancedVolcano)
knitr::knit_engines$set(python = reticulate::eng_python)
reticulate::py_available(TRUE)

## [1] TRUE

# bug in rstudio/reticulate:
matplotlib <- import("matplotlib")
matplotlib$use("Agg", force = TRUE)
registerDoMC(cores=16)

import scanpy as sc

## /home/sturm/scratch/projects/2020/kortekaas2020_paper/work/conda/vanderburg_oropharyngeal_cancer-ca9006a72cd4c32ca83708ebea1b2975/lib/python3.7/site-packages/dask/config.py:161: YAMLLoadWarning: calling yaml.load() without Loader=... is deprecated, as the default Loader is unsafe. Please read https://msg.pyyaml.org/load for full details.
##   data = yaml.load(f.read()) or {}

import matplotlib.pyplot as plt
import matplotlib
import sys
import seaborn as sns

sys.path.extend(("lib", "../lib"))
from jupytertools import *
fix_logging(sc.settings)
matplotlib.rcParams.update({"figure.autolayout": True, "figure.max_open_warning": 0})

adata = sc.read_h5ad(r.params['input_file'])

obs = py$adata$obs

get_path_cd39 = function(file_name) {file.path(params$input_de_res_dir_cd39, file_name)}

de_pairwise = read_tsv(get_path_cd39("cd39_status.rda.res.tsv"))

writexl::write_xlsx(de_pairwise, "pairwise_comparison.xlsx")

DE-analysis

Differntial expression analysis (pairwise comparison of clusters)

How to read the plot: cd4_cd8 shows genes that are differentially expressed between cells in the CD4 and CD8 cluster. Blue: Enriched in CD8, under-represented in CD4; Orange: Enriched in CD4, uner-represented in CD8. Analogous for the other two comparisons.

## Joining, by = "cluster"

DE as volcano plot

## [[1]]

## 
## [[2]]

## 
## [[3]]

Violin Plots

sc.pl.dotplot(adata_cd39, groupby="cd39_status", 
                     var_names=genes)

## GridSpec(2, 5, height_ratios=[0, 10.5], width_ratios=[5.6, 0, 0.2, 0.5, 0.25])
## 
## /data/scratch/sturm/scratch/projects/2020/kortekaas2020_paper/work/conda/vanderburg_oropharyngeal_cancer-ca9006a72cd4c32ca83708ebea1b2975/lib/R/library/reticulate/python/rpytools/call.py:13: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
##   res = rpycall.call_r_function(f, *args, **kwargs)

sc.pl.stacked_violin(adata_cd39, groupby="cd39_status", 
                     var_names=genes, 
                                figsize=(6, 3))

## [<matplotlib.axes._subplots.AxesSubplot object at 0x2b704d346fd0>, <matplotlib.axes._subplots.AxesSubplot object at 0x2b704d2eee48>, <matplotlib.axes._subplots.AxesSubplot object at 0x2b704d40d898>]

sc.pl.stacked_violin(adata_cd39, groupby="cd39_status", 
                     var_names=genes, swap_axes=True, figsize=(5, 12), use_raw=True, save="_wide.svg")

## saving figure to file figures/stacked_violin_wide.svg
## [<matplotlib.axes._subplots.AxesSubplot object at 0x2b704d525c50>, <matplotlib.axes._subplots.AxesSubplot object at 0x2b704d7c28d0>, <matplotlib.axes._subplots.AxesSubplot object at 0x2b704d8809b0>, <matplotlib.axes._subplots.AxesSubplot object at 0x2b704d89f748>, <matplotlib.axes._subplots.AxesSubplot object at 0x2b704d8d8828>, <matplotlib.axes._subplots.AxesSubplot object at 0x2b704d912860>, <matplotlib.axes._subplots.AxesSubplot object at 0x2b704d94a5f8>, <matplotlib.axes._subplots.AxesSubplot object at 0x2b704d985780>, <matplotlib.axes._subplots.AxesSubplot object at 0x2b704d9a5a58>, <matplotlib.axes._subplots.AxesSubplot object at 0x2b704d9f95f8>, <matplotlib.axes._subplots.AxesSubplot object at 0x2b704da1d978>, <matplotlib.axes._subplots.AxesSubplot object at 0x2b704da75518>, <matplotlib.axes._subplots.AxesSubplot object at 0x2b704daaaf60>, <matplotlib.axes._subplots.AxesSubplot object at 0x2b704db18780>, <matplotlib.axes._subplots.AxesSubplot object at 0x2b704dad6ac8>, <matplotlib.axes._subplots.AxesSubplot object at 0x2b704db569e8>]