Set-Up

Load Libraries

library(Seurat)
## Loading required package: ggplot2
## Loading required package: cowplot
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggplot2':
## 
##     ggsave
## Loading required package: Matrix
library(monocle)
## Loading required package: Biobase
## Loading required package: BiocGenerics
## Loading required package: parallel
## 
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB
## The following objects are masked from 'package:Matrix':
## 
##     colMeans, colSums, rowMeans, rowSums, which
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, basename, cbind,
##     colMeans, colnames, colSums, dirname, do.call, duplicated,
##     eval, evalq, Filter, Find, get, grep, grepl, intersect,
##     is.unsorted, lapply, lengths, Map, mapply, match, mget, order,
##     paste, pmax, pmax.int, pmin, pmin.int, Position, rank, rbind,
##     Reduce, rowMeans, rownames, rowSums, sapply, setdiff, sort,
##     table, tapply, union, unique, unsplit, which, which.max,
##     which.min
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## Loading required package: VGAM
## Loading required package: stats4
## Loading required package: splines
## Loading required package: DDRTree
## Loading required package: irlba
library(reshape2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:Biobase':
## 
##     combine
## The following objects are masked from 'package:BiocGenerics':
## 
##     combine, intersect, setdiff, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(pheatmap)
library(SingleR)
library(GSEABase)
## Loading required package: annotate
## Loading required package: AnnotationDbi
## Loading required package: IRanges
## Loading required package: S4Vectors
## 
## Attaching package: 'S4Vectors'
## The following objects are masked from 'package:dplyr':
## 
##     first, rename
## The following object is masked from 'package:Matrix':
## 
##     expand
## The following object is masked from 'package:base':
## 
##     expand.grid
## 
## Attaching package: 'IRanges'
## The following objects are masked from 'package:dplyr':
## 
##     collapse, desc, slice
## 
## Attaching package: 'AnnotationDbi'
## The following object is masked from 'package:dplyr':
## 
##     select
## Loading required package: XML
## Loading required package: graph
## 
## Attaching package: 'graph'
## The following object is masked from 'package:XML':
## 
##     addNode

Load Data

data <- readRDS("~/Downloads/All_patients_norm_final_v3.rds")
dim(data@raw.data)
## [1] 18398 24992

Here we are loading normalized data from all 3 patients, which have undergone canonical correlation analysis as described by the Sajita lab here. From this dataset, you can see we have 24,992 individual immune cells. Although we hope to characterize all the immune cells, this analysis will only look at CD8+ T cells.

CD8 Subseting

Before we start, I want to take a second to attach the annotation data to the cells, so I am quick loading a table I made of the cell type and final annotation (which I misspelled as Final.annotion cause I am not smart). We will add this data to the Seurat object using AddMetaData.

ID_final <- read.delim("~/Documents/GitHub/SC_Immune_Renal_Project/data//ID_final.txt")

meta <- as.data.frame(data@ident)
colnames(meta)[1] <- "clusters"
meta$names <- row.names(data@meta.data)
meta <- merge(meta, ID_final, by.x="clusters", by.y="Cluster.ID")
rownames(meta) <- meta$names
meta <- meta[,-2]
meta <- meta[,-1]

data <- AddMetaData(object = data, metadata = meta)

Now we can subset the Seurat object using the “CellType” variable in the meta data. I will also quickly add a new meta variable of “orig.cluster” for the sake of downstream analysis.

CD8 <- SubsetData(object =data,  subset.name = "CellType", accept.value = c("CD8"), subset.raw = T)
CD8 <-AddMetaData(CD8, CD8@ident, col.name = "orig.cluster")

Clustering Figures

All CD8 clusters included by the Type (either peripheral blood or tumor-infiltrating)

TSNEPlot(CD8, do.return = T, pt.size = 0.3, do.label = T, no.legend=T, colors.use = c("#009e73", "#99d8c7", "#bebebe", "#99d8c7", "black"))

ggsave("CD8only_TSNE.pdf", height=3, width=4)



cells.use <- data@cell.names[which(data@meta.data$CellType == "CD8")]
DimPlot(data, reduction.use = "tsne", pt.size = 0.3, do.label = F, cells.highlight =cells.use, cols.highlight = "#009e73", cols.use="#bebebe")

ggsave("CD8_highlighted.pdf", height=4, width=4.5)

Before moving on to a more comprehensive analysis, we are going to remove the data object from the environment because it is so large. This will help users that do not have much RAM.

rm(data)

Immune Signatures

Examining the immune signatures built into the SingleR package.

Manipulating the GSEA signatures to make ggplot graph.

signatures <- singler$signatures
rownames(signatures) <- singler$singler[[2]]$SingleR.single$cell.names
Sign <- CD8@meta.data

Sign <- merge(Sign, signatures, by = "row.names")
Sign2 <- Sign[,(ncol(Sign)-4):ncol(Sign)]
melted <- melt(Sign2, id.vars = c("orig.cluster"))
melted$variable  <- factor(melted$variable, levels = c("Naive", "Exhaustion", "Cytotoxicity", "Cell_Cycle"))
melted$orig.cluster  <- factor(melted$orig.cluster, levels = c("0", "3", "6", "10", "17"))
ggplot(melted, aes(x=orig.cluster, y=value)) + 
  geom_boxplot(outlier.size=0.5, lwd=0.5, aes(fill=orig.cluster)) + 
  facet_wrap(~variable) +
  guides(fill=F) + 
  scale_fill_manual(values=c("#009E73", "#99d8c7", "#bebebe","#99d8c7", "black")) +
  theme(axis.title = element_blank())

ggsave("CD8_scores_v3.pdf", height=3, width=4)

Monocle Analysis of orignial clusters

Here we are going to remove cluster 17 because it represents and activated state that is transcriptionally different. We will adjust for cell cycle later, but for now, we are looking at the basic transcriptionals structures.

CD8_v2 <- SubsetData(object =CD8,  subset.name = "orig.cluster", accept.value = c("0","3","6","10"), subset.raw = T)
monocleCD8T <- importCDS(CD8_v2, import_all = TRUE)

monocleCD8T <- estimateSizeFactors(monocleCD8T)
monocleCD8T <- suppressWarnings(estimateDispersions(monocleCD8T))
## Removing 165 outliers
monocleCD8T <- detectGenes(monocleCD8T, min_expr = 0.1)
expressed_genes2 <- row.names(subset(fData(monocleCD8T), num_cells_expressed >= 10))
diff_test_res3 <- differentialGeneTest(monocleCD8T[expressed_genes2,], fullModelFormulaStr = "~orig.cluster")
ordering_genes2 <- row.names (subset(diff_test_res3, qval < 0.01))
monocleCD8T <- setOrderingFilter(monocleCD8T, ordering_genes2)
monocleCD8T <- reduceDimension(monocleCD8T, method = 'DDRTree')
monocleCD8T <- orderCells(monocleCD8T)

plot_cell_trajectory(monocleCD8T, color_by = "State")

plot_cell_trajectory(monocleCD8T, color_by = "Type", cell_size = 0.75) + 
  scale_color_manual(values=c("#772D24","#69C3C2"))

ggsave("CD8_noGamma_byCluster_Trajectory_Type.pdf", height=4, width=4)
plot_cell_trajectory(monocleCD8T, color_by = "Final.annotion", cell_size = 0.75) + 
  scale_color_manual(values=c("#bebebe","#009E73", "#99d8c7"))

ggsave("CD8_noGamma_byCluster_Trajectory_Annot_v2.pdf", height=4, width=4)
plot_cell_trajectory(monocleCD8T, color_by = "Patient")

ggsave("CD8_noGamma_byCluster_Trajectory_Patient.pdf", height=4, width=4)
plot_cell_trajectory(monocleCD8T, color_by = "orig.cluster", cell_size = 0.75)

ggsave("CD8_noGamma_byCluster_Trajectory_cluster.pdf", height=4, width=4)
plot_cell_trajectory(monocleCD8T, color_by = "State")

Differential gene expression by the ordering of the cell trajectory.

diff_test_res_pt <- differentialGeneTest(monocleCD8T[expressed_genes2,], fullModelFormulaStr = "~sm.ns(Pseudotime)")
sig_gene_names_pt <- rownames(diff_test_res_pt[head(order(diff_test_res_pt$qval),100),])
write.table(diff_test_res_pt, file="Ajay_CD8_differentialGenes_byPseudoTime.txt", sep="\t",append=F)

Se the root state for peripheral blood.

monocleCD8T <- orderCells(monocleCD8T, root_state = 2)

Branched Analysis for Effector versus Central Memory

BEAM_res <- BEAM(monocleCD8T, branch_point = 1, cores = 2)
## Warning in if (progenitor_method == "duplicate") {: the condition has
## length > 1 and only the first element will be used
## Warning in if (progenitor_method == "sequential_split") {: the condition
## has length > 1 and only the first element will be used
BEAM_res <- BEAM_res[order(BEAM_res$qval),]

Filter out TCR genes

BEAM_res_subset <- dplyr::filter(BEAM_res, !grepl("TRBV",gene_short_name))
BEAM_res_subset <- dplyr::filter(BEAM_res_subset, !grepl("TRAV",gene_short_name))
rownames(BEAM_res_subset) <- BEAM_res_subset$gene_short_name
write.table(BEAM_res_subset, file="Ajay_CD8_BEAM.txt", sep="\t",append=F)

Define the color pallette to use.

library("colorRamps")
library("RColorBrewer")
hmcols<-colorRampPalette(rev(brewer.pal(n = 7, name ="RdYlBu")))(65)

Looking at expression of genes across the cell trajectory

plot_cell_trajectory(monocleCD8T, markers = c("LMNA","NR4A1", "REL", "FOSB", "HSPD1"), use_color_gradient=T, show_branch_points=F) + scale_color_gradient2(low="#4575B4", mid="#F9FCC1",
                     high="#D7191C", midpoint=0.25)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave("Ajay_CD8_clusterC5Trajectory_EMV2.eps", heigh=4.5, width=6)
plot_cell_trajectory(monocleCD8T, markers = c("PDCD1", "LAG3", "TNFRSF9", "HAVCR2",  "GZMK"), use_color_gradient=T, show_branch_points=F) + scale_color_gradient2(low="#4575B4", mid="#F9FCC1",
                     high="#D7191C", midpoint=0.25)
## Scale for 'colour' is already present. Adding another scale for
## 'colour', which will replace the existing scale.

ggsave("Ajay_CD8_cluster2Trajectory_CMV2.eps", heigh=4.5, width=6)

Cell Cycle Regression

We previously identified our CD8 cells using Spearman correlational analysis on the Blue+Encode and HPCA datasets. We found cluster 17 to be comprised of “gamma-delta-like” T cells, further investigation actually shows these CD8+ T cells are enriched for cell cycle progression. To reduce the affects of cell cycle on the clustering analysis, we will need to assign cell cycle state. The process and files required are again available at the Sajita Laboratory located here.

cc.genes <- readLines(con = "~/Documents/GitHub/SC_Immune_Renal_Project/data/regev_lab_cell_cycle_genes.txt")

# We can segregate this list into markers of G2/M phase and markers of S
# phase
s.genes <- cc.genes[1:43]
g2m.genes <- cc.genes[44:97]

Here is the assignment step for cell cycle state.

cellCycle <- CellCycleScoring(object = CD8, s.genes = s.genes, g2m.genes = g2m.genes, 
    set.ident = TRUE)

To visualize the difference between clusters, we are going to take the meta data variables, “Phase” and “orig.cluster” to make a proportion table and plug that into ggplots.

freq_table <- prop.table(table(cellCycle@meta.data[, "Phase"], cellCycle@meta.data[, "orig.cluster"]), 
    margin = 2)
freq_table <- as.data.frame(freq_table)

freq_table$Var1 <- factor(freq_table$Var1, levels = c("G1", "S", "G2M"))
ggplot(freq_table, aes(x=Var2, y=Freq, fill=Var1)) + 
  geom_bar(stat="identity", position="fill", color="black", lwd=0.25) + 
  theme(axis.title.x = element_blank())+ 
  guides(fill=F)+
   scale_fill_manual(values=c("#FDE723", "#22958B", "#430154"))

ggsave("CD8_cellCycleDistribution.pdf", height = 2, width=3)

We can see that Cluster 17, as suspected, has the greatest percentage of cells in S and G2M phases, implying replicating cells. In order to remove the effects of cell cycle, we are going to scale the data to regress out the difference in cell cycles.

cellCycle@meta.data$CC.Difference <- cellCycle@meta.data$S.Score - cellCycle@meta.data$G2M.Score
cellCycle <- ScaleData(object = cellCycle, vars.to.regress = "CC.Difference", display.progress = FALSE)

Now we can re-cluster the CD8 T cells and examine the changes of our regression. Here we need to:

1. Find Variable Genes

cellCycle <- FindVariableGenes(cellCycle, mean.function = ExpMean, dispersion.function = LogVMR,
                               x.low.cutoff = 0.0125, x.high.cutoff = 3, y.cutoff = 0.5)

2. Find Principal Compoenets

cellCycle <- RunPCA(object = cellCycle, pc.genes = cellCycle@var.genes, do.print = F, pcs.print = 1:5,
                    genes.print = 5)

After this we can visualize the PC using a number of different methods using Seurat

VizPCA(cellCycle, pcs.use = 10)

PCAPlot(cellCycle, dim.1 = 1, dim.2 = 3)

PCHeatmap(cellCycle, pc.use = 1:6, cells.use = 500, do.balanced = TRUE, label.columns = FALSE)

PCElbowPlot(cellCycle)

We can also perform a JackStraw analysis to attempt to determine what PC we should use.

cellCycle <- JackStraw(object = cellCycle, num.replicate = 100, display.progress = FALSE)
JackStrawPlot(object = cellCycle, PCs = 1:12)
## Warning: Removed 10896 rows containing missing values (geom_point).

## An object of class seurat in project 10X_GU0700_aggr 
##  18398 genes across 6529 samples.
ggsave("jackstraw.pdf", height=5, width=5)
## Warning: Removed 10896 rows containing missing values (geom_point).

3. Find Clusters and Run tSNE

Here we are going to use the top 10 principal components.

cellCycle <- FindClusters(cellCycle, reduction.type = 'pca', dims.use = 1:10,
                          resolution = 0.3, print.output = 0, force.recalc=T)

cellCycle <- RunTSNE(cellCycle, dims.use=1:10, check_duplicates = FALSE)
cellCycle <- AddMetaData(cellCycle, cellCycle@ident, col.name = "cc.cluster")

4. Visualize

The tSNE plots by new cluster ID and original cluster IDs.

TSNEPlot(cellCycle, pt.size= 0.3, alpha=0.75, do.label=TRUE)

ggsave("CD8_newclusters_TSNE.pdf", height=3, width=4)
TSNEPlot(cellCycle, do.label=TRUE, group.by = "orig.cluster", do.return = TRUE, pt.size= 0.3, alpha=0.75)

ggsave("CD8cc_TSNE.pdf", height=3, width=4)

In subsequent analysis we believe cluster 7 may be myeloid contaminant and eliminated it.

cellCycle <- SubsetData(object =cellCycle,  subset.name = "cc.cluster", accept.value = c("0", "1", "2", "3", "4", "5", "6"), subset.raw = T)
TSNEPlot(cellCycle, pt.size= 0.3, alpha=0.75, do.label=TRUE)

ggsave("CD8_newclusters_TSNE_no7.pdf", height=3, width=4)
TSNEPlot(cellCycle, do.label=TRUE, group.by = "orig.cluster", do.return = TRUE, pt.size= 0.3, alpha=0.75)

ggsave("CD8cc_TSNE_no7.pdf", height=3, width=4)

Relative percentage of old clusters that form the new clusters

freq_table <- prop.table(table(cellCycle@ident, cellCycle@meta.data[, "orig.cluster"]), 
    margin = 2)
freq_table <- as.data.frame(freq_table)

freq_table$Var2 <- factor(freq_table$Var2, levels = c("6", "3", "10", "0", "17"))
ggplot(freq_table, aes(x=reorder(Var1, -Var1), y=Freq, fill=Var2)) + 
  geom_bar(stat="identity", position="fill", color="black", lwd=0.25) + 
  theme(axis.title.x = element_blank()) +
  scale_fill_manual(values=c("#bebebe", "#99d8c7", "#99d8c7","#009e73", "black")) + 
  guides(fill=F)
## Warning in Ops.factor(Var1): '-' not meaningful for factors

## Warning in Ops.factor(Var1): '-' not meaningful for factors

ggsave("CD8_newclusterBreakdown.pdf", height=3, width=2)
## Warning in Ops.factor(Var1): '-' not meaningful for factors

## Warning in Ops.factor(Var1): '-' not meaningful for factors

5. New cluster IDentification and Built in Pathway Analysis

Cluster cell type identification using HPCA

SingleR.DrawHeatmap(singler$singler[[1]]$SingleR.single, top.n = 15, clusters=singler$meta.data$cc.cluster, order.by.clusters = T)

Using Blue+ENCODE datasets

SingleR.DrawHeatmap(singler$singler[[2]]$SingleR.single, top.n = 15, clusters=singler$meta.data$cc.cluster, order.by.clusters = T)

Here we can see that cluster 7 may be a contimate of myeloid cells. Which may account for the distinct transcriptional state of these cells.

Built-in signatures for individual cells

signatures <- singler$signatures
rownames(signatures) <- singler$singler[[1]]$SingleR.single$cell.names
Sign <- cellCycle@meta.data
Sign <- merge(Sign, signatures, by = "row.names")
Sign2 <- Sign[,42:46] #selecting the signatures
melted <- melt(Sign2, id.vars = c("cc.cluster"))
melted$variable  <- factor(melted$variable, levels = c("Naive", "Exhaustion", "Cytotoxicity", "Cell_Cycle"))
melted$cc.cluster  <- factor(melted$cc.cluster, levels = c("0", "1", "2", "3", "4", "5", "6", "7"))
ggplot(melted, aes(x=cc.cluster, y=value)) + 
  geom_boxplot(outlier.shape = NA, lwd=0.5, aes(fill=cc.cluster)) + 
  facet_wrap(~variable) + 
  guides(fill=F) + 
  theme(axis.title = element_blank())


Immune Signatures

Downloaded the gene signatures from Aziz, A. and Carr, A.J., et al work located here. I had to update the list manually to remove protein names and add major histocompatibility complex genes.In the end, the data frame is 20 rows of gene lists.

signature <- read.delim("~/Documents/GitHub/SC_Immune_Renal_Project/data/Immune_sign.txt")
names <- colnames(signature)

Looping through the data frame of gene lists to form Genesets for the GSEAbase R package. Because the lists are factors of unequal lengths, I had to reassign them as characters, remove empty strings, and use unique() to remove possible repeats.

Lastly before converting the lists into gene sets, I converted the genes to lower case. This is important as the SingleR R package only works with lower case genes. Running this loop is going to add 20 objects to your environment.

for (i in 1:length(names)) {
  out <- signature[,i]
  out <- as.character(out)
  out <- out[out != ""] 
  out <- unique(out)
  out <- sapply(out, tolower)
  y <- GeneSet(out, setName=paste(names[i]))
  assign(noquote(paste(names[i])),  y)
}

Making the Gene Set collection from all the new Gene Sets, listing the ones I am specifically interesting in looking at for CD8 T cells.

GeneSet <- GeneSetCollection(Antinflammatory, Proinflammatory, Lipid_mediators, Glycolysis, TCA_cycle, PPS, Glycogen_Metabolism, Glucose_Deprivation, Cytolytic, T1_Interferon, T2_Interferon, Hypoxia, CD8_Activation, Treg, Anergy, T_Cell_Terminal_Differentiation)

Preparing matrix for correlation analysis

counts=as.matrix(cellCycle@raw.data)
sc.data <- ReadSingleCellData(cellCycle@raw.data, annot=NULL)
N = colSums(counts>0)
sc.data.gl = counts
rownames(sc.data.gl) = tolower(rownames(sc.data.gl))

This is the actual calculation of the signatures across each cell. Be careful before running, this will break the calculation into 1000 gene clusters and iterate until it is finished. On my late-2015 iMac i7 4GHz station, it takes about 30 minutes, but I had to throttle down the number of cores this function used (numCores) because it would freeze my computer.

The output will me a dataframe with individual scores per cell for each of the signatures, this will process the output and attach it to the meta data of the Seurat object making a really versatile data frame known as “Sign”

signatures <- output
colnames(signatures) <- c("Antinflammatory", "Proinflammatory", "Lipid_mediators", "Glycolysis", "TCA_cycle", "PPS", "Glycogen_Metabolism", "Glucose_Deprivation", "Cytolytic", "T1_Interferon", "T2_Interferon", "Hypoxia", "CD8_Activation", "Treg", "Anergy", "T_Cell_Terminal_Differentiation")
rownames(signatures) <- colnames(sc.data.gl)
Sign <- cellCycle@meta.data
ident <- as.data.frame(cellCycle@ident)
Sign <- merge(Sign, signatures, by = "row.names")
Sign <-  merge(Sign, ident, by.x = "Row.names", by.y="row.names")

cellCycle <- AddMetaData(object = cellCycle, metadata = signatures)
pathway <- c("Antinflammatory", "Proinflammatory", "Lipid_mediators", "Glycolysis", "TCA_cycle", "PPS", "Glycogen_Metabolism", "Glucose_Deprivation", "Cytolytic", "T1_Interferon", "T2_Interferon", "Hypoxia", "CD8_Activation", "Treg", "Anergy", "T_Cell_Terminal_Differentiation")

Before we graph, we can examine the p-values for the pathway analysis using a loop of one-way ANOVAs.

out3 <- NULL
for (i in seq_along(pathway)) {
  aov <- aov(Sign[,pathway[i]] ~ cc.cluster, data=Sign)
  overall <- summary(aov)[[1]]$'Pr(>F)'[[1]]
  tukey <- TukeyHSD(aov)
  mc <- tukey$cc.cluster[,4]
  df <- c(overall, mc)
  out3 <- cbind(out3, df)
  df <- NULL
}

colnames(out3) <- pathway
rownames(out3)[1] <- "ANOVAp"

out3
##        Antinflammatory Proinflammatory Lipid_mediators    Glycolysis
## ANOVAp   9.932847e-252   5.718907e-179    1.908292e-08 5.038827e-249
## 1-0       0.000000e+00    8.740042e-01    4.395104e-01  0.000000e+00
## 2-0       3.195693e-03    0.000000e+00    9.999197e-01  0.000000e+00
## 3-0       0.000000e+00    0.000000e+00    2.233007e-02  0.000000e+00
## 4-0       7.290556e-03    8.884854e-01    6.899727e-01  8.000896e-03
## 5-0       9.671093e-06    8.558579e-07    3.299382e-04  0.000000e+00
## 6-0       1.396484e-03    0.000000e+00    8.698088e-01  9.998471e-01
## 2-1       0.000000e+00    0.000000e+00    7.675802e-01  0.000000e+00
## 3-1       2.730126e-10    0.000000e+00    7.795485e-01  3.624368e-02
## 4-1       0.000000e+00    3.916400e-01    9.999982e-01  0.000000e+00
## 5-1       0.000000e+00    2.070748e-08    1.158485e-06  0.000000e+00
## 6-1       0.000000e+00    0.000000e+00    1.903130e-01  0.000000e+00
## 3-2       0.000000e+00    0.000000e+00    9.995158e-02  0.000000e+00
## 4-2       9.942894e-01    0.000000e+00    8.651412e-01  1.642859e-01
## 5-2       2.498502e-01    0.000000e+00    3.446361e-04  0.000000e+00
## 6-2       1.713531e-08    4.797397e-01    8.011651e-01  1.035530e-06
## 4-3       0.000000e+00    0.000000e+00    9.585427e-01  0.000000e+00
## 5-3       0.000000e+00    0.000000e+00    1.688978e-08  0.000000e+00
## 6-3       0.000000e+00    0.000000e+00    1.734826e-02  0.000000e+00
## 5-4       7.813366e-01    6.347691e-03    4.265773e-05  0.000000e+00
## 6-4       7.427948e-08    0.000000e+00    2.876313e-01  6.398257e-02
## 6-5       5.757672e-11    0.000000e+00    2.756365e-01  0.000000e+00
##            TCA_cycle           PPS Glycogen_Metabolism Glucose_Deprivation
## ANOVAp 7.743935e-312 1.976487e-162         0.001627663        0.000000e+00
## 1-0     0.000000e+00  2.015193e-01         0.843340719        0.000000e+00
## 2-0     7.833204e-01  1.202077e-02         0.005772089        0.000000e+00
## 3-0     0.000000e+00  2.600988e-07         0.997290764        0.000000e+00
## 4-0     9.238655e-01  9.901154e-01         0.016183117        1.205071e-03
## 5-0     0.000000e+00  0.000000e+00         0.363870217        0.000000e+00
## 6-0     9.823085e-01  5.904576e-01         0.966539541        0.000000e+00
## 2-1     1.772360e-12  2.181112e-06         0.300005823        0.000000e+00
## 3-1     9.801773e-01  1.039941e-02         0.998506135        0.000000e+00
## 4-1     2.822135e-07  9.644932e-01         0.258289529        0.000000e+00
## 5-1     0.000000e+00  0.000000e+00         0.909889482        0.000000e+00
## 6-1     3.480721e-10  2.770934e-02         0.999999884        9.990443e-01
## 3-2     2.581824e-12  0.000000e+00         0.178294143        3.286780e-05
## 4-2     1.000000e+00  3.127673e-02         0.997789472        0.000000e+00
## 5-2     0.000000e+00  0.000000e+00         0.999581757        0.000000e+00
## 6-2     6.248863e-01  9.991978e-01         0.821048544        0.000000e+00
## 4-3     6.580841e-08  6.196673e-03         0.153449529        0.000000e+00
## 5-3     0.000000e+00  0.000000e+00         0.762192493        0.000000e+00
## 6-3     8.578727e-11  1.578276e-06         0.999218043        0.000000e+00
## 5-4     0.000000e+00  0.000000e+00         0.983124234        2.792422e-05
## 6-4     7.418694e-01  4.179644e-01         0.661430295        0.000000e+00
## 6-5     0.000000e+00  0.000000e+00         0.983473504        1.959323e-07
##           Cytolytic T1_Interferon T2_Interferon      Hypoxia
## ANOVAp 0.000000e+00  9.882478e-68  0.000000e+00 0.000000e+00
## 1-0    0.000000e+00  5.064523e-07  0.000000e+00 0.000000e+00
## 2-0    6.531159e-04  0.000000e+00  0.000000e+00 0.000000e+00
## 3-0    5.330627e-01  0.000000e+00  0.000000e+00 0.000000e+00
## 4-0    0.000000e+00  6.818919e-05  0.000000e+00 0.000000e+00
## 5-0    0.000000e+00  8.559925e-05  0.000000e+00 0.000000e+00
## 6-0    4.627416e-01  1.485748e-02  0.000000e+00 0.000000e+00
## 2-1    0.000000e+00  0.000000e+00  0.000000e+00 5.123496e-01
## 3-1    0.000000e+00  1.323221e-07  2.268058e-03 0.000000e+00
## 4-1    0.000000e+00  9.968052e-01  0.000000e+00 0.000000e+00
## 5-1    0.000000e+00  9.815973e-01  0.000000e+00 6.831628e-02
## 6-1    0.000000e+00  9.999982e-01  0.000000e+00 9.999994e-01
## 3-2    5.724703e-01  6.172010e-04  0.000000e+00 0.000000e+00
## 4-2    0.000000e+00  0.000000e+00  2.650596e-03 0.000000e+00
## 5-2    0.000000e+00  3.747080e-11  0.000000e+00 8.205876e-04
## 6-2    9.883495e-01  0.000000e+00  4.596837e-01 8.153759e-01
## 4-3    0.000000e+00  1.151413e-03  0.000000e+00 0.000000e+00
## 5-3    0.000000e+00  7.455422e-03  1.384643e-09 2.259812e-03
## 6-3    9.981033e-01  4.381684e-04  0.000000e+00 3.431635e-08
## 5-4    2.954464e-05  9.999877e-01  0.000000e+00 0.000000e+00
## 6-4    0.000000e+00  9.968419e-01  8.698547e-01 0.000000e+00
## 6-5    8.668288e-12  9.862164e-01  0.000000e+00 3.524181e-01
##        CD8_Activation         Treg       Anergy
## ANOVAp     0.00000000 2.838255e-99 0.000000e+00
## 1-0        0.00000000 0.000000e+00 0.000000e+00
## 2-0        0.00000000 0.000000e+00 0.000000e+00
## 3-0        0.00000000 0.000000e+00 0.000000e+00
## 4-0        0.07971477 6.132064e-07 0.000000e+00
## 5-0        0.00000000 0.000000e+00 0.000000e+00
## 6-0        0.00000000 9.760215e-01 0.000000e+00
## 2-1        0.00000000 4.429020e-02 0.000000e+00
## 3-1        0.00000000 2.767606e-02 0.000000e+00
## 4-1        0.00000000 0.000000e+00 0.000000e+00
## 5-1        0.00000000 9.998943e-01 9.966387e-01
## 6-1        0.00000000 4.514056e-12 1.563235e-01
## 3-2        0.00000000 2.555883e-07 1.109031e-06
## 4-2        0.00000000 0.000000e+00 0.000000e+00
## 5-2        0.13902613 5.493547e-01 0.000000e+00
## 6-2        0.80411523 4.496681e-06 0.000000e+00
## 4-3        0.00000000 0.000000e+00 0.000000e+00
## 5-3        0.00000000 1.131090e-01 0.000000e+00
## 6-3        0.00000000 0.000000e+00 0.000000e+00
## 5-4        0.00000000 0.000000e+00 3.761289e-09
## 6-4        0.00000000 3.590316e-05 1.862054e-04
## 6-5        0.97868866 1.441710e-07 6.785814e-01
##        T_Cell_Terminal_Differentiation
## ANOVAp                    0.000000e+00
## 1-0                       0.000000e+00
## 2-0                       7.032708e-12
## 3-0                       0.000000e+00
## 4-0                       4.808879e-01
## 5-0                       0.000000e+00
## 6-0                       1.748380e-02
## 2-1                       0.000000e+00
## 3-1                       0.000000e+00
## 4-1                       0.000000e+00
## 5-1                       0.000000e+00
## 6-1                       0.000000e+00
## 3-2                       0.000000e+00
## 4-2                       1.372021e-02
## 5-2                       6.604700e-04
## 6-2                       0.000000e+00
## 4-3                       0.000000e+00
## 5-3                       0.000000e+00
## 6-3                       0.000000e+00
## 5-4                       6.069863e-09
## 6-4                       6.958497e-04
## 6-5                       0.000000e+00

We can do the same staistical testing for the original clusters as well

out4 <- NULL
for (i in seq_along(pathway)) {
  aov <- aov(Sign[,pathway[i]] ~ orig.cluster, data=Sign)
  overall <- summary(aov)[[1]]$'Pr(>F)'[[1]]
  tukey <- TukeyHSD(aov)
  mc <- tukey$orig.cluster[,4]
  df <- c(overall, mc)
  out4 <- cbind(out4, df)
  df <- NULL
}

colnames(out4) <- pathway
rownames(out4)[1] <- "ANOVAp"

out4
##        Antinflammatory Proinflammatory Lipid_mediators    Glycolysis
## ANOVAp   1.877335e-258   6.842585e-128    4.379477e-19 2.304101e-237
## 3-0       4.387934e-12    4.387934e-12    8.601462e-01  4.387934e-12
## 6-0       4.387934e-12    4.387934e-12    3.677099e-02  4.387934e-12
## 10-0      4.387934e-12    6.383949e-11    8.085630e-10  4.387934e-12
## 17-0      4.387934e-12    4.387934e-12    1.482281e-11  4.387934e-12
## 6-3       2.377215e-08    9.999126e-01    1.353848e-02  7.878567e-10
## 10-3      4.436562e-12    4.387934e-12    5.443705e-09  7.016084e-08
## 17-3      9.975307e-01    4.484777e-05    2.104017e-11  4.387934e-12
## 10-6      8.231260e-06    4.387934e-12    1.713006e-03  9.993691e-01
## 17-6      2.233171e-04    4.801864e-06    2.353649e-06  4.387934e-12
## 17-10     5.288769e-12    4.387934e-12    1.429198e-01  4.387934e-12
##            TCA_cycle           PPS Glycogen_Metabolism Glucose_Deprivation
## ANOVAp 1.059573e-305 6.330583e-169         0.001529774        0.000000e+00
## 3-0     8.361564e-10  2.789049e-03         0.108367591        4.404477e-12
## 6-0     4.434786e-12  3.453727e-03         0.757570374        4.387934e-12
## 10-0    4.435452e-12  4.450884e-12         0.999312688        7.175747e-01
## 17-0    4.387934e-12  4.387934e-12         0.079958452        4.387934e-12
## 6-3     1.197474e-01  9.758248e-01         0.011420630        4.387934e-12
## 10-3    2.170112e-02  2.989549e-03         0.146646385        4.435896e-12
## 17-3    4.387934e-12  4.387934e-12         0.964438015        3.726345e-08
## 10-6    8.857964e-01  2.579970e-05         0.942421618        4.387934e-12
## 17-6    4.387934e-12  4.387934e-12         0.012884701        4.424794e-12
## 17-10   4.387934e-12  4.387934e-12         0.091823025        4.387934e-12
##            Cytolytic T1_Interferon T2_Interferon      Hypoxia
## ANOVAp 1.462819e-271  5.153466e-37  0.000000e+00 0.000000e+00
## 3-0     9.696675e-01  3.114796e-07  4.387934e-12 4.431233e-12
## 6-0     4.387934e-12  4.430567e-12  4.387934e-12 4.387934e-12
## 10-0    4.387934e-12  9.277891e-01  4.387934e-12 4.387934e-12
## 17-0    2.898982e-03  8.025805e-01  4.431677e-12 9.991666e-01
## 6-3     4.387934e-12  4.387934e-12  4.387934e-12 4.387934e-12
## 10-3    4.387934e-12  4.858709e-07  2.797405e-05 2.284135e-01
## 17-3    2.341553e-03  2.462954e-05  4.387934e-12 6.347284e-09
## 10-6    4.426570e-12  2.060541e-11  4.387934e-12 4.387934e-12
## 17-6    4.403922e-12  2.829817e-05  4.387934e-12 4.387934e-12
## 17-10   4.387934e-12  9.884305e-01  4.387934e-12 4.472755e-12
##        CD8_Activation         Treg       Anergy
## ANOVAp   0.000000e+00 4.788788e-72 0.000000e+00
## 3-0      4.387934e-12 4.387934e-12 4.446887e-12
## 6-0      4.387934e-12 4.387934e-12 4.387934e-12
## 10-0     4.387934e-12 4.387934e-12 4.429679e-12
## 17-0     4.387934e-12 3.572309e-01 1.521473e-08
## 6-3      7.684270e-10 3.405985e-01 4.387934e-12
## 10-3     6.991185e-12 8.280573e-02 9.400368e-01
## 17-3     4.494294e-12 4.494627e-12 9.999826e-01
## 10-6     4.387934e-12 8.828635e-01 4.387934e-12
## 17-6     4.387934e-12 3.566492e-11 4.387934e-12
## 17-10    1.825004e-01 2.190010e-08 9.557217e-01
##        T_Cell_Terminal_Differentiation
## ANOVAp                    0.000000e+00
## 3-0                       4.387934e-12
## 6-0                       4.387934e-12
## 10-0                      4.387934e-12
## 17-0                      4.387934e-12
## 6-3                       2.036967e-02
## 10-3                      2.574562e-03
## 17-3                      1.158961e-08
## 10-6                      8.757135e-01
## 17-6                      4.424683e-12
## 17-10                     4.445666e-12

Outputs

Loop for ggplot to make box plots for each gene signature by cluster ID.

for (i in seq_along(pathway)) {
plot <- ggplot(Sign, aes(x=`cellCycle@ident`, y=as.numeric(Sign[,pathway[i]]))) + 
  geom_boxplot(aes(fill=`cellCycle@ident`)) +
  theme(axis.title = element_blank())
plot
ggsave(paste("Boxplot_CD8_newCCCluster_", pathway[i], ".pdf", sep=""), plot, height=3, width=5)
}

Loop for ggplot to make box plots for each gene signature by Final Annotation.

for (i in seq_along(pathway)) {
plot <- ggplot(Sign, aes(x=Final.annotion, y=as.numeric(Sign[,pathway[i]]))) + 
  geom_boxplot(aes(fill=Final.annotion)) +
  theme(axis.title = element_blank())+
guides(fill=F)
plot
ggsave(paste("Boxplot_CD8_Celltype_", pathway[i], ".pdf", sep=""), plot, height=3, width=5)
}

Forming a heatmap of the gene signatures, takes some processing, first step is to remove the access data in the Sign object.

x <- Sign[,c(1, 43:59)]
colnames(x)[18] <- "cluster"
melted <- melt(x, id.vars = c("Row.names", "cluster"))

Secondly, we are going to take the mean value for each signature by cluster.

meanvalues <- melted %>%
  group_by(cluster, variable) %>%
  summarise(mean(value))

The last step in processing is going to cast the melted data frame into a matrix.

matrix <- dcast(meanvalues, cluster ~ variable)
## Using mean(value) as value column: use value.var to override.
rownames(matrix) <- matrix[,1]
matrix <- matrix[,-1]

Finally the output, with calculating the z-score across rows.

pdf("CD8_CC_PathwayHM_v2.pdf", height=4, width = 6)
pheatmap(as.matrix(t(matrix)), scale = "row", cluster_cols=F)
dev.off()
## pdf 
##   3

Monocle Analysis

Monocle is another popular single-cell pipeline that was developed by the Trapnell Laboratory.. Instead of clustering, this pipeline allows us to order single-cells into Monocle and perform the ordering.

monocle <- importCDS(cellCycle, import_all = TRUE)
monocle <- estimateSizeFactors(monocle)
monocle <- suppressWarnings(estimateDispersions(monocle))
## Removing 23 outliers
monocle <- detectGenes(monocle, min_expr = 0.1)
fData(monocle)$use_for_ordering <- fData(monocle)$num_cells_expressed > 0.05 * ncol(monocle)

monocle <- reduceDimension(monocle,
                              max_components = 2,
                              norm_method = 'log',
                              num_dim = 15,
                              reduction_method = 'tSNE',
                              verbose = T)
## Remove noise by PCA ...
## Reduce dimension by tSNE ...
monocle <- clusterCells(monocle, verbose = F)
## Distance cutoff calculated to 4.001019
HSMM_expressed_genes <-  row.names(subset(fData(monocle), num_cells_expressed >= 10))
clustering_DEG_genes <-
    differentialGeneTest(monocle[HSMM_expressed_genes,],
          fullModelFormulaStr = '~Cluster',
          cores = 2)

HSMM_ordering_genes <- row.names(subset(clustering_DEG_genes, qval < 0.01))

monocle <-
    setOrderingFilter(monocle,
        ordering_genes = HSMM_ordering_genes)

monocle <-
    reduceDimension(monocle, method = 'DDRTree')

monocle <-
    orderCells(monocle)

We can visualize the trajectory using the function plot_cell_trajectory below.

plot_cell_trajectory(monocle, color_by = "orig.cluster", cell_size = 0.75)+ 
  scale_color_manual(values=c("#009E73", "#99d8c7", "#bebebe", "#99d8c7", "black"))

ggsave("CD8_Trajectory_newCluster_byNewCluste_v3.pdf", width=4, height=4)



plot_cell_trajectory(monocle, color_by = "State", cell_size = 0.75)

ggsave("CD8_Trajectory__newCluster_byState_v2.pdf", width=4, height=4)
plot_cell_trajectory(monocle, color_by = "Final.annotion", cell_size = 0.75)

ggsave("CD8_Trajectory_newCluster_byState_v2.pdf", width=4, height=4)
plot_cell_trajectory(monocle, color_by = "CellType", cell_size = 0.75) 

plot_cell_trajectory(monocle, color_by = "cc.cluster", cell_size = 0.75) 

ggsave("CD8_Trajectory_newCluster_newCCcluster_v2.pdf", width=4, height=4.25)

Pseudotime Analysis

Using the structure we can test for significant genes as a function of pseudo-time, or the order we set in place above.

diff_test_res_pt <- differentialGeneTest(monocle[HSMM_expressed_genes,], fullModelFormulaStr = "~sm.ns(Pseudotime)")
sig_gene_names_pt <- rownames(diff_test_res_pt[head(order(diff_test_res_pt$qval),100),])
write.table(diff_test_res_pt, file="CD8_differentialGenes_byPseudoTime.txt", sep="\t",append=F)

BEAM Analysis

The last set of analyses we can perform on this CD8 T cell data using the monocle package is branched expression analysis modeling (BEAM). This let’s us look at differential gene expression between two branches within our cell trajectory.

Defining the root state to peripheral blood cells

monocle <- orderCells(monocle, root_state = 4)

Performing the differential analysis

BEAM_res <- BEAM(monocle, branch_point = 1, cores = 2)
## Warning in if (progenitor_method == "duplicate") {: the condition has
## length > 1 and only the first element will be used
## Warning in if (progenitor_method == "sequential_split") {: the condition
## has length > 1 and only the first element will be used
BEAM_res <- BEAM_res[order(BEAM_res$qval),]

Let’s filter out genes for the TCR that might be affecting results

BEAM_res_subset <- dplyr::filter(BEAM_res, !grepl("TRBV",gene_short_name))
BEAM_res_subset <- dplyr::filter(BEAM_res_subset, !grepl("TRAV",gene_short_name))
rownames(BEAM_res_subset) <- BEAM_res_subset$gene_short_name
write.table(BEAM_res_subset, file="CD8_BEAM.txt", sep="\t",append=F)

Let’s define the color scheme to match the default for pheatmap, just for the sake of my OCD.

library("colorRamps")
library("RColorBrewer")
hmcols<-colorRampPalette(rev(brewer.pal(n = 7, name ="RdYlBu")))(65)
hmcols2<-colorRampPalette(rev(brewer.pal(n = 5, name ="RdYlBu")))(65)

Now let’s draw the heatmap, here we use the filtered BEAM results, but also add the additional stipulation of p-values < 0.000001 and display genes used for ordering. The latter is looking at the role of the individual gene in the production of the cell trajectory and significantly deferentially-expressed between the branches we are comparing.

CD8Heatmap<- plot_genes_branched_heatmap(monocle[rownames(subset(BEAM_res_subset,
                                          qval < 1e-6 & use_for_ordering == TRUE)),],
                                          branch_point = 1,
                                          num_clusters = 5,
                                          cores = 2, 
                            hmcols = hmcols,
                            show_rownames=T, 
                            branch_colors = c("#bebebe","#009E73", "#99d8c7"), return_heatmap = F)

pdf("CD8_BEAM_Heatmap.pdf", height=11, width=8.5)
CD8Heatmap
## NULL
dev.off()
## quartz_off_screen 
##                 2

Session Info

sessionInfo()
## R version 3.5.1 (2018-07-02)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Sierra 10.12.6
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
##  [1] splines   stats4    parallel  stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] RColorBrewer_1.1-2   colorRamps_2.3       bindrcpp_0.2.2      
##  [4] GSEABase_1.42.0      graph_1.58.0         annotate_1.58.0     
##  [7] XML_3.98-1.16        AnnotationDbi_1.44.0 IRanges_2.14.10     
## [10] S4Vectors_0.18.3     SingleR_0.2.0        pheatmap_1.0.10     
## [13] dplyr_0.7.6          reshape2_1.4.3       monocle_2.8.0       
## [16] DDRTree_0.1.5        irlba_2.3.2          VGAM_1.0-5          
## [19] Biobase_2.40.0       BiocGenerics_0.26.0  Seurat_2.3.4        
## [22] Matrix_1.2-14        cowplot_0.9.3        ggplot2_3.0.0       
## 
## loaded via a namespace (and not attached):
##   [1] snow_0.4-2             backports_1.1.2        Hmisc_4.1-1           
##   [4] plyr_1.8.4             igraph_1.2.1           lazyeval_0.2.1        
##   [7] listenv_0.7.0          densityClust_0.3       fastICA_1.2-1         
##  [10] digest_0.6.15          foreach_1.4.4          htmltools_0.3.6       
##  [13] viridis_0.5.1          lars_1.2               gdata_2.18.0          
##  [16] magrittr_1.5           checkmate_1.8.5        memoise_1.1.0         
##  [19] cluster_2.0.7-1        mixtools_1.1.0         ROCR_1.0-7            
##  [22] limma_3.36.5           globals_0.12.1         matrixStats_0.54.0    
##  [25] R.utils_2.6.0          docopt_0.4.5           colorspace_1.3-2      
##  [28] blob_1.1.1             ggrepel_0.8.0          sparsesvd_0.1-4       
##  [31] crayon_1.3.4           RCurl_1.95-4.11        jsonlite_1.5          
##  [34] bindr_0.1.1            survival_2.42-6        zoo_1.8-3             
##  [37] iterators_1.0.10       ape_5.2                glue_1.3.0            
##  [40] gtable_0.2.0           kernlab_0.9-27         prabclus_2.2-6        
##  [43] DEoptimR_1.0-8         scales_1.0.0           mvtnorm_1.0-8         
##  [46] DBI_1.0.0              Rcpp_1.0.0             metap_0.9             
##  [49] dtw_1.20-1             viridisLite_0.3.0      xtable_1.8-2          
##  [52] htmlTable_1.12         reticulate_1.9         foreign_0.8-70        
##  [55] bit_1.1-14             proxy_0.4-22           mclust_5.4.1          
##  [58] SDMTools_1.1-221       Formula_1.2-3          GSVA_1.28.0           
##  [61] tsne_0.1-3             htmlwidgets_1.3        httr_1.3.1            
##  [64] FNN_1.1                gplots_3.0.1           fpc_2.1-11.1          
##  [67] acepack_1.4.1          modeltools_0.2-22      ica_1.0-2             
##  [70] pkgconfig_2.0.1        R.methodsS3_1.7.1      flexmix_2.3-14        
##  [73] nnet_7.3-12            labeling_0.3           later_0.7.3           
##  [76] tidyselect_0.2.4       rlang_0.2.1            pbmcapply_1.2.5       
##  [79] munsell_0.5.0          tools_3.5.1            RSQLite_2.1.1         
##  [82] ggridges_0.5.0         evaluate_0.11          stringr_1.3.1         
##  [85] yaml_2.2.0             outliers_0.14          knitr_1.20            
##  [88] bit64_0.9-7            fitdistrplus_1.0-9     robustbase_0.93-3     
##  [91] caTools_1.17.1         purrr_0.2.5            RANN_2.6              
##  [94] future_1.9.0           pbapply_1.3-4          nlme_3.1-137          
##  [97] mime_0.5               slam_0.1-43            R.oo_1.22.0           
## [100] shinythemes_1.1.1      hdf5r_1.0.0            compiler_3.5.1        
## [103] rstudioapi_0.7         png_0.1-7              geneplotter_1.58.0    
## [106] tibble_1.4.2           stringi_1.2.3          lattice_0.20-35       
## [109] trimcluster_0.1-2.1    HSMMSingleCell_0.114.0 pillar_1.3.0          
## [112] combinat_0.0-8         lmtest_0.9-36          data.table_1.11.4     
## [115] bitops_1.0-6           httpuv_1.4.5           R6_2.3.0              
## [118] latticeExtra_0.6-28    promises_1.0.1         KernSmooth_2.23-15    
## [121] gridExtra_2.3          codetools_0.2-15       MASS_7.3-50           
## [124] gtools_3.8.1           assertthat_0.2.0       rprojroot_1.3-2       
## [127] withr_2.1.2            qlcMatrix_0.9.7        diptest_0.75-7        
## [130] doSNOW_1.0.16          grid_3.5.1             rpart_4.1-13          
## [133] tidyr_0.8.1            class_7.3-14           rmarkdown_1.10        
## [136] segmented_0.5-3.0      Rtsne_0.13             shiny_1.1.0           
## [139] base64enc_0.1-3