<jupyter_start><jupyter_text>06a HSC sexual dimorpism 2 replicatesdifferential expression analysis using MAST to discover sexual dimorpism accross two replicates of unperturbed HSCs.**HSCs**Run this model:`zlmCond_all <- zlm(formula = ~female + n_genes + leiden, sca=sca)`Comparisons:compare both replicates of old and newdone with this docker image:docker run --rm -d --name scanpy -p 8883:8888 -e JUPYTER_ENABLE_LAB=YES -v /Users/efast/Documents/:/home/jovyan/work r_scanpy:vs5<jupyter_code>import scanpy as sc
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
from gprofiler import GProfiler

import rpy2.rinterface_lib.callbacks
import logging

from rpy2.robjects import pandas2ri
import anndata2ri
# Ignore R warning messages
#Note: this can be commented out to get more verbose R output
rpy2.rinterface_lib.callbacks.logger.setLevel(logging.ERROR)

# Automatically convert rpy2 outputs to pandas dataframes
pandas2ri.activate()
anndata2ri.activate()
%load_ext rpy2.ipython

plt.rcParams['figure.figsize']=(8,8) #rescale figures
sc.settings.verbosity = 3
#sc.set_figure_params(dpi=200, dpi_save=300)
sc.logging.print_versions()
%%R
# Load libraries from correct lib Paths for my environment - ignore this!
.libPaths(.libPaths()[c(3,2,1)])

# Load all the R libraries we will be using in the notebook
library(scran)
library(ggplot2)
library(plyr)
library(MAST)<jupyter_output><empty_output><jupyter_text>Replicate 1<jupyter_code># load data

adata = sc.read('./sc_objects/old_LT_preprocessed.h5ad', cache = True)
#Create new Anndata object for use in MAST with non-batch corrected data as before
adata_raw = adata.copy()
adata_raw = sc.AnnData(X=adata.raw.X, obs=adata.obs, var=adata.raw.var)
adata_raw.obs['n_genes'] = (adata_raw.X > 0).sum(1) # recompute number of genes expressed per cell
adata = None
adata_raw.obs.head()<jupyter_output><empty_output><jupyter_text>Run MAST on total cells - Select genes expressed in >5% of cells (no adaptive thresholding)<jupyter_code>%%R -i adata_raw

#Convert SingleCellExperiment to SingleCellAssay type as required by MAST
sca <- SceToSingleCellAssay(adata_raw, class = "SingleCellAssay")

#Scale Gene detection rate
colData(sca)$n_genes = scale(colData(sca)$n_genes)

