17a. Plot genes of interest on the trajectory
gene_list <- c("Foxd1", "Ren1", "Cnn1", "Acta2", "Myh11", "Tagln", "Rbpj",
"Mef2a", "Mef2b", "Mef2c", "Mef2d", "Bach1", "Bach2", "Nfix",
"Smarcc1", "Snai2", "Fos", "Jun", "Junb", "Jund", "Mesp1",
"Mesp2", "Tcf4", "Tcf12", "Zeb1", "Stat1",
"Stat3","Tead4","Fosl1", "Foxc1", "Foxc2", "Klf6",
"Nfkb1", "Ebf1", "Akr1b7", "Nfic", "Yy1", "Usf1", "Usf2",
"Rara", "Rxrb", "Rxra", "E4f1", "Smad1", "Wt1", "Fosb",
"Nfya", "Crem", "Creb1", "Pbx1", "Pknox1", "Sp1", "Sp3",
"Hoxa9", "Hoxa10", "Hoxb7", "Srf", "Sdc1", "Gata3", "Nr3c2",
"Noxa1", "Sfrp1", "Kl", "Hnf4a", "Cdh6", "Maf", "Vamp2",
"Sema3f", "Kiss1")
# Klotho = Kl
for (gene in gene_list) {
if (!file.exists(paste0(
sample_name, "_JG-Trajectory-DotPlot_", gene, "_", Sys.Date(), ".svg"))) {
p_JG <- plotTrajectory(proj6, trajectory = "JG",
embedding = "UMAP_LSI_Harmony",
colorBy = "GeneScoreMatrix_no_ribosomal",
name = gene, continuousSet = "horizonExtra")
p2_JG <- plotTrajectory(proj6,
embedding = "UMAP_LSI_Harmony",
trajectory = "JG",
colorBy = "GeneIntegrationMatrix_LSI_Harmony",
name = gene, continuousSet = "blueYellow")
if (!file.exists(paste0(sample_name, "_JG-Trajectory_", gene,
"_", Sys.Date(), ".svg"))) {
svglite(paste0(sample_name, "_JG-Trajectory_",
gene, "_", Sys.Date(), ".svg"))
ggAlignPlots(p_JG[[1]], p2_JG[[1]], type = "h")
dev.off()
}
if (!file.exists(paste0(sample_name, "_JG-Trajectory-DotPlot_", gene,
"_", Sys.Date(), ".svg"))) {
svglite(paste0(sample_name, "_JG-Trajectory-DotPlot_",
gene, "_", Sys.Date(), ".svg"))
ggAlignPlots(p_JG[[2]], p2_JG[[2]], type = "h")
dev.off()
}
}
}
save.image(paste0("sc_kidney_trajectory_plots_", Sys.Date(), ".RData"))
17b. Pairwise testing between groups
Helper function:
getPairwiseComparisons <- function (ArchRProj,
fg_group,
bg_group,
FDR_cutoff = 0.1,
F2C_cutoff = 1.0,
use_clusters,
motif_name) {
markerTest <- getMarkerFeatures(
ArchRProj = ArchRProj,
useMatrix = "PeakMatrix",
groupBy = use_clusters,
testMethod = "wilcoxon",
bias = c("TSSEnrichment", "log10(nFrags)"),
useGroups = fg_group,
bgdGroups = bg_group
)
pma <- plotMarkers(seMarker = markerTest, name = fg_group,
cutOff = paste0("FDR <= ", FDR_cutoff, " & abs(Log2FC) >= ", F2C_cutoff),
plotAs = "MA")
svglite(paste0(sample_name, "_", fg_group, "-vs-", bg_group, "_markerPlot_",
Sys.Date(), ".svg"))
print(pma)
dev.off()
# Motif enrichment in differential peaks
motifsUp <- peakAnnoEnrichment(
seMarker = markerTest,
ArchRProj = ArchRProj,
peakAnnotation = motif_name,
cutOff = paste0("FDR <= ", FDR_cutoff, " & Log2FC >= ", F2C_cutoff)
#cutOff = "FDR <= 0.1 & Log2FC >= 0.5"
)
motifsUp
dfUp <- data.frame(TF = rownames(motifsUp),
mlog10Padj = assay(motifsUp)[,1])
dfUp <- dfUp[order(dfUp$mlog10Padj, decreasing = TRUE),]
dfUp$rank <- seq_len(nrow(dfUp))
fwrite(dfUp,
paste0(sample_name, "_", fg_group, "-vs-", bg_group, "_motifsUp_",
Sys.Date(), ".csv"),
sep=",")
ggUp <- ggplot(dfUp, aes(rank, mlog10Padj, color = mlog10Padj)) +
geom_point(size = 1) +
ggrepel::geom_label_repel(
data = dfUp[rev(seq_len(30)), ],
aes(x = rank, y = mlog10Padj, label = TF),
size = 1.5,
nudge_x = 5,
max.overlaps=30,
color = "black"
) + theme_ArchR() +
ylab("-log10(P-adj) Motif Enrichment") +
xlab("Rank Sorted TFs Enriched") +
scale_color_gradientn(colors = paletteContinuous(set = "comet"))
svglite(paste0(sample_name, "_", fg_group, "-vs-", bg_group, "_motifsUp_",
Sys.Date(), ".svg"))
print(ggUp)
dev.off()
# Down motifs
motifsDo <- peakAnnoEnrichment(
seMarker = markerTest,
ArchRProj = ArchRProj,
peakAnnotation = motif_name,
cutOff = paste0("FDR <= ", FDR_cutoff, " & Log2FC <= -", F2C_cutoff)
#cutOff = "FDR <= 0.1 & Log2FC <= -0.5"
)
dfDo <- data.frame(TF = rownames(motifsDo),
mlog10Padj = assay(motifsDo)[,1])
dfDo <- dfDo[order(dfDo$mlog10Padj, decreasing = TRUE),]
dfDo$rank <- seq_len(nrow(dfDo))
fwrite(dfDo,
paste0(sample_name, "_", fg_group, "-vs-", bg_group, "_motifsDown_",
Sys.Date(), ".csv"),
sep=",")
ggDo <- ggplot(dfDo, aes(rank, mlog10Padj, color = mlog10Padj)) +
geom_point(size = 1) +
ggrepel::geom_label_repel(
data = dfDo[rev(seq_len(30)), ],
aes(x = rank, y = mlog10Padj, label = TF),
size = 1.5,
nudge_x = 5,
max.overlaps = 30,
color = "black"
) + theme_ArchR() +
ylab("-log10(FDR) Motif Enrichment") +
xlab("Rank Sorted TFs Enriched") +
scale_color_gradientn(colors = paletteContinuous(set = "comet"))
svglite(paste0(sample_name, "_", fg_group, "-vs-", bg_group,
"_motifsDown_", Sys.Date(), ".svg"))
print(ggDo)
dev.off()
}
Compare each group along our trajectory to JG cells.
cell_groups <- c(trajectory)
for (i in 1:(length(cell_groups)-1)) {
fg_group <- cell_groups[i]
bg_group <- cell_groups[i+1]
getPairwiseComparisons(proj6,
fg_group = fg_group,
bg_group = bg_group,
use_clusters = "ManualClustersJG",
motif_name = "cisbp_Motif_JG")
}
getPairwiseComparisons(proj6,
fg_group = "late_JGs",
bg_group = "early_JGs",
use_clusters = "ManualClustersJG",
motif_name = "cisbp_Motif_JG")
getPairwiseComparisons(proj6,
fg_group = "late_JGs",
bg_group = "early_VSMCs",
use_clusters = "ManualClustersJG",
motif_name = "cisbp_Motif_JG")
getPairwiseComparisons(proj6,
fg_group = "late_VSMCs",
bg_group = "early_JGs",
use_clusters = "ManualClustersJG",
motif_name = "cisbp_Motif_JG")
getPairwiseComparisons(proj6,
fg_group = "late_JGs",
bg_group = "early_VSMCs",
use_clusters = "ManualClustersJG",
motif_name = "cisbp_Motif_JG")
getPairwiseComparisons(proj6,
fg_group = "late_JGs",
bg_group = "late_VSMCs",
use_clusters = "ManualClustersJG",
motif_name = "cisbp_Motif_JG")
getPairwiseComparisons(proj6,
fg_group = "early_JGs",
bg_group = "early_VSMCs",
use_clusters = "ManualClustersJG",
motif_name = "cisbp_Motif_JG",
F2C_cutoff = 0.5)
getPairwiseComparisons(proj6,
fg_group = "early_VSMCs",
bg_group = "PCs",
use_clusters = "ManualClustersJG",
motif_name = "cisbp_Motif_JG",
F2C_cutoff = 0.5)
getPairwiseComparisons(proj6,
fg_group = "late_JGs",
bg_group = "MCs",
use_clusters = "ManualClustersJG",
motif_name = "cisbp_Motif_JG",
F2C_cutoff = 0.5)
17c. Plot chromVar deviation score plots
# "Gm4881" == "Erfl"
# Cot2 == Nr2f2
gene_list <- c("Foxd1", "Cnn1", "Acta2", "Myh11", "Tagln", "Rbpj",
"Mef2a", "Mef2b", "Mef2c", "Mef2d", "Bach1", "Bach2", "Nfix",
"Smarcc1", "Snai2", "Fos", "Jun", "Junb", "Jund", "Mesp1",
"Mesp2", "Tcf4", "Tcf12", "Zeb1", "Stat1", "Lef1",
"Stat3","Tead4","Fosl1", "Foxc1", "Foxc2", "Klf6",
"Nfkb1", "Ebf1", "Akr1b7", "Nfic", "Yy1", "Usf1", "Usf2",
"Rara", "Rxrb", "Rxra", "Lmx1b", "Gm4881", "E2f1", "Figla",
"Nfya", "Crem", "Creb1", "Pbx1", "Pknox1", "Sp1", "Sp3",
"Hoxa9", "Hoxa10", "Hoxb7", "Nr2f6", "Rfx1", "Rfx2", "Rfx3",
"Zfa", "Zfy1", "Patz1", "Klf2", "Klf5", "Klf15", "Wt1",
"Creb3l3", "Creb3", "Srf", "Maz", "Smad1", "Smad5", "Tcf21",
"Tead4", "Pou2f1", "Hoxc5", "Rela", "Klf7", "Pax5", "Gbx1",
"Atf2", "Nr2f2", "Rfx8", "Rfx4", "Rfx7", "Rfx6", "Rfx5",
"Rbpsuhrs3", "Lcor", "Nfib", "Nfia", "Pgr", "Thrb", "Thra",
"Rarb", "Tcfap2a", "Tcfap2b", "Tcfap2c", "Tcfap2e",
"Tcfap2d", "Arid3b", "Arid3a", "Arid5b", "Arid5a", "Arid3c",
"Arid2", "Setbp1", "Ahctf1", "Hmga2", "Phf21a", "Zfp957",
"Tcfe3", "Mnt", "Npas1", "Tcfap4", "Hic1", "Hic2",
"Sp2", "Foxq1", "Foxn3", "Gm5294", "Foxl1", "Hnf4g", "Zfp637",
"Foxo6", "Foxj3", "Foxn2", "Foxl2", "ENSMUSG00000090020",
"Foxg1", "Foxo4", "Foxd4", "Hoxb9", "Hoxa11", "Cdx2", "Cdx1",
"Cdx4", "Hoxc12", "Hoxc13", "Sox4", "Hmga1", "Hmga1rs1",
"Runx1", "Runx3", "Gata4", "Gata1", "Gata2", "Gata3", "Gata5",
"Hoxd11", "Myog", "Mbd2", "Sp6", "Sp5", "E2f4", "Sp4",
"Zfp161", "Ctcfl", "Mbd1", "Zfp148", "Klf14", "Zfp263", "Egr1",
"Elf1")
positiveTFs <- toupper(gene_list)
positiveTFsMotifs <- getFeatures(proj6,
select = paste(positiveTFs, collapse="|"),
useMatrix = "cisbp_Motif_JGMatrix")
# Get just z scores
positiveTFsMotifs <- grep("z:", positiveTFsMotifs, value = TRUE)
positiveTFsMotifs
# Have to use all clusters present (set order of trajectory first)
p <- plotGroups(ArchRProj = proj6,
groupBy = "ManualClustersJG",
useGroups = trajectory,
groupOrder = c("MCs",
"late_VSMCs",
"late_JGs",
"early_VSMCs",
"early_JGs",
"PCs",
"MCs_and_PCs",
"proliferating_mesenchyme_and_PODs",
"proliferating_MED_stroma_and_FBs",
"proliferating_cortical_stroma_and_EPIs",
"proliferating_stroma_and_PODs",
"MED_and_cortical_stroma_and_FBs",
"PODs",
"PCTs",
"proliferating_CD",
"mesenchymal_progenitors",
"CD_stroma_and_FBs",
"MED_stroma_and_FBs",
"ITSs_and_FBs",
"ureteric_epithelium"),
colorBy = "cisbp_Motif_JGMatrix",
name = positiveTFsMotifs,
imputeWeights = getImputeWeights(proj6)
)
dir.create(file.path(getwd(), "chromVarDeviations"))
for (name in names(p)) {
new_name <- gsub("z:", "", name)
svglite(paste0(file.path(getwd(), "chromVarDeviations"), "/",
new_name, "_chromVarDeviation_", Sys.Date(), ".svg"))
print(p[[name]])
dev.off()
}
save.image(paste0("sc_kidney_chromVarPlots_",Sys.Date(),".RData"))
17d. Generate frequency distribution plots
table(proj6$ManualClustersJG)
CD_principal_cells CD_stroma
195 621
early_JGs early_VSMCs
460 150
glomerular_MCs ICCs_MS_and_FBs
209 238
late_JGs MC_progenitors
267 401
MNPs MNPs_and_glomerular_stroma
670 534
PCs PCs_and_VSMCs
203 98
PCs_VSMCs_and_PODs PODs
87 721
progenitor_PCs proliferating_ICCs
310 681
proliferating_ITSs proliferating_MCs_PCs_and_VSMCs
55 213
proliferating_MS_and_FBs PT
980 49
PT_and_MNPs VSMCs
71 103
library(forcats)
library(ggalluvial)
all_cells[, cluster := proj6$ManualClustersJG]
cell_distribution <- as.data.table(table(all_cells$timepoint,
all_cells$cluster))
colnames(cell_distribution) <- c("timepoint", "cluster", "N")
cell_distribution$timepoint <- factor(cell_distribution$timepoint,
levels = c("E12", "E18", "P5", "P30"))
fwrite(cell_distribution, "proj6_cell_distribution.csv", sep=",")
pal <- paletteDiscrete(values = cell_distribution$cluster)
p <- ggplot(cell_distribution, aes(x = timepoint, y = N, fill = cluster)) +
geom_bar(position = "fill", stat = "identity", color="black") +
scale_fill_manual(values = pal) +
labs(
x = "Mouse kidney developmental timepoint",
y = "Frequency"
) +
custom_theme_ArchR() +
theme(legend.position = "right")
svglite(width=20,
paste0("Cell_distribution_frequency_stacked_barplot_", Sys.Date() ,".svg"))
print(p)
dev.off()
cell_distribution <- cell_distribution[order(timepoint, N),]
cell_pct <- cell_distribution %>%
group_by(timepoint, cluster) %>%
summarise(count = N) %>%
mutate(perc = count/sum(count))
cell_pct$timepoint <- ordered(cell_pct$timepoint,
levels = c("E12", "E18", "P5", "P30"))
connected_alluvial_plot <- ggplot(cell_pct,
aes(y = perc, x = factor(timepoint),
fill = fct_reorder2(cluster, timepoint, perc, .desc = FALSE))) +
geom_alluvium(aes(alluvium = fct_reorder2(cluster, timepoint,
perc, .desc = FALSE)),
alpha=1, color = "black") +
scale_fill_manual(NULL, breaks = 17:0, values=pal)+
scale_y_continuous(expand = c(0,0)) +
custom_theme_ArchR() +
theme(legend.position = "right") +
labs(x = "Mouse kidney developmental timepoint",
y = "Frequency (%)",
fill = "cluster")
svglite(width=20,
paste0("Cell_distribution_alluvial_plot_", Sys.Date() ,".svg"))
print(connected_alluvial_plot)
dev.off()
17e. Plot additional genes of interest along trajectory
gene_list <- c("Sdc1", "Gata3", "Nr3c2", "Noxa1", "Sfrp1", "Kl",
"Hnf4a", "Cdh6", "Maf", "Vamp2", "Sema3f", "Kiss1",
"Mef2a", "Mef2b", "Mef2c", "Mef2d", "Smarcc1", "Nfix",
"Nfic", "Nfib", "Nfia", "Srf", "Stat3", "Pgr", "Bach1", "Bach2",
"Creb3", "Creb3l3", "Creb1", "Rara", "Ebf1", "Elf1", "Elf3",
"Elf5", "Erg", "Rarb", "Klf5", "Wt1", "Zfa", "Zfy1", "Smad1",
"Smad5", "Maz", "E2f1", "Sp1", "Sp3", "Tcfap2d", "Patz1",
"Zfp219", "Zfp148", "Jun", "Junb", "Jund", "Fos", "Crem",
"Esrrg", "Rxra", "Rxrg", "Tcf21", "Snai2", "Usf1", "Usf2",
"Foxd1", "Ren1", "Akr1b7", "Rbpj",
"Acta2", "Smtn", "Tagln", "Myh11", "Cnn1")
for (gene in gene_list) {
if (!file.exists(paste0(
sample_name, "_JG-Trajectory-DotPlot_", gene, "_", Sys.Date(), ".svg"))) {
p_JG <- plotTrajectory(proj6, trajectory = "JG",
embedding = "UMAP_LSI_Harmony",
colorBy = "GeneScoreMatrix_no_ribosomal",
name = gene, continuousSet = "horizonExtra")
p2_JG <- plotTrajectory(proj6,
embedding = "UMAP_LSI_Harmony",
trajectory = "JG",
colorBy = "GeneIntegrationMatrix_LSI_Harmony",
name = gene, continuousSet = "blueYellow")
if (!file.exists(
paste0(sample_name, "_JG-Trajectory_",
gene, "_", Sys.Date(), ".svg"))) {
svglite(paste0(sample_name, "_JG-Trajectory_",
gene, "_", Sys.Date(), ".svg"))
ggAlignPlots(p_JG[[1]], p2_JG[[1]], type = "h")
dev.off()
}
if (!file.exists(
paste0(sample_name, "_JG-Trajectory-DotPlot_",
gene, "_", Sys.Date(), ".svg"))) {
svglite(paste0(sample_name, "_JG-Trajectory-DotPlot_",
gene, "_", Sys.Date(), ".svg"))
ggAlignPlots(p_JG[[2]], p2_JG[[2]], type = "h")
dev.off()
}
}
}