Comparing Igor model

comparing_igor_model.R

library(data.table)
library(readr)
library(magrittr)
library(dplyr)
library(plotly)
library(lambda.tools)

read_model_data = function(model_file) {
    # @params model_file the Igor model files 

    model_file_lines = readLines(model_file)

    stat_model = model_file_lines %>% grep("^%", ., value = T) %>% sub("%", "", .) %>% strsplit(",") %>% unlist %>% as.numeric

    events_name = model_file_lines %>% grep("@",., value = T)
    events_num = model_file_lines %>% grep("Dim",., value = T) %>% sub("\\$Dim\\[", "", .) %>% sub("\\]", "", .)  %>% strsplit(",") %>% sapply(function(x){fold(as.integer(x), function(a, b) a * b, 1)})
    events_v = rep(events_name, events_num)

    d = list(p_model = stat_model, event = events_v)

    d

}

comparing_two_model = function(model1, model2, model_names = c('model1', 'model2'), exclude_events = c(), fig = T) {
    # @params model1, model2: the output of read_model_data() function
    # @params model_names: the model names to display in the ouput figure
    # @params exclude_events: the events to exclude in the statistics and visualization
    # @params fig: if the visualization will be carried out or not

    if (any(model1$event != model2$event)) {
        warning("Two model are different, they can not be compared.")
    }

    d = data.table(p_model1 = -log10(model1$p_model + 0.000001), p_model2 = -log10(model2$p_model + 0.000001), event = model1$event)
    d = d[!(event %in% exclude_events)]

    if (fig) {
        g = ggplot(d) + aes(x = p_model1, y = p_model2, color = event)
        p = g + geom_point() + theme_bw() + labs(x = model_names[1], y = model_names[2])
        ggplotly(p) %>% print
    }
    l = lm(p_model1 ~ p_model2 + 0, data = d)
    cat("=================================================\n")
    cat("Spearman Correlation:\n")
    with(d, cor.test(p_model1, p_model2, method = "spearman")) %>% print

    cat("=================================================\n")
    cat("Pearson Correlation:\n")
    with(d, cor.test(p_model1, p_model2, method = "spearman")) %>% print

    cat("=================================================\n")
    cat("Linear regression:\n")
    summary(l) %>% print
}

## example:
m1 = read_model_data("../data/Anne-Marie_20200111/igor_sample_outputs_reads/meghan_bcmouse_igor_inference/final_marginals.txt")
m2 = read_model_data("../data/Anne-Marie_20200111/igor_sample_outputs_reads_Emilie/emilie_bcmouse_igor_inference/final_marginals.txt")

comparing_two_model(m1, m2, model_names = c("Meghan's data", "Emilie's data"))
comparing_two_model(m1, m2, model_names = c("Meghan's data", "Emilie's data"), exclude_events = c("@d_3_del"))