<jupyter_start><jupyter_text>Fig. 5: Benchmark FateID----In this notebook, we extract FateID's fate biases and smoothed expression.  Preliminaries  Import packages<jupyter_code># import standard packages
import numpy as np
from rpy2.robjects import ListVector

import os
import sys

# import single-cell packages
import cellrank as cr
import anndata2ri

# set verbosity levels
cr.settings.verbosity = 2

anndata2ri.activate()
%load_ext rpy2.ipython
%%R
library(RaceID)
library(FateID)
library(destiny)<jupyter_output>R[write to console]: 
Attaching package: ‘FateID’


R[write to console]: The following objects are masked from ‘package:RaceID’:

    diffexpnb, plotdiffgenesnb<jupyter_text>Print package versions for reproducibility If you want to exactly reproduce the results shown here, please make sure that your package versions match what is printed below.<jupyter_code>cr.logging.print_versions()
%%R
print(packageVersion("RaceID"))
print(packageVersion("FateID"))<jupyter_output>[1] ‘0.2.2’
[1] ‘0.1.9’<jupyter_text>Set up paths Define the paths to load data, cache results and write figure panels.<jupyter_code>sys.path.insert(0, "../../../../")  # this depends on the notebook depth and must be adapted per notebook

from paths import DATA_DIR, CACHE_DIR, FIG_DIR<jupyter_output><empty_output><jupyter_text>Set global parameters<jupyter_code>root = str(DATA_DIR / "benchmarking" / "fateid")<jupyter_output><empty_output><jupyter_text>Define utility functions<jupyter_code>%%R
# copy of plotexpressionProfile which extracts the LOESSed expression
extractData <- function (x, y, g, n, logsc = FALSE, col = NULL, name = NULL, cluster = FALSE, alpha = 0.5, lwd = 1, ylim = NULL){
    
    if ( logsc ){
        if ( min(x) == 0 ) x <- x + .1
        x <- x/apply(x,1,sum)
        x <- log2(x)
    }else{
        x <- x/apply(x,1,sum)
    }
    
    cl <- unique(y[n])
    set.seed(111111)
    if (is.null(col)) col <- sample(rainbow(length(n)))
    xlim <- c(1, length(n))
    z  <- x[g, n]
    zc <- z
    u <- 1:length(n)
       
    for ( i in 1:nrow(z) ){
        v <- as.vector(t(z[i,]))
        k <- predict(loess(v ~ u, span = alpha))
        if ( ! logsc ) k[k < 0] <- 0
        zc[i,] <- k
    }
    
    return(zc)
}

dptTraj <- function(x,y,fb,root_idxs,trthr=NULL,distance="euclidean",sigma=1000,...){
  trc <- list()
  pts <- list()  
  for ( j in colnames(fb$probs) ){
    if ( ! is.null(trthr) ){
      probs <- fb$probs
      n  <- rownames(probs)[probs[,j] > trthr]
    }else{
      votes <- fb$votes
      b <- bias(votes)
      n  <- rownames(votes)[b$bias[,j] > 1 & b$pv < .05]
    }
    dm <- destiny::DiffusionMap(as.matrix(t(x[,n])),distance=distance,sigma=sigma,...)
    root_idx <- root_idxs[[j]]
    pt <- destiny::DPT(dm, root_idx)
    pto <- pt[root_idx, ]
    
    #b <- pt@branch[, 1]
    #tip_idx <- which(b==1 & !is.na(b) & pt@tips[, 1])
    #pto <- pt[tip_idx, ]
    
    od <- order(pto,decreasing=FALSE)
    n <- n[od]

    #ts <- Transitions(as.matrix(t(x[,n])),distance=distance,sigma=sigma,...)
    #pt <- dpt::dpt(ts, branching = FALSE)
    #n <- n[order(pt$DPT,decreasing=FALSE)]
  
    if ( median((1:length(n))[y[n] == sub("t","",j)]) < median((1:length(n))[y[n] != sub("t","",j)]) ) {
        n <- rev(n)
        od <- rev(od)
    }
    trc[[j]] <- n
    pts[[j]] <- pto[od]
  }
  return(list(trc=trc,pts=pts))
}

