---
title: "Figures | DE + heatmaps/PCA"
author: "Lauren Wolfe"
date: "`r Sys.Date()`"
output: 
  html_document:
    code_folding: hide
    theme: flatly
    highlight: tango
    toc: true
    toc_depth: 2
    toc_float: true
    fig_width: 12
    fig_height: 8
    fig_caption: true
---

## Goals  
  
* Generate figures for RNA seq paper  
* Annotate code  
* Consider figure captions / methods throughout

***
## DE Analyses Included  

**Bulk RNA Seq**  
  
* HH13 Invasive Front vs HH15 Invasive Front  
* HH13 Stream vs HH15 Stream  
* HH13 Invasive Front vs HH13 Stream  
* HH15 Invasive Front vs HH15 Stream  
  
** scRNA-Seq all subpops**  
  
* HH11 Invasive Front vs HH13 Invasive Front  
* HH11 Invasive Front vs HH15 Invasive Front  
* HH11 Invasive Front vs HH13 Leader  
* HH11 Invasive Front vs HH15 Leader  
* HH11 Invasive Front vs HH13 Trailer  
* HH11 Invasive Front vs HH15 Trailer  
  
* HH13 Invasive Front vs HH15 Invasive Front  
* HH13 Invasive Front vs HH13 Leader  
* HH13 Invasive Front vs HH15 Leader  
* HH13 Invasive Front vs HH13 Trailer  
* HH13 Invasive Front vs HH15 Trailer  
  
* HH15 Invasive Front vs HH13 Leader  
* HH15 Invasive Front vs HH15 Leader  
* HH15 Invasive Front vs HH13 Trailer  
* HH15 Invasive Front vs HH15 Trailer  
  
* HH13 Leader vs HH15 Leader  
* HH13 Leader vs HH13 Trailer  
* HH13 Leader vs HH15 Trailer  
  
* HH15 Leader vs HH13 Trailer  

* HH15 Leader vs HH15 Trailer  
  
* HH13 Trailer vs HH15 Trailer  
  
* NT Culture vs HH11 Invasive Front  
  
* NT Culture vs HH13 Invasive Front  
* NT Culture vs HH13 Leader  
* NT Culture vs HH13 Trailer  
  
* NT Culture vs HH15 Invasive Front  
* NT Culture vs HH15 Leader  
* NT Culture vs HH15 Trailer 

** scRNA-seq IF vs St **  
  
* HH15 Invasive Front vs HH15 Stream  
* HH13 Invasive Front vs HH13 Stream  
  
* HH15 and HH13 Invasive Fronts vs HH15 and HH13 Stream  

** Dend clusters **  
  
* Alpha - epsilon vs all other IF  
    
* Alpha - epsilon vs all other in vivo cells  

***

```{r global opts, include = FALSE}
knitr::opts_chunk$set( warning=FALSE, message=FALSE)
```
## Packages Used
```{r setup, message=FALSE}
library(RColorBrewer)
library(NMF)
library(ggplot2)
library(ggbiplot)
library(edgeR)
library(pheatmap)
library(VennDiagram)
library(wesanderson)
library(reshape2)
library(dplyr)
```

## Set Color Pallet
Keep colors the same throughout analysis
```{r}
col1 <-brewer.pal(8,"Dark2")
col2 <- brewer.pal(9, "Set1")
cols <- c(col1[c(4,5,6)], col2[c(4,1,2,7,9)])

set.seed(002) # just to make it reproducible
cols <- sample(cols)

IF11_col <- cols[1]
IF13_col <- cols[2]
L13_col <- cols[3]
Tr13_col <- cols[4]
IF15_col <- cols[5]
L15_col <- cols[6]
Tr15_col <- cols[7]
NT_col <- cols[8]

group_col <- cols
names(group_col) <- c("HH11_Invasive_Front", "HH13_Invasive_Front", "HH13_Leader", "HH13_Trailer", "HH15_Invasive_Front", "HH15_Leader", "HH15_Trailer", "NT_culture")
time_col <- brewer.pal(11, "Spectral")
time_col <- time_col[c(1,3,9,11)]
names(time_col) <- c("HH11", "HH13", "HH15", "NT_culture")
celltype_col <- brewer.pal(5,"Set2")
celltype_col <- celltype_col[2:6]
names(celltype_col) <- c("Invasive_Front", "Leader", "Trailer", "NT_culture")

ann_color <- list(group_col, time_col, celltype_col)

time.iv <- c(FALSE, TRUE, FALSE)
group.iv <- c(TRUE, FALSE, FALSE)

```


## Global Data Used
```{r data, warning=FALSE, message=FALSE}
# genome data (Used to match gene names and ensemble IDs)
genedata<-read.table("/n/data1/genomes/bowtie-index/galGal4/Ens_80/galGal4.Ens_80.genedata.txt",sep='\t',header=T,quote="",comment.char="")
```

## Figure 1 | Unique IF Signature as Determined by Bulk RNA-seq

**Data from bulk seq**  
  
- Differentially expressed == abs(lfc) > 2, padj < 0.05  
  
```{r Figure 1 prep}

# DE script can be found in /n/core/Bioinformatics/analysis/Kulesa/jason_morrison/cbio.jam.102/code/edgeR.DE.R
# Analysis done by ....
# Read in Bulk RNA seq data

Stream_13_15_Bulk_DE <- 
  read.table("/n/projects/lwo/singlecell/txt_files/bulk_rnaseq_DE/Stream_13_15_Bulk.txt", 
             header = T, sep = "\t", quote = "")
Invasive_Front_13_15_Bulk_DE <- 
  read.table("/n/projects/lwo/singlecell/txt_files/bulk_rnaseq_DE/Trailblazer_13_15_Bulk.txt", 
             header = T, sep = "\t", quote = "")
Invasive_Front_Stream_15_Bulk_DE <-
  read.table("/n/projects/lwo/singlecell/txt_files/bulk_rnaseq_DE/Trailblazer_Stream_15_Bulk.txt",
             header = T, sep = "\t", quote = "")
Invasive_Front_Stream_13_Bulk_DE <-
  read.table("/n/projects/lwo/singlecell/txt_files/bulk_rnaseq_DE/Trailblazer_Stream_13_Bulk.txt", 
             header = T, sep = "\t", quote = "")


# change Trailblazer to invasive front
names(edger.rpkm.bulk) <- gsub("Trailblazer_", "Invasive_Front_", names(edger.rpkm.bulk))

# DE from bulk RNA seq
# subset out DE genes
Invasive_Front_Stream_15_sigs<-subset(Invasive_Front_Stream_15_Bulk_DE, FDR<=0.05 & abs(logFC) > 2) # 91
Invasive_Front_Stream_13_sigs<-subset(Invasive_Front_Stream_13_Bulk_DE, FDR<=0.05 & abs(logFC) > 2)  # 316
Stream_13_15_sigs<-subset(Stream_13_15_Bulk_DE, FDR<=0.05 & abs(logFC) > 2)  # 1
Invasive_Front_13_15_sigs<-subset(Invasive_Front_13_15_Bulk_DE, FDR<=0.05 & abs(logFC) > 2)  # 283

# merge gene list
bulk_DE_genelist <- union(Invasive_Front_Stream_13_sigs$geneId, Invasive_Front_Stream_15_sigs$geneId) # 365
bulk_DE_genelist <- union(bulk_DE_genelist, Stream_13_15_sigs$geneId) # 365
bulk_DE_genelist <- union(bulk_DE_genelist, Invasive_Front_13_15_sigs$geneId) # Total == 477


```

## Heatmap of 477 DE genes

```{r Figure 1B}

# logical vector to subsets sigs list
genes.iv <- edger.rpkm.bulk$geneId %in% bulk_DE_genelist

# get rpkm 
temp <- edger.rpkm.bulk[genes.iv,]
# gene symbol to ensid
row.names(temp) <- paste(temp$geneId, "  ", temp$Name)

temp <- temp[,2:13]
temp<-log2(temp+.1)
temp<-scale(t(temp), center = T, scale= T) # z-score
temp <- data.frame(t(temp))
temp[temp>3]<-3
temp[temp< -3]<- -3
df <- temp[complete.cases(temp),]

# draw heatmap

#pdf("/n/projects/lwo/singlecell/R_out/Finalized_figures/bulkseq_heatmap_pearson.pdf", width = 12, height = 15)
aheatmap(df, color = rev(brewer.pal(11,"RdBu")), scale = "none", breaks = 0, treeheight = 35, distfun = "pearson")
#dev.off()

```

## Averaged Bulk Replicates

```{r}
Stream_13.iv <- grepl("Stream_13", names(edger.rpkm.bulk))
Stream_15.iv <- grepl("Stream_15", names(edger.rpkm.bulk))
Invasive_Front_15.iv <- grepl("Invasive_Front_15", names(edger.rpkm.bulk))
Invasive_Front_13.iv <- grepl("Invasive_Front_13", names(edger.rpkm.bulk))

temp <- edger.rpkm.bulk[genes.iv,]

means <- data.frame(ensid = temp$geneId, 
                    st_13_mean = rowMeans(temp[,Stream_13.iv]),
                    st_15_mean = rowMeans(temp[,Stream_15.iv]),
                    if_13_mean = rowMeans(temp[,Invasive_Front_13.iv]),
                    if_15_mean = rowMeans(temp[,Invasive_Front_15.iv]))

row.names(temp) <- paste(temp$geneId, "  ", temp$Name)

temp <- means[,-1]
temp<-log2(temp+.1)
temp<-scale(t(temp), center = T, scale= T) # z-score
temp <- data.frame(t(temp))
temp[temp>3]<-3
temp[temp< -3]<- -3
df <- temp[complete.cases(temp),]

#pdf("/n/projects/lwo/singlecell/R_out/Finalized_figures/bulkseq_heatmap_pearson_averaged.pdf", width = 12, height = 15)
aheatmap(df, color = rev(brewer.pal(11,"RdBu")), scale = "none", breaks = 0, treeheight = 35, distfun = "pearson")
#dev.off()

```


## Data for single-cell analysis
```{r}

# not subsetted by QC yet
ucount_df <- read.table('/n/core/Bioinformatics/analysis/Kulesa/jason_morrison/cbio.jam.104/data/ucounts.txt', header = T, sep = '\t')
targets <- read.table('/n/projects/lwo/singlecell/OLD_R/Finalized_txt/targets.sub.txt', header = T, sep = '\t')
edger.rpkm.sc <- read.table("/n/projects/lwo/singlecell/OLD_R/Finalized_txt/edger.rpkm.sub.txt", header = T, sep = "\t")
## change TB to IF
targets$alias <- gsub("Trailblazer", "Invasive_Front",targets$alias)
targets$group <- gsub("Trailblazer", "Invasive_Front",targets$group)

# subset ucount_df
# remove low quality cells
names(ucount_df) <- gsub("Trailblazer", "Invasive_Front", names(ucount_df))
ensid <- ucount_df[,1]
ucount_df <- ucount_df[,names(ucount_df) %in% targets$alias] 
ucount_df <- data.frame(ensid, ucount_df)

```


## DEG lists for all single-cell subpopulation comparisons  
  

```{r Differential Expression Single Cell}

# ucount_sub == remove ensid
ucount_sub <- ucount_df[2:length(ucount_df)]
y <- DGEList(counts=ucount_sub, group=targets$group)
keep<-apply(cpm(y),1,function(x) length(x[x>5])>=2)
y <- y[keep, , keep.lib.sizes=FALSE]
y <- calcNormFactors(y)
y <- estimateCommonDisp(y)
y <- estimateTagwiseDisp(y)

# Note that the first group listed in the pair is the baseline for the comparison—so if the pair is c("A","B") then the comparison is B - A, so genes with positive log-fold change are up-regulated in group B compared with group A (and vice versa for genes with negative log-fold change).

res<-list()
res[[1]] <- exactTest(y, pair=c("HH13_Invasive_Front", "HH11_Invasive_Front"))
res[[2]] <- exactTest(y, pair = c("HH15_Invasive_Front", "HH11_Invasive_Front"))
res[[3]] <- exactTest(y, pair = c("HH13_Leader", "HH11_Invasive_Front"))
res[[4]] <- exactTest(y, pair =c("HH15_Leader", "HH11_Invasive_Front"))
res[[5]] <- exactTest(y, pair = c("HH13_Trailer", "HH11_Invasive_Front"))
res[[6]] <- exactTest(y, pair = c("HH15_Trailer", "HH11_Invasive_Front"))

res[[7]] <- exactTest(y, pair=c("HH15_Invasive_Front", "HH13_Invasive_Front"))
res[[8]] <- exactTest(y, pair = c("HH13_Leader", "HH13_Invasive_Front"))
res[[9]] <- exactTest(y, pair =c("HH15_Leader", "HH13_Invasive_Front"))
res[[10]]<- exactTest(y, pair = c("HH13_Trailer", "HH13_Invasive_Front"))
res[[11]] <- exactTest(y, pair = c("HH15_Trailer", "HH13_Invasive_Front"))

res[[12]] <- exactTest(y, pair = c("HH13_Leader", "HH15_Invasive_Front"))
res[[13]] <- exactTest(y, pair =c("HH15_Leader", "HH15_Invasive_Front"))
res[[14]]<- exactTest(y, pair = c("HH13_Trailer", "HH15_Invasive_Front"))
res[[15]] <- exactTest(y, pair = c("HH15_Trailer", "HH15_Invasive_Front"))

res[[16]] <- exactTest(y, pair =c("HH15_Leader", "HH13_Leader"))
res[[17]]<- exactTest(y, pair = c("HH13_Trailer", "HH13_Leader"))
res[[18]] <- exactTest(y, pair = c("HH15_Trailer", "HH13_Leader"))

res[[19]]<- exactTest(y, pair = c("HH13_Trailer", "HH15_Leader"))
res[[20]] <- exactTest(y, pair = c("HH15_Trailer", "HH15_Leader"))

res[[21]] <- exactTest(y, pair = c("HH15_Trailer", "HH13_Trailer"))


# List of conts
# X_v_Y == X relative to Y
conts <- c("HH11_Invasive_Front_v_HH13_Invasive_Front", 
           "HH11_Invasive_Front_v_HH15_Invasive_Front", 
           "HH11_Invasive_Front_v_HH13_Leader",
           "HH11_Invasive_Front_v_HH15_Leader",
           "HH11_Invasive_Front_v_HH13_Trailer",  
           "HH11_Invasive_Front_v_HH15_Trailer", 
           "HH13_Invasive_Front_v_HH15_Invasive_Front",
           "HH13_Invasive_Front_v_HH13_Leader",  
           "HH13_Invasive_Front_v_HH15_Leader",
           "HH13_Invasive_Front_v_HH13_Trailer",
           "HH13_Invasive_Front_v_HH15_Trailer", 
           "HH15_Invasive_Front_v_HH13_Leader", 
           "HH15_Invasive_Front_v_HH15_Leader",
           "HH15_Invasive_Front_v_HH13_Trailer", 
           "HH15_Invasive_Front_v_HH15_Trailer",
           "HH13_Leader_v_HH15_Leader",
           "HH13_Leader_v_HH13_Trailer", 
           "HH13_Leader_v_HH15_Trailer", 
           "HH15_Leader_v_HH13_Trailer", 
           "HH15_Leader_v_HH15_Trailer", 
           "HH13_Trailer_v_HH15_Trailer")
# Add Padj
for(i in 1:length(conts))
{
  res[[i]]$table$padj<-p.adjust(res[[i]]$table$PValue,method="BH")
}

# Summary table
summarytable <- data.frame(ucount_df[which(keep),1],res[[1]]$table,res[[2]]$table,res[[3]]$table,res[[4]]$table, res[[5]]$table, res[[6]]$table, res[[7]]$table,res[[8]]$table,res[[9]]$table,res[[10]]$table, res[[11]]$table, res[[12]]$table, res[[13]]$table,res[[14]]$table,res[[15]]$table,res[[16]]$table, res[[17]]$table, res[[18]]$table, res[[19]]$table,res[[20]]$table,res[[21]]$table)

# Column lables
colnames(summarytable)[1]<-"ensid"
colnames(summarytable)[2:length(summarytable[1,])]<-paste(c(paste(rep(conts,each=4),rep(c(".log2fc",".log2cpm",".pval", ".padj"),times=1),sep='')))

# add annotation to summary of DE analysis
summarytable_ann_invivoDE<-data.frame(ensid=summarytable[,1],genedata[match(summarytable[,1],genedata[,1]),c("Name", "Chr", "Start", "End", "Strand", "Biotype","Description")],summarytable[,2:length(summarytable[1,])])
rownames(summarytable_ann_invivoDE)<-summarytable_ann_invivoDE$ensid

# write.table(summarytable_ann_invivoDE, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/summary_tables/summarytable_invivoDE.txt", row.names = FALSE, sep = "\t")
```
  
