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

# This notebook has commands to recreate results discussed in Results-section-3: Gene length-associated detection bias in single-nuclei RNA-sequencing 
```{r setup}
knitr::opts_knit$set(root.dir = '~/Desktop/scRNA-seq_snRNA-seq_adipocyte_lineages/')
```

```{r}
# Load libraries
library(Seurat)
library(VISION)
library(Matrix)
size=20

# 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
load(file = "rdata/umi.count.white.nuclei.intron.only.RData")
i.wnuclei=mat.intronic.subset
i.cells=Read10X(data.dir = "umi_counts/intronic_preadipocytes/")
colnames(i.cells)=paste(colnames(i.cells),"1",sep = "-")

# 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
white.cells <- read.table("additional_files/barcodes/singlets.high.quality.white.cells.txt", quote="\"", comment.char="", stringsAsFactors=FALSE)
white.cells=white.cells$V1
i.cells=i.cells[,white.cells]

e.wnuclei=Read10X(data.dir = "umi_counts/exonic_white_nuclei/")
e.cells=Read10X(data.dir = "umi_counts/exonic_preadipocytes/")
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,]

# Merge intronic and exonic matrices for nuclei
i.wnuclei=i.wnuclei[gene_name,]
i.wnuclei.norm=i.wnuclei
i.cells=i.cells[gene_name,]
i.cells.norm=i.cells
wnuclei=i.wnuclei.norm+e.wnuclei
cells=i.cells.norm+e.cells

# Create Seurat objects for cells and nuclei and then merge them together
cells=CreateSeuratObject(cells,project = "cells")
wnuclei=CreateSeuratObject(wnuclei,project = "nuclei")
s.merge.inex <- merge(cells, y = wnuclei, add.cell.ids = c("c", "wn"))

#Normalize and perform DE between cells and nuclei
s.merge.inex=NormalizeData(s.merge.inex)
markers.nuclei=FindMarkers(object = s.merge.inex,ident.1 = "nuclei",assay = "RNA",slot = "data",logfc.threshold = 1,verbose = FALSE,only.pos = TRUE)
markers.nuclei=subset.data.frame(markers.nuclei,subset = p_val_adj<.05)
save(markers.nuclei,file = "rdata/mnuclei.inex.RData")

markers.cells=FindMarkers(object = s.merge.inex,ident.1 = "cells",assay = "RNA",slot = "data",logfc.threshold = 1,verbose = FALSE,only.pos = TRUE)
markers.cells=subset.data.frame(markers.cells,subset = p_val_adj<.05)
save(markers.cells,file = "rdata/mcells.inex.RData")

# Extracting ensemble IDs for DE genes
# features.3.0.0.RData contains gene-length information for all the genes in cellranger reference v3.0 (November 2018)
load("rdata/features.3.0.0.RData")
length.cells=vector(mode = "double",length = nrow(markers.cells))
length.nuclei=vector(mode = "double",length = nrow(markers.nuclei))
for (i in 1:nrow(markers.cells)) {
  length.cells[i]=features$length[features$symbol==rownames(markers.cells)[i]]
}
for (i in 1:nrow(markers.nuclei)) {
  length.nuclei[i]=features$length[features$symbol==rownames(markers.nuclei)[i]]
}

# Recreating panel (B) from Fig. 3
median(length.nuclei)/median(length.cells)s
df.length=data.frame(type=c(replicate(n = length(length.nuclei),"nuclei"),replicate(n = length(length.cells),"cells")), length=log(c(length.nuclei/1000,length.cells/1000),base = 10))
library(ggplot2)
ggplot(df.length,aes(x=type,y=length,color=type))+geom_boxplot()+scale_color_manual(values = c("navyblue","goldenrod"))+theme_classic(base_size = 20)+xlab("Genes Enriched In") +ylab("Length (log-scaled Kbp)")
```