# reason:
# Error in bias(votes) : could not find function "bias"
bias <- function(tvn){
  bias <- tvn/apply(tvn,1,function(x) x[order(x,decreasing=TRUE)][2] + 1e-3)
  pv <- apply(tvn,1,function(x){ h <- x[order(x,decreasing=TRUE)][1];  l <- x[order(x,decreasing=TRUE)][2]; binom.test( c(h, l), alternative="t" )$p.value})
  return(list(bias=bias,pv=pv))
}<jupyter_output><empty_output><jupyter_text>Load the data Load the AnnData object from the CellRank software package.<jupyter_code>adata = cr.datasets.pancreas(DATA_DIR / "pancreas" / "pancreas.h5ad")
del adata.uns["neighbors"]  # creshed anndata2ri, we don't need it
adata<jupyter_output><empty_output><jupyter_text>Analysis  Convert AnnData to SCseq<jupyter_code>%%R -i adata -i root
out_dir = root
X <- adata@assays@data[['X']]  # X, spliced, unspliced
rownames(X) <- rownames(adata)
colnames(X) <- colnames(adata)
# The raw expression data matrix with cells as columns and genes as rows in sparse matrix format.
sc <- SCseq(X)<jupyter_output><empty_output><jupyter_text>Filter, cluster and compute the embedding<jupyter_code>%%R
sc <- filterdata(sc, mintotal=10)
sc <- compdist(sc, metric="logpearson")
sc <- clustexp(sc, samp=1000)
sc <- findoutliers(sc)
sc <- comptsne(sc)<jupyter_output>Clustering k = 1,2,..., K.max (= 30): .. 
k = 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30 
done.
subset 1 
subset 2 
subset 3 
subset 4 
subset 5 
subset 6 
subset 7 
subset 8 
subset 9 
subset 10 
subset 11 
subset 12 
subset 13 
subset 14 
subset 15 
subset 16 
subset 17 
subset 18 
subset 19 
subset 20 
subset 21 
subset 22 
subset 23 
subset 24 
subset 25 
subset 26 
subset 27 
subset 28 
subset 29 
subset 30 
subset 31 
subset 32 
subset 33 
subset 34 
subset 35 
subset 36 
subset 37 
subset 38 
subset 39 
subset 40 
subset 41 
subset 42 
subset 43 
subset 44 
subset 45 
subset 46 
subset 47 
subset 48 
subset 49 
subset 50 
find outliers in cluster 1  2  3  4  5  6  7  8  9  10  11 
merging outliers 1  2  3  4  5  6  7  8  9  10  11 
determine final clustering partition 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22<jupyter_text>Extract data needed for fateBias<jupyter_code>%%R
x <- as.matrix(getfdata(sc)[sc@cluster$features,])
print(dim(x))
genes <- c("Pax4", "Pdx1", "Arx", "Peg10", "Irs4" ,"Ghrl", "Hhex", "Cd24a")
for (gene in genes) {
    if (!gene %in% rownames(x)) {
        print(gene)
        rnames <- rownames(x)
        x <- rbind(x, getfdata(sc)[gene,])
        rownames(x) <- c(rnames, gene)
    }
}
for (gene in genes) {
    if (!gene %in% rownames(x)) {
        print(gene)
    }
}
print(dim(x))

markers <- list(Beta=c("Ins1"), Alpha=c("Gcg"), Epsilon=c("Ghrl"), Delta=c("Sst"))

# cell type partition based, n - positive integer number. For each component of FMarker the expression
# of all genes is aggregated in every cell and the n top-expressing cells are extracted.
pa <- getPart(x, markers, n=50)

# A vector with a partitioning, i. e. cluster assignment for each cell
clustering <- pa$part
# A vector with the numbers of target (end point) clusters.
# Cluster 1 comprises all cells with no enrichment of marker genes
endpoints <- pa$tar
z <- sc@distances<jupyter_output>[1]  526 2531
[1] "Irs4"
[1]  527 2531<jupyter_text>Run fateBias<jupyter_code>%%R
# fateBias(x, y, tar, z=NULL, minnr=5, minnrh=10, adapt=TRUE, confidence=0.75, nbfactor=5, use.dist=FALSE, seed=12345, nbtree=NULL)
# minr - step size of the algo (test size - number of cells to be classified per end point, recommended 5)
# minnrh - cells used for training (can be inf to use all previous cells, recommended: 20)
# x - data representation (e.g. could be only HVGs, PCA, ...)

