Counting RT mutations from illumina sequencing data.

1 Calculating error rates.

This document is intended to step through the process followed in these estimates of RT error rates. This R package is a companion to the similarly named ‘errrt’ Perl package, which is responsible for preprocessing the raw sequencing data.

2 Introduction

An outline of the experimental process is represented by:

Experimental Design

The errrt package is responsible for the bottom and right side of the image. It takes the raw fastq data, merges the reads using flash, writes out the index and identifier for all reads identical to the template, writes out the sequences for the ones not identical to the template, aligns them with fasta36 against the template, interprets the fasta36 output, and writes a table of the mismatches, insertions, deletions from it.

This package, creatively named ‘errRt’ takes a sample sheet describing the samples, reads the tables of identical reads/mismatches, and creates matrices of mutation rates for various classes of mutations.

3 Installation

The perl package may be installed via:

git clone https://github.com/abelew/errrt
cd errrt
perl Build.PL
perl Build install

The R package may be installed via:

devtools::install_github("abelew/errRt")

4 Preprocessing Steps

I split up the preprocessing into discrete steps so that I could check on each phase and try to ensure correctness. The errrt Perl package was written so that pretty much everything may be set via command line parameters, but all the defaults were written to suite this data set. Thus I have a pretty simple bash for loop which handled the work of preprocessing the data.

The following fragment of bash was run inside my directory preprocessing/, which contained 3 subdirectories: s1/, s2/, and s3/ which contain the raw reads for the 3 samples. The excel spreadsheet in inst/extdata/all_samples.xlsx describes them, but s1 is the control, s2 has low magnesium, and s3 has high.

#!/usr/bin/env bash

## The default options should take care of most of what these require.
start=$(pwd)
for sample in $(/bin/ls -d s*); do
    cd ${sample}
    ## Merge the forward and reverse reads.
    ../../errrt/script/merge_reads.pl --read1 *_1.fastq.gz --read2 *_2.fastq.gz
    ## Extract the reads with the template sequence into a fasta file.
    ## This assumes input named 'step1.extendedFrags.fastq.xz' and output named 'step2.fasta'
    ## This output is not compressed because it is easier to invoke fasta36 on uncompressed data.
    ## Thus it defaults to step2.fasta as output.
    ../../errrt/script/extract_fasta.pl
    ## Run fasta against the template sequence.
    ## This assumes step2.fasta as input and provides step3.txt.xz as output.
    ../../errrt/script/run_fasta.pl
    ## Boil the fasta output into something we can read.
    ## This assumes step3.txt.xz as input and step4.txt.xz as output.
    ../../errrt/script/parse_fasta.pl
    cd ${start}
done

As the above shell fragment states, the final outputs of errrt are step4.txt.xz and step2_identical_reads.txt.xz. These filenames are therefore recorded in the sample sheet in the columns ‘ident_table’ and ‘mutation_table’ for each sample.

I decided for the moment at least to copy these files into inst/extdata and modify the sample sheet in inst/extdata to refer to them.

5 Processing the data

I wrapped the various steps performed by errRt into a primary main function, create_matrices(). The following is the invocation I used for this data, along with some explanation of what is happening.

5.1 Getting Started

library(errRt)
basedir <- system.file("extdata", package="errRt")
metadata_file <- file.path(basedir, "all_samples.xlsx")
metadata_file

## [1] "/mnt/sshfs/major/local/R/3.6.1/3.10/errRt/extdata/all_samples.xlsx"

setwd(basedir)
meta <- hpgltools::extract_metadata(metadata_file)

Here is the resulting experimental design. The create_matrices() function will use this to read the tables of identical reads and mutated reads.

knitr::kable(meta)

	sampleid	condition	batch	identtable	mutationtable
s1	s1	control	b201910	s1/step2_identical_reads.txt.xz	s1/step4.txt.xz
s2	s2	low	b201910	s2/step2_identical_reads.txt.xz	s2/step4.txt.xz
s3	s3	high	b201910	s3/step2_identical_reads.txt.xz	s3/step4.txt.xz

