In jonathonthill/sangerseqR: Tools for Sanger Sequencing Data in R

  BiocStyle::markdown()

library(knitr, quietly=TRUE)
library(sangerseqR, quietly=TRUE)
library(Biostrings, quietly=TRUE)
opts_chunk$set(tidy=TRUE, tidy.opts=list(width.cutoff=70))

Introduction

The sangerseqR package provides basic functions for importing and working with sanger sequencing data files. It currently functions with Scf and ABIF files. The Scf file specification is an open source, although somewhat limited, data file type. Several tools designed to view and or edit chromatogram data can convert file types to Scf. The ABIF file specification is a proprietary data storage file specification for sequencing data generated by Applied Biosystems machines. More information on each filetype can be found at the following sites:

• ABIF: https://projects.nfstc.org/workshops/resources/articles/ABIF_File_Format.pdf
• SCF: http://staden.sourceforge.net/manual/formats_unix_2.html

The objects and functions included in this package were developed as part of the Poly Peak Parser web application (http://yost.genetics.utah.edu/software.php), which automates the process of seperating ambiguous double peaks from Sanger sequencing individuals containing heterozygous indels. This package contains a complete working copy of the Poly Peak Parser web application that can be run locally using the PolyPeakParser() function. In addition to the web program, this package also provides general objects and functions for working with Sanger sequencing data.

This vignette will walk you through a typical workflow using two sequence files: 1) homozygous.scf and 2) heterozygous.ab1. As their names indicate, the first example contains results typical of sequencing from a PCR product of a homozygous individual or from a plasmid. The second example contains results from sequencing the same region in an individual with a small indel.

Loading Data

The first step of most workflows will be to upload data from a sequencing results file. This can be done using one of three included functions: read.abif(), read.scf() and readsangerseq(). The first two functions directly import all of the fields into abif and scf class objects, respectively. These classes are meant as intermediate classes and exist to allow the user to inspect the file contents, as file contents may vary between basecallers and sequencing machines. Users should generally use the readsangerseq() function. This function automatically detects and reads in the file type and then extracts the fields necessary to create a sangerseq class object, which is used by all of the other functions in this package.

read.abif

read.abif takes a single argument for the filename of the abif file to be read. The resulting object contains three major parts. The header, containing information on the file structure, the directory, containing information on each of the data fields included in the file, and the data fields. Here is an example:

hetab1 <- read.abif(system.file("extdata", 
                                "heterozygous.ab1", 
                                package="sangerseqR"))
str(hetab1, list.len=20)

As you can see, the file is very long and contains a lot of Data fields (130 in this example). However, most of these contain run information and only a few are directly relevant to data analysis:

| Data Field | Description | |:--------- |:--------- | |DATA.9-DATA.12 | Vectors containing the signal intensities for each channel. | |FWO.1 | A string containing the base corresponding to each channel. For example, if it is "ACGT", then DATA.9 = A, DATA.10 = C, DATA.11 = G and DATA.12 = T. | |PLOC.2 | Peak locations as an index of the trace vectors. | |PBAS.1, PBAS.2 | Primary basecalls. PBAS.1 may contain bases edited in the original basecaller, while PBAS.2 always contains the basecaller's calls. | |P1AM.1 | Amplitude of primary basecall peaks. | |P2BA.1 | (optional) Contains the secondary basecalls. | |P2AM.1 | (optional) Amplitude of the secondary basecall peaks. |

read.scf

Like read.abif, read.scf takes a single argument with the filename. However, the data structure of the resulting scf object is far less complicated, containing only a header with file structure information, a matrix of the trace data (@sample_points), a matrix of relative probabilities of each base at each position (@sequence_probs), basecall positions (@basecall_positions), basecalls (@basecalls) and optionally a comments sections with the run data (@comments). The last slot (@private) is rarely used and impossible to interpret without knowing how it was created.

homoscf <- read.scf(system.file("extdata", 
                                "homozygous.scf", 
                                package="sangerseqR"))
str(homoscf)

readsangerseq

The readsangerseq function is a convenience function equivalent to sangerseq(read.abif(file)) or sangerseq(read.scf(file)). It should generally be used when the contents of the file do not need to be directly accessed because it returns a sangerseq object, described below.

Sangerseq Class Objects

The sangerseq class is the backbone of the sangerseqR package and contains the chromatogram data necesary to perform all other functions. It can be created in two ways: from an abif or scf object using the sangerseq method or directly from an abif or scf file using readsangerseq.

#from a sequence file object
homosangerseq <- sangerseq(homoscf)

#directly from the file
hetsangerseq <- readsangerseq(system.file("extdata", 
                                          "heterozygous.ab1", 
                                          package="sangerseqR"))
str(hetsangerseq)

The slots are as follows:

| Slot | Description | |:--------- |:--------- | | primarySeqID | Identification of the primary Basecalls. | | primarySeq | The primary Basecalls formatted as a DNAString object. | | secondarySeqID | Identification of the secondary Basecalls. | | secondarySeq | The secondary Basecalls formatted as a DNAString object. | | traceMatrix | A numerical matrix containing 4 columns corresponding to the normalized signal values for the chromatogram traces. Column order = A,C,G,T. | | peakPosMatrix | A numerical matrix containing the position of the maximum peak values for each base within each Basecall window. If no peak was detected for a given base in a given window, then "NA". Column order = A,C,G,T. | | peakAmpMatrix | A numerical matrix containing the maximum peak amplitudes for each base within each Basecall window. If no peak was detected for a given base in a given window, then 0. Column order = A,C,G,T. |

