SPLASH simulation

Overview

I implemented the SPLASH Chaung et al. (n.d.) data analysis to understand it. I use a subset of two SARS-CoV2 assembled genomes and simulate sample reads, some with read errors (noise), from that. From sample reads alone, we want to recover sequence variants. The SPLASH (née NOMAD), method tiles kmers–actually, a “reading frame” of an anchor and target. Then you can count up unique anchor-target pairs per sample (here, only 2).

Code

library(Biostrings)
library(rentrez)
library(tidyverse)
library(here)

Download

Code

cov2_du_oct23 <- entrez_fetch(db="nucleotide", id="OY751671.1", rettype="fasta", retmode="text")
cov2_uk_oct23 <- entrez_fetch(db="nucleotide", id="OY753298.1", rettype="fasta", retmode="text")
# cov2_wu_jul20 <- entrez_fetch(db="nucleotide", id="NC_045512.2", rettype="fasta", retmode="text")

# Write to FASTA files
write(cov2_du_oct23, file = "data/cov2_du_oct23.fasta")
write(cov2_uk_oct23, file = "data/cov2_uk_oct23.fasta")
# write(cov2_wu_jul20, file = "data/cov2_wu_jul20.fasta")

Align

Code

genome_du <- readDNAStringSet(here("data", "cov2_du_oct23.fasta"))
genome_uk <- readDNAStringSet(here("data", "cov2_uk_oct23.fasta"))
# genome_wu <- readDNAStringSet(here("data", "cov2_wu_jul20.fasta"))

cov2_align <- Biostrings::pairwiseAlignment(genome_uk, genome_du, type = "overlap")

cov2_align

Overlap PairwiseAlignmentsSingleSubject (1 of 1)
pattern: [53] --------------------------------...NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
subject:  [1] TATACCTTCCCAGGTAACAAACCAACCAACTT...---------------------------NNNN
score: 58027.25 

Code

missum <- cov2_align |> mismatchSummary()

missum$subject |> 
  head() |> 
  knitr::kable()

SubjectPosition	Subject	Pattern	Count	Probability
44	T	N	1	1
208	C	T	1	1
282	C	T	1	1
395	A	G	1	1
796	G	A	1	1
1921	A	C	1	1

Code

missum$subject |> 
  tail() |> 
  knitr::kable()

	SubjectPosition	Subject	Pattern	Count	Probability
63	28328	G	C	1	1
64	28380	G	T	1	1
65	28638	T	G	1	1
66	29663	G	T	1	1
67	29712	T	G	1	1
68	29794	A	N	1	1

Subset

Code

sub_du <- genome_du |> substr(200, 1200) |> as.character()
sub_uk <- genome_uk |> substr(200, 1200) |> as.character()

Define functions

Code

extract_anc_tar <- function(input_string, kmer_length = 27, step = 5) {
  L <- nchar(input_string)
  
  # space
  R <- max(0, (L - 2 * kmer_length) / 2)
  
  # Calculate the starting positions for each kmer
  start_positions <- seq(1, L - R - kmer_length + 1, by = step)
  
  # Extract kmers based on the starting positions
  pairs <- map(start_positions, function(pos) {
    list(
      anc = substr(input_string, pos, pos + kmer_length - 1),
      tar = substr(input_string, pos + R, pos + kmer_length + R - 1)
    )
  })
  
  return(pairs)
}

rand_mut <- function(input_string, max_mut = 2) {
  
  n <- nchar(input_string)
  
  # Number of characters to change (randomly between 0 and 2)
  num_changes <- sample(0:max_mut, 1)
  
  # Generate random positions to change
  positions_to_change <- sample(1:n, num_changes, replace = FALSE)
  
  # Generate random characters to replace
  new_chars <- sample(c("G", "T","A", "C"), num_changes, replace = TRUE)
  
  # Convert string to a list of characters
  char_list <- strsplit(input_string, "")[[1]]
  
  # Perform the changes
  char_list[positions_to_change] <- new_chars
  
  # Convert back to a string and return
  return(paste(char_list, collapse = ""))
}


