Splash Results

Author

Matt Broerman

Code 
library(tidyverse)

Rationale

For the presentation about SPLASH and the statistical tests reported, it is helpful to know what the actual outputs of the SPLASH2 package are.

Summary

Running SPLASH

The authors provide an example run.

I tried to run SPLASH2 on my M1 Mac, enticed by the README saying this might be possible. Inside the Docker container that the authors provide, it was not. However, I tracked down that the error had to do with the resources available on my machine.

I then ran on the CRC cluster. This was not immediately possible, but I got related errors, so with a highly provisioned node (and not the default), I was able to run it.

Code 
crc-interactive -s -b 16 -c 6

The example script involves four clinical samples of SARS-CoV-2 co-infection from this experiment: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA817870

Code 
readLines("data/sample_name_to_id.mapping.txt")
[1] "SRX14565338_SRR18431620_1 0" "SRX14565342_SRR18431616_1 1"
[3] "SRX14565346_SRR18431612_1 2" "SRX14565361_SRR18431597_1 3"

There are two sets of data to consider, the results and the “dumps”.

Four file results

Code 
files <- list.files("data", "^result")

files
[1] "result.after_correction.all_anchors.tsv"              
[2] "result.after_correction.scores.top_target_entropy.tsv"
[3] "result.after_correction.scores.tsv"                   
[4] "result.compactors.tsv"

The dumps are the intermediate representation of the data, not provided by default. This is what the tests operate on.

Code 
splash  --dump_sample_anchor_target_count_binary --dump_sample_anchor_target_count_txt --bin_path .. input.txt

Here I look at only four pieces.

Code 
list.files("data", "^bin", full.names = TRUE)
[1] "data/bin0.satc.dump" "data/bin1.satc.dump" "data/bin2.satc.dump"
[4] "data/bin3.satc.dump"

Exploration

Code 
res_anc <- read_tsv("data/result.after_correction.all_anchors.tsv")
res_scr_ent <- read_tsv("data/result.after_correction.scores.top_target_entropy.tsv")
res_scr <- read_tsv("data/result.after_correction.scores.tsv")
res_com <- read_tsv("data/result.compactors.tsv")


res_scr |> 
  glimpse()
Rows: 10,227
Columns: 24
$ anchor                              <chr> "AAAAAAAGAACAAAGACCATTGAGTAC", "AA…
$ pval_opt                            <dbl> 2.53003e-06, 1.49488e-04, 4.67710e…
$ effect_size_bin                     <dbl> 0.1164390, 0.1005560, 0.8712200, 0…
$ pval_base                           <dbl> 2.24687e-03, 4.23241e-02, 3.82822e…
$ M                                   <dbl> 3649, 3422, 142, 4424, 700, 1023, …
$ anch_uniqTargs                      <dbl> 70, 59, 11, 70, 30, 44, 61, 120, 6…
$ number_nonzero_samples              <dbl> 4, 4, 2, 4, 4, 4, 4, 4, 4, 2, 4, 4…
$ target_entropy                      <dbl> 0.840874, 0.704373, 1.597170, 1.46…
$ avg_no_homopolymer_targets          <dbl> 0.000822143, 0.000876680, 0.000000…
$ avg_hamming_distance_max_target     <dbl> 1.090710, 0.808007, 1.429580, 0.68…
$ avg_hamming_distance_all_pairs      <dbl> 2.102380, 1.570760, 2.496550, 1.22…
$ avg_edit_distance_max_target        <dbl> 0.909290, 0.673875, 1.204230, 0.60…
$ avg_edit_distance_all_pairs         <dbl> 1.739520, 1.302940, 2.096690, 1.06…
$ anchor_2mer_seq_entropy             <dbl> 2.97759, 3.15180, 2.98393, 2.92103…
$ anchor_3mer_seq_entropy             <dbl> 3.86347, 4.01947, 3.68386, 3.73366…
$ most_freq_target_1                  <chr> "TCTGGACTAAAACTAAAAGTGAAGTCA", "CT…
$ cnt_most_freq_target_1              <dbl> 3319, 3172, 94, 3101, 440, 599, 32…
$ most_freq_target_1_2mer_seq_entropy <dbl> 3.22327, 3.13294, 3.35353, 3.53641…
$ most_freq_target_1_3mer_seq_entropy <dbl> 3.92386, 3.92386, 4.32386, 4.24386…
$ most_freq_target_2                  <chr> "CTGTCTCTTATACACATCTCCGAGCCC", "CT…
$ cnt_most_freq_target_2              <dbl> 60, 45, 31, 1010, 193, 315, 45, 18…
$ most_freq_target_2_2mer_seq_entropy <dbl> 3.59448, 3.19969, 3.58430, 3.48853…
$ most_freq_target_2_3mer_seq_entropy <dbl> 4.29366, 3.92386, 4.56386, 4.32386…
$ pval_opt_corrected                  <dbl> 3.430984e-04, 1.680252e-02, 1.0087…
Code 
dump <- 
  list.files("data", "^bin", full.names = TRUE) |> 
  read_tsv(col_names = c("fastq", "anchor", "trg", "n"), id = "src")

