HMEC_ECMs.Rmd

---
title: HMEC ECM Protein Analysis 
author: Mark Dane

output: 
  flexdashboard::flex_dashboard:
    storyboard: true
    source_code: embed
---


```{r "setup", include=FALSE}
knitr::opts_chunk$set(fig.width=12, fig.height=8, echo=FALSE, warning=FALSE, message=FALSE, cache=FALSE)
knitr::opts_knit$set(root.dir = "/data/share/dane/HMEC_ECMs")

```


```{r LoadAndCleanData, cache=FALSE}
#Author: Mark Dane, copyright 2017
library(ggplot2)
library(scales)
library(data.table)
library(MEMA)
library(grid)
library(knitr)
library(gplots)
library(RColorBrewer)
library(DT)
library(d3heatmap)
library(plotly)
library(corrplot)
library(Rtsne)
library(synapseClient)
library(NMF)
library(stringr)

#Setup colors for Barcode and text in all heatmaps
selDark2 <- colorRampPalette(brewer.pal(8,"Dark2"))
plateCol <- selDark2(8)
hmcols<-colorRampPalette(c("blue","white","red"))(16)

normToCOL1Median <- function(dt){
  m <- numericMedian(unlist(dt[dt$ECMp=="COL1",1]))
  if(m==0)stop(paste("Can't normalize when median equals 0",unique(dt$Barcode),colnames(dt)))
  r <- dt[,1]/m
}

#Load the latest level 3 SS4 HMEC datasets from Synapse
synapseLogin()
HMEC122L_l3 <- "syn8565101"
HMEC240L_l3 <- "syn8606270"
#Get and clean the datasets
l3L <- lapply(c(HMEC240L_l3,HMEC122L_l3), function(x){
  synGet(x) %>%
    getFileLocation() %>%
    fread(.,showProgress=FALSE) %>%
    filter(grepl("FBS",Ligand)) #Filter to only the FBS wells
})

#Combine into one data.table
l3 <- rbindlist(l3L)

#Normalize all log and logit signalsto their Collagen 1values without doing their SE values
sigNames <- grep("_SE",grep("Norm$", colnames(l3),value=TRUE),value=TRUE,invert=TRUE)
for(sigName in sigNames){
l3 <- l3[,sigNameNorm := normToCOL1Median(.SD), .SDcols=c(sigName,"ECMp","Barcode"),by=c("Barcode")]
setnames(l3,"sigNameNorm",paste0(sigName,"COL1MedianNorm"))
}

#Calculate wilcoxon rank tests within a cell line
wilcoxpValuesWithinACellLine <- function(ECMp,signalValues) {
  #Get the COL1 values for the signal
  COL1Values <- signalValues[ECMp=="COL1"]
  res <- lapply(unique(ECMp[!ECMp=="COL1"]), function(e){
    p <- wilcox.test(x = signalValues[ECMp==e], y=COL1Values)$p.value
    names(p) <- e
    p
  })
  data.table(ECMp=names(unlist(res)), wilcoxPValue=unlist(res))
}

#Summarize the signals at each spot within a cell line
normMedians <- unique(l3[,.(Spot_PA_SpotCellCountNormCOL1MedianNorm=numericMedian(Spot_PA_SpotCellCountNormCOL1MedianNorm),
                               Nuclei_PA_Gated_EdUPositiveProportionNormCOL1MedianNorm=numericMedian(Nuclei_PA_Gated_EdUPositiveProportionNormCOL1MedianNorm),
                               Cytoplasm_PA_Gated_KRT19PositiveProportionNormCOL1MedianNorm = numericMedian(Cytoplasm_PA_Gated_KRT19PositiveProportionNormCOL1MedianNorm),
                             ECMp),by=c("CellLine","Spot")]) 

#Calculate wilcoxon rank sum p values within each cell line comparing spot cell counts for each ECMp to COL1 
dtp122L <- normMedians[grepl("HMEC122L",CellLine), wilcoxpValuesWithinACellLine(ECMp, Spot_PA_SpotCellCountNormCOL1MedianNorm)]
dtp122L$CellLine <- "HMEC122L"
dtp122L <- dtp122L[,pAdjustedWithinCellLineSCC := p.adjust(wilcoxPValue,method="fdr")]

dtp240L <- normMedians[grepl("HMEC240L",CellLine), wilcoxpValuesWithinACellLine(ECMp, Spot_PA_SpotCellCountNormCOL1MedianNorm)]
dtp240L$CellLine <- "HMEC240L"
dtp240L <- dtp240L[,pAdjustedWithinCellLineSCC := p.adjust(wilcoxPValue,method="fdr")]

#Add the p values that compare each SpotCellCount to the COL1 values
l3 <- merge(l3,rbind(dtp122L,dtp240L),by=c("ECMp","CellLine"), all=TRUE)

#Calculate wilcoxon rank sum p values within each cell line comparing EdU+ Proportion for each ECMp to COL1 
dtp122L <- normMedians[grepl("HMEC122L",CellLine), wilcoxpValuesWithinACellLine(ECMp, Nuclei_PA_Gated_EdUPositiveProportionNormCOL1MedianNorm)]
dtp122L$CellLine <- "HMEC122L"
dtp122L <- dtp122L[,pAdjustedWithinCellLineEdUPP := p.adjust(wilcoxPValue,method="fdr")]

dtp240L <- normMedians[grepl("HMEC240L",CellLine), wilcoxpValuesWithinACellLine(ECMp, Nuclei_PA_Gated_EdUPositiveProportionNormCOL1MedianNorm)]
dtp240L$CellLine <- "HMEC240L"
dtp240L <- dtp240L[,pAdjustedWithinCellLineEdUPP := p.adjust(wilcoxPValue,method="fdr")]

#Add the p values that compare each EdU positive proportion to the COL1 values
l3 <- merge(l3,rbind(dtp122L,dtp240L),by=c("ECMp","CellLine"), all=TRUE)

#Calculate wilcoxon rank sum p values within each cell line comparing KRT19 High Proportion for each ECMp to COL1 
dtp122L <- normMedians[grepl("HMEC122L",CellLine), wilcoxpValuesWithinACellLine(ECMp, Cytoplasm_PA_Gated_KRT19PositiveProportionNormCOL1MedianNorm)]
dtp122L$CellLine <- "HMEC122L"
dtp122L <- dtp122L[,pAdjustedWithinCellLineKRT19PP := p.adjust(wilcoxPValue,method="fdr")]

dtp240L <- normMedians[grepl("HMEC240L",CellLine), wilcoxpValuesWithinACellLine(ECMp, Cytoplasm_PA_Gated_KRT19PositiveProportionNormCOL1MedianNorm)]
dtp240L$CellLine <- "HMEC240L"
dtp240L <- dtp240L[,pAdjustedWithinCellLineKRT19PP := p.adjust(wilcoxPValue,method="fdr")]

#Add the p values that compare each KRT19 High Proportion to the COL1 values
l3 <- merge(l3,rbind(dtp122L,dtp240L),by=c("ECMp","CellLine"), all=TRUE)

#Calculate wilcoxon rank tests across cell lines
wilcoxpValuesAcrossCellLines <- function(DT){
  sigName <- colnames(DT)[2]
  tmp <- wilcox.test(DT[[sigName]][DT$CellLine==unique(DT$CellLine)[1]],
                DT[[sigName]][DT$CellLine==unique(DT$CellLine)[2]])$p.value
}

dtp <- normMedians[,.(wilcoxPValueAcrossCellLines=wilcoxpValuesAcrossCellLines(.SD)),by="ECMp",.SDcols=c("CellLine","Spot_PA_SpotCellCountNormCOL1MedianNorm")]
dtp <- dtp[,pAdjustedAcrosCellLinesSCC := p.adjust(wilcoxPValueAcrossCellLines,method="fdr")]
#Add the p values that compare each SpotCellCount across cell lines
l3 <- merge(l3,dtp,by="ECMp")

dtp <- normMedians[,.(wilcoxPValueAcrossCellLines=wilcoxpValuesAcrossCellLines(.SD)),by="ECMp",.SDcols=c("CellLine","Nuclei_PA_Gated_EdUPositiveProportionNormCOL1MedianNorm")]
dtp <- dtp[,pAdjustedAcrosCellLinesEdUPP := p.adjust(wilcoxPValueAcrossCellLines,method="fdr")]
#Add the p values that compare each SpotCellCount across cell lines
l3 <- merge(l3,dtp,by="ECMp")

dtp <- normMedians[,.(wilcoxPValueAcrossCellLines=wilcoxpValuesAcrossCellLines(.SD)),by="ECMp",.SDcols=c("CellLine","Cytoplasm_PA_Gated_KRT19PositiveProportionNormCOL1MedianNorm")]
dtp <- dtp[,pAdjustedAcrosCellLinesKRT19PP := p.adjust(wilcoxPValueAcrossCellLines,method="fdr")]
#Add the p values that compare each SpotCellCount across cell lines
l3 <- merge(l3,dtp,by="ECMp")

```


