## ----options, include = FALSE-------------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    warning = FALSE,
    message = FALSE,
    fig.width = 10,
    fig.height = 10,
    out.width = "100%"
)

## ----installation, eval=FALSE-------------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE)) {
#     install.packages("BiocManager")
# }
# 
# BiocManager::install("TraianProt")

## ----library, eval=TRUE-------------------------------------------------------
library(TraianProt)
library(dplyr)
library(knitr)

## ----run-app, eval=FALSE------------------------------------------------------
# # Application structure
# app <- runTraianProt()
# 
# # Execution whether it is an interactive session
# if (interactive()) {
#     shiny::runApp(app)
# }

## ----parameters, eval=TRUE----------------------------------------------------
# 1. Define file paths for example data included in the package
file_path_raw <- system.file("extdata", "proteinGroups.txt", package = "TraianProt")
file_path_metadata <- system.file("extdata", "metadata_MaxQuant.tsv", package = "TraianProt")
# file_path_report #In case we are using DIANN data for reporting where the report.tsv file is located.
output_directory <- tempdir() # Use a temporary directory for outputs in this vignette

# 2. Define Experimental Parameters
platform_id <- 1 # 1=MaxQuant, 2=MSFragger, 3=DIA-NN, 4=PD, 5=ProteoScape
organism_id <- 1 # 1= Candida albicans model organism, 2= Other organism (e.g Homo sapiens/Mus musculus etc)
selected_conditions <- c("WT_H2O2", "WT") # Condition id defined in group column from metadata file.
conditions <- "WT_H2O2 vs WT" # Title for plots.
id_target <- "ENSG" # Gene ID type
organismo <- "calbicans" # Organism for annotation, more at https://biit.cs.ut.ee/gprofiler/page/organism-list
threshold <- 0.05 # Significance threshold for enrichment analyis.
font_size <- 12 # Term size for plots in the enrichment analysis.
terms_number <- 12 # Number of terms for enrichment plots in  the enrichment analysis.
psms <- TRUE # Whether PSMs correction is applied in differential analysis.

## ----metadata, eval=TRUE------------------------------------------------------
metadata <- read.delim(file_path_metadata, header = TRUE, row.names = NULL, stringsAsFactors = FALSE)

if (platform_id == 1 | platform_id == 2 | platform_id == 5) {
    metadata <- metadata %>%
        mutate(
            raw_name = intensity_sample_name,
            intensity_sample_name = raw_name,
            log2_col = sub("Intensity", "LOG2", raw_name)
        ) %>%
        mutate(
            raw_name = intensity_sample_name,
            unique_peptides_col = sub("Intensity", "Unique.peptides", raw_name)
        ) %>%
        select(intensity_sample_name, group, sample_name, log2_col, unique_peptides_col)
} else if (platform_id == 3) {
    metadata <- metadata %>%
        mutate(
            raw_name = basename(intensity_sample_name),
            intensity_sample_name = raw_name,
            log2_col = sub("\\.(d|raw)$", ".LOG2", raw_name, ignore.case = TRUE),
            unique_peptides_col = paste0("Unique peptides ", sub("\\.(d|raw)$", "", raw_name, ignore.case = TRUE))
        ) %>%
        select(intensity_sample_name, group, sample_name, log2_col, unique_peptides_col)
} else if (platform_id == 4) {
    metadata <- metadata %>%
        mutate(
            raw_name = intensity_sample_name,
            intensity_sample_name = raw_name,
            log2_col = sub("Abundance:", "LOG2", raw_name)
        ) %>%
        mutate(
            raw_name = intensity_sample_name,
            unique_peptides_col = sub("Abundance:", "Unique.peptides", raw_name)
        ) %>%
        select(intensity_sample_name, group, sample_name, log2_col, unique_peptides_col)
}

filtered_metadata <- metadata[metadata$group %in% selected_conditions, ]

## ----dataset, eval=TRUE-------------------------------------------------------
if (platform_id == 1 | platform_id == 2 | platform_id == 5) {
    raw <- read.delim(file_path_raw, sep = "\t", stringsAsFactors = FALSE, colClasses = "character")
    df <- TraianProt::quick_filtering(raw, platform_id, organism_id, filtered_metadata, selected_conditions)
} else if (platform_id == 3) {
    raw <- read.delim(file_path_raw, sep = "\t", stringsAsFactors = FALSE, colClasses = "character", check.names = FALSE)
    df <- quick_filtering(raw, platform_id, organism_id, filtered_metadata, selected_conditions, file_path_report)
    df <- as.data.frame(df)
} else if (platform_id == 4) {
    raw <- as.data.frame(readxl::read_xlsx(file_path_raw))

    df <- quick_filtering(raw, platform_id, organism_id, filtered_metadata, selected_conditions)
}

# LOG 2 Intensity

LOG2.names <- obtain_LOG.names(df)

## ----unique-proteins, eval = TRUE, out.width="100%"---------------------------
# Unique Proteins
unique_proteins <- TraianProt::obtain_unique_proteins(df, filtered_metadata, selected_conditions)
unique_proteins_control <- as.data.frame(unique_proteins[1], check.names = FALSE)
unique_proteins_treatment <- as.data.frame(unique_proteins[2], check.names = FALSE)

