MEP-LINCS_SS1_Analysis.Rmd

---
output: html_document
---

```{r global_options, include=FALSE}
knitr::opts_chunk$set(fig.width=12, fig.height=8,
                      echo=FALSE, warning=FALSE, message=FALSE)
```


```{r }
#Author: Mark Dane, copyright 2015

source("MEPLINCSFunctions.R")

#Set the staining set to be analyzed (SS1|SS2|SS3)
ss <- "SS1"
#Set the cell line to be analyzed (PC3|MCF7|YAPC)
cellLine <- "MCF7"
```



```{r setup}
library("ggplot2")
library("data.table")
library("MEMA")
library("grid")
library(knitr)
library("gplots")
library("RColorBrewer")

  
#Setup colors for Barcode and text in all heatmaps
selDark2 <- colorRampPalette(brewer.pal(8,"Dark2"))
plateCol = selDark2(8)

l1 <- fread(paste0("./",ss,"/AnnotatedData/",cellLine,"_",ss,"_Level1.txt"), showProgress = FALSE)
l2 <- fread(paste0("./",ss,"/AnnotatedData/",cellLine,"_",ss,"_Level2.txt"), showProgress = FALSE)
l3 <- fread(paste0("./",ss,"/AnnotatedData/",cellLine,"_",ss,"_Level3.txt"), showProgress = FALSE)
l4 <- fread(paste0("./",ss,"/AnnotatedData/",cellLine,"_",ss,"_Level4.txt"), showProgress = FALSE)

#Work around issue in fread about logicals read as characters
l1$Spot_PA_Perimeter <- as.logical(gsub(" ","",l1$Spot_PA_Perimeter))
l1$Spot_PA_Sparse <- as.logical(gsub(" ","",l1$Spot_PA_Sparse))
l1$Spot_PA_OuterCell <- as.logical(gsub(" ","",l1$Spot_PA_OuterCell))

barcodes <- sort(unique(l3$Barcode))

#Set a threshold for filtering wells on their QA score
wellQAThresh <- 0.7

#TODO: Read this from Level 3 data
lthresh <- 0.6

#Number of PCS components to use
nrPCs <- 9

```
---
title: "MEP-LINCS `r cellLine` `r ss` Pilot Analysis"
date: "`r Sys.Date()`"
output: pdf_document
---

##Summary
The MEP-LINCS `r cellLine` `r ss` datasets include four levels of high content imaging data from on Microenvironment Microarrays (MEMAs). After QA filtering, there are `r nrow(l4)` Microenvironment Perturbations (MEPs) that are pairwise combinations of `r length(unique(l4$ECMpAnnotID))` printed ECM proteins and `r length(unique(l4$LigandAnnotID))-1` ligands or growth factors.  

##Introduction
The LINCS Pilot `r cellLine` `r ss` experiment was performed with cells grown in eight 8-well plates. The `r ss` staining set includes, DAPI, `r unique(l4$Endpoint488)` (488nm), `r unique(l4$Endpoint555)` (555nm) and `r unique(l4$Endpoint647)` (647nm). Four color images of the cells at each spot were gathered on an Nikon automated microscope. All data for this staining set comes from the nuclei as defined by the DAPI staining.

Intensity, position and a limited set of morphology data are gatherd for each cell, merged with the experiment metadata, normalized and aggregated. The dataset is organized to the four LINCS imaging categories as follows:

Level 1 - Raw data  
Level 2 - Normalized data  
Level 3 - Normalized data aggregated to the spot level  
Level 4 - Normalized data aggregated to the replicate (MEP) level  

The data merging and analysis is done in R using open source software. 

\newpage

##QA Scoring of the dataset
Each well is scored for even cell seeding according to the count of the DAPI-stained nuclei. A detailed explanation of the QA method is in the supplemental material. In brief, the level 2 and 3 data have cell counts at the spot level and locally-averaged cell counts at the neighborhood level. Both of these parameters are used to score the wells and filter the dataset. QA Scores range from 0 to 1 and represent the proportion of the spots that have at least one cell and are not in low cell count neighborhoods.