###Experiment Overview


```{r Content Layout, fig.height=5, fig.width=5.5}

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

p <- ggplot(data=dt[!grepl("COL1",dt$ECMp),],aes(x = ArrayColumn, y = ArrayRow, fill=ECMp))+
  geom_point(shape=21, size = 2.7)+
  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("MEMA ECM Protein Layout")+
  theme(legend.text=element_text(size = rel(.5)),legend.title=element_text(size = rel(.8)))
p <- p+geom_point(data = dt[grepl("COL1",dt$ECMp),], shape=24, fill="black")
print(p)

#Well level COL1 Spot Cell Count Boxplots
dt <- l3[grepl("COL1",ECMp),]
#Melt the spot cell count value to show effects of normalization
dtm <- suppressWarnings(melt(dt, id.vars = c("Barcode"),measure.vars = c("Spot_PA_SpotCellCount","Spot_PA_SpotCellCountNorm","Spot_PA_SpotCellCountNormCOL1MedianNorm"), variable.name = "ProcessLevel",value.name = "Spot_Cell_Count", factorsAsstrings=FALSE))
levels(dtm$ProcessLevel) <- c("Raw","RUVLoessResidual","RUVLoessResidualandCol1")

p <- ggplot(dtm, aes(x=Barcode, y=Spot_Cell_Count, colour=ProcessLevel))+
  geom_boxplot(outlier.size = .5)+
  scale_y_log10()+
  xlab("Well Plate Barcode")+ylab("Spot Cell Count")+
  ggtitle("Normalization Results of Spot Cell Count")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0, hjust=0.5, size=rel(1)), 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(.8)),legend.title=element_text(size = rel(.8)))

print(p)
```

