---
title: "Reproducing results discussed in Results-Section-4"
author: Anushka Gupta
date: February 28th 2021
---

# This notebook has commands to recreate results discussed in Results-section-4: Nuclear transcriptome is enriched for long non-coding RNAs that regulate adipogenesis and drive cell-type differences 
```{r setup}
knitr::opts_knit$set(root.dir = '~/Desktop/scRNA-seq_snRNA-seq_adipocyte_lineages/')
```

```{r}
# Load libraries
library(Seurat)
library(Matrix)
library(ggplot2)
library(RColorBrewer)
library(viridis)
library(tibble)
library(clusterCrit)
library(dplyr)
library(ComplexHeatmap)
library("circlize")
library(DropletUtils)
size=20

# loading all genes in ref-2020A. This object also has information regarding gene length and gene-type. 
# the notebook to create this R object can be found in generating_rds folder (features.annotate.ref2020A.Rmd)
load("rdata/features.ref2020A.RData")
gene_name=features.new$V2
gene_length=features.new$length

# For this analysis, cellranger's ref2020-A reference was used as it contains lncRNAs
# Performing DE between white cells and white nuclei using both intronic and exonic reads for UMI quantification
# We separately generated count-matrices for only intronic reads and only exonic reads, and sum the two up to get the final count-matrices
# R Notebooks for generating intronic-only UMI-count matrix can be found in generating_rds folder
# Load intronic only and exonic only matrix for white cells and nuclei
i.wnuclei=Read10X(data.dir = "umi_counts/intronic_white_nuclei_ref-2020-A/")
i.cells=Read10X(data.dir = "umi_counts/intronic_preadipocytes_ref-2020-A/")
# Intronic only matrix for scRNA-seq data has brown preadipocytes as well.
# Therefore, reading in high quality barcodes for white preadipocytes from scRNA-seq dataset
b_cells <- read.table("additional_files/barcodes/singlets.high.quality.white.cells.txt", quote="\"", comment.char="")
b_cells=b_cells$V1
i.cells=i.cells[,intersect(colnames(i.cells),b_cells)]

e.wnuclei=Read10X(data.dir = "umi_counts/exonic_white_nuclei_ref-2020-A/")
e.cells=Read10X(data.dir = "umi_counts/exonic_preadipocytes_ref-2020-A/")
e.wnuclei=e.wnuclei[,colnames(i.wnuclei)]
e.wnuclei=e.wnuclei[gene_name,]
e.cells=e.cells[,colnames(i.cells)]
e.cells=e.cells[gene_name,]

# Normalize and merge intronic and exonic matrices for white cells and white nuclei
pd=0.07
i.wnuclei=i.wnuclei[gene_name,]
i.wnuclei.norm=i.wnuclei/(pd*gene_length/1000)
i.cells=i.cells[gene_name,]
i.cells.norm=i.cells/(pd*gene_length/1000)
wnuclei=i.wnuclei.norm+e.wnuclei
cells=i.cells.norm+e.cells

# Create seurat objects and merge them
s.cells=CreateSeuratObject(cells,project = "cells")
s.nuclei=CreateSeuratObject(wnuclei,project = "nuclei")
s.white=merge(s.cells,s.nuclei,add.cell.ids = c("c", "wn"))
s.white=NormalizeData(s.white,verbose = FALSE)
FindMarkers(object = s.white,features = c("MEG3","NEAT1","MIR31HG","PVT1"),ident.1="cells",ident.2="nuclei",slot = "data",logfc.threshold = .5,verbose = FALSE,only.pos = FALSE)

# Recreating panel (A) in Fig. 4
s.white@meta.data$mir31hg=s.white@assays$RNA@data["MIR31HG",rownames(s.white@meta.data)]
ggplot(s.white@meta.data,aes(x=orig.ident,y=mir31hg,color=orig.ident))+geom_boxplot(size=1.75)+ylab("Normalized Expression Value")+theme_classic(base_size = 40)+theme(legend.position = "none",axis.title.x = element_blank())+scale_colour_manual(values=c("navyblue","goldenrod"))

s.white@meta.data$pvt1=s.white@assays$RNA@data["PVT1",rownames(s.white@meta.data)]
ggplot(s.white@meta.data,aes(x=orig.ident,y=pvt1,color=orig.ident))+geom_boxplot(size=1.75)+ylab("Normalized Expression Value")+theme_classic(base_size = 40)+theme(legend.position = "none",axis.title.x = element_blank())+scale_colour_manual(values=c("navyblue","goldenrod"))+theme(axis.title.y = element_blank())

s.white@meta.data$meg3=s.white@assays$RNA@data["MEG3",rownames(s.white@meta.data)]
ggplot(s.white@meta.data,aes(x=orig.ident,y=meg3,color=orig.ident))+geom_boxplot(size=1.75)+ylab("Normalized Expression Value")+theme_classic(base_size = 40)+theme(legend.position = "none",axis.title.x = element_blank())+scale_colour_manual(values=c("navyblue","goldenrod"))+theme(axis.title.y  = element_blank())

s.white@meta.data$neat1=s.white@assays$RNA@data["NEAT1",rownames(s.white@meta.data)]
ggplot(s.white@meta.data,aes(x=orig.ident,y=neat1,color=orig.ident))+geom_boxplot(size=1.75)+ylab("Normalized Expression Value")+theme_classic(base_size = 40)+theme(legend.position = "none",axis.title.x = element_blank())+scale_colour_manual(values=c("navyblue","goldenrod"))+theme(axis.title.y = element_blank())

# Examining # of lncRNAs DE between white cells and white-nuclei
markers=FindMarkers(object = s.white,ident.1="cells",ident.2="nuclei",slot = "data",logfc.threshold = .5,verbose = FALSE,only.pos = FALSE,features=features.new$V2[features.new$biotype=="lncRNA"])
markers=subset.data.frame(markers,subset = p_val_adj<.05)
markers_cells=markers[markers$avg_logFC>0,]
markers_nuclei=markers[markers$avg_logFC<0,]

save(s.white,file = "rdata/seurat.lncRNA.adipogenic.white.RData")
```