## ----pre-processing, eval = TRUE, out.width="100%"----------------------------
# 5. Preprocessing
# ------------------------------------------------------------------------------

# 5.1) UNIQUE PEPTIDES FILTERING

df <- unique_peptides_filter(df, filtered_metadata, number = 1, min_fraction = 0.5)

# 5.2) Quantification filtering

df.F <- TraianProt::filter_valids(df, filtered_metadata, unique_proteins, min_prop = 0.5, at_least_one <- FALSE, labeltype = 1)


# 5.3) Normalization

df.F <- median_centering(df.F, LOG2.names)
df.F <- normalization_func(df.F, LOG2.names, method = "mean")


# 5.4) Imputation

# df.FNI <- impute_KNN_data(as.data.frame(df.F), LOG2.names, k = 5)
df.FNI <- impute_data(as.data.frame(df.F), LOG2.names)

df.FNI <- df.F # No imputation preferred action by developers.


total <- bind_rows(df.FNI, as.data.frame(unique_proteins[1], check.names = FALSE))
total_dataset <- bind_rows(total, as.data.frame(unique_proteins[2], check.names = FALSE))

## ----venn, eval=FALSE---------------------------------------------------------
# library(VennDiagram)
# venn_diagram_plot <- grid.draw(venn_diagram(df.F, unique_proteins, color1 = "blue", color2 = "maroon"))
# plot(venn_diagram_plot)

## ----quantified-prots, eval=TRUE----------------------------------------------
identify_proteins(df, filtered_metadata, platform_id, selected_conditions)

## ----boxplot, eval=TRUE-------------------------------------------------------
boxplot <- boxplot_function(df.FNI, filtered_metadata, selected_conditions)
print(boxplot)

## ----preimputation, eval=TRUE-------------------------------------------------
preimputation_state(df.F, filtered_metadata$log2_col)

## ----correlation-plot, eval=TRUE----------------------------------------------
corrplot_function(df.FNI[filtered_metadata$log2_col], filtered_metadata, display = "shade")

## ----pca, eval=TRUE-----------------------------------------------------------
pca(df.FNI, filtered_metadata, selected_conditions)

## ----tsne, eval=TRUE----------------------------------------------------------
tsne(df.FNI, filtered_metadata, perplexity_num = 2, selected_conditions)

## ----differential-analysis, eval = TRUE, out.width="100%"---------------------
limma_1 <- TraianProt::statistical_analysis(df.FNI, test = 2, paired = FALSE, filtered_metadata, logfc = 1, sig = 0.05, adjval = "fdr", statval = 1, unique_proteins, way = 2, psms = TRUE, platform_id, selected_conditions, diann_dir = if (platform_id == 3) file_path_report else NULL)

limma <- limma_1[-6]
limma <- merge(total_dataset, limma, by = "Protein", check.names = FALSE)

if (psms == TRUE) {
    limma <- limma %>%
        dplyr::select("Protein", "Protein_description", "logFC", "sca.P.Value", "sca.adj.pval", "expression", everything())
} else if (psms == FALSE) {
    limma <- limma %>%
        select("Protein", "Protein_description", "logFC", "p.value", "adj.P.Val", "expression", everything())
}

limma <- limma[order(limma$expression), ]
row.names(limma) <- limma$Protein

# write.xlsx(limma, file = "AD_.xlsx", rowNames = FALSE)

## ----volcano-plot, eval=TRUE--------------------------------------------------
volcano <- volcano_plot_tiff(limma, title = "Treatment vs Control", label = 3, statval = 2, psms)
print(volcano)

## ----heatmap, eval=FALSE------------------------------------------------------
# heatmap <- my_heatmap_differential(limma, df.FNI, filtered_metadata$log2_col, title = "Treatment vs Control")

## ----functional-analysis, eval = TRUE, out.width="100%"-----------------------
Go_terms <- Goterms_finder(limma, df, target = "ENSG", numeric_ns = "", mthreshold = Inf, filter_na = TRUE, organismo = "calbicans", custombg = FALSE, platform_id, user_threshold = 0.05, multi_query = FALSE, evcodes = TRUE, sources = c("GO", "KEGG", "WP", "REAC"))

# write.xlsx(Go_terms, file = "AF_.xlsx", rowNames = FALSE)

## ----dot-barplot,  eval = TRUE, echo=FALSE, out.width="100%"------------------
dotplot_func(Go_terms, x = "GeneRatio", title = conditions, split = "Conditions", font.size = 10, showCategory = 10, color = "adj.P.Val")

barplot_func(Go_terms, 10, conditions = conditions, font.size = 10)

## ----interaction-analysis, echo=FALSE, out.width="100%"-----------------------
interactions <- interactions_up(limma, taxonid = 237561, score = 400)

interactions_down <- interactions_down(limma, taxonid = 237561, score = 400)


graph_analysis <- igraph_analysis(interactions, taxonid = 237561, score = 400)

## ----session-info-------------------------------------------------------------
sessionInfo()