***
This project is a data mining exploration of the signals modulated by the ECM (Extracellular Matrix) proteins in the two MEMA (Microenvironment Microarray) HMEC cell line datasets from Oregon Health and Science University's MEP-LINCS project https://www.synapse.org/#!Synapse:syn2862345/wiki/394513.

There are eight replicate MEMAs in each cellline dataset, each in a different plate. 
All eight replicates of a cell line were processed in the wet lab as a batch.

This analysis is based on immunofluorescent cell-level data that has been summarized to the spot and ECM protein levels. 

#####MEMA Layout
All MEMAs in the experiment are in separate wells and have the same design of 48 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. There is a higher number of COL1 spots (shown as black triangles) throughout the array that are control spots. 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.

#####Boxplots
This dataset has raw, RUVLoessResidual normalized and RUVLoessResidual normalized that has then been normalized to the collagen 1 spots in the same MEMA. These figures show results of these three states of processing. The RUVLoessResiduals normalization removes the variance between the replicates that are in the same location in the MEMA.  

###Spot Cell Counts by ECM Protein and CellLine

```{r, fig.height=5, fig.width=14}

#ECMp Spot cell counts Boxplots
p <- ggplot(l3, aes(x=factor(ECMp), y=Spot_PA_SpotCellCountNormCOL1MedianNorm, colour=ECMp=="COL1", fill=pAdjustedWithinCellLineSCC<.2))+
  geom_boxplot(outlier.size = .5)+
    scale_fill_manual(values=c("transparent","cornflowerblue"))+
  scale_colour_manual(values=c("black","blue"))+
  guides(colour=FALSE)+
  labs(x="ECM Protein",
       y="Normalized Spot Cell Count",
       fill="Wilcoxon Rank Sum\n Adjusted pValue < .2",
       title="Normalized Spot Cell Count by ECM Protein")+
  facet_wrap(~CellLine)+
  theme(axis.text.x = element_text(angle = 90, vjust = 0, hjust=0.5, size=rel(1)), 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(.8)),legend.title=element_text(size = rel(.8)))
p <- p + geom_hline(yintercept = 1, colour="blue", alpha=.5)
p

#ECMp Spot cell counts Boxplots
p <- ggplot(l3, aes(x=factor(ECMp), y=Spot_PA_SpotCellCountNormCOL1MedianNorm, colour=CellLine, fill=pAdjustedWithinCellLineSCC<.2))+
  geom_boxplot(outlier.size = .5)+
    scale_fill_manual(values=c("transparent","cornflowerblue"))+
  scale_colour_manual(values=c("black","coral3"))+
  labs(x="ECM Protein",
       y="Normalized Spot Cell Count",
       fill="Wilcoxon Rank Sum\n Adjusted pValue < .2",
       title="Normalized Spot Cell Count by ECM Protein")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0, hjust=0.5, size=rel(1)), 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(.8)),legend.title=element_text(size = rel(.8)))
p <- p + geom_hline(yintercept = 1, colour="blue", alpha=.5)
p

p <- ggplot(l3, aes(x=factor(ECMp), y=Spot_PA_SpotCellCountNormCOL1MedianNorm, fill=pAdjustedAcrosCellLinesSCC<.2, colour=CellLine))+
  geom_boxplot(outlier.size = .5)+
    scale_fill_manual(values=c("transparent","cornflowerblue"))+
  scale_colour_manual(values=c("black","coral3"))+
    labs(x="ECM Protein",
       y="Normalized Spot Cell Count",
       colour="Cell Line",
       fill="Wilcoxon Rank Sum\n Adjusted pValue < .2",
       title="Normalized Spot Cell Count by ECM Protein")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0, hjust=0.5, size=rel(1)), 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(.8)),legend.title=element_text(size = rel(.8)))
p <- p + geom_hline(yintercept = 1, colour="blue", alpha=.5)
p

```