The following plots are pseudoimages each MEMA's spot cell count and a histogram of the loess model used for QA scoring. 


```{r Heatmaps_QAScores, echo=FALSE, fig.width=3.7,fig.height=4, eval=TRUE}

plotSCCHeatmapsQAHistograms(l3, barcodes)

```



```{r Filtering, echo=FALSE}

#Remove failed QA wells
l3F <- l3[!l3$QA_LowWellQA]

#Remove the fiducial and blank data
setkey(l3F,ECMp)
l3F <- l3F[!"fiducial"]
l3F <- l3F[!"blank"]

#After filtering, summarize spot level data to MEP level by taking 
#the medians of the parameters

mepNames<-grep("AreaShape|Children|Intensity|Spot_PA_SpotCellCount|Loess|EdUPositiveProportion|Population|Nuclei_PA_AreaShape_Neighbors|Spot_PA_ReplicateCount|Cytoplasm_PA_Intensity_LineageRatio|LigandAnnotID|ECMpAnnotID", x=names(l3F),value=TRUE)

mepKeep<-l3F[,mepNames,with=FALSE]
l4F<-mepKeep[,lapply(.SD,numericMedian),keyby="LigandAnnotID,ECMpAnnotID"]
  
#Merge back in the replicate metadata
mDT <- l3F[,list(CellLine,Ligand,Barcode,Endpoint488,Endpoint555,Endpoint647,EndpointDAPI,ECMp,MEP,Spot_PA_ReplicateCount),keyby="LigandAnnotID,ECMpAnnotID"]
l4F <- mDT[l4F, mult="first"]

#Add Robust Z Scores of the normalized Spot Cell Counts across
#the entire staining set
l4F <- l4F[,Spot_PA_SpotCellCount_RZSNorm_RobustZ := RZScore(Spot_PA_SpotCellCount_RZSNorm)]

if (ss == "SS2"){
  #Add Robust Z Score of the Edu Signal
  l4F <- l4F[,EduPositiveProportion_RZSNorm_RobustZ := RZScore(EduPositiveProportion_RZSNorm)]
}

if (ss == "SS3"){
  #Add Robust Z Scores of the normalized lineage ratios
  l4F <- l4F[, Cytoplasm_PA_Intensity_LineageRatio_RZSNorm_RobustZ:= RZScore(Cytoplasm_PA_Intensity_LineageRatio_RZSNorm)]
  
}

```

\newpage

##Filtering
Wells with QA scores below `r wellQAThresh ` and the FBS control wells are removed from further analysis of the dataset. After filtering on the well QA score there are `r length(unique(l4F$LigandAnnotID))-1` ligands in the dataset. 

Each spot represents a MEP that is a pairwise combination of the ECM protein printed at the spot and the ligand in the solution of the well. The number of replicate MEPs after removing low-quality wells are shown in the supplemental material.

##Normalization
The raw data was normalized to the control well in each plate as follows:

For each feature, the median and median absolute deviation (MAD) of each plate's control well were calculated.
Feature values for all cells in a plate were normalized by subtracting the plate's control well median and dividing by the plate's control well MAD*1.48.

##Spot Cell Count Analysis
The spot cell count analysis uses robust Z scores to identify MEPs with extreme population sizes. The normalized spot cell counts are summarized by the median of their replicates. The median and mad of the distribution of normalized and summarized values are used to convert to robust Z scores and are shown below. The blue lines at +/- 2 show thresholds for selecting MEPs of interest. Below the distribution plot are plots with Z scores stratified by ligand and ECM protein. A listing of the MEPs outside of the blue lines is in the supplemental material.

```{r SCCRobustZScores, echo=FALSE, fig.width=8, fig.height=3.5, eval=TRUE}

plotSCCRobustZScores(l4F)
```


