In transcriptomics applications, one of the most utilized exploratory plots is the multi-dimensional scaling (MDS) plot or a principal component analysis (PCA) plot. Such plots show similarities between samples measured in an unsupervised way and give a sense of how much differential expression can be detected before conducting any formal tests. An MDS plot can be generated with the plotMDS function from the limma package. In transcriptomics, distances between samples are calculated based on the expression of the top varying genes. We propose a similar plot for HDCyto data using median marker expression over all cells to calculate dissimilarities between samples (other aggregations are also possible, and one could reduce the number of top varying markers to include in the calculation). Ideally, samples should cluster well within the same condition, although this depends on the magnitude of the difference between experimental conditions. With this diagnostic, one can identify the outlier samples and eliminate them if the circumstances warrant it. In our MDS plot on median marker expression values (see Figure 3), we can see that the first dimension (MDS1) separates the unstimulated and stimulated samples reasonably well. The second dimension (MDS2) represents, to some degree, differences between patients. Most of the samples that originate from the same patient are placed at a similar point along the y-axis, for example, samples from patients 7, 5, and 8 are at the top of the plot, samples from patient 4 are located at the bottom of the plot. This also indicates that the marker expression of individual patients is driving similarity and perhaps should be formally accounted for in the downstream statistical modeling.

library (dplyr) # Get the median marker expression per sample expr_median_sample_tbl <- data.frame ( sample_id = sample_ids, expr) %>% group_by (sample_id) %>% summarize_all ( funs (median)) expr_median_sample <- t (expr_median_sample_tbl[, - 1 ]) colnames(expr_median_sample) <- expr_median_sample_tbl $ sample_id library (limma) mds <- plotMDS (expr_median_sample, plot = FALSE ) library (ggrepel) ggdf <- data.frame( MDS1 = mds $ x, MDS2 = mds $ y, sample_id = colnames (expr_median_sample)) mm <- match (ggdf $ sample_id, md $ sample_id) ggdf $condition <- md $ condition[mm] ggplot (ggdf, aes ( x = MDS1, y = MDS2, color = condition)) + geom_point( size = 2, alpha = 0.8) + geom_label_repel( aes( label = sample_id)) + theme_bw () + scale_color_manual ( values = color_conditions) + coord_fixed ()

In contrast to genomic applications, the number of variables measured for each sample is much lower in HDCyto data. In the former, thousands of genes are surveyed, whereas in the latter, ~20–50 antigens are typically targeted. Similar to the MDS plot above, a heatmap of the same data also gives insight into the structure of the data. The heatmap shows median marker intensities with clustered columns (samples) and rows (markers). We have used hierarchical clustering with average linkage and euclidean distance, but also Ward’s linkage could be used ( Bruggner et al., 2014), and in CyTOF applications, a cosine distance shows good performance ( Bendall et al., 2014). In this plot, we can see which markers drive the observed clustering of samples (see Figure 4).

As with the MDS plot, the dendrogram separates the reference and stimulated samples very well. Also, similar groupings of patients within experimental conditions are observed (patients 1-2, 5-7-8 and 3-4 are together in both conditions).

library (RColorBrewer) library (pheatmap) # Column annotation for the heatmap mm <- match ( colnames (expr_median_sample), md $ sample_id) annotation_col <- data.frame ( condition = md $ condition[mm], row.names = colnames (expr_median_sample)) annotation_colors <- list ( condition = color_conditions) # Colors for the heatmap color <- colorRampPalette ( brewer.pal ( n = 9 , name = "YlGnBu" ))( 100 ) pheatmap (expr_median_sample, color = color, display_numbers = TRUE , number_color = "black" , fontsize_number = 5 , annotation_col = annotation_col, annotation_colors = annotation_colors, clustering_method = "average")

One of the most popular plots for representing single cell data are t-SNE plots, where each cell is represented in a lower, usually two-dimensional, space computed using t-stochastic neighbor embedding (t-SNE) ( Van Der Maaten & Hinton, 2008). More generally, dimensionality reduction techniques represent the similarity of points in 2 or 3 dimensions, such that similar objects in high dimensional space are also similar in lower dimensional space. Mathematically, there are a myriad of ways to define this similarity. For example, principal components analysis (PCA) uses linear combinations of the original features to find orthogonal dimensions that show the highest levels of variability; the top 2 or 3 principal components can then be visualized.

Nevertheless, there are few notes of caution when using t-SNE or any other dimensionality reduction technique. Since they are based on preserving similarities between cells, those that are similar in the original space will be close in the 2D/3D representation, but the opposite does not always hold. In our experience, t-SNE with default parameters for HDCyto data is often suitable; for more guidance on the specifics of t-SNE, see How to Use t-SNE Effectively ( Wattenberg et al., 2016). Due to the stochastic nature of t-SNE optimization, rerunning the method will result in different lower dimensional projections, thus it is advisable to run it a few times to identify the common trends and get a feeling about the variability of the results. As with other methods, to be sure that the analysis is reproducible, the user can define the random seed.

t-SNE is a method that requires significant computational time to process the data even for tens of thousands of cells. CyTOF datasets are usually much larger and thus to keep running times reasonable, one may use a subset of cells; for example, here we use 2000 cells from each sample. The t-SNE map below is colored according to the expression level of the CD4 marker, highlighting that the CD4+ T-cells are placed to the left side of the plot (see Figure 9). In this way, one can use a collection of markers to highlight where cell types of interest are located on the map.

Instead of t-SNE, one could also use other dimension reduction techniques, such as PCA, diffusion maps, SIMLR ( Wang et al., 2017) or isomaps, some of which are conveniently available via the cytof_dimReduction function from the cytofkit package ( Chen et al., 2016). To speed up the t-SNE analysis, one could use a multicore version that is available via the Rtsne.multicore package. Alternative algorithms, such as largeVis ( Tang et al., 2016) (available via the largeVis package), can be used for dimensionality reduction of very large datasets without downsampling. Alternatively, the dimensionality reduction can be performed on the codes of the SOM, at a resolution (size of the SOM) specified by the user (see Figure 13).

## Find and skip duplicates dups <- which ( ! duplicated (expr[, lineage_markers])) ## Data subsampling: create indices by sample inds <- split ( 1 : length (sample_ids), sample_ids) ## How many cells to downsample per-sample tsne_ncells <- pmin ( table (sample_ids), 2000 ) ## Get subsampled indices set.seed ( 1234 ) tsne_inds <- lapply ( names (inds), function (i){ s <- sample (inds[[i]], tsne_ncells[i], replace = FALSE ) intersect (s, dups) }) tsne_inds <- unlist (tsne_inds) tsne_expr <- expr[tsne_inds, lineage_markers]

## Run t-SNE library (Rtsne) set.seed ( 1234 ) tsne_out <- Rtsne (tsne_expr, check_duplicates = FALSE , pca = FALSE )