5.2 Create matrices

The create_matrices() function has a few options to adjust the stringency of the various ways to define an error. It functions to count the various error types with quantify_parsed(), gather them into matrices by sample along with the perceived errors by the sequencer, and calculate aggregate error rates for the various classes of error. In other words, it does pretty much everything in one function call.

error_data <- create_matrices(sample_sheet=metadata_file, ident_column="identtable",
                              mut_column="mutationtable", min_reads=3, min_indexes=3,
                              min_sequencer=10, min_position=24, prune_n=TRUE,
                              verbose=TRUE)

## Starting sample: 1.

## Reading the file containing mutations: s1/step4.txt.xz

## Error: 's1/step4.txt.xz' does not exist in current working directory ('/mnt/sshfs/cbcbsub/fs/cbcb-lab/nelsayed/scratch/atb/dnaseq/rt_erates_destefano_2019/errRt/vignettes').

The options are explained in the help documentation, but here is an overview of the ones which might not be self-evident.

• min_reads: How many reads/index are required for a given index to be considered usable? This is the first toggle of stringency and is affected primarily by the sequencing depth.
• min_indexes: How many separate indexes are expected to agree on a given mutation in order for it to be considered real? This toggle of stringency is primarily affected by the success in getting a large number of different RT primers to bind and extend during the early stages of creating the sequencing libraries.
• min_sequencer: In order for a given mutation to be deemed ‘sequencer-specific’, 1 and only 1 read out of ‘n’ for a given index must have the mutation. This toggles ‘n’, higher numbers are thus more stringent and require greater sequencing depth. If this is set too low, a pool of false positives will result; the default of 10 is completely arbitrary.
• min_position: In the actual experiment, it appears that the primer used for the experimental samples may have had an error, or it was differently cut out of the PAGE purification gel, or something else happened which resulted in a tremendous percentage of the reads presenting themselves as either a deletion at ~ position 20, or a A to C mutation at position 22 or somesuch. This parameter just ignores any mutations occuring before this position. We observed that this strange trend stopped at position 22, and so set it to 24 to be safe.
• prune_n: When true, this simply drops entries in the mutation table which have Ns.

create_matrices() works in 3 main stages:

1. For each sample in the metadata, invoke quantify_parsed(). This function is responsible for reading the tables and collating the various error types. This is the most interesting portion of code and most likely to have problems.
2. Collect the various tables from step 1 and merge them into 1 matrix per table where the rows are mutation types and columns are samples.
3. Calculate error rates by read/index and write new versions of the matrices from step 2 accordingly.

While create_matrices() is running, it should print some information about what it is doing and give some idea about how much data is lost due to the stringency parameters. Its result is a list:

• samples: The result from quantify_parsed() for each sample in the metadata.
• reads_per_sample: The sum of the reads observed in the mutation data and the identity data, e.g. how many reads survived.
• indexes_per_sample: The sum of the indexes observed. This is our primary normalization factor, but along with reads_per_sample should help define a healthy vs. problematic sample.
• matrices: The ‘raw’ matrices for each table.
• matrices_by_counts: The matrices dividied by the indexes.
• normalized: The ‘raw’ matrices as cpm, thus normalized by library sizes.
• normalized_by_counts: The matrices as cpm and dividied by indexes.

Thus, as one might imagine, there are a lot of matrices returned, here are some examples of the data after normalization, in each case I will first print the data for the RT mutations followed by the sequencer:

1. Mismatches by position as a heatmap: This is quite curious, as we see apparently fewer variants at the end of the sequence.

gplots::heatmap.2(error_data[["normalized_by_counts"]][["miss_indexes_by_position"]],
                  trace="none",
                  Rowv=FALSE)

## Error in gplots::heatmap.2(error_data[["normalized_by_counts"]][["miss_indexes_by_position"]], : object 'error_data' not found

