## ----set, echo = FALSE--------------------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
)

## ----setup, message = FALSE, warning = FALSE, eval = FALSE--------------------
#  if (!require("BiocManager", quietly = TRUE))
#      install.packages("BiocManager")
#  BiocManager::install("MICSQTL")

## ----lib, message = FALSE, warning = FALSE------------------------------------
library(MICSQTL)

## ----lib2, message = FALSE, warning = FALSE, eval = FALSE---------------------
#  library(reshape2)
#  library(GGally)
#  library(ggplot2)

## ----obj----------------------------------------------------------------------
data(se)

## ----eg, eval = FALSE---------------------------------------------------------
#  se <- SummarizedExperiment(
#      assays = list(protein = your_protein_data),
#      rowData = your_anno_protein
#  )
#  metadata(se) <- list(
#      gene_data = your_gene_data
#  )

## ----cross, message = FALSE, warning = FALSE, results = FALSE, eval = FALSE----
#  se <- ajive_decomp(se, use_marker = FALSE, refactor_loading = TRUE)
#  se <- deconv(se, source = "cross", method = "JNMF",
#               Step = c(10^(-9), 10^(-7)),
#               use_refactor = 1000,
#               pinit = se@metadata$prop_gene,
#               ref_pnl = se@metadata$ref_gene)

## ----plot2, echo = FALSE, warning = FALSE, eval = FALSE-----------------------
#  ggplot(
#      cbind(data.frame(melt(metadata(se)$prop), metadata(se)$meta)),
#      aes(x = Var2, y = value, fill = Var2)
#  ) +
#      geom_point(
#          position = position_jitterdodge(
#              jitter.width = 0.1,
#              dodge.width = 0.7
#          ),
#          aes(fill = Var2, color = Var2),
#          pch = 21, alpha = 0.5
#      ) +
#      geom_boxplot(lwd = 0.7, outlier.shape = NA) +
#      theme_classic() +
#      facet_wrap(~disease) +
#      xlab("Cell type") +
#      ylab("Estimated proportion") +
#      theme(legend.position = "none")

## ----ajive, eval = FALSE------------------------------------------------------
#  se <- ajive_decomp(se, plot = TRUE,
#                     group_var = "disease",
#                     scatter = TRUE, scatter_x = "cns_1", scatter_y = "cns_2")
#  metadata(se)$cns_plot

## ----pca----------------------------------------------------------------------
pca_res <- prcomp(t(assay(se)), rank. = 3, scale. = FALSE)
pca_res_protein <- data.frame(pca_res[["x"]])
pca_res_protein <- cbind(pca_res_protein, metadata(se)$meta$disease)
colnames(pca_res_protein)[4] <- "disease"

## ----pcaplot, eval = FALSE----------------------------------------------------
#  ggpairs(pca_res_protein,
#          columns = seq_len(3), aes(color = disease, alpha = 0.5),
#          upper = list(continuous = "points")
#  ) + theme_classic()
#  
#  
#  pca_res <- prcomp(t(metadata(se)$gene_data), rank. = 3, scale. = FALSE)
#  pca_res_gene <- data.frame(pca_res[["x"]])
#  pca_res_gene <- cbind(pca_res_gene, metadata(se)$meta$disease)
#  colnames(pca_res_gene)[4] <- "disease"
#  ggpairs(pca_res_gene,
#          columns = seq_len(3), aes(color = disease, alpha = 0.5),
#          upper = list(continuous = "points")
#  ) + theme_classic()

## ----filter1------------------------------------------------------------------
head(rowData(se))

## ----filter2------------------------------------------------------------------
head(metadata(se)$anno_SNP)

## ----filter3------------------------------------------------------------------
target_protein <- rowData(se)[rowData(se)$Chr == 9, ][seq_len(3), "Symbol"]

## ----filter 4-----------------------------------------------------------------
se <- feature_filter(se,
    target_protein = target_protein,
    filter_method = c("allele", "distance"),
    filter_allele = 0.15,
    filter_geno = 0.05,
    ref_position = "TSS"
)

## ----filter5------------------------------------------------------------------
unlist(lapply(metadata(se)$choose_SNP_list, length))

## ----csQTL1, eval = FALSE-----------------------------------------------------
#  system.time(se <- csQTL(se))

## ----csQTL2, eval = FALSE-----------------------------------------------------
#  res <- metadata(se)$TOAST_output[[2]]
#  head(res[order(apply(res, 1, min)), ])

## ----sessionInfo, echo=FALSE--------------------------------------------------
sessionInfo()