***
The first figure compares Spot Cell Count values from each ECM protein to the collagen 1 spots of the same cell line. Each boxplot shows the distribution of ~14 replicate spots from 8 different plates. Significant differences between the response are based on Wilcoxon Rank Sum tests between each ECM protein and collagen 1. Boxplots with fdr adjusted q-values less than 0.2 are filled blue.  


###Proliferation Signal by ECM Protein and CellLine


```{r, fig.height=5, fig.width=14}

#ECMp EdU Boxplots
p <- ggplot(l3, aes(x=factor(ECMp), y=Nuclei_PA_Gated_EdUPositiveProportionNormCOL1MedianNorm, colour=ECMp=="COL1", fill=pAdjustedWithinCellLineEdUPP<.2))+
  geom_boxplot(outlier.size = .5)+
    scale_fill_manual(values=c("transparent","cornflowerblue"))+
  scale_colour_manual(values=c("black","blue"))+
  guides(colour=FALSE)+
  labs(x="ECM Protein",
       y="Normalized EdU",
       fill="Wilcoxon Rank Sum\n Adjusted pValue < .2",
       title="Normalized EdU by ECM Protein")+
  facet_wrap(~CellLine)+
  theme(axis.text.x = element_text(angle = 90, vjust = 0, hjust=0.5, size=rel(1)), 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(.8)),legend.title=element_text(size = rel(.8)))
p <- p + geom_hline(yintercept = 1, colour="blue", alpha=.5)
p

#ECMp EdU Boxplots
p <- ggplot(l3, aes(x=factor(ECMp), y=Nuclei_PA_Gated_EdUPositiveProportionNormCOL1MedianNorm, colour=CellLine, fill=pAdjustedWithinCellLineEdUPP<.2))+
  geom_boxplot(outlier.size = .5)+
    scale_fill_manual(values=c("transparent","cornflowerblue"))+
  scale_colour_manual(values=c("black","coral3"))+
  labs(x="ECM Protein",
       y="Normalized EdU",
       fill="Wilcoxon Rank Sum\n Adjusted pValue < .2",
       title="Normalized EdU Positive Proportion by ECM Protein")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0, hjust=0.5, size=rel(1)), 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(.8)),legend.title=element_text(size = rel(.8)))
p <- p + geom_hline(yintercept = 1, colour="blue", alpha=.5)
p

p <- ggplot(l3, aes(x=factor(ECMp), y=Nuclei_PA_Gated_EdUPositiveProportionNormCOL1MedianNorm, fill=pAdjustedAcrosCellLinesEdUPP<.2, colour=CellLine))+
  geom_boxplot(outlier.size = .5)+
    scale_fill_manual(values=c("transparent","cornflowerblue"))+
  scale_colour_manual(values=c("black","coral3"))+
    labs(x="ECM Protein",
       y="Normalized EdU Positive Proportion",
       colour="Cell Line",
       fill="Wilcoxon Rank Sum\n Adjusted pValue < .2",
       title="Normalized EdU Positive Proportion by ECM Protein")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0, hjust=0.5, size=rel(1)), 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(.8)),legend.title=element_text(size = rel(.8)))
p <- p + geom_hline(yintercept = 1, colour="blue", alpha=.5)
p
```