## Make DE tables

```{r make DE tables}
HH11_Invasive_Front_v_HH13_Invasive_Front_DE.up <- subset(summarytable_ann_invivoDE, HH11_Invasive_Front_v_HH13_Invasive_Front.log2fc > 2 & HH11_Invasive_Front_v_HH13_Invasive_Front.padj < 0.05)

HH11_Invasive_Front_v_HH13_Invasive_Front_DE.dn <- subset(summarytable_ann_invivoDE, HH11_Invasive_Front_v_HH13_Invasive_Front.log2fc < -2 & HH11_Invasive_Front_v_HH13_Invasive_Front.padj < 0.05)

HH11_Invasive_Front_v_HH15_Invasive_Front_DE.up <- subset(summarytable_ann_invivoDE, HH11_Invasive_Front_v_HH15_Invasive_Front.log2fc > 2 & HH11_Invasive_Front_v_HH15_Invasive_Front.padj < 0.05)

HH11_Invasive_Front_v_HH15_Invasive_Front_DE.dn <- subset(summarytable_ann_invivoDE, HH11_Invasive_Front_v_HH15_Invasive_Front.log2fc < -2 & HH11_Invasive_Front_v_HH15_Invasive_Front.padj < 0.05)

HH11_Invasive_Front_v_HH13_Leader_DE.up <- subset(summarytable_ann_invivoDE, HH11_Invasive_Front_v_HH13_Leader.log2fc > 2 & HH11_Invasive_Front_v_HH13_Leader.padj < 0.05)

HH11_Invasive_Front_v_HH13_Leader_DE.dn <- subset(summarytable_ann_invivoDE, HH11_Invasive_Front_v_HH13_Leader.log2fc < -2 & HH11_Invasive_Front_v_HH13_Leader.padj < 0.05)

HH11_Invasive_Front_v_HH15_Leader_DE.up <- subset(summarytable_ann_invivoDE, HH11_Invasive_Front_v_HH15_Leader.log2fc > 2 & HH11_Invasive_Front_v_HH15_Leader.padj < 0.05)

HH11_Invasive_Front_v_HH15_Leader_DE.dn <- subset(summarytable_ann_invivoDE, HH11_Invasive_Front_v_HH15_Leader.log2fc < -2 & HH11_Invasive_Front_v_HH15_Leader.padj < 0.05)

HH11_Invasive_Front_v_HH13_Trailer_DE.up <- subset(summarytable_ann_invivoDE, HH11_Invasive_Front_v_HH13_Trailer.log2fc > 2 & HH11_Invasive_Front_v_HH13_Trailer.padj < 0.05)

HH11_Invasive_Front_v_HH13_Trailer_DE.dn <- subset(summarytable_ann_invivoDE, HH11_Invasive_Front_v_HH13_Trailer.log2fc < -2 & HH11_Invasive_Front_v_HH13_Trailer.padj < 0.05)

HH11_Invasive_Front_v_HH15_Trailer_DE.up <- subset(summarytable_ann_invivoDE, HH11_Invasive_Front_v_HH15_Trailer.log2fc > 2 & HH11_Invasive_Front_v_HH15_Trailer.padj < 0.05)

HH11_Invasive_Front_v_HH15_Trailer_DE.dn <- subset(summarytable_ann_invivoDE, HH11_Invasive_Front_v_HH15_Trailer.log2fc < -2 & HH11_Invasive_Front_v_HH15_Trailer.padj < 0.05)

HH13_Invasive_Front_v_HH15_Invasive_Front_DE.up <- subset(summarytable_ann_invivoDE, HH13_Invasive_Front_v_HH15_Invasive_Front.log2fc > 2 & HH13_Invasive_Front_v_HH15_Invasive_Front.padj < 0.05)

HH13_Invasive_Front_v_HH15_Invasive_Front_DE.dn <- subset(summarytable_ann_invivoDE, HH13_Invasive_Front_v_HH15_Invasive_Front.log2fc < -2 & HH13_Invasive_Front_v_HH15_Invasive_Front.padj < 0.05)

HH13_Invasive_Front_v_HH13_Leader_DE.up <- subset(summarytable_ann_invivoDE, HH13_Invasive_Front_v_HH13_Leader.log2fc > 2 & HH13_Invasive_Front_v_HH13_Leader.padj < 0.05)

HH13_Invasive_Front_v_HH13_Leader_DE.dn <- subset(summarytable_ann_invivoDE, HH13_Invasive_Front_v_HH13_Leader.log2fc < -2 & HH13_Invasive_Front_v_HH13_Leader.padj < 0.05)

HH13_Invasive_Front_v_HH15_Leader_DE.up <- subset(summarytable_ann_invivoDE, HH13_Invasive_Front_v_HH15_Leader.log2fc > 2 & HH13_Invasive_Front_v_HH15_Leader.padj < 0.05)

HH13_Invasive_Front_v_HH15_Leader_DE.dn <- subset(summarytable_ann_invivoDE, HH13_Invasive_Front_v_HH15_Leader.log2fc < -2 & HH13_Invasive_Front_v_HH15_Leader.padj < 0.05)

HH13_Invasive_Front_v_HH13_Trailer_DE.up <- subset(summarytable_ann_invivoDE, HH13_Invasive_Front_v_HH13_Trailer.log2fc > 2 & HH13_Invasive_Front_v_HH13_Trailer.padj < 0.05)

HH13_Invasive_Front_v_HH13_Trailer_DE.dn <- subset(summarytable_ann_invivoDE, HH13_Invasive_Front_v_HH13_Trailer.log2fc < -2 & HH13_Invasive_Front_v_HH13_Trailer.padj < 0.05)

HH13_Invasive_Front_v_HH15_Trailer_DE.up <- subset(summarytable_ann_invivoDE, HH13_Invasive_Front_v_HH15_Trailer.log2fc > 2 & HH13_Invasive_Front_v_HH15_Trailer.padj < 0.05)

HH13_Invasive_Front_v_HH15_Trailer_DE.dn <- subset(summarytable_ann_invivoDE, HH13_Invasive_Front_v_HH15_Trailer.log2fc < -2 & HH13_Invasive_Front_v_HH15_Trailer.padj < 0.05)

HH15_Invasive_Front_v_HH13_Leader_DE.up <- subset(summarytable_ann_invivoDE, HH15_Invasive_Front_v_HH13_Leader.log2fc > 2 & HH15_Invasive_Front_v_HH13_Leader.padj < 0.05)

HH15_Invasive_Front_v_HH13_Leader_DE.dn <- subset(summarytable_ann_invivoDE, HH15_Invasive_Front_v_HH13_Leader.log2fc < -2 & HH15_Invasive_Front_v_HH13_Leader.padj < 0.05)
           
HH15_Invasive_Front_v_HH15_Leader_DE.up <- subset(summarytable_ann_invivoDE, HH15_Invasive_Front_v_HH15_Leader.log2fc > 2 & HH15_Invasive_Front_v_HH15_Leader.padj < 0.05)

HH15_Invasive_Front_v_HH15_Leader_DE.dn <- subset(summarytable_ann_invivoDE, HH15_Invasive_Front_v_HH15_Leader.log2fc < -2 & HH15_Invasive_Front_v_HH15_Leader.padj < 0.05)

HH15_Invasive_Front_v_HH13_Trailer_DE.up <- subset(summarytable_ann_invivoDE,HH15_Invasive_Front_v_HH13_Trailer.log2fc > 2 & HH15_Invasive_Front_v_HH13_Trailer.padj < 0.05)

HH15_Invasive_Front_v_HH13_Trailer_DE.dn <- subset(summarytable_ann_invivoDE,HH15_Invasive_Front_v_HH13_Trailer.log2fc < -2 & HH15_Invasive_Front_v_HH13_Trailer.padj < 0.05)

HH15_Invasive_Front_v_HH15_Trailer_DE.up <- subset(summarytable_ann_invivoDE,HH15_Invasive_Front_v_HH15_Trailer.log2fc > 2 & HH15_Invasive_Front_v_HH15_Trailer.padj < 0.05)

HH15_Invasive_Front_v_HH15_Trailer_DE.dn <- subset(summarytable_ann_invivoDE,HH15_Invasive_Front_v_HH15_Trailer.log2fc < -2 & HH15_Invasive_Front_v_HH15_Trailer.padj < 0.05)

HH13_Leader_v_HH15_Leader_DE.up <- subset(summarytable_ann_invivoDE, HH13_Leader_v_HH15_Leader.log2fc > 2 & HH13_Leader_v_HH15_Leader.padj < 0.05)

HH13_Leader_v_HH15_Leader_DE.dn <- subset(summarytable_ann_invivoDE, HH13_Leader_v_HH15_Leader.log2fc < -2 & HH13_Leader_v_HH15_Leader.padj < 0.05)

HH13_Leader_v_HH13_Trailer_DE.up <- subset(summarytable_ann_invivoDE, HH13_Leader_v_HH13_Trailer.log2fc > 2 & HH13_Leader_v_HH13_Trailer.padj < 0.05)

HH13_Leader_v_HH13_Trailer_DE.dn <- subset(summarytable_ann_invivoDE, HH13_Leader_v_HH13_Trailer.log2fc < -2 & HH13_Leader_v_HH13_Trailer.padj < 0.05)

HH13_Leader_v_HH15_Trailer_DE.up <- subset(summarytable_ann_invivoDE, HH13_Leader_v_HH15_Trailer.log2fc > 2 & HH13_Leader_v_HH15_Trailer.padj < 0.05)

HH13_Leader_v_HH15_Trailer_DE.dn <- subset(summarytable_ann_invivoDE, HH13_Leader_v_HH15_Trailer.log2fc < -2 & HH13_Leader_v_HH15_Trailer.padj < 0.05)

HH15_Leader_v_HH13_Trailer_DE.up <- subset(summarytable_ann_invivoDE, HH15_Leader_v_HH13_Trailer.log2fc > 2 & HH15_Leader_v_HH13_Trailer.padj < 0.05)

HH15_Leader_v_HH13_Trailer_DE.dn <- subset(summarytable_ann_invivoDE, HH15_Leader_v_HH13_Trailer.log2fc < -2 & HH15_Leader_v_HH13_Trailer.padj < 0.05)

HH15_Leader_v_HH15_Trailer_DE.up <- subset(summarytable_ann_invivoDE, HH15_Leader_v_HH15_Trailer.log2fc > 2 & HH15_Leader_v_HH15_Trailer.padj < 0.05)

HH15_Leader_v_HH15_Trailer_DE.dn <- subset(summarytable_ann_invivoDE, HH15_Leader_v_HH15_Trailer.log2fc < -2 & HH15_Leader_v_HH15_Trailer.padj < 0.05)

HH13_Trailer_v_HH15_Trailer_DE.up <- subset(summarytable_ann_invivoDE, HH13_Trailer_v_HH15_Trailer.log2fc > 2 & HH13_Trailer_v_HH15_Trailer.padj < 0.05)

HH13_Trailer_v_HH15_Trailer_DE.dn <- subset(summarytable_ann_invivoDE, HH13_Trailer_v_HH15_Trailer.log2fc < -2 & HH13_Trailer_v_HH15_Trailer.padj < 0.05)

```

## Write DEG list

```{r, eval = FALSE}
# write.table(HH11_Invasive_Front_v_HH13_Invasive_Front_DE.up, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH11_Invasive_Front_v_HH13_Invasive_Front_DE.up.txt", row.names = FALSE, sep = "\t")
# write.table(HH11_Invasive_Front_v_HH13_Invasive_Front_DE.dn, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH11_Invasive_Front_v_HH13_Invasive_Front_DE.dn.txt", row.names = FALSE, sep = "\t")
# write.table(HH11_Invasive_Front_v_HH15_Invasive_Front_DE.up, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH11_Invasive_Front_v_HH15_Invasive_Front_DE.up.txt", row.names = FALSE, sep = "\t")
# write.table(HH11_Invasive_Front_v_HH15_Invasive_Front_DE.dn, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH11_Invasive_Front_v_HH15_Invasive_Front_DE.dn.txt", row.names = FALSE, sep = "\t")
# write.table(HH11_Invasive_Front_v_HH13_Leader_DE.up, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH11_Invasive_Front_v_HH13_Leader_DE.up.txt", row.names = FALSE, sep = "\t")
# write.table(HH11_Invasive_Front_v_HH13_Leader_DE.dn, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH11_Invasive_Front_v_HH13_Leader_DE.dn.txt", row.names = FALSE, sep = "\t")
# write.table(HH11_Invasive_Front_v_HH15_Leader_DE.up, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH11_Invasive_Front_v_HH15_Leader_DE.up.txt", row.names = FALSE, sep = "\t")
# write.table(HH11_Invasive_Front_v_HH15_Leader_DE.dn, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH11_Invasive_Front_v_HH15_Leader_DE.dn.txt", row.names = FALSE, sep = "\t")
# write.table(HH11_Invasive_Front_v_HH13_Trailer_DE.up, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH11_Invasive_Front_v_HH13_Trailer_DE.up.txt", row.names = FALSE, sep = "\t")
# write.table(HH11_Invasive_Front_v_HH13_Trailer_DE.dn, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH11_Invasive_Front_v_HH13_Trailer_DE.dn.txt", row.names = FALSE, sep = "\t")
# write.table(HH11_Invasive_Front_v_HH15_Trailer_DE.up, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH11_Invasive_Front_v_HH15_Trailer_DE.up.txt", row.names = FALSE, sep = "\t")
# write.table(HH11_Invasive_Front_v_HH15_Trailer_DE.dn, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH11_Invasive_Front_v_HH15_Trailer_DE.dn.txt", row.names = FALSE, sep = "\t")
# write.table(HH13_Invasive_Front_v_HH15_Invasive_Front_DE.up, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH13_Invasive_Front_v_HH15_Invasive_Front_DE.up.txt", row.names = FALSE, sep = "\t")
# write.table(HH13_Invasive_Front_v_HH15_Invasive_Front_DE.dn, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH13_Invasive_Front_v_HH15_Invasive_Front_DE.dn.txt", row.names = FALSE, sep = "\t")
# write.table(HH13_Invasive_Front_v_HH13_Leader_DE.up, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH13_Invasive_Front_v_HH13_Leader_DE.up.txt", row.names = FALSE, sep = "\t")
# write.table(HH13_Invasive_Front_v_HH13_Leader_DE.dn, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH13_Invasive_Front_v_HH13_Leader_DE.dn.txt", row.names = FALSE, sep = "\t")
# write.table(HH13_Invasive_Front_v_HH15_Leader_DE.up, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH13_Invasive_Front_v_HH15_Leader_DE.up.txt", row.names = FALSE, sep = "\t")
# write.table(HH13_Invasive_Front_v_HH15_Leader_DE.dn, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH13_Invasive_Front_v_HH15_Leader_DE.dn.txt", row.names = FALSE, sep = "\t")
# write.table(HH13_Invasive_Front_v_HH13_Trailer_DE.up, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH13_Invasive_Front_v_HH13_Trailer_DE.up.txt", row.names = FALSE, sep = "\t")
# write.table(HH13_Invasive_Front_v_HH13_Trailer_DE.dn, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH13_Invasive_Front_v_HH13_Trailer_DE.dn.txt", row.names = FALSE, sep = "\t")
# write.table(HH13_Invasive_Front_v_HH15_Trailer_DE.up, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH13_Invasive_Front_v_HH15_Trailer_DE.up.txt", row.names = FALSE, sep = "\t")
# write.table(HH13_Invasive_Front_v_HH15_Trailer_DE.dn, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH13_Invasive_Front_v_HH15_Trailer_DE.dn.txt", row.names = FALSE, sep = "\t")
# write.table(HH15_Invasive_Front_v_HH13_Leader_DE.up, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH15_Invasive_Front_v_HH13_Leader_DE.up.txt", row.names = FALSE, sep = "\t")
# write.table(HH15_Invasive_Front_v_HH13_Leader_DE.dn, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH15_Invasive_Front_v_HH13_Leader_DE.dn.txt", row.names = FALSE, sep = "\t")
# write.table(HH15_Invasive_Front_v_HH15_Leader_DE.up, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH15_Invasive_Front_v_HH15_Leader_DE.up.txt", row.names = FALSE, sep = "\t")
# write.table(HH15_Invasive_Front_v_HH15_Leader_DE.dn, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH15_Invasive_Front_v_HH15_Leader_DE.dn.txt", row.names = FALSE, sep = "\t")
# write.table(HH15_Invasive_Front_v_HH13_Trailer_DE.up, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH15_Invasive_Front_v_HH13_Trailer_DE.up.txt", row.names = FALSE, sep = "\t")
# write.table(HH15_Invasive_Front_v_HH13_Trailer_DE.dn, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH15_Invasive_Front_v_HH13_Trailer_DE.dn.txt", row.names = FALSE, sep = "\t")
# write.table(HH15_Invasive_Front_v_HH15_Trailer_DE.up, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH15_Invasive_Front_v_HH15_Trailer_DE.up.txt", row.names = FALSE, sep = "\t")
# write.table(HH15_Invasive_Front_v_HH15_Trailer_DE.dn, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH15_Invasive_Front_v_HH15_Trailer_DE.dn.txt", row.names = FALSE, sep = "\t")
# write.table(HH13_Leader_v_HH15_Leader_DE.up, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH13_Leader_v_HH15_Leader_DE.up.txt", row.names = FALSE, sep = "\t")
# write.table(HH13_Leader_v_HH15_Leader_DE.dn, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH13_Leader_v_HH15_Leader_DE.dn.txt", row.names = FALSE, sep = "\t")
# write.table(HH13_Leader_v_HH13_Trailer_DE.up, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH13_Leader_v_HH13_Trailer_DE.up.txt", row.names = FALSE, sep = "\t")
# write.table(HH13_Leader_v_HH13_Trailer_DE.dn, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH13_Leader_v_HH13_Trailer_DE.dn.txt", row.names = FALSE, sep = "\t")
# write.table(HH13_Leader_v_HH15_Trailer_DE.up, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH13_Leader_v_HH15_Trailer_DE.up.txt", row.names = FALSE, sep = "\t")
# write.table(HH13_Leader_v_HH15_Trailer_DE.dn, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH13_Leader_v_HH15_Trailer_DE.dn.txt", row.names = FALSE, sep = "\t")
# write.table(HH15_Leader_v_HH13_Trailer_DE.up, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH15_Leader_v_HH13_Trailer_DE.up.txt", row.names = FALSE, sep = "\t")
# write.table(HH15_Leader_v_HH13_Trailer_DE.dn, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH15_Leader_v_HH13_Trailer_DE.dn.txt", row.names = FALSE, sep = "\t")
# write.table(HH15_Leader_v_HH15_Trailer_DE.up, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH15_Leader_v_HH15_Trailer_DE.up.txt", row.names = FALSE, sep = "\t")
# write.table(HH15_Leader_v_HH15_Trailer_DE.dn, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH15_Leader_v_HH15_Trailer_DE.dn.txt", row.names = FALSE, sep = "\t")
# write.table(HH13_Trailer_v_HH15_Trailer_DE.up, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH13_Trailer_v_HH15_Trailer_DE.up.txt", row.names = FALSE, sep = "\t")
# write.table(HH13_Trailer_v_HH15_Trailer_DE.dn, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/HH13_Trailer_v_HH15_Trailer_DE.dn.txt", row.names = FALSE, sep = "\t")
```

