The extreme diversity of retinal amacrine cells has deep evolutionary roots.

Code ↔ Paper

30 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 30 matches

[1] § RESULTS › Evolutionarily grounded annotation of oACs ↔ src/5_tf_analysis.Rmd, lines 1831–1919 · score 0.96 · oAC39, oAC14, oAC7, oAC37, oAC13, oAC19
[2] § RESULTS › A conserved regulatory logic underlies AC heterogeneity ↔ src/5_tf_analysis.Rmd, lines 1831–1919 · score 0.96 · ZNF804A, oAC4, oAC3, oAC7, oAC6, oAC2
[3] § RESULTS › Evolutionarily grounded annotation of oACs ↔ src/utils/objects.R, lines 3–12 · score 0.96 · oAC39, oAC14, oAC7, oAC37, oAC13, oAC19
[4] § RESULTS › Conservation of AC diversity in basal vertebrates ↔ src/utils/objects.R, lines 3–12 · score 0.91 · oAC37, oAC8, oAC6, oAC2, oAC36, oAC1
[5] § RESULTS › Evolutionarily grounded annotation of oACs ↔ src/8_figures.Rmd, lines 3708–3741 · score 0.89 · oAC20, oAC40, oAC14, oAC18, oAC38, oAC13
[6] § RESULTS › A conserved regulatory logic underlies AC heterogeneity ↔ src/2_orthotype_analysis.Rmd, lines 4058–4131 · score 0.82 · ZNF804A, TFAP2B, EBF1, NEUROD2, BNC1, BNC2
[7] § RESULTS › Neurotransmitter identity of oACs ↔ src/1_species_clustering.Rmd, lines 715–734 · score 0.80 · Slc6a11, Slc6a13, Slc6a5, Slc6a9, GAD2, GABAergic
[8] § RESULTS › Coordinated evolution of AC and RGC diversity ↔ src/8_figures.Rmd, lines 2662–2742 · score 0.76 · nearest neighbor, SAMap alignment, TFAP2B, lamprey clusters, NN, RBPMS
[9] § RESULTS › Integration across amniotes identifies 42 oACs ↔ src/8_figures.Rmd, lines 1008–1076 · score 0.70 · confidence interval, ferret SACs, SLC5A7, precision, recall, discrete
[10] § MATERIALS AND METHODS › Regulatory tree construction by MP ↔ src/utils/wrappers.R, lines 3032–3148 · score 0.68 · equally parsimonious trees, consensus tree, pratchet, binary, MP, seq
[11] § MATERIALS AND METHODS › Integration of nonmammalian atlases ↔ src/8_figures.Rmd, lines 2662–2742 · score 0.67 · nearest neighbor classification, lamprey cluster, unmapped, alignment, MG, predict
[12] § MATERIALS AND METHODS › Statistical analysis › Enrichment of protein-protein interactions ↔ src/7_sac_analysis.Rmd, lines 724–770 · score 0.66 · protein interactions, interaction score, STRINGdb, v12, downloaded, genes
[13] § RESULTS › Cross-species histological validation of select AC orthotypes ↔ src/8_figures.Rmd, lines 1604–1649 · score 0.65 · VGluT3, SLC17A8, nNOS, GAD, NOS1, NXPH1
[14] § RESULTS › Shared evolutionary origin of GABAergic ACs and RGCs ↔ src/8_figures.Rmd, lines 2848–2880 · score 0.63 · TFAP2B, ambiguous cluster, NEFL, NEFM, OPN4, RBPMS2
[15] § RESULTS › Integration across amniotes identifies 42 oACs ↔ src/2_orthotype_analysis.Rmd, lines 839–865 · score 0.59 · integrated LISI, iLISI, cLISI, VG3, Seurat, SAC
[16] § MATERIALS AND METHODS › Regulatory tree construction by MP ↔ src/utils/wrappers.R, lines 3032–3148 · score 0.58 · equally parsimonious, Branch lengths, ACCTRAN, MP, tree
[17] § MATERIALS AND METHODS › Statistical analysis › Enrichment of protein-protein interactions ↔ src/5_tf_analysis.Rmd, lines 2007–2062 · score 0.57 · STRINGdb, protein interactions, shuffled, Enrichment, tree, score
[18] § MATERIALS AND METHODS › Assessing robustness of amniote orthotypes ↔ src/2_orthotype_analysis.Rmd, lines 839–865 · score 0.57 · iLISI, cLISI, metrics, precision, recall, VG3
[19] § MATERIALS AND METHODS › Dimensionality reduction and clustering within each species ↔ src/utils/documented-functions.R/clustering.R, lines 76–193 · score 0.57 · principal components, PCs, dimensionality, Seurat, variable, log
[20] § MATERIALS AND METHODS › Assessing robustness of amniote orthotypes ↔ src/8_figures.Rmd, lines 799–843 · score 0.56 · cLISI, iLISI, mixing, amniote, VG3, SAC
[21] § RESULTS › Shared evolutionary origin of GABAergic ACs and RGCs ↔ src/utils/objects.R, lines 839–905 · score 0.54 · TFAP2B, NEFL, NEFM, OPN4, RBPMS2, NOS1
[22] § RESULTS › AC atlases across 24 vertebrates ↔ src/utils/objects.R, lines 1040–1085 · score 0.54 · snRNA, mouse lemur, sc, newt, axolotl, killifish
[23] § RESULTS › Integration across amniotes identifies 42 oACs ↔ src/8_figures.Rmd, lines 799–843 · score 0.53 · cLISI, iLISI, mixed, shuffled, amniotes, VG3
[24] § MATERIALS AND METHODS › Integration of nonmammalian atlases ↔ src/8_figures.Rmd, lines 3018–3119 · score 0.52 · BC integration, Leiden clustering, cell class, nonmammals, Jaccard, heatmaps
[25] § MATERIALS AND METHODS › Assessing robustness of amniote orthotypes ↔ src/utils/documented-functions.R/clustering.R, lines 1–73 · score 0.52 · nearest neighbors, clustering resolution, identity
[26] § RESULTS › Conservation of AC diversity across amniotes ↔ src/8_figures.Rmd, lines 2055–2143 · score 0.51 · tree shrew, species clusters, curve, abundance, quantify, opossum
[27] § RESULTS › Conservation of AC diversity in basal vertebrates ↔ src/utils/wrappers.R, lines 3448–3561 · score 0.51 · genome duplication, cell classes, retinal, PR, goldfish, orthology
[28] § MATERIALS AND METHODS › Regulatory tree construction by MP ↔ src/5_tf_analysis.Rmd, lines 1243–1277 · score 0.51 · anc_pars, ACCTRAN, reconstructing, TF, tree
[29] § MATERIALS AND METHODS › Statistical analysis › Divergence analysis ↔ src/utils/wrappers.R, lines 547–608 · score 0.51 · pseudo bulked, PC, space, distances, BC, RGC
[30] § MATERIALS AND METHODS › Integration of nonmammalian atlases ↔ src/utils/wrappers.R, lines 3448–3561 · score 0.50 · clustering resolutions, cell class, downstream, atlases, Seurat, goldfish

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

R Markdown · 4,412 lines · 195 KB · MIT · 10 matches

