ClusterMap - Compare multiple Single cell RNA-seq profiling

First, use a list of marker genes of each sub-group in each sample to map sub-groups across samples. The marker gene file can be the direct output of "FindAllMarkers" of Seruat package, or can be a table with columns of "cluster" and "gene".
The "edge_cutoff" can be adjust to decide if a singleton should be merged with a paired groups or not. See "2.2 More options"

library(ClusterMap)

marker_file_list <- c(p7d = 'pre_analysis/p7d.markers.csv', p7e = 'pre_analysis/p7e.markers.csv')

res <- cluster_map(marker_file_list, edge_cutoff = 0.1, output = 'p7')
res

By calling "cluster_map", the mapped sub-groups information can be generate as the following figures and table. They are auto-saved.

Cluster match

Hierachical clustering:

The sub-groups were matched by hierarchical clustering. The distance between the sub-groups across the samples were determined by the existence of the marker genes in each group. The union of all marker genes were considered together. The sub-groups of each sample were clustered by the averaged binary distance of the existence of the marker genes.
The similarity of the matched groups is defined as 1 minus the height of the merging node of the matched groups in the dendragram.

Tree cut:

We search through all the nodes in the dendragram from the bottom. If all the offspring nodes of a node are from the same sample, the node is considered as a pure node. For a given node N with two direct sub-nodes n1 and n2, N will be kept only under three conditions:

• n1 and n2 are both singleton or pure.
  
• n1 is singleton or pure, AND n2 is not removed, AND the edge between N and n2 is smaller than the cutoff.
  
• none of n1 or n2 is cut, AND none of the edge between N and n1 or n2 is bigger than the cutoff, AND S1 and S2 are not included by each other. S1 is the unique sample list that all the offspring nodes of n1 comes from, similar for S2.

Cluster match results were saved as file p7.cluster.map.csv . The contents are

Further more, besides the marker gene list, if the Seurat object of sinlge dataset and combined dataset are provided, the full analysis can be done.

marker_file_list <- c(p7d = 'pre_analysis/p7d.markers.csv', p7e = 'pre_analysis/p7e.markers.csv')

fList <- c(p7d = 'pre_analysis/p7d.RDS', p7e = 'pre_analysis/p7e.RDS', p7 = 'pre_analysis/p7.RDS')
objList <- lapply(fList, readRDS)

single_obj_list <- c(p7d = objList$p7d, p7e = objList$p7e)

res <- cluster_map(marker_file_list, edge_cutoff = 0.1, output = 'p7', single_obj_list = single_obj_list, comb_obj = objList$p7)
res

By calling "cluster_map" this way, Circos plot, re-colored t-SNE plot will be generated as following. As well as the cell percentage change and sample separability of each new group.

Circos plot and population size

To visualized the results, we use Circos plot to summarize.

• The cords of the Circos plot indicate the linkage between the sub-groups.
  
• The percentage of cell numbers in each sample is represented by the width of the black sectors. The sub-population size change is intuitive by comparing the sector size of matched (cord linked) groups.
  
• The color of the cords indicates different new groups, which matches the recolored t-SNE plot.
  
• The transparency of the cord color indicates the similarity of the matched groups.

Recolor

For easy comparison, each sample and combined sample are recolored by their new groups after matching under t-SNE coordinates in the pre-analysis. Cells belonging to the same new group in different samples are recolored in the same color. Thus, some samples may miss one or more colors, which means no such group in the sample. New groups will be saved as a list into RDS file, which can be used later.

Separability

After matching and regrouping, we estimate the condition effect on each paired groups besides the percentage change. To quantify the difference between paired groups, separability was calculated for each new group. Separability will be generated pairwise if more than two samples.

For each new group, the separability is defined as the median difference of inna- and inter- sample distance of each cell in the combined t-SNE coordinates. The inna-sample distance is defined as the median distance of K-nearest cells of the same sample. Defualt K=5. The inter-sample distance is the same but cells of different sample.