## Compile genes from all in vivo DE
  
## Collect unique genes and subset rpkm table
```{r unique genes?}

# Capture 42 data frames in a list
df_list <- lapply(ls(pattern="_DE.up|_DE.dn"), function(x) get(x))
names(df_list) <- seq(from =1, to =42, by = 1)

# Collect ensids
gene_list <- vector()
for(i in 1:length(df_list))
{
  gene_list <- c(gene_list, as.character(df_list[[i]][[1]]))
}

# Get rid of repeats
sc_invivo_DE_genelist <- unique(gene_list) 


```

## Get DE genes for NT_culture

```{r}

y <- DGEList(counts=ucount_df[,-1], group=targets$group)
keep<-apply(cpm(y),1,function(x) length(x[x>5])>=2)
y <- y[keep, , keep.lib.sizes=FALSE]
y <- calcNormFactors(y)
y <- estimateCommonDisp(y)
y <- estimateTagwiseDisp(y)

# Note that the first group listed in the pair is the baseline for the comparison—so if the pair is c("A","B") then the comparison is B - A, so genes with positive log-fold change are up-regulated in group B compared with group A (and vice versa for genes with negative log-fold change).

res<-list()
res[[1]] <- exactTest(y, pair=c("HH11_Invasive_Front", "NT_culture"))
res[[2]] <- exactTest(y, pair = c("HH13_Invasive_Front", "NT_culture"))
res[[3]] <- exactTest(y, pair = c("HH15_Invasive_Front", "NT_culture"))
res[[4]] <- exactTest(y, pair =c("HH13_Trailer", "NT_culture"))
res[[5]] <- exactTest(y, pair = c("HH15_Trailer", "NT_culture"))
res[[6]] <- exactTest(y, pair = c("HH13_Leader", "NT_culture"))
res[[7]] <- exactTest(y, pair = c("HH15_Leader", "NT_culture"))


# List of conts
# X_v_Y == X relative to Y
conts <- c("NT_v_HH11_Invasive_Front", "NT_v_HH13_Invasive_Front", "NT_v_HH15_Invasive_Front", "NT_v_HH13_Trailer", "NT_v_HH15_Trailer", "NT_v_HH13_Leader", "NT_v_HH15_Leader") 

# Add Padj
for(i in 1:length(conts))
{
  res[[i]]$table$padj<-p.adjust(res[[i]]$table$PValue,method="BH")
}

# Summary table
summarytable <- data.frame(ucount_df[which(keep),1],res[[1]]$table,res[[2]]$table,res[[3]]$table,res[[4]]$table,res[[5]]$table, res[[6]]$table,res[[7]]$table)
# Column lables
colnames(summarytable)[1]<-"ensid"
colnames(summarytable)[2:length(summarytable[1,])]<-paste(c(paste(rep(conts,each=4),rep(c(".log2fc",".log2cpm",".pval", ".padj"),times=1),sep='')))

# add annotation to summary of DE analysis
summarytable_ann_invitroDE<-data.frame(ensid=summarytable[,1],genedata[match(summarytable[,1],genedata[,1]),c("Name", "Chr", "Start", "End", "Strand", "Biotype","Description")],summarytable[,2:length(summarytable[1,])])
rownames(summarytable_ann_invitroDE)<-summarytable_ann_invitroDE$ensid

# write.table(summarytable_ann_invitroDE, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/summary_tables/summarytable_invitroDE.txt", row.names = FALSE, sep = "\t")
```

## Make DE tables

```{r}
NT_v_HH11_Invasive_Front <- subset(summarytable_ann_invitroDE, abs(summarytable_ann_invitroDE$NT_v_HH11_Invasive_Front.log2fc) > 2 & summarytable_ann_invitroDE$NT_v_HH11_Invasive_Front.padj < 0.05)

NT_v_HH13_Invasive_Front <- subset(summarytable_ann_invitroDE, abs(summarytable_ann_invitroDE$NT_v_HH13_Invasive_Front.log2fc) > 2 & summarytable_ann_invitroDE$NT_v_HH13_Invasive_Front.padj < 0.05)

NT_v_HH15_Invasive_Front <- subset(summarytable_ann_invitroDE, abs(summarytable_ann_invitroDE$NT_v_HH15_Invasive_Front.log2fc) > 2 & summarytable_ann_invitroDE$NT_v_HH15_Invasive_Front.padj < 0.05)

NT_v_HH13_Trailer <- subset(summarytable_ann_invitroDE, abs(summarytable_ann_invitroDE$NT_v_HH13_Trailer.log2fc) > 2 & summarytable_ann_invitroDE$NT_v_HH13_Trailer.padj < 0.05)

NT_v_HH15_Trailer <- subset(summarytable_ann_invitroDE, abs(summarytable_ann_invitroDE$NT_v_HH15_Trailer.log2fc) > 2 & summarytable_ann_invitroDE$NT_v_HH15_Trailer.padj < 0.05)

NT_v_HH13_Leader <- subset(summarytable_ann_invitroDE, abs(summarytable_ann_invitroDE$NT_v_HH13_Leader.log2fc) > 2 & summarytable_ann_invitroDE$NT_v_HH13_Leader.padj < 0.05)

NT_v_HH15_Leader <- subset(summarytable_ann_invitroDE, abs(summarytable_ann_invitroDE$NT_v_HH15_Leader.log2fc) > 2 & summarytable_ann_invitroDE$NT_v_HH15_Leader.padj < 0.05)


sc_invitro_DE_genelist <- as.vector(c(as.character(NT_v_HH11_Invasive_Front$ensid), as.character(NT_v_HH13_Invasive_Front$ensid), as.character(NT_v_HH15_Invasive_Front$ensid), as.character(NT_v_HH13_Trailer$ensid),  as.character(NT_v_HH15_Trailer$ensid), as.character(NT_v_HH13_Leader$ensid), as.character(NT_v_HH15_Leader$ensid))) #8235

#write.table(NT_v_HH11_Invasive_Front, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/NT_v_HH11_Invasive_Front.txt", sep = "\t", row.names = FALSE)
#write.table(NT_v_HH13_Invasive_Front, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/NT_v_HH13_Invasive_Front.txt", sep = "\t", row.names = FALSE)
#write.table(NT_v_HH15_Invasive_Front, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/NT_v_HH15_Invasive_Front.txt", sep = "\t", row.names = FALSE)
#write.table(NT_v_HH13_Leader, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/NT_v_NT_v_HH13_Leader.txt", sep = "\t", row.names = FALSE)
#write.table(NT_v_HH13_Trailer, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/NT_v_HH13_Trailer.txt", sep = "\t", row.names = FALSE)
#write.table(NT_v_HH15_Leader, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/NT_v_HH15_Leader.txt", sep = "\t", row.names = FALSE)
#write.table(NT_v_HH15_Trailer, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/NT_v_HH15_Trailer.txt", sep = "\t", row.names = FALSE)


```

## Combine in vitro and in vivo DE gene lists

```{r}

sc_DE_genelist <- as.vector(c(sc_invitro_DE_genelist, sc_invivo_DE_genelist)) #10925
sc_DE_genelist <- unique(sc_DE_genelist) #3723

```

## Unclustered single cell heatmap including in vitro cells

```{r}

# subset edger.rpkm.sc to pull out unique genes
unique_genes.iv <- row.names(edger.rpkm.sc) %in% sc_DE_genelist

# add time and celltype categories to targets for annotation
time<-gsub("_.*","",targets$group)
celltype<-gsub("HH[1-9]+_","",targets$group)
# add to targets df
targets$time <- time
targets$celltype <- celltype

# Fix NT culture 

targets$time[with(targets, group == "NT_culture")] <- "NT_culture"
targets$celltype[with(targets, group == "NT_culture")] <- "NT_culture"

# manipulate edger.sub for heatmap (log2, z-score, remove inf/NA)
temp <- edger.rpkm.sc[unique_genes.iv,]
temp<-log2(temp+.1)
temp<-scale(t(temp), scale = T, center = T)
temp <- data.frame(t(temp))
temp[temp>3]<- 3
temp[temp< -3]<- -3
df <- temp[complete.cases(temp),]

## Unclustered
unique_genes.iv

#pdf("/n/projects/lwo/singlecell/R_out/Finalized_figures/scRNAseq_heatmap_allcells_unclustered.pdf", height = 15, width = 12)
aheatmap(df, breaks = 0, color = rev(brewer.pal(11,"RdBu")), annColors = ann_color, scale="none",annCol=targets[,c("group","time","celltype")], distfun = "pearson", Colv = NA)
#dev.off()

```

## Z-score filter single cell heat map

```{r}

zscore_count <- data.frame(number_greater_than_0 = rowSums(df > 0))

## Histogram of number of cells expression

hist(zscore_count$number_greater_than_0)

## Set z-score threshold at 30
## 30 cells must have genes expressed above the mean for that row for it to be included

zscore_count <- subset(zscore_count, number_greater_than_0 > 60) #760

hist(zscore_count$number_greater_than_0)

zscore.iv <- row.names(df) %in% row.names(zscore_count)

#pdf("/n/projects/lwo/singlecell/R_out/Finalized_figures/scRNAseq_heatmap_allcells_unclustered_z_score_60.pdf", height = 10, width = 15)
aheatmap(df[zscore.iv,], breaks = 0, color = rev(brewer.pal(11,"RdBu")), annColors = ann_color, scale="none",annCol=targets[,c("group","time","celltype")], distfun = "pearson", Colv = NA)
#dev.off()
```


## Correlation matrix

Use top and bottom 2% of genes from PC1 and 2, overlap with unique DE list

```{r}
cor_genes <- intersect(collected2pct, sc_DE_genelist)

temp <- edger.rpkm.sc[row.names(edger.rpkm.sc) %in% cor_genes,]
temp <- log2(temp + 0.1)
cor_mat <- cor(temp)

#pdf("/n/projects/lwo/singlecell/R_out/Figures/Figures_7_12_2016/pearson_corr_1pct.pdf", height= 12, width = 12)
aheatmap(cor_mat, color = "-RdBu:100", cellwidth = 1, cellheight = 1, breaks = 0.5, Colv = NA, Rowv= NA,  annRow=targets[,c("group","time","celltype")], annCol=targets[,c("group","time","celltype")], annColors = ann_color)
#dev.off()

```


## Unique genes (in vivo only) sc heatmap
```{r}

# subset edger.rpkm.sc to pull out unique genes
unique_genes_vivo.iv <- row.names(edger.rpkm.sc) %in% sc_invivo_DE_genelist

# Do not include NT in heat map
rest.iv <- !grepl("NT", names(edger.rpkm.sc))

# manipulate edger.sub for heatmap (log2, z-score, remove inf/NA)
temp <- edger.rpkm.sc[unique_genes_vivo.iv,rest.iv]
temp<-log2(temp+.1)
temp<-scale(t(temp), scale = T, center = T)
temp <- data.frame(t(temp))
temp[temp>3]<-3
temp[temp< -3]<- -3
df <- temp[complete.cases(temp),]

## Clustered

# Pearson

#pdf("/n/projects/lwo/singlecell/R_out/Finalized_figures/scRNAseq_clustered_heatmap_invivo_pearson.pdf", height = 15, width = 12)
#aheatmap(df, breaks = 0, color = rev(brewer.pal(11,"RdBu")), annColors = ann_color, scale="none",annCol=targets[rest.iv,c("group","time","celltype")], distfun = "pearson")
#dev.off()

# Euclidean 

#pdf("/n/projects/lwo/singlecell/R_out/Finalized_figures/scRNAseq_clustered_heatmap_invivo_euclid.pdf", height = 15, width = 12)
#aheatmap(df, breaks = 0, color = rev(brewer.pal(11,"RdBu")), annColors = ann_color, scale="none",annCol=targets[rest.iv,c("group","time","celltype")])
#dev.off()

## Unclustered

#pdf("/n/projects/lwo/singlecell/R_out/Finalized_figures/scRNAseq_unclustered_heatmap_invivo.pdf", height = 15, width = 12)
aheatmap(df, breaks = 0, color = rev(brewer.pal(11,"RdBu")), annColors = ann_color, scale="none",annCol=targets[rest.iv,c("group","time","celltype")], Colv = NA,distfun = "pearson")
#dev.off()
```

