Package 'DECIPHER'

Description

DECIPHER is a software toolset that can be used for deciphering and managing biological sequences efficiently using the R statistical programming language. The program is designed to be used with non-destructive workflows for importing, maintaining, analyzing, manipulating, and exporting a massive amount of sequences.

Details

Index:

AA_REDUCED               Reduced amino acid alphabets
Add2DB                   Add Data to a Database
AdjustAlignment          Improve An Existing Alignment By Adjusting Gap
                         Placements
AlignDB                  Align Two Sets of Aligned Sequences in a Sequence
                         Database
AlignPairs               Align pairs of sequences
AlignProfiles            Align Two Sets of Aligned Sequences
AlignSeqs                Align a Set of Unaligned Sequences
AlignSynteny             Pairwise Aligns Syntenic Blocks
AlignTranslation         Align Sequences By Their Amino Acid Translation
AmplifyDNA               Simulate Amplification of DNA by PCR
Array2Matrix             Create a Matrix Representation of a Microarray
BLOSUM                   BLOSUM Amino Acid Substitution Matrices
BrowseDB                 View a Database Table in a Web Browser
BrowseSeqs               View Sequences in a Web Browser
CalculateEfficiencyArray Predict the Hybridization Efficiency of
                         Probe/Target Sequence Pairs
CalculateEfficiencyFISH  Predict Thermodynamic Parameters of Probe/Target
                         Sequence Pairs
CalculateEfficiencyPCR   Predict Amplification Efficiency of Primer Sequences
Clusterize               Cluster Sequences By Distance
Codec                    Compression/Decompression of Character Vectors
ConsensusSequence        Create a Consensus Sequence
Cophenetic               Compute cophenetic distances on dendrogram objects
CorrectFrameshifts       Corrects Frameshift Errors In Protein Coding
                         Sequences
CreateChimeras           Create Artificial Chimeras
DB2Seqs                  Export Database Sequences to a FASTA or FASTQ File
deltaGrules              Free Energy of Hybridization of Probe/Target
                         Quadruplets
deltaHrules              Change in Enthalpy of Hybridization of DNA/DNA
                         Quadruplets in Solution
deltaHrulesRNA           Change in Enthalpy of Hybridization of RNA/RNA
                         Quadruplets in Solution
deltaSrules              Change in Entropy of Hybridization of DNA/DNA
                         Quadruplets in Solution
deltaSrulesRNA           Change in Entropy of Hybridization of RNA/RNA
                         Quadruplets in Solution
DesignArray              Design a Set of DNA Microarray Probes for Detecting
                         Sequences
DesignPrimers            Design Primers Targeting a Specific Group of
                         Sequences
DesignProbes             Design FISH Probes Targeting a Specific Group of
                         Sequences
DesignSignatures         Design PCR Primers for Amplifying Group-Specific
                         Signatures
DetectRepeats            Detect Repeats in a Sequence
DigestDNA                Simulate Restriction Digestion of DNA
Disambiguate             Expand Ambiguities into All Permutations of a
                         DNAStringSet
DistanceMatrix           Calculate the Distance Between Sequences
ExtractGenes             Extract Predicted Genes from a Genome
FindChimeras             Find Chimeras in a Sequence Database
FindGenes                Find Genes in a Genome
FindNonCoding            Find Non-Coding RNAs in a Genome
FindSynteny              Finds Synteny in a Sequence Database
FormGroups               Forms Groups By Rank
Genes-class              Genes objects and accessors
HEC_MI                   Mutual Information for Protein Secondary Structure
                         Prediction
IdConsensus              Create Consensus Sequences by Groups
IdentifyByRank           Identify By Taxonomic Rank
IdLengths                Determine the Number of Characters in Each Sequence
                         of Each Sequence
IdTaxa                   Assign Sequences a Taxonomic Classification
IndexSeqs                Build an inverted index
InvertedIndex-class      InvertedIndex objects
LearnNonCoding           Learn a Non-Coding RNA Model
LearnTaxa                Train a Classifier for Assigning Taxonomy
MapCharacters            Map Changes in Ancestral Character States
MaskAlignment            Mask Highly Variable Regions of An Alignment
MeltDNA                  Simulate Melting of DNA
MIQS                     MIQS Amino Acid Substitution Matrix
MMLSUM                   MMLSUM Amino Acid Substitution Matrices
MODELS                   Available Models of Sequence Evolution
NNLS                     Sequential Coordinate-wise Algorithm for the
                         Non-negative Least Squares Problem
NonCoding                NonCoding Models for Common Non-Coding RNA Families
NonCoding-class          NonCoding Objects and Methods
OrientNucleotides        Orient Nucleotide Sequences
PAM                      PAM Amino Acid Substitution Matrices
PFASUM                   PFASUM Amino Acid Substitution Matrices
PredictDBN               Predict RNA Secondary Structure in Dot-Bracket
                         Notation
PredictHEC               Predict Protein Secondary Structure
Read Dendrogram          Read a Dendrogram from a Newick Formatted File
RemoveGaps               Remove Gap Characters in Sequences
RESTRICTION_ENZYMES      Common Restriction Enzyme's Cut Sites
ScoreAlignment           Score a Multiple Sequence Alignment
SearchDB                 Obtain Specific Sequences from a Database
SearchIndex              Search an inverted index
Seqs2DB                  Add Sequences from Text File to Database
StaggerAlignment         Produce a Staggered Alignment
Synteny-class            Synteny blocks and hits
Taxa-class               Taxa training and testing objects
TerminalChar             Determine the Number of Terminal Characters
TileSeqs                 Form a Set of Tiles for Each Group of Sequences
TrainingSet_16S          Training Set for Classification of 16S rRNA Gene
                         Sequences
Treeline                 Construct a Phylogenetic Tree
TrimDNA                  Trims DNA Sequences to the High Quality Region
                         Between Patterns
WriteDendrogram          Write a Dendrogram to Newick Format
WriteGenes               Write Genes to a File

Author(s)

Erik Wright

Maintainer: Erik Wright <[email protected]>

Description

The AA_REDUCED list contains reductions of the standard amino acid alphabet (AA_STANDARD).

Usage

AA_REDUCED

Details

Reduced amino alphabets can sometimes improve sensitivity and specificity of finding homologous matches between amino acid sequences. The first 12 AA_REDUCED alphabets were optimized for finding synteny between genomic sequences. The next 113 alphabets are from a review of published amino acid alphabets (Solis, 2015). The following 17 alphabets were optimized for amino acid classification. The subsequent 18 alphabets are progressive mergers based on average similarities in PFASUM. The following 25 alphabets were optimized for protein search. The final 8 alphabets were optimized for clustering.

References

Solis, A. (2015). Amino acid alphabet reduction preserves fold information contained in contact interactions in proteins. Proteins: Structure, Function, and Genetics, 83(12), 2198-2216.

See Also

FindSynteny, LearnTaxa, PFASUM

Examples

str(AA_REDUCED)
AA_REDUCED[[1]]

Description

Adds a data.frame to a database table by its row.names.

Usage

Add2DB(myData,
       dbFile,
       tblName = "Seqs",
       clause = "",
       verbose = TRUE)

Arguments

Details

Data contained in myData will be added to the tblName by its respective row.names.

Value

Returns TRUE if the data was added successfully, or FALSE otherwise.

Author(s)

Erik Wright [email protected]

References

ES Wright (2016) "Using DECIPHER v2.0 to Analyze Big Biological Sequence Data in R". The R Journal, 8(1), 352-359.

See Also

Seqs2DB, SearchDB, BrowseDB

Examples

if (require("RSQLite", quietly=TRUE)) {
	# Create a sequence database
	gen <- system.file("extdata", "Bacteria_175seqs.gen", package="DECIPHER")
	dbConn <- dbConnect(dbDriver("SQLite"), ":memory:")
	Seqs2DB(gen, "GenBank", dbConn, "Bacteria")
	
	# Identify the sequence lengths
	l <- IdLengths(dbConn)
	
	# Add lengths to the database
	Add2DB(l, dbConn)
	
	# View the added lengths
	BrowseDB(dbConn)
	
	# Change the value of existing columns
	ids <- data.frame(identifier=rep("Bacteroidetes", 18), stringsAsFactors=FALSE)
	rownames(ids) <- 10:27
	Add2DB(ids, dbConn)
	BrowseDB(dbConn)
	
	# Add data to a subset of rows using a clause
	ids[[1]][] <- "Changed"
	nrow(ids) # 18 rows
	Add2DB(ids, dbConn, clause="accession like 'EU808318%'")
	BrowseDB(dbConn) # only 1 row effected
	
	dbDisconnect(dbConn)
}

Description

Makes small adjustments by shifting groups of gaps left and right to find their optimal positioning in a multiple sequence alignment.

Usage

AdjustAlignment(myXStringSet,
                perfectMatch = 2,
                misMatch = -1,
                gapLetter = -3,
                gapOpening = -0.1,
                gapExtension = 0,
                substitutionMatrix = NULL,
                shiftPenalty = -0.2,
                threshold = 0.1,
                weight = 1,
                processors = 1)

Arguments

Details

The process of multiple sequence alignment often results in the integration of small imperfections into the final alignment. Some of these errors are obvious by-eye, which encourages manual refinement of automatically generated alignments. However, the manual refinement process is inherently subjective and time consuming. AdjustAlignment refines an existing alignment in a process similar to that which might be applied manually, but in a repeatable and must faster fashion. This function shifts all of the gaps in an alignment to the left and right to find their optimal positioning. The optimal position is defined as the position that maximizes the alignment “score”, which is determined by the input parameters. The resulting alignment will be similar to the input alignment but with many imperfections eliminated. Note that the affine gap penalties here are different from the more flexible penalties used in AlignProfiles, and have been optimized independently.

Value

An XStringSet of aligned sequences.

Author(s)

Erik Wright [email protected]

References

Wright, E. S. (2015). DECIPHER: harnessing local sequence context to improve protein multiple sequence alignment. BMC Bioinformatics, 16, 322. http://doi.org/10.1186/s12859-015-0749-z

Wright, E. S. (2020). RNAconTest: comparing tools for noncoding RNA multiple sequence alignment based on structural consistency. RNA 2020, 26, 531-540.

See Also

AlignSeqs, AlignTranslation, PFASUM, StaggerAlignment

Examples

# a trivial example
aa <- AAStringSet(c("ARN-PK", "ARRP-K"))
aa # input alignment
AdjustAlignment(aa) # output alignment

# specifying an alternative substitution matrix
AdjustAlignment(aa, substitutionMatrix="BLOSUM62")

# a real example
fas <- system.file("extdata", "Streptomyces_ITS_aligned.fas", package="DECIPHER")
dna <- readDNAStringSet(fas)
dna # input alignment
adjustedDNA <- AdjustAlignment(dna) # output alignment
BrowseSeqs(adjustedDNA, highlight=1)
adjustedDNA==dna # most sequences were adjusted (those marked FALSE)

Description

Merges the two separate sequence alignments in a database. The aligned sequences must have separate identifiers in the same table or be located in different database tables.

Usage

AlignDB(dbFile,
        tblName = "Seqs",
        identifier = "",
        type = "DNAStringSet",
        add2tbl = "Seqs",
        batchSize = 10000,
        perfectMatch = 2,
        misMatch = -1,
        gapOpening = -12,
        gapExtension = -3,
        gapPower = -1,
        terminalGap = -4,
        normPower = c(1, 0),
        standardize = TRUE,
        substitutionMatrix = NULL,
        processors = 1,
        verbose = TRUE,
        ...)

Arguments

Details

Sometimes it is useful to align two large sets of sequences, where each set of sequences is already aligned but the two sets are not aligned to each other. AlignDB first builds a profile of each sequence set in increments of batchSize so that the entire sequence set is not required to fit in memory. Next the two profiles are aligned using dynamic programming. Finally, the new alignment is applied to all the sequences as they are incrementally added to the add2tbl.

Two identifiers or tblNames must be provided, indicating the two sets of sequences to align. The sequences corresponding to the first identifier and tblName will be aligned to those of the second identifier or tblName. The aligned sequences are added to add2tbl under a new identifier formed from the concatenation of the two identifiers or tblNames. (See examples section below.)

Value

Returns the number of newly aligned sequences added to the database.

Author(s)

Erik Wright [email protected]

References

Wright, E. S. (2015). DECIPHER: harnessing local sequence context to improve protein multiple sequence alignment. BMC Bioinformatics, 16, 322. http://doi.org/10.1186/s12859-015-0749-z

Wright, E. S. (2020). RNAconTest: comparing tools for noncoding RNA multiple sequence alignment based on structural consistency. RNA 2020, 26, 531-540.

See Also

AlignProfiles, AlignSeqs, AlignTranslation, PFASUM

Examples

