Set chunk options

Load prepped data

# seur.sub <- readRDS(here("data/fgf_filtered_nuclei.RDS"))
# ```
# 
# # Food Intake and BG of mice
# ```{r FI_BG}
# #Read in excel file of FI
# fi_v <- readxl::read_xlsx(here("data/mouse_data/fig1/191120_FIBG.xlsx"), range = "C7:N24")
# colnames(fi_v) <- c(paste0("V_D", seq_len(6)), paste0("F_D", seq_len(6)))
# melt(fi_v, id.vars = NULL) %>%
#   na.omit() %>%
#   mutate(value = as.numeric(value)) %>%
#   dplyr::group_by(variable) %>%
#   dplyr::summarise(
#     mean = mean(value, na.rm = T),
#     sd = sd(value, na.rm = T),
#     se = sd / sqrt(length(variable))
#   ) %>%
#   mutate(day = rep(c(0:5), 2)) %>%
#   separate(variable, sep = "_", into = "trt") %>%
#   ggplot(aes(x = day, y = mean, color = trt)) + geom_point(size = 1) + 
#   geom_line(size = 1) +
#   geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = .2) + 
#   ggpubr::theme_pubr() +
#   scale_color_manual(name = NULL, labels = c(expression(icv ~ hFGF1 ~ 3 * mu * g), "icv Vehicle pair-fed"), 
#                      values = c("gray30", "gray80")) +
#   scale_x_continuous(breaks = c(0, 1, 2, 3, 4, 5)) +
#   ylab("Daily food intake (g)") + 
#   xlab("Days") + ylim(c(0, 10)) +
#   theme(legend.position = "none", legend.background = element_blank()) + 
#   theme_figure -> fi_fig1
# 
# #Read in excel file of BG
# bg <- readxl::read_xlsx(here("data/mouse_data/fig1/191120_FIBG.xlsx"), range = "B27:O44")
# colnames(bg) <- c(paste0("V_D", seq_len(7)), paste0("F_D", seq_len(7)))
# melt(bg, id.vars = NULL) %>%
#   na.omit() %>%
#   mutate(value = as.numeric(value)) %>%
#   dplyr::group_by(variable) %>%
#   dplyr::summarise(
#     mean = mean(value, na.rm = T),
#     sd = sd(value, na.rm = T),
#     se = sd / sqrt(length(variable))
#   ) %>%
#   mutate(day = rep(c(0:6), 2)) %>%
#   separate(variable, sep = "_", into = "trt") %>%
#   ggplot(aes(x = day, y = mean, color = trt)) + 
#   geom_point(size = 1) + geom_line(size = 1) +
#   geom_errorbar(aes(ymin = mean - se, ymax = mean + se), width = .2) + 
#   ggpubr::theme_pubr() +
#   scale_color_manual(name = NULL, labels = c(expression(icv ~ hFGF1 ~ 3 * mu * g), 
#                                              "icv Vehicle pair-fed"), values = c("gray30", "gray80")) +
#   scale_x_continuous(breaks = c(0, 1, 2, 3, 4, 5, 6),
#     label = c("Pre", "0", "1", "2", "3", "4", "5")) +
#   ylab("BG levels (mg/dl)") + xlab("Days") +
#   theme(legend.direction = "vertical", legend.position = "none", legend.background = element_blank()) + 
#   theme_figure -> bg_fig1
# 
# blank <- plot_grid("")
# bg_fi <- plot_grid(fi_fig1, bg_fig1)
# ```
# 
# # Publication figures of clustering (all nuclei)
# ```{r fig1 panel b}
# seur.sub <- RenameIdents(seur.sub, "Agrp" = "Neuron")
# seur.sub <- RenameIdents(seur.sub, "Hist" = "Neuron")
# seur.sub <- RenameIdents(seur.sub, "Neur" = "Neuron")
# seur.sub <- RenameIdents(seur.sub, "Glia" = "Astro/Tany/Epend")
# seur.sub <- RenameIdents(seur.sub, "COP" = "OPC/COP")
# seur.sub <- RenameIdents(seur.sub, "Endo" = "Endothelial")
# seur.sub <- RenameIdents(seur.sub, "Micro" = "Microglia")
# 
# data.frame(Embeddings(seur.sub, reduction = "umap")) %>% 
#   mutate(group = seur.sub$group) %>%
#   mutate(celltype = Idents(seur.sub)) %>%
#   .[sample(nrow(.)),] %>% 
#   mutate(group = replace(group, group == "FGF_Day-5", "FGF_d5")) %>% 
#   mutate(group = replace(group, group == "FGF_Day-1", "FGF_d1")) %>% 
#   mutate(group = replace(group, group == "PF_Day-1", "Veh_d1")) %>%
#   mutate(group = replace(group, group == "PF_Day-5", "Veh_d5")) -> umap_embed
# 
# label.df <- data.frame(cluster=levels(umap_embed$celltype),label=levels(umap_embed$celltype))
# label.df_2 <- umap_embed %>% 
#   dplyr::group_by(celltype) %>% 
#   dplyr::summarize(x = median(UMAP_1), y = median(UMAP_2)) 
# 
# prop_seur_byclus <- ggplot(umap_embed, aes(x=UMAP_1, y=UMAP_2, color=celltype)) + 
#   geom_point(size=1, alpha=0.5)  + 
#   geom_text_repel(data = label.df_2, aes(label = celltype, x=x, y=y), 
#                   size=3, fontface="bold", inherit.aes = F, bg.colour="white") +
#   xlab("UMAP1") + ylab("UMAP2") +
#   ggpubr::theme_pubr() + ggsci::scale_color_igv() + theme_figure + theme(legend.position = "none", legend.title=NULL)
# 
# prop_seur_group <- ggplot(umap_embed, aes(x=UMAP_1, y=UMAP_2, color=group)) + 
#   geom_point(size=1, alpha=0.5)  + guides(color = guide_legend(override.aes = list(size = 3))) +
#   ggpubr::theme_pubr() + ggsci::scale_color_jco(name = "Treatment Group") +
#   xlab("UMAP1") + ylab("UMAP2") +
#   theme_figure + theme(legend.position = c(.9,.3), legend.title = element_blank(), 
#                        legend.text = element_text(size=12, face="bold"))
# plot_grid(prop_seur_group, prop_seur_byclus)
# ```
# 
# ```{r recluster neurons, fig.width=10, fig.height=4}
# subset(seur.sub, ident="Neuron") %>% 
#   reprocess_subset(., dims = 30, resolution = 0.1) -> fgf.neur
# DimPlot(fgf.neur, label = T)
# ```
# 
# # Load Datasets
# ```{r load mapping datasets}
# load("/projects/mludwig/Dataset_alignment/Data_preprocessing/Campbell_neurons_preprocessed.RData")
# campbell <- UpdateSeuratObject(campbell)
# Idents(campbell) <- "cell_type"
# subset(campbell, idents = c("n34.unassigned(2)", "n33.unassigned(1)"), invert = T) %>%
#   NormalizeData(verbose = FALSE) %>%
#   FindVariableFeatures(
#     selection.method = "vst",
#     nfeatures = 2000, verbose = FALSE
#   ) %>%
#   ScaleData(verbose = FALSE) %>%
#   RunPCA(npcs = 30, verbose = FALSE) -> campbell
# 
# # Load Chen Data
# load("/projects/mludwig/Dataset_alignment/Data_preprocessing/Chen_neurons_preprocessed.RData")
# chen <- UpdateSeuratObject(chen)
# Idents(chen) <- "cell_type"
# chen %>%
#   NormalizeData(verbose = FALSE) %>%
#   FindVariableFeatures(
#     selection.method = "vst",
#     nfeatures = 2000, verbose = FALSE
#   ) %>%
#   ScaleData(verbose = FALSE) %>%
#   RunPCA(npcs = 30, verbose = FALSE) -> chen
# 
# # Load Romanov Data
# load("/projects/mludwig/Dataset_alignment/Data_preprocessing/Romanov_neurons_preprocessed.RData")
# romanov <- UpdateSeuratObject(romanov)
# Idents(romanov) <- "cell_type"
# romanov %>%
#   NormalizeData(verbose = FALSE) %>%
#   FindVariableFeatures(
#     selection.method = "vst",
#     nfeatures = 2000, verbose = FALSE
#   ) %>%
#   ScaleData(verbose = FALSE) %>%
#   RunPCA(npcs = 30, verbose = FALSE) -> romanov
# ```
# 
# #Propagate Labels
# ```{r iterative label propagation, warning=FALSE, message=FALSE}
# #propagate Campbell labels
# fgf.neur <- prop_function(reference = campbell, query = fgf.neur)
# hist(fgf.neur$pred_CONF)
# fgf.neur.camp <- subset(fgf.neur, pred_CONF > 0.75)
# 
# #propagate Chen labels
# subset(fgf.neur, pred_CONF < 0.75) %>%
#   FindVariableFeatures(selection.method = "vst", nfeatures = 2000, verbose = FALSE) %>%
#   prop_function(reference = chen, query = .) -> fgf.neur.relab_chen
# hist(fgf.neur.relab_chen$pred_CONF)
# fgf.neur.chen <- subset(fgf.neur.relab_chen, pred_CONF > 0.75)
# 
# #propagate Romanov labels
# subset(fgf.neur.relab_chen, pred_CONF < 0.75) %>%
#   FindVariableFeatures(selection.method = "vst", nfeatures = 2000, verbose = FALSE) %>%
#   prop_function(reference = romanov, query = .) -> fgf.neur.relab_rom
# hist(fgf.neur.relab_rom$pred_CONF)
# fgf.neur.rom <- subset(fgf.neur.relab_rom, pred_CONF > 0.75)
# 
# #transfer labels to fgf object
# prop_lab <- data.frame(
#   cell.names = c(colnames(fgf.neur.camp), colnames(fgf.neur.chen), colnames(fgf.neur.rom)),
#   labels = c(fgf.neur.camp$pred_ID, fgf.neur.chen$pred_ID, fgf.neur.rom$pred_ID))
# fgf.neur$ref <- as.character(prop_lab[match(colnames(fgf.neur), prop_lab$cell.names), "labels"])
# fgf.neur$ref[is.na(fgf.neur$ref)] <- "unmap"
# 
# #save object with labels
# saveRDS(fgf.neur, file = here("data/neuron/fgf_neur_mappingscores.RDS"))
# ```
# 
# #Filter and recluster dataset
# ```{r reclustering of mapped neurons}
# fgf.neur.prop <- subset(fgf.neur, ref != "unmap")
# fgf.neur.prop <- reprocess_subset(fgf.neur.prop, dims = 30, resolution = 0.1)
# DefaultAssay(fgf.neur.prop) <- "SCT"
# lab.mark <- FindAllMarkers(fgf.neur.prop, only.pos = T, logfc.threshold = 0.5)
# table(fgf.neur.prop@active.ident, fgf.neur.prop$ref) %>%
#   as.data.frame() %>%
#   group_by(Var1) %>%
#   top_n(1, Freq) %>%
#   select(Var1, Var2) -> label_mapping
# lab.mark$prop_label <- as.character(pull(label_mapping[match(lab.mark$cluster, label_mapping$Var1), "Var2"]))
# write_csv(x = lab.mark, here("neuron_clusters.csv"))
# ```
# 
# #Rename clusters
# ```{r find specific markers for both Glu5 populations}
# rename <- c("Hcrt (Chen_Glu12)","Tcf7l2 (Chen_Glu4)","Grm7/Foxb1 (Chen_Glu5)","Nxph1/Foxb1 (Chen_Glu5)",
#             "Tac1/Htr2c (Chen_Glu6)","Hist (Camp_n1)","Tac1/Gad2 (Camp_n2)", 
#             "Avp/Oxt (Camp_n6)","Agrp (Camp_n13)",
#             "Vip (Camp_n16)","Tac2 (Camp_n20)","Sst (Camp_n23)","Hs3st4/Nr5a1 (Camp_n29)",
#             "Fam19a2 (Camp_n31)","Pmch (Roma_Pmch)")
# rename <- rename
# names(rename) <- as.character(label_mapping$Var1)
# ```
# 
# 
# # Neuron specific clustering
# ```{r fig1 panel c}
# fgf.neur.prop[["recluster_0.1"]] <- Idents(object = fgf.neur.prop)
# fgf.neur.prop <- RenameIdents(fgf.neur.prop, rename)
# 
# data.frame(Embeddings(fgf.neur.prop, reduction = "umap")) %>% 
#   mutate(group = fgf.neur.prop$group) %>%
#   mutate(celltype = Idents(fgf.neur.prop)) %>%
#   sample_frac(1L) -> umap_embed
# colnames(umap_embed)[1:2] <- c("UMAP 1", "UMAP 2")
# 
# label.df <- data.frame(cluster=levels(umap_embed$celltype),label=levels(umap_embed$celltype))
# label.df_2 <- umap_embed %>% 
#   group_by(celltype) %>% 
#   dplyr::summarize(x = median(`UMAP 1`), y = median(`UMAP 2`)) 
# 
# prop_neur_byclus <- ggplot(umap_embed, aes(x=`UMAP 1`, y=`UMAP 2`, color=celltype)) + 
#   geom_point(size=0.5, alpha=0.5)  + 
#   geom_text_repel(data = label.df_2, aes(label = celltype, x=x, y=y), 
#                   size=2, 
#                   inherit.aes = F, bg.colour="white", fontface="bold",
#                   force=1, min.segment.length = unit(0, 'lines')) +
#   xlab("UMAP1") + ylab("UMAP2") +
#   ggpubr::theme_pubr(legend="none") + ggsci::scale_color_igv() + theme_figure 
# ```
# 
# #Extract color scheme
# ```{r}
# g <- ggplot_build(prop_neur_byclus)
# cols<-data.frame(colours = as.character(unique(g$data[[1]]$colour)), 
#              label = as.character(unique(g$plot$data[, g$plot$labels$colour])))
# colvec<-as.character(cols$colours)
# names(colvec)<-as.character(cols$label)
# ```
# 
# #Resampling DEG
# ```{r pseudobulk resampling, message=F, warning=FALSE}
# #Generate matrices
# split_mats<-splitbysamp(fgf.neur.prop, split_by="sample")
# names(split_mats)<-unique(Idents(fgf.neur.prop))
# pb<-replicate(100, gen_pseudo_counts(split_mats, ncells=10)) 
# names(pb)<-paste0(rep(names(split_mats)),"_",rep(1:100, each=length(names(split_mats))))
# 
# # Generate DESeq2 Objects
# res<-rundeseq(pb)
# ```
# 
# # Identify neuronal populations with most DE genes at 24 hr
# ```{r plot resampling, message=F, warning=FALSE}
# degenes<-lapply(res, function(x) {
#   tryCatch({
#     y<-x[[2]]
#     y<-na.omit(y)
#     data.frame(y)%>%filter(padj<0.1)%>%nrow()}, 
#     error=function(err) {NA})
# })
# 
# boxplot<-lapply(unique(Idents(fgf.neur.prop)), function(x) {
#   z<-unlist(degenes[grep(as.character(x), names(degenes), fixed = T)])
# })
# 
# names(boxplot)<-unique(Idents(fgf.neur.prop))
# boxplot<-t(as.data.frame(do.call(rbind, boxplot)))
# rownames(boxplot)<-1:100
# genenum<-melt(boxplot)
# write_csv(genenum, path = here("output/neuron/genenum.csv"))