## Unique genes sc heatmap AVERAGED
```{r}

# use all DE genes from in vivo and in vitro DE
scDE.iv <- row.names(edger.rpkm.sc) %in% sc_DE_genelist
temp <- edger.rpkm.sc[scDE.iv,]


# Create logical vectors to subset data and get averages
IF11.iv <- grepl("HH11_Invasive_Front", names(edger.rpkm.sc))
IF13.iv <- grepl("HH13_Invasive_Front", names(edger.rpkm.sc))
IF15.iv <- grepl("HH15_Invasive_Front", names(edger.rpkm.sc))
L13.iv <- grepl("HH13_Leader", names(edger.rpkm.sc))
L15.iv <- grepl("HH15_Leader", names(edger.rpkm.sc))
Tr13.iv <- grepl("HH13_Trailer", names(edger.rpkm.sc))
Tr15.iv <- grepl("HH15_Trailer", names(edger.rpkm.sc))
NT.iv <- grepl("NT", names(edger.rpkm.sc))

# Take averages for each subgroup
IF11_mean <- rowMeans(temp[,IF11.iv])
IF13_mean <- rowMeans(temp[,IF13.iv])
IF15_mean <- rowMeans(temp[,IF15.iv])
L13_mean <- rowMeans(temp[,L13.iv])
L15_mean <- rowMeans(temp[,L15.iv])
Tr13_mean <- rowMeans(temp[,Tr13.iv])
Tr15_mean <- rowMeans(temp[,Tr15.iv])
NT_mean <- rowMeans(temp[,NT.iv])


# Collect means in Data Frame
sc_means <- data.frame(IF11_mean, IF13_mean, L13_mean, Tr13_mean, IF15_mean, L15_mean, Tr15_mean, NT_mean)

temp<-log2(sc_means+.1)
temp<-scale(t(temp), center = T, scale = T)
temp <- data.frame(t(temp))
temp[temp>3]<-3
temp[temp< -3]<- -3
df <- temp[complete.cases(temp),]

#pdf("/n/projects/lwo/singlecell/R_out/Finalized_figures/scRNAseq_average_heatmap_pearson.pdf", height = 15, width = 12)
aheatmap(df, breaks = 0, color = rev(brewer.pal(11,"RdBu")), scale="none", distfun = "pearson")
#dev.off()

#pdf("/n/projects/lwo/singlecell/R_out/Finalized_figures/scRNAseq_average_heatmap_euclid.pdf", height = 15, width = 12)
aheatmap(df, breaks = 0, color = rev(brewer.pal(11,"RdBu")), scale="none")
#dev.off()
```

## DE vivo vs vitro

```{r}
# Create in vivo vs in vitor group
targets$vivo_vitro <- ifelse(targets$group == "NT_culture", "in_vitro", "in_vivo")


# Note that the first group listed in the pair is the baseline for the comparison—so if the pair is c("A","B") then the comparison is B - A, so genes with positive log-fold change are up-regulated in group B compared with group A (and vice versa for genes with negative log-fold change).


ucount_sub <- ucount_df[2:length(ucount_df)]
y <- DGEList(counts=ucount_sub, group=targets$vivo_vitro)
keep<-apply(cpm(y),1,function(x) length(x[x>5])>=2)
y <- y[keep, , keep.lib.sizes=FALSE]
y <- calcNormFactors(y)
y <- estimateCommonDisp(y)
y <- estimateTagwiseDisp(y)

res<-list()
res[[1]] <- exactTest(y, pair=c("in_vitro", "in_vivo"))

conts <- c("in_vivo_vs_in_vitro") 

# Add Padj
for(i in 1:length(conts))
{
  res[[i]]$table$padj<-p.adjust(res[[i]]$table$PValue,method="BH")
}

# Summary table
summarytable <- data.frame(ucount_df[which(keep),1],res[[1]]$table)
# Column lables
colnames(summarytable)[1]<-"ensid"
colnames(summarytable)[2:length(summarytable[1,])]<-paste(c(paste(rep(conts,each=4),rep(c(".log2fc",".log2cpm",".pval", ".padj"),times=1),sep='')))

# add annotation to summary of DE analysis
summarytable_ann_invitro_v_invivo<-data.frame(ensid=summarytable[,1],genedata[match(summarytable[,1],genedata[,1]),c("Name", "Chr", "Start", "End", "Strand", "Biotype","Description")],summarytable[,2:length(summarytable[1,])])
rownames(summarytable_ann_invitro_v_invivo)<-summarytable_ann_invitro_v_invivo$ensid

#### Get sig genes ####

vivo_v_vitro.up <- subset(summarytable_ann_invitro_v_invivo, summarytable_ann_invitro_v_invivo$in_vivo_vs_in_vitro.log2fc > 2 & summarytable_ann_invitro_v_invivo$in_vivo_vs_in_vitro.padj < 0.05)

vivo_v_vitro.dn <- subset(summarytable_ann_invitro_v_invivo, summarytable_ann_invitro_v_invivo$in_vivo_vs_in_vitro.log2fc < -2 & summarytable_ann_invitro_v_invivo$in_vivo_vs_in_vitro.padj < 0.05)

#write.table(vivo_v_vitro.up, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/vivo_v_vitro_up.txt", row.names = FALSE, sep = "\t")
#write.table(vivo_v_vitro.dn, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/vivo_v_vitro_dn.txt", row.names = FALSE, sep = "\t")

```



## PCA of in vitro/in vivo

```{r}
vivo_v_vitro_genelist <- union(vivo_v_vitro.up$ensid, vivo_v_vitro.dn$ensid)

# PCA with all genes
temp <- t(edger.rpkm.sc[row.names(edger.rpkm.sc) %in% sc_DE_genelist,])
temp <- log2(temp +0.01)
rpkm.pca<-prcomp(temp,center=T,scale=T)

#pdf("/n/projects/lwo/singlecell/R_out/Finalized_figures/biplots_DEgenes_vivo_vitro.pdf", height = 15, width= 15)
print(ggbiplot(rpkm.pca,groups=targets$vivo_vitro,ellipse=F,circle=F,var.axes=F,alpha=.7, size = 3) + theme_minimal())
print(ggbiplot(rpkm.pca,choices=2:3,groups=targets$vivo_vitro,ellipse=F,circle=F,var.axes=F,alpha=.7)+ theme_minimal())
print(ggbiplot(rpkm.pca,choices=3:4,groups=targets$vivo_vitro,ellipse=F,circle=F,var.axes=F,alpha=.7)+ theme_minimal())
#dev.off()

nobs.factor <- sqrt(nrow(rpkm.pca$x) - 1)
d <- rpkm.pca$sdev
v <- rpkm.pca$rotation
u <- sweep(rpkm.pca$x, 2, 1 / (d * nobs.factor), FUN = '*')

df.u <- as.data.frame(sweep(u[,c(1:4)], 2, d[1:4]^0, FUN='*'))

df.u <- df.u * nobs.factor

df <- data.frame(PC1 = df.u$PC1, PC2 = df.u$PC2, PC3 = df.u$PC3, row.names = row.names(rpkm.pca$x))



pdf("/n/projects/lwo/singlecell/Finalized_folder/pdf_out/vivo_vitro_PCA_12.pdf", height = 10, width = 13)
ggplot(df, aes(x = PC1, y = PC2)) + 
  geom_point(alpha = 0.8, size = 5, aes(color = targets$vivo_vitro))+
  scale_color_manual(values = c("slateblue4", "goldenrod3"))+
  ylab("standardized PC2 (2.5% explained var.)") +
  xlab("standardized PC1 (3.0% explained var.)") +
  theme_classic()
dev.off()
```

## Vivo vs vitro Pearson plot averaged DE genes

```{r}
temp <- edger.rpkm.sc[row.names(edger.rpkm.sc) %in% sc_DE_genelist,]

IF11_mean <- rowMeans(temp[,IF11.iv])
IF13_mean <- rowMeans(temp[,IF13.iv])
IF15_mean <- rowMeans(temp[,IF15.iv])
L13_mean <- rowMeans(temp[,L13.iv])
L15_mean <- rowMeans(temp[,L15.iv])
Tr13_mean <- rowMeans(temp[,Tr13.iv])
Tr15_mean <- rowMeans(temp[,Tr15.iv])
NT_mean <- rowMeans(temp[,NT.iv])


# Collect means in Data Frame
sc_means <- data.frame(IF11_mean, IF13_mean, L13_mean, Tr13_mean, IF15_mean, L15_mean, Tr15_mean, NT_mean)

temp <- sc_means
temp <- log2(temp + 0.1)
cor_mat <- cor(temp)

pdf("/n/projects/lwo/singlecell/Finalized_folder/pdf_out/averaged_pearson_DEgenes_allsubpops.pdf", height = 10, width = 12)
aheatmap(cor_mat, color = "-RdBu:100", cellwidth = 10, cellheight = 10, breaks = 0.7, distfun = "pearson")
dev.off()

```



## DE gene list of HH13 and HH15 IF vs HH13 and HH15 Leaders/Trailers

```{r}
# create Targets group that includes stream as a subgroup

targets$group2 <- gsub("Leader", "Stream", targets$group)
targets$group2 <- gsub("Trailer", "Stream", targets$group2)
targets$group3 <- gsub("HH[0-9]+_", "", targets$group2)

ucount_sub <- ucount_df[,-1]
y <- DGEList(counts=ucount_sub, group=targets$group3)
keep<-apply(cpm(y),1,function(x) length(x[x>5])>=2)
y <- y[keep, , keep.lib.sizes=FALSE]
y <- calcNormFactors(y)
y <- estimateCommonDisp(y)
y <- estimateTagwiseDisp(y)

res <- list()
res[[1]] <- exactTest(y, pair = c("Stream", "Invasive_Front"))


conts <- c("Invasive_Front_v_Stream")

# Add Padj
for(i in 1:length(conts))
{
  res[[i]]$table$padj<-p.adjust(res[[i]]$table$PValue,method="BH")
}

# Summary table
summarytable <- data.frame(ucount_df[which(keep),1],res[[1]]$table)

# Column lables
colnames(summarytable)[1]<-"ensid"
colnames(summarytable)[2:length(summarytable[1,])]<-paste(c(paste(rep(conts,each=4),rep(c(".log2fc",".log2cpm",".pval", ".padj"),times=1),sep='')))

# add annotation to summary of DE analysis
summarytable_ann<-data.frame(ensid=summarytable[,1],genedata[match(summarytable[,1],genedata[,1]),c("Name", "Chr", "Start", "End", "Strand", "Biotype","Description")],summarytable[,2:length(summarytable[1,])])
rownames(summarytable_ann)<-summarytable_ann$ensid

Invasive_Front_v_Stream_singlecell.up <- subset(summarytable_ann, Invasive_Front_v_Stream.log2fc > 2 & Invasive_Front_v_Stream.padj < 0.05)
Invasive_Front_v_Stream_singlecell.dn <- subset(summarytable_ann, Invasive_Front_v_Stream.log2fc < 2 & Invasive_Front_v_Stream.padj < 0.05)

#write.table(Invasive_Front_v_Stream_singlecell.dn, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/Invasive_Front_v_Stream_combined13adn15_dn.txt", row.names = F, sep = "\t")

#write.table(Invasive_Front_v_Stream_singlecell.up, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/Invasive_Front_v_Stream_combined13and15_up.txt", row.names = F, sep = "\t")

```

## Combine Invasive Front gene lists
```{r}
# All IF vs each other
IF_gene_list <- c(as.character(HH11_Invasive_Front_v_HH15_Invasive_Front_DE.up$ensid), as.character(HH11_Invasive_Front_v_HH15_Invasive_Front_DE.dn$ensid), as.character(HH11_Invasive_Front_v_HH13_Invasive_Front_DE.dn$ensid), as.character(HH11_Invasive_Front_v_HH13_Invasive_Front_DE.up$ensid), as.character(HH13_Invasive_Front_v_HH15_Invasive_Front_DE.dn$ensid), as.character(HH13_Invasive_Front_v_HH15_Invasive_Front_DE.dn$ensid))

IF_gene_list <- unique(IF_gene_list) #1433

IF_gene_list <- c(IF_gene_list, as.character(Invasive_Front_v_Stream_singlecell.up$ensid), as.character(Invasive_Front_v_Stream_singlecell.dn$ensid))

IF_gene_list <- unique(IF_gene_list) # 2206

```

## Invasive_Front heat map
```{r}
IF.iv <- grepl("Invasive_Front", names(edger.rpkm.sc))

temp <- edger.rpkm.sc[row.names(edger.rpkm.sc) %in% IF_gene_list, IF.iv] 

row.names(temp) <- paste(row.names(temp), '  ', genedata[match(row.names(temp), genedata[,1]),2], sep = '')

temp<-log2(temp+.1)
temp<-scale(t(temp), scale =T, center = T)
temp<- data.frame(t(temp))
temp[temp>3]<-3
temp[temp< -3]<- -3
df<- as.data.frame(temp[complete.cases(temp),])


#pdf("/n/projects/lwo/singlecell/R_out/Finalized_figures/IF_scRNAseq_heatmap_euclid.pdf", height = 30, width = 20)
aheatmap(df, breaks = 0, color = rev(brewer.pal(11,"RdBu")), scale="none", labCol = targets$alias[IF.iv], annColors = ann_color[group.iv], annCol=targets[IF.iv,c("group")])
#dev.off()

#pdf("/n/projects/lwo/singlecell/R_out/Finalized_figures/IF_scRNAseq_heatmap_pearson.pdf", height = 30, width = 20)
aheatmap(df, distfun = "pearson", breaks = 0, color = rev(brewer.pal(11,"RdBu")), scale="none", labCol = targets$alias[IF.iv], annColors = ann_color[group.iv], annCol=targets[IF.iv,c("group")])
#dev.off()
```

## Remove genes that have low variation

For example, one or two cells express very high

```{r}
temp <- edger.rpkm.sc[row.names(edger.rpkm.sc) %in% IF_gene_list, IF.iv] 
temp<-log2(temp+.1)
temp<-scale(t(temp), scale =T, center = T)
temp<- data.frame(t(temp))
temp[temp>3]<-3
temp[temp< -3]<- -3
df<- as.data.frame(temp[complete.cases(temp),])

IF_zscore_count <- data.frame(number_greater_than_0 = rowSums(df > 0))

## Histogram of number of cells expression

hist(IF_zscore_count$number_greater_than_0)
```

Due to technical/biological variablity in scRNA-seq there are many genes with only a few cells highly expressing. We are interested in broader trends than a few cells. 

## Z-score filter of IF DE genes

```{r}

## Set z-score threshold at 30
## 30 cells must have genes expressed above the mean for that row for it to be included

IF_zscore_count <- subset(IF_zscore_count, number_greater_than_0 > 30) #760

hist(IF_zscore_count$number_greater_than_0)

zscorefilt.iv <- row.names(df) %in% row.names(IF_zscore_count)

df <- df[zscorefilt.iv,]

#pdf("/n/projects/lwo/singlecell/R_out/Finalized_figures/IF_scRNAseq_heatmap_zfiltered_760_thresh30_pearson.pdf", height = 30, width = 20)
#png("/n/projects/lwo/singlecell/R_out/YISR/IF_scRNAseq_heatmap_zfiltered_760_thresh30.png", height = 20, width = 30, res = 600, units = "in")
aheatmap(df, breaks = 0, color = rev(brewer.pal(11,"RdBu")), scale="none", labCol = targets$alias[IF.iv], annColors = ann_color[group.iv], annCol=targets[IF.iv,c("group")], distfun = "pearson")
#dev.off()

#pdf("/n/projects/lwo/singlecell/R_out/Finalized_figures/IF_scRNAseq_heatmap_zfiltered_760_thresh30_euclid.pdf", height = 30, width = 20)
#png("/n/projects/lwo/singlecell/R_out/YISR/IF_scRNAseq_heatmap_zfiltered_760_thresh30.png", height = 20, width = 30, res = 600, units = "in")
aheatmap(df, breaks = 0, color = rev(brewer.pal(11,"RdBu")), scale="none", labCol = targets$alias[IF.iv], annColors = ann_color[group.iv], annCol=targets[IF.iv,c("group")], distfun = "euclidean")
#dev.off()

```

## Get dendrogram from heat map
```{r}

heatmap.obj <- aheatmap(df, breaks = 0, color = rev(brewer.pal(11,"RdBu")), scale="none", labCol = targets$alias[IF.iv], annColors = ann_color[group.iv], annCol=targets[IF.iv,c("group")], distfun = "pearson")

dend <- heatmap.obj$Colv

#pdf("/n/projects/lwo/singlecell/R_out/dend_pearson.pdf", width = 15, height = 10)
par(mar=c(10,2,2,2), cex = 0.5)
plot(dend)
#dev.off()
```

*****






## CLUSTERS


## DE clusters vs all other in vivo