## Plot t-SNE colored by CD4 expression dr <- data.frame ( tSNE1 = tsne_out $ Y[, 1 ], tSNE2 = tsne_out $ Y[, 2 ], expr[tsne_inds, lineage_markers]) ggplot (dr, aes ( x = tSNE1, y = tSNE2, color = CD4)) + geom_point ( size = 0.8 ) + theme_bw () + scale_color_gradientn ( "CD4" , colours = colorRampPalette ( rev ( brewer.pal ( n = 11 , name = "Spectral" )))( 50 ))

We can color the cells by cluster. Ideally, cells of the same color should be close to each other (see Figure 10). When the plots are further stratified by sample (see Figure 11), we can verify whether similar cell populations are present in all replicates, which can help in identifying outlying samples. Optionally, stratification can be done by condition (see Figure 12). With such a spot-check plot, we can inspect whether differences in cell abundance are strong between conditions, and we can visualize and identify distinguishing clusters before applying formal statistical testing. A similar approach of data exploration was proposed in studies of treatment-specific differences of polyfunctional antigen-specific T-cells ( Lin et al., 2015).

dr $ sample_id <- sample_ids[tsne_inds] mm <- match (dr $ sample_id, md $ sample_id) dr $ condition <- md $ condition[mm] dr $ cell_clustering1 <- factor (cell_clustering1[tsne_inds], levels = 1 : nmc) ## Plot t-SNE colored by clusters ggp <- ggplot (dr, aes ( x = tSNE1, y = tSNE2, color = cell_clustering1)) + geom_point ( size = 0.8 ) + theme_bw () + scale_color_manual ( values = color_clusters) + guides ( color = guide_legend ( override.aes = list ( size = 4 ), ncol = 2 )) ggp

## Facet per sample ggp + facet_wrap ( ~ sample_id)

## Facet per condition ggp + facet_wrap ( ~ condition)

The SOM codes represent characteristics of the 100 (by default) clusters generated in the first step of the FlowSOM pipeline. Their visualization can also be helpful in understanding the cell population structure and determining the number of clusters. Ultimately, the metaclustering step uses the codes and not the original cells. We treat the codes as new representative cells and apply the t-SNE dimension reduction to visualize them in 2D (see Figure 13). The size of the points corresponds to the number of cells that were assigned to a given code. The points are colored according to the results of metaclustering. Since we have only 100 data points, the t-SNE analysis is fast.

As there are multiple ways to mathematically define similarity in high dimensional space, it is always good practice visualizing projections from other methods to see how consistent the observed patterns are. For instance, we also represent the FlowSOM codes via the first two principal components (see Figure 13).

## Get code sizes; sometimes not all the codes have mapped cells so they will have size 0 code_sizes <- table ( factor (som $ map $ mapping[, 1 ], levels = 1 : nrow (codes))) code_sizes <- as.numeric (code_sizes)

## Run t-SNE on codes set.seed ( 1234 ) tsne_out <- Rtsne (codes, perplexity = 5, pca = FALSE ) ## Run PCA on codes pca_out <- prcomp (codes, center = TRUE , scale. = FALSE )

codes_dr <- data.frame ( tSNE1 = tsne_out $ Y[, 1 ], tSNE2 = tsne_out $ Y[, 2 ], PCA1 = pca_out $ x[, 1 ], PCA2 = pca_out $ x[, 2 ]) codes_dr $ code_clustering1 <- factor (code_clustering1) codes_dr $ size <- code_sizes ## Plot t-SNE on codes gg_tsne_codes <- ggplot (codes_dr, aes ( x = tSNE1, y = tSNE2, color = code_clustering1, size = size)) + geom_point ( alpha = 0.9 ) + theme_bw () + scale_color_manual ( values = color_clusters) + guides ( color = guide_legend ( override.aes = list ( size = 4 ), ncol = 2 ))

## Plot PCA on codes gg_pca_codes <- ggplot (codes_dr, aes ( x = PCA1, y = PCA2, color = code_clustering1, size = size)) + geom_point ( alpha = 0.9 ) + theme_bw () + scale_color_manual ( values = color_clusters) + guides ( color = guide_legend ( override.aes = list ( size = 4 ), ncol = 2 )) + theme ( legend.position = "right" , legend.box = "vertical" ) library (cowplot) legend <- get_legend (gg_tsne_codes) ggdraw () + draw_plot (gg_tsne_codes + theme ( legend.position = "none" ), 0 , . 5 , . 7 , . 5 ) + draw_plot (gg_pca_codes + theme ( legend.position = "none" ), 0 , 0 , . 7 , . 5 ) + draw_plot (legend, . 7 , . 0 , . 2 , 1 ) + draw_plot_label ( c ( "A" , "B" , "" ), c ( 0 , 0 , . 7 ), c ( 1 , . 5 , 1 ), size = 15 )

Using heatmaps, we can also visualize median marker expression in the 100 SOM codes as in Figure 14. Of note, the clustering presented with the dendrogram does not completely agree with the clustering depicted by the 20 colors because the first one is based on the hierarchical clustering with average linkage and Euclidean distance, while the second one results from the consensus clustering.

plot_clustering_heatmap_wrapper2 ( expr = expr, expr01 = expr01, lineage_markers = lineage_markers, functional_markers = "pS6" , sample_ids = sample_ids, cell_clustering = cell_clustering_som, color_clusters = rep (color_clusters, length.out = 100 ), cluster_merging = data.frame( original_cluster = 1 : 100 , new_cluster = code_clustering1), plot_cluster_annotation = FALSE )

In our experience, manual merging of clusters leads to slightly different results compared to an algorithm with a specified number of clusters. In order to detect somewhat rare populations, some level of over-clustering is necessary so that the more subtle populations become separated from the main populations. In addition, merging can always follow an over-clustering step, but splitting of existing clusters is not generally feasible. In our setup, over-clustering is also useful when the interest is identifying the “natural” number of clusters present in the data. In addition to the t-SNE plots, one could investigate the delta area plot from the ConsensusClusterPlus package and the hierarchical clustering dendrogram of the over-clustered subpopulations, as shown in Figure 16 and Figure 18, respectively.

In our example, we expect around 6 specific cell types, and we have performed FlowSOM clustering into 20 groups as a reasonable over-estimate. After analyzing the heatmaps ( Figure 6) and t-SNE plots ( Figure 10), we can clearly see that stratification of the data into 20 clusters may be too strong. In the t-SNE map, many clusters are placed very close to each other, indicating that they could be merged together. The same can be deduced from the heatmaps, highlighting that marker expression patterns for some neighboring clusters are very similar. Cluster merging and annotating is somewhat manual, based partially on visual inspection of t-SNE plots and heatmaps and thus, benefits from expert knowledge of the cell types.