```{r}
# Expression of lncRNAs in brown nuclei over brown cells
# Performing DE between brown cells and brown nuclei using both intronic and exonic reads for UMI quantification
# We separately generated count-matrices for only intronic reads and only exonic reads, and sum the two up to get the final count-matrices
# R Notebooks for generating intronic-only UMI-count matrix can be found in generating_rds folder
# Load intronic only and exonic only matrix for brown cells and nuclei
i.bnuclei=Read10X(data.dir = "umi_counts/intronic_brown_nuclei_ref-2020-A/")
i.cells=Read10X(data.dir = "umi_counts/intronic_preadipocytes_ref-2020-A/")
# Intronic only matrix for scRNA-seq data has white preadipocytes as well.
# Therefore, reading in high quality barcodes for brown preadipocytes from scRNA-seq dataset
b_cells <- read.table("additional_files/barcodes/singlets.high.quality.brown.cells.txt", quote="\"", comment.char="")
b_cells=b_cells$V1
i.cells=i.cells[,intersect(colnames(i.cells),b_cells)]
e.bnuclei=Read10X(data.dir = "umi_counts/exonic_brown_nuclei_ref-2020-A/")
e.cells=Read10X(data.dir = "umi_counts/exonic_preadipocytes_ref-2020-A/")
e.bnuclei=e.bnuclei[,colnames(i.bnuclei)]
e.bnuclei=e.bnuclei[gene_name,]
e.cells=e.cells[,colnames(i.cells)]
e.cells=e.cells[gene_name,]

# Normalize and merge intronic and exonic matrices for brown cells and brown nuclei
pd=0.07
i.bnuclei=i.bnuclei[gene_name,]
i.bnuclei.norm=i.bnuclei/(pd*gene_length/1000)
i.cells=i.cells[gene_name,]
i.cells.norm=i.cells/(pd*gene_length/1000)
bnuclei=i.bnuclei.norm+e.bnuclei
cells=i.cells.norm+e.cells

# Create seurat objects and merge them
s.cells=CreateSeuratObject(cells,project = "cells")
s.nuclei=CreateSeuratObject(bnuclei,project = "nuclei")
s.brown=merge(s.cells,s.nuclei,add.cell.ids = c("c", "bn"))
s.brown=NormalizeData(s.brown,verbose = FALSE)
FindMarkers(object = s.brown,features = c("MEG3","NEAT1","MIR31HG","PVT1"),ident.1="cells",ident.2="nuclei",slot = "data",logfc.threshold = .25,verbose = FALSE,only.pos = FALSE)
```

