---
title: "Report: PGS Catalog Calculator Pipeline"
output: 
  html_document:
    toc: true
    toc_depth: 5
    toc_float: 
      collapsed: false
date: "`r Sys.time()`"
params:
  img: PGS_Logo.png
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(httr)
library(jsonlite)
library(dplyr)
library(ggplot2)

log_paths <- list.files(pattern = "*_summary.csv", 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 <- Reduce(dplyr::bind_rows, lapply(log_paths, read_log))
log$sampleset <- gsub("_", " ", log$sampleset)  # page breaking issues
```

<style>
#TOC {
  background: url("https://github.com/PGScatalog/pgsc_calc/blob/main/assets/PGS_Logo.png?raw=true");
  background-size: 33%;
  background-position: top;
  padding-top: 100px !important;
  background-repeat: no-repeat;
}
</style>

## Pipeline command 

```{r engine='bash', comment='', echo=FALSE}
cat command.txt
```

## Scoring file metadata 

<div class="alert alert-info", role="alert">
<p>Additional <a href="https://pgsc-calc.readthedocs.io/en/dev/explanation/output.html">documentation is available</a> that explains some of the terms used this report in more detail</p>
</div>

### Scoring files

```{r, echo = FALSE}
json_scorefiles <- read_json('log_combined.json', simplifyVector=TRUE)

parse_json2rows <- function(score_header){
  result <- list(pgs_id=NA, Publication = NA, Trait=NA, n_variants = as.integer(score_header[['variants_number']]), genome_build=NA)
  # pgs_id
  pgs_id <- ''
  if(is.null(score_header[['pgs_id']]) == FALSE){
    pgs_id <- paste0('<a href="https://www.pgscatalog.org/score/', score_header[['pgs_id']], '">', score_header[['pgs_id']], '</a>')
  }
  # Add name
  if(is.null(score_header[['pgs_name']]) == FALSE){
    if(pgs_id == ''){
      pgs_id <- score_header[['pgs_name']]
    }else{
      pgs_id <- paste0(pgs_id, '<br><small>(', score_header[['pgs_name']], ')</small>')
    }
  }
  result[['pgs_id']] <- pgs_id
  
  # pgp_id
  pgp_id <- ''
  if(is.null(score_header[['pgp_id']]) == FALSE){
    pgp_id <- paste0('<a href="https://www.pgscatalog.org/publication/', score_header[['pgp_id']], '">', score_header[['pgp_id']], '</a>')
  }
  if(is.null(score_header[['citation']]) == FALSE){
    if(pgp_id == ''){
      pgp_id <- paste0('<small>', score_header[['citation']], '</small>')
    }else{
      pgp_id <- paste(pgp_id, paste0('<br><small>', score_header[['citation']], '</small>'))
    }
  }
  if(pgp_id != ''){
    result[['Publication']] <- pgp_id
  }

  # Trait
  trait_mapped <- NA
  if(length(score_header[['trait_efo']]) > 0){
    efo <- data.frame(label = score_header[['trait_mapped']], id = score_header[['trait_efo']])
    efo$display <- paste0('<a href="http://www.ebi.ac.uk/efo/', efo$id, '">', efo$label,'</a>')
    trait_mapped <- paste0('<u>Mapped trait</u>: ', paste(efo$display, collapse=" "))
  }
  if(is.null(score_header[['trait_reported']])){
    result[['Trait']] <- trait_mapped
  }else if(is.na(trait_mapped)){
    result[['Trait']] <- paste0('<u>Reported trait</u>: ', score_header[['trait_reported']])
  }else{
    result[['Trait']] <- paste0('<u>Reported trait</u>: ', score_header[['trait_reported']], "<br>", trait_mapped)
  }
  
  # Genome build
  build_info <- paste0('<u>Reported</u>: ', score_header[['genome_build']])
  if(score_header[['use_harmonised']] == TRUE){
    build_info <- paste0(build_info, '<br><u>Harmonized Build</u>: ', score_header[['HmPOS_build']])
  }
  result[['genome_build']] <- build_info
  
  return(result)
} 
df_scorefiles <- do.call(rbind.data.frame, lapply(json_scorefiles, parse_json2rows))
df_scorefiles$accession <- rownames(df_scorefiles)

df_scorefiles %>%
  select(accession, pgs_id, Publication, Trait, n_variants, genome_build) %>% 
  DT::datatable(
    rownames = FALSE,
    escape = FALSE,
    colnames = c(
      "Scoring file" = "accession",
      "Polygenic Score ID" = "pgs_id",
      "Number of variants" = "n_variants",
      "Genome build" = "genome_build"
    ),
    extensions = 'Buttons',
    options = list(dom = 'Bfrtip',
                   buttons = c('csv'))
  )
```

### Variant matching

<hr>

#### Parameters

```{bash match_params, echo = FALSE}
cat "params.txt"
```

```{r match_summary, echo = FALSE, warning=FALSE}
log %>%
  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
```

<hr>

#### Summary

```{r match_table, echo = FALSE}
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')))
```

<hr>
#### Detailed results

```{r match_table_detailed, echo = FALSE, warning=FALSE}
log %>%
  group_by(sampleset, accession) %>%
  count(match_status, ambiguous, is_multiallelic, 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) %>%
  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",
      "%" = "percent"
    ))