```{r}
# Recreating panel (C) and (D) from Fig. 3

load("rdata/seurat.nuclei.preadipocytes.post.QC.RData")
load("rdata/seurat.preadipocytes.post.demultiplexingAndQC.RData")
dataset=s_nuclei
Idents(dataset.subset)="type"

cells=Read10X(data.dir = "~/Desktop/scRNA-seq_snRNA-seq_adipocyte_lineages/umi_counts/exonic_preadipocytes/")
Idents(dataset.subset)="type"
wcells.exonic=cells[,WhichCells(dataset.subset,idents = 1)]
bcells.exonic=cells[,WhichCells(dataset.subset,idents = 0)]
wcells.exonic=CreateSeuratObject(counts = wcells.exonic,project = "1")
bcells.exonic=CreateSeuratObject(counts = bcells.exonic,project = "0")
dataset.subset=merge(wcells.exonic,y = bcells.exonic,add.cell.ids = c("wc","bc"))
dataset.subset=NormalizeData(dataset.subset)

# 1 is WAT in scRNA-seq
#markers.white=FindMarkers(object = dataset.subset,ident.1 = 1,assay = "SCT",slot = "data",logfc.threshold = .5,verbose = FALSE,only.pos = TRUE)
markers.white=FindMarkers(object = dataset.subset,ident.1 = "1",assay = "RNA",slot = "data",logfc.threshold = .5,verbose = FALSE,only.pos = TRUE)
markers.white=subset.data.frame(markers.white,subset = p_val_adj<.05)

wnuclei.exonic=Read10X(data.dir = "~/Desktop/scRNA-seq_snRNA-seq_adipocyte_lineages/umi_counts/exonic_white_nuclei/")
bnuclei.exonic=Read10X(data.dir = "~/Desktop/scRNA-seq_snRNA-seq_adipocyte_lineages/umi_counts/exonic_brown_nuclei/")

#Subset exonic matrix to high quality barcodes
# 0 is WAT in scRNA-seq
Idents(s_nuclei)="batch.indices"
white.nuclei.barcodes=substr(WhichCells(s_nuclei,idents = 0),start = 1,stop = 18)
brown.nuclei.barcodes=substr(WhichCells(s_nuclei,idents = 1),start = 1,stop = 18)
wnuclei.exonic=wnuclei.exonic[,white.nuclei.barcodes]
bnuclei.exonic=bnuclei.exonic[,intersect(brown.nuclei.barcodes,colnames(bnuclei.exonic))]

wnuclei.exonic=CreateSeuratObject(counts = wnuclei.exonic,project = "wnuclei")
bnuclei.exonic=CreateSeuratObject(counts = bnuclei.exonic,project = "bnuclei")
s.nuclei.ex <- merge(wnuclei.exonic, y = bnuclei.exonic, add.cell.ids = c("wn", "bn"))
#Normalize and perform DE between cells and nuclei
s.nuclei.ex=NormalizeData(s.nuclei.ex)

markers.white.nuclei.exonic=FindMarkers(s.nuclei.ex,ident.1 = "wnuclei",assay = "RNA",slot = "data",features = rownames(markers.white),logfc.threshold = 0,only.pos = FALSE,min.pct = -1)
library(Matrix)
mean.umi=rowMeans(s.nuclei.ex@assays$RNA@counts[rownames(markers.white),WhichCells(s.nuclei.ex,idents = "wnuclei")])
library(ggplot2)
markers.white.nuclei.exonic[names(mean.umi),"mean.umi"]=mean.umi
markers.white.nuclei.exonic$label=""
markers.white.nuclei.exonic["ID1","label"]="ID1"
markers.white.nuclei.exonic["BAALC","label"]="BAALC"
#markers.white.nuclei.exonic["LRMDA","label"]="LRMDA"
#markers.white.nuclei.exonic["PGF","label"]="PGF"
markers.white.nuclei.exonic[rownames(markers.white),"avg_logFC_cells"]=markers.white$avg_logFC
#markers.white.nuclei$FDR=as.integer(markers.white.nuclei$p_val_adj<0.05)
ggplot(markers.white.nuclei.exonic,aes(x=log(mean.umi,base = 10),y=avg_logFC,color=avg_logFC_cells,label=label))+geom_point(size=3)+geom_hline(yintercept = 0.5,linetype='dotted', col = 'red',size=2)+theme_classic(base_size = 20)+xlab("log10UMI (nuclei, exonic)")+geom_text(color="red",size=4,aes(fontface=2))+ylab("logFC (nuclei, exonic)")+geom_vline(xintercept = log(markers.white.nuclei.exonic["ID1","mean.umi"],base=10),linetype='dotted', col = 'red',size=2)+scale_color_viridis_c()+labs(color="logFC_cells")
load("/Users/anushka/Desktop/scRNA-seq_snRNA-seq_adipocyte_lineages/rdata/features.3.0.0.RData")
rownames(features)=features$symbol
markers.white.nuclei.exonic$length=features[rownames(markers.white.nuclei.exonic),"length"]

x=markers.white.nuclei.exonic[markers.white.nuclei.exonic$mean.umi<=2.02,]
x=x[x$avg_logFC>0.5,]
markers.white.nuclei.exonic$label=""
markers.white.nuclei.exonic[rownames(x),"label"]=rownames(x)
markers.white.nuclei.exonic["ID1","label"]="ID1"
ggplot(markers.white.nuclei.exonic,aes(x=log(mean.umi,base = 10),y=avg_logFC,color=avg_logFC_cells,label=label))+geom_point(size=3)+geom_hline(yintercept = 0.5,linetype='dotted', col = 'red',size=2)+theme_classic(base_size = 20)+xlab("log10UMI (exonic reads)")+geom_text(color="red",size=4,aes(fontface=2))+ylab("logFC (exonic reads)")+geom_vline(xintercept = log(markers.white.nuclei.exonic["ID1","mean.umi"],base=10),linetype='dotted', col = 'red',size=2)+scale_color_viridis_c()+labs(color="logFC_cells")+geom_vline(xintercept = log(markers.white.nuclei.exonic["KCNMA1","mean.umi"],base=10),linetype='dotted', col = 'blue',size=2)
```

