---
title: PGS Catalog Calculator (`pgsc_calc`) report
author: PGS Catalog Calculator (`pgsc_calc`)
date: today
title-block-banner: "#383434"
toc: true
format:
  html:
    code-fold: true
    anchor-sections: true
engine: knitr
params:
  log_pattern: "*_summary.csv"
  log_scorefiles: "log_scorefiles.json"
  score_path: ""
  sampleset: ""
  run_ancestry: false
  reference_panel_name: ""
  version: ""
css: logo.css
---

::: {.callout-note}
See the online [documentation](https://pgsc-calc.readthedocs.io/en/latest/explanation/output.html) for additional
explanation of the terms and data presented in this report.
:::

```{r setup, echo=FALSE, warning=FALSE, message=FALSE}
library(jsonlite)
library(dplyr)
library(tidyr)
library(stringr)
library(purrr)
library(ggplot2)
library(DT)
library(tibble)
library(forcats)
library(readr)

# prevent plots with small sample sets
MINIMUM_N_SAMPLES <- 50
LOW_SAMPLE_SIZE <- TRUE
```

```{r setup_logs, echo=FALSE}
log_paths <- list.files(pattern = params$log_pattern, full.names = TRUE)
read_log <- function(path) {
  log <- read.csv(path, stringsAsFactors = FALSE)
  
  return(
    log %>%
      mutate(
        is_multiallelic =  factor(is_multiallelic, levels = c("false", "true")),
        ambiguous = factor(ambiguous, levels = c("false", "true"))
      ) %>%
      rename(sampleset = dataset)
  )
}
log_df <- Reduce(dplyr::bind_rows, lapply(log_paths, read_log))
log_df$sampleset <- gsub("_", " ", log_df$sampleset)  # page breaking issues
```

# Workflow metadata

## Command

```{bash}
cat command.txt | fold -w 80 -s | awk -F ' ' 'NR==1 { print "$", $0} NR>1 { print "    " $0}' | sed 's/$/\\/' | sed '$ s/.$//' 
```

```{asis, echo = grepl("-profile test", readLines("command.txt"))}
:::{.callout-tip}
* If you're using the test profile, this report and these results are not biologically meaningful 
* The test profile is only used to check that all software is installed and working correctly 
* If you're reading this message, then that means everything is OK and you're ready to use your own data!
:::
```

## Version

```{r, echo=FALSE}
message(params$version)
``` 

# Scoring file metadata

## Scoring file summary

```{r load_scorefiles}
json_list <- jsonlite::fromJSON(params$log_scorefiles, simplifyVector = FALSE)

link_traits <- function(trait_efo, mapped) {
  if (length(trait_efo) == 0) {
    return("")
  } else {
    return(purrr::map2_chr(trait_efo, mapped, ~ stringr::str_glue('<a href="http://www.ebi.ac.uk/efo/{.x}">{.y}</a>')))
  }
}

extract_traits <- function(x) {
  trait_efo <- purrr::map(x, ~ extract_chr_handle_null(.x$header, "trait_efo"))
  mapped <- purrr::map(x, ~ extract_chr_handle_null(.x$header, "trait_mapped"))
  trait_display <- purrr::map2(trait_efo, mapped, link_traits)
  mapped_trait_links <- purrr::map_chr(trait_display, ~ paste(.x, collapse = "<br />"))
  reported_traits <- purrr::map(x, ~ extract_chr_handle_null(.x, "trait_reported"))
  purrr::map2_chr(reported_traits, mapped_trait_links, ~ {
    stringr::str_glue("<u>Reported trait:</u> {.x} <br /> <u>Mapped trait(s):</u> {.y}")
  })
}

extract_chr_handle_null <- function(x, field) {
  return(replace(x[[field]], is.null(x[[field]]), ""))
}

link_pgscatalog <- function(id, link_type) {
  if (id != "") {
    return(stringr::str_glue('<a href="https://www.pgscatalog.org/{link_type}/{id}">{id}</a>'))
  } else {
    return(id)
  }
}

add_note <- function(id, note) {
  if (id != "") {
    return(stringr::str_glue("{id} <br /> <small>{note}</small>"))
  } else {
    return(id)
  }
}

annotate_genome_build <- function(original_build, harmonised_build) {
  return(stringr::str_glue("<u>Original build:</u> {original_build} <br /> <u>Harmonised build:</u> {harmonised_build}"))
}

# extract fields from json list
tibble(
  pgs_id = map_chr(json_list, "pgs_id"),
  pgs_name = map_chr(json_list, ~ extract_chr_handle_null(.x$header, "pgs_name")),
  pgp_id = map_chr(json_list, ~ extract_chr_handle_null(.x$header, "pgp_id")),
  citation = map_chr(json_list, ~ extract_chr_handle_null(.x$header, "citation")),
  trait_display = extract_traits(json_list),
  genome_build = purrr::map_chr(json_list, ~ extract_chr_handle_null(.x$header, "genome_build")),
  harmonised_build = purrr::map_chr(json_list, ~ extract_chr_handle_null(.x$header, "HmPOS_build")),
  n_variants = purrr::map_chr(json_list, ~ extract_chr_handle_null(.x$header, "variants_number")),
  compatible_effect_type = map_lgl(json_list, "compatible_effect_type"),
  has_complex_alleles = map_lgl(json_list, "has_complex_alleles")) %>%
  # add links to pgs catalog identifiers
  mutate(pgs_id = purrr::map_chr(pgs_id, ~ link_pgscatalog(.x, "score")),
         pgp_id = purrr::map_chr(pgp_id, ~ link_pgscatalog(.x, "publication"))) %>%
  # add notes
  mutate(pgp_id = purrr::map2_chr(pgp_id, citation, ~ add_note(.x, .y)),
         pgs_id = purrr::map2_chr(pgs_id, pgs_name, ~ add_note(.x, .y)),
         genome_build = purrr::map2_chr(genome_build, harmonised_build, ~ annotate_genome_build(.x, .y))) %>% 
  # pick columns
  select(pgs_id, pgp_id, trait_display, n_variants, genome_build, has_complex_alleles, compatible_effect_type) -> scorefile_metadata
```

:::{.column-body-outset}

```{r, echo=FALSE}
tooltip_text <- c(
  "Polygenic Score ID" = "Unique identifier for the polygenic score.",
  "Publication" = "Reference publication for the score.",
  "Traits" = "Traits associated with the score.",
  "Number of variants" = "Total number of genetic variants (defined in the header)",
  "Genome build" = "The genome assembly version used.",
  "Complex alleles present?" = "Describes if complex non-SNP alleles included in the scoring file, e.g. APOE/HLA. These variants are excluded from the PGS calculation in the current version",
  "Effect types compatible?" = "Describes if the scoring file is compatible with the Calculator. Scores with dosage-specific weights are removed."
)

DT::datatable(
  scorefile_metadata,
  rownames = FALSE,
  escape = FALSE,
  colnames = setNames(
    paste0('<span title="', tooltip_text, '">', names(tooltip_text), '</span>'), 
    NULL
  ),
  extensions = 'Buttons',
  options = list(dom = 'Bfrtip',
                 buttons = c('csv'))
) %>%
    DT::formatStyle('has_complex_alleles', 
                  backgroundColor = DT::styleEqual(c(FALSE, TRUE), c('#a6dba0', '#c2a5cf'))) %>% 
    DT::formatStyle('compatible_effect_type', 
                  backgroundColor = DT::styleEqual(c(FALSE, TRUE), c('#c2a5cf', '#a6dba0')))

```

:::

```{asis, echo = any(!scorefile_metadata$compatible_effect_type)}
::: {.callout-warning title="Incompatible effect types detected"}
* Some scoring files contain variants with dosage dependent effect weights (for example, [PGS002253](https://www.pgscatalog.org/score/PGS002253/))
* Scores with dosage-specific weights are removed from calculation 
* Scores that contain variants with recessive, dominant, or additive effect types are supported
:::
```

```{asis, echo = any(scorefile_metadata$has_complex_alleles)}
::: {.callout-warning title="Complex alleles detected"}
* Some scoring files contain complex alleles (e.g. APOE / HLA / CYP)
* These variants are excluded from the PGS calculation in the current version
* Please check [Appendix A - Curation of PGS including complex alleles](https://www.pgscatalog.org/docs/curation) for more detailed information
:::
```

# Variant matching 

## Parameters

```{bash}
cat params.txt
```

```{asis, echo = params$run_ancestry}
## Reference matching summary
```

```{r, echo = FALSE, message = FALSE, eval=params$run_ancestry}
intersect_stats_files <- list.files(pattern = "intersect_counts*")
intersect_stats <- lapply(intersect_stats_files, function(x) scan(x, sep = "\n", what=integer()))

n_target <- purrr::map_int(intersect_stats, ~ .x[[1]])
n_ref <- purrr::map_int(intersect_stats, ~ .x[[2]])
n_match <- purrr::map_int(intersect_stats, ~ .x[[3]])

data.frame(reference = params$reference_panel_name, n_target = n_target, n_match = n_match, n_ref = n_ref) %>%
  group_by(reference) %>%
  summarise(n_target = sum(n_target), n_ref = sum(n_ref), n_match = sum(n_match)) %>%
  mutate("% matched" = round(n_match / n_ref * 100, 2)) %>%
  rename("n target" = n_target, "n (matched)" = n_match, "N variants in panel" = n_ref) %>%
  DT::datatable()
```

## Summary

```{r setup_matches}
log_df %>%
  mutate(match_status = forcats::fct_collapse(match_status, matched = "matched", other_level = "unmatched")) %>%
  group_by(sampleset, accession, match_status, score_pass) %>%
  count(wt = count) %>%
  group_by(sampleset, accession) %>%
  mutate(percent = round(n / sum(n) * 100, 1), n_variants = sum(n)) %>%
  arrange(accession, desc(percent)) %>%
  tidyr::pivot_wider(names_from = match_status, values_from = c(n, percent)) %>%
  replace(is.na(.), 0) -> compat
```

```{r match_table}
if (!"n_unmatched" %in% colnames(compat)) {
  # handle missing column if all PGS matches perfectly (e.g. no unmatched or excluded variants)
  compat <- compat %>%
    mutate(n_unmatched = 0) 
}

compat %>%
  select(sampleset, accession, n_variants, score_pass, percent_matched,
         n_matched, n_unmatched) %>%
  mutate(score_pass = as.logical(score_pass)) %>%
  DT::datatable(rownames = FALSE,
                extensions = 'Buttons',
    options = list(dom = 'Bfrtip',
                   buttons = c('csv')),
    colnames = c(
      "Sampleset" = "sampleset",
      "Scoring file" = "accession",
      "Number of variants" = "n_variants",
      "Passed matching" = "score_pass",
      "Match %" = "percent_matched",
      "Total matched" = "n_matched",
      "Total unmatched" = "n_unmatched"
    )) %>%
  DT::formatStyle('Scoring file', 
                  valueColumns = 'Passed matching',
                  backgroundColor = DT::styleEqual(c(FALSE, TRUE), c('#c2a5cf', '#a6dba0')))
```

## Detailed log

:::{.column-body-outset}

```{r match_table_detailed, echo = FALSE, warning=FALSE}
if(params$run_ancestry == TRUE){
    # Include match_IDs in the groupby to account for
    log_df %>%
      group_by(sampleset, accession) %>%
      count(match_status, match_IDs, ambiguous, is_multiallelic, match_flipped, duplicate_best_match, duplicate_ID, wt = count) %>%
      rename(is_ambiguous = ambiguous) %>%
      mutate(percent = round(n / sum(n) * 100, 2),
             match_status = forcats::fct_relevel(match_status, "matched", "excluded", "unmatched")) %>%
      arrange(accession, match_status) %>%
      mutate(accession = stringr::str_replace_all(accession, "_", " ")) %>%
      DT::datatable(rownames=FALSE,
                    extensions = 'Buttons',
        options = list(dom = 'Bfrtip',
                       buttons = c('csv')),
        colnames = c(
          "Sampleset" = "sampleset",
          "Scoring file" = "accession",
          "Match type" = "match_status",
          "Variant in reference panel" = "match_IDs",
          "Multiple potential matches" = "duplicate_best_match",
          "Duplicated matched variants" = "duplicate_ID",
          "Ambiguous" = "is_ambiguous",
          "Multiallelic" = "is_multiallelic",
          "Matches strand flip" = "match_flipped",
          "%" = "percent"
        ))
} else{
    log_df %>%
      group_by(sampleset, accession) %>%
      count(match_status, ambiguous, is_multiallelic,match_flipped, duplicate_best_match, duplicate_ID, wt = count) %>%
      rename(is_ambiguous = ambiguous) %>%
      mutate(percent = round(n / sum(n) * 100, 2),
             match_status = forcats::fct_relevel(match_status, "matched", "excluded", "unmatched")) %>%
      arrange(accession, match_status) %>%
      mutate(accession = stringr::str_replace_all(accession, "_", " ")) %>%
      DT::datatable(rownames=FALSE,
                    extensions = 'Buttons',
        options = list(dom = 'Bfrtip',
                       buttons = c('csv')),
        colnames = c(
          "Sampleset" = "sampleset",
          "Scoring file" = "accession",
          "Match type" = "match_status",
          "Multiple potential matches" = "duplicate_best_match",
          "Duplicated matched variants" = "duplicate_ID",
          "Ambiguous" = "is_ambiguous",
          "Multiallelic" = "is_multiallelic",
          "Matches strand flip" = "match_flipped",
          "%" = "percent"
        ))
}
```

:::

```{asis, echo = params$run_ancestry}
# Genetic Ancestry
```

```{r colour_palette, echo = FALSE, eval=params$run_ancestry}
# source: https://github.com/PGScatalog/PGS_Catalog/blob/master/catalog/static/catalog/pgs.scss#L2493-L2520
# $ancestry_colours
if({params$reference_panel_name} == '1000G'){
  thousand_genomes_colours <- c("#FFD900", "#E41A1C", "#B15928", "#4DAF4A",
                                "#377EB8", "#00CED1", "#984EA3", "#A6CEE3",
                                "#FF7F00", "#BBB", "#999")
  names(thousand_genomes_colours) <- c("AFR", "AMR", "ASN", "EAS",
                                       "EUR", "GME", "SAS", "MAE",
                                       "MAO", "NR", "OTH")
  current_population_palette <- scale_colour_manual(name = "Populations", values = thousand_genomes_colours)
} else if({params$reference_panel_name} == 'HGDP+1kGP'){
  gnomAD_pop_colours <- c("#97519d", "#e42523", "#f67e1e", "#48b24b",
                          "#3280bb", "#a65528", "#9a9c9b")
  names(gnomAD_pop_colours) <- c("AFR", "AMR", "CSA", "EAS",
                                 "EUR", "MID", "OCE")
  current_population_palette <- scale_colour_manual(name = "Populations", values = gnomAD_pop_colours)
} else{
  current_population_palette <- scale_colour_brewer(palette="Set3")
}
```

```{r, echo = FALSE, message = FALSE, eval=params$run_ancestry}
popsim <- readr::read_tsv(gsub("_pgs.", "_popsimilarity.", params$score_path))
head(popsim)

new_label_target = paste({params$sampleset}, '(Most Similar Population)')
new_label_ref = paste0('reference (', {params$reference_panel_name}, ' populations)')

popsim$slabel <- new_label_ref
popsim[popsim$sampleset == {params$sampleset}, 'slabel'] <- new_label_target

map_shapes <- c(1, 16)
map_shapes <- setNames(map_shapes, c(new_label_target, new_label_ref))

# Placeholder figure: needs better legend labelling
for(pc in seq.int(1,5,2)){
  pcY = paste('PC', pc, sep='')
  pcX = paste('PC', pc+1, sep='')
  if (pcX %in% colnames(popsim)){
    p_pca <- ggplot(popsim[popsim$REFERENCE == TRUE,], aes(x=!!sym(pcX), y=!!sym(pcY))) + geom_point(aes(colour=SuperPop, shape=slabel), alpha=0.25)
    p_pca <- p_pca + geom_point(data=popsim[popsim$REFERENCE != TRUE,], aes(color=MostSimilarPop, shape=slabel))
    p_pca <- p_pca + theme_bw() + current_population_palette + scale_shape_manual(values=map_shapes, name='sampleset')
    print(p_pca)
  }
}
```

```{asis, echo = params$run_ancestry}
## Population similarity summary
```

```{r, echo = FALSE, message = FALSE, eval=params$run_ancestry}
popsim %>% 
  filter(sampleset == "reference") %>%
  group_by(sampleset) %>%
  count(`Most similar population` = SuperPop ) %>%
  mutate(percent = round(n / sum(n) * 100, 2)) -> ref_count

popsim %>%
  filter(sampleset != "reference") %>%
  group_by(sampleset) %>%
  count(`Most similar population` = MostSimilarPop) %>%
  mutate(percent = round(n / sum(n) * 100, 2)) %>%
  bind_rows(ref_count) %>%
  mutate(count = stringr::str_glue("{n} ({percent}%)")) %>%
  select(-n, -percent) %>%
  tidyr::pivot_wider(names_from = sampleset, values_from = count) -> pop_summary

write.csv(pop_summary, "pop_summary.csv", quote=FALSE, row.names = FALSE)

pop_summary %>%
  DT::datatable()
```

# Scores

```{r, echo = FALSE, message = FALSE}
scores <- readr::read_tsv(params$score_path) 
n_scores <- length(unique(scores$PGS))
n_samples <- length(unique(scores$IID))
if (n_samples < MINIMUM_N_SAMPLES) {
  LOW_SAMPLE_SIZE <- TRUE
} else {
  LOW_SAMPLE_SIZE <- FALSE
}
```


```{asis, echo = (LOW_SAMPLE_SIZE && !params$run_ancestry)}

::: {.callout-important title="Warning: small sampleset size (n < 50) detected"}
* plink2 uses allele frequency data to [mean impute](https://www.cog-genomics.org/plink/2.0/score) the dosages of missing genotypes
* Currently `pgsc_calc` disables mean-imputation in these small sample sets to make sure that the calculated PGS is as consistent with the genotype data as possible
* With a small sample size, the resulting score sums may be inconsistent between samples
* The average `([scorename]_AVG)` may be more applicable as it calculates an average weighting over all genotypes present

It's recommended to use `--run_ancestry` with small samplesets to ensure consistent calculation of score sums (e.g. 1000G Genomes).
:::

```

```{asis, echo = (nrow(compat) == n_scores)}
:::{.callout-tip title="Success"}
* All requested scores were calculated successfully
:::
```

```{asis, echo = (nrow(compat) != n_scores)}
:::{.callout-caution}
* Some requested scores could not be calculated
* Please check the variant matching summary table to understand why
:::
```

`r n_scores` scores for `r n_samples` samples processed.

### Score data 

#### Density plot(s)

```{asis, echo = !LOW_SAMPLE_SIZE}
::: {.callout-note}
The summary density plots show up to six scoring files
:::
```

```{asis, echo = LOW_SAMPLE_SIZE}
::: {.callout-warning}
Density plots are disabled for low sample sizes
:::
```

```{r density_ancestry, echo=FALSE, message=FALSE, warning=FALSE, eval=(!LOW_SAMPLE_SIZE & params$run_ancestry)}
# Select which PGS to plot
uscores <- unique(scores$PGS)
uscores_plot <- uscores[1:min(length(uscores), 6)] # plot max 6 PGS

# Plot multiple adjustment methods at once per PGS
for(current_pgs in uscores_plot){
  long_scores <- scores %>% select(!percentile_MostSimilarPop) %>% filter(PGS == current_pgs) %>% gather(Method, score, -sampleset, -FID, -IID, -PGS)
  long_scores %>%
  ggplot(aes(x = score, fill = sampleset)) +
    geom_density(alpha = 0.3) +
    facet_wrap(~Method, scales = "free", nrow=1) +
    theme_bw() +
    labs(x = 'PGS', y = "Density", title = paste(current_pgs, '(adjusted distributions)')) -> p_dist
  print(p_dist)
}
```

```{r, echo = FALSE, message=FALSE, warning=FALSE, eval=(!LOW_SAMPLE_SIZE & !params$run_ancestry)}
scores %>%
  ungroup() %>%
  select(IID, sampleset, PGS, SUM) %>%
  mutate(PGS = ifelse(
    stringr::str_detect(PGS, "^PGS"),
    stringr::str_extract(PGS, "^PGS[0-9]{6}"),
    PGS)) %>%
  group_by(sampleset, PGS) -> long_scores
group_keys(long_scores) %>%
  slice(1:6) -> score_keys
long_scores %>%
  inner_join(score_keys) %>%
  ggplot(aes(x = SUM, fill = sampleset)) +
      geom_density(alpha = 0.3) +
      facet_wrap(~PGS, ncol = 2, scales = "free") +
      theme_bw() +
      labs(x = "PGS (SUM)", y = "Density", title = "PGS Distribution(s)")
```

### Get all scores

All scores can be found in the results directory, at:

```{r, eval=params$run_ancestry, echo=FALSE}
stringr::str_glue("{params$sampleset}/score/{params$sampleset}_pgs.txt.gz")
```

```{r, eval=!params$run_ancestry, echo=FALSE}
stringr::str_glue("{params$sampleset}/score/aggregated_scores.txt.gz")
```

# Citation

> Samuel A. Lambert, Benjamin Wingfield, Joel T. Gibson, Laurent Gil, Santhi Ramachandran, Florent Yvon, Shirin Saverimuttu, Emily Tinsley, Elizabeth Lewis, Scott C. Ritchie, Jingqin Wu, Rodrigo Canovas, Aoife McMahon, Laura W. Harris, Helen Parkinson, Michael Inouye.
Enhancing the Polygenic Score Catalog with tools for score calculation and ancestry normalization.
Nature Genetics (2024) | doi: [10.1038/s41588-024-01937-x](https://doi.org/10.1038/s41588-024-01937-x)

::: {.callout-important}
For scores from the PGS Catalog, please remember to cite the original publications from which they came (these are listed in the metadata table).
:::

## Score licenses

::: {.callout-tip}
* Scores deposited in the PGS Catalog may have specific license terms
* It's important to follow the license terms when you reuse scoring files
* Please check below for a summary of license terms
* License terms for custom scoring files aren't reported here, please check how the creators of the scoring file licensed their data
:::

```{r}
# as of 2023-12-12 only non-default licenses are recorded in the scoring file header
default_ebi_terms <- "PGS obtained from the Catalog should be cited appropriately, and used in accordance with any licensing restrictions set by the authors. See EBI Terms of Use (https://www.ebi.ac.uk/about/terms-of-use/) for additional details."

tibble(
    pgs_id = map_chr(json_list, "pgs_id"),
    license_text = map_chr(json_list, ~ extract_chr_handle_null(.x$header, "license"))) %>%
  mutate(license_text = ifelse(license_text == "", default_ebi_terms, license_text)) %>%
  # display license terms for files in the PGS Catalog only (with a PGS ID)
  filter(startsWith(pgs_id, "PGS")) %>%
  DT::datatable(., colnames = c(
      "PGS ID" = "pgs_id",
      "License text" = "license_text"
    ))
```

