Here we’ll use a recent dataset on chromaffin cell differentiation to look at different visualization and pseudotime techniques.

Loading data

We generally start with count matrices. These are obtained from alignment results and protocol-specific demultiplexing routines. Chromaffin dataset was measured using SMART-seq2 protocol, and the reads for each cell are then reported in a separate bam file. We use featureCounts to extract the read counts per gene:

for name in ``find ./ -name \"*.bam\" -print``; do sem -j 50 featureCounts -p -a mm10.genes.gtf -o ${name%.bam}.counts $name; done;

data.path <- '~/workshop_materials/transcriptomics/chromaffin';
counts.dir <- paste(data.path,'counts',sep='/')
str(list.files(path=counts.dir))
 chr [1:384] "A10_unique.counts" "A11_unique.counts" "A12_unique.counts" "A13_unique.counts" "A14_unique.counts" ...

We’ll define a quick function to read them in in parallel:

t.read.in.counts <- function(path,n.cores) {
  require(parallel)
  fs <- list.files(path=path,pattern=".*_unique.counts$",recursive=T)[1:6]
  names(fs) <- gsub("_unique.counts","",fs)
  names(fs) <- gsub(".+?/","",names(fs))
  dat <- do.call(cbind,mclapply(fs,function(file) { df <- read.delim(paste(path,file,sep="/"),header=F,stringsAsFactors=F,skip=2); x <- df[,7]; names(x) <- df[,1]; x;},mc.cores=n.cores))
  colnames(dat) <- paste(path,colnames(dat),sep=".");
  dat
}
dat <- t.read.in.counts(counts.dir,2)
Loading required package: parallel
str(dat)
 int [1:23420, 1:6] 0 0 0 43 26 5 0 18 0 0 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:23420] "Xkr4" "Rp1" "Sox17" "Mrpl15" ...
  ..$ : chr [1:6] "~/workshop_materials/transcriptomics/chromaffin/counts.A10" "~/workshop_materials/transcriptomics/chromaffin/counts.A11" "~/workshop_materials/transcriptomics/chromaffin/counts.A12" "~/workshop_materials/transcriptomics/chromaffin/counts.A13" ...
colnames(dat) <- gsub("~/workshop_materials/transcriptomics/chromaffin/counts.","",colnames(dat))
str(dat)
 int [1:23420, 1:6] 0 0 0 43 26 5 0 18 0 0 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:23420] "Xkr4" "Rp1" "Sox17" "Mrpl15" ...
  ..$ : chr [1:6] "A10" "A11" "A12" "A13" ...
rm(dat);

The data was the processed using scde/PAGODA:

source("pagoda.helpers.r")
library(Cairo); library(parallel)
n.cores <- 30;
min.cell.genes <- 3e3; min.cell.reads <- 1e5; min.gene.reads <- 10; min.gene.cells <- 5; min.nonfailed <- 8;
n.groups <- 10; trim <- 3;
res <- t.process.dataset(cd,nam)
res <- c(res,t.postprocess(res,name="E12_5_c",port = NULL,perplexity = 10,include.aspects = TRUE,distance.threshold = 0.9, n.clusters = 5, top.aspects = 7,return.details=T))

But let’s read in results of the processing:

suppressMessages(require(scde))
res <- readRDS(paste(data.path,"processed.rds",sep='/'))
names(res)
 [1] "cd"               "knn"              "prior"            "varinfo"          "pwpca"            "clpca"           
 [7] "batch"            "app"              "tam"              "hc"               "cl"               "clcol"           
[13] "colcols"          "tamr"             "emb"              "tamr2"            "cols"             "ct"              
[19] "ctcol"            "mat"              "matw"             "vc"               "igw"              "wgm"             
[25] "wgwm"             "ag"               "bpvt"             "btvg"             "gene.df"          "background.genes"
show.app(res$app,name='chromaffin',browse=F)
Loading required package: Rook
[1] "app loaded at: http://127.0.0.1:8787/custom/chromaffin/index.html"

Let’s recapitulate the t-SNE embedding plot:

par(mfrow=c(1,1), mar = c(0.5,0.5,0.5,0.5), mgp = c(2,0.65,0), cex = 0.9);
plot(res$emb,col=1,bg=adjustcolor(res$cols,alpha=0.6),cex=2,pch=21,panel.first=grid(),axes=F); box();