```{r}
mean.umi.inex=rowMeans(s_nuclei@assays$RNA@counts[rownames(markers.white),WhichCells(s_nuclei,idents = 0)])
markers.white.nuclei.exonic[names(mean.umi.inex),"mean.umi.inex"]=mean.umi.inex

markers.white.nuclei.inex=FindMarkers(s_nuclei,ident.1 = 0,assay = "SCT",slot = "data",features = rownames(markers.white),logfc.threshold = 0,only.pos = FALSE,min.pct = -1)
markers.white.nuclei.exonic[rownames(markers.white.nuclei.inex),"avg_logFC_inex"]=markers.white.nuclei.inex$avg_logFC
ggplot(markers.white.nuclei.exonic,aes(x=log(mean.umi.inex,base = 10),y=avg_logFC_inex,color=log(length/1000,base = 10),label=label))+geom_point(size=3)+geom_hline(yintercept = 0.5,linetype='dotted', col = 'red',size=1)+theme_classic(base_size = 20)+xlab("log10UMI")+geom_text(color="red",size=4,aes(fontface=2))+labs(col="length")+scale_color_viridis_c()+xlab("log10UMI (nuclei,inex)")+ylab("logFC (nuclei,inex)")+geom_vline(xintercept = log(mean(s_nuclei@assays$RNA@counts["ID1",WhichCells(s_nuclei,idents = 0)]),base = 10),linetype='dotted', col = 'blue',size=1)

#ggplot(markers.white.nuclei,aes(x=log(mean.umi.nuclei,base = 10),y=avg_logFC,color=log(length/1000,base = 10),label=label))+geom_point(size=3)+geom_hline(yintercept = 0.5,linetype='dotted', col = 'red',size=1)+theme_classic(base_size = 20)+xlab("log10UMI")+geom_text(color="red",size=6,aes(fontface=2))+labs(col="length")+scale_color_viridis_c()+xlab("log10UMI (nuclei)")+ylab("logFC (nuclei)")+geom_vline(xintercept = log(mean(s_nuclei@assays$RNA@counts[gene,WhichCells(s_nuclei,idents = 0)]),base = 10),linetype='dotted', col = 'blue',size=1)
```
```{r}
x=markers.white.nuclei.exonic[markers.white.nuclei.exonic$avg_logFC<0.5,]
x=x[x$avg_logFC_inex>0.5,]
x=x[x$length>50000,]
head(x)

y=markers.white.nuclei.exonic[markers.white.nuclei.exonic$avg_logFC<0.5,]
y=y[y$avg_logFC_inex<0.5,]
y=y[y$length<10000,]
y=y[y$mean.umi<1.2,]
head(y)
```

