---
title: "Analysis of vgll3 expression (PLOS Genetics)"
author: "J-P Verta"
date: "8/12/2020"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(reshape2)
library(ggplot2)
library(data.table)
library(plyr)
library(plotly)
library(dplyr)
library(broom)
library(tidyr)
library(ggExtra)
library(car)
library(ggpubr)
library(nlme)
library(ggbeeswarm)
library(mclust)
library(factoextra)
library(emmeans)
library(ggsignif)
library(cowplot)
```

## Comparison of vgll3 expression acros tissues 

Load locally saved data files.
```{r}
load("ddPCRdataFiles.RData")
```


Assign shorted tissue labels.
```{r}
design$TissueOrRep[which(design$TissueOrRep == "heart")] = "Hr"
design$TissueOrRep[which(design$TissueOrRep == "liver")] = "Liv"
design$TissueOrRep[which(design$TissueOrRep == "muscle")] = "Mc"
design$TissueOrRep[which(design$TissueOrRep == "adipose")] = "Ad"
design$TissueOrRep[which(design$TissueOrRep == "gonad")] = "Gn"
design$TissueOrRep[which(design$TissueOrRep == "brain")] = "Br"
design$TissueOrRep[which(design$TissueOrRep == "sperm")] = "Sp"
design$TissueOrRep[which(design$TissueOrRep == "stripped_gonad")] = "SG"
head(design)
```


Compile the results tables from QuantStudio. Join results table from different runs (reverse the order to get the first timepoint first). The "micro" in the column names is non-ASCII character that we ignore using check.names=F and then convert to "u" using *iconv*. Also the column names are shifted one towards the right, so this is corrected for.
```{r}
ddPCRres = rbind(ddPCRres1, ddPCRres2, ddPCRres3, ddPCRres4, ddPCRres5, ddPCRres6, ddPCRres7, ddPCRres8, ddPCRres9, ddPCRres10, ddPCRres12, ddPCRres13, ddPCRres14)
colnames(ddPCRres) = iconv(colnames(ddPCRres), to = "ASCII", sub = "u")
colnames(ddPCRres) = c(colnames(ddPCRres)[-2],"Empty")
```


Compile table containing the input concentrations of RNA for these samples. Paste *vgll3 genotype*, *position/tank*,  *tissue/rep* and *collection date* to make up sample ID's as above. 
```{r}
conc$Sample = 'Na'

# match by sample names given in ddPCR plate (this may vary from run to run) 

# whole fry samples from 19022018
conc$Sample[which(conc$CollectionDate == '19022018')] = paste(conc$TankOrPosition[which(conc$CollectionDate == '19022018')], conc$TissueOrRep[which(conc$CollectionDate == '19022018')], sep='-')

# samples with vgll3-tank-tissue-date IDs
conc$Sample[which(conc$CollectionDate %in% c("20032018","20042018","17052018","11062018","11072018"))] = 
  paste(conc$VGLL3[which(conc$CollectionDate %in% c("20032018","20042018","17052018","11062018","11072018"))],
        conc$TankOrPosition[which(conc$CollectionDate %in% c("20032018","20042018","17052018","11062018","11072018"))], 
        conc$TissueOrRep[which(conc$CollectionDate %in% c("20032018","20042018","17052018","11062018","11072018"))], 
        conc$CollectionDate[which(conc$CollectionDate %in% c("20032018","20042018","17052018","11062018","11072018"))], sep='-')

# heterozygote samples from 11072018 with individual # IDs
conc$Sample[which(conc$CollectionDate == '11072018' & conc$VGLL3 %in% c("EL","LE"))] = 
  paste(conc$Individual[which(conc$CollectionDate == '11072018' & conc$VGLL3 %in% c("EL","LE"))], 
        conc$TissueOrRep[which(conc$CollectionDate == '11072018' & conc$VGLL3 %in% c("EL","LE"))], sep='-')

# homozygote/heterozygote samples with individual # IDs
conc$Sample[which(conc$CollectionDate %in% c("6082018","8082018","9082018","13082018","14082018", "15082018","16082018","17082018","27082018","28082018","29082018",
                                             "20092018","9102018","10102018","11102018","12102018",
                                                       "5112018","6112018","7112018","8112018","9112018","12112018",
                                                       "23112018","26112018","27112018","28112018","29112018","30112018","3122018","4122018","5122018","10012019"))] =
  paste(conc$Individual[which(conc$CollectionDate %in% c("6082018","8082018","9082018","13082018","14082018", "15082018","16082018","17082018","27082018","28082018","29082018",
                                                          "20092018","9102018","10102018","11102018","12102018",
                                                       "5112018","6112018","7112018","8112018","9112018","12112018",
                                                       "23112018","26112018","27112018","28112018","29112018","30112018","3122018","4122018","5122018","10012019"))],
        conc$TissueOrRep[which(conc$CollectionDate %in% c("6082018","8082018","9082018","13082018","14082018", "15082018","16082018","17082018","27082018","28082018","29082018",
                                                           "20092018","9102018","10102018","11102018","12102018",
                                                       "5112018","6112018","7112018","8112018","9112018","12112018",
                                                       "23112018","26112018","27112018","28112018","29112018","30112018","3122018","4122018","5122018","10012019"))],sep='-')

head(conc)
tail(conc)
```

### Duplicate concentration measurements

Second concentration reading for gonad samples that have been analyzed.
```{r}
conc2 = na.omit(conc2)
```

Merge the two concentration measurements and compare by plotting 1st concentration readings versus 2nd.
```{r}
aveConc = merge(conc,conc2,by.x="Sample",by.y="Sample")
summary(lm(aveConc$RNAconc.y ~ 0 + aveConc$RNAconc.x))
plot(aveConc$RNAconc.x,aveConc$RNAconc.y, ylab="RNA concentration measurement 1", xlab="RNA concentration measurement 2")
abline(lm(aveConc$RNAconc.y ~ 0 + aveConc$RNAconc.x))
```

Below we identify problematic samples where the 2nd concentration is much lower than the 1st (degradation of RNA more likely than false concentration on the 1st reading).
```{r}
aveConc[which(aveConc$RNAconc.x/aveConc$RNAconc.y > 2)]
# turn the second concentration reading for these samples to NA's
omitConcSamples = aveConc[which(aveConc$RNAconc.x/aveConc$RNAconc.y > 2),"Sample"]
conc2$RNAconc[which(conc2$Sample %in% omitConcSamples$Sample)] = NA
# verify
plot(conc$RNAconc[match(conc2$Sample,conc$Sample)],conc2$RNAconc, ylab="RNA concentration measurement 1", xlab="RNA concentration measurement 2")
```

Merge the two measurements again and take the average of the two measurements and use this as the value for all analyses.
```{r}
aveConc = merge(conc,conc2,by.x="Sample",by.y="Sample")
aveConc$aveRNAconc = NA 
aveConc$aveRNAconc = apply(aveConc[,c("RNAconc.x","RNAconc.y")],1,mean,na.rm=T)
aveConc
```

### Calculation of total concentrations

Intersect design table, results table and concentrations to make a master table to be used in downstream analyses.
```{r}
# Initiate empty table
countTable = setNames(data.frame(matrix(ncol = 18, nrow = 0)), c("VGLL3","FamID","TankOrPosition","Individual","TissueOrRep","Temp","CollectionDate","Sex","Sample","Ch1Conc","Ch2Conc","Ch1Copiesper20ulWell","Ch2Copiesper20ulWell","FractionalAbundance","Ch1Pos","Ch2Pos","Ch1Neg","Ch2Neg"))

# Split ddPCR results into two based on channels
ddPCRch1 = ddPCRres[which(ddPCRres$TargetType == "Ch1Unknown"),]
ddPCRch2 = ddPCRres[which(ddPCRres$TargetType == "Ch2Unknown"),]

#################################################################################
# Add ddPCR channel concentrations by mathing sample name to information on design table 