***
The first figure compares EdU Positive Proportion values from each ECM protein to the collagen 1 spots of the same cell line. Each boxplot shows the distribution of ~14 replicate spots from 8 different plates. Significant differences between the response are based of Wilcoxon Rank Sum tests between each ECM protein and collagen 1. Boxplots with fdr adjusted p-values less than .2 are blue.  


###Lineage Signal by ECM Protein and CellLine



```{r, fig.height=5, fig.width=14}

#ECMp Lineage Boxplots
p <- ggplot(l3, aes(x=factor(ECMp), y=Cytoplasm_PA_Gated_KRT19PositiveProportionNormCOL1MedianNorm, colour=ECMp=="COL1", fill=pAdjustedWithinCellLineKRT19PP<.2))+
  geom_boxplot(outlier.size = .5)+
    scale_fill_manual(values=c("transparent","cornflowerblue"))+
  scale_colour_manual(values=c("black","blue"))+
  guides(colour=FALSE)+
  labs(x="ECM Protein",
       y="Normalized KRT19 Positive Proportion",
       fill="Wilcoxon Rank Sum\n Adjusted pValue < .2",
       title="Normalized KRT19 Positive Proportion by ECM Protein")+
  facet_wrap(~CellLine)+
  theme(axis.text.x = element_text(angle = 90, vjust = 0, hjust=0.5, size=rel(1)), 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(.8)),legend.title=element_text(size = rel(.8)))
p <- p + geom_hline(yintercept = 1, colour="blue", alpha=.5)
p

#ECMp lineage Boxplots
p <- ggplot(l3, aes(x=factor(ECMp), y=Cytoplasm_PA_Gated_KRT19PositiveProportionNormCOL1MedianNorm, colour=CellLine, fill=pAdjustedWithinCellLineKRT19PP<.2))+
  geom_boxplot(outlier.size = .5)+
    scale_fill_manual(values=c("transparent","cornflowerblue"))+
  scale_colour_manual(values=c("black","coral3"))+
  labs(x="ECM Protein",
       y="Normalized KRT19 Positive Proportion",
       fill="Wilcoxon Rank Sum\n Adjusted pValue < .2",
       title="Normalized KRT19 Positive Proportion by ECM Protein")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0, hjust=0.5, size=rel(1)), 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(.8)),legend.title=element_text(size = rel(.8)))
p <- p + geom_hline(yintercept = 1, colour="blue", alpha=.5)
p

p <- ggplot(l3, aes(x=factor(ECMp), y=Cytoplasm_PA_Gated_KRT19PositiveProportionNormCOL1MedianNorm, fill=pAdjustedAcrosCellLinesKRT19PP<.2, colour=CellLine))+
  geom_boxplot(outlier.size = .5)+
    scale_fill_manual(values=c("transparent","cornflowerblue"))+
  scale_colour_manual(values=c("black","coral3"))+
    labs(x="ECM Protein",
       y="Normalized KRT19 Positive Proportion",
       colour="Cell Line",
       fill="Wilcoxon Rank Sum\n Adjusted pValue < .2",
       title="Normalized KRT19 Positive Proportion by ECM Protein")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0, hjust=0.5, size=rel(1)), 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(.8)),legend.title=element_text(size = rel(.8)))
p <- p + geom_hline(yintercept = 1, colour="blue", alpha=.5)
p

```

***
The first figure compares KRT19 Positive Proportion values from each ECM protein to the collagen 1 spots of the same cell line. Each boxplot shows the distribution of ~14 replicate spots from 8 different plates. Significant differences between the response are based of Wilcoxon Rank Sum tests between each ECM protein and collagen 1. Boxplots with fdr adjusted p-values less than .2 are blue.  


