# this script runs mashr on prepped eQTL results to find response and cell-type-specific eQTLs
# 1/28/2020 JR
# based on mashr vignette https://stephenslab.github.io/mashr/articles/intro_mash.html
# last editted and ran 10/7/2020 JR: cleaned up

# prereq: ./mashr_Tcell_prep.R

library(ashr)
library(mashr)
library(data.table)
library(dplyr)
library(doParallel)
cores <- as.integer(Sys.getenv("SLURM_STEP_TASKS_PER_NODE"))
registerDoParallel(cores = cores)
library(RhpcBLASctl)
blas_set_num_threads(cores)

folder = "SCAIP1-6"

# load the mashr input data:
load("Tcell_mash-ready-data.Rd")

# Step 3: fit the model
# this fits a mixture model to the data, estimating the mixture proportions
m.ed = mash(data, U.ed) #ISSUE: not done after two weeks

#### END 10/7/2020 (never got past this step)
m   = mash(data, c(U.c,U.ed))
save(m,m.ed,file="m_m.ed_Tcell.Rd")
# Step 4: Extract Posterior Summaries:
# posterior mean:
head(get_lfsr(m.ed))
# posterior standard deviation:
head(get_pm(m.ed))
# local false sign rate (measure of significance)
head(get_psd(m.ed))


sessionInfo()