Manual cluster merging and annotating based on heatmaps. In our experience, the main reference for manual merging of clusters is the heatmap of marker characteristics across metaclusters (e.g. Figure 6), with dendrograms showing the hierarchy of similarities. Such plots aggregate information over many cells and thus show average marker expression for each cluster. Together with dimensionality reduction, these plots give good insight into the relationships between clusters and the marker levels within each cluster. Given expert knowledge of the cell types and markers, it is then left to the researcher to decide how exactly to merge clusters (e.g. with higher weight given to some markers).

The dendrogram highlights the similarity between the metaclusters and can be used explicitly for the merging. However, there are reasons why we would not always follow the dendrogram exactly. In general, when it comes to clustering, blindly following the hierarchy of codes will lead to identification of populations of similar cells, but it does not necessarily mean that they are of biological interest. The distances between metaclusters are calculated across all the markers, and it may be that some markers carry higher weight for certain cell types. In addition, different linkage methods may lead to different hierarchy, especially when clusters are not fully distinct. Another aspect to consider in cluster merging is the cluster size, represented in the parentheses next to the cluster label in our plots. If the cluster size is very small, but the cluster seems relevant and distinct, one can keep it as separate. However, if it is small and different from the neighboring clustering only in a somewhat unimportant marker, it could be merged. And, if some of the metaclusters do not represent any specific cell types, they could be dropped out of the downstream analysis instead of being merged. However, in case an automated solution for cluster merging is truly needed, one could use the cutree() function applied to the dendrogram.

Based on the seed that was set, cluster merging of the 20 metaclusters is defined in the PBMC8_cluster_merging1.xlsx file on the Robinson Lab server with the IDs of the original clusters and new cluster names, and we save it as a cluster_merging1 data frame. The expert has annotated 8 cell populations: CD8 T-cells, CD4 T-cells, B-cells IgM-, B-cells IgM+, NK cells, dendritic cells (DCs), monocytes and surface negative cells; monocytes could be further subdivided based on HLA-DR into high, medium and low subtypes.

cluster_merging1_filename <- "PBMC8_cluster_merging1.xlsx" download.file ( paste0 (url, "/" , cluster_merging1_filename), destfile = cluster_merging1_filename, mode = "wb" ) cluster_merging1 <- read_excel (cluster_merging1_filename) data.frame (cluster_merging1)

## original_cluster new_cluster ## 1 1 B-cells IgM+ ## 2 2 surface- ## 3 3 NK cells ## 4 4 CD8 T-cells ## 5 5 B-cells IgM- ## 6 6 monocytes ## 7 7 monocytes ## 8 8 CD8 T-cells ## 9 9 CD8 T-cells ## 10 10 monocytes ## 11 11 monocytes ## 12 12 CD4 T-cells ## 13 13 DC ## 14 14 CD8 T-cells ## 15 15 CD4 T-cells ## 16 16 DC ## 17 17 CD4 T-cells ## 18 18 CD4 T-cells ## 19 19 CD4 T-cells ## 20 20 CD4 T-cells

## Convert to factor with merged clusters in desired order levels_clusters_merged <- c ( "B-cells IgM+" , "B-cells IgM-" , "CD4 T-cells" , "CD8 T-cells" , "DC" , "NK cells" , "monocytes" , "surface-" ) cluster_merging1 $ new_cluster <- factor (cluster_merging1 $ new_cluster, levels = levels_clusters_merged) ## New clustering1m mm <- match (cell_clustering1, cluster_merging1 $ original_cluster) cell_clustering1m <- cluster_merging1 $ new_cluster[mm] mm <- match (code_clustering1, cluster_merging1 $ original_cluster) code_clustering1m <- cluster_merging1 $ new_cluster[mm]

We update the t-SNE plot with the new annotated cell populations, Figure 15.

dr $ cell_clustering1m <- cell_clustering1m[tsne_inds] ggplot (dr, aes ( x = tSNE1, y = tSNE2, color = cell_clustering1m)) + geom_point ( size = 0.8 ) + theme_bw () + scale_color_manual ( values = color_clusters) + guides ( color = guide_legend ( override.aes = list ( size = 4 )))

One of the useful representations of merging is a heatmap of median marker expression in each of the original clusters, which are labeled according to the proposed merging, Figure 16.

plot_clustering_heatmap_wrapper ( expr = expr[, lineage_markers_ord], expr01 = expr01[, lineage_markers_ord], cell_clustering = cell_clustering1, color_clusters = color_clusters, cluster_merging = cluster_merging1)

To get a final summary of the annotated cell types, one can plot a heatmap of median marker expression, calculated based on the cells in each of the annotated populations, Figure 17.

plot_clustering_heatmap_wrapper ( expr = expr[, lineage_markers_ord], expr01 = expr01[, lineage_markers_ord], cell_clustering = cell_clustering1m, color_clusters = color_clusters)

Reducing the number of clusters in ConsensusClusterPlus. The ConsensusClusterPlus package provides visualizations that can help to understand the metaclustering process and the characteristics of the analyzed data. For example, the delta area plot (see Figure 18) highlights the amount of extra cluster stability gained when clustering into k groups as compared to k-1 groups (from k=2 to k=20). It can be expected that high stability of clusters can be reached when clustering into the number of groups that best fits the data. Thus, using the delta area plot could help finding the “natural” number of clusters present in the data, which would correspond to the value of k where there is no appreciable increase in stability. This strategy can be referred as the “elbow criterion”. For more details about the meaning of this plot, the user can refer to the original description of the consensus clustering method ( Monti et al., 2003).

The elbow criterion is quite subjective since the “appreciable” increase is not defined exactly. For example, in the delta plot below, we could say that the last point before plateau is for k=6, or for k=5, or for k=3, depending on our perception of sufficient decrease of the delta area score. Moreover, the exact point where a plateau is reached may vary for runs with different random seeds, the drop may not always be so sharp and and the function is not guaranteed to be decreasing. It is advisable to investigate more of those plots and the resulting t-SNE and heatmaps before drawing any conclusions about the final number of “natural” clusters.

Manual merging of up to 20 clusters can be laborious. To simplify this task, one could reduce the strength of over-clustering and allow the metaclustering method to do a part of the merging, which then can be completed manually. Analyzing the delta plot from the right side, we can see how much we can reduce the strength of over-clustering while still obtaining stable clusters. In parallel, one should check the heatmaps to see whether the less stringent stratification is able to capture cell populations of interest.

As an example, we chose to reduce the strength of metaclustering to 12 groups.

## Get cluster ids for each cell nmc2 <- 12 code_clustering2 <- mc[[nmc2]] $ consensusClass cell_clustering2 <- code_clustering2[som $ map $ mapping[, 1 ]]

In the t-SNE plot (see Figure 19), we can see that many small clusters obtained when stratifying data into 20 groups are now merged together, which should simplify the new cluster annotation.