###Scatterplots of the RUVLoessResiduals normalized signals


```{r ScatterplotsOfReplicatesRUV, fig.height=4, fig.width=7}

compScatterplots <- function(dt,sig,title,limits,sdInt=.25) {
  dtc <- dcast(dt,ECMp~CellLine,fun.aggregate = median, value.var = sig, fill=NA)
  p <- ggplot(dtc, aes(x=HMEC122L, y=HMEC240L, colour=ECMp))+
    geom_point(size = rel(1), alpha = .6) +
    guides(colour=FALSE)+
    labs(colour = "ECM Protein")+
    ggtitle(title)+
    coord_cartesian(xlim=limits, ylim=limits)+
    theme(legend.text=element_text(size = 5),
          plot.title=element_text(size = rel(.7)))
  p <- p + geom_abline(slope = 1,intercept = 0)
  p <- p + geom_abline(slope = 1,intercept = sdInt, linetype=2)
  p <- p + geom_abline(slope = 1,intercept = -sdInt, linetype=2)
  print(p)
}

dt <- l3

#Select the COL1+RUV normalized values
COL1NormFVNames <- grep("NormCOL1MedianNorm",colnames(l3), value=TRUE)
#Median summarize the signals to the MEP/ECMp level within each cell line
mc <- l3[!grepl("NID1|ELN",ECMp),lapply(.SD,numericMedian), by=c("ECMp","CellLine"),.SDcols=COL1NormFVNames]
#Calculate the difference between the cell lines
MEPs <- mc[,lapply(.SD,diff), by=c("ECMp"),.SDcols=COL1NormFVNames]

compScatterplots(dt, sig = "Spot_PA_SpotCellCountNormCOL1MedianNorm",title = "Comparison of RUV+COL1 Normalized Spot Cell Counts\n in the HMEC Cell Lines",limits=c(0,2), sdInt = sd(MEPs$Spot_PA_SpotCellCountNormCOL1MedianNorm))
compScatterplots(dt,sig = "Nuclei_PA_Gated_EdUPositiveProportionNormCOL1MedianNorm",title = "Comparison of RUV+COL1 Normalized Proliferation\n in the HMEC Cell Lines",limits=c(0,5), sdInt = sd(MEPs$Nuclei_PA_Gated_EdUPositiveProportionNormCOL1MedianNorm))
compScatterplots(dt, sig = "Cytoplasm_PA_Gated_KRT19PositiveProportionNormCOL1MedianNorm",title = "Comparison of RUV+COL1 Normalized Luminal Signal\n  in the HMEC Cell Lines",limits=c(0,2.5), sdInt = sd(MEPs$Cytoplasm_PA_Gated_KRT19PositiveProportionNormCOL1MedianNorm))

```

***
These scatterplots are of the RUV+COL1 normalized cell count, proliferation and lineage signals. The diagonal lines have slope=1 and are for reference.

###PCA of all normalized signals

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

createPCAModel <- function(dt,fvNames,verbose=FALSE){
  #delete any non-varying signals that cannot be scaled to unit variance
  stdDevs <- apply(dt[,fvNames, with=FALSE], 2, sd)
  if(any(stdDevs==0)&verbose)  cat(paste("Deleting",fvNames[stdDevs==0], "due to 0 variance\n"))
  fvNames <- fvNames[!stdDevs==0]
  PCAModel <- prcomp(as.formula(paste(" ~ ", paste(fvNames, collapse = "+"))), data = dt, scale.= TRUE)
  rownames(PCAModel$x)<-paste(dt$MEP,gsub("LI8X00","",dt$Barcode), sep="_")
  return(PCAModel)
}

plotScree <- function(m,nrPCs=9,main  = paste("PCA Scree Plot")){
  var <- m$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=main)
}


#Get a PCA model of the RUV+COL1 normalized data
COL1NormPCAModel <- createPCAModel(dt=MEPs, fvNames = COL1NormFVNames)
dt <- cbind(data.table(COL1NormPCAModel$x), ECMp=MEPs$ECMp)

plotScree(COL1NormPCAModel)