Let’s try different visualizations of the same dataset.

PCA

# get an FPM matrix as log10(FPM+1)
#fpm <- log10(exp(scde.expression.magnitude(res$knn, res$cd))+1)
lib.size <- colSums(res$cd);
fpm <- log10(t(t(res$cd)/lib.size)+1)

Run PCA:

suppressMessages(require(pcaMethods))
base.pca <- pca(t(fpm),nPcs = 2)

And visualize:

par(mfrow=c(1,1), mar = c(0.5,0.5,0.5,0.5), mgp = c(2,0.65,0), cex = 0.9);
plot(scores(base.pca),col=1,bg=adjustcolor(res$cols,alpha=0.6),cex=2,pch=21,panel.first=grid(),axes=F); box();

Let’s use PAGODA library size estimation instead:

# get an FPM matrix as log10(FPM+1)
fpm <- log10(exp(scde.expression.magnitude(res$knn, res$cd))+1)
base.pca <- pca(t(fpm),nPcs = 2)
par(mfrow=c(1,1), mar = c(0.5,0.5,0.5,0.5), mgp = c(2,0.65,0), cex = 0.9);
plot(scores(base.pca),col=1,bg=adjustcolor(res$cols,alpha=0.6),cex=2,pch=21,panel.first=grid(),axes=F); box();

We can run PCA on just a set of overdispersed genes to make things faster:

od.genes <- names(which(res$varinfo$arv>1.2))
length(od.genes)
[1] 2412
od.pca <- pca(t(fpm[od.genes,]),nPcs = 2)
par(mfrow=c(1,1), mar = c(0.5,0.5,0.5,0.5), mgp = c(2,0.65,0), cex = 0.9);
plot(scores(od.pca),col=1,bg=adjustcolor(res$cols,alpha=0.6),cex=2,pch=21,panel.first=grid(),axes=F); box();

tSNE

By default, tSNE will measure cell similarity using Euclidean distance in PCA space:

suppressMessages(require(Rtsne))
tsne <- Rtsne(t(fpm[od.genes,]), perplexity=50, verbose = FALSE, num_threads=2) 
par(mfrow=c(1,1), mar = c(0.5,0.5,0.5,0.5), mgp = c(2,0.65,0), cex = 0.9);
plot(tsne$Y,col=1,bg=adjustcolor(res$cols,alpha=0.6),cex=2,pch=21,panel.first=grid(),axes=F); box();

Another good alternative is correlation distance (on od genes, or on PCs):

d <- as.dist(1-cor(fpm[od.genes,]))
tsne2 <- Rtsne(d, is_distance = TRUE, perplexity=50, verbose = F, num_threads=2) 
par(mfrow=c(1,1), mar = c(0.5,0.5,0.5,0.5), mgp = c(2,0.65,0), cex = 0.9);
plot(tsne2$Y,col=1,bg=adjustcolor(res$cols,alpha=0.6),cex=2,pch=21,panel.first=grid(),axes=F); box();

We can see what increasing or decreasing perplexity does:

tsne3 <- Rtsne(d, is_distance = TRUE, perplexity=3, verbose = F, num_threads=2) 
par(mfrow=c(1,1), mar = c(0.5,0.5,0.5,0.5), mgp = c(2,0.65,0), cex = 0.9);
plot(tsne3$Y,col=1,bg=adjustcolor(res$cols,alpha=0.6),cex=2,pch=21,panel.first=grid(),axes=F); box();

Diffusion maps

Calculate diffusion map using overdispersed genes:

require("destiny")
Loading required package: destiny
dm <- DiffusionMap(t(fpm[od.genes,]), verbose = FALSE,k=15,sigma=100)

visualize:

par(mfrow=c(1,1), mar = c(0.5,0.5,0.5,0.5), mgp = c(2,0.65,0), cex = 0.9);
plot(dm@eigenvectors[,c(1,2)],col=1,bg=adjustcolor(res$cols,alpha=0.6),cex=2,pch=21,panel.first=grid(),axes=F); box();

Monocle2

Monocle aims to reconstruct a tree and determine embedding that best represents that tree. First we have to represent a dataset in a different way:

require("monocle")
Loading required package: monocle
Loading required package: Matrix
Loading required package: ggplot2
Loading required package: VGAM
Loading required package: stats4
Loading required package: splines
Loading required package: DDRTree
Loading required package: irlba
fpm <- res$cd[,colnames(res$mat)];
pd <- as.data.frame(cbind(res$cols ));
rownames(pd) <- colnames(fpm); colnames(pd) <- c("color")
pd = new("AnnotatedDataFrame", data = pd )
fd <- data.frame(ifelse(rownames(fpm) %in% od.genes,1,0),rownames(fpm));
rownames(fd) <- rownames(fpm); colnames(fd) <- c("od","gene_short_name")
fd = new("AnnotatedDataFrame", data = fd)
fpm1 <- newCellDataSet(fpm, phenoData = pd, featureData = fd)
fpm1 <- estimateSizeFactors(fpm1)
fpm1 <- estimateDispersions(fpm1)
Deprecated, use tibble::rownames_to_column() instead.Removing 368 outliers
fpm1 <- setOrderingFilter(fpm1,od.genes)
mncl <- reduceDimension(fpm1)
suppressWarnings(mncl <- orderCells(mncl))

Visualize the resulting tree

plot_cell_trajectory(mncl)

Check correspondence with pre-defined clusters:

plot(mncl@reducedDimS[1,],mncl@reducedDimS[2,],col=as.character(mncl@phenoData@data$col))