dr $ cell_clustering2 <- factor (cell_clustering2[tsne_inds], levels = 1 : nmc2) ggplot (dr, aes ( x = tSNE1, y = tSNE2, color = cell_clustering2)) + geom_point ( size = 0.8 ) + theme_bw () + scale_color_manual ( values = color_clusters) + guides ( color = guide_legend ( override.aes = list ( size = 4 ), ncol = 2 ))

plot_clustering_heatmap_wrapper ( expr = expr[, lineage_markers_ord], expr01 = expr01[, lineage_markers_ord], cell_clustering = cell_clustering2, color_clusters = color_clusters)

Over-clustering into as few as 12 groups still allows us to identify the same 8 cell populations as when merging 20 clusters (see Figure 20– Figure 23), but it is simpler to define since fewer profiles need to be manually scanned. The expert-based merging is saved in the PBMC8_cluster_merging2.xlsx file on the Robinson Lab server.

cluster_merging2_filename <- "PBMC8_cluster_merging2.xlsx" download.file ( paste0 (url, "/" , cluster_merging2_filename), destfile = cluster_merging2_filename, mode = "wb" ) cluster_merging2 <- read_excel (cluster_merging2_filename) data.frame (cluster_merging2)

## original_cluster new_cluster ## 1 1 B-cells IgM+ ## 2 2 surface- ## 3 3 NK cells ## 4 4 CD8 T-cells ## 5 5 B-cells IgM- ## 6 6 monocytes ## 7 7 CD8 T-cells ## 8 8 CD8 T-cells ## 9 9 monocytes ## 10 10 CD4 T-cells ## 11 11 DC ## 12 12 CD4 T-cells

## Convert to factor with merged clusters in correct order cluster_merging2 $ new_cluster <- factor (cluster_merging2 $ new_cluster, levels = levels_clusters_merged) ## New clustering2m mm <- match (cell_clustering2, cluster_merging2 $ original_cluster) cell_clustering2m <- cluster_merging2 $ new_cluster[mm]

dr $ cell_clustering2m <- cell_clustering2m[tsne_inds] gg_tsne_cl2m <- ggplot (dr, aes ( x = tSNE1, y = tSNE2, color = cell_clustering2m)) + geom_point ( size = 0.8 ) + theme_bw () + scale_color_manual ( values = color_clusters) + guides ( color = guide_legend ( override.aes = list ( size = 4 ))) gg_tsne_cl2m

plot_clustering_heatmap_wrapper ( expr = expr[, lineage_markers_ord], expr01 = expr01[, lineage_markers_ord], cell_clustering = cell_clustering2, color_clusters = color_clusters, cluster_merging = cluster_merging2)

plot_clustering_heatmap_wrapper ( expr = expr[, lineage_markers_ord], expr01 = expr01[, lineage_markers_ord], cell_clustering = cell_clustering2m, color_clusters = color_clusters)

Comparison of automated and manual merging. The manual merging of 20 (or 12) clusters by an expert resulted in identification of 8 cell populations. To highlight the impact of manual merging versus algorithm-defined subpopulations, we compare to the results of an automated cluster merging that is set to stratify the data also into 8 clusters. We extract the result from the ConsensusClusterPlus output. Out of interest, we can see which clusters are split by tabulating the cell labels.

## Get cluster ids for each cell nmc3 <- 8 code_clustering3 <- mc[[nmc3]] $ consensusClass cell_clustering3 <- code_clustering3[som $ map $ mapping[, 1 ]]

# tabular comparison of cell_clustering3 and cell_clustering2m table ( algorithm= cell_clustering3, manual= cell_clustering2m)

## manual ## algorithm B-cells IgM+ B-cells IgM- CD4 T-cells CD8 T-cells DC NK cells ## 1 6651 0 0 0 0 0 ## 2 0 0 0 0 0 0 ## 3 0 0 0 32112 0 24518 ## 4 0 3265 0 0 0 0 ## 5 0 0 0 0 0 0 ## 6 0 0 2603 19038 0 0 ## 7 0 0 60287 0 0 0 ## 8 0 0 0 0 1980 0 ## manual ## algorithm monocytes surface- ## 1 0 0 ## 2 0 3901 ## 3 0 0 ## 4 0 0 ## 5 18436 0 ## 6 0 0 ## 7 0 0 ## 8 0 0

In the t-SNE map (see Figure 24), we can see that part of the new cell populations (cluster 7, 1 and 4, 2, 5 and 8) overlap substantially with populations obtained by the means of manual merging (CD4 T-cells, B-cells, surface-, monocytes and DC). However, cells that belong to clusters 3 and 6 are subdivided in a different manner according to the manual merging. Cluster 3 consists of CD8 T-cells and NK cells, and the latter cannot be identified anymore based on the heatmap corresponding to clustering into 8 groups (see Figure 25).

dr $ cell_clustering3 <- factor (cell_clustering3[tsne_inds], levels = 1 : nmc3) gg_tsne_cl3 <- ggplot (dr, aes ( x = tSNE1, y = tSNE2, color = cell_clustering3)) + geom_point ( size = 0.8 ) + theme_bw () + scale_color_manual ( values = color_clusters) + guides ( color = guide_legend ( override.aes = list ( size = 4 ))) plot_grid (gg_tsne_cl2m, gg_tsne_cl3, ncol = 1, labels = c ( "A" , "B" ))

plot_clustering_heatmap_wrapper ( expr = expr[, lineage_markers_ord], expr01 = expr01[, lineage_markers_ord], cell_clustering = cell_clustering3, color_clusters = color_clusters)

The example above highlights the difference between automatic clustering and manual merging of algorithm-generated clusters in the search for biologically meaningful cell populations. Automated and manual merging may give different weight to marker importance and thus result in different populations being detected. However, in our view, the manual merging done here in a reproducible fashion results in a more biologically meaningful cell stratification.

Differential analysis of cell population abundance compares the proportions of cell types across experimental conditions and aims to highlight populations that are present at different ratios. First, we calculate two tables: one that contains cell counts for each sample and population and one with the corresponding proportions of cell types by sample. The proportions are used only for plotting, since the statistical modeling takes the cell counts by cluster and sample as input.

counts_table <- table (cell_clustering1m, sample_ids) props_table <- t ( t (counts_table) / colSums (counts_table)) * 100

counts <- as.data.frame.matrix (counts_table) props <- as.data.frame.matrix (props_table)

For each sample, we plot its PBMC cell type composition represented with colored bars, where the size of a given stripe reflects the proportion of the corresponding cell type in a given sample (see Figure 26).

ggdf <- melt ( data.frame ( cluster = rownames (props), props), id.vars = "cluster" , value.name = "proportion" , variable.name = "sample_id" ) ggdf $ cluster <- factor (ggdf $ cluster, levels = levels_clusters_merged) ## Add condition info mm <- match (ggdf $ sample_id, md $ sample_id) ggdf $ condition <- factor (md $ condition[mm]) ggplot (ggdf, aes ( x = sample_id, y = proportion, fill = cluster)) + geom_bar ( stat = "identity" ) + facet_wrap ( ∼ condition, scales = "free_x" ) + theme_bw () + theme ( axis.text.x = element_text ( angle = 90 , hjust = 1 )) + scale_fill_manual ( values = color_clusters)