Accessor functions also exist for each slot in the sangerseq object. Most of the accessors return the data in its native format, but the primarySeq and secondarySeq accessors can optionally return the data as a character string or a DNAString class object from the Biostrings package by setting string=TRUE or string=FALSE, respectively. The DNAString class contains several convenient functions for manipulating the sequence, including generating the reverse compliment and performing alignments. The Biostrings package is automatically loaded with the sangerseq package, so all methods should be available.

#default is to return a DNAString object
Seq1 <- primarySeq(homosangerseq)
reverseComplement(Seq1)

#can return as string
primarySeq(homosangerseq, string=TRUE)

Creating Chromatograms

Basic chromatogram plots can be made using the chromatogram function. These plots are optimized for printing, so they contain several rows to plot all of the data simultaneously. The downside of this is that it can give an error if the graphics device dimensions are not large enough. If this occurs, we suggest you provide a filename in the command to save it to a pdf automatically sized to fit everything. Several parameters can also be set to affect how the plot appears. These are documented in the chromatogram help file.

chromatogram(hetsangerseq, width=80, height=3, trim5=50, trim3=100, 
             showcalls='both')

Making Basecalls

As shown in the chromatogram, secondary basecalls are sometimes provided in ab1 files (Scf files are unable to show them). However, the exact nature of these calls is inconsistent. In the heterozygous.ab1 file used here, it is any peak near the primary peak, no matter how small. For example, base 100 (first base on the second line) has a primary call of "C" and a secondary call of "A", even though the A peak is very small and likely noise. In homozygous sequencing results, these calls should simply be ignored and are hidden in the chromatogram by default (showcalls="primary"). When heterozygous regions of the sequence are present, the makeBaseCalls can be used to determine whether a particular peak is homozygous or heterozygous and call the appropriate bases.

Let's use the chromatogram we created in the previous section as an example. The chromatogram contains a homozygous region from bases 1 to approximately 160, but then breaks down into a series of double peaks for the remainder of the chromatogram. This is due to an indel in one allele of the sequenced region. makeBaseCalls can be used to show this more clearly or to add the secondary basecalls if the data file does not contain them. The function essentially divides the sequence into a series of basecall windows and identifies the tallest peak for each fluorescence channel within the window. These peaks are converted to signal ratio to the tallest peak. A cutoff ratio is then applied to determine if a peak is signal or noise. Peaks below this ratio are ignored. Remaining peaks in each window are used to make primary and secondary basecalls.

hetcalls <- makeBaseCalls(hetsangerseq, ratio=0.33)
hetcalls

The resulting file now contains the maximum peak in the @primarySeq slot and the second tallest peak, if it is above the cutoff, in the @secondarySeq slot. If only one peak is above the cutoff ratio, then this call matches the primary basecall. If three peaks were above the cutoff ratio, then the peak with the maximum amplitude is the primary basecall and an ambiguous base code is used as the secondary basecall. The resulting chromatogram also shows this:

chromatogram(hetcalls, width=80, height=3, trim5=50, trim3=100, 
             showcalls='both')

Chromatogram of heterozygous sequencing results after making basecalls. Primary and secondary basecalls now match for homozygous peaks

Parsing Alleles

Although makeBaseCalls has fixed the primary and secondary peak calls. It still does not tell us anything about the nature of the mutation. For this, we need to set the allele phase using a reference base sequence from an online source or from another sequencing run on a homozygous sample. The examples used in this vignette are from heterozygous and homozygous siblings, so we will use the primary basecalls from the homozygous sibling (loaded earlier) as our reference. The beginnings and ends of these sequences do not need to match, but the reference should ideally encompass the sequenced region. setAllelePhase will then separate the primary and secondary basecalls into reference and non-reference bases at each position and set (@primarySeq) to the reference and @secondarySeq to the non-reference allele.

ref <- subseq(primarySeq(homosangerseq, string=TRUE), start=30, width=500)
hetseqalleles <- setAllelePhase(hetcalls, ref, trim5=50, trim3=300)
hetseqalleles

At this point, we could plot the chromatogram again, but it is more informative to align the resulting sequences to see how the alleles differ. Since sangerseqR depends on Biostrings, pairwiseAlignment can be used.

pa <- pairwiseAlignment(primarySeq(hetseqalleles)[1:400], 
                        secondarySeq(hetseqalleles)[1:400], 
                        type="global-local")
writePairwiseAlignments(pa)

Conclusion

In this vignette, we have walked you through the basic functions in the sangerseqR package. This work is a work in progress and we hope to improve its functionality. For example, improving the base calling algorithm and adding an interactive chromatogram function. If you have any suggestions or requested features, please email Jonathon Hill at [e-mail].

jonathonthill/sangerseqR documentation built on July 1, 2023, 4:55 p.m.