p12 <- ggplot(dt, aes(x = PC1, y = PC2, colour = ECMp)) +
      geom_point(size = rel(2), alpha = .8) +
      ggtitle("PCA of The Cell Line Differences of All RUV+COL1 Normalized Signals")+
      #coord_cartesian(xlim=xlim, ylim=ylim)+
  guides(colour=FALSE)+
      theme(legend.text=element_text(size = 6),
            plot.title=element_text(size = rel(.7)))
p12

p13 <- ggplot(dt, aes(x = PC1, y = PC3, colour = ECMp)) +
      geom_point(size = rel(2), alpha = .8) +
      ggtitle("PCA of The Cell Line Differences of All RUV+COL1 Normalized Signals")+
      #coord_cartesian(xlim=xlim, ylim=ylim)+
  guides(colour=FALSE)+
      theme(legend.text=element_text(size = 6),
            plot.title=element_text(size = rel(.7)))
p13

p23 <- ggplot(dt, aes(x = PC2, y = PC3, colour = ECMp)) +
      geom_point(size = rel(2), alpha = .8) +
      ggtitle("PCA of The Cell Line Differences of All RUV+COL1 Normalized Signals")+
      #coord_cartesian(xlim=xlim, ylim=ylim)+
  guides(colour=FALSE)+
      theme(legend.text=element_text(size = 6),
            plot.title=element_text(size = rel(.7)))
p23

p14 <- ggplot(dt, aes(x = PC1, y = PC4, colour = ECMp)) +
      geom_point(size = rel(2), alpha = .8) +
    guides(colour=FALSE)+
      ggtitle("PCA of The Cell Line Differences of All RUV+COL1 Normalized Signals")+
      #coord_cartesian(xlim=xlim, ylim=ylim)+
      theme(legend.text=element_text(size = 6),
            plot.title=element_text(size = rel(.7)))
p14

```

###Interactive PCA of all RUV+Col1 Normalized Signals

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

p12 <- ggplot(dt, aes(x = PC1, y = PC2, colour = ECMp)) +
      geom_point(size = rel(2), alpha = .8) +
      labs(colour = "ECM Protein")+
      ggtitle("PCA of The Cell Line Differences of All RUV+COL1 Normalized Signals")+
      #coord_cartesian(xlim=xlim, ylim=ylim)+
      theme(legend.text=element_text(size = 6),
            plot.title=element_text(size = rel(.7)))
ggplotly(p12)

```

***

####Method
The responses in all normalized signals are first median summerized across the replicates.  
Next, the differences between the cell lines are calculated.  
PCA is run on the differences to reduce from `r length(COL1NormFVNames)` signals.  
The scatterplot shows the first 2 principal components of all normalized signals, colored by ECM protein. Hovering over a spot shows its ECM protein. 


###tSNE analysis on all normalized features

```{r createtSNE}
dt <- MEPs
tSNEParameters <- grep("Norm",colnames(dt), value=TRUE)
stdDevs <- apply(dt[,tSNEParameters, with=FALSE], 2, sd)
tSNEParameters <- tSNEParameters[!stdDevs==0]
data_matrix <- scale(as.matrix(dt[,tSNEParameters, with=FALSE]))

set.seed(1234)
data_tsne <- Rtsne(data_matrix, perplexity = 15)
dtt <- cbind(dt,data_tsne$Y)
```


```{r plottSNE}

p <- ggplot(dtt, aes_string(x="V1", y="V2", text="ECMp", colour = "ECMp"))+
  geom_point(size=rel(2), alpha=.8)+
  ggtitle(paste("tSNE plot of All N"))+
  #scale_colour_gradient(low="black",high="red",oob = scales::squish)+
  guides(colour =FALSE)
ggplotly(p)

```

***

tSNE dimension reduction on collagen 1 MEP's paramaters, colored by ligand.


###Interactive Scatterplot of RUV+Col1 Normalized Spot Cell Count

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


compScatterplotly <- function(dt,sig,title,limits,sdInt=.25) {
  dtc <- dcast(dt,ECMp~CellLine,fun.aggregate = median, value.var = sig, fill=NA)
  p <- ggplot(dtc, aes(x=HMEC122L, y=HMEC240L, colour=ECMp))+
    geom_point(size = rel(2), alpha = .6) +
    labs(colour = "ECM Protein")+
    ggtitle(title)+
    coord_cartesian(xlim=limits, ylim=limits)+
    theme(legend.text=element_text(size = 5),
          plot.title=element_text(size = rel(.7)))
  p <- p + geom_abline(slope = 1,intercept = 0)
  p <- p + geom_abline(slope = 1,intercept = sdInt, linetype=2)
  p <- p + geom_abline(slope = 1,intercept = -sdInt, linetype=2)
  ggplotly(p)
}