```{r SCC_response, fig.width=8,  fig.height=6, echo=FALSE}

p <- ggplot(l4F, aes(x = reorder(Ligand, Spot_PA_SpotCellCount_RZSNorm_RobustZ, FUN=median), y = Spot_PA_SpotCellCount_RZSNorm_RobustZ, colour = Barcode))+geom_boxplot()+
  ggtitle(paste("\n\nMEP Normalized Spot Cell Count Robust Z Scores by Ligand"))+
  xlab("")+ylab("")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=rel(.8)), axis.text.y = element_text(angle = 0, vjust = 0.5, hjust=1, size=rel(1)), plot.title = element_text(size = rel(1)),legend.text=element_text(size = rel(.3)),legend.title=element_text(size = rel(.3)))
suppressWarnings(print(p))
```

```{r, fig.width=8,  fig.height=3.5, echo=FALSE}
p <- ggplot(l4F, aes(x = reorder(ECMp, Spot_PA_SpotCellCount_RZSNorm_RobustZ, FUN=median), y = Spot_PA_SpotCellCount_RZSNorm_RobustZ))+geom_boxplot()+
  ggtitle(paste("\n\nMEP Normalized Spot Cell Count Robust Z Scores by ECM Protein"))+
  xlab("")+ylab("")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=rel(.8)), axis.text.y = element_text(angle = 0, vjust = 0.5, hjust=1, size=rel(1)), plot.title = element_text(size = rel(1)),legend.text=element_text(size = rel(.3)),legend.title=element_text(size = rel(.3)))
suppressWarnings(print(p))

```


\newpage

##Unbiased Analysis

The unbiased analysis approach is:  
  Normalize the cell level data using the median and MAD of each plate's control well
  Median summarize the normalized data to the spot level
  Median summarize the spot level data to the MEP level  
  Select a curated feature vector of MEP intensities and morphologies  
  Calculate the euclidean distance of each MEP's feature vector to the control feature vector  
  Perform hierarchical clustering on the MEPs 
  
The following heatmap shows the normalized MEP responses based on the curated feature vector normalized dataset without NID1 (an outlier ECM protein). The row color bars and text colors identify the 8 different experiment plates and help show plate batch effects.
  
```{r, fig.height=16}

lowQALigands <- unique(l3$Ligand[l3$QA_LowWellQA])
fvDT <- filterl4(l4, lowQALigands)

#Shorten the feature names
setnames(fvDT,grep("Barcode|MEP",colnames(fvDT), value = TRUE, invert = TRUE),gsub(".*_","",gsub("_RZSNorm","",grep("Barcode|MEP",colnames(fvDT), value = TRUE, invert = TRUE))))
#Remove when not needed
if("2NProportion" %in% colnames(fvDT)) setnames(fvDT,"2NProportion","DNA2NProportion")

```

The next heatmap shows the same analysis run on the dataset after removing NID1 and using data normalized to each plate's control well.
  
```{r, fig.height=16, eval = TRUE}

#Plot the heat map without NID1 in the dataset
activeMEPs <- ubHeatmap(fvDT[!grepl("NID1",fvDT$MEP)], title = "Normalized MEP Feature Vectors", cols = plateCol, activeThresh = 0 )

```


###PCA Unbiased Analysis  
The unbiased analysis is extended by using Principal Component Analysis (PCA) to reduce the dimensions of the feature vector.

The unbiased PCA analysis method is:
  Use PCA to transform the entire normalized feature set to a much smaller number of principal components (PCFV)
  Calculate the euclidean distance of each MEP's PCFV to the FBS control PCFV 
  Exclude inactive MEPs by selecting for distances above a selected percentile  
  Perform hierarchical clustering of the first `r nrPCs` PCFVs of the active MEPs 
  

