Report breakdown by amplicon sequence
In amplican: Automated analysis of CRISPR experiments

library(amplican)
library(ggplot2)
alignments <- data.table::fread(params$alignments)
data.table::setDF(alignments)
config <- data.frame(data.table::fread(params$config_summary))
config$AmpliconUpper <- toupper(config$Amplicon)
uniqueAmlicons <- unique(config$AmpliconUpper)
labels <- sapply(uniqueAmlicons, function(x) {
  toString(config$ID[config$AmpliconUpper == x])
})
height <- amplican::plot_height(length(uniqueAmlicons))

Description

Read distribution plot - plot shows number of reads assigned during read grouping
Filtered Reads - plot shows percentage of assigned reads that have been recognized as PRIMER DIMERS or filtered based on low alignment score
Edit rates - plot gives overview of percentage of reads (not filtered as PRIMER DIMER) that have edits
Frameshift - plot shows what percentage of reads that have frameshift
Read heterogeneity plot - shows what is the share of each of the unique reads in total count of all reads. The more yellow each row, the less heterogeneity in the reads, more black means reads don't repeat often and are unique
Deletions plot - shows summary of deletions detected after alignments with distinction for forward (top plot) and reverse (bottom) reads, blue dotted lines represent primers as black dotted line represents cut site box, for deletions overlapping with cut site box there is distinction in color
Mismatches plot - shows summary of mismatches detected after alignments split by forward (top plot) and reverse (bottom) reads, mismatches are colored in the same manner as amplicon
Insertions plot - shows summary of insertions detected after alignments split by forward (top plot) and reverse (bottom) reads, insertion is shown as right-angled triangle where size of the insertion corresponds to the width of the triangle, size and transparency of triangle reflect on the frequency of the insertion
Cuts plot - shows summary of cuts in the amplicon (insertions that are overlapping with specified box of uppercase letters in the amplicon, and are supported by both forward and reverse reads)

Amplicon Summary

Amplicon groups

library(knitr)
ampliconDF <- data.frame(group = 1:length(uniqueAmlicons), IDs = unname(labels))
kable(ampliconDF)
ampliconDF$Amplicon <- uniqueAmlicons #for sorting
config$group <- match(config$AmpliconUpper, ampliconDF$Amplicon)

Read distribution

ampliconTable <- aggregate(cbind(Reads, Reads_Filtered, PRIMER_DIMER, Low_Score,
                                 Reads_Edited, Reads_Frameshifted) ~ group, 
                           data = config, sum)

ggplot(data = ampliconTable, aes(x = as.factor(group), y = log10(Reads + 1)), order = group) +
  geom_bar(stat = 'identity') +
  ylab('Number of reads + 1, log10 scaled')  +
  xlab('Amplicon Group') +
  theme(legend.position = 'none',
        axis.text = element_text(size = 12),
        axis.title = element_text(size = 14, face = 'bold')) +
  coord_flip() +
  geom_text(aes(x = as.factor(group), y = log10(Reads + 1), label = Reads), hjust = -1)

Filtered Reads

ampliconTable$F_percentage <- 
  (ampliconTable$PRIMER_DIMER + ampliconTable$Low_Score) * 100/ampliconTable$Reads
ampliconTable$F_percentage[is.nan(ampliconTable$PD_percentage)] <- 0 

ggplot(data = ampliconTable, aes(x = as.factor(group), y = F_percentage, order = group)) +
  geom_bar(stat='identity') +
  ylab('Percentage of filtered reads')  +
  xlab('Amplicon Group') +
  theme(axis.text = element_text(size=12),
        axis.title = element_text(size=14, face = 'bold')) +
  ylim(0, 100) +
  coord_flip() +
  geom_text(aes(x = as.factor(group), y = F_percentage, 
                label = PRIMER_DIMER), hjust = -1)

Edit rates

ampliconTable$edit_percentage <- 
  ampliconTable$Reads_Edited * 100/ampliconTable$Reads_Filtered
