---
title: "Single Cell Functional Genomics Analysis in R (clean)"
author: "Matthew Lalli"
date: "4/11/2020"
output: html_document
---


```{r setup, include=FALSE}
rm(list = ls())
knitr::opts_chunk$set(echo = TRUE)
library(Seurat)
library(Matrix)
library(ggthemes)
library(cellrangerRkit)
library(edgeR)
library(ComplexHeatmap)
library(RColorBrewer)
library(tidyr)
library(viridis)
library(dplyr)
library(monocle)
require(gridExtra)
library(readxl)
library(gplots)
library(ggrepel)
library(RRHO)
library(outliers)
library(ggbeeswarm)

library('ggpubr') ### ggboxplot with statistics
 

coul = brewer.pal(12, "Set3")
myPal2 <- magma(10)[c(2,3,4,6,7,9,10)]  ### main color palette throughout 
geneSet <- c("SETD5","PTEN","POGZ","MYT1L","MECP2","HDAC5","DYRK1A","CTNND2","CHD8","CHD2","ASH1L","ARID1B","ADNP")
SFARI_genes <- sort(c('RELN',geneSet))

setwd('/Users/Matthew/Documents/')

## Load custom functions 
source('Lalli_et_al_Rfunctions.R')

todays_history <- paste('Rhistory_',format(Sys.time(), "%b_%d_%Y"),".Rhistory",sep="")

#savehistory(paste(getwd(), todays_history, sep="/"))  ## save history at the end of your session 

##  Load pre-processed data.  To repeat pre-processing from scratch, follow code pre-processing section. 
#SFARI_ALL_with_guide <- readRDS('SFARI_ALL.rds')   ## Seurat output
#SFARI_all_mon_pt3 <- readRDS('SFARI_pseudotime_paper.rds')  ## Monocle output
#timeCourse <- read.table('Timecourse_LUHMES_counts.txt', sep='\t', head=T) # Timecourse bulk RNA-seq data (counts)

```


```{r, Analysis of bulk time course data, include = T}

## Standard bulk-RNA-seq edgeR differential expression analysis 

groups <- as.factor(c('8','8','6','6','4','4','2','2','0','0'))
rownames(timeCourse) = make.names(timeCourse[,3], unique=TRUE)
timeCourseDGE <- DGEList(timeCourse[,c(-1,-2,-3,-4,-5)], group=groups)

x <- timeCourseDGE
## x$samples$group <- relevel(x$samples$group, ref="2") This will set day4 to reference so DE genes will be wrt to this point 
selr1 <- rowSums(cpm(x$counts)>20)>=2  ###  count > 20 in 2 or more samples
x <- x[selr1,]
x <- calcNormFactors(x)

design <- model.matrix(~x$samples$group)
xglm <- estimateGLMCommonDisp(x, design, verbose=T)
fit <- glmFit(xglm, design)
lrt <- glmLRT(fit, coef=5)  ###  DE from day 0 to 5 = genes in neural stem cells and neuron diff 
#lrt <- glmLRT(fit, coef=4)  ###  DE from day 0 to 5 = genes in neural stem cells and neuron diff 

deGenes <- topTags(lrt, n=Inf)
deGenes <- deGenes$table
deSig <- subset(deGenes, deGenes$FDR < 0.00001 & abs(deGenes$logFC) > 1.2)
summary(deSig$logFC > 0)  ##  number of up & down regulated genes each ~ 2000 well balanced up and down regulated genes during differentiation.  

#write.table(lrt$fitted.values[rownames(deSig),c(1,10)], file="Documents/D8_DO_DE.txt")

### correlation heatmap 

M <- cor(x$counts)^2 ## R^2

my_dend <- as.dendrogram(hclust(dist(cor(x$counts))))
my_dend <- reorder(my_dend, 10:1)

Heatmap(M, col=magma(10)[c(-1,-2)], column_dend_reorder = F, cluster_columns = my_dend, cluster_rows = my_dend, row_dend_reorder = F, show_row_dend = F, cell_fun = function(j, i, x, y, width, height, fill) { grid.text(sprintf("%.2f", M[i, j]), x, y, gp = gpar(fontsize = 20)) } )
# or
#Heatmap(M, col=magma(10)[c(-1,-2,-3,-4)], column_dend_reorder = F, cluster_columns = my_dend, cluster_rows = my_dend, row_dend_reorder = F, show_row_dend = F, cell_fun = function(j, i, x, y, width, height, fill) { grid.text(sprintf("%.2f", M[i, j]), x, y, gp = gpar(fontsize = 20)) } )

### Figure 1A -- heatmap of top 1000 variable genes in LUHMES differentiation 

hm_data <- log(cpm(fit$fitted.values[row.names(topTags(lrt,n=1000)),])+1)
colnames(hm_data) <- c('day6','day6','day8','day8','day4','day4','day2','day2','day0','day0')
clusters <- cutree(hclust(dist(hm_data)), 4 )
nofclust.height <-  length(unique(as.vector(clusters)))
selcol <- colorRampPalette(brewer.pal(12,"Set3"))
clustcol.height = selcol(nofclust.height)
labelz <- row.names(hm_data)
label_these <- c('SOX2','CCND2','CCNB1','TP53','MKI67','DCX','MAP2','RELN','SYN1', 'NEUROG2','NEFM', 'NCAM1', 'CDK4')
labelz[-which(row.names(hm_data)%in%label_these)] <- ""
heatmap.2(hm_data, trace="none", dendrogram="column", col= magma, cexRow = 1, density.info = "none", keysize=0.8, key.title='Expression',  labRow=labelz)


plotGeneSetDplyrFacet(geneSet, lrt)  ## Figure 1D


### relevel and repeat above code, changing ceoff to find genes variable on any day of diff 

#day0_genes <- deSig[deSig$logFC > 0,]  ### relevel, 4, coef 2
#day2_genes <- deSig[deSig$logFC > 0,]  ### relevel 4, coef 3
#day4_genes <- deSig[deSig$logFC > 0,]  ### level 0, coef 3
#day6_genes <- deSig[deSig$logFC > 0,]  ### level 0 coef 4
#day8_genes <- deSig[deSig$logFC > 0,]  ### level 0 coef 5

## Temporal genes are used to organize single-cells in to pseudotime 

#temporal_genes <- unique(c(row.names(day0genes), row.names(day2_genes),  row.names(day4_genes),  row.names(day6_genes),  row.names(day8_genes)))

#write.table(temporal_genes, file="Temporal_gene_list_from_bulk.txt", quote=F, row.names=F, col.names=F)

## Plot genes in bulk-rna seq, using dplyr -- Figure 3C
plotGeneSetDplyr(c('MAP2','DCX','TP53','CDK4'), lrt)

## Excitatory genes

ExcitatoryGenes <- c('CPLX1','CAMK2A','WFS1','GRIN2D','CUX2','FOXP1','PROX1','POU3F2','SLC17A7','SLC17A6','SYT1','ZBTB20','NEUROD2','RORB','')
InhibitoryGenesAll <- c('SLC6A1', 'SLC6A12','GAD1','SLC32A1','GAD2')
Dopaminergic <- c('EN','TH','DRD2','NR4A2','SLC18A2','PITX3','DRD1','LMX1B','SLC6A3','LMX1A')

EID_genes <- c(ExcitatoryGenes, InhibitoryGenesAll, Dopaminergic)

my_dend <- rev(as.dendrogram(hclust(dist(t(log(x[EID_genes[EID_genes%in%row.names(x)],]$counts[,1:4]+1))))))

x <- timeCourseDGE  ## unfiltered count data to retain the low or non-expressed genes for E/I/D markers

Heatmap(log(x[EID_genes[EID_genes%in%row.names(x)],]$counts[,1:4]+1), cluster_rows = F, col=magma(10), column_dend_reorder=F, cluster_columns = my_dend, rect_gp = gpar(col = "white", lwd = 2))

```