```{r}
markers.white.nuclei.exonic$label=""
markers.white.nuclei.exonic["AK5","label"]="AK5"
markers.white.nuclei.exonic["BAALC","label"]="BAALC"
markers.white.nuclei.exonic["KCNMA1","label"]="KCNMA1"
markers.white.nuclei.exonic["CRABP2","label"]="CRABP2"
markers.white.nuclei.exonic["TM4SF1","label"]="TM4SF1"
markers.white.nuclei.exonic["EVA1B","label"]="EVA1B"
ggplot(markers.white.nuclei.exonic,aes(x=log(mean.umi,base = 10),y=avg_logFC,label=label))+geom_point(size=3)+geom_hline(yintercept = 0.5,linetype='dotted', col = 'red',size=2)+theme_classic(base_size = 20)+xlab("log10UMI (nuclei, exonic)")+geom_text(color="red",size=4,aes(fontface=2))+ylab("logFC (nuclei, exonic)")+geom_vline(xintercept = log(markers.white.nuclei.exonic["KCNMA1","mean.umi"],base=10),linetype='dotted', col = 'blue',size=2)

ggplot(markers.white.nuclei.exonic,aes(x=log(mean.umi.inex,base = 10),y=avg_logFC_inex,color=log(length/1000,base = 10),label=label))+geom_point(size=3)+geom_hline(yintercept = 0.5,linetype='dotted', col = 'red',size=1)+theme_classic(base_size = 20)+xlab("log10UMI")+geom_text(color="red",size=4,aes(fontface=2))+labs(col="length")+scale_color_viridis_c()+xlab("log10UMI (nuclei,inex)")+ylab("logFC (nuclei,inex)")
```


