In farrellja/URD: URD

\fontsize{8}{18}

library("knitr")
opts_chunk$set(tidy.opts=list(width.cutoff=80),tidy=TRUE,dev="png",dpi=150)

library(URD)
library(rgl)
library(gridExtra) # For arranging plots
library(RColorBrewer) # For color palettes

# Set up knitr to capture rgl output
rgl::setupKnitr()

Load previous saved object

object <- readRDS("obj/object_6_tree.rds")

Color Palettes

Define some color palettes to use for figures. Most of these are gentle modifications of RColorBrewer palettes.

# Colors to use for stage
stage.colors <- c("#CCCCCC", RColorBrewer::brewer.pal(9, "Set1")[9], RColorBrewer::brewer.pal(12, "Paired")[c(9,10,7,8,5,6,3,4,1,2)])

# Preference colors for the preference plots
pref.colors <- c("#CECECE", "#CBDAC2", RColorBrewer::brewer.pal(9, "YlGnBu")[3:9])

# Red-orange color scheme for gene expression
fire.with.grey <- c("#CECECE", "#DDC998", RColorBrewer::brewer.pal(9, "YlOrRd")[3:9])

# Grey-blue-green color scheme for module and gene expression
pond.with.grey <- c("#CECECE", "#CBDAC2", RColorBrewer::brewer.pal(9, "YlGnBu")[3:9])

branch.colors <- c("#CECECE", "#E6298B")

plotDim

The plotDim command allows plotting cells according to various dimensionality reductions that have been performed (such as tSNE, PCA, or the diffusion map.)

tSNE - Stage

plotDim(object, "stage.nice", reduction.use = "tSNE", discrete.colors=stage.colors, plot.title = "Stage")

tSNE - Clustering

plotDim(object, "ZF6S-Cluster", reduction.use = "tSNE", plot.title="6-somite clusters")

tSNE - Gene Expression

Individual markers can be plotted with the default blue-to-red color scheme, or a custom palette.

plotDim(object, "PRDX5", reduction.use = "tSNE", plot.title="PRDX5 (An endoderm marker)")
plotDim(object, "TAGLN2", reduction.use = "tSNE", plot.title="TAGLN2 (Pharyngeal Endoderm and Epidermal Marker", colors = fire.with.grey)

Two markers can be simultaneously plotted using a red-green color scheme.

plotDimDual(object, label.red="PRDX5", label.green="TAGLN2", reduction.use = "tSNE", plot.title="PRDX5 (Red) vs. TAGLN2 (Green)", legend.offset.x = 5)

Diffusion Map

Pairs of components from the diffusion map (or PCA, if reduction.use="pca") can also be plotted with any metadata or gene expression.

plotDim(object, reduction.use = "dm", dim.x=3, dim.y=4, label = "n.Genes", plot.title = "Number of Genes")

Additionally, pairs of labels can be plotted in a red-green color scheme.

plotDimDual(object, reduction.use="dm", dim.x=3, dim.y=4, label.red = "GSC", label.green="ICN")

Furthermore, any plotDim command can be run against many pairs of components using plotDimArray -- all parameters are passed to plotDim, and the x and y dimensions are taken in pairs from the dims.to.plot vector.

plotDimArray(object, reduction.use="dm", dims.to.plot = 1:8, label="NOTO", outer.title="NOTO expression", plot.title="")

Tree Dendrogram

URD's recovered tree dendrogram can also be decorated with any gene expression or metadata.

plotTree(object, "NOTO", title = "NOTO expression")

plotTree(object, "n.Trans", title="")

plotTree(object, "Louvain-15")

Force-directed Layout

Similarly, the force-directed layout can be decorated with any clustering, gene expression, or metadata.

plotTreeForce(object, "stage.nice", title="STAGE", title.line=1, discrete.colors = stage.colors, alpha=0.4)

plotTreeForce(object, "WNT8A", title="WNT8A expression", title.line=1)