#Figure 1 Panel C

genenum <- read_csv(here("output/neuron/genenum.csv"))
ggplot(genenum,aes(x=reorder(Var2, -value), y=value, fill=factor(Var2))) + 
  geom_boxplot(notch = T, alpha=0.75) + 
  #scale_fill_manual(values = colvec) +
  ggpubr::theme_pubr() +
  theme(axis.text.x = element_text(angle=45, hjust=1), legend.position = "none") + 
  ylab("Number DEG") + xlab(NULL) + theme_figure

#deboxplot

Session information

sessionInfo()

R version 3.5.3 (2019-03-11)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Storage

Matrix products: default
BLAS/LAPACK: /usr/lib64/libopenblas-r0.3.3.so

locale:
 [1] LC_CTYPE=en_DK.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_DK.UTF-8        LC_COLLATE=en_DK.UTF-8    
 [5] LC_MONETARY=en_DK.UTF-8    LC_MESSAGES=en_DK.UTF-8   
 [7] LC_PAPER=en_DK.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_DK.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] ggsignif_0.5.0              gProfileR_0.6.7            
 [3] reshape2_1.4.3              future.apply_1.3.0         
 [5] ggrepel_0.8.0.9000          cowplot_1.0.0              
 [7] parallelDist_0.2.4          cluster_2.1.0              
 [9] future_1.14.0               here_0.1                   
