R/plot_is_bridgeable.R

In OlinkAnalyze: Facilitate Analysis of Proteomic Data from Olink

#' Plots for each bridgeable assays between two products.
#'
#' @author
#'   Amrita Kar
#'   Klev Diamanti
#'
#' Generates a combined plot per assay containing a violin and boxplot for IQR
#' ranges; correlation plot of NPX values; a median count bar plot and KS plots
#' from the 2 products.
#'
#' @param df A tibble containing the cross-product bridge normalized dataset
#' generated by \code{\link{olink_normalization}}.
#' @param check_log A named list returned by [`check_npx()`]. If `NULL`,
#' [`check_npx()`] will be run internally using `df`.
#' @param olink_id Character vector of Olink assay identifiers \var{OlinkID} for
#' which bridgeability plots will be created. If null, plots for all assays in
#' \var{data} will be created. (default = NULL)
#' @param median_counts_threshold Threshold indicating the minimum median counts
#' for each product (default = 150).
#' @param min_count Threshold indicating the minimum number of counts per
#' data point (default = 10). Data below \var{min_count} are excluded.
#'
#' @return An object of class "ggplot" containing 4 plots for each assay.
#'
#' @export
#'
#' @examples
#' \donttest{
#' if (rlang::is_installed(pkg = c("ggpubr"))) {
#'   npx_ht <- OlinkAnalyze:::data_ht_small |>
#'     dplyr::filter(
#'       .data[["SampleType"]] == "SAMPLE"
#'     )
#'
#'   npx_3072 <- OlinkAnalyze:::data_3k_small |>
#'     dplyr::filter(
#'       .data[["SampleType"]] == "SAMPLE"
#'     )
#'
#'   overlapping_samples <- intersect(
#'     x = npx_ht$SampleID,
#'     y = npx_3072$SampleID
#'   )
#'
#'   data_norm <- OlinkAnalyze::olink_normalization(
#'     df1 = npx_ht,
#'     df2 = npx_3072,
#'     overlapping_samples_df1 = overlapping_samples,
#'     df1_project_nr = "Explore HT",
#'     df2_project_nr = "Explore 3072",
#'     reference_project = "Explore HT",
#'     df1_check_log = check_npx(df = npx_ht) |>
#'       suppressMessages() |>
#'       suppressWarnings(),
#'     df2_check_log = check_npx(df = npx_3072) |>
#'       suppressMessages() |>
#'       suppressWarnings()
#'   )
#'
#'   data_norm_bridge_p <- OlinkAnalyze::olink_bridgeability_plot(
#'     df = data_norm,
#'     check_log = check_npx(df = data_norm) |>
#'       suppressMessages() |>
#'       suppressWarnings(),
#'     olink_id = c("OID40770", "OID40835"),
#'     median_counts_threshold = 150L,
#'     min_count = 10L
#'   )
#' }
#' }
#'
olink_bridgeability_plot <- function(df,
                                     check_log = NULL,
                                     olink_id = NULL,
                                     median_counts_threshold = 150L,
                                     min_count = 10L) {

  # Check if all required libraries for this function are installed
  rlang::check_installed(
    pkg = c("ggpubr"),
    call = rlang::caller_env()
  )

  # set seed
  set.seed(seed = 1L)

  # Check data format
  check_log <- run_check_npx(df = df, check_log = check_log)

  # check that check_log$col_names$count exists
  check_log_colname(check_log = check_log, col_key = "count")

  # check that columns not pre-cleared by check_log are present
  check_columns(
    df = df,
    col_list = list("Project", "BridgingRecommendation")
  )

  # check OlinkID ----

  olink_id_set <- df |>
    dplyr::pull(
      .data[[check_log$col_names$olink_id]]
    ) |>
    unique()

  if (is.null(olink_id)) {

    # if olink_id is NULL, then use all OlinkID
    olink_id <- olink_id_set

  } else if (!all(olink_id %in% olink_id_set)) {

    # check that all OlinkID are present in the df
    non_overlap_oid <- olink_id[!(olink_id %in% olink_id_set)] # nolint: object_usage_linter
    cli::cli_warn(
      c(
        "{cli::qty(non_overlap_oid)}{.val {length(non_overlap_oid)}} assay{?s}
        in {.arg olink_id} {?is/are} not present in the dataset {.arg df}:
        {.val {non_overlap_oid}}.",
        "i" = "Keeping only Olink assay identifiers present in the dataset:
        {.val {olink_id[olink_id %in% olink_id_set]}}."
      )
    )
    olink_id <- olink_id[olink_id %in% olink_id_set]
  }

  # keep only rows with OlinkIDs to be plotted
  df <- df |>
    dplyr::filter(
      .data[[check_log$col_names$olink_id]] %in% .env[["olink_id"]]
    )

  # Clean up df to assays with sufficient counts for plotting ----

  if (nrow(df) > 0L) {
    nrow_df <- nrow(df)
    df <- bridgeability_prep_data(
      df = df,
      check_log = check_log,
      min_count = min_count
    )
    if (nrow(df) != nrow_df) {
      cli::cli_inform(
        "{cli::qty(nrow_df - nrow(df))}Removed {.val {nrow_df - nrow(df)}}
        row{?s} with less than {.val {min_count}} counts from dataset
        {.arg df}!"
      )
    }
    rm(nrow_df)
  }

  # Check that all BridgingRecommendation are valid ----

  if (nrow(df) > 0L) {
    bridge_recommend_all <- unique(df[["BridgingRecommendation"]])

    if (!all(bridge_recommend_all %in% bridge_recommendations)) {
      invalid_recommendation <- setdiff(x = bridge_recommend_all, # nolint: object_usage_linter
                                        y = bridge_recommendations)

      cli::cli_abort(
        c(
          "x" = "{cli::qty(invalid_recommendation)}Identified invalid bridging
        recommendation{?s} in column {.arg {\"BridgingRecommendation\"}}:
        {.val {invalid_recommendation}}.",
          "i" = "Expected values are {.val {bridge_recommendations}}."
        )
      )
    }
  }

  # Check that each assay has a unique BridgingRecommendation ----

  if (nrow(df) > 0L) {
    df_br <- df |>
      dplyr::distinct(
        .data[[check_log$col_names$olink_id]],
        .data[["BridgingRecommendation"]]
      ) |>
      dplyr::count(
        .data[[check_log$col_names$olink_id]]
      ) |>
      dplyr::filter(
        .data[["n"]] > 1L
      ) |>
      dplyr::pull(
        .data[[check_log$col_names$olink_id]]
      )

    if (length(df_br) > 0L) {
      cli::cli_abort(
        c(
          "x" = "{cli::qty(df_br)}Identified {.val {length(df_br)}} assay{?s}
          with multiple bridging recommendations in column
          {.arg {\"BridgingRecommendation\"}}: {.val {df_br}}.",
          "i" = "Each assay should have a unique bridging recommendation!"
        )
      )
    }
  }

  # Drop BridgingRecommendation not relevant for plotting (NotOverlapping) ----

  accepted_br <- unname(bridge_recommendations[ # nolint: object_usage_linter
    !(names(bridge_recommendations) %in% c("not_overlapping"))
  ])
  non_accepted_br <- unname(bridge_recommendations[
    names(bridge_recommendations) %in% c("not_overlapping")
  ])

  if (nrow(df) > 0L &&
        any(unique(df[["BridgingRecommendation"]]) %in%
              bridge_recommendations[c("not_overlapping")])) {

    df_br_no_overlap <- df |> # nolint: object_usage_linter
      dplyr::filter(
        .data[["BridgingRecommendation"]] %in% .env[["non_accepted_br"]]
      ) |>
      dplyr::pull(
        .data[[check_log$col_names$olink_id]]
      ) |>
      unique()

    cli::cli_warn(
      c(
        "Identified {.val {length(df_br_no_overlap)}} assay{?s} with
        {.arg {\"BridgingRecommendation\"}} equal to {.val {non_accepted_br}}:
        {.val {df_br_no_overlap}}.",
        "i" = "Only assays with bridging recommendations {.val {accepted_br}}
        will be plotted!"
      )
    )

    df <- df |>
      dplyr::filter(
        !(.data[[check_log$col_names$olink_id]] %in% .env[["df_br_no_overlap"]])
      )
  }

  # Check there are exactly 2 projects for each assay ----

  if (nrow(df) > 0L) {
    df_proj <- df |>
      dplyr::distinct(
        .data[[check_log$col_names$olink_id]],
        .data[["Project"]]
      ) |>
      dplyr::count(
        .data[[check_log$col_names$olink_id]]
      ) |>
      dplyr::filter(
        .data[["n"]] != 2L
      ) |>
      dplyr::pull(
        .data[[check_log$col_names$olink_id]]
      )

    if (length(df_proj) > 0L) {
      cli::cli_abort(
        c(
          "x" = "{cli::qty(df_proj)}Identified {.val {length(df_proj)}}
          assay{?s} not belonging to exactly 2 projects: {.val {df_proj}}.",
          "i" = "Each assay should appear in exactly 2 projects in {.arg {df}}!"
        )
      )
    }
  }

  # Check that there are no duplicate samples ----

  if (nrow(df) > 0L) {
    df_dups <- df |>
      dplyr::count(
        .data[[check_log$col_names$sample_id]],
        .data[[check_log$col_names$olink_id]],
        .data[["Project"]]
      ) |>
      dplyr::filter(
        .data[["n"]] > 1L
      ) |>
      dplyr::pull(
        .data[[check_log$col_names$sample_id]]
      ) |>
      unique()

    if (length(df_dups) > 0L) {
      cli::cli_abort(
        c(
          "x" = "{cli::qty(df_dups)}Identified {.val {length(df_dups)}}
          duplicate sample{?s} in dataset {.arg df}: {.val {df_dups}}.",
          "i" = "There should be exactly one sample per combination of sample
          identifier, assay identifier and project!"
        )
      )
    }
  }

  # check if any rows are remaining ----

  if (nrow(df) > 0L) {
    olink_id <- olink_id[olink_id %in% df[[check_log$col_names$olink_id]]]
  } else {
    cli::cli_abort(
      c(
        "x" = "Dataset {.arg df} has {.val {0}} rows left!",
        "i" = "No plots can be generated!"
      )
    )
  }

  # Bridgeable plot ----

  out_plts <- lapply(
    olink_id,
    function(oid) {
      data_tmp <- df |>
        dplyr::filter(
          .data[[check_log$col_names$olink_id]] %in% .env[["oid"]]
        ) |>
        dplyr::mutate(
          oid_assay = paste(.data[[check_log$col_names$assay]],
                            .data[[check_log$col_names$olink_id]],
                            sep = " - ")
        )

      # unique bridging recommendation for assay
      bridge_suggest <- unique(data_tmp[["BridgingRecommendation"]])

      # iqr plot
      iqr_p <- bridgeability_iqr_range_plt(
        df = data_tmp,
        check_log = check_log
      )
      # r2 plot
      r2_p <- bridgeability_r2_plt(
        df = data_tmp,
        check_log = check_log
      )
      # counts plot
      counts_p <- bridgeability_counts_plt(
        df = data_tmp,
        median_counts_threshold = median_counts_threshold,
        check_log = check_log
      )
      # ks plot
      ks_p <- bridgeability_ks_plt(
        df = data_tmp,
        check_log = check_log
      )

      # combine plots

      out_plot <- bridgeability_combine_plots(
        iqr = iqr_p,
        r2 = r2_p,
        counts = counts_p,
        ks = ks_p,
        title = ifelse(
          bridge_suggest %in%
            bridge_recommendations[
              names(bridge_recommendations) %in%
                c("median_centering", "quantile_smoothing")
            ],
          paste0(unique(data_tmp[["oid_assay"]]), " (bridgeable)"),
          paste0(unique(data_tmp[["oid_assay"]]), " (non bridgeable)")
        )
      )

      return(out_plot)
    }
  )
  names(out_plts) <- olink_id

  return(out_plts)
}

bridgeability_prep_data <- function(df,
                                    check_log,
                                    min_count = 10L) {
  df <- df |>
    dplyr::filter(
      .data[[check_log$col_names$count]] > .env[["min_count"]]
    )
  return(df)
}

bridgeability_iqr_range_plt <- function(df,
                                        check_log) {

  iqr_range_plt <- ggplot2::ggplot(
    data = df,
    mapping = ggplot2::aes(
      x = .data[["oid_assay"]],
      y = .data[[check_log$col_names$quant]],
      fill = .data[["Project"]]
    )
  ) +
    ggplot2::geom_violin(
      alpha = 0.4,
      position = ggplot2::position_nudge(x = 0),
      width = 0.4
    ) +
    ggplot2::geom_boxplot(
      width = .1,
      outlier.shape = NA,
      alpha = 0.4
    ) +
    ggplot2::labs(
      x = "Assay",
      y = "NPX Distribution",
      fill = "Platform: "
    ) +
    ggplot2::guides(
      fill = ggplot2::guide_legend(
        nrow = 1L,
        byrow = TRUE
      )
    ) +
    OlinkAnalyze::set_plot_theme() +
    OlinkAnalyze::olink_fill_discrete() +
    OlinkAnalyze::olink_color_discrete() +
    ggplot2::theme(
      axis.text.x = ggplot2::element_blank()
    )

  return(iqr_range_plt)
}

bridgeability_r2_plt <- function(df,
                                 check_log) {

  projects <- unique(df[["Project"]])

  df <- df |>
    dplyr::select(
      dplyr::all_of(
        c(check_log$col_names$sample_id,
          "oid_assay",
          check_log$col_names$quant,
          "Project")
      )
    ) |>
    tidyr::pivot_wider(
      names_from = dplyr::all_of("Project"),
      values_from = dplyr::all_of(check_log$col_names$quant)
    ) |>
    tidyr::drop_na()

  r2_lm <- stats::cor(
    x = dplyr::pull(df, .data[[projects[1L]]]),
    y = dplyr::pull(df, .data[[projects[2L]]]),
    use = "everything",
    method = "pearson"
  ) |>
    (\(.) . ^ 2L)() |>
    signif(2L)

  caps <- ifelse(r2_lm < 0.2,
                 "Possibly bridging background to background",
                 "")

  r2_plt <- ggplot2::ggplot(
    data = df,
    mapping = ggplot2::aes(
      x = .data[[projects[1L]]],
      y = .data[[projects[2L]]]
    )
  ) +
    ggplot2::geom_point(
      color = "blue",
      alpha = 0.4
    ) +
    ggplot2::geom_smooth(
      method = "lm",
      formula = "y ~ x",
      color = "black"
    ) +
    ggpubr::stat_cor(
      method = "pearson",
      ggplot2::aes(
        label = ggplot2::after_stat(x = .data[["rr.label"]])
      ),
      geom = "label"
    ) +
    OlinkAnalyze::set_plot_theme() +
    OlinkAnalyze::olink_color_discrete() +
    ggplot2::labs(
      caption = caps
    )

  return(r2_plt)
}

bridgeability_counts_plt <- function(df,
                                     median_counts_threshold,
                                     check_log) {

  counts_plt <- df |>
    dplyr::group_by(
      dplyr::across(
        dplyr::all_of(
          c(check_log$col_names$olink_id, "Project")
        )
      )
    ) |>
    dplyr::mutate(
      median_count = stats::median(x = .data[[check_log$col_names$count]],
                                   na.rm = TRUE)
    ) |>
    dplyr::ungroup() |>
    dplyr::distinct(
      .data[["oid_assay"]], .data[["Project"]], .data[["median_count"]]
    ) |>
    ggplot2::ggplot(
      mapping = ggplot2::aes(
        x = .data[["oid_assay"]],
        fill = .data[["Project"]],
        y = .data[["median_count"]],
        label = .data[["median_count"]]
      )
    ) +
    ggplot2::geom_col(
      width = 0.5,
      position = "dodge"
    ) +
    ggplot2::geom_text(
      position = ggplot2::position_dodge(0.5),
      vjust = 0.5
    ) +
    ggplot2::geom_hline(
      yintercept = median_counts_threshold,
      color = "black",
      linetype = "dashed",
      linewidth = 0.7
    ) +
    ggplot2::labs(
      x = "Assay",
      y = "Median Count",
      fill = "Platform:"
    ) +
    OlinkAnalyze::set_plot_theme() +
    OlinkAnalyze::olink_fill_discrete() +
    ggplot2::theme(
      axis.text.x = ggplot2::element_blank()
    )

  return(counts_plt)
}

bridgeability_ks_plt <- function(df,
                                 check_log) {

  projects <- unique(df[["Project"]])

  data_wide <- df |>
    dplyr::select(
      dplyr::all_of(
        c(check_log$col_names$sample_id,
          "oid_assay",
          check_log$col_names$quant,
          "Project")
      )
    ) |>
    tidyr::pivot_wider(
      names_from = dplyr::all_of("Project"),
      values_from = dplyr::all_of(check_log$col_names$quant)
    ) |>
    tidyr::drop_na()

  # Calculate empirical cumulative distribution function (ECDF) per platform
  ecdf_p1 <- dplyr::pull(data_wide, .data[[projects[1L]]]) |> stats::ecdf()
  ecdf_p2 <- dplyr::pull(data_wide, .data[[projects[2L]]]) |> stats::ecdf()

  # Find min/max statistics to draw line between points of greatest distance
  min_max <- seq(
    from = dplyr::select(data_wide, dplyr::all_of(projects)) |> min(),
    to = dplyr::select(data_wide, dplyr::all_of(projects)) |> max(),
    length.out = nrow(data_wide)
  )
  x0 <- min_max[which(abs(ecdf_p1(min_max) - ecdf_p2(min_max)) ==
                        max(abs(ecdf_p1(min_max) - ecdf_p2(min_max))))]
  y0 <- ecdf_p1(x0)
  y1 <- ecdf_p2(x0)

  ks_result <- stats::ks.test(
    x = dplyr::pull(data_wide, .data[[projects[1L]]]),
    y = dplyr::pull(data_wide, .data[[projects[2L]]]),
    alternative = "two.sided",
    exact = NULL,
    simulate.p.value = FALSE,
    B = 2000L
  ) |>
    (\(.) signif(x = .$statistic, digits = 2L))() |>
    (\(.) paste("D =", .))()

  # Main KS plot construction
  ks_plt <- ggplot2::ggplot(
    data = df,
    mapping = ggplot2::aes(
      x = .data[[check_log$col_names$quant]],
      group = .data[["Project"]],
      color = .data[["Project"]]
    )
  ) +
    ggplot2::stat_ecdf(
      linewidth = 1L
    ) +
    ggplot2::geom_segment(
      mapping = ggplot2::aes(
        x = x0[1L],
        y = y0[1L],
        xend = x0[1L],
        yend = y1[1L]
      ),
      linetype = "dashed",
      color = "black"
    ) +
    ggplot2::geom_point(
      mapping = ggplot2::aes(
        x = x0[1L],
        y = y0[1L]
      ),
      color = "black",
      size = 2L
    ) +
    ggplot2::geom_point(
      mapping = ggplot2::aes(
        x = x0[1L],
        y = y1[1L]
      ),
      color = "black",
      size = 2L
    ) +
    ggplot2::annotate(
      geom = "text",
      x = Inf,
      y = 0.1,
      hjust = 1,
      cex = 4,
      label = ks_result
    ) +
    ggplot2::labs(
      x = check_log$col_names$quant,
      y = "ECDF",
      color = "Platform:"
    ) +
    OlinkAnalyze::set_plot_theme() +
    OlinkAnalyze::olink_color_discrete() +
    ggplot2::theme(
      legend.position = "top",
      legend.title = ggplot2::element_blank()
    )

  return(ks_plt)
}

bridgeability_combine_plots <-  function(iqr, r2, counts, ks, title) {
  out_plot <- ggpubr::ggarrange(
    ggpubr::ggarrange(
      r2,
      ggpubr::ggarrange(
        iqr, counts,
        ncol = 2L,
        widths = c(1L, 1L),
        common.legend = TRUE
      ),
      ncol = 2L
    ),
    ks,
    nrow = 2L,
    heights = c(1L, 1L)
  )

  out_plot <- ggpubr::annotate_figure(
    p = out_plot,
    top = ggpubr::text_grob(
      label = title,
      size = 14L,
      just = "centre",
      face = "plain"
    )
  )

  return(out_plot)
}