---
title: "Reproducing results discussed in Results-Section-1"
author: Anushka Gupta
date: November 2 2021
---

# This notebook has commands to recreate results discussed in Results-section-1: scRNA-seq reveals transcriptional landscape of white and brown preadipocytes 
```{r setup}
knitr::opts_knit$set(root.dir = '~/Desktop/scRNA-seq_snRNA-seq_adipocyte_lineages/')
```

```{r}
# load libraries
library(Seurat)
library(ggplot2)
library(RColorBrewer)
library(viridis)
library(dplyr)
library(limma)
library(org.Hs.eg.db)
library(ComplexHeatmap)
size=20

# Before analysis of scRNA-seq data for Figure 1, we performed demultiplexing, followed by QC. The steps for both analyses can be found in the following notebooks in pre-processing folder:
# Demultiplexing: demultiplexing.htos.Rmd
# QC: QC.scRNAseq.preadipocytes.post.demultiplexing.Rmd

# load seurat object after demultiplexing and QC
# Users can generate this seurat object by using raw FASTQ files, followed by cellranger analysis, demultiplexing, and QC. 
load(file = "rdata/seurat.preadipocytes.post.demultiplexingAndQC.RData")
UMAPPlot(dataset.subset,group.by='cell_type',pt.size=.5)+theme_void()+theme(legend.position = "none")+scale_color_manual(values=c("brown", "grey"))+theme(axis.text.x=element_blank(),axis.ticks.x=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank())
UMAPPlot(dataset.subset,group.by="seurat_clusters",pt.size=.5,cols=brewer.pal(n = 3,name = "Dark2"),label=TRUE,label.size=5)+theme_void()+theme(legend.position = "none")+theme(axis.text.x=element_blank(),axis.ticks.x=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank())

# Expression profile of endothelial (CD31) and hematopoietic marker genes (CD45)
FeaturePlot(dataset.subset,"PECAM1",pt.size=.5)+scale_color_viridis(option = "B")+theme_void(base_size = size)+theme( plot.title = element_text(size = 0))
FeaturePlot(dataset.subset,"PTPRC",pt.size=.5)+scale_color_viridis(option = "B")+theme_void(base_size = size)+theme( plot.title = element_text(size = 0))

# Expression of mesenchymal stem cell markers
FeaturePlot(dataset.subset,"ATXN1",pt.size=.5)+scale_color_viridis(option = "B")+theme_void(base_size = size)+theme( plot.title = element_text(size = 0))
FeaturePlot(dataset.subset,"ITGB1",pt.size=.5)+scale_color_viridis(option = "B")+theme_void(base_size = size)+theme( plot.title = element_text(size = 0))
FeaturePlot(dataset.subset,"THY1",pt.size=.5)+scale_color_viridis(option = "B")+theme_void(base_size = size)+theme( plot.title = element_text(size = 0))
FeaturePlot(dataset.subset,"CD44",pt.size=.5)+scale_color_viridis(option = "B")+theme_void(base_size = size)+theme( plot.title = element_text(size = 0))
FeaturePlot(dataset.subset,"ENG",pt.size=.5)+scale_color_viridis(option = "B")+theme_void(base_size = size)+theme( plot.title = element_text(size = 0))

# Expression of adipogenesis regulators
FeaturePlot(dataset.subset,"CEBPB",pt.size=.5)+scale_color_viridis(option = "B")+theme_void(base_size = size)+theme( plot.title = element_text(size = 0))
FeaturePlot(dataset.subset,"PPARG",pt.size=.5)+scale_color_viridis(option = "B")+theme_void(base_size = size)+theme( plot.title = element_text(size = 0))
FeaturePlot(dataset.subset,"ZEB1",pt.size=.5)+scale_color_viridis(option = "B")+theme_void(base_size = size)+theme( plot.title = element_text(size = 0))


#### Investigating differences between white and brown preadipocytes
#### For this, I am not working with clusters but with manual annotation of cell types 
#### 0 is BAT and 1 is WAT

Idents(dataset.subset)="type"
FindMarkers(object = dataset.subset,features = c("TCF21","PAX3","PDGFRA"),ident.1 = 1,assay = "SCT",slot = "data",logfc.threshold = 0,verbose = FALSE,only.pos = TRUE)

markers.white=FindMarkers(object = dataset.subset,ident.1 = 1,assay = "SCT",slot = "data",logfc.threshold = .5,verbose = FALSE,only.pos = TRUE)
markers.white=subset.data.frame(markers.white,subset = p_val_adj<.05)
markers.brown=FindMarkers(object = dataset.subset,ident.1 = 0,assay = "SCT",slot = "data",logfc.threshold = .5,verbose = FALSE,only.pos = TRUE)
markers.brown=subset.data.frame(markers.brown,subset = p_val_adj<.05)

# Recreating Figure 1 panel (B)
Idents(dataset.subset)="cell_type"
markers=FindAllMarkers(dataset.subset, only.pos = TRUE,assay = "SCT",slot = "data",logfc.threshold = .5,verbose = FALSE)
markers=subset.data.frame(markers,subset = p_val_adj<.05)
top <- markers %>% group_by(cluster) %>% top_n(n = 20, wt = avg_logFC)
DoHeatmap(dataset.subset, features = top$gene,group.colors = c("brown","grey"),raster = FALSE)+scale_fill_viridis(option = "B")+theme(legend.position = "none")

# Expression of myogenic marker genes
FindMarkers(object = dataset.subset,features = c("ACTA2","TAGLN","MYL9","CNN1"),ident.1 = "BAT",assay = "SCT",slot = "data",verbose = FALSE,only.pos = FALSE,logfc.threshold = 0)

#### Classifying genes as markers of white preadipocytes
m_w_cells_df=subset.data.frame(markers.white,subset = avg_logFC>1)
m_w_cells=rownames(m_w_cells_df)[order(m_w_cells_df$pct.2,decreasing = FALSE)]
top_m_w_cells=m_w_cells[1:7]
FeaturePlot(dataset.subset,features = top_m_w_cells)
write.table(x = top_m_w_cells,file = "additional_files/markers/white.cells.txt",sep = "\n",row.names = FALSE,col.names = FALSE,quote = FALSE)

#### Classifying genes as markers of brown preadipocytes
m_b_cells_df=subset.data.frame(markers.brown,subset = avg_logFC>1 & pct.1>.6)
m_b_cells=rownames(m_b_cells_df)[order(m_b_cells_df$pct.2,decreasing = FALSE)]
top_m_b_cells=m_b_cells[1:11]
# Remove marker genes that are differentially expressed between 2 brown clusters
Idents(dataset.subset)="seurat_clusters"
x=rownames(FindMarkers(object = dataset.subset,features = top_m_b_cells,ident.1 = 1,ident.2 = 2,assay = "SCT",slot = "data",logfc.threshold = .5,verbose = FALSE,only.pos = FALSE))
top_m_b_cells=setdiff(top_m_b_cells,x)
Idents(dataset.subset)="type"
FeaturePlot(dataset.subset,features = top_m_b_cells)
write.table(x = top_m_b_cells,file = "additional_files/markers/brown.cells.txt",sep = "\n",row.names = FALSE,col.names = FALSE,quote = FALSE)
```