```{r}

# Take ordered_dend_list
# Apply to aheatmap ordered list and match clusters to cells "sections of interest" on heatmap
# Done in Excel
# read in excel table

annotated_cluster <- read.table("/n/projects/lwo/singlecell/Finalized_txt/annotated_IFcluster_pearson.txt", sep = "\t", header=  TRUE)

noNT.iv <- !grepl("NT|ensid", names(ucount_df))


# Create groups
group <- data.frame(cell = names(ucount_df[noNT.iv]))
group$cell <- as.character(group$cell)
group$all <- as.character(annotated_cluster[match(group[,1], annotated_cluster$Cell),"Cluster"])
group[is.na(group)] <- "other"
group$Alpha_v_other <- ifelse(group$all == "Alpha", "Alpha", "other")
group$Beta_v_other <- ifelse(group$all == "Beta", "Beta", "other")
group$Gamma_v_other <- ifelse(group$all == "Gamma", "Gamma", "other")
group$Delta_v_other <- ifelse(group$all == "Delta", "Delta", "other")
group$Epsilon_v_other <- ifelse(group$all == "Epsilon", "Epsilon", "other")

y <- DGEList(counts=ucount_df[,noNT.iv], group=group$Beta_v_other)
keep<-apply(cpm(y),1,function(x) length(x[x>5])>=2)
y <- y[keep, , keep.lib.sizes=FALSE]
y <- calcNormFactors(y)
y <- estimateCommonDisp(y)
y <- estimateTagwiseDisp(y)

# Note that the first group listed in the pair is the baseline for the comparison—so if the pair is c("A","B") then the comparison is B - A, so genes with positive log-fold change are up-regulated in group B compared with group A (and vice versa for genes with negative log-fold change).

res<-list()
res[[1]] <- exactTest(y, pair=c("other", "Beta"))

y <- DGEList(counts=ucount_df[,noNT.iv], group=group$Alpha_v_other)
keep<-apply(cpm(y),1,function(x) length(x[x>5])>=2)
y <- y[keep, , keep.lib.sizes=FALSE]
y <- calcNormFactors(y)
y <- estimateCommonDisp(y)
y <- estimateTagwiseDisp(y)


res[[2]] <- exactTest(y, pair=c("other", "Alpha"))

y <- DGEList(counts=ucount_df[,noNT.iv], group=group$Gamma_v_other)
keep<-apply(cpm(y),1,function(x) length(x[x>5])>=2)
y <- y[keep, , keep.lib.sizes=FALSE]
y <- calcNormFactors(y)
y <- estimateCommonDisp(y)
y <- estimateTagwiseDisp(y)


res[[3]] <- exactTest(y, pair=c("other", "Gamma"))

y <- DGEList(counts=ucount_df[,noNT.iv], group=group$Delta_v_other)
keep<-apply(cpm(y),1,function(x) length(x[x>5])>=2)
y <- y[keep, , keep.lib.sizes=FALSE]
y <- calcNormFactors(y)
y <- estimateCommonDisp(y)
y <- estimateTagwiseDisp(y)

res[[4]] <- exactTest(y, pair=c("other", "Delta"))

y <- DGEList(counts=ucount_df[,noNT.iv], group=group$Epsilon_v_other)
keep<-apply(cpm(y),1,function(x) length(x[x>5])>=2)
y <- y[keep, , keep.lib.sizes=FALSE]
y <- calcNormFactors(y)
y <- estimateCommonDisp(y)
y <- estimateTagwiseDisp(y)


res[[5]] <- exactTest(y, pair=c("other", "Epsilon"))


# List of conts
# X_v_Y == X relative to Y
conts <- c("Beta_v_other","Alpha_v_other","Gamma_v_other","Delta_v_other","Epsilon_v_other") 

# Add Padj
for(i in 1:length(conts))
{
  res[[i]]$table$padj<-p.adjust(res[[i]]$table$PValue,method="BH")
}

# Summary table
summarytable <- data.frame(ucount_df[which(keep),1],res[[1]]$table,res[[2]]$table,res[[3]]$table,res[[4]]$table, res[[5]]$table)
# Column lables
colnames(summarytable)[1]<-"ensid"
colnames(summarytable)[2:length(summarytable[1,])]<-paste(c(paste(rep(conts,each=4),rep(c(".log2fc",".log2cpm",".pval", ".padj"),times=1),sep='')))

# add annotation to summary of DE analysis
summarytable_ann_clusters_v_allinvivo<-data.frame(ensid=summarytable[,1],genedata[match(summarytable[,1],genedata[,1]),c("Name", "Chr", "Start", "End", "Strand", "Biotype","Description")],summarytable[,2:length(summarytable[1,])])
rownames(summarytable_ann_clusters_v_allinvivo)<-summarytable_ann_clusters_v_allinvivo$ensid

#write.table(summarytable_ann_clusters_v_allinvivo, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/summary_tables/summary_clusters_v_invivo_DE.txt", sep = "\t", row.names = FALSE)
```

## Characterize DEG from cluster vs all other in vivo
For each DEG include %expression in subgroup and avg level of expression
```{r}

Alpha_cells <- as.character(annotated_cluster[grepl("Alpha", annotated_cluster$Cluster),1])
Beta_cells <- as.character(annotated_cluster[grepl("Beta", annotated_cluster$Cluster),1])
Gamma_cells <- as.character(annotated_cluster[grepl("Gamma", annotated_cluster$Cluster),1])
Delta_cells <- as.character(annotated_cluster[grepl("Delta", annotated_cluster$Cluster),1])
Epsilon_cells <- as.character(annotated_cluster[grepl("Epsilon", annotated_cluster$Cluster),1])

# Get up DE genes for each cluster
Alpha_v_other_DE.up <- subset(summarytable_ann_clusters_v_allinvivo, summarytable_ann_clusters_v_allinvivo$Alpha_v_other.log2fc > 2 & summarytable_ann_clusters_v_allinvivo$Alpha_v_other.padj < 0.05)
Beta_v_other_DE.up <- subset(summarytable_ann_clusters_v_allinvivo, summarytable_ann_clusters_v_allinvivo$Beta_v_other.log2fc > 2 & summarytable_ann_clusters_v_allinvivo$Beta_v_other.padj < 0.05)
Gamma_v_other_DE.up <- subset(summarytable_ann_clusters_v_allinvivo, summarytable_ann_clusters_v_allinvivo$Gamma_v_other.log2fc > 2 & summarytable_ann_clusters_v_allinvivo$Gamma_v_other.padj < 0.05)
Delta_v_other_DE.up <- subset(summarytable_ann_clusters_v_allinvivo, summarytable_ann_clusters_v_allinvivo$Delta_v_other.log2fc > 2 & summarytable_ann_clusters_v_allinvivo$Delta_v_other.padj < 0.05)
Epsilon_v_other_DE.up <- subset(summarytable_ann_clusters_v_allinvivo, summarytable_ann_clusters_v_allinvivo$Epsilon_v_other.log2fc > 2 & summarytable_ann_clusters_v_allinvivo$Epsilon_v_other.padj < 0.05)

# Get down DE
Alpha_v_other_DE.dn <- subset(summarytable_ann_clusters_v_allinvivo, summarytable_ann_clusters_v_allinvivo$Alpha_v_other.log2fc < -2 & summarytable_ann_clusters_v_allinvivo$Alpha_v_other.padj < 0.05)
Beta_v_other_DE.dn <- subset(summarytable_ann_clusters_v_allinvivo, summarytable_ann_clusters_v_allinvivo$Beta_v_other.log2fc < -2 & summarytable_ann_clusters_v_allinvivo$Beta_v_other.padj < 0.05)
Gamma_v_other_DE.dn <- subset(summarytable_ann_clusters_v_allinvivo, summarytable_ann_clusters_v_allinvivo$Gamma_v_other.log2fc < -2 & summarytable_ann_clusters_v_allinvivo$Gamma_v_other.padj < 0.05)
Delta_v_other_DE.dn <- subset(summarytable_ann_clusters_v_allinvivo, summarytable_ann_clusters_v_allinvivo$Delta_v_other.log2fc < -2 & summarytable_ann_clusters_v_allinvivo$Delta_v_other.padj < 0.05)
Epsilon_v_other_DE.dn <- subset(summarytable_ann_clusters_v_allinvivo, summarytable_ann_clusters_v_allinvivo$Epsilon_v_other.log2fc < -2 & summarytable_ann_clusters_v_allinvivo$Epsilon_v_other.padj < 0.05)

# logical vectors to pull out cells by cluster
# remove .rpkm from edger.rpkm.sc names so we can match to DEG_dfs
names(edger.rpkm) <- gsub("rpkm.", "", names(edger.rpkm.sc))
Alpha.iv <- names(edger.rpkm) %in% Alpha_cells
Beta.iv <- names(edger.rpkm) %in% Beta_cells
Gamma.iv <- names(edger.rpkm) %in% Gamma_cells
Delta.iv <- names(edger.rpkm) %in% Delta_cells
Epsilon.iv <- names(edger.rpkm) %in% Epsilon_cells


#### Alpha #####


# Put DE genes and symbols in df
Alpha_DEG_df <- data.frame(Alpha_v_other_DE.up[,c(1,2)], Alpha_v_other_DE.up$Alpha_v_other.log2fc, Alpha_v_other_DE.up$Alpha_v_other.padj)
#A_DEG_df <- A_DEG_df[A_DEG_df$ensid %in% genelist,]


# Calculate percentage expressing AND mean expressio

Alpha_DEG_df$mean_rpkm <- rowMeans(edger.rpkm.sc[row.names(edger.rpkm.sc) %in% Alpha_DEG_df$ensid, Alpha.iv])

# pct above 0 rpkm
Alpha_DEG_df$pct_above_0_rpkm <- rowSums(edger.rpkm.sc[row.names(edger.rpkm.sc) %in% Alpha_DEG_df$ensid, Alpha.iv] > 0)/length(Alpha_cells)

##### Beta #####


# Put DE genes and symbols in df
Beta_DEG_df <- data.frame(Beta_v_other_DE.up[,c(1,2)], Beta_v_other_DE.up$Beta_v_other.log2fc, Beta_v_other_DE.up$Beta_v_other.padj)

# Calculate percentage expressing AND mean expression
Beta_DEG_df$mean_rpkm <- rowMeans(edger.rpkm.sc[row.names(edger.rpkm.sc) %in% Beta_DEG_df$ensid, Beta.iv])

# pct above 0 rpkm
Beta_DEG_df$pct_above_0_rpkm <- rowSums(edger.rpkm.sc[row.names(edger.rpkm.sc) %in% Beta_DEG_df$ensid, Beta.iv] > 0)/length(Beta_cells)


##### Gamma #####


# Put DE genes and symCols in df
Gamma_DEG_df <- data.frame(Gamma_v_other_DE.up[,c(1,2)], Gamma_v_other_DE.up$Gamma_v_other.log2fc, Gamma_v_other_DE.up$Gamma_v_other.padj)

# Calculate percentage expressing AND mean expression
Gamma_DEG_df$mean_rpkm <- rowMeans(edger.rpkm.sc[row.names(edger.rpkm.sc) %in% Gamma_DEG_df$ensid, Gamma.iv])

# pct above 0 rpkm
Gamma_DEG_df$pct_above_0_rpkm <- rowSums(edger.rpkm.sc[row.names(edger.rpkm.sc) %in% Gamma_DEG_df$ensid, Gamma.iv] > 0)/length(Gamma_cells)


##### D #####


Delta_DEG_df <- data.frame(Delta_v_other_DE.up[,c(1,2)], Delta_v_other_DE.up$Delta_v_other.log2fc, Delta_v_other_DE.up$Delta_v_other.padj)

# Calculate percentage expressing AND mean expression
Delta_DEG_df$mean_rpkm <- rowMeans(edger.rpkm.sc[row.names(edger.rpkm.sc) %in% Delta_DEG_df$ensid, Delta.iv])

# pct above 0 rpkm
Delta_DEG_df$pct_above_0_rpkm <- rowSums(edger.rpkm.sc[row.names(edger.rpkm.sc) %in% Delta_DEG_df$ensid, Delta.iv] > 0)/length(Delta_cells)


##### E #####


Epsilon_DEG_df <- data.frame(Epsilon_v_other_DE.up[,c(1,2)], Epsilon_v_other_DE.up$Epsilon_v_other.log2fc, Epsilon_v_other_DE.up$Epsilon_v_other.padj)

# Calculate percentage expressing AND mean expression
Epsilon_DEG_df$mean_rpkm <- rowMeans(edger.rpkm.sc[row.names(edger.rpkm.sc) %in% Epsilon_DEG_df$ensid, Epsilon.iv])


# pct above 0 rpkm
Epsilon_DEG_df$pct_above_0_rpkm <- rowMeans(edger.rpkm.sc[row.names(edger.rpkm.sc) %in% Epsilon_DEG_df$ensid, Epsilon.iv] > 0/length(Epsilon_cells))

#write.table(Alpha_DEG_df, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/DEG_info_Alpha_up.txt", row.names = FALSE, sep = "\t")
#write.table(Beta_DEG_df, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/DEG_info_Beta_up.txt", row.names = FALSE, sep = "\t")
#write.table(Gamma_DEG_df, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/DEG_info_Gamma_up.txt", row.names = FALSE, sep = "\t")
#write.table(Delta_DEG_df, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/DEG_info_Delta_up.txt", row.names = FALSE, sep = "\t")
#write.table(Epsilon_DEG_df, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/DEG_info_Epsilon_up.txt", row.names = FALSE, sep = "\t")

###################
# DOWN GENES
#### A #####


# Put DE genes and symbols in df
Alpha_DEG_df <- data.frame(Alpha_v_other_DE.dn[,c(1,2)], Alpha_v_other_DE.dn$Alpha_v_other.log2fc, Alpha_v_other_DE.dn$Alpha_v_other.padj)

# Calculate percentage expressing AND mean expressio

Alpha_DEG_df$mean_rpkm <- rowMeans(edger.rpkm.sc[row.names(edger.rpkm.sc) %in% Alpha_DEG_df$ensid, Alpha.iv])

# pct above 0 rpkm
Alpha_DEG_df$pct_above_0_rpkm <- rowSums(edger.rpkm.sc[row.names(edger.rpkm.sc) %in% Alpha_DEG_df$ensid, Alpha.iv] > 0)/length(Alpha_cells)


##### B #####


# Put DE genes and symbols in df
Beta_DEG_df <- data.frame(Beta_v_other_DE.dn[,c(1,2)], Beta_v_other_DE.dn$Beta_v_other.log2fc, Beta_v_other_DE.dn$Beta_v_other.padj)

# Calculate percentage expressing AND mean expression
Beta_DEG_df$mean_rpkm <- rowMeans(edger.rpkm.sc[row.names(edger.rpkm.sc) %in% Beta_DEG_df$ensid, Beta.iv])

# pct above 0 rpkm
Beta_DEG_df$pct_above_0_rpkm <- rowSums(edger.rpkm.sc[row.names(edger.rpkm.sc) %in% Beta_DEG_df$ensid, Beta.iv] > 0)/length(Beta_cells)


##### C #####


# Put DE genes and symCols in df
Gamma_DEG_df <- data.frame(Gamma_v_other_DE.dn[,c(1,2)], Gamma_v_other_DE.dn$Gamma_v_other.log2fc, Gamma_v_other_DE.dn$Gamma_v_other.padj)

# Calculate percentage expressing AND mean expression
Gamma_DEG_df$mean_rpkm <- rowMeans(edger.rpkm.sc[row.names(edger.rpkm.sc) %in% Gamma_DEG_df$ensid, Gamma.iv])

# pct above 0 rpkm
Gamma_DEG_df$pct_above_0_rpkm <- rowSums(edger.rpkm.sc[row.names(edger.rpkm.sc) %in% Gamma_DEG_df$ensid, Gamma.iv] > 0)/length(Gamma_cells)

##### D #####

Delta_DEG_df <- data.frame(Delta_v_other_DE.dn[,c(1,2)], Delta_v_other_DE.dn$Delta_v_other.log2fc, Delta_v_other_DE.dn$Delta_v_other.padj)

# Calculate percentage expressing AND mean expression
Delta_DEG_df$mean_rpkm <- rowMeans(edger.rpkm.sc[row.names(edger.rpkm.sc) %in% Delta_DEG_df$ensid, Delta.iv])

# pct above 0 rpkm
Delta_DEG_df$pct_above_0_rpkm <- rowSums(edger.rpkm.sc[row.names(edger.rpkm.sc) %in% Delta_DEG_df$ensid, Delta.iv] > 0)/length(Delta_cells)


##### E #####

Epsilon_DEG_df <- data.frame(Epsilon_v_other_DE.dn[,c(1,2)], Epsilon_v_other_DE.dn$Epsilon_v_other.log2fc, Epsilon_v_other_DE.dn$Epsilon_v_other.padj)

# Calculate percentage expressing AND mean expression
Epsilon_DEG_df$mean_rpkm <- rowMeans(edger.rpkm.sc[row.names(edger.rpkm.sc) %in% Epsilon_DEG_df$ensid, Epsilon.iv])

# pct above 0 rpkm
Epsilon_DEG_df$pct_above_0_rpkm <- rowSums(edger.rpkm.sc[row.names(edger.rpkm.sc) %in% Epsilon_DEG_df$ensid, Epsilon.iv] > 0)/length(Epsilon_cells)


#write.table(Alpha_DEG_df, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/DEG_info_Alpha_dn.txt", row.names = FALSE, sep = "\t")
#write.table(Beta_DEG_df, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/DEG_info_Beta_dn.txt", row.names = FALSE, sep = "\t")
#write.table(Gamma_DEG_df, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/DEG_info_Gamma_dn.txt", row.names = FALSE, sep = "\t")
#write.table(Delta_DEG_df, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/DEG_info_Delta_dn.txt", row.names = FALSE, sep = "\t")
#write.table(Epsilon_DEG_df, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/DEG_lists/DEG_info_Epsilon_dn.txt", row.names = FALSE, sep = "\t")

```