It may be quite hard to see the differences in cluster abundances in the plot above, especially for clusters with very low frequency. And, since boxplots cannot represent multimodal distributions, we show boxplots with jittered points of the sample-level cluster proportions overlaid (see Figure 27). The y-axes are scaled to the range of data plotted for each cluster, to better visualize the differences in frequency of lower abundance clusters. For this experiment, it may be interesting to additionally depict the patient information. We do this by plotting a different point shape for each patient. Indeed, we can see that often the direction of abundance changes between the two conditions are concordant among the patients.

ggdf $ patient_id <- factor (md $ patient_id[mm]) ggplot (ggdf) + geom_boxplot ( aes ( x = condition, y = proportion, color = condition, fill = condition), position = position_dodge (), alpha = 0.5, outlier.color = NA ) + geom_point ( aes ( x = condition, y = proportion, color = condition, shape = patient_id), alpha = 0.8, position = position_jitterdodge ()) + facet_wrap ( ∼ cluster, scales = "free" , nrow = 2 ) + theme_bw () + theme ( axis.text.x = element_blank (), axis.ticks.x = element_blank (), axis.title.x = element_blank (), strip.text = element_text ( size = 6 )) + scale_color_manual ( values = color_conditions) + scale_fill_manual ( values = color_conditions) + scale_shape_manual ( values = c ( 16 , 17 , 8 , 3 , 12 , 0 , 1 , 2 ))

As our goal is to compare proportions, one could take these values, transform them (e.g. logit) and use them as a dependent variable in a linear model. However, this approach does not take into account the uncertainty of proportion estimates, which is higher when ratios are calculated for samples with lower total cell counts. A distribution that naturally accounts for such uncertainty is the binomial distribution (i.e. logistic regression), which takes the cell counts as input (relative to the total for each sample). Nevertheless, as in the genomic data analysis, the pure logistic regression is not able to capture the overdispersion that is present in HDCyto data. A natural extension to model the extra variation is the generalized linear mixed model (GLMM), where the random effect is defined by the sample ID ( Jia et al., 2014; Zhao et al., 2013). Additionally, in our example the patient pairing could be accounted in the model by incorporating a random intercept for each patient. Thus, we present two GLMMs. Both of them comprise a random effect defined by the sample ID to model the overdispersion in proportions. The second model includes a random effect defined by the patient ID to account for the experiment pairing.

In our model, the blocking variable is patient ID i = 1, ..., n, where n = 8. For each patient, there are n _i samples measured, and j = 1, ..., n _i indicates the sample ID. Here, n _i = 2 for all i (one from reference and one from BCR/FcR-XL stimulated).

We assume that for a given cell population, the cell counts Y _ij follow a binomial distribution Y _ij ∼ Bin( m _ij , π _ij ), where m _ij is a total number of cells in a sample corresponding to patient i and condition j. The generalized linear mixed model with observation-level random effects ξ _ij is defined as follows: E ( Y i j | β 0 , β 1 , ξ i j ) = l o g i t − 1 ( β 0 + β 1 x i j + ξ i j ) ,

where ξ i j ∼ N ( 0 , σ ξ 2 ) and x _ij corresponds to the conditionBCRXL column in the design matrix and indicates whether a sample ij belongs to the reference ( x _ij = 0) or treated condition ( x _ij = 1). Since E( Y _ij|β ₀, β ₁, ξ _ij ) = π _ij , the above formula can be written as follows: l o g i t ( π i j ) = β 0 + β 1 x i j + ξ i j .

The generalized linear mixed model that furthermore accounts for the patient pairing incorporates additionally a random intercept for each patient i: E ( Y i j | β 0 , β 1 , γ i , ξ i j ) = l o g i t − 1 ( β 0 + β 1 x i j + γ i + ξ i j ) , where γ i ∼ N ( 0 , σ γ 2 ) .

formula_glmer_binomial1 <- y / total ∼ condition + ( 1 | sample_id) formula_glmer_binomial2 <- y / total ∼ condition + ( 1 | patient_id) + ( 1 | sample_id)

The wrapper function below takes as input a data frame with cell counts (each row is a population, each column is a sample), the metadata table, and the formula, and performs the differential analysis specified with contrast K for each population separately, returning a table with non-adjusted and adjusted p-values.