```{r}
# Recreating panel (B) from Fig. 4
# fastq files were sub-sampled to the following number of transcriptome-mapped reads for both preadipocyte scRNA-seq and snRNA_seq dataset
nreads=5
reads=c(8029,23529, 30647, 45429, 50351)

# Initiatng empty vector for number of lncRNAs detected
median_all_cells=vector(mode = "double",length = nreads)
iqr_all_cells=vector(mode = "list",length = nreads)
median_all_nuclei=vector(mode = "double",length = nreads)
iqr_all_nuclei=vector(mode = "list",length = nreads)

# scRNA-seq and snRNA-seq datasest were aggregated together to have the same number of transcriptome-mapped reads, for varying read depths
# The output of this cellranger aggregation can be found in umi_counts/aggregated_ref2020A folder
# The output is merged UMI-count for preadipocytes from both scRNA-seq and snRNA-seq datasets at varying read depths
dir="umi_counts/aggregated_ref2020A/"
files=list.files(path=dir,full.names = TRUE)

bc_brown <- read.table("additional_files/barcodes/singlets.high.quality.brown.cells.txt", quote="\"", comment.char="")
bc_white <- read.table("additional_files/barcodes/singlets.high.quality.white.cells.txt", quote="\"", comment.char="")
bn_brown <- read.table("additional_files/barcodes/singlets.high.quality.brown.nuclei.txt", quote="\"", comment.char="")
bn_white <- read.table("additional_files/barcodes/singlets.high.quality.white.nuclei.txt", quote="\"", comment.char="")

# After aggergation, cellranger adds a suffic to barcodes based on the order of aggregation
bn_white=substr(bn_white$V1,start = 1,stop = 16)
bn_white=paste(bn_white,2,sep = "-")
bn_brown=substr(bn_brown$V1,start = 1,stop = 16)
bn_brown=paste(bn_brown,3,sep = "-")
bc_white=bc_white$V1
bc_brown=bc_brown$V1

# Initiating a matrix for storing the Mann-Whitney U-test statistic when examining significance in # of lncRNAs detected between scRNA-seq and snRNA-seq at a given read-depth
p.values=matrix(nrow=length(reads),ncol=1)
rownames(p.values)=as.character(reads)
colnames(p.values)="P-Value"

# Only quantify # of lncRNAs detected
genes=features.new$V1[features.new$biotype=="lncRNA"]

# Count # of lncRNAs in scRNA-seq and snRNA-seq as a function of sequencing depth
for (i in 1:nreads) {
  m_aggr=Read10X(data.dir = files[i],gene.column = 1)
  m_aggr=m_aggr[genes,]
  x=colSums(m_aggr[,intersect(colnames(m_aggr),c(bc_white,bc_brown))]>0)
  y=colSums(m_aggr[,intersect(colnames(m_aggr),c(bn_white,bn_brown))]>0)
  median_all_cells[i]=median(x)
  iqr_all_cells[[i]]=quantile(x,probs = c(.25,.75))
  median_all_nuclei[i]=median(y)
  iqr_all_nuclei[[i]]=quantile(y,probs = c(.25,.75))
  x=wilcox.test(x,y)
  p.values[i,1]=x$p.value
}

# Create a tibble for plotting
summary.tibble=tibble(sample=c(replicate(nreads,'cells'),replicate(nreads,'nuclei')),Depth=c(reads,reads),median.umi=c(median_all_cells,median_all_nuclei),quant.25=c(unlist(iqr_all_cells)[seq(from=1,to=9,by=2)],unlist(iqr_all_nuclei)[seq(from=1,to=9,by=2)]),quant.75=c(unlist(iqr_all_cells)[seq(from=2,to=10,by=2)],unlist(iqr_all_nuclei)[seq(from=2,to=10,by=2)]))
formatC(p.values[,1],format="e")
options(scipen = 999)
save(summary.tibble,file = "rdata/summary.tibble.RData")

# Recreating panel (B) from Fig. 4
ggplot(summary.tibble,aes(x=Depth,y=median.umi,color=sample,shape=sample)) +geom_point(size=5)+ geom_line(size=1)+ theme_classic(base_size=20)+scale_color_manual(values=c("navyblue","goldenrod"))+ geom_errorbar(aes(ymin=quant.25, ymax=quant.75,color=sample), width=1e3) + xlab('Transcriptome mapped reads per sample') + ylab('Median lncRNA Detected') + theme(legend.position = "top")
```

