---
title: "R Notebook"
output: html_notebook
---
```{r}
library(forestplot)
library(readr)

data <- read_delim("~/Documents/Roberts Lab/geoduck transcriptome/Geoduck-transcriptome_v3_bigtable_v08.tab.txt", "\t", escape_double = FALSE, trim_ws = TRUE)

GOslim_bin <- read_delim("~/Documents/Roberts Lab/geoduck transcriptome/GOslim-bin.txt", 
    "\t", escape_double = FALSE, trim_ws = TRUE)
```


```{r}

data_merged <- merge(data, GOslim_bin, by.x = "ContigID", by.y = "Column1", all = FALSE)
```

```{r}
#install.packages("mixtools")
#install.packages("mclust")
library("mixtools")
library("mclust")

data_cpg <- data_merged[,4]
data_cpg <- data_cpg[data_cpg >= 0.001 & data_cpg <= 1.5]

model <- normalmixEM(data_cpg, k = 2)

plot(model, which = 2, col2 = c("blue", "red"), xlab2 = " ", main2 = "Panopea generosa", font.main = 3)

lines(density(data_cpg), lty=2, lwd=2)


model.mclust <- Mclust(data_cpg)
summary(model.mclust)

model.mclust.1 <- Mclust(data_cpg, G = 1)
unimodal <- logLik(model.mclust.1)
bimodal <- logLik(model.mclust)
difference <- bimodal[1] - unimodal[1]

p.val <- 1- pchisq(difference, df = 3)
data_merged$GOSlim_bin <- as.factor(data_merged$GOSlim_bin)
GOSlim_types <- levels(data_merged$GOSlim_bin)

```

```{r}

for(i in 1:length(GOSlim_types))   {
  hist(data_cpg[data_merged$GOSlim_bin == GOSlim_types[i]], main = GOSlim_types[i], xlab = "CpG O/E", xlim = range(0.01, 1.5), breaks = 15)
}


```