``` {r, Pre-processing the data, importing into Seurat, include = F }
## download the GEO data and read that in here
SFARI_run1.data <- Read10X(data.dir = "GEO/Run1_filtered_gene_bc_matrices/hg38_transgenes_guides/")
SFARI_run2.data <-  Read10X(data.dir = "GEO/Run2_filtered_gene_bc_matrices/hg38_transgenes_guides/")
SFARI_run3.data <-  Read10X(data.dir = "GEO/Run3_filtered_gene_bc_matrices/hg38_transgenes_guides/")

colnames(SFARI_run1.data) <- paste(colnames(SFARI_run1.data),"-1", sep="") ## some cellbarcodes are shared across runs so add run number to cell name 
colnames(SFARI_run2.data) <- paste(colnames(SFARI_run2.data),"-2", sep="")
colnames(SFARI_run3.data) <- paste(colnames(SFARI_run3.data),"-3", sep="")


#SFARI_run3_transgene.data <- SFARI_run3.data[rownames(SFARI_run3.data)[grep("gene", rownames(SFARI_run3.data), ignore.case = T)],]  ## these are all the custom transgenes, guides / cas9 / gfp etc. Leaving them in doesn't affect global clustering but gRNAs will be found differentially expressed as a sanity check

SFARI <- CreateSeuratObject(SFARI_run1.data, min.cells = 50, project = "SFARI")
SFARI2 <- CreateSeuratObject(SFARI_run2.data, min.cells = 50, project = "SFARI")
SFARI3 <- CreateSeuratObject(SFARI_run3.data, min.cells = 50, project = "SFARI")

bctable <- read.csv('Run1_singlets.csv', head=T)  # These lookup tables are provided in Supplementary Code
names(bctable) <- c("V1","V2")
bctable$V1 <- paste(bctable$V1,"-1",sep="")

bctable2 <- read.csv('Run2_singlets.csv', head=T)
names(bctable2) <- c("V1","V2")
bctable2$V1 <- paste(bctable2$V1,"-2",sep="")

bctable3 <- read.csv('Run3_singlets.csv', head=T)
names(bctable3) <- c("V1","V2")
bctable3$V1 <- paste(bctable3$V1,"-3",sep="")

all_singlets <- rbind(bctable, bctable2[,1:2], bctable3[,1:2])

## Although neurons are post-mitotic, added cell-cycle regression for generalizability 
cc.genes <- readLines('regev_lab_cell_cycle_genes.txt')
s.genes <- cc.genes[1:43]
g2m.genes <- cc.genes[44:97]

```

```{r, Analyzing Seurat perturbseq data: Step 1, add sgRNA to datasets, include = F  } 

SFARI_ALL <- add_sgRNA(SFARI_ALL, all_singlets)

batch <- rep("None", length(names(SFARI_ALL$orig.ident)))
batch[grep("-1", names(SFARI_ALL$orig.ident))] <- "1"
batch[grep("-2", names(SFARI_ALL$orig.ident))] <- "2"
batch[grep("-3", names(SFARI_ALL$orig.ident))] <- "3"
SFARI_ALL$batch <- batch 

# subset of cells that have a gRNA
SFARI_ALL_with_guide <- subset(x = SFARI_ALL, subset = gene_level_numeric > 1)

```

``` {r, Quality Control and Visualizing the Data in Seurat }

SFARI_ALL_with_guide <- mito_qc(SFARI_ALL_with_guide)
VlnPlot(object = SFARI_ALL_with_guide, features = c("nCount_RNA", "nFeature_RNA", "percent.mito"), group.by = "batch")  
SFARI_ALL_with_guide <- QC_Filter_Seurat(SFARI_ALL_with_guide)

# I save this file so I don't need to repeat this step in re-analysis 
#saveRDS(SFARI_ALL_with_guide, file = "~/Documents/SFARI_ALL.rds")

##############################################################################################
### If you have saved the Monocle and Seurat processed R datasets, you can start here for future re-analysis

SFARI_ALL_with_guide <- PCA_TSNE(SFARI_ALL_with_guide)

## Note PC1 is enriched for Neuronal maturation genes (NEFM, STMN2, NEFL, MAP1B, already indicating variability in neuronal maturation in single-cell data

## Fixing typos
SFARI_ALL_with_guide$gene_level <- relevel(factor(SFARI_ALL_with_guide$gene_level), "Nontargeting")
SFARI_ALL_with_guide$sgRNA <- revalue(SFARI_ALL_with_guide$sgRNA, c("PTEN_G2_" = "PTEN_G2"))
SFARI_ALL_with_guide$sgRNA <- revalue(SFARI_ALL_with_guide$sgRNA, c("ARID1B_g1" = "ARID1B_G2"))

### For paper, use SAFRI_ALL_with_guide, UMAP = 6 PCA dims, PCA resolution = 0.25.  Res = 1 to show cluster enrichment, res 0.25 to show few clusters ... Supplemental UMAP Figures  

Idents(SFARI_ALL_with_guide) <- SFARI_ALL_with_guide$RNA_snn_res.0.25

# DimPlot(object = SFARI_ALL_with_guide, reduction = "umap", pt.size = 2, cols=myPal) ## UMAP plot 


#DimPlot(object = SFARI_ALL_with_guide, reduction = "umap", pt.size = 2) + scale_color_economist() + theme_classic(base_size=20)
#DimPlot(object = SFARI_ALL_with_guide, reduction = "umap", pt.size = 2) + scale_color_viridis(option="magma", discrete="T") + theme_classic(base_size=20)

## UMAP state enrichment testing, no sgRNA was really significantly enriched / depleted for broad UMAP states 

a <- sapply(split(SFARI_ALL_with_guide$RNA_snn_res.0.25, SFARI_ALL_with_guide$gene_level), function(x){table(x)/sum(table(x))})
plot(hclust(dist(t(a))))  ### sort the columns to make the figure prettier 
b<- hclust(dist(t(a)), method="ward.D2")
barplot(a[,rev(b$order)], col=myPal)

## Removed ribosomal genes for differential gene expression analysis due to high variability across cells 
ribo_gene.indexes1 <- grep(pattern = "^RPL|^RPS", rownames(SFARI_ALL_with_guide), ignore.case=TRUE) 
SFARI_ALL_with_guide <- SFARI_ALL_with_guide[-ribo_gene.indexes1]

### Visualization of sgRNA efficiency using Seurat DotPlot 

geneSet <- c("SETD5","PTEN","POGZ","MYT1L","MECP2","HDAC5","DYRK1A","CTNND2","CHD8","CHD2","ASH1L","ARID1B","ADNP")
Idents(SFARI_ALL_with_guide) <- SFARI_ALL_with_guide$gene_level
b <- DotPlot(SFARI_ALL_with_guide, features = geneSet, dot.scale=20, scale.by='size') + scale_color_viridis(option="magma")
b$data$id <- factor(b$data$id, rev(sort(levels(b$data$id))))
b$data$id <- relevel(b$data$id, "Nontargeting")
b + coord_flip() + scale_x_discrete(limits = rev(sort(levels(b$data$id)))) + scale_y_discrete(limits = sort(levels(b$data$id)))  ## X = sgRNA, Y = Gene Expression of indicated Gene -- basically a heatmap of candidate genes in cells grouped by sgRNA. Diagonal represents the gene expression of 'On-target' gene.  Size and colors are meaningful within rows. 

```


```{r, Pseudotime Ordering analysis in Monocle, include = T}

#if not running on a server, perform garbage collection to free up space before running Monocle2
gc()

## import Seurat batch corrected dataset into Monocle 
SFARI_all <- newCellDataSet(as.matrix(SFARI_ALL_with_guide$RNA@scale.data), lowerDetectionLimit = 0.5, expressionFamily = VGAM::negbinomial.size())

SFARI_all <- detect_expression(SFARI_all)

some_mon <- SFARI_all
expressed_genes <- row.names(subset(fData(some_mon), num_cells_expressed >= 60))
some_mon <- some_mon[expressed_genes]
some_mon <- some_mon[,pData(some_mon)$num_genes_expressed > 500 & pData(some_mon)$num_genes_expressed < 5000]
some_mon <- estimateSizeFactors(some_mon)
SFARI_all <- some_mon

SFARI_all <- add_sgRNA_monocle(SFARI_all, all_singlets)  ## Add identified gRNA to cells 
SFARI_all_mon <- SFARI_all[,pData(SFARI_all)$gene_level > 0]   ### SFARI_with_guide 


#ordering_genes <- VariableFeatures(object = SFARI_ALL_with_guide)  ##  Used genes from bulk-time-course time-point-DE-genes but all pseduotime orderings based on highly variable genes gave broadly similar results because the variable genes ARE the neuronal differentiation genes 

temporal_genes <- read.table('Temporal_gene_list_from_bulk.txt')  ## These genes are differentially expressed across timepoints in bulk RNA-seq 
ordering_genes <- as.character(temporal_genes$V1)

SFARI_all_mon <- setOrderingFilter(SFARI_all_mon, ordering_genes)

SFARI_all_mon_pt2 <- reduceDimension(SFARI_all_mon, max_components = 2, norm_method="none")
gc()

## Next command takes 30 minutes to 1 hour, and eats all my RAM
SFARI_all_mon_pt3 <- orderCells(SFARI_all_mon_pt2)
gc()

## Re-number PT states from left to right 

SFARI_all_mon_pt3$State[which(SFARI_all_mon_pt3$State == 2)] <- rep(3, length(which(SFARI_all_mon_pt3$State == 2)))
SFARI_all_mon_pt3$State[which(SFARI_all_mon_pt3$State == 7)] <- rep(2, length(which(SFARI_all_mon_pt3$State == 7)))
levels(SFARI_all_mon_pt3$State) <- c(levels(SFARI_all_mon_pt3$State),10) ## add dummy level for swap
SFARI_all_mon_pt3$State[which(SFARI_all_mon_pt3$State == 6)] <- rep(10, length(which(SFARI_all_mon_pt3$State == 6)))
SFARI_all_mon_pt3$State[which(SFARI_all_mon_pt3$State == 5)] <- rep(6, length(which(SFARI_all_mon_pt3$State == 5)))
SFARI_all_mon_pt3$State[which(SFARI_all_mon_pt3$State == 10)] <- rep(5, length(which(SFARI_all_mon_pt3$State == 10)))
SFARI_all_mon_pt3$State <- droplevels(SFARI_all_mon_pt3$State) ### drop empty levels 

saveRDS(SFARI_all_mon_pt3, file = "SFARI_pseudotime_processed.rds") # save this file for future re-analysis

### Plot cells in Pseduotime   Figure 3A
plot_cell_trajectory(SFARI_all_mon_pt3, theta=180, cell_size=2, show_branch_points = F) + scale_color_manual(values = myPal2)  + theme_classic(base_size=20)

state1genes <- getSeuratPseudoTimeDE(SFARI_all_mon_pt3, SFARI_ALL_with_guide, 1) # identifying positive marker genes for cells in a given pseduotime state
plotGeneSetBulk(row.names(head(state1genes, n=20)), lrt)

plotGeneSetDplyr(c('HES7','VIM','TP53', 'TCF7L1','BAZ1A','CCND2','NCOR2'), lrt) ## state1 genes
plotGeneSetDplyr(c('AGPAT2','CD63','EMP3','GNG5','GPX1','NACA','RAN', 'CHRNA1','TSPO'), lrt) ## state2
plotGeneSetDplyr(c('CALM1','CCER2','GNG5','PDLIM7','SNCG','TMSB4X', 'BSG','CTNNB1'), lrt) ## state3
plotGeneSetDplyr(c('CALM1','GAP43','RTN1','STMN2','TUBB2B','SYT1', 'RAB3A'), lrt) ## state4
plotGeneSetDplyr(c('CXXC4','FXYD7','PEG10','SOX4','TLE1','CELF4','EBF3','MAP2','KMT2A'), lrt) ## state5
plotGeneSetDplyr(c('DCX','MAP1B','NEFL','NEFM','ANK3','NCAM1','CD24','NRXN1'), lrt) ### state6

```