if (require("RSQLite", quietly=TRUE)) {
	gen <- system.file("extdata", "Bacteria_175seqs.gen", package="DECIPHER")
	fas <- system.file("extdata", "Bacteria_175seqs.fas", package="DECIPHER")
	
	# Align two tables and place result into a third
	dbConn <- dbConnect(dbDriver("SQLite"), ":memory:")
	Seqs2DB(gen, "GenBank", dbConn, "Seqs1", tblName="Set1")
	Seqs2DB(fas, "FASTA", dbConn, "Seqs2", tblName="Set2")
	AlignDB(dbConn, tblName=c("Set1", "Set2"), add2tbl="AlignedSets")
	l <- IdLengths(dbConn, "AlignedSets", add2tbl=TRUE)
	BrowseDB(dbConn, tblName="AlignedSets") # all sequences have the same width
	dbDisconnect(dbConn)
	
	# Align two identifiers and place the result in the same table
	dbConn <- dbConnect(dbDriver("SQLite"), ":memory:")
	Seqs2DB(gen, "GenBank", dbConn, "Seqs1")
	Seqs2DB(fas, "FASTA", dbConn, "Seqs2")
	AlignDB(dbConn, identifier=c("Seqs1", "Seqs2"))
	l <- IdLengths(dbConn, add2tbl=TRUE)
	BrowseDB(dbConn) # note the sequences with a new identifier
	dbDisconnect(dbConn)
}

Description

Aligns pairs of sequences globally or locally using anchored adaptive banding.

Usage

AlignPairs(pattern,
           subject,
           pairs = NULL,
           type = "values",
           perfectMatch = 2,
           misMatch = -1,
           gapOpening = -16,
           gapExtension = -1.2,
           substitutionMatrix = NULL,
           bandWidth = 50,
           dropScore = -100,
           processors = 1,
           verbose = TRUE)

Arguments

Details

Performs pairwise alignment of pattern and subject sequences, either in their respective pairs or as specified by pairs. Uses adaptive banding and (optionally) anchoring to accelerate the alignment. Unlike AlignProfiles, AlignPairs can perform local alignment around anchor positions to align local regions of the pattern and/or subject sequences. In the absence of anchors, AlignPairs performs global alignment without terminal gap penalties or with terminal gap penalties when virtual anchors are provided immediately outside the bounds of the sequences. (See examples section below.)

AlignPairs is designed to maximize speed, and provides slightly less accuracy than using AlignProfiles for pairwise alignment. Adaptive banding is applied with a fixed bandWidth to reduce memory consumption, rather than the dynamic band width used by AlignProfiles. There are a few other differences: AlignPairs applies affine rather than Zipfian gap penalties, there is no option to incorporate secondary structures, scores are not standardized by length, gap penalties are not modulated around specific residues, and any anchors must be supplied in pairs rather than determined automatically. For very dissimilar sequences, it is preferable to use AlignTranslation (best for coding sequences), AlignSeqs (best for nucleotide/protein sequences), or AlignProfiles.

AlignPairs does not directly output the pairwise alignments. Instead, it outputs statistics about the alignment and the position(s) of gaps in each sequence. This makes alignment more efficient because no sequences are copied. For many applications only the percent identity or number of gaps is needed, which can be calculated directly from the returned data.frame. However, the aligned sequences can also easily be obtained from the output if desired. (See examples section below.)

Value

If type is "values" (the default), a data.frame is returned with one alignment per input pattern or row of pairs if not NULL. Columns are defined as the sequences' index in pattern (Pattern), start position in the pattern sequence (PatternStart), end position in the pattern sequence (PatternEnd), index in the subject sequence (Subject), start position in the subject sequence (SubjectStart), end position in the subject sequence (SubjectEnd), number of matching positions in the alignment (Matches), number of mismatched positions in the alignment (Mismatches), total number of positions in the alignment (AlignmentLength), alignment score (Score), and the position/length of gaps in the pattern and subject (i.e., PatternGapPosition, PatternGapLength, SubjectGapPosition, & SubjectGapLength).

If type is "sequences", a list containing two components: the aligned pattern and subject.

If type is "both", a list with three components: the values, aligned pattern, and aligned subject.

Author(s)

Erik Wright [email protected]

See Also

AlignProfiles, IndexSeqs, SearchIndex

Examples

# import target sequences and build an inverted index
fas <- system.file("extdata", "PlanctobacteriaNamedGenes.fas.gz", package="DECIPHER")
target <- readAAStringSet(fas)
index <- IndexSeqs(target, K=6L)
index

# import query sequences and search the index
fas <- system.file("extdata", "50S_ribosomal_protein_L2.fas", package="DECIPHER")
dna <- readDNAStringSet(fas)
query <- translate(dna)
hits <- SearchIndex(query, index, scoreOnly=FALSE)
head(hits)
nrow(hits) # number of query/target pairs

# locally pairwise align the query/target pairs representing each hit
aligned <- AlignPairs(query, target, hits) # local alignment around hits
head(aligned)

# plot the percent identity of each alignment versus score
PIDs <- aligned$Matches/aligned$AlignmentLength
plot(PIDs, hits$Score) # versus the hit score
plot(PIDs, aligned$Score) # versus the alignment score

# plot the number of pattern gaps versus subject gaps per alignment
subjectGaps <- sapply(aligned$SubjectGapLength, sum)
patternGaps <- sapply(aligned$PatternGapLength, sum)
plot(subjectGaps, patternGaps, pch=16, col="#00000011")
abline(a=0, b=1) # identity line (y = x)

# extract the pairwise aligned regions
alignments <- AlignPairs(query, target, hits, type="sequences")
BrowseSeqs(c(alignments[[1]][1], alignments[[2]][1])) # view the first pair

# perform global pairwise alignment by creating virtual endpoint anchors
# virtual anchors are immediately out-of-bounds (positions 0 and width + 1)
# this causes gap opening/extension penalties to be applied at each terminus
hits$Position <- mapply(function(x, y, z)
		cbind(matrix(0L, 4), x, matrix(c(y, y, z, z), 4)),
	hits$Position,
	width(query)[hits$Pattern] + 1L,
	width(target)[hits$Subject] + 1L,
	SIMPLIFY=FALSE)
aligned <- AlignPairs(query, target, hits) # penalizes terminal gaps
head(aligned) # all alignments now span start to end of the sequences

# perform global pairwise alignment of sequences without anchors (approach 1)
pattern <- query[rep(1, length(query))] # first sequence repeated
subject <- query
aligned1 <- AlignPairs(pattern, subject) # no terminal gap penalties
head(aligned1)
# perform global pairwise alignment of sequences without anchors (approach 2)
pairs <- data.frame(Pattern=1L, Subject=seq_along(query))
aligned2 <- AlignPairs(query, query, pairs) # no terminal gap penalties
head(aligned2) # note the Pattern column is always 1 (first sequence)

Description

Aligns two sets of one or more aligned sequences by first generating representative profiles, then aligning the profiles with dynamic programming, and finally merging the two aligned sequence sets.

Usage

AlignProfiles(pattern,
              subject,
              p.weight = 1,
              s.weight = 1,
              p.struct = NULL,
              s.struct = NULL,
              perfectMatch = 2,
              misMatch = -1,
              gapOpening = -12,
              gapExtension = -3,
              gapPower = -1,
              terminalGap = -4,
              restrict = c(-1000, 2, 10),
              anchor = 0.7,
              normPower = c(1, 0),
              standardize = TRUE,
              substitutionMatrix = NULL,
              structureMatrix = NULL,
              processors = 1)

Arguments

Details

Profiles are aligned using dynamic programming, a variation of the Needleman-Wunsch algorithm for global alignment. The dynamic programming method requires order N*M time and memory space where N and M are the width of the pattern and subject. This method works by filling in a matrix of the possible “alignment space” by considering all matches, insertions, and deletions between two sequence profiles. The highest scoring alignment is then used to add gaps to each of the input sequence sets.

Heuristics can be useful to improve performance on long input sequences. The restrict parameter can be used to dynamically constrain the possible “alignment space” to only paths that will likely include the final alignment, which in the best case can improve the speed from quadratic time to nearly linear time. The degree of restriction is important, and the default value of restrict is reasonable in the vast majority of cases. It is also possible to prevent restriction by setting restrict to such extreme values that these requirements will never be met (e.g., c(-1e10, 1e10, 1e10)).

The argument anchor can be used to split the global alignment into multiple sub-alignments. This can greatly decrease the memory requirement for long sequences when appropriate anchor points can be found. Anchors are 15-mer (for DNA/RNA) or 7-mer (for AA) subsequences that are shared between at least anchor fraction of pattern(s) and subject(s). Anchored ranges are extended along the length of each sequence in a manner designed to split the alignment into sub-alignments that can be separately solved. For most input sequences, the default anchoring has no effect on accuracy, but anchoring can be disabled by setting anchor=NA.

Alternatively, anchor can be a matrix with 4 rows and one column per anchor. The first two rows correspond to the anchor start and end positions in the pattern sequence(s), and the second two rows are the equivalent anchor region in the subject sequence(s). Anchors specified in this manner must be strictly increasing (non-overlapping) in both sequences, and have an anchor width of at least two positions. Note that the anchors do not have to be equal length, in which case the anchor regions will also be aligned. Manually splitting the alignment into more subtasks can sometimes make it more efficient, but typically automatic anchoring is effective.

The argument normPower determines how the distribution of information is treated during alignment. Higher values of normPower encourage alignment between columns with higher occupancy (1 - fraction of gaps), and de-emphasize the alignment of columns containing many gaps. A normPower of 0 will treat all columns equally regardless of occupancy, which can be useful when the pattern or subject contain many incomplete (fragment) sequences. For example, normPower should be set to 0 when aligning many short reads to a longer reference sequence.

The arguments p.struct and s.struct may be used to provide secondary structure probabilities in the form of a list containing matrices or a single matrix. If the input is a list, then each list element must contain a matrix with dimensions q*w, where q is the number of possible secondary structure states, and w is the width of the unaligned pattern sequence. Values in each matrix represent the probability of the given state at that position in the sequence. Alternatively, a single matrix can be used as input if w is the width of the entire pattern or subject alignment. A structureMatrix must be supplied along with the structures. The functions PredictHEC and PredictDBN can be used to predict secondary structure probabilities in the format required by AlignProfiles (for amino acid and RNA sequences, respectively).

The gap cost function is based on the observation that gap lengths are best approximated by a Zipfian distribution (Chang & Benner, 2004). The cost of inserting a gap of length L is equal to: gapOpening + gapExtension*sum(seq_len(L - 1)^gapPower) when L > 1, and gapOpen when L = 1. This function effectively penalizes shorter gaps significantly more than longer gaps when gapPower < 0, and is equivalent to the affine gap penalty when gapPower is 0.

Value

An XStringSet of aligned sequences.

Author(s)

Erik Wright [email protected]

References

Chang, M. S. S., & Benner, S. A. (2004). Empirical Analysis of Protein Insertions and Deletions Determining Parameters for the Correct Placement of Gaps in Protein Sequence Alignments. Journal of Molecular Biology, 341(2), 617-631.

Needleman, S., & Wunsch, C. (1970). A general method applicable to the search for similarities in the amino acid sequence of two proteins. Journal of Molecular Biology, 48(3), 443-453.

Wright, E. S. (2015). DECIPHER: harnessing local sequence context to improve protein multiple sequence alignment. BMC Bioinformatics, 16, 322. http://doi.org/10.1186/s12859-015-0749-z

Wright, E. S. (2020). RNAconTest: comparing tools for noncoding RNA multiple sequence alignment based on structural consistency. RNA 2020, 26, 531-540.

Yu, Y.-K., et al. (2015). Log-odds sequence logos. Bioinformatics, 31(3), 324-331. http://doi.org/10.1093/bioinformatics/btu634

See Also

AlignDB, AlignSeqs, AlignSynteny, AlignTranslation, PFASUM, MIQS

Examples

# align two sets of sequences
fas <- system.file("extdata", "Bacteria_175seqs.fas", package="DECIPHER")
dna1 <- readDNAStringSet(fas, n=100) # the first 100 sequences
dna2 <- readDNAStringSet(fas, n=100, skip=100) # the rest
dna1 <- RemoveGaps(dna1, "common")
dna2 <- RemoveGaps(dna2, "common")
alignedDNA <- AlignProfiles(dna1, dna2)
BrowseSeqs(alignedDNA, highlight=1)

# specify a DNA substitution matrix
subMatrix <- matrix(0,
                    nrow=4, ncol=4,
                    dimnames=list(DNA_BASES, DNA_BASES))
diag(subMatrix) <- 5 # perfectMatch
alignedDNA.defaultSubM <- AlignProfiles(dna1, dna2, substitutionMatrix=subMatrix)
all(alignedDNA.defaultSubM==alignedDNA)

# specify a different DNA substitution matrix
subMatrix2 <- matrix(c(12, 3, 5, 3, 3, 12, 3, 6, 5, 3, 11, 3, 3, 6, 3, 9),
                    nrow=4, ncol=4,
                    dimnames=list(DNA_BASES, DNA_BASES))
alignedDNA.alterSubM <- AlignProfiles(dna1, dna2, substitutionMatrix=subMatrix2)
all(alignedDNA.alterSubM==alignedDNA)

# anchors are found automatically by default, but it is also
# possible to specify anchor regions between the sequences
anchors <- matrix(c(774, 788, 752, 766), nrow=4)
anchors
subseq(dna1, anchors[1, 1], anchors[2, 1])
subseq(dna2, anchors[3, 1], anchors[4, 1])
alignedDNA2 <- AlignProfiles(dna1, dna2, anchor=anchors)

Description

Performs profile-to-profile alignment of multiple unaligned sequences following a guide tree.