dump
# A tibble: 144,508 × 5
   src                 fastq                     anchor              trg       n
   <chr>               <chr>                     <chr>               <chr> <dbl>
 1 data/bin0.satc.dump SRX14565361_SRR18431597_1 AAAACAGTTACACAACAA… CAGA…     1
 2 data/bin0.satc.dump SRX14565361_SRR18431597_1 AAAACAGTTACACAACAA… CAGC…     1
 3 data/bin0.satc.dump SRX14565361_SRR18431597_1 AAAACAGTTACACAACAA… CAGG…     2
 4 data/bin0.satc.dump SRX14565346_SRR18431612_1 AAAACAGTTACACAACAA… CAGT…     1
 5 data/bin0.satc.dump SRX14565346_SRR18431612_1 AAAACAGTTACACAACAA… CAGT…     1
 6 data/bin0.satc.dump SRX14565361_SRR18431597_1 AAAACAGTTACACAACAA… CAGT…     1
 7 data/bin0.satc.dump SRX14565361_SRR18431597_1 AAAACAGTTACACAACAA… CAGT…     3
 8 data/bin0.satc.dump SRX14565342_SRR18431616_1 AAAACAGTTACACAACAA… CAGT…     1
 9 data/bin0.satc.dump SRX14565361_SRR18431597_1 AAAACAGTTACACAACAA… CAGT…     1
10 data/bin0.satc.dump SRX14565361_SRR18431597_1 AAAACAGTTACACAACAA… CAGT…     1
# ℹ 144,498 more rows

Note that in this example

  1. there are different numbers of rows per facet
  2. the counts are possibly incomplete (I am only using the first four “dump” files here, I suspect all are need)
Code 
dump |> 
  dplyr::count(src)
# A tibble: 4 × 2
  src                     n
  <chr>               <int>
1 data/bin0.satc.dump 36128
2 data/bin1.satc.dump 35981
3 data/bin2.satc.dump 35282
4 data/bin3.satc.dump 37117
Code 
bot_dump <- dump |> 
  distinct(anchor) |> 
  inner_join(res_scr) |> 
  arrange(pval_base) |> 
  select(anchor, pval_base, anch_uniqTargs, M) |> 
  tail(16)
Joining with `by = join_by(anchor)`
Code 
top_dump <- dump |> 
  distinct(anchor) |> 
  inner_join(res_scr) |> 
  arrange(pval_base) |> 
  select(anchor, pval_base, anch_uniqTargs, M) |> 
  head(16)
Joining with `by = join_by(anchor)`
Code 
top_dump
# A tibble: 16 × 4
   anchor                      pval_base anch_uniqTargs     M
   <chr>                           <dbl>          <dbl> <dbl>
 1 AGAGACTAGTGGCAATAAAACAAGAAA 0                    165 16998
 2 TGTCAAGCAGCAGCAAAGCAAGAGCAG 7.83e-265            189 22834
 3 TTACTGCCAGTTGAATCTGAGGGTCCA 1.15e-260            103  8209
 4 CGGGTGCCAATGTGATCTTTTGGTGTA 3.43e-220            130  8223
 5 GTATGGGAACTTCTCCTGCTAGAATGG 4.22e-198            136 14727
 6 CACCATCATCATACACAGTTCTTGCTG 2.05e-196             95  6436
 7 GGGTGCCAATGTGATCTTTTGGTGTAT 6.48e-168            131  8222
 8 ACAGTTCAAGAAATTCAACTCCAGGCA 2.78e-167            173 11463
 9 TGCGGGTGCCAATGTGATCTTTTGGTG 3.13e-154            122  8333