sim_pairs <- function(genome, n_reads = 200, read_len = 300, max_mut = 10) {
  # Generate a random starting position
  start_positions <-
    sample(1:(nchar(genome) - (read_len - 1)), n_reads)
  
  start_positions |>
    # Extract a 300-character substring
    map( ~ substr(genome, start = .x, stop = .x + (read_len - 1))) |>
    map(rand_mut, max_mut) |>
    map(as.character) |>
    map(extract_anc_tar)
}

# Function to convert each sublist to a tibble row
frame_to_row <- function(item, read, frame) {
  tibble(read = read, frame = frame, anc = item$anc, tar = item$tar)
}

read_to_rows <- function(read){
  map2_dfr(read, seq_along(read), ~frame_to_row(.x, 1, .y))
}

samp_to_rows <- function(samp){
  samp |> 
    map(read_to_rows) |> 
    bind_rows()
}

Simulation

Code

set.seed(123)

# EXAMPLE PARAMETERS, many are defaults
# n_reads <- 300
# 
# L <- 300
# k = 27
# step = 5
# R <- max(0, (L - 2 * k) / 2)
# 
# read_len <- 300
# 
# 
# max_mut = 30



tbl_pairs_du <- 
  sub_du |> 
  sim_pairs(max_mut = 0) |> 
  samp_to_rows()

tbl_pairs_uk <- 
  sub_uk |> 
  sim_pairs(max_mut = 0) |> 
  samp_to_rows()

tbl_pairs_du_noise <- 
  sub_du |> 
  sim_pairs(max_mut = 10) |> 
  samp_to_rows()

tbl_pairs_uk_noise <- 
  sub_uk |> 
  sim_pairs(max_mut = 10) |> 
  samp_to_rows()

result <- bind_rows(uk = tbl_pairs_uk, du = tbl_pairs_du, .id = "samp")
result_noise <- bind_rows(uk = tbl_pairs_uk_noise, du = tbl_pairs_du_noise, .id = "samp")


cnt_res <- 
  result |> 
  count(anc, tar) |> 
  filter(n() > 1, .by = c(anc)) |> 
  arrange(anc) |> 
  # TODO: this is awkward
  left_join(result |> distinct(samp, anc, tar)) |> 
  arrange(anc) |> 
  pivot_wider(names_from = "samp", values_from = "n", values_fill = 0)

Joining with `by = join_by(anc, tar)`

Code

cnt_res_noise <- 
  result_noise |> 
  count(anc, tar) |> 
  filter(n() > 1, .by = c(anc)) |> 
  arrange(anc) |> 
  # TODO: this is awkward
  left_join(result_noise |> distinct(samp, anc, tar)) |> 
  arrange(anc) |> 
  pivot_wider(names_from = "samp", values_from = "n", values_fill = 0)

Joining with `by = join_by(anc, tar)`

Examine

Code

cnt_res

# A tibble: 72 × 4
   anc                         tar                            du    uk
   <chr>                       <chr>                       <int> <int>
 1 AAAGGTAAGATGGAGAGCCTTGTCCCT TTATCAGAGGCACGTCAACATCTTAAA     2     0
 2 AAAGGTAAGATGGAGAGCCTTGTCCCT TTATCAGAGGCACGTCAACATCTTAGA     0     4
 3 AAGATGGAGAGCCTTGTCCCTGGTTTC GAGGCACGTCAACATCTTAAAGATGGC     4     0
 4 AAGATGGAGAGCCTTGTCCCTGGTTTC GAGGCACGTCAACATCTTAGAGATGGC     0     5
 5 AAGGTAAGATGGAGAGCCTTGTCCCTG TATCAGAGGCACGTCAACATCTTAAAG     3     0
 6 AAGGTAAGATGGAGAGCCTTGTCCCTG TATCAGAGGCACGTCAACATCTTAGAG     0     4
 7 ACGGCGCCGATCTAAAGTCATTTGACT TTAACGGAGGGACATACACTCGCTATG     0    10
 8 ACGGCGCCGATCTAAAGTCATTTGACT TTAACGGAGGGGCATACACTCGCTATG     7     0
 9 AGATGGAGAGCCTTGTCCCTGGTTTCA AGGCACGTCAACATCTTAAAGATGGCA     5     0