## Create new Targets group designating clusters

```{r}

targets$cluster <- as.character(annotated_cluster[match(targets$alias, annotated_cluster$Cell),"Cluster"])
targets[is.na(targets)] <- "other"

```

## Regenerate pearson IF zscore filtered heatmap w/ A-E
```{r}
temp <- edger.rpkm.sc[row.names(edger.rpkm.sc) %in% IF_gene_list, IF.iv] 
temp<-log2(temp+.1)
temp<-scale(t(temp), scale =T, center = T)
temp<- data.frame(t(temp))
temp[temp>3]<-3
temp[temp< -3]<- -3
df<- as.data.frame(temp[complete.cases(temp),])

## Set z-score threshold at 30
## 30 cells must have genes expressed above the mean for that row for it to be included

IF_zscore_count <- subset(IF_zscore_count, number_greater_than_0 > 30) #760

hist(IF_zscore_count$number_greater_than_0)

zscorefilt.iv <- row.names(df) %in% row.names(IF_zscore_count)

df <- df[zscorefilt.iv,]

# create colors for clusters
cols <- wes_palette("Rushmore")[3:4]
clust_col <- c(cols, wes_palette("FantasticFox")[3:5])
set.seed(004)
clust_col <- sample(clust_col)
names(clust_col) <- c("Alpha", "Beta", "Gamma", "Delta", "Epsilon")

# put into list with group for aheatmap argument
list <- list(group_col, clust_col)


#pdf("/n/projects/lwo/singlecell/R_out/Finalized_figures/IF_scRNAseq_heatmap_zfiltered_clusters.pdf", height = 30, width = 20)
#png("/n/projects/lwo/singlecell/R_out/figs_labmeeting/IF_cluster_heatmap.png", width = 13, height = 10, units = "in", res = 300)
aheatmap(df, breaks = 0, color = rev(brewer.pal(11,"RdBu")), scale="none", labCol = targets$alias[IF.iv], annColors = list, annCol=targets[IF.iv,c("group", "cluster")], distfun = "pearson")
#dev.off()



```

## PCA of Trailblazers

Use IF_DE gene list

```{r}

# PCA with all genes
temp <- edger.rpkm.sc[row.names(edger.rpkm.sc) %in% IF_gene_list,IF.iv]
temp <- log2(temp +0.01)
temp1 <- data.frame(temp[apply(temp, 1, var, na.rm=TRUE) != 0,])
common <- intersect(row.names(temp), row.names(temp1))
common.iv <- row.names(edger.rpkm.sc) %in% common

rpkm.pca<-prcomp(t(temp1),center=T,scale=T)

df <- rpkm.pca$x

nobs.factor <- sqrt(nrow(rpkm.pca$x) - 1)
d <- rpkm.pca$sdev
v <- rpkm.pca$rotation
u <- sweep(rpkm.pca$x, 2, 1 / (d * nobs.factor), FUN = '*')

df.u <- as.data.frame(sweep(u[,c(1:4)], 2, d[1:4]^0, FUN='*'))

df.u <- df.u * nobs.factor

df <- data.frame(PC1 = df.u$PC1, PC2 = df.u$PC2, PC3 = df.u$PC3, row.names = row.names(rpkm.pca$x))

table(targets$cluster)

#pdf("/n/projects/lwo/singlecell/R_out/Finalized_figures/Trailblazer_clusters_PCA.pdf", height = 10, width = 13)
#png("/n/projects/lwo/singlecell/R_out/figs_labmeeting/Trailblazer_clusters_PCA.png", res = 300, units = "in", height = 10, width = 13)
ggplot(df, aes(x = PC1, y = PC2)) + 
  geom_point(alpha = 0.8, size = 3, aes(color = targets$cluster[IF.iv], shape = targets$group[IF.iv]))+
  scale_color_manual(values = clust_col)+
  ylab("standardized PC2 (3.4% explained var.)") +
  xlab("standardized PC1 (5.0% explained var.)") +
  theme_minimal()
#dev.off()

pdf("/n/projects/lwo/singlecell/Reworked_figures/Trailblazer_clusters_PCA.pdf", height = 10, width = 13)
ggplot(df, aes(x = PC1, y = PC2)) + 
  geom_point(alpha = 0.8, size = 5, aes(color = targets$cluster[IF.iv]))+
  scale_color_manual(values = clust_col)+
  ylab("standardized PC2 (3.4% explained var.)") +
  xlab("standardized PC1 (5.0% explained var.)") +
  theme_minimal()
dev.off()

#ggbiplot(rpkm.pca, var.axes = F)
```

## Table of who is in Trailblazers

```{r}


temp <- data.frame(table(gsub("_[0-9]+", "", names(edger.rpkm.sc))))
row.names(temp)<-temp[,1]
names(temp) <- c("subgroup", "total_freq")

alpha <- data.frame(table(gsub("_[0-9]+", "", Alpha_cells)))
names(alpha) <- c("subgroup", "alpha_frequency")

temp <- merge(temp, alpha, all = F)


#write.table(temp, "/n/projects/lwo/singlecell/txt_files/Finalized_txt/pct_tables/alpha_IF_pcts.txt", sep = "\t", row.names = FALSE)

```


#####################################################################################
GO TERMS



```{r setup, include=FALSE}
genedata<-read.table("/n/data1/genomes/bowtie-index/galGal4/Ens_80/galGal4.Ens_80.genedata.txt",sep='\t',header=T,quote="",comment.char="")
alpha.dn <- read.table("/n/projects/lwo/singlecell/Finalized_txt/DEG_lists/DEG_info_Alpha_dn.txt", header = T, sep = "\t")
alpha.up <- read.table("/n/projects/lwo/singlecell/Finalized_txt/DEG_lists/DEG_info_Alpha_up.txt", header = T, sep = "\t")
```
  
subset alpha.up by % expression (40% cutoff)
  
```{r}
alpha40 <- subset(alpha.up, alpha.up$pct_above_0_rpkm >= .40)
```
  

Alpha GO

```{r}
library(clusterProfiler)

gene <- as.character(alpha40$ensid)
gene.df <- bitr(gene, fromType = "ENSEMBL",
        toType = c("ENTREZID", "SYMBOL"),
        OrgDb = org.Gg.eg.db)

## Try enrichGO

ego <- enrichGO(gene          = gene,
                universe      = genedata$GeneID,
                OrgDb         = org.Gg.eg.db,
                ont           = "BP",
                keytype = "ENSEMBL",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.05,
                qvalueCutoff  = 0.1,
                readable      = TRUE)

result.ego <- data.frame(ego@result)

ego.simplify <- simplify(ego, cutoff=0.5, by="p.adjust", select_fun=min)

result.ego.simplify <- data.frame(ego.simplify@result)

write.table(result.ego, "/n/projects/lwo/singlecell/GO_final/alpha_40_up_ego.txt", row.names = FALSE, sep = "\t")
write.table(result.ego.simplify, "/n/projects/lwo/singlecell/GO_final/alpha_40_up_ego_simplify.txt", row.names = FALSE, sep = "\t")



gene <- as.character(alpha.dn$ensid)
gene.df <- bitr(gene, fromType = "ENSEMBL",
        toType = c("ENTREZID", "SYMBOL"),
        OrgDb = org.Gg.eg.db)

## Try enrichGO

ego <- enrichGO(gene          = gene,
                universe      = genedata$GeneID,
                OrgDb         = org.Gg.eg.db,
                ont           = "BP",
                keytype = "ENSEMBL",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.05,
                qvalueCutoff  = 0.1,
                readable      = TRUE)

result.ego <- data.frame(ego@result)

ego.simplify <- simplify(ego, cutoff=0.5, by="p.adjust", select_fun=min)

result.ego.simplify <- data.frame(ego.simplify@result)

write.table(result.ego, "/n/projects/lwo/singlecell/GO_final/alpha_dn_ego.txt", row.names = FALSE, sep = "\t")
write.table(result.ego.simplify, "/n/projects/lwo/singlecell/GO_final/alpha_dn_ego_simplify.txt", row.names = FALSE, sep = "\t")



```


GO In vivo/vitro

```{r}
vivo_vitro_dn <- read.table("/n/projects/lwo/singlecell/Finalized_txt/DEG_lists/vivo_v_vitro_dn.txt", header = T, sep = "\t")
vivo_vitro_up <- read.table("/n/projects/lwo/singlecell/Finalized_txt/DEG_lists/vivo_v_vitro_up.txt", header = T, sep = "\t")



gene <- as.character(vivo_vitro_up$ensid)
gene.df <- bitr(gene, fromType = "ENSEMBL",
        toType = c("ENTREZID", "SYMBOL"),
        OrgDb = org.Gg.eg.db)

## Try enrichGO

ego <- enrichGO(gene          = gene,
                universe      = genedata$GeneID,
                OrgDb         = org.Gg.eg.db,
                ont           = "BP",
                keytype = "ENSEMBL",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.05,
                qvalueCutoff  = 0.1,
                readable      = TRUE)

result.ego <- data.frame(ego@result)

ego.simplify <- simplify(ego, cutoff=0.5, by="p.adjust", select_fun=min)

result.ego.simplify <- data.frame(ego.simplify@result)


write.table(result.ego, "/n/projects/lwo/singlecell/GO_final/vivo_vitro_up_ego.txt", row.names = FALSE, sep = "\t")
write.table(result.ego.simplify, "/n/projects/lwo/singlecell/GO_final/vivo_vitro_up_ego_simplify.txt", row.names = FALSE, sep = "\t")

################ DOWN ########################

gene <- as.character(vivo_vitro_dn$ensid)
gene.df <- bitr(gene, fromType = "ENSEMBL",
        toType = c("ENTREZID", "SYMBOL"),
        OrgDb = org.Gg.eg.db)

## Try enrichGO

ego <- enrichGO(gene          = gene,
                universe      = genedata$GeneID,
                OrgDb         = org.Gg.eg.db,
                ont           = "BP",
                keytype = "ENSEMBL",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.15,
                qvalueCutoff  = 0.15,
                readable      = TRUE)

result.ego <- data.frame(ego@result)

ego.simplify <- simplify(ego, cutoff=0.5, by="p.adjust", select_fun=min)

result.ego.simplify <- data.frame(ego.simplify@result)

### DID NOT WRITE TABLES BC ITS NOT SIGNIFICANT ###

write.table(result.ego, "/n/projects/lwo/singlecell/GO_final/vivo_vitro_dn_ego.txt", row.names = FALSE, sep = "\t")
write.table(result.ego.simplify, "/n/projects/lwo/singlecell/GO_final/vivo_vitro_dn_ego_simplify.txt", row.names = FALSE, sep = "\t")

```

Delta

```{r}
delta_up <- read.table("/n/projects/lwo/singlecell/Finalized_txt/DEG_lists/DEG_info_Delta_up.txt", header =T, sep = "\t")
delta_dn <- read.table("/n/projects/lwo/singlecell/Finalized_txt/DEG_lists/DEG_info_Delta_dn.txt", header =T, sep = "\t")

delta40 <- subset(delta_up, delta_up$pct_above_0_rpkm > 0.4)

gene <- as.character(delta40$ensid)
gene.df <- bitr(gene, fromType = "ENSEMBL",
        toType = c("ENTREZID", "SYMBOL"),
        OrgDb = org.Gg.eg.db)

## Try enrichGO

ego <- enrichGO(gene          = gene,
                universe      = genedata$GeneID,
                OrgDb         = org.Gg.eg.db,
                ont           = "BP",
                keytype = "ENSEMBL",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.05,
                qvalueCutoff  = 0.1,
                readable      = TRUE)

result.ego <- data.frame(ego@result)

ego.simplify <- simplify(ego, cutoff=0.5, by="p.adjust", select_fun=min)

result.ego.simplify <- data.frame(ego.simplify@result)

write.table(result.ego, "/n/projects/lwo/singlecell/GO_final/Delta_40_up_ego.txt", row.names = FALSE, sep = "\t")
write.table(result.ego.simplify, "/n/projects/lwo/singlecell/GO_final/Delta_40_up_ego_simplify.txt", row.names = FALSE, sep = "\t")



gene <- as.character(delta_dn$ensid)
gene.df <- bitr(gene, fromType = "ENSEMBL",
        toType = c("ENTREZID", "SYMBOL"),
        OrgDb = org.Gg.eg.db)

## Try enrichGO

ego <- enrichGO(gene          = gene,
                universe      = genedata$GeneID,
                OrgDb         = org.Gg.eg.db,
                ont           = "BP",
                keytype = "ENSEMBL",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.05,
                qvalueCutoff  = 0.1,
                readable      = TRUE)

result.ego <- data.frame(ego@result)

ego.simplify <- simplify(ego, cutoff=0.5, by="p.adjust", select_fun=min)

result.ego.simplify <- data.frame(ego.simplify@result)


```


Epsilon

```{r}
Epsilon_up <- read.table("/n/projects/lwo/singlecell/Finalized_txt/DEG_lists/DEG_info_Epsilon_up.txt", header =T, sep = "\t")
Epsilon_dn <- read.table("/n/projects/lwo/singlecell/Finalized_txt/DEG_lists/DEG_info_Epsilon_dn.txt", header =T, sep = "\t")

Epsilon40 <- subset(Epsilon_up, Epsilon_up$pct_above_0_rpkm > 0.4)

gene <- as.character(Epsilon40$ensid)
gene.df <- bitr(gene, fromType = "ENSEMBL",
        toType = c("ENTREZID", "SYMBOL"),
        OrgDb = org.Gg.eg.db)

## Try enrichGO

ego <- enrichGO(gene          = gene,
                universe      = genedata$GeneID,
                OrgDb         = org.Gg.eg.db,
                ont           = "BP",
                keytype = "ENSEMBL",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.05,
                qvalueCutoff  = 0.1,
                readable      = TRUE)

result.ego <- data.frame(ego@result)

ego.simplify <- simplify(ego, cutoff=0.5, by="p.adjust", select_fun=min)

result.ego.simplify <- data.frame(ego.simplify@result)


gene <- as.character(Epsilon_dn$ensid)
gene.df <- bitr(gene, fromType = "ENSEMBL",
        toType = c("ENTREZID", "SYMBOL"),
        OrgDb = org.Gg.eg.db)

## Try enrichGO

ego <- enrichGO(gene          = gene,
                universe      = genedata$GeneID,
                OrgDb         = org.Gg.eg.db,
                ont           = "BP",
                keytype = "ENSEMBL",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.05,
                qvalueCutoff  = 0.1,
                readable      = TRUE)

result.ego <- data.frame(ego@result)

ego.simplify <- simplify(ego, cutoff=0.5, by="p.adjust", select_fun=min)

result.ego.simplify <- data.frame(ego.simplify@result)

write.table(result.ego, "/n/projects/lwo/singlecell/GO_final/Epsilon_dn_ego.txt", row.names = FALSE, sep = "\t")
write.table(result.ego.simplify, "/n/projects/lwo/singlecell/GO_final/Epsilon_dn_ego_simplify.txt", row.names = FALSE, sep = "\t")

```

Leaders and Trailers combined