# filter genes based on hard cutoff (have to be expressed in at least 5% of all cells)
freq_expressed <- 0.05
expressed_genes <- freq(sca) > freq_expressed
sca <- sca[expressed_genes,]<jupyter_output><empty_output><jupyter_text>everythingbackground:  `zlmCond_all <- zlm(formula = ~Female + n_genes + leiden, sca=sca)` this runs the model`a formula with the measurement variable (gene expression) on the LHS (left hand side) and predictors present in colData on the RHSexpression of genes controlling for cluster, condition, sex + n_genesquestions I can ask:sex differences controlling for treatmentssex differences controlling for clusters - not necessary analyze all the clustersoverall gene expression changes in treatment<jupyter_code>%%R 
#Define & run hurdle model 
zlmCond_all <- zlm(formula = ~Female + n_genes + leiden, sca=sca) # this runs the model
summaryCond_all <- summary(zlmCond_all, doLRT=TRUE) # extracts the data, gives datatable with summary of fit, doLRT=TRUE extracts likelihood ratio test p-value
summaryDt_all <- summaryCond_all$datatable # reformats into a table
%%R
head(summaryDt_all)
%%R -o female_all

# reformat for female
result_all_Female <- merge(summaryDt_all[contrast=='FemaleTRUE' & component=='H',.(primerid, `Pr(>Chisq)`)], #P-vals
                  summaryDt_all[contrast=='FemaleTRUE' & component=='logFC', .(primerid, coef)],
                  by='primerid') #logFC coefficients
#Correct for multiple testing (FDR correction) and filtering
result_all_Female[,FDR:=p.adjust(`Pr(>Chisq)`, 'fdr')] # create column named FDR - probably that p.adjust function
female_all = result_all_Female[result_all_Female$FDR<0.01,, drop=F] # create new table where rows with FDR<0.01 are droped
female_all = female_all[order(female_all$FDR),] # sorts the table
%%R -o MAST_raw_all

MAST_raw_all <- summaryDt_all
# save files as .csvs

MAST_raw_all.to_csv('./write/MAST_raw_old_LT.csv')
female_all.to_csv('./write/MAST_female_old_LT.csv')
%%R
# remove previous variables

rm(zlmCond_all)
rm(summaryDt_all)
rm(summaryCond_all)
rm(MAST_raw_all)<jupyter_output><empty_output><jupyter_text>Replicate 2<jupyter_code># load data
adata = sc.read('./sc_objects/new_ct_LT_preprocessed.h5ad', cache = True)
#Create new Anndata object for use in MAST with non-batch corrected data as before
adata_raw = adata.copy()
adata_raw = sc.AnnData(X=adata.raw.X, obs=adata.obs, var=adata.raw.var)
adata_raw.obs['n_genes'] = (adata_raw.X > 0).sum(1) # recompute number of genes expressed per cell
adata = None
adata_raw.obs.head()<jupyter_output><empty_output><jupyter_text>Run MAST on total cells - Select genes expressed in >5% of cells (no adaptive thresholding)<jupyter_code>%%R -i adata_raw

#Convert SingleCellExperiment to SingleCellAssay type as required by MAST
sca <- SceToSingleCellAssay(adata_raw, class = "SingleCellAssay")

#Scale Gene detection rate
colData(sca)$n_genes = scale(colData(sca)$n_genes)

# filter genes based on hard cutoff (have to be expressed in at least 5% of all cells)
freq_expressed <- 0.05
expressed_genes <- freq(sca) > freq_expressed
sca <- sca[expressed_genes,]<jupyter_output><empty_output><jupyter_text>everythingbackground:  `zlmCond_all <- zlm(formula = ~Female + n_genes + leiden, sca=sca)` this runs the model`a formula with the measurement variable (gene expression) on the LHS (left hand side) and predictors present in colData on the RHSexpression of genes controlling for cluster, condition, sex + n_genesquestions I can ask:sex differences controlling for treatmentssex differences controlling for clusters - not necessary analyze all the clustersoverall gene expression changes in treatment<jupyter_code>%%R 
#Define & run hurdle model 
zlmCond_all <- zlm(formula = ~Female + n_genes + leiden, sca=sca) # this runs the model
summaryCond_all <- summary(zlmCond_all, doLRT=TRUE) # extracts the data, gives datatable with summary of fit, doLRT=TRUE extracts likelihood ratio test p-value
summaryDt_all <- summaryCond_all$datatable # reformats into a table
%%R
head(summaryDt_all)
%%R -o female_all

# reformat for female
result_all_Female <- merge(summaryDt_all[contrast=='FemaleTRUE' & component=='H',.(primerid, `Pr(>Chisq)`)], #P-vals
                  summaryDt_all[contrast=='FemaleTRUE' & component=='logFC', .(primerid, coef)],
                  by='primerid') #logFC coefficients
#Correct for multiple testing (FDR correction) and filtering
result_all_Female[,FDR:=p.adjust(`Pr(>Chisq)`, 'fdr')] # create column named FDR - probably that p.adjust function
female_all = result_all_Female[result_all_Female$FDR<0.01,, drop=F] # create new table where rows with FDR<0.01 are droped
female_all = female_all[order(female_all$FDR),] # sorts the table
%%R -o MAST_raw_all

MAST_raw_all <- summaryDt_all
# save files as .csvs

MAST_raw_all.to_csv('./write/MAST_raw_new_LT.csv')
female_all.to_csv('./write/MAST_female_new_LT.csv')
%%R
# remove previous variables

rm(zlmCond_all)
rm(summaryDt_all)
rm(summaryCond_all)
rm(MAST_raw_all)
!pip list<jupyter_output>Package                Version            
---------------------- -------------------
alembic                1.3.0              
anndata                0.7.1              
anndata2ri             1.0.2              
async-generator        1.10               
attrs                  19.3.0             
backcall               0.1.0              
bleach                 3.1.0              
blinker                1.4                
certifi                2019.11.28         
certipy                0.1.3              
cffi                   1.13.2             
chardet                3.0.4              
conda                  4.7.12             
conda-package-handling 1.6.0              
cryptography           2.8                
cycler                 0.10.0             
decorator              4.4.1              
defusedxml             0.6.0              
entrypoints            0.3                
get-version            2.1                
gprofiler-official     1.0.0              
h5py       [...]
{'fraction_non_alphanumeric': 0.0988747905195116, 'fraction_numerical': 0.016040220253770648, 'mean_word_length': 3.693820224719101, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 38, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 1}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '288969', 'license': 'MIT', 'n_tokens_mistral': 2866, 'n_tokens_neox': 2647, 'n_words': 749, 'path': '06a_HSC_sexual_dimorphism_2_replicates_MAST.ipynb', 'repo_name': 'evafast/scrnaseq_paper', 'size': '8354'}