10 TCTTGTAGATCTGTTCTCTAAACGAAC 1.37e-148            145  6850
11 AAAGCTTTCGCTAGCATTTCAGTAGTG 3.54e-145             87  6904
12 TTCGTTCTATGAAGACTTTTTAGAGTA 2.79e-135             56  3803
13 TCAAAGCTTTCGCTAGCATTTCAGTAG 5.63e-132             91  7070
14 GTTAATTGGAACGCCTTGTCCTCGAGG 3.62e-129             94  5542
15 CTGTGTTGCTGATTATTCTGTCCTATA 7.28e-121            105  7119
16 TGTTGGCACTTTTCTCAAAGCTTTCGC 1.03e-120             73  6932

Now we can look at only anchorhors in a set with many nonzero targets. The methods sets M total counts as a threshold for the test.

Code 
anchor_cont_tbl_top <-
  dump |> 
  filter(anchor %in% top_dump$anchor) |> 
  separate(fastq, into = c("srx", "srr", "read")) |> 
  select(anchor, srr, trg, n) |> 
  pivot_wider(names_from = "srr", values_from = "n", values_fill = 0)

anchor_cont_tbl_bot <-
  dump |> 
  filter(anchor %in% bot_dump$anchor) |> 
  separate(fastq, into = c("srx", "srr", "read")) |> 
  select(anchor, srr, trg, n) |> 
  pivot_wider(names_from = "srr", values_from = "n", values_fill = 0)

anchor_cont_tbl_top_log <-
  dump |> 
  filter(anchor %in% top_dump$anchor) |> 
  separate(fastq, into = c("srx", "srr", "read")) |> 
  select(anchor, srr, trg, n) |> 
  mutate(n = log(n)) |>
  pivot_wider(names_from = "srr", values_from = "n", values_fill = 0)

anchor_cont_tbl_bot_log <-
  dump |> 
  filter(anchor %in% bot_dump$anchor) |> 
  separate(fastq, into = c("srx", "srr", "read")) |> 
  select(anchor, srr, trg, n) |> 
  mutate(n = log(n)) |>
  pivot_wider(names_from = "srr", values_from = "n", values_fill = 0)
Code 
anchor_cont_tbl_top_log |> 
  pivot_longer(!c(trg, anchor), names_to = "sample", values_to = "log_n") |> 
  ggplot(aes(sample, trg, fill = log_n)) +
  geom_raster() +
  facet_wrap(vars(anchor), scales = "free_y") +
  scale_fill_viridis_c() +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    axis.text.x = element_text(angle = -90)
  ) +
  ylab("Targets") +
  labs(title = "Log Counts", subtitle = "Best anchors")

Code 
anchor_cont_tbl_bot_log |> 
  pivot_longer(!c(trg, anchor), names_to = "sample", values_to = "log_n") |> 
  ggplot(aes(sample, trg, fill = log_n)) +
  geom_raster() +
  facet_wrap(vars(anchor), scales = "free_y") +
  scale_fill_viridis_c() +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    axis.text.x = element_text(angle = -90)
  ) +
  ylab("Targets") +
  labs(title = "Log Counts", subtitle = "Worst anchors")

Code 
anchor_cont_tbl <- anchor_cont_tbl_top

Comment

  • There are small row sums in most cases, probably due to measurement error
  • One sample seems to have higher counts of targets overall
  • At least one row has high counts

That is to say, across samples, within individual reads but at some fixed distance from an “anchor”, these are the counts of the target. There is at least one very common target. But in some cases, the a target exists in one set of samples, but not in another. This is a target

Each facet will be tested and in the base case, multiple times, each according to a random rowwise and columnwise partition.

Statistic

Below is a partial implementation of the SPLASH random statistic. TODO: use Hoeffding’s inequality to convert statistics to valid P-values