```{r}

Leader_up <- read.table("/n/projects/lwo/singlecell/Finalized_folder/text_out/DEG_lists/LeadersCombined_v_invivo_up.txt", 
                        header =T, sep = "\t")
Leader_dn <- read.table("/n/projects/lwo/singlecell/Finalized_folder/text_out/DEG_lists/LeadersCombined_v_invivo_dn.txt", 
                        header =T, sep = "\t")
table(targets$comboLeader)
Leader.iv <- grepl("Leader", names(edger.rpkm))
targets$comboLeader == "Leader"

pct_above_0_rpkm <- 
  rowSums(edger.rpkm[row.names(edger.rpkm) %in% Leader_up$ensid, Leader.iv] > 0)/length(names(edger.rpkm)[Leader.iv])

max(pct_above_0_rpkm)

gene <- as.character(Epsilon40$ensid)
gene.df <- bitr(gene, fromType = "ENSEMBL",
        toType = c("ENTREZID", "SYMBOL"),
        OrgDb = org.Gg.eg.db)

## Try enrichGO

ego <- enrichGO(gene          = gene,
                universe      = genedata$GeneID,
                OrgDb         = org.Gg.eg.db,
                ont           = "BP",
                keytype = "ENSEMBL",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.05,
                qvalueCutoff  = 0.1,
                readable      = TRUE)

result.ego <- data.frame(ego@result)

ego.simplify <- simplify(ego, cutoff=0.5, by="p.adjust", select_fun=min)

result.ego.simplify <- data.frame(ego.simplify@result)


gene <- as.character(Leader_dn$ensid)
gene.df <- bitr(gene, fromType = "ENSEMBL",
        toType = c("ENTREZID", "SYMBOL"),
        OrgDb = org.Gg.eg.db)

## Try enrichGO

ego <- enrichGO(gene          = gene,
                universe      = genedata$GeneID,
                OrgDb         = org.Gg.eg.db,
                ont           = "BP",
                keytype = "ENSEMBL",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.05,
                qvalueCutoff  = 0.1,
                readable      = TRUE)

result.ego <- data.frame(ego@result)

ego.simplify <- simplify(ego, cutoff=0.5, by="p.adjust", select_fun=min)

result.ego.simplify <- data.frame(ego.simplify@result)

write.table(result.ego, "/n/projects/lwo/singlecell/Finalized_folder/text_out/GO/Leader_dn_ego.txt", row.names = FALSE, sep = "\t")
write.table(result.ego.simplify, "/n/projects/lwo/singlecell/Finalized_folder/text_out/GO/Leader_dn_ego_simplify.txt", row.names = FALSE, sep = "\t")




```

```{r}

Trailer_up <- read.table("/n/projects/lwo/singlecell/Finalized_folder/text_out/DEG_lists/TrailersCombined_v_invivo_up.txt", 
                        header =T, sep = "\t")
Trailer_dn <- read.table("/n/projects/lwo/singlecell/Finalized_folder/text_out/DEG_lists/TrailersCombined_v_invivo_dn.txt", 
                        header =T, sep = "\t")

Trailer.iv <- grepl("Trailer", names(edger.rpkm))

pct_above_0_rpkm <- 
  rowSums(edger.rpkm[row.names(edger.rpkm) %in% Trailer_up$ensid, Trailer.iv] > 0)/length(names(edger.rpkm)[Trailer.iv])

max(pct_above_0_rpkm)

gene <- as.character(Epsilon40$ensid)
gene.df <- bitr(gene, fromType = "ENSEMBL",
        toType = c("ENTREZID", "SYMBOL"),
        OrgDb = org.Gg.eg.db)

## Try enrichGO

ego <- enrichGO(gene          = gene,
                universe      = genedata$GeneID,
                OrgDb         = org.Gg.eg.db,
                ont           = "BP",
                keytype = "ENSEMBL",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.05,
                qvalueCutoff  = 0.1,
                readable      = TRUE)

result.ego <- data.frame(ego@result)

ego.simplify <- simplify(ego, cutoff=0.5, by="p.adjust", select_fun=min)

result.ego.simplify <- data.frame(ego.simplify@result)


gene <- as.character(Trailer_dn$ensid)
gene.df <- bitr(gene, fromType = "ENSEMBL",
        toType = c("ENTREZID", "SYMBOL"),
        OrgDb = org.Gg.eg.db)

## Try enrichGO

ego <- enrichGO(gene          = gene,
                universe      = genedata$GeneID,
                OrgDb         = org.Gg.eg.db,
                ont           = "BP",
                keytype = "ENSEMBL",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.05,
                qvalueCutoff  = 0.1,
                readable      = TRUE)

result.ego <- data.frame(ego@result)

ego.simplify <- simplify(ego, cutoff=0.5, by="p.adjust", select_fun=min)

result.ego.simplify <- data.frame(ego.simplify@result)

# write.table(result.ego, "/n/projects/lwo/singlecell/Finalized_folder/text_out/GO/Trailer_dn_ego.txt", row.names = FALSE, sep = "\t")
# write.table(result.ego.simplify, "/n/projects/lwo/singlecell/Finalized_folder/text_out/GO/Trailer_dn_ego_simplify.txt", row.names = FALSE, sep = "\t")


```

## try uncombined vs in vivo

```{r}

summary <- read.table("/n/projects/lwo/singlecell/Finalized_folder/text_out/summarytable_subpop_v_all_invivo.txt", 
           header = T, sep = "\t")

L_13_DE <- subset(summary, 
                   summary$HH13_Leader_v_other.padj < 0.05 & abs(summary$HH13_Leader_v_other.log2fc) >2)
L_15_DE <- subset(summary, 
                   summary$HH15_Leader_v_other.padj < 0.05 & abs(summary$HH15_Leader_v_other.log2fc) >2)

Tr_13_DE <- subset(summary, 
                   summary$HH13_Trailer_v_other.padj < 0.05 & abs(summary$HH13_Trailer_v_other.log2fc) >2)
Tr_15_DE <- subset(summary, 
                   summary$HH15_Trailer_v_other.padj < 0.05 & abs(summary$HH15_Trailer_v_other.log2fc) >2)

## logical vectors
L15.iv <- grepl("HH15_Leader", names(edger.rpkm))
L13.iv <- grepl("HH13_Leader", names(edger.rpkm))
Tr15.iv <- grepl("HH15_Trailer", names(edger.rpkm))
Tr13.iv <- grepl("HH13_Trailer", names(edger.rpkm))

## Check percentages
df <- subset(L_13_DE, L_13_DE$HH13_Leader_v_other.log2fc > 2)
pct_above_0_rpkm <- 
  rowSums(edger.rpkm[row.names(edger.rpkm) %in% df$ensid, L13.iv] > 0)/length(names(edger.rpkm)[L13.iv])

max(pct_above_0_rpkm) #0.48
hist(pct_above_0_rpkm)


df <- subset(L_15_DE, L_15_DE$HH15_Leader_v_other.log2fc > 2)
pct_above_0_rpkm <- 
  rowSums(edger.rpkm[row.names(edger.rpkm) %in% df$ensid, L15.iv] > 0)/length(names(edger.rpkm)[L15.iv])

max(pct_above_0_rpkm) # 0.52
hist(pct_above_0_rpkm)

df <- subset(Tr_13_DE, Tr_13_DE$HH13_Trailer_v_other.log2fc > 2)
pct_above_0_rpkm <- 
  rowSums(edger.rpkm[row.names(edger.rpkm) %in% df$ensid, Tr13.iv] > 0)/length(names(edger.rpkm)[Tr13.iv])

max(pct_above_0_rpkm) #0.3
hist(pct_above_0_rpkm)

df <- subset(Tr_15_DE, Tr_15_DE$HH15_Trailer_v_other.log2fc > 2)
pct_above_0_rpkm <- 
  rowSums(edger.rpkm[row.names(edger.rpkm) %in% df$ensid,Tr15.iv] > 0)/length(names(edger.rpkm)[Tr15.iv])

max(pct_above_0_rpkm) #0.91
hist(pct_above_0_rpkm)
```

```{r}

# write.table(L_13_DE, "/n/projects/lwo/singlecell/Finalized_folder/text_out/DEG_lists/HH13_leaders_v_all_other_invivo.txt", sep = "\t", row.names = F)
# write.table(L_15_DE, "/n/projects/lwo/singlecell/Finalized_folder/text_out/DEG_lists/HH15_leaders_v_all_other_invivo.txt", sep = "\t", row.names = F)
# write.table(Tr_13_DE, "/n/projects/lwo/singlecell/Finalized_folder/text_out/DEG_lists/HH13_trailers_v_all_other_invivo.txt", sep = "\t", row.names = F)
# write.table(Tr_15_DE, "/n/projects/lwo/singlecell/Finalized_folder/text_out/DEG_lists/HH15_trailers_v_all_other_invivo.txt", sep = "\t", row.names = F)

```

```{r}
interest.up <- subset(L_15_DE, L_15_DE$HH15_Leader_v_other.log2fc > 0)
interest.dn <- subset(L_15_DE, L_15_DE$HH15_Leader_v_other.log2fc < 0)

L15.iv <- grepl("HH15_Leader", names(edger.rpkm))

pct_above_0_rpkm <- 
  rowSums(edger.rpkm[row.names(edger.rpkm) %in% interest.up$ensid, L15.iv] > 0)/length(names(edger.rpkm)[L15.iv])

max(pct_above_0_rpkm)
hist(pct_above_0_rpkm)

df<-data.frame(pct_above_0_rpkm)
df <- subset(df, df$pct_above_0_rpkm >= 0.4)

gene <- as.character(interest.up$ensid)
gene.df <- bitr(gene, fromType = "ENSEMBL",
        toType = c("ENTREZID", "SYMBOL"),
        OrgDb = org.Gg.eg.db)

## enrichGO

ego <- enrichGO(gene          = gene,
                universe      = genedata$GeneID,
                OrgDb         = org.Gg.eg.db,
                ont           = "BP",
                keytype = "ENSEMBL",
                pAdjustMethod = "BH",
                pvalueCutoff  = .05,
                qvalueCutoff  = 0.1,
                readable      = TRUE)

result.ego <- data.frame(ego@result)
 
ego.simplify <- simplify(ego, cutoff=0.5, by="p.adjust", select_fun=min)
 
result.ego.simplify <- data.frame(ego.simplify@result)
 
write.table(result.ego,
            "/n/projects/lwo/singlecell/Finalized_folder/text_out/GO/HH15_Leader_ego_up.txt",
            row.names = FALSE, sep = "\t")

write.table(result.ego.simplify,
            "/n/projects/lwo/singlecell/Finalized_folder/text_out/GO/HH15_Leader_ego_simplified_up.txt",
            row.names = FALSE, sep = "\t")

gene <- as.character(interest.dn$ensid)
gene.df <- bitr(gene, fromType = "ENSEMBL",
        toType = c("ENTREZID", "SYMBOL"),
        OrgDb = org.Gg.eg.db)

## enrichGO

ego <- enrichGO(gene          = gene,
                universe      = genedata$GeneID,
                OrgDb         = org.Gg.eg.db,
                ont           = "BP",
                keytype = "ENSEMBL",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.05,
                qvalueCutoff  = 0.1,
                readable      = TRUE)

result.ego <- data.frame(ego@result)

ego.simplify <- simplify(ego, cutoff=0.5, by="p.adjust", select_fun=min)

result.ego.simplify <- data.frame(ego.simplify@result)

write.table(result.ego,
            "/n/projects/lwo/singlecell/Finalized_folder/text_out/GO/HH15_Leader_ego_dn.txt",
            row.names = FALSE, sep = "\t")

write.table(result.ego.simplify,
            "/n/projects/lwo/singlecell/Finalized_folder/text_out/GO/HH15_Leader_ego_simplified_dn.txt",
            row.names = FALSE, sep = "\t")

```


############################################################################################################
VOLCANO PLOTS

## Set up

```{r setup, include=FALSE}

library("ggplot2")
library("dplyr")
library("RColorBrewer")
library("plotly")
library("htmlwidgets")

summary_clusters <- read.table("/n/projects/lwo/singlecell/R_out/Finalized_txt/summary_tables/summary_clusters_v_invivo_DE.txt", sep = "\t", header = T)

```

## Pull out Trailblazer DEG

```{r}
summary_clusters <- tbl_df(summary_clusters)

trailblazer <- summary_clusters %>% 
  dplyr::select(ensid, Name, Alpha_v_other.log2fc, Alpha_v_other.padj) %>% 
  filter(abs(Alpha_v_other.log2fc) > 2 & Alpha_v_other.padj < 0.05) %>%
  mutate(log_padj = -log(Alpha_v_other.padj))

trailblazer$sign <- ifelse(trailblazer$Alpha_v_other.log2fc > 0, "pos", "neg")

trailblazer$log_padj <- ifelse(trailblazer$log_padj < 8, 8, trailblazer$log_padj)
trailblazer$log_padj <- ifelse(trailblazer$log_padj > 45, 45, trailblazer$log_padj)
trailblazer$Alpha_v_other.log2fc <- ifelse(trailblazer$Alpha_v_other.log2fc > 7, 7, trailblazer$Alpha_v_other.log2fc)
trailblazer$Alpha_v_other.log2fc <- ifelse(trailblazer$Alpha_v_other.log2fc < -7, -7, trailblazer$Alpha_v_other.log2fc)

cols <- brewer.pal(3, "Set1")
cols <- rev(cols[1:2])

#pdf("/n/projects/lwo/singlecell/Finalized_folder/pdf_out/trailblazer_logfc_v_logpadj.pdf", height = 10, width = 14)
ggplot(trailblazer, aes(x = Alpha_v_other.log2fc, y = log_padj, color = sign)) +
  geom_point(alpha = 0.5, size = 5) +
  scale_color_manual(values = cols) +
  theme_classic()
#dev.off()

trailblazer <- mutate(trailblazer, name_ensid = paste(trailblazer$Name, trailblazer$ensid, sep = "   "))

g <- ggplot(trailblazer, aes(x = Alpha_v_other.log2fc, y = log_padj, color = sign)) +
  geom_point(alpha = 0.5, size = 3, aes(text = name_ensid)) +
  scale_color_manual(values = cols) +
  theme_classic()

gg <- ggplotly(g)

saveWidget(gg, "/n/projects/lwo/singlecell/Finalized_folder/pdf_out/trailblazer_logfc_v_logpadj.html")
```

```{r}

epsilon <- summary_clusters %>% 
  dplyr::select(ensid, Name, Epsilon_v_other.log2fc, Epsilon_v_other.padj) %>% 
  filter(abs(Epsilon_v_other.log2fc) > 2 & Epsilon_v_other.padj < 0.05) %>%
  mutate(log_padj = -log(Epsilon_v_other.padj))

epsilon$sign <- ifelse(epsilon$Epsilon_v_other.log2fc > 0, "pos", "neg")

length(epsilon[epsilon$Epsilon_v_other.log2fc < 0,1])

epsilon$log_padj <- ifelse(epsilon$log_padj < 8, 8, epsilon$log_padj)
epsilon$log_padj <- ifelse(epsilon$log_padj > 45, 45, epsilon$log_padj)
epsilon$Epsilon_v_other.log2fc <- ifelse(epsilon$Epsilon_v_other.log2fc > 7, 7, epsilon$Epsilon_v_other.log2fc)
epsilon$Epsilon_v_other.log2fc <- ifelse(epsilon$Epsilon_v_other.log2fc < -7, -7, epsilon$Epsilon_v_other.log2fc)

cols <- brewer.pal(3, "Set1")
cols <- rev(cols[1:2])

pdf("/n/projects/lwo/singlecell/Finalized_folder/pdf_out/epsilon_logfc_v_logpadj.pdf", height = 10, width = 14)
ggplot(epsilon, aes(x = Epsilon_v_other.log2fc, y = log_padj, color = sign)) +
  geom_point(alpha = 0.5, size = 5) +
  scale_color_manual(values = cols) +
  theme_classic()
dev.off()

epsilon <- mutate(epsilon, name_ensid = paste(epsilon$Name, epsilon$ensid, sep = "   "))

g <- ggplot(epsilon, aes(x = Epsilon_v_other.log2fc, y = log_padj, color = sign)) +
  geom_point(alpha = 0.5, size = 3, aes(text = name_ensid)) +
  scale_color_manual(values = cols) +
  theme_classic()

gg <- ggplotly(g)

saveWidget(gg, "/n/projects/lwo/singlecell/Finalized_folder/pdf_out/epsilon_logfc_v_logpadj.html")

```