``` {r, Back to Seurat for more analysis, include =T }

### Transfer monocle labels to Seurat dataset (slow function)
SFARI_ALL_with_guide <- add_pseudotime_state_to_Seurat(SFARI_all_mon_pt3, SFARI_ALL_with_guide)

## UMAP highlighting small cluster of cells at the bottom 

DimPlot(SFARI_ALL_with_guide, cells.highlight=colnames(SFARI_ALL_with_guide)[SFARI_ALL_with_guide$RNA_snn_res.0.25 %in% c(4)], cols.highlight = "purple4", cols="lightgrey", pt.size=1)

### Adding dichotomized pseudotime labels 

PT_binary <- rep('early', length(colnames(SFARI_ALL_with_guide)))
PT_binary[which(SFARI_ALL_with_guide$pt_state > 3)] <- "late"
SFARI_ALL_with_guide$pt_binary <- PT_binary
SFARI_ALL_with_guide$pt_guide <- paste(SFARI_ALL_with_guide$gene_level, SFARI_ALL_with_guide$pt_binary)

## Recluster extremes from pseduotime 
Idents(SFARI_ALL_with_guide) <- SFARI_ALL_with_guide$pt_state
extreme_PT <- subset(SFARI_ALL_with_guide, idents=c("1","6"))
extreme_PT <- PCA_TSNE(extreme_PT)
Idents(extreme_PT) <- extreme_PT$pt_state

### Pseudotime States 1 and 6 are clearly distinct 'Cell States' by UMAP, defined by progressoin of neuronal maturation. Again, note PC1 genes (STMN2, NEFM, MAP1B, DCX)
### Supplementary UMAP figures 
DimPlot(object = extreme_PT, reduction = "umap", pt.size = 1, cols=c("lightgrey","purple4"))
FeaturePlot(object = extreme_PT, features = "STMN2", pt.size = 1, cols=c("lightgrey","lightgrey","lightgrey","purple4"))
FeaturePlot(object = extreme_PT, features = "NEUROD1", pt.size = 1, cols=c("lightgrey","lightgrey","purple4"))

### Figure 3B

Idents(SFARI_ALL_with_guide) <- factor(SFARI_ALL_with_guide$pt_state)
VlnPlot(object = SFARI_ALL_with_guide, features = c('MAP2','DCX','TP53','CDK4'), cols=myPal2)  

### Figure 3D 

a <- sapply(split(SFARI_all_mon_pt3$State, SFARI_all_mon_pt3$gene_level), function(x){table(x)/sum(table(x))})
plot(hclust(dist(t(a))))  ### sort the columns to make the figure prettier 
b<- hclust(dist(t(a[-4,])), method="ward.D2")
barplot(a[,rev(b$order)], col=myPal2, cex.names=1.1, cex.axis = 1.5, las=2 )  

## Also double check pseudotime labels correctly added to Seurat data: 
a <- sapply(split(SFARI_ALL_with_guide$pt_state, SFARI_ALL_with_guide$gene_level), function(x){table(x)/sum(table(x))})
plot(hclust(dist(t(a))))  ### sort the columns to make the figure prettier 
b<- hclust(dist(t(a[-4,])), method="ward.D2")
barplot(a[,rev(b$order)], col=myPal2, cex.names=1.1, cex.axis = 1.5)  

### Separating by sgRNA level 

a <- sapply(split(SFARI_all_mon_pt3$State, SFARI_all_mon_pt3$sgRNA), function(x){table(x)/sum(table(x))})
barplot(a, col=myPal2, cex.names=1.1, cex.axis = 1.5, las=2)

## Plotting Pseudotime state mean +/- SEM 

### Pre-revision Figure 3E -- average pseduotime state 
plotData <- as.data.frame(cbind(SFARI_all_mon_pt3$State, SFARI_all_mon_pt3$gene_level))
colnames(plotData) <- c("State","Guide")
pd2 <- plotData %>% group_by(Guide) %>% summarise(mean = mean(as.numeric(State)), sem = sd(as.numeric(State))/sqrt(n()))
pd3 <- subset(pd2, pd2$Guide %in% c('ARID1B','ASH1L','CHD2','PTEN','CHD8','DYRK1A','Nontargeting'))
g <- ggplot(pd3, aes(reorder(Guide, mean), mean))
g +  geom_point(size=5) + geom_errorbar(aes(ymin=mean-sem, ymax=mean+sem, width=0.3)) +  ylim(1,6) + geom_hline(yintercept=3.94, linetype="dashed", color = "black", size=1) + theme_tufte(base_size=30)


### Figure 3E - individual sgRNA estimates of average pseudotime-  point estimates came from indiv_sgRNA_pt_enrich_test(SFARI_all_mon_pt3, GOI)

guide_level_PT_point_estimates <- read_excel("guide_level_PT_point_estimates.xlsx")
pd4 <- gather(guide_level_PT_point_estimates, key=Level, value=Expression, Guide1:All)
pd4$Gene <- factor(pd4$Gene, levels=c('CHD2','ASH1L','ARID1B','DYRK1A','Nontargeting','CHD8','PTEN'))
pd4$Group1 <- pd4$Level == "All"
g <- ggplot(pd4, aes(Gene, Expression, label=Level, group=Group1, color=Group1, shape=Group1))
g + geom_point(size=5, position=position_dodge(.4)) + geom_text(aes(label=Level),hjust=2, vjust=0) + ylim(2,5) + geom_hline(yintercept=3.94, linetype="dashed", color = "black", size=1) + theme_grey(base_size = 30) + scale_color_grey(start = .7, end = 0) 

### Cells in early / late pseduotime by sgRNA

myPal <- viridis(15)[c(4,10)]  ## Two colors for batch, colors from viridis so not to be confused with pseudotime colors

SFARI_batch1 <- subset(SFARI_ALL_with_guide, batch == "1")
SFARI_batch2 <- subset(SFARI_ALL_with_guide, batch == "2" | batch =="3")
SFARI_ALL_with_guide$batch2 <- SFARI_ALL_with_guide$batch
SFARI_ALL_with_guide$batch2[SFARI_ALL_with_guide$batch == "3"] <- "2"

a <- sapply(split(SFARI_ALL_with_guide$batch2, SFARI_ALL_with_guide$gene_level), function(x){table(x)/sum(table(x))})
plot(hclust(dist(t(a))))  ### sort the columns to make the figure prettier 
b<- hclust(dist(t(a[-4,])), method="ward.D2")
barplot(a[,rev(b$order)], col=myPal, cex.names=1.1, cex.axis = 1.5)  

a <- sapply(split(SFARI_batch1$pt_state, SFARI_batch1$gene_level), function(x){table(x)/sum(table(x))})
barplot(a, col=myPal2, cex.names=.9, cex.axis = 1.5, las=2)

## Number of cells in each pseudotime state for each batch 
taba <- table(SFARI_batch1$pt_state, SFARI_batch1$gene_level)
tabb <- table(SFARI_batch2$pt_state, SFARI_batch2$gene_level)
barplot(t(rbind(rowSums(taba)/sum(rowSums(taba)), rowSums(tabb)/sum(rowSums(tabb)))), col=myPal2)

```