```{r, fig.height=12}
l4F <- l4F[!grepl("NID1",l4F$ECMp)]

PCFVNames <- grep("RZSNorm",colnames(l4F),value=TRUE)
PCFVNames <- grep("Euler|MaximumRadius|_SE|_RobustZ",PCFVNames,invert = TRUE, value = TRUE)

#Get the PCFVs
l4PCAModel <- prcomp(as.formula(paste(" ~ ", paste(PCFVNames, collapse = "+"))), data = l4F, scale= TRUE)

#Extract the first PCs of each MEP into a matrix
l4pcvDT <- data.table(l4PCAModel$x[,1:nrPCs],Barcode =l4F$Barcode, MEP = l4F$MEP)

#Plot the heat map of the entire dataset
activeMEPs <- ubHeatmap(l4pcvDT, title = "MEP PCA Vectors", cols = plateCol, activeThresh = 0 )

```

```{r, fig.height=3, fig.width=4}

var <- l4PCAModel$sd[1:nrPCs]^2
var.percent <- var/sum(var) * 100
barplot(var.percent, xlab="PC", ylab="Percent Variance", names.arg=1:length(var.percent), las=1, ylim=c(0,max(var.percent)), col="gray", main  = paste("PCA Scree Plot for",cellLine,ss))

p <- ggplot(data.frame(l4PCAModel$x), aes(x = PC1, y = PC2, col = l4F$Barcode)) +
  geom_point(size = rel(.8), alpha = .8) +
  labs(colour = "Barcode")+
  theme(legend.text=element_text(size = 6))+
  guides(colour = guide_legend(override.aes = list(size=6)))
print(p)

p <- ggplot(data.frame(l4PCAModel$x), aes(x = PC1, y = PC3, col = l4F$Barcode)) +
  geom_point(size = rel(.8), alpha = .8) +
  labs(colour = "Barcode")+
  theme(legend.text=element_text(size = 6))+
  guides(colour = guide_legend(override.aes = list(size=6)))
print(p)

p <- ggplot(data.frame(l4PCAModel$x), aes(x = PC2, y = PC3, col = l4F$Barcode)) +
  geom_point(size = rel(.8), alpha = .8) +
  labs(colour = "Barcode")+
  theme(legend.text=element_text(size = 6))+
  guides(colour = guide_legend(override.aes = list(size=6)))
print(p)

p <- ggplot(data.frame(l4PCAModel$x), aes(x = PC1, y = PC2, col = l4F$Cytoplasm_CP_Intensity_IntegratedIntensity_MitoTracker_RZSNorm)) +
  geom_point(size = rel(.8), alpha = .8) +
  labs(colour = "Total MitoTracker")+
  theme(legend.text=element_text(size = 6))+
  guides(colour = guide_legend(override.aes = list(size=6)))
print(p)


p <- ggplot(data.frame(l4PCAModel$x), aes(x = PC1, y = PC2, col = l4F$Cells_CP_AreaShape_Perimeter_RZSNorm)) +
  geom_point(size = rel(.8), alpha = .8) +
  labs(colour = "Perimeter ProportionStatus")+
  theme(legend.text=element_text(size = 6))+
  guides(colour = guide_legend(override.aes = list(size=6)))
print(p)

kable(l4PCAModel$rotation[,1:nrPCs], digits = 2)

#Extract the selectedPCs of each MEP into a matrix
l4pcv <- l4PCAModel$x[,1:nrPCs]

```


\newpage

#Supplemental Material


##MEMA Layout
All MEMAs in the experiment are in separate wells and have the same design of 46 ECM proteins spotted in 35 rows and 20 columns. The proteins are randomly assigned to spots in the top 30 rows. Rows 31-35 are replicates of rows 1-5. The upper left and bottom right corners of each MEMA are image fiducials in the 488nm channel and there are four blank spots for checking orientation in all channels.   