10 AGATGGAGAGCCTTGTCCCTGGTTTCA AGGCACGTCAACATCTTAGAGATGGCA     0     4
# ℹ 62 more rows

Code

cnt_res_noise

# A tibble: 3,105 × 4
   anc                         tar                            uk    du
   <chr>                       <chr>                       <int> <int>
 1 AAAAAGGCGTTTTGCCTCAACTTGAAC GTCGTAGTGGTGAGACACTTGGTGTCC    12    12
 2 AAAAAGGCGTTTTGCCTCAACTTGAAC GTCGTAGTGGTGATACACTTGGTGTCC     1     0
 3 AAAACACACGTCCAACTCAGTTTGCCT GGCTTAGTAGAAGTTCAAAAAGGCGTT     1     0
 4 AAAACACACGTCCAACTCAGTTTGCCT GGCTTAGTAGAAGTTGAAAAAGGCGTT     8     8
 5 AAAACACACGTCCAACTCAGTTTGCCT TGCTTAGTAGAAGTTGAAAAAGGTGTT     1     0
 6 AAAACTGGAACACTAAACATAGCAGTG GCATGAAAGACCTTCTAGCACGTGCTG     0     1
 7 AAAACTGGAACACTAAACATAGCAGTG GCATTAAAGACCTTCTAGCACGTGCTG     5     5
 8 AAAACTGGAACACTAAACATAGCAGTG GCATTAAAGACCTTCTAGCACTTGCTG     1     0
 9 AAAACTGGAACACTAAACATAGCAGTG GCATTAAAGACTTTCTAGCACGTGCTG     0     1
10 AAAACTGGAACACTAAACATAGCAGTG GCATTAAGGACCTTCTAGCACGTGCTG     0     1
# ℹ 3,095 more rows

Code

cnt_res |>
  count(anc) |> 
  ggplot(aes(n)) +
  geom_histogram(binwidth = 1) +
  scale_x_continuous(breaks = 1:10) +
  xlab("Rows of contingency table") +
  labs(title = "Many contingency tables")

Code

cnt_res_noise |>
  count(anc) |> 
  ggplot(aes(n)) +
  geom_histogram(binwidth = 1) +
  scale_x_continuous(breaks = 1:10) +
  xlab("Rows of contingency table") +
  labs(title = "Many contingency tables")

Code

cnt_res |> 
  summarise(du_sum = sum(du),  uk_sum = sum(uk), .by = anc) |> 
  mutate(M = du_sum + uk_sum) |> 
  ggplot(aes(M)) +
  geom_histogram(binwidth = 1) +
  labs(title = "Total obs M per contingency table, no noise")

Code

cnt_res_noise |> 
  summarise(du_sum = sum(du),  uk_sum = sum(uk), .by = anc) |> 
  mutate(M = du_sum + uk_sum) |> 
  ggplot(aes(M)) +
  geom_histogram(binwidth = 1) +
  labs(title = "Total obs M per contingency table, with noise")

References

Chaung, Kaitlin, Tavor Z. Baharav, George Henderson, Ivan N. Zheludev, Peter L. Wang, and Julia Salzman. n.d. “SPLASH: A Statistical, Reference-Free Genomic Algorithm Unifies Biological Discovery.” https://doi.org/10.1101/2022.06.24.497555.