[11] DESeq2_1.22.2               SummarizedExperiment_1.12.0
[13] DelayedArray_0.8.0          BiocParallel_1.16.6        
[15] matrixStats_0.54.0          Biobase_2.42.0             
[17] GenomicRanges_1.34.0        GenomeInfoDb_1.18.2        
[19] IRanges_2.16.0              S4Vectors_0.20.1           
[21] BiocGenerics_0.28.0         forcats_0.4.0              
[23] stringr_1.4.0               dplyr_0.8.3                
[25] purrr_0.3.2                 readr_1.3.1.9000           
[27] tidyr_0.8.3                 tibble_2.1.3               
[29] ggplot2_3.2.1               tidyverse_1.2.1            
[31] Seurat_3.0.3.9036          

loaded via a namespace (and not attached):
  [1] reticulate_1.13        R.utils_2.9.0          tidyselect_0.2.5      
  [4] RSQLite_2.1.1          AnnotationDbi_1.44.0   htmlwidgets_1.3       
  [7] grid_3.5.3             Rtsne_0.15             munsell_0.5.0         
 [10] codetools_0.2-16       ica_1.0-2              withr_2.1.2           
 [13] colorspace_1.4-1       highr_0.8              knitr_1.23            
 [16] rstudioapi_0.10        ROCR_1.0-7             gbRd_0.4-11           
 [19] listenv_0.7.0          Rdpack_0.11-0          git2r_0.25.2          
 [22] GenomeInfoDbData_1.2.0 bit64_0.9-7            rprojroot_1.3-2       
 [25] vctrs_0.2.0            generics_0.0.2         xfun_0.8              
 [28] R6_2.4.0               rsvd_1.0.2             locfit_1.5-9.1        
 [31] bitops_1.0-6           assertthat_0.2.1       SDMTools_1.1-221.1    
 [34] scales_1.0.0           nnet_7.3-12            gtable_0.3.0          
 [37] npsurv_0.4-0           globals_0.12.4         workflowr_1.4.0       
 [40] rlang_0.4.0            zeallot_0.1.0          genefilter_1.64.0     
 [43] splines_3.5.3          lazyeval_0.2.2         acepack_1.4.1         
 [46] broom_0.5.2            checkmate_1.9.4        yaml_2.2.0            
 [49] modelr_0.1.4           backports_1.1.4        Hmisc_4.2-0           
 [52] tools_3.5.3            gplots_3.0.1.1         RColorBrewer_1.1-2    
 [55] ggridges_0.5.1         Rcpp_1.0.2             plyr_1.8.4            
 [58] base64enc_0.1-3        zlibbioc_1.28.0        RCurl_1.95-4.12       
 [61] rpart_4.1-15           pbapply_1.4-1          zoo_1.8-6             
 [64] haven_2.1.0            fs_1.3.1               magrittr_1.5          
 [67] data.table_1.12.2      lmtest_0.9-37          RANN_2.6.1            
 [70] whisker_0.3-2          fitdistrplus_1.0-14    hms_0.5.0             
 [73] lsei_1.2-0             evaluate_0.14          xtable_1.8-4          
 [76] XML_3.98-1.20          readxl_1.3.1           gridExtra_2.3         
 [79] compiler_3.5.3         KernSmooth_2.23-15     crayon_1.3.4          
 [82] R.oo_1.22.0            htmltools_0.3.6        Formula_1.2-3         
 [85] geneplotter_1.60.0     RcppParallel_4.4.3     lubridate_1.7.4       
 [88] DBI_1.0.0              MASS_7.3-51.4          Matrix_1.2-17         
 [91] cli_1.1.0              R.methodsS3_1.7.1      gdata_2.18.0          
 [94] metap_1.1              igraph_1.2.4.1         pkgconfig_2.0.2       
 [97] foreign_0.8-71         plotly_4.9.0           xml2_1.2.0            
[100] annotate_1.60.1        XVector_0.22.0         bibtex_0.4.2          
[103] rvest_0.3.4            digest_0.6.20          sctransform_0.2.0     
[106] RcppAnnoy_0.0.12       tsne_0.1-3             rmarkdown_1.13        
[109] cellranger_1.1.0       leiden_0.3.1           htmlTable_1.13.1      
[112] uwot_0.1.3             gtools_3.8.1           nlme_3.1-140          
[115] jsonlite_1.6           viridisLite_0.3.0      pillar_1.4.2          
[118] lattice_0.20-38        httr_1.4.1             survival_2.44-1.1     
[121] glue_1.3.1             png_0.1-7              bit_1.1-14            
[124] stringi_1.4.3          blob_1.1.1             latticeExtra_0.6-28   
[127] caTools_1.17.1.2       memoise_1.1.0          irlba_2.3.3           
[130] ape_5.3