Trajectory analysis of flow cytometry data

Explanation

The first concept to define are partitions. Monocle3 partitions act as supergroups of cells, or metaclusters if you will. They are meant to identify larger groups of cells that are still somewhat connected despite being organised in multiple clusters, but to separate from other groups of cells with a larger distance gap. A good example would be, in biological terms, CD4 and CD8 T cells. CD4 T cells would be part of one partition, and CD8 T cells would be part of another. They both contain naive, effector and memory subsets, but differ enough to be separated.

By default, the trajectory analysis will try to find trajectories across all cells, but very often that doesn’t make much biological sense. There is no reason for CD4 T cells to differentiate into CD8 T cells in the periphery. Partitioning helps you set up the analysis to generate parallel trajectories within each partition.

However, partitioning is not a perfect capture all of each group of cells. Sometimes, you are looking at related entities (say Naive CD4 T cells) and activated effector cells, but partitioning split them into separate partitions because the gap between them is too large: this is something you see when looking at a very potent stimulating agent like a TLR agonist, or a MHC tetramer.

Partitioning is meant to guide you in how to analyze your data, but ultimately, it’s your knowledge of the biology and the hypotheses that you’re testing that will be the final judge.

Code

By default, all cells are set to partition "1" automatically. By running SetPartitions(), you can either let monocle3 detect partitions automatically (automatic = T), or supply it with a data.frame that contains 1. row.names as cellID and 2. partition labels in a column called partition. The former can be more practical if you want to have your partitions named instead of generic numbers.

th_trajectory <- SetPartitions(th_trajectory, automatic = T)


plot_umap(th_trajectory, group.by = "Condition") + plot_umap(th_trajectory, group.by = "partition")

Here, you see that cells towards the top (Th2; Partition 1) are split into a separate partition from the bottom (Partition 1) despite having a small group of cells connecting them to Partition 2. We know that biologically , those Th0 and Th1 cells are related, as the former differentiates into the latter. The partitioning here suggests that either 1. there are not enough events in that “bridge”, or that Th2 cells are too different from the rest phenotypically. Indeed, you can see they separate in the UMAP plot quite drastically.

Cluster subsampling

It is not always advisable to generate a trajectory analysis on the intact dataset, especially if you subsampled more heavily as mentioned in the Before you start section. While this approach makes it easier to find transition states, it also makes it computationally heavy.

The solution for this, is to perform cluster-based subsampling. This subsampling method randomly selects n events per cluster to equalize the data and make it more computationally manageable.

flow_obj_th <- ClusterSubsampling(th_trajectory, n = 2000)

Export flow_object to cds

Monocle3 works on cell_data_set objects, not flow_object. Unsupflowhelper provides a wrapper function to generate this object. However, you can provide either the intact dataset (if it’s computationally feasible) or the subsampled one.

Furthermore, you can use the SubsetFlowObject function to select the partition(s) you wish to generate a trajectory analysis for.

cds <- flow_to_cds(SubsetFlowObject(flow_obj_th, subset = partition %in% c("1", "2")))

### monocle3 has some plotting functions of it's own. 
plot_cells(cds, color_cells_by = "cluster", show_trajectory_graph = F)

Learn graph

This is where monocle3 will attempt to find a trajectory across the cells. By default, monocle3 will use embedded partitions to find trajectories within each. If you wish to perform a trajectory analysis within partitions, set use_partition = FALSE.

Also, different parameters will influence if the trajectory gets closed in a loop or open-ended, or the minimal branch length of secondary branches. More details in monocle3’s help documentation.

cds <- learn_graph(cds, close_loop = F, use_partition = F, learn_graph_control = list(minimal_branch_len = 15))


plot_cells(cds, color_cells_by = "cluster", show_trajectory_graph = T)

cds umap plot

Pick graph segments

Having a graph is a great first start, but now we have to give it some directionality and meaning.

First step would be to infer where the tree starts. If your dataset has multiple timepoints, you could compare them using plot_density_umap() to figure out if one tree end matches with a region that is enriched in the earliest timepoint. Another approach is to use your knowledge of the biological system: in this case, we know that the bottom left clusters, i.e. nodes 1,8 and 9, are located in the region where Th0 samples are located. Given that uncommitted Th0 cells are known to differentiate into Th1 or Th2 cells, this seems like a good starting point.

Second, we observe that the graph splits off in 2 trajectories at the node 7. Thus we are looking at 2 possible trajectories. Each of those trajectories will have a parallel analysis.

cds_th1 <- choose_graph_segments(cds, clear_cds = F)
cds_th1 <- order_cells(cds_th1)

cds_th2 <- choose_graph_segments(cds, clear_cds = F)
cds_th2 <- order_cells(cds_th2)

Pseudotime imputation

If you did use Cluster subsampling, you have lineages and markers, but it will be difficult to compare enrichment of cells witin or across trajectories in different metadata groups, as you potentially subsampled more heavily one group than another.

This issue is then solved by imputing pseudotime values from nearest neighbors within a trajectory.

th_trajectory <- impute_pseudotime(th_trajectory, trajectory = "Th1 trajectory")

Imputation

Pseudotime UMAP

Now that trajectories and pseudotime have been added to your flow object, you have access to more visuals.

You can visualize pseudotime values for a given trajectory on the UMAP plot with the plot_umap_pseudotime function.

plot_umap_pseudotime(th_trajectory, trajectory = "Th1 trajectory", dot_size = 1)

Moreover, you can plot the trajectories on other UMAP plots, such as plot_umap, plot_density_umap or plot_umap_heat with the parameter show_trajectories = T.

plot_umap(th_trajectory, show_trajectories = T, group.by = "louvain")

plot_umap_heat(th_trajectory, markers = c("Tbet","GATA3"), show_trajectories = T )

Pseudotime density

You may want to compare sample condition distributions

plot_grouped_pseudotime(th_trajectory, group.by = "Condition", trajectory = "Th1 trajectory", adjust = 2)

Trajectory markers

You will likely wish to see the expression of markers in relation to a given trajectory.

By default, the plot_pseudotime_markers function will produce a heatmap for all included markers.

plot_pseudotime_markers(th_trajectory, trajectory = "Th1 trajectory")

Marker heatmap

However, as you can see, some markers (i.e. GATA3) do not vary in relation to pseudotime. You can visually appreciate that, perhaps select some biologically relevant markers, or get markers from FindTrajectoryMarkers (above).

plot_pseudotime_markers(th_trajectory, trajectory = "Th1 trajectory", markers = c("Tbet", "IFNg", "CD45RO", "CD27", "CD45RA"))

Selected markers heatmap