```{r}

delta <- summary_clusters %>% 
  dplyr::select(ensid, Name, Delta_v_other.log2fc, Delta_v_other.padj) %>% 
  filter(abs(Delta_v_other.log2fc) > 2 & Delta_v_other.padj < 0.05) %>%
  mutate(log_padj = -log(Delta_v_other.padj))

delta$sign <- ifelse(delta$Delta_v_other.log2fc > 0, "pos", "neg")

length(delta[delta$Delta_v_other.log2fc < 0,1])

delta$log_padj <- ifelse(delta$log_padj < 8, 8, delta$log_padj)
delta$log_padj <- ifelse(delta$log_padj > 45, 45, delta$log_padj)
delta$Delta_v_other.log2fc <- ifelse(delta$Delta_v_other.log2fc > 7, 7, delta$Delta_v_other.log2fc)
delta$Delta_v_other.log2fc <- ifelse(delta$Delta_v_other.log2fc < -7, -7, delta$Delta_v_other.log2fc)

cols <- brewer.pal(3, "Set1")
cols <- rev(cols[1:2])

pdf("/n/projects/lwo/singlecell/Finalized_folder/pdf_out/delta_logfc_v_logpadj.pdf", height = 10, width = 14)
ggplot(delta, aes(x = Delta_v_other.log2fc, y = log_padj, color = sign)) +
  geom_point(alpha = 0.5, size = 5) +
  scale_color_manual(values = cols) +
  theme_classic()
dev.off()

delta <- mutate(delta, name_ensid = paste(delta$Name, delta$ensid, sep = "   "))

g <- ggplot(delta, aes(x = Delta_v_other.log2fc, y = log_padj, color = sign)) +
  geom_point(alpha = 0.5, size = 3, aes(text = name_ensid)) +
  scale_color_manual(values = cols) +
  theme_classic()

gg <- ggplotly(g)

saveWidget(gg, "/n/projects/lwo/singlecell/Finalized_folder/pdf_out/delta_logfc_v_logpadj.html")

```

```{r}
summary_combo <-
  read.table("/n/projects/lwo/singlecell/Finalized_folder/text_out/summary_comboLeader_comboTrailer_v_invivo_DE.txt",
           header = T, sep = "\t")


Leader <- summary_combo %>% 
  dplyr::select(ensid, Name, Leader_v_other.log2fc, Leader_v_other.padj) %>% 
  filter(abs(Leader_v_other.log2fc) > 2 & Leader_v_other.padj < 0.05) %>%
  mutate(log_padj = -log(Leader_v_other.padj))

Leader$sign <- ifelse(Leader$Leader_v_other.log2fc > 0, "pos", "neg")

Leader$log_padj <- ifelse(Leader$log_padj < 5, 5, Leader$log_padj)
Leader$log_padj <- ifelse(Leader$log_padj > 45, 45, Leader$log_padj)
Leader$Leader_v_other.log2fc <- ifelse(Leader$Leader_v_other.log2fc > 7, 7, Leader$Leader_v_other.log2fc)
Leader$Leader_v_other.log2fc <- ifelse(Leader$Leader_v_other.log2fc < -7, -7, Leader$Leader_v_other.log2fc)

cols <- brewer.pal(3, "Set1")
cols <- rev(cols[1:2])

pdf("/n/projects/lwo/singlecell/Finalized_folder/pdf_out/Leadercombo_logfc_v_logpadj.pdf", height = 10, width = 14)
ggplot(Leader, aes(x = Leader_v_other.log2fc, y = log_padj, color = sign)) +
  geom_point(alpha = 0.5, size = 5) +
  scale_color_manual(values = cols) +
  theme_classic()
dev.off()

Leader <- mutate(Leader, name_ensid = paste(Leader$Name, Leader$ensid, sep = "   "))

g <- ggplot(Leader, aes(x = Leader_v_other.log2fc, y = log_padj, color = sign)) +
  geom_point(alpha = 0.5, size = 3, aes(text = name_ensid)) +
  scale_color_manual(values = cols) +
  theme_classic()

gg <- ggplotly(g)

saveWidget(gg, "/n/projects/lwo/singlecell/Finalized_folder/pdf_out/Leadercombo_logfc_v_logpadj.html")


```


```{r}

Trailer <- summary_combo %>% 
  dplyr::select(ensid, Name, Trailer_v_other.log2fc, Trailer_v_other.padj) %>% 
  filter(abs(Trailer_v_other.log2fc) > 2 & Trailer_v_other.padj < 0.05) %>%
  mutate(log_padj = -log(Trailer_v_other.padj))

Trailer$sign <- ifelse(Trailer$Trailer_v_other.log2fc > 0, "pos", "neg")

Trailer$log_padj <- ifelse(Trailer$log_padj < 5, 5, Trailer$log_padj)
Trailer$log_padj <- ifelse(Trailer$log_padj > 45, 45, Trailer$log_padj)
Trailer$Trailer_v_other.log2fc <- ifelse(Trailer$Trailer_v_other.log2fc > 7, 7, Trailer$Trailer_v_other.log2fc)
Trailer$Trailer_v_other.log2fc <- ifelse(Trailer$Trailer_v_other.log2fc < -7, -7, Trailer$Trailer_v_other.log2fc)

cols <- brewer.pal(3, "Set1")
cols <- rev(cols[1:2])

pdf("/n/projects/lwo/singlecell/Finalized_folder/pdf_out/Trailercombo_logfc_v_logpadj.pdf", height = 10, width = 14)
ggplot(Trailer, aes(x = Trailer_v_other.log2fc, y = log_padj, color = sign)) +
  geom_point(alpha = 0.5, size = 5) +
  scale_color_manual(values = cols) +
  theme_classic()
dev.off()

Trailer <- mutate(Trailer, name_ensid = paste(Trailer$Name, Trailer$ensid, sep = "   "))

g <- ggplot(Trailer, aes(x = Trailer_v_other.log2fc, y = log_padj, color = sign)) +
  geom_point(alpha = 0.5, size = 3, aes(text = name_ensid)) +
  scale_color_manual(values = cols) +
  theme_classic()

gg <- ggplotly(g)

saveWidget(gg, "/n/projects/lwo/singlecell/Finalized_folder/pdf_out/Trailercombo_logfc_v_logpadj.html")


```

## Subpopulations vs all other in vivo

```{r}
df <- read.table("/n/projects/lwo/singlecell/Finalized_folder/text_out/DEG_lists/HH15_trailers_v_all_other_invivo.txt", header = T, sep = "\t")



HH15_trailer <- df %>% 
  dplyr::select(ensid, Name, HH15_Trailer_v_other.log2fc, HH15_Trailer_v_other.padj) %>% 
  filter(abs(HH15_Trailer_v_other.log2fc) > 2 & HH15_Trailer_v_other.padj < 0.05) %>%
  mutate(log_padj = -log(HH15_Trailer_v_other.padj))

HH15_trailer$sign <- ifelse(HH15_trailer$HH15_Trailer_v_other.log2fc > 0, "pos", "neg")

length(HH15_trailer[HH15_trailer$HH15_Trailer_v_other.log2fc < 0,1])

HH15_trailer$log_padj <- ifelse(HH15_trailer$log_padj < 5, 5, HH15_trailer$log_padj)
HH15_trailer$log_padj <- ifelse(HH15_trailer$log_padj > 45, 45, HH15_trailer$log_padj)
HH15_trailer$HH15_Trailer_v_other.log2fc <- ifelse(HH15_trailer$HH15_Trailer_v_other.log2fc > 7, 7, HH15_trailer$HH15_Trailer_v_other.log2fc)
HH15_trailer$HH15_Trailer_v_other.log2fc <- ifelse(HH15_trailer$HH15_Trailer_v_other.log2fc < -7, -7, HH15_trailer$HH15_Trailer_v_other.log2fc)

cols <- brewer.pal(3, "Set1")
cols <- rev(cols[1:2])

pdf("/n/projects/lwo/singlecell/Finalized_folder/pdf_out/HH15_trailers_v_all_other_invivo_logfc_v_logpadj.pdf", height = 10, width = 14)
ggplot(HH15_trailer, aes(x = HH15_Trailer_v_other.log2fc, y = log_padj, color = sign)) +
  geom_point(alpha = 0.5, size = 5) +
  scale_color_manual(values = cols) +
  theme_classic()
dev.off()

HH15_trailer <- mutate(HH15_trailer, name_ensid = paste(HH15_trailer$Name, HH15_trailer$ensid, sep = "   "))

g <- ggplot(HH15_trailer, aes(x = HH15_Trailer_v_other.log2fc, y = log_padj, color = sign)) +
  geom_point(alpha = 0.5, size = 3, aes(text = name_ensid)) +
  scale_color_manual(values = cols) +
  theme_classic()

gg <- ggplotly(g)

saveWidget(gg, "/n/projects/lwo/singlecell/Finalized_folder/pdf_out/HH15_trailers_v_all_other_invivo_logfc_v_logpadj.html")


```


## Subpops vs all other invivo
```{r}
df <- read.table("/n/projects/lwo/singlecell/Finalized_folder/text_out/DEG_lists/HH13_trailers_v_all_other_invivo.txt", header = T, sep = "\t")



HH13_trailer <- df %>% 
  dplyr::select(ensid, Name, HH13_Trailer_v_other.log2fc, HH13_Trailer_v_other.padj) %>% 
  filter(abs(HH13_Trailer_v_other.log2fc) > 2 & HH13_Trailer_v_other.padj < 0.05) %>%
  mutate(log_padj = -log(HH13_Trailer_v_other.padj))

HH13_trailer$sign <- ifelse(HH13_trailer$HH13_Trailer_v_other.log2fc > 0, "pos", "neg")

length(HH13_trailer[HH13_trailer$HH13_Trailer_v_other.log2fc < 0,1])

HH13_trailer$log_padj <- ifelse(HH13_trailer$log_padj < 5, 5, HH13_trailer$log_padj)
HH13_trailer$log_padj <- ifelse(HH13_trailer$log_padj > 45, 45, HH13_trailer$log_padj)
HH13_trailer$HH13_Trailer_v_other.log2fc <- ifelse(HH13_trailer$HH13_Trailer_v_other.log2fc > 7, 7, HH13_trailer$HH13_Trailer_v_other.log2fc)
HH13_trailer$HH13_Trailer_v_other.log2fc <- ifelse(HH13_trailer$HH13_Trailer_v_other.log2fc < -7, -7, HH13_trailer$HH13_Trailer_v_other.log2fc)

cols <- brewer.pal(3, "Set1")
cols <- rev(cols[1:2])

pdf("/n/projects/lwo/singlecell/Finalized_folder/pdf_out/HH13_trailers_v_all_other_logfc_v_logpadj.pdf", height = 10, width = 14)
ggplot(HH13_trailer, aes(x = HH13_Trailer_v_other.log2fc, y = log_padj, color = sign)) +
  geom_point(alpha = 0.5, size = 5) +
  scale_color_manual(values = cols) +
  theme_classic()
dev.off()

HH13_trailer <- mutate(HH13_trailer, name_ensid = paste(HH13_trailer$Name, HH13_trailer$ensid, sep = "   "))

g <- ggplot(HH13_trailer, aes(x = HH13_Trailer_v_other.log2fc, y = log_padj, color = sign)) +
  geom_point(alpha = 0.5, size = 3, aes(text = name_ensid)) +
  scale_color_manual(values = cols) +
  theme_classic()

gg <- ggplotly(g)

saveWidget(gg, "/n/projects/lwo/singlecell/Finalized_folder/pdf_out/HH13_trailers_v_all_other_logfc_v_logpadj.html")


```

```{r}
df <- read.table("/n/projects/lwo/singlecell/Finalized_folder/text_out/DEG_lists/HH15_leaders_v_all_other_invivo.txt", header = T, sep = "\t")



HH15_Leader<- df %>% 
  dplyr::select(ensid, Name, HH15_Leader_v_other.log2fc, HH15_Leader_v_other.padj) %>% 
  filter(abs(HH15_Leader_v_other.log2fc) > 2 & HH15_Leader_v_other.padj < 0.05) %>%
  mutate(log_padj = -log(HH15_Leader_v_other.padj))

HH15_Leader$sign <- ifelse(HH15_Leader$HH15_Leader_v_other.log2fc > 0, "pos", "neg")

length(HH15_Leader[HH15_Leader$HH15_Leader_v_other.log2fc < 0,1])

HH15_Leader$log_padj <- ifelse(HH15_Leader$log_padj < 5, 5, HH15_Leader$log_padj)
HH15_Leader$log_padj <- ifelse(HH15_Leader$log_padj > 45, 45, HH15_Leader$log_padj)
HH15_Leader$HH15_Leader_v_other.log2fc <- ifelse(HH15_Leader$HH15_Leader_v_other.log2fc > 7, 7, HH15_Leader$HH15_Leader_v_other.log2fc)
HH15_Leader$HH15_Leader_v_other.log2fc <- ifelse(HH15_Leader$HH15_Leader_v_other.log2fc < -7, -7, HH15_Leader$HH15_Leader_v_other.log2fc)

cols <- brewer.pal(3, "Set1")
cols <- rev(cols[1:2])

pdf("/n/projects/lwo/singlecell/Finalized_folder/pdf_out/HH15_leaders_v_all_other_invivo_logfc_v_logpadj.pdf", height = 10, width = 14)
ggplot(HH15_Leader, aes(x = HH15_Leader_v_other.log2fc, y = log_padj, color = sign)) +
  geom_point(alpha = 0.5, size = 5) +
  scale_color_manual(values = cols) +
  theme_classic()
dev.off()

HH15_Leader <- mutate(HH15_Leader, name_ensid = paste(HH15_Leader$Name, HH15_Leader$ensid, sep = "   "))

g <- ggplot(HH15_Leader, aes(x = HH15_Leader_v_other.log2fc, y = log_padj, color = sign)) +
  geom_point(alpha = 0.5, size = 3, aes(text = name_ensid)) +
  scale_color_manual(values = cols) +
  theme_classic()

gg <- ggplotly(g)

saveWidget(gg, "/n/projects/lwo/singlecell/Finalized_folder/pdf_out/HH15_leaders_v_all_other_invivo_logfc_v_logpadj.html")


```

```{r}
df <- read.table("/n/projects/lwo/singlecell/Finalized_folder/text_out/DEG_lists/HH13_leaders_v_all_other_invivo.txt", header = T, sep = "\t")



HH13_Leader<- df %>% 
  dplyr::select(ensid, Name, HH13_Leader_v_other.log2fc, HH13_Leader_v_other.padj) %>% 
  filter(abs(HH13_Leader_v_other.log2fc) > 2 & HH13_Leader_v_other.padj < 0.05) %>%
  mutate(log_padj = -log(HH13_Leader_v_other.padj))

HH13_Leader$sign <- ifelse(HH13_Leader$HH13_Leader_v_other.log2fc > 0, "pos", "neg")

length(HH13_Leader[HH13_Leader$HH13_Leader_v_other.log2fc < 0,1])

HH13_Leader$log_padj <- ifelse(HH13_Leader$log_padj < 5, 5, HH13_Leader$log_padj)
HH13_Leader$log_padj <- ifelse(HH13_Leader$log_padj > 45, 45, HH13_Leader$log_padj)
HH13_Leader$HH13_Leader_v_other.log2fc <- ifelse(HH13_Leader$HH13_Leader_v_other.log2fc > 7, 7, HH13_Leader$HH13_Leader_v_other.log2fc)
HH13_Leader$HH13_Leader_v_other.log2fc <- ifelse(HH13_Leader$HH13_Leader_v_other.log2fc < -7, -7, HH13_Leader$HH13_Leader_v_other.log2fc)

cols <- brewer.pal(3, "Set1")
cols <- rev(cols[1:2])

pdf("/n/projects/lwo/singlecell/Finalized_folder/pdf_out/HH13_leaders_v_all_other_invivo_logfc_v_logpadj.pdf", height = 10, width = 14)
ggplot(HH13_Leader, aes(x = HH13_Leader_v_other.log2fc, y = log_padj, color = sign)) +
  geom_point(alpha = 0.5, size = 5) +
  scale_color_manual(values = cols) +
  theme_classic()
dev.off()

HH13_Leader <- mutate(HH13_Leader, name_ensid = paste(HH13_Leader$Name, HH13_Leader$ensid, sep = "   "))

g <- ggplot(HH13_Leader, aes(x = HH13_Leader_v_other.log2fc, y = log_padj, color = sign)) +
  geom_point(alpha = 0.5, size = 3, aes(text = name_ensid)) +
  scale_color_manual(values = cols) +
  theme_classic()

gg <- ggplotly(g)

saveWidget(gg, "/n/projects/lwo/singlecell/Finalized_folder/pdf_out/HH13_leaders_v_all_other_invivo_logfc_v_logpadj.html")

```