```{r Content Layout,echo=FALSE, message=FALSE, warnings=FALSE, fig.width=6}

#Select the A row wells and delete the blanks
setkey(l1,Well)
DT <- unique(l1[grep("A",unique(l1$Well),value=TRUE),list(ArrayRow,ArrayColumn,ECMp)])
setkey(DT,ECMp)
DT <- DT[!"blank"]

p <- ggplot(DT,aes(x = ArrayColumn, y = ArrayRow, fill=ECMp))+
  geom_point(shape=21, size = 2.2)+
  guides(fill=guide_legend(ncol = 4))+
  theme(legend.text = element_text(size = rel(.5)),legend.title=element_text(size = rel(.5)),plot.title=element_text(size = rel(.8)))+
  scale_y_reverse()+
  xlab("")+ylab("")+
  ggtitle(" \n\nLINCS MEMA A Row Layout")
print(p)

```


##Replicate Count
The MEMAs have an average of 15 replicates with a range from 13 to 19.  

```{r Layout Replicate Count,echo=FALSE, message=FALSE, warnings=FALSE, fig.width=6.5, fig.height=3}
#Remove the fiducial and blank entries
setkey(DT,ECMp)
DT <- DT[!"fiducial"]
DT <- DT[!"blank"]


p <- ggplot(DT, aes(x=ECMp))+
  geom_bar(width=.8)+geom_hline(yintercept = mean(table(DT$ECMp)), colour="blue")+
  ggtitle(" \n\nCount of Replicate ECM Proteins In Each MEMA")+
  xlab("Printed ECM Protein")+ylab("Number of spots")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=rel(.8)),axis.title.x = element_text(size=rel(.8)),axis.title.y = element_text(size=rel(.8)),plot.title = element_text(size = rel(1)),legend.text=element_text(size = rel(.3)),legend.title=element_text(size = rel(.5)))

print(p)

```

##Quality Analysis
The variance of the signal in MEMA data comes from biological and technical factors. The technical factors create regions of low cell counts per spot and uneven staining across the array. The goal of the QA pipeline is to quantify the technical factors to identify wells or plates that need to be removed from downstream processing and/or be replaced by wells from a new experiment.

The hypothesis for the MEMA QA process is that the biological signal comes from individual spots while the technical variations come from regions of low signal. A bivariate loess model can be used to quantify the number of spots in low signal regions, leading to a MEMA QA score. 
\newpage

###Loess Model Explanation
The loess model of a MEMA is the mean value of a weighted version of each spot's region or neighborhood. In a 700 spot array, a loess span value of 0.1 sets the size of the neighborhood to be the nearest 70 points (within approximately 5 spots in all directions). The weights are a tricubic function of the euclidean distance between the spot being modeled and the neighborhood spots. These weights vary from 1 to 0 as distances increase from the nearest to the farthest neighbor. In other words, each spot in the model takes on the mean value of its 70 nearest neighbors with the closest neighbors having the largest impact. Therefore, the loess model is dominated by the technical regional factors as opposed to individual biological responses.

A MEMA's QA score is derived from the loess model of the control-well-normalized values by calculating the proportion of spots in low signal regions(LSR). A threshold for classifying spots as LSR is based on the median of each plate's control well. To have higher scores reflect increasing quality, the MEMA QA score is defined as the proportion of non-LSR spots to total spots. This value will be 1 for MEMAs with no low signal regions and approach 0 as the number of LSR spots increases.

Below are plots showing data from well B01 from plate LI8X00110 from LINCS staining set 2. The LSR spots are those to the left of the blue vertical line at the threshold value of `r lthresh ` in the histogram. 