gplots::heatmap.2(error_data[["normalized_by_counts"]][["miss_sequencer_by_position"]],
                  trace="none",
                  Rowv=FALSE)

## Error in gplots::heatmap.2(error_data[["normalized_by_counts"]][["miss_sequencer_by_position"]], : object 'error_data' not found

2. Mismatches by reference nucleotide

knitr::kable(error_data[["normalized_by_counts"]][["miss_indexes_by_ref_nt"]])

## Error in knitr::kable(error_data[["normalized_by_counts"]][["miss_indexes_by_ref_nt"]]): object 'error_data' not found

knitr::kable(error_data[["normalized_by_counts"]][["miss_sequencer_by_ref_nt"]])

## Error in knitr::kable(error_data[["normalized_by_counts"]][["miss_sequencer_by_ref_nt"]]): object 'error_data' not found

3. Mismatches by mutated nucleotide

knitr::kable(error_data[["normalized_by_counts"]][["miss_indexes_by_hit_nt"]])

## Error in knitr::kable(error_data[["normalized_by_counts"]][["miss_indexes_by_hit_nt"]]): object 'error_data' not found

knitr::kable(error_data[["normalized_by_counts"]][["miss_sequencer_by_hit_nt"]])

## Error in knitr::kable(error_data[["normalized_by_counts"]][["miss_sequencer_by_hit_nt"]]): object 'error_data' not found

4. Mismatches by type: A to C, G to A etc etc

knitr::kable(error_data[["normalized_by_counts"]][["miss_indexes_by_type"]])

## Error in knitr::kable(error_data[["normalized_by_counts"]][["miss_indexes_by_type"]]): object 'error_data' not found

knitr::kable(error_data[["normalized_by_counts"]][["miss_sequencer_by_type"]])

## Error in knitr::kable(error_data[["normalized_by_counts"]][["miss_sequencer_by_type"]]): object 'error_data' not found

5. Mismatches by transition/transversion

knitr::kable(error_data[["normalized_by_counts"]][["miss_indexes_by_trans"]])

## Error in knitr::kable(error_data[["normalized_by_counts"]][["miss_indexes_by_trans"]]): object 'error_data' not found

knitr::kable(error_data[["normalized_by_counts"]][["miss_sequencer_by_trans"]])

## Error in knitr::kable(error_data[["normalized_by_counts"]][["miss_sequencer_by_trans"]]): object 'error_data' not found

6. Mismatches by strong/weak

knitr::kable(error_data[["normalized_by_counts"]][["miss_indexes_by_strength"]])

## Error in knitr::kable(error_data[["normalized_by_counts"]][["miss_indexes_by_strength"]]): object 'error_data' not found

knitr::kable(error_data[["normalized_by_counts"]][["miss_sequencer_by_strength"]])

## Error in knitr::kable(error_data[["normalized_by_counts"]][["miss_sequencer_by_strength"]]): object 'error_data' not found

7. Insertions by position

knitr::kable(error_data[["normalized_by_counts"]][["insert_indexes_by_position"]])

## Error in knitr::kable(error_data[["normalized_by_counts"]][["insert_indexes_by_position"]]): object 'error_data' not found

knitr::kable(error_data[["normalized_by_counts"]][["insert_sequencer_by_position"]])

## Error in knitr::kable(error_data[["normalized_by_counts"]][["insert_sequencer_by_position"]]): object 'error_data' not found

8. Insertions by nt

knitr::kable(error_data[["normalized_by_counts"]][["insert_indexes_by_nt"]])

## Error in knitr::kable(error_data[["normalized_by_counts"]][["insert_indexes_by_nt"]]): object 'error_data' not found

knitr::kable(error_data[["normalized_by_counts"]][["insert_sequencer_by_nt"]])

## Error in knitr::kable(error_data[["normalized_by_counts"]][["insert_sequencer_by_nt"]]): object 'error_data' not found

9. Deletions by position

There was insufficient data to quantify these.