Using median instead of mean will reduce the variation due to the outliners. The higher the more separable.

Map results

From the following results, we see that new group 2 and 4 are one to multiple match. There is no matched groups for the new group 8, 10, 11, 12 at current cutoff. Both circos plot and recolored t-SNE plot showed the cell percentage change of the group 2, 4 and 9. The separability of each new group further emerged the most affected sub-types, such as group 5, 7 and 4, besides the no matched groups.

The full results were saved as in file p7.results.csv

Core basal genes: Krt14, Acta2, Myl9, Sparc, Mylk, Cxcl14 – group 3, 4, 7, 8
Luminal genes: Elf5, Prlr, Areg, Ly6d, Stc2, Krt19 – group 1, 11, 12, 9
ML cell-specific genes: Prlr, Cited1, Esrrb, and Cxcl15 – group 2, 5, 6, 10

library(Seurat)

FeaturePlot(object = objList$p7d, features.plot = c("Krt14", "Acta2", "Elf5", "Prlr", "Krt19","Pgr"), cols.use = c("grey", "blue"), reduction.use = "tsne")
savePlot('feature.gene.p7d.png')

FeaturePlot(object = objList$p7e, features.plot = c("Krt14", "Acta2", "Elf5", "Prlr", "Krt19","Pgr"), cols.use = c("grey", "blue"), reduction.use = "tsne")
savePlot('feature.gene.p7e.png')

P7D

P7E

The cutoff is the edge length cutoff for deciding if a singleton should be merged or not. We suggest the cutoff is better to be less than 0.2.

marker_file_list <- c(p7d = 'pre_analysis/p7d.markers.csv', p7e = 'pre_analysis/p7e.markers.csv')
## cutoff 0.1
mapRes <- cluster_map_by_marker(marker_file_list, cutoff = 0.1, output = 'p7')
## cutoff 0.15
mapRes <- cluster_map_by_marker(marker_file_list, cutoff = 0.15, output = 'p7_cutoff0.15')
mapRes ## results didn't change.
## cutoff 0.2
mapRes <- cluster_map_by_marker(marker_file_list, cutoff = 0.2, output = 'p7_cutoff0.2')
mapRes ## p7d_9 was grouped together with p7d_10 and p7e_8.

cutoff = 0.2

The cell numbers for each sub-group can be pulled out from Seruat object or provided by the user. cell_num_list is a list of vectors for each sample. Each vector contains the cell numbers for each sub-group of the sample.

fList <- c(p7d = 'pre_analysis/p7d.RDS', p7e = 'pre_analysis/p7e.RDS')
objList <- lapply(fList, readRDS)
single_obj_list <- c(p7d = objList$p7d, p7e = objList$p7e)

cell_num_list <- lapply(single_obj_list, function(obj) summary(obj@ident))

mapRes <- read.csv('p7.results.csv') ## This file was auto-saved when cluster_map_by_marker function was run.

circos_map(mapRes, cell_num_list, output = 'p7')
res <- add_perc(mapRes, cell_num_list)
res

Number of samples in circos plot

If there are more than 3 samples, the circos plot will become complicated. We can choose a subset of samples to plot. To achieve this, just remove extra samples from the mapRes and cell_num_list variables. Then use circos_map to plot a subset.

p7 <- readRDS('pre_analysis/p7.RDS')
tsne_coord <- as.data.frame(p7@dr$tsne@cell.embeddings)

new_group_list <- readRDS('p7.new.group.list.RDS')
group <- new_group_list$comb
head(group)

sample_label <- as.factor(sub('-.*', '', names(group)))
## k=5
sepa_k5 <- separability_pairwise(tsne_coord, group, sample_label, k = 5)
sepa_k5
## k=10
sepa_k10 <- separability_pairwise(tsne_coord, group, sample_label, k = 10)
sepa_k10

k = 5

k = 10