```{r, Differential Gene Expression Analysis }

## Differential expression by target gene 

Idents(SFARI_ALL_with_guide) <- SFARI_ALL_with_guide$gene_level

for (GOI in SFARI_genes[-c(10,11,13,14)])   ## myt1l and reln had no DE, pogz + setd5 had many due to false positives (few cells)  
 { 
de_GOI <- FindMarkers(SFARI_ALL_with_guide, ident.1 = GOI, logfc.threshold = 0.1)
write.table(de_GOI, file=paste("~/Documents/2_19_2019_LFC1_de_in",GOI,sep="_"), quote=F, row.names=T, col.names=T)
}


## Differentiatial expression by Target and Pseudotime Staage 

## Early Pseudotime 
Idents(SFARI_ALL_with_guide) <- SFARI_ALL_with_guide$pt_guide   

for (GOI in SFARI_genes[-c(11,14)])  ## Exclude POGZ and SETD5
{ 
de_GOI <- FindMarkers(SFARI_ALL_with_guide, ident.1 = paste(GOI,"early"), ident.2="Nontargeting early", logfc.threshold = 0.1)
write.table(de_GOI, file=paste("~/Documents/2_18_2019_SFARI_all_cells_EARLY_DE_in",GOI,sep="_"), quote=F, row.names=T, col.names=T)
}

## Late Pseudotime 

for (GOI in SFARI_genes[-c(11,13)])  ## Exclude POGZ and RELN
{ 
de_GOI <- FindMarkers(SFARI_ALL_with_guide, ident.1 = paste(GOI,"late"), ident.2="Nontargeting late", logfc.threshold = 0.1)
write.table(de_GOI, file=paste("~/Documents/2_18_2019_SFARI_all_cells_LATE_DE_in",GOI,sep="_"), quote=F, row.names=T, col.names=T)
}


## Grouping Average Transcriptomes by Certain Featuers (Tareget gene, pseduotime state, target gene by pseudotime stage etc)

Idents(SFARI_ALL_with_guide) <- SFARI_ALL_with_guide$gene_level 
GroupedByGeneSeurat <- AverageExpression(subset(SFARI_ALL_with_guide, idents="Nontargeting"))$RNA
for (GOI in SFARI_genes) { 
all_cells_with_guide <- AverageExpression(subset(SFARI_ALL_with_guide, idents=GOI))$RNA
GroupedByGeneSeurat <- cbind(GroupedByGeneSeurat, all_cells_with_guide)
}
GroupedByGeneSeurat <- data.matrix(GroupedByGeneSeurat)

Idents(SFARI_ALL_with_guide) <- SFARI_ALL_with_guide$pt_state 
GroupedByPTSeurat <- AverageExpression(subset(SFARI_ALL_with_guide, idents="1"))$RNA
for (GOI in 2:6) { 
all_cells_with_guide <- AverageExpression(subset(SFARI_ALL_with_guide, idents=GOI))$RNA
GroupedByPTSeurat <- cbind(GroupedByPTSeurat, all_cells_with_guide)
}
GroupedByPTSeurat <- data.matrix(GroupedByPTSeurat)

PT_state_markers <- c('HES7','VIM','TP53', 'TCF7L1','BAZ1A','CCND2','NCOR2', 'AGPAT2','CD63','EMP3','GNG5','GPX1','NACA','RAN', 'CHRNA1','TSPO','CALM1','CCER2','GNG5','PDLIM7','SNCG','BSG','CTNNB1','GAP43','RTN1','STMN2','TUBB2B','SYT1', 'RAB3A','CXXC4','FXYD7','PEG10','SOX4','TLE1','CELF4','EBF3','MAP2','KMT2A','DCX','MAP1B','NEFL','NEFM','ANK3','NCAM1','CD24','NRXN1')


## Visualizing sets of cells grouped by pseudotime 

heatData <- GroupedByPTSeurat[PT_state_markers,]
heatmap(heatData, scale="row")
my_dend2 <- as.dendrogram(hclust(dist(t(heatData))))
my_dend2 <- rev(reorder(my_dend2, c(100,10,100,-5,1,1,1,1,0)))
a <- Heatmap(t(scale(t(heatData))), column_dend_reorder = F, cluster_columns = my_dend2, show_row_dend = F, km=6, col=magma(10)) 


## time specific dot-plot 

Idents(SFARI_ALL_with_guide) <- SFARI_ALL_with_guide$pt_guide
### Colors scaled to own show repression, so any circle means some degree of repression.  Repression at the early or late pseudotime stage 

low_col <- magma(12)[3]
mid_col <- magma(12)[7]

## Early 
b <- DotPlot(subset(SFARI_ALL_with_guide, idents=paste(geneSet, "early")), features = geneSet, dot.scale=20, scale.by='size') +  scale_color_viridis(option="magma")
b$data$id <- factor(b$data$id, rev(sort(levels(b$data$id))))
b + scale_color_gradient2(low = low_col, high="white")  + theme_classic(base_size=20)

## Late
b <- DotPlot(subset(SFARI_ALL_with_guide, idents=paste(geneSet, "late")), features = geneSet, dot.scale=20, scale.by='size') +  scale_color_viridis(option="magma")
b$data$id <- factor(b$data$id, rev(sort(levels(b$data$id))))
b + scale_color_gradient2(low = low_col, high="white")  + theme_classic(base_size=20)

                    
Idents(SFARI_ALL_with_guide) <- SFARI_ALL_with_guide$pt_guide

GroupedByGenePTSeurat <- AverageExpression(subset(SFARI_ALL_with_guide, idents="Nontargeting early"))$RNA
for (GOI in SFARI_genes[-13]) { 
all_cells_with_guide <- AverageExpression(subset(SFARI_ALL_with_guide, idents=paste(GOI, "early")))$RNA
GroupedByGenePTSeurat <- cbind(GroupedByGenePTSeurat, all_cells_with_guide)
all_cells_with_guide <- AverageExpression(subset(SFARI_ALL_with_guide, idents=paste(GOI, "late")))$RNA
GroupedByGenePTSeurat <- cbind(GroupedByGenePTSeurat, all_cells_with_guide)
}
all_cells_with_guide <- AverageExpression(subset(SFARI_ALL_with_guide, idents=paste("Nontargeting", "late")))$RNA
GroupedByGenePTSeurat <- cbind(GroupedByGenePTSeurat, all_cells_with_guide)
GroupedByGenePTSeurat <- data.matrix(GroupedByGenePTSeurat)
```

```{r, PCA plot of cells grouped by pseduotime stage, include = T}
## PCA plot of pseduo-time stratified samples, grouped by target gene 
myPal <- viridis(15)[c(4,10)]

pca1 <- prcomp(t(log2(GroupedByGenePTSeurat[,c(-22,-23,-26,-27)]+1)))
scores = as.data.frame(pca1$x)
percentage <- round(pca1$sdev / sum(pca1$sdev) * 100, 2)
percentage <- paste( colnames(pca1$x), "(", paste( as.character(percentage), "%", ")", sep="") )
scores$State <- sapply(strsplit(row.names(scores), " "), "[", 2)
scores$Gene <- sapply(strsplit(row.names(scores), " "), "[", 1)

ggplot(data = scores, aes(x = PC1, y = PC2, col = State, shape = State, label = c(Gene))) + geom_text_repel(size=6) + theme_gray(base_size = 20) + geom_point(size=5) +  xlab(percentage[1]) + ylab(percentage[2]) + scale_color_manual(values = myPal)

```