```{r Loess_Model_explanation , echo=FALSE, fig.width=2.5,fig.height=4}
setkey(l3,Barcode,Well)
DT <-l3[.(barcodes[1],"A01")]
#Remove the fiducial entries
setkey(DT,ECMp)
DT <- DT[!"fiducial"]
DT <- DT[!"blank"]

p <- ggplot(DT, aes(x=ArrayColumn, y=ArrayRow,colour=Spot_PA_SpotCellCount))+
  geom_point(size=1.8)+
  scale_y_reverse()+   scale_x_continuous(breaks= c(min(DT$ArrayColumn),round(mean(c(min(DT$ArrayColumn),max(DT$ArrayColumn)))),max(DT$ArrayColumn)))+
  scale_colour_gradient(low = "white", high = "red")+
  guides(colour = guide_legend("Spot Cell\nCount", keywidth = .5, keyheight = .5))+
  ggtitle(paste("\n\n","Spot Cell Count for",unique(DT$CellLine), "cells \nin plate",unique(DT$Barcode)))+
  xlab("")+ylab("")+theme(axis.text.x = element_text(angle = 0, vjust = 0.5, hjust=1),plot.title = element_text(size = rel(.5)),legend.text=element_text(size = rel(.3)),legend.title=element_text(size = rel(.3)))
suppressWarnings(print(p))

p <- ggplot(DT, aes(x=ArrayColumn, y=ArrayRow,colour=Spot_PA_LoessSCC))+
  geom_point(size=1.8)+
  scale_y_reverse()+   scale_x_continuous(breaks= c(min(DT$ArrayColumn),round(mean(c(min(DT$ArrayColumn),max(DT$ArrayColumn)))),max(DT$ArrayColumn)))+
  scale_colour_gradient(low = "white", high = "red")+
  guides(colour = guide_legend("Normalized \nSpot Cell \nCount", keywidth = .5, keyheight = .5))+
  ggtitle(paste("\n\n","Loess Model of Spot Cell Count \nfor",unique(DT$CellLine), "cells in plate",unique(DT$Barcode)))+
  xlab("")+ylab("")+theme(axis.text.x = element_text(angle = 0, vjust = 0.5, hjust=1),plot.title = element_text(size = rel(.5)),legend.text=element_text(size = rel(.3)),legend.title=element_text(size = rel(.3)))
suppressWarnings(print(p))

  DT <- DT[,QAScore := calcQAScore(.SD,threshold=lthresh,value="Spot_PA_LoessSCC"),by="Well"]
  wellScores <- unique(DT[,list(Well,QAScore=sprintf("%.2f",QAScore))])

  p <- ggplot(DT, aes(x=Spot_PA_LoessSCC))+
  geom_histogram(binwidth=.02)+
  geom_vline(xintercept=lthresh, colour="blue")+
  geom_text(data=wellScores, aes(label=paste0("QA\n",QAScore)), x = .9, y = 30, size = rel(5), colour="red")+
  ggtitle(paste("\n\n","Loess Model of Spot Cell Count \nfor",unique(DT$CellLine), "cells in plate",unique(DT$Barcode)))+xlab("Spot Cell Count")+
  theme(axis.text.x = element_text(angle = 0, vjust = 0.5, hjust=1, size=rel(.5)), axis.text.y = element_text(angle = 0, vjust = 0.5, hjust=1, size=rel(1)), plot.title = element_text(size = rel(.5)),legend.text=element_text(size = rel(.3)),legend.title=element_text(size = rel(.3)))
suppressWarnings(print(p))


```

\newpage

##Replicates after filtering

```{r FilteredReplicateCount,echo=FALSE, message=FALSE, warnings=FALSE, fig.width=6.5, fig.height=5}

p <- ggplot(l4F[!grepl("FBS",l4F$Ligand)], aes(x = ECMp, y=Spot_PA_ReplicateCount))+
  geom_boxplot()+
  ggtitle(" \n\nCount of MEP Replicates by ECM Protein")+
  xlab("")+ylab("Replicate Count")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=rel(.8)),axis.title.x = element_text(size=rel(.8)),axis.title.y = element_text(size=rel(.8)),plot.title = element_text(size = rel(1)),legend.text=element_text(size = rel(.3)),legend.title=element_text(size = rel(.5)))

print(p)

p <- ggplot(l4F[!grepl("FBS",l4F$Ligand)], aes(x = Ligand, y=Spot_PA_ReplicateCount))+
  geom_boxplot()+
  ggtitle(" \n\nCount of MEP Replicates by Ligand")+
  xlab("")+ylab("Replicate Count")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=rel(.8)),axis.title.x = element_text(size=rel(.8)),axis.title.y = element_text(size=rel(.8)),plot.title = element_text(size = rel(1)),legend.text=element_text(size = rel(.3)),legend.title=element_text(size = rel(.5)))

print(p)

```