differential_abundance_wrapper <- function (counts, md, formula, K){ ## Fit the GLMM for each cluster separately ntot <- colSums (counts) fit_binomial <- lapply( 1 : nrow (counts), function (i){ data_tmp <- data.frame( y = as.numeric (counts[i, md $ sample_id]), total = ntot[md $ sample_id], md) fit_tmp <- glmer (formula, weights = total, family = binomial, data = data_tmp) ## Fit contrasts one by one out <- apply(K, 1, function (k){ contr_tmp <- glht (fit_tmp, linfct = matrix (k, 1 )) summ_tmp <- summary (contr_tmp) pval <- summ_tmp $ test $ pvalues return (pval) }) return (out) }) pvals <- do.call (rbind, fit_binomial) colnames (pvals) <- paste0( "pval_" , contrast_names) rownames (pvals) <- rownames (counts) ## Adjust the p-values adjp <- apply (pvals, 2 , p.adjust, method = "BH" ) colnames (adjp) <- paste0 ( "adjp_" , contrast_names) return ( list ( pvals = pvals, adjp = adjp)) }

We fit both of the GLMMs specified above. We can see that accounting for the patient pairing increases the sensitivity to detect differentially abundant cell populations.

da_out1 <- differential_abundance_wrapper (counts, md = md, formula = formula_glmer_binomial1, K = K) apply (da_out1 $ adjp < FDR_cutoff, 2 , table)

## adjp_BCRXLvsRef ## FALSE 5 ## TRUE 3

da_out2 <- differential_abundance_wrapper (counts, md = md, formula = formula_glmer_binomial2, K = K) apply (da_out2 $ adjp < FDR_cutoff, 2 , table)

## adjp_BCRXLvsRef ## FALSE 2 ## TRUE 6

An output table containing the observed cell population proportions in each sample and p-values can be assembled (and optionally written to a file).

da_output2 <- data.frame ( cluster = rownames (props), props, da_out2 $ pvals, da_out2 $ adjp, row.names = NULL ) print ( head (da_output2), digits = 2 )

## cluster BCRXL1 BCRXL2 BCRXL3 BCRXL4 BCRXL5 BCRXL6 BCRXL7 BCRXL8 ## 1 B-cells IgM+ 3.95 1.43 4.1 3.8 3.9 4.0 2.6 2.7 ## 2 B-cells IgM- 1.09 1.01 2.0 1.1 1.5 1.2 2.1 1.7 ## 3 CD4 T-cells 26.81 35.78 32.3 31.8 36.7 44.1 26.8 31.5 ## 4 CD8 T-cells 46.05 40.87 26.5 25.9 28.2 24.6 34.3 39.7 ## 5 DC 0.18 0.83 1.5 1.4 2.1 1.4 1.2 1.2 ## 6 NK cells 15.64 11.69 16.6 18.3 12.6 6.8 25.5 12.9 ## Ref1 Ref2 Ref3 Ref4 Ref5 Ref6 Ref7 Ref8 pval_BCRXLvsRef ## 1 4.82 2.8 8.3 4.7 4.4 5.68 4.34 3.82 3.5e-08 ## 2 1.90 1.3 3.3 1.4 2.5 2.34 2.79 2.19 2.2e-11 ## 3 44.72 49.1 39.7 32.4 38.4 47.33 28.16 36.94 1.9e-03 ## 4 23.66 23.8 15.5 17.6 26.0 25.31 33.49 34.21 1.2e-03 ## 5 0.22 0.9 1.2 1.2 1.6 0.86 0.93 0.89 7.1e-05 ## 6 14.31 9.7 15.1 14.5 10.2 6.67 22.54 10.99 4.5e-13 ## adjp_BCRXLvsRef ## 1 9.2e-08 ## 2 8.8e-11 ## 3 2.5e-03 ## 4 1.9e-03 ## 5 1.4e-04 ## 6 3.6e-12

We use a heatmap to report the differential cell populations (see Figure 28). Proportions are first scaled with the arcsine-square-root transformation (as an alternative to logit that does not return NAs when ratios are equal to zero or one). Then, Z-score normalization is applied to each population to better highlight the relative differences between compared conditions. We created two wrapper functions: normalization_wrapper performs the normalization of submitted expression expr to mean 0 and standard deviation 1, and plot_differential_heatmap_wrapper generates a heatmap of submitted expression expr_norm, where samples are grouped by condition, indicated with a color bar on top of the plot. Additionally, labels of clusters contain the adjusted p-values in parenthesis.

normalization_wrapper <- function (expr, th = 2.5 ){ expr_norm <- apply (expr, 1 , function (x){ sdx <- sd (x, na.rm = TRUE ) if (sdx == 0 ){ x <- (x - mean (x, na.rm = TRUE )) } else { x <- (x - mean (x, na.rm = TRUE )) / sdx } x[x > th] <- th x[x < - th] <- - th return (x) }) expr_norm <- t (expr_norm) }

plot_differential_heatmap_wrapper <- function(expr_norm, sign_adjp, condition, color_conditions, th = 2.5){ ## Order samples by condition oo <- order (condition) condition <- condition[oo] expr_norm <- expr_norm[, oo, drop = FALSE] ## Create the row labels with adj p-values and other objects for pheatmap labels_row <- paste0 ( rownames (expr_norm), " (" , sprintf ( "%.02e" , sign_adjp), ")" ) labels_col <- colnames (expr_norm) annotation_col <- data.frame ( condition = factor (condition)) rownames (annotation_col) <- colnames (expr_norm) annotation_colors <- list ( condition = color_conditions) color <- colorRampPalette ( c ( "#87CEFA" , "#56B4E9" , "#56B4E9" , "#0072B2" , "#000000" , "#D55E00" , "#E69F00" , "#E69F00" , "#FFD700" ))( 100 ) breaks = seq ( from = - th, to = th, length.out = 101 ) legend_breaks = seq ( from = - round (th), to = round (th), by = 1 ) gaps_col <- as.numeric ( table (annotation_col $ condition)) pheatmap (expr_norm, color = color, breaks = breaks, legend_breaks = legend_breaks, cluster_cols = FALSE , cluster_rows = FALSE , labels_col = labels_col, labels_row = labels_row, gaps_col = gaps_col, annotation_col = annotation_col, annotation_colors = annotation_colors, fontsize = 8 ) }

## Apply the arcsine-square-root transformation to the proportions asin_table <- asin ( sqrt (( t ( t (counts_table) / colSums (counts_table))))) asin <- as.data.frame.matrix (asin_table) ## Get significant clusters and sort them by significance sign_clusters <- names ( which ( sort (da_out2 $ adjp[, "adjp_BCRXLvsRef" ]) < FDR_cutoff)) ## Get the adjusted p-values for the significant clusters sign_adjp <- da_out2 $ adjp[sign_clusters , "adjp_BCRXLvsRef" , drop= FALSE ] ## Normalize the transformed proportions to mean = 0 and sd = 1 asin_norm <- normalization_wrapper (asin[sign_clusters, ])

mm <- match( colnames (asin_norm), md $ sample_id) plot_differential_heatmap_wrapper ( expr_norm = asin_norm, sign_adjp = sign_adjp, condition = md $ condition[mm], color_conditions = color_conditions)

For this part of the analysis, we calculate the median expression of the 14 signaling markers in each cell population (merged cluster) and sample. These will be used as the response variable Y _ij in the linear model (LM) or linear mixed model (LMM), for which we assume that the median marker expression follows a Gaussian distribution (on the already arcsinh-transformed scale). The linear model is formulated as follows: Y i j = β 0 + β 1 x i j + ϵ i j ,

where ∈ _ij ∼ N(0, σ ²), and the mixed model includes a random intercept for each patient: Y i j = β 0 + β 1 x i j + γ i + ϵ i j ,

where γ i ∼ N ( 0 , σ γ 2 ) . In the current experiment, we have an intercept (basal level) and a single covariate, x _ij , which is represented as a binary (stimulated/unstimulated) variable. For more complicated designs or batch effects, additional columns of a design matrix can be used.

One drawback of summarizing the protein marker intensity with a median over cells is that all the other characteristics of the distribution, such as bimodality, skewness and variance, are ignored. On the other hand, it results in a simple, easy to interpret approach, which in many cases will be able to detect interesting changes. Another issue that arises from using a summary statistic is the level of uncertainty, which increases as the number of cells used to calculate it decreases. In the statistical modeling, this problem could be partially handled by assigning observation weights (number of cells) to each cluster and sample (parameter weights in the lm and lmer functions). However, since each cluster is tested separately, these weights do not account for the differences in size between clusters.

There might be instances of small cell populations for which no cells are observed in some samples or where the number of cells is very low. For clusters absent from a sample (e.g. due to biological variance or insufficient sampling), NAs are introduced because no median expression can be calculated; in the case of few cells, the median may be quite variable. Thus, we apply a filter to remove samples that have fewer than 5 cells. We also remove cases where marker expression is equal to zero in all the samples, as this leads to an error during model fitting.

## Get median marker expression per sample and cluster expr_median_sample_cluster_tbl <- data.frame (expr[, functional_markers], sample_id = sample_ids, cluster = cell_clustering1m) %>% group_by (sample_id, cluster) %>% summarize_all ( funs (median)) ## Melt expr_median_sample_cluster_melt <- melt (expr_median_sample_cluster_tbl, id.vars = c ( "sample_id" , "cluster" ), value.name = "median_expression" , variable.name = "antigen" ) ## Rearange so the rows represent clusters and markers expr_median_sample_cluster <- dcast (expr_median_sample_cluster_melt, cluster + antigen ∼ sample_id, value.var = "median_expression" ) rownames (expr_median_sample_cluster) <- paste0 (expr_median_sample_cluster $ cluster, "_" , expr_median_sample_cluster $ antigen) ## Eliminate clusters with low frequency clusters_keep <- names ( which (( rowSums (counts < 5 ) == 0 ))) keepLF <- expr_median_sample_cluster $ cluster %in% clusters_keep expr_median_sample_cluster <- expr_median_sample_cluster[keepLF, ] ## Eliminate cases with zero expression in all samples keep0 <- rowSums(expr_median_sample_cluster[, md $sample_id]) > 0 expr_median_sample_cluster <- expr_median_sample_cluster[keep0, ]

It is helpful to plot the median expression of all the markers in each cluster for each sample colored by condition, to get a rough image of how strong the differences might be (see Figure 29). We do this by combining boxplots and jitter.

ggdf <- expr_median_sample_cluster_melt[expr_median_sample_cluster_melt $ cluster %in% clusters_keep, ] ## Add info about samples mm <- match (ggdf $ sample_id, md $ sample_id) ggdf $ condition <- factor (md $ condition[mm]) ggdf $ patient_id <- factor (md $ patient_id[mm]) ggplot (ggdf) + geom_boxplot ( aes ( x = antigen, y = median_expression, color = condition, fill = condition), position = position_dodge (), alpha = 0.5 , outlier.color = NA ) + geom_point ( aes ( x = antigen, y = median_expression, color = condition, shape = patient_id), alpha = 0.8 , position = position_jitterdodge (), size = 0.7 ) + facet_wrap ( ∼ cluster, scales = "free_y" , ncol= 2 ) + theme_bw () + theme ( axis.text.x = element_text ( angle = 90 , vjust = 0.5, hjust = 1 )) + scale_color_manual ( values = color_conditions) + scale_fill_manual ( values = color_conditions) + scale_shape_manual ( values = c ( 16 , 17 , 8 , 3 , 12 , 0 , 1 , 2 )) + guides ( shape = guide_legend ( override.aes = list ( size = 2 )))

We created a wrapper function differential_expression_wrapper that performs the differential analysis of marker expression. The user needs to specify a data frame expr_median with marker expression, where each column corresponds to a sample and each row to a cluster/marker combination. One can choose between fitting a regular linear model model = "lm" or a linear mixed model model = "lmer". The formula parameter must be adjusted adequately to the model choice. The wrapper function returns the non-adjusted and adjusted p-values for each of the specified contrasts K for each cluster/marker combination.

differential_expression_wrapper <- function (expr_median, md, model = "lmer" , formula, K){ ## Fit LMM or LM for each marker separately fit_gaussian <- lapply ( 1 : nrow (expr_median), function (i){ data_tmp <- data.frame ( y = as.numeric (expr_median[i, md $ sample_id]), md) switch (model, lmer = { fit_tmp <- lmer (formula, data = data_tmp) }, lm = { fit_tmp <- lm (formula, data = data_tmp) }) ## Fit contrasts one by one out <- apply (K, 1 , function (k){ contr_tmp <- glht (fit_tmp, linfct = matrix (k, 1 )) summ_tmp <- summary (contr_tmp) pval <- summ_tmp $ test $ pvalues return (pval) }) return (out) }) pvals <- do.call (rbind, fit_gaussian) colnames (pvals) <- paste0 ( "pval_" , contrast_names) rownames (pvals) <- rownames (expr_median) ## Adjust the p-values adjp <- apply (pvals, 2 , p.adjust, method = "BH" ) colnames (adjp) <- paste0 ( "adjp_" , contrast_names) return ( list ( pvals = pvals, adjp = adjp)) }

To present how accounting for the within patient variability with the mixed model increases sensitivity, we also fit a regular linear model. The linear mixed model has a random intercept for each patient.

formula_lm <- y ∼ condition formula_lmer <- y ∼ condition + ( 1 | patient_id)

By accounting for the patient effect, we detect almost twice as many cases of differential signaling compared to the regular linear model.

de_out1 <- differential_expression_wrapper ( expr_median = expr_median_sample_cluster, md = md, model = "lm" , formula = formula_lm, K = K) apply (de_out1 $ adjp < FDR_cutoff, 2 , table)

## adjp_BCRXLvsRef ## FALSE 51 ## TRUE 42

de_out2 <- differential_expression_wrapper ( expr_median = expr_median_sample_cluster, md = md, model = "lmer" , formula = formula_lmer, K = K) apply (de_out2 $ adjp < FDR_cutoff, 2 , table)

## adjp_BCRXLvsRef ## FALSE 23 ## TRUE 70

One can assemble together an output table with the information about median marker expression in each cluster and sample, and the obtained p-values.

de_output2 <- data.frame (expr_median_sample_cluster, de_out2 $ pvals, de_out2 $ adjp, row.names = NULL ) print ( head (de_output2), digits = 2 )

## cluster antigen BCRXL1 BCRXL2 BCRXL3 BCRXL4 BCRXL5 BCRXL6 BCRXL7 ## 1 B-cells IgM+ pNFkB 1.179 0.880 0.808 1.47 1.361 1.725 1.436 ## 2 B-cells IgM+ pp38 0.109 -0.012 0.044 0.24 -0.046 0.083 -0.039 ## 3 B-cells IgM+ pAkt 3.247 2.960 2.951 3.26 2.382 3.184 2.762 ## 4 B-cells IgM+ pStat1 0.343 0.126 0.242 0.33 -0.010 0.616 -0.050 ## 5 B-cells IgM+ pZap70 0.317 0.287 0.351 0.40 0.132 0.604 0.267 ## 6 B-cells IgM+ pStat3 -0.047 -0.059 0.451 0.35 -0.058 -0.026 0.534 ## BCRXL8 Ref1 Ref2 Ref3 Ref4 Ref5 Ref6 Ref7 Ref8 ## 1 1.5747 1.9639 1.869 1.7726 2.1833 1.861 1.953 1.915 1.979 ## 2 -0.0055 0.8891 1.113 0.8534 0.6424 0.126 0.210 0.128 0.126 ## 3 3.1439 2.3195 2.310 2.2688 3.0858 1.729 2.024 2.145 2.603 ## 4 0.3795 -0.0058 0.064 0.0079 0.5151 -0.047 0.030 -0.034 0.191 ## 5 0.3202 -0.0198 -0.033 -0.0336 -0.0056 -0.061 -0.060 -0.032 -0.017 ## 6 0.3092 -0.0479 -0.082 0.2652 0.1567 -0.060 -0.066 0.275 0.381 ## pval_BCRXLvsRef adjp_BCRXLvsRef ## 1 6.1e-11 2.7e-10 ## 2 7.5e-04 1.6e-03 ## 3 2.6e-11 1.3e-10 ## 4 6.2e-02 7.5e-02 ## 5 1.6e-14 1.0e-13 ## 6 5.6e-02 7.1e-02

To report the significant results, we use a heatmap (see Figure 30). Instead of plotting the absolute expression, we display the normalized expression, which better highlights the direction of marker changes. Additionally, we order the cluster-marker instances by their significance and group them by cell type (cluster).

## Keep the significant markers, sort them by significance and group by cluster sign_clusters_markers <- names ( which (de_out2 $ adjp[, "adjp_BCRXLvsRef" ] < FDR_cutoff)) oo <- order (expr_median_sample_cluster[sign_clusters_markers, "cluster" ], de_out2 $ adjp[sign_clusters_markers, "adjp_BCRXLvsRef" ]) sign_clusters_markers <- sign_clusters_markers[oo]

## Get the significant adjusted p-values sign_adjp <- de_out2 $ adjp[sign_clusters_markers , "adjp_BCRXLvsRef" ]

## Normalize expression to mean = 0 and sd = 1 expr_s <- expr_median_sample_cluster[sign_clusters_markers,md $ sample_id] expr_median_sample_cluster_norm <- normalization_wrapper (expr_s)

mm <- match ( colnames (expr_median_sample_cluster_norm), md $ sample_id) plot_differential_heatmap_wrapper ( expr_norm = expr_median_sample_cluster_norm, sign_adjp = sign_adjp, condition = md $ condition[mm], color_conditions = color_conditions)

The analysis of overall expression is analogous to the previous section, except that median marker expression is aggregated from all the cells in a given sample, Figure 31.

ggdf <- melt ( data.frame (expr_median_sample[functional_markers, ], antigen = functional_markers), id.vars = "antigen" , value.name = "median_expression" , variable.name = "sample_id" ) ## Add condition info mm <- match (ggdf $ sample_id, md $ sample_id) ggdf $ condition <- factor (md $ condition[mm]) ggdf $ patient_id <- factor (md $ patient_id[mm]) ggplot (ggdf) + geom_boxplot ( aes ( x = condition, y = median_expression, color = condition, fill = condition), position = position_dodge (), alpha = 0.5 , outlier.color = NA ) + geom_point ( aes ( x = condition, y = median_expression, color = condition, shape = patient_id), alpha = 0.8 , position = position_jitterdodge ()) + facet_wrap ( ∼ antigen, scales = "free" , nrow = 5 ) + theme_bw () + theme ( axis.text.x = element_blank (), axis.ticks.x = element_blank ()) + scale_color_manual ( values = color_conditions) + scale_fill_manual ( values = color_conditions) + scale_shape_manual ( values = c ( 16 , 17 , 8 , 3 , 12 , 0 , 1 , 2 ))

Similar to the analysis above, we identify more markers being differentially expressed with the LMM, which accounts for the within patient variability.

## Fit a linear model de_out3 <- differential_expression_wrapper ( expr_median = expr_median_sample[functional_markers, ], md = md, model = "lm" , formula = formula_lm, K = K) apply (de_out3 $ adjp < FDR_cutoff, 2 , table)

## adjp_BCRXLvsRef ## FALSE 9 ## TRUE 5

## Fit a linear mixed model with patient ID as a random effect de_out4 <- differential_expression_wrapper ( expr_median = expr_median_sample[functional_markers, ], md = md, model = "lmer" , formula = formula_lmer, K = K) apply (de_out4 $ adjp < FDR_cutoff, 2 , table)

## adjp_BCRXLvsRef ## FALSE 3 ## TRUE 11

As before, we create an output table with the median marker expression calculated in each sample and the p-values, and we plot a heatmap with the significant markers sorted by their statistical significance ( Figure 32).

de_output4 <- data.frame ( antigen = functional_markers, expr_median_sample[functional_markers, ], de_out4 $ pvals, de_out4 $ adjp) print ( head (de_output4), digits= 2 )

## antigen BCRXL1 BCRXL2 BCRXL3 BCRXL4 BCRXL5 BCRXL6 BCRXL7 BCRXL8 ## pNFkB pNFkB 1.070 0.520 1.144 1.7397 1.143 1.6951 1.195 1.245 ## pp38 pp38 -0.052 -0.057 -0.024 -0.0034 -0.061 -0.0006 -0.064 -0.053 ## pStat5 pStat5 -0.043 -0.058 -0.041 -0.0103 -0.075 -0.0404 -0.067 -0.060 ## pAkt pAkt 3.053 2.727 2.876 3.2424 1.958 2.6068 2.075 2.416 ## pStat1 pStat1 1.356 1.096 1.504 1.6960 0.535 1.9823 0.526 0.566 ## pSHP2 pSHP2 -0.053 -0.057 -0.048 -0.0352 -0.073 -0.0435 -0.068 -0.065 ## Ref1 Ref2 Ref3 Ref4 Ref5 Ref6 Ref7 Ref8 ## pNFkB 2.392 2.469 2.670 2.940 1.979 2.025 1.980 1.985 ## pp38 1.280 1.474 1.467 0.939 0.091 0.218 0.062 0.122 ## pStat5 -0.049 -0.049 -0.040 0.100 -0.065 -0.063 -0.059 -0.052 ## pAkt 2.416 2.122 2.300 3.273 1.443 1.573 1.680 2.114 ## pStat1 0.889 0.930 1.474 2.056 0.493 1.097 0.416 0.549 ## pSHP2 -0.071 -0.068 -0.059 -0.038 -0.080 -0.069 -0.073 -0.066 ## pval_BCRXLvsRef adjp_BCRXLvsRef ## pNFkB 2.3e-09 1.1e-08 ## pp38 1.0e-03 1.8e-03 ## pStat5 3.0e-01 3.0e-01 ## pAkt 3.0e-06 8.3e-06 ## pStat1 1.9e-01 2.1e-01 ## pSHP2 5.4e-04 1.1e-03

## Keep the significant markers and sort them by significance sign_markers <- names ( which ( sort (de_out4 $ adjp[, "adjp_BCRXLvsRef" ]) < FDR_cutoff)) ## Get the adjusted p-values sign_adjp <- de_out4 $ adjp[sign_markers , "adjp_BCRXLvsRef" ] ## Normalize expression to mean = 0 and sd = 1 expr_median_sample_norm <- normalization_wrapper (expr_median_sample[sign_markers, ]) mm <- match ( colnames (expr_median_sample_norm), md $ sample_id) plot_differential_heatmap_wrapper ( expr_norm = expr_median_sample_norm, sign_adjp = sign_adjp, condition = md $ condition[mm], color_conditions = color_conditions)