``` {r, Hierarchical clustering plots, include = T }

de_219_log2_1ormore <- read.table('DE_LFC2_uniq', head=F)  ## log2 FC > 2 in 2 or more knockdown conditions 

exclude <- c(-22,-23,-26,-27)
exclude_early <- c(-12,-14,-15)


## Figure S6, panel A - recurrently dyresgulated genes, heatmap 

exclude <- c(-12,-14,-15)  ## Exclude variable samples (too few cells) 
heatData <- GroupedByGeneSeurat[which(row.names(GroupedByGeneSeurat)%in%as.character(de_219_log1_2ormore$V1)), exclude]
heatmap(heatData, col=myPal2)
my_dend2 <- as.dendrogram(hclust(dist(t(heatData))))
my_dend2 <- rev(reorder(my_dend2, c(100,0,-2,-5,1,1,1,1,0,0,0,0)))
Heatmap(t(scale(t(heatData))), column_dend_reorder = F, cluster_columns = my_dend2, show_row_dend = F, km=5, col=magma(10))

#### 

PT_early <- grep('early', colnames(GroupedByGenePTSeurat))
PT_late <- grep('late', colnames(GroupedByGenePTSeurat))

## Figure 4E clustering by neuronal projection genes 

neuron_projection <- c('ACAP3','APBB1','CCDC88A','CDH2','CNTN2','CRABP2','CRMP1','DPYSL3','EFNA5','FKBP4','GSK3B','HSPA5','ITM2C','MAP1B','MAPT','MEF2C','NEFL','PTPRS','SERPINF1','STYXL1','UBE2V2','ZEB2','ZMYND8')
heatData <- GroupedByGenePTSeurat[neuron_projection, PT_late]
my_dend2 <- as.dendrogram(hclust(dist(scale(t(heatData[,exclude_late])))))
my_dend_neuron_projection <- rev(reorder(my_dend2, c(200,0,0,0,0,0,0,0,0,0,100,0)))
Heatmap(t(scale(t(heatData[,exclude_late]))), col=magma(10), show_row_dend = F, split=3, column_dend_reorder = F, cluster_columns = my_dend_neuron_projection)


## These files are just the gene names of any DE gene from any target knock-down condition 

ALL_early_genes <- read.table('Early_ALL_DE_genes', head=F, row.names=1)
ALL_late_genes <- read.table('Late_ALL_DE_genes', head=F, row.names=1)
ALL_earlylate_genes <- read.table('DE_Early_or_late', head=F, row.names=1)

### Making the binary matrix of DE yes/no for each target -- read in your DE lists for each gene 

for (GOI in SFARI_genes[c(-10,-11,-13,-14)]) {

early <- rep(0,length(row.names(ALL_early_genes)))
early_genes <- read.table( paste('~/Documents/SFARI_Analysis/Differential Expression/EARLY', GOI, sep="_"), head=F)
early[which(row.names(ALL_early_genes)%in%early_genes$V1)] <- 1

ALL_early_genes <- cbind(ALL_early_genes, early)
}

## Clustering samples at Early Pseudotime 

names(ALL_early_genes) <- SFARI_genes[c(-10,-11,-13,-14)]
heatmap(data.matrix(ALL_early_genes[which(rowSums(ALL_early_genes) > 3 & rowSums(ALL_early_genes) < 9),]), scale="none", col=myPal2[c(-7)])

early_cluster_genes <-row.names(ALL_early_genes[which(rowSums(ALL_early_genes) > 3 & rowSums(ALL_early_genes) < 9),])

## Binary matrix for all the DE genes from the late stage 
for (GOI in SFARI_genes[c(-10,-11,-13,-14)]) {

late <- rep(0,length(row.names(ALL_late_genes)))
late_genes <- read.table( paste('~/Documents/SFARI_Analysis/Differential Expression/LATE', GOI, sep="_"), head=F)
late[which(row.names(ALL_late_genes)%in%late_genes$V1)] <- 1

ALL_late_genes <- cbind(ALL_late_genes, late)
}

## Clustering LATE stage samples 
names(ALL_late_genes) <- SFARI_genes[c(-10,-11,-13,-14)]
heatmap(data.matrix(ALL_late_genes[which(rowSums(ALL_late_genes) > 2 & rowSums(ALL_late_genes) < 8),]), scale="none", col=myPal2[c(-7)])
late_cluster_genes <- row.names(data.matrix(ALL_late_genes[which(rowSums(ALL_late_genes) > 2 & rowSums(ALL_late_genes) < 9),]))


### Clustering all Samples 

## Binary matrix for DE at early first 
for (GOI in SFARI_genes[c(-10,-11,-13,-14)]) {

early_genes <- read.table( paste('~/Documents/SFARI_Analysis/Differential Expression/EARLY', GOI, sep="_"), head=F)
early <- rep(0,length(row.names(ALL_earlylate_genes)))
early[which(row.names(ALL_earlylate_genes)%in%early_genes$V1)] <- 1

ALL_earlylate_genes <- cbind(ALL_earlylate_genes, early)
}

## Again for late stage genes

for (GOI in SFARI_genes[c(-10,-11,-13,-14)]) {

late_genes <- read.table( paste('~/Documents/SFARI_Analysis/Differential Expression/LATE', GOI, sep="_"), head=F)
late <- rep(0,length(row.names(ALL_earlylate_genes)))
late[which(row.names(ALL_earlylate_genes)%in%late_genes$V1)] <- 1

ALL_earlylate_genes <- cbind(ALL_earlylate_genes, late)
}

colnames(ALL_earlylate_genes)[1:10] <-  paste(SFARI_genes[c(-10,-11,-13,-14)],"Early")
colnames(ALL_earlylate_genes)[11:20] <-  paste(SFARI_genes[c(-10,-11,-13,-14)],"Late")


# Heatmap(data.matrix(ALL_earlylate_genes[which(rowSums(ALL_earlylate_genes) > 3 & rowSums(ALL_earlylate_genes) <= 12),]), col=myPal2[c(1,5)], km=8) 

clusterThese <- unique(c(early_cluster_genes, late_cluster_genes))
clusterThese2 <- unique(c(clusterThese, row.names(data.matrix(ALL_earlylate_genes[which(rowSums(ALL_earlylate_genes) > 3 & rowSums(ALL_earlylate_genes) <= 12),]))))


#Heatmap(data.matrix(ALL_earlylate_genes[clusterThese2,]), col=myPal2[c(1,5)], km=8)

my_dend <- as.dendrogram(hclust(dist(t(data.matrix(ALL_earlylate_genes[clusterThese2,])))))
my_dend <- rev(reorder(my_dend, c(1000,100,200,0,100,0,0,0,0,1000,0,-1000,-1000,-1,-1,-1,-1,-1000,-1000,100)))

## Figure 4A 

Heatmap(data.matrix(ALL_earlylate_genes[clusterThese2,]), col=myPal2[c(1,5)], km=8, column_dend_reorder = F, cluster_columns = my_dend, show_row_dend = F) # row clustering is somewhat stochastic, this does not change enrichment analyses or column clustering 

```



``` {r, Pseudotime Chi-squared enrichment testing }

geneSet <- unique(SFARI_all_mon_pt3$gene_level)

sink('GuideRNA_Pseudotime_Enrichment')  
for(i in 1:length(geneSet))
{
cat('________________________________________________________\n')
GOI <- geneSet[i] 
cat(GOI)
sgRNA_pt_enrich_test(SFARI_all_mon_pt3, GOI)
}
sink()

geneSet <- paste(unique(SFARI_all_mon_pt3$sgRNA))

sink('Individual_Guide_Pseudotime_Enrich')
for(i in 1:length(geneSet))
{
cat('________________________________________________________\n')
GOI <- geneSet[i] 
cat(GOI)
indiv_sgRNA_pt_enrich_test(SFARI_all_mon_pt3, GOI)
}
sink()

```



``` {r, Visualizing sgRNA efficiency (not used), include = T  }

Idents(SFARI_ALL_with_guide) <- SFARI_ALL_with_guide$gene_level
GOI <- "ASH1L"

ctrl_mean <- mean(GetAssayData(SFARI_ALL_with_guide, slots=scale.data)[GOI, WhichCells(SFARI_ALL_with_guide, idents="Nontargeting")])

pct_untarg <- table(GetAssayData(SFARI_ALL_with_guide, slots=scale.data)[GOI, WhichCells(SFARI_ALL_with_guide, idents="Nontargeting")] == 0)[1] / length(GetAssayData(SFARI_ALL_with_guide, slots=scale.data)[GOI, WhichCells(SFARI_ALL_with_guide, idents="Nontargeting")] == 0)  ### percent expressing in controls 

##targeted_mean
targeted_mean <- mean(GetAssayData(SFARI_ALL_with_guide, slots=scale.data)[GOI, WhichCells(SFARI_ALL_with_guide, idents=GOI)])

## percent expressing in targeted
pct_targ <- table(GetAssayData(SFARI_ALL_with_guide, slots=scale.data)[GOI, WhichCells(SFARI_ALL_with_guide, idents=GOI)] == 0)[1] / length(GetAssayData(SFARI_ALL_with_guide, slots=scale.data)[GOI, WhichCells(SFARI_ALL_with_guide, idents=GOI)] == 0) 

kd_eff <- matrix(c(ctrl_mean, targeted_mean, targeted_mean/ctrl_mean, pct_untarg, pct_targ, pct_targ/pct_untarg),ncol=3,byrow=TRUE)
colnames(kd_eff) <- c("Nontargeting",paste("Targeting\n",GOI),"Effect Size\n (Targeted/Nontargeted)")
rownames(kd_eff) <- c("Mean     \nExpression ", "Percent of    \n Cells Expressing")

a <- VlnPlot(SFARI_ALL_with_guide, features=c(GOI), idents=c(GOI,"Nontargeting"), log=T, adjust=2, pt.size = 0.1, cols=c("lightgrey","purple4"))
a$labels$x <- ""
b <- tableGrob(round(kd_eff, 3))
grid.arrange(a, b, nrow = 1)

```