ampliconTable$edit_percentage[is.nan(ampliconTable$edit_percentage)] <- 0  

ggplot(data = ampliconTable, aes(x = as.factor(group), 
                                 y = edit_percentage, order = group)) +
  geom_bar(stat = 'identity') +
  ylab('Percentage of reads (not filtered) that have edits')  +
  xlab('Amplicon Group') +
  theme(axis.text = element_text(size=12),
        axis.title = element_text(size=14, face = 'bold')) +
  ylim(0,100) +
  coord_flip() +
  geom_text(aes(x = as.factor(group), y = edit_percentage, 
                label = Reads_Edited), hjust = -1)

Frameshift

ampliconTable$frameshift_percentage <- 
  ampliconTable$Reads_Frameshifted * 100/ampliconTable$Reads_Filtered
ampliconTable$frameshift_percentage[is.nan(ampliconTable$frameshift_percentage)] <- 0 

ggplot(data = ampliconTable, aes(x = as.factor(group), 
                                 y = frameshift_percentage, order = group)) +
  geom_bar(position = 'stack', stat = 'identity') +
  ylab('Percentage of reads (not filtered) that have frameshift')  +
  xlab('Amplicon Group') +
  theme(axis.text = element_text(size=12),
        axis.title = element_text(size=14,face = 'bold')) +
  ylim(0, 100) +
  coord_flip() +
  geom_text(aes(x = as.factor(group), y = frameshift_percentage, 
                label = Reads_Frameshifted), hjust = -1)

Heterogeneity of reads

plot_heterogeneity(alignments, config, level = 'group')

Alignments plots

alignments <- alignments[alignments$consensus & alignments$overlaps, ]
alignments$strand <- "+" # strand does not matter after consensus filtering
src = sapply(seq_along(uniqueAmlicons), function(i) {
  ID = config$ID[config$AmpliconUpper == uniqueAmlicons[i]]
  ID = paste0("c('", paste0(ID, collapse = "','"), "')", collapse = "")
  knitr::knit_expand(text = c(
    "## Group {{i}}  \n", 
    "### Deletions  \n", 
    paste('```r}, echo = F, results = "asis", ',
          'fig.width=25, message=F, warning=F}', collapse = ''), 
    paste('p <- amplican::plot_deletions(alignments, config, {{ID}},',
          ' params$cut_buffer, params$xlab_spacing)', collapse = ''), 
    '```\n',
    "### Insertions  \n", 
    paste('```r}, echo = F, results = "asis", ',
          'fig.width=25, message=F, warning=F}', collapse = ''), 
    paste('p <- amplican::plot_insertions(alignments, config, {{ID}},',
          ' params$cut_buffer, params$xlab_spacing)', collapse = ''), 
    '```\n', 
    "### Mismatches  \n", 
    paste('```r}, echo = F, results = "asis", ',
          'fig.width=25, message=F, warning=F}', collapse = ''), 
    paste('p <- amplican::plot_mismatches(alignments, config, {{ID}},',
          ' params$cut_buffer, params$xlab_spacing)', collapse = ''), 
    '```\n', 
    "### Cuts  \n", 
    paste('```r}, echo = F, results = "asis", ',
          'fig.width=25, message=F, warning=F}', collapse = ''), 
    paste('amplican::plot_cuts(alignments, config, {{ID}},',
          ' params$cut_buffer, params$xlab_spacing)', collapse = ''), 
    '```\n',
    "### Variants  \n", 
    paste('```r}, echo = F, message=F, results = "asis", ',
          'message=F, warning=F}', collapse = ''),
    paste('p <- amplican::plot_variants(alignments, config, {{ID}}, ',
          ' params$cut_buffer)', collapse = ''),
    '```\n'))
})
# knit the source
res = knitr::knit_child(text = src, quiet = TRUE)
cat(res, sep = '\n')