```{r}
# Recreating panel (D) from Fig. 4. Panel (C) code chunk is below. 

# Reading high-quality barcodes in brown cluster 1 and cluster 2 in scRNA-seq and snRNA-seq dataset
bc_brown1 <- read.table("additional_files/barcodes/singlets.high.quality.brown.cells.cluster.1.txt", quote="\"", comment.char="")
bc_brown1=bc_brown1$V1
bc_brown2 <- read.table("additional_files/barcodes/singlets.high.quality.brown.cells.cluster.2.txt", quote="\"", comment.char="")
bc_brown2=bc_brown2$V1
bn_brown1 <- read.table("additional_files/barcodes/singlets.high.quality.brown.nuclei.cluster.1.txt", quote="\"", comment.char="")
bn_brown1=substr(bn_brown1$V1,start = 1,stop = 16)
bn_brown1=paste(bn_brown1,3,sep = "-")
bn_brown2 <- read.table("additional_files/barcodes/singlets.high.quality.brown.nuclei.cluster.2.txt", quote="\"", comment.char="")
bn_brown2=substr(bn_brown2$V1,start = 1,stop = 16)
bn_brown2=paste(bn_brown2,3,sep = "-")

# Initiating empty vector to store cluster-separation indices
sc_cells=vector(mode = "double",length = length(nreads))
sc_nuclei=vector(mode = "double",length = length(nreads))
ch_cells=vector(mode = "double",length = length(nreads))
ch_nuclei=vector(mode = "double",length = length(nreads))
db_cells=vector(mode = "double",length = length(nreads))
db_nuclei=vector(mode = "double",length = length(nreads))


for (i in 1:nreads) { 
  print(i)
  m_aggr=Read10X(data.dir = files[i])
  mc_brown=m_aggr[,intersect(colnames(m_aggr),c(bc_brown1,bc_brown2))]
  mn_brown=m_aggr[,intersect(colnames(m_aggr),c(bn_brown1,bn_brown2))]
  
  # Sub-sampling to have the same number of barcodes in scRNA-seq and snRNA-seq
  mn_brown=mn_brown[,sample(x = colnames(mn_brown),size = min(ncol(mc_brown),ncol(mn_brown)))]
  mc_brown=mc_brown[,sample(x = colnames(mc_brown),size = min(ncol(mc_brown),ncol(mn_brown)))]
  
  # Manually adding cluster annotations identified in the analysis of scRNA-seq (fig.1.Rmd) and snRNA-seq (fig.2.Rmd) notebooks
  s.nuclei=CreateSeuratObject(mn_brown)
  s.nuclei@meta.data$cluster=1
  s.nuclei@meta.data[intersect(rownames(s.nuclei@meta.data),bn_brown2),"cluster"]=2
  Idents(s.nuclei)="cluster"
  s.cells=CreateSeuratObject(mc_brown)
  s.cells@meta.data$cluster=1
  s.cells@meta.data[intersect(rownames(s.cells@meta.data),bc_brown2),"cluster"]=2
  Idents(s.cells)="cluster"
  
  # Calculating cluster-separtion indices in PCA space generated using only lncRNAs
  s.nuclei=NormalizeData(s.nuclei,verbose = FALSE)
  s.nuclei=ScaleData(s.nuclei,features = features.new$V2[features.new$biotype=="lncRNA"],verbose = FALSE)
  s.nuclei=RunPCA(s.nuclei,features = features.new$V2[features.new$biotype=="lncRNA"],verbose = FALSE)
  x=intCriteria(s.nuclei@reductions$pca@cell.embeddings[,1:20], as.integer(s.nuclei@meta.data$cluster), crit=c("Silhouette","Calinski_Harabasz","Davies_Bouldin"))
  sc_nuclei[i]=x[[1]]
  ch_nuclei[i]=x[[2]]
  db_nuclei[i]=x[[3]]
  
  s.cells=NormalizeData(s.cells,verbose = FALSE)
  s.cells=ScaleData(s.cells,features = features.new$V2[features.new$biotype=="lncRNA"],verbose = FALSE)
  s.cells=RunPCA(s.cells,features = features.new$V2[features.new$biotype=="lncRNA"],verbose = FALSE)
  x=intCriteria(s.cells@reductions$pca@cell.embeddings[,1:20], as.integer(s.cells@meta.data$cluster), crit=c("Silhouette","Calinski_Harabasz","Davies_Bouldin"))
  sc_cells[i]=x[[1]]
  ch_cells[i]=x[[2]]
  db_cells[i]=x[[3]]
}

# Generating a dataframe for plotting
# Recreating panel (D) from Fig. 4
df=data.frame(reads=reads,type=c(replicate(n = nreads,"cells"),replicate(n = nreads,"nuclei")),silhouette=c(sc_cells,sc_nuclei),Calinski_Harabasz=c(ch_cells,ch_nuclei),Davies_Bouldin=c(db_cells,db_nuclei))
save(df,file = "rdata/separation.coefficients.lncrna.Rdata")

ggplot(df,aes(x=reads,y=silhouette,color=type))+geom_point()+geom_point(size=3)+xlab("Sequencing Depth (x 1e3)") + ylab("Silhouette Coefficient")+theme_classic(base_size = size)+theme(legend.position = "bottom",legend.title = element_blank())+scale_color_manual(values = c("navyblue","goldenrod"))+scale_x_continuous(breaks=seq(from=0,to=50000,by=10000),labels = seq(from=0,to=50,by=10))
```