``` {r, Gene Ontology Enrichment plots, include = T}

goMF <- read.csv('~/Documents/LUHMES_GO_MF.csv', head=F)
goDIS <- read.csv('~/Documents/LUHMES_GO_DISEASE.csv', head=F)
goSC <- read.csv('~/Documents/single_cell_recurrent_enrich.csv', head=F)

upCortexandLuhmes <- read_xlsx('GO LUHMES Cortex UP UP.xlsx')
dd <- read_xlsx('GO_LUHMES_cortex_DOWN.xlsx')
ud <- read_xlsx('GO_LUHMES_cortex_up_luhmes_down.xlsx')
uu <- read_xlsx('GO_LUHMES_UP_CORTEX_UP.xlsx')
du <- read_xlsx('GO_LUHMES_increase_cortex_decrease.xlsx')

dd <- du
dd$Description <- factor(dd$Description, dd$Description[order(as.numeric(dd$`P Value`))])

ggplot(dd[1:10,], aes(Description, -log10(as.numeric(`P Value`)))) + geom_bar(stat="identity", width=0.4, fill=magma(10)[6]) + coord_flip() + scale_x_discrete(limits = rev(levels(dd$Description)[1:10]))

goEarly <- read.csv('~/Documents/Early3_GOBP.csv', head=F)
goLate <- read.csv('~/Documents/Late2_GOBP.csv', head=F)

recDE <- goEarly
recDE2 <- goLate
#recDE <- rbind(recDE, recDE2)
recDE <- recDE2

colnames(recDE) <- c("GO_name","Pathway", "DontCare","AlsoDontCare","Ditto", "Enrichment", "FDR")
recDE$Enrichment <- -log10(recDE$Enrichment)

#recDE <- goSC
#colnames(recDE) <- c("Pathway","Enrichment")

recDE$Pathway <- factor(recDE$Pathway, recDE$Pathway[order(recDE$Enrichment)])

ggplot(recDE, aes(Pathway, Enrichment), fill=2, width=1)  +
  geom_bar(stat="identity", width=.4) + 
  coord_flip() + theme_grey(base_size=20) 

```





``` {r, Analyze AllenBrain data, include = T }

### Allen Brain data downloaded from:  https://www.brainspan.org/api/v2/well_known_file_download/267666525 
AllenBrain <- read.csv('~/Downloads/AllenBrain/expression_matrix.csv', head=F)
AllenRow <- read.csv('~/Downloads/AllenBrain/rows_metadata.csv', head=T)
AllenCol <- read.csv('~/Downloads/AllenBrain/columns_metadata.csv', head=T)
colnames(AllenBrain) <- paste(AllenCol$age,AllenCol$structure_acronym, AllenCol$donor_id)
rownames(AllenBrain) <- make.unique(as.character(AllenRow$gene_symbol))

AllenBrainSubset <- AllenBrain[rownames(AllenBrain)%in%rownames(lrt),]  ### genes that are in my LUHMES dataset

corticalSamples <- colnames(AllenBrain)[grep('FC', colnames(AllenBrain))]
CorticalSubset <- AllenBrainSubset[,corticalSamples]

fetalSamples <- colnames(AllenBrain)[grep('pcw', colnames(AllenBrain))]
FetalSubset <- AllenBrainSubset[,fetalSamples[-c(1:75)]]  ### exclude E9,10,12 due to progenitors, and 37 pcw due to outlier on MDS
FetalSubset <- FetalSubset[,-grep(37, colnames(FetalSubset))]

region <- sapply(strsplit(colnames(FetalSubset), " "), "[", 3)
FetalDGE <- DGEList(FetalSubset, group=factor(region))
plotMDS(FetalDGE)

groupLabel <- rep('NA',length(CorticalSubset))
groupLabel[grep("^8 pcw|^9 pcw", names(CorticalSubset))] <- "EarlyEightTen"
groupLabel[grep("^12 pcw|^13 pcw", names(CorticalSubset))] <- "EarlyTenThirteen"
groupLabel[grep("^16 pcw", names(CorticalSubset))] <- "EarlySixteen"
groupLabel[grep("^17 pcw|^19 pcw", names(CorticalSubset))] <- "EarlyMidSeventeenNineteen"
groupLabel[grep("^21 pcw|^24 pcw|^26 pcw|^35 pcw|^37 pcw", names(CorticalSubset))] <- "Late2438"
groupLabel[grep("^4 mos", names(CorticalSubset))] <- "BirthSix"
groupLabel[grep("^10 mos|^1 yrs", names(CorticalSubset))] <- "SixMonthsOneYear"
groupLabel[grep("^2 yrs|^3 yrs|^4 yrs", names(CorticalSubset))] <- "OneYearSixYear"
groupLabel <- factor(groupLabel,levels=c('EarlyEightTen','EarlyTenThirteen','EarlySixteen','EarlyMidSeventeenNineteen','Late2438','BirthSix','SixMonthsOneYear','OneYearSixYear'))

myDatainAllen <- lrt[rownames(lrt)%in%rownames(AllenBrainSubset),]  ### LUHMES genes in Allen 

### Allen Brain / LUHMES overlap re-analysis 
## LUHMES DE paramaters (coef 5, lenient FC)

timeCourse <- read.table('Timecourse_LUHMES_counts.txt', sep='\t', head=T)
groups <- as.factor(c('8','8','6','6','4','4','2','2','0','0'))
rownames(timeCourse) = make.names(timeCourse[,3], unique=TRUE)
timeCourseDGE <- DGEList(timeCourse[,c(-1,-2,-3,-4,-5)], group=groups)
x <- timeCourseDGE
## x$samples$group <- relevel(x$samples$group, ref="0") This will set day4 to reference so DE genes will be wrt to this point
selr1 <- rowSums(cpm(x$counts)>20)>=2  ###  count > 20 in 2 or more samples
x <- x[selr1,]
x <- calcNormFactors(x)
design <- model.matrix(~x$samples$group)
xglm <- estimateGLMCommonDisp(x, design, verbose=T)
fit <- glmFit(xglm, design)
lrt <- glmLRT(fit, coef=5)  ###  DE from day 0 to 5 = genes in neural stem cells and neuron diff
deGenes <- topTags(lrt, n=Inf)
deGenes <- deGenes$table
deSig <- subset(deGenes, deGenes$FDR < 0.0001 & abs(deGenes$logFC) > .5)  ## 3215, and 3115


myDatainAllen <- lrt[rownames(lrt)%in%rownames(AllenBrainSubset),]  ### LUHMES genes in Allen
AllenDGE <- DGEList(CorticalSubset[,1:88], group=groupLabel[1:88])
x2 <- AllenDGE
selr1 <- rowSums(cpm(x2$counts)>20)>=5  ###  count > 20 in 2 or more samples
x2 <- x2[selr1,]
x2 <- calcNormFactors(x2)
x2$samples$group <- relevel(x2$samples$group, ref='BirthSix')
design <- model.matrix(~x2$samples$group)
xglm2 <- estimateGLMCommonDisp(x2, design, verbose=T)
contMat <- makeContrasts(EarlyTenThirteen-EarlyEightTen, levels=levels(x2$samples$group))
fit2 <- glmFit(xglm2, design)
lrt2 <- glmLRT(fit2, contrast = contMat)
deGenes <- topTags(lrt2, n=2500)
deGenes <- deGenes$table
oerlap1 <- deSig[which(row.names(deSig)%in%row.names(deGenes)),1:2]   ### deSig is all DE from LUHMES d0-d8, deGenes = top 2500 DE from E8_10v10_13

### LUHMES derived from 8‐week‐old human ventral mesencephalic tissue, and compares well E8_10 v E10_13
overlap2 <- deGenes[rownames(oerlap1),1:2]   
overlaps <- cbind(oerlap1, overlap2)
plot(overlaps[,c(1,3)])
cor(overlaps[,c(1,3)])
cor.test(overlaps[,3],overlaps[,1])
dat <- overlaps[,c(1,3)]
dat$density <- get_density(dat$logFC, dat$logFC.1, n = 100)
ggplot(dat, aes(x=logFC, y=logFC.1)) + geom_point(aes(logFC, logFC.1, color = density)) + scale_color_viridis(option='magma')

table(dat$logFC >0 & dat$logFC.1 > 0)
table(dat$logFC < 0 & dat$logFC.1 > 0)
design
cor.test(overlaps[,3],overlaps[,1])
cor(overlaps[,c(1,3)])
ggplot(dat, aes(x=logFC, y=logFC.1)) + geom_point(aes(logFC, logFC.1, color = density)) + scale_color_viridis(option='magma') + theme_grey(base_size=20)

## Up in LUHMES, down in Cortex
subset(dat, dat$logFC > 0 & dat$logFC.1 < 0)

BothUP <- subset(dat, dat$logFC >0 & dat$logFC.1 > 0)
#write.table(BothUP, file='LUHMES_cortex_UP.txt')

BothDOWN <- subset(dat, dat$logFC < 0 & dat$logFC.1 < 0)
#write.table(BothDOWN, file='LUHMES_cortex_DOWN.txt')

LUHMES_decrease_up_cortex <- subset(dat, dat$logFC < 0 & dat$logFC.1 > 0)
#write.table(LUHMES_decrease_up_cortex, file='LUHMES_decreased_up_cortex.txt')

Cortex_decrease_up_LUHMES <- subset(dat, dat$logFC > 0 & dat$logFC.1 < 0)
#write.table(Cortex_decrease_up_LUHMES, file='Cortex_decrease_up_LUHMES.txt')

```