Usage

AlignSeqs(myXStringSet,
         guideTree = NULL,
         iterations = 2,
         refinements = 1,
         gapOpening = c(-18, -10),
         gapExtension = -3,
         useStructures = TRUE,
         structures = NULL,
         FUN = AdjustAlignment,
         levels = c(0.9, 0.7, 0.7, 0.4, 10, 5, 5, 2),
         alphabet = AA_REDUCED[[1]],
         processors = 1,
         verbose = TRUE,
         ...)

Arguments

Details

The profile-to-profile method aligns a sequence set by merging profiles along a guide tree until all the input sequences are aligned. This process has three main steps: (1) If guideTree=NULL, an initial single-linkage guide tree is constructed based on a distance matrix of shared k-mers. Alternatively, a dendrogram can be provided as the initial guideTree. (2) If iterations is greater than zero, then a UPGMA guide tree is built based on the initial alignment and the sequences are re-aligned along this tree. This process repeated iterations times or until convergence. (3) If refinements is greater than zero, then subsets of the alignment are re-aligned to the remainder of the alignment. This process generates two alignments, the best of which is chosen based on its sum-of-pairs score. This refinement process is repeated refinements times, or until convergence.

The purpose of levels is to speed-up the alignment process by not running time consuming processes when they are unlikely to change the outcome. The first four levels control when refinements occur and the function FUN is run on the alignment. The default levels specify that these events should happen when above 0.9 (AA; levels[1]) or 0.7 (DNA/RNA; levels[3]) average dissimilarity on the initial tree, when above 0.7 (AA; levels[2]) or 0.4 (DNA/RNA; levels[4]) average dissimilarity on the iterative tree(s), and after every tenth improvement made during refinement. The sixth element of levels (levels[6]) prevents FUN from being applied at any point to less than 5 sequences.

The FUN function is always applied before returning the alignment so long as there are at least levels[6] sequences. The default FUN is AdjustAlignment, but FUN can be any function that takes in an XStringSet as its first argument, as well as weights, processors, and substitutionMatrix as optional arguments. For example, the default FUN could be altered to not perform any changes by setting it equal to function(x, ...) return(x), where x is an XStringSet.

Secondary structures are automatically computed for amino acid and RNA sequences unless structures are provided or useStructures is FALSE. Use of structures generally provides a moderate improvement in average accuracy on difficult-to-align sequences. The default structureMatrix is used unless an alternative is provided. For RNA sequences, consensus secondary structures are only computed when the total length of the initial guide tree is at least 5 (levels[7]) or the length of subsequent trees is at least 2 (levels[8]). Note that input RNAStringSets are assumed to be structured non-coding RNAs. Largely unstructured RNAs should be aligned with useStructures set to FALSE or, ideally, aligned with AlignTranslation if coding sequences (i.e., mRNAs).

Value

An XStringSet of aligned sequences.

Author(s)

Erik Wright [email protected]

References

Wright, E. S. (2015). DECIPHER: harnessing local sequence context to improve protein multiple sequence alignment. BMC Bioinformatics, 16, 322. http://doi.org/10.1186/s12859-015-0749-z

Wright, E. S. (2020). RNAconTest: comparing tools for noncoding RNA multiple sequence alignment based on structural consistency. RNA 2020, 26, 531-540.

See Also

AdjustAlignment, AlignDB, AlignProfiles, AlignSynteny, AlignTranslation, ReadDendrogram, Treeline, StaggerAlignment

Examples

fas <- system.file("extdata", "Bacteria_175seqs.fas", package="DECIPHER")
dna <- readDNAStringSet(fas)
dna <- RemoveGaps(dna)
alignedDNA <- AlignSeqs(dna)
BrowseSeqs(alignedDNA, highlight=1)

# use secondary structure with RNA sequences
alignedRNA <- AlignSeqs(RNAStringSet(dna))
BrowseSeqs(alignedRNA, highlight=1)

Description

Performs pairwise alignment of all blocks of synteny between sets of sequences.

Usage

AlignSynteny(synteny,
             dbFile,
             tblName = "Seqs",
             identifier = "",
             processors = 1,
             verbose = TRUE,
             ...)

Arguments

Details

AlignSynteny will extract all sequence regions belonging to syntenic blocks in synteny, and perform pairwise alignment with AlignProfiles. Hits are used to anchor the alignment such that only the regions between anchors are aligned.

Value

A list with elements for each pair of identifiers in synteny. Each list element contains a DNAStringSetList one pairwise alignment per syntenic block.

Author(s)

Erik Wright [email protected]

See Also

FindSynteny, Synteny-class

Examples

if (require("RSQLite", quietly=TRUE)) {
	db <- system.file("extdata", "Influenza.sqlite", package="DECIPHER")
	synteny <- FindSynteny(db, minScore=50)
	DNA <- AlignSynteny(synteny, db)
	names(DNA)
	DNA[[1]] # the first set of pairwise alignments
	DNA[[1]][[1]] # the first block of synteny between H9N2 & H5N1
	unlist(DNA[[2]]) # a DNAStringSet of synteny between H9N2 & H2N2
}

Description

Performs alignment of a set of DNA or RNA sequences by aligning their corresponding amino acid sequences.

Usage

AlignTranslation(myXStringSet,
                 sense = "+",
                 direction = "5' to 3'",
                 readingFrame = NA,
                 type = class(myXStringSet),
                 geneticCode = GENETIC_CODE,
                 ...)

Arguments

Details

Alignment of proteins is often more accurate than alignment of their coding nucleic acid sequences. This function aligns the input nucleic acid sequences via aligning their translated amino acid sequences. First, the input sequences are translated according to the specified sense, direction, and readingFrame. The resulting amino acid sequences are aligned using AlignSeqs, and this alignment is (conceptually) reverse translated into the original sequence type, sense, and direction. Not only is alignment of protein sequences generally more accurate, but aligning translations will ensure that the reading frame is maintained in the nucleotide sequences.

If the readingFrame is NA (the default) then an attempt is made to guess the reading frame of each sequence based on the number of stop codons in the translated amino acids. For each sequence, the first reading frame will be chosen (either 1, 2, or 3) without stop codons, except in the final position. If the number of stop codons is inconclusive for a sequence then the reading frame will default to 1. The entire length of each sequence is translated in spite of any stop codons identified. Note that this method is only constructive in circumstances where there is a substantially long coding sequence with at most a single stop codon expected in the final position, and therefore it is preferable to specify the reading frame of each sequence if it is known.

Value

An XStringSet of the class specified by type, or a list of two components (nucleotides and amino acids) if type is "both". Note that incomplete starting and ending codons will be translated into the mask character ("+") if the result includes an AAStringSet.

Author(s)

Erik Wright [email protected]

References

Wright, E. S. (2015). DECIPHER: harnessing local sequence context to improve protein multiple sequence alignment. BMC Bioinformatics, 16, 322. http://doi.org/10.1186/s12859-015-0749-z

See Also

AlignDB, AlignProfiles, AlignSeqs, AlignSynteny, CorrectFrameshifts

Examples

# first three sequences translate to  MGFITP*
# and the last sequence translates as MGF-TP*
rna <- RNAStringSet(c("AUGGGUUUCAUCACCCCCUAA", "AUGGGCUUCAUAACUCCUUGA",
	"AUGGGAUUCAUUACACCGUAG", "AUGGGGUUUACCCCAUAA"))
RNA <- AlignSeqs(rna, verbose=FALSE)
RNA

RNA <- AlignTranslation(rna, verbose=FALSE)
RNA

AA <- AlignTranslation(rna, type="AAStringSet", verbose=FALSE)
AA

BOTH <- AlignTranslation(rna, type="both", verbose=FALSE)
BOTH

# example of aligning many protein coding sequences:
fas <- system.file("extdata", "50S_ribosomal_protein_L2.fas", package="DECIPHER")
dna <- readDNAStringSet(fas)
DNA <- AlignTranslation(dna) # align the translation then reverse translate
DNA

# using a mixture of standard and non-standard genetic codes
gC1 <- getGeneticCode(id_or_name2="1", full.search=FALSE, as.data.frame=FALSE)
# Mollicutes use an alternative genetic code
gC2 <- getGeneticCode(id_or_name2="4", full.search=FALSE, as.data.frame=FALSE)
w <- grep("Mycoplasma|Ureaplasma", names(dna))
gC <- vector("list", length(dna))
gC[-w] <- list(gC1)
gC[w] <- list(gC2)
AA <- AlignTranslation(dna, geneticCode=gC, type="AAStringSet")
BrowseSeqs(AA)

Description

Predicts the amplification efficiency of theoretical PCR products (amplicons) given one or more primer sequences.

Usage

AmplifyDNA(primers,
        myDNAStringSet,
        maxProductSize,
        annealingTemp,
        P,
        ions = 0.2,
        includePrimers=TRUE,
        minEfficiency = 0.001,
        ...)

Arguments

Details