plotTreeForce(object, "NOTO", title="NOTO expression", title.line=1, colors = fire.with.grey)

object <- groupFromCells(object, group.id="lineage_Tailbud", cells=cellsAlongLineage(object, "Tailbud", remove.root=F))
plotTreeForce(object, "lineage_Tailbud", title="Tailbud", title.line=1, discrete.colors=branch.colors, alpha=0.4)

Gene Expression

We also include several functions for exploring gene expression within lineages or populations. Here we use a violin plot and dot plot to illustrate markers of the clusters used as hindbrain tips.

plotViolin(object, labels.plot = c("EGR2B", "FGF8A", "MAFBA"), clustering="ZF6S-Cluster", clusters=c("Hindbrain R3", "Hindbrain R4", "Hindbrain R5+6"))

plotDot(object, genes=c("EGR2B", "IRX3A", "HOXA2B", "HOXB2A", "MAFBA", "FGF8A", "HOXA3A", "HOXB3A"), clustering = "ZF6S-Cluster", clusters.use = c("Hindbrain R3", "Hindbrain R4", "Hindbrain R5+6"))

Markers of a specific lineage

In the manuscript, we illustrated the expression of the determined markers for a single cascade on the force-directed layout (Figure S5 shows markers of the prechordal plate). We took all genes that were part of the prechordal plate gene cascade and applied an additional layer of specificity -- we required that they were ~4 times better than a random precision-recall classifier when compared globally between the prechordal plate cascade and the rest of the embryo. We then determined whether they had already been annotated as having expression in the prechordal plate (or one of its synonyms) in ZFIN (the Zebrafish Information Network). Finally, we plotted the expression of each gene, colored according to whether it had been previously annotated or not.

# Load gene cascade
pcp.cascade <- readRDS("cascades/impulse/casc_Prechordal Plate.rds")
pcp.markers <- rownames(pcp.cascade$scaled.expression)

# Determine which genes are also global markers
pcp.axial.cells <- cellsInCluster(object, "segment", c("29","79"))
pcp.markers.global <- markersAUCPR(object, cells.1 = pcp.axial.cells, genes.use=pcp.markers)
marker.thresh <- aucprThreshold(cells.1=pcp.axial.cells, cells.2=setdiff(unlist(object@tree$cells.in.segment), pcp.axial.cells), factor=2.5, max.auc = Inf)
pcp.de.markers <- pcp.markers.global[pcp.markers.global$AUCPR >= marker.thresh,]

# Who is annotated already in ZFIN?
# Search terms: anterior axial hypoblast, prechordal plate, polster, hatching gland
zfin.pcp.markers <- read.csv(file="data/ZFIN_annotated_Markers.csv", header=F)
new.markers <- setdiff(rownames(pcp.de.markers), toupper(zfin.pcp.markers$V1))

# Still had to hand verify these in ZFIN, since sometimes genes have been renamed or have a weird but equivalent anatomy term (like "dorsal marginal blastomeres")
renamed.in.zfin <- c("HE1A", "HE1B", "LGALS3L", "SHISA2", "CHD", "OTX1A", "OTX1B", "ATPIF1B", "FBXO2")
known.marker <- c("RIPPLY1", "NDR1", "LHX1A", "MIXL1", "ISM1", "FSCN1A", "CTH1", "DHRS3B", "SND1") # Genes with semi-random, but equivalent anatomy terms

# Final new/old markers list.
new.markers <- setdiff(new.markers, c(renamed.in.zfin, known.marker)) # 49
old.markers <- setdiff(rownames(pcp.de.markers), new.markers) # 71
pcp.markers <- c(new.markers, old.markers) # 120

# Order markers according to gene cascade
timing <- pcp.cascade$timing[pcp.markers,]
timing[intersect(which(is.na(timing$time.on)), which(is.infinite(timing$time.off))), "time.on"] <- Inf
ordered.markers <- pcp.markers[order(timing$time.on, timing$time.off, na.last=F)]
ordered.markers.new <- ordered.markers %in% new.markers