dt <- l3[!grepl("NID1|ELN",ECMp),]

#Select the COL1+RUV normalized values
COL1NormFVNames <- grep("NormCOL1MedianNorm",colnames(l3), value=TRUE)
#Median summarize the signals to the MEP/ECMp level within each cell line
mc <- l3[!grepl("NID1|ELN",ECMp),lapply(.SD,numericMedian), by=c("ECMp","CellLine"),.SDcols=COL1NormFVNames]
#Calculate the difference between the cell lines
MEPs <- mc[,lapply(.SD,diff), by=c("ECMp"),.SDcols=COL1NormFVNames]

compScatterplotly(dt, sig = "Spot_PA_SpotCellCountNormCOL1MedianNorm",title = "Comparison of RUV+COL1 Normalized Spot Cell Counts\n in the HMEC Cell Lines",limits=c(0,2), sdInt = sd(MEPs$Spot_PA_SpotCellCountNormCOL1MedianNorm))


```

***

###Interactive Scatterplot of RUV+Col1 Normalized Proliferation

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

compScatterplotly(dt, sig = "Nuclei_PA_Gated_EdUPositiveProportionNormCOL1MedianNorm",title = "Comparison of RUV+COL1 Normalized Proliferation\n in the HMEC Cell Lines",limits=c(0,5), sdInt = sd(MEPs$Nuclei_PA_Gated_EdUPositiveProportionNormCOL1MedianNorm))

```

***


###Interactive Scatterplot of RUV+Col1 Normalized Lineage

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

compScatterplotly(dt, sig = "Cytoplasm_PA_Gated_KRT19PositiveProportionNormCOL1MedianNorm",title = "Comparison of RUV+COL1 Normalized Lineage\n in the HMEC Cell Lines",limits=c(0,2.5), sdInt = sd(MEPs$Cytoplasm_PA_Gated_KRT19PositiveProportionNormCOL1MedianNorm))

```

***

###QA Figures
```{r}
#Show the row scc values and evaluate for filtering
dt <- l3[,ArrayRowSCC := numericMedian(Spot_PA_SpotCellCount), by=c("ArrayRow","Barcode")]
dt <- l3[,ArrayRowSCCNorm := numericMedian(Spot_PA_SpotCellCountNorm), by=c("ArrayRow","Barcode")]

p <- ggplot(dt, aes(x=ArrayRow, y=ArrayRowSCC, colour=CellLine))+
  geom_point()
p
p <- ggplot(dt, aes(x=ArrayRow, y=ArrayRowSCCNorm, colour=CellLine))+
  geom_point()
p
```

***

###Normalization Checks

```{r}
dt <- l3[l3$CellLine=="HMEC240L",]
dt <- dt[dt$ECMp=="COL1",]
dt <- dt[dt$Well=="A03",]
dt <- dt[,SCCSpotMedian := numericMedian(Spot_PA_SpotCellCount),by=c("Spot")]
dt <- dt[,SCCECMpMedian := numericMedian(Spot_PA_SpotCellCount),by=c("Barcode","ECMp")]
dtc <- dcast(dt, Spot~Barcode, value.var = "Spot_PA_SpotCellCountNorm")

p <- ggplot(dtc, aes(x=LI8X00651, y=LI8X00652))+
  geom_point()+
  labs(x="Normalized Cell Count for LI8X00651",
       y="Normalized Cell Count for LI8X00652",
       title="Normalized values of COL1 spots at the same position in different wells")
p

p <- ggplot(dt, aes(x=Spot_PA_SpotCellCountNorm, y=SCCECMpMedian))+
  geom_point(alpha=.3)+
  facet_wrap(~Barcode)+
  labs(x="Normalized Spot Cell Count",
       y="Median COL1 SCC for the well",
       title="Median COL1_FBS SCC vs Normalized Value")
p

```

***

These plots provide evidence of the normalization results. Spots with the same conditions and locations in different plates are normalized until they have similar/identical values. Spots within a plate that have the same conditions but different locations have different values.