Exponential amplification in PCR requires the annealing and elongation of two primers from target sites on opposing strands of the template DNA. If the template DNA sequence (e.g., chromosome) is known then predictions of theoretical amplicons can be obtained from in silico simulations of amplification. AmplifyDNA first searches for primer target sites on the template DNA, and then calculates an amplification efficiency from each target site using CalculateEfficiencyPCR. Ambiguity codes (IUPAC_CODE_MAP) are supported in the primers, but not in myDNAStringSet to prevent trivial matches (e.g., runs of N's).

If taqEfficiency is TRUE (the default), the amplification efficiency of each primer is defined as the product of hybridization efficiency and elongation efficiency. Amplification efficiency must be at least minEfficiency for a primer to be amplified in silico. Overall amplification efficiency of the PCR product is then calculated as the geometric mean of the two (i.e., forward and reverse) primers' efficiencies. Finally, amplicons are generated if the two primers are within maxProductSize nucleotides downstream of each other.

Potential PCR products are returned, either with or without including the primer sequences in the amplicon. The default (includePrimers=TRUE) is to incorporate the primer sequences as would normally occur during amplification. The alternative is to return the complete template sequence including the target sites, which may not exactly match the primer sequences. Note that amplicons may be duplicated when different input primers can amplify the same region of DNA.

Value

A DNAStringSet object containing potential PCR products, sorted from highest-to-lowest amplification efficiency. The sequences are named by their predicted amplification efficiency followed by the index of each primer in primers and the name (or index if names are missing) of the amplified sequence in myDNAStringSet. (See examples section below.)

Note

The program OligoArrayAux (http://www.unafold.org/Dinamelt/software/oligoarrayaux.php) must be installed in a location accessible by the system. For example, the following code should print the installed OligoArrayAux version when executed from the R console:

system("hybrid-min -V")

Author(s)

Erik Wright [email protected]

References

ES Wright et al. (2013) "Exploiting Extension Bias in PCR to Improve Primer Specificity in Ensembles of Nearly Identical DNA Templates." Environmental Microbiology, doi:10.1111/1462-2920.12259.

See Also

CalculateEfficiencyPCR, DesignPrimers, DesignSignatures, MeltDNA

Examples

data(yeastSEQCHR1)

# not run (must have OligoArrayAux installed first):

# match a single primer that acts as both the forward and reverse
primer1 <- "TGGAAGCTGAAACG"
## Not run: AmplifyDNA(primer1, yeastSEQCHR1, annealingTemp=55, P=4e-7, maxProductSize=500)

# perform a typical amplification with two primer sequences:
primer2 <- c("GGCTGTTGTTGGTGTT", "TGTCATCAGAACACCAA")
## Not run: AmplifyDNA(primer2, yeastSEQCHR1, annealingTemp=55, P=4e-7, maxProductSize=500)

# perform a multiplex PCR amplification with multiple primers:
primers <- c(primer1, primer2)
## Not run: AmplifyDNA(primers, yeastSEQCHR1, annealingTemp=55, P=4e-7, maxProductSize=500)

Description

Converts the output of DesignArray into the sparse matrix format used by NNLS.

Usage

Array2Matrix(probes,
             verbose = TRUE)

Arguments

Details

A microarray can be represented by a matrix of hybridization efficiencies, where the rows represent each of the probes and the columns represent each the possible templates. This matrix is sparse since microarray probes are designed to only target a small subset of the possible templates.

Value

A list specifying the hybridization efficiency of each probe to its potential templates.

Author(s)

Erik Wright [email protected]

References

ES Wright et al. (2013) Identification of Bacterial and Archaeal Communities From Source to Tap. Water Research Foundation, Denver, CO.

DR Noguera, et al. (2014). Mathematical tools to optimize the design of oligonucleotide probes and primers. Applied Microbiology and Biotechnology. doi:10.1007/s00253-014-6165-x.

See Also

DesignArray, NNLS

Examples

fas <- system.file("extdata", "Bacteria_175seqs.fas", package="DECIPHER")
dna <- readDNAStringSet(fas)
names(dna) <- 1:length(dna)
probes <- DesignArray(dna)
A <- Array2Matrix(probes)

Description

The BLOSUM amino acid substitution matrices defined by Henikoff, S., & Henikoff, J. (1992).

Usage

data("BLOSUM")

Format

The format is: num [1:24, 1:24, 1:15] 2.4 -0.6 0 0 -1.8 0.6 0 0 -1.2 0 ... - attr(*, "dimnames")=List of 3 ..$ : chr [1:24] "A" "R" "N" "D" ... ..$ : chr [1:24] "A" "R" "N" "D" ... ..$ : chr [1:15] "30" "35" "40" "45" ...

Details

Substitution matrix values represent the log-odds of observing an aligned pair of amino acids versus the likelihood of finding the pair by chance. The PFASUM substitution matrices are stored as an array named by each sub-matrix's similarity threshold. (See examples section below.) In all cases values are in units of third-bits ($log(odds\ ratio)*3/log(2)$).

Source

Henikoff, S., & Henikoff, J. (1992). Amino Acid Substitution Matrices from Protein Blocks. PNAS, 89(22), 10915-10919.

Examples

data(BLOSUM)
BLOSUM62 <- BLOSUM[,, "62"] # the BLOSUM62 matrix
BLOSUM62["A", "R"] # score for A/R pairing

data(PFASUM)
plot(PFASUM[1:20, 1:20, "62"], BLOSUM62[1:20, 1:20])
abline(a=0, b=1)

Description

Opens an html file in a web browser to show the contents of a table in a database.

Usage

BrowseDB(dbFile,
         htmlFile = tempfile(fileext=".html"),
         openURL = interactive(),
         tblName = "Seqs",
         identifier = "",
         limit = -1,
         orderBy = "row_names",
         maxChars = 50,
         title = "",
         clause="")

Arguments

Value

Creates an html table containing all the fields of the database table and (if openURL is TRUE) opens it in the web browser for viewing.

Returns htmlFile if the html file was written successfully.

Note

If viewing a table containing sequences, the sequences are purposefully not shown in the output.

Author(s)

Erik Wright [email protected]

References

ES Wright (2016) "Using DECIPHER v2.0 to Analyze Big Biological Sequence Data in R". The R Journal, 8(1), 352-359.

See Also

BrowseSeqs

Examples

if (require("RSQLite", quietly=TRUE)) {
	db <- system.file("extdata", "Bacteria_175seqs.sqlite", package="DECIPHER")
	BrowseDB(db)
}

Description

Opens an html file in a web browser to show the sequences in an XStringSet.

Usage

BrowseSeqs(myXStringSet,
           htmlFile = tempfile(fileext=".html"),
           openURL = interactive(),
           colorPatterns = TRUE,
           highlight = NA,
           patterns = c("-", alphabet(myXStringSet, baseOnly=TRUE)),
           colors = substring(rainbow(length(patterns),
                              v=0.8, start=0.9, end=0.7), 1, 7),
           colWidth = Inf,
           title = "",
           ...)

Arguments

Details

BrowseSeqs converts an XStringSet into html format for viewing in a web browser. The sequences are colored in accordance with the patterns that are provided, or left uncolored if colorPatterns is FALSE or patterns is NULL. Character or XStringSet patterns are matched as regular expressions and colored according to colors. If patterns is a list of matrices, then it must contain one element per sequence. Each matrix is interpreted as providing the fraction red, blue, and green for each letter in the sequence. Thus, colors is ignored when patterns is a list. (See examples section below.)

Patterns are not matched across column breaks, so multi-character patterns should be carefully considered when colWidth is less than the maximum sequence length. Patterns are matched sequentially in the order provided, so it is feasible to use nested patterns such as c("ACCTG", "CC"). In this case the “CC” could be colored differently inside the previously colored “ACCTG”. Note that patterns overlapping the boundaries of a previously matched pattern will not be matched. For example, “ACCTG” would not be matched if patterns=c("CC", "ACCTG").

Some web browsers cannot quickly display a large amount colored text, so it is recommended to use colorPatterns = FALSE or to highlight a sequence when viewing a large XStringSet. Highlighting will only show all of the characters in the highlighted sequence, and convert all matching positions in the other sequences into dots without color. Also, note that some web browsers display small shifts between fixed-width characters that may become noticeable as color offsets between the ends of long sequences.

Value

Creates an html file containing sequence data and (if openURL is TRUE) opens it in a web browser for viewing. The layout has the sequence name on the left, position legend on the top, cumulative number of nucleotides on the right, and consensus sequence on the bottom.

Returns htmlFile if the html file was written successfully.

Note

Some web browsers do not display colored characters with equal widths. If positions do not align across sequences then try opening the htmlFile with a different web browser.

Author(s)

Erik Wright [email protected]

References

ES Wright (2016) "Using DECIPHER v2.0 to Analyze Big Biological Sequence Data in R". The R Journal, 8(1), 352-359. Kunzmann P., et al. (2020) "Substitution matrix based color schemes for sequence alignment visualization". BMC Bioinformatics, 21(1):209.

See Also

BrowseDB, ConsensusSequence

Examples

# load the example DNA sequences
fas <- system.file("extdata", "Bacteria_175seqs.fas", package="DECIPHER")
dna <- readDNAStringSet(fas) # non-coding ribosomal RNA gene sequences

# example of using the defaults with DNA sequences
BrowseSeqs(dna) # view the XStringSet

# color only "ACTG" and "CSC" patterns (where S is C or G)
BrowseSeqs(dna, patterns=DNAStringSet(c("ACTG", "CSC")))

# highlight (i.e., only fully-color) the first sequence
BrowseSeqs(dna, highlight=1) # other sequences are dots where matching

# highlight the consensus sequence at the bottom
BrowseSeqs(dna, highlight=0) # other sequences are dots where matching

# split the wide view into multiple vertical pages (for printing)
BrowseSeqs(dna, colWidth=100, highlight=1)

# specify an alternative color scheme for -, A, C, G, T
BrowseSeqs(dna, colors=c("#1E90FF", "#32CD32", "#9400D3", "black", "#EE3300"))

# only color the positions within certain positional ranges (100-200 & 250-500)
BrowseSeqs(dna, colorPatterns=c(100, 200, 250, 500))

# example of calling attention to letters by coloring gaps black
BrowseSeqs(dna, patterns="-", colors="black")

# add a title to the top of the page
BrowseSeqs(dna, title="Title")

## Not run: 
# color according to base-pairing by supplying the fraction RGB for every position
dbn <- PredictDBN(dna, type="structures") # calculate the secondary structures
# dbn now contains the scores for whether a base is paired (left/right) or unpaired
dbn[[1]][, 1] # the scores for the first position in the first sequence
dbn[[2]][, 10] # the scores for the tenth position in the second sequence
# these positional scores can be used as shades of red, green, and blue:
BrowseSeqs(dna, patterns=dbn) # red = unpaired, green = left-pairing, blue = right
# positions in black are not part of the consensus secondary structure

## End(Not run)

# color all restriction sites
data(RESTRICTION_ENZYMES) # load dataset containing restriction enzyme sequences
sites <- RESTRICTION_ENZYMES
sites <- gsub("[^A-Z]", "", sites) # remove non-letters
sites <- DNAStringSet(sites) # convert the character vector to a DNAStringSet
rc_sites <- reverseComplement(DNAStringSet(sites))
w <- which(sites != rc_sites) # find non-palindromic restriction sites
sites <- c(sites, rc_sites[w]) # append their reverse complement
sites <- sites[order(nchar(sites))] # match shorter sites first
BrowseSeqs(dna, patterns=sites)

# color bases by quality score
fastq <- system.file("extdata", "s_1_sequence.txt", package="Biostrings")
reads <- readQualityScaledDNAStringSet(fastq, quality.scoring="solexa")
colors <- colorRampPalette(c("red", "yellow", "green"))(42)
colors <- col2rgb(colors)/255
quals <- as(quality(reads), "IntegerList")
quals <- lapply(quals, function(x) colors[, x])
BrowseSeqs(DNAStringSet(reads), patterns=quals) # green = high quality, red = low quality

# load the example protein coding sequences
fas <- system.file("extdata", "50S_ribosomal_protein_L2.fas", package="DECIPHER")
dna <- readDNAStringSet(fas)

# example of using the defaults with amino acid sequences
aa <- unique(translate(dna)) # the unique amino acid sequences
BrowseSeqs(aa)

# example of highlighting the consensus amino acid sequence
AA <- AlignSeqs(aa)
BrowseSeqs(AA, highlight=0)

# example of highlighting positions that differ from the majority consensus
BrowseSeqs(AA, highlight=0, threshold=0.5)

# specify an alternative color scheme for amino acids (from Kunzmann et al.)
colors <- c(`-`="#000000", `A`="#BDB1E8", `R`="#EFA2C5", `N`="#F6602F",
    `D`="#FD5559", `C`="#12C7FE", `Q`="#DDACB4", `E`="#FEA097", `G`="#F46802",
    `H`="#FCA708", `I`="#369BD9", `L`="#2E95EC", `K`="#CF7690", `M`="#4B8EFE",
    `F`="#76997D", `P`="#FD2AE3", `S`="#A08A9A", `T`="#9A84D5", `W`="#74C80D",
    `Y`="#9BB896", `V`="#89B9F9")
BrowseSeqs(AA, colors=colors, patterns=names(colors))

# example of coloring in a reduced amino acid alphabet
alpha <- AA_REDUCED[[15]]
alpha # clustering of amino acids based on similarity
BrowseSeqs(AA, patterns=c("-", paste("[", alpha, "]", sep="")))

# color amino acids according to their predicted secondary structure
hec <- PredictHEC(AA, type="probabilities") # calculate the secondary structures
# hec now contains the probability that a base is in an alpha-helix or beta-sheet
hec[[3]][, 18] # the 18th position in sequence 3 is likely part of a beta-sheet (E)
# the positional probabilities can be used as shades of red, green, and blue:
BrowseSeqs(AA, patterns=hec) # red = alpha-helix, green = beta-sheet, blue = coil

# color codons according to their corresponding amino acid
DNA <- AlignTranslation(dna) # align the translation then reverse translate
colors <- rainbow(21, v=0.8, start=0.9, end=0.7) # codon colors
m <- match(GENETIC_CODE, unique(GENETIC_CODE)) # corresponding amino acid
codonBounds <- matrix(c(seq(1, width(DNA)[1], 3), # start of codons
	seq(3, width(DNA)[1], 3)), # end of codons
	nrow=2,
	byrow=TRUE)
BrowseSeqs(DNA,
	colorPatterns=codonBounds,
	patterns=c("---", names(GENETIC_CODE)), # codons to color
	colors=c("black", substring(colors[m], 1, 7)))

Description

Calculates the Gibbs free energy and hybridization efficiency of probe/target pairs at varying concentrations of the denaturant formamide.

Usage

CalculateEfficiencyArray(probe,
                         target,
                         FA = 0,
                         dGini = 1.96,
                         Po = 10^-2.0021,
                         m = 0.1731,
                         temp = 42,
                         deltaGrules = NULL)

Arguments

Details

This function calculates the free energy and hybridization efficiency (HE) for a given formamide concentration ([FA]) using the linear free energy model given by:

$HE = Po*exp[-(dG_0 + m*FA)/RT]/(1+Po*exp[-(dG_0 + m*FA)/RT])$

The probe and target input sequences must be aligned in pairs, such that the first probe is aligned to the first target, second-to-second, and so on. Ambiguity codes in the IUPAC_CODE_MAP are accepted in probe and target sequences. Any ambiguities will default to perfect match pairings by inheriting the nucleotide in the same position on the opposite sequence whenever possible. If the ambiguity results in a mismatch then “T”, “G”, “C”, and “A” are substituted, in that order. For example, if a probe nucleotide is “S” (“C” or “G”) then it will be considered a “C” if the target nucleotide in the same position is a “C”, otherwise the ambiguity will be interpreted as a “G”.

If deltaGrules is NULL then the rules defined in data(deltaGrules) will be used. Note that deltaGrules of the same format may be customized for any application and specified as an input.

Value

A matrix with the predicted Gibbs free energy (dG) and hybridization efficiency (HE) at each concentration of formamide ([FA]).

Author(s)

Erik Wright [email protected]

References

Yilmaz LS, Loy A, Wright ES, Wagner M, Noguera DR (2012) Modeling Formamide Denaturation of Probe-Target Hybrids for Improved Microarray Probe Design in Microbial Diagnostics. PLoS ONE 7(8): e43862. doi:10.1371/journal.pone.0043862.

See Also

deltaGrules

Examples

probes <- c("AAAAACGGGGAGCGGGGGGATACTG", "AAAAACTCAACCCGAGGAGCGGGGG")
targets <- c("CAACCCGGGGAGCGGGGGGATACTG", "TCGGGCTCAACCCGAGGAGCGGGGG")
result <- CalculateEfficiencyArray(probes, targets, FA=0:40)
dG0 <- result[, "dG_0"]
HE0 <- result[, "HybEff_0"]
plot(result[1, 1:40], xlab="[FA]", ylab="HE", main="Probe/Target # 1", type="l")

Description

Calculates the Gibbs free energy, formamide melt point, and hybridization efficiency of probe/target (DNA/RNA) pairs.

Usage

CalculateEfficiencyFISH(probe,
                        target,
                        temp,
                        P,
                        ions,
                        FA,
                        batchSize = 1000)

Arguments

Details

Hybridization of pairwise probe/target (DNA/RNA) pairs is simulated in silico. Gibbs free energies are obtained from system calls to OligoArrayAux, which must be properly installed (see the Notes section below). Probe/target pairs are sent to OligoArrayAux in batches of batchSize, which prevents systems calls from being too many characters. Note that OligoArrayAux does not support degeneracy codes (non-base letters), although they are accepted without error. Any sequences with ambiguity should be expanded into multiple permutations with Disambiguate before input.

Value

A matrix of predicted hybridization efficiency (HybEff), formamide melt point (FAm), and free energy (ddG1 and dG1) for each probe/target pair of sequences.

Note

The program OligoArrayAux (http://www.unafold.org/Dinamelt/software/oligoarrayaux.php) must be installed in a location accessible by the system. For example, the following code should print the installed OligoArrayAux version when executed from the R console:

system("hybrid-min -V")

Author(s)

Erik Wright [email protected]

References

ES Wright et al. (2014) "Automated Design of Probes for rRNA-Targeted Fluorescence In Situ Hybridization Reveals the Advantages of Using Dual Probes for Accurate Identification." Applied and Environmental Microbiology, doi:10.1128/AEM.01685-14.

See Also

DesignProbes, TileSeqs

Examples

probe <- c("GGGCTTTCACATCAGACTTAAGAAACC", "CCCCACGCTTTCGCGCC")
target <- reverseComplement(DNAStringSet(probe))
# not run (must have OligoArrayAux installed first):
## Not run: CalculateEfficiencyFISH(probe, target, temp=46, P=250e-9, ions=1, FA=35)

Description

Calculates the amplification efficiency of primers from their hybridization efficiency and elongation efficiency at the target site.

Usage

CalculateEfficiencyPCR(primer,
                       target,
                       temp,
                       P,
                       ions,
                       batchSize = 1000,
                       taqEfficiency = TRUE,
                       maxDistance = 0.4,
                       maxGaps = 2,
                       processors = 1)

Arguments

Details

Amplification of pairwise primer/target pairs is simulated in silico. A complex model of hybridization is employed that takes into account the side reactions resulting from probe-folding, target-folding, and primer-dimer formation. The resulting hybridization efficiency is multiplied by the elongation efficiency to predict the overall efficiency of amplification.

Free energy is obtained from system calls to OligoArrayAux, which must be properly installed (see the Notes section below). Primer/target pairs are sent to OligoArrayAux in batches of batchSize, which prevents systems calls from being too many characters. Note that OligoArrayAux does not support degeneracy codes (non-base letters), although they are accepted without error. Any sequences with ambiguity should be expanded into multiple permutations with Disambiguate before input.

Value

A vector of predicted efficiencies for amplifying each primer/target pair of sequences.

Note

The program OligoArrayAux (http://www.unafold.org/Dinamelt/software/oligoarrayaux.php) must be installed in a location accessible by the system. For example, the following code should print the installed OligoArrayAux version when executed from the R console:

system("hybrid-min -V")

Author(s)

Erik Wright [email protected]

References

ES Wright et al. (2013) "Exploiting Extension Bias in PCR to Improve Primer Specificity in Ensembles of Nearly Identical DNA Templates." Environmental Microbiology, doi:10.1111/1462-2920.12259.

See Also

AmplifyDNA, DesignPrimers, DesignSignatures

Examples

primers <- c("AAAAACGGGGAGCGGGGGG", "AAAAACTCAACCCGAGGAGCGCGT")
targets <- reverseComplement(DNAStringSet(primers))
# not run (must have OligoArrayAux installed first):
## Not run: CalculateEfficiencyPCR(primers, targets, temp=75, P=4e-7, ions=0.225)

Description

Groups the sequences into approximate clusters of similarity.

Usage

Clusterize(myXStringSet,
           cutoff = 0,
           method = "overlap",
           includeTerminalGaps = FALSE,
           penalizeGapLetterMatches = NA,
           minCoverage = 0.5,
           maxPhase1 = 2e4,
           maxPhase2 = 2e3,
           maxPhase3 = 2e3,
           maxAlignments = 200,
           rareKmers = 50,
           probability = 0.99,
           invertCenters = FALSE,
           singleLinkage = FALSE,
           maskRepeats = FALSE,
           maskLCRs = FALSE,
           alphabet = AA_REDUCED[[186]],
           processors = 1,
           verbose = TRUE)

Arguments

Details

Clusterize groups the input sequences into approximate clusters using a heuristic algorithm with linear time and memory complexity. In phase 1, the sequences are partitioned into groups of similarity. In phase 2, the sequences are ordered by k-mer similarity by relatedness sorting. In phase 3, the sequences are iteratively clustered in this order by their similarity to surrounding sequences in the sorting. That is, the first sequence becomes the representative of cluster #1. If the second sequence is within cutoff distance then it is added to the cluster, otherwise it becomes a new cluster representative. The remaining sequences are matched to cluster representatives in a similar fashion until all sequences belong to a cluster. In the majority of cases, this process results in clusters with members separated by less than cutoff distance, and all cluster members must be within cutoff distance of their cluster representative.

The calculation of distance can be controlled in multiple ways, with each parameterization of distance having advantages and disadvantages. By default, distance is the fraction of positions that are different, including gaps, within the overlapping region in a pairwise alignment. The defaults will handle partial-length sequences well, but also cluster sequences with high similarity between their opposite ends. For this reason, it is important to set minimumCoverage such that distances are based off of considerable overlap between sequences in the pairwise alignment. The default (0.5) requires sequences to share at least 50% of their positions with the cluster representative. This distance parameterization works well, but there are reasonable alternatives.

If penalizeGapLetterMatches is FALSE, the distance will exclude gap regions. If includeTerminalGaps is TRUE, the calculation of distance will use the entire (global) alignment. If method is "shortest" and includeTerminalGaps is set to TRUE, then the distance is calculated for the region encompassed by the shorter sequence in each pair, which is the common definition of distance used by other clustering programs. This common definition of distance will sometimes separate partly overlapping sequences, which is why it is not the default. Another option is to set minCoverage to a negative fraction. This requires sequences to have substantial overlap with the cluster representative, which is the longest sequence in the cluster. For example, setting minCoverage to -0.5 would require every clustered sequence to share at least 50% of positions with the cluster representative.

The algorithm requires time proportional to the number of input sequences in myXStringSet. The phase 1, up to maxPhase1 sequences sharing a k-mer are tabulated while partitioning each sequence. In phase 2, the sequences are compared with up to maxPhase2 passes that each take linear time. Ordering of the sequences is performed in linear time using radix sorting. In phase 3, each sequence is compared with up to maxPhase3 previous cluster representatives of sequences sharing rareKmers or nearby sequences in the relatedness ordering. This is possible because the sequences are sorted by relatedness, such that more recent cluster representatives are more similar. Hence, the complete algorithm scales in linear time asymptotically and returns clusters of sequences within cutoff distance of their center sequence.

Multiple cutoffs can be provided in sorted order, which saves time because phases 1 and 2 only need to be performed once. If the cutoffs are provided in descending order then clustering at each new value of cutoff is continued within the prior cutoff's clusters. In this way clusters at lower values of cutoff are completely contained within their “umbrella” clusters at higher values of cutoff. This slightly accelerates the clustering process, because each subsequent group is only clustered within the previous group. If multiple cutoffs are provided in ascending order then clustering at each level of cutoff is independent of the prior level.

Note, the clustering algorithm is stochastic. Hence, clusters can vary from run-to-run unless the random number seed is set for repeatability (i.e., with set.seed). Also, invertCenters can be used to determine the center sequence of each cluster from the output. Since identical sequences will always be assigned the same cluster numbers, it is possible for more than one input sequence in myXStringSet to be assigned as the center of a cluster if they are identical.

Value

A data.frame is returned with dimensions $N*M$, where each one of $N$ sequences is assigned to a cluster at the $M$-level of cutoff. The row.names of the data.frame correspond to the names of myXStringSet.

Author(s)

Erik Wright [email protected]

References

Wright, E. S. (2024). Accurately clustering biological sequences in linear time by relatedness sorting. Nature Communications, 15, 1-13. http://doi.org/10.1038/s41467-024-47371-9

See Also

AA_REDUCED, DistanceMatrix, Treeline

Examples

fas <- system.file("extdata", "50S_ribosomal_protein_L2.fas", package="DECIPHER")
dna <- readDNAStringSet(fas)
aa <- translate(dna)

# typical usage (e.g., clustering at >= 90 percent similarity)
clusters <- Clusterize(aa, cutoff=0.1) # set processors = NULL for max speed
head(clusters)

# typical usage (e.g., obtaining cluster representatives)
clusters <- Clusterize(aa, cutoff=0.1, invertCenters=TRUE)
aa[clusters[[1]] < 0]

# cluster each cutoff within the previous cluster (nested groups)
clusters <- Clusterize(aa, cutoff=seq(0.7, 0, -0.1))
head(clusters)
apply(clusters, 2, max) # number of clusters per cutoff

# cluster each cutoff independently (possibly fewer clusters per cutoff)
clusters <- Clusterize(aa, cutoff=seq(0, 0.7, 0.1))
head(clusters)
apply(clusters, 2, max) # number of clusters per cutoff

# make cluster center(s) negative for tracking
clusters <- Clusterize(aa, cutoff=0.5, invertCenters=TRUE)
head(clusters)
clusters[clusters$cluster < 0,, drop=FALSE]
unique(aa[clusters$cluster < 0]) # unique cluster centers
apply(clusters, 2, function(x) max(abs(x))) # number of clusters

# cluster nucleotide sequences
clusters <- Clusterize(dna, cutoff=0.2, invertCenters=TRUE)
head(clusters)
apply(clusters, 2, function(x) max(abs(x))) # number of clusters

Description

Compresses character vectors into raw vectors, or decompresses raw vectors into character vectors using a variety of codecs.

Usage

Codec(x,
      compression,
      compressRepeats = FALSE,
      processors = 1)

Arguments

Details

Codec can be used to compress/decompress character vectors using different algorithms. The "nbit" and "qbit" methods are tailored specifically to nucleotides and quality scores, respectively. These two methods will store the data as plain text ("ASCII" format) when it is incompressible. In such cases, a second compression method can be given to use in lieu of plain text. For example compression = c("nbit", "gzip") will use "gzip" compression when "nbit" compression is inappropriate.

When performing the reverse operation, decompression, the type of compression is automatically detected based on the unique signature ("magic number") added by each compression algorithm.

Value

If x is a character vector to be compressed, the output is a list with one element containing a raw vector per character string. If x is a list of raw vectors to be decompressed, then the output is a character vector with one string per list element.

Author(s)

Erik Wright [email protected]

Examples

fas <- system.file("extdata", "Bacteria_175seqs.fas", package="DECIPHER")
dna <- as.character(readDNAStringSet(fas)) # aligned sequences
object.size(dna)

# compression
system.time(x <- Codec(dna, compression="nbit"))
object.size(x)/sum(nchar(dna)) # bytes per position

system.time(g <- Codec(dna, compression="gzip"))
object.size(g)/sum(nchar(dna)) # bytes per position

# decompression
system.time(y <- Codec(x))
stopifnot(dna==y)

system.time(z <- Codec(g))
stopifnot(dna==z)

Description

Forms a consensus sequence representing a set of sequences.

Usage

ConsensusSequence(myXStringSet,
                  threshold = 0.05,
                  ambiguity = TRUE,
                  noConsensusChar = "+",
                  minInformation = 1 - threshold,
                  includeNonLetters = FALSE,
                  includeTerminalGaps = FALSE)

Arguments

Details

ConsensusSequence removes the least frequent characters at each position, so long as they represent less than threshold fraction of the sequences in total. If necessary, ConsensusSequence represents the remaining characters using a degeneracy code from the IUPAC_CODE_MAP. Degeneracy codes are always used in cases where multiple characters are equally abundant.

Two key parameters control the degree of consensus: threshold and minInformation. The default threshold (0.05) means that at less than 5% of sequences will not be represented by the consensus sequence at any given position. The default minInformation (1 - 0.05) specifies that at least 95% of sequences must contain the information in the consensus, otherwise the noConsensusChar is used. This enables an alternative character (e.g., "+") to be substituted at positions that would otherwise yield an ambiguity code.

If ambiguity = TRUE (the default) then degeneracy codes in myXStringSet are split between their respective bases according to the IUPAC_CODE_MAP for DNA/RNA and AMINO_ACID_CODE for AA. For example, an “R” in a DNAStringSet would count as half an “A” and half a “G”. If ambiguity = FALSE then degeneracy codes are not considered in forming the consensus. For an AAStringSet input, the lack of degeneracy codes generally results in “X” at positions with mismatches, unless the threshold is set to a higher value than the default.

If includeNonLetters = TRUE (the default) then gap ("-"), mask ("+"), and unknown (".") characters are counted towards the consensus, otherwise they are omitted from calculation of the consensus. Note that gap ("-") and unknown (".") characters are treated interchangeably as gaps when forming the consensus sequence. For this reason, the consensus of a position with all unknown (".") characters will be a gap ("-"). Also, note that if consensus is formed between different length sequences then it will represent only the longest sequences at the end. For this reason the consensus sequence is generally based on a sequence alignment such that all of the sequences have equal lengths.

Value

An XStringSet with a single consensus sequence matching the input type.

Author(s)

Erik Wright [email protected]

See Also

Disambiguate, IdConsensus

Examples

fas <- system.file("extdata", "Bacteria_175seqs.fas", package="DECIPHER")
dna <- readDNAStringSet(fas, n=10) # limit 10 sequences for visability
BrowseSeqs(dna) # consensus at bottom
BrowseSeqs(dna, threshold=0.5) # consensus at bottom

# controlling the degree of consensus
AAAT <- DNAStringSet(c("A", "A", "A", "T"))
ConsensusSequence(AAAT) # "W"
ConsensusSequence(AAAT, threshold=0.3) # "A"
ConsensusSequence(AAAT, threshold=0.3, minInformation=0.8) # "+"
ConsensusSequence(AAAT, threshold=0.3, minInformation=0.8, noConsensusChar="N") # "N"

# switch between degenerate-based and majority-based consensus
majority <- DNAStringSet(c("GTT", "GAA", "CTG"))
ConsensusSequence(majority) # degenerate-based
ConsensusSequence(majority, threshold=0.5) # majority-based
ConsensusSequence(majority, threshold=0.5, minInformation=0.75)

# behavior in the case of a tie
ConsensusSequence(DNAStringSet(c("A", "T"))) # "W"
ConsensusSequence(DNAStringSet(c("A", "T")), threshold=0.5) # "W"
ConsensusSequence(AAStringSet(c("A", "T"))) # "X"
ConsensusSequence(AAStringSet(c("A", "T")), threshold=0.5) # "X"
ConsensusSequence(AAStringSet(c("I", "L"))) # "J"
ConsensusSequence(AAStringSet(c("I", "L")), threshold=0.5) # "J"

# handling terminal gaps
dna <- DNAStringSet(c("ANGCT-","-ACCT-"))
ConsensusSequence(dna) # "ANSCT-"
ConsensusSequence(dna, includeTerminalGaps=TRUE) # "+NSCT-"

# the "." character is treated is a "-"
aa <- AAStringSet(c("ANQIH-", "ADELW."))
ConsensusSequence(aa) # "ABZJX-"

# internal non-bases are included by default
ConsensusSequence(DNAStringSet(c("A-+.A", "AAAAA")), noConsensusChar="N") # "ANNNA"
ConsensusSequence(DNAStringSet(c("A-+.A", "AAAAA")), includeNonLetters=TRUE) # "AAAAA"

# degeneracy codes in the input are considered by default
ConsensusSequence(DNAStringSet(c("AWNDA", "AAAAA"))) # "AWNDA"
ConsensusSequence(DNAStringSet(c("AWNDA", "AAAAA")), ambiguity=FALSE) # "AAAAA"

Description

Calculates the matrix of cophenetic distances represented by a dendrogram object.

Usage

Cophenetic(x)

Arguments

Details

The cophenetic distance between two observations is defined here as the branch length separating them on a dendrogram, also known as the patristic distance. This function differs from the cophenetic function in that it does not assume the tree is ultrametric and outputs the branch length separating pairs of observations rather than the height of their merger. A dendrogram that better preserves a distance matrix will show higher correlation between the distance matrix and its cophenetic distances.

Value

An object of class 'dist'.

Author(s)

Erik Wright [email protected]

See Also

Treeline

Examples

fas <- system.file("extdata", "Bacteria_175seqs.fas", package="DECIPHER")
dna <- readDNAStringSet(fas)
d1 <- DistanceMatrix(dna, type="dist")
dend <- Treeline(myDistMatrix=d1)
d2 <- Cophenetic(dend)
cor(d1, d2)

Description

Corrects the reading frame to mitigate the impact of frameshift errors caused by insertions or deletions in unaligned nucleotide sequences.

Usage

CorrectFrameshifts(myXStringSet,
                   myAAStringSet,
                   type = "indels",
                   acceptDistance = 0.01,
                   rejectDistance = 0.60,
                   maxComparisons = 10,
                   gapOpening = -13,
                   gapExtension = -1,
                   frameShift = -15,
                   geneticCode = GENETIC_CODE,
                   substitutionMatrix = "PFASUM50",
                   verbose = TRUE,
                   processors = 1)

Arguments

Details

Accurate translation of protein coding sequences can be greatly disrupted by one or two nucleotide phase shifts that occasionally occur during DNA sequencing. These frameshift errors can potentially be corrected through comparison with other unshifted protein sequences. This function uses a set of reference amino acid sequences (AAStringSet) to find and correct frameshift errors in a set of nucleotide sequences (myXStringSet). First, three frame translation of the nucleotide sequences is performed, and the nearest reference sequence is selected. Then the optimal reading frame at each position is determined based on a variation of the Guan & Uberbacher (1996) method. Putative insertions and/or deletions (indels) are returned in the result, typically with close proximity to the true indel locations. For a comparison of this method to others, see Wang et al. (2013).

If type is "sequences" or "both", then frameshifts are corrected by adding N's and/or removing nucleotides. Note that this changes the nucleotide sequence, and the new sequence often has minor errors because the exact location of the indel(s) cannot be determined. However, the original frameshifts that disrupted the entire downstream sequence are reduced to local perturbations. All of the returned nucleotide sequences will have a reading frame starting from the first position. This allows direct translation, and in practice works well if there is a similar reference myAAStringSet with the correct reading frame. Hence it is more important that myAAStringSet contain a wide variety of sequences than it is that it contain a lot of sequences.

Multiple inputs control the time required for frameshift correction. The number of sequences in the reference set (myAAStringSet) will affect the speed of the first search for similar sequences. Assessing frameshifts in the second step requires order N*M time, where N and M are the lengths of the query (myXStringSet) and reference sequences. Two parameters control the number of assessments that are made for each sequence: (1) maxComparisons determines the maximum number of reference sequences to compare to each query sequence, and (2) acceptDist defines the maximum distance between a query and reference that is acceptable before continuing to the next query sequence. A lower value for maxComparisons or a higher value for acceptDist will accelerate frameshift correction, potentially at the expense of some accuracy.

Value

If type is "indels" then the returned object is a list with the same length as myXStringSet. Each element is a list with four components:

Note that positions in insertions and deletions are sometimes repeated to indicate that the same position needs to be shifted successively more than once to correct the reading frame.

If type is "sequences" then the returned object is an XStringSet of the same type as the input (myXStringSet). Nucleotides are added or deleted as necessary to correct for frameshifts. The returned sequences all have a reading frame starting from position 1, so that they can be translated directly.

If type is "both" then the returned object is a list with two components: one for the "indels" and the other for the "sequences".

Author(s)

Erik Wright [email protected]

References

Guan, X., & Uberbacher, E. C. (1996). Alignments of DNA and protein sequences containing frameshift errors. Computer Applications in the Biosciences : CABIOS, 12(1), 31-40.

Wang, Q., et al. (2013). Ecological Patterns of nifH Genes in Four Terrestrial Climatic Zones Explored with Targeted Metagenomics Using FrameBot, a New Informatics Tool. mBio, 4(5), e00592-13-e00592-13.

See Also

AlignTranslation, OrientNucleotides, PFASUM

Examples

fas <- system.file("extdata", "50S_ribosomal_protein_L2.fas", package="DECIPHER")
dna <- readDNAStringSet(fas)

# introduce artificial indels
n_ins <- 2 # insertions per sequence
shifted <- replaceAt(dna,
	lapply(width(dna),
		sample,
		n_ins),
	sample(DNA_BASES,
		n_ins,
		replace=TRUE))
n_dels <- 1 # deletions per sequence
shifted <- replaceAt(shifted,
	as(lapply(width(shifted),
		function(x) {
			IRanges(sample(x,
					n_dels),
				width=1)
		}), "IRangesList"))

# to make frameshift correction more challenging,
# only supply 20 reference amino acid sequences
s <- sample(length(dna), 20)
x <- CorrectFrameshifts(shifted,
	translate(dna[s]),
	type="both")

# there was a wide range of distances
# to the nearest reference sequence
quantile(unlist(lapply(x[[1]], `[`, "distance")))

# none of the sequences were > rejectDistance
# from the nearest reference sequence
length(which(unlist(lapply(x[[1]], `[`, "index"))==0))

# the number of indels was generally correct
table(unlist(lapply(x[[1]], function(x) {
	length(x$insertions)})))/length(shifted)
table(unlist(lapply(x[[1]], function(x) {
	length(x$deletions)})))/length(shifted)

# align and display the translations
AA <- AlignTranslation(x$sequences,
	readingFrame=1,
	type="AAStringSet")
BrowseSeqs(AA)

Description

Creates artificial random chimeras from a set of sequences.

Usage

CreateChimeras(myDNAStringSet,
               numChimeras = 10,
               numParts = 2,
               minLength = 80,
               maxLength = Inf,
               minChimericRegionLength = 30,
               randomLengths = TRUE,
               includeParents = TRUE,
               processors = 1,
               verbose = TRUE)

Arguments

Details

Forms a set of random chimeras from the input set of (typically good quality) sequences. The chimeras are created by merging random sequences at random breakpoints. These chimeras can be used for testing the accuracy of the FindChimeras or other chimera finding functions.

Value

A DNAStringSet object containing chimeras. The names of the chimeras are specified as "parent #1 name [chimeric region] (distance from parent to chimera), ...".

If includeParents = TRUE then the parents of the chimeras are included at the end of the result. The parents are trimmed to the same length as the chimera if randomLengths = TRUE. The names of the parents are specified as "parent #1 name [region] (distance to parent #2, ...)".

Author(s)

Erik Wright [email protected]

See Also

FindChimeras, Seqs2DB

Examples

fas <- system.file("extdata", "Bacteria_175seqs.fas", package="DECIPHER")
dna <- readDNAStringSet(fas)
chims <- CreateChimeras(dna)
BrowseSeqs(chims)

Description

Exports a database containing sequences to a FASTA or FASTQ formatted file of sequence records.

Usage

DB2Seqs(file,
         dbFile,
         tblName = "Seqs",
         identifier = "",
         type = "BStringSet",
         limit = -1,
         replaceChar = NA,
         nameBy = "description",
         orderBy = "row_names",
         removeGaps = "none",
         append = FALSE,
         width = 80,
         compress = FALSE,
         chunkSize = 1e5,
         sep = "::",
         clause = "",
         processors = 1,
         verbose = TRUE)

Arguments

Details

Sequences are exported into either a FASTA or FASTQ file as determined by the type of sequences. If type is an XStringSet then sequences are exported to FASTA format. Quality information for QualityScaledXStringSets are interpreted as PredQuality scores before export to FASTQ format.

If type is "BStringSet" (the default) then sequences are exported to a FASTA file exactly the same as they were when imported. If type is "DNAStringSet" then all U's are converted to T's before export, and vise-versa if type is "RNAStringSet". All remaining characters not in the XStringSet's alphabet are converted to replaceChar or removed if replaceChar is "". Note that if replaceChar is NA (the default), it will result in an error when an unexpected character is found.

Value

Writes a FASTA or FASTQ formatted file containing the sequence records in the database.

Returns the number of sequence records written to the file.

Author(s)

Erik Wright [email protected]

References

ES Wright (2016) "Using DECIPHER v2.0 to Analyze Big Biological Sequence Data in R". The R Journal, 8(1), 352-359.

Examples

if (require("RSQLite", quietly=TRUE)) {
	db <- system.file("extdata", "Bacteria_175seqs.sqlite", package="DECIPHER")
	tf <- tempfile()
	DB2Seqs(tf, db, limit=10)
	file.show(tf) # press 'q' to exit
	unlink(tf)
}

Description

An 8D array with four adjacent base pairs of the probe and target sequences at a time. Each dimension has five elements defining the nucleotide at that position ("A", "C", "G", "T", or "-"). The array contains the standard Gibbs free energy change of probe binding (dG, [kcal/mol]) for every quadruple base pairing.

Usage

data(deltaGrules)

Format

The format is: num [1:5, 1:5, 1:5, 1:5, 1:5, 1:5, 1:5, 1:5] -0.141 0 0 0 0 ... - attr(*, "dimnames")=List of 8 ..$ : chr [1:5] "A" "C" "G" "T" ... ..$ : chr [1:5] "A" "C" "G" "T" ... ..$ : chr [1:5] "A" "C" "G" "T" ... ..$ : chr [1:5] "A" "C" "G" "T" ... ..$ : chr [1:5] "A" "C" "G" "T" ... ..$ : chr [1:5] "A" "C" "G" "T" ... ..$ : chr [1:5] "A" "C" "G" "T" ... ..$ : chr [1:5] "A" "C" "G" "T" ...

Details

The first four dimensions correspond to the four probe positions from 5' to 3'. The fifth to eighth dimensions correspond to the four positions from 5' to 3' of the target sequence.

Source

Data obtained using NimbleGen microarrays and a Linear Free Energy Model developed by Yilmaz et al.

References

Yilmaz LS, Loy A, Wright ES, Wagner M, Noguera DR (2012) Modeling Formamide Denaturation of Probe-Target Hybrids for Improved Microarray Probe Design in Microbial Diagnostics. PLoS ONE 7(8): e43862. doi:10.1371/journal.pone.0043862.

Examples

data(deltaGrules)
# dG of probe = AGCT / target = A-CT pairing
deltaGrules["A", "G", "C", "T", "A", "-", "C", "T"]

Description

An 8D array with four adjacent base pairs of the RNA duplex. Each dimension has five elements defining the nucleotide at that position ("A", "C", "G", "U", or "-"). The array contains the pseudoenergy of duplex formation for every quadruple base pairing.

Usage

data("deltaGrulesRNA")

Format

The format is: num [1:5, 1:5, 1:5, 1:5, 1:5, 1:5, 1:5, 1:5] -0.776 -0.197 -0.197 -0.291 0 ... - attr(*, "dimnames")=List of 8 ..$ : chr [1:5] "A" "C" "G" "U" ... ..$ : chr [1:5] "A" "C" "G" "U" ... ..$ : chr [1:5] "A" "C" "G" "U" ... ..$ : chr [1:5] "A" "C" "G" "U" ... ..$ : chr [1:5] "A" "C" "G" "U" ... ..$ : chr [1:5] "A" "C" "G" "U" ... ..$ : chr [1:5] "A" "C" "G" "U" ... ..$ : chr [1:5] "A" "C" "G" "U" ...

Details

The first four dimensions correspond to the four top strand positions from 5' to 3'. The fifth to eighth dimensions correspond to the four bottom strand positions from 5' to 3'.

Source

Psuedoenergy values of ungapped quadruplets are inferred from their log-odds of being found in palindromes within hairpin regions across thousands of non-coding RNA alignments. Each value represents the log-odds of in vivo folding relative to chance.

Examples

data(deltaGrulesRNA)
# dG of the duplex AGCU / ACCU pairing (1 mismatch)
deltaGrulesRNA["A", "G", "C", "U", "A", "C", "C", "U"]

Description

An 8D array with four adjacent base pairs of the DNA duplex. Each dimension has five elements defining the nucleotide at that position ("A", "C", "G", "T", or "-"). The array contains the standard enthalpy change of probe binding (dH, [kcal/mol]) for every quadruple base pairing.

Usage

data(deltaHrules)

Format

The format is: num [1:5, 1:5, 1:5, 1:5, 1:5, 1:5, 1:5, 1:5] -7.97 0 0 0 0 ... - attr(*, "dimnames")=List of 8 ..$ : chr [1:5] "A" "C" "G" "T" ... ..$ : chr [1:5] "A" "C" "G" "T" ... ..$ : chr [1:5] "A" "C" "G" "T" ... ..$ : chr [1:5] "A" "C" "G" "T" ... ..$ : chr [1:5] "A" "C" "G" "T" ... ..$ : chr [1:5] "A" "C" "G" "T" ... ..$ : chr [1:5] "A" "C" "G" "T" ... ..$ : chr [1:5] "A" "C" "G" "T" ...

Details

The first four dimensions correspond to the four top strand positions from 5' to 3'. The fifth to eighth dimensions correspond to the four bottom strand positions from 5' to 3'.

Source

Data from a variety of publications by SantaLucia et al.

References

SantaLucia, J., Jr., & Hicks, D. (2004) The Thermodynamics of DNA Structural Motifs. Annual Review of Biophysics and Biomolecular Structure, 33(1), 415-440. doi:10.1146/annurev.biophys.32.110601.141800.

Examples

data(deltaHrules)
# dH of the duplex AGCT / A-CT pairing
deltaHrules["A", "G", "C", "T", "A", "-", "C", "T"]

Description

An 8D array with four adjacent base pairs of the RNA duplex. Each dimension has five elements defining the nucleotide at that position ("A", "C", "G", "U", or "-"). The array contains the standard enthalpy change of probe binding (dH, [kcal/mol]) for every quadruple base pairing.

Usage

data(deltaHrulesRNA)

Format

The format is: num [1:5, 1:5, 1:5, 1:5, 1:5, 1:5, 1:5, 1:5] -6.55 0 0 0 0 ... - attr(*, "dimnames")=List of 8 ..$ : chr [1:5] "A" "C" "G" "U" ... ..$ : chr [1:5] "A" "C" "G" "U" ... ..$ : chr [1:5] "A" "C" "G" "U" ... ..$ : chr [1:5] "A" "C" "G" "U" ... ..$ : chr [1:5] "A" "C" "G" "U" ... ..$ : chr [1:5] "A" "C" "G" "U" ... ..$ : chr [1:5] "A" "C" "G" "U" ... ..$ : chr [1:5] "A" "C" "G" "U" ...

Details

The first four dimensions correspond to the four top strand positions from 5' to 3'. The fifth to eighth dimensions correspond to the four bottom strand positions from 5' to 3'.

Source

Data from a variety of publications by SantaLucia et al.

References

SantaLucia, J., Jr., & Hicks, D. (2004) The Thermodynamics of DNA Structural Motifs. Annual Review of Biophysics and Biomolecular Structure, 33(1), 415-440. doi:10.1146/annurev.biophys.32.110601.141800.

Examples

data(deltaHrulesRNA)
# dH of the duplex AGCU / A-CU pairing
deltaHrulesRNA["A", "G", "C", "U", "A", "-", "C", "U"]

Description

An 8D array with four adjacent base pairs of the DNA duplex. Each dimension has five elements defining the nucleotide at that position ("A", "C", "G", "T", or "-"). The array contains the standard entropy change of probe binding (dS, [kcal/mol]) for every quadruple base pairing.

Usage

data(deltaSrules)

Format

The format is: num [1:5, 1:5, 1:5, 1:5, 1:5, 1:5, 1:5, 1:5] -0.0226 0 0 0 0 ... - attr(*, "dimnames")=List of 8 ..$ : chr [1:5] "A" "C" "G" "T" ... ..$ : chr [1:5] "A" "C" "G" "T" ... ..$ : chr [1:5] "A" "C" "G" "T" ... ..$ : chr [1:5] "A" "C" "G" "T" ... ..$ : chr [1:5] "A" "C" "G" "T" ... ..$ : chr [1:5] "A" "C" "G" "T" ... ..$ : chr [1:5] "A" "C" "G" "T" ... ..$ : chr [1:5] "A" "C" "G" "T" ...

Details

The first four dimensions correspond to the four top strand positions from 5' to 3'. The fifth to eighth dimensions correspond to the four bottom strand positions from 5' to 3'.

Source

Data from a variety of publications by SantaLucia et al.

References

SantaLucia, J., Jr., & Hicks, D. (2004) The Thermodynamics of DNA Structural Motifs. Annual Review of Biophysics and Biomolecular Structure, 33(1), 415-440. doi:10.1146/annurev.biophys.32.110601.141800.

Examples

data(deltaSrules)
# dS of the duplex AGCT / A-CT pairing
deltaSrules["A", "G", "C", "T", "A", "-", "C", "T"]

Description

An 8D array with four adjacent base pairs of the RNA duplex. Each dimension has five elements defining the nucleotide at that position ("A", "C", "G", "T", or "-"). The array contains the standard entropy change of probe binding (dS, [kcal/mol]) for every quadruple base pairing.

Usage

data(deltaSrulesRNA)

Format

The format is: num [1:5, 1:5, 1:5, 1:5, 1:5, 1:5, 1:5, 1:5] -0.0182 0 0 0 0 ... - attr(*, "dimnames")=List of 8 ..$ : chr [1:5] "A" "C" "G" "U" ... ..$ : chr [1:5] "A" "C" "G" "U" ... ..$ : chr [1:5] "A" "C" "G" "U" ... ..$ : chr [1:5] "A" "C" "G" "U" ... ..$ : chr [1:5] "A" "C" "G" "U" ... ..$ : chr [1:5] "A" "C" "G" "U" ... ..$ : chr [1:5] "A" "C" "G" "U" ... ..$ : chr [1:5] "A" "C" "G" "U" ...

Details

The first four dimensions correspond to the four top strand positions from 5' to 3'. The fifth to eighth dimensions correspond to the four bottom strand positions from 5' to 3'.

Source

Data from a variety of publications by SantaLucia et al.

References

SantaLucia, J., Jr., & Hicks, D. (2004) The Thermodynamics of DNA Structural Motifs. Annual Review of Biophysics and Biomolecular Structure, 33(1), 415-440. doi:10.1146/annurev.biophys.32.110601.141800.

Examples

data(deltaSrulesRNA)
# dS of the duplex AGCU / A-CU pairing
deltaSrulesRNA["A", "G", "C", "U", "A", "-", "C", "U"]

Description

Chooses the set of microarray probes maximizing sensitivity and specificity to each target consensus sequence.

Usage

DesignArray(myDNAStringSet,
            maxProbeLength = 24,
            minProbeLength = 20,
            maxPermutations = 4,
            numRecordedMismatches = 500,
            numProbes = 10,
            start = 1,
            end = NULL,
            maxOverlap = 5,
            hybridizationFormamide = 10,
            minMeltingFormamide = 15,
            maxMeltingFormamide = 20,
            minScore = -1e+12,
            processors = 1,
            verbose = TRUE)

Arguments

Details

The algorithm begins by determining the optimal length of probes required to meet the input constraints while maximizing sensitivity to the target consensus sequence at the specified hybridization formamide concentration. This set of potential target sites is then scored based on the possibility of cross-hybridizing to the other non-target sequences. The set of probes is returned with the minimum possibility of cross-hybridizing.

Value

A data.frame with the optimal set of probes matching the specified constraints. Each row lists the probe's target sequence (name), start position, length in nucleotides, start and end position in the sequence alignment, number of permutations, score, melt point in percent formamide at 42 degrees Celsius, hybridization efficiency (hyb_eff), target site, and probe(s). Probes are designed such that the stringency is determined by the equilibrium hybridization conditions and not subsequent washing steps.

Author(s)

Erik Wright [email protected]

References

ES Wright et al. (2013) Identification of Bacterial and Archaeal Communities From Source to Tap. Water Research Foundation, Denver, CO.

DR Noguera, et al. (2014). Mathematical tools to optimize the design of oligonucleotide probes and primers. Applied Microbiology and Biotechnology. doi:10.1007/s00253-014-6165-x.

See Also

Array2Matrix, NNLS

Examples

fas <- system.file("extdata", "Bacteria_175seqs.fas", package="DECIPHER")
dna <- readDNAStringSet(fas)
names(dna) <- 1:length(dna)
probes <- DesignArray(dna)
probes[1,]

Description

Assists in the design of primer sets targeting a specific group of sequences while minimizing the potential to cross-amplify other groups of sequences.

Usage

DesignPrimers(tiles,
              identifier = "",
              start = 1,
              end = NULL,
              minLength = 17,
              maxLength = 26,
              maxPermutations = 4,
              minCoverage = 0.9,
              minGroupCoverage = 0.2,
              annealingTemp = 64,
              P = 4e-07,
              monovalent = 0.07,
              divalent = 0.003,
              dNTPs = 8e-04,
              minEfficiency = 0.8,
              worstScore = -Inf,
              numPrimerSets = 0,
              minProductSize = 75,
              maxProductSize = 1200,
              maxSearchSize = 1500,
              batchSize = 1000,
              maxDistance = 0.4,
              primerDimer = 1e-07,
              ragged5Prime = TRUE,
              taqEfficiency = TRUE,
              induceMismatch = FALSE,
              processors = 1,
              verbose = TRUE)

Arguments

Details

Primers are designed for use with Taq DNA Polymerase to maximize sensitivity and specificity for the target group of sequences. The design makes use of Taq's bias against certain 3' terminal mismatch types in order to increase specificity further than can be achieve with hybridization efficiency alone.

Primers are designed from a set of tiles to target each identifier while minimizing affinity for all other tiled groups. Arguments provide constraints that ensure the designed primer sets meet the specified criteria as well as being optimized for the particular experimental conditions. A search is conducted through all tiles in the same alignment position to estimate the chance of cross-amplification with a non-target group.

If numPrimers is greater than or equal to one then the set of forward and reverse primers that minimizes potential false positive overlap is returned. This will also initiate a thorough search through all target sites upstream and downstream of the expected binding sites to ensure that the primers do not bind to nearby positions. Lowering the maxSearchSize will speed up the thorough search at the expense of potentially missing an unexpected target site. The number of possible primer sets assessed is increased with the size of numPrimers.

Value

A different data.frame will be returned depending on number of primer sets requested. If no primer sets are required then columns contain the forward and reverse primers for every possible position scored by their potential to amplify other identified groups. If one or more primer sets are requested then columns contain information for the optimal set of forward and reverse primers that could be used in combination to give the fewest potential cross-amplifications.

Note

The program OligoArrayAux (http://www.unafold.org/Dinamelt/software/oligoarrayaux.php) must be installed in a location accessible by the system. For example, the following code should print the installed OligoArrayAux version when executed from the R console:

system("hybrid-min -V")

To install OligoArrayAux from the downloaded source folder on Unix-like platforms, open the shell (or Terminal on Mac OS) and type:

cd oligoarrayaux # change directory to the correct folder name

./configure

make

sudo make install

Author(s)

Erik Wright [email protected]

References

ES Wright et al. (2013) "Exploiting Extension Bias in PCR to Improve Primer Specificity in Ensembles of Nearly Identical DNA Templates." Environmental Microbiology, doi:10.1111/1462-2920.12259.

See Also

AmplifyDNA, CalculateEfficiencyPCR, DesignSignatures, TileSeqs

Examples

if (require("RSQLite", quietly=TRUE)) {
	db <- system.file("extdata", "Bacteria_175seqs.sqlite", package="DECIPHER")
	# not run (must have OligoArrayAux installed first):
	## Not run: tiles <- TileSeqs(db,
		identifier=c("Rhizobiales", "Sphingomonadales"))
## End(Not run)
	## Not run: primers <- DesignPrimers(tiles, identifier="Rhizobiales",
		start=280, end=420, minProductSize=50, numPrimerSets=1)
## End(Not run)
}

Description

Assists in the design of single or dual probes targeting a specific group of sequences while minimizing the potential to cross-hybridize with other groups of sequences.

Usage

DesignProbes(tiles,
             identifier = "",
             start = 1,
             end = NULL,
             minLength = 17,
             maxLength = 26,
             maxPermutations = 4,
             minCoverage = 0.9,
             minGroupCoverage = 0.2,
             hybTemp = 46,
             P = 2.5e-07,
             Na = 1,
             FA = 35,
             minEfficiency = 0.5,
             worstScore = -Inf,
             numProbeSets = 0,
             batchSize = 1000,
             target = "SSU",
             verbose = TRUE)

Arguments

Details

Probes are designed to maximize sensitivity and specificity to the target group(s) (identifier(s)). If numProbeSets > 0 then that many pairs of probes with minimal cross-hybridization overlap are returned, enabling increased specificity with a dual-color approach.

Probes are designed from a set of tiles to target each identifier while minimizing affinity for all other tiled groups. Arguments provide constraints that ensure the designed probes meet the specified criteria as well as being optimized for the particular experimental conditions. A search is conducted through all tiles in the same alignment position to estimate the chance of cross-hybridization with a non-target group.

Two models are used in design, both of which were experimentally calibrated using denaturation profiles from 5 organisms belonging to all three domains of life. Probe lengths are chosen to meet the minEfficiency using a fast model of probe-target hybridization. Candidate probes are then confirmed using a slower model that also takes into account probe-folding and target-folding. Finally, probes are scored for their inability to cross-hybridize with non-target groups by using the fast model and taking into account any mismatches.

Value

A different data.frame will be returned depending on number of primer sets requested. If no probe sets are required then columns contain the designed probes for every possible position scored by their potential to cross-hybridize with other identified groups. If one or more probe sets are requested then columns contain information for the optimal set of probes (probe one and probe two) that could be used in combination to give the fewest potential cross-hybridizations.

Note

The program OligoArrayAux (http://www.unafold.org/Dinamelt/software/oligoarrayaux.php) must be installed in a location accessible by the system. For example, the following code should print the installed OligoArrayAux version when executed from the R console:

system("hybrid-min -V")

To install OligoArrayAux from the downloaded source folder on Unix-like platforms, open the shell (or Terminal on Mac OS) and type:

cd oligoarrayaux # change directory to the correct folder name

./configure

make

sudo make install

Author(s)

Erik Wright [email protected]

References

ES Wright et al. (2014) "Automated Design of Probes for rRNA-Targeted Fluorescence In Situ Hybridization Reveals the Advantages of Using Dual Probes for Accurate Identification." Applied and Environmental Microbiology, doi:10.1128/AEM.01685-14.

See Also

CalculateEfficiencyFISH, TileSeqs

Examples

if (require("RSQLite", quietly=TRUE)) {
	db <- system.file("extdata", "Bacteria_175seqs.sqlite", package="DECIPHER")
	# not run (must have OligoArrayAux installed first):
	## Not run: tiles <- TileSeqs(db,
		identifier=c("Rhizobiales", "Sphingomonadales"))
## End(Not run)
	## Not run: probes <- DesignProbes(tiles, identifier="Rhizobiales",
		start=280, end=420)
## End(Not run)
}

Description

Aids the design of pairs of primers for amplifying a unique “signature” from each group of sequences. Signatures are distinct PCR products that can be differentiated by their length, melt temperature, or sequence.

Usage

DesignSignatures(dbFile,
                 tblName = "Seqs",
                 identifier = "",
                 focusID = NA,
                 type = "melt",
                 resolution = 0.5,
                 levels = 10,
                 enzymes = NULL,
                 minLength = 17,
                 maxLength = 26,
                 maxPermutations = 4,
                 annealingTemp = 64,
                 P = 4e-07,
                 monovalent = 0.07,
                 divalent = 0.003,
                 dNTPs = 8e-04,
                 minEfficiency = 0.8,
                 ampEfficiency = 0.5,
                 numPrimerSets = 100,
                 minProductSize = 70,
                 maxProductSize = 400,
                 kmerSize = 8,
                 searchPrimers = 500,
                 maxDictionary = 20000,
                 primerDimer = 1e-07,
                 pNorm = 1,
                 taqEfficiency = TRUE,
                 processors = 1,
                 verbose = TRUE)

Arguments

Details

Signatures are group-specific PCR products that can be differentiated by either their melt temperature profile, length, or sequence. DesignSignatures assists in finding the optimal pair of forward and reverse primers for obtaining a distinguishable signature from each group of sequences. Groups are delineated by their unique identifier in the database. The algorithm works by progressively narrowing the search for optimal primers: (1) the most frequent k-mers are found; (2) these are used to design primers initially matching the focusID group; (3) the most common forward and reverse primers are selected based on all of the groups, and ambiguity is added up to maxPermutations; (4) a final search is performed to find the optimal forward and reverse primer. Pairs of primers are scored by the distance between the signatures generated for each group, which depends on the type of experiment.

The arguments resolution and levels control the theoretical resolving power of the experiment. The signature for a group is discretized or grouped into “bins” each with a certain magnitude of the signal. Here resolution determines the separation between distinguishable “bins”, and levels controls the range of values in each bin. A high-accuracy experiment would have many bins and/or many levels. While levels is interpreted similarly for every type of experiment, resolution is treated differently depending on type. If type is "melt", then resolution can be either a vector of different melt temperatures, or a single number giving the change in temperatures that can be differentiated. A high-resolution melt (HRM) assay would typically have a resolution between 0.25 and 1 degree Celsius. If type is "length" then resolution is either the number of bins between the minProductSize and maxProductSize, or the bin boundaries. For example, resolution can be lower (wider bins) at long lengths, and higher (narrower bins) at shorter lengths. If type is "sequence" then resolution sets the k-mer size used in differentiating amplicons. Oftentimes, 4 to 6-mers are used for the classification of amplicons.

The signatures can be diversified by using a restriction enzyme to digest the PCR products when type is "melt" or "length". If enzymes are supplied then the an additional search is made to find the best enzyme to use with each pair of primers. In this case, the output includes all of the primer pairs, as well as any enzymes that will digest the PCR products of that primer pair. The output is re-scored to rank the top primer pair and enzyme combination. Note that enzymes is inapplicable when type is "sequence" because restriction enzymes do not alter the sequence of the DNA. Also, it is recommended that only a subset of the available RESTRICTION_ENZYMES are used as input enzymes in order to accelerate the search for the best enzyme.

Value

A data.frame with the top-scoring pairs of forward and reverse primers, their score, the total number of PCR products, and associated columns for the restriction enzyme (if enzyme is not NULL).

Author(s)

Erik Wright [email protected]

References

Wright, E.S. & Vetsigian, K.H. (2016) "DesignSignatures: a tool for designing primers that yields amplicons with distinct signatures." Bioinformatics, doi:10.1093/bioinformatics/btw047.

See Also

AmplifyDNA, CalculateEfficiencyPCR, DesignPrimers, DigestDNA, Disambiguate, MeltDNA, RESTRICTION_ENZYMES

Examples

if (require("RSQLite", quietly=TRUE)) {
	# below are suggested inputs for different types of experiments
	db <- system.file("extdata", "Bacteria_175seqs.sqlite", package="DECIPHER")
	
	## Not run: 
	# High Resolution Melt (HRM) assay:
	primers <- DesignSignatures(db,
	           resolution=seq(75, 100, 0.25), # degrees Celsius
	           minProductSize=55, # base pairs
	           maxProductSize=400)
	
	# Primers for next-generation sequencing:
	primers <- DesignSignatures(db,
	           type="sequence",
	           minProductSize=300, # base pairs
	           maxProductSize=700,
	           resolution=5, # 5-mers
	           levels=5)
	
	# Primers for community fingerprinting:
	primers <- DesignSignatures(db,
	           type="length",
	           levels=2, # presence/absence
	           minProductSize=200, # base pairs
	           maxProductSize=1400,
	           resolution=c(seq(200, 700, 3),
	                        seq(705, 1000, 5),
	                        seq(1010, 1400, 10)))
	
	# Primers for restriction fragment length polymorphism (RFLP):
	data(RESTRICTION_ENZYMES)
	myEnzymes <- RESTRICTION_ENZYMES[c("EcoRI", "HinfI", "SalI")]
	primers <- DesignSignatures(db,
	           type="length",
	           levels=2, # presence/absence
	           minProductSize=200, # base pairs
	           maxProductSize=600,
	           resolution=c(seq(50, 100, 3),
	                        seq(105, 200, 5),
	                        seq(210, 600, 10)),
	           enzymes=myEnzymes)
	
## End(Not run)
}

Description

Detects approximate copies of sequence patterns that likely arose from duplication events and therefore share a common ancestor.

Usage

DetectRepeats(myXStringSet,
              type = "tandem",
              minScore = 10,
              allScores = FALSE,
              maxCopies = 1000,
              maxPeriod = 1000,
              maxFailures = 3,
              maxShifts = 5,
              alphabet = AA_REDUCED[[47]],
              useEmpirical = TRUE,
              correctBackground = TRUE,
              processors = 1,
              verbose = TRUE,
              ...)

Arguments

Details

Many sequences are composed of a substantial fraction of repetitive sequence. Two main forms of repetition are tandem repeats and interspersed repeats, which can be caused by duplication events followed by divergence. Detecting duplications is challenging because of variability in repeat length and composition due to evolution. The significance of a repeat can be quantified by its time since divergence from a common ancestor. DetectRepeats uses a seed-and-extend approach to identify candidate repeats, and tests whether a set of repeats is statistically significant using a background-corrected substitution matrix and gap (i.e., insertion and deletion) penalties. A higher score implies the repeats are more conserved and, therefore, are more likely to have diverged within a finite amount of time from a common ancestor. When myXStringSet is an AAStringSet, similarity includes agreement among predicted secondary structures.

Two possible types of repeats are detectable: * type is "tandem" (the default) or "both" Contiguous approximate copies of a nucleotide or amino acid sequence. First, repeated k-mers are identified along the length of the input sequence(s). Once a k-mer seed is identified, repeated attempts are made to extend the repeat left and right, as well as optimize the beginning and ending positions. * type is "interspersed" or "both" Dispersed approximate copies of a nucleotide sequence on the same strand or opposite strands. These are identified with FindSynteny, aligned with AlignSynteny, and then scored using the same statistical framework as tandem repeats.

The highest scoring repeat in each region is returned, unless allScores is TRUE, in which case overlapping repeats are permitted in the result.

Value

If type is "tandem", a data.frame giving the "Index" of the sequence in myXStringSet, "Begin" and "End" positions of tandem repeats, "Left" and "Right" positions of each repeat, and its "Score".

If type is "interspersed", a data.frame similar to the matrix in the lower diagonal of Synteny objects (see Synteny-class).

If type is "both", a list with the above two elements.

Author(s)

Erik Wright [email protected]

References

Schaper, E., et al. (2012). Repeat or not repeat?-Statistical validation of tandem repeat prediction in genomic sequences. Nucleic Acids Research, 40(20), 10005-17.

See Also

ScoreAlignment

Examples

fas <- system.file("extdata", "Human_huntingtin_cds.fas.gz", package="DECIPHER")
dna <- readDNAStringSet(fas)

x <- DetectRepeats(dna)
x

# number of tandem repeats
lengths(x[, "Left"])

# average periodicity of tandem repeats
per <- mapply(function(a, b) b - a + 1,
	x[, "Left"],
	x[, "Right"],
	SIMPLIFY=FALSE)
sapply(per, mean)

# extract a tandem repeat
i <- 1
reps <- extractAt(dna[[x[i, "Index"]]],
	IRanges(x[[i, "Left"]], x[[i, "Right"]]))
reps
reps <- AlignSeqs(reps, verbose=FALSE) # align the repeats
reps
BrowseSeqs(reps)

# highlight tandem repeats in the sequence
colors <- c("deeppink", "deepskyblue")
colors <- lapply(colors, function(x) col2rgb(x)/255)
cols <- vector("list", length(dna))
for (i in seq_along(cols)) {
	cols[[i]] <- matrix(0, nrow=3, ncol=width(dna)[i])
	for (j in which(x[, "Index"] == i)) {
		left <- x[[j, "Left"]]
		right <- x[[j, "Right"]]
		n <- 0
		for (k in seq_along(left)) {
			r <- left[k]:right[k]
			n <- n + 1
			if (n > length(colors))
				n <- 1
			cols[[i]][, r] <- colors[[n]]
		}
	}
}
BrowseSeqs(dna, patterns=cols)

# find interspersed (rather than tandem) repeats
data(yeastSEQCHR1)
chr1 <- DNAStringSet(yeastSEQCHR1)

if (require("RSQLite", quietly=TRUE)) {
	z <- DetectRepeats(chr1, type="interspersed")
	z
}