```{r}
# Recreating panel (E) from Fig. 3

# Performing DE between white cells and white nuclei using only exonic reads for UMI quantification
# Read count matrix for cells and nuclei
cells=Read10X(data.dir = "umi_counts/exonic_preadipocytes/")
wnuclei=Read10X(data.dir = "umi_counts/exonic_white_nuclei/")

# Subset to only have high quality barcodes identified after performing QC
cells=cells[,white.cells]
wnuclei=wnuclei[,intersect(colnames(wnuclei),white.nuclei)]

# Create Seurat objects for cells and nuclei and then merge them together
cells=CreateSeuratObject(cells,project = "cells")
wnuclei=CreateSeuratObject(wnuclei,project = "nuclei")
s.merge.ex <- merge(cells, y = wnuclei, add.cell.ids = c("c", "wn"))

#Normalize and perform DE between cells and nuclei
s.merge.ex=NormalizeData(s.merge.ex)
markers.nuclei=FindMarkers(object = s.merge.ex,ident.1 = "nuclei",assay = "RNA",slot = "data",logfc.threshold = 0.25,verbose = FALSE,only.pos = TRUE)
markers.nuclei=subset.data.frame(markers.nuclei,subset = p_val_adj<.05)

# Find the logFC values for genes that are DE above using only exonic reads, upon inclusion of intronic reads
# Essentially, identify the logFC values for nuclear-enriched genes when only exonic reads are used, and when both exonic and intronic reads are used
x=FindMarkers(object = s.merge,ident.1 = "nuclei",assay = "RNA",slot = "data",logfc.threshold = 0,verbose = FALSE,only.pos = T,features = rownames(markers.nuclei))
x=subset.data.frame(x,subset = p_val_adj<.05)
load("/Users/anushka/Desktop/nuclei_v_cells/RData/features.3.0.0.RData")

# Add gene length information
x$gl=NA
for (i in 1:nrow(x)) {
  x$gl[i]=features$length[features$symbol==rownames(x)[i]]
}

df_inex=data.frame(genes=rownames(x),ex=markers.nuclei[rownames(x),"avg_logFC"],inex=x$avg_logFC,length=x$gl)
save(df_inex,file = "rdata/df.inex.RData")

# Recreating panel (E) from Fig. 3
ggplot(df_inex,aes(x=ex,y=inex))+geom_point(color="black",pch=21,fill="white",size=0.75,alpha=0.2)+coord_fixed()+theme_classic(base_size = size)+xlab("logFC Exonic")+ylab("logFC Intronic + Exonic")+geom_abline(slope = 1,linetype="dashed",color="red",intercept = 0)+xlim(0,max(df_inex$ex))+ylim(0,max(df_inex$ex))
ggplot(df_inex,aes(x=gl/1000,y=inex/ex))+geom_point(size=.5,alpha=0.2)+xlab("Gene length (Kbp)")+ylab("logFC InEx/Ex Reads")+theme_classic(base_size = size)+ylim(0,max(df_inex$inex/df_inex$ex))
```

```{r}
# Recreating panel (F) from Fig. 3
# For normalization, we divide intronic-only count matrix with an estimate of A{15} density
# Normalize and merge intronic and exonic matrices for nuclei
pd=0.07
i.wnuclei.norm=i.wnuclei/(pd*gene_length/1000)
i.cells.norm=i.cells/(pd*gene_length/1000)
wnuclei=i.wnuclei.norm+e.wnuclei
cells=i.cells.norm+e.cells

# Create Seurat objects for cells and nuclei and then merge them together
cells=CreateSeuratObject(cells,project = "cells")
wnuclei=CreateSeuratObject(wnuclei,project = "nuclei")
s.merge.inex.norm <- merge(cells, y = wnuclei, add.cell.ids = c("c", "wn"))
s.merge.inex.norm=NormalizeData(s.merge.inex.norm)

# Find the logFC values for genes that are DE using only exonic reads, upon inclusion of intronic reads, after normalization
# Essentially, identify the logFC values for nuclear-enriched genes when only exonic reads are used, and when both exonic and intronic reads are used, after normalization
y=FindMarkers(object = s.merge.inex.norm,ident.1 = "nuclei",assay = "RNA",slot = "data",logfc.threshold = 0,verbose = FALSE,only.pos = T,features = rownames(markers.nuclei))
y=subset.data.frame(y,subset = p_val_adj<.05)

# Add gene length information
y$gl=NA
for (i in 1:nrow(y)) {
  y$gl[i]=features$length[features$symbol==rownames(y)[i]]
}

df_inex_norm=data.frame(genes=rownames(y),ex=markers.nuclei[rownames(y),"avg_logFC"],inex=y$avg_logFC,length=y$gl)
save(df_inex_norm,file = "rdata/df.inex.norm.RData")
  
# Recreating panel (F) from Fig. 3
ggplot(df_inex_norm,aes(x=ex,y=inex))+geom_point(color="black",pch=21,fill="white",size=0.5,alpha=0.2)+coord_fixed()+theme_classic(base_size = size)+xlab("logFC Exonic")+ylab("logFC Intronic + Exonic")+geom_abline(slope = 1,linetype="dashed",color="red",intercept = 0)+xlim(0,max(df_inex$ex))+ylim(0,max(df_inex$ex))
ggplot(df_inex_norm,aes(x=length/1000,y=inex/ex))+geom_point(size=.5,alpha=0.2)+xlab("Gene length (Kbp)")+ylab("logFC InEx/Ex Reads")+theme_classic(base_size = size)+ylim(0,max(df_inex$inex/df_inex$ex))
```