# whole fry samples from 19022018
sel = function(x) unlist(strsplit(x,"-"))[1] == "W1" 
tempConc1 = ddPCRch1[unlist(lapply(ddPCRch1$Sample,sel)),]
tempConc2 = ddPCRch2[unlist(lapply(ddPCRch2$Sample,sel)),]
tempDesign = design[unlist(lapply(design$TankOrPosition,sel)),]
tempDesign$Sample = paste(tempDesign$TankOrPosition,tempDesign$TissueOrRep,sep="-")
tempDesign$Ch1Copiesper20ulWell = tempConc1$CopiesPer20uLWell[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$Ch2Copiesper20ulWell = tempConc2$CopiesPer20uLWell[match(tempDesign$Sample,tempConc2$Sample)]
tempDesign$Ch1Conc = tempConc1$Concentration[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$Ch2Conc = tempConc2$Concentration[match(tempDesign$Sample,tempConc2$Sample)]
tempDesign$FractionalAbundance = tempConc1$FractionalAbundance[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$Ratio = tempConc1$Ratio[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$PoissonRatioMax = tempConc1$PoissonRatioMax[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$PoissonRatioMin = tempConc1$PoissonRatioMin[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$Ch1Pos = tempConc1$"Ch1+Ch2-"[match(tempDesign$Sample,tempConc1$Sample)] + tempConc1$"Ch1+Ch2+"[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$Ch2Pos = tempConc2$"Ch1-Ch2+"[match(tempDesign$Sample,tempConc2$Sample)] + tempConc2$"Ch1+Ch2+"[match(tempDesign$Sample,tempConc2$Sample)]
tempDesign$Ch1Neg = tempConc1$"Ch1-Ch2+"[match(tempDesign$Sample,tempConc1$Sample)] + tempConc1$"Ch1-Ch2-"[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$Ch2Neg = tempConc2$"Ch1+Ch2-"[match(tempDesign$Sample,tempConc2$Sample)] + tempConc2$"Ch1-Ch2-"[match(tempDesign$Sample,tempConc2$Sample)]
countTable = rbind(countTable,tempDesign)
rm(sel,tempConc1,tempConc2,tempDesign)

# samples with vgll3-tank-tissue-date IDs
sel = function(x) unlist(strsplit(x,"-"))[4] %in% c("20032018","20042018","17052018","11062018","11072018") & unlist(strsplit(x,"-"))[1] %in% c("EL","LE") == F
tempConc1 = ddPCRch1[unlist(lapply(ddPCRch1$Sample,sel)),]
tempConc2 = ddPCRch2[unlist(lapply(ddPCRch2$Sample,sel)),]
tempDesign = design[which(design$CollectionDate %in% c("20032018","20042018","17052018","11062018","11072018") & design$VGLL3 %in% c("EL","LE") == F),]
tempDesign$Sample = paste(tempDesign$VGLL3,tempDesign$TankOrPosition,tempDesign$TissueOrRep,tempDesign$CollectionDate, sep="-")
tempDesign$Ch1Copiesper20ulWell = tempConc1$CopiesPer20uLWell[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$Ch2Copiesper20ulWell = tempConc2$CopiesPer20uLWell[match(tempDesign$Sample,tempConc2$Sample)]
tempDesign$Ch1Conc = tempConc1$Concentration[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$Ch2Conc = tempConc2$Concentration[match(tempDesign$Sample,tempConc2$Sample)]
tempDesign$FractionalAbundance = tempConc1$FractionalAbundance[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$Ratio = tempConc1$Ratio[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$PoissonRatioMax = tempConc1$PoissonRatioMax[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$PoissonRatioMin = tempConc1$PoissonRatioMin[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$Ch1Pos = tempConc1$"Ch1+Ch2-"[match(tempDesign$Sample,tempConc1$Sample)] + tempConc1$"Ch1+Ch2+"[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$Ch2Pos = tempConc2$"Ch1-Ch2+"[match(tempDesign$Sample,tempConc2$Sample)] + tempConc2$"Ch1+Ch2+"[match(tempDesign$Sample,tempConc2$Sample)]
tempDesign$Ch1Neg = tempConc1$"Ch1-Ch2+"[match(tempDesign$Sample,tempConc1$Sample)] + tempConc1$"Ch1-Ch2-"[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$Ch2Neg = tempConc2$"Ch1+Ch2-"[match(tempDesign$Sample,tempConc2$Sample)] + tempConc2$"Ch1-Ch2-"[match(tempDesign$Sample,tempConc2$Sample)]
countTable = rbind(countTable,tempDesign)
rm(sel,tempConc1,tempConc2,tempDesign)

# heterozygote samples from 11072018 with individual # IDs
tempDesign = design[which(design$CollectionDate == '11072018' & design$VGLL3 %in% c("EL","LE")),]
tempDesign$Sample = paste(tempDesign$Individual,tempDesign$TissueOrRep, sep="-")
tempConc1 = ddPCRch1[which(ddPCRch1$Sample %in% tempDesign$Sample),]
tempConc2 = ddPCRch2[which(ddPCRch2$Sample %in% tempDesign$Sample),]
tempDesign$Ch1Copiesper20ulWell = tempConc1$CopiesPer20uLWell[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$Ch2Copiesper20ulWell = tempConc2$CopiesPer20uLWell[match(tempDesign$Sample,tempConc2$Sample)]
tempDesign$Ch1Conc = tempConc1$Concentration[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$Ch2Conc = tempConc2$Concentration[match(tempDesign$Sample,tempConc2$Sample)]
tempDesign$FractionalAbundance = tempConc1$FractionalAbundance[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$Ratio = tempConc1$Ratio[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$PoissonRatioMax = tempConc1$PoissonRatioMax[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$PoissonRatioMin = tempConc1$PoissonRatioMin[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$Ch1Pos = tempConc1$"Ch1+Ch2-"[match(tempDesign$Sample,tempConc1$Sample)] + tempConc1$"Ch1+Ch2+"[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$Ch2Pos = tempConc2$"Ch1-Ch2+"[match(tempDesign$Sample,tempConc2$Sample)] + tempConc2$"Ch1+Ch2+"[match(tempDesign$Sample,tempConc2$Sample)]
tempDesign$Ch1Neg = tempConc1$"Ch1-Ch2+"[match(tempDesign$Sample,tempConc1$Sample)] + tempConc1$"Ch1-Ch2-"[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$Ch2Neg = tempConc2$"Ch1+Ch2-"[match(tempDesign$Sample,tempConc2$Sample)] + tempConc2$"Ch1-Ch2-"[match(tempDesign$Sample,tempConc2$Sample)]
countTable = rbind(countTable,tempDesign)
rm(tempConc1,tempConc2,tempDesign)

# samples with individual # IDs
tempDesign = design[which(design$CollectionDate %in% c("6082018","8082018","9082018","13082018","14082018", "15082018","16082018","17082018","27082018","28082018","29082018",
                                                       "20092018","9102018","10102018","11102018","12102018",
                                                       "5112018","6112018","7112018","8112018","9112018","12112018",
                                                       "23112018","26112018","27112018","28112018","29112018","30112018","3122018","4122018","5122018","10012019")),]
tempDesign$Sample = paste(tempDesign$Individual,tempDesign$TissueOrRep, sep="-")
tempConc1 = ddPCRch1[which(ddPCRch1$Sample %in% tempDesign$Sample),]
tempConc2 = ddPCRch2[which(ddPCRch2$Sample %in% tempDesign$Sample),]
tempDesign$Ch1Copiesper20ulWell = tempConc1$CopiesPer20uLWell[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$Ch2Copiesper20ulWell = tempConc2$CopiesPer20uLWell[match(tempDesign$Sample,tempConc2$Sample)]
tempDesign$Ch1Conc = tempConc1$Concentration[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$Ch2Conc = tempConc2$Concentration[match(tempDesign$Sample,tempConc2$Sample)]
tempDesign$FractionalAbundance = tempConc1$FractionalAbundance[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$Ratio = tempConc1$Ratio[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$PoissonRatioMax = tempConc1$PoissonRatioMax[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$PoissonRatioMin = tempConc1$PoissonRatioMin[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$Ch1Pos = tempConc1$"Ch1+Ch2-"[match(tempDesign$Sample,tempConc1$Sample)] + tempConc1$"Ch1+Ch2+"[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$Ch2Pos = tempConc2$"Ch1-Ch2+"[match(tempDesign$Sample,tempConc2$Sample)] + tempConc2$"Ch1+Ch2+"[match(tempDesign$Sample,tempConc2$Sample)]
tempDesign$Ch1Neg = tempConc1$"Ch1-Ch2+"[match(tempDesign$Sample,tempConc1$Sample)] + tempConc1$"Ch1-Ch2-"[match(tempDesign$Sample,tempConc1$Sample)]
tempDesign$Ch2Neg = tempConc2$"Ch1+Ch2-"[match(tempDesign$Sample,tempConc2$Sample)] + tempConc2$"Ch1-Ch2-"[match(tempDesign$Sample,tempConc2$Sample)]
countTable = rbind(countTable,tempDesign)
rm(tempConc1,tempConc2,tempDesign)

######################################
# add starting total RNA concentration
countTable$RNAconc = conc$RNAconc[match(countTable$Sample,conc$Sample)]
#head(countTable)
#tail(countTable)

# replace with average concentrations for samples analyzed in duplicates
countTable$RNAconc[match(aveConc$Sample,countTable$Sample)] = aveConc$aveRNAconc

######################################
# remove samples that failed the ddPCR
failedSamples = c()
failed = c(which(apply(countTable[,c("Ch1Pos","Ch2Pos","Ch1Neg","Ch2Neg")],1,sum) < 10000),
                 which(countTable$Sample %in% failedSamples))
countTable = countTable[-failed,]
```

We wish to normalize the number of vgll3 copies to the amount of input RNA in each sample. Input RNA concentration is ng/ul of template (ng/ul_H2O) and ddPCR concentration is copies/ul of PCR reaction (cp/ul_PCR). 

We first calculate the concentration of template RNA in the PCR reaction: we know the original template concentration, X ng/ul_H2O. We mixed 2 ul of template with 20 ul of PCR reaction, meaning we put 2•X ng into 22 ul of PCR reaction and our concentration is 2•X/22 = N ng/ul_PCR.

Quantasoft reports target concentration as Y copies/ul_PCR. So now we have Y copies/ul_PCR and we know there was N ng/ul_PCR of template, meaning we can calculate our normalized value of expression as Y copies/N ng of template by dividing Y with N:

$\frac{Y copies}{1 ul}$ : $\frac{N ng}{1 ul}$ = $\frac{Y copies}{1 ul}$ • $\frac{1 ul}{N ng}$ = $\frac{Y copies}{N ng}$

Therefore, we need to calculate the concentration of the total RNA per PCR reaction (N = [RNAconc*2]/22) and divide ddPCR concentration with this value to get copies/ng of input RNA.  
```{r}
# correct RNA concentration for one sample with pipetting error (40ul of RT-ddPCR MM instead of 20ul)
#countTable[which(grepl("582-Ad",countTable$Sample)),"RNAconc"] = countTable[which(grepl("582-Ad",countTable$Sample)),"RNAconc"]/2

# concentration measurement for individual 918 gonad is false
#countTable[which(grepl("918-Gn",countTable$Sample)),"RNAconc"] = NA

# calculation
# countTable$Ch1copies = as.numeric(countTable$Ch1Conc)/((countTable$RNAconc*2)/22)
# countTable$Ch2copies = as.numeric(countTable$Ch2Conc)/((countTable$RNAconc*2)/22)
# countTable$TotCopies = countTable$Ch1copies + countTable$Ch2copies 
# head(countTable)

# alternative way of calculating
# # calculation of target concentration in template (2 ul)
# countTable$Ch1TemplateConc = as.numeric(countTable$Ch1Conc)*20/2
# countTable$Ch2TemplateConc = as.numeric(countTable$Ch2Conc)*20/2
# countTable$TotTemplateConc = countTable$Ch1TemplateConc+countTable$Ch2TemplateConc
# 
# # calculation of template copies per ng of input RNA
# countTable$Ch1copies = countTable$Ch1TemplateConc / countTable$RNAconc
# countTable$Ch2copies = countTable$Ch2TemplateConc / countTable$RNAconc
# countTable$TotCopies = countTable$Ch1copies + countTable$Ch2copies 
# head(countTable)
```

Alternative way of normalizing to input RNA amount. QuantaSoft also provides calculated number of template molecules per 20 ul well. Each 20 ul well contains 2 ul of template at concentration X ng/ul_H2O, meaning 2•X ng/ul_H2O total RNA per well. Copies per ng total input RNA is then simply $\frac{copies per well}{2*RNAconc}$
```{r}
# correct RNA concentration for one sample with pipetting error (40ul of RT-ddPCR MM instead of 20ul)
countTable[which(grepl("582-Ad",countTable$Sample)),"RNAconc"] = countTable[which(grepl("582-Ad",countTable$Sample)),"RNAconc"]/2

# calculation
countTable$Ch1copies = as.numeric(countTable$Ch1Copiesper20ulWell)/(countTable$RNAconc*2)
countTable$Ch2copies = as.numeric(countTable$Ch2Copiesper20ulWell)/(countTable$RNAconc*2)
countTable$TotCopies = countTable$Ch1copies + countTable$Ch2copies 
head(countTable)
```

### Analysis of total expression in tissues

Group data by sampling time-point (months post-hatching).
```{r}
#+++++++++++++++++++++++++
# Function to calculate the mean, the standard deviation and the standard error of the mean for each group
#+++++++++++++++++++++++++
# data : a data frame
# varname : the name of a column containing the variable
  #to be summariezed
# groupnames : vector of column names to be used as
  # grouping variables
data_summary <- function(data, varname, groupnames){
  require(plyr)
  summary_func <- function(x, col){
    c(mean = mean(x[[col]], na.rm=TRUE),
      sd = sd(x[[col]], na.rm=TRUE),
      sem = sd(x[[col]], na.rm=TRUE)/sqrt(length(x[[col]])),
      n = length(x[[col]]))
  }
  data_sum<-ddply(data, groupnames, .fun=summary_func,varname)
  data_sum <- plyr::rename(data_sum, c("mean" = varname))
 return(data_sum)
}

# duplicate main results table
plotData = countTable

# group all whole fry samples together (timepoint 2)
plotData$TissueOrRep = revalue(plotData$TissueOrRep, c("1"="Fry","2"="Fry","3"="Fry"))

# re-order tissues for plotting
plotData$TissueOrRep = factor(plotData$TissueOrRep,levels=c("Fry","Hr","Liv","Mc","Ad","Gn","Br","SG","Sp"))

# re-order vgll3 genotypes so that EE is plotted first
#plotData$VGLL3 = factor(plotData$VGLL3,levels=c("EE","EL","LE","LL"))

# assign gonads as factors
plotData$Gonads = factor(plotData$Gonads,levels=c("m","s","l","xl"))

# remove infinite and missing values
plotData = plotData[which(is.infinite(plotData$TotCopies) == F),]
plotData = plotData[which(is.na(plotData$TotCopies) == F),]

# add information of timepoints
plotData$Timepoint = "Na"
plotData$Timepoint[which(plotData$CollectionDate == "19022018")] = "2"
plotData$Timepoint[which(plotData$CollectionDate == "20032018")] = "3"
plotData$Timepoint[which(plotData$CollectionDate == "20042018")] = "4"
plotData$Timepoint[which(plotData$CollectionDate == "17052018")] = "5"
plotData$Timepoint[which(plotData$CollectionDate == "11062018")] = "6"
plotData$Timepoint[which(plotData$CollectionDate == "11072018")] = "7"
plotData$Timepoint[which(plotData$CollectionDate %in% c("6082018","8082018","9082018","13082018","14082018", "15082018","16082018","17082018"))] = "8"
plotData$Timepoint[which(plotData$CollectionDate %in% c("27082018","28082018","29082018","20092018"))] = "9"
plotData$Timepoint[which(plotData$CollectionDate %in% c("8102018","9102018","10102018","11102018","12102018"))] = "10"
plotData$Timepoint[which(plotData$CollectionDate %in% c("5112018","6112018","7112018","8112018","9112018","12112018"))] = "11"
plotData$Timepoint[which(plotData$CollectionDate %in% c("23112018","26112018","27112018","28112018","29112018","30112018","3122018","4122018","5122018"))] = "12"
plotData$Timepoint[which(plotData$CollectionDate %in% c("10012019"))] = "13"

# copy data to a new data frame to be modified
plotData2 = plotData
```

*Sample exclusions* Below we verify that genotypes and RT-ddPCR signal channel match (by indentifying individuals that e.g. are recorded as LL but only contain signal from the EE probes), and that recorder morphological sex matches genotypic sex. We also exclude data from heterozygotes showing monoallelic vgll3 expression, likely representing sampling or labeling error.
```{r}
# LL genotypes should express LL probe higher compared to EE probe
LLsamples = plotData2$Sample[which(plotData2$VGLL3 == "LL")]
LLinconsistent = LLsamples[which(plotData2$Ch1copies[which(plotData2$VGLL3 == "LL")] - plotData2$Ch2copies[which(plotData2$VGLL3 == "LL")] < 0)]

# EE genotypes should express EE probe higher compared to LL probe
EEsamples = plotData2$Sample[which(plotData2$VGLL3 == "EE")]
EEinconsistent = EEsamples[which(plotData2$Ch1copies[which(plotData2$VGLL3 == "EE")] - plotData2$Ch2copies[which(plotData2$VGLL3 == "EE")] > 0)]

# remove these individuals from the data
plotData2 = plotData2[which(plotData2$Sample %in% c(EEinconsistent,LLinconsistent) == F),]

# remove individuals where genotypic sex and morphology do not match
inconsistentSex = plotData2[which(plotData2$SexGeno == "male" & plotData2$Gonads == "m" | (plotData2$SexGeno == "female" & plotData2$Gonads %in% c("s","l","xl"))),]
plotData2 = plotData2[which(plotData2$Individual %in% inconsistentSex$Individual == F),]

# remove heterozygotes with monoallelic vgll3 expression
monoallelic = plotData2[which(plotData2$VGLL3 %in% c("EL","LE") & plotData2$FractionalAbundance > 90),]
plotData2 = plotData2[which(plotData2$Individual %in% monoallelic$Individual == F),]
```


Impute missing annotations on early samples (morphological differences between gonads were unrecorded for first dissectable gonads because of their very small size) and assign gonad types based on gonad morphology. Create a new column "TissueOrGonadType" where gonads are assigned to ovaries, immature testes, pubertal testes and mature testes based on their visual morphology (recorded in "Gonads" column).
```{r}
# assign missing gonad types by genotypic sex in early samples (all males before timepoint 8 are type "s" and females type "m")
missingGonadInfo = plotData2[which(is.na(plotData2$Gonads) & is.na(plotData2$SexGeno)==F), "Individual"]

for (i in missingGonadInfo){
  if (unique(plotData2$SexGeno[which(plotData2$Individual == i)]) == "female"){
    plotData2[which(plotData2$Individual == i),"Gonads"] = "m"
  }
  if (unique(plotData2$SexGeno[which(plotData2$Individual == i)]) == "male"){
    plotData2[which(plotData2$Individual == i),"Gonads"] = "s"
  }
}

# verify that assignment is correct
table(plotData2[which(plotData2$Individual%in% missingGonadInfo),c("SexGeno","Gonads")])

# assign ovaries, immature testes and mature testes by morphology 
# here l gonads considered "pubertal" because they have clearly initiated pubertal cell type changes
plotData2$TissueOrGonadType = NA
plotData2$TissueOrGonadType = as.character(plotData2$TissueOrRep)
plotData2[which(plotData2$TissueOrRep == "Gn" & plotData2$Gonads == "m"),"TissueOrGonadType"] = "Ov" 
plotData2[which(plotData2$TissueOrRep == "Gn" & plotData2$Gonads %in% c("s")),"TissueOrGonadType"] = "IT"
plotData2[which(plotData2$TissueOrRep == "Gn" & plotData2$Gonads %in% c("l")),"TissueOrGonadType"] = "PT"
plotData2[which(plotData2$TissueOrRep == "Gn" & plotData2$Gonads %in% c("xl")),"TissueOrGonadType"] = "MT"
```


Average expression across tissues and time points (with number of samples in each time-point indicated by number in grey).
```{r}
# calculate sd and sem for plotting
plotStats = data_summary(plotData2, varname="TotCopies", groupnames=c("Timepoint","TissueOrGonadType"))

# sort by timepoint
plotStats$Timepoint = factor(plotStats$Timepoint,levels=c("2","3","4","5","6","7","8","9","10","11","12","13","14"))

# order tissues
plotStats$TissueOrGonadType = factor(plotStats$TissueOrGonadType, levels=c("Fry","Hr","Liv","Mc","Ad","Br","Ov","IT","PT","MT","SG","Sp"))

# remove rows where we don't know the gonad morphology
plotStats = plotStats[is.na(plotStats$TissueOrGonadType)==F,]

# labels
tissueLabels = c('Fry'='Whole fry','Liv'="Liver",'Mc'="Muscle","Hr"="Heart","Ad"="Adipose","Gn"="Gonad", "Br"="Brain", "Ov" = "Ovary","IT"="Immature testis","PT"="Pubertal testes","MT"="Mature testis","SG"="Stripped testis","Sp"="Ejaculate")

ggplot(plotStats,aes(Timepoint, TotCopies)) +
 geom_errorbar(aes(ymin=TotCopies-sem, ymax=TotCopies+sem), width=.1,position=position_dodge(0.05)) +
 geom_line() +
 geom_point() +
 theme_classic() +
 theme(axis.title.y=element_text(size=16),
       axis.text.y=element_text(size=16),
       axis.title.x=element_text(size=16),
       axis.text.x=element_text(size=16),
       strip.text.x=element_text(size=16),
       strip.background = element_blank()) +
 geom_text(aes(y=TotCopies-15,x=as.numeric(Timepoint)+0.45,label=n),vjust=0.5,hjust=0.4,color="darkgrey", size=3) +
 facet_wrap(~TissueOrGonadType,labeller=labeller(TissueOrGonadType=tissueLabels)) +
 ylab(expression(~italic(vgll3)~'copies / ng RNA')) +
 xlab("Months post hatch") +
 #ylim(-50,250) +
 scale_x_discrete(breaks=seq(2,14,2)) 
```

Contrast immature testes expression with expression in mature testes, stripped testes and ejaculate.
```{r}
# summary stats for timepoints where we have mature gonads (for panel b)
plotStats2 = data_summary(plotData2[which(plotData2$Timepoint %in% plotData2$Timepoint[which(plotData2$TissueOrGonadType=="MT")]),], varname="TotCopies", groupnames=c("TissueOrGonadType"))

# order tissues
plotStats2$TissueOrGonadType = factor(plotStats2$TissueOrGonadType, levels=c("Fry","Hr","Liv","Mc","Ad","Br","Ov","IT","PT","MT","SG","Sp"))

ggplot(plotStats2[which(plotStats2$TissueOrGonadType %in% c("IT","MT", "Sp","SG")),],aes(TissueOrGonadType, TotCopies, color=TissueOrGonadType)) +
         geom_errorbar(aes(ymin=TotCopies-sem, ymax=TotCopies+sem), width=.1,position=position_dodge(0.05)) +
         geom_point(size=4) +
         theme_classic() +
         scale_color_manual(values=c("black","black", "black","black")) + 
         ylab(expression(~italic(vgll3)~'copies / ng RNA')) +
         theme(axis.title.x=element_blank(),
               axis.text.x=element_text(hjust=1, vjust=1, size=14,color="black",angle=65),
               axis.text.y=element_text(size=14,color="black"),
               axis.title.y=element_text(size=14,color="black"),
               legend.position="none")+
        geom_text(aes(y=TotCopies-4,x=TissueOrGonadType,label=n),vjust=1,hjust=-0.25,color="darkgrey", size=3) +
         scale_x_discrete(labels=c("Immature testes","Mature testes","Stripped testes","Ejaculate"))

# t-tests
# mature testes versus immature testes
t.test(plotData2[which(plotData2$Timepoint %in% plotData2$Timepoint[which(plotData2$TissueOrGonadType=="MT")] & plotData2$TissueOrGonadType=="MT"),"TotCopies"],
       plotData2[which(plotData2$Timepoint %in% plotData2$Timepoint[which(plotData2$TissueOrGonadType=="MT")] & plotData2$TissueOrGonadType=="IT"),"TotCopies"])

# stripped testes versus immature testes
t.test(plotData2[which(plotData2$Timepoint %in% plotData2$Timepoint[which(plotData2$TissueOrGonadType=="SG")] & plotData2$TissueOrGonadType=="SG"),"TotCopies"],
       plotData2[which(plotData2$Timepoint %in% plotData2$Timepoint[which(plotData2$TissueOrGonadType=="SG")] & plotData2$TissueOrGonadType=="IT"),"TotCopies"])

# ejaculate versus immature testes
t.test(plotData2[which(plotData2$Timepoint %in% plotData2$Timepoint[which(plotData2$TissueOrGonadType=="Sp")] & plotData2$TissueOrGonadType=="Sp"),"TotCopies"],
       plotData2[which(plotData2$Timepoint %in% plotData2$Timepoint[which(plotData2$TissueOrGonadType=="Sp")] & plotData2$TissueOrGonadType=="IT"),"TotCopies"])

# percent reduction in expression:
ITmean = mean(plotData2[which(plotData2$Timepoint %in% plotData2$Timepoint[which(plotData2$TissueOrGonadType=="MT")] & plotData2$TissueOrGonadType=="IT"),"TotCopies"],na.rm=T)
MTmean = mean(plotData2[which(plotData2$Timepoint %in% plotData2$Timepoint[which(plotData2$TissueOrGonadType=="MT")] & plotData2$TissueOrGonadType=="MT"),"TotCopies"],na.rm=T)
(ITmean-MTmean)/ITmean*100

# contrast p-value
my_comparisons5  = data.frame(start=c("IT"), 
                            end=c("MT"), 
                            y=c(120,120),
                            label=c("P=1e-06"))
my_comparisons6  = data.frame(start=c("IT"), 
                            end=c("SG"), 
                            y=c(130,130),
                            label=c("P=0.002"))
my_comparisons7  = data.frame(start=c("IT"), 
                            end=c("Sp"), 
                            y=c(140,140),
                            label=c("P=0.005"))

testesFig2 = ggplot(plotStats2[which(plotStats2$TissueOrGonadType %in% c("IT","MT", "Sp","SG")),]) +
         geom_errorbar(aes(x=TissueOrGonadType, ymin=TotCopies-sem, ymax=TotCopies+sem), width=.1,position=position_dodge(0.05)) +
         geom_point(aes(TissueOrGonadType, TotCopies, color=TissueOrGonadType),size=2) +
         theme_classic() +
         geom_signif(data=my_comparisons5,
                    aes(xmin=start, xmax=end, annotations=label, y_position=y),
                    textsize = 3, size=0.5, tip_length = 0.01,
                    manual=TRUE) +
         geom_signif(data=my_comparisons6,
                    aes(xmin=start, xmax=end, annotations=label, y_position=y),
                    textsize = 3, size=0.5, tip_length = 0.01,
                    manual=TRUE) +
         geom_signif(data=my_comparisons7,
                    aes(xmin=start, xmax=end, annotations=label, y_position=y),
                    textsize = 3, size=0.5, tip_length = 0.01,
                    manual=TRUE) +
         scale_color_manual(values=c("black","black", "black","black")) + 
         ylab(expression(~italic(vgll3)~'copies / ng RNA')) +
         theme(axis.title.x=element_blank(),
               axis.text.x=element_text(hjust=1, vjust=1, size=10,color="black",angle=65),
               axis.text.y=element_text(size=10,color="black"),
               axis.title.y=element_text(size=10,color="black"),
               legend.position="none")+
        #geom_text(aes(y=TotCopies-4,x=TissueOrGonadType,label=n),vjust=1,hjust=-0.25,color="darkgrey", size=3) +
         scale_x_discrete(labels=c("Immature testes","Mature testes","Stripped testes","Ejaculate"))

testesFig2
```

Testes expression colored by genotype.
```{r}
# collapse heterozygotes
plotData2$VGLL3RECPRO = revalue(plotData2$VGLL3,c("EL"="EL","LE"="EL"))

# summary stats for timepoints where we have mature gonads (for panel b)
plotStats2_2 = data_summary(plotData2[which(plotData2$Timepoint %in% plotData2$Timepoint[which(plotData2$TissueOrGonadType=="MT")]),], varname="TotCopies", groupnames=c("TissueOrGonadType","VGLL3RECPRO"))

# order tissues
plotStats2_2$TissueOrGonadType = factor(plotStats2_2$TissueOrGonadType, levels=c("Fry","Hr","Liv","Mc","Ad","Br","Ov","IT","PT","MT","SG","Sp"))

ggplot(plotStats2_2[which(plotStats2_2$TissueOrGonadType %in% c("IT","MT")),], aes(TissueOrGonadType, TotCopies, color=VGLL3RECPRO)) +
         geom_point(size=3, position=position_dodge(width=0.05)) +
         geom_errorbar(aes(ymin=TotCopies-sem, ymax=TotCopies+sem), width=.1, position=position_dodge(width=0.05)) +
         theme_classic() +
         scale_color_manual(name=expression(~italic("vgll3")),values=c("blue","darkgrey","red")) +
         ylab(expression(~italic(vgll3)~'copies / ng RNA')) +
         theme(axis.title.x=element_blank(),
               axis.text.x=element_text(hjust=1, vjust=1, size=14,color="black",angle=65),
               axis.text.y=element_text(size=14,color="black"),
               axis.title.y=element_text(size=14,color="black"),
               legend.position="right",
               legend.title = element_text(size=14))+
        #geom_text(aes(y=TotCopies-4,x=TissueOrGonadType,label=n),vjust=1,hjust=-0.25,color="darkgrey", size=3) +
         scale_x_discrete(labels=c("Immature testes","Mature testes","Stripped testes","Ejaculate"))

```


## Analysis in context of Amh/Igf3 expression ratio and clustering

Combine result tables from QuantStudio output.
Channel 1 (FAM) Amh
channel 2 (VIC) Igf3
```{r}
# load results
ddPCRresAmhIgf1 = rbind(ddPCRresAmhIgf1_1,ddPCRresAmhIgf1_2)

colnames(ddPCRresAmhIgf1) = iconv(colnames(ddPCRresAmhIgf1), to = "ASCII", sub = "u")
colnames(ddPCRresAmhIgf1) = c(colnames(ddPCRresAmhIgf1)[-2],"Empty")

# calculate efficiency factor for amh (FAM) vs igf3 (VIC) probes based on genomic DNA
amhPos = ddPCRresAmhIgf1$Positives[which(ddPCRresAmhIgf1$Sample=="e28-DNA" & ddPCRresAmhIgf1$TargetType=="Ch1Positive Control")]
amhNeg = ddPCRresAmhIgf1$Negatives[which(ddPCRresAmhIgf1$Sample=="e28-DNA" & ddPCRresAmhIgf1$TargetType=="Ch1Positive Control")]
igf3Pos = ddPCRresAmhIgf1$Positives[which(ddPCRresAmhIgf1$Sample=="e28-DNA" & ddPCRresAmhIgf1$TargetType=="Ch2Positive Control")]
igf3Neg = ddPCRresAmhIgf1$Negatives[which(ddPCRresAmhIgf1$Sample=="e28-DNA" & ddPCRresAmhIgf1$TargetType=="Ch2Positive Control")]
effFactor2 = mean(apply(cbind(amhPos,igf3Pos,amhNeg,igf3Neg), 1, function(x) poisson.test(c(x[1],x[2]),c(x[3],x[4]))[["estimate"]]))
```

Flag and remove amh/igf3 expression data for samples where the RT-ddPCR reaction was saturated with positive droplets (more than 98% of all droplets are positive). In the figures below, saturated samples are plotted in red and non-saturated samples in black.
```{r}
# can't use reactions where positives are saturated
plot(ddPCRresAmhIgf1$Positives,
     ddPCRresAmhIgf1$Positives+ddPCRresAmhIgf1$Negatives,
     col=as.factor(ddPCRresAmhIgf1$Positives/c(ddPCRresAmhIgf1$Negatives+ddPCRresAmhIgf1$Positives)>0.98),ylab="Total number of droplets", xlab="Number of positive droplets")

# add a flag to samples that are not saturated
ddPCRresAmhIgf1$AbundanceOK = ddPCRresAmhIgf1$Positives/c(ddPCRresAmhIgf1$Negatives+ddPCRresAmhIgf1$Positives)<0.98

# clean for samples where Amh:Igf1 reaction has OK negative:positive droplet numbers (on both channels)
ddPCRresAmhIgf1$Remove = NA
for (i in unique(ddPCRresAmhIgf1$Sample[which(ddPCRresAmhIgf1$Sample != "H2O")])){
  tmp = ddPCRresAmhIgf1[which(ddPCRresAmhIgf1$Sample == i),]
  if (all(tmp$AbundanceOK)){
    ddPCRresAmhIgf1[which(ddPCRresAmhIgf1$Sample == i),"Remove"] = FALSE
  }
  else {
    ddPCRresAmhIgf1[which(ddPCRresAmhIgf1$Sample == i),"Remove"] = TRUE
  }
}
ddPCRresAmhIgf1$Remove[which(ddPCRresAmhIgf1$Sample == "H2O")] = TRUE

ddPCRresAmhIgf1Clean = data.frame()
ddPCRresAmhIgf1Clean = ddPCRresAmhIgf1[which(ddPCRresAmhIgf1$Remove == F),]

# verify that abundances are now ok
plot(ddPCRresAmhIgf1Clean$Positives,
     ddPCRresAmhIgf1Clean$Positives+ddPCRresAmhIgf1Clean$Negatives,
     col=as.factor(ddPCRresAmhIgf1Clean$Positives/c(ddPCRresAmhIgf1Clean$Negatives+ddPCRresAmhIgf1Clean$Positives)>0.98),ylab="Total number of droplets", xlab="Number of positive droplets")
```

Clean for any samples not analyzed for vgll3.
```{r}
# clean results from any samples not analyzed for vgll3
ddPCRresAmhIgf1Clean = ddPCRresAmhIgf1Clean[which(ddPCRresAmhIgf1Clean$Sample %in% plotData2$Sample),]
```

Split the ddPCR output file into Channel 1 (Amh) and Channel 2 (Igf1) to avoid bugs.
```{r}
ddPCRresAmh = ddPCRresAmhIgf1Clean[which(ddPCRresAmhIgf1Clean$TargetType == "Ch1Unknown"),]
ddPCRresIgf1 = ddPCRresAmhIgf1Clean[which(ddPCRresAmhIgf1Clean$TargetType == "Ch2Unknown"),]
```

Add columns to existing data frame.
```{r}
# add Amh copies
plotData2$AmhConc = NA
plotData2$AmhConc = replace(plotData2$AmhConc,
                              match(ddPCRresAmh$Sample,
                                    plotData2$Sample),
                              ddPCRresAmh$CopiesPer20uLWell[which(ddPCRresAmh$TargetType == "Ch1Unknown")])

# add Igf1 copies
plotData2$Igf1Conc = NA
plotData2$Igf1Conc = replace(plotData2$Igf1Conc,
                              match(ddPCRresIgf1$Sample,
                                    plotData2$Sample),
                              ddPCRresIgf1$CopiesPer20uLWell[which(ddPCRresIgf1$TargetType == "Ch2Unknown")])

# Amh:Igf1 ratio of concentrations
plotData2$AmhIgf1ConcRatio = NA
plotData2$AmhIgf1ConcRatio = replace(plotData2$AmhIgf1ConcRatio,
                              match(ddPCRresAmh$Sample,
                                    plotData2$Sample),
                              ddPCRresAmh$Ratio[which(ddPCRresAmh$TargetType == "Ch1Unknown")])

# Amh:Igf1 fractional abundance
plotData2$AmhIgf1FA = NA
plotData2$AmhIgf1FA = replace(plotData2$AmhIgf1FA,
                              match(ddPCRresAmh$Sample,
                                    plotData2$Sample),
                              ddPCRresAmh$FractionalAbundance[which(ddPCRresAmh$TargetType == "Ch1Unknown")])

# Amh:igf3 Poisson rate ratio (add speudocount +1 to avoid Inf values)
plotData2$AmhPos = NA
plotData2$AmhPos = replace(plotData2$AmhPos,
                              match(ddPCRresAmh$Sample,
                                    plotData2$Sample),
                              ddPCRresAmh$"Ch1+Ch2-"[which(ddPCRresAmh$TargetType == "Ch1Unknown")])
plotData2$AmhPos = plotData2$AmhPos + 1
  
plotData2$Igf1Pos = NA
plotData2$Igf1Pos = replace(plotData2$Igf1Pos,
                              match(ddPCRresIgf1$Sample,
                                    plotData2$Sample),
                              ddPCRresIgf1$"Ch1-Ch2+"[which(ddPCRresIgf1$TargetType == "Ch2Unknown")])
plotData2$Igf1Pos = plotData2$Igf1Pos + 1

plotData2$AmhNeg = NA
plotData2$AmhNeg = replace(plotData2$AmhNeg,
                              match(ddPCRresAmh$Sample,
                                    plotData2$Sample),
                              c(ddPCRresAmh$"Ch1-Ch2-"[which(ddPCRresAmh$TargetType == "Ch1Unknown")]+ddPCRresAmh$"Ch1-Ch2+"[which(ddPCRresAmh$TargetType == "Ch1Unknown")]))
plotData2$AmhNeg = plotData2$AmhNeg + 1

plotData2$Igf1Neg = NA
plotData2$Igf1Neg = replace(plotData2$Igf1Neg,
                              match(ddPCRresIgf1$Sample,
                                    plotData2$Sample),
                              c(ddPCRresIgf1$"Ch1-Ch2-"[which(ddPCRresIgf1$TargetType == "Ch2Unknown")]+ddPCRresIgf1$"Ch1+Ch2-"[which(ddPCRresIgf1$TargetType == "Ch2Unknown")]))
plotData2$Igf1Neg = plotData2$Igf1Neg + 1
```

Subset for males with Amh and Igf1 expression data
```{r}
plotData3 = plotData2[which(is.na(plotData2$AmhIgf1ConcRatio) == F & plotData2$Gonads != "m"),]

plotData3$AmhIgf3PoissonRateRatio = apply(plotData3[,c("AmhPos","Igf1Pos","AmhNeg","Igf1Neg")], 1, function(x) poisson.test(c(x[1],x[2]),c(x[3],x[4]),r=effFactor2)[["estimate"]])
```


Amh/igf3 log- Poisson rate ratio in mature, immature and pubertal gonads.
```{r}
# tissue labels
tissueLabels = c('IT'="Immature","PT"="Pubertal","MT"="Mature")

ggplot(plotData3)+
  geom_quasirandom(aes(TissueOrGonadType,log(AmhIgf3PoissonRateRatio, base=2),color=TissueOrGonadType),size=5,alpha=0.7) +
  scale_color_manual(values=c("darkgrey","black","grey")) + 
  theme_classic() +
  facet_wrap(~VGLL3) +
  scale_x_discrete(breaks=c("IT","PT","MT"),labels=c("Immature","Pubertal","Mature")) + 
  theme(axis.text.x=element_text(size=14,color="black"),
       axis.title.x=element_blank(),
       axis.ticks.x=element_blank(),
       axis.text.y=element_text(size=14,color="black"),
       axis.title.y=element_text(size=14,color="black"),
       axis.line.x = element_blank(),
       legend.position="none") +
  ylab("log2 amh/igf3 rate ratio") + 
  xlab("Gonad maturity")
```

Unsupervised clustering of samples according to amh/igf3 ratio and vgll3 expression using Gaussian mixture modeling as implemented in the package Mclust.
```{r}
clustTable = as.data.frame(cbind(plotData3$TotCopies,log(plotData3$AmhIgf3PoissonRateRatio, base=2)))

BIC = mclustBIC(clustTable)
summary(BIC)
mod = Mclust(clustTable,x=BIC)
summary(mod,parameter=T)

# plots of best model/classification
fviz_mclust(mod, "BIC", palette = "jco")
fviz_mclust(mod, "classification", geom = "point", pointsize = 1.5, palette = "jco", stand=F, ellipse.level=0.68)
fviz_mclust(mod, "uncertainty", palette = "jco", stand=F)
plot(mod,what="density")
plot(mod, what = "density", type = "persp")

LRT = mclustBootstrapLRT(clustTable, nboot=999, modelName = "VEI")
LRT
```

Plot of Mclust best classification (vs vgll3 genotype).
```{r}
# visualize Mclust clustering confidence interval by adding a layer of this ggplot object
ellipsePlot = fviz_mclust(mod, "classification", geom = "point", stand=F, ellipse.level=0.95) 

ggplot(plotData3)+
  geom_point(aes(TotCopies,log(AmhIgf3PoissonRateRatio, base=2),shape=TissueOrGonadType,fill=VGLL3),size=4,alpha=0.7)+
  #scale_color_manual(values=c("black","darkgrey")) + 
  scale_fill_manual(values=c("grey","grey", "blue","darkgreen","orange","red")) +
  #scale_fill_manual(values=c("slategrey","grey")) +
  scale_shape_manual(values=c(21,22,23)) +
  ellipsePlot$layers[[2]] +
  geom_label(aes(x=60,y=-4), label="Mature cluster") +
  geom_label(aes(x=120,y=22), label="Immature cluster") +
  theme_classic() +
  ylab("log2 amh/igf3 rate ratio") + 
  xlab("vgll3 copies / ng RNA") +
  theme(axis.text=element_text(size=14,color="black"),
        axis.title=element_text(size=14,color="black"),
        legend.position="top")
```

Plot of Mclust best model classification. We observe that "pubertal" samples show a large range of amh/igf3 expression ratio and vgll3 expression, indicating that the class contains samples that were either in ht early stages of transitioning towards maturity, or samples that had reached maturity while resembling immature gonads more than mature gonads in morphology.
```{r}
plotData3$clust = mod$classification

bxFig = ggplot(plotData3)+
          geom_quasirandom(aes(TissueOrGonadType,log(AmhIgf3PoissonRateRatio, base=2),shape=TissueOrGonadType),color="black",fill="darkgrey",size=5,alpha=0.7) +
          scale_fill_manual(values=c("darkgrey","darkgrey")) + 
          scale_shape_manual(values=c(21,22,23)) +
          theme_classic() +
          scale_x_discrete(breaks=c("IT","MT","PT"),labels=c("Immature","Mature","Pubertal")) + 
          theme(axis.text.x=element_text(size=14,color="black"),
               axis.title.x=element_blank(),
               axis.ticks.x=element_blank(),
               axis.text.y=element_text(size=14,color="black"),
               axis.title.y=element_text(size=14,color="black"),
               axis.line.x = element_blank(),
               legend.position="none") +
          ylab("log2 amh/igf3 rate ratio") + 
          xlab("Gonad maturity")

scatterFig = ggplot(plotData3)+
  geom_point(aes(TotCopies,log(AmhIgf3PoissonRateRatio, base=2),shape=TissueOrGonadType,fill=as.factor(clust)),size=4,alpha=0.7)+
  scale_color_manual(values=c("black","darkgrey")) + 
  scale_fill_manual(values=c("slategrey","grey")) +
  scale_shape_manual(values=c(21,22,23)) +
  ellipsePlot$layers[[2]] +
  geom_label(aes(x=60,y=-4), label="Mature cluster") +
  geom_label(aes(x=120,y=22), label="Immature cluster") +
  theme_classic() +
  ylab("log2 amh/igf3 rate ratio") + 
  xlab("vgll3 copies / ng RNA") +
  theme(axis.text=element_text(size=14,color="black"),
        axis.title=element_text(size=14,color="black"),
        legend.position="none")

ggarrange(bxFig, scatterFig, 
          labels = c("a","b"),
          ncol = 2, nrow = 1,
          widths = c(0.7, 1))

scatterFig2 = ggplot(plotData3)+
  ellipsePlot$layers[[2]] +
  geom_point(aes(TotCopies,log(AmhIgf3PoissonRateRatio, base=2),shape=TissueOrGonadType,fill=as.factor(clust)),size=3,alpha=0.7)+
  scale_color_manual(values=c("darkgrey","black")) + 
  scale_fill_manual(values=c("grey", "slategrey")) +
  scale_shape_manual(values=c(21,22,23)) +
  geom_label(aes(x=60,y=-4), label="Mature cluster",size=3) +
  geom_label(aes(x=120,y=22), label="Immature cluster",size=3) +
  theme_classic() +
  ylab(expression('log2'~italic(amh/igf3)~'rate ratio')) +
  xlab(expression(~italic(vgll3)~'copies / ng RNA')) +
  theme(axis.text=element_text(size=10,color="black"),
        axis.title=element_text(size=10,color="black"),
        legend.position="none")

scatterFig2
```


## Modeling EE versus LL and immature versus mature gonad expression

Below we calculate age for samples based on date when 90% of eggs had hatched.
```{r}
convertToDate = function(x){
  x = as.character(x)
  if (strsplit(x,"") %>% unlist() %>% length() == 8){
      day = unlist(strsplit(x,""))[1:2] %>% as.numeric() %>% as.matrix() %>% t() %>% as.data.frame() %>% unite("d",c("V1","V2"),sep="")
      month = unlist(strsplit(x,""))[3:4] %>% as.numeric() %>% as.matrix() %>% t() %>% as.data.frame() %>% unite("m",c("V1","V2"),sep="")
      year = unlist(strsplit(x,""))[5:8] %>% as.numeric() %>% as.matrix() %>% t() %>% as.data.frame() %>% unite("y",c("V1","V2","V3","V4"),sep="")
      return(paste(month,day,year,sep="/"))
  }
  else {
      day = paste(0,unlist(strsplit(x,""))[1],sep="") 
      month = unlist(strsplit(x,""))[2:3] %>% as.numeric() %>% as.matrix() %>% t() %>% as.data.frame() %>% unite("m",c("V1","V2"),sep="")
      year = unlist(strsplit(x,""))[4:7] %>% as.numeric() %>% as.matrix() %>% t() %>% as.data.frame() %>% unite("y",c("V1","V2","V3","V4"),sep="")
      return(paste(month,day,year,sep="/"))
  }
}

daysPostHatch = function(x){
  datetimes = strptime(c("12/29/2017",x), format = "%m/%d/%Y")
  return(as.numeric(difftime(datetimes[2], datetimes[1], units = "days")))
}

# OUL warm 90% hatch date 12/29/2017
plotData2$Date = unlist(lapply(plotData2$CollectionDate,convertToDate))
plotData2$DPH = unlist(lapply(plotData2$Date,daysPostHatch))
plotData3$Date = unlist(lapply(plotData3$CollectionDate,convertToDate))
plotData3$DPH = unlist(lapply(plotData3$Date,daysPostHatch))
```

Plot of tissue expression levels as above but where age is converted from months post-hatch to days post-hatch.
```{r}
plotData2$TissueOrGonadType = factor(plotData2$TissueOrGonadType, levels=c("Fry","Hr","Liv","Mc","Ad","Br","Ov","IT","PT", "MT","SG","Sp"))

tissueLabels = c('Fry'='Whole fry','Liv'="Liver",'Mc'="Muscle","Hr"="Heart","Ad"="Adipose","Gn"="Gonad", "Br"="Brain", "Ov" = "Immature ovary","IT"="Immature testes","MT"="Mature testes","SG"="Stripped testes","Sp"="Ejaculate")

allFigDPH = ggplot(plotData2[which(plotData2$TissueOrGonadType %in% c("Hr","Liv","Mc","Ad","Ov","IT","MT","SG","Sp")),],aes(DPH, TotCopies)) +
 geom_smooth(data = plotData2[which(plotData2$TissueOrGonadType %in% c("Hr","Liv","Mc","Ad","Ov","IT","MT")),],fill="grey70",colour="aquamarine4") +
 geom_point(colour="black",fill="black",alpha=0.25) +
 theme_classic() +
 theme(axis.text=element_text(size=10,color="black"),
       axis.title=element_text(size=10,color="black"),
       strip.text.x=element_text(size=10),
       legend.position="none",
       strip.background = element_blank()) +
# scale_x_discrete(breaks=seq(2,14,2)) +
 ylim(0,350) +
 facet_wrap(~TissueOrGonadType,labeller=labeller(TissueOrGonadType=tissueLabels)) +
 ylab(expression(~italic(vgll3)~'copies / ng RNA')) +
 xlab("Age (days post-hatch)") 
allFigDPH

# color points according to vgll3 genotypes
allFigDPHcolor = ggplot(plotData2[which(plotData2$TissueOrGonadType %in% c("Hr","Liv","Mc","Ad","Ov","IT")),],aes(DPH, TotCopies,colour=VGLL3RECPRO)) +
 geom_smooth(fill="grey70",colour="aquamarine4") +
 geom_point(fill="black",alpha=0.25) +
 scale_color_manual(name=expression(~italic("vgll3")), values=c( "blue","darkgrey","red")) +
 theme_classic() +
 theme(axis.text=element_text(size=14,color="black"),
       axis.title=element_text(size=14,color="black"),
       strip.text.x=element_text(size=14),
       legend.position="right",
       strip.background = element_blank()) +
# scale_x_discrete(breaks=seq(2,14,2)) +
 facet_wrap(~TissueOrGonadType,labeller=labeller(TissueOrGonadType=tissueLabels)) +
 ylab(expression(~italic(vgll3)~'copies / ng RNA')) +
 xlab("Age (days post-hatch)") 
allFigDPHcolor
```

Reordered panels for presentation.
```{r}
plotData2$TissueOrGonadType = factor(plotData2$TissueOrGonadType, levels=c("Ad","Ov","IT","Hr","Liv","Mc","Br","Fry","PT", "SG","MT","Sp"))

tissueLabels = c('Fry'='Whole fry','Liv'="Liver",'Mc'="Muscle","Hr"="Heart","Ad"="Adipose","Gn"="Gonad", "Br"="Brain", "Ov" = "Immature ovary","IT"="Immature testes","MT"="Mature testes","SG"="Stripped testes","Sp"="Ejaculate")

ggplot(plotData2[which(plotData2$TissueOrGonadType %in% c("Hr","Liv","Mc","Ad","Ov","IT")),],aes(DPH, TotCopies)) +
 geom_smooth(fill="grey70",colour="aquamarine4") +
 geom_point(colour="black",fill="black",alpha=0.25) +
 theme_classic() +
 theme(axis.text=element_text(size=10,color="black"),
       axis.title=element_text(size=10,color="black"),
       strip.text.x=element_text(size=10),
       legend.position="none",
       strip.background = element_blank()) +
# scale_x_discrete(breaks=seq(2,14,2)) +
 facet_wrap(~TissueOrGonadType,labeller=labeller(TissueOrGonadType=tissueLabels)) +
 ylab(expression(~italic(vgll3)~'copies / ng RNA')) +
 xlab("Age (days post-hatch)") 
```

```{r}
ggarrange(allFigDPH,                                                 
          ggarrange(testesFig2,scatterFig2, ncol = 2, labels = c("B", "C"), widths = c(0.7, 1)),
          nrow = 2, heights=c(1.5,1),
          labels = "A")
```

Plot all tissues for supplemental figure.
```{r}
# re-order tissues
plotData2$TissueOrGonadType = factor(plotData2$TissueOrGonadType, levels=c("Fry","Hr","Liv","Mc","Ad","Br","Ov","IT","PT", "SG","MT","Sp"))
tissueLabels = c('Fry'='Whole fry','Liv'="Liver",'Mc'="Muscle","Hr"="Heart","Ad"="Adipose","Gn"="Gonad", "Br"="Brain", "Ov" = "Immature ovary","IT"="Immature testes","PT"="Pubertal testes", "MT"="Mature testes","SG"="Stripped testes","Sp"="Ejaculate")

# color points according to vgll3 genotypes
allTissuesGenotypes = ggplot(plotData2,aes(DPH, TotCopies,colour=VGLL3RECPRO)) +
 #geom_smooth(fill="grey70",colour="aquamarine4") +
 geom_point(fill="black",alpha=0.25) +
 scale_color_manual(name=expression(~italic("vgll3")), values=c( "blue","darkgrey","red")) +
 theme_classic() +
 theme(axis.text=element_text(size=14,color="black"),
       axis.title=element_text(size=14,color="black"),
       strip.text.x=element_text(size=14),
       legend.position="right",
       strip.background = element_blank()) +
# scale_x_discrete(breaks=seq(2,14,2)) +
 facet_wrap(~TissueOrGonadType,labeller=labeller(TissueOrGonadType=tissueLabels)) +
 ylab(expression(~italic(vgll3)~'copies / ng RNA')) +
 xlab("Age (days post-hatch)")
allTissuesGenotypes
```


Individual plots
```{r}
tissueLabels = c('Fry'='Whole fry','Liv'="Liver",'Mc'="Muscle","Hr"="Heart","Ad"="Adipose","Gn"="Gonad", "Br"="Brain", "Ov" = "Immature ovary","IT"="Immature testes","MT"="Mature testes","SG"="Stripped testes","Sp"="Ejaculate")

ggplot(plotData2[which(plotData2$TissueOrGonadType %in% c("Hr")),],aes(DPH, TotCopies)) +
 geom_smooth(fill="grey70",colour="aquamarine4") +
 geom_point(colour="black",fill="black",alpha=0.25,size=4) +
 theme_classic() +
 theme(axis.text=element_text(size=16,color="black"),
       axis.title=element_text(size=16,color="black"),
       strip.text.x=element_text(size=16),
       legend.position="none",
       strip.background = element_blank()) +
# scale_x_discrete(breaks=seq(2,14,2)) +
 facet_wrap(~TissueOrGonadType,labeller=labeller(TissueOrGonadType=tissueLabels)) +
 ylab(expression(~italic(vgll3)~'copies / ng RNA')) +
 xlab("Age (days post-hatch)") 

ggplot(plotData2[which(plotData2$TissueOrGonadType %in% c("Liv")),],aes(DPH, TotCopies, color=VGLL3)) +
 geom_smooth(fill="grey70",colour="aquamarine4") +
 geom_point(colour="black",fill="black",alpha=0.25,size=4) +
 theme_classic() +
 theme(axis.text=element_text(size=16,color="black"),
       axis.title=element_text(size=16,color="black"),
       strip.text.x=element_text(size=16),
       legend.position="none",
       strip.background = element_blank()) +
# scale_x_discrete(breaks=seq(2,14,2)) +
 facet_wrap(~TissueOrGonadType,labeller=labeller(TissueOrGonadType=tissueLabels)) +
 ylab(expression(~italic(vgll3)~'copies / ng RNA')) +
 xlab("Age (days post-hatch)") 

ggplot(plotData2[which(plotData2$TissueOrGonadType %in% c("Mc")),],aes(DPH, TotCopies)) +
 geom_smooth(fill="grey70",colour="aquamarine4") +
 geom_point(colour="black",fill="black",alpha=0.25,size=4) +
 theme_classic() +
 theme(axis.text=element_text(size=16,color="black"),
       axis.title=element_text(size=16,color="black"),
       strip.text.x=element_text(size=16),
       legend.position="none",
       strip.background = element_blank()) +
# scale_x_discrete(breaks=seq(2,14,2)) +
 facet_wrap(~TissueOrGonadType,labeller=labeller(TissueOrGonadType=tissueLabels)) +
 ylab(expression(~italic(vgll3)~'copies / ng RNA')) +
 xlab("Age (days post-hatch)")

ggplot(plotData2[which(plotData2$TissueOrGonadType %in% c("Ad")),],aes(DPH, TotCopies)) +
 geom_smooth(fill="grey70",colour="aquamarine4") +
 geom_point(colour="black",fill="black",alpha=0.25,size=4) +
 theme_classic() +
 theme(axis.text=element_text(size=16,color="black"),
       axis.title=element_text(size=16,color="black"),
       strip.text.x=element_text(size=16),
       legend.position="none",
       strip.background = element_blank()) +
# scale_x_discrete(breaks=seq(2,14,2)) +
 facet_wrap(~TissueOrGonadType,labeller=labeller(TissueOrGonadType=tissueLabels)) +
 ylab(expression(~italic(vgll3)~'copies / ng RNA')) +
 xlab("Age (days post-hatch)") 

ggplot(plotData2[which(plotData2$TissueOrGonadType %in% c("Ov")),],aes(DPH, TotCopies)) +
 geom_smooth(fill="grey70",colour="aquamarine4") +
 geom_point(colour="black",fill="black",alpha=0.25,size=4) +
 theme_classic() +
 theme(axis.text=element_text(size=16,color="black"),
       axis.title=element_text(size=16,color="black"),
       strip.text.x=element_text(size=16),
       legend.position="none",
       strip.background = element_blank()) +
# scale_x_discrete(breaks=seq(2,14,2)) +
 facet_wrap(~TissueOrGonadType,labeller=labeller(TissueOrGonadType=tissueLabels)) +
 ylab(expression(~italic(vgll3)~'copies / ng RNA')) +
 xlab("Age (days post-hatch)") 


ggplot(plotData2[which(plotData2$TissueOrGonadType %in% c("IT")),],aes(DPH, TotCopies)) +
 geom_smooth(fill="grey70",colour="aquamarine4") +
 geom_point(colour="black",fill="black",alpha=0.25,size=4) +
 theme_classic() +
 theme(axis.text=element_text(size=16,color="black"),
       axis.title=element_text(size=16,color="black"),
       strip.text.x=element_text(size=16),
       legend.position="none",
       strip.background = element_blank()) +
# scale_x_discrete(breaks=seq(2,14,2)) +
 facet_wrap(~TissueOrGonadType,labeller=labeller(TissueOrGonadType=tissueLabels)) +
 ylab(expression(~italic(vgll3)~'copies / ng RNA')) +
 xlab("Age (days post-hatch)") 
```


Average plots.
```{r}
sexLabels = c('m'="Female","s"="Immature male","xl"="Mature male")

ggplot(plotData2[which(plotData2$TissueOrGonadType %in% c("Hr")),],aes(VGLL3, TotCopies,color=VGLL3)) +
 geom_boxplot(aes(VGLL3,TotCopies,color=VGLL3),outlier.alpha = 0) + 
 stat_summary(fun.y=mean, geom="point", shape=17, size=4, color="black", fill="black",alpha=0.5) +
 geom_quasirandom(alpha=0.25,size=4) +
 scale_color_manual(values=c("blue","red")) +
 theme_classic() +
 theme(axis.text=element_text(size=16,color="black"),
       axis.title=element_text(size=16,color="black"),
       strip.text.x=element_text(size=16),
       legend.position="none",
       strip.background = element_blank()) +
 facet_wrap(~Gonads,labeller=labeller(Gonads=sexLabels)) +
 ylab(expression(~italic(vgll3)~'copies / ng RNA')) +
 ggtitle("Heart") +
 xlab("")

ggplot(plotData2[which(plotData2$TissueOrGonadType %in% c("Ad")),],aes(VGLL3, TotCopies,color=VGLL3)) +
 geom_boxplot(aes(VGLL3,TotCopies,color=VGLL3),outlier.alpha = 0) +
 stat_summary(fun.y=mean, geom="point", shape=17, size=4, color="black", fill="black",alpha=0.5) +
 geom_quasirandom(alpha=0.25,size=4) +
 scale_color_manual(values=c("blue","darkgreen", "darkorange","red")) +
 theme_classic() +
 theme(axis.text=element_text(size=16,color="black"),
       axis.title=element_text(size=16,color="black"),
       strip.text.x=element_text(size=16),
       legend.position="none",
       strip.background = element_blank()) +
 facet_wrap(~Gonads,labeller=labeller(Gonads=sexLabels)) +
 ylab(expression(~italic(vgll3)~'copies / ng RNA')) +
 ggtitle("Adipose") +
 xlab("")

#t.test(plotData2$TotCopies[which(plotData2$VGLL3=="EE" & plotData2$TissueOrGonadType %in% c("Ad"))],plotData2$TotCopies[which(plotData2$VGLL3=="LL" & plotData2$TissueOrGonadType %in% c("Ad"))])

ggplot(plotData2[which(plotData2$TissueOrGonadType %in% c("Mc")),],aes(VGLL3, TotCopies,color=VGLL3)) +
 geom_boxplot(aes(VGLL3,TotCopies,color=VGLL3),outlier.alpha = 0) + 
 stat_summary(fun.y=mean, geom="point", shape=17, size=4, color="black", fill="black",alpha=0.5) +
 geom_quasirandom(alpha=0.25,size=4) +
 scale_color_manual(values=c("blue","red")) +
 theme_classic() +
 theme(axis.text=element_text(size=16,color="black"),
       axis.title=element_text(size=16,color="black"),
       strip.text.x=element_text(size=16),
       legend.position="none",
       strip.background = element_blank()) +
 facet_wrap(~Gonads,labeller=labeller(Gonads=sexLabels)) +
 ylab(expression(~italic(vgll3)~'copies / ng RNA')) +
 ggtitle("Muscle") +
  xlab("")

ggplot(plotData2[which(plotData2$TissueOrGonadType %in% c("Ov")),],aes(VGLL3, TotCopies,color=VGLL3)) +
 geom_boxplot(aes(VGLL3,TotCopies,color=VGLL3),outlier.alpha = 0) + 
 stat_summary(fun.y=mean, geom="point", shape=17, size=4, color="black", fill="black",alpha=0.5) +
 geom_quasirandom(alpha=0.25,size=4) +
 scale_color_manual(values=c("blue","darkgreen", "darkorange","red")) +
 theme_classic() +
 theme(axis.text=element_text(size=16,color="black"),
       axis.title=element_text(size=16,color="black"),
       strip.text.x=element_text(size=16),
       legend.position="none",
       strip.background = element_blank()) +
 facet_wrap(~Gonads,labeller=labeller(Gonads=sexLabels)) +
 ylab(expression(~italic(vgll3)~'copies / ng RNA')) +
 ggtitle("Ovaries") +
  xlab("")

ggplot(plotData2[which(plotData2$TissueOrGonadType %in% c("IT")),],aes(VGLL3, TotCopies,color=VGLL3)) +
 geom_boxplot(aes(VGLL3,TotCopies,color=VGLL3),outlier.alpha = 0) + 
 stat_summary(fun.y=mean, geom="point", shape=17, size=4, color="black", fill="black",alpha=0.5) +
 geom_quasirandom(alpha=0.25,size=4) +
 scale_color_manual(values=c("blue","darkgreen", "darkorange","red")) +
 theme_classic() +
 theme(axis.text=element_text(size=16,color="black"),
       axis.title=element_text(size=16,color="black"),
       strip.text.x=element_text(size=16),
       legend.position="none",
       strip.background = element_blank()) +
 facet_wrap(~Gonads,labeller=labeller(Gonads=sexLabels)) +
 ylab(expression(~italic(vgll3)~'copies / ng RNA')) +
 ggtitle("Immature testes") +
 xlab("")

```


Plot of EE & LL vgll3 expression time-series (immature gonads). Time-point 13 is exlcuded because there's no immature EE data.
```{r}
# subset for male gonads (immature) and only EE & LL genotypes - exclude TP 13 because only LL data
plotData4 = plotData2[which(plotData2$TissueOrRep=="Gn" & plotData2$TissueOrGonadType %in% c("IT") & plotData2$VGLL3 %in% c("EE","LL") & plotData2$Timepoint!=13),]

# fit polynomial regression - test models with and without interaction term 
fitSimple <- lm(plotData4$TotCopies ~ plotData4$VGLL3 + plotData4$DPH)
fitFull <- lm(plotData4$TotCopies ~  plotData4$VGLL3 + plotData4$VGLL3 * plotData4$DPH)
fitVgll3 = lm(plotData4$TotCopies ~ plotData4$VGLL3)
fitAge = lm(plotData4$TotCopies ~ plotData4$DPH)

# select for the best model
AIC(fitSimple,fitFull)
BIC(fitSimple,fitFull)

# significance of the models
summary(fitSimple)
summary(fitFull)
summary(fitVgll3)
summary(fitAge)

# are the slopes of the linear model different for EE and LL?
# option 1: look at interaction term
anova(lm(plotData4$TotCopies ~ plotData4$VGLL3 + plotData4$VGLL3 * plotData4$DPH))

# calculate fitted values for means and 95% confidence interval
plotData4$fittedSimple = predict(fitSimple)
plotData4$fittedSimpleLwr = predict(fitSimple,interval="confidence")[,'lwr']
plotData4$fittedSimpleUpr = predict(fitSimple,interval="confidence")[,'upr']

# confidence interval for EE-LL contrast
marginal= emmeans(fitSimple, ~ VGLL3, weights="proportional")
pairs(marginal)

# contrast p-value
my_comparisons1  = data.frame(start=c("EE"), 
                            end=c("LL"), 
                            y=c(145,145),
                            label=c("*"))

# plot EE & LL
expPlot = ggplot(plotData4, aes(DPH,TotCopies,colour=VGLL3,group=VGLL3)) +
  geom_jitter(size=4,alpha=0.7) +
  theme_classic() + 
    scale_colour_manual(values=c("blue","red")) + 
  geom_smooth(aes(DPH,fittedSimple),method="loess") +
  #geom_ribbon(aes(ymin=fittedSimpleLwr, ymax=fittedSimpleUpr), alpha=0.1, color=NA) +
  xlab("Age (days post-hatch)") + 
  ylab(expression(~italic(vgll3)~'copies / ng RNA')) +
  #geom_text(aes(label=Sample),position=position_jitter(width=0.2),hjust=0, vjust=0,size=3) + 
  ylim(0,250) + 
  theme(axis.text=element_text(size=14,color="black"),
        axis.title=element_text(size=14,color="black")
        ,plot.title = element_text(size=14, face="bold",hjust=0.5),
        legend.position = "none")

insert = ggplot(as.data.frame(marginal)) + 
  geom_pointrange(aes(y=emmean, ymin=lower.CL, ymax=upper.CL, x=VGLL3,color=VGLL3),size=1) + 
  scale_color_manual(values=c("blue","red")) +
  geom_signif(data=my_comparisons1,
              aes(xmin=start, xmax=end, annotations=label, y_position=y),
              textsize = 3, size=0.5, tip_length = 0.01,
              manual=TRUE) +
  ylab("EMM") + 
  xlab(expression(~italic(vgll3)~'')) +
  theme_classic() +
  theme(axis.text=element_text(size=14,color="black"),
        axis.title=element_text(size=14,color="black")
        ,plot.title = element_text(size=14, face="bold",hjust=0.5),
        legend.position = "none")

ggdraw() +
  draw_plot(expPlot) +
  draw_plot(insert, x = 0.17, y = .65, width = .2, height = .3)
```

Model female gonad EE and LL expression.
```{r}
# subset for female gonads (both immature and mature) and only EE & LL genotypes
plotData8 = plotData2[which(plotData2$TissueOrRep=="Gn" & plotData2$TissueOrGonadType %in% c("Ov") & plotData2$VGLL3 %in% c("EE","LL")),]

# fit linear regression - alternatively (use poly() to perform the polynomial regressions)
fitSimple <- lm(plotData8$TotCopies ~ plotData8$VGLL3 + plotData8$DPH)
Anova(fitSimple)

# significance of the model
summary(fitSimple)

# confidence interval for EE-LL contrast
marginal= emmeans(fitSimple, ~ VGLL3, weights="proportional")
pairs(marginal)

# calculate fitted values
plotData8$fittedSimple = predict(fitSimple)

# plot EE & LL
ggplot(plotData8, aes(DPH,TotCopies,color=VGLL3,group=VGLL3)) +
  geom_jitter(size=4,alpha=0.7) +
  theme_classic() + 
  scale_color_manual(values=c("blue","red")) + 
  geom_smooth(aes(DPH,fittedSimple),method="loess") +
  xlab("Age (days post-hatch)") + 
  ylab(expression(~italic(vgll3)~'copies / ng RNA')) +
  theme(axis.text=element_text(size=14,color="black"),
        axis.title=element_text(size=14,color="black"),
        plot.title = element_text(size=14, face="bold",hjust=0.5),
        legend.position = "none")
```

Plot **heart** expression for males and females. Compile plot data and fitted values from lm into one data frame
```{r}
# split by gonad type
plotData10 = plotData2[which(plotData2$TissueOrRep=="Hr" & plotData2$Gonads %in% c("s") & plotData2$VGLL3 %in% c("EE","LL")),]
plotData11 = plotData2[which(plotData2$TissueOrRep=="Hr" & plotData2$Gonads %in% c("m") & plotData2$VGLL3 %in% c("EE","LL")),]
plotData12 = plotData2[which(plotData2$TissueOrRep=="Hr" & plotData2$Gonads %in% c("xl") & plotData2$VGLL3 %in% c("EE","LL")),]

# predict fitted values
fitSimple1 <- lm(plotData10$TotCopies ~ plotData10$VGLL3 + plotData10$DPH)
Anova(fitSimple1)
plotData10$fittedSimple = predict(fitSimple1)

fitSimple2 <- lm(plotData11$TotCopies ~ plotData11$VGLL3 + plotData11$DPH)
Anova(fitSimple2)
plotData11$fittedSimple = predict(fitSimple2)

fitSimple3 <- lm(plotData12$TotCopies ~ plotData12$VGLL3 + plotData12$DPH)
Anova(fitSimple3)
plotData12$fittedSimple = predict(fitSimple3)

# combine into one data frame
plotDataHeart = rbind(plotData10,plotData11,plotData12)

# plot
tissueLabels = c("m"="Female","s"="Immature male","xl"="Mature male")
ggplot(plotDataHeart,aes(DPH,TotCopies,group=VGLL3, color=VGLL3)) +
 geom_point(aes(DPH, TotCopies),size=3,alpha=0.7) +
 facet_wrap(~Gonads,labeller=labeller(Gonads=tissueLabels)) +
 geom_line(aes(y=fittedSimple)) +
 theme_classic() +
 scale_colour_manual(values=c("blue","red")) +
 theme(axis.title.y=element_text(size=16),
       axis.text.y=element_text(size=16, color="black"),
       axis.title.x=element_text(size=16),
       axis.text.x=element_text(hjust=1, vjust=0.5, size=16, color="black", angle=90),
       strip.text.x=element_text(size=16),
       strip.background = element_blank()) +
 ylab(expression(~italic(vgll3)~'copies / ng RNA')) +
 xlab("Age (days post-hatch)")

# heart expression difference between immature and mature maless - means
t.test(
plotData2$TotCopies[which(plotData2$TissueOrRep=='Hr' & plotData2$Gonads == "xl")],
plotData2$TotCopies[which(plotData2$TissueOrRep=='Hr' & plotData2$Gonads == "s")])

```

Plot **adipose** expression for males and females.
```{r}
# split by gonad type
plotData13 = plotData2[which(plotData2$TissueOrRep=="Ad" & plotData2$Gonads %in% c("s") & plotData2$VGLL3 %in% c("EE","LL")),]
plotData14 = plotData2[which(plotData2$TissueOrRep=="Ad" & plotData2$Gonads %in% c("m") & plotData2$VGLL3 %in% c("EE","LL")),]
plotData15 = plotData2[which(plotData2$TissueOrRep=="Ad" & plotData2$Gonads %in% c("xl") & plotData2$VGLL3 %in% c("EE","LL")),]

# predict fitted values
fitSimple4 <- lm(plotData13$TotCopies ~ plotData13$VGLL3 + plotData13$DPH)
Anova(fitSimple4)
plotData13$fittedSimple = predict(fitSimple4)

fitSimple5 <- lm(plotData14$TotCopies ~ plotData14$VGLL3 + plotData14$DPH)
Anova(fitSimple5)
plotData14$fittedSimple = predict(fitSimple5)

fitSimple6 <- lm(plotData15$TotCopies ~ plotData15$VGLL3 + plotData15$DPH)
Anova(fitSimple6)
plotData15$fittedSimple = predict(fitSimple6)

# combine into one data frame
plotDataAdipose = rbind(plotData13,plotData14,plotData15)

# plot
tissueLabels = c("m"="Female","s"="Immature male","xl"="Mature male")
ggplot(plotDataAdipose,aes(DPH,TotCopies,group=VGLL3, color=VGLL3)) +
 geom_point(aes(DPH, TotCopies),size=3,alpha=0.7) +
 facet_wrap(~Gonads,labeller=labeller(Gonads=tissueLabels)) +
 geom_line(aes(y=fittedSimple)) +
 theme_classic() +
 scale_colour_manual(values=c("blue","red")) +
 theme(axis.title.y=element_text(size=16),
       axis.text.y=element_text(size=16, color="black"),
       axis.title.x=element_text(size=16),
       axis.text.x=element_text(hjust=1, vjust=0.5, size=16, color="black", angle=90),
       strip.text.x=element_text(size=16),
       strip.background = element_blank()) +
 ylab(expression(~italic(vgll3)~'copies / ng RNA')) +
 xlab("Age (days post-hatch)")
```



## Comparison of vgll3 exon expression

Data from averaaged replicates (N=2).
```{r}
colnames(scamp1Means)[2]="TotCopies"
colnames(vgll3Exon1Means)[2]="TotCopies"

scamp1Means$Assay = "scamp1"
vgll3Mis1Means$Assay = "vgll3Mis1"
vgll3Mis2Means$Assay = "vgll3Mis2"
vgll3Exon1Means$Assay = "vgll3Exon1"

head(scamp1Means)
head(vgll3Mis1Means)
head(vgll3Mis2Means)
head(vgll3Exon1Means)

combinedAssays = rbind(scamp1Means,vgll3Exon1Means,vgll3Mis1Means,vgll3Mis2Means)
head(combinedAssays)
```

Plot expression levels in different exons.
```{r}
assayLabels = c("vgll3Exon1"="Exon 1-2 boundary","vgll3Mis1"="Exon 2","vgll3Mis2"="Exon 3")

ggplot(combinedAssays[which(combinedAssays$Assay != "scamp1"),],aes(VGLL3,TotCopies,color=VGLL3))+
  #geom_boxplot(outlier.alpha = 0) + 
  geom_quasirandom(size=5,alpha=0.3) +
  stat_summary(fun=mean, geom="point", shape=17, size=4, color="black", fill="black",alpha=0.5) +
  scale_color_manual(values=c("blue","red")) + 
  facet_wrap(~Assay, labeller=labeller(Assay=assayLabels)) +
  theme_classic() +
  theme(axis.text.x=element_text(size=14,color="black"),
       axis.title.x=element_blank(),
       axis.ticks.x=element_blank(),
       axis.text.y=element_text(size=14,color="black"),
       axis.title.y=element_text(size=14,color="black"),
       axis.line.x = element_blank(),
       strip.background=element_blank(),
       strip.text=element_text(size=14,color="black"),
       legend.position="none") +
  xlab("Age (days post-hatch)") + 
  ylab("copies / ng RNA") 

t.test(combinedAssays$TotCopies[which(combinedAssays$Assay == "vgll3Exon1" & combinedAssays$VGLL3 == "EE")],
       combinedAssays$TotCopies[which(combinedAssays$Assay == "vgll3Exon1" & combinedAssays$VGLL3 == "LL")])

t.test(combinedAssays$TotCopies[which(combinedAssays$Assay == "vgll3Mis1" & combinedAssays$VGLL3 == "EE")],
       combinedAssays$TotCopies[which(combinedAssays$Assay == "vgll3Mis1" & combinedAssays$VGLL3 == "LL")])

t.test(combinedAssays$TotCopies[which(combinedAssays$Assay == "vgll3Mis2" & combinedAssays$VGLL3 == "EE")],
       combinedAssays$TotCopies[which(combinedAssays$Assay == "vgll3Mis2" & combinedAssays$VGLL3 == "LL")])


```


## ASE in vgll3 mis-sense SNP 1
**Expression levels of the alleles in heterozygotes.** The column FractionalAbundance gives Ch1 versus Ch2 percentage of counts relative to the sum of counts. This is within samples so we don't need to take into account the starting amount of RNA. If there is a cis-regulatory difference that makes E-allele lower expressed it should be seen as ASE towards the L-allele in heterozygotes (keeping in mind the possibly higher efficiency of the L-allele probe). The analysis takes into account the possibly higher efficiency of the L (C)-probe using 1-effFactor correction.

Channel 1 corresponds to FAM dye and the C-allele (late maturation).
Channel 2 corresponds to VIC dye and the T-allele (early maturation).

Below we analyze if Ch1 and Ch2 give the same number of copies in the control DNA. Any difference from 50:50 ratios is due to one probe being lower efficiency (giving lower number of counts). We are calculating a correction factor to account for differences in efficiency based on the mean of all e28 DNA controls we are including in the results table.
```{r}
controlResMis1 = ddPCRres[which(ddPCRres$Sample == "e28-DNA"),]

# for poisson.test --> based on rate ratio
LnegCountsCtr = controlResMis1[which(controlResMis1$TargetType == "Ch1Unknown"),"Negatives"]
EnegCountsCtr = controlResMis1[which(controlResMis1$TargetType == "Ch2Unknown"),"Negatives"]
LcountsCtr = controlResMis1[which(controlResMis1$TargetType == "Ch1Unknown"),"Positives"]
EcountsCtr = controlResMis1[which(controlResMis1$TargetType == "Ch2Unknown"),"Positives"]
effFactor = apply(cbind(LcountsCtr,EcountsCtr,LnegCountsCtr,EnegCountsCtr), 1, function(x) poisson.test(c(x[1],x[2]),c(x[3],x[4]))[["estimate"]])
effFactor
effFactor = mean(effFactor)
paste("Mean of all control reps:",effFactor)
```

Below we test the difference between L and E allele counts within each heterozygous samples using a comparison of two Poisson rates. There is no consensus in the literature which test to use. For example in RNAseq people typically use a simple binomial test, but there it's hard to know the "negative" rate required to model the process using a Poisson distribution. Now we have this - it's the number of negative droplets - so using Poisson modeling may be more appropriate (this is what the QuantaSoft software from Bio-Rad uses). For this we need the number of positive AND negative counts. See more about comparing Poisson rates here: https://cran.r-project.org/web/packages/rateratio.test/vignettes/rateratio.test.pdf We are using the value stored as effFactor (relative freq of channel 1 positives to channel 2 positives based on genomic DNA from a heterozygote) as a null hypothesis.
```{r}
LnegCountsHet = plotData2$Ch1Neg[which(plotData2$VGLL3 %in% c("LE","EL") & plotData2$TissueOrRep == 'Gn' )]
EnegCountsHet = plotData2$Ch2Neg[which(plotData2$VGLL3 %in% c("LE","EL") & plotData2$TissueOrRep == 'Gn')]
LcountsHet = plotData2$Ch1Pos[which(plotData2$VGLL3 %in% c("LE","EL") & plotData2$TissueOrRep == 'Gn' )]
EcountsHet = plotData2$Ch2Pos[which(plotData2$VGLL3 %in% c("LE","EL") & plotData2$TissueOrRep == 'Gn' )]

# Initiate table for results of Poisson test
poissonRes = as.data.frame(matrix(data=NA,ncol=10,nrow=length(LcountsHet),dimnames=list(c(plotData2$Sample[which(plotData2$VGLL3 %in% c("LE","EL") & plotData2$TissueOrRep == 'Gn' )]),c("Sample","vgll3","sex","TotCopies","Lcounts","LcountsNeg","Ecounts","EcountsNeg","poissonPvalue","Timepoint"))))
poissonRes$Sample = rownames(poissonRes)
poissonRes$vgll3 = plotData2$VGLL3[which(plotData2$VGLL3 %in% c("LE","EL") & plotData2$TissueOrRep == 'Gn' )]
poissonRes$SexPheno = plotData2$SexPheno[which(plotData2$VGLL3 %in% c("LE","EL") & plotData2$TissueOrRep == 'Gn' )]
poissonRes$Gonads = plotData2$Gonads[which(plotData2$VGLL3 %in% c("LE","EL") & plotData2$TissueOrRep == 'Gn' )]
poissonRes$TotCopies = plotData2$TotCopies[which(plotData2$VGLL3 %in% c("LE","EL") & plotData2$TissueOrRep == 'Gn')]
poissonRes$Lcounts = LcountsHet
poissonRes$Ecounts = EcountsHet
poissonRes$LcountsNeg = LnegCountsHet
poissonRes$EcountsNeg = EnegCountsHet
poissonRes$Timepoint = plotData2$Timepoint[which(plotData2$VGLL3 %in% c("LE","EL") & plotData2$TissueOrRep == 'Gn')]

# Remove na's (individuals that have been collected but that not yet been analyzed).
poissonRes = poissonRes[which(is.na(poissonRes$TotCopies) == F),]

# Test ASE using Poisson rate ratio
poissonRes$poissonPvalue = apply(poissonRes[,c("Lcounts","Ecounts","LcountsNeg","EcountsNeg")], 1, function(x) poisson.test(c(x[1],x[2]),c(x[3],x[4]), r=effFactor)[["p.value"]])
poissonRes$poissonRateRatio = apply(poissonRes[,c("Lcounts","Ecounts","LcountsNeg","EcountsNeg")], 1, function(x) poisson.test(c(x[1],x[2]),c(x[3],x[4]), r=effFactor)[["estimate"]])
poissonRes$rateRatioCI95low = apply(poissonRes[,c("Lcounts","Ecounts","LcountsNeg","EcountsNeg")], 1, function(x) poisson.test(c(x[1],x[2]),c(x[3],x[4]), r=effFactor)[["conf.int"]][1])
poissonRes$rateRatioCI95high = apply(poissonRes[,c("Lcounts","Ecounts","LcountsNeg","EcountsNeg")], 1, function(x) poisson.test(c(x[1],x[2]),c(x[3],x[4]), r=effFactor)[["conf.int"]][2])

# Transfer DPH from previous data frames
poissonRes$DPH = plotData2$DPH[match(poissonRes$Sample,plotData2$Sample)]
head(poissonRes)

# sort according to Poisson rate ratio
poissonRes$Sample = factor(poissonRes$Sample, levels = poissonRes$Sample[order(poissonRes$poissonRateRatio)])
```

Plot Poisson rate ratio of L allele expression versus E allele expression. Null hypothesis of 1:1 expression rate is calculated based on heterozygous DNA as template. Points are colored according to wether the 95% confidence interval of rate ratio includes zero.
```{r}
plotData9 = plotData2[which(plotData2$VGLL3 %in% c("EL","LE") & plotData2$CollectionDate %in% c("19022018","20032018","20042018","17052018") == F & plotData2$TissueOrGonadType == 'IT'),]
plotData9$poissonPvalue = poissonRes$poissonPvalue[match(plotData9$Sample,poissonRes$Sample)]
plotData9$poissonRateRatio = poissonRes$poissonRateRatio[match(plotData9$Sample,poissonRes$Sample)]
plotData9$poissonRateRatio = poissonRes$poissonRateRatio[match(plotData9$Sample,poissonRes$Sample)]
plotData9$rateRatioCI95low = poissonRes$rateRatioCI95low[match(plotData9$Sample,poissonRes$Sample)]
plotData9$rateRatioCI95high = poissonRes$rateRatioCI95high[match(plotData9$Sample,poissonRes$Sample)]

# order by poisson rate ratio
plotData9$Sample = factor(plotData9$Sample, levels = plotData9$Sample[order(plotData9$poissonRateRatio)])

g = ggplot(plotData9,aes(x=Sample,log(poissonRateRatio,base=2),ymin=log(rateRatioCI95low,base=2),ymax=log(rateRatioCI95high,base=2),color=poissonPvalue<0.05)) +
 geom_pointrange(alpha=0.7) +
 geom_hline(yintercept = log(effFactor,base=2), linetype=2) +
 theme_classic() +
 theme(axis.title.x=element_blank(),
       axis.text.x=element_blank(),
       axis.ticks.x=element_blank(),
       axis.line.x = element_blank(),
       axis.text.y=element_text(size=14,color="black"),
       axis.title.y=element_text(size=14,color="black")) +
 scale_color_manual(values=c("darkgrey","black")) +
 theme(legend.position = "none") +
 ylab("Log2 L-rate/E-rate") + 
 ylim(-1.1,1.1)
rateFig = ggMarginal(g, type = "density", fill="darkgrey", margins="y")
rateFig

ggplot(plotData9) +
         geom_quasirandom(aes(1, FractionalAbundance,color=poissonPvalue),dodge.width=0.1,groupOnX=T,size=5,alpha=0.7) +
         geom_segment(y=mean(plotData9$FractionalAbundance),yend=mean(plotData9$FractionalAbundance),x=0.8,xend=1.2,size=1.5,lineend="round") +
         geom_hline(yintercept=mean(controlResMis1$FractionalAbundance,na.rm=T),color="darkgrey",linetype="dotted") +
         theme_classic() +
         ylab("L/(L+E) allele abundance %") +
         theme(axis.title.x=element_blank(),
               axis.text.x=element_blank(),
               axis.ticks.x=element_blank(),
               axis.text.y=element_text(size=14,color="black"),
               axis.title.y=element_text(size=14,color="black"),
               axis.line.x = element_blank(),
               legend.position="bottom",
               legend.title = element_text(hjust=1,vjust=0.75,size=10,color="black")) +
          #guides(colour=guide_legend(title.position = "top")) + 
          scale_color_gradient2(name = "ASE P-value", high="green", low="black", mid="green4", midpoint=0.4)
```

Can ASE account for EE-LL difference in exon 2 assay? Compare number of expressed molecules from each allele.
```{r}
# average molecule difference between LL and EE
aveHomozDiffExon2 = mean(combinedAssays$TotCopies[which(combinedAssays$Assay == "vgll3Mis1" & combinedAssays$VGLL3 == "EE")],na.rm=T) -
       mean(combinedAssays$TotCopies[which(combinedAssays$Assay == "vgll3Mis1" & combinedAssays$VGLL3 == "LL")],na.rm=T)

# average molecule difference between L and E
aveAlleleDiffExon2 = mean(plotData9$Ch2copies[which(plotData9$Sample %in% c("e28-DNA","G-gBlock","C-gBlock") == F)] ,na.rm=T) - mean(plotData9$Ch1copies[which(plotData9$Sample %in% c("e28-DNA","G-gBlock","C-gBlock") == F)],na.rm=T)

# how much % does 2*average L-E difference account from LL-EE difference (in molecules)?
2*(aveAlleleDiffExon2) / (aveHomozDiffExon2) * 100
```

## ASE in vgll3 mis-sense SNP 2

```{r}
colnames(ddPCRresMis2) = iconv(colnames(ddPCRresMis2), to = "ASCII", sub = "u")
colnames(ddPCRresMis2) = c(colnames(ddPCRresMis2)[-2],"Empty")
colnames(ddPCRresMis2)[1] = "Well"

# assign falsely labeled well
#ddPCRres$Sample[which(ddPCRres$Well == "H04")] = "e28-DNA"
ddPCRresMis2$Sample[which(ddPCRresMis2$Well == "D03")] = "437-Gn"
ddPCRresMis2$Sample[which(ddPCRresMis2$Well == "D09")] = "437-Gn"

# remove failed wells
ddPCRresMis2 = ddPCRresMis2[which(ddPCRresMis2$AcceptedDroplets > 10000),]

# remove samples that could not be re-analyzed from non-diluted RNA samples
ddPCRresMis2 = ddPCRresMis2[which(ddPCRresMis2$Sample %in% c("289-Gn","322-Gn","326-Gn","337-Gn","338-Gn","358-Gn","370-Gn","373-Gn","385-Gn","389-Gn","390-Gn") == F),]
```

Additional control wells.
```{r}
colnames(ddPCRresAdd) = iconv(colnames(ddPCRresAdd), to = "ASCII", sub = "u")
colnames(ddPCRresAdd) = c(colnames(ddPCRresAdd)[-2],"Empty")
colnames(ddPCRresAdd)[1] = "Well"

# remove failed wells
ddPCRresAdd = ddPCRresAdd[which(ddPCRresAdd$AcceptedDroplets > 10000),]
```


```{r}
controlResMis2 = rbind(ddPCRresMis2[which(ddPCRresMis2$Sample == "e28-DNA"),],ddPCRresAdd[which(ddPCRresAdd$Sample == "e28-DNA"),])

# for poisson.test --> based on rate ratio
LnegCountsCtr = controlResMis2[which(controlResMis2$TargetType == "Ch1Unknown"),"Negatives"]
EnegCountsCtr = controlResMis2[which(controlResMis2$TargetType == "Ch2Unknown"),"Negatives"]
LcountsCtr = controlResMis2[which(controlResMis2$TargetType == "Ch1Unknown"),"Positives"]
EcountsCtr = controlResMis2[which(controlResMis2$TargetType == "Ch2Unknown"),"Positives"]
effFactor = apply(cbind(LcountsCtr,EcountsCtr,LnegCountsCtr,EnegCountsCtr), 1, function(x) poisson.test(c(x[1],x[2]),c(x[3],x[4]))[["estimate"]])
effFactor
effFactor = mean(effFactor)
paste("Mean of all control reps:",effFactor)
```

```{r}
plotData16 = ddPCRresMis2

LnegCountsHet = plotData16$Negatives[which(plotData16$TargetType == "Ch1Unknown")]
EnegCountsHet = plotData16$Negatives[which(plotData16$TargetType == "Ch2Unknown")]
LcountsHet = plotData16$Positives[which(plotData16$TargetType == "Ch1Unknown")]
EcountsHet = plotData16$Positives[which(plotData16$TargetType == "Ch2Unknown")]

# Initiate table for results of Poisson test
poissonRes = as.data.frame(matrix(data=NA,ncol=6,nrow=length(LcountsHet),
                                  dimnames=list(plotData16$Positives[which(plotData16$TargetType == "Ch1Unknown")],c("Sample","Lcounts","LcountsNeg","Ecounts","EcountsNeg","poissonPvalue"))))
poissonRes$Sample = plotData16$Sample[which(plotData16$TargetType == "Ch1Unknown")]
poissonRes$Lcounts = LcountsHet
poissonRes$Ecounts = EcountsHet
poissonRes$LcountsNeg = LnegCountsHet
poissonRes$EcountsNeg = EnegCountsHet

# Test ASE using Poisson rate ratio
poissonRes$poissonPvalue = apply(poissonRes[,c("Lcounts","Ecounts","LcountsNeg","EcountsNeg")], 1, function(x) poisson.test(c(x[1],x[2]),c(x[3],x[4]), r=effFactor)[["p.value"]])
poissonRes$poissonRateRatio = apply(poissonRes[,c("Lcounts","Ecounts","LcountsNeg","EcountsNeg")], 1, function(x) poisson.test(c(x[1],x[2]),c(x[3],x[4]), r=effFactor)[["estimate"]])
poissonRes$rateRatioCI95low = apply(poissonRes[,c("Lcounts","Ecounts","LcountsNeg","EcountsNeg")], 1, function(x) poisson.test(c(x[1],x[2]),c(x[3],x[4]), r=effFactor)[["conf.int"]][1])
poissonRes$rateRatioCI95high = apply(poissonRes[,c("Lcounts","Ecounts","LcountsNeg","EcountsNeg")], 1, function(x) poisson.test(c(x[1],x[2]),c(x[3],x[4]), r=effFactor)[["conf.int"]][2])
```


```{r}
plotData17 = plotData16[which(plotData16$TargetType == "Ch1Unknown" & plotData16$Sample %in% c("e28-DNA","G-gBlock","C-gBlock") == F & is.na(plotData16$FractionalAbundance) == F),c("Sample","FractionalAbundance")]

plotData17$poissonPvalue = poissonRes$poissonPvalue[match(plotData17$Sample,poissonRes$Sample)]
plotData17$poissonRateRatio = poissonRes$poissonRateRatio[match(plotData17$Sample,poissonRes$Sample)]
plotData17$rateRatioCI95low = poissonRes$rateRatioCI95low[match(plotData17$Sample,poissonRes$Sample)]
plotData17$rateRatioCI95high = poissonRes$rateRatioCI95high[match(plotData17$Sample,poissonRes$Sample)]

ggplot(plotData17) +
         geom_quasirandom(aes(1, FractionalAbundance,color=poissonPvalue),dodge.width=0.1,groupOnX=T,size=5,alpha=0.7) +
         geom_segment(y=mean(plotData17$FractionalAbundance),yend=mean(plotData17$FractionalAbundance),x=0.8,xend=1.2,size=1.5,lineend="round") +
         geom_hline(yintercept=mean(controlResMis2$FractionalAbundance,na.rm=T),color="darkgrey",linetype="dotted") +
         theme_classic() +
         ylab("L/(L+E) allele abundance %") +
         theme(axis.title.x=element_blank(),
               axis.text.x=element_blank(),
               axis.ticks.x=element_blank(),
               axis.text.y=element_text(size=14,color="black"),
               axis.title.y=element_text(size=14,color="black"),
               axis.line.x = element_blank(),
               legend.position="bottom",
               legend.title = element_text(hjust=1,vjust=0.75,size=10,color="black")) +
          #guides(colour=guide_legend(title.position = "top")) + 
          scale_color_gradient2(name = "ASE P-value", high="green", low="black", mid="green4", midpoint=0.4)

nrow(plotData17)
nrow(plotData17[which(plotData17$poissonPvalue<0.05 & plotData17$FractionalAbundance>49.4),])
nrow(plotData17[which(plotData17$poissonPvalue<0.05 & plotData17$FractionalAbundance<49.4),])
mean(plotData17$FractionalAbundance)
```


Attact concentrations fo mis-sense SNP 2.
```{r}
# add starting total RNA concentration
ddPCRresMis2$RNAconc = RNAconcMis2$Concentration[match(ddPCRresMis2$Sample,RNAconcMis2$Sample)]
```


Normalizing to input RNA amount. QuantaSoft also provides calculated number of template molecules per 20 ul well. Each 20 ul well contains 2 ul of template at concentration X ng/ul_H2O, meaning 2•X ng/ul_H2O total RNA per well. Copies per ng total input RNA is then simply $\frac{copies per well}{2*RNAconc}$
```{r}
# Split ddPCR results into two based on channels
ddPCRch1 = ddPCRresMis2[which(ddPCRresMis2$TargetType == "Ch1Unknown"),]
ddPCRch2 = ddPCRresMis2[which(ddPCRresMis2$TargetType == "Ch2Unknown"),]

# calculation
ddPCRch1$Copies = as.numeric(ddPCRch1$CopiesPer20uLWell)/(ddPCRch1$RNAconc*2)
ddPCRch2$Copies = as.numeric(ddPCRch2$CopiesPer20uLWell)/(ddPCRch2$RNAconc*2)

# attach back to main data frame
plotData17$Ch1copies = ddPCRch1$Copies[match(plotData17$Sample,ddPCRch1$Sample)]
plotData17$Ch2copies = ddPCRch2$Copies[match(plotData17$Sample,ddPCRch2$Sample)]
plotData17$TotCopies = plotData17$Ch1copies + plotData17$Ch2copies 
head(plotData17)
```


Can ASE account for EE-LL difference? Compare number of expressed molecules from each allele.
```{r}
# average molecule difference between LL and EE
aveHomozDiff = mean(combinedAssays$TotCopies[which(combinedAssays$Assay == "vgll3Mis2" & combinedAssays$VGLL3 == "EE")],na.rm=T) -
       mean(combinedAssays$TotCopies[which(combinedAssays$Assay == "vgll3Mis2" & combinedAssays$VGLL3 == "LL")],na.rm=T)

# average molecule difference between L and E
aveAlleleDiff = mean(plotData17$Ch2copies[which(plotData17$Sample %in% c("e28-DNA","G-gBlock","C-gBlock") == F)] ,na.rm=T) - mean(plotData17$Ch1copies[which(plotData17$Sample %in% c("e28-DNA","G-gBlock","C-gBlock") == F)],na.rm=T)

# how much % does 2*average L-E difference account from LL-EE difference (in molecules)?
2*(aveAlleleDiff) / (aveHomozDiff) * 100
```



## ASE in alternative 5' UTR SNP
```{r}
colnames(ddPCRresUTR) = iconv(colnames(ddPCRresUTR), to = "ASCII", sub = "u")
colnames(ddPCRresUTR) = c(colnames(ddPCRresUTR)[-2],"Empty")
colnames(ddPCRresUTR)[1] = "Well"

# assign falsely labeled well
ddPCRresUTR$Sample[which(ddPCRresUTR$Well == "D03")] = "437-Gn"
ddPCRresUTR$Sample[which(ddPCRresUTR$Well == "D09")] = "437-Gn"

# remove failed wells
ddPCRresUTR = ddPCRresUTR[which(ddPCRresUTR$AcceptedDroplets > 10000),]

# remove samples that could not be re-analyzed from non-diluted RNA samples
ddPCRresUTR = ddPCRresUTR[which(ddPCRresUTR$Sample %in% c("289-Gn","322-Gn","326-Gn","337-Gn","338-Gn","358-Gn","370-Gn","373-Gn","385-Gn","389-Gn","390-Gn") == F),]
```

```{r}
plotData18 = ddPCRresUTR

LnegCountsHet = plotData18$Negatives[which(plotData18$TargetType == "Ch1Unknown")]
EnegCountsHet = plotData18$Negatives[which(plotData18$TargetType == "Ch2Unknown")]
LcountsHet = plotData18$Positives[which(plotData18$TargetType == "Ch1Unknown")]
EcountsHet = plotData18$Positives[which(plotData18$TargetType == "Ch2Unknown")]

# Initiate table for results of Poisson test
poissonResUTR = as.data.frame(matrix(data=NA,ncol=6,nrow=length(LcountsHet),
                                  dimnames=list(plotData18$Positives[which(plotData18$TargetType == "Ch1Unknown")],c("Sample","Lcounts","LcountsNeg","Ecounts","EcountsNeg","poissonPvalue"))))
poissonResUTR$Sample = plotData18$Sample[which(plotData18$TargetType == "Ch1Unknown")]
poissonResUTR$Lcounts = LcountsHet
poissonResUTR$Ecounts = EcountsHet
poissonResUTR$LcountsNeg = LnegCountsHet
poissonResUTR$EcountsNeg = EnegCountsHet

# Test ASE using Poisson rate ratio
poissonResUTR$poissonPvalue = apply(poissonResUTR[,c("Lcounts","Ecounts","LcountsNeg","EcountsNeg")], 1, function(x) poisson.test(c(x[1],x[2]),c(x[3],x[4]))[["p.value"]])
poissonResUTR$poissonRateRatio = apply(poissonResUTR[,c("Lcounts","Ecounts","LcountsNeg","EcountsNeg")], 1, function(x) poisson.test(c(x[1],x[2]),c(x[3],x[4]))[["estimate"]])
poissonResUTR$rateRatioCI95low = apply(poissonResUTR[,c("Lcounts","Ecounts","LcountsNeg","EcountsNeg")], 1, function(x) poisson.test(c(x[1],x[2]),c(x[3],x[4]))[["conf.int"]][1])
poissonResUTR$rateRatioCI95high = apply(poissonResUTR[,c("Lcounts","Ecounts","LcountsNeg","EcountsNeg")], 1, function(x) poisson.test(c(x[1],x[2]),c(x[3],x[4]))[["conf.int"]][2])
```

```{r}
plotData19 = plotData18[which(plotData18$TargetType == "Ch1Unknown",plotData18$Sample %in% c("EL-gBlock","NTC") == F & is.na(plotData18$InverseFractionalAbundance) == F),c("Sample","InverseFractionalAbundance")]

plotData19$poissonPvalue = poissonResUTR$poissonPvalue[match(plotData19$Sample,poissonResUTR$Sample)]
plotData19$poissonRateRatio = poissonResUTR$poissonRateRatio[match(plotData19$Sample,poissonResUTR$Sample)]
plotData19$rateRatioCI95low = poissonResUTR$rateRatioCI95low[match(plotData19$Sample,poissonResUTR$Sample)]
plotData19$rateRatioCI95high = poissonResUTR$rateRatioCI95high[match(plotData19$Sample,poissonResUTR$Sample)]

ggplot(plotData19) +
         geom_quasirandom(aes(1, InverseFractionalAbundance,color=poissonPvalue),dodge.width=0.1,groupOnX=T,size=5,alpha=0.7) +
         geom_segment(y=mean(plotData19$InverseFractionalAbundance),yend=mean(plotData19$InverseFractionalAbundance),x=0.8,xend=1.2,size=1.5,lineend="round") +
         geom_hline(yintercept=50,color="darkgrey",linetype="dotted") +
         theme_classic() +
         ylab("L/(L+E) allele abundance %") +
         theme(axis.title.x=element_blank(),
               axis.text.x=element_blank(),
               axis.ticks.x=element_blank(),
               axis.text.y=element_text(size=14,color="black"),
               axis.title.y=element_text(size=14,color="black"),
               axis.line.x = element_blank(),
               legend.position="bottom",
               legend.title = element_text(hjust=1,vjust=0.75,size=10,color="black")) +
          #guides(colour=guide_legend(title.position = "top")) + 
          scale_color_gradient2(name = "ASE P-value", high="green", low="black", mid="green4", midpoint=0.4)

nrow(plotData19)
nrow(plotData19[which(plotData19$poissonPvalue<0.05 & plotData19$InverseFractionalAbundance>50),])
nrow(plotData19[which(plotData19$poissonPvalue<0.05 & plotData19$InverseFractionalAbundance<50),])
mean(plotData19$InverseFractionalAbundance)

```

Combine ddPCR ASE results from URT and mis-sense SNPs 1 and 2 to the same figure.
```{r}
plotData19$Assay = "UTR"
plotData19$FractionalAbundance = plotData19$InverseFractionalAbundance # dummy fractional abundance column
plotData17$Assay = "Mis2"
plotData9$Assay = "Mis1"

plotData20 = rbind(plotData19[,c("Sample","FractionalAbundance","poissonPvalue","poissonRateRatio","rateRatioCI95low","rateRatioCI95high","Assay")], plotData17[,c("Sample","FractionalAbundance","poissonPvalue","poissonRateRatio","rateRatioCI95low","rateRatioCI95high","Assay")], plotData9[,c("Sample","FractionalAbundance","poissonPvalue","poissonRateRatio","rateRatioCI95low","rateRatioCI95high","Assay")])

plotDataControl = data.frame(Assay = c("UTR","Mis1", "Mis2"), mean = c(50,mean(controlResMis1$FractionalAbundance,na.rm=T),mean(controlResMis2$FractionalAbundance,na.rm=T)))
plotDataControl$Assay = factor(plotDataControl$Assay,levels=c("UTR","Mis1","Mis2"))

assayLabels = c("UTR"="5' UTR", "Mis1"="Exon 2","Mis2"="Exon 3")

plotData20$Assay = factor(plotData20$Assay,levels=c("UTR","Mis1","Mis2"))

plotData20$ColorPvalue=1-plotData20$poissonPvalue
plotData20$ColorPvalue[which(plotData20$Assay == "UTR" & plotData20$FractionalAbundance<50)] = -(plotData20$ColorPvalue[which(plotData20$Assay == "UTR" & plotData20$FractionalAbundance<50)])
plotData20$ColorPvalue[which(plotData20$Assay == "Mis1" & plotData20$FractionalAbundance<50.44)] = -(plotData20$ColorPvalue[which(plotData20$Assay == "Mis1" & plotData20$FractionalAbundance<50.44)])
plotData20$ColorPvalue[which(plotData20$Assay == "Mis2" & plotData20$FractionalAbundance<49.4)] = -(plotData20$ColorPvalue[which(plotData20$Assay == "Mis2" & plotData20$FractionalAbundance<49.4)])

ggplot(plotData20) +
         geom_quasirandom(aes(1, FractionalAbundance,color=ColorPvalue),dodge.width=0.1,groupOnX=T,size=5,alpha=0.7) +
         #geom_segment(y=mean(plotData4$FractionalAbundance),yend=mean(plotData4$FractionalAbundance),x=0.8,xend=1.2,size=1.5,lineend="round") +
         geom_hline(data=plotDataControl,aes(yintercept=mean),color="darkgrey",linetype="dotted") +
         facet_wrap(~Assay,labeller=labeller(Assay=assayLabels)) +
         theme_classic() +
         ylab("L/(L+E) allele abundance %") +
         theme(axis.title.x=element_blank(),
               axis.text.x=element_blank(),
               axis.ticks.x=element_blank(),
               axis.text.y=element_text(size=14,color="black"),
               axis.title.y=element_text(size=14,color="black"),
               axis.line.x = element_blank(),
               strip.background=element_blank(),
               strip.text=element_text(size=14,color="black"),
               legend.position="bottom",
               legend.title = element_text(hjust=1,vjust=0.75,size=10,color="black")) +
          #guides(colour=guide_legend(title.position = "top")) + 
          scale_color_gradient2(name = "ASE P-value", high="red", low="blue", mid="grey", midpoint=0,labels=c(0,0.5,1,0.5,0))

```