---
title: "Amacrine cell manuscript figures"
output:
BiocStyle::html_document:
toc: true
---
Load necessary libraries for analysis
```{r setup}
source("../../utils/dario_functions.R")
LoadLibraries(load.lisi = FALSE)
SourceFiles()
JS.THRESHOLD = 0.10
```
Add new orthotype labels
```{r}
ac.ortho = LoadACOrtho()
# order = levels(ac.ortho$type)
meta = Metadata2(ac.ortho, c('type', 'type_no', 'lit_type', 'displaced'))
old.values = levels(ac.ortho$type)
# annotation = ExtractString(old.values, before = '_')
# new.levels = gsub(' \\[\\]', '', paste0('oAC', 1:42, ' [', annotation, ']'))
number = paste0('oAC', 1:42, ifelse(meta$displaced, '*', ''))
new.levels = gsub(' \\[NA\\]', '', paste0(number, ' [', meta$lit_type, ']'))
ac.ortho$orthotype = factor(convert_values(ac.ortho$type, key_table = data.frame(old.names = old.values,
new.names = new.levels)),
levels = new.levels)
table(ac.ortho$orthotype)
saveRDS(ac.ortho, '../../Ortho_Objects/vertebrateAC_BIGSEURAT.rds')
```
# Species figures
## Raster solution
```{r}
library(magick)
dir.create('../../Evolution/images/PhyloPic/colored2/')
colors = species_palette3
input_dir <- "../../Evolution/images/PhyloPic"
files <- list.files(input_dir, pattern = "\\.svg$", full.names = TRUE)
for (f in files) {
nm <- tools::file_path_sans_ext(basename(f))
if (!nm %in% names(colors)) next
img <- image_read_svg(f)
img <- image_colorize(img, opacity = 100, color = colors[nm])
image_write(img, paste0('../../Evolution/images/PhyloPic/colored2/', gsub('.svg', '.png', basename(f))), format = "png")
image_write(img, paste0('../../Evolution/images/PhyloPic/colored2/', basename(f)), format = "svg")
image_write(img, paste0('../../Evolution/images/PhyloPic/colored2/', gsub('.svg', '.pdf', basename(f))), format = "pdf")
}
```
## Raster solution (gray)
```{r}
library(magick)
dir.create('../../Evolution/images/PhyloPic/gray/')
colors = species_palette3
input_dir <- "../../Evolution/images/PhyloPic"
files <- list.files(input_dir, pattern = "\\.svg$", full.names = TRUE)
for (f in files) {
nm <- tools::file_path_sans_ext(basename(f))
if (!nm %in% names(colors)) next
img <- image_read_svg(f)
img <- image_colorize(img, opacity = 100, color = 'grey15')
image_write(img, paste0('../../Evolution/images/PhyloPic/gray/', basename(f)), format = "svg")
image_write(img, paste0('../../Evolution/images/PhyloPic/gray/', gsub('.svg', '.pdf', basename(f))), format = "pdf")
}
```
## SVG solution (fails for half)
```{r}
library(rsvg)
outdir <- "../../Evolution/images/PhyloPic/colored3/"
dir.create(outdir, showWarnings = FALSE)
colors <- species_palette3
input_dir <- "../../Evolution/images/PhyloPic"
files <- list.files(input_dir, pattern = "\\.svg$", full.names = TRUE)
for (f in files) {
nm <- tools::file_path_sans_ext(basename(f))
if (!nm %in% names(colors)) next
color <- colors[nm]
# Read SVG as text
svg_text <- readLines(f, warn = FALSE)
svg_text <- paste(svg_text, collapse = "\n")
# Replace all black fills (both style and direct attribute)
new_svg_text <- svg_text
new_svg_text <- gsub("fill:#000000", paste0("fill:", color), new_svg_text, ignore.case = TRUE)
new_svg_text <- gsub('fill="#000000"', paste0('fill="', color, '"'), new_svg_text, ignore.case = TRUE)
if(identical(svg_text, new_svg_text)){
new_svg_text <- gsub(
'<path([^>]*?)(/?>)',
paste0('<path\\1 fill="', color, '"\\2'),
svg_text
)
}
if (identical(svg_text, new_svg_text)) {
message("No black fill found in ", nm, "; skipping PDF.")
next
}
# Save to temporary SVG
tmp_svg <- tempfile(fileext = ".svg")
writeLines(new_svg_text, tmp_svg)
# Output PDF
out_pdf <- file.path(outdir, paste0(nm, ".pdf"))
if (file.exists(out_pdf)) file.remove(out_pdf)
rsvg::rsvg_pdf(tmp_svg, out_pdf)
message("Saved PDF for ", nm, " with fill color ", color)
}
```
# Legends
```{r}
ggarrange(get_legend(JSHeatmap2(ac.umap$orthotypes, ac.umap$leiden_clusters, col.high = 'chartreuse4') + theme(legend.position = "bottom")),
get_legend(JSHeatmap2(ac.umap$orthotypes, ac.umap$leiden_clusters) + theme(legend.position = "bottom")),
get_legend(JSHeatmap2(ac.umap$orthotypes, ac.umap$leiden_clusters, col.high = 'cyan4') + theme(legend.position = "bottom")))
ggsave('../../figures/my_figs/jaccard-legends.pdf')
```
```{r}
plot_grid(get_legend(PrettyUmap2(ac.ortho.no.tetrapods,
group.by = "orthotype.no",
color.by = "species",
cols = species_palette2,
show.legend = TRUE,
geom.label = geom_label,
color = 'black', fill = 'white',
label.size = 0, alpha = 0.5,
label.padding = unit(0, "lines"),
label.r = unit(0.2, "lines"),
label = TRUE)))
ggsave('figures/mammals-legend.png')
```
# Figure 1: intro
## Phylogeny
```{r, fig.height=5, fig.width=10}
key = data.table::fread("../../Evolution/phylogeny/species_common_latin_24.txt")
write.table(key$LatinName, "../../Evolution/phylogeny/species_latin_24.txt", quote = FALSE, row.names = FALSE, col.names = FALSE)
# Make newick tree using timetree (https://timetree.org)
tree = ape::read.tree("../../Evolution/phylogeny/species_latin_24.nwk")
tree$tip.label = key$CommonName[match(gsub("_", " ", tree$tip.label), key$LatinName)]
tree = ape::rotate(tree, node = c(25))
plot(tree)
saveRDS(tree, '../../Ortho_Objects/phylo.hs.pm.rds')
# Pruned
dendro = phylogram::as.dendrogram.phylo(tree)
sorted.dendro = reorder(dendro, 1:length(tree$tip.label))
dendro.hs.liz = dendextend::prune(sorted.dendro, c('Lamprey', 'Cat shark', 'Killifish', "Zebrafish", 'Goldfish', 'Newt', 'Axolotl'))
dendro.hs.liz
labels(dendro.hs.liz)
plot(as.phylo(dendro.hs.liz))
saveRDS(as.phylo(dendro.hs.liz), '../../Ortho_Objects/phylo.hs.liz.rds')
```
```{r, fig.height=2.5, fig.width=6}
tree = readRDS('../../Ortho_Objects/phylo.hs.pm.rds')
tree_fixed = tree
tree_fixed$edge.length <- tree$edge.length / 2
dendro = as.dendrogram(tree_fixed)
labels(dendro) = species_names_pretty((labels(dendro)))
# Sort by species order
dendro = dendextend::rotate(dendro, rev(species.names.pretty))
plot(dendro)
# saveRDS(dendro, '../../Ortho_Objects/phylo.hs.pm.sorted.dendro.rds')
library(ggdendro)
ggdendrogram(dendro) +
# coord_flip() +
# scale_y_reverse(expand = c(0.2, 0))+
# theme_classic()+
ylab("Evolutionary distance (MYA)")+
# theme(axis.text.x = element_text(size = 13, hjust = 1, vjust = 1, angle = 45), axis.ticks.length.y = unit(.25, "cm"))+
scale_y_continuous(breaks = seq(0, 500, len = 6), expand = expansion(mult = c(0.01,0.01)))+ #expand = expansion(add = c(0,0.1)))+
scale_x_discrete(expand = expansion(add = 0.6))+
theme_cowplot()+
theme(axis.line.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.x=element_blank(),
axis.title.y = element_text(angle = -90, hjust = 0.5),
axis.text.y = element_text(angle = -90, hjust = 0.5),
panel.grid.minor.x=element_blank(),
panel.grid.major.x=element_blank())
ggsave('../../figures/my_figs/fig-phylo-24.pdf', height=2.5, width=6)
```
```{r, fig.height=2.5, fig.width=6}
dendro = readRDS('../../Ortho_Objects/phylo.hs.pm.sorted.dendro.rds')
dd <- dendro_data(dendro)
ggplot() +
geom_segment(
data = dd$segments,
aes(x = x, y = y, xend = xend, yend = yend),
linewidth = .75
) +
ylab("Evolutionary distance (MYA)") +
scale_y_continuous(
breaks = seq(0, 500, length.out = 6),
expand = expansion(mult = c(0.01, 0.01))
) +
scale_x_discrete(
expand = expansion(add = 0.6)
) +
theme_cowplot() +
theme(
axis.line.y = element_line(linewidth = 0.75),
axis.line.x = element_blank(),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.y = element_text(angle = -90, hjust = 0.5),
axis.text.y = element_text(angle = -90, hjust = 0.5),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank()
)
ggsave('../../figures/my_figs/fig-phylo-24.pdf', height=2.5, width=6)
```
## Quantification
```{r}
# Save Catshark
temp = readRDS('../../Full_Objects/Shark_full_v5.rds')
temp$type = temp$annotated
temp$animal = temp$orig.ident
saveRDS(subset(temp, cell_class == 'AC'), '../../Species_Objects/SharkAC_v6.rds')
# Save Axolotl
temp = readRDS('../../Full_Objects/Axolotl_full_v3.rds')
temp$type = temp$annotated
temp$animal = temp$animal_group
temp$species = 'Axolotl'
saveRDS(subset(temp, cell_class == 'AC'), '../../Species_Objects/AxolotlAC_v6.rds')
# Save Newt
temp = readRDS('../../Full_Objects/Newt_full_v3.rds')
temp$type = temp$annotated
temp$animal = 'N1'
temp$species = 'Newt'
saveRDS(subset(temp, cell_class == 'AC'), '../../Species_Objects/NewtAC_v6.rds')
# Save Killifish
temp = readRDS('../../Species_Objects/Killifish_ncbiAC_v6.rds')
temp$animal = temp$orig.file
temp$classification = ifelse(temp$cell_class2 == 'gabaAC', 'GABA', 'Gly')
saveRDS(temp, '../../Species_Objects/Killifish_ncbiAC_v6.rds')
```
```{r}
objectList = ReadAmacrineData2()
objectList$Zebrafish$animal = objectList$Zebrafish$sample
objectList$Mouse$animal = objectList$Mouse$batch
# If there is no enrichment info, set enrichment to none
objectList = lapply(objectList, function(object) {
if(!'enrichment' %in% colnames([email hidden])){
object$enrichment = "NONE"
}
object
})
# Count total cells
sum(sapply(objectList, function(object) ncol(object)))
max(sapply(objectList, function(object) ncol(object)))
```
```{r, fig.height=6, fig.width=7}
# nTypes.df = readRDS('metadata/nTypes.df.v2.rds')
# PrettyBarplot(nTypes.df, x = 'species', y = 'nCells') + coord_flip()
line.width = 0.5
bar.width = 0.85
nTypes.df = data.frame(species = factor(names(objectList), levels = names(objectList)),
nClusters = sapply(objectList, function(object) length(unique(object$type))),
nCells = sapply(objectList, function(object) ncol(object)),
nReplicates = sapply(objectList, function(object) length(unique(object$animal))),
nUMI = sapply(objectList, function(object) mean(object$nFeature_RNA)))
# Use Yan et al. number of clusters for mouse
nTypes.df$nClusters[nTypes.df$species == 'Mouse'] = 63
nTypes.df$nClusters[nTypes.df$species == 'Macaque'] = 44
nTypes.df$species2 = sapply(seq_along(nTypes.df$species), function(x) {
if(x %% 2 == 0) paste0('', nTypes.df$species[[x]]) else paste0('', nTypes.df$species[[x]])
}) %>% species_labels
nTypes.df$species2 = factor(nTypes.df$species2, levels = unique(nTypes.df$species2))
# coord_flip_discrete(ggbarplot(nTypes.df, x = 'species', y = 'nCells'), 'species')
PrettyBarplot(flip_var(nTypes.df, 'species'), x = 'species', y = 'nCells', fill = 'species', color = 'black', width = bar.width) +
# scale_y_continuous(name = '# of cells', labels = label_number(scale = 1e-3, suffix = "k"), expand = expansion(mult = c(0,0)))+
# scale_y_continuous(name = '# of cells', transform = 'log10', expand = expansion(mult = c(0,0)))+
# scale_y_log10(name = 'Number of\ncells', breaks = 10^(0:4),
# # labels = comma,
# labels = scales::trans_format("log10", scales::math_format(10^.x)),
# minor_breaks = rep(1:9, times = 4) * 10^(rep(0:3, each = 9)),
# expand = expansion(mult = c(0,0))) +
scale_y_log10(
name = 'Number of\ncells',
# limits = c(10^(1), 10^(5)),
# limits = c(10, NA),
breaks = 10^(1:4),
labels = scales::trans_format("log10", scales::math_format(10^.x)),
minor_breaks = rep(1:9, times = 3) * 10^(rep(1:3, each = 9)),
expand = expansion(mult = c(0, 0))
)+
# coord_cartesian(ylim = c(10, Inf))+
annotation_logticks(sides = "b", outside = TRUE,
short = unit(0.6,"mm"), mid = unit(0.5,"mm"),long = unit(1,"mm")) +
# coord_cartesian(clip = 'off')+
scale_fill_manual(values = species_palette3)+
theme(axis.ticks.y = element_blank(),
axis.text.y = element_blank(),
axis.line = element_line(linewidth = line.width),
# axis.text.x = element_text(angle = 90),
# axis.title.x = element_text(angle = 180),
axis.title.y = element_blank()
)+
coord_flip(clip = 'off')+
NoLegend() |
PrettyBarplot(flip_var(nTypes.df, 'species'), x = 'species', y = 'nUMI', fill = 'species', color = 'black', width = bar.width) +
# scale_y_continuous(name = 'Number of cells', labels = label_number(scale = 1e-3, suffix = "k"), expand = expansion(mult = c(0,0)))+
scale_fill_manual(values = species_palette3)+
# ylab('Number of\nbiological\nreplicates')+
ylab('Number of\ngenes detected')+
theme(axis.ticks.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_line(linewidth = line.width),
axis.line = element_line(linewidth = line.width),
axis.title.y = element_blank(),
# axis.text.x = element_text(angle = 90),
# axis.title.x = element_text(angle = 180)
)+
coord_flip()+
NoLegend() |
PrettyBarplot(flip_var(nTypes.df, 'species2'), x = 'species2', y = 'nClusters', fill = 'species', color = 'black', width = bar.width) +
# scale_y_continuous(name = 'Number of cells', labels = label_number(scale = 1e-3, suffix = "k"), expand = expansion(mult = c(0,0)))+
scale_fill_manual(values = species_palette3)+
scale_x_discrete(position = "top")+ # labels = sapply(1:24, function(x) ))+
# scale_y_
ylab('Number of\nnominal\nclusters')+
geom_hline(yintercept = 0, linewidth = line.width)+
theme(axis.title.y = element_blank(),
axis.line = element_line(linewidth = line.width),
# axis.text.x = element_text(angle = 90),
axis.ticks = element_line(linewidth = line.width),
axis.text.y.right = element_text(size = 12, margin = margin(l = 3), color = darken(rev(species_palette3), 0.1)),
# axis.title.x = element_text(angle = 180)
)+
# axis.text.y.right = element_text(angle = 135, hjust = 1, vjust = 1))+# = element_text(angle = 135, hjust = 0, vjust = 0.5))+
coord_flip(clip = 'off')+
NoLegend()
ggsave('../../figures/my_figs/species-cell-cluster-quant-tight.pdf', height = 6, width = 5)
```
```{r, fig.height=6, fig.width=7}
# objectList = mclapply(speciesList, function(species) readRDS(paste0("../../Species_Objects/", species, "AC_v6.rds")), mc.cores = 1)
# names(objectList) = speciesList
merged.metadata = Reduce(rbind, lapply(objectList, function(object) [email hidden][,c("classification", "enrichment", "seurat_clusters", "type", 'animal')])) #"lit_type",
merged.metadata$species = factor(rep(speciesList, sapply(objectList, ncol)), levels = speciesList)
merged.metadata[,"enrichment"] = toupper(merged.metadata[,"enrichment"])
merged.metadata$`Enrichment group` = gsub('\n', '', merged.metadata$enrichment)
merged.metadata$`Enrichment group` = gsub('NONE', 'No enrichment', merged.metadata$`Enrichment group`)
merged.metadata$`Enrichment group` = factor(merged.metadata$`Enrichment group`,
levels = c('NEUN+', 'NEUN+CHX10-',
'CHX10+',
'NEUN-', 'NEUN-CHX10-', 'NEUN-CHX10+',
'CD90+',
'CD73-',
'No enrichment'))
coord_flip_discrete(stackedBarGraph2(merged.metadata, x = 'species', y = 'Enrichment group', normalize = FALSE, border.color = 'black'), discrete_axis = 'Var1') +
scale_fill_manual(values = ClusterPalette(c(3,3,1,4,4,4,5,2,6), vary.by = 0.2,
colors = c('deeppink', 'orange', 'chartreuse1', 'cyan2', 'violet', 'antiquewhite2'), dark.first = FALSE)) +
theme(axis.text.y = element_blank(), axis.title.y = element_blank(), axis.title.x = element_text()) +
ArialFont()+
scale_y_continuous(name = 'Number of cells', labels = label_number(scale = 1e-3, suffix = "k"), expand = expansion(mult = c(0,0)))
ggsave('../../figures/my_figs/species-enrich-quant.pdf', height = 6, width = 7)
```
Number of cells and number of nominal clusters
```{r, fig.height=6, fig.width=7}
nTypes.df = readRDS('metadata/nTypes.df.v2.rds')
nReplicates = apply(table(merged.metadata$animal, merged.metadata$species), 2, function(col) length(which(col > 0)))
nTypes.df$nReplicates = nReplicates[match(nTypes.df$species, names(nReplicates))]
saveRDS(nTypes.df, 'metadata/nTypes.df.v2.rds')
```
```{r, fig.height=5, fig.width=7}
nTypes.df = readRDS('metadata/nTypes.df.v2.rds')
# PrettyBarplot(nTypes.df, x = 'species', y = 'nCells') + coord_flip()
library(ggplot2)
library(scales)
library(patchwork)
log10_reverse <- trans_new(
name = "log10-reverse",
transform = function(x) -log10(x),
inverse = function(x) 10^(-x),
breaks = log_breaks(base = 10),
format = math_format(10^.x)
)
species_order <- rev(unique(nTypes.df$species)) # reverse order for downward bars
p1 <- PrettyBarplot(flip_var(nTypes.df, 'species'),
x = 'species', y = 'nCells', fill = 'species',
color = NA, width = 0.9) +
# scale_y_log10(name = 'Number of\ncells', breaks = 10^(0:4),
# labels = trans_format("log10", math_format(10^.x)),
# minor_breaks = rep(1:9, times = 4) * 10^(rep(0:3, each = 9)),
# expand = expansion(mult = c(0,0))) +
scale_y_continuous(
name = 'Number of\ncells',
trans = log10_reverse,
breaks = 10^(0:4),
minor_breaks = rep(1:9, times = 4) * 10^(rep(0:3, each = 9)),
expand = expansion(mult = c(0, 0))
) +
annotation_logticks(sides = "l", outside = TRUE,
short = unit(0.6,"mm"), mid = unit(0.5,"mm"), long = unit(1,"mm")) +
scale_fill_manual(values = species_palette2) +
scale_x_discrete(limits = species_order) +
theme(#axis.ticks.y = element_blank(),
axis.ticks.x = element_blank(),
axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_text(),
axis.line.y = element_line(),
axis.title.y = element_text()) +
# coord_flip(clip = 'off') +
NoLegend()
p2 <- PrettyBarplot(flip_var(nTypes.df, 'species'),
x = 'species', y = 'nReplicates', fill = 'species',
color = NA, width = 0.9) +
scale_fill_manual(values = species_palette2) +
ylab('Number of\nbiological\nreplicates') +
scale_x_discrete(limits = species_order) +
scale_y_continuous(trans = "reverse")+
theme(#axis.ticks.y = element_blank(),
axis.ticks.x = element_blank(),
axis.text.y = element_text(),
axis.line.y = element_line(),
axis.title.y = element_text(),
axis.text.x = element_blank(),
axis.title.x = element_blank()) +
# coord_flip() +
NoLegend()
p3 <- PrettyBarplot(flip_var(nTypes.df, 'species'),
x = 'species', y = 'nClusters', fill = 'species',
color = NA, width = 0.9) +
scale_fill_manual(values = species_palette2) +
scale_x_discrete(limits = species_order, position = "bottom") +
scale_y_continuous(trans = "reverse")+
ylab('Number of\nnominal\nclusters') +
# geom_hline(yintercept = -0.5, linewidth = 1) +
theme(axis.title.y = element_text(),
axis.title.x = element_blank()) +
# coord_flip() +
NoLegend()
# Combine plots vertically
p1 / p2 / p3
# ggsave('../../figures/my_figs/species-cell-cluster-quant.pdf', height = 5.4, width = 5)
```
## Manual annotation of select ACs
```{r, fig.height=8, fig.width=6}
DotPlot3(objectList$Killifish, features = c('LOC107381579', 'slc5a7a', 'isl1a', 'slc17a8', 'prox1a'), group.by = 'annotated')
# DotPlot3(objectList$Goldfish, features = toupper(c('slc5a7a', 'slc17a8', 'prox1a', 'dab1')), group.by = 'type')
DotPlot3(objectList$CatShark, features = c('chata', 'slc5a7a', 'slc17a8', 'prox1a'), group.by = 'annotated')
DotPlot3(objectList$Axolotl, features = c('CHAT', 'SLC5A7', 'SLC17A8', 'GJD2', 'PROX1', 'NFIA'), group.by = 'annotated')
DotPlot3(objectList$Newt, features = c('CHAT', 'SLC5A7', 'SLC17A8', 'GJD2', 'PROX1', 'NFIA'), group.by = 'annotated')
objectList$Axolotl = AssignAnnotations(objectList$Axolotl, list(SAC = 'gabaAC-1', VG3 = 'glyAC-12', A2 = 'glyAC-3'),
use = 'annotated')
objectList$Newt = AssignAnnotations(objectList$Newt, list(SAC = 'gabaAC-1', VG3 = 'glyAC-7', A2 = 'glyAC-21'),
use = 'annotated')
objectList$Killifish = AssignAnnotations(objectList$Killifish, list(SAC = 'AC-1'),
use = 'annotated')
objectList$CatShark = AssignAnnotations(objectList$CatShark, list(SAC = 'gabaAC-3', VG3 = 'glyAC-13'),
use = 'annotated')
# None found in goldfish by manual annotation
objectList2 = objectList[!names(objectList) %in% c('Goldfish')]
# Remove zebrafish types found by SAMap integration
objectList2$Zebrafish$lit_type[objectList2$Zebrafish$lit_type == 'A2'] = NA
# Remove rhabdomys A2 which doesn't express markers
objectList2$Rhabdomys$lit_type[objectList2$Rhabdomys$lit_type == 'A2'] = NA
# Keep lamprey SAC
objectList2$Lamprey$lit_type[objectList2$Lamprey$lit_type == 'VG3'] = NA
# marker.df = data.frame(matrix(rep(c('CHAT', 'SLC5A7', 'ISL1', 'SOX2', 'SLC18A3', 'GJD2', 'PROX1', 'NFIA', 'DAB1', 'SLC17A8', 'SDK2'), 20), ncol = 20)) %>% setNames(speciesList)
marker.df = data.frame(matrix(rep(c('CHAT', 'SLC5A7', 'ISL1',
'GJD2', 'PROX1', 'NFIA',
'SLC17A8', 'SDK2', 'NXPH1'), length(objectList2)),
ncol = length(objectList2))) %>% setNames(names(objectList2))
marker.df$row_annotation = c(rep('SAC', 3), rep('A2', 3), rep('VG3', 3))
marker.df
# Modify symbols for killifish and catshark, which were not converted
marker.df$Killifish = c('LOC107381579', 'slc5a7a', 'isl1a', 'gjd2a', 'prox1a', 'nfia', 'slc17a8', 'sdk2', 'nxph1')
marker.df$CatShark = c('chata', 'slc5a7a', 'isl1a', 'gjd2a', 'prox1a', 'nfia', 'slc17a8', 'sdk2', 'nxph1')
# Subset to annotated cells
objectList2 = lapply(objectList2, function(object) {
object$annotated = as.character(object$type)
object$annotated[object$lit_type %in% c('SAC', 'A2', 'VG3')] = object$lit_type[object$lit_type %in% c('SAC', 'A2', 'VG3')]
object
})
# species_palette = colorRampPalette(as.character(paletteer::paletteer_d("lisa::OskarSchlemmer", 5)))(19) %>% setNames(speciesList)
# species_palette = colorspace::darken(rainbow(19), 0.2) %>% setNames(speciesList)
pdf('../../figures/my_figs/manual-annotation-3.pdf', width=5, height=7, family = 'ArialMT')
SAMapHeatmap2(objectList2,
marker.df,
type_palette = rep('white', 3) %>% setNames(c('SAC', 'A2', 'VG3')),
#colorspace::lighten(AcPalette2(3) %>% setNames(c('SAC', 'A2', 'VG3')), 0.2),
#colorspace::lighten(c('blue', 'orange2', 'deeppink3'), 0.5) %>% setNames(c('SAC', 'A2', 'VG3')),
species_palette = species_palette3, #colorspace::darken(AcPalette2(length(speciesList)), 0) %>% setNames(speciesList),
species.use = names(objectList2), #names(species_ids),
types.use = c('SAC', 'A2', 'VG3'),
# font.size = 12,
show_annotation_legend = FALSE,
show_row_names = FALSE,
show_column_names = TRUE,
# pseudocount = 10,
font.family = 'ArialMT',
column_names_gp = gpar(fontface = "italic", fontsize = 10),
row_title_gp = gpar(fontsize = 30),
row_title_rot = 45,
column_names_rot = 45,
names.show = 'Human',
# width = unit(4, "in"),
rotate = TRUE,
dotplot = TRUE,
dot.scale.factor = 0.18)
dev.off()
```
## Amniote integration
```{r}
types.to.consider = c('A2', 'SAC', 'VG3')
vertebrateAC_unintegrated = readRDS("../../Ortho_Objects/vertebrateAC_unintegrated.rds")
vertebrateAC_unintegrated$species_littype = str_split_fixed(vertebrateAC_unintegrated$species_type, "_", 2)[,2]
vertebrateAC_unintegrated$species_littype[!vertebrateAC_unintegrated$species_littype %in% types.to.consider] = "Other"
vertebrateAC_unintegrated$species_littype = factor(vertebrateAC_unintegrated$species_littype, levels = c('SAC', 'A2', 'VG3', 'Other'))
vertebrateAC_integrated = readRDS("../../Ortho_Objects/vertebrateAC_seurat.rds")
vertebrateAC_integrated$species_littype = str_split_fixed(vertebrateAC_integrated$species_type, "_", 2)[,2]
vertebrateAC_integrated$species_littype[!vertebrateAC_integrated$species_littype %in% types.to.consider] = "Other"
vertebrateAC_integrated$species_littype = factor(vertebrateAC_integrated$species_littype, levels = c('SAC', 'A2', 'VG3', 'Other'))
# Remove cluster 0 (poorly integrated cells)
vertebrateAC_integrated = subset(vertebrateAC_integrated, seurat_clusters == 0, invert = TRUE)
vertebrateAC_shuffled = readRDS("../../Ortho_Objects/vertebrateAC_shuffled_row_entries.rds")
vertebrateAC_shuffled$species_littype = str_split_fixed(vertebrateAC_shuffled$species_type, "_", 2)[,2]
vertebrateAC_shuffled$species_littype[!vertebrateAC_shuffled$species_littype %in% types.to.consider] = "Other"
vertebrateAC_shuffled$species_littype = factor(vertebrateAC_shuffled$species_littype, levels = c('SAC', 'A2', 'VG3', 'Other'))
ortho_path = '../../Ortho_Objects/'
vertebrateAC_iLISI = readRDS(paste0(ortho_path, 'vertebrateAC_iLISI.rds'))
vertebrateAC_seurat_iLISI = readRDS(paste0(ortho_path, 'vertebrateAC_seurat_iLISI.rds'))
shuffledAC_iLISI = readRDS(paste0(ortho_path, 'shuffledAC_iLISI.rds'))
vertebrateAC_cLISI = readRDS(paste0(ortho_path, 'vertebrateAC_cLISI.rds'))
vertebrateAC_seurat_cLISI = readRDS(paste0(ortho_path, 'vertebrateAC_seurat_cLISI.rds'))
shuffledAC_cLISI = readRDS(paste0(ortho_path, 'shuffledAC_cLISI.rds'))
```
### Simple for main figure
```{r, fig.height=4, fig.width=12}
# types.to.consider = c('A2', 'SAC', 'VG3')
# vertebrateAC_unintegrated = readRDS("../../Ortho_Objects/vertebrateAC_unintegrated.rds")
# vertebrateAC_unintegrated$species_littype = str_split_fixed(vertebrateAC_unintegrated$species_type, "_", 2)[,2]
# vertebrateAC_unintegrated$species_littype[!vertebrateAC_unintegrated$species_littype %in% types.to.consider] = "Other"
# vertebrateAC_unintegrated$species_littype = factor(vertebrateAC_unintegrated$species_littype, levels = c('SAC', 'A2', 'VG3', 'Other'))
#
# vertebrateAC_integrated = readRDS("../../Ortho_Objects/vertebrateAC_seurat.rds")
# vertebrateAC_integrated$species_littype = str_split_fixed(vertebrateAC_integrated$species_type, "_", 2)[,2]
# vertebrateAC_integrated$species_littype[!vertebrateAC_integrated$species_littype %in% types.to.consider] = "Other"
# vertebrateAC_integrated$species_littype = factor(vertebrateAC_integrated$species_littype, levels = c('SAC', 'A2', 'VG3', 'Other'))
#
# ortho_path = '../../Ortho_Objects/'
# shuffledAC_cLISI = readRDS(paste0(ortho_path, 'shuffledAC_cLISI.rds'))
# shuffledAC_cLISI$species_littype = factor(shuffledAC_cLISI$species_littype, levels = c('SAC', 'A2', 'VG3', 'Other'))
umap.params = list(shuffle = TRUE, raster.dpi = 1000, pad = 0.1, pt.size = 0.001, pt.alpha = 1, nbreaks = 30, rasterise = TRUE)
theme_umap(
do.call(PrettyUmap2, c(list(vertebrateAC_unintegrated,
group.by = "species",
color.by = 'species_littype',
title = 'Unintegrated',
cols = lighten(c(AcPalette3(3), 'grey'), 0.3),
label = TRUE),
umap.params))
) + NoLegend() |
theme_umap(
do.call(PrettyUmap2, c(list(vertebrateAC_integrated,
group.by = "species_littype",
# color.by = "species_littype",
title = 'Integrated',
cols = lighten(c(AcPalette3(3), 'grey'), 0.3),
label = TRUE,
geom.label = geom_text_repel,
box.padding = 1), umap.params))
) + NoLegend() |
theme_umap(
do.call(PrettyUmap2, c(list(vertebrateAC_shuffled,
group.by = "species_littype",
# color.by = "species_littype",
title = 'Shuffled + integrated',
cols = lighten(c(AcPalette3(3), 'grey'), 0.3),
label = FALSE), umap.params))
) + NoLegend()
# ggsave("../../figures/my_figs/fig-tetrapod-integration-simple.pdf", width=8, height=4)
# ggsave("../../figures/my_figs/fig-tetrapod-integration-simple.png", width=8, height=4)
ggsave("../../figures/my_figs/fig-tetrapod-integration-supp.png", width=11, height=4)
```
### Just mammals
```{r, fig.height=4, fig.width=8}
ac.ortho = LoadACOrtho()
ac.ortho$orthotype.no = gsub('oAC|\\*', '', ExtractString(ac.ortho$orthotype, after = ' '))
# ac.ortho.no.tetrapods = subset(ac.ortho, species %in% c("Chicken", "Lizard"), invert = TRUE)
types.to.consider = c('A2', 'SAC', 'VG3')
ac.ortho$species_littype[!ac.ortho$species_littype %in% types.to.consider] = "Other"
ac.ortho$species_littype = factor(ac.ortho$species_littype, levels = c('SAC', 'A2', 'VG3', 'Other'))
umap.params = list(shuffle = TRUE, raster.dpi = 1000, pad = 0.1,
pt.size = 0.001, pt.alpha = 1, nbreaks = 30, rasterise = TRUE)
theme_umap(
do.call(PrettyUmap2, c(list(ac.ortho,
group.by = "orthotype.no",
color.by = "species",
cols = species_palette2,
geom.label = geom_label,
color = 'black', fill = 'white',
label.size = 0, alpha = 0.5,
label.padding = unit(0, "lines"),
label.r = unit(0.2, "lines"),
label = TRUE), umap.params))
) + NoLegend() + theme(plot.title = element_blank()) |
theme_umap(
do.call(PrettyUmap2, c(list(ac.ortho,
group.by = "species_littype",
# color.by = "species_littype",
cols = lighten(c(AcPalette3(3), 'grey'), 0.3),
label = TRUE), umap.params))
) + NoLegend() + theme(plot.title = element_blank())
# ggsave("../../figures/my_figs/fig-tetrapod-integration-mammals.pdf", width=8, height=4)
ggsave("../../figures/my_figs/fig-tetrapod-integration-amniotes.pdf", width=8, height=4)
```
### Without shuffled control
```{r, fig.width=9, fig.height=7.5}
# raster.dpi = 1024
separator = 0
# nbreaks = 30
umap.params = list(shuffle = TRUE, raster.dpi = 1000, pad = 0.1, pt.size = 0.001, pt.alpha = 1, nbreaks = 30, rasterise = TRUE)
LISI_umap <- plot_grid(
theme_umap(
do.call(PrettyUmap2, c(list(vertebrateAC_unintegrated, group.by = "species", label = TRUE, cols = colorspace::lighten(species_palette2, 0.3)), umap.params)),
title = "Unintegrated"
) + NoLegend(), #+ theme(plot.title = element_blank()),
NULL,
theme_umap(
do.call(PrettyUmap2, c(list(vertebrateAC_integrated, group.by = "species", label = FALSE, cols = colorspace::lighten(species_palette2, 0.3)), umap.params)),
title = "Integrated"
) + NoLegend(), #+ theme(plot.title = element_blank()),
theme_umap(
do.call(PrettyUmap2, c(list(vertebrateAC_unintegrated, group.by = "species_littype", cols = lighten(c(AcPalette3(3), 'grey'), 0.3), label = FALSE), umap.params))
) + NoLegend() + theme(plot.title = element_blank()),
NULL,
theme_umap(
do.call(PrettyUmap2, c(list(vertebrateAC_integrated, group.by = "species_littype", cols = lighten(c(AcPalette3(3), 'grey'), 0.3), label = TRUE), umap.params))
) + NoLegend() + theme(plot.title = element_blank()),
nrow = 2, ncol = 3,
rel_heights = c(1.1, 1),
rel_widths = c(1, separator, 1, separator, 1),
align = 'v', axis = 'lr'
)
# LISI_umap
iLISI.df = data.frame(unintegrated = vertebrateAC_iLISI$normalized.iLISI,
integrated = vertebrateAC_seurat_iLISI$normalized.iLISI,
shuffled = shuffledAC_iLISI$normalized.iLISI)
cLISI.df = data.frame(unintegrated = subset(vertebrateAC_cLISI, species_littype %in% c("A2", "SAC", "VG3"))$normalized.cLISI,
integrated = subset(vertebrateAC_seurat_cLISI, species_littype %in% c("A2", "SAC", "VG3"))$normalized.cLISI,
shuffled = subset(shuffledAC_cLISI, species_littype %in% c("A2", "SAC", "VG3"))$normalized.cLISI)
boxplots = plot_grid(
NULL,
ggplot(reshape2::melt(iLISI.df), aes(x = variable, y = value, color = variable))+
geom_boxplot(linetype = "dashed", outlier.shape = NA) +
stat_boxplot(aes(ymin = ..lower.., ymax = ..upper..), outlier.shape = NA) +
stat_boxplot(geom = "errorbar", aes(ymin = ..ymax..), width = 0.5) +
stat_boxplot(geom = "errorbar", aes(ymax = ..ymin..), width = 0.5) +
ylab("Species mixing (iLISI)")+
# ggtitle("Mixing of species")+
theme_cowplot()+
theme(axis.text.x = element_blank(), legend.position = "none", axis.title.x = element_blank()),
NULL,
ggplot(reshape2::melt(cLISI.df), aes(x = variable, y = value, color = variable))+
geom_boxplot(linetype = "dashed", outlier.shape = NA) +
stat_boxplot(aes(ymin = ..lower.., ymax = ..upper..), outlier.shape = NA) +
stat_boxplot(geom = "errorbar", aes(ymin = ..ymax..), width = 0.5) +
stat_boxplot(geom = "errorbar", aes(ymax = ..ymin..), width = 0.5) +
ylab("Cell type mixing (cLISI)")+
xlab(NULL)+
# ggtitle("Mixing of cell types")+
theme_cowplot()+
theme(legend.position = "none")+
RotatedAxis(),
nrow = 4, rel_heights = c(0.2, 1, 0.2, 1.25), align = "v"
)
lisi_summary = plot_grid(
LISI_umap,
# plot_grid(iLISI_umap, plot_grid(cLISI_umap, NULL, nrow = 1, rel_widths = c(1,0.03)), nrow = 2, rel_heights = c(1,1)),
NULL,
boxplots,
ncol = 3,
rel_widths = c(4, 0.1, 1))
lisi_summary
# ggsave("../../figures/my_figs/LISI-summary-cairo.pdf", width=14.5, height=7.5, device = cairo_pdf)
ggsave("../../figures/my_figs/LISI-summary-small.pdf", width=11, height=7.5)
```
### Just boxplots
```{r, fig.height=4, fig.width=4}
#' Run amniote integration chunk to load files
iLISI.df = data.frame(Unintegrated = vertebrateAC_iLISI$normalized.iLISI,
Integrated = vertebrateAC_seurat_iLISI$normalized.iLISI,
Shuffled = shuffledAC_iLISI$normalized.iLISI)
cLISI.df = data.frame(Unintegrated = subset(vertebrateAC_cLISI, species_littype %in% c("A2", "SAC", "VG3"))$normalized.cLISI,
Integrated = subset(vertebrateAC_seurat_cLISI, species_littype %in% c("A2", "SAC", "VG3"))$normalized.cLISI,
Shuffled = subset(shuffledAC_cLISI, species_littype %in% c("A2", "SAC", "VG3"))$normalized.cLISI)
plot_grid(
ggplot(reshape2::melt(iLISI.df), aes(x = factor(variable, levels = rev(c('Unintegrated', 'Integrated', 'Shuffled'))), y = value, color = variable))+
geom_boxplot(linetype = "dashed", outlier.shape = NA) +
stat_boxplot(aes(ymin = ..lower.., ymax = ..upper..), outlier.shape = NA) +
stat_boxplot(geom = "errorbar", aes(ymin = ..ymax..), width = 0.5) +
stat_boxplot(geom = "errorbar", aes(ymax = ..ymin..), width = 0.5) +
scale_color_manual(values = c(Shuffled = '#4646FD', Integrated = '#A63FE4', Unintegrated = 'darkgrey'))+
ylab("Species mixing (iLISI)")+
# ggtitle("Mixing of species")+
theme_cowplot()+
coord_flip()+
RotatedAxis()+
theme(legend.position = "none", axis.title.y = element_blank()),
ggplot(reshape2::melt(cLISI.df), aes(x = factor(variable, levels = rev(c('Unintegrated', 'Integrated', 'Shuffled'))), y = value, color = variable))+
geom_boxplot(linetype = "dashed", outlier.shape = NA) +
stat_boxplot(aes(ymin = ..lower.., ymax = ..upper..), outlier.shape = NA) +
stat_boxplot(geom = "errorbar", aes(ymin = ..ymax..), width = 0.5) +
stat_boxplot(geom = "errorbar", aes(ymax = ..ymin..), width = 0.5) +
scale_color_manual(values = c(Shuffled = '#4646FD', Integrated = '#A63FE4', Unintegrated = 'darkgrey'))+
ylab("Cell type mixing (cLISI)")+
xlab(NULL)+
# ggtitle("Mixing of cell types")+
theme_cowplot()+
theme(legend.position = "none")+
coord_flip()+
RotatedAxis(),
ncol = 1,
align = "v", axis = 'lr'
)
ggsave("../../figures/my_figs/LISI-boxplots.pdf", width=3.5, height=4)
```
### With shuffled control
```{r, fig.width=14.5, fig.height=7.5}
# raster.dpi = 1024
separator = -0.3
# nbreaks = 30
umap.params = list(shuffle = TRUE, raster.dpi = 1000, pad = 0.1, pt.size = 0.001, pt.alpha = 1, nbreaks = 30, label = FALSE, rasterise = TRUE)
LISI_umap <- plot_grid(
theme_umap(
do.call(PrettyUmap2, c(list(vertebrateAC_unintegrated, group.by = "species"), umap.params)),
title = "Unintegrated"
) + NoLegend(), #+ theme(plot.title = element_blank()),
NULL,
theme_umap(
do.call(PrettyUmap2, c(list(vertebrateAC_integrated, group.by = "species"), umap.params)),
title = "Integrated"
) + NoLegend(), #+ theme(plot.title = element_blank()),
NULL,
theme_umap(
do.call(PrettyUmap2, c(list(vertebrateAC_shuffled, group.by = "species", show.legend = TRUE), umap.params)),
title = "Integrated & shuffled"
) + theme(legend.key.height = unit(1, "lines")) +
guides(color = guide_legend(override.aes = list(size = 3))),
theme_umap(
do.call(PrettyUmap2, c(list(vertebrateAC_unintegrated, group.by = "species_littype", cols = c(AcPalette2(4)[1:3], 'grey')), umap.params))
) + NoLegend() + theme(plot.title = element_blank()),
NULL,
theme_umap(
do.call(PrettyUmap2, c(list(vertebrateAC_integrated, group.by = "species_littype", cols = c(AcPalette2(4)[1:3], 'grey')), umap.params))
) + NoLegend() + theme(plot.title = element_blank()),
NULL,
theme_umap(
do.call(PrettyUmap2, c(list(vertebrateAC_shuffled, group.by = "species_littype", show.legend = TRUE, cols = c(AcPalette2(4)[1:3], 'grey')), umap.params))
) + theme(plot.title = element_blank(), legend.key.height = unit(1, "lines")) +
guides(color = guide_legend(override.aes = list(size = 3))),
nrow = 2, ncol = 5,
rel_heights = c(1.1, 1),
rel_widths = c(1, separator, 1, separator, 1),
align = 'v', axis = 'lr'
)
# LISI_umap
iLISI.df = data.frame(unintegrated = vertebrateAC_iLISI$normalized.iLISI,
integrated = vertebrateAC_seurat_iLISI$normalized.iLISI,
shuffled = shuffledAC_iLISI$normalized.iLISI)
cLISI.df = data.frame(unintegrated = subset(vertebrateAC_cLISI, species_littype %in% c("A2", "SAC", "VG3"))$normalized.cLISI,
integrated = subset(vertebrateAC_seurat_cLISI, species_littype %in% c("A2", "SAC", "VG3"))$normalized.cLISI,
shuffled = subset(shuffledAC_cLISI, species_littype %in% c("A2", "SAC", "VG3"))$normalized.cLISI)
boxplots = plot_grid(
NULL,
ggplot(reshape2::melt(iLISI.df), aes(x = variable, y = value, color = variable))+
geom_boxplot(linetype = "dashed", outlier.shape = NA) +
stat_boxplot(aes(ymin = ..lower.., ymax = ..upper..), outlier.shape = NA) +
stat_boxplot(geom = "errorbar", aes(ymin = ..ymax..), width = 0.5) +
stat_boxplot(geom = "errorbar", aes(ymax = ..ymin..), width = 0.5) +
ylab("Species mixing (iLISI)")+
# ggtitle("Mixing of species")+
theme_cowplot()+
theme(axis.text.x = element_blank(), legend.position = "none", axis.title.x = element_blank()),
NULL,
ggplot(reshape2::melt(cLISI.df), aes(x = variable, y = value, color = variable))+
geom_boxplot(linetype = "dashed", outlier.shape = NA) +
stat_boxplot(aes(ymin = ..lower.., ymax = ..upper..), outlier.shape = NA) +
stat_boxplot(geom = "errorbar", aes(ymin = ..ymax..), width = 0.5) +
stat_boxplot(geom = "errorbar", aes(ymax = ..ymin..), width = 0.5) +
ylab("Cell type mixing (cLISI)")+
xlab(NULL)+
# ggtitle("Mixing of cell types")+
theme_cowplot()+
theme(legend.position = "none")+
RotatedAxis(),
nrow = 4, rel_heights = c(0.2, 1, 0.2, 1.25), align = "v"
)
lisi_summary = plot_grid(
LISI_umap,
# plot_grid(iLISI_umap, plot_grid(cLISI_umap, NULL, nrow = 1, rel_widths = c(1,0.03)), nrow = 2, rel_heights = c(1,1)),
NULL,
boxplots,
ncol = 3,
rel_widths = c(6, 0.1, 1))
# lisi_summary
# ggsave("../../figures/my_figs/LISI-summary-cairo.pdf", width=14.5, height=7.5, device = cairo_pdf)
ggsave("../../figures/my_figs/LISI-summary.pdf", width=14.5, height=7.5)
```
```{r LISI-summary, fig.width=15, fig.height=8.5}
ggsave("../../figures/my_figs/LISI-summary-cairo.pdf", width=14.5, height=7.5, device = cairo_pdf)
ggsave("../../figures/my_figs/LISI-summary.pdf", width=14.5, height=7.5)
```
## AC orthotypes
```{r}
ac.ortho = LoadACOrtho()
```
```{r, fig.height=5, fig.width=5}
# pdf('../../figures/my_figs/umap-orthotypes.pdf', height = 4, width = 4.5) # family="ArialMT")
TitlePlot(PrettyUmap2(ac.ortho, group.by = 'type_no', nbreaks = 30, pt.size = 0.2,
raster.dpi = 1000, cols = type_cols2.no,
geom.label = geom_label, color = 'black', fill = 'white',
label.size = 0, alpha = 0.2, label.padding = unit(0, "lines"),
label.r = unit(0.2, "lines")), #label.size = NA, label.r = unit(0, "lines"))
title = '42 unsupervised clusters')
# dev.off()
```
## Species composition
```{r, fig.height=5, fig.width=9}
# stackedBarGraph2(ac.ortho, 'type', 'species', border.color = 'black') +
# theme(text = element_text(family = 'ArialMT'), axis.text.x = element_text(size = 8)) +
# ArialFont()+
# scale_fill_manual(name = 'Species', values = lighten(species_palette2, 0.2))
# ggsave('../../figures/my_figs/figS_species-composition-1.pdf', height=2.5, width=9)
```
```{r, fig.height=6, fig.width=12}
# ac.ortho$orthotype_no = factor(as.character((gsub('\\*', '', ExtractString(ac.ortho$orthotype, after = ' ')))),
# levels = paste0('oAC', 1:42))
# coord_flip_discrete(
# TitlePlot(stackedBarGraph2(ac.ortho, 'orthotype_no', 'species', border.color = 'black'), 'Proportion') +
# # scale_x_discrete(labels = function(x) str_split_fixed(gsub('*', '', x), '_', 2)[,1])+
# theme(text = element_text(family = 'ArialMT'), axis.text.y = element_text(size = 10)) +
# ArialFont()+
# xlab(NULL) + ylab(NULL)+
# scale_fill_manual(name = 'Species', values = lighten(species_palette2, 0.2)), discrete_axis = 'Var1')
# stackedBarGraph2(ac.ortho, 'type', 'species', border.color = NA) +
# # scale_x_discrete(labels = function(x) str_split_fixed(gsub('*', '', x), '_', 2)[,1])+
# theme(text = element_text(family = 'ArialMT'), axis.text.y = element_text(size = 10)) +
# ArialFont()+
# xlab(NULL) + ylab('y')+
# scale_fill_manual(name = 'Species', values = lighten(species_palette2, 0.2)) + coord_flip()
TitlePlot(stackedBarGraph2(ac.ortho, 'orthotype', 'species', border.color = 'black'), 'Species composition') +
# scale_x_discrete(labels = function(x) str_split_fixed(gsub('*', '', x), '_', 2)[,1])+
theme(text = element_text(family = 'ArialMT'), axis.text.y = element_text(size = 10)) +
ArialFont()+
xlab(NULL) +
ylab('Percent')+
theme(legend.position = 'bottom')+
scale_fill_manual(name = 'Species', values = lighten(species_palette2, 0.2))
ggsave('../../figures/my_figs/figS_species-composition-1.pdf')
```
## Precision and recall across species
```{r, fig.height=6, fig.width=6}
# ac.ortho = readRDS('../../Ortho_Objects/vertebrateAC_BIGSEURAT.rds')
vertebrateAC = ac.ortho
vertebrateAC = AssignAnnotations(ac.ortho, list(SAC = 8, A2 = 22, VG3 = 27), use = 'type_no')
# Ferret SACs post hoc
[email hidden][WhichCells(PositiveCell(subset(vertebrateAC, species == 'Ferret'), features = 'SLC5A7', return.object = TRUE), expr = positive),'species_littype'] = 'SAC'
raw_tables = lapply(speciesList.liz, function(this.species) {
lapply(KNOWNTYPES, function(type) {
# if(this.species == 'Ferret' & type == 'SAC') {
# matrix(NA, nrow = 2, ncol = 2)
# } else {
message(paste0('working on ', this.species, ' and ', type))
lit_types = vertebrateAC$lit_type[vertebrateAC$species == this.species]
species_littype = vertebrateAC$species_littype[vertebrateAC$species == this.species]
lit_types[is.na(lit_types)] = 'other'
raw_table = table(lit_types == type, species_littype == type)[2:1,2:1]
raw_table
# }
})
})
res = lapply(unlist(raw_tables, recursive = FALSE), EvaluateModel)
res
# make a matrix
precision = matrix(sapply(res, function(x) x$precision),
nrow = length(speciesList.liz),
ncol = 3,
dimnames = list(speciesList.liz, KNOWNTYPES),
byrow = TRUE)
recall = matrix(sapply(res, function(x) x$recall),
nrow = length(speciesList.liz),
ncol = 3,
dimnames = list(speciesList.liz, KNOWNTYPES),
byrow = TRUE)
# barplot(t(precision))
# barplot(t(recall))
metrics.df = rbind(reshape2::melt(precision) %>% mutate(group = 'precision'),
reshape2::melt(recall) %>% mutate(group = 'recall')) %>%
setNames(c('Species', 'Type', 'Value', 'Metric'))
# coord_flip_discrete(ggbarplot(subset(metrics.df, Metric == 'precision'), x = 'Species', y = 'Value', fill = 'Type', position = position_dodge()), discrete_axis = 'Species')
metrics.df$Species = factor(factor(metrics.df$Species, levels = rev(phylogenetic_order)))
# Subset to mammals
# metrics.df = subset(metrics.df, Species %in% species_metadata$Species[1:15]) # leave all species in
TitlePlot(ggbarplot(subset(metrics.df, Metric == 'precision'),
x = 'Species', y = 'Value', fill = 'Type', color = 'black',
position = position_dodge(width = 0.7)), title = 'Precision') +
coord_flip() + NoLegend() + #RotatedAxis() +
ylab(NULL) + xlab(NULL) +
scale_y_continuous(expand = expansion(mult = c(0,0.1)), labels = c(0, '', '', '', 1)) |
TitlePlot(ggbarplot(subset(metrics.df, Metric == 'recall'), x = 'Species', y = 'Value', fill = 'Type',
position = position_dodge(width = 0.7), legend = "right", color = 'black'), 'Recall') +
coord_flip() + #RotatedAxis() +
ylab(NULL) + scale_y_continuous(expand = expansion(mult = c(0,0.1)), labels = c(0, '', '', '', 1)) +
ArialFont() + theme(axis.text.y = element_blank(), axis.title.y = element_blank()) #+NoLegend()
ggsave('../../figures/my_figs/fig-precision-recall.pdf', height = 6, width = 4.5)
# Confidence intervals
paste0(round(t.test(subset(metrics.df, Metric == 'precision')$Value)$conf.int, 2), collapse = '-')
paste0(round(t.test(subset(metrics.df, Metric == 'recall')$Value)$conf.int, 2), collapse = '-')
```
Averaged across species
```{r}
vertebrateAC = AssignAnnotations(vertebrateAC, list(SAC = 8, A2 = 22, VG3 = 27), use = 'type_no')
KNOWNTYPES = c('SAC', 'A2', 'VG3')
# Raw tabulation
lit_types = vertebrateAC$lit_type
lit_types[is.na(lit_types)] = 'other'
raw_tables = lapply(KNOWNTYPES, function(type) {
raw_table = table(lit_types == type, vertebrateAC$species_littype == type)[2:1,2:1]
return(raw_table)
})
lapply(raw_tables, EvaluateModel)
```
Averaged across types
```{r}
raw_tables = lapply(speciesList.liz, function(this.species) {
# lapply(KNOWNTYPES, function(type) {
message(paste0('working on ', this.species))
lit_types = vertebrateAC$lit_type[vertebrateAC$species == this.species]
species_littype = vertebrateAC$species_littype[vertebrateAC$species == this.species]
lit_types[is.na(lit_types)] = 'other'
raw_table = table(lit_types == type, species_littype == type)[2:1,2:1]
raw_table
# })
})
```
## Conserved genes in primates, rodents, and laurasiatherans
```{r, fig.height=3, fig.width=12}
primate.ortho = DownsampleSeurat(subset(ac.ortho, species %in% subset(species_metadata, Clade == 'Primate')$Species),
group.by = 'orthotype', size = 50)
rodent.ortho = DownsampleSeurat(subset(ac.ortho, species %in% subset(species_metadata, Clade == 'Rodent')$Species),
group.by = 'orthotype', size = 50)
laurasia.ortho = DownsampleSeurat(subset(ac.ortho, species %in% subset(species_metadata, Clade == 'Laurasiatheria')$Species),
group.by = 'orthotype', size = 50)
ac.ortho.sub = DownsampleSeurat(ac.ortho, group.by = 'orthotype', size = 500)
de_all = TopNDEGs(ac.ortho.sub, group.by = 'orthotype', n = 100, sort.by = )
primate.ortho = ScaleData(primate.ortho, features = unique(de_all$gene))
rodent.ortho = ScaleData(rodent.ortho, features = unique(de_all$gene))
laurasia.ortho = ScaleData(laurasia.ortho, features = unique(de_all$gene))
SeuratHeatmap(ShuffleObject(primate.ortho, group.by = 'orthotype'), features = de_all$gene, group.by = 'orthotype',
show_row_names = FALSE, show_column_names = FALSE, color = c('white', 'white', 'deeppink'))
SeuratHeatmap(rodent.ortho, features = de_all$gene, group.by = 'orthotype',
show_row_names = FALSE, show_column_names = FALSE, color = c('white', 'white', 'orange2')) +
SeuratHeatmap(laurasia.ortho, features = de_all$gene, group.by = 'orthotype',
show_row_names = FALSE, show_column_names = FALSE, color = c('white', 'white', 'chartreuse3'))
```
# Figure 2: classification and annotation of ACs
## GABA gly scores
```{r, fig.height=2.5, fig.width=2.5}
BIGSEURAT = readRDS('../../Ortho_Objects/vertebrateAC_BIGSEURAT.rds')
# Check OT presence in species
table(BIGSEURAT$species, BIGSEURAT$type)
BIGSEURAT$species_type2 = paste0(BIGSEURAT$species, '-', BIGSEURAT$species_type)
X = as.matrix(ConfusionMatrix(BIGSEURAT$type, BIGSEURAT$species_type2, plot = FALSE))
dim(X)
colnames(X)
class.meta = readRDS('metadata/GabaGlyMetadata.rds')
class.meta = dapply(seq_along(class.meta), function(i) {
class.meta[[i]]$species_type2 = paste0(names(class.meta)[[i]], '-', class.meta[[i]]$type)
return(class.meta[[i]])
})
# names(class.meta) = names(SPECIESFILES)
Vdata = do.call(rbind, class.meta[!names(class.meta) %in% c('Zebrafish', 'Goldfish', 'Lamprey')])
dim(Vdata)
data.frame(c(colnames(X), rep(0, length(Vdata$type) - length(colnames(X)))), Vdata$species_type2)
# Two lizard clusters are missing as expected (HC contamination)
VennDiagram(colnames(X), Vdata$species_type2)
# V = as.matrix(data.frame(pGly = Vdata$Gly_posterior[match(colnames(X), Vdata$species_type2)],
# pGaba = Vdata$GABA_posterior[match(colnames(X), Vdata$species_type2)]))
V = as.matrix(data.frame(`Gly_score` = Vdata$Gly_score[match(colnames(X), Vdata$species_type2)],
`GABA_score` = Vdata$GABA_score[match(colnames(X), Vdata$species_type2)]))
# Check dimensions
dim(X)
dim(V)
U = as.data.frame(X %*% V)
U$cluster = rownames(U)
U$Subclass = 'GABAergic'
U$Subclass[U$GABA_score < 0.2] = 'Glycinergic'
U$cluster_no = convert_values(U$cluster, Metadata(ac.ortho, 'type', 'type_no') %>% setNames(c('old.names', 'new.names')))
U$orthotype.no = orthotype_no_labels(U$cluster_no)
ScatterPlot(U, 'GABA_score', 'Gly_score', fill = 'Subclass', r_pval = FALSE, lm = FALSE, labels = 'orthotype.no',
max.overlaps = 10, xlab = 'GABAergic score', ylab = 'Glycinergic score') +
scale_fill_manual(values = c('cyan', 'cyan4') %>% setNames(c('GABAergic', 'Glycinergic'))) +
LegendTopRight()+
# LegendLowerRight(x = 0.99, y = 0.53) +
geom_vline(xintercept = 0.2, linetype = 'dashed', color = 'black', alpha = 0.2) +
geom_hline(yintercept = 0.2, linetype = 'dashed', color = 'black', alpha = 0.2)
ggsave('../../figures/my_figs/gaba-gly-scores.pdf', width = 2.5, height = 2.5)
```
## nGnGs from mouse have lowest scores
```{r, fig.height=3, fig.width=3}
# c('16_NNgaba*', '2_NNgly')
# vec1 = as.character(ac.ortho$type)
# vec1[!vec1 %in% c('16_NNgaba*', '2_NNgly')] = 'Other'
# vec2 = as.character(ac.ortho$mouse.annotated)
# vec2[!vec2 %in% c('36', '10_CCK', '24', '30')] = 'Other'
# nn_matrix = JSMatrix(table(vec1, vec2))
nn_matrix = JSMatrix(table(ifelse(ac.ortho$type == '16_NNgaba*', '16_NNgaba*', ifelse(ac.ortho$type == '2_NNgly', '2_NNgly', 'Other')), ac.ortho$mouse.annotated))[, c('36', '10_CCK', '24', '30')]
colnames(nn_matrix) = ExtractString(colnames(nn_matrix), after = '_')
rownames(nn_matrix) = c('oAC36', 'oAC8', 'Other') #as.character((rownames(nn_matrix)))
rownames(nn_matrix)[is.na(rownames(nn_matrix))] = 'Other'
JSHeatmap(t(nn_matrix), title = 'nGnG oACs') +
ylab('Yan et al. 2020')+
ArialFont()+
theme(plot.margin = unit(c(1,1,1,1.3),"cm"))+
xlab(NULL)
ggsave('../../figures/my_figs/fig-ngng-yan-comparison.pdf', height=2.8, width=3)
```
```{r, fig.height=4, fig.width=9}
dendro = readRDS('../../Ortho_Objects/dendro.hs.liz.rds')
gaba.scores = readRDS('../../Ortho_Objects/gaba-scores.rds')
gly.scores = readRDS('../../Ortho_Objects/gly-scores.rds')
bar.col = 'black'
gly.df = data.frame(`NNgly` = gly.scores[,'2_NNgly'],
`NNgaba` = gly.scores[,'16_NNgaba*']) %>% rownames_to_column('Species')
gaba.df = data.frame(`NNgly` = gaba.scores[,'2_NNgly'],
`NNgaba` = gaba.scores[,'16_NNgaba*']) %>% rownames_to_column('Species')
# Define individual plots
p1_title <- TitlePlot(dendro, orthotype_labels('2_NNgly')) +
theme(panel.background = element_rect(fill = "transparent", color = NA),
plot.background = element_rect(fill = "transparent", color = NA))
p1_top <- PrettyBarplot(gly.df, x = 'Species', y = 'NNgly',
fill = 'cyan4', color = bar.col) +
# geom_hline(yintercept = 0.2, linetype = 'dashed', color = 'grey') +
scale_y_reverse(name = 'Glycinergic', limits = c(1, 0)) +
theme(axis.text.x = element_blank(),
axis.title.x = element_blank(),
axis.line.x = element_blank(),
panel.background = element_rect(fill = "transparent", color = NA),
plot.background = element_rect(fill = "transparent", color = NA),
axis.ticks.x = element_blank())
p1_bottom <- PrettyBarplot(gaba.df, x = 'Species', y = 'NNgly',
fill = 'cyan', color = bar.col) +
# geom_hline(yintercept = 0.2, linetype = 'dashed', color = 'grey') +
RotatedAxis() +
scale_y_reverse(name = 'GABAergic', limits = c(1, 0))
# Left column: Glycinergic score
left_col <- plot_grid(p1_title + scale_x_reverse(),
NULL,
p1_top + theme(plot.margin = unit(c(0,0,0,0), 'cm'),
panel.background = element_rect(fill = "transparent", color = NA),
plot.background = element_rect(fill = "transparent", color = NA)),
p1_bottom + theme(axis.title.x = element_blank(),
panel.background = element_rect(fill = "transparent", color = NA),
plot.background = element_rect(fill = "transparent", color = NA)),
ncol = 1,
align = "v",
rel_heights = c(0.4, -0.11, 0.3, 0.7))
# Right column: GABAergic score
p2_title <- TitlePlot(dendro, orthotype_labels('16_NNgaba*'))+
theme(panel.background = element_rect(fill = "transparent", color = NA),
plot.background = element_rect(fill = "transparent", color = NA))
p2_top <- PrettyBarplot(gly.df, x = 'Species', y = 'NNgaba',
fill = 'cyan4', color = bar.col) +
# geom_hline(yintercept = 0.2, linetype = 'dashed', color = 'grey') +
scale_y_reverse(limits = c(1, 0)) +
theme(axis.text.x = element_blank(),
axis.title.x = element_blank(),
panel.background = element_rect(fill = "transparent", color = NA),
plot.background = element_rect(fill = "transparent", color = NA),
axis.line.x = element_blank(),
axis.ticks.x = element_blank())
p2_bottom <- add_rectangle_annotation(
PrettyBarplot(gaba.df, x = 'Species', y = 'NNgaba', fill = c('cyan'), color = bar.col), index = 7, axis = 'row') +
RotatedAxis()+
# geom_hline(yintercept = 0.2, linetype = 'dashed', color = 'grey') +
scale_y_reverse(limits = c(1, 0))
# Left column: Glycinergic score
right_col <- plot_grid(p2_title + scale_x_reverse(),
NULL,
p2_top + theme(plot.margin = unit(c(0,0,0,0), 'cm'),
panel.background = element_rect(fill = "transparent", color = NA),
plot.background = element_rect(fill = "transparent", color = NA),
axis.text.y = element_blank(),
axis.title.y = element_blank()),
p2_bottom + theme(
panel.background = element_rect(fill = "transparent", color = NA),
plot.background = element_rect(fill = "transparent", color = NA),
axis.text.y = element_blank(),
axis.title.y = element_blank(),
axis.title.x = element_blank()),
ncol = 1,
align = "v",
rel_heights = c(0.4, -0.11, 0.3, 0.7))
# Combine both columns
final_plot <- plot_grid(left_col, right_col, ncol = 2, rel_widths = c(1,0.9))
final_plot
ggsave('../../figures/my_figs/fig-ngng-scores.pdf', height=3.5, width=8)
```
Dendrogram
```{r phylogeny-ncells-barplot, fig.height=6, fig.width=8}
library(ggdendro)
# Convert to common name
key = data.table::fread("../../Evolution/phylogeny/species_common_latin.txt")
tree = ape::read.tree("../../Evolution/phylogeny/20_species.nwk")
tree$tip.label = key$CommonName[match(gsub("_", " ", tree$tip.label), key$LatinName)]
dendro = phylogram::as.dendrogram.phylo(tree)
sorted.dendro = reorder(dendro, 1:20)
dendro = ggdendrogram(dendextend::prune(sorted.dendro, c("Zebrafish", 'Lamprey', 'Goldfish')))+
theme(axis.text.y = element_blank(), axis.text.x = element_blank())
dendro
saveRDS(dendextend::prune(sorted.dendro, c("Zebrafish", 'Lamprey', 'Goldfish')), '../../Ortho_Objects/phylo.hs.liz.rds')
saveRDS(dendro, '../../Ortho_Objects/dendro.hs.liz.rds')
ggsave('../../figures/my_figs/17-species-phylo.pdf', height = 1, width = 5)
# ylab("Evolutionary distance (MYA)")
# theme(axis.text.x = element_text(size = 13, hjust = 1, vjust = 1, angle = 45), axis.ticks.length.y = unit(.25, "cm"))+
# scale_y_continuous(breaks = seq(0, 500, len = 6))+ #expand = expansion(add = c(0,0.1)))+
# scale_x_discrete(expand = expansion(add = 0.6))+
# theme_cowplot()
# theme(axis.line.x=element_blank(),
# # axis.text.x=element_blank(),
# axis.ticks.x=element_blank(),
# axis.title.x=element_blank(),
# panel.grid.minor.x=element_blank(),
# panel.grid.major.x=element_blank())
```
New phylogeny
```{r, fig.height=5, fig.width=4}
library(ggdendro)
# Convert to common name
key = data.table::fread("../../Evolution/phylogeny/species_common_latin.txt")
tree = ape::read.tree("../../Evolution/phylogeny/22_species.nwk")
tree$tip.label = key$CommonName[match(gsub("_", " ", tree$tip.label), key$LatinName)]
dendro = phylogram::as.dendrogram.phylo(tree)
sorted.dendro = reorder(dendro, 1:20)
sorted.dendro.filt = sorted.dendro#dendextend::prune(sorted.dendro, c("Zebrafish", 'Lamprey', 'Goldfish'))
# dendro = ggdendrogram(sorted.dendro.filt)+
# # scale_y_reverse()
# coord_flip()
# theme(axis.text.y = element_blank(), axis.text.x = element_blank())
dendro
library(ape)
phylo = ape::as.phylo(sorted.dendro.filt)
library(ggtree)
ggtree(phylo) +
geom_tiplab() +
coord_cartesian(clip = 'off')+
# theme_tree2() +
theme(plot.margin = unit(c(1, 3, 1, 1), 'cm')) # xlim(0, max(nodeHeights(phylo)) + 1)
ggsave('figures/phylogeny-22-species.pdf', height=5, width=4.5)
```
## Manhattan marker plot
```{r, fig.width=20}
BIGSEURAT = smartReadRDS("../../Ortho_Objects/vertebrateAC_BIGSEURAT.rds")
gene_type = fread('../../reference_files/gene-type.csv')
markers = gene_type$gene
scores.matrix = readRDS("../../Ortho_Objects/scores.matrix.rds")
scores.matrix.shuffled = readRDS("../../Ortho_Objects/scores.matrix.shuffled.rds")
conservation.summary = readRDS("../../Ortho_Objects/vertebrateAC_conservation.rds")
# Normalization
scores.matrix.melt = reshape2::melt(scores.matrix) %>% setNames(c("gene", "OT", "score"))
scores.matrix.melt$rand.score = reshape2::melt(scores.matrix.shuffled)[,3]
scores.matrix.melt$normalized.score = (scores.matrix.melt$score - mean(scores.matrix.melt$rand.score))/(0 - mean(scores.matrix.melt$rand.score))
scores.matrix.melt$normalized.rand.score = (scores.matrix.melt$rand.score - mean(scores.matrix.melt$rand.score))/(0 - mean(scores.matrix.melt$rand.score))
conservation.summary$best.normalized.score = apply(reshape2::dcast(scores.matrix.melt, gene ~ OT, value.var = 'normalized.score') %>% select(-gene), 1, function(x) max(x))
saveRDS(conservation.summary, "../../Ortho_Objects/vertebrateAC_conservation.rds")
```
```{r, fig.width=20}
BIGSEURAT = LoadACOrtho()
gene_type = fread('../../reference_files/gene-type.csv')
markers = gene_type$gene
scores.matrix = readRDS("../../Ortho_Objects/scores.matrix.rds")
scores.matrix.shuffled = readRDS("../../Ortho_Objects/scores.matrix.shuffled.rds")
conservation.summary = readRDS("../../Ortho_Objects/vertebrateAC_conservation.rds")
# Score normalization
scores.matrix.melt = reshape2::melt(scores.matrix) %>% setNames(c("gene", "OT", "score"))
scores.matrix.melt$rand.score = reshape2::melt(scores.matrix.shuffled)[,3]
scores.matrix.melt$normalized.score = (scores.matrix.melt$score - mean(scores.matrix.melt$rand.score))/(0 - mean(scores.matrix.melt$rand.score))
scores.matrix.melt$normalized.rand.score = (scores.matrix.melt$rand.score - mean(scores.matrix.melt$rand.score))/(0 - mean(scores.matrix.melt$rand.score))
# Threshold selection
specificity.threshold.low = min(scores.matrix.melt$normalized.rand.score)
specificity.threshold.high = max(scores.matrix.melt$normalized.rand.score)
# Marker selection
scores.matrix.melt$marker = scores.matrix.melt$gene %in% markers
scores.matrix.melt$significant = scores.matrix.melt$normalized.score > specificity.threshold.high
scores.matrix.melt$label_all = paste0("italic('", as.character(scores.matrix.melt$gene), "')")
scores.matrix.melt$label_marker = ifelse(scores.matrix.melt$marker & scores.matrix.melt$significant, paste0("italic('", as.character(scores.matrix.melt$gene), "')"), NA)
scores.matrix.melt$label_signif = ifelse(scores.matrix.melt$significant, paste0("italic('", as.character(scores.matrix.melt$gene), "')"), NA)
scores.matrix.melt$label_high = ifelse(scores.matrix.melt$normalized.score > 0.5, paste0("italic('", as.character(scores.matrix.melt$gene), "')"), NA)
# Top 2
scores.matrix.melt %>%
mutate(RowIndex = row_number()) %>%
group_by(OT) %>%
arrange(-normalized.score) %>%
slice_head(n = 2) %>%
ungroup() %>%
pull(RowIndex) -> top2_indices
scores.matrix.melt$label_top2 = NA
scores.matrix.melt$label_top2[top2_indices] = scores.matrix.melt$label_all[top2_indices]
scores.matrix.melt$top2 = FALSE
scores.matrix.melt$top2[top2_indices] = TRUE
scores.matrix.melt$show.label = scores.matrix.melt$top2 | scores.matrix.melt$significant & scores.matrix.melt$marker
# Final labels
scores.matrix.melt$label = ifelse(scores.matrix.melt$show.label, paste0("italic('", as.character(scores.matrix.melt$gene), "')"), NA)
# Odd and even for manhattan plot
scores.matrix.melt$odd = as.logical(scores.matrix.melt$OT %% 2)
# Colors
scores.matrix.melt$color = ifelse(scores.matrix.melt$significant & scores.matrix.melt$marker,
'signif_marker',
ifelse(scores.matrix.melt$significant,
'signif', 'other'))
# scores.matrix.melt$fill = ifelse(scores.matrix.melt$top2, 'signif', 'other'))
ots.with.markers = subset(scores.matrix.melt, marker & normalized.score > specificity.threshold.high)$OT %>% unique
message('Number of OTs with specific and conserved known markers: ', length(ots.with.markers))
ots.with.anymarkers = subset(scores.matrix.melt, normalized.score > specificity.threshold.high)$OT %>% unique
message('Number of OTs with specific and conserved any markers: ', length(ots.with.anymarkers))
scores.matrix.melt$ot.with.markers = scores.matrix.melt$OT %in% ots.with.markers
scores.matrix.melt$ot.with.anymarkers = scores.matrix.melt$OT %in% ots.with.anymarkers
xmin = 1
xmax = nClusters(BIGSEURAT)
ymin = min(scores.matrix.melt$normalized.score)
ymax = max(scores.matrix.melt$normalized.score)
scores.matrix.melt$OT = factor(scores.matrix.melt$OT, levels = levels(BIGSEURAT$dendro.order.type))
scores.matrix.melt = scores.matrix.melt %>% arrange(show.label)
# scores.matrix.melt$type = Metadata(BIGSEURAT, 'type_no', 'type')$type[match(scores.matrix.melt$OT, Metadata(BIGSEURAT, 'type_no', 'type')$type_no)]
scores.matrix.melt$orthotype = Metadata(BIGSEURAT, 'type_no', 'orthotype')$orthotype[match(scores.matrix.melt$OT, Metadata(BIGSEURAT, 'type_no', 'orthotype')$type_no)]
obs = scores.matrix.melt$normalized.score
null = scores.matrix.melt$normalized.rand.score # length 718746 (42 oAC x 17k genes)
# PrettyHistogram2(scores.matrix.melt$normalized.score,
# scores.matrix.melt$normalized.rand.score,
# bins = 100, logY = F,
# vline = max(null)) +
# scale_x_continuous(limits = c(0.1,1))
# coord_cartesian(xlim = c(0.1,1), expand = 0)
# qqplot(
# quantile(null, probs = ppoints(length(obs))),
# sort(obs),
# xlab = "Null distribution quantiles",
# ylab = "Observed quantiles"
# )
# abline(0, 1, col = "red")
p.val <- (length(null) - findInterval(obs, sort(null)) + 1) /
(length(null) + 1)
# Adjustment
p.val.adj = p.adjust(p.val, method = 'fdr')
head(sort(unique(p.val.adj))) # Current cutoff translates to an adjusted p-value of 0.002456999
scores.matrix.melt$p.val = p.val
scores.matrix.melt$p.adj = p.val.adj
# Save top 50 markers for each oAC
# top50_markers = TopN(scores.matrix.melt, 'orthotype', 'normalized.score', 50)
# library(openxlsx)
# wb <- loadWorkbook("../../manuscript/AC/supplement/TableS3-orthotype_markers.xlsx")
# addWorksheet(wb, "Conserved markers 2") # give your new sheet a name
# writeData(wb, sheet = "Conserved markers 2", top50_markers[,c('gene', 'orthotype', 'normalized.score', 'p.val', 'p.adj', 'significant')])
# saveWorkbook(wb, "../../manuscript/AC/supplement/TableS3-orthotype_markers.xlsx", overwrite = TRUE)
```
```{r, fig.width=11, fig.height=4, dev='cairo_pdf''}
p2 = ggplot(scores.matrix.melt,
aes(x = orthotype, y = normalized.score, label = label, color = show.label, fill = color, alpha = show.label))+
annotate("rect",
xmin = seq(xmin-0.5, xmax-0.5, by = 1),
xmax = seq(xmin+0.5, xmax+0.5, by = 1),
ymin = 0, #min(scores.matrix.melt$normalized.score),
ymax = 1, #max(scores.matrix.melt$normalized.score),
alpha = .2, fill = rep(c('white', 'lightgrey'), nClusters(BIGSEURAT)/2))+
# ggtitle('Orthotype markers Known, specific markers Other specific markers Non-specific markers')+
# scale_shape_manual(values = c(21, 21, 21))+
# geom_bin2d()+
# scale_fill_continuous(type = "viridis") +
xlab('Orthotype')+
ylab('Specificity score')+
scale_y_continuous(expand = expansion(mult = c(0, 0)), limits = c(0, 1))+
scale_x_discrete(expand = expansion(mult = c(0, 0)))+ #breaks = seq(1,nClusters(BIGSEURAT),1),
ggrastr::rasterise(geom_jitter(data = subset(scores.matrix.melt, !show.label), position = position_jitter(seed = 1), shape = 21), dpi = 300)+
geom_jitter(data = subset(scores.matrix.melt, show.label), position = position_jitter(seed = 1), shape = 21)+
# geom_jitter(data = subset(scores.matrix.melt, !show.label), position = position_jitter(seed = 1), shape = 21)+
geom_hline(yintercept = c(specificity.threshold.high), linetype = 'dashed', color = fcols[2], linewidth = 0.25)+
# geom_jitter(data = subset(scores.matrix.melt, !significant), position = position_jitter(seed = 1), shape = 21, alpha = 0.1)+
# geom_jitter(data = subset(scores.matrix.melt, significant), position = position_jitter(seed = 1), shape = 21)+
scale_color_manual('', values = c('grey', 'black'), labels = c('Novel', 'Known'))+
scale_fill_manual('', values = c('grey', 'grey', fcols[2]), labels = c('Novel', 'Novel', 'Known'))+
scale_alpha_manual(values = c(1, 1))+
guides(color = 'none', alpha = 'none', label = 'none')+
geom_text_repel(data = subset(scores.matrix.melt, show.label), max.overlaps = Inf, min.segment.length = 0.05, position = position_jitter(seed = 1), size = 2, parse = TRUE)+
theme_dario() +
theme(axis.title.x = element_blank())+
RotatedAxis() +
NoLegend() +
ArialFont()
# theme(legend.position = 'right')
p2
# ggsave('../../figures/my_figs/specificity-scores.pdf', width = 11, height = 4, device = cairo_pdf)
ggsave('../../figures/my_figs/specificity-scores.pdf', width = 11, height = 4.2)
```
Only displaced ACs
```{r}
p3 = ggplot(scores.matrix.melt[grepl('\\*', scores.matrix.melt$type),],
aes(x = type, y = normalized.score, label = label, color = show.label, fill = color, alpha = show.label))+
annotate("rect",
xmin = seq(1-0.5, 11-0.5, by = 1),
xmax = seq(1+0.5, 11+0.5, by = 1),
ymin = 0, #min(scores.matrix.melt$normalized.score),
ymax = 1, #max(scores.matrix.melt$normalized.score),
alpha = .2, fill = rep(c('white', 'lightgrey'), 12/2)[1:11])+
# ggtitle('Orthotype markers Known, specific markers Other specific markers Non-specific markers')+
# scale_shape_manual(values = c(21, 21, 21))+
# geom_bin2d()+
# scale_fill_continuous(type = "viridis") +
xlab('Orthotype')+
ylab('Specificity score')+
scale_y_continuous(expand = expansion(mult = c(0, 0)), limits = c(0, 1))+
scale_x_discrete(expand = expansion(mult = c(0, 0)))+ #breaks = seq(1,nClusters(BIGSEURAT),1),
ggrastr::rasterise(geom_jitter(data = subset(scores.matrix.melt[grepl('\\*', scores.matrix.melt$type),], !show.label),
position = position_jitter(seed = 1), shape = 21), dpi = 300)+
geom_jitter(data = subset(scores.matrix.melt[grepl('\\*', scores.matrix.melt$type),], show.label),
position = position_jitter(seed = 1), shape = 21)+
# geom_jitter(data = subset(scores.matrix.melt, !show.label), position = position_jitter(seed = 1), shape = 21)+
geom_hline(yintercept = c(specificity.threshold.high), linetype = 'dashed', color = fcols[2], linewidth = 0.25)+
# geom_jitter(data = subset(scores.matrix.melt, !significant), position = position_jitter(seed = 1), shape = 21, alpha = 0.1)+
# geom_jitter(data = subset(scores.matrix.melt, significant), position = position_jitter(seed = 1), shape = 21)+
scale_color_manual('', values = c('grey', 'black'), labels = c('Novel', 'Known'))+
scale_fill_manual('', values = c('grey', 'grey', fcols[2]), labels = c('Novel', 'Novel', 'Known'))+
scale_alpha_manual(values = c(1, 1))+
guides(color = 'none', alpha = 'none', label = 'none')+
geom_text_repel(data = subset(scores.matrix.melt[grepl('\\*', scores.matrix.melt$type),], show.label),
max.overlaps = Inf, min.segment.length = 0.05, position = position_jitter(seed = 1), size = 2, parse = TRUE)+
theme_dario() +
RotatedAxis() +
NoLegend() +
ArialFont()
# theme(legend.position = 'right')
p3
# ggsave('../../figures/my_figs/specificity-scores.pdf', width = 11, height = 4, device = cairo_pdf)
ggsave('../../figures/my_figs/specificity-scores-dACs.pdf', width = 4, height = 4)
```
## Conserved markers example
```{r, fig.width=4, fig.height=2.5}
# coord_flip_discrete(geneDotPlotFast(ac.ortho, gene = 'PDGFRA', mini = TRUE), discrete_axis = 'OT')
# MultigeneDotPlot(ac.ortho, c('PDGFRA', 'EOMES'), mini = TRUE)
geneDotPlotFast(ac.ortho, gene = 'VIP', mini = TRUE) +
coord_cartesian(clip = 'off') +
theme_void() + theme_dario(linewidth = 0.6) + NoLegend() + theme(plot.title = element_blank())
ggsave('../../figures/my_figs/fig3-vip-example.pdf', width = 5.5, height = 3)
```
# Figure 3: histological validation
```{r, fig.height=3, fig.width=2}
TwoGeneDotPlot(ac.ortho, '27_VG3', c('SLC17A8', 'NXPH1'), font.size = 12, legend = TRUE) +
scale_color_manual(values = c('magenta3', 'cyan3'), guide = "none") +
theme(plot.margin = unit(c(0,0.5,0,0.5), 'cm')) + ArialFont()
ggsave('../../figures/my_figs/vglut3-nxph1-dot-legend.pdf', height=3, width=3)
width = 1.9
# VGluT3
TwoGeneDotPlot(ac.ortho, '27_VG3', c('SLC17A8', 'NXPH1'), font.size = 12) +
scale_color_manual(values = c('magenta3', 'cyan3'), guide = "none") +
theme(plot.margin = unit(c(0,0.5,0,0.5), 'cm'))
ggsave('../../figures/my_figs/vglut3-nxph1-dot.pdf', height=3, width=width)
# PDGFRA
TwoGeneDotPlot(ac.ortho, '21_PDGFRA', c('PDGFRA', 'GAD1')) +
scale_color_manual(values = c('magenta3', 'cyan3'), guide = "none") +
theme(plot.margin = unit(c(0,0.5,0,0.5), 'cm'))
ggsave('../../figures/my_figs/pdgfra-gad1-dot.pdf', height=3, width=width)
# CART SLC35D3
TwoGeneDotPlot(ac.ortho, '29_SLC35D3', c('SLC35D3', 'CARTPT')) +
scale_color_manual(values = c('magenta3', 'cyan3'), guide = "none") +
theme(plot.margin = unit(c(0,0.5,0,0.5), 'cm'))
ggsave('../../figures/my_figs/slc35d3-cartpt-dot.pdf', height=3, width=width)
# MAF GAD
TwoGeneDotPlot(ac.ortho, '1_MAF*', c('MAF', 'GAD1')) +
scale_color_manual(values = c('magenta3', 'cyan3'), guide = "none") +
theme(plot.margin = unit(c(0,0.5,0,0.5), 'cm'))
ggsave('../../figures/my_figs/maf1-gad-dot.pdf', height=3, width=width)
TwoGeneDotPlot(ac.ortho, '10_MAF*', c('MAF', 'GAD1')) +
scale_color_manual(values = c('magenta3', 'cyan3'), guide = "none") +
theme(plot.margin = unit(c(0,0.5,0,0.5), 'cm'))
ggsave('../../figures/my_figs/maf10-gad-dot.pdf', height=3, width=width)
# NOS/EOMES
TwoGeneDotPlot(ac.ortho, '12_nNOS*', c('NOS1', 'EOMES')) +
scale_color_manual(values = c('magenta3', 'cyan3'), guide = "none") +
theme(plot.margin = unit(c(0,0.5,0,0.5), 'cm'))
ggsave('../../figures/my_figs/nos-eomes-dot.pdf', height=3, width=width)
```
## MAF proportion displaced quantification
```{r, fig.height=3.5, fig.width=2.6}
YanAC = readRDS('../../Species_Reference/YanAC_v4.rds')
table(YanAC$cluster_no)
prop.displaced.merfish = readRDS('../../reference_files/Choi_fig4d_automeris.rds')
prop.displaced.merfish
# 1_MAF is 2, 10_MAF is 32
ac2.prop = prop.displaced.merfish$proportion[prop.displaced.merfish$type == 'AC-2']
ac32.prop = prop.displaced.merfish$proportion[prop.displaced.merfish$type == 'AC-32']
weighted.average = (ac2.prop*table(YanAC$cluster_no)[[2]] + ac32.prop*table(YanAC$cluster_no)[[32]])/
(table(YanAC$cluster_no)[[2]] + table(YanAC$cluster_no)[[32]])
# Mouse replicates: C538A, C580, C588
# Macaque replicates: C285, C271, C292
prop.displaced.df = data.frame(Method = c('Mouse MERFISH', 'Mouse IHC 1', 'Mouse IHC 2 ', 'Mouse IHC 3', 'Macaque IHC 1', 'Macaque IHC 2', 'Macaque IHC 3'),
Sample = c('Mouse MERFISH', 'C538', 'C580', 'C588', 'C285', 'C271', 'C292'),
GCL = c(304, 13, 10, 18, 4, 7, 16),
INL = c(1000-304, 37-13, (23), (34), 12-4, (12), (27)),
proportion.displaced = c(weighted.average, 13/37, 10/(23+10), 18/(34+18), 4/12, 7/(12+7), 16/(27+16))) # 196 MAF+ cells counted
# prop.displaced.df$Method = factor(prop.displaced.df$Method, levels = rev(c('Mouse MERFISH', 'Mouse IHC', 'Macaque IHC')))
# prop.displaced.df$prop.displaced = prop.displaced.df$n.cells.GCL/(prop.displaced.df$n.cells.GCL + prop.displaced.df$n.cells.INL)
# TitlePlot(
# PrettyBarplot(prop.displaced.df, x = 'Method', y = 'proportion.displaced', add = c('mean_sd'), fill = 'white', width = 0.6) +
# geom_jitter(color = 'darkgrey')+
# ylab('Proportion\nin GCL') +
# xlab(NULL)+
# # RotatedAxis()+
# scale_x_discrete(labels = function(x) gsub(' ', '\n', x))+
# RotatedXAxis(angle = 0, hjust = 0.5)+
# # scale_y_continuous(labels = c(0,'', '', '', 0.4))
# # scale_y_continuous(breaks = function(x) range(x))+
# # scale_y_continuous(
# # breaks = c(0,0.4)
# # # breaks = function(x) {
# # # b <- scales::pretty_breaks()(x)
# # # c(min(b), max(b))
# # # }
# # )
# scale_ticks_first_last()+
# coord_flip()
# # theme(plot.margin = unit(c(0,0,0,0.5), 'cm'), axis.title.x = element_text(size = 16), axis.text.y = element_text(size = 14))
# # ylab(NULL)+
# # theme(axis.text.x = element_text(hjust = 0.5, angle = 0, size = 9))
# # 'Proportion\ndisplaced')
temp = reshape2::melt(prop.displaced.df[,c('Method', 'INL', 'GCL')])
ggplot(temp, aes(fill=variable, y=value, x=Method)) +
geom_bar(position="fill", stat="identity", color = "black") +
scale_x_discrete(labels = function(x) gsub(' ', '\n', gsub(" [0-9]+", "", x)))+
theme_minimal()+
# ggtitle(title)+
# xlab(feature.2)+
# labs(fill=feature.1)+
theme_cowplot(font_family = 'ArialMT')+
scale_fill_manual(name = '', values = c(INL = 'cyan3', GCL = 'deeppink2'))+
coord_flip()+
scale_ticks_first_last()+
# scale_x_continuous(labels = function(x) gsub('0\\.', '', x))+
ylab('Proportion')+
xlab(NULL)+
theme(legend.position = 'top')
# theme(axis.text.x = element_text(hjust = 1, angle = 45), plot.title = element_text(hjust = 0.5))
# ggsave('../../figures/my_figs/fig-hist-maf-prop.pdf', height = 3.5, width = 2.6)
```
## oAC30 [PDGFRA] abundance
```{r, fig.width = 4.5, fig.height = 2.6}
freq.data2 = readRDS('../../Ortho_Objects/freq.data2.rds')
# PDGFRA is enriched in higher primates
PrettyBarplot(subset(freq.data2, Type == '21_PDGFRA' & jaccard > 0.25), x = 'Species', y = 'Frequency', add = c('mean_sd'), fill = 'white') +
geom_jitter(aes(color = Region), width = 0.2)+
scale_color_manual(values = c(Periphery = 'cyan3', Fovea = 'deeppink2'), na.value = "darkgrey")+
ylab('oAC30 [PDGFRA] frequency')+
xlab(NULL)+
theme(axis.title.y = element_text(hjust = 1, size = 13))+
NoLegend()+
RotatedAxis()
ggsave('../../figures/my_figs/fig-pdgfra-freq.pdf', width = 4.5, height = 2.6)
```
# Figure 4: conservation of orthotypes
```{r, fig.height=9, fig.width=10}
ac.ortho = LoadACOrtho()
js.matrices = lapply(speciesList.liz, function(this.species){
object = subset(ac.ortho, species == this.species)
JSMatrix(table(object$type_unordered, object$species_type2))
# addmargins(table(object$type, object$species_type2))
})
best.jaccard = do.call(cbind, lapply(js.matrices, function(matrix){
apply(matrix, 1, max)
}))
saveRDS(best.jaccard, '../../Ortho_Objects/best.jaccard.matrix.rds')
```
## Panel A: jaccard matrices
```{r orthotype-speciestype-tight, fig.height=3, fig.width=20}
nTypes.df = readRDS('metadata/nTypes.df.rds')
tableList = readRDS('../../Ortho_Objects/tableList.rds')
ncells = nTypes.df$nCells[match(c(speciesList, 'Zebrafish', 'Lamprey'), nTypes.df$species)] %>% setNames(c(speciesList, 'Zebrafish', 'Lamprey'))
heatmapList2 = lapply(seq_along(tableList), function(index) {
JSHeatmap(JSMatrix(tableList[[index]]),
heatmap = TRUE,
border.col = NA,
title = paste0(unique(BIGSEURAT$species)[[index]], '\n(', ncells[[index]], ')'),
row.order = levels(ac.ortho$type),
stagger.threshold = 0.20) + theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.title = element_blank(),
plot.margin=grid::unit(c(0,0,0,0), "mm"),
axis.ticks = element_blank()) + NoLegend()
})
nclusters = sapply(speciesACList, function(object) length(unique(object$species_type)))
ggarrange(plotlist = heatmapList2,
ncol = (length(heatmapList2)),
nrow = 1,
widths = nclusters)
```
```{r orthotype-speciestype-tight, fig.height=4, fig.width=20}
nTypes.df = readRDS('metadata/nTypes.df.rds')
tableList = readRDS('../../Ortho_Objects/tableList.rds')
ncells = nTypes.df$nCells[match(c(speciesList, 'Zebrafish', 'Lamprey'), nTypes.df$species)] %>% setNames(c(speciesList, 'Zebrafish', 'Lamprey'))
heatmapList2 = lapply(seq_along(tableList), function(index) {
JSHeatmap(JSMatrix(tableList[[index]]),
heatmap = TRUE,
border.col = NA,
title = paste0(unique(BIGSEURAT$species)[[index]], '\n(', comma(ncells[[index]]), ')'),
row.order = levels(ac.ortho$type),
stagger.threshold = 0.20) + theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.title = element_blank(),
plot.margin=grid::unit(c(0,0,0,0), "mm"),
axis.ticks = element_blank()) + NoLegend() + ArialFont()
})
nclusters = sapply(speciesACList, function(object) length(unique(object$species_type)))
ggarrange(plotlist = heatmapList2,
ncol = (length(heatmapList2)),
nrow = 1,
widths = nclusters)
ggsave('../../figures/my_figs/orthotype-speciestype-supertight.pdf')
```
Split by clades
```{r, fig.height=7, fig.width=14}
# plot_grid(
# ggarrange(plotlist = heatmapList2[1:5],
# ncol = 5,
# nrow = 1,
# widths = nclusters[1:5]),
# ggarrange(plotlist = heatmapList2[6:10],
# ncol = 5,
# nrow = 1,
# widths = nclusters[6:10]),
# ggarrange(plotlist = heatmapList2[11:14],
# ncol = 5,
# nrow = 1,
# widths = nclusters[11:14]),
# ggarrange(plotlist = heatmapList2[15],
# ncol = 5,
# nrow = 1,
# widths = nclusters[[15]]),
# ggarrange(plotlist = heatmapList2[16:17],
# ncol = 5,
# nrow = 1,
# widths = nclusters[16:17]),
# ncol = 1, nrow = 2, labels = LETTERS, label_size = 20
# )
plot_grid(
ggarrange(plotlist = heatmapList2[1:8],
ncol = 8,
nrow = 1,
widths = nclusters[1:8]),
NULL,
ggarrange(plotlist = heatmapList2[9:17],
ncol = 9,
nrow = 1,
widths = nclusters[9:17]),
ncol = 1, nrow = 3, rel_heights = c(1,0.1,1) #, labels = LETTERS, label_size = 20
) + theme(plot.margin = unit(c(1,1,1,1), units = 'cm'))
ggsave('../../figures/my_figs/orthotype-speciestype-two-row.pdf')
```
## Cluster merges plot
```{r}
bestMatchDf = readRDS('../../Ortho_Objects/best-match.rds')[,speciesList.liz]
bestJaccardDf = readRDS('../../Ortho_Objects/best-jaccard.rds')[,speciesList.liz]
summary = reshape2::melt(as.matrix(bestMatchDf[,speciesList.liz])) %>% setNames(c('type', 'species', 'label'))
summary$score = reshape2::melt(as.matrix(bestJaccardDf[,speciesList.liz]))$value
# For each label in each species, find its best type (that will be its color!)
summary$species_type = paste0(summary$species, '-', summary$label)
best.ot = dapply(unique(summary$species_type), function(x) {
this = subset(summary, species_type == x)
as.character(this$type[which.max(this$score)])
})
summary$best.ot = best.ot[summary$species_type]
cor.mat = reshape2::dcast(summary, species ~ type, value.var = "best.ot")[,-1]
cor.mat <- apply(cor.mat, 2, as.character)
jac.mat = t(bestJaccardDf)
cor.mat[jac.mat < JS.THRESHOLD] = NA
rownames(cor.mat) = speciesList.liz
colnames(cor.mat) = orthotype_labels(colnames(cor.mat))
# pdf('../../figures/my_figs/merged-clusters-plot.pdf', height = 4, width = 12)
# Heatmap(as.matrix(cor.mat),
# cell_fun = function(j, i, x, y, width, height, fill) {
# grid.text(round(jac.mat[i, j], 2), x, y, gp = gpar(fontsize = 6))
# },
# col = type_cols2,
# cluster_rows = FALSE,
# cluster_columns = FALSE)
# dev.off()
```
```{r, fig.height=5.5, fig.width=11}
# Make colors a bit more distinct
type_cols4
# Quantification
ac.enriched = readRDS('../../Ortho_Objects/AC-frequency-all.rds')
ac.enriched.df = reshape2::melt(ac.enriched) %>% setNames(c('Species', 'Orthotype', 'Frequency'))
ac.enriched.df$Orthotype2 = orthotype_labels(gsub('\\*', '', ac.enriched.df$Orthotype), asterisk = F)
plt1 = PrettyBoxplot2(ac.enriched.df, x = 'Orthotype2', y = 'Frequency', color = 'Orthotype2', jitter = F) +
ylab('Observed\nfrequency')+
scale_color_manual(values = darken(type_cols4, 0.1))+
theme(axis.line.x = element_blank(), axis.ticks.x = element_blank(), axis.text.x = element_blank(),
plot.margin = unit(c(0.5,0,0,0), 'cm')) +
scale_y_continuous(expand = expansion(mult = c(0,0)))+
coord_cartesian(ylim = c(0, 0.1501))
# Snake
rownames(cor.mat) = speciesList.liz
cor.mat2 = matrix(orthotype_labels(cor.mat), nrow = nrow(cor.mat), ncol = ncol(cor.mat), dimnames = dimnames(cor.mat))
plt2 = SnakePlot(cor.mat2,
pt.size = 4,
lwd1 = 5,
lwd2 = NA,
lty2 = 'solid',
grid.pt.size = 0.2,
cols = lighten(type_cols4, 0.1))
plot_grid(plt1, plt2, nrow = 2, align = 'v', axis = 'lr', rel_heights = c(1,3.5))
ggsave('../../figures/my_figs/fig3-snake-plot.pdf', height=5.5, width=11)
```
```{r, fig.height=3.8, fig.width=10}
cor.mat.mammal = cor.mat[rownames(cor.mat) %in% speciesList.mammals,]
SnakePlot(cor.mat.mammal, pt.size = 4, lwd1 = 5, grid.pt.size = 0.2)
ggsave('../../figures/my_figs/fig3-snake-plot-mammal.pdf')
```
Investigate human NN-gly cluster to see if they can be subclustered
```{r}
human = readRDS('../../Species_Objects/HumanAC_v6.rds')
BIGSEURAT = readRDS('../../Ortho_Objects/vertebrateAC-allcells-counts.rds')
human.10 = subset(BIGSEURAT, species_type == 'Human-10_SEG')
human.10 = MapBack(human.10, human, group.by = 'donor')
human.10.harmony = Harmonize(human.10, batch = 'donor')
DimPlotLabeled(human.10.harmony, 'type')
DimPlotLabeled(human.10.harmony, 'donor')
human.10.int = ClusterSeurat(human.10, integrate.by = 'donor')
DimPlotLabeled(human.10.int, 'type')
DimPlotLabeled(human.10.int, 'donor')
stackedBarGraph2(human.10.int, 'type', 'donor')
FeaturePlot(human.10.int, 'PROM1')
FeaturePlot(human.10.int, 'SATB2')
FeaturePlot(human.10.int, 'NFIB')
# Challenging to split these three types
# Human and macaque integration
hum_mac = subset(BIGSEURAT, species %in% c('Human', 'Macaque'))
# hum_mac = ClusterSeurat(hum_mac, integrate.by = 'species')
```
Also interested in oAC splits, not just merges
```{r}
jaccardList = readRDS('../../Ortho_Objects/jaccardList.rds')
jaccardList
# How many species types map to this orthotype above some threshold?
type.depth = lapply(jaccardList, function(this.species) apply(this.species, 1, function(row) length(which(row > js.threshold))))
type.depth.mat = do.call(rbind, type.depth)
Heatmap2(type.depth.mat)
```
Co-clustering frequency
```{r}
tableList = readRDS('../../Ortho_Objects/tableList.rds')
# JSHeatmap(t(jaccardList$Human) %*% jaccardList$Mouse)
tableList = lapply(names(tableList), function(species) {
colnames(tableList[[species]]) = paste0(species, '-', colnames(tableList[[species]]))
tableList[[species]]
})
names(tableList) = speciesList.liz
ccf = row_norm(t(tableList$Human)) %*% row_norm(tableList$Rat)
JSHeatmap(ccf, stagger.threshold = 0.2)
ccf = row_norm(t(tableList$Chicken)) %*% row_norm(tableList$Rat)
JSHeatmap(ccf, stagger.threshold = 0.2)
```
```{r, fig.height=6, fig.width=15}
library(ggplot2)
library(dplyr)
library(tidyr)
# Example character matrix (3 rows x 4 cols)
# mat <- matrix(
# c("A", "A", "B", "C",
# "B", "B", "B", "C",
# "C", "C", "A", "A"),
# nrow = 3, byrow = TRUE,
# dimnames = list(paste0("row", 1:3), paste0("col", 1:4))
# )
mat = cor.mat
rownames(mat) = speciesList.liz
# Convert to long
df <- as.data.frame(mat) %>%
tibble::rownames_to_column("row") %>%
pivot_longer(-row, names_to = "col", values_to = "value")
# Create numeric coordinates once (consistent for points & segments)
# Keep the column order as in mat and row order as in mat
col_levels <- colnames(mat)
row_levels <- rownames(mat)
df <- df %>%
mutate(
col_num = as.integer(factor(col, levels = col_levels)),
row_num = as.integer(factor(row, levels = row_levels))
)
# Build segments: connect adjacent columns within the same row when values match
segments <- df %>%
arrange(row_num, col_num) %>%
group_by(row) %>%
mutate(next_value = lead(value), next_col_num = lead(col_num)) %>%
filter(value == next_value) %>%
transmute(
row, value,
x = col_num,
xend = next_col_num,
y = row_num,
yend = row_num
) %>%
ungroup()
# Plot: use numeric coords, convert axes to show labels
ggplot() +
# segments first so points draw on top
geom_segment(data = segments,
aes(x = x, xend = xend, y = y, yend = yend, color = value),
linewidth = 1.2, lineend = "round", show.legend = FALSE) +
geom_point(data = df,
aes(x = col_num, y = row_num, color = value),
size = 5, show.legend = FALSE) +
scale_x_continuous(breaks = seq_along(col_levels), labels = col_levels, expand = expansion(0.1)) +
# reverse y so row1 is on top (optional, UpSet-like)
scale_y_reverse(breaks = seq_along(row_levels), labels = row_levels, expand = expansion(0.1)) +
theme_minimal() +
RotatedAxis()+
theme(panel.grid = element_blank(),
axis.title = element_blank())
```
## Conservation level quantification and comparison
```{r, fig.height=3, fig.width=3.2}
# Quantification of AC conservation in mammals; we have no way of disentangling presence and conservation
# bestMatchDf = readRDS('../../Ortho_Objects/best-match.rds')[,speciesList.liz]
# bestJaccardDf = readRDS('../../Ortho_Objects/best-jaccard.rds')[,speciesList.liz]
# jac.mat = t(bestJaccardDf[,speciesList.mammals])
# Area under the presence cutoff curve
AUPCC = function(jac.mat, return.values = FALSE){
library(pracma)
x = seq(0, 1, length.out = 101)
y = sapply(x, function(i) length(which(jac.mat > i)) / prod(dim(jac.mat)))
plot(x, y)
area <- trapz(x, y)
auc_norm <- area / (max(x) - min(x))
if(return.values) {
y
} else {
auc_norm
}
}
# With ferret
ac.ortho = LoadACOrtho()
ac.match = MakeBestMatchDF(ac.ortho)
# AUPCC(t(ac.match[[2]]))
# Removing some species to make sure they have the same set
species.use = c('Human', 'Macaque', 'Marmoset', 'TreeShrew',
'Mouse', 'Rhabdomys', 'Peromyscus', 'Squirrel',
'Sheep', 'Pig','Opossum')
# Remove some lowly abundant clusters
ac.match = MakeBestMatchDF(ac.ortho)
ac.match2 = ac.match[[2]]
# ac.match2 = ac.match2[,!colnames(ac.match2) %in% c('Ferret', 'Chicken', 'Lizard')]
ac.match2 = ac.match2[!rownames(ac.match2) %in% c('38', '39', '40*', '42'), colnames(ac.match2) %in% species.use]
AUPCC(t(ac.match2))
# BCs
bc.ortho = readRDS('../../Ortho_Objects/vertebrateBC_BIGSEURAT.rds')
bc.match = MakeBestMatchDF(bc.ortho, species_cluster = 'species_type', orthotype = 'NOG')
AUPCC(t(bc.match[[2]] %>% dplyr::select(species.use)))
# RGCs
rgc.ortho = readRDS('../../Ortho_Objects/vertebrateRGC_BIGSEURAT.rds')
rgc.match = MakeBestMatchDF(rgc.ortho, species_cluster = 'species_type', orthotype = 'NOG')
AUPCC(t(rgc.match[[2]] %>% dplyr::select(species.use)))
# Collect the values
auc.df = data.frame(cutoff = rep(seq(0, 1, length.out = 101), 3),
AUC = factor_unique(c(rep(paste0(c('BC', round(AUPCC(t(bc.match[[2]] %>% dplyr::select(species.use))), 2)), collapse = ' = '), 101),
rep(paste0(c('AC', round(AUPCC(t(ac.match[[2]] %>% dplyr::select(species.use))), 2)), collapse = ' = '), 101),
rep(paste0(c('RGC', round(AUPCC(t(rgc.match[[2]] %>% dplyr::select(species.use))), 2)), collapse = ' = '), 101))),
fraction = c(AUPCC(t(bc.match[[2]] %>% dplyr::select(species.use)), return.values = T),
AUPCC(t(ac.match[[2]] %>% dplyr::select(species.use)), return.values = T),
AUPCC(t(rgc.match[[2]] %>% dplyr::select(species.use)), return.values = T)))
ggplot(auc.df, aes(x = cutoff, y = fraction * 100, color = AUC)) +
LegendLowerLeft(background.col = 'white', x = 0.01, y = 0.04, background.alpha = 1)+
geom_line() +
theme_dario(linewidth = 1) +
ylab('% orthotypes found')+
xlab('Jaccard index')+
scale_x_continuous(
breaks = seq(0, 1, by = 0.25),
labels = function(b) {
labs <- rep("", length(b))
labs[1] <- b[1]
labs[length(b)] <- b[length(b)]
labs
})+
scale_y_continuous(
breaks = seq(0, 100, by = 25),
labels = function(b) {
labs <- rep("", length(b))
labs[1] <- b[1]
labs[length(b)] <- b[length(b)]
labs
})+
theme(axis.title = element_text(size = 12))+
ArialFont()+
geom_vline(xintercept = 0.1, linetype = 'dashed', color = 'grey')
ggsave('../../figures/my_figs/fig-aupcc.pdf', height=3, width=3.2)
```
## SAC UMAPs
```{r, fig.height=8, fig.width=10}
objectList = list()
objectList$Primates = readRDS('../../Ortho_Objects/primateSAC.rds')
objectList$TreeShrew = subset(readRDS('../../Ortho_Objects/otherSAC.rds'), species == 'TreeShrew')
objectList$Rodents = readRDS('../../Ortho_Objects/rodentSAC.rds')
objectList$Laurasiatherians = readRDS('../../Ortho_Objects/lauraSAC.rds')
# objectList$Other = readRDS('../../Ortho_Objects/otherSAC.rds')
objectList$Opossum = subset(readRDS('../../Ortho_Objects/otherSAC.rds'), species == 'Opossum')
objectList$Sauropsids = readRDS('../../Ortho_Objects/sauroSAC.rds')
objectList$Amphibians = readRDS('../../Ortho_Objects/amph.SAC.rds')
objectList$Zebrafish = readRDS('../../Ortho_Objects/zeSAC.rds')
objectList$Lamprey = readRDS('../../Ortho_Objects/SAC.pm.clean.rds')
pm2ch = readRDS('../cones/orthology_graphs/pm2ch.rds')
objectList$Lamprey = ConvertGeneSymbols2(objectList$Lamprey, pm2ch)
umap_params = list(group.by = 'species', cols = species_palette3, label = FALSE, show.legend = FALSE, pad = 0.4, pt.size = 2, pt.alpha = 0.6)
# angle.list = c(0, pi, pi/1.1, 0, pi, pi, pi/1.5, 0)
angle.list = c(0, pi+0.1, pi, -0.3, pi+0.1, pi, pi, pi/1.3, 0.3)
plt.list = lapply(seq_along(objectList), function(i){
theme_umap(
do.call(PrettyUmap2, c(objectList[[i]], umap_params, angle = angle.list[[i]])) +
theme(plot.background = element_rect(fill = 'transparent')),
remove.axes = T) + theme(plot.margin = unit(c(2,10,10,10), 'mm'))
# LegendLowerLeft(background.col = 'white')
})
plot_grid(plotlist = plt.list)
ggsave('../../figures/my_figs/fig-sac-umap-2.pdf', width = 10, height = 7)
# FeaturePlot(objectList$Amphibians, features = 'TENM3')
# FeaturePlot(objectList$Zebrafish, features = 'tenm3')
# ON SAC is on the left
```
```{r}
umap_params = list(group.by = 'lit_subtype', label = FALSE, show.legend = TRUE, pad = 0.4, pt.size = 2, pt.alpha = 0.6)
angle.list = c(0, pi, 0, pi, pi, 0, 0, 0)
plt.list = lapply(seq_along(objectList), function(i){
do.call(PrettyUmap2, c(objectList[[i]], umap_params, angle = angle.list[[i]])) +
LegendLowerLeft()
})
plot_grid(plotlist = plt.list)
```
```{r, fig.height=4, fig.width=15}
lamprey.full = readRDS('../../Full_Objects/Lamprey_Wang_full_v2.rds')
lamprey.full = ConvertGeneSymbols2(lamprey.full, pm2ch)
DotPlot3(lamprey.full, features = c('NOS1', 'OPN4'), group.by = 'annotated')
FeaturePlot(objectList$Lamprey, 'LOC116944204') | DimPlotLabeled(objectList$Lamprey, 'annotated')
DotPlot3(objectList$Lamprey, features = c('TENM3', 'FEZF1', 'FEZF2', 'LOC116944204'), group.by = 'annotated')
de = TopNDEGs(objectList$Lamprey, group.by = 'annotated')
obj = load('../../Species_Reference/Wang_Lamprey/seu.lamprey.SAC.types.rda')
wang_lamprey <- get(obj)
DotPlot3(wang_lamprey, features = c('TENM3', 'FEZF2'), group.by = 'subtype')
objectList$Lamprey = MapBack(objectList$Lamprey, wang_lamprey, group.by = 'subtype')
table(objectList$Lamprey$annotated, objectList$Lamprey$subtype)
```
## SAC clustering
```{r}
# Need to redo some graph because I removed a cluster
objectList$Amphibians = ClusterSeurat(objectList$Amphibians, integrate.by = 'species')
objectList$Zebrafish = ClusterSeurat(objectList$Zebrafish)
# DimPlotLabeled(objectList$Amphibians)
DimPlotLabeled(objectList$Zebrafish)
objectList$Lamprey = Harmonize(objectList$Lamprey, batch = 'animal')
objectList2 = dapply(names(objectList), function(group){
message('Trying ', group, '...')
object = objectList[[group]]
if(!group %in% c('Zebrafish', 'Lamprey')) DefaultAssay(object) = 'integrated'
ClusterUntil(object, k = 2)
# DimPlotLabeled(objectList$Primates)
})
umap_params = list(group.by = 'seurat_clusters', label = FALSE, show.legend = FALSE, pad = 0.4, pt.size = 2, pt.alpha = 0.6)
angle.list = c(0, pi, pi/1.1, 0, pi, pi, pi/1.5, 0)
plt.list = lapply(seq_along(objectList2), function(i){
theme_umap(
do.call(PrettyUmap2, c(objectList2[[i]], umap_params, angle = angle.list[[i]])) +
theme(plot.background = element_rect(fill = 'transparent')),
remove.axes = T) + theme(plot.margin = unit(c(2,10,10,10), 'mm'))
# LegendLowerLeft(background.col = 'white')
})
plot_grid(plotlist = plt.list)
# Fix the sauropsids
objectList2$Sauropsids = FindClusters(objectList2$Sauropsids, resolution = 0.2)
DimPlotLabeled(objectList2$Sauropsids)
objectList2$Sauropsids = MergeClusters(objectList2$Sauropsids, c(1,2), refactor = TRUE)
DimPlotLabeled(objectList2$Sauropsids)
saveRDS(objectList2, '../../Ortho_Objects/sac-ortho-list.rds')
```
## SAC violin
```{r, fig.height=7, fig.width=6}
objectList2 = readRDS('../../Ortho_Objects/sac-ortho-list.rds')
```
```{r, fig.height=6, fig.width=10}
# Annotate on and off sacs
anno_list = list(Primates = list('SAC' = c(0,1)),
Rodents = list('ON' = c(1), 'OFF' = c(0)),
Laurasiatherians = list('ON' = c(0), 'OFF' = c(1)),
Other = list('ON' = c(1), 'OFF' = c(0)),
Sauropsids = list('ON' = c(1), 'OFF' = c(0)),
Amphibians = list('ON' = c(1), 'OFF' = c(0)),
Zebrafish = list('ON' = c(0), 'OFF' = c(1)),
Lamprey = list('ON' = c(0), 'OFF' = c(1)))
objectList2 = dapply(names(objectList2), function(group){
AssignAnnotations(objectList2[[group]], anno_list[[group]])
})
umap_params = list(group.by = 'lit_type', label = FALSE, show.legend = TRUE, pad = 0.4, pt.size = 2, pt.alpha = 0.6)
angle.list = c(0, pi, -0.3, pi, pi, pi, pi/1.3, 0.3)
plt.list = lapply(seq_along(objectList2), function(i){
theme_umap(
do.call(PrettyUmap2, c(objectList2[[i]], umap_params, angle = angle.list[[i]])) +
theme(plot.background = element_rect(fill = 'transparent')),
remove.axes = T) + theme(plot.margin = unit(c(2,10,10,10), 'mm'))
# LegendLowerLeft(background.col = 'white')
})
plot_grid(plotlist = plt.list)
```
```{r, fig.height=7, fig.width=6}
# Convert some gene symbols
objectList2$Zebrafish = RenameFeatures(objectList2$Zebrafish, gsub('tenm3', 'TENM3',
gsub('zfhx3', 'ZFHX3',
gsub('rnd3a', 'RND3',
gsub('fezf1', 'FEZF1',
gsub('fezf2', 'FEZF2',
rownames(objectList2$Zebrafish)))))))
# Add FEZF2 from Wang et al. (replacing a random unexpressed gene)
obj = load('../../Species_Reference/Wang_Lamprey/seu.lamprey.SAC.types.rda')
wang_sac <- get(obj)
objectList2$Lamprey@assays$RNA@data['LOC116942270',] = wang_sac@assays$RNA@data['FEZF2',Cells(objectList2$Lamprey)]
objectList2$Lamprey = RenameFeatures(objectList2$Lamprey, gsub('LOC116944204', 'TENM3',
gsub('TENM3', 'TENM4',
gsub('LOC116942270', 'FEZF2',
rownames(objectList2$Lamprey)))))
merged = merge(objectList2[[1]], objectList2[2:length(objectList2)])
merged # 11.5k SACs
merged$species = factor(merged$species, levels = rev(phylogenetic_order2))
# Switch to RNA assay
DefaultAssay(merged) = 'RNA'
# VlnPlot(merged, features = c('FEZF1', 'FEZF2', 'TENM3', 'ZFHX3', 'RND3'),
# split.by = 'lit_type', group.by = 'species', stack = TRUE, flip = FALSE,
# cols = c(ON = col1, OFF = col2, SAC = 'grey'))
# FeaturePlot(objectList2$Lamprey, 'TENM3')
# FeaturePlot(objectList2$Sauropsids, 'RND3')
# VlnPlot(objectList2$Lamprey, group.by = 'lit_type', features = 'TENM3')
col1 = 'deeppink2'
col2 = 'chartreuse2'
# Linear combination
merged$FEZF1_FEZF2 = merged@assays$RNA@data['FEZF1',] + merged@assays$RNA@data['FEZF2',]
merged$TENM3_ZFHX3_RND3 = merged@assays$RNA@data['TENM3',] + merged@assays$RNA@data['ZFHX3',] + merged@assays$RNA@data['RND3',]
merged$Subtype = merged$lit_type
merged$`Off.score` = (merged@assays$RNA@data['TENM3',] + merged@assays$RNA@data['ZFHX3',] + merged@assays$RNA@data['RND3',] -
merged@assays$RNA@data['FEZF1',] - merged@assays$RNA@data['FEZF2',]) / 5
# VlnPlot(merged, features = c('FEZF1_FEZF2', 'TENM3_ZFHX3_RND3'),
# split.by = 'lit_type', group.by = 'species', stack = TRUE, flip = FALSE,
# cols = c(ON = 'deeppink3', OFF = 'cyan3', SAC = 'grey'))
# plot_grid(VlnPlot(merged,
# features = c("FEZF1_FEZF2", "TENM3_ZFHX3_RND3"),
# group.by = "species",
# # cols = c(ON = col1, OFF = col2, SAC = 'grey'),
# split.by = "Subtype",
# split.plot = FALSE,
# pt.size = 100,
# stack = TRUE) +
# scale_fill_manual(values = c(ON = col1, OFF = col2, SAC = 'grey'))+
# NoLegend() +
# theme(axis.title.y = element_blank()),
# stackedBarGraph(merged, feature.2 = "species", feature.1 = "Subtype") +
# labs(y = 'Proportion')+
# scale_fill_manual(values = c(ON = col1, OFF = col2, SAC = 'grey'))+
# guides(fill = guide_legend(title = "Subtype"))+
# # scale_fill_manual(values = c(params$color1, params$color2))+
# coord_flip() +
# scale_y_continuous(expand = expansion(mult = c(0,0)))+
# guides(fill = guide_legend(ncol = 1))+
# theme(axis.text.y = element_blank(),
# legend.position = 'top',
# axis.title.y = element_blank(),
# axis.ticks.y = element_blank()),
# ncol = 2, rel_widths = c(2,1), align = "h", axis = "bt")
# ggsave('../../figures/my_figs/sac-violn-1.pdf', height = 8, width = 4.5)
# plt.list = VlnPlot(merged,
# features = c("FEZF1_FEZF2", "TENM3_ZFHX3_RND3"),
# group.by = "species",
# # cols = c(ON = col1, OFF = col2, SAC = 'grey'),
# split.by = "Subtype",
# split.plot = FALSE,
# pt.size = 0.01,
# stack = FALSE,
# combine = FALSE)
#
# wrap_plots(lapply(plt.list, function(x) x + coord_flip()))
#
# plt.list[[1]] + coord_flip()
col1 = 'orange2'
col2 = 'cyan3'
merged2 = subset(merged, species %in% c('Human', 'Macaque', 'Marmoset', 'MouseLemur'), invert = T)
merged2$species_pretty = factor(factor(species_names_pretty(merged2$species), levels = rev(phylogenetic_order_pretty)))
PrettyBarplot([email hidden], x = 'species_pretty', y = 'FEZF1_FEZF2', fill = 'Subtype',
add = c('mean_se'), position = position_dodge(width = 0.8), width = 0.8) +
scale_fill_manual(values = c(ON = col1, OFF = col2, SAC = 'grey')) +
ylab('FEZF1+\nFEZF2') +
scale_ticks_first_last()+
NoRotatedXAxis()+
theme(axis.title.x = element_text(face = 'italic'), axis.title.y = element_blank())+
coord_flip() + NoLegend() |
PrettyBarplot([email hidden], x = 'species_pretty', y = 'TENM3_ZFHX3_RND3', fill = 'Subtype',
add = c('mean_se'), position = position_dodge(width = 0.8), width = 0.8) +
scale_fill_manual(values = c(ON = col1, OFF = col2, SAC = 'grey')) +
ylab('TENM3+\nZFXH3+\nRND3') +
scale_ticks_first_last()+
NoRotatedXAxis()+
theme(axis.title.x = element_text(face = 'italic'), axis.title.y = element_blank())+
coord_flip() +
theme(axis.text.y = element_blank(), axis.title.y = element_blank(), legend.position = 'top') |
stackedBarGraph(merged2, feature.2 = "species_pretty", feature.1 = "Subtype", percent = F) +
labs(y = 'Percent')+
scale_fill_manual(values = c(ON = col1, OFF = col2, SAC = 'grey'))+
guides(fill = guide_legend(title = "Subtype"))+
# scale_fill_manual(values = c(params$color1, params$color2))+
coord_flip() +
NoLegend()+
NoRotatedXAxis()+
# scale_y_continuous(expand = expansion(mult = c(0,0)), labels = function(x) x*100)+
scale_ticks_first_last(do = function(x) x*100)+
# guides(fill = guide_legend(ncol = 1))+
theme(axis.text.y = element_blank(),
axis.title.y = element_blank(),
axis.ticks.y = element_blank())
ggsave('../../figures/my_figs/sac-violn-1.pdf', height = 6, width = 4.5)
```
# Figure 6: non-mammal integration
## Read in samap data
```{r, fig.height=8, fig.width=10}
#' Prepare umap metadata
class.umap = as.data.frame(fread('samap/pm_sh_ca_kf_ze_ne_am_ch_li-class-v2-umap.csv'))
# colnames(class.umap) = gsub('UMAP', 'UMAP_', colnames(class.umap))
class.umap$species_full = convert_values(class.umap$species, index$species %>% setNames(index$ident))
class.umap$species_full = factor(factor(class.umap$species_full, levels = index$species))
# ac.umap = as.data.frame(fread('samap/pm_sh_ca_kf_ze_am_ch_li-ac-v1-umap.csv')) # v1 used Wang annotations
# ac.umap = as.data.frame(fread('samap/pm_sh_ca_kf_ze_am_ch_li-ac-v2-umap.csv'))
ac.umap = as.data.frame(fread('samap/pm_sh_ca_kf_ze_ne_am_ch_li-ac-v3-umap.csv'))
# colnames(ac.umap) = gsub('UMAP', 'UMAP_', colnames(ac.umap))
ac.umap$species_full = convert_values(ac.umap$species, index$species %>% setNames(index$ident))
ac.umap$species_full = factor(factor(ac.umap$species_full, levels = index$species))
# rgc.umap = as.data.frame(fread('samap/pm_sh_ca_kf_ze_am_ch_li-class-v1-umap.csv'))
# colnames(rgc.umap) = gsub('UMAP', 'UMAP_', colnames(rgc.umap))
bc.umap = as.data.frame(fread('samap/pm_sh_ca_kf_ze_ne_am_ch_li-bc-v2-umap.csv'))
# colnames(bc.umap) = gsub('UMAP', 'UMAP_', colnames(bc.umap))
bc.umap$species_full = convert_values(bc.umap$species, index$species %>% setNames(index$ident))
bc.umap$species_full = factor(factor(bc.umap$species_full, levels = index$species))
prefix.key = index$species %>% setNames(index$ident)
```
## Cell class confusion matrix
```{r, fig.height=4.5, fig.width=16}
# Assign majority class
class.umap$leiden_clusters = class.umap$leiden_clusters_5
# Remove leiden cluster 33 (poorly integrated cells)
class.umap = subset(class.umap, leiden_clusters != 31)
class.umap$leiden_clusters = as.numeric(as.character(RenumberClustering(class.umap$leiden_clusters)))
cm = ConfusionMatrix(subset(class.umap, species_full != "Lamprey")$leiden_clusters,
subset(class.umap, species_full != "Lamprey")$cell_class, plot = FALSE)
majority.class = factor(apply(cm, 1, function(x) colnames(cm)[which.max(x)]),
levels = c('PR', 'onBC', 'offBC', 'HC', 'glyAC', 'gabaAC', 'AC', 'RGC', 'MG'))
cluster.order = seq(0, length(majority.class)-1)[rev(order(majority.class))]
bg = stackedBarGraph2(subset(class.umap, species_full != 'Lamprey') %>%
mutate(leiden_clusters = factor(as.character(leiden_clusters),
levels = (cluster.order))),
'leiden_clusters', 'cell_class', as.factor = TRUE) +
scale_y_continuous(expand = c(0, 0),
breaks = c(0, 0.25, 0.5, 0.75, 1), # only first and last
labels = c('0', '', '', '', '1'),
minor_breaks = waiver() # keeps automatic minor ticks
)+
coord_flip() +
scale_fill_manual(values = cell_class3_colors)+
theme_cowplot()+
theme(axis.text.y = element_blank(),
axis.title.y = element_blank(),
axis.ticks.y = element_blank(),
axis.title.x = element_text(size = 12),
axis.text.x = element_text(angle = 0, hjust = 0.5))
ident.2 = 'annotated'
ident.1 = 'leiden_clusters'
objectList = split(class.umap, class.umap$species_full)
heatmapList = lapply(seq_along(objectList), function(index) {
JSHeatmap2(objectList[[index]][[ident.1]],
objectList[[index]][[ident.2]],
title = names(objectList)[[index]],
row.order = rev(cluster.order),
border.col = NA,
max.value = 0.5,
stagger.threshold = 0.1) +
theme(axis.ticks = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
plot.margin = margin(0, 0, 0, 0))
#seq(0, max(class.umap$leiden_clusters)))
})
names(heatmapList) = names(objectList)
ncol = length(objectList)+2
nrow = 1
# ggarrange(plotlist = c(list(NULL),
# ggarrange(plotlist = heatmapList,
# widths = sapply(objectList, function(x) length(unique(x$annotated))),
# common.legend = TRUE), list(bg)),
# ncol = ncol,
# nrow = nrow,
# # common.legend = TRUE,
# widths = c(10, 1000, 50),
# legend = "none",
# align = 'hv',
# axis = 'btlr')
# ggarrange(plotlist = c(list(NULL), heatmapList, list(bg)))
# ggsave('../../figures/my_figs/nonmammal-integration.pdf', height=4, width=12)
# JSHeatmap2(objectList[[3]][[ident.1]],
# objectList[[3]][[ident.2]],
# title = names(objectList)[[3]],
# row.order = cluster.order)
ggarrange(plotlist = c(heatmapList, list(bg)),
widths = c(sapply(objectList, function(x) length(unique(x$annotated))), 40),
common.legend = TRUE,
align = 'h',
nrow = 1)
ggsave('../../figures/my_figs/nonmammal-class-jaccard.pdf', height=4.5, width=16)
```
## Cell class UMAPs
```{r, fig.height=4, fig.width=15}
PrettyUmap2(class.umap, group.by = 'species_full',
label = FALSE, geom.label = geom_text_repel,
box.padding = 0.7, show.legend = TRUE,
title = 'Non-mammal integration (species)', pt.alpha = 0.1,
cols = species_palette3, rasterise = TRUE,
nbreaks = 30, remove.times = 3) |
PrettyUmap2(subset(class.umap, species_full != 'Lamprey'), group.by = 'cell_class2',
label = TRUE, geom.label = geom_text_repel,
box.padding = 0.7,
title = 'Manually annotated cell classes (excluding lamprey)', pt.alpha = 0.1,
cols = cell_class3_colors, rasterise = TRUE,
nbreaks = 30, remove.times = 3) |
PrettyUmap2(subset(class.umap, species_full == 'Lamprey'), group.by = 'cell_class2',
label = TRUE, geom.label = geom_text_repel,
box.padding = 0.7,
title = 'Manually annotated cell classes (lamprey)', pt.alpha = 0.2,
cols = cell_class3_colors, rasterise = TRUE,
nbreaks = 30, remove.times = 3)
ggsave('../../figures/my_figs/nonmammal-class-umap.pdf', height=4, width=15)
```
```{r}
#' Tabulate errors
class.umap$cell_class2_inferred = MatchClusters(class.umap$leiden_clusters, class.umap$cell_class2)
accuracies = dapply(levels(class.umap$species_full), function(this.species){
sub = subset(class.umap, species_full == this.species)
mat = table(sub$cell_class2_inferred, sub$cell_class2)
diag(mat) = 0
# Print accuracy
(1 - sum(mat) / nrow(sub)) * 100
})
paste0(round(quantile(unlist(accuracies[-1]), 0.025), 2), '-', round(quantile(unlist(accuracies[-1]), 0.975), 2))
dapply(levels(class.umap$species_full), function(this.species){
sub = subset(class.umap, species_full == this.species)
table(sub$cell_class2_inferred, sub$cell_class2)
})
```
```{r, fig.height=10, fig.width=10}
# Cluster 33 is suspicious
class.umap$leiden_clusters = as.character(class.umap$leiden_clusters)
PrettyUmap2(class.umap,
group.by = 'leiden_clusters',
# color.by = 'cell_class2',
label = TRUE,
geom.label = geom_text_repel,
# nudge_x = 2, nudge_y = ,
title = 'Non-mammal integration', pt.alpha = 0.1,
# cols = cell_class3_colors,
cols = c(`33` = 'red'),
rasterise = FALSE, nbreaks = 30, remove.times = 3)
# Which species are in 33
table(subset(class.umap, leiden_clusters == 33)$species) # pretty even
sort(table(subset(class.umap, species == 'ze')$leiden_clusters,
subset(class.umap, species == 'ze')$annotated)['33',]) # ze_glyAC-43
sort(table(subset(class.umap, species == 'pm')$leiden_clusters,
subset(class.umap, species == 'pm')$annotated)['33',]) # pm_MG-4
sort(table(subset(class.umap, species == 'ne')$leiden_clusters,
subset(class.umap, species == 'ne')$annotated)['33',]) # multiple
sort(table(subset(class.umap, species == 'ca')$leiden_clusters,
subset(class.umap, species == 'ca')$annotated)['33',]) # multiple
# Seems like a cluster of poorly integrated cells more than contamination
# Shark HCs? 15, 21
JSHeatmap2(subset(class.umap, species == 'sh')$leiden_clusters,
subset(class.umap, species == 'sh')$annotated)
cm[as.character(cluster.order),]
# Cluster 39 (RGC) has some ACs
sort(table(subset(class.umap, species == 'ze')$leiden_clusters,
subset(class.umap, species == 'ze')$annotated)['39',]) # no cells
sort(table(subset(class.umap, species == 'ne')$leiden_clusters,
subset(class.umap, species == 'ne')$annotated)['39',]) # RGC
sort(table(subset(class.umap, species == 'sh')$leiden_clusters,
subset(class.umap, species == 'sh')$annotated)['39',]) # no cells
sort(table(subset(class.umap, species == 'li')$leiden_clusters,
subset(class.umap, species == 'li')$annotated)['39',]) # RGC
sort(table(subset(class.umap, species == 'ch')$leiden_clusters,
subset(class.umap, species == 'ch')$annotated)['39',]) # ch_gabaAC-34; no contamination tho
sort(table(subset(class.umap, leiden_clusters == '39')$annotated))
```
## Alignment matrix: class
```{r, fig.height=8, fig.width=9}
samap_alignment = as.matrix(read.csv("samap/pm_sh_ca_kf_ze_am_ch_li-class-v1-mappingtable.csv", row.names = 1))
pdf('../../figures/my_figs/nonmammal-integration-matrix.pdf', height=5, width=6)
sm.class = SAMapAlignmentHeatmap2(samap_alignment,
cell_class3_colors,
nm_palette2,
# width = unit(19, "in"),
# height = unit(19, "in"),
show_annotation_legend = TRUE,
show_row_names = FALSE,
show_column_names = FALSE,
species.use = c('pm', 'sh', 'kf', 'ca', 'ze', 'am', 'ch', 'li'),
type.order = CELL_CLASSES3,
max.value = 0.4,
cols = c('white', 'deeppink', 'deeppink4'),
rect_gp = gpar(col = "grey", lwd = 0),
use_raster = TRUE,
raster_quality = 2,
raster_by_magick = TRUE)
dev.off()
cnames = rownames([email hidden])
pdf('../../figures/my_figs/nonmammal-integration-matrix-class.pdf', height=50, width=50)
sm.class = SAMapAlignmentHeatmap2(samap_alignment,
cell_class3_colors,
nm_palette2,
show_annotation_legend = TRUE,
show_row_names = TRUE,
show_column_names = TRUE,
species.use = c('pm', 'sh', 'kf', 'ca', 'ze', 'am', 'ch', 'li'),
type.order = CELL_CLASSES3,
max.value = 0.4,
cols = c('white', 'deeppink', 'deeppink4'),
rect_gp = gpar(col = "grey", lwd = 0),
use_raster = TRUE,
raster_quality = 2,
raster_by_magick = TRUE)
dev.off()
```
## One nearest neighbor classification for lamprey
```{r, fig.height=5, fig.width=14}
# samap_alignment = (read.csv("samap/pm_sh_ca_kf_ze_am_ch_li-class-v1-mappingtable.csv", row.names = 1))
samap_alignment = (read.csv("samap/pm_sh_ca_kf_ze_ne_am_ch_li-class-v2-mappingtable.csv", row.names = 1))
species = str_split_fixed(rownames(samap_alignment), "_", 2)[,1]
types = str_split_fixed(str_split_fixed(rownames(samap_alignment), "_", 2)[,2], '-', 2)[,1]
idents = str_split_fixed(str_split_fixed(rownames(samap_alignment), "_", 2)[,2], '-', 2)[,2]
smatrix = SmartMatrix(as.matrix(samap_alignment), row.data = data.frame(species, types, idents), col.data = data.frame(species, types, idents))
# Set threshold of alignment score
# smatrix@matrix[smatrix@matrix < 0.1] = 0
score.threshold = 0.05
# For each lamprey cluster, find its nearest neighbor in each species
matrix.list = list(Shark = smatrix[[email hidden]$species == 'pm', [email hidden]$species == 'sh'],
Killifish = smatrix[[email hidden]$species == 'pm', [email hidden]$species == 'kf'],
Goldfish = smatrix[[email hidden]$species == 'pm', [email hidden]$species == 'ca'],
Zebrafish = smatrix[[email hidden]$species == 'pm', [email hidden]$species == 'ze'],
Newt = smatrix[[email hidden]$species == 'pm', [email hidden]$species == 'ne'],
Axolotl = smatrix[[email hidden]$species == 'pm', [email hidden]$species == 'am'],
Chicken = smatrix[[email hidden]$species == 'pm', [email hidden]$species == 'ch'],
Lizard = smatrix[[email hidden]$species == 'pm', [email hidden]$species == 'li'])
nn.res = do.call(cbind, lapply(matrix.list, function(smat)
apply(smat@matrix, 1, function(x) {
this.max = max(x)
if(this.max > score.threshold)
[email hidden]$type[which.max(x)]
else
'Unmapped'
})))
nn.res
tabulation = do.call(gtools::smartbind, apply(nn.res, 1, function(x) t(as.data.frame.vector(table(x)))))
tabulation[is.na(tabulation)] = 0
melted = reshape2::melt(as.matrix(tabulation)) %>% setNames(c('Cluster', 'Prediction', 'nSpecies'))
melted = subset(melted, Prediction != 'Unmapped')
max.value = 8; col.high = "#584B9FFF"; col.low = 'white'; legend_name = '# species';
if(!exists('objectList')) objectList = ReadSAMapObjects()
lamprey = objectList$Lamprey
# lamprey = OrderAnnotated(readRDS('../../Full_Objects/Lamprey_Wang_full_v1.rds'))
# Make dotplot first
pm2ch = readRDS('../cones/orthology_graphs/pm2ch.rds')
lamprey = ConvertGeneSymbols2(lamprey, pm2ch)
plt2 = AnnotatedUmap(lamprey,
annotation = subset(major_annotation_lamprey_on_off,
genes = c('PDE6H', 'RHO', 'VSX2', 'GRIK1', 'ONECUT3', 'SLC6A9', 'SLC6A5', 'PAX6',
'TFAP2B', 'GAD1', 'RBPMS', 'RBPMS2', 'SLC17A6', 'POU4F3', 'NEFL', 'NEFM', 'RLBP1')),
group.by = 'annotated',
color.clusters.by = 'cell_class2',
umap.legend = TRUE,
plot.umap = FALSE,
title = 'Lamprey')
# Make NN plot after, same order as dotplot
melted = subset(melted, Cluster %in% paste0('pm_', lamprey$annotated))
melted$Cluster = factor(melted$Cluster, levels = paste0('pm_', levels(plt2$data$id)))
melted$Prediction = factor(melted$Prediction, levels = rev(c('PR', 'onBC', 'offBC', 'HC', 'glyAC', 'gabaAC', 'RGC', 'MG')))
plt1 = ggplot(melted, aes(x = Cluster, y = Prediction))+
geom_point(aes(colour = nSpecies, size=nSpecies))+
scale_color_gradient(legend_name, low=col.low, high = col.high, limits=c(0, max.value), na.value = 'grey') +
scale_radius(legend_name, limits=c(0, max.value)) +
# scale_size(range = c(1, max.size), limits = c(0, max.perc))+
theme_bw() +
RotatedAxis() +
# ylab('True')+
# xlab('Predicted')+
theme(axis.title = element_text(color = 'black'),
axis.text = element_text(color = 'black'),
plot.title = element_text(hjust = 0.5))
plot_grid(plt1 + theme(axis.text.x = element_blank(), axis.title.x = element_blank()),
plt2, nrow = 2, align = 'v', axis = 'lr', rel_heights = c(1,2))
# ggsave('figures/lamprey-1NN-predictions-05.pdf', height=5, width=14)
```
Permissive: 9, 27, 30, 34, 35, 37, 39, 41
Stringent: 27, 34
```{r, fig.height=3, fig.width=10}
# Best alignment score within each species... 3 entries x 382 rows
best.scores = do.call(rbind, lapply(1:nrow(sm.class@matrix), function(i){
species = unique([email hidden]$species)
this.species = [email hidden][i,'species']
do.call(cbind, lapply(species, function(other.species){
max(sm.class@matrix[i,[email hidden]$species == other.species])
}) %>% setNames(species))
}))
best.scores.df = reshape2::melt(best.scores) %>% setNames(c('id', 'target', 'Score'))
best.scores.df$species = rep([email hidden]$species, 4)
best.scores.df$type = rep([email hidden]$types, 4)
ggboxplot(subset(best.scores.df, type %in% c('PR', 'HC', 'MG')), x = 'type', y = 'Score', color = 'species', legend = 'none') |
ggboxplot(subset(best.scores.df, !type %in% c('PR', 'HC', 'MG')), x = 'type', y = 'Score', color = 'species', legend = 'right')
```
### Heatmap variant for SFN
```{r, fig.width=4, fig.height=6}
# melted2$Cluster = factor(melted2$Cluster,
# levels = unique(arrange(melted2, Prediction, nSpecies)$Cluster))
tabulation2 = tabulation[,colnames(tabulation) != 'Unmapped']
tabulation2$BC = tabulation2$onBC + tabulation2$offBC
tabulation2 = tabulation2[,!colnames(tabulation2) %in% c('onBC', 'offBC')]
best.match = apply(tabulation2, 1, function(x) colnames(tabulation2)[which.max(x)])
rgc.fraction = apply(tabulation2, 1, function(x) x['RGC'] / 7) #max(x)/7)
ac.fraction = apply(tabulation2, 1, function(x) 1 - x['gabaAC'] / 7)
sort.df = data.frame(Cluster = rownames(tabulation2),
best.match = best.match,
rgc.fraction = rgc.fraction,
ac.fraction = ac.fraction)
melted2 = reshape2::melt(as.matrix(tabulation2)) %>% setNames(c('Cluster', 'Prediction', 'nSpecies'))
# melted2 = subset(melted, Prediction != 'Unmapped')
melted2$Prediction = gsub('gabaAC', 'gaAC', gsub('glyAC', 'glAC', melted2$Prediction))
melted2$Prediction = factor(melted2$Prediction, levels = (c('PR', 'BC', 'HC', 'glAC', 'gaAC', 'RGC', 'MG')))
# melted2$Cluster = factor(melted2$Cluster, levels = rev(arrange(sort.df,
# factor(best.match,
# levels = (c('PR', 'BC', 'HC', 'glyAC', 'gabaAC', 'RGC', 'MG'))),
# rgc.fraction,
# ac.fraction)$Cluster))
pm_order = c("pm_MG-3", "pm_MG-2", "pm_MG-1", "pm_RGC-34", "pm_RGC-30", "pm_RGC-37", "pm_RGC-27",
"pm_RGC-9", "pm_RGC-39", "pm_RGC-35", "pm_RGC-3", "pm_RGC-41",
"pm_RGC-1", "pm_RGC-33", "pm_MG-4", "pm_RGC-26", "pm_gabaAC-3", "pm_RGC-22",
"pm_RGC-23", "pm_RGC-7", "pm_RGC-31", "pm_RGC-20", "pm_RGC-2", "pm_RGC-18",
"pm_RGC-11", "pm_RGC-8", "pm_RGC-19", "pm_RGC-4", "pm_RGC-36", "pm_RGC-28",
"pm_RGC-5", "pm_RGC-21", "pm_RGC-15", "pm_gabaAC-4", "pm_gabaAC-2",
"pm_gabaAC-1", "pm_RGC-25", "pm_RGC-24", "pm_RGC-12", "pm_RGC-6",
"pm_RGC-40", "pm_RGC-38", "pm_RGC-32", "pm_RGC-29", "pm_RGC-17",
"pm_RGC-16", "pm_RGC-14", "pm_RGC-10", "pm_RGC-13", "pm_glyAC-4", "pm_offBC-6",
"pm_glyAC-9", "pm_glyAC-8", "pm_glyAC-7", "pm_glyAC-6", "pm_glyAC-5",
"pm_glyAC-3", "pm_glyAC-2", "pm_glyAC-11", "pm_glyAC-10", "pm_glyAC-1",
"pm_HC-4", "pm_HC-3", "pm_HC-2", "pm_HC-1", "pm_onBC-5",
"pm_onBC-3", "pm_onBC-1", "pm_offBC-8", "pm_offBC-7", "pm_offBC-4",
"pm_offBC-2", "pm_PR-2", "pm_PR-1")
melted2$Cluster <- factor(melted2$Cluster,
levels = pm_order
)
# Remove suspicious clusters
sus.clusters = c('pm_offBC-6', 'pm_MG-4', 'pm_gabaAC-3', 'pm_gabaAC-4')
melted2 = subset(melted2, !Cluster %in% sus.clusters)
# Plot
ggplot(melted2, aes(y = Cluster, x = Prediction))+
geom_tile(aes(fill = nSpecies))+
scale_fill_gradient(legend_name, low='white', high = 'deeppink2', limits=c(0, max.value), na.value = 'grey',
guide = guide_colorbar(frame.colour = "black", ticks.colour = "black")) +
# scale_radius(legend_name, limits=c(0, max.value)) +
# scale_size(range = c(1, max.size), limits = c(0, max.perc))+
theme_dario() +
scale_x_discrete(expand = expansion(mult = c(0,0)))+
RotatedAxis() +
ArialFont()+
# ylab('True')+
# xlab('Predicted')+
theme(axis.title = element_text(color = 'black'),
axis.text = element_text(color = 'black'),
plot.title = element_text(hjust = 0.5),
axis.ticks.y = element_blank(),
axis.text.x = element_text(size = 12))
ggsave('../../figures/my_figs/lamprey-1NN-heatmap-v2.pdf', width=3.9, height=5)
# ggsave('figures/lamprey-1NN-heatmap.png')
```
Check for doublet clusters in lamprey
```{r, fig.height=8, fig.width=15}
objectList$Lamprey = FindDoublets(objectList$Lamprey, channels = 'animal')
DoubletAnalysis(objectList$Lamprey, group.by = 'annotated')
BrowseSeurat(objectList$Lamprey)
VlnPlot(objectList$Lamprey, features = 'nCount_RNA', group.by = 'annotated')
```
Check ambiguous clusters
```{r, fig.height=4.8, fig.width=4.2}
objectList = ReadSAMapObjects(new.lamprey = TRUE)
pm2ch = readRDS('../cones/orthology_graphs/pm2ch.rds')
objectList$Lamprey = ConvertGeneSymbols2(objectList$Lamprey, pm2ch)
# ClusterFeaturePlot(objectList$Lamprey, group.by = 'annotated', features = c('TFAP2B', 'RBPMS2'), cols = c('white', 'deeppink2'))
# clusters.of.interest = gsub('pm_', '', c("pm_RGC-34", "pm_RGC-30", "pm_RGC-27", "pm_RGC-37", "pm_RGC-9", "pm_RGC-35", "pm_RGC-41", "pm_RGC-3", "pm_RGC-39", "pm_RGC-1", "pm_RGC-33", "pm_RGC-26", "pm_RGC-31", "pm_RGC-23", "pm_gabaAC-4", "pm_RGC-7"))
clusters.of.interest = gsub('pm_', '', c("pm_RGC-27", "pm_RGC-9", "pm_RGC-39", "pm_RGC-35", "pm_RGC-3",
"pm_RGC-41", "pm_RGC-1", "pm_RGC-33", "pm_RGC-26",
"pm_RGC-22", "pm_RGC-23", "pm_RGC-7", "pm_RGC-31"))
# pm.ambig = subset(objectList$Lamprey, annotated %in% clusters.of.interest)
# pm.ambig$annotated = factor(pm.ambig$annotated, levels = rev(clusters.of.interest))
# objectList$Lamprey$annotated = factor(objectList$Lamprey$annotated, levels = gsub('pm_', '', rev(c(tail(pm_order, -6), pm_order[1:6]))))
objectList$Lamprey$annotated = factor(objectList$Lamprey$annotated,
levels = rev(unique(c(clusters.of.interest, levels(objectList$Lamprey$annotated)))))
sus.clusters = c('pm_offBC-6', 'pm_MG-4', 'pm_gabaAC-3', 'pm_gabaAC-4')
DotPlot3(subset(objectList$Lamprey, annotated %in% gsub('pm_', '', sus.clusters), invert = T),
features = c('TFAP2B', 'GAD1', 'RBPMS', 'RBPMS2', 'NEFL', 'NEFM',
'OPN4', 'LOC116943951', 'NOS1', 'NPY'),
# idents = clusters.of.interest,
show = clusters.of.interest,
coord.flip = FALSE,
group.by = 'annotated',
col.low = 'grey90',
col.high = 'deeppink3')
ggsave('figures/fig-lamprey-amiguous-dotplots.pdf', height=4.8, width=4.2)
```
```{r, fig.height=5, fig.width=20}
DotPlot3(subset(objectList$Lamprey, annotated %in% gsub('pm_', '', sus.clusters), invert = T),
features = c('TFAP2B', 'GAD1', 'RBPMS', 'RBPMS2', 'NEFL', 'NEFM',
'OPN4', 'LOC116943951', 'NOS1', 'NPY'),
# idents = clusters.of.interest,
# show = clusters.of.interest,
coord.flip = TRUE,
group.by = 'annotated',
col.low = 'grey90',
col.high = 'deeppink3')
```
## Gene heatmaps
```{r, fig.height=5, fig.width=10}
nmList = ReadSAMapObjects()
nmList = lapply(nmList, DownsampleSeurat, group.by = 'annotated', size = 50)
nmList = sapply(names(nmList), function(species){
nmList[[species]]$annotated = paste0(index$ident[index$species == species], '_', nmList[[species]]$annotated)
nmList[[species]]
}, USE.NAMES = TRUE)
nmList$Zebrafish = NormalizeData(nmList$Zebrafish)
# nmList$Killifish$annotated = factor(nmList$Killifish$annotated)
```
```{r, fig.height=5, fig.width=10}
gene_pairs = as.data.frame(fread('samap/pm_sh_ca_kf_ze_am_ch_li-class-v1-genepairs.csv')[,-1])
markers.dfs = lapply(setdiff(CELL_CLASSES3, 'Other'), function(class){
df_onecol = pivot_longer(head(gene_pairs[,grepl(class, colnames(gene_pairs)) &
!grepl('pval', colnames(gene_pairs))], 600),
everything(),
values_to = "value")
if(class %in% c('onBC', 'offBC', 'gabaAC', 'glyAC', 'RGC'))
df_onecol = df_onecol[df_onecol$value %in% names(which(table(df_onecol$value) >= 2)), ]
df_onecol$species1 = ExtractString(ExtractString(df_onecol$value, after = ';'), after = '_')
df_onecol$gene1 = ExtractString(ExtractString(df_onecol$value, after = ';'), before = '_')
df_onecol$species2 = ExtractString(ExtractString(df_onecol$value, before = ';'), after = '_')
df_onecol$gene2 = ExtractString(ExtractString(df_onecol$value, before = ';'), before = '_')
# browser()
if(class == 'RGC') browser()
Reduce(function(dtf1,dtf2) merge(dtf1, dtf2, by="Chicken", all = TRUE), list(
subset(df_onecol, species1 == 'ch' & species2 == 'pm')[,c('gene1', 'gene2')] %>% setNames(c('Chicken', 'Lamprey')),
subset(df_onecol, species1 == 'ch' & species2 == 'sh')[,c('gene1', 'gene2')] %>% setNames(c('Chicken', 'Shark')),
subset(df_onecol, species1 == 'ch' & species2 == 'kf')[,c('gene1', 'gene2')] %>% setNames(c('Chicken', 'Killifish')),
subset(df_onecol, species1 == 'ch' & species2 == 'ca')[,c('gene1', 'gene2')] %>% setNames(c('Chicken', 'Goldfish')),
subset(df_onecol, species1 == 'ch' & species2 == 'ze')[,c('gene1', 'gene2')] %>% setNames(c('Chicken', 'Zebrafish')),
subset(df_onecol, species1 == 'ch' & species2 == 'am')[,c('gene1', 'gene2')] %>% setNames(c('Chicken', 'Axolotl')),
subset(df_onecol, species1 == 'ch' & species2 == 'li')[,c('gene1', 'gene2')] %>% setNames(c('Chicken', 'Lizard'))
)) %>% unique
})
markers.df = do.call(rbind, lapply(markers.dfs, function(x) {
# Each gene can appear no more than n times
x = x %>%
group_by(Chicken) %>%
slice_head(n = 10) %>%
ungroup() %>% as.data.frame()
# if(nrow(x) >= 1000) x[sample(seq_len(nrow(x)), 1000),] else x[sample(seq_len(nrow(x))),]
# if(nrow(x) >= 1000) x[1:1000,] else x
# Scramble order to make prettier
x[sample(seq_len(nrow(x))),]
}))
markers.df$row_annotation = ExtractString(rownames(markers.df), after = '\\.')
table(markers.df$row_annotation)
types.use = gsub('PR-L-M', 'PR-L/M', gsub('\\.', '-', cnames))
types.order = paste0(prefix.key[ExtractString(types.use, after = '_')], ' ', types.use)
# Check types
stopifnot(all(types.use %in% unlist(lapply(nmList, function(x) unique(x$annotated)))))
setdiff(types.use, unlist(lapply(nmList, function(x) unique(x$annotated))))
setdiff(unlist(lapply(nmList, function(x) unique(x$annotated))), types.use)
# Check convergence (10:1)
length(unique(markers.df$Chicken))
length(unique(markers.df$Lizard))
length(unique(markers.df$Zebrafish))
length(unique(markers.df$Killifish))
# pdf('../../figures/my_figs/nonmammal-integration-genes.pdf', height=5, width=9)
SAMapHeatmap3(nmList,
markers.df,
type_palette = major_annotation_palette2,
species_palette = species_palette2,
species.use = c("Chicken", "Lizard", 'Killifish', "Zebrafish"),
types.use = types.use, #rownames(samap_alignment),
show_heatmap_legend = FALSE,
show_row_names = FALSE,
show_column_names = FALSE,
min.z.score = -1,
max.z.score = 2,
col_fun = circlize::colorRamp2(c(-1, 0, 2), c("white", "white", "deeppink")),
# width = unit(10, "in"),
rotate = TRUE,
types.order = types.order, #rownames(samap_alignment),
use_raster = TRUE)
# dev.off()
```
Troubleshooting
```{r}
SAMapHeatmap3(nmList,
head(subset(markers.df, row_annotation == 'gabaAC'), 20),
type_palette = major_annotation_palette2,
species_palette = species_palette2,
species.use = c("Chicken", "Lizard", 'Killifish', "Zebrafish"),
types.use = types.use, #rownames(samap_alignment),
show_heatmap_legend = FALSE,
show_row_names = FALSE,
show_column_names = TRUE,
min.z.score = -2,
max.z.score = 2,
col_fun = circlize::colorRamp2(c(-2, 0, 2), c("white", "white", "deeppink")),
# width = unit(10, "in"),
rotate = TRUE,
types.order = types.order, #rownames(samap_alignment),
use_raster = TRUE)
```
## Clustered heatmap
```{r, fig.height=50, fig.width=50}
pdf('figures/nonmammals-clustered-heatmap.pdf', height=50, width=50)
Heatmap(sm.class@matrix)
Heatmap(sm.bc@matrix)
dev.off()
```
## Confusion matrix: BC
```{r, fig.height=4, fig.width=15}
# Assign majority class
bc.umap$leiden_clusters = bc.umap$leiden_clusters_2
# bc.umap$leiden_clusters[bc.umap$leiden_clusters == 10] = 6
cm = ConfusionMatrix(bc.umap$leiden_clusters, bc.umap$cell_class, plot = FALSE)
# majority.class = factor(apply(cm, 1, function(x) colnames(cm)[which.max(x)]),
# levels = c('PR', 'onBC', 'offBC', 'HC', 'glyAC', 'gabaAC', 'AC', 'RGC', 'MG'))
# How many onBC clusters?
n.clusters = length(which(cm[,'onBC'] > 0.5))
cluster.order = rev(rownames(cm %>% arrange(offBC)))
bg = stackedBarGraph2(bc.umap %>%
mutate(leiden_clusters = factor(as.character(leiden_clusters),
levels = rev(cluster.order))),
'leiden_clusters', 'cell_class', as.factor = TRUE) +
scale_y_continuous(expand = c(0, 0),
breaks = c(0, 0.25, 0.5, 0.75, 1), # only first and last
labels = c('0', '', '', '', '1'),
minor_breaks = waiver() # keeps automatic minor ticks
)+
coord_flip() +
scale_fill_manual(values = cell_class3_colors)+
theme_cowplot()+
NoLegend()+
theme(axis.text.y = element_blank(),
axis.title.y = element_blank(),
axis.ticks.y = element_blank(),
axis.title.x = element_text(size = 12),
axis.text.x = element_text(angle = 0, hjust = 0.5))
ident.2 = 'annotated'
ident.1 = 'leiden_clusters'
ob.list = split(bc.umap, bc.umap$species_full)
heatmapList = lapply(seq_along(ob.list), function(index) {
JSHeatmap2(ob.list[[index]][[ident.1]],
ob.list[[index]][[ident.2]],
title = names(ob.list)[[index]],
row.order = cluster.order,
border.col = NA,
# col.high = 'chartreuse4',
max.value = 0.5,
stagger.threshold = 0.15) +
geom_hline(yintercept = n.clusters+0.5, linetype = 'dashed', color = 'grey')+
theme(axis.ticks = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
plot.margin = margin(0, 0, 0, 0),
plot.background = element_rect(fill = "transparent", colour = NA)) +
NoLegend()
#seq(0, max(class.umap$leiden_clusters)))
})
names(heatmapList) = names(ob.list)
ncol = length(ob.list)+2
nrow = 1
# plt.umap = PrettyUmap2(bc.umap,
# group.by = 'species_full',
# label = FALSE,
# geom.label = geom_text_repel,
# cols = species_palette3, #major_annotation_palette2,
# show.legend = FALSE,
# title = 'BC integration',
# pt.alpha = 0.3,
# rasterise = FALSE, pad = 0) +
# LegendLowerRight(x = 0.99, y = 0.1)
#
# plt.umap2 = PrettyUmap2(bc.umap,
# group.by = 'cell_class2',
# label = TRUE,
# geom.label = geom_text_repel,
# cols = cell_class3_colors,
# show.legend = FALSE,
# title = 'BC subclass',
# pt.alpha = 0.3,
# rasterise = FALSE,
# pad = 0) +
# LegendLowerRight(x = 0.7, y = 0.6) +
# theme(axis.title.y = element_blank(),
# plot.margin = unit(c(0,0,0,0), 'cm'))
#
# ggarrange(plotlist = c(list(plt.umap),
# list(plt.umap2),
# list(NULL), heatmapList, list(bg)),
# ncol = ncol,
# nrow = nrow,
# # common.legend = TRUE,
# widths = c(40, 35, 8, sapply(ob.list, function(x) length(unique(x$annotated))), 13),
# # legend = "none",
# align = 'h')
ggarrange(plotlist = c(list(NULL), heatmapList, list(bg)),
ncol = ncol,
nrow = nrow,
common.legend = TRUE,
widths = c(8, sapply(ob.list, function(x) length(unique(x$annotated))), 13),
# legend = "none",
align = 'h')
ggsave('../../figures/my_figs/nonmammal-bc-jaccard.pdf', height=3.5, width=15)
```
Find the errors for each species
```{r, fig.height=5, fig.width=20}
bc.umap$cell_class2_inferred = MatchClusters(bc.umap$leiden_clusters, bc.umap$cell_class2)
EvaluateModel(table(bc.umap$cell_class2_inferred, bc.umap$cell_class2))
dapply(levels(bc.umap$species_full), function(this.species){
sub = subset(bc.umap, species_full == this.species)
EvaluateModel(table(sub$cell_class2_inferred, sub$cell_class2))
})
# Goldfish has most errors
# What are these clusters being inferred as
subset(bc.umap, species_full == 'Goldfish')[,c('annotated', 'cell_class2', 'cell_class2_inferred')] %>% unique
table(subset(bc.umap, species_full == 'Goldfish')$cell_class2,
subset(bc.umap, species_full == 'Goldfish')$annotated)
table(subset(bc.umap, species_full == 'Goldfish')$cell_class2_inferred,
subset(bc.umap, species_full == 'Goldfish')$annotated)
# ca_onBC-13, ca_onBC-12 is flipped, ca_offBC-10 is mixed (might be low quality or doublets)
# 12 and 13 seem OFF
goldfish = readRDS('../../Full_Objects/Goldfish_full_v2.rds')
DotPlot3(subset(goldfish, cell_class == 'BC'), features = c('isl1', 'grm6a', 'grm6b', 'grik1a', 'grik1b'), group.by = 'annotated')
VlnPlot(goldfish, features = 'nFeature_RNA', group.by = 'annotated', pt.size = 0) + NoLegend()
# Check for doublets
goldfish = FindDoublets(goldfish, 'orig.ident')
DoubletAnalysis(goldfish, group.by = 'annotated')
# Other classes
class.umap$cell_class2_inferred = MatchClusters(class.umap$leiden_clusters, class.umap$cell_class2)
dapply(levels(bc.umap$species_full), function(this.species){
sub = subset(class.umap, species_full == this.species)
table(sub$cell_class2_inferred, sub$cell_class2)
})
table(subset(class.umap, species_full == 'Goldfish')$cell_class2_inferred,
subset(class.umap, species_full == 'Goldfish')$cell_class2)
# Updated to goldfish v3
# Check newt; ne_offBC-17 (GRIK1+), onBC-18, 20 as ISL1 low GRIK high so moving to offBC
table(subset(class.umap, species_full == 'Newt')$cell_class2_inferred,
subset(class.umap, species_full == 'Newt')$annotated)
table(subset(bc.umap, species_full == 'Newt')$cell_class2_inferred,
subset(bc.umap, species_full == 'Newt')$annotated)
# Check killifish; kf_offBC-17 (mixed), kf_onBC-19 (off, not on), kf_onBC-7 (mixed)
table(subset(bc.umap, species_full == 'Killifish')$cell_class2_inferred,
subset(bc.umap, species_full == 'Killifish')$annotated)
# Check kf for doublets
kf = readRDS('../../Species_Objects/Killifish_ncbi_initial_v5.rds')
kf = FindDoublets(kf, 'orig.file')
DoubletAnalysis(kf, group.by = 'annotated')
# Not many doublets found, keeping v5 as is
```
```{r}
ob.list = split(bc.umap, bc.umap$species_full)
lapply(seq_along(ob.list), function(index) {
JSHeatmap2(ob.list[[index]][[ident.1]],
ob.list[[index]][[ident.2]],
title = names(ob.list)[[index]],
row.order = cluster.order,
# border.col = NA,
col.high = 'chartreuse4',
max.value = 0.5,
stagger.threshold = 0.15) +
theme(axis.ticks = element_blank(),
# axis.text.x = element_blank(),
# axis.text.y = element_blank(),
plot.margin = margin(0, 0, 0, 0)) +
NoLegend()
})
```
RBC:
pm onBC-5
sh onBC-3
ze onBC-14
```{r}
PrettyUmap2(bc.umap %>% mutate(leiden_clusters = as.character(leiden_clusters)), group.by = 'leiden_clusters')
PrettyUmap2(subset(bc.umap, species == 'ch'), group.by = 'annotated')
PrettyUmap2(subset(bc.umap, species == 'li'), group.by = 'annotated')
PrettyUmap2(subset(bc.umap, species == 'ze'), group.by = 'annotated')
PrettyUmap2(subset(bc.umap, species == 'kf'), group.by = 'annotated')
```
```{r, fig.width=10, fig.height=3}
DotPlot3(subset(nmList$Chicken, cell_class == 'BC'), features = c('PRKCA'))
DotPlot3(nmList$Lizard, features = c('PRKCA'))
DotPlot3(nmList$Zebrafish, features = c('gramd1bb', 'gramd1ba', 'prkcaa', 'isl1'))
# LOC107381579 (megf11) LOC107387807 (chat)
DotPlot3(nmList$Killifish, features = c('gramd1bb', 'gramd1ba', 'prkcaa', 'isl1a','syt5b', 's100a10b'), group.by = 'annotated')
```
Addition of mammals
```{r}
bc.umap2 = as.data.frame(fread('samap/kf_ze_ch_li_mf-BC-v1-umap.csv'))
colnames(bc.umap2) = gsub('UMAP', 'UMAP_', colnames(bc.umap2))
bc.umap2$leiden_clusters = as.character(bc.umap2$leiden_clusters_1.5)
bc.umap2$leiden_clusters[bc.umap2$leiden_clusters == '5'] = '0'
PrettyUmap2(bc.umap2, group.by = 'leiden_clusters')
ob.list = split(bc.umap2, bc.umap2$species)
lapply(seq_along(ob.list), function(index) {
JSHeatmap2(ob.list[[index]][['leiden_clusters']],
ob.list[[index]][['annotated']],
title = names(ob.list)[[index]],
row.order = sort(unique(bc.umap2$leiden_clusters)),
# border.col = NA,
col.high = 'chartreuse4',
max.value = 0.5,
stagger.threshold = 0.15) +
theme(axis.ticks = element_blank(),
# axis.text.x = element_blank(),
# axis.text.y = element_blank(),
plot.margin = margin(0, 0, 0, 0)) +
NoLegend()
})
```
Validate SAMap signatures
```{r, fig.height=10, fig.width=15}
gene_pairs = read.csv('samap/kf_ze_ch_li_mf-BC-v1-genepairs.csv', row.names = 1, check.names = FALSE)
gene_pairs = gene_pairs[,!grepl('pval', colnames(gene_pairs))]
pdf('samap/kf_ze_ch_li_mf-BC-v1-genepairs.pdf', width = 15, height = 15)
lapply(colnames(gene_pairs), function(cname){
# x = gsub('BC.', 'BC-', cname)
x = cname
t1 = ExtractString(x, after = ';')
t2 = ExtractString(x, before = ';')
s1 = ExtractString(t1, after = '_')
s2 = ExtractString(t2, after = '_')
species1 = prefix.key[s1]
species2 = prefix.key[s2]
object1 = nmList[[species1]]
object2 = nmList[[species2]]
features1 = setdiff(ExtractString(ExtractString(gene_pairs[,cname], after = ';'), before = '_'), '')
features2 = setdiff(ExtractString(ExtractString(gene_pairs[,cname], before = ';'), before = '_'), '')
print(paste0(species1, '-', species2))
TitlePlot(DotPlot3(subset(object1, cell_class == 'BC'), features = features1), paste0(species1, '-', t1)) + NoLegend() |
TitlePlot(DotPlot3(subset(object2, cell_class == 'BC'), features = features2), paste0(species2, '-', t2))
})
dev.off()
DotPlot3(subset(nmList$Chicken, cell_class == 'BC'), features = 'PRKCA')
DotPlot3(subset(nmList$Lizard, cell_class == 'BC'), features = 'PRKCA')
```
## UMAP: BCs
```{r, fig.height=4, fig.width=15}
PrettyUmap2(bc.umap, group.by = 'species_full',
label = FALSE, geom.label = geom_text_repel,
box.padding = 0.7, show.legend = TRUE,
title = 'BC integration (species)', pt.alpha = 0.2,
cols = species_palette3, rasterise = TRUE,
nbreaks = 30, remove.times = 3) |
PrettyUmap2(subset(bc.umap, species_full != 'Lamprey'), group.by = 'cell_class2',
label = FALSE, geom.label = geom_text_repel,
box.padding = 2, show.legend = TRUE,
title = 'Manually annotated BC subclasses (excluding lamprey)', pt.alpha = 0.2,
cols = cell_class3_colors, rasterise = TRUE,
nbreaks = 30, remove.times = 3) + LegendLowerLeft() |
PrettyUmap2(subset(bc.umap, species_full == 'Lamprey'), group.by = 'cell_class2',
label = FALSE, geom.label = geom_text_repel,
box.padding = 2, show.legend = TRUE,
title = 'Manually annotated BC subclasses (lamprey)', pt.alpha = 0.4,
cols = cell_class3_colors, rasterise = TRUE,
nbreaks = 30, remove.times = 3) + LegendLowerLeft()
ggsave('../../figures/my_figs/nonmammal-bc-umap.pdf', height=4, width=15)
```
## Confusion matrix: AC
```{r, fig.height=8, fig.width=10}
# Assign majority class
ac.umap$leiden_clusters = factor(as.numeric(factor(ac.umap$leiden_clusters_4.5)))
message('Found ', length(unique(ac.umap$leiden_clusters)), ' clusters!')
# Remove tiny clusters
table(ac.umap$leiden_clusters)
ac.umap = subset(ac.umap, !leiden_clusters %in% names(which(table(ac.umap$leiden_clusters) < 5)))
# Merge similar clusters
# Heatmap(CoclusteringMatrix(factor(subset(ac.umap, species == 'ch')$leiden_clusters),
# subset(ac.umap, species == 'ch')$annotated))
# Heatmap(CoclusteringMatrix(factor(subset(ac.umap, species == 'ze')$leiden_clusters),
# subset(ac.umap, species == 'ze')$annotated))
# Averged across all species
Heatmap(CoclusteringMatrix(factor(ac.umap$leiden_clusters),
ac.umap$annotated))
# clust = hclust(dist(1-CoclusteringMatrix(factor(subset(ac.umap, species == 'ch')$leiden_clusters),
# subset(ac.umap, species == 'ch')$annotated)), method = 'average')
# plot(clust)
# Merge 12 and 26 with resolution 4
# ac.umap$leiden_clusters[ac.umap$leiden_clusters == 26] = 12
# Merge 35/7, 20/2, 12/4, 5/6, 24/13 with resolution 4.5
# ac.umap$leiden_clusters[ac.umap$leiden_clusters == 35] = 7
# ac.umap$leiden_clusters[ac.umap$leiden_clusters == 20] = 2
# ac.umap$leiden_clusters[ac.umap$leiden_clusters == 12] = 4
# ac.umap$leiden_clusters[ac.umap$leiden_clusters == 6] = 5
# ac.umap$leiden_clusters[ac.umap$leiden_clusters == 24] = 13
# ac.umap$leiden_clusters = factor(ac.umap$leiden_clusters)
# Merge 32 and 12 with resolution 5
# Merge 35 and 11 with resolution 4.5
ac.umap$leiden_clusters[ac.umap$leiden_clusters == 35] = 11
ac.umap$leiden_clusters = RenumberClustering(ac.umap$leiden_clusters, zero.base = FALSE)
table(ac.umap$leiden_clusters)
```
Supplementary figure
```{r}
# ac.umap = as.data.frame(fread('samap/pm_sh_ca_kf_ze_ne_am_ch_li-ac-v3-umap.csv'))
# ac.umap$species_full = convert_values(ac.umap$species, index$species %>% setNames(index$ident))
# ac.umap$species_full = factor(factor(ac.umap$species_full, levels = index$species))
# prefix.key = index$species %>% setNames(index$ident)
ac.umap$cell_class2 = gsub('gabaAC', 'gaAC', gsub('glyAC', 'glAC', ac.umap$cell_class2))
major_annotation_palette3 = major_annotation_palette2
names(major_annotation_palette3) = gsub('gabaAC', 'gaAC', gsub('glyAC', 'glAC', names(major_annotation_palette3)))
PrettyUmap2(ac.umap,
group.by = 'cell_class2',
label = FALSE,
geom.label = geom_text_repel,
cols = major_annotation_palette3,
show.legend = TRUE,
title = 'Non-mammalian AC integration (subclass)',
pt.alpha = 0.2,
remove.times = 3,
rasterise = FALSE,
pad = 0) +
LegendTopLeft() +
theme(axis.title.y = element_blank(),
plot.margin = unit(c(0,0,0,0), units = 'cm')) |
PrettyUmap2(ac.umap,
group.by = 'leiden_clusters',
label = TRUE,
geom.label = geom_text_repel,
# cols = species_palette2,
show.legend = FALSE,
title = 'Non-mammalian oACs',
pt.alpha = 0.2,
remove.times = 3,
rasterise = FALSE,
min.segment.length = 0.1,
size = 2.5,
pad = 0)
ggsave('../../figures/my_figs/nonmammalian-oacs-umap.pdf', height = 3.2, width = 8)
```
```{r, fig.height=3.5, fig.width=18}
#' Assign majority class, plot umap, and plot jaccard matrices
# ac.umap$leiden_clusters = factor(as.numeric(factor(ac.umap$leiden_clusters_4.5)))
message('Found ', length(unique(ac.umap$leiden_clusters)), ' clusters!')
cm = ConfusionMatrix(ac.umap$leiden_clusters, ac.umap$cell_class2, plot = FALSE)
cluster.order = rev(rownames(cm %>% arrange(glyAC)))
majority.class = factor(apply(cm, 1, function(x) colnames(cm)[which.max(x)]), levels = c('PR', 'onBC', 'offBC', 'HC', 'glyAC', 'gabaAC', 'AC', 'RGC', 'MG'))
# apply(cm[cluster.order, ], 1, max)
bg = stackedBarGraph2(ac.umap %>%
mutate(leiden_clusters = factor(as.character(leiden_clusters),
levels = rev(cluster.order))),
'leiden_clusters', 'cell_class2', as.factor = TRUE) +
scale_y_continuous(expand = c(0, 0),
breaks = c(0, 0.25, 0.5, 0.75, 1), # only first and last
labels = c('0', '', '', '', '1'),
minor_breaks = waiver() # keeps automatic minor ticks
)+
coord_flip() +
scale_fill_manual(values = major_annotation_palette2)+
theme_cowplot()+
NoLegend()+
theme(axis.text.y = element_blank(),
axis.title.y = element_blank(),
axis.ticks.y = element_blank(),
axis.title.x = element_text(size = 12),
axis.text.x = element_text(angle = 0, hjust = 0.5, size = 10))
ident.2 = 'annotated'
ident.1 = 'leiden_clusters'
ob.list = split(ac.umap, ac.umap$species_full)
heatmapList = lapply(seq_along(ob.list), function(index) {
JSHeatmap2(ob.list[[index]][[ident.1]],
ob.list[[index]][[ident.2]],
title = names(ob.list)[[index]],
row.order = cluster.order,
border.col = NA,
# col.high = 'cyan4',
max.value = 0.5,
stagger.threshold = 0.10) +
# geom_hline(yintercept = difference(majority.class[rev(cluster.order)]) - 0.5, linetype = 'solid', color = 'grey', linewidth = 0.4, alpha = 0.4)+
theme(axis.ticks = element_blank(),
axis.text.x = element_blank(),
axis.text.y = element_blank(),
plot.margin = margin(0, 0, 0, 0),
plot.background = element_rect(fill = "transparent", colour = NA)) +
NoLegend()
#seq(0, max(class.umap$leiden_clusters)))
})
names(heatmapList) = names(ob.list)
ncol = length(ob.list)+3
nrow = 1
ac.umap$species_full2 = as.character(ac.umap$species_full)
plt.umap = PrettyUmap2(ac.umap,
group.by = 'species_full',
label = FALSE,
geom.label = geom_text_repel,
cols = species_palette2,
show.legend = TRUE,
title = 'Non-mammal integration',
pt.alpha = 0.2,
remove.times = 3,
rasterise = FALSE,
pad = 0) +
LegendLowerRight(size = 8, line.spacing = 0.7, x = 1.01)
plt.umap2 = PrettyUmap2(ac.umap,
group.by = 'cell_class2',
label = FALSE,
geom.label = geom_text_repel,
cols = major_annotation_palette2,
show.legend = TRUE,
title = 'AC subclasses',
pt.alpha = 0.2,
rasterise = FALSE,
pad = 0) +
LegendLowerLeft() +
theme(axis.title.y = element_blank(),
plot.margin = unit(c(0,0,0,0), units = 'cm'))
ggarrange(plotlist = c(list(plt.umap),
# list(plt.umap2),
list(NULL), heatmapList, list(bg)),
ncol = ncol,
nrow = nrow,
# common.legend = TRUE,
widths = c(100, 10, sapply(ob.list, function(x) length(unique(x$annotated))), 22),
# legend = "none",
align = 'h')
# ggsave('../../figures/my_figs/nonmammal-integration-AC-noline.pdf', height=3.5, width=18)
```
```{r, fig.height=5, fig.width=40}
# ac.umap$leiden_clusters = factor(as.numeric(factor(ac.umap$leiden_clusters_4.5)))
# ac.umap$leiden_clusters[ac.umap$leiden_clusters == 10] = 6
cm = ConfusionMatrix(ac.umap$leiden_clusters, ac.umap$cell_class, plot = FALSE)
# majority.class = factor(apply(cm, 1, function(x) colnames(cm)[which.max(x)]),
# levels = c('PR', 'onBC', 'offBC', 'HC', 'glyAC', 'gabaAC', 'AC', 'RGC', 'MG'))
cluster.order = rev(rownames(cm %>% arrange(glyAC)))
majority.class = factor(apply(cm, 1, function(x) colnames(cm)[which.max(x)]), levels = c('PR', 'onBC', 'offBC', 'HC', 'glyAC', 'gabaAC', 'AC', 'RGC', 'MG'))
ident.2 = 'annotated'
ident.1 = 'leiden_clusters'
ob.list = split(ac.umap, ac.umap$species_full)
heatmapList = lapply(seq_along(ob.list), function(index) {
JSHeatmap2(ob.list[[index]][[ident.1]],
ob.list[[index]][[ident.2]],
title = names(ob.list)[[index]],
row.order = cluster.order,
border.col = 'black',
# col.high = 'cyan4',
max.value = 0.5,
stagger.threshold = 0.10) +
geom_hline(yintercept = difference(majority.class[rev(cluster.order)]) - 0.5, linetype = 'dashed', color = 'grey')+
theme(axis.ticks = element_blank(),
#axis.text.x = element_blank(),
#axis.text.y = element_blank(),
plot.margin = margin(0, 0, 0, 0)) +
NoLegend()
#seq(0, max(class.umap$leiden_clusters)))
})
plot_grid(plotlist = heatmapList, nrow = 1, align = 'h', axis = 'bt', rel_widths = sapply(ob.list, function(x) length(unique(x$annotated))))
ggsave('figures/nonmammal-ac-jaccard-v3.pdf')
```
```{r, fig.height=5, fig.width=30}
#' Find the errors
ac.umap$cell_class2_inferred = MatchClusters(ac.umap$leiden_clusters, ac.umap$cell_class2)
EvaluateModel(table(ac.umap$cell_class2_inferred, ac.umap$cell_class2))
dapply(levels(ac.umap$species_full), function(this.species){
sub = subset(ac.umap, species_full == this.species)
message(this.species)
EvaluateModel(table(sub$cell_class2_inferred, sub$cell_class2))
})
# errors are rare
# Zebrafish has one cluster ze_glyAC-18 but it is solidly slc6a9+ and gad1/2-
table(subset(ac.umap, species == 'ze')$cell_class2_inferred,
subset(ac.umap, species == 'ze')$annotated)
ze = readRDS('../../Full_Objects/Zebrafish_full_v3.rds')
DotPlot3(ze, group.by = 'annotated', features = c('gad1a', 'gad1b', 'gad2', 'slc6a9', 'slc6a5'))
```
pm_RGC-13 seems glycinergic, expresses GlyT2
pm_RGC-4
pm_RGC-19
Near perfect 1:1 correspondence between chicken and zebrafish
```{r}
PrettyUmap2(subset(ac.umap, species_full == 'Lamprey'),
group.by = 'cell_class2',
label = FALSE,
geom.label = geom_text_repel,
cols = major_annotation_palette2,
show.legend = FALSE,
title = 'AC subclass',
pt.alpha = 1,
rasterise = FALSE,
pad = 0)
PrettyUmap2(subset(ac.umap, species_full == 'Killifish'),
group.by = 'cell_class2',
label = FALSE,
geom.label = geom_text_repel,
cols = major_annotation_palette2,
show.legend = FALSE,
title = 'AC subclass',
pt.alpha = 1,
rasterise = FALSE,
pad = 0)
```
```{r, fig.height=8, fig.width=30}
PrettyUmap2(ac.umap,
group.by = 'leiden_clusters',
label = TRUE,
geom.label = geom_text_repel,
# cols = major_annotation_palette2,
show.legend = FALSE,
title = 'AC subclass',
pt.alpha = 0.2,
rasterise = FALSE,
pad = 0)
# ConfusionMatrix(ac.umap$leiden_clusters, ac.umap$annotated, plot = FALSE)['27',] %>% unlist %>% sort %>% tail # oAC37 ortholog
# ConfusionMatrix(ac.umap$leiden_clusters, ac.umap$annotated, plot = FALSE)['23',] %>% unlist %>% sort %>% tail # VG3 ortholog
# ConfusionMatrix(ac.umap$leiden_clusters, ac.umap$annotated, plot = FALSE)['29',] %>% unlist %>% sort %>% tail # VG3 ortholog
plot_grid(
plotlist = lapply(seq_along(ob.list), function(index) {
JSHeatmap2(ob.list[[index]][[ident.1]],
ob.list[[index]][[ident.2]],
title = names(ob.list)[[index]],
row.order = cluster.order,
border.col = NA,
col.high = 'cyan4',
max.value = 0.5,
stagger.threshold = 0.15) +
theme(axis.ticks = element_blank(),
# axis.text.x = element_blank(),
# axis.text.y = element_blank(),
plot.margin = margin(0, 0, 0, 0)) +
NoLegend()
#seq(0, max(class.umap$leiden_clusters)))
}),
align = 'h',
rel_widths = sapply(ob.list, function(x) length(unique(x$annotated))),
nrow = 1)
```
## Alignment matrix: AC
```{r, fig.height=8, fig.width=9}
ac_alignment = as.matrix(read.csv("samap/pm_sh_ca_kf_ze_ne_am_ch_li-ac-v2-mappingtable.csv", row.names = 1))
pdf('../../figures/my_figs/nonmammal-integration-matrix-ac.pdf', height=45, width=46)
sm.class = SAMapAlignmentHeatmap2(ac_alignment,
cell_class3_colors,
nm_palette2,
# width = unit(19, "in"),
# height = unit(19, "in"),
show_annotation_legend = TRUE,
show_row_names = TRUE,
show_column_names = TRUE,
species.use = c('pm', 'sh', 'kf', 'ca', 'ze', 'ne', 'am', 'ch', 'li'),
type.order = CELL_CLASSES3,
max.value = 0.4,
cols = c('white', 'deeppink', 'deeppink4'),
rect_gp = gpar(col = "grey", lwd = 0),
use_raster = TRUE,
raster_quality = 2,
raster_by_magick = TRUE)
dev.off()
cnames = rownames([email hidden])
```
## Orthotype correspondence
```{r, fig.height=5, fig.width=6}
if(!exists('ac.ortho')) ac.ortho = LoadACOrtho()
ac.umap$orthotypes_li = ac.ortho$type[paste0(ac.umap$species_full, '_', apply(str_split_fixed(ac.umap$V1, '-', 3)[,1:2], 1, paste, collapse = '-'))]
ac.umap$orthotypes_ch = ac.ortho$type[paste0(ac.umap$species_full, '_', ExtractString(ac.umap$V1, after = '-'))]
ac.umap$orthotypes = ifelse(is.na(ac.umap$orthotypes_li), as.character(ac.umap$orthotypes_ch), as.character(ac.umap$orthotypes_li))
ac.umap$orthotypes = orthotype_labels(ac.umap$orthotypes, as.factor = FALSE)
# Tabulation and counting
js.mat = JSMatrix(table(factor(ac.umap$orthotypes, ORTHOTYPES), ac.umap$leiden_clusters))
apply(js.mat, 1, function(x) max(x) > 0.25)
over.stats = OverlapStatistics(table(factor(ac.umap$orthotypes, ORTHOTYPES), ac.umap$leiden_clusters))
unique(subset(over.stats, padj < 0.05)$ident1)
# ac.umap$leiden_clusters = ac.umap$leiden_clusters_4
# JSHeatmap2(subset(ac.umap, species == 'ch')$orthotypes_ch, subset(ac.umap, species == 'ch')$leiden_clusters, stagger.threshold = 0.1)
# JSHeatmap2(subset(ac.umap, species == 'ch')$orthotypes_ch, subset(ac.umap, species == 'ch')$annotated, stagger.threshold = 0.1)
# JSHeatmap2(subset(ac.umap, species == 'li')$orthotypes_li, subset(ac.umap, species == 'li')$leiden_clusters, stagger.threshold = 0.1)
# JSHeatmap2(ac.umap$orthotypes, ac.umap$leiden_clusters, stagger.threshold = 0.1, row.order = levels(ac.ortho$orthotype))
ggarrange(
JSHeatmap2(factor(ac.umap$orthotypes, ORTHOTYPES),
ac.umap$leiden_clusters,
stagger.threshold = 0.1,
row.order = ORTHOTYPES,
border.col = NA,
title = '') +
# coord_flip()+
ArialFont()+
theme(axis.text.y = element_text(size = 5, color = 'black',
# margin = margin(b = -200)
# margin = unit(c(0,0,0,0), 'cm')
),
axis.ticks.length.y = unit(1, "pt"),
axis.ticks = element_blank(),
axis.text.x = element_blank(), #element_text(size = 5, color = 'black'),
plot.background = element_rect(fill = "transparent", colour = NA)) +
scale_y_discrete(position = "right"),
common.legend = TRUE,
legend = 'none')
ggsave('../../figures/my_figs/orthotype-correspondence.pdf', height = 3.5, width = 3.5)
js.dat = JSHeatmap2(ac.umap$orthotypes,
ac.umap$leiden_clusters,
stagger.threshold = 0.1,
# row.order = levels(ac.ortho$type),
border.col = NA,
title = '')
saveRDS(js.dat, '../../Ortho_Objects/js.dat.rds')
overlap.list = OverlapStatistics(table(ac.umap$orthotypes, ac.umap$leiden_clusters))
lapply(unique(overlap.list$ident1), function(i) CorrespondingCluster(overlap.list, cluster = i))
key = MatchClusters(ac.umap$orthotypes, ac.umap$leiden_clusters, return.key = TRUE)
table(MatchClusters(subset(ac.umap, species %in% c('ch', 'li'))$orthotypes,
subset(ac.umap, species %in% c('ch', 'li'))$leiden_clusters, return.key = FALSE))
# Bidirectional best hits (25)
bbh = BidirectionalBestHits(JSMatrix(table(ac.umap$orthotypes, ac.umap$leiden_clusters)))
bbh = subset(cbind(bbh, do.call(rbind, apply(bbh, 1, function(x) subset(overlap.list, ident1 == x['row'] & ident2 == x['col'])))), overlap > 9)
saveRDS(bbh, '../../Ortho_Objects/bbh_v2.rds')
```
Unordered for SFN
```{r, fig.height=4, fig.width=4}
custom.order = c(
"oAC1 [A2]", "oAC2 [VG3]", "oAC3 [TRHDE]", "oAC4 [TRHDE]", "oAC5",
"oAC6 [A8]", "oAC7 [SEG]", "oAC8 [NNgly]", "oAC9", "oAC10 [ROBO3]",
"oAC11 [SLC35D3]", "oAC12", "oAC13 [A17]", "oAC14 [A17]", "oAC15", "oAC20 [NPY]",
"oAC16* [MAF]", "oAC21* [MAF]", "oAC17*", "oAC29", "oAC18 [nNOS]", "oAC19 [CA2]",
"oAC22*", "oAC23*", "oAC24", "oAC25*", "oAC26", "oAC27",
"oAC28", "oAC30 [PDGFRA]", "oAC31", "oAC32", "oAC33 [VIP]",
"oAC34 [RXRG]", "oAC35", "oAC36* [NNgaba]", "oAC37* [CRH]",
"oAC38* [nNOS]", "oAC39* [CA1]", "oAC40 [nNOS]", "oAC41",
"oAC42* [SAC]"
)
ggarrange(
JSHeatmap2(ac.umap$orthotypes,
ac.umap$leiden_clusters,
stagger.threshold = 0.1,
row.order = custom.order,
border.col = NA) +
# coord_flip()+
ArialFont()+
theme(axis.text.y = element_text(size = 5, color = 'black'),
axis.ticks = element_blank(),
plot.title = element_blank(),
axis.text.x = element_blank(), #element_text(size = 5, color = 'black'),
plot.background = element_rect(fill = "transparent", colour = NA)) +
scale_y_discrete(position = "right"),
common.legend = TRUE,
legend = 'none')
ggsave('../../figures/my_figs/orthotype-correspondence-unordered.pdf', height = 3.5, width = 3.5)
```
Need to run cells above (confusion AC, ac alignment, and bbh before this)
```{r, fig.height=20, fig.width=20}
# Split by species
cm.ac = do.call(cbind, lapply(seq_along(ob.list), function(index) JSMatrix(table(ob.list[[index]]$leiden_clusters, ob.list[[index]]$annotated))))
# cm.ac = ConfusionMatrix(ac.umap$leiden_clusters, ac.umap$annotated, plot = FALSE, row.norm = TRUE)
ac.groups = lapply(seq_len(nrow(cm.ac)), function(i) colnames(cm.ac)[which(cm.ac[i,] > 0.15)]) # jaccard threshold
ac.groups.sorted = ac.groups[as.numeric((levels(js.dat$data$col)))]
names(ac.groups.sorted) = levels(js.dat$data$col)
saveRDS(ac.groups.sorted, 'samap/ac.groups.sorted.rds')
# Remove singletons? Removing any cluster who's max alignment score was less than some threshold
singletons = rownames(ac_alignment)[apply(ac_alignment, 1, function(x) max(x) < 0.20)]
ac.groups.sorted = lapply(ac.groups.sorted, function(x) setdiff(x, singletons))
ac.groups = lapply(ac.groups, function(x) setdiff(x, singletons))
saveRDS(ac.groups, 'samap/ac.groups.rds')
```
```{r, fig.height=4, fig.width=20}
# NNgaba
sm.ac@matrix['li_gabaAC-36',] %>% sort
sm.ac@matrix['li_gabaAC-32',] %>% sort
sm.ac@matrix['ch_gabaAC-62',] %>% sort
# BHLHE22+
sm.ac@matrix['ze_glyAC-18',] %>% sort
sm.ac@matrix['ze_gabaAC-39',] %>% sort
DotPlot3(nmList$Chicken, features = c('BHLHE22', 'CHAT', 'SLC18A3', 'CADPS2', 'KCNH5', 'MAP7D2', 'CARMIL1', 'PCDH11X', 'ARHGAP21'), group.by = 'annotated')
DotPlot3(nmList$Lizard, features = c('BHLHE22', 'CHAT', 'SLC18A3', 'CADPS2', 'KCNH5', 'MAP7D2', 'CARMIL1', 'PCDH11X', 'ARHGAP21'))
DotPlot3(nmList$Zebrafish, features = c('bhlhe22', 'chata', 'slc18a3a', 'cadps2', 'kcnh5a', 'map7d2a', 'foxp4', 'carmil3', 'pcdh11', 'arhgap23b'))
DotPlot3(nmList$Killifish, features = c('bhlhe22', 'chata', 'slc18a3a', 'cadps2', 'kcnh5a', 'map7d2a', 'foxp4', 'carmil3', 'pcdh11', 'arhgap23b'), group.by = 'annotated')
DotPlot3(nmList$Goldfish, features = c('bhlhe22', 'chata', 'slc18a3a', 'cadps2', 'kcnh5a', 'map7d2a', 'foxp4', 'carmil3', 'pcdh11', 'arhgap23b'), group.by = 'annotated')
# PDGFRA
sm.ac@matrix['li_gabaAC-35',] %>% sort
sm.ac@matrix['ch_gabaAC-20',] %>% sort
sm.ac@matrix['ze_gabaAC-50_nNOS',] %>% sort
sm.ac@matrix['ze_gabaAC-30',] %>% sort
sm.class@matrix['ze_gabaAC-50_nNOS',] %>% sort
sm.class@matrix['ze_gabaAC-30',] %>% sort
DotPlot3(nmList$Killifish, features = c('pdgfra', 'ddc'))
DotPlot3(nmList$Zebrafish, features = c('pdgfra', 'ddc'))
DotPlot3(nmList$Chicken, features = 'PDGFRA')
DotPlot3(nmList$Lizard, features = 'PDGFRA')
kfish.ac = Harmonize(subset(nmList$Killifish, cell_class == 'AC'), 'orig.file')
FeaturePlot(kfish.ac, features = 'pdgfra', order = TRUE)
# li_CADPS2;ze_cadps2
# li_KCNH5;ze_kcnh5a
# li_MAP7D2;ze_map7d2a
# li_FOXP1;ze_foxp4
# li_CARMIL1;ze_carmil3
# li_PCDH11X;ze_pcdh11
# li_ARHGAP21;ze_arhgap23b
```
### Figure S11
```{r, fig.height=4, fig.width=15}
ac.umap = readRDS('../../Ortho_Objects/ac.umap.final.rds')
ac.umap$orthotypes_li = ac.ortho$type[paste0(ac.umap$species_full, '_', apply(str_split_fixed(ac.umap$V1, '-', 3)[,1:2], 1, paste, collapse = '-'))]
ac.umap$orthotypes_ch = ac.ortho$type[paste0(ac.umap$species_full, '_', ExtractString(ac.umap$V1, after = '-'))]
ac.umap$orthotypes = ifelse(is.na(ac.umap$orthotypes_li), as.character(ac.umap$orthotypes_ch), as.character(ac.umap$orthotypes_li))
ac.umap$orthotypes = orthotype_labels(ac.umap$orthotypes, as.factor = FALSE)
# Tabulation and counting
js.mat = JSMatrix(table(factor(ac.umap$orthotypes, ORTHOTYPES), ac.umap$leiden_clusters))
p1 = JSHeatmap2(factor(ac.umap$orthotypes, ORTHOTYPES),
ac.umap$leiden_clusters,
stagger.threshold = 0.1,
row.order = ORTHOTYPES,
border.col = NA,
title = paste0('Bridge: chicken and lizard', '\n ARI: ',
round(adj.rand.index(factor(ac.umap$orthotypes, ORTHOTYPES), ac.umap$leiden_clusters), 2))) +
# coord_flip()+
ArialFont()+
theme(axis.text.y = element_text(size = 6, color = 'black',
# margin = margin(b = -200)
# margin = unit(c(0,0,0,0), 'cm')
),
axis.ticks.length.y = unit(1, "pt"),
axis.ticks = element_blank(),
axis.text.x = element_blank(), #element_text(size = 5, color = 'black'),
plot.background = element_rect(fill = "transparent", colour = NA)) +
scale_y_discrete(position = "right") + NoLegend()
# res.chli = process_ortho_heatmap('samap/pm_sh_ca_kf_ze_ne_am_ch_li-ac-v1-umap.csv',
# stagger.threshold = 0.1,
# cluster_col = "leiden_clusters_4.5")
res.ch = process_ortho_heatmap('samap/pm_sh_ca_kf_ze_ne_am_ch-ac-v1-umap.csv',
stagger.threshold = 0.1,
cluster_col = "leiden_clusters_5",
title_label = 'Bridge: chicken')
res.chop = process_ortho_heatmap('samap/pm_sh_ca_kf_ze_ne_am_ch_op-ac-v1-umap.csv',
stagger.threshold = 0.1,
cluster_col = "leiden_clusters_4",
title_label = 'Bridge: chicken and opossum')
res.chliop = process_ortho_heatmap('samap/pm_sh_ca_kf_ze_ne_am_ch_li_op-ac-v1-umap.csv',
stagger.threshold = 0.15,
cluster_col = "leiden_clusters_4.5",
title_label = 'Bridge: chicken, lizard, opossum')
res.mfms = process_ortho_heatmap('samap/pm_sh_kf_ze_am_ch_li_op_mm_mf-ac-v1-umap.csv',
stagger.threshold = 0.15,
cluster_col = "leiden_clusters_4",
title_label = 'Bridge: chicken, lizard, \nopossum, mouse, macaque')
p1 | res.ch$plot | res.chop$plot | res.chliop$plot | res.mfms$plot
ggsave('../../figures/my_figs/samap-robustness.pdf')
# res.ch$plot | p2
```
```{r}
ac.umap.chliop = as.data.frame(fread('samap/pm_sh_ca_kf_ze_ne_am_ch-ac-v1-umap.csv')) # chicken and opossum
ac.umap.chliop$leiden_clusters = ac.umap.chliop$leiden_clusters_5
ac.umap.chliop$species_full = convert_values(ac.umap.chliop$species, index$species %>% setNames(index$ident))
ac.umap.chliop$species_full = factor(factor(ac.umap.chliop$species_full, levels = index$species))
ac.umap.chliop$orthotypes_li = ac.ortho$type[paste0(ac.umap.chliop$species_full, '_', apply(str_split_fixed(ac.umap.chliop$V1, '-', 3)[,1:2], 1, paste, collapse = '-'))]
ac.umap.chliop$orthotypes_ch = ac.ortho$type[paste0(ac.umap.chliop$species_full, '_', ExtractString(ac.umap.chliop$V1, after = '-'))]
ac.umap.chliop$orthotypes = ifelse(is.na(ac.umap.chliop$orthotypes_li), as.character(ac.umap.chliop$orthotypes_ch), as.character(ac.umap.chliop$orthotypes_li))
ac.umap.chliop$orthotypes = orthotype_labels(ac.umap.chliop$orthotypes, as.factor = FALSE)
# Tabulation and counting
js.mat.chliop = JSMatrix(table(factor(ac.umap.chliop$orthotypes, ORTHOTYPES), ac.umap.chliop$leiden_clusters))
# over.stats = OverlapStatistics(table(factor(ac.umap.noliz$orthotypes, ORTHOTYPES), ac.umap.noliz$leiden_clusters))
ggarrange(
JSHeatmap2(factor(ac.umap.chliop$orthotypes, ORTHOTYPES),
ac.umap.chliop$leiden_clusters,
stagger.threshold = 0.15,
row.order = ORTHOTYPES,
border.col = NA,
title = 'Bridge: chicken, lizard, and opossum') +
# coord_flip()+
ArialFont()+
theme(axis.text.y = element_text(size = 6, color = 'black',
# margin = margin(b = -200)
# margin = unit(c(0,0,0,0), 'cm')
),
axis.ticks.length.y = unit(1, "pt"),
axis.ticks = element_blank(),
axis.text.x = element_blank(), #element_text(size = 5, color = 'black'),
plot.background = element_rect(fill = "transparent", colour = NA)) +
scale_y_discrete(position = "right"),
common.legend = TRUE,
legend = 'none')
```
```{r, fig.height= 4, fig.width=16}
nrow(BidirectionalBestHits(js.mat, verbose = F))
nrow(BidirectionalBestHits(res.ch$matrix, verbose = F))
nrow(BidirectionalBestHits(res.chop$matrix, verbose = F))
nrow(BidirectionalBestHits(res.chliop$matrix, verbose = F))
nrow(BidirectionalBestHits(res.mfms$matrix, verbose = F))
best.jac.df = data.frame(orthotype = get_numbers(rownames(js.mat)),
best.chli = apply(js.mat, 1, max),
best.ch = apply(res.ch$matrix, 1, max),
best.chop = apply(res.chop$matrix, 1, max),
best.chliop = apply(res.chliop$matrix, 1, max),
best.chliopmmmf = apply(res.mfms$matrix, 1, max)
)
ScatterPlot(best.jac.df, x = 'best.ch', y = 'best.chli', y_equals_x = T, labels = 'orthotype', max.overlaps = 10,
ylab = 'Best jaccard score\n(chicken and lizard)', xlab = 'Best jaccard score\n(chicken)', r = T) |
ScatterPlot(best.jac.df, x = 'best.chop', y = 'best.chli', y_equals_x = T, labels = 'orthotype', max.overlaps = 10,
ylab = 'Best jaccard score\n(chicken and lizard)', xlab = 'Best jaccard score\n(chicken and opossum)', r = T) |
ScatterPlot(best.jac.df, x = 'best.chliop', y = 'best.chli', y_equals_x = T, labels = 'orthotype', max.overlaps = 10,
ylab = 'Best jaccard score\n(chicken and lizard)', xlab = 'Best jaccard score\n(chicken, lizard, and opossum)', r = T) |
ScatterPlot(best.jac.df, x = 'best.chliopmmmf', y = 'best.chli', y_equals_x = T, labels = 'orthotype', max.overlaps = 10,
ylab = 'Best jaccard score\n(chicken and lizard)', xlab = 'Best jaccard score \n(chicken, lizard, opossum, mouse, macaque)', r = T)
ggsave('../../figures/my_figs/conservation-robustness.pdf')
```
```{r, fig.height=8, fig.width=14}
res.mfms$data$leiden_clusters = as.character(res.mfms$data$leiden_clusters)
(TitlePlot((PrettyUmap(ac.ortho, group.by = 'species', alpha = 0.6, pt.size = 0.1, label = F, show.legend = T) +
guides(color = guide_legend(byrow = TRUE)) + theme(legend.key.height = unit(12, "pt"))), 'Seurat CCA (17 amniotes)') |
TitlePlot(PrettyUmap(ac.ortho, group.by = 'classification', alpha = 0.6, pt.size = 0.1), 'Subclass') |
TitlePlot(PrettyUmap(ac.ortho, group.by = 'orthotype.no', alpha = 0.6, pt.size = 0.1), 'Mean silhouette: 0.6')) /
(TitlePlot((PrettyUmap(res.mfms$data, group.by = 'species_full', alpha = 0.6, pt.size = 0.1, label = F, show.legend = T) +
guides(color = guide_legend(byrow = TRUE)) + theme(legend.key.height = unit(12, "pt"))), 'SAMap (10 vertebrates)') |
TitlePlot(PrettyUmap(res.mfms$data, group.by = 'cell_class2', alpha = 0.6, pt.size = 0.1), 'Subclass') |
TitlePlot(PrettyUmap(res.mfms$data, group.by = 'leiden_clusters', alpha = 0.6, pt.size = 0.1), 'Mean silhouette: 0.2'))
ggsave('../../figures/my_figs/samap-lower-res.pdf')
```
```{r, fig.height=6, fig.width=7}
JSHeatmap2(res.mfms$data$orthotypes, res.mfms$data$leiden_clusters, row.order = ORTHOTYPES,
xlab = 'SAMap vertebrate leiden clusters', ylab = 'Seurat amniote orthotypes', ari = T)
```
```{r, fig.height=5, fig.width=10}
seurat = readRDS('../../Ortho_Objects/vertebrateAC_seurat.rds')
seurat = MapBack(seurat, ac.ortho, group.by = 'orthotype.no')
seurat.ds = DownsampleSeurat(seurat[,!is.na(seurat$orthotype.no)], group.by = 'orthotype.no', size = 1000)
# seurat.ds = seurat[,!is.na(seurat$orthotype.no)]
# dist = 1-seurat@graphs$integrated_snn[Cells(seurat.ds), Cells(seurat.ds)]
# Computing silhouette in UMAP space because neither PCA nor graph space are comparable between Seurat and SAMap
silhouttes = FindSilhouette(seurat.ds,
group.by = 'orthotype.no',
# dist = dist,
reduction = 'umap',
method = 'Euclidean',
average = T
)
silhouttes
paste0('Mean cluster silhouette Seurat: ', mean(silhouttes$sil_width))
res.mfms
res.mfms.ds = res.mfms$data[sample(1:nrow(res.mfms$data), 10000),]
sil.samap = FindSilhouette(res.mfms.ds,
group.by = 'leiden_clusters',
dist = dist(res.mfms.ds[,c('UMAP_1', 'UMAP_2')]),
# reduction = 'umap',
# method = 'Euclidean',
average = T
)
paste0('Mean cluster silhouette SAMap: ', mean(sil.samap$sil_width))
```
## Checkpoint: Heatmap 1
```{r, fig.height=3, fig.width=10, dpi=300}
nmList = ReadSAMapObjects()
nmList = lapply(nmList, DownsampleSeurat, group.by = 'annotated', size = 50)
nmList = sapply(names(nmList), function(species){
nmList[[species]]$annotated = paste0(index$ident[index$species == species], '_', nmList[[species]]$annotated)
nmList[[species]]
}, USE.NAMES = TRUE)
nmList$Zebrafish = NormalizeData(nmList$Zebrafish)
# Newest version of files
gene_pairs = fread('samap/pm_sh_ca_kf_ze_ne_am_ch_li-ac-v3-genepairs.csv')[,-1] %>% data.frame()
# gene_pairs = fread('samap/kf_ze_ch_li-AC-v5-genepairs.csv')[,-1] %>% data.frame()
# gene_pairs = fread('samap/pm_sh_ca_kf_ze_am_ch_li-ac-v2-genepairs.csv')[,-1] %>% data.frame()
bbh = readRDS('../../Ortho_Objects/bbh_v2.rds')
ac.groups = readRDS('samap/ac.groups.rds')
ac.groups.sorted = readRDS('samap/ac.groups.sorted.rds')
# Get markers for each group of cell types!
marker.dfs = lapply(seq_along(ac.groups), function(i){
this.group = ac.groups[[i]]
message('Working on group ', i)
# print(this.group)
this.species = unique(ExtractString(this.group, after = '_'))
if(length(this.group) < 2) return(NULL)
col.names = gsub('-', '\\.', unlist(Iterate(sort(this.group), function(x,y) paste0(x, '.', y), return.matrix = FALSE)))
col.names = intersect(unlist(col.names), colnames(gene_pairs))
if(length(col.names) == 0) return(NULL)
df_onecol = pivot_longer(gene_pairs[,col.names,drop = FALSE],
everything(),
values_to = "value")
df_onecol$species1 = ExtractString(ExtractString(df_onecol$value, after = ';'), after = '_')
df_onecol$gene1 = ExtractString(ExtractString(df_onecol$value, after = ';'), before = '_')
df_onecol$species2 = ExtractString(ExtractString(df_onecol$value, before = ';'), after = '_')
df_onecol$gene2 = ExtractString(ExtractString(df_onecol$value, before = ';'), before = '_')
df_onecol = subset(df_onecol, value != '') # Remove rows with empty values
# Chicken symbols and lizard symbols
chicken.symbols = unique(c(subset(df_onecol, species1 == 'ch')$gene1, subset(df_onecol, species2 == 'ch')$gene2))
lizard.symbols = unique(c(subset(df_onecol, species1 == 'li')$gene1, subset(df_onecol, species2 == 'li')$gene2))
# Projection onto chicken
marker.df.ch = Reduce(function(dtf1,dtf2) full_join(dtf1, dtf2, by="Chicken", relationship = "many-to-many"), list(
unique(subset(df_onecol, species1 == 'ch' & species2 == 'pm')[,c('gene1', 'gene2')] %>% setNames(c('Chicken', 'Lamprey'))),
unique(subset(df_onecol, species1 == 'ch' & species2 == 'sh')[,c('gene1', 'gene2')] %>% setNames(c('Chicken', 'Shark'))),
unique(subset(df_onecol, species1 == 'ch' & species2 == 'kf')[,c('gene1', 'gene2')] %>% setNames(c('Chicken', 'Killifish'))),
unique(subset(df_onecol, species1 == 'ca' & species2 == 'ch')[,c('gene1', 'gene2')] %>% setNames(c('Goldfish', 'Chicken'))),
unique(subset(df_onecol, species1 == 'ch' & species2 == 'ze')[,c('gene1', 'gene2')] %>% setNames(c('Chicken', 'Zebrafish'))),
unique(subset(df_onecol, species1 == 'am' & species2 == 'ch')[,c('gene1', 'gene2')] %>% setNames(c('Axolotl', 'Chicken'))),
unique(subset(df_onecol, species1 == 'ch' & species2 == 'ne')[,c('gene1', 'gene2')] %>% setNames(c('Chicken', 'Newt'))),
unique(subset(df_onecol, species1 == 'ch' & species2 == 'li')[,c('gene1', 'gene2')] %>% setNames(c('Chicken', 'Lizard')))
)) %>% unique
# Projection onto lizard
marker.df.li = Reduce(function(dtf1,dtf2) full_join(dtf1, dtf2, by="Lizard", relationship = "many-to-many"), list(
unique(subset(df_onecol, species1 == 'li' & species2 == 'pm')[,c('gene1', 'gene2')] %>% setNames(c('Lizard', 'Lamprey'))),
unique(subset(df_onecol, species1 == 'li' & species2 == 'sh')[,c('gene1', 'gene2')] %>% setNames(c('Lizard', 'Shark'))),
unique(subset(df_onecol, species1 == 'kf' & species2 == 'li')[,c('gene1', 'gene2')] %>% setNames(c('Killifish', 'Lizard'))),
unique(subset(df_onecol, species1 == 'ca' & species2 == 'li')[,c('gene1', 'gene2')] %>% setNames(c('Goldfish', 'Lizard'))),
unique(subset(df_onecol, species1 == 'li' & species2 == 'ze')[,c('gene1', 'gene2')] %>% setNames(c('Lizard', 'Zebrafish'))),
unique(subset(df_onecol, species1 == 'am' & species2 == 'li')[,c('gene1', 'gene2')] %>% setNames(c('Axolotl', 'Lizard'))),
unique(subset(df_onecol, species1 == 'li' & species2 == 'ne')[,c('gene1', 'gene2')] %>% setNames(c('Lizard', 'Newt'))),
unique(subset(df_onecol, species1 == 'ch' & species2 == 'li')[,c('gene1', 'gene2')] %>% setNames(c('Chicken', 'Lizard')))
)) %>% unique
# Found in chicken, but not in lizard? or vice versa
# Retreive genes missed from chicken projection
missed.by.ch = subset(marker.df.li, Lizard %in% setdiff(lizard.symbols, marker.df.ch$Lizard))
# Retrieve genes missed from lizard projection
missed.by.li = subset(marker.df.ch, Chicken %in% setdiff(chicken.symbols, marker.df.li$Chicken))
marker.df1 = rbind(marker.df.li, missed.by.li)
marker.df2 = rbind(marker.df.ch, missed.by.ch)
if(nrow(marker.df1) > nrow(marker.df2)){
message('Using lizard as primary reference...')
marker.df = marker.df1
} else {
message('Using chicken as primary reference...')
marker.df = marker.df2
}
# Add annotation
marker.df$row_annotation = as.character(i)
marker.df$fraction = apply(marker.df, 1, function(x) length(which(!is.na(x[1:8]))) / length(this.species))
if(nrow(marker.df) == 0) return(NULL) else as.data.frame(unique(marker.df))
})
# Check redundancy factor
lapply(marker.dfs, function(x) nrow(x) / length(unique(x[,1])))
# Mostly good, only cluster 1 (NNgly) has high value
# Which clusters did we not find genes for?
missing.clusters = which(sapply(marker.dfs, is.null))
# nmListModified = nmList[NM_SPECIES]
markers.df = Reduce(gtools::smartbind, marker.dfs[-missing.clusters])[,c(NM_SPECIES, 'row_annotation', 'fraction')] #lapply(marker.dfs, function(x) x[,NM_SPECIES]))
js.dat = readRDS('../../Ortho_Objects/js.dat.rds')
markers.df.sorted = markers.df %>% arrange(factor(row_annotation, levels = (levels(js.dat$data$col))))
dir.create('samap/genes', showWarnings = F)
saveRDS(markers.df.sorted, 'samap/genes/markers.df.sorted.20251230.rds')
# pdf('../../figures/my_figs/ac-fish-types-heatmap-v6.pdf', height=5, width=10)
# SAMapHeatmap4(nmListModified,
# markers.df.sorted[,c(NM_SPECIES, 'row_annotation')],
# # type_palette = colorspace::lighten(c('SAC' = 'orangered', PDGFRA = 'chartreuse', A2 = 'deepskyblue', VG3 = 'pink'), 0.2),
# species_palette = species_palette2,
# species.use = NM_SPECIES, #c("Chicken", "Lizard", "Zebrafish", "Killifish"), # names(species_ids),
# types.use = unlist(ac.groups.sorted),
# show_heatmap_legend = FALSE,
# show_row_names = FALSE,
# show_column_names = FALSE,
# col_fun = circlize::colorRamp2(c(-1, 0, 2), c("white", "white", "#584B9FFF")),
# types.order = unlist(ac.groups.sorted),
# # min.z.score = -1,
# # max.z.score = 2,
# # width = unit(10, "in"),
# rotate = TRUE,
# dotplot = FALSE,
# use_raster = TRUE,
# add.breaks = FALSE)
# dev.off()
```
## Heatmap 2: controlled for multiplicity
```{r}
markers.df.sorted = readRDS('samap/genes/markers.df.sorted.20251230.rds')
markers.df.sorted.trimmed = ControlMultiplicity(markers.df.sorted, multiplicity = 2, sort = TRUE)
bbh = readRDS('../../Ortho_Objects/bbh_v2.rds')
markers.df.sorted.trimmed.filtered = subset(markers.df.sorted.trimmed, row_annotation %in% bbh$ident2)
ha.metadata = data.frame(
name = rep(names(ac.groups.sorted), lengths(ac.groups.sorted)),
value = unlist(ac.groups.sorted)
)
ha.metadata$species = ExtractString(ha.metadata$value, after = '_')
ha.metadata$type = as.character(bbh$ident1[match(ha.metadata$name, bbh$ident2)])
pdf('../../figures/my_figs/ac-fish-types-heatmap-v6.pdf', height=5, width=10)
SAMapHeatmap4(nmList,
markers.df.sorted.trimmed.filtered, # subset(markers.df.sorted, row_annotation %in% c(5))[1:20,],
# type_palette = colorspace::lighten(c('SAC' = 'orangered', PDGFRA = 'chartreuse', A2 = 'deepskyblue', VG3 = 'pink'), 0.2),
species_palette = species_palette3,
species.use = NM_SPECIES, #c("Chicken", "Lizard", "Zebrafish", "Killifish"), #names(species_ids),
types.use = unlist(ac.groups.sorted),
show_heatmap_legend = TRUE,
show_row_names = FALSE,
show_column_names = FALSE,
col_fun = circlize::colorRamp2(c(-1, 0, 2), c("white", "white", "#584B9FFF")),
types.order = (unlist(ac.groups.sorted)),
show_annotation_legend = TRUE,
# min.z.score = -1,
# max.z.score = 2,
# width = unit(10, "in"),
rotate = TRUE,
dotplot = FALSE,
use_raster = TRUE,
raster_quality = 5,
add.breaks = FALSE,
na_col = "grey90",
mybreaks = ha.metadata$name,
left_annotation = rowAnnotation(type = ha.metadata$type,
species = factor(ha.metadata$species, levels = names(prefix.key)),
col = list(type = type_cols3, species = nm_palette2),
border = TRUE,
show_legend = TRUE)
)
dev.off()
unique(ha.metadata$type)
```
## Fig. 5D: SAMap gene combo heatmap
```{r}
ha.metadata = data.frame(
leiden_clusters = rep(names(ac.groups.sorted), lengths(ac.groups.sorted)),
species_type = unlist(ac.groups.sorted)
)
ha.metadata$species = ExtractString(ha.metadata$species_type, after = '_')
ha.metadata$type = as.character(bbh$ident1[match(ha.metadata$leiden_clusters, bbh$ident2)])
write.csv(bbh, 'figures/bbh_v2.csv')
write.csv(ha.metadata, 'figures/ha.metadata.v7.csv')
# write.csv(ha.metadata.clean, 'figures/ha.metadata.clean.v5.csv')
# Now its correct:
# subset(ha.metadata, species_type == 'ch_gabaAC-60')
# leiden_clusters species_type species type
# 48 35 ch_gabaAC-60 ch oAC41
```
```{r}
# ha.metadata.clean = read.csv('figures/ha.metadata.clean.csv')
# ha.metadata.clean = subset(ha.metadata.clean, !leiden_clusters %in% c('17', '22', '11', '20'))
ha.metadata = read.csv('figures/ha.metadata.v7.csv')
ha.metadata.clean = ha.metadata[!is.na(ha.metadata$type),]
ha.metadata.clean = subset(ha.metadata.clean, !type %in% c('oAC27'))
# Order by oAC
custom.order = levels(orthotype_labels(c("22_A2", "27_VG3", "18_TRHDE", "41_TRHDE", "34", "26_A8", "24_SEG", "2_NNgly", "15",
"17_ROBO3", "29_SLC35D3", "35", "30_A17", "32_A17", "9", "1_MAF*", "3*", "20", "4_nNOS",
"14_CA2", "6_NPY", "10_MAF*", "11*", "25*", "13", "40*", "42", "38", "7", "21_PDGFRA",
"28", "5", "36_VIP", "23_RXRG", "39", "16_NNgaba*", "19_CRH*", "12_nNOS*", "33_CA1*", "31_nNOS", "37", "8_SAC*")))
ha.metadata.clean = ha.metadata.clean %>% arrange(factor(type, levels = custom.order))
markers.clean = subset(markers.df.sorted.trimmed, row_annotation %in% ha.metadata.clean$leiden_clusters) %>%
arrange(factor(row_annotation, levels = unique(ha.metadata.clean$leiden_clusters)))
saveRDS(markers.clean, 'samap/genes/markers.clean.v2.rds')
stopifnot(all(ha.metadata.clean$species_type %in% unlist(lapply(nmList, function(x) unique(x$annotated)))))
pdf('../../figures/my_figs/ac-fish-types-heatmap-v3.pdf', height=3.5, width=7.6)
SAMapHeatmap4(nmList,
subset(markers.clean, fraction > 0.3), # subset(markers.df.sorted, row_annotation %in% c(5))[1:20,],
# type_palette = colorspace::lighten(c('SAC' = 'orangered', PDGFRA = 'chartreuse', A2 = 'deepskyblue', VG3 = 'pink'), 0.2),
species_palette = species_palette3,
species.use = NM_SPECIES, #names(species_ids),
types.use = ha.metadata.clean$species_type,
show_heatmap_legend = TRUE,
show_row_names = FALSE,
show_column_names = FALSE,
col_fun = circlize::colorRamp2(c(-1, 0, 2), c("white", "white", "#584B9FFF")),
types.order = ha.metadata.clean$species_type,
show_annotation_legend = TRUE,
na_col = "grey90",
# min.z.score = -1,
# max.z.score = 2,
# width = unit(10, "in"),
rotate = TRUE,
dotplot = FALSE,
use_raster = TRUE,
raster_quality = 5,
add.rectangles = TRUE,
row.breaks = ha.metadata.clean$type,
col.breaks = subset(markers.clean, fraction > 0.3)$row_annotation,
left_annotation = rowAnnotation(type = factor(ha.metadata.clean$type, levels = unique(ha.metadata.clean$type)),
species = ha.metadata.clean$species,
col = list(type = type_cols3, species = nm_palette2),
border = TRUE,
show_legend = TRUE),
lty = 2,
lwd = 0.75
)
dev.off()
unique(ha.metadata.clean$type)
```
## Figure 5E: Custom dotplots for select types
SAC BHLEH22 PDGFRA A8 VG3
Might consider also
37 very nice but not famous
14_CA2 no teleosts
33_CA1 not very clean markers
10_MAF try this next
```{r, fig.height=15, fig.width=16}
# Load nmList
# Load ha.metadata.clean
ha.metadata = read.csv('figures/ha.metadata.v7.csv')
ha.metadata.clean = ha.metadata[!is.na(ha.metadata$type),]
markers.clean = readRDS('samap/genes/markers.clean.v2.rds')
markers.clean.1 = ControlMultiplicity(markers.clean, multiplicity = 1, sort = FALSE)
# lc.use = c(5, 3, 12, 28, 36, 11, 19, 1, 10, 2)
orthotype.use = c('oAC42* [SAC]', 'oAC17*', 'oAC30 [PDGFRA]', 'oAC41',' oAC19 [CA2]', 'oAC39* [CA1]', 'oAC2 [VG3]', 'oAC1 [A2]', 'oAC6 [A8]', 'oAC7 [SEG]')
types.use = subset(ha.metadata.clean, type %in% orthotype.use) %>%
# arrange(factor(leiden_clusters, levels = lc.use)) %>%
arrange(factor(leiden_clusters, levels = orthotype.use)) %>%
pull(species_type)
n_genes = 30
# markers.clean.subset = do.call(rbind, lapply(lc.use, function(lc) head(subset(markers.clean.1, row_annotation == lc), n_genes)))
# pdf('figures/nonmammal-ac-dotplots-v3.pdf', height=30, width=20)
# SAMapHeatmap4(nmList,
# markers.clean.subset,
# # type_palette = colorspace::lighten(c('SAC' = 'orangered', PDGFRA = 'chartreuse', A2 = 'deepskyblue', VG3 = 'pink'), 0.2),
# species_palette = species_palette2,
# species.use = names(nmList), #setdiff(names(nmList), 'Goldfish'),
# types.use = types.use,
# show_heatmap_legend = TRUE,
# show_row_names = TRUE,
# show_column_names = TRUE,
# col_fun = circlize::colorRamp2(c(-1, 0, 2), c("white", "white", "#584B9FFF")),
# types.order = types.use,
# show_annotation_legend = TRUE,
# font.size = 8,
# rotate = FALSE,
# dot.scale.factor = 0.2,
# dotplot = TRUE,
# lty = 2
# )
# dev.off()
```
```{r, fig.width=15, fig.width=15}
markers.clean.1 = ControlMultiplicity(markers.clean, multiplicity = 1, sort = FALSE)
# Additional filtration; doing it manually
# singletons_0.4 = rownames(ac_alignment)[apply(ac_alignment, 1, function(x) max(x) < 0.4)]
# pdf('../../figures/my_figs/nonmammal-integration-matrix-ac-subset.pdf', height=35, width=36)
sm = MakeSmartMatrix(ac_alignment[unique(types.use),unique(gsub('-', '\\.', types.use))])
SmartHeatmap2(sm)
# dev.off()
```
```{r, fig.height=8, fig.width=12}
markers.clean.1 = ControlMultiplicity(markers.clean, multiplicity = 1, sort = FALSE)
orthotype.use = c('oAC42* [SAC]', 'oAC17*', 'oAC30 [PDGFRA]', 'oAC41','oAC2 [VG3]', 'oAC1 [A2]', 'oAC6 [A8]', 'oAC7 [SEG]')
types.use2 = c(
# SAC oAC42* [SAC]
"pm_gabaAC-1", #"pm_gabaAC-2", "pm_gabaAC-3", "pm_gabaAC-4",
"sh_gabaAC-3",
"kf_gabaAC-1",
"ze_gabaAC-0_SAC",
"ne_gabaAC-1", #'ne_gabaAC-14',
"am_gabaAC-1", #'am_gabaAC-17',
"ch_gabaAC-2_SAC", "li_gabaAC-14_SAC", #"li_gabaAC-24_SAC",
# BHLHE22 (oAC17*)
"kf_gabaAC-9",
"ze_glyAC-18", # "ze_gabaAC-39",
# "am_gabaAC-5", 'am_gabaAC-16',
'am_gabaAC-28',
"ch_gabaAC-13", #"ch_gabaAC-57",
'li_gabaAC-5', #"li_gabaAC-0", 'li_gabaAC-2',
# PDGFRA (oAC30 [PDGFRA])
"pm_RGC-5",
# 'sh_gabaAC-8',
"kf_gabaAC-10", #"kf_gabaAC-3",
"ze_gabaAC-50_nNOS",
'ne_gabaAC-23',
"am_gabaAC-13",
"ch_gabaAC-20", #"ch_gabaAC-61",
"li_gabaAC-35",
# 37 (oAC41)
"kf_gabaAC-19",
"ze_gabaAC-37",
'ne_gabaAC-28',
"am_gabaAC-54",
"ch_gabaAC-60",
"li_gabaAC-23",
# 14_CA2
# "kf_gabaAC-19", "ca_gabaAC-2", "ze_gabaAC-37", "am_gabaAC-35", "ch_gabaAC-60", "li_gabaAC-23",
# VG3
"pm_glyAC-1",
"sh_glyAC-12",
"ca_glyAC-3",
"ze_glyAC-9",
'ne_glyAC-7',
"am_glyAC-12",
"ch_glyAC-40_VG3",
"li_glyAC-21_VG3",
# A2
"pm_glyAC-8",
"sh_glyAC-10",
"kf_glyAC-8",
"ze_glyAC-13_A2",
'ne_glyAC-21',
"am_glyAC-3",
"ch_glyAC-1_A2",
"li_glyAC-17_A2",
# A8
"pm_glyAC-5",
"sh_glyAC-4",
"kf_glyAC-7",
"ca_glyAC-7",
"ze_glyAC-16",
# "ze_glyAC-44_SEG",
# 'ne_glyAC-24',
'ne_glyAC-34',
# 'am_glyAC-27',
'am_glyAC-8',
# "ch_glyAC-0",
"ch_glyAC-14_SEG",
# "ch_glyAC-33",
"li_glyAC-19"
# nGnG-GLY
# "pm_RGC-5", "kf_gabaAC-10", "kf_gabaAC-3",
# "ze_gabaAC-50_nNOS", "am_gabaAC-21", "ch_gabaAC-20", "ch_gabaAC-61",
# "li_gabaAC-35",
)
# Get corresponding leiden clusters
# ha.metadata$leiden_clusters[match(orthotype.use, ha.metadata$type)]
markers.clean.subset2 = rbind(subset(markers.clean.1, row_annotation %in% c(1) & (Chicken %in% c('ISL1', 'SLC5A7', 'MEF10', 'MEGF11', 'CHAT') |
Lizard %in% c('ISL1', 'SLC5A7', 'MEF10', 'MEGF11'))),
c('SOX2', 'SOX2', 'SOX2', 'SOX2', 'sox2', 'sox2', NA, 'sox2', 'LOC116954613', 1, 0.75), # SOX2 probe not working in chicken
subset(markers.clean.1, row_annotation %in% c(5) & (Chicken %in% c('BHLHE22', 'BHLHE23', 'CHAT', 'SYT7', 'SEMA6B') |
Lizard %in% c('BHLHE22', 'BHLHE23', 'CHAT'))),
# Recover some genes for oAC30
c('GJC2', 'GJC2', 'GJC1', 'LOC138259626', 'cx47.1', 'gjc1', NA, NA, 'LOC116937600', 21, 0.75),
c('TJP1', 'TJP1', 'TJP1', 'TJP1', 'tjp2b', 'tjp1a', NA, NA, 'LOC116939131', 21, 0.75),
c('PDGFRA', 'PDGFRA', 'PDGFRA', 'PDGFRA', 'pdgfra', NA, NA, NA, 'FLT1', 21, 0.444),
c('CBLN1', NA, 'CBLN1', 'CBLN1', 'cbln1', 'cbln1', 'cbln1', 'cbln2b', 'LOC116954123', 21, 0.889),
# subset(markers.clean.1, row_annotation %in% c(21) & (Chicken %in% c('PDGFRA', 'GJC1', 'GJC2', 'TJP1', 'CBLN1') |
# Lizard %in% c('PDGFRA', 'GJC1', 'GJC2', 'TJP1', 'CBLN1'))),
subset(markers.clean.1, row_annotation %in% c(35) & (Chicken %in% c('FSTL5', 'CNTNAP5', 'NR2F1', 'EPHA5', 'NAV2', 'ELAVL2'))),
# Lizard %in% c('PDGFRA', 'GJC1', 'GJC2', 'TJP1', 'CBLN1'))),
subset(markers.clean.1, row_annotation %in% c(19) & (Chicken %in% c('SLC17A8', 'ERBB4', 'NXPH1', 'BNC2') |
Lizard %in% c('SLC17A8', 'ERBB4', 'NXPH1', 'BNC2'))),
subset(markers.clean.1, row_annotation %in% c(22) & (Chicken %in% c('NFIA', 'PROX1', 'ESRRG') | #'CNTFR'
Lizard %in% c('NFIA', 'PROX1', 'ESRRG'))),
c('PTPRZ1', 'PTPRZ1', 'LOC138497686',

8_figures.Rmd at commit 560dae2, under MIT · at the source

Overview

Authors: Dario Tommasini^1,2, Aboozar Monavarfeshani³, Vishruth Dinesh⁴, Joshua Hahn⁴, Jared Tangeman⁵, Olivier Marre^6,7, Seth Blackshaw⁵, Teresa Puthussery^8,9, Joshua R. Sanes³, Karthik Shekhar^{2,4,9,10,11,12}

ORCID iDs: Dario Tommasini, Aboozar Monavarfeshani, Joshua Hahn, Jared Tangeman, Olivier Marre, Seth Blackshaw, Teresa Puthussery, Joshua R. Sanes, Karthik Shekhar

Department of Neuroscience, University of California, Berkeley, Berkeley, CA, USA
Helen Wills Neuroscience Institute, University of California, Berkeley, Berkeley, CA, USA
Department of Molecular and Cellular Biology and Center for Brain Science, Harvard University, Cambridge, MA, USA
Department of Chemical and Biomolecular Engineering, University of California, Berkeley, Berkeley, CA, USA
Solomon H. Snyder Department of Neuroscience, Johns Hopkins University, Baltimore, MD, USA
Institut de la Vision, Sorbonne Université, INSERM, CNRS, Paris, France
Rothschild Foundation Hospital, Paris, France
Herbert Wertheim School of Optometry and Vision Science, University of California, Berkeley, Berkeley, CA, USA
Vision Sciences Graduate Program, University of California, Berkeley, Berkeley, CA, USA
Center for Computational Biology, University of California, Berkeley, CA, USA
Biophysics Graduate Group, University of California, Berkeley, Berkeley, CA, USA
Biological Systems Division, Lawrence Berkeley National Laboratory, Berkeley, CA, USA

Institutions: University of California, Berkeley (United States); Harvard University (United States); Johns Hopkins University (United States); Centre National de la Recherche Scientifique (France); Institut de la Vision (France); Inserm (France); Sorbonne Université (France); Hôpital Rothschild (France); Fondation de Rothschild (France); Lawrence Berkeley National Laboratory (United States)

Journal: Science advances, volume 12, issue 33, article eaeg3223

Dates: received 9 February 2026; accepted 7 July 2026; published online 14 August 2026; in print August 2026

Type: Research article · Language: English

License: CC BY-NC

Identifiers: DOI 10.1126/sciadv.aeg3223 · PMID 42600003 · PMCID PMC13475493 · OpenAlex W7134841878

Open access: gold, a free copy (OpenAlex)

Preprint: doi:10.64898/2026.03.07.710289

Status: code verified

Categories: genetics / omics (modality), human (organism), cellular / molecular (subfield)

Methods: Connectivity, Statistics, Smoothing, state filtering, decompositions, Machine learning

MeSH: Amacrine Cells*, Biological Evolution*, Retina*, Animals, Evolution, Molecular, Humans, Phylogeny, Retinal Ganglion Cells, Transcriptome, Vertebrates (* major topic)

Journal subjects: Neuroscience, Evolutionary Biology

Topic: Retinal Development and Disorders (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex

Funding: National Institutes of Health (EY028625, EY024265, F31EY038101, EY028633, MH105960, U01NS136405); Research Corporation for Science Advancement (SA-MND-2023-031); Glaucoma Research Foundation; McKnight Foundation; National Science Foundation CRCNS (2309039); BrightFocus Foundation (BFF) (National Glaucoma Award); Howard Hughes Medical Institute Hanna H. Gray; Burroughs Wellcome Fund Postdoctoral Diversity Enrichment Program (1287193); European Research Council DEEPRETINA (101045253)

Citations: not cited yet (Europe PMC); 94 references in the paper

Research resources: RRID:SCR_017852

Abstract

Amacrine cells (ACs) comprise a heterogeneous class of inhibitory neurons in the vertebrate retina, exhibiting morphological and functional complexity rivaling that of cortical interneurons. Here, we integrate single-cell and single-nucleus transcriptomic atlases from 24 vertebrate species to reconstruct the evolutionary origins of this extreme diversity. We identify 42 orthologous AC types, most of which exhibit a one-to-one correspondence across amniotes and, in many cases, across vertebrates. While core molecular identities are conserved, AC types vary in abundance and gene expression across species, likely reflecting adaptations to distinct visual ecologies. AC diversity scales with that of retinal ganglion cells (RGCs), indicative of coevolution. Last, we suggest that ACs arose from an AC-RGC hybrid precursor, with glycinergic ACs diverging early in vertebrate evolution, followed by a bifurcation between RGCs and GABAergic ACs. Together, these findings establish a unified evolutionary framework for understanding the diversity, development, and function of a class of inhibitory neurons across vertebrates.

Reproduced under the paper's license (CC BY-NC), from the paper cited above.

Repositories

Its files are read in the Code ↔ Paper reader above, with 30 matches between paragraphs and lines of code.

shekharlab/AmacrineCells

License: MIT

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Commit: 560dae2585e4d5a6115e5f1a7e03f5add2ae4133, 2 June 2026

Languages: R (20), Python (2), Jupyter (1), Shell (1)

Size: 28 files, 24 scripts

Software Heritage: not archived

Found in: “Data, code, and materials availability:”

Holds: README, license file, 8 notebooks

Not found: CITATION.cff, environment file, tests, continuous integration, documentation

Tools: reshape2 (9 files), tidyverse (6 files), XGBoost (4 files), data.table (3 files), ggplot2 (3 files), SciPy (3 files), anndata (2 files), circlize (2 files), Matplotlib (2 files), NumPy (2 files), pandas (2 files), Scanpy (2 files), scikit-learn (2 files), seaborn (2 files), WGCNA (2 files), clusterProfiler (1 file), ComplexHeatmap (1 file), cowplot (1 file), DESeq2 (1 file), glmnet (1 file), igraph (1 file), patchwork (1 file), Plotly (1 file), psych (1 file), reticulate (1 file), rstatix (1 file), Seurat (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers

27 September 2026: the link answers

26 files

src/0_run.R — R, 385 lines
src/1_species_clustering.Rmd — R, 968 lines, 1 match
src/2_orthotype_analysis.Rmd — R, 4,859 lines, 3 matches
src/3_literature_proportions.Rmd — R, 236 lines
src/4_cluster_reproducibility.Rmd — R, 677 lines
src/5_tf_analysis.Rmd — R, 2,138 lines, 4 matches
src/6_samap_analysis.ipynb — Jupyter, 655 lines
src/7_sac_analysis.Rmd — R, 1,375 lines, 1 match
src/8_figures.Rmd — R, 4,412 lines, 10 matches
src/utils/DarioFunctions.py — Python, 1,027 lines
src/utils/KS_fxns.R — R, 552 lines
src/utils/MappingExample.R — R, 37 lines
src/utils/SamapFunctions.py — Python, 1,034 lines
src/utils/SmartMatrix.R — R, 149 lines
src/utils/dario_functions.R — R, 5,677 lines
src/utils/documented-functions.R/clustering.R — R, 193 lines, 2 matches
src/utils/objects.R — R, 1,187 lines, 4 matches
src/utils/plotCellTypePhylo.R — R, 580 lines
src/utils/plottingFxns.R — R, 581 lines
src/utils/rename.sh — Shell, 27 lines
src/utils/utilFxns.R — R, 457 lines
src/utils/wrappers.R — R, 3,949 lines, 5 matches
src/utils/xgboost_train.R — R, 238 lines
src/utils/xgboost_train_DT.R — R, 393 lines
LICENSE — License, 21 lines
README.md — Text, 45 lines

At the source: github.com/shekharlab/AmacrineCells

Zenodo 20451319

License: MIT

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Size: 1 file

Software Heritage: not checked

Found in: “Data, code, and materials availability:”

Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation

Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)

27 September 2026: the link answers (HTTP 200)

26 files

src/0_run.R — R, 385 lines
src/1_species_clustering.Rmd — R, 968 lines
src/2_orthotype_analysis.Rmd — R, 4,859 lines
src/3_literature_proportions.Rmd — R, 236 lines
src/4_cluster_reproducibility.Rmd — R, 677 lines
src/5_tf_analysis.Rmd — R, 2,138 lines
src/6_samap_analysis.ipynb — Jupyter, 655 lines
src/7_sac_analysis.Rmd — R, 1,375 lines
src/8_figures.Rmd — R, 4,412 lines
src/utils/DarioFunctions.py — Python, 1,027 lines
src/utils/KS_fxns.R — R, 552 lines
src/utils/MappingExample.R — R, 37 lines
src/utils/SamapFunctions.py — Python, 1,034 lines
src/utils/SmartMatrix.R — R, 149 lines
src/utils/dario_functions.R — R, 5,677 lines
src/utils/documented-functions.R/clustering.R — R, 193 lines
src/utils/objects.R — R, 1,187 lines
src/utils/plotCellTypePhylo.R — R, 580 lines
src/utils/plottingFxns.R — R, 581 lines
src/utils/rename.sh — Shell, 27 lines
src/utils/utilFxns.R — R, 457 lines
src/utils/wrappers.R — R, 3,949 lines
src/utils/xgboost_train.R — R, 238 lines
src/utils/xgboost_train_DT.R — R, 393 lines
LICENSE — License, 21 lines
README.md — Text, 45 lines

At the source: doi.org/10.5281/zenodo.20451319

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

2 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
48 scripts, each with its path and the digest of its content;
30 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data, code, and materials availability

All data and code needed to evaluate and reproduce the results in the paper are present in the paper and/or the Supplementary Materials. This study did not generate new materials. The raw and processed sequencing data produced in this work are available via the Gene Expression Omnibus (GEO): GSE332621 (mouse lemur), GSE331097 (rat), GSE335334 (axolotl), and GSE335331 (Iberian ribbed newt). Species reference genomes are described in table S2 and are available on Ensembl (jun2026.archive.ensembl.org) or NCBI (https://ncbi.nlm.nih.gov/datasets/genome/). Code to reproduce the analyses is available on GitHub (https://github.com/shekharlab/AmacrineCells) and Zenodo (https://doi.org/10.5281/zenodo.20451319). Interactive visualization of orthotype atlases is available on the Single Cell Portal (https://singlecell.broadinstitute.org/single_cell): SCP3480 (BC, AC, and RGC orthotypes), SCP3462 (mouse lemur), SCP3446 (rat), SCP3463 (axolotl), and SCP3454 (newt).

Reproduced under the paper's license (CC BY-NC), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 10 authors, 10 MeSH terms, 9 funders, 91 references, 1 RRID.

Cite

This paper

Tommasini, D., Monavarfeshani, A., Dinesh, V., Hahn, J., Tangeman, J., Marre, O., Blackshaw, S., Puthussery, T., Sanes, J. R., & Shekhar, K. (2026). The extreme diversity of retinal amacrine cells has deep evolutionary roots. Science advances, 12(33), eaeg3223. https://doi.org/10.1126/sciadv.aeg3223

BibTeX

@article{tommasini2026extreme, author = {Tommasini, Dario and Monavarfeshani, Aboozar and Dinesh, Vishruth and Hahn, Joshua and Tangeman, Jared and Marre, Olivier and Blackshaw, Seth and Puthussery, Teresa and Sanes, Joshua R. and Shekhar, Karthik}, title = {{The extreme diversity of retinal amacrine cells has deep evolutionary roots}}, journal = {Science advances}, year = {2026}, month = aug, volume = {12}, number = {33}, pages = {eaeg3223}, publisher = {American Association for the Advancement of Science}, issn = {2375-2548}, doi = {10.1126/sciadv.aeg3223}, url = {https://doi.org/10.1126/sciadv.aeg3223}, pmid = {42600003}, pmcid = {PMC13475493} }

RIS

TY - JOUR AU - Tommasini, Dario AU - Monavarfeshani, Aboozar AU - Dinesh, Vishruth AU - Hahn, Joshua AU - Tangeman, Jared AU - Marre, Olivier AU - Blackshaw, Seth AU - Puthussery, Teresa AU - Sanes, Joshua R. AU - Shekhar, Karthik TI - The extreme diversity of retinal amacrine cells has deep evolutionary roots T2 - Science advances J2 - Sci Adv PY - 2026 DA - 2026/08/14 VL - 12 IS - 33 SP - eaeg3223 SN - 2375-2548 PB - American Association for the Advancement of Science DO - 10.1126/sciadv.aeg3223 UR - https://doi.org/10.1126/sciadv.aeg3223 LA - en ER -

CSL-JSON

{ "id": "10.1126/sciadv.aeg3223", "type": "article-journal", "title": "The extreme diversity of retinal amacrine cells has deep evolutionary roots", "container-title": "Science advances", "author": [ { "family": "Tommasini", "given": "Dario" }, { "family": "Monavarfeshani", "given": "Aboozar" }, { "family": "Dinesh", "given": "Vishruth" }, { "family": "Hahn", "given": "Joshua" }, { "family": "Tangeman", "given": "Jared" }, { "family": "Marre", "given": "Olivier" }, { "family": "Blackshaw", "given": "Seth" }, { "family": "Puthussery", "given": "Teresa" }, { "family": "Sanes", "given": "Joshua R." }, { "family": "Shekhar", "given": "Karthik" } ], "container-title-short": "Sci Adv", "volume": "12", "issue": "33", "page": "eaeg3223", "DOI": "10.1126/sciadv.aeg3223", "PMID": "42600003", "PMCID": "PMC13475493", "ISSN": "2375-2548", "publisher": "American Association for the Advancement of Science", "URL": "https://doi.org/10.1126/sciadv.aeg3223", "language": "en", "issued": { "date-parts": [ [ 2026, 8, 14 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41586-026-10629-x [code]: Whole-genome duplication shaped cell-type evolution in the vertebrate brain.
Journal: Nature
In common: reticulate, rstatix, anndata, 17 other tools, genetics / omics, cellular / molecular, 3 references
[2] doi:10.1002/imt2.70163 [code]: Spatial multi-omics unveils sphingolipid metabolic reprogramming within the retinal pathological niche.
Journal: iMeta
In common: rstatix, anndata, DESeq2, 19 other tools, genetics / omics, cellular / molecular, 1 reference
[3] doi:10.1038/s42003-026-10957-8 [code]: Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: glmnet, WGCNA, XGBoost, 19 other tools, cellular / molecular
[4] doi:10.1093/bioinformatics/btag592 [code]: Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: WGCNA, reticulate, rstatix, 18 other tools, genetics / omics, 1 reference
[5] doi:10.1186/s13059-026-04177-w [code]: Genomic sequence evolution underlying human neocortical interareal diversification.
Journal: Genome biology
In common: reticulate, rstatix, anndata, 16 other tools, genetics / omics, cellular / molecular, 1 reference
[6] doi:10.3389/fnmol.2026.1844705 [code]: Risperidone regulates the expression of schizophrenia-related genes in the forebrain of adult male mice.
Journal: Frontiers in molecular neuroscience
In common: psych, anndata, DESeq2, 17 other tools, genetics / omics, cellular / molecular
[7] doi:10.3390/ijms27104466 [code]: Uncovering the Key Circuit FOSL2/FOS/EGR3/EGR1, Contributing to the Hyperexcitability of Excitatory Neurons in the Epileptic Temporal Cortex and Hippocampus.
Journal: International journal of molecular sciences
In common: glmnet, reticulate, anndata, 14 other tools, genetics / omics
[8] doi:10.1038/s41593-026-02300-5 [code]: Integrated single-cell and spatial transcriptomic profiling in ALS uncovers peripheral-to-central immune infiltration and reprogramming.
Journal: Nature neuroscience
In common: anndata, DESeq2, igraph, 16 other tools, genetics / omics, cellular / molecular
[9] doi:10.1038/s41467-026-76675-1 [code]: Long-read proteogenomic atlas of human neuronal differentiation reveals isoform diversity informing neurodevelopmental risk mechanisms.
Journal: Nature communications
In common: WGCNA, reticulate, rstatix, 16 other tools, genetics / omics
[10] doi:10.1038/s41380-026-03629-w [code]: Maternal fasting during early gestation induces epigenetic alterations and schizophrenia-related phenotypes.
Journal: Molecular psychiatry
In common: WGCNA, anndata, DESeq2, 14 other tools, genetics / omics, cellular / molecular, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.