```{r}
# Normalize and perform DE between cells and nuclei after normalization
markers.nuclei=FindMarkers(object = s.merge.inex.norm,ident.1 = "nuclei",assay = "RNA",slot = "data",logfc.threshold = 1,verbose = FALSE,only.pos = TRUE)
markers.nuclei=subset.data.frame(markers.nuclei,subset = p_val_adj<.05)
markers.cells=FindMarkers(object = s.merge.inex.norm,ident.1 = "cells",assay = "RNA",slot = "data",logfc.threshold = 1,verbose = FALSE,only.pos = TRUE)
markers.cells=subset.data.frame(markers.cells,subset = p_val_adj<.05)

# Extracting ensemble IDs for DE genes
length.cells=vector(mode = "double",length = nrow(markers.cells))
length.nuclei=vector(mode = "double",length = nrow(markers.nuclei))
for (i in 1:nrow(markers.cells)) {
  length.cells[i]=features$length[features$symbol==rownames(markers.cells)[i]]
}
for (i in 1:nrow(markers.nuclei)) {
  length.nuclei[i]=features$length[features$symbol==rownames(markers.nuclei)[i]]
}

median(length.nuclei)/median(length.cells)s
df.length=data.frame(type=c(replicate(n = length(length.nuclei),"nuclei"),replicate(n = length(length.cells),"cells")), length=log(c(length.nuclei/1000,length.cells/1000),base = 10))
ggplot(df.length,aes(x=type,y=length,color=type))+geom_boxplot()+scale_color_manual(values = c("navyblue","goldenrod"))+theme_classic(base_size = 20)+xlab("Genes Enriched In") +ylab("Length (log-scaled Kbp)")
```

```{r}
# Gene ontology results using rownames in markers.cells above
# Showing results from geneontology.org
GO <- read.delim("additional_files/gene ontology/results/white.cells.vs.nuclei.after.norm.txt", comment.char="#")
GO$FDR.Plot=-log10(GO$upload_1..FDR)
GO$GO.biological.process.complete=substr(GO$GO.biological.process.complete,start = 1,stop = nchar(GO$GO.biological.process.complete)-13)
ggplot(GO[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())
```
```{r}
# Perform transcriptional signature analysis using Vision with genes enriched in cells after normalization defind as the signature. This analysis takes time, hence you could directly load the vision object using the command in line # 256
# Read in expression counts (Genes X Cells)
load("rdata/seurat.preadipocytes.post.demultiplexingAndQC.RData")

# Read in expression counts (Genes X Cells)
counts <- dataset.subset@assays$RNA@counts
# Scale counts within a sample
n.umi <- colSums(counts)
scaled_counts <- t(t(counts) / n.umi) * median(n.umi)
# Adding UMAP
projection <- dataset.subset@reductions$umap@cell.embeddings
# Read in meta data (Cells x Vars)
meta = dataset.subset@meta.data
# Read in signatures
# Signatures were defined as the top 100 genes in markers.cells object above after normalization,based on logFC values
sig="~/Desktop/nuclei_v_cells/Signatures/ribosomal_proteins/top100.white.cells.gmt"
vis <- Vision(scaled_counts,signatures = sig,meta = meta)
vis <- addProjection(vis, "UMAP", projection)
vis <- analyze(vis)
save(vis,file = "rdata/vision.preadipocytes.scrnaseq.ribsomal.signature.RData")
# Or directly load Vision object
#load(file = "rdata/vision.preadipocytes.scrnaseq.ribsomal.signature.RData")
viewResults(vis)
```

