R Notebook

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)

Parsed with column specification:
cols(
  .default = col_character(),
  length = col_integer(),
  GC_content = col_double(),
  `CpG_o/e` = col_double(),
  male_unique_count = col_integer(),
  female_unique_count = col_integer(),
  Evalue = col_double(),
  `evalue-Gigaton` = col_double(),
  `evalue-Ruphibase` = col_double(),
  `evalue-Sigenae` = col_double(),
  Dheilly_Cluster = col_integer()
)
See spec(...) for full column specifications.

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

Parsed with column specification:
cols(
  Column1 = col_character(),
  GOSlim_bin = col_character()
)

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

#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)

number of iterations= 325 

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)

----------------------------------------------------
Gaussian finite mixture model fitted by EM algorithm 
----------------------------------------------------

Mclust V (univariate, unequal variance) model with 2 components:

 log.likelihood     n df      BIC       ICL
       2034.472 19442  5 4019.568 -154.6012

Clustering table:
    1     2 
14565  4877 

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)

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)
}