``` {r, Parikshak et al data SFARI genes, include = T }

## ASD candidate gene module assignments from Parikshak et al 2013 Cell 
Parikshak <- read.csv('SFARI_Parikshak.csv', head=T)
plotThis <- as.data.frame(Parikshak[,c(2,10:18)])
row.names(plotThis) <- plotThis[,1]
colnames(plotThis) <- c("Gene","PCW8","PCW12","PCW13","PCW16","PCW17.22","PCW25.26","FourMonths","TenMonths","TwelveMonths")
plotThis <- gather(plotThis, key = Timepoint, value = Expression, PCW8:TwelveMonths)
plotThis$TimePoint <- factor(plotThis$Timepoint, levels= c("PCW8","PCW12","PCW13","PCW16","PCW17.22","PCW25.26","FourMonths","TenMonths","TwelveMonths"))
ggplot(plotThis, aes(y=log2(Expression), x=TimePoint, group=Gene, color=Gene)) + geom_line(lwd=2)

```

``` {r, Eichler ASD_ID dataset, include = T }
## Data from Stessman et al , 2017 Nature Genetics 

Eichler <- read.csv('Eichler_ID_ASD.csv', head=T)
row.names(Eichler) <- Eichler$Gene

my_genes <- c('ADNP','ARID1B','ASH1L','DYRK1A','PTEN','CHD2','CHD8')
my_data <- Eichler[my_genes, -1]

my_d3 <- as.dendrogram(hclust(dist(scale(Eichler[my_genes,-1]))))
my_d3 <- reorder(my_d3, c(5,50,-10,-10,-10,0))

Heatmap(t(scale(Eichler[my_genes,-1])), column_dend_reorder = F, cluster_columns = my_d3, col=myPal2, cell_fun = function(j, i, x, y, width, height, fill) {
    grid.text(sprintf("%.0f%%", my_data[j, i]), x, y, gp = gpar(fontsize = 24)) }, split=2) 

## Our functional phenotype data from LUHMES 
my_pheno <- read.csv('My_Pheno.csv', head=T)
row.names(my_pheno) <- my_pheno$Gene
my_genes <- c('ADNP','ARID1B','ASH1L','DYRK1A','PTEN','CHD2','CHD8')
my_data <- my_pheno[my_genes, c(-1,-2)]

my_d4 <- as.dendrogram(hclust(dist(scale(my_pheno[my_genes,c(-1,-2)]))))
my_d4 <- reorder(my_d4, c(5,50,-10,-10,-10,0))

Heatmap(t(scale(my_pheno[my_genes,c(-1,-2)])), column_dend_reorder = F, cluster_columns = my_d4, col=myPal2, cell_fun = function(j, i, x, y, width, height, fill) {
    grid.text(sprintf("%.2f", my_data[j, i]), x, y, gp = gpar(fontsize = 24)) }, split=2) 

```



``` {r, Ranked sgRNA curve plot, include = TRUE}

myDF <- as.data.frame(log2(sort(summary(SFARI_ALL_with_guide$sgRNA), decreasing=T)+1))
sample <- row.names(myDF)
colnames(myDF) <- 'Cells'

plotThis <- myDF

ggplot(plotThis) + geom_point(aes(x=reorder(sample, Cells, rev=T), y=Cells)) + theme_grey(base_size=24) + coord_flip() + geom_vline(xintercept = 4.3, col='red')
```


``` {r, Interpret MIMOSCA results, include = T }
GeneBeta <- read.csv('MIMOSCA_BETA_genelevel.csv', head=T) ## output from MIMOSCA pipeline merging 3 gRNAs at the target gene level
GB2 <- GeneBeta[,-c(12,15)]
row.names(GB2) <- GB2$X
GB3 <- GB2[-1]
#  data.matrix(GB3) 
GB4 <- matrix(unlist(GB3), nrow=13)
row.names(GB4) <- row.names(GB3)
colnames(GB4) <- colnames(GB3)

Heatmap(scale(GB4), cluster_rows = F, cluster_columns = F)

GB5 <- gather(GB3, key = TargetGene, value = BetaCoef)
GB5$sgRNA <- row.names(GB3)

#ggplot2::qplot(sgRNA, BetaCoef, data = GB5, geom='beeswarm')
GB5$Color <- GB5$TargetGene == GB5$sgRNA
ggplot(GB5, aes(x=sgRNA, y=BetaCoef)) + geom_beeswarm(aes(col=Color)) + scale_color_manual(values = c('black','red'))

## Alternatively, plot Z score 
ggplot(GB5, aes(x=sgRNA, y=scale(BetaCoef))) + geom_beeswarm(aes(col=Color)) + scale_color_manual(values = c('black','red'))

###  By TargetGene 
ggplot(GB5, aes(x=TargetGene, y=scale(BetaCoef))) + geom_beeswarm(aes(col=Color, size=5)) + scale_color_manual(values = c('black','red')) + geom_text(aes(y=scale(BetaCoef), label=sgRNA)) + theme_grey(base_size = 30)


### By sgRNA - for paper Fig 2D
ggplot(GB5, aes(x=sgRNA, y=scale(BetaCoef))) + geom_beeswarm(aes(col=Color, size=5)) + scale_color_manual(values = c('black','red')) + theme_grey(base_size = 30)

GB5$Color <- factor(GB5$Color, levels = c("TRUE","FALSE"))
ggboxplot(GB5, x = 'Color', y = 'BetaCoef',  add = "jitter", palette = c('red','black'), color='Color', width = 0.3, add.params = list(alpha=.5)) + stat_compare_means(method="t.test")  

## t-test and Grubb's outlier test for each gene: 

GOI <- 'SETD5'
a <- subset(GB5, GB5$TargetGene == GOI & GB5$sgRNA != GOI)$BetaCoef
b <- subset(GB5, GB5$TargetGene == GOI & GB5$sgRNA == GOI)$BetaCoef

Z <- (mean(a)-b)/sd(a)
Z
pvalue <- pt(-Z,1)
pvalue 

grubbs.test(subset(GB5, GB5$TargetGene == GOI)$BetaCoef, type=10)


GuideBeta <- read.csv('MIMOSCA_sgRNA_rep1and2.csv', head=T) ### MIMOSCA beta estimates for inividual gRNAs across replicate experiments
ggscatter(GuideBeta, x="Rep1", y="Rep2", label="X",    label.select = list(criteria = "`x` > 0 & `y` > 0")) + geom_hline(yintercept = 0) + geom_vline(xintercept = 0) + stat_cor(method="pearson")

ggscatter(GuideBeta, x="Rep1", y="Rep2", label="X") + geom_hline(yintercept = 0) + geom_vline(xintercept = 0) + stat_cor(method="pearson")


AllGuides <- read.csv('~/Downloads/on_target_sgRNA_level.csv', head=T) ### Subset of MIMOSCA results. Expression of only our candidate genes 
ag <- data.matrix(AllGuides)
row.names(ag) <- AllGuides$X
ag2 <- as.data.frame(ag)
ag3 <- gather(ag2[-1], key = sgRNA, value = BetaCoef)
ag3$Gene <- row.names(ag)
ag3$TargetGene <- sapply(strsplit(as.character(ag3$sgRNA), "_"), "[[", 1)

ag3$Color <- ag3$TargetGene == ag3$Gene

#ag3$Color <- factor(ag3$Color, levels = c("TRUE","FALSE"))
ag3$Color <- factor(ag3$Color, levels=c('TRUE','FALSE','OTHER'))
ag3$Color[ag3$TargetGene == 'Nontargeting'] <- 'OTHER'  ## special color for non-targeting sgRNAs


ggplot(ag3, aes(x=sgRNA, y=BetaCoef)) + geom_beeswarm(aes(col=Color)) + scale_color_manual(values = c('red','black','blue'))
ggplot(ag3, aes(x=Gene, y=BetaCoef)) + geom_beeswarm(aes(col=Color)) + scale_color_manual(values = c('red','black', 'blue'))


# ggpboxplot vs geom_boxplot  

ag3$GuideNumber <- sapply(strsplit(as.character(ag3$sgRNA), "_"), "[", 2)

ggboxplot(subset(ag3, !Gene%in% c('none','Nontargeting') & Color !='FALSE'), x='Gene', y='BetaCoef', col='Color') + geom_point(aes(col=Color), position = position_dodge(width=0.8)) + scale_color_manual(values = c('red','black','blue')) + geom_text(data=subset(ag3, Color == T), aes(label=GuideNumber), nudge_x = -.2) + stat_compare_means(aes(group = Color), method='t.test'  ) 


### 
ggplot(subset(ag3, !Gene%in% c('none','Nontargeting') & Color !='FALSE'), aes(x=Gene, y=scale(BetaCoef))) + geom_boxplot(outlier.shape = NA, aes(col=Color)) + geom_point(aes(col=Color), position = position_dodge(width=0.75)) + scale_color_manual(values = c('red','black','blue')) + geom_text(data=subset(ag3, Color==T), aes(label=sgRNA)) + stat_compare_means(aes(group = Color), method='t.test') 

ggplot(subset(ag3, !Gene%in% c('none','Nontargeting') & Color !='FALSE'), aes(x=Gene, y=BetaCoef)) + geom_boxplot(aes(col=Color)) + scale_color_manual(values = c('red','black','blue')) + geom_text(data=subset(ag3, Color==T), aes(label=sgRNA))


ggboxplot(ag3, x = 'Color', y = 'BetaCoef',  add = "jitter", palette = c('red','black','blue'), color='Color', width = 0.3, add.params = list(alpha=.5)) + stat_compare_means() 

ag3$GuideNum <- sapply(strsplit(as.character(ag3$sgRNA), "_"), "[", 2)
ag4 <- subset(ag3, ag3$GuideNum %in% c('G1','G2','G3'))

ggboxplot(ag4, x = 'GuideNum', y = 'BetaCoef',  add = "jitter", palette = c('red','black'), color='Color', width = 0.3, add.params = list(alpha=.5))  + stat_compare_means(aes(group = Color)) #method="t.test")

## guide number and repression => no effect of guide number 
ag5 <- subset(ag4, ag4$Color=="TRUE")
ggboxplot(ag5, x = 'GuideNum', y = 'BetaCoef',  add = "jitter", palette = c('red','black'), color='Color', width = 0.3, add.params = list(alpha=.5))  + stat_compare_means(aes(group = GuideNum))

```