Intuitively, the “non-optimized” statistic is on a random partition of one “target table” (such as that above) where the partition f is on j rows and the partition c is on k columns. The random partition is conducted K times for f, and L times for c (assuming no metadata). Once the statistics are converted to p-values (and with Bonforoni correction and Benjamini-Yekutieli (BY) correction), the minimal p-value is reported.

From the paper, I wasn’t sure if proportions of the partition was fixed. I think no, but for the case of the columns (samples), when this is small, the test wouldn’t work. So for now, I use a fixed 50/50 partition.

To get a sense of this, let’s start with taking the mean of each partition, first rowwise.

Code 
part_f <- function(tbl, fn = NULL){
  tbl |> 
    mutate(v_f = sample(
      c(rep(0, floor(n()/2)), rep(1, n() - floor(n()/2))),
      size = n(),
      replace = FALSE
    )) |> 
    summarize(n_f = n(), across(starts_with("SRR"), mean), .by = v_f)
}
Code 
anchor_cont_tbl |> 
  filter(anchor == "CGGGTGCCAATGTGATCTTTTGGTGTA") |> 
  part_f()
# A tibble: 2 × 6
    v_f   n_f SRR18431616 SRR18431612 SRR18431597 SRR18431620
  <dbl> <int>       <dbl>       <dbl>       <dbl>       <dbl>
1     1    65        21.8        5.97        1.17       0.6  
2     0    65        23.2       19.6        53.4        0.831

Then column wise

Code 
part_c <- function(f_part){
  
  f_part |> 
    select(-n_f) |> 
    column_to_rownames(var = "v_f") |> 
    t() |> 
    as_tibble() |> 
    mutate(v_c = sample(
      c(rep(0, floor(n()/2)), rep(1, n() - floor(n()/2))),
      size = n(),
      replace = FALSE
    )) |> 
    summarise(n_c = n(), across(everything(), mean), .by = v_c)
}
Code 
anchor_cont_tbl |> 
  filter(anchor == "CGGGTGCCAATGTGATCTTTTGGTGTA") |> 
  part_f() |> 
  part_c()
# A tibble: 2 × 4
    v_c   n_c   `0`   `1`
  <dbl> <dbl> <dbl> <dbl>
1     0     2  7.38  15.8
2     1     2  2.13  37.9

Now this is done many times

Code 
repl_parts <- function(trg_tbl, n_rep_f, n_rep_c){
  
  rep_part_f <- replicate(n = n_rep_f, part_f(trg_tbl), simplify = FALSE)
  rep_part_c <- replicate(n = n_rep_c, map(rep_part_f, ~ part_c(.x)), simplify = FALSE)
  
  rep_part_c |> 
    map(~ bind_rows(.x, .id = "c_rep")) |> 
    bind_rows(.id = "f_rep") |> 
    select(-n_c)
}

anchor_cont_tbl |> 
  filter(anchor == "CGGGTGCCAATGTGATCTTTTGGTGTA") |> 
  repl_parts(n_rep_c = 2, n_rep_f = 3)
# A tibble: 12 × 5
   f_rep c_rep   v_c   `0`   `1`
   <chr> <chr> <dbl> <dbl> <dbl>
 1 1     1         1  46.4 3.38 
 2 1     1         0  12.6 0.931
 3 1     2         0  17.6 5.58 
 4 1     2         1  37.9 2.16 
 5 1     3         1  18.0 5.16 
 6 1     3         0  37.1 2.92 
 7 2     1         1  46.4 3.38 
 8 2     1         0  12.6 0.931
 9 2     2         0  28.5 6.73 
10 2     2         1  27.0 1.01 
11 2     3         1  43.8 5.95 
12 2     3         0  11.4 2.13

These statistics are converted to p-values, and we assign the minimum, after correction, to that anchor.

The actual splash statistic S is as follows. After partitioning with f:

Code 
S_j <- function(f1_mu_hat, f1_col_j) {
    sqrt(sum(f1_col_j)) * (mean(f1_col_j) - f1_mu_hat)
}

where f1_mu_hat is the total count of the partition (f1) over the total count of the target table.

So within partitioning:

Code 
part_f_S_j <- function(trg_tbl){
  
  trg_num <- trg_tbl[-1:-2] # remove bookkeeping
  M <- sum(trg_num)

  S_j <- function(mu_hat, f1_col_j) {
    sqrt(sum(f1_col_j)) * (mean(f1_col_j) - mu_hat)
  }

  trg_part <- trg_tbl |>
    mutate(v_f = sample(
      c(rep(0, floor(n()/2)), rep(1, n() - floor(n()/2))),
      size = n(),
      replace = FALSE
    )) |> 
    filter(v_f == 1)
  
  mu_hat <- trg_part |> 
    select(starts_with("SRR")) |> 
    sum()/M
  
  trg_part |>
    summarize(n_f = n(), across(starts_with("SRR"), ~ S_j(mu_hat, .x)), .by = v_f)
}
Code 
anchor_cont_tbl |> 
  filter(anchor == "CGGGTGCCAATGTGATCTTTTGGTGTA") |> 
  part_f_S_j()
# A tibble: 1 × 6
    v_f   n_f SRR18431616 SRR18431612 SRR18431597 SRR18431620
  <dbl> <int>       <dbl>       <dbl>       <dbl>       <dbl>
1     1    65        373.        45.7        14.0       0.915

Then for columnwise partitioning, we just do the weighted sum:

Code 
S <- function(v_c, S_j){
  sum(v_c * S_j)
}
Code 
part_c_S <- function(f_part){
  
  S <- function(v_c, S_j) {
    sum(v_c * S_j)
  }
  
  f_part |> 
    select(-n_f) |> 
    column_to_rownames(var = "v_f") |> 
    t() |> 
    as_tibble() |> 
    mutate(v_c = sample(
      c(rep(0, floor(n()/2)), rep(1, n() - floor(n()/2))),
      size = n(),
      replace = FALSE
    )) |> 
    summarise(S = S(v_c = v_c, S_j = `1`))
}

anchor_cont_tbl |> 
  filter(anchor == "CGGGTGCCAATGTGATCTTTTGGTGTA") |> 
  part_f_S_j() |> 
  part_c_S()
# A tibble: 1 × 1
      S
  <dbl>
1  692.
Code 
repl_parts_S <- function(trg_tbl, n_rep_f, n_rep_c){
  
  rep_part_f <- replicate(n = n_rep_f, part_f_S_j(trg_tbl), simplify = FALSE)
  rep_part_c <- replicate(n = n_rep_c, map(rep_part_f, ~ part_c_S(.x)), simplify = FALSE)
  
  rep_part_c |> 
    map(~ bind_rows(.x, .id = "c_rep")) |> 
    bind_rows(.id = "f_rep")
}
Code 
mk_noise_tbl <- function() {
  # matrix(sample(0:4, 130 * 4, replace = TRUE), nrow = 130) |>
  matrix(rpois(n = 130 * 4, lambda = 4), nrow = 130) |>  
    as_tibble() |>
    set_names(paste0("SRR", 1:4)) |>
    mutate(anchor = "AAA") |>
    mutate(trg = sample(1:1000, n())) |>
    relocate(anchor, trg)
}

noise_tbl <- mk_noise_tbl()
Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
`.name_repair` is omitted as of tibble 2.0.0.
ℹ Using compatibility `.name_repair`.
Code 
repl_s_noise <- 
  noise_tbl |> 
  repl_parts_S(n_rep_c = 10, n_rep_f = 10) |> 
  dplyr::rename(K = f_rep, L = c_rep)

repl_s_top <- anchor_cont_tbl |> 
  filter(anchor == "CGGGTGCCAATGTGATCTTTTGGTGTA") |> 
  repl_parts_S(n_rep_c = 10, n_rep_f = 10) |> 
  dplyr::rename(K = f_rep, L = c_rep)

repl_s_bot <- 
  anchor_cont_tbl_bot |> 
  filter(anchor == "AACTGCTTTTGGCTTAGTTGCAGAGTG") |> 
  repl_parts_S(n_rep_c = 10, n_rep_f = 10) |> 
  dplyr::rename(K = f_rep, L = c_rep)

bind_rows(
  "top_target" = repl_s_top,
  "bottom_target" = repl_s_bot,
  "noise" = repl_s_noise,
  .id = "repl_type"
) |> 
  ggplot(aes(S)) +
  geom_histogram() +
  facet_wrap(vars(repl_type))
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

These are

The choice of f and c that minimize the p-value are reported, and are used for computing the p-value of this anchor.