```{r}
# Recreating panel (C) from Fig. 4
# Panel (C) uses the same sets of barcodes for scRNA-seq and snRNA-seq as in panel (D) above for  ~50,000 reads per cell or nuclei

s.nuclei=CreateSeuratObject(mn_brown)
s.nuclei=NormalizeData(s.nuclei,verbose = FALSE)
s.nuclei@meta.data$cluster=1
s.nuclei@meta.data[intersect(rownames(s.nuclei@meta.data),bn_brown2),"cluster"]=2
Idents(s.nuclei)="cluster"
markers=FindAllMarkers(s.nuclei, only.pos = TRUE,assay = "RNA",slot = "data",logfc.threshold = .25,verbose = FALSE,features=features.new$V2[features.new$biotype=="lncRNA"])
markers=subset.data.frame(markers,subset = p_val_adj<.05)
top <- markers %>% group_by(cluster) %>% top_n(n = 20, wt = avg_logFC)
s.nuclei=ScaleData(s.nuclei,features = top$gene)

# Choosing 100 random barcodes
# This may generate a different looking plot then panel (D) since selection of 100 barcodes is random
y=floor(runif(100, min=1, max=ncol(s.nuclei)))
col_fun = colorRamp2(c(-0.8,0,0,.3), c("darkgrey","ivory", "#BB3754FF","seagreen"))
x=s.nuclei@assays$RNA@scale.data[top$gene,y]
ht=Heatmap(x,column_labels = s.nuclei@meta.data$cluster[y],row_labels = top$gene,show_column_names = TRUE,show_row_names = TRUE,column_km = 2,row_order = top$ensg.id,col=col_fun,column_names_gp = gpar(fontsize = 6))
ht=draw(ht)
```