```

## Scores 

```{r, echo = FALSE, message = FALSE}
scores <- readr::read_tsv("aggregated_scores.txt.gz") 
```

```{asis, echo = any(table(scores$sampleset) < 50)}
<div class="alert alert-warning", role="alert">
  <strong>⚠️ Warning: small sampleset size (n < 50)</strong>
  <hr>
  <ul>
  <li>plink2 uses allele frequency data to <a href="https://www.cog-genomics.org/plink/2.0/score">mean-impute the dosages of missing genotypes</a></li>
  <li>Currently the pipeline disables mean-imputation in these small sample sets to make sure that the calculated PGS is as consistent with the genotype data as possible</li>
  <li>With a small sample size, the resulting score sums may be inconsistent between samples</li>
  <li>The average ([scorename]_AVG) may be more applicable as it calculates an average weighting over all genotypes present</li>
  </ul> 
  <hr>
  <p>
  In the future mean-imputation will be supported in small samplesets using ancestry-matched reference samplesets to ensure consistent calculation of score sums (e.g. 1000G Genomes).
  </p>
</div>
```

<div class="alert alert-success", role="alert">
<p>`r sum(grepl("*_SUM$", colnames(scores)))` scores for `r nrow(scores)` samples processed</p>
</div>


```{r missing_scores, echo = FALSE}
log %>% 
  filter(score_pass == "false") %>%
  group_by(sampleset, accession) %>% 
  count(score_pass) -> missing_scores
```

```{asis, echo = nrow(missing_scores) > 0}
<div class="alert alert-danger", role="alert">
<p>Some scores are missing because they failed matching, check passed matching column 
in summary table</p>
</div>
```

### Score data {.tabset}

#### Score extract

<div class="alert alert-info", role="alert">
<p>Below is a summary of the aggregated scores, which might be useful for 
debugging.</p>
</div>


```{r, echo = FALSE}
scores %>%
  select(sampleset, IID, ends_with("SUM")) %>%
  tibble::as_tibble(.)
```

<div class="alert alert-info", role="alert">
<p>See here for an explanation of <a href="https://www.cog-genomics.org/plink/2.0/formats#sscore">plink2 column names</a>
</p>
</div>

#### Density plot

<div class="alert alert-info", role="alert">
<p>The summary density plots show up to six scoring files</p>
</div>

```{r, echo = FALSE, message=FALSE, warning=FALSE}
scores %>%
  ungroup() %>%
  select(IID, sampleset, ends_with("SUM")) %>%
  group_by(IID, sampleset) %>%
  tidyr::pivot_longer(cols = -group_cols()) %>%
  ungroup() %>%
  filter(name != "DENOM_SUM") %>%
  mutate(name = ifelse(
    stringr::str_detect(name, "^PGS"),
    stringr::str_extract(name, "^PGS[0-9]{6}"),
    name)) %>%
  group_by(sampleset, name) -> long_scores

group_keys(long_scores) %>%
  slice(1:6) -> score_keys

long_scores %>%
  inner_join(score_keys) %>%
  ggplot(aes(x = value, fill = sampleset)) +
      geom_density(alpha = 0.3) +
      facet_wrap(~name, ncol = 2, scales = "free") +
      theme_bw() +
      labs(x = "PGS (SUM)", y = "Density", title = "PGS Distribution(s)")
```

### Get scores

All scores can be found in "aggregated_scores.txt.gz", in the results folder output by the pipeline.

## Citations

<div class="alert alert-info", role="alert">
<p>For scores from the PGS Catalog, please remember to cite the original publications from which they came (these are listed in the metadata table.) </p>
</div>

> PGS Catalog Calculator (in development). PGS Catalog Team. https://github.com/PGScatalog/pgsc_calc

> Lambert et al. (2021) The Polygenic Score Catalog as an open database for reproducibility and systematic evaluation. Nature Genetics. 53:420–425 doi:10.1038/s41588-021-00783-5.