# Save the list of markers for later.
write(ordered.markers, "cascades/pcp_markers_FigS5.txt")

We plot all of the markers to a folder for building the supplemental figure...

# Plot each marker, colored based on whether it was annotated in ZFIN previously.
for (i in 1:length(ordered.markers)) {
  m <- ordered.markers[i]
  if (ordered.markers.new[i]) colors.use <- pond.with.grey else colors.use <- fire.with.grey
  plotTreeForce(object, m, alpha=0.7, alpha.fade=0.08, size=10, density.alpha=T, label.tips=F, colors=colors.use, view = "figure1")
  text3d(x=-8, y=4.2, z=100, m, cex=4)
  Sys.sleep(0.2)
  rgl.snapshot(file=paste0("cascades/pcp_markers/", sprintf("%03d", i), "-", m, ".png"))
  rgl.close()
}

... but include here a couple of examples of a new marker (pnocb) and a very classic marker of the prechordal plate (ctslb/hgg1):

plotTreeForce(object, ordered.markers[71], title=ordered.markers[71], alpha=0.7, alpha.fade=0.08, size=10, density.alpha=T, label.tips=F, colors=pond.with.grey, view = "figure1")

plotTreeForce(object, label=ordered.markers[72], title=ordered.markers[72], alpha=0.7, alpha.fade=0.08, size=10, density.alpha=T, label.tips=F, colors=fire.with.grey, view = "figure1")

Preference plot at a branchpoint

In the manuscript, we described that the axial mesoderm branchpoint (between the notochord and prechordal plate) has cells that expressed genes characteristic of both downstream populations (Figure 6). We illustrated this using preference plots, colored by gene expression.

Define the preference layout for the branchpoint

First, the preference plot layout is defined. Cells that were visited by random walks started in the two tips in question are included in the layout. Cells are ordered along the y-axis according to pseudotime. Cells are placed along the x-axis based on their preference, which is based on their ratio of visits by the random walks from each tip.

# Define layout for plots
np.layout <- branchpointPreferenceLayout(object, pseudotime = "pseudotime", lineages.1 = "29", lineages.2 = "32", parent.of.lineages = "79", opposite.parent = c("72","78"), min.visit = 1)

Plot stage on the preference plot

We found that the intermediate cells were prevalent at mid-gastrulation, from 60%-90% epiboly.

plotBranchpoint(object, np.layout, label="stage.nice", point.alpha=0.5, populations = c("P", "N"), pt.lim=c(0.7,0.1), xlab="", ylab="", legend=T, axis.lines = F, fade.low=0, discrete.colors = stage.colors[c(1,3:12)], title="Stage")

Plot gene expression on the branchpoint.

We found that the intermediate cells no longer expressed progenitor markers (nanog, mex3b), expressed early markers of both cell types (ta, noto, gsc, and frzb), and expressed late markers of the notochord (ntd5 and shha), but did not express late markers of the prechordal plate (prdm1a, icn). By the time that differentiation genes (col8a1a and he1a) were expressed, there were no longer intermediate cells between the two trajectories.

# Define genes to plot
axial.genes.plot <- c("NANOG", "TA", "NOTO", "NTD5", "SHHA", "COL8A1A", "MEX3B", "GSC", "FRZB", "PRDM1A", "ICN", "HE1A")

# Plot gene expression on the branchpoint preference plot
axial.branchpoint.plots <- lapply(axial.genes.plot, function(gene) plotBranchpoint(object, np.layout, label=gene, point.alpha=1, populations = c("P", "N"), pt.lim=c(0.7,0.11), color.scale = pref.colors, xlab="", ylab="", title=gene, legend=F, axis.lines = F, fade.low=0.66))

grid.arrange(grobs=axial.branchpoint.plots, ncol=6)

farrellja/URD documentation built on June 17, 2020, 4:48 a.m.