```{r}
# This code chunk analyzes differences between two brown clusters
Idents(dataset.subset)="seurat_clusters"
markers.1.vs.2.SCT=FindMarkers(object = dataset.subset,ident.1 = 1,ident.2=2,assay = "SCT",slot = "data",logfc.threshold = .5,verbose = FALSE,only.pos = TRUE)
markers.1.vs.2.SCT=subset.data.frame(markers.1.vs.2.SCT,subset = p_val_adj<.05)
markers.2.vs.1.SCT=FindMarkers(object = dataset.subset,ident.1 = 2,ident.2=1,assay = "SCT",slot = "data",logfc.threshold = .5,verbose = FALSE,only.pos = TRUE)
markers.2.vs.1.SCT=subset.data.frame(markers.2.vs.1.SCT,subset = p_val_adj<.05)
# Write these markers and perform GO analysis
write.table(x = rownames(markers.1.vs.2.SCT),file = "additional_files/de_genes/1.vs.2.cells.txt",sep = "\n",row.names = FALSE,col.names = FALSE,quote = FALSE)
write.table(x = rownames(markers.2.vs.1.SCT),file = "additional_files/de_genes/2.vs.1.cells.txt",sep = "\n",row.names = FALSE,col.names = FALSE,quote = FALSE)

# GO analysis was performed using two tools, GO.db and at geneontology.org
# Results from geneontology.org can be found in additional files along with their version use for GO and date when database was accessed
# Showing results from geneontology.org
GO.1.vs.2.cells <- read.delim("additional_files/gene ontology/results/1.vs.2.cells.txt", comment.char="#")
GO.1.vs.2.cells$FDR.Plot=-log10(GO.1.vs.2.cells$upload_1..FDR)
GO.1.vs.2.cells$GO.biological.process.complete=substr(GO.1.vs.2.cells$GO.biological.process.complete,start = 1,stop = nchar(GO.1.vs.2.cells$GO.biological.process.complete)-13)
GO.1.vs.2.cells$GO.biological.process.complete[4]="response to decre. O2"
GO.1.vs.2.cells$GO.biological.process.complete[6]="anat. struc. morphogenesis"
GO.1.vs.2.cells$GO.biological.process.complete[8]="pos. reg. cell migration"
GO.1.vs.2.cells$GO.biological.process.complete[10]="pos. reg. cell motility"
ggplot(GO.1.vs.2.cells[1:10,],mapping = aes(x = reorder(GO.biological.process.complete, FDR.Plot), FDR.Plot))+geom_bar(stat = "identity",fill="#D95F02",color="black")+ coord_flip()+theme_classic(base_size = size)+ylab("-log10(FDR)")+theme(axis.title.y = element_blank())

GO.2.vs.1.cells <- read.delim("additional_files/gene ontology/results/2.vs.1.cells.txt", comment.char="#")
GO.2.vs.1.cells$FDR.Plot=-log10(GO.2.vs.1.cells$upload_1..FDR)
GO.2.vs.1.cells$GO.biological.process.complete=substr(GO.2.vs.1.cells$GO.biological.process.complete,start = 1,stop = nchar(GO.2.vs.1.cells$GO.biological.process.complete)-13)
GO.2.vs.1.cells$GO.biological.process.complete[1]="neg. reg. cellular comp. movement"
GO.2.vs.1.cells$GO.biological.process.complete[10]="neg. reg. multicellular organismal pro."
ggplot(GO.2.vs.1.cells[1:10,],mapping = aes(x = reorder(GO.biological.process.complete, FDR.Plot), FDR.Plot))+geom_bar(stat = "identity",fill="#7570B3",color="black")+ coord_flip()+theme_classic(base_size = size)+ylab("-log10(FDR)")+theme(axis.title.y = element_blank())

# showing results from go.db
entrez.ids <- mapIds(org.Hs.eg.db, keys=rownames(markers.1.vs.2.SCT), column="ENTREZID", keytype="SYMBOL")
entrez.ids.universe=mapIds(org.Hs.eg.db, keys=rownames(dataset.subset), column="ENTREZID", keytype="SYMBOL")
go.out <- goana(unique(entrez.ids), species="Hs",universe=unique(entrez.ids.universe))
go.out <- go.out[order(go.out$P.DE),]
go.useful <- go.out[go.out$Ont=="BP" & go.out$N <= 600,]
head(go.useful, 20)

entrez.ids <- mapIds(org.Hs.eg.db, keys=rownames(markers.2.vs.1.SCT), column="ENTREZID", keytype="SYMBOL")
go.out <- goana(unique(entrez.ids), species="Hs",universe=unique(entrez.ids.universe))
go.out <- go.out[order(go.out$P.DE),]
go.useful <- go.out[go.out$Ont=="BP" & go.out$N <= 600,]
head(go.useful, 20)

# Expression of TGFB superfamily genes and regulators of adipose tissue development
FindMarkers(object = dataset.subset,features = c("BMP4","TGFB2","MYL9","SEMA5A","NPPB","MEST","FST"),ident.1 = 1,ident.2=2,assay = "SCT",slot = "data",verbose = FALSE,only.pos = FALSE)

# TFEA analysis was performed using CheA3 using DE genes in additional_files/de_genes folder as input

# Expression of extracellular matrix (ECM) components 
FindMarkers(object = dataset.subset,features = c("MFAP5","ECM1","COL6A2","ACAN"),ident.1 = 2,ident.2=1,assay = "SCT",slot = "data",verbose = FALSE,only.pos = FALSE)

# Expression of markers identified in Vijay et al. study
FindMarkers(object = dataset.subset,features = c("FBN1","S100A4","IGFBP6","PI16","MFAP5"),ident.1 = 2,ident.2 = 1,assay = "SCT",slot = "data",verbose = FALSE,only.pos = FALSE)

# recreating panel c from Figure 1
Idents(dataset.subset)="cell_type"
markers=FindAllMarkers(dataset.subset, only.pos = TRUE,assay = "SCT",slot = "data",logfc.threshold = .5,verbose = FALSE)
markers=subset.data.frame(markers,subset = p_val_adj<.05)
top <- markers %>% group_by(cluster) %>% top_n(n = 20, wt = avg_logFC)
Idents(dataset.subset)="seurat_clusters"
dataset=subset(dataset.subset,cells=rownames(dataset.subset@meta.data)[dataset.subset@meta.data$seurat_clusters!=0])
markers=FindAllMarkers(dataset, only.pos = TRUE,assay = "SCT",slot = "data",logfc.threshold = .5,verbose = FALSE)
markers=subset.data.frame(markers,subset = p_val_adj<.05)
top.brown <- markers %>% group_by(cluster) %>% top_n(n = 20, wt = avg_logFC)

y=floor(runif(200, min=1, max=ncol(dataset.subset)))
dataset.subset=ScaleData(dataset.subset,assay = "SCT",features = union(top$gene,top.brown$gene))
x=dataset.subset@assays$SCT@scale.data[union(top$gene,top.brown$gene),y]
col_fun = colorRamp2(c(-1,0,2), c("#781B6D","#C33B4E", "#FBFAA2"))
ht=Heatmap(x,column_labels = dataset.subset$seurat_clusters[y],row_labels = union(top$gene,top.brown$gene),show_column_names = TRUE,show_row_names = TRUE,col=col_fun,column_names_gp = gpar(fontsize = 6),row_names_gp = gpar(fontsize = 12,font=3))
ht=draw(ht)

# Classifying genes as markers of brown cluster 1
m_b1_cells_df=subset.data.frame(markers.1.vs.2.SCT,subset = avg_logFC>.8 & pct.1>.6)
top_m_b1_cells=m_b1_cells_df[order(m_b1_cells_df$pct.2,decreasing = FALSE),]
top_m_b1_cells=rownames(top_m_b1_cells)[top_m_b1_cells$pct.2<.4]
write.table(x = top_m_b1_cells,file = "additional_files/markers/1.cells.txt",sep = "\n",row.names = FALSE,col.names = FALSE,quote = FALSE)

# Classifying genes as markers of brown cluster 2
m_b2_cells_df=subset.data.frame(markers.2.vs.1.SCT,subset = avg_logFC>.8 & pct.1>.6)
top_m_b2_cells=m_b2_cells_df[order(m_b2_cells_df$pct.2,decreasing = FALSE),]
top_m_b2_cells=rownames(top_m_b2_cells)[top_m_b2_cells$pct.2<.4]
write.table(x = top_m_b2_cells,file = "additional_files/markers/2.cells.txt",sep = "\n",row.names = FALSE,col.names = FALSE,quote = FALSE)
```