```{r}
# Performing cluster-separation analysis as a function of UMI 
m_cells=Read10X(data.dir = "umi_counts/intronic_exonic_preadipocytes_ref-2020-A/",gene.column = 1)
mc_brown=m_cells[,c(bc_brown1,bc_brown2)]
mn_brown=Read10X(data.dir = "umi_counts/intronic_exonic_brown_nuclei_ref-2020-A/",gene.column = 1)
mn_brown=mn_brown[,c(bn_brown1,bn_brown2)]

# Sub-sample matrices to have the same number of barcodes in cells and nuclei
mc_brown=mc_brown[,sample(x = colnames(mc_brown),size = ncol(mn_brown))]

# Sub-sample matrices to have mean UMIs stored in depths variable
depths=seq(from=2000,to=23000,by=5000)
# Sub-sample scRNA-seq dataset
fumic_brown=mean(colSums(mc_brown))
dmc_brown=vector(mode = "list",length = length(depths))
for (i in 1:length(depths)) {
  prop=depths[i]/fumic_brown
  dmc_brown[[i]]=downsampleMatrix(x = mc_brown,prop = prop ,bycol = FALSE)
}
# Sub-sample snRNA-seq dataset
fumin_brown=mean(colSums(mn_brown))
dmn_brown=vector(mode = "list",length = length(depths))
for (i in 1:length(depths)) {
  prop=depths[i]/fumin_brown
  dmn_brown[[i]]=downsampleMatrix(x = mn_brown,prop = prop ,bycol = FALSE)
}

# Initiating empty vector to store cluster-separation indices
sc_cells=vector(mode = "double",length = length(depths))
sc_nuclei=vector(mode = "double",length = length(depths))
ch_cells=vector(mode = "double",length = length(depths))
ch_nuclei=vector(mode = "double",length = length(depths))
db_cells=vector(mode = "double",length = length(depths))
db_nuclei=vector(mode = "double",length = length(depths))

# calcultaing separation indices for scRNA-seq
for (i in 1:length(depths)) { 
  print(i)
  s.cells=CreateSeuratObject(dmc_brown[[i]], project = "silhouette")
  s.cells=NormalizeData(s.cells,verbose = FALSE)
  s.cells=ScaleData(s.cells,features = features.new$V1[features.new$biotype=="lncRNA"],verbose = FALSE) 
  s.cells=RunPCA(s.cells,features = features.new$V1[features.new$biotype=="lncRNA"],verbose = FALSE)
  s.cells@meta.data$cluster=1
  s.cells@meta.data[bc_brown2,"cluster"]=2
  x=intCriteria(s.cells@reductions$pca@cell.embeddings[,1:20], as.integer(s.cells@meta.data$cluster), crit=c("Silhouette","Calinski_Harabasz","Davies_Bouldin"))
  sc_cells[i]=x[[1]]
  ch_cells[i]=x[[2]]
  db_cells[i]=x[[3]]
}

# calcultaing separation indices for snRNA-seq
for (i in 1: length(depths)) { 
  print(i)
  s.nuclei=CreateSeuratObject(dmn_brown[[i]], project = "silhouette")
  s.nuclei=NormalizeData(s.nuclei,verbose = FALSE)
  s.nuclei=ScaleData(s.nuclei,features = features.new$V1[features.new$biotype=="lncRNA"],verbose = FALSE)
  s.nuclei=RunPCA(s.nuclei,features = features.new$V1[features.new$biotype=="lncRNA"],verbose = FALSE)
  s.nuclei@meta.data$cluster=1
  s.nuclei@meta.data[intersect(rownames(s.nuclei@meta.data),bn_brown2),"cluster"]=2
  x=intCriteria(s.nuclei@reductions$pca@cell.embeddings[,1:20], as.integer(s.nuclei@meta.data$cluster), crit=c("Silhouette","Calinski_Harabasz","Davies_Bouldin"))
  sc_nuclei[i]=x[[1]]
  ch_nuclei[i]=x[[2]]
  db_nuclei[i]=x[[3]]
}

# Generating a dataframe for plotting
df=data.frame(mean.umi.count=c(depths,depths),type=c(replicate(n = length(depths),"cells"),replicate(n = length(depths),"nuclei")),silhouette=c(sc_cells,sc_nuclei),Calinski_Harabasz=c(ch_cells,ch_nuclei),Davies_Bouldin=c(db_cells,db_nuclei))
```

