\fontsize{8}{18}
library("knitr") opts_chunk$set(tidy.opts=list(width.cutoff=80),tidy=TRUE,dev="png",dpi=150)
library(URD) library(rgl) # Set up knitr to capture rgl output rgl::setupKnitr()
setwd("~/Dropbox/Jeff-Yiqun/URD-walkthrough/")
object <- readRDS("obj/object_5_withWalks.rds")
A few clusters that were run as separate tips are totally intermixed in the diffusion map. In these cases, it's often best to combine the two tips before starting to build the tree. (This averages their visitation frequency, according to the number of cells in each tip, and avoids having to re-run the random walks.)
# Load tip cells object <- loadTipCells(object, tips="ZF6S-Cluster-Num") # Combine a few sets of tips where you walked from two groups of cells # that probably should be considered one, based on the fact that they # are intermixed in the diffusion map. # Diencephalon object <- combineTipVisitation(object, "14", "27", "14") object <- combineTipVisitation(object, "14", "19", "19") # Optic Cup object <- combineTipVisitation(object, "2", "9", "2") # Combine epidermis and integument object <- combineTipVisitation(object, "1", "36", "1") # Tailbud object <- combineTipVisitation(object, "3", "13", "3")
Additionally, only "tips" that are actually terminal populations should be used in construction of the tree if possible. We ran random walks from all 6-somite clusters, and so we exclude several of them here, based on prior knowledge or based on their position in the diffusion map.
# Decide on the tips to use in the tree construction tips.to.exclude <- c("4", "6", "7", "9", "11", "13", "14", "15", "20", "23", "24", "25", "27", "28", "30", "35", "36", "37", "39", "41", "42", "44") tips.to.use <- setdiff(as.character(1:53), tips.to.exclude)
# Build the tree object.built <- buildTree(object = object, pseudotime="pseudotime", divergence.method = "ks", tips.use=tips.to.use, weighted.fusion = T, use.only.original.tips = T, cells.per.pseudotime.bin=80, bins.per.pseudotime.window = 5, minimum.visits = 1, visit.threshold = 0.7, p.thresh = 0.025, save.breakpoint.plots = NULL, dendro.node.size = 100, min.cells.per.segment = 10, min.pseudotime.per.segment = .01, verbose = F)
First, just use the names of the clusters to inspect the tree structure.
tip.names <- unique(object@group.ids[,c("ZF6S-Cluster", "ZF6S-Cluster-Num")]) tip.names <- tip.names[complete.cases(tip.names),] object.built <- nameSegments(object.built, segments=tip.names$`ZF6S-Cluster-Num`, segment.names=tip.names$`ZF6S-Cluster`)
plotTree(object.built)
And then after inspecting the tree, choose a set of refined names to use going forward, including short names that will look better in the force-directed layout.
# Descriptive names that will be used on dendrogram new.seg.names <- c("Spinal Cord", "Diencephalon", "Optic Cup", "Midbrain+Neural Crest", "Hindbrain R3", "Hindbrain R4+5+6", "Telencephalon", "Epidermis", "Neural Plate Border", "Placode Adeno.+Lens+Trigeminal", "Placode Epibranchial+Otic", "Placode Olfactory", "Tailbud", "Adaxial Cells", "Somites", "Hematopoietic (ICM)", "Hematopoietic (RBI)+Pronephros", "Endoderm Pharyngeal", "Endoderm Pancreatic+Intestinal", "Heart Primordium", "Cephalic Mesoderm", "Prechordal Plate", "Notochord", "Primordial Germ Cells", "EVL/Periderm") # Short names / Abbreviations for use on force-directed layout new.short.names <- c("SC", "Di", "Optic", "MB+NC", "HB R3", "HB R4-6", "Tel", "Epi", "NPB", "P(A+L+T)", "P(E+Ot)", "P(Olf)", "TB", "Adax", "Som", "Hem(ICM)", "Hem(RBI)+Pro", "Endo Phar", "Endo Pan+Int", "Heart", "CM", "PCP", "Noto", "PGC", "EVL") # Segment numbers segs.to.name <- c("8", "19", "2", "59", "50", "56", "16", "1", "10", "55", "57", "53", "3", "34", "12", "52", "58", "17", "26", "5", "18", "29", "32", "40", "38") # Run the naming object.built <- nameSegments(object.built, segments = segs.to.name, segment.names=new.seg.names, short.names=new.short.names)
plotTree(object.built, label.segments=T)
genes.plot <- c("SOX17", "NOTO", "TA", "GSC", "MEOX1", "GATA2A", "NANOS3", "KRT4", "FSTA", "WNT8A", "CRABP2A", "EGR2B", "ENG2B", "TAL1", "MAFBA", "EMX3") for (gene in genes.plot) { plot(plotTree(object.built, gene)) }
We generate a force-directed layout as a nice visualization of the data. It is constructed by generating a weighted k-nearest neighbor network, based on euclidean distance in visitation space (using the frequency of visitation of each cell by the walks from each tip). The nearest neighbor network is optionally refined based on cells' distance in the dendrogram (in terms of which segments are connected). Cells then push and pull against their neighbors, as the amount of freedom they have to move is slowly decreased until cells are locked into place. This produces a two-dimensional layout, and we then add pseudotime as a third dimension.
We find that the force directed layout works best if the most poorly visited (and thus likely poorly connected) cells are excluded.
# Data frame to measure cell visitation visitation <- data.frame( cell=rownames(object.built@diff.data), seg=object.built@diff.data$segment, stringsAsFactors=F, row.names=rownames(object.built@diff.data) ) visitation$visit <- log10(apply(visitation, 1, function(cr) object.built@diff.data[as.character(cr['cell']), paste0("visitfreq.raw.", as.character(cr['seg']))])+1) # Choose those cells that were well visited robustly.visited.cells <- visitation[visitation$visit >= 0.5, "cell"] # Since some tips of the tree were combined in their entirety, get the terminal segments to use as the tips of the force-directed layout. final.tips <- segTerminal(object.built)
It can be important to try several sets of parameters here. Varying the number of nearest neighbors (num.nn) and the amount of refinement based on the dendrogram (cut.unconnected.segments) affects the layout significantly.
# Generate the force-directed layout object.built <- treeForceDirectedLayout(object.built, num.nn=120, pseudotime="pseudotime", method = "fr", dim = 2, cells.to.do = robustly.visited.cells, tips=final.tips, cut.unconnected.segments = 2, min.final.neighbors=4, verbose=T)
Once calculated and plotted, the plot can be rotated to a desired view, and you can save the view using the function plotTreeForceStore3DView. Thus, many plots with a comparable orientation can then be produced.
Since the force-directed layout is not deterministic, we instead load a saved layout and orientation to finish the tutorial.
# Load view object.built@tree$force.view.list <- readRDS("fdls/force.view.list.rds") object.built@tree$force.view.default <- "figure1" # Load layout precalc.fdl <- readRDS("fdls/layout.51.cells05.nn120.mn4.cut2.rds") object.built@tree$walks.force.layout <- precalc.fdl$walks.force.layout
plotTreeForce(object.built, "segment", alpha=0.2, view="figure1")
For optimal 2D presentation in the paper, we further tuned the layout by hand to reduce overlaps and ensure as much of the tree is visible simultaneously as possible. This was done by increasing the angular distance at the first branchpoint, and moving the two completely disconnected populations (EVL & PGCs).
## ECTODERM # Rotate the ectoderm branch around the z-axis to optimize the # orientation of the neural and non-neural ectoderm later. object.built <- treeForceRotateCoords(object.built, seg="72", angle = -3.3, axis="z", around.cell = 10, throw.out.cells=1000, pseudotime="pseudotime") # Curve the ectoderm outwards and forwards, to prevent it # overlapping the mesoderm, so that the branching structure # within the two domains is easily viewed. for (throw.out in c(0,100,250,500,1000)) { object.built <- treeForceRotateCoords(object.built, seg="72", angle = pi/20, axis="y", around.cell = 10, throw.out.cells=throw.out, pseudotime="pseudotime") } for (throw.out in c(0,100,250,500,1000)) { object.built <- treeForceRotateCoords(object.built, seg="72", angle = pi/24, axis="z", around.cell = 10, throw.out.cells=throw.out, pseudotime="pseudotime") } ## AXIAL MESODERM # Rotate the axial mesoderm a bit to the right to make # some extra space for the remainder of the mesendoderm. object.built <- treeForceRotateCoords(object.built, seg="79", angle = -pi/10, axis="y", around.cell = 10, throw.out.cells=0, pseudotime="pseudotime") object.built <- treeForceRotateCoords(object.built, seg="79", angle = -pi/8, axis="x", around.cell = 10, throw.out.cells=0, pseudotime="pseudotime") ## REMAINDER OF THE MESENDODERM # Rotate the rest of the mesoderm a bit to the right into the # empty space between the ectoderm and axial mesoderm. object.built <- treeForceRotateCoords(object.built, seg="78", angle = pi/4, axis="z", around.cell = 10, throw.out.cells=0, pseudotime="pseudotime") object.built <- treeForceRotateCoords(object.built, seg="78", angle = -pi/5, axis="y", around.cell = 10, throw.out.cells=0, pseudotime="pseudotime") ## PGCs # The PGCs are disconnected from the tree totally, # so just rotate and move them into place. object.built <- treeForceRotateCoords(object.built, seg="40", angle = -pi/2, axis="y", around.cell = 1, throw.out.cells=0, pseudotime="pseudotime") object.built <- treeForceTranslateCoords(object.built, seg="40", x=0, y=0, z=-3) ## EVL / Periderm # The EVL/Periderm is also nearly completely disconnected. # Rotate these cells, and also close the enormous gap in them # somewhat so that they fit neatly next to the ectoderm. # Determine the EVL cells evl.cells <- intersect(cellsInCluster(object.built, "segment", "38"), rownames(object.built@tree$walks.force.layout)) evl.cells.move.1 <- evl.cells[which(object.built@tree$walks.force.layout[evl.cells, "telescope.pt"] > 15)] evl.cells.move.2 <- evl.cells[which(object.built@tree$walks.force.layout[evl.cells, "telescope.pt"] > 30)] # Close the gaps a bit in this lineage to shorten it so that it doesn't overlap with the ectoderm object.built <- treeForceTranslateCoords(object.built, cells=evl.cells.move.1, x=0, y=0, z=-10) object.built <- treeForceTranslateCoords(object.built, cells=evl.cells.move.2, x=0, y=0, z=-15) # Rotate the EVL cells in next to the ectoderm object.built <- treeForceRotateCoords(object.built, seg="38", angle = 1.35, axis="y", around.cell = 1, throw.out.cells=0, pseudotime="pseudotime")
plotTreeForce(object.built, "segment", alpha=0.2, view="figure1")
plotTreeForce(object.built, "SOX17", alpha=0.2, view="figure1")
plotTreeForce(object.built, "GSC", alpha=0.2, view="figure1")
plotTreeForce(object.built, "NOTO", alpha=0.2, view="figure1")
plotTreeForce(object.built, "TAL1", alpha=0.2, view="figure1")
plotTreeForce(object.built, "GATA2A", alpha=0.2, view="figure1")
plotTreeForce(object.built, "SOX19A", alpha=0.2, view="figure1")
saveRDS(object.built, file="obj/object_6_tree.rds")
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.