\newpage

##Stain Pseudoimages
The pseudoimages of each well's raw signals are shown in the plots below. Wells that could not be sucessfully imaged due to focus issues are missing from the pseudoimages.

```{r Pseudoimages_all_stains, echo=FALSE, fig.width=3.7,fig.height=4, eval=TRUE}


for (barcode in barcodes){
  DT <-l3[l3$Barcode==barcode,]
  #Remove the fiducial entries
  setkey(DT,ECMp)
  DT <- DT[!"fiducial"]
  DT <- DT[!"blank"]
  
  p <- create8WellPseudoImage(DT, pr = "Nuclei_CP_Intensity_MedianIntensity_Dapi", prDisplay = "Median DAPI")
  suppressWarnings(print(p))
  
  p <- create8WellPseudoImage(DT, pr = "Cytoplasm_CP_Intensity_MedianIntensity_Actin", prDisplay = unique(DT$Endpoint488))
  suppressWarnings(print(p))
  
    p <- create8WellPseudoImage(DT, pr = "Cytoplasm_CP_Intensity_MedianIntensity_CellMask",prDisplay = unique(DT$Endpoint555))
  suppressWarnings(print(p))
  
  p <- create8WellPseudoImage(DT, pr = "Cytoplasm_CP_Intensity_MedianIntensity_MitoTracker",prDisplay = unique(DT$Endpoint647))
  suppressWarnings(print(p))
  
}

```

\newpage

```{r SCCFBS_response, eval = TRUE}

p <- ggplot(l3[l3$Ligand == "FBS"], aes(x = reorder(Ligand, Spot_PA_SpotCellCount, FUN=median), y = Spot_PA_SpotCellCount, colour = Barcode))+geom_boxplot()+
  ggtitle(paste("\n\nSpot Cell Count in High Serum Wells"))+
  xlab("")+ylab("")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size=rel(.8)), axis.text.y = element_text(angle = 0, vjust = 0.5, hjust=1, size=rel(1)), plot.title = element_text(size = rel(1)),legend.text=element_text(size = rel(.3)),legend.title=element_text(size = rel(.3)))
suppressWarnings(print(p))

```

\newpage

###Extreme Spot Cell Count MEPs

```{r SCCMeps,echo=FALSE, fig.width=8, fig.height=5, eval=TRUE}

setkey(l4F, Spot_PA_SpotCellCount_RZSNorm_RobustZ)
kable(l4F[l4F$Spot_PA_SpotCellCount_RZSNorm_RobustZ >=2| 
            l4F$Spot_PA_SpotCellCount_RZSNorm_RobustZ <=-2, list(Ligand,ECMp,Spot_PA_SpotCellCount_RZSNorm_RobustZ, Spot_PA_SpotCellCount)],digits = 2)
```

\newpage

###Active MEPS from Unbiased Analysis
The table below shows the first and last active MEPs and their PC values.

```{r ActiveMeps, echo = FALSE}
kable(head(activeMEPs, n=20), digits = 2)
kable(tail(activeMEPs, n=20), digits = 2)

```

\newpage

##Cell Cycle Plots
Cell Cycle plots include univariate plots of the total DAPI signal.

```{r Cell_cycle_plots,echo=FALSE, fig.width=8, fig.height=5, eval=TRUE}

plotTotalDAPI(l1, barcodes)

if (ss == "SS1"){
  
  
} else if (ss == "SS2"){
  
} 
```