LS0tCnRpdGxlOiAiSGV0ZXJvZ2VuZWl0eSBhbmFseXNpcyBJSSIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKSGVyZSB3ZSdsbCB1c2UgYSByZWNlbnQgZGF0YXNldCBvbiBjaHJvbWFmZmluIGNlbGwgZGlmZmVyZW50aWF0aW9uIHRvIGxvb2sgYXQgZGlmZmVyZW50IHZpc3VhbGl6YXRpb24gYW5kIHBzZXVkb3RpbWUgdGVjaG5pcXVlcy4KCiMjIExvYWRpbmcgZGF0YQpXZSBnZW5lcmFsbHkgc3RhcnQgd2l0aCBjb3VudCBtYXRyaWNlcy4gVGhlc2UgYXJlIG9idGFpbmVkIGZyb20gYWxpZ25tZW50IHJlc3VsdHMgYW5kIHByb3RvY29sLXNwZWNpZmljIGRlbXVsdGlwbGV4aW5nIHJvdXRpbmVzLiBDaHJvbWFmZmluIGRhdGFzZXQgd2FzIG1lYXN1cmVkIHVzaW5nIFNNQVJULXNlcTIgcHJvdG9jb2wsIGFuZCB0aGUgcmVhZHMgZm9yIGVhY2ggY2VsbCBhcmUgdGhlbiByZXBvcnRlZCBpbiBhIHNlcGFyYXRlIGJhbSBmaWxlLiBXZSB1c2UgZmVhdHVyZUNvdW50cyB0byBleHRyYWN0IHRoZSByZWFkIGNvdW50cyBwZXIgZ2VuZToKCmBmb3IgbmFtZSBpbiBgYGZpbmQgLi8gLW5hbWUgXCIqLmJhbVwiIC1wcmludGBgOyBkbyBzZW0gLWogNTAgZmVhdHVyZUNvdW50cyAtcCAtYSBtbTEwLmdlbmVzLmd0ZiAtbyAke25hbWUlLmJhbX0uY291bnRzICRuYW1lOyBkb25lO2AKCgpgYGB7cn0KZGF0YS5wYXRoIDwtICd+L3dvcmtzaG9wX21hdGVyaWFscy90cmFuc2NyaXB0b21pY3MvY2hyb21hZmZpbic7CmNvdW50cy5kaXIgPC0gcGFzdGUoZGF0YS5wYXRoLCdjb3VudHMnLHNlcD0nLycpCnN0cihsaXN0LmZpbGVzKHBhdGg9Y291bnRzLmRpcikpCmBgYAoKV2UnbGwgZGVmaW5lIGEgcXVpY2sgZnVuY3Rpb24gdG8gcmVhZCB0aGVtIGluIGluIHBhcmFsbGVsOgpgYGB7cn0KdC5yZWFkLmluLmNvdW50cyA8LSBmdW5jdGlvbihwYXRoLG4uY29yZXMpIHsKICByZXF1aXJlKHBhcmFsbGVsKQogIGZzIDwtIGxpc3QuZmlsZXMocGF0aD1wYXRoLHBhdHRlcm49Ii4qX3VuaXF1ZS5jb3VudHMkIixyZWN1cnNpdmU9VClbMTo2XQogIG5hbWVzKGZzKSA8LSBnc3ViKCJfdW5pcXVlLmNvdW50cyIsIiIsZnMpCiAgbmFtZXMoZnMpIDwtIGdzdWIoIi4rPy8iLCIiLG5hbWVzKGZzKSkKICBkYXQgPC0gZG8uY2FsbChjYmluZCxtY2xhcHBseShmcyxmdW5jdGlvbihmaWxlKSB7IGRmIDwtIHJlYWQuZGVsaW0ocGFzdGUocGF0aCxmaWxlLHNlcD0iLyIpLGhlYWRlcj1GLHN0cmluZ3NBc0ZhY3RvcnM9Rixza2lwPTIpOyB4IDwtIGRmWyw3XTsgbmFtZXMoeCkgPC0gZGZbLDFdOyB4O30sbWMuY29yZXM9bi5jb3JlcykpCiAgY29sbmFtZXMoZGF0KSA8LSBwYXN0ZShwYXRoLGNvbG5hbWVzKGRhdCksc2VwPSIuIik7CiAgZGF0Cn0KYGBgCgpgYGB7cn0KZGF0IDwtIHQucmVhZC5pbi5jb3VudHMoY291bnRzLmRpciwyKQpzdHIoZGF0KQpgYGAKCmBgYHtyfQpjb2xuYW1lcyhkYXQpIDwtIGdzdWIoIn4vd29ya3Nob3BfbWF0ZXJpYWxzL3RyYW5zY3JpcHRvbWljcy9jaHJvbWFmZmluL2NvdW50cy4iLCIiLGNvbG5hbWVzKGRhdCkpCnN0cihkYXQpCnJtKGRhdCk7CmBgYAoKVGhlIGRhdGEgd2FzIHRoZSBwcm9jZXNzZWQgdXNpbmcgc2NkZS9QQUdPREE6CmBgYHtyIGV2YWw9RkFMU0V9CnNvdXJjZSgicGFnb2RhLmhlbHBlcnMuciIpCmxpYnJhcnkoQ2Fpcm8pOyBsaWJyYXJ5KHBhcmFsbGVsKQpuLmNvcmVzIDwtIDMwOwptaW4uY2VsbC5nZW5lcyA8LSAzZTM7IG1pbi5jZWxsLnJlYWRzIDwtIDFlNTsgbWluLmdlbmUucmVhZHMgPC0gMTA7IG1pbi5nZW5lLmNlbGxzIDwtIDU7IG1pbi5ub25mYWlsZWQgPC0gODsKbi5ncm91cHMgPC0gMTA7IHRyaW0gPC0gMzsKcmVzIDwtIHQucHJvY2Vzcy5kYXRhc2V0KGNkLG5hbSkKcmVzIDwtIGMocmVzLHQucG9zdHByb2Nlc3MocmVzLG5hbWU9IkUxMl81X2MiLHBvcnQgPSBOVUxMLHBlcnBsZXhpdHkgPSAxMCxpbmNsdWRlLmFzcGVjdHMgPSBUUlVFLGRpc3RhbmNlLnRocmVzaG9sZCA9IDAuOSwgbi5jbHVzdGVycyA9IDUsIHRvcC5hc3BlY3RzID0gNyxyZXR1cm4uZGV0YWlscz1UKSkKYGBgCgoKQnV0IGxldCdzIHJlYWQgaW4gcmVzdWx0cyBvZiB0aGUgcHJvY2Vzc2luZzoKYGBge3J9CnN1cHByZXNzTWVzc2FnZXMocmVxdWlyZShzY2RlKSkKcmVzIDwtIHJlYWRSRFMocGFzdGUoZGF0YS5wYXRoLCJwcm9jZXNzZWQucmRzIixzZXA9Jy8nKSkKbmFtZXMocmVzKQpgYGAKCmBgYHtyfQpzaG93LmFwcChyZXMkYXBwLG5hbWU9J2Nocm9tYWZmaW4nLGJyb3dzZT1GKQpgYGAKCkxldCdzIHJlY2FwaXR1bGF0ZSB0aGUgdC1TTkUgZW1iZWRkaW5nIHBsb3Q6CmBgYHtyfQpwYXIobWZyb3c9YygxLDEpLCBtYXIgPSBjKDAuNSwwLjUsMC41LDAuNSksIG1ncCA9IGMoMiwwLjY1LDApLCBjZXggPSAwLjkpOwpwbG90KHJlcyRlbWIsY29sPTEsYmc9YWRqdXN0Y29sb3IocmVzJGNvbHMsYWxwaGE9MC42KSxjZXg9MixwY2g9MjEscGFuZWwuZmlyc3Q9Z3JpZCgpLGF4ZXM9Rik7IGJveCgpOwpgYGAKCkxldCdzIHRyeSBkaWZmZXJlbnQgdmlzdWFsaXphdGlvbnMgb2YgdGhlIHNhbWUgZGF0YXNldC4gCgoKIyMgUENBCgpgYGB7cn0KIyBnZXQgYW4gRlBNIG1hdHJpeCBhcyBsb2cxMChGUE0rMSkKI2ZwbSA8LSBsb2cxMChleHAoc2NkZS5leHByZXNzaW9uLm1hZ25pdHVkZShyZXMka25uLCByZXMkY2QpKSsxKQpsaWIuc2l6ZSA8LSBjb2xTdW1zKHJlcyRjZCk7CmZwbSA8LSBsb2cxMCh0KHQocmVzJGNkKS9saWIuc2l6ZSkrMSkKYGBgCgpSdW4gUENBOgpgYGB7cn0Kc3VwcHJlc3NNZXNzYWdlcyhyZXF1aXJlKHBjYU1ldGhvZHMpKQpiYXNlLnBjYSA8LSBwY2EodChmcG0pLG5QY3MgPSAyKQpgYGAKCkFuZCB2aXN1YWxpemU6CmBgYHtyfQpwYXIobWZyb3c9YygxLDEpLCBtYXIgPSBjKDAuNSwwLjUsMC41LDAuNSksIG1ncCA9IGMoMiwwLjY1LDApLCBjZXggPSAwLjkpOwpwbG90KHNjb3JlcyhiYXNlLnBjYSksY29sPTEsYmc9YWRqdXN0Y29sb3IocmVzJGNvbHMsYWxwaGE9MC42KSxjZXg9MixwY2g9MjEscGFuZWwuZmlyc3Q9Z3JpZCgpLGF4ZXM9Rik7IGJveCgpOwpgYGAKCkxldCdzIHVzZSBQQUdPREEgbGlicmFyeSBzaXplIGVzdGltYXRpb24gaW5zdGVhZDoKYGBge3J9CiMgZ2V0IGFuIEZQTSBtYXRyaXggYXMgbG9nMTAoRlBNKzEpCmZwbSA8LSBsb2cxMChleHAoc2NkZS5leHByZXNzaW9uLm1hZ25pdHVkZShyZXMka25uLCByZXMkY2QpKSsxKQpiYXNlLnBjYSA8LSBwY2EodChmcG0pLG5QY3MgPSAyKQpwYXIobWZyb3c9YygxLDEpLCBtYXIgPSBjKDAuNSwwLjUsMC41LDAuNSksIG1ncCA9IGMoMiwwLjY1LDApLCBjZXggPSAwLjkpOwpwbG90KHNjb3JlcyhiYXNlLnBjYSksY29sPTEsYmc9YWRqdXN0Y29sb3IocmVzJGNvbHMsYWxwaGE9MC42KSxjZXg9MixwY2g9MjEscGFuZWwuZmlyc3Q9Z3JpZCgpLGF4ZXM9Rik7IGJveCgpOwpgYGAKCldlIGNhbiBydW4gUENBIG9uIGp1c3QgYSBzZXQgb2Ygb3ZlcmRpc3BlcnNlZCBnZW5lcyB0byBtYWtlIHRoaW5ncyBmYXN0ZXI6CmBgYHtyfQpvZC5nZW5lcyA8LSBuYW1lcyh3aGljaChyZXMkdmFyaW5mbyRhcnY+MS4yKSkKbGVuZ3RoKG9kLmdlbmVzKQpvZC5wY2EgPC0gcGNhKHQoZnBtW29kLmdlbmVzLF0pLG5QY3MgPSAyKQpwYXIobWZyb3c9YygxLDEpLCBtYXIgPSBjKDAuNSwwLjUsMC41LDAuNSksIG1ncCA9IGMoMiwwLjY1LDApLCBjZXggPSAwLjkpOwpwbG90KHNjb3JlcyhvZC5wY2EpLGNvbD0xLGJnPWFkanVzdGNvbG9yKHJlcyRjb2xzLGFscGhhPTAuNiksY2V4PTIscGNoPTIxLHBhbmVsLmZpcnN0PWdyaWQoKSxheGVzPUYpOyBib3goKTsKYGBgCgoKIyMgdFNORSAKCkJ5IGRlZmF1bHQsIHRTTkUgd2lsbCBtZWFzdXJlIGNlbGwgc2ltaWxhcml0eSB1c2luZyBFdWNsaWRlYW4gZGlzdGFuY2UgaW4gUENBIHNwYWNlOgpgYGB7cn0Kc3VwcHJlc3NNZXNzYWdlcyhyZXF1aXJlKFJ0c25lKSkKdHNuZSA8LSBSdHNuZSh0KGZwbVtvZC5nZW5lcyxdKSwgcGVycGxleGl0eT01MCwgdmVyYm9zZSA9IEZBTFNFLCBudW1fdGhyZWFkcz0yKSAKYGBgCgpgYGB7cn0KcGFyKG1mcm93PWMoMSwxKSwgbWFyID0gYygwLjUsMC41LDAuNSwwLjUpLCBtZ3AgPSBjKDIsMC42NSwwKSwgY2V4ID0gMC45KTsKcGxvdCh0c25lJFksY29sPTEsYmc9YWRqdXN0Y29sb3IocmVzJGNvbHMsYWxwaGE9MC42KSxjZXg9MixwY2g9MjEscGFuZWwuZmlyc3Q9Z3JpZCgpLGF4ZXM9Rik7IGJveCgpOwpgYGAKCkFub3RoZXIgZ29vZCBhbHRlcm5hdGl2ZSBpcyBjb3JyZWxhdGlvbiBkaXN0YW5jZSAob24gb2QgZ2VuZXMsIG9yIG9uIFBDcyk6CmBgYHtyfQpkIDwtIGFzLmRpc3QoMS1jb3IoZnBtW29kLmdlbmVzLF0pKQpgYGAKCmBgYHtyfQp0c25lMiA8LSBSdHNuZShkLCBpc19kaXN0YW5jZSA9IFRSVUUsIHBlcnBsZXhpdHk9NTAsIHZlcmJvc2UgPSBGLCBudW1fdGhyZWFkcz0yKSAKYGBgCgpgYGB7cn0KcGFyKG1mcm93PWMoMSwxKSwgbWFyID0gYygwLjUsMC41LDAuNSwwLjUpLCBtZ3AgPSBjKDIsMC42NSwwKSwgY2V4ID0gMC45KTsKcGxvdCh0c25lMiRZLGNvbD0xLGJnPWFkanVzdGNvbG9yKHJlcyRjb2xzLGFscGhhPTAuNiksY2V4PTIscGNoPTIxLHBhbmVsLmZpcnN0PWdyaWQoKSxheGVzPUYpOyBib3goKTsKYGBgCgpXZSBjYW4gc2VlIHdoYXQgaW5jcmVhc2luZyBvciBkZWNyZWFzaW5nIHBlcnBsZXhpdHkgZG9lczoKYGBge3J9CnRzbmUzIDwtIFJ0c25lKGQsIGlzX2Rpc3RhbmNlID0gVFJVRSwgcGVycGxleGl0eT0zLCB2ZXJib3NlID0gRiwgbnVtX3RocmVhZHM9MikgCmBgYAoKYGBge3J9CnBhcihtZnJvdz1jKDEsMSksIG1hciA9IGMoMC41LDAuNSwwLjUsMC41KSwgbWdwID0gYygyLDAuNjUsMCksIGNleCA9IDAuOSk7CnBsb3QodHNuZTMkWSxjb2w9MSxiZz1hZGp1c3Rjb2xvcihyZXMkY29scyxhbHBoYT0wLjYpLGNleD0yLHBjaD0yMSxwYW5lbC5maXJzdD1ncmlkKCksYXhlcz1GKTsgYm94KCk7CmBgYAoKIyMgRGlmZnVzaW9uIG1hcHMKCkNhbGN1bGF0ZSBkaWZmdXNpb24gbWFwIHVzaW5nIG92ZXJkaXNwZXJzZWQgZ2VuZXM6CmBgYHtyfQpyZXF1aXJlKCJkZXN0aW55IikKZG0gPC0gRGlmZnVzaW9uTWFwKHQoZnBtW29kLmdlbmVzLF0pLCB2ZXJib3NlID0gRkFMU0Usaz0xNSxzaWdtYT0xMDApCmBgYAoKdmlzdWFsaXplOgpgYGB7cn0KcGFyKG1mcm93PWMoMSwxKSwgbWFyID0gYygwLjUsMC41LDAuNSwwLjUpLCBtZ3AgPSBjKDIsMC42NSwwKSwgY2V4ID0gMC45KTsKcGxvdChkbUBlaWdlbnZlY3RvcnNbLGMoMSwyKV0sY29sPTEsYmc9YWRqdXN0Y29sb3IocmVzJGNvbHMsYWxwaGE9MC42KSxjZXg9MixwY2g9MjEscGFuZWwuZmlyc3Q9Z3JpZCgpLGF4ZXM9Rik7IGJveCgpOwpgYGAKCgojIyBNb25vY2xlMgoKTW9ub2NsZSBhaW1zIHRvIHJlY29uc3RydWN0IGEgdHJlZSBhbmQgZGV0ZXJtaW5lIGVtYmVkZGluZyB0aGF0IGJlc3QgcmVwcmVzZW50cyB0aGF0IHRyZWUuCkZpcnN0IHdlIGhhdmUgdG8gcmVwcmVzZW50IGEgZGF0YXNldCBpbiBhIGRpZmZlcmVudCB3YXk6CmBgYHtyfQpyZXF1aXJlKCJtb25vY2xlIikKZnBtIDwtIHJlcyRjZFssY29sbmFtZXMocmVzJG1hdCldOwpwZCA8LSBhcy5kYXRhLmZyYW1lKGNiaW5kKHJlcyRjb2xzICkpOwpyb3duYW1lcyhwZCkgPC0gY29sbmFtZXMoZnBtKTsgY29sbmFtZXMocGQpIDwtIGMoImNvbG9yIikKcGQgPSBuZXcoIkFubm90YXRlZERhdGFGcmFtZSIsIGRhdGEgPSBwZCApCgpmZCA8LSBkYXRhLmZyYW1lKGlmZWxzZShyb3duYW1lcyhmcG0pICVpbiUgb2QuZ2VuZXMsMSwwKSxyb3duYW1lcyhmcG0pKTsKcm93bmFtZXMoZmQpIDwtIHJvd25hbWVzKGZwbSk7IGNvbG5hbWVzKGZkKSA8LSBjKCJvZCIsImdlbmVfc2hvcnRfbmFtZSIpCmZkID0gbmV3KCJBbm5vdGF0ZWREYXRhRnJhbWUiLCBkYXRhID0gZmQpCgpmcG0xIDwtIG5ld0NlbGxEYXRhU2V0KGZwbSwgcGhlbm9EYXRhID0gcGQsIGZlYXR1cmVEYXRhID0gZmQpCmBgYAoKYGBge3J9CmZwbTEgPC0gZXN0aW1hdGVTaXplRmFjdG9ycyhmcG0xKQpmcG0xIDwtIGVzdGltYXRlRGlzcGVyc2lvbnMoZnBtMSkKZnBtMSA8LSBzZXRPcmRlcmluZ0ZpbHRlcihmcG0xLG9kLmdlbmVzKQptbmNsIDwtIHJlZHVjZURpbWVuc2lvbihmcG0xKQpzdXBwcmVzc1dhcm5pbmdzKG1uY2wgPC0gb3JkZXJDZWxscyhtbmNsKSkKYGBgCgpWaXN1YWxpemUgdGhlIHJlc3VsdGluZyB0cmVlCmBgYHtyfQpwbG90X2NlbGxfdHJhamVjdG9yeShtbmNsKQpgYGAKCkNoZWNrIGNvcnJlc3BvbmRlbmNlIHdpdGggcHJlLWRlZmluZWQgY2x1c3RlcnM6CmBgYHtyfQpwbG90KG1uY2xAcmVkdWNlZERpbVNbMSxdLG1uY2xAcmVkdWNlZERpbVNbMixdLGNvbD1hcy5jaGFyYWN0ZXIobW5jbEBwaGVub0RhdGFAZGF0YSRjb2wpKQpgYGAKCg==