# as.matrix(x) necessary if use.dist=FALSE
fb  <- fateBias(as.matrix(x), clustering, endpoints, z=z, minnr=5, minnrh=20, seed=123, use.dist=FALSE)<jupyter_output>minnr: 5 
minnrh: 20 
test set size iteration 1 : 5 5 5 5 
randomforest iteration 1 of 20 cells
test set size iteration 2 : 5 3 3 2 
randomforest iteration 2 of 13 cells
test set size iteration 3 : 5 2.5 3.75 3.75 
randomforest iteration 3 of 15 cells
test set size iteration 4 : 5 2.5 1.666667 3.333333 
randomforest iteration 4 of 12 cells
test set size iteration 5 : 5 1.666667 0.8333333 1.666667 
randomforest iteration 5 of 10 cells
test set size iteration 6 : 5 1 1 2 
randomforest iteration 6 of 9 cells
test set size iteration 7 : 5 2.5 2.5 2.5 
randomforest iteration 7 of 11 cells
test set size iteration 8 : 5 2 1 3 
randomforest iteration 8 of 11 cells
test set size iteration 9 : 3.333333 3.333333 3.333333 5 
randomforest iteration 9 of 14 cells
test set size iteration 10 : 5 3.333333 3.333333 5 
randomforest iteration 10 of 16 cells
test set size iteration 11 : 5 1 4 3 
randomforest iteration 11 of 13 cells
test set size iteration 12 : 5 2.5 5 5 
randomforest iteration 12 of 17 ce[...]<jupyter_text>Extract the values<jupyter_code>%%R -o res
trthr = 0.15  # fate bias threshold
res = list()
for ( j in colnames(fb$probs) ){
  probs <- fb$probs
  n  <- rownames(probs)[probs[,j] > trthr]
  res[[j]] = n
}<jupyter_output><empty_output><jupyter_text>Identify root cells witin each lineage<jupyter_code>adata_preprocessed = cr.datasets.pancreas_preprocessed(DATA_DIR / "pancreas" / "pancreas_preprocessed.h5ad")
root_names = []

for i in range(2, 6):
    indices = res.rx2(f"t{i}")
    tmp = adata_preprocessed[indices].obs['dpt_pseudotime']
    ix = tmp.argmin()
    root_names.append(tmp.index[ix])
    
print(root_names)

root_ixs = []
for root_name, key in zip(root_names, [f"t{i}" for i in range(2, 6)]):
    # +1 because R indexes from 1
    root_ixs.append(np.where(res.rx2(key) == root_name)[0][0] + 1)
root_ixs = dict(zip([f"t{i}" for i in range(2, 6)], root_ixs))
print(root_ixs)
root_ixs = ListVector(root_ixs)<jupyter_output>['CATCGGGTCTTCGGTC-1-3', 'ACATACGGTCCCGACA-1-3', 'AAACGGGTCAGCTCTC-1-3', 'TCTGAGATCTGAGTGT-1-3']
{'t2': 143, 't3': 69, 't4': 7, 't5': 1701}<jupyter_text>Run DPT<jupyter_code>%%R -i root_ixs
trc <- dptTraj(as.matrix(x), clustering, fb, root_ixs, trthr=trthr, distance="euclidean")
pts <- trc[['pts']]
trc <- trc[['trc']]
names(trc) <- c('Beta', 'Alpha', 'Epsilon', 'Delta')
names(pts) <- c('Beta', 'Alpha', 'Epsilon', 'Delta')<jupyter_output><empty_output><jupyter_text>Extract and save the data for plotting  Save the fate biases, smoothed expression values and DPT<jupyter_code>%%R
for (lineage in names(trc)) {
    n <- trc[[lineage]]
    fs <- filterset(x,n=n, minexpr=0)
    for (gene in genes) {
        stopifnot(gene %in% rownames(fs))
    }
    # this is the same as plotexpressionProfile, but we return the LOESS'd data
    tmp <- extractData(as.matrix(x), clustering, g=genes, n=n, cluster=FALSE, logsc=FALSE)
    
    write.csv(t(tmp), paste0(out_dir, "/", lineage, ".csv"))
    write.csv(pts[[lineage]], paste0(out_dir, "/", lineage, "_pt.csv"))
    write.csv(trc[[lineage]], paste0(out_dir, "/", lineage, "_ixs.csv"))
}

write.csv(fb$probs, paste0(out_dir, "/", "probs.csv"))<jupyter_output><empty_output><jupyter_text>Visualize smoothed trends for the Delta lineage<jupyter_code>%%R
print(lineage)
for (gene in genes) {
    stopifnot(gene %in% rownames(x))
}
plotexpressionProfile(x, clustering, n=n, g=genes, cluster=FALSE, logsc=FALSE)<jupyter_output>[1] "Delta"
{'fraction_non_alphanumeric': 0.14017205909855995, 'fraction_numerical': 0.05030858425285207, 'mean_word_length': 3.798115746971736, 'pattern_counts': {'":': 0, '<': 131, '<?xml version=': 0, '>': 57, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '226216', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4630, 'n_tokens_neox': 4041, 'n_words': 1167, 'path': 'notebooks/fig_5_benchmarking/fateid/analysis_notebooks/MK_2020-10-17_extract_fate_bias.ipynb', 'repo_name': 'theislab/cellrank_reproducibility', 'size': '10694'}