```{r BRB-seq analysis, error=TRUE, include=TRUE, echo=F, warning=F}

BRB_counts <-read.table('BRB_seq_counts', head=F, row.names = 1)

a <- Matrix(unlist(BRB_counts[-1]), nrow=length(row.names(BRB_counts)), sparse=T)

myGroups <- c("ADNP","PTEN","ADNP","DYRK1A","DYRK1A","NonTargeting","CHD8","CHD8","ASH1L","ASH1L","ASH1L", "ARID1B","NonTargeting","PTEN","DYRK1A","CHD2","ARID1B","NonTargeting","ARID1B","CHD2","CHD2","ADNP","PTEN","CHD8")

colnames(a) <- make.unique(as.character(myGroups))
row.names(a) <-  make.unique(as.character(raw_counts[,1]))

Run1_NoUMIS <- CreateSeuratObject(a, min.cells=1, min.features=50)
Run1_NoUMIS <- mito_qc(Run1_NoUMIS)
VlnPlot(object = Run1_NoUMIS, features = c("nCount_RNA", "nFeature_RNA", "percent.mito"))
BRB_Seurat <- Run1_NoUMIS


### Start differential expression analysis of BRB-seq data

ribo_gene.indexes1 <- grep(pattern = "^RPL|^RPS", rownames(BRB_Seurat), ignore.case=TRUE)
BRB_Seurat <- BRB_Seurat[-ribo_gene.indexes1]

raw.counts <- as.data.frame(GetAssayData(BRB_Seurat, slot='counts'))
BRB_DGE <- DGEList(raw.counts, group=myGroups)   ### edgeR analysis of count data 

x <- BRB_DGE  

selr1 <- rowSums(cpm(x$counts)>20)>=2  ###  count > 20 in 2 or more samples
x <- x[selr1,]
x <- calcNormFactors(x)
x$samples$group <- relevel(x$samples$group, ref="NonTargeting") 
design <- model.matrix(~x$samples$group)
xglm <- estimateGLMCommonDisp(x, design, verbose=T)
fit <- glmFit(xglm, design)


### DE for each gene versus Non-targeting gRNA condition 
for (i in 2:8)  { 
lrt <- glmLRT(fit, coef=i)  
deGenes <- topTags(lrt, n=Inf)
deGenes <- deGenes$table
deSig <- subset(deGenes, deGenes$FDR < 0.01 & abs(deGenes$logFC) > .5)
write.table(deSig, file=paste("~/Documents/BRB_seq_DE",as.character(x$samples$group[i])
,sep="_"), quote=F, row.names=T, col.names=T)
}


## Hierarchical clustering of samples by top 2000 variable genes 
top2000 <- VariableFeatures(object = BRB_Seurat)
my_d3 <- as.dendrogram(hclust(dist(t(log(cpm(x[which(row.names(x)%in%top2000),]$counts)+1)))))
my_d3 <- reorder(my_d3, c(rep(10000,12), rep(-20000,6), rep(-10000,6)))
plot(my_d3)

### Plot expression of knock-dwon conditions vs nontargeting cells and off-target gRNAs 
x <- BRB_DGE
selr1 <- rowSums(cpm(x$counts)>5)>=2  ###  count > 5 in 2 or more samples to see expression of lowly expressed genes 
x <- x[selr1,]
x <- calcNormFactors(x)
x$samples$group <- relevel(x$samples$group, ref="NonTargeting") 
design <- model.matrix(~x$samples$group)
xglm <- estimateGLMCommonDisp(x, design, verbose=T)
fit <- glmFit(xglm, design)

geneSet <- c('ADNP','ARID1B','ASH1L','CHD2','CHD8','DYRK1A','PTEN')

plotValues <- ''
nonTarg <- grep('NonTargeting', x$samples$group)

for(i in 1:length(geneSet))
{
GOI <- geneSet[i] 
samples <- grep(GOI, x$samples$group)
sample_Values <- fit$fitted.values[GOI, samples]
names(sample_Values) <- rep(paste(GOI,GOI, sep="_"), length(sample_Values))
nonTarg_Values <- fit$fitted.values[GOI, nonTarg]
names(nonTarg_Values) <-  rep(paste(GOI,'Nontargeting', sep = "_"), length(nonTarg_Values))
offTarg <- fit$fitted.values[GOI,-c(samples, nonTarg)]
names(offTarg) <-  rep(paste(GOI,'Offtarget', sep = "_"), length(offTarg))
newPV <- c(sample_Values, nonTarg_Values, offTarg)
plotValues <- c(plotValues, newPV)
} 

plotDF <- as.data.frame(plotValues[-1])
plotDF$gene <- names(plotValues[-1])
colnames(plotDF)[1] <- 'Expression'

plotDF$TargetGene <- sapply(strsplit(as.character(plotDF$gene), "_"), "[[", 1)
plotDF$Guide <- sapply(strsplit(as.character(plotDF$gene), "_"), "[[", 2)
plotDF$Guide <- factor(plotDF$Guide, levels = c(sort(geneSet), 'Nontargeting','Offtarget'))


ggplot(plotDF, aes(y=as.numeric(paste(Expression)), x=TargetGene)) + geom_boxplot(aes(fill=Guide),outlier.shape=NA, position='dodge') + scale_fill_manual(values = rev(magma(13)[c(-1,-2,-4)]))  + theme_grey(base_size=30) + ylim(0,200) + ylab('Expression (logCPM)')

## Plot PCA of samples after Seurat normalization 

BRB_Seurat <- Run1_NoUMIS
#VlnPlot(object = BRB_Seurat, features = c("nCount_RNA", "nFeature_RNA", "percent.mito"))

BRB_Seurat$gRNA <- sapply(strsplit(names(BRB_Seurat$orig.ident), "[.]"), "[", 1)

BRB_Seurat <- FindVariableFeatures(BRB_Seurat)
BRB_Seurat <- NormalizeData(BRB_Seurat)
BRB_Seurat <- ScaleData(BRB_Seurat) #vars.to.regress = c("Adapter","Lane"))
BRB_Seurat <- RunPCA(BRB_Seurat, npcs=10, verbose=F)
BRB_Seurat <- FindNeighbors(BRB_Seurat, dims=1:10, k=10)
BRB_Seurat <- RunUMAP(BRB_Seurat, reduction="pca",dims=1:10)

Idents(BRB_Seurat) <- BRB_Seurat$gRNA
DimPlot(BRB_Seurat, reduction = "pca", pt.size=5) + scale_color_tableau()

geneSet <- c("PTEN","DYRK1A","CHD8","CHD2","ASH1L","ARID1B","ADNP")
b <- DotPlot(BRB_Seurat, features = geneSet, dot.scale=20, scale.by='size') + scale_color_viridis(option="magma")
b$data$id <- factor(b$data$id, sort(levels(b$data$id)))
b$data$id <- relevel(b$data$id, "NonTargeting")
b + coord_flip() + scale_x_discrete(limits = rev(sort(levels(b$data$id))))

```

```{r, session info}
sessionInfo(package = NULL)
```

