pRESTO - The REpertoire Sequencing TOolkit¶

presto.Annotation.annotationConsensus(seq_iter, field, delimiter=('|', '=', ','))¶

Calculate a consensus annotation for a set of sequences

Parameters

seq_iter – an iterator or list of SeqRecord objects
field – the annotation field to take a consensus of
delimiter – a tuple of delimiters for (annotations, field/values, value lists)

Returns

Dictionary with keys: set containing a list of unique annotation values, count containing annotation counts, cons containing the consensus annotation, freq containing the majority annotation frequency

Return type

presto.Annotation.collapseAnnotation(ann_dict, action, fields=None, delimiter=('|', '=', ','))¶

Collapses multiple annotations into new single annotations for each field

Parameters

ann_dict – dictionary of field/value pairs
action – collapse action to take; one of {min, max, sum, first, last, set, cat}
fields – subset of ann_dict to _collapse; if None, collapse all but the ID field
delimiter – Tuple of delimiters for (fields, values, value lists)

Returns

Modified field dictionary

Return type

OrderedDict

presto.Annotation.collapseHeader(header, fields, actions, delimiter=('|', '=', ','))¶

Collapses a sequence header

Parameters

header – an annotation dictionary returned by parseAnnotation.
fields – the list of fields to collapse.
actions – the list of collapse action take; one of (max, min, sum, first, last, set, cat) for each field.
delimiter – a tuple of delimiters for (fields, values, value lists).

Returns

modified header dictionary.

Return type

presto.Annotation.convert454Header(desc)¶

Parses 454 headers into the pRESTO format

Parameters: desc (str) – a sequence description string.
Returns: a dictionary of header field and value pairs.
Return type: dict

Examples

New style 454 header:

@<accession> <length=##>
@GXGJ56Z01AE06X length=222

Old style 454 header:

@<rank_x_y> <length=##> <uaccno=accession>
@000034_0199_0169 length=437 uaccno=GNDG01201ARRCR

presto.Annotation.convertGenbankHeader(desc, delimiter=('|', '=', ','))¶

Converts GenBank and RefSeq headers into the pRESTO format

Parameters

desc (str) – a sequence description string.
delimiter (tuple) – a tuple of delimiters for (fields, values, value lists).

Returns

a dictionary of header field and value pairs.

Return type

Examples

New style GenBank header:

<accession>.<version> <description>
>CM000663.2 Homo sapiens chromosome 1, GRCh38 reference primary assembly

Old style GenBank header:

gi|<GI record number>|<dbsrc>|<accession>.<version>|<description>
>gi|568336023|gb|CM000663.2| Homo sapiens chromosome 1, GRCh38 reference primary assembly

presto.Annotation.convertGenericHeader(desc, delimiter=('|', '=', ','))¶

Converts any header to the pRESTO format

Parameters

desc (str) – a sequence description string.
delimiter (tuple) – a tuple of delimiters for (fields, values, value lists).

Returns

a dictionary of header field and value pairs.

Return type

presto.Annotation.convertIMGTHeader(desc, simple=False)¶

Converts germline headers from IMGT/GENE-DB into the pRESTO format

Parameters

desc (str) – a sequence description string.
simple (bool) – if True then the header will be converted to only the allele name.

Returns

a dictionary of header field and value pairs.

Return type

Examples

IMGT header:

>X60503|IGHV1-18*02|Homo sapiens|F|V-REGION|142..417|276 nt|1| | | | |276+24=300|partial in 3'| |

Header contains 15 fields separated by | (http://imgt.org/genedb):

IMGT/LIGM-DB accession number(s).

Gene and allele name.

Species.

Functionality.

Exon(s), region name(s), or extracted label(s).

Start and end positions in the IMGT/LIGM-DB accession number(s).

Number of nucleotides in the IMGT/LIGM-DB accession number(s).

Codon start, or ‘NR’ (not relevant) for non coding labels and out-of-frame pseudogenes.

Number of nucleotides added in 5' compared to the corresponding label extracted from IMGT/LIGM-DB.

Number of nucleotides added or removed in 3' compared to the corresponding label extracted from IMGT/LIGM-DB.

Number of added, deleted, and/or substituted nucleotides to correct sequencing errors, or ‘not corrected’ if non corrected sequencing errors.

Number of amino acids (AA). This field indicates that the sequence is in amino acids.

Number of characters in the sequence. Nucleotides (or AA) plus IMGT gaps.

Partial (if it is).

Reverse complementary (if it is).

presto.Annotation.convertIlluminaHeader(desc)¶

Converts Illumina headers into the pRESTO format

Parameters: desc (str) – a sequence description string.
Returns: a dictionary of header field and value pairs.
Return type: dict

Examples

New style Illumina header:

@<instrument>:<run number>:<flowcell ID>:<lane>:<tile>:<x-pos>:<y-pos> <read number>:<is filtered>:<control number>:<index sequence>
@MISEQ:132:000000000-A2F3U:1:1101:14340:1555 2:N:0:ATCACG

Old style Illumina header:

@<instrument>:<flowcell lane>:<tile>:<x-pos>:<y-pos>#<index sequence>/<read number>
@HWI-EAS209_0006_FC706VJ:5:58:5894:21141#ATCACG/1
@MS6_33112:1:1101:18371:1066/1

presto.Annotation.convertMIGECHeader(desc)¶

Parses headers from the MIGEC tool into the pRESTO format

Parameters: desc (str) – a sequence description string.
Returns: a dictionary of header field and value pairs.
Return type: dict

Examples

MIGEC header:

@MIG UMI:<UMI sequence>:<consensus read count>
@MIG UMI:TCGGCCAACAAA:8

presto.Annotation.convertSRAHeader(desc)¶

Parses NCBI SRA or EMBL-EBI ENA headers into the pRESTO format

Parameters: desc (str) – a sequence description string.
Returns: a dictionary of header field and value pairs.
Return type: dict

Examples

Header from fastq-dump --split-files:

@<accession>.<spot> <original sequence description> <length=#>
@SRR001666.1 071112_SLXA-EAS1_s_7:5:1:817:345 length=36
@SRR1383326.1 1 length=250

Header from fastq-dump --split-files -I:

@<accession>.<spot>.<read number> <original sequence description> <length=#>
@SRR1383326.1.1 1 length=250

Header from ENA:

@<accession>.<spot> <original sequence description>
@ERR220397.1 HKSQ1MM01DXT2W/3
@ERR346596.1 BS-DSFCONTROL04:4:000000000-A3F0Y:1:1101:12758:1640/1
@ERR346596.1 BS-DSFCONTROL04:4:000000000-A3F0Y:1:1101:12758:1640/2

presto.Annotation.copyHeader(header, fields, names, actions=None, delimiter=('|', '=', ','))¶

Copies fields in a sequence header

Parameters

header – an annotation dictionary returned by parseAnnotation.
fields – a list of the field names to copy.
names – a list of the new field names.
actions – the list of collapse action take after the copy; one of (max, min, sum, first, last, set, cat) for each field.
delimiter – a tuple of delimiters for (fields, values, value lists).

Returns

modified header dictionary.

Return type

presto.Annotation.deleteHeader(header, fields, delimiter=('|', '=', ','))¶

Deletes fields from a sequence header

Parameters

header – an annotation dictionary returned by parseAnnotation.
fields – the list of fields to delete.
delimiter – a tuple of delimiters for (fields, values, value lists).

Returns

modified header dictionary

Return type

presto.Annotation.expandHeader(header, fields, separator=',', delimiter=('|', '=', ','))¶

Splits and annotation value into separate fields in a sequence header

Parameters

header – an annotation dictionary returned by parseAnnotation.
fields – the field to split.
separator – the delimiter to split the values by.
delimiter – a tuple of delimiters for (fields, values, value lists).

Returns

modified header dictionary.

Return type

presto.Annotation.flattenAnnotation(ann_dict, delimiter=('|', '=', ','))¶

Converts annotations from a dictionary to a FASTA/FASTQ sequence description

Parameters

ann_dict – Dictionary of field/value pairs
delimiter – Tuple of delimiters for (fields, values, value lists)

Returns

Formatted sequence description string

Return type

str

presto.Annotation.getAnnotationValues(seq_iter, field, unique=False, delimiter=('|', '=', ','))¶

Gets the set of unique annotation values in a sequence set

Parameters

seq_iter – Iterator or list of SeqRecord objects
field – Annotation field to retrieve values for
unique – If True return a list of only the unique values; if False return a list of all values
delimiter – Tuple of delimiters for (fields, values, value lists)

Returns

List of values for the field

Return type

presto.Annotation.getCoordKey(header, coord_type='presto', delimiter=('|', '=', ','))¶

Return the coordinate identifier for a sequence description

Parameters

header (str) – Sequence header string
coord_type (str) – Sequence header format; one of ‘illumina’, ‘solexa’, ‘sra’, ‘ena’, ‘454’, or ‘presto’; if unrecognized type or None, then return the input header.
delimiter (tuple) – Tuple of delimiters for (fields, values, value lists)

Returns

Coordinate identifier as a string.

Return type

str

presto.Annotation.mergeAnnotation(ann_dict_1, ann_dict_2, prepend=False, delimiter=('|', '=', ','))¶

Merges non-ID field annotations from one field dictionary into another

Parameters

ann_dict_1 – Dictionary of field/value pairs to append to
ann_dict_2 – Dictionary of field/value pairs to merge with ann_dict_2
prepend – If True then add ann_dict_2 values to the front of any ann_dict_1 values that are already present, rather than the default behavior of appending ann_dict_2 values.
delimiter – Tuple of delimiters for (fields, values, value lists)

Returns

Modified ann_dict_1 dictonary of field/value pairs

Return type

OrderedDict

presto.Annotation.mergeHeader(header, fields, name, action=None, delete=False, delimiter=('|', '=', ','))¶

Merges fields in a sequence header

Parameters

header – an annotation dictionary returned by parseAnnotation.
fields – a list of the field names to merge.
name – the name of the new field.
delete – if True delete the merged fields.
actions – the list of collapse action take after the merge one of (max, min, sum, first, last, set, cat).
delimiter – a tuple of delimiters for (fields, values, value lists)

Returns

modified header dictionary.

Return type

presto.Annotation.parseAnnotation(record, fields=None, delimiter=('|', '=', ','))¶

Extracts annotations from a FASTA/FASTQ sequence description

Parameters

record – Description string to extract annotations from
fields – List of fields to subset the return dictionary to; if None return all fields
delimiter – a tuple of delimiters for (fields, values, value lists)

Returns

An OrderedDict of field/value pairs

Return type

OrderedDict

presto.Annotation.parseLog(record)¶

Parses an pRESTO log record

Parameters: record (str) – a string of lines representing a log record including newline characters.
Returns: parsed log contain field and values pairs as a dictionary.
Return type: collections.OrderedDict

presto.Annotation.renameAnnotation(ann_dict, old_field, new_field, delimiter=('|', '=', ','))¶

Renames an annotation and merges annotations if the new name already exists

Parameters

ann_dict – Dictionary of field/value pairs
old_field – Old field name
new_field – New field name
delimiter – Tuple of delimiters for (fields, values, value lists)

Returns

Modified fields dictonary

Return type

OrderedDict

presto.Annotation.renameHeader(header, fields, names, actions=None, delimiter=('|', '=', ','))¶

Renames fields in a sequence header

Parameters

header – an annotation dictionary returned by parseAnnotation.
fields – a list of the current field names.
names – a list of the new field names.
actions – the list of collapse action take after the rename; one of (max, min, sum, first, last, set, cat) for each field.
delimiter – a tuple of delimiters for (fields, values, value lists).

Returns

modified header dictionary.

Return type

presto.Applications¶

External application wrappers

presto.Applications.makeBlastnDb(ref_file, db_exec='makeblastdb')¶

Makes a blastn database file

Parameters

ref_file – the path to the reference database file
db_exec – the path to the makeblastdb executable

Returns

(name and location of the database, handle of the tempfile.TemporaryDirectory)

Return type

tuple

presto.Applications.makeUBlastDb(ref_file, db_exec='usearch')¶

Makes a ublast database file

Parameters

ref_file – path to the reference database file.
db_exec – path to the usearch executable.

Returns

(location of the database, handle of the tempfile.NamedTemporaryFile)

Return type

tuple

presto.Applications.runBlastn(seq, database, evalue=1e-05, max_hits=100, aligner_exec='blastn')¶

Aligns a sequence against a reference database using BLASTN

Parameters

seq – a list of SeqRecord objects to align.
database – the path and name of the blastn database.
evalue – the E-value cut-off.
maxhits – the maximum number of hits returned.
aligner_exec – the path to the blastn executable.

Returns

Alignment results.

Return type

pandas.DataFrame

presto.Applications.runCDHit(seq_list, ident=0.9, length_ratio=0.0, seq_start=0, seq_end=None, max_memory=3000, threads=1, cluster_exec='cd-hit-est', min_word_match=12)¶

Cluster a set of sequences using CD-HIT

Parameters

seq_list (list) – a list of SeqRecord objects to align.
ident (float) – the sequence identity cutoff to be passed to cd-hit-est.
length_ratio (float) – cd-hit-est parameter defining the minimum short/long length ratio allowed within a cluster.
seq_start (int) – the start position to trim sequences at before clustering.
seq_end (int) – the end position to trim sequences at before clustering.
max_memory (int) – cd-hit-est max memory limit (Mb)
threads (int) – number of threads for cd-hit-est.
cluster_exec (str) – the path to the cd-hit-est executable.
min_word_match (int) – (not used for this function, see runUClust)

Returns

{cluster id: list of sequence ids}.

Return type

presto.Applications.runMuscle(seq_list, aligner_exec='muscle')¶

Multiple aligns a set of sequences using MUSCLE

Parameters

seq_list – a list of SeqRecord objects to align
aligner_exec – the MUSCLE executable

Returns

Multiple alignment results.

Return type

Bio.Align.MultipleSeqAlignment

presto.Applications.runUBlast(seq, database, evalue=1e-05, max_hits=100, aligner_exec='usearch')¶

Aligns a sequence against a reference database using the usearch_local algorithm of USEARCH

Parameters

seq – a list of SeqRecord objects to align.
database – the path to the ublast database or a fasta file.
evalue – the E-value cut-off.
maxhits – the maximum number of hits returned.
aligner_exec – the path to the usearch executable.

Returns

Alignment results.

Return type

pandas.DataFrame

presto.Applications.runUClust(seq_list, ident=0.9, length_ratio=0.0, seq_start=0, seq_end=None, max_memory=3000, threads=1, cluster_exec='usearch', min_word_match=12)¶

Cluster a set of sequences using the UCLUST algorithm from USEARCH

Parameters

seq_list (list) – a list of SeqRecord objects to align.
ident (float) – the sequence identity cutoff to be passed to usearch.
length_ratio (float) – usearch parameter defining the minimum short/long length ratio allowed within a cluster.
seq_start (int) – the start position to trim sequences at before clustering.
seq_end (int) – the end position to trim sequences at before clustering.
max_memory (int) – currently ignored.
threads (int) – number of threads for usearch.
cluster_exec (str) – the path to the usearch executable.
min_word_match (int) – minimum number of words that must match to be clustered.

Returns

{cluster id: list of sequence ids}.

Return type

presto.Commandline¶

Commandline interface functions

class presto.Commandline.CommonHelpFormatter(prog, indent_increment=2, max_help_position=24, width=None)¶

Bases: argparse.RawDescriptionHelpFormatter, argparse.ArgumentDefaultsHelpFormatter

Custom argparse.HelpFormatter

presto.Commandline.checkArgs(parser)¶

Checks that arguments have been provided and prints help if they have not.

Parameters: parser – An argparse.ArgumentParser defining the commandline arguments.
Returns: True if arguments are present. Prints help and exits if not.
Return type: boolean

presto.Commandline.getCommonArgParser(seq_in=True, seq_out=True, seq_paired=False, db_in=False, db_out=False, out_file=True, failed=True, log=True, annotation=True, multiproc=False, add_help=True)¶

Defines an ArgumentParser object with common pRESTO arguments

Parameters

seq_in (bool) – if True include sequence input arguments.
seq_out (bool) – if True include sequence output arguments.
seq_paired (bool) – if True defined paired-end sequence input and output arguments.
db_in (bool) – if True include tab-delimited database input arguments.
db_out (bool) – if True include tab-delimited database output arguments.
out_file (bool) – if True add explicit output file name arguments.
failed (bool) – If True include arguments for output of failed results.
log (bool) – If True include log arguments.
annotation (bool) – If True include annotation arguments.
multiproc (bool) – If True include multiprocessing arguments.
add_help (bool) – If True add help and version arguments.

Returns

ArgumentParser object.

Return type

argparse.ArgumentParser

presto.Commandline.parseCommonArgs(args, in_arg=None, in_types=None)¶

Checks common arguments from getCommonArgParser and transforms output options to a dictionary

Parameters

args – Argument Namespace defined by ArgumentParser.parse_args
in_arg – String defining a non-standard input file argument to verify; by default [‘db_files’, ‘seq_files’, ‘seq_files_1’, ‘seq_files_2’, ‘primer_file’] are supported in that order
in_types – List of types (file extensions as strings) to allow for files in file_arg if None do not check type

Returns

Dictionary copy of args with output arguments embedded in the dictionary out_args

Return type

presto.IO¶

File I/O and logging functions

presto.IO.countSeqFile(seq_file)¶

Counts the records in FASTA/FASTQ files

Parameters: seq_file – FASTA or FASTQ file containing sample sequences
Returns: Count of records in the sequence file
Return type: int

presto.IO.countSeqSets(seq_file, field='BARCODE', delimiter=('|', '=', ','))¶

Identifies sets of sequences with the same ID field

Parameters

seq_file – FASTA or FASTQ file containing sample sequences
field – Annotation field containing set IDs
delimiter – Tuple of delimiters for (fields, values, value lists)

Returns

Count of unit set IDs in the sequence file

Return type

presto.IO.getFileType(filename)¶

Determines the type of a file by file extension

Parameters: filename – Filename
Returns: String defining the sequence type for SeqIO operations
Return type: str

presto.IO.getOutputHandle(in_file, out_label=None, out_dir=None, out_name=None, out_type=None)¶

Opens an output file handle

Parameters

in_file – Input filename
out_label – Text to be inserted before the file extension; if None do not add a label
out_type – the file extension of the output file; if None use input file extension
out_dir – the output directory; if None use directory of input file
out_name – the short filename to use for the output file; if None use input file short name

Returns

File handle

Return type

file

presto.IO.printCount(current, step, start_time=None, task=None, end=False)¶

Prints a progress bar to standard out

Parameters

current (int) – count of completed tasks.
step (int) – an int defining the progress increment to print at.
start_time (time.time) – task start time returned by time.time(); if None do not add run time to progress
task (str) – name of task to display.
end (bool) – if True print final log (add newline).

presto.IO.printDebug(message, debug=True, newline=False)¶

Prints a debug message to standard error

Parameters

message (str) – message.
debug (bool) – if True print the message.
newline (bool) – prefix with a newline if True.

presto.IO.printError(message, exit=True, newline=False)¶

Prints an error to standard error and exits

Parameters

message (str) – error message.
exit (bool) – if True exit after the message.
newline (bool) – prefix with a newline if True.

presto.IO.printLog(record, handle=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='UTF-8'>, inset=None)¶

Formats a dictionary into a log string

Parameters

record – a dict or OrderedDict of field names mapping to values.
handle – the file handle to write the log to; if None do not write to file.
inset – minimum field name inset; if None automatically space field names.

Returns

Formatted multi-line string in the log format.

Return type

str

presto.IO.printMessage(message, start_time=None, end=False, width=25)¶

Prints a progress message to standard out

Parameters

message – Current task message
start_time – task start time returned by time.time(); if None do not add run time to progress
end – If True print final message (add newline)
width – Maximum number of characters for messages

presto.IO.printProgress(current, total, step, start_time=None, task=None, end=False)¶

Prints a progress bar to standard out

Parameters

current (int) – count of completed tasks.
total (int) – total task count.
step (float) – float defining the fractional progress increment to print at.
start_time (float) – task start time returned by time.time(); if None do not add run time to progress
task (str) – name of task to display.
end (bool) – if True print final log (add newline).

presto.IO.printWarning(message, newline=False)¶

Prints a warning to standard error

Parameters

message (str) – warning message.
newline (bool) – prefix with a newline if True.

presto.IO.readPrimerFile(primer_file, replace_special=True)¶

Processes primer sequences from file

Parameters

primer_file (str) – name of the FASTA file containing primer sequences.
replace_special (bool) – if True replace whitespace and pRESTO annotation delimiters with an underscore.

Returns

Dictionary mapping primer ID to sequence.

Return type

presto.IO.readReferenceFile(ref_file)¶

Create a dictionary of cleaned and ungapped reference sequences.

Parameters

ref_file – reference sequences in fasta format.

Returns

cleaned and ungapped reference sequences;: with the key as the sequence ID and value as a Bio.SeqRecord for each reference sequence.

Return type

presto.IO.readSeqFile(seq_file, index=False, key_func=None)¶

Reads FASTA/FASTQ files

Parameters

seq_file – FASTA or FASTQ file containing sample sequences
index – If True return a dictionary from SeqIO.index(); if False return an iterator from SeqIO.parse()
key_func – the key_function argument to pass to SeqIO.index if index=True

Returns

an interator of SeqRecords if index=False. A dict if True.

Return type

iter

presto.Multiprocessing¶

Multiprocessing functions

class presto.Multiprocessing.SeqData(id, data)¶

Bases: object

Class defining sequence data objects for worker processes

id¶: unique identifier

data¶: single object or a list of data objects.

valid¶: if True data is suitable for processing.

__bool__()¶

Boolean evaluation

Returns: True if the valid attribute is True
Return type: bool

__len__()¶

Length evaluation

Returns: number of objects in the data attribute.
Return type: int

class presto.Multiprocessing.SeqResult(id, data)¶

Bases: object

Class defining sequence result objects for collector processes

id¶: unique identifier

data¶: single unprocessed object or a list of unprocessed data objects.

results¶: single processed object or a list of processed data objects.

valid¶: if True processing was successful.

log¶: dictionary containing the processing log.

__bool__()¶

Boolean evaluation

Returns: True if the valid attribute is True.
Return type: bool

__len__()¶

Length evaluation

Returns: number of objects in the results attribute.
Return type: int

property data_count¶

Data length

Returns: number of objects in the data attribute.
Return type: int

presto.Multiprocessing.collectPairQueue(alive, result_queue, collect_queue, seq_file_1, seq_file_2, label, out_file=None, out_args={'delimiter': ('|', '=', ','), 'failed': True, 'log_file': None, 'out_dir': None, 'out_name': None, 'out_type': None, 'separator': ','})¶

Pulls from results queue, assembles results and manages log and file IO

Parameters

alive – a multiprocessing.Value boolean controlling whether processing continues; when False function returns.
result_queue – a multiprocessing.Queue holding worker results.
collect_queue – a multiprocessing.Queue holding collector return values.
seq_file_1 – the first sequence file name.
seq_file_2 – the second sequence file name.
label – task label used to tag the output files.
out_file – output file name. Automatically generated from the input file if None.
out_args – common output argument dictionary from parseCommonArgs.

Returns

adds a dictionary of {log: log object, out_files: output file names} to collect_queue.

Return type

None

presto.Multiprocessing.collectSeqQueue(alive, result_queue, collect_queue, seq_file, label, index_field=None, out_file=None, out_args={'delimiter': ('|', '=', ','), 'failed': True, 'log_file': None, 'out_dir': None, 'out_name': None, 'out_type': None, 'separator': ','})¶

Pulls from results queue, assembles results and manages log and file IO

Parameters

alive – a multiprocessing.Value boolean controlling whether processing continues; when False function returns.
result_queue – Multiprocessing.Queue holding worker results.
collect_queue – Multiprocessing.Queue to store collector return values.
seq_file – sample sequence file name.
label – task label used to tag the output files.
out_file – output file name. Automatically generated from the input file if None.
out_args – Common output argument dictionary from parseCommonArgs.
index_field – Field defining set membership for sequence sets if None data queue contained individual records.

Returns

Adds a dictionary with key value pairs to collect_queue containing: ’log’ defining a log object, ‘out_files’ defining the output file names

Return type

None

presto.Multiprocessing.feedPairQueue(alive, data_queue, seq_file_1, seq_file_2, coord_type='presto', delimiter=('|', '=', ','))¶

Feeds the data queue with sequence pairs for processQueue processes

Parameters

alive – a multiprocessing.Value boolean controlling whether processing continues; when False function returns
data_queue – an multiprocessing.Queue to hold data for processing
seq_file_1 – the name of sequence file 1
seq_file_2 – the name of sequence file 2
coord_type – the sequence header format
delimiter – a tuple of delimiters for (fields, values, value lists)

Returns

None

presto.Multiprocessing.feedSeqQueue(alive, data_queue, seq_file, index_func=None, index_args={})¶

Feeds the data queue with SeqRecord objects

Parameters

alive – multiprocessing.Value boolean controlling whether processing continues; when False function returns
data_queue – multiprocessing.Queue to hold data for processing
seq_file – Sequence file to read input from
index_func – Function to use to define sequence sets if None do not index sets and feed individual records
index_args – Dictionary of arguments to pass to index_func

Returns

None

presto.Multiprocessing.manageProcesses(feed_func, work_func, collect_func, feed_args={}, work_args={}, collect_args={}, nproc=None, queue_size=None)¶

Manages feeder, worker and collector processes

Parameters

feed_func (function) – Data Queue feeder function.
work_func (function) – Worker function.
collect_func (function) – Result Queue collector function.
feed_args (dict) – Dictionary of arguments to pass to feed_func.
work_args (dict) – Dictionary of arguments to pass to work_func.
collect_args (dict) – Dictionary of arguments to pass to collect_func.
nproc (int) – Number of processQueue processes; if None defaults to the number of CPUs
queue_size (int) – Maximum size of the argument queue; if None defaults to 2*nproc

Returns

Dictionary of collector results

Return type

presto.Multiprocessing.processSeqQueue(alive, data_queue, result_queue, process_func, process_args={})¶

Pulls from data queue, performs calculations, and feeds results queue

Parameters

alive – multiprocessing.Value boolean controlling whether processing continues; when False function returns
data_queue – multiprocessing.Queue holding data to process
result_queue – multiprocessing.Queue to hold processed results
process_func – function to use for processing sequences
process_args – Dictionary of arguments to pass to process_func

Returns

None

presto.Sequence¶

Sequence processing functions

class presto.Sequence.AssemblyRecord(seq=None)¶

Bases: object

A class defining a paired-end assembly result

property overlap¶

class presto.Sequence.AssemblyStats(n)¶

Bases: object

Class containing p-value and z-score matrices for scoring assemblies

class presto.Sequence.PrimerAlignment(seq=None)¶

Bases: object

A class defining a primer alignment result

__bool__()¶

Boolean evaluation of the alignment

Returns: evaluates to the value of the valid attribute
Return type: int

__len__()¶

Length of alignment

Returns: length of align_seq attribute
Return type: int

presto.Sequence.alignAssembly(head_seq, tail_seq, alpha=1e-05, max_error=0.3, min_len=8, max_len=1000, scan_reverse=False, assembly_stats=None, score_dict={('-', '-'): 0, ('-', '.'): 0, ('-', 'A'): 0, ('-', 'B'): 0, ('-', 'C'): 0, ('-', 'D'): 0, ('-', 'G'): 0, ('-', 'H'): 0, ('-', 'K'): 0, ('-', 'M'): 0, ('-', 'N'): 0, ('-', 'R'): 0, ('-', 'S'): 0, ('-', 'T'): 0, ('-', 'U'): 0, ('-', 'V'): 0, ('-', 'W'): 0, ('-', 'Y'): 0, ('.', '-'): 0, ('.', '.'): 0, ('.', 'A'): 0, ('.', 'B'): 0, ('.', 'C'): 0, ('.', 'D'): 0, ('.', 'G'): 0, ('.', 'H'): 0, ('.', 'K'): 0, ('.', 'M'): 0, ('.', 'N'): 0, ('.', 'R'): 0, ('.', 'S'): 0, ('.', 'T'): 0, ('.', 'U'): 0, ('.', 'V'): 0, ('.', 'W'): 0, ('.', 'Y'): 0, ('A', '-'): 0, ('A', '.'): 0, ('A', 'A'): 1, ('A', 'B'): 0, ('A', 'C'): 0, ('A', 'D'): 1, ('A', 'G'): 0, ('A', 'H'): 1, ('A', 'K'): 0, ('A', 'M'): 1, ('A', 'N'): 1, ('A', 'R'): 1, ('A', 'S'): 0, ('A', 'T'): 0, ('A', 'U'): 0, ('A', 'V'): 1, ('A', 'W'): 1, ('A', 'Y'): 0, ('B', '-'): 0, ('B', '.'): 0, ('B', 'A'): 0, ('B', 'B'): 1, ('B', 'C'): 1, ('B', 'D'): 1, ('B', 'G'): 1, ('B', 'H'): 1, ('B', 'K'): 1, ('B', 'M'): 1, ('B', 'N'): 1, ('B', 'R'): 1, ('B', 'S'): 1, ('B', 'T'): 1, ('B', 'U'): 1, ('B', 'V'): 1, ('B', 'W'): 1, ('B', 'Y'): 1, ('C', '-'): 0, ('C', '.'): 0, ('C', 'A'): 0, ('C', 'B'): 1, ('C', 'C'): 1, ('C', 'D'): 0, ('C', 'G'): 0, ('C', 'H'): 1, ('C', 'K'): 0, ('C', 'M'): 1, ('C', 'N'): 1, ('C', 'R'): 0, ('C', 'S'): 1, ('C', 'T'): 0, ('C', 'U'): 0, ('C', 'V'): 1, ('C', 'W'): 0, ('C', 'Y'): 1, ('D', '-'): 0, ('D', '.'): 0, ('D', 'A'): 1, ('D', 'B'): 1, ('D', 'C'): 0, ('D', 'D'): 1, ('D', 'G'): 1, ('D', 'H'): 1, ('D', 'K'): 1, ('D', 'M'): 1, ('D', 'N'): 1, ('D', 'R'): 1, ('D', 'S'): 1, ('D', 'T'): 1, ('D', 'U'): 1, ('D', 'V'): 1, ('D', 'W'): 1, ('D', 'Y'): 1, ('G', '-'): 0, ('G', '.'): 0, ('G', 'A'): 0, ('G', 'B'): 1, ('G', 'C'): 0, ('G', 'D'): 1, ('G', 'G'): 1, ('G', 'H'): 0, ('G', 'K'): 1, ('G', 'M'): 0, ('G', 'N'): 1, ('G', 'R'): 1, ('G', 'S'): 1, ('G', 'T'): 0, ('G', 'U'): 0, ('G', 'V'): 1, ('G', 'W'): 0, ('G', 'Y'): 0, ('H', '-'): 0, ('H', '.'): 0, ('H', 'A'): 1, ('H', 'B'): 1, ('H', 'C'): 1, ('H', 'D'): 1, ('H', 'G'): 0, ('H', 'H'): 1, ('H', 'K'): 1, ('H', 'M'): 1, ('H', 'N'): 1, ('H', 'R'): 1, ('H', 'S'): 1, ('H', 'T'): 1, ('H', 'U'): 1, ('H', 'V'): 1, ('H', 'W'): 1, ('H', 'Y'): 1, ('K', '-'): 0, ('K', '.'): 0, ('K', 'A'): 0, ('K', 'B'): 1, ('K', 'C'): 0, ('K', 'D'): 1, ('K', 'G'): 1, ('K', 'H'): 1, ('K', 'K'): 1, ('K', 'M'): 0, ('K', 'N'): 1, ('K', 'R'): 1, ('K', 'S'): 1, ('K', 'T'): 1, ('K', 'U'): 1, ('K', 'V'): 1, ('K', 'W'): 1, ('K', 'Y'): 1, ('M', '-'): 0, ('M', '.'): 0, ('M', 'A'): 1, ('M', 'B'): 1, ('M', 'C'): 1, ('M', 'D'): 1, ('M', 'G'): 0, ('M', 'H'): 1, ('M', 'K'): 0, ('M', 'M'): 1, ('M', 'N'): 1, ('M', 'R'): 1, ('M', 'S'): 1, ('M', 'T'): 0, ('M', 'U'): 0, ('M', 'V'): 1, ('M', 'W'): 1, ('M', 'Y'): 1, ('N', '-'): 1, ('N', '.'): 1, ('N', 'A'): 1, ('N', 'B'): 1, ('N', 'C'): 1, ('N', 'D'): 1, ('N', 'G'): 1, ('N', 'H'): 1, ('N', 'K'): 1, ('N', 'M'): 1, ('N', 'N'): 1, ('N', 'R'): 1, ('N', 'S'): 1, ('N', 'T'): 1, ('N', 'U'): 1, ('N', 'V'): 1, ('N', 'W'): 1, ('N', 'Y'): 1, ('R', '-'): 0, ('R', '.'): 0, ('R', 'A'): 1, ('R', 'B'): 1, ('R', 'C'): 0, ('R', 'D'): 1, ('R', 'G'): 1, ('R', 'H'): 1, ('R', 'K'): 1, ('R', 'M'): 1, ('R', 'N'): 1, ('R', 'R'): 1, ('R', 'S'): 1, ('R', 'T'): 0, ('R', 'U'): 0, ('R', 'V'): 1, ('R', 'W'): 1, ('R', 'Y'): 0, ('S', '-'): 0, ('S', '.'): 0, ('S', 'A'): 0, ('S', 'B'): 1, ('S', 'C'): 1, ('S', 'D'): 1, ('S', 'G'): 1, ('S', 'H'): 1, ('S', 'K'): 1, ('S', 'M'): 1, ('S', 'N'): 1, ('S', 'R'): 1, ('S', 'S'): 1, ('S', 'T'): 0, ('S', 'U'): 0, ('S', 'V'): 1, ('S', 'W'): 0, ('S', 'Y'): 1, ('T', '-'): 0, ('T', '.'): 0, ('T', 'A'): 0, ('T', 'B'): 1, ('T', 'C'): 0, ('T', 'D'): 1, ('T', 'G'): 0, ('T', 'H'): 1, ('T', 'K'): 1, ('T', 'M'): 0, ('T', 'N'): 1, ('T', 'R'): 0, ('T', 'S'): 0, ('T', 'T'): 1, ('T', 'U'): 1, ('T', 'V'): 0, ('T', 'W'): 1, ('T', 'Y'): 1, ('U', '-'): 0, ('U', '.'): 0, ('U', 'A'): 0, ('U', 'B'): 1, ('U', 'C'): 0, ('U', 'D'): 1, ('U', 'G'): 0, ('U', 'H'): 1, ('U', 'K'): 1, ('U', 'M'): 0, ('U', 'N'): 1, ('U', 'R'): 0, ('U', 'S'): 0, ('U', 'T'): 1, ('U', 'U'): 1, ('U', 'V'): 0, ('U', 'W'): 1, ('U', 'Y'): 1, ('V', '-'): 0, ('V', '.'): 0, ('V', 'A'): 1, ('V', 'B'): 1, ('V', 'C'): 1, ('V', 'D'): 1, ('V', 'G'): 1, ('V', 'H'): 1, ('V', 'K'): 1, ('V', 'M'): 1, ('V', 'N'): 1, ('V', 'R'): 1, ('V', 'S'): 1, ('V', 'T'): 0, ('V', 'U'): 0, ('V', 'V'): 1, ('V', 'W'): 1, ('V', 'Y'): 1, ('W', '-'): 0, ('W', '.'): 0, ('W', 'A'): 1, ('W', 'B'): 1, ('W', 'C'): 0, ('W', 'D'): 1, ('W', 'G'): 0, ('W', 'H'): 1, ('W', 'K'): 1, ('W', 'M'): 1, ('W', 'N'): 1, ('W', 'R'): 1, ('W', 'S'): 0, ('W', 'T'): 1, ('W', 'U'): 1, ('W', 'V'): 1, ('W', 'W'): 1, ('W', 'Y'): 1, ('Y', '-'): 0, ('Y', '.'): 0, ('Y', 'A'): 0, ('Y', 'B'): 1, ('Y', 'C'): 1, ('Y', 'D'): 1, ('Y', 'G'): 0, ('Y', 'H'): 1, ('Y', 'K'): 1, ('Y', 'M'): 1, ('Y', 'N'): 1, ('Y', 'R'): 0, ('Y', 'S'): 1, ('Y', 'T'): 1, ('Y', 'U'): 1, ('Y', 'V'): 1, ('Y', 'W'): 1, ('Y', 'Y'): 1})¶

Stitches two sequences together by aligning the ends

Parameters

head_seq – the head SeqRecord.
head_seq – the tail SeqRecord.
alpha – the minimum p-value for a valid assembly.
max_error – the maximum error rate for a valid assembly.
min_len – minimum length of overlap to test.
max_len – maximum length of overlap to test.
scan_reverse – if True allow the head sequence to overhang the end of the tail sequence if False end alignment scan at end of tail sequence or start of head sequence.
assembly_stats – optional successes by trials numpy.array of p-values.
score_dict – optional dictionary of character scores in the . form {(char1, char2): score}.

Returns

assembled sequence object.

Return type

presto.Sequence.calculateDiversity(seq_list, score_dict=getDNAScoreDict())¶

Determine the average pairwise error rate for a list of sequences

Parameters

seq_list – List of SeqRecord objects to score
score_dict – Optional dictionary of alignment scores as {(char1, char2): score}

Returns

Average pairwise error rate for the list of sequences

Return type

float

presto.Sequence.calculateSetError(seq_list, ref_seq, ignore_chars=['n', 'N'], score_dict=getDNAScoreDict())¶

Counts the occurrence of nucleotide mismatches from a reference in a set of sequences

Parameters

seq_list – list of SeqRecord objects with aligned sequences.
ref_seq – SeqRecord object containing the reference sequence to match against.
ignore_chars – list of characters to exclude from mismatch counts.
score_dict – optional dictionary of alignment scores as {(char1, char2): score}.

Returns

error rate for the set.

Return type

float

presto.Sequence.checkSeqEqual(seq1, seq2, ignore_chars={'-', '.', 'N', 'n'})¶

Determine if two sequences are equal, excluding missing positions

Parameters

seq1 – SeqRecord object
seq2 – SeqRecord object
ignore_chars – Set of characters to ignore

Returns

True if the sequences are equal

Return type

bool

presto.Sequence.compilePrimers(primers)¶

Translates IUPAC Ambiguous Nucleotide characters to regular expressions and compiles them

Parameters: key – Dictionary of sequences to translate
Returns: Dictionary of compiled regular expressions
Return type: dict

presto.Sequence.consensusUnify(data, field, delimiter=('|', '=', ','))¶

Reassigns all annotations to the consensus annotation in group

Parameters

data – SeqData object contain sequences to process.
field – field containing annotations to collapse.
delimiter – a tuple of delimiters for (annotations, field/values, value lists).

Returns

modified sequences.

Return type

presto.Sequence.deleteSeqPositions(seq, positions)¶

Deletes a list of positions from a SeqRecord

Parameters

seq – SeqRecord objects
positions – Set of positions (indices) to delete

Returns

Modified SeqRecord with the specified positions removed

Return type

SeqRecord

presto.Sequence.deletionUnify(data, field, delimiter=('|', '=', ','))¶

Removes all sequences with differing annotations in a group

Parameters

data – SeqData object contain sequences to process.
field – field containing annotations to collapse.
delimiter – a tuple of delimiters for (annotations, field/values, value lists).

Returns

modified sequences.

Return type

presto.Sequence.PrimerAlignment

presto.Sequence.extractAlignment(seq_record, start, length, rev_primer=False)¶

Extracts a subsequence from sequence

Parameters

data – SeqRecord to process.
start – position where subsequence starts.
length – the length of the subsequence to extract.
rev_primer – if True extract relative to the tail end of the sequence.

Returns

extraction results as an alignment object

Return type

presto.Sequence.filterLength(data, min_length=250, inner=True, missing_chars='-nN.')¶

Filters sequences by length

Parameters

data (SeqData) – a SeqData object with a single SeqRecord to process.
min_length (int) – the minimum length allowed.
inner (bool) – if True exclude outer missing characters from calculation.
missing_chars (str) – a string of missing character values.

Returns

SeqResult object.

Return type

presto.Sequence.filterMissing(data, max_missing=10, inner=True, missing_chars='-nN.')¶

Filters sequences by number of missing nucleotides

Parameters

data (SeqData) – SeqData object with a single SeqRecord to process.
max_missing (int) – the maximum number of allowed ambiguous characters.
inner (bool) – if True exclude outer missing characters from calculation.
missing_chars (str) – a string of missing character values.

Returns

SeqResult object.

Return type

presto.Sequence.filterQuality(data, min_qual=0, inner=True, missing_chars='-nN.')¶

Filters sequences by quality score

Parameters

data (SeqData) – a SeqData object with a single SeqRecord to process.
min_qual (int) – minimum mean quality score for retained sequences.
inner (bool) – if True exclude outer missing characters from calculation.
missing_chars (str) – a string of missing character values.

Returns

SeqResult object.

Return type

presto.Sequence.filterRepeats(data, max_repeat=15, include_missing=False, inner=True, missing_chars='-nN.')¶

Filters sequences by fraction of ambiguous nucleotides

Parameters

data (SeqData) – a SeqData object with a single SeqRecord to process.
max_repeat (int) – the maximum number of allowed repeating characters.
include_missing (int) – if True count ambiguous character repeats; if False do not consider ambiguous character repeats.
inner (int) – if True exclude outer missing characters from calculation.
missing_chars (str) – a string of missing character values.

Returns

SeqResult object.

Return type

presto.Sequence.findGapPositions(seq_list, max_gap, gap_chars={'-', '.'})¶

Finds positions in a set of aligned sequences with a high number of gap characters.

Parameters

seq_list – List of SeqRecord objects with aligned sequences
max_gap – Float of the maximum gap frequency to consider a position as non-gapped
gap_chars – Set of characters to consider as gaps

Returns

Positions (indices) with gap frequency greater than max_gap

Return type

presto.Sequence.frequencyConsensus(seq_list, min_freq=0.6, ignore_chars={'-', '.', 'N', 'n'})¶

Builds a consensus sequence from a set of sequences

Parameters

set_seq – List of SeqRecord objects
min_freq – Frequency cutoff to assign a base
ignore_chars – Set of characters to exclude when building a consensus sequence

Returns

Consensus SeqRecord object

Return type

SeqRecord

presto.Sequence.getAAScoreDict(mask_score=None, gap_score=None)¶

Generates a score dictionary

Parameters

mask_score – Tuple of length two defining scores for all matches against an X character for (a, b), with the score for character (a) taking precedence; if None score symmetrically according to IUPAC character identity
gap_score – Tuple of length two defining score for all matches against a [-, .] character for (a, b), with the score for character (a) taking precedence; if None score symmetrically according to IUPAC character identity

Returns

Score dictionary with keys (char1, char2) mapping to scores

Return type

presto.Sequence.getDNAScoreDict(mask_score=None, gap_score=None)¶

Generates a score dictionary

Parameters

mask_score – Tuple of length two defining scores for all matches against an N character for (a, b), with the score for character (a) taking precedence; if None score symmetrically according to IUPAC character identity
gap_score – Tuple of length two defining score for all matches against a [-, .] character for (a, b), with the score for character (a) taking precedence; if None score symmetrically according to IUPAC character identity

Returns

Score dictionary with keys (char1, char2) mapping to scores

Return type

presto.Sequence.indexSeqSets(seq_dict, field='BARCODE', delimiter=('|', '=', ','))¶

Identifies sets of sequences with the same ID field

Parameters

seq_dict – a dictionary index of sequences returned from SeqIO.index()
field – the annotation field containing set IDs
delimiter – a tuple of delimiters for (fields, values, value lists)

Returns

Dictionary mapping set name to a list of record names

Return type

presto.Sequence.joinAssembly(head_seq, tail_seq, gap=0, insert_seq=None)¶

Concatenates two sequences

Parameters

head_seq – the head SeqRecord.
tail_seq – the tail SeqRecord.
gap – number of gap characters to insert between head and tail ignored if insert_seq is not None.
insert_seq – a string or Bio.Seq.Seq object, to insert between the head and tail; if None insert with N characters.

Returns

assembled sequence object.

Return type

presto.Sequence.PrimerAlignment

presto.Sequence.localAlignment(seq_record, primers, primers_regex=None, max_error=0.3, max_len=1000, rev_primer=False, skip_rc=False, gap_penalty=(1, 1), score_dict={('-', '-'): 0, ('-', '.'): 0, ('-', 'A'): 0, ('-', 'B'): 0, ('-', 'C'): 0, ('-', 'D'): 0, ('-', 'G'): 0, ('-', 'H'): 0, ('-', 'K'): 0, ('-', 'M'): 0, ('-', 'N'): 0, ('-', 'R'): 0, ('-', 'S'): 0, ('-', 'T'): 0, ('-', 'U'): 0, ('-', 'V'): 0, ('-', 'W'): 0, ('-', 'Y'): 0, ('.', '-'): 0, ('.', '.'): 0, ('.', 'A'): 0, ('.', 'B'): 0, ('.', 'C'): 0, ('.', 'D'): 0, ('.', 'G'): 0, ('.', 'H'): 0, ('.', 'K'): 0, ('.', 'M'): 0, ('.', 'N'): 0, ('.', 'R'): 0, ('.', 'S'): 0, ('.', 'T'): 0, ('.', 'U'): 0, ('.', 'V'): 0, ('.', 'W'): 0, ('.', 'Y'): 0, ('A', '-'): 0, ('A', '.'): 0, ('A', 'A'): 1, ('A', 'B'): 0, ('A', 'C'): 0, ('A', 'D'): 1, ('A', 'G'): 0, ('A', 'H'): 1, ('A', 'K'): 0, ('A', 'M'): 1, ('A', 'N'): 1, ('A', 'R'): 1, ('A', 'S'): 0, ('A', 'T'): 0, ('A', 'U'): 0, ('A', 'V'): 1, ('A', 'W'): 1, ('A', 'Y'): 0, ('B', '-'): 0, ('B', '.'): 0, ('B', 'A'): 0, ('B', 'B'): 1, ('B', 'C'): 1, ('B', 'D'): 1, ('B', 'G'): 1, ('B', 'H'): 1, ('B', 'K'): 1, ('B', 'M'): 1, ('B', 'N'): 1, ('B', 'R'): 1, ('B', 'S'): 1, ('B', 'T'): 1, ('B', 'U'): 1, ('B', 'V'): 1, ('B', 'W'): 1, ('B', 'Y'): 1, ('C', '-'): 0, ('C', '.'): 0, ('C', 'A'): 0, ('C', 'B'): 1, ('C', 'C'): 1, ('C', 'D'): 0, ('C', 'G'): 0, ('C', 'H'): 1, ('C', 'K'): 0, ('C', 'M'): 1, ('C', 'N'): 1, ('C', 'R'): 0, ('C', 'S'): 1, ('C', 'T'): 0, ('C', 'U'): 0, ('C', 'V'): 1, ('C', 'W'): 0, ('C', 'Y'): 1, ('D', '-'): 0, ('D', '.'): 0, ('D', 'A'): 1, ('D', 'B'): 1, ('D', 'C'): 0, ('D', 'D'): 1, ('D', 'G'): 1, ('D', 'H'): 1, ('D', 'K'): 1, ('D', 'M'): 1, ('D', 'N'): 1, ('D', 'R'): 1, ('D', 'S'): 1, ('D', 'T'): 1, ('D', 'U'): 1, ('D', 'V'): 1, ('D', 'W'): 1, ('D', 'Y'): 1, ('G', '-'): 0, ('G', '.'): 0, ('G', 'A'): 0, ('G', 'B'): 1, ('G', 'C'): 0, ('G', 'D'): 1, ('G', 'G'): 1, ('G', 'H'): 0, ('G', 'K'): 1, ('G', 'M'): 0, ('G', 'N'): 1, ('G', 'R'): 1, ('G', 'S'): 1, ('G', 'T'): 0, ('G', 'U'): 0, ('G', 'V'): 1, ('G', 'W'): 0, ('G', 'Y'): 0, ('H', '-'): 0, ('H', '.'): 0, ('H', 'A'): 1, ('H', 'B'): 1, ('H', 'C'): 1, ('H', 'D'): 1, ('H', 'G'): 0, ('H', 'H'): 1, ('H', 'K'): 1, ('H', 'M'): 1, ('H', 'N'): 1, ('H', 'R'): 1, ('H', 'S'): 1, ('H', 'T'): 1, ('H', 'U'): 1, ('H', 'V'): 1, ('H', 'W'): 1, ('H', 'Y'): 1, ('K', '-'): 0, ('K', '.'): 0, ('K', 'A'): 0, ('K', 'B'): 1, ('K', 'C'): 0, ('K', 'D'): 1, ('K', 'G'): 1, ('K', 'H'): 1, ('K', 'K'): 1, ('K', 'M'): 0, ('K', 'N'): 1, ('K', 'R'): 1, ('K', 'S'): 1, ('K', 'T'): 1, ('K', 'U'): 1, ('K', 'V'): 1, ('K', 'W'): 1, ('K', 'Y'): 1, ('M', '-'): 0, ('M', '.'): 0, ('M', 'A'): 1, ('M', 'B'): 1, ('M', 'C'): 1, ('M', 'D'): 1, ('M', 'G'): 0, ('M', 'H'): 1, ('M', 'K'): 0, ('M', 'M'): 1, ('M', 'N'): 1, ('M', 'R'): 1, ('M', 'S'): 1, ('M', 'T'): 0, ('M', 'U'): 0, ('M', 'V'): 1, ('M', 'W'): 1, ('M', 'Y'): 1, ('N', '-'): 0, ('N', '.'): 0, ('N', 'A'): 0, ('N', 'B'): 0, ('N', 'C'): 0, ('N', 'D'): 0, ('N', 'G'): 0, ('N', 'H'): 0, ('N', 'K'): 0, ('N', 'M'): 0, ('N', 'N'): 0, ('N', 'R'): 0, ('N', 'S'): 0, ('N', 'T'): 0, ('N', 'U'): 0, ('N', 'V'): 0, ('N', 'W'): 0, ('N', 'Y'): 0, ('R', '-'): 0, ('R', '.'): 0, ('R', 'A'): 1, ('R', 'B'): 1, ('R', 'C'): 0, ('R', 'D'): 1, ('R', 'G'): 1, ('R', 'H'): 1, ('R', 'K'): 1, ('R', 'M'): 1, ('R', 'N'): 1, ('R', 'R'): 1, ('R', 'S'): 1, ('R', 'T'): 0, ('R', 'U'): 0, ('R', 'V'): 1, ('R', 'W'): 1, ('R', 'Y'): 0, ('S', '-'): 0, ('S', '.'): 0, ('S', 'A'): 0, ('S', 'B'): 1, ('S', 'C'): 1, ('S', 'D'): 1, ('S', 'G'): 1, ('S', 'H'): 1, ('S', 'K'): 1, ('S', 'M'): 1, ('S', 'N'): 1, ('S', 'R'): 1, ('S', 'S'): 1, ('S', 'T'): 0, ('S', 'U'): 0, ('S', 'V'): 1, ('S', 'W'): 0, ('S', 'Y'): 1, ('T', '-'): 0, ('T', '.'): 0, ('T', 'A'): 0, ('T', 'B'): 1, ('T', 'C'): 0, ('T', 'D'): 1, ('T', 'G'): 0, ('T', 'H'): 1, ('T', 'K'): 1, ('T', 'M'): 0, ('T', 'N'): 1, ('T', 'R'): 0, ('T', 'S'): 0, ('T', 'T'): 1, ('T', 'U'): 1, ('T', 'V'): 0, ('T', 'W'): 1, ('T', 'Y'): 1, ('U', '-'): 0, ('U', '.'): 0, ('U', 'A'): 0, ('U', 'B'): 1, ('U', 'C'): 0, ('U', 'D'): 1, ('U', 'G'): 0, ('U', 'H'): 1, ('U', 'K'): 1, ('U', 'M'): 0, ('U', 'N'): 1, ('U', 'R'): 0, ('U', 'S'): 0, ('U', 'T'): 1, ('U', 'U'): 1, ('U', 'V'): 0, ('U', 'W'): 1, ('U', 'Y'): 1, ('V', '-'): 0, ('V', '.'): 0, ('V', 'A'): 1, ('V', 'B'): 1, ('V', 'C'): 1, ('V', 'D'): 1, ('V', 'G'): 1, ('V', 'H'): 1, ('V', 'K'): 1, ('V', 'M'): 1, ('V', 'N'): 1, ('V', 'R'): 1, ('V', 'S'): 1, ('V', 'T'): 0, ('V', 'U'): 0, ('V', 'V'): 1, ('V', 'W'): 1, ('V', 'Y'): 1, ('W', '-'): 0, ('W', '.'): 0, ('W', 'A'): 1, ('W', 'B'): 1, ('W', 'C'): 0, ('W', 'D'): 1, ('W', 'G'): 0, ('W', 'H'): 1, ('W', 'K'): 1, ('W', 'M'): 1, ('W', 'N'): 1, ('W', 'R'): 1, ('W', 'S'): 0, ('W', 'T'): 1, ('W', 'U'): 1, ('W', 'V'): 1, ('W', 'W'): 1, ('W', 'Y'): 1, ('Y', '-'): 0, ('Y', '.'): 0, ('Y', 'A'): 0, ('Y', 'B'): 1, ('Y', 'C'): 1, ('Y', 'D'): 1, ('Y', 'G'): 0, ('Y', 'H'): 1, ('Y', 'K'): 1, ('Y', 'M'): 1, ('Y', 'N'): 1, ('Y', 'R'): 0, ('Y', 'S'): 1, ('Y', 'T'): 1, ('Y', 'U'): 1, ('Y', 'V'): 1, ('Y', 'W'): 1, ('Y', 'Y'): 1})¶

Performs pairwise local alignment of a list of short sequences against a long sequence

Parameters

seq_record – a SeqRecord object to align primers against
primers – dictionary of {names: short IUPAC ambiguous sequence strings}
primers_regex – optional dictionary of {names: compiled primer regular expressions}
max_error – maximum acceptable error rate before aligning reverse complement
max_len – maximum length of sample sequence to align
rev_primer – if True align with the tail end of the sequence
skip_rc – if True do not check reverse complement sequences
gap_penalty – a tuple of positive (gap open, gap extend) penalties
score_dict – optional dictionary of alignment scores as {(char1, char2): score}

Returns

primer alignment result object

Return type

presto.Sequence.maskQuality(data, min_qual=0)¶

Masks characters by in sequence by quality score

Parameters

data (SeqData) – a SeqData object with a single SeqRecord to process.
min_qual (int) – minimum quality for retained characters.

Returns

SeqResult object.

Return type

presto.Sequence.maskSeq(align, mode='mask', barcode=False, barcode_field='BARCODE', primer_field='PRIMER', delimiter=('|', '=', ','))¶

Create an output sequence with primers masked or cut

Parameters

align – a PrimerAlignment object.
mode – defines the action taken; one of ‘cut’, ‘mask’, ‘tag’ or ‘trim’.
barcode – if True add sequence preceding primer to description.
barcode_field – name of the output barcode annotation.
primer_field – name of the output primer annotation.
delimiter – a tuple of delimiters for (annotations, field/values, value lists).

Returns

masked sequence.

Return type

Bio.SeqRecord.SeqRecord

presto.Sequence.meanQuality(qual, prob=(1.0, 0.7943282347242815, 0.6309573444801932, 0.5011872336272722, 0.3981071705534972, 0.31622776601683794, 0.251188643150958, 0.19952623149688797, 0.15848931924611134, 0.12589254117941673, 0.1, 0.07943282347242814, 0.06309573444801933, 0.05011872336272722, 0.039810717055349734, 0.03162277660168379, 0.025118864315095794, 0.0199526231496888, 0.015848931924611134, 0.012589254117941675, 0.01, 0.007943282347242814, 0.00630957344480193, 0.005011872336272725, 0.003981071705534973, 0.0031622776601683794, 0.0025118864315095794, 0.001995262314968879, 0.001584893192461114, 0.0012589254117941675, 0.001, 0.0007943282347242813, 0.000630957344480193, 0.0005011872336272725, 0.00039810717055349735, 0.00031622776601683794, 0.00025118864315095795, 0.00019952623149688788, 0.00015848931924611142, 0.00012589254117941674, 0.0001, 7.943282347242822e-05, 6.309573444801929e-05, 5.011872336272725e-05, 3.9810717055349695e-05, 3.1622776601683795e-05, 2.5118864315095822e-05, 1.9952623149688786e-05, 1.584893192461114e-05, 1.2589254117941661e-05, 1e-05, 7.943282347242822e-06, 6.30957344480193e-06, 5.011872336272725e-06, 3.981071705534969e-06, 3.162277660168379e-06, 2.5118864315095823e-06, 1.9952623149688787e-06, 1.584893192461114e-06, 1.2589254117941661e-06, 1e-06, 7.943282347242822e-07, 6.30957344480193e-07, 5.011872336272725e-07, 3.981071705534969e-07, 3.162277660168379e-07, 2.5118864315095823e-07, 1.9952623149688787e-07, 1.584893192461114e-07, 1.2589254117941662e-07, 1e-07, 7.943282347242822e-08, 6.30957344480193e-08, 5.011872336272725e-08, 3.981071705534969e-08, 3.162277660168379e-08, 2.511886431509582e-08, 1.9952623149688786e-08, 1.5848931924611143e-08, 1.2589254117941661e-08, 1e-08, 7.943282347242822e-09, 6.309573444801943e-09, 5.011872336272715e-09, 3.981071705534969e-09, 3.1622776601683795e-09, 2.511886431509582e-09, 1.9952623149688828e-09, 1.584893192461111e-09, 1.2589254117941663e-09, 1e-09, 7.943282347242822e-10, 6.309573444801942e-10, 5.011872336272714e-10, 3.9810717055349694e-10, 3.1622776601683795e-10, 2.511886431509582e-10, 1.9952623149688828e-10, 1.584893192461111e-10, 1.2589254117941662e-10, 1e-10, 7.943282347242822e-11, 6.309573444801942e-11, 5.011872336272715e-11, 3.9810717055349695e-11, 3.1622776601683794e-11, 2.5118864315095823e-11, 1.9952623149688828e-11, 1.5848931924611107e-11, 1.2589254117941662e-11, 1e-11, 7.943282347242821e-12, 6.309573444801943e-12, 5.011872336272715e-12, 3.9810717055349695e-12, 3.1622776601683794e-12, 2.5118864315095823e-12, 1.9952623149688827e-12, 1.584893192461111e-12, 1.258925411794166e-12, 1e-12, 7.943282347242822e-13, 6.309573444801942e-13, 5.011872336272715e-13, 3.981071705534969e-13, 3.162277660168379e-13, 2.511886431509582e-13, 1.9952623149688827e-13))¶

Calculate mean quality score

Parameters

qual (list) – numeric Phred quality scores.
prob (list) – mapping of Phred score (index) to probability values

Returns

floor of the mean Phred quality score.

Return type

presto.Sequence.overlapConsensus(head_seq, tail_seq, ignore_chars={'-', '.', 'N', 'n'})¶

Creates a consensus overlap sequences from two segments

Parameters

head_seq – the overlap head SeqRecord.
tail_seq – the overlap tail SeqRecord.
ignore_chars – list of characters which do not contribute to consensus.

Returns

A SeqRecord object with consensus characters and quality scores.

Return type

SeqRecord

presto.Sequence.qualityConsensus(seq_list, min_qual=0, min_freq=0.6, dependent=False, ignore_chars={'-', '.', 'N', 'n'})¶

Builds a consensus sequence from a set of sequences

Parameters

seq_list – List of SeqRecord objects
min_qual – Quality cutoff to assign a base
min_freq – Frequency cutoff to assign a base
dependent – If False assume sequences are independent for quality calculation
ignore_chars – Set of characters to exclude when building a consensus sequence

Returns

Consensus SeqRecord object

Return type

SeqRecord

presto.Sequence.referenceAssembly(head_seq, tail_seq, ref_dict, ref_db, min_ident=0.5, evalue=1e-05, max_hits=100, fill=False, aligner='usearch', aligner_exec='usearch', score_dict={('-', '-'): 0, ('-', '.'): 0, ('-', 'A'): 0, ('-', 'B'): 0, ('-', 'C'): 0, ('-', 'D'): 0, ('-', 'G'): 0, ('-', 'H'): 0, ('-', 'K'): 0, ('-', 'M'): 0, ('-', 'N'): 0, ('-', 'R'): 0, ('-', 'S'): 0, ('-', 'T'): 0, ('-', 'U'): 0, ('-', 'V'): 0, ('-', 'W'): 0, ('-', 'Y'): 0, ('.', '-'): 0, ('.', '.'): 0, ('.', 'A'): 0, ('.', 'B'): 0, ('.', 'C'): 0, ('.', 'D'): 0, ('.', 'G'): 0, ('.', 'H'): 0, ('.', 'K'): 0, ('.', 'M'): 0, ('.', 'N'): 0, ('.', 'R'): 0, ('.', 'S'): 0, ('.', 'T'): 0, ('.', 'U'): 0, ('.', 'V'): 0, ('.', 'W'): 0, ('.', 'Y'): 0, ('A', '-'): 0, ('A', '.'): 0, ('A', 'A'): 1, ('A', 'B'): 0, ('A', 'C'): 0, ('A', 'D'): 1, ('A', 'G'): 0, ('A', 'H'): 1, ('A', 'K'): 0, ('A', 'M'): 1, ('A', 'N'): 1, ('A', 'R'): 1, ('A', 'S'): 0, ('A', 'T'): 0, ('A', 'U'): 0, ('A', 'V'): 1, ('A', 'W'): 1, ('A', 'Y'): 0, ('B', '-'): 0, ('B', '.'): 0, ('B', 'A'): 0, ('B', 'B'): 1, ('B', 'C'): 1, ('B', 'D'): 1, ('B', 'G'): 1, ('B', 'H'): 1, ('B', 'K'): 1, ('B', 'M'): 1, ('B', 'N'): 1, ('B', 'R'): 1, ('B', 'S'): 1, ('B', 'T'): 1, ('B', 'U'): 1, ('B', 'V'): 1, ('B', 'W'): 1, ('B', 'Y'): 1, ('C', '-'): 0, ('C', '.'): 0, ('C', 'A'): 0, ('C', 'B'): 1, ('C', 'C'): 1, ('C', 'D'): 0, ('C', 'G'): 0, ('C', 'H'): 1, ('C', 'K'): 0, ('C', 'M'): 1, ('C', 'N'): 1, ('C', 'R'): 0, ('C', 'S'): 1, ('C', 'T'): 0, ('C', 'U'): 0, ('C', 'V'): 1, ('C', 'W'): 0, ('C', 'Y'): 1, ('D', '-'): 0, ('D', '.'): 0, ('D', 'A'): 1, ('D', 'B'): 1, ('D', 'C'): 0, ('D', 'D'): 1, ('D', 'G'): 1, ('D', 'H'): 1, ('D', 'K'): 1, ('D', 'M'): 1, ('D', 'N'): 1, ('D', 'R'): 1, ('D', 'S'): 1, ('D', 'T'): 1, ('D', 'U'): 1, ('D', 'V'): 1, ('D', 'W'): 1, ('D', 'Y'): 1, ('G', '-'): 0, ('G', '.'): 0, ('G', 'A'): 0, ('G', 'B'): 1, ('G', 'C'): 0, ('G', 'D'): 1, ('G', 'G'): 1, ('G', 'H'): 0, ('G', 'K'): 1, ('G', 'M'): 0, ('G', 'N'): 1, ('G', 'R'): 1, ('G', 'S'): 1, ('G', 'T'): 0, ('G', 'U'): 0, ('G', 'V'): 1, ('G', 'W'): 0, ('G', 'Y'): 0, ('H', '-'): 0, ('H', '.'): 0, ('H', 'A'): 1, ('H', 'B'): 1, ('H', 'C'): 1, ('H', 'D'): 1, ('H', 'G'): 0, ('H', 'H'): 1, ('H', 'K'): 1, ('H', 'M'): 1, ('H', 'N'): 1, ('H', 'R'): 1, ('H', 'S'): 1, ('H', 'T'): 1, ('H', 'U'): 1, ('H', 'V'): 1, ('H', 'W'): 1, ('H', 'Y'): 1, ('K', '-'): 0, ('K', '.'): 0, ('K', 'A'): 0, ('K', 'B'): 1, ('K', 'C'): 0, ('K', 'D'): 1, ('K', 'G'): 1, ('K', 'H'): 1, ('K', 'K'): 1, ('K', 'M'): 0, ('K', 'N'): 1, ('K', 'R'): 1, ('K', 'S'): 1, ('K', 'T'): 1, ('K', 'U'): 1, ('K', 'V'): 1, ('K', 'W'): 1, ('K', 'Y'): 1, ('M', '-'): 0, ('M', '.'): 0, ('M', 'A'): 1, ('M', 'B'): 1, ('M', 'C'): 1, ('M', 'D'): 1, ('M', 'G'): 0, ('M', 'H'): 1, ('M', 'K'): 0, ('M', 'M'): 1, ('M', 'N'): 1, ('M', 'R'): 1, ('M', 'S'): 1, ('M', 'T'): 0, ('M', 'U'): 0, ('M', 'V'): 1, ('M', 'W'): 1, ('M', 'Y'): 1, ('N', '-'): 1, ('N', '.'): 1, ('N', 'A'): 1, ('N', 'B'): 1, ('N', 'C'): 1, ('N', 'D'): 1, ('N', 'G'): 1, ('N', 'H'): 1, ('N', 'K'): 1, ('N', 'M'): 1, ('N', 'N'): 1, ('N', 'R'): 1, ('N', 'S'): 1, ('N', 'T'): 1, ('N', 'U'): 1, ('N', 'V'): 1, ('N', 'W'): 1, ('N', 'Y'): 1, ('R', '-'): 0, ('R', '.'): 0, ('R', 'A'): 1, ('R', 'B'): 1, ('R', 'C'): 0, ('R', 'D'): 1, ('R', 'G'): 1, ('R', 'H'): 1, ('R', 'K'): 1, ('R', 'M'): 1, ('R', 'N'): 1, ('R', 'R'): 1, ('R', 'S'): 1, ('R', 'T'): 0, ('R', 'U'): 0, ('R', 'V'): 1, ('R', 'W'): 1, ('R', 'Y'): 0, ('S', '-'): 0, ('S', '.'): 0, ('S', 'A'): 0, ('S', 'B'): 1, ('S', 'C'): 1, ('S', 'D'): 1, ('S', 'G'): 1, ('S', 'H'): 1, ('S', 'K'): 1, ('S', 'M'): 1, ('S', 'N'): 1, ('S', 'R'): 1, ('S', 'S'): 1, ('S', 'T'): 0, ('S', 'U'): 0, ('S', 'V'): 1, ('S', 'W'): 0, ('S', 'Y'): 1, ('T', '-'): 0, ('T', '.'): 0, ('T', 'A'): 0, ('T', 'B'): 1, ('T', 'C'): 0, ('T', 'D'): 1, ('T', 'G'): 0, ('T', 'H'): 1, ('T', 'K'): 1, ('T', 'M'): 0, ('T', 'N'): 1, ('T', 'R'): 0, ('T', 'S'): 0, ('T', 'T'): 1, ('T', 'U'): 1, ('T', 'V'): 0, ('T', 'W'): 1, ('T', 'Y'): 1, ('U', '-'): 0, ('U', '.'): 0, ('U', 'A'): 0, ('U', 'B'): 1, ('U', 'C'): 0, ('U', 'D'): 1, ('U', 'G'): 0, ('U', 'H'): 1, ('U', 'K'): 1, ('U', 'M'): 0, ('U', 'N'): 1, ('U', 'R'): 0, ('U', 'S'): 0, ('U', 'T'): 1, ('U', 'U'): 1, ('U', 'V'): 0, ('U', 'W'): 1, ('U', 'Y'): 1, ('V', '-'): 0, ('V', '.'): 0, ('V', 'A'): 1, ('V', 'B'): 1, ('V', 'C'): 1, ('V', 'D'): 1, ('V', 'G'): 1, ('V', 'H'): 1, ('V', 'K'): 1, ('V', 'M'): 1, ('V', 'N'): 1, ('V', 'R'): 1, ('V', 'S'): 1, ('V', 'T'): 0, ('V', 'U'): 0, ('V', 'V'): 1, ('V', 'W'): 1, ('V', 'Y'): 1, ('W', '-'): 0, ('W', '.'): 0, ('W', 'A'): 1, ('W', 'B'): 1, ('W', 'C'): 0, ('W', 'D'): 1, ('W', 'G'): 0, ('W', 'H'): 1, ('W', 'K'): 1, ('W', 'M'): 1, ('W', 'N'): 1, ('W', 'R'): 1, ('W', 'S'): 0, ('W', 'T'): 1, ('W', 'U'): 1, ('W', 'V'): 1, ('W', 'W'): 1, ('W', 'Y'): 1, ('Y', '-'): 0, ('Y', '.'): 0, ('Y', 'A'): 0, ('Y', 'B'): 1, ('Y', 'C'): 1, ('Y', 'D'): 1, ('Y', 'G'): 0, ('Y', 'H'): 1, ('Y', 'K'): 1, ('Y', 'M'): 1, ('Y', 'N'): 1, ('Y', 'R'): 0, ('Y', 'S'): 1, ('Y', 'T'): 1, ('Y', 'U'): 1, ('Y', 'V'): 1, ('Y', 'W'): 1, ('Y', 'Y'): 1})¶

Stitches two sequences together by aligning against a reference database

Parameters

head_seq – the head SeqRecord.
head_seq – the tail SeqRecord.
ref_dict – a dictionary of reference SeqRecord objects.
ref_db – the path and name of the reference database.
min_ident – the minimum identity for a valid assembly.
evalue – the E-value cut-off for ublast.
max_hits – the maxhits output limit for ublast.
fill – if False non-overlapping regions will be assigned Ns; if True non-overlapping regions will be filled with the reference sequence.
aligner – the alignment tool; one of ‘blastn’ or ‘usearch’.
aligner_exec – the path to the alignment tool executable.
score_dict – optional dictionary of character scores in the form {(char1, char2): score}.

Returns

assembled sequence object.

Return type

presto.Sequence.reverseComplement(seq)¶

Takes the reverse complement of a sequence

Parameters: seq – a SeqRecord object, Seq object or string to reverse complement
Returns: Object of the same type as the input with the reverse complement sequence
Return type: Seq

presto.Sequence.scoreAA(a, b, mask_score=None, gap_score=None)¶

Returns the score for a pair of IUPAC Extended Protein characters

Parameters

a – First character
b – Second character
mask_score – Tuple of length two defining scores for all matches against an X character for (a, b), with the score for character (a) taking precedence; if None score symmetrically according to IUPAC character identity
gap_score – Tuple of length two defining score for all matches against a gap (-, .) character for (a, b), with the score for character (a) taking precedence; if None score symmetrically according to IUPAC character identity

Returns

Score for the character pair

Return type

presto.Sequence.PrimerAlignment

presto.Sequence.scoreAlignment(seq_record, primers, start=0, rev_primer=False, score_dict={('-', '-'): 0, ('-', '.'): 0, ('-', 'A'): 0, ('-', 'B'): 0, ('-', 'C'): 0, ('-', 'D'): 0, ('-', 'G'): 0, ('-', 'H'): 0, ('-', 'K'): 0, ('-', 'M'): 0, ('-', 'N'): 0, ('-', 'R'): 0, ('-', 'S'): 0, ('-', 'T'): 0, ('-', 'U'): 0, ('-', 'V'): 0, ('-', 'W'): 0, ('-', 'Y'): 0, ('.', '-'): 0, ('.', '.'): 0, ('.', 'A'): 0, ('.', 'B'): 0, ('.', 'C'): 0, ('.', 'D'): 0, ('.', 'G'): 0, ('.', 'H'): 0, ('.', 'K'): 0, ('.', 'M'): 0, ('.', 'N'): 0, ('.', 'R'): 0, ('.', 'S'): 0, ('.', 'T'): 0, ('.', 'U'): 0, ('.', 'V'): 0, ('.', 'W'): 0, ('.', 'Y'): 0, ('A', '-'): 0, ('A', '.'): 0, ('A', 'A'): 1, ('A', 'B'): 0, ('A', 'C'): 0, ('A', 'D'): 1, ('A', 'G'): 0, ('A', 'H'): 1, ('A', 'K'): 0, ('A', 'M'): 1, ('A', 'N'): 1, ('A', 'R'): 1, ('A', 'S'): 0, ('A', 'T'): 0, ('A', 'U'): 0, ('A', 'V'): 1, ('A', 'W'): 1, ('A', 'Y'): 0, ('B', '-'): 0, ('B', '.'): 0, ('B', 'A'): 0, ('B', 'B'): 1, ('B', 'C'): 1, ('B', 'D'): 1, ('B', 'G'): 1, ('B', 'H'): 1, ('B', 'K'): 1, ('B', 'M'): 1, ('B', 'N'): 1, ('B', 'R'): 1, ('B', 'S'): 1, ('B', 'T'): 1, ('B', 'U'): 1, ('B', 'V'): 1, ('B', 'W'): 1, ('B', 'Y'): 1, ('C', '-'): 0, ('C', '.'): 0, ('C', 'A'): 0, ('C', 'B'): 1, ('C', 'C'): 1, ('C', 'D'): 0, ('C', 'G'): 0, ('C', 'H'): 1, ('C', 'K'): 0, ('C', 'M'): 1, ('C', 'N'): 1, ('C', 'R'): 0, ('C', 'S'): 1, ('C', 'T'): 0, ('C', 'U'): 0, ('C', 'V'): 1, ('C', 'W'): 0, ('C', 'Y'): 1, ('D', '-'): 0, ('D', '.'): 0, ('D', 'A'): 1, ('D', 'B'): 1, ('D', 'C'): 0, ('D', 'D'): 1, ('D', 'G'): 1, ('D', 'H'): 1, ('D', 'K'): 1, ('D', 'M'): 1, ('D', 'N'): 1, ('D', 'R'): 1, ('D', 'S'): 1, ('D', 'T'): 1, ('D', 'U'): 1, ('D', 'V'): 1, ('D', 'W'): 1, ('D', 'Y'): 1, ('G', '-'): 0, ('G', '.'): 0, ('G', 'A'): 0, ('G', 'B'): 1, ('G', 'C'): 0, ('G', 'D'): 1, ('G', 'G'): 1, ('G', 'H'): 0, ('G', 'K'): 1, ('G', 'M'): 0, ('G', 'N'): 1, ('G', 'R'): 1, ('G', 'S'): 1, ('G', 'T'): 0, ('G', 'U'): 0, ('G', 'V'): 1, ('G', 'W'): 0, ('G', 'Y'): 0, ('H', '-'): 0, ('H', '.'): 0, ('H', 'A'): 1, ('H', 'B'): 1, ('H', 'C'): 1, ('H', 'D'): 1, ('H', 'G'): 0, ('H', 'H'): 1, ('H', 'K'): 1, ('H', 'M'): 1, ('H', 'N'): 1, ('H', 'R'): 1, ('H', 'S'): 1, ('H', 'T'): 1, ('H', 'U'): 1, ('H', 'V'): 1, ('H', 'W'): 1, ('H', 'Y'): 1, ('K', '-'): 0, ('K', '.'): 0, ('K', 'A'): 0, ('K', 'B'): 1, ('K', 'C'): 0, ('K', 'D'): 1, ('K', 'G'): 1, ('K', 'H'): 1, ('K', 'K'): 1, ('K', 'M'): 0, ('K', 'N'): 1, ('K', 'R'): 1, ('K', 'S'): 1, ('K', 'T'): 1, ('K', 'U'): 1, ('K', 'V'): 1, ('K', 'W'): 1, ('K', 'Y'): 1, ('M', '-'): 0, ('M', '.'): 0, ('M', 'A'): 1, ('M', 'B'): 1, ('M', 'C'): 1, ('M', 'D'): 1, ('M', 'G'): 0, ('M', 'H'): 1, ('M', 'K'): 0, ('M', 'M'): 1, ('M', 'N'): 1, ('M', 'R'): 1, ('M', 'S'): 1, ('M', 'T'): 0, ('M', 'U'): 0, ('M', 'V'): 1, ('M', 'W'): 1, ('M', 'Y'): 1, ('N', '-'): 0, ('N', '.'): 0, ('N', 'A'): 0, ('N', 'B'): 0, ('N', 'C'): 0, ('N', 'D'): 0, ('N', 'G'): 0, ('N', 'H'): 0, ('N', 'K'): 0, ('N', 'M'): 0, ('N', 'N'): 0, ('N', 'R'): 0, ('N', 'S'): 0, ('N', 'T'): 0, ('N', 'U'): 0, ('N', 'V'): 0, ('N', 'W'): 0, ('N', 'Y'): 0, ('R', '-'): 0, ('R', '.'): 0, ('R', 'A'): 1, ('R', 'B'): 1, ('R', 'C'): 0, ('R', 'D'): 1, ('R', 'G'): 1, ('R', 'H'): 1, ('R', 'K'): 1, ('R', 'M'): 1, ('R', 'N'): 1, ('R', 'R'): 1, ('R', 'S'): 1, ('R', 'T'): 0, ('R', 'U'): 0, ('R', 'V'): 1, ('R', 'W'): 1, ('R', 'Y'): 0, ('S', '-'): 0, ('S', '.'): 0, ('S', 'A'): 0, ('S', 'B'): 1, ('S', 'C'): 1, ('S', 'D'): 1, ('S', 'G'): 1, ('S', 'H'): 1, ('S', 'K'): 1, ('S', 'M'): 1, ('S', 'N'): 1, ('S', 'R'): 1, ('S', 'S'): 1, ('S', 'T'): 0, ('S', 'U'): 0, ('S', 'V'): 1, ('S', 'W'): 0, ('S', 'Y'): 1, ('T', '-'): 0, ('T', '.'): 0, ('T', 'A'): 0, ('T', 'B'): 1, ('T', 'C'): 0, ('T', 'D'): 1, ('T', 'G'): 0, ('T', 'H'): 1, ('T', 'K'): 1, ('T', 'M'): 0, ('T', 'N'): 1, ('T', 'R'): 0, ('T', 'S'): 0, ('T', 'T'): 1, ('T', 'U'): 1, ('T', 'V'): 0, ('T', 'W'): 1, ('T', 'Y'): 1, ('U', '-'): 0, ('U', '.'): 0, ('U', 'A'): 0, ('U', 'B'): 1, ('U', 'C'): 0, ('U', 'D'): 1, ('U', 'G'): 0, ('U', 'H'): 1, ('U', 'K'): 1, ('U', 'M'): 0, ('U', 'N'): 1, ('U', 'R'): 0, ('U', 'S'): 0, ('U', 'T'): 1, ('U', 'U'): 1, ('U', 'V'): 0, ('U', 'W'): 1, ('U', 'Y'): 1, ('V', '-'): 0, ('V', '.'): 0, ('V', 'A'): 1, ('V', 'B'): 1, ('V', 'C'): 1, ('V', 'D'): 1, ('V', 'G'): 1, ('V', 'H'): 1, ('V', 'K'): 1, ('V', 'M'): 1, ('V', 'N'): 1, ('V', 'R'): 1, ('V', 'S'): 1, ('V', 'T'): 0, ('V', 'U'): 0, ('V', 'V'): 1, ('V', 'W'): 1, ('V', 'Y'): 1, ('W', '-'): 0, ('W', '.'): 0, ('W', 'A'): 1, ('W', 'B'): 1, ('W', 'C'): 0, ('W', 'D'): 1, ('W', 'G'): 0, ('W', 'H'): 1, ('W', 'K'): 1, ('W', 'M'): 1, ('W', 'N'): 1, ('W', 'R'): 1, ('W', 'S'): 0, ('W', 'T'): 1, ('W', 'U'): 1, ('W', 'V'): 1, ('W', 'W'): 1, ('W', 'Y'): 1, ('Y', '-'): 0, ('Y', '.'): 0, ('Y', 'A'): 0, ('Y', 'B'): 1, ('Y', 'C'): 1, ('Y', 'D'): 1, ('Y', 'G'): 0, ('Y', 'H'): 1, ('Y', 'K'): 1, ('Y', 'M'): 1, ('Y', 'N'): 1, ('Y', 'R'): 0, ('Y', 'S'): 1, ('Y', 'T'): 1, ('Y', 'U'): 1, ('Y', 'V'): 1, ('Y', 'W'): 1, ('Y', 'Y'): 1})¶

Performs a simple fixed position alignment of primers

Parameters

seq_record – a SeqRecord object to align primers against
primers – dictionary of {names: short IUPAC ambiguous sequence strings}
start – position where primer alignment starts
rev_primer – if True align with the tail end of the sequence
score_dict – optional dictionary of {(char1, char2): score} alignment scores

Returns

primer alignment result object

Return type

presto.Sequence.scoreDNA(a, b, mask_score=None, gap_score=None)¶

Returns the score for a pair of IUPAC ambiguous nucleotide characters

Parameters

a – First character
b – Second character
n_score – Tuple of length two defining scores for all matches against an N character for (a, b), with the score for character (a) taking precedence; if None score symmetrically according to IUPAC character identity
gap_score – Tuple of length two defining score for all matches against a gap (-, .) character for (a, b), with the score for character (a) taking precedence; if None score symmetrically according to IUPAC character identity

Returns

Score for the character pair

Return type

presto.Sequence.scoreSeqPair(seq1, seq2, ignore_chars={}, score_dict={('-', '-'): 1, ('-', '.'): 1, ('-', 'A'): 0, ('-', 'B'): 0, ('-', 'C'): 0, ('-', 'D'): 0, ('-', 'G'): 0, ('-', 'H'): 0, ('-', 'K'): 0, ('-', 'M'): 0, ('-', 'N'): 0, ('-', 'R'): 0, ('-', 'S'): 0, ('-', 'T'): 0, ('-', 'U'): 0, ('-', 'V'): 0, ('-', 'W'): 0, ('-', 'Y'): 0, ('.', '-'): 1, ('.', '.'): 1, ('.', 'A'): 0, ('.', 'B'): 0, ('.', 'C'): 0, ('.', 'D'): 0, ('.', 'G'): 0, ('.', 'H'): 0, ('.', 'K'): 0, ('.', 'M'): 0, ('.', 'N'): 0, ('.', 'R'): 0, ('.', 'S'): 0, ('.', 'T'): 0, ('.', 'U'): 0, ('.', 'V'): 0, ('.', 'W'): 0, ('.', 'Y'): 0, ('A', '-'): 0, ('A', '.'): 0, ('A', 'A'): 1, ('A', 'B'): 0, ('A', 'C'): 0, ('A', 'D'): 1, ('A', 'G'): 0, ('A', 'H'): 1, ('A', 'K'): 0, ('A', 'M'): 1, ('A', 'N'): 1, ('A', 'R'): 1, ('A', 'S'): 0, ('A', 'T'): 0, ('A', 'U'): 0, ('A', 'V'): 1, ('A', 'W'): 1, ('A', 'Y'): 0, ('B', '-'): 0, ('B', '.'): 0, ('B', 'A'): 0, ('B', 'B'): 1, ('B', 'C'): 1, ('B', 'D'): 1, ('B', 'G'): 1, ('B', 'H'): 1, ('B', 'K'): 1, ('B', 'M'): 1, ('B', 'N'): 1, ('B', 'R'): 1, ('B', 'S'): 1, ('B', 'T'): 1, ('B', 'U'): 1, ('B', 'V'): 1, ('B', 'W'): 1, ('B', 'Y'): 1, ('C', '-'): 0, ('C', '.'): 0, ('C', 'A'): 0, ('C', 'B'): 1, ('C', 'C'): 1, ('C', 'D'): 0, ('C', 'G'): 0, ('C', 'H'): 1, ('C', 'K'): 0, ('C', 'M'): 1, ('C', 'N'): 1, ('C', 'R'): 0, ('C', 'S'): 1, ('C', 'T'): 0, ('C', 'U'): 0, ('C', 'V'): 1, ('C', 'W'): 0, ('C', 'Y'): 1, ('D', '-'): 0, ('D', '.'): 0, ('D', 'A'): 1, ('D', 'B'): 1, ('D', 'C'): 0, ('D', 'D'): 1, ('D', 'G'): 1, ('D', 'H'): 1, ('D', 'K'): 1, ('D', 'M'): 1, ('D', 'N'): 1, ('D', 'R'): 1, ('D', 'S'): 1, ('D', 'T'): 1, ('D', 'U'): 1, ('D', 'V'): 1, ('D', 'W'): 1, ('D', 'Y'): 1, ('G', '-'): 0, ('G', '.'): 0, ('G', 'A'): 0, ('G', 'B'): 1, ('G', 'C'): 0, ('G', 'D'): 1, ('G', 'G'): 1, ('G', 'H'): 0, ('G', 'K'): 1, ('G', 'M'): 0, ('G', 'N'): 1, ('G', 'R'): 1, ('G', 'S'): 1, ('G', 'T'): 0, ('G', 'U'): 0, ('G', 'V'): 1, ('G', 'W'): 0, ('G', 'Y'): 0, ('H', '-'): 0, ('H', '.'): 0, ('H', 'A'): 1, ('H', 'B'): 1, ('H', 'C'): 1, ('H', 'D'): 1, ('H', 'G'): 0, ('H', 'H'): 1, ('H', 'K'): 1, ('H', 'M'): 1, ('H', 'N'): 1, ('H', 'R'): 1, ('H', 'S'): 1, ('H', 'T'): 1, ('H', 'U'): 1, ('H', 'V'): 1, ('H', 'W'): 1, ('H', 'Y'): 1, ('K', '-'): 0, ('K', '.'): 0, ('K', 'A'): 0, ('K', 'B'): 1, ('K', 'C'): 0, ('K', 'D'): 1, ('K', 'G'): 1, ('K', 'H'): 1, ('K', 'K'): 1, ('K', 'M'): 0, ('K', 'N'): 1, ('K', 'R'): 1, ('K', 'S'): 1, ('K', 'T'): 1, ('K', 'U'): 1, ('K', 'V'): 1, ('K', 'W'): 1, ('K', 'Y'): 1, ('M', '-'): 0, ('M', '.'): 0, ('M', 'A'): 1, ('M', 'B'): 1, ('M', 'C'): 1, ('M', 'D'): 1, ('M', 'G'): 0, ('M', 'H'): 1, ('M', 'K'): 0, ('M', 'M'): 1, ('M', 'N'): 1, ('M', 'R'): 1, ('M', 'S'): 1, ('M', 'T'): 0, ('M', 'U'): 0, ('M', 'V'): 1, ('M', 'W'): 1, ('M', 'Y'): 1, ('N', '-'): 0, ('N', '.'): 0, ('N', 'A'): 1, ('N', 'B'): 1, ('N', 'C'): 1, ('N', 'D'): 1, ('N', 'G'): 1, ('N', 'H'): 1, ('N', 'K'): 1, ('N', 'M'): 1, ('N', 'N'): 1, ('N', 'R'): 1, ('N', 'S'): 1, ('N', 'T'): 1, ('N', 'U'): 1, ('N', 'V'): 1, ('N', 'W'): 1, ('N', 'Y'): 1, ('R', '-'): 0, ('R', '.'): 0, ('R', 'A'): 1, ('R', 'B'): 1, ('R', 'C'): 0, ('R', 'D'): 1, ('R', 'G'): 1, ('R', 'H'): 1, ('R', 'K'): 1, ('R', 'M'): 1, ('R', 'N'): 1, ('R', 'R'): 1, ('R', 'S'): 1, ('R', 'T'): 0, ('R', 'U'): 0, ('R', 'V'): 1, ('R', 'W'): 1, ('R', 'Y'): 0, ('S', '-'): 0, ('S', '.'): 0, ('S', 'A'): 0, ('S', 'B'): 1, ('S', 'C'): 1, ('S', 'D'): 1, ('S', 'G'): 1, ('S', 'H'): 1, ('S', 'K'): 1, ('S', 'M'): 1, ('S', 'N'): 1, ('S', 'R'): 1, ('S', 'S'): 1, ('S', 'T'): 0, ('S', 'U'): 0, ('S', 'V'): 1, ('S', 'W'): 0, ('S', 'Y'): 1, ('T', '-'): 0, ('T', '.'): 0, ('T', 'A'): 0, ('T', 'B'): 1, ('T', 'C'): 0, ('T', 'D'): 1, ('T', 'G'): 0, ('T', 'H'): 1, ('T', 'K'): 1, ('T', 'M'): 0, ('T', 'N'): 1, ('T', 'R'): 0, ('T', 'S'): 0, ('T', 'T'): 1, ('T', 'U'): 1, ('T', 'V'): 0, ('T', 'W'): 1, ('T', 'Y'): 1, ('U', '-'): 0, ('U', '.'): 0, ('U', 'A'): 0, ('U', 'B'): 1, ('U', 'C'): 0, ('U', 'D'): 1, ('U', 'G'): 0, ('U', 'H'): 1, ('U', 'K'): 1, ('U', 'M'): 0, ('U', 'N'): 1, ('U', 'R'): 0, ('U', 'S'): 0, ('U', 'T'): 1, ('U', 'U'): 1, ('U', 'V'): 0, ('U', 'W'): 1, ('U', 'Y'): 1, ('V', '-'): 0, ('V', '.'): 0, ('V', 'A'): 1, ('V', 'B'): 1, ('V', 'C'): 1, ('V', 'D'): 1, ('V', 'G'): 1, ('V', 'H'): 1, ('V', 'K'): 1, ('V', 'M'): 1, ('V', 'N'): 1, ('V', 'R'): 1, ('V', 'S'): 1, ('V', 'T'): 0, ('V', 'U'): 0, ('V', 'V'): 1, ('V', 'W'): 1, ('V', 'Y'): 1, ('W', '-'): 0, ('W', '.'): 0, ('W', 'A'): 1, ('W', 'B'): 1, ('W', 'C'): 0, ('W', 'D'): 1, ('W', 'G'): 0, ('W', 'H'): 1, ('W', 'K'): 1, ('W', 'M'): 1, ('W', 'N'): 1, ('W', 'R'): 1, ('W', 'S'): 0, ('W', 'T'): 1, ('W', 'U'): 1, ('W', 'V'): 1, ('W', 'W'): 1, ('W', 'Y'): 1, ('Y', '-'): 0, ('Y', '.'): 0, ('Y', 'A'): 0, ('Y', 'B'): 1, ('Y', 'C'): 1, ('Y', 'D'): 1, ('Y', 'G'): 0, ('Y', 'H'): 1, ('Y', 'K'): 1, ('Y', 'M'): 1, ('Y', 'N'): 1, ('Y', 'R'): 0, ('Y', 'S'): 1, ('Y', 'T'): 1, ('Y', 'U'): 1, ('Y', 'V'): 1, ('Y', 'W'): 1, ('Y', 'Y'): 1})¶

Determine the error rate for a pair of sequences

Parameters

seq1 – SeqRecord object
seq2 – SeqRecord object
ignore_chars – Set of characters to ignore when scoring and counting the weight
score_dict – Optional dictionary of alignment scores

Returns

Tuple of the (score, minimum weight, error rate) for the pair of sequences

Return type

Tuple

presto.Sequence.sequentialAssembly(head_seq, tail_seq, ref_dict, ref_db, alpha=1e-05, max_error=0.3, min_len=8, max_len=1000, scan_reverse=False, min_ident=0.5, evalue=1e-05, max_hits=100, fill=False, aligner='usearch', aligner_exec='usearch', assembly_stats=None, score_dict={('-', '-'): 0, ('-', '.'): 0, ('-', 'A'): 0, ('-', 'B'): 0, ('-', 'C'): 0, ('-', 'D'): 0, ('-', 'G'): 0, ('-', 'H'): 0, ('-', 'K'): 0, ('-', 'M'): 0, ('-', 'N'): 0, ('-', 'R'): 0, ('-', 'S'): 0, ('-', 'T'): 0, ('-', 'U'): 0, ('-', 'V'): 0, ('-', 'W'): 0, ('-', 'Y'): 0, ('.', '-'): 0, ('.', '.'): 0, ('.', 'A'): 0, ('.', 'B'): 0, ('.', 'C'): 0, ('.', 'D'): 0, ('.', 'G'): 0, ('.', 'H'): 0, ('.', 'K'): 0, ('.', 'M'): 0, ('.', 'N'): 0, ('.', 'R'): 0, ('.', 'S'): 0, ('.', 'T'): 0, ('.', 'U'): 0, ('.', 'V'): 0, ('.', 'W'): 0, ('.', 'Y'): 0, ('A', '-'): 0, ('A', '.'): 0, ('A', 'A'): 1, ('A', 'B'): 0, ('A', 'C'): 0, ('A', 'D'): 1, ('A', 'G'): 0, ('A', 'H'): 1, ('A', 'K'): 0, ('A', 'M'): 1, ('A', 'N'): 1, ('A', 'R'): 1, ('A', 'S'): 0, ('A', 'T'): 0, ('A', 'U'): 0, ('A', 'V'): 1, ('A', 'W'): 1, ('A', 'Y'): 0, ('B', '-'): 0, ('B', '.'): 0, ('B', 'A'): 0, ('B', 'B'): 1, ('B', 'C'): 1, ('B', 'D'): 1, ('B', 'G'): 1, ('B', 'H'): 1, ('B', 'K'): 1, ('B', 'M'): 1, ('B', 'N'): 1, ('B', 'R'): 1, ('B', 'S'): 1, ('B', 'T'): 1, ('B', 'U'): 1, ('B', 'V'): 1, ('B', 'W'): 1, ('B', 'Y'): 1, ('C', '-'): 0, ('C', '.'): 0, ('C', 'A'): 0, ('C', 'B'): 1, ('C', 'C'): 1, ('C', 'D'): 0, ('C', 'G'): 0, ('C', 'H'): 1, ('C', 'K'): 0, ('C', 'M'): 1, ('C', 'N'): 1, ('C', 'R'): 0, ('C', 'S'): 1, ('C', 'T'): 0, ('C', 'U'): 0, ('C', 'V'): 1, ('C', 'W'): 0, ('C', 'Y'): 1, ('D', '-'): 0, ('D', '.'): 0, ('D', 'A'): 1, ('D', 'B'): 1, ('D', 'C'): 0, ('D', 'D'): 1, ('D', 'G'): 1, ('D', 'H'): 1, ('D', 'K'): 1, ('D', 'M'): 1, ('D', 'N'): 1, ('D', 'R'): 1, ('D', 'S'): 1, ('D', 'T'): 1, ('D', 'U'): 1, ('D', 'V'): 1, ('D', 'W'): 1, ('D', 'Y'): 1, ('G', '-'): 0, ('G', '.'): 0, ('G', 'A'): 0, ('G', 'B'): 1, ('G', 'C'): 0, ('G', 'D'): 1, ('G', 'G'): 1, ('G', 'H'): 0, ('G', 'K'): 1, ('G', 'M'): 0, ('G', 'N'): 1, ('G', 'R'): 1, ('G', 'S'): 1, ('G', 'T'): 0, ('G', 'U'): 0, ('G', 'V'): 1, ('G', 'W'): 0, ('G', 'Y'): 0, ('H', '-'): 0, ('H', '.'): 0, ('H', 'A'): 1, ('H', 'B'): 1, ('H', 'C'): 1, ('H', 'D'): 1, ('H', 'G'): 0, ('H', 'H'): 1, ('H', 'K'): 1, ('H', 'M'): 1, ('H', 'N'): 1, ('H', 'R'): 1, ('H', 'S'): 1, ('H', 'T'): 1, ('H', 'U'): 1, ('H', 'V'): 1, ('H', 'W'): 1, ('H', 'Y'): 1, ('K', '-'): 0, ('K', '.'): 0, ('K', 'A'): 0, ('K', 'B'): 1, ('K', 'C'): 0, ('K', 'D'): 1, ('K', 'G'): 1, ('K', 'H'): 1, ('K', 'K'): 1, ('K', 'M'): 0, ('K', 'N'): 1, ('K', 'R'): 1, ('K', 'S'): 1, ('K', 'T'): 1, ('K', 'U'): 1, ('K', 'V'): 1, ('K', 'W'): 1, ('K', 'Y'): 1, ('M', '-'): 0, ('M', '.'): 0, ('M', 'A'): 1, ('M', 'B'): 1, ('M', 'C'): 1, ('M', 'D'): 1, ('M', 'G'): 0, ('M', 'H'): 1, ('M', 'K'): 0, ('M', 'M'): 1, ('M', 'N'): 1, ('M', 'R'): 1, ('M', 'S'): 1, ('M', 'T'): 0, ('M', 'U'): 0, ('M', 'V'): 1, ('M', 'W'): 1, ('M', 'Y'): 1, ('N', '-'): 1, ('N', '.'): 1, ('N', 'A'): 1, ('N', 'B'): 1, ('N', 'C'): 1, ('N', 'D'): 1, ('N', 'G'): 1, ('N', 'H'): 1, ('N', 'K'): 1, ('N', 'M'): 1, ('N', 'N'): 1, ('N', 'R'): 1, ('N', 'S'): 1, ('N', 'T'): 1, ('N', 'U'): 1, ('N', 'V'): 1, ('N', 'W'): 1, ('N', 'Y'): 1, ('R', '-'): 0, ('R', '.'): 0, ('R', 'A'): 1, ('R', 'B'): 1, ('R', 'C'): 0, ('R', 'D'): 1, ('R', 'G'): 1, ('R', 'H'): 1, ('R', 'K'): 1, ('R', 'M'): 1, ('R', 'N'): 1, ('R', 'R'): 1, ('R', 'S'): 1, ('R', 'T'): 0, ('R', 'U'): 0, ('R', 'V'): 1, ('R', 'W'): 1, ('R', 'Y'): 0, ('S', '-'): 0, ('S', '.'): 0, ('S', 'A'): 0, ('S', 'B'): 1, ('S', 'C'): 1, ('S', 'D'): 1, ('S', 'G'): 1, ('S', 'H'): 1, ('S', 'K'): 1, ('S', 'M'): 1, ('S', 'N'): 1, ('S', 'R'): 1, ('S', 'S'): 1, ('S', 'T'): 0, ('S', 'U'): 0, ('S', 'V'): 1, ('S', 'W'): 0, ('S', 'Y'): 1, ('T', '-'): 0, ('T', '.'): 0, ('T', 'A'): 0, ('T', 'B'): 1, ('T', 'C'): 0, ('T', 'D'): 1, ('T', 'G'): 0, ('T', 'H'): 1, ('T', 'K'): 1, ('T', 'M'): 0, ('T', 'N'): 1, ('T', 'R'): 0, ('T', 'S'): 0, ('T', 'T'): 1, ('T', 'U'): 1, ('T', 'V'): 0, ('T', 'W'): 1, ('T', 'Y'): 1, ('U', '-'): 0, ('U', '.'): 0, ('U', 'A'): 0, ('U', 'B'): 1, ('U', 'C'): 0, ('U', 'D'): 1, ('U', 'G'): 0, ('U', 'H'): 1, ('U', 'K'): 1, ('U', 'M'): 0, ('U', 'N'): 1, ('U', 'R'): 0, ('U', 'S'): 0, ('U', 'T'): 1, ('U', 'U'): 1, ('U', 'V'): 0, ('U', 'W'): 1, ('U', 'Y'): 1, ('V', '-'): 0, ('V', '.'): 0, ('V', 'A'): 1, ('V', 'B'): 1, ('V', 'C'): 1, ('V', 'D'): 1, ('V', 'G'): 1, ('V', 'H'): 1, ('V', 'K'): 1, ('V', 'M'): 1, ('V', 'N'): 1, ('V', 'R'): 1, ('V', 'S'): 1, ('V', 'T'): 0, ('V', 'U'): 0, ('V', 'V'): 1, ('V', 'W'): 1, ('V', 'Y'): 1, ('W', '-'): 0, ('W', '.'): 0, ('W', 'A'): 1, ('W', 'B'): 1, ('W', 'C'): 0, ('W', 'D'): 1, ('W', 'G'): 0, ('W', 'H'): 1, ('W', 'K'): 1, ('W', 'M'): 1, ('W', 'N'): 1, ('W', 'R'): 1, ('W', 'S'): 0, ('W', 'T'): 1, ('W', 'U'): 1, ('W', 'V'): 1, ('W', 'W'): 1, ('W', 'Y'): 1, ('Y', '-'): 0, ('Y', '.'): 0, ('Y', 'A'): 0, ('Y', 'B'): 1, ('Y', 'C'): 1, ('Y', 'D'): 1, ('Y', 'G'): 0, ('Y', 'H'): 1, ('Y', 'K'): 1, ('Y', 'M'): 1, ('Y', 'N'): 1, ('Y', 'R'): 0, ('Y', 'S'): 1, ('Y', 'T'): 1, ('Y', 'U'): 1, ('Y', 'V'): 1, ('Y', 'W'): 1, ('Y', 'Y'): 1})¶

Stitches sequences together by first attempting de novo assembly then falling back to reference guided assembly

Parameters

head_seq – the head SeqRecord
head_seq – the tail SeqRecord
ref_dict – a dictionary of reference SeqRecord objects
ref_db – the path and name of the reference database
alpha – the minimum p-value for a valid de novo assembly
max_error – the maximum error rate for a valid de novo assembly
min_len – minimum length of overlap to test for de novo assembly
max_len – maximum length of overlap to test for de novo assembly
scan_reverse – if True allow the head sequence to overhang the end of the tail sequence in de novo assembly if False end alignment scan at end of tail sequence or start of head sequence
min_ident – the minimum identity for a valid reference guided assembly
evalue – the E-value cut-off for reference guided assembly
max_hits – the maxhits output limit for reference guided assembly
fill – if False non-overlapping regions will be assigned Ns in reference guided assembly; if True non-overlapping regions will be filled with the reference sequence.
aligner – the alignment tool; one of ‘blastn’ or ‘usearch’
aligner_exec – the path to the alignment tool executable
assembly_stats – optional successes by trials numpy.array of p-values
score_dict – optional dictionary of character scores in the form {(char1, char2): score}.

Returns

assembled sequence object.

Return type

presto.Sequence.subsetSeqIndex(seq_dict, field, values, delimiter=('|', '=', ','))¶

Subsets a sequence set by annotation value

Parameters

seq_dict – Dictionary index of sequences returned from SeqIO.index()
field – Annotation field to select keys by
values – List of annotation values that define the retained keys
delimiter – Tuple of delimiters for (annotations, field/values, value lists)

Returns

List of keys

Return type

presto.Sequence.subsetSeqSet(seq_iter, field, values, delimiter=('|', '=', ','))¶

Subsets a sequence set by annotation value

Parameters

seq_iter – Iterator or list of SeqRecord objects
field – Annotation field to select by
values – List of annotation values that define the retained sequences
delimiter – Tuple of delimiters for (annotations, field/values, value lists)

Returns

Modified list of SeqRecord objects

Return type

presto.Sequence.translateAmbigDNA(key)¶

Translates IUPAC ambiguous nucleotide characters to or from character sets

Parameters: key – String or re.search object containing the character set to translate
Returns: Character translation
Return type: str

presto.Sequence.trimQuality(data, min_qual=0, window=10, reverse=False)¶

Cuts sequences using a moving mean quality score

Parameters

data (SeqData) – a SeqData object with a single SeqRecord to process.
min_qual (int) – minimum mean quality to define a cut point.
window (int) – nucleotide window size.
reverse (bool) – if True cut the head of the sequence; if False cut the tail of the sequence

Returns

SeqResult object.

Return type

presto.Sequence.weightSeq(seq, ignore_chars={})¶

Returns the length of a sequencing excluding ignored characters

Parameters

seq – SeqRecord or Seq object
ignore_chars – Set of characters to ignore when counting sequence length

Returns

Sum of the character scores for the sequence

Return type

Jiang, He and Weinstein et al, 2013 Example Files

Roche 454 BCR mRNA with Multiplexed Samples¶

Overview of Experimental Data¶

This example uses publicly available data from:

Lineage structure of the human antibody repertoire in response to influenza vaccination.

Jiang N, He J, and Weinstein JA, et al.

Sci Transl Med. 2013. 5(171):171ra19. doi:10.1126/scitranslmed.3004794.

Which may be downloaded from the NCBI Sequence Read Archive under accession ID: SRX190717. For this example, we will use the first 50,000 sequences of sample 43 (accession: SRR765688), which may downloaded downloaded using fastq-dump from the SRA Toolkit:

fastq-dump -X 50000 SRR765688

Primer and sample barcode (referred to as MID in Jiang, He and Weinstein et al, 2013) sequences are available in the published manuscript. This example assumes that the sample barcodes, forward primers (V-region), and reverse primers (C-region) have been extracted from the manuscript and placed into three corresponding FASTA files.

Read Configuration¶

Schematic of the Roche 454 read configuration. The start of each sequences is labeled with the sample barcode. The following sequence may either be in the forward (top) orientation, proceeding 5’ to 3’ in the direction of the V(D)J reading frame, or the reverse complement orientation (bottom), proceeding in the opposite direction.¶

Example Data¶

We have hosted a small subset of the data (Accession: SRR765688) on the pRESTO website in FASTQ format, with accompanying primer and sample barcode (MID) files. The sample data set and workflow script may be downloaded from here:

Overview of the Workflow¶

The example that follows performs all processing steps to arrive at high-quality unique sequences using this example data set. The workflow is derived into three high level tasks:

Quality control

Sample barcode and primer identification

Deduplication and filtering

A graphical representation of the workflow along with the corresponding sequence of pRESTO commands is shown below.

Flowchart¶

Flowchart of processing steps. Each pRESTO tool is shown as a colored box. The workflow is divided into three primary tasks: (orange) quality control, (green) sample barcode and primer identification, (blue) deduplication and filtering of the repertoire. Grey boxes indicate the initial and final data files. The intermediate files output by each tool are not shown for the sake of brevity.¶

Commands¶

#!/usr/bin/env bash
FilterSeq.py length -s SRR765688.fastq -n 300 --outname S43 --log FSL.log
FilterSeq.py quality -s S43_length-pass.fastq -q 20 --outname S43 --log FSQ.log
MaskPrimers.py score -s S43_quality-pass.fastq -p SRR765688_MIDs.fasta \
    --start 0 --maxerror 0.1 --mode cut --pf MID \
    --outname S43-MID --log MPM.log
MaskPrimers.py align -s S43-MID_primers-pass.fastq -p SRX190717_VPrimers.fasta \
    --maxlen 50 --maxerror 0.3 --mode mask --pf VPRIMER \
    --outname S43-FWD --log MPV.log
MaskPrimers.py align -s S43-FWD_primers-pass.fastq -p SRX190717_CPrimers.fasta \
    --maxlen 50 --maxerror 0.3 --mode cut --pf CPRIMER --revpr --skiprc \
    --outname S43-REV --log MPC.log
CollapseSeq.py -s S43-REV_primers-pass.fastq -n 20 --inner --uf MID CPRIMER \
    --cf VPRIMER --act set --outname S43
SplitSeq.py group -s S43_collapse-unique.fastq -f DUPCOUNT --num 2 --outname S43
ParseHeaders.py table -s S43_atleast-2.fastq -f ID DUPCOUNT MID CPRIMER VPRIMER
ParseLog.py -l FSL.log -f ID LENGTH
ParseLog.py -l FSQ.log -f ID QUALITY
ParseLog.py -l MPM.log MPV.log MPC.log -f ID PRSTART PRIMER ERROR

Greiff et al, 2014 Example Files

Quality control¶

The initial stage of the workflow involves two executions of the FilterSeq.py tool. First, the length subcommand is used to filter reads which are too short to yield full V(D)J sequences using a liberal minimum length requirement of 300 bp (-n 300):

2	FilterSeq.py length -s SRR765688.fastq -n 300 --outname S43 --log FSL.log

Next, the quality subcommand removes sequences having a mean Phred quality score below 20 (-q 20).

FilterSeq.py quality -s S43_length-pass.fastq -q 20 --outname S43 --log FSQ.log

The ParseLog.py tool is then used to extract the results from the two FilterSeq.py log files:

17	ParseLog.py -l FSL.log -f ID LENGTH

To create two tab-delimited files containing the following results for each read:

Field	Description
ID	Sequence name
LENGTH	Sequence length
QUALITY	Quality score

Sample barcode and primer identification¶

Annotation of sample barcodes¶

Following the initial filtering steps, additional filtering is performed with three iterations of the MaskPrimers.py tool based upon the presence of recognized sample barcode (MID), forward primer, and reverse primer sequences. As the orientation and position of the sample barcode is known, the first pass through MaskPrimers uses the faster score subcommand which requires a fixed start position (--start 0) and a low allowable error rate (--maxerror 0.1). The name of the MID will be annotated in the sequence header as the MID field (--pf MID):

MaskPrimers.py score -s S43_quality-pass.fastq -p SRR765688_MIDs.fasta \
    --start 0 --maxerror 0.1 --mode cut --pf MID \
    --outname S43-MID --log MPM.log

Primer masking and annotation¶

The next MaskPrimers.py task uses the align subcommand to identify both the start position of the V-segment primer and correct the orientation of the sequence such that all reads are now oriented in the direction of the V(D)J reading frame), as determined by the orientation of the V-segment primer match. The name of the V-segment primer will be annotated in the sequence header as the VPRIMER field (--pf VPRIMER):

MaskPrimers.py align -s S43-MID_primers-pass.fastq -p SRX190717_VPrimers.fasta \
    --maxlen 50 --maxerror 0.3 --mode mask --pf VPRIMER \
    --outname S43-FWD --log MPV.log

The final MaskPrimers task locates the C-region primer, which is used for isotype assignment of each read. As all sequences are assumed to have been properly oriented by the second MaskPrimers task, the additional arguments MaskPrimers align --revpr and MaskPrimers align --skiprc are added to the third execution. The MaskPrimers align --revpr argument informs the tool that primers sequences should be reverse complemented prior to alignment, and that a match should be searched for (and cut from) the tail end of the sequence. The MaskPrimers align --skiprc argument tells the tool to align against only the forward sequence; meaning, it will not check primer matches against the reverse complement sequence and it will not reorient sequences. The name of the C-region primer will be annotated in the sequence header as the CPRIMER field (--pf CPRIMER):

MaskPrimers.py align -s S43-FWD_primers-pass.fastq -p SRX190717_CPrimers.fasta \
    --maxlen 50 --maxerror 0.3 --mode cut --pf CPRIMER --revpr --skiprc \
    --outname S43-REV --log MPC.log

At this stage, a table of primers and alignment error rates may be generated by executing ParseLog on the log file of each MaskPrimers.py tasks:

19	ParseLog.py -l MPM.log MPV.log MPC.log -f ID PRSTART PRIMER ERROR

Which will contain the following information for each log file:

Field	Description
ID	Sequence name
PRIMER	Primer or sample barcode name
ERROR	Primer match error rate

Deduplication and filtering¶

Removal of duplicate sequences¶

The final stage of the workflow involves two filtering steps to yield unique sequences for each sample barcode. First, the set of unique sequences is identified using the CollapseSeq.py tool, allowing for up to 20 interior N-valued positions (-n 20 and --inner), and requiring that all reads considered duplicated share the same isotype and sample barcode tag (--uf MID CPRIMER). Additionally, the V-region primer annotations of the set of duplicate reads are propagated into the annotation of each retained unique sequence (--cf VPRIMER and --act set):

CollapseSeq.py -s S43-REV_primers-pass.fastq -n 20 --inner --uf MID CPRIMER \
    --cf VPRIMER --act set --outname S43

Filtering to repeated sequences¶

CollapseSeq stores the count of duplicate reads for each sequence in the DUPCOUNT annotation. Following duplicate removal, the data is subset to only those unique sequence with at least two representative reads by using the group subcommand of SplitSeq.py on the count field (-f DUPCOUNT) and specifying a numeric threshold (--num 2):

SplitSeq.py group -s S43_collapse-unique.fastq -f DUPCOUNT --num 2 --outname S43

Creating an annotation table¶

Finally, the annotations, including the sample barcode (MID), duplicate read count (DUPCOUNT), isotype (CPRIMER) and V-region primer (VPRIMER), of the final repertoire are then extracted from the SplitSeq.py output into a tab-delimited file using the table subcommand of ParseHeaders.py:

ParseHeaders.py table -s S43_atleast-2.fastq -f ID DUPCOUNT MID CPRIMER VPRIMER

Note

Optionally, you may split each sample into separate files using the MID annotation and an alternate invocation of SplitSeq.py. The group subcommand may be used to split files on a categorical field, rather than a numerical field, by skipping the --num argument:

SplitSeq.py group -s M1_collapse-unique.fastq -f MID

Will split the unique sequence file into a set of separate files according the the valud in the MID field (-f MID), such that each file will contain sequences from only one sample.

Output files¶

The final set of sequences, which serve as input to a V(D)J reference aligner (Eg, IMGT/HighV-QUEST or IgBLAST), and tables that can be plotted for quality control are:

File	Description
S43_collapse-unique.fastq	Total unique sequences
S43_atleast-2.fastq	Unique sequences represented by at least 2 reads
S43_atleast-2_headers.tab	Annotation table of the atleast-2 file
FSL_table.tab	Table of the FilterSeq-length log
FSQ_table.tab	Table of the FilterSeq-quality log
MPM_table.tab	Table of the MID MaskPrimers log
MPV_table.tab	Table of the V-segment MaskPrimers log
MPC_table.tab	Table of the C-region MaskPrimers log

A number of other intermediate and log files are generated during the workflow, which allows easy tracking/reversion of processing steps. These files are not listed in the table above.

Performance¶

Example performance statistics for a comparable, but larger, 454 workflow are presented below. Performance was measured on a 64-core system with 2.3GHz AMD Opteron(TM) 6276 processors and 512GB of RAM, with memory usage measured at peak utilization. The data set contained 1,346,039 raw reads, and required matching of 11 sample barcodes, 11 V-segment primers, and 5 C-region primers.

Line	Tool	Reads	Cores	MB	Minutes
01	FilterSeq.py length	1,346,039	10	669	12.2
02	FilterSeq.py quality	1,184,905	10	619	11.5
03	MaskPrimers.py score	1,184,219	10	621	18.5
04	MaskPrimers.py align	1,177,245	10	618	263.1
05	MaskPrimers.py align	959,920	10	563	122.9
07	ParseHeaders.py expand	548,658	1	329	3.6
08	ParseHeaders.py rename	548,658	1	329	4.1
09	CollapseSeq.py	548,658	1	1,509	5.6

Illumina MiSeq 2x250 BCR mRNA¶

Overview of Experimental Data¶

This example uses publicly available data from:

Quantitative assessment of the robustness of next-generation sequencing of antibody variable gene repertoires from immunized mice.

Greiff, V. et al.

BMC Immunol. 2014. 15(1):40. doi:10.1186/s12865-014-0040-5.

Which may be downloaded from the EBI European Nucleotide Archive under accession ID: ERP003950. For this example, we will use the first 25,000 sequences of sample Replicate-1-M1 (accession: ERR346600), which may downloaded using fastq-dump from the SRA Toolkit:

fastq-dump --split-files -X 25000 ERR346600

Primers sequences are available in additional file 1 of the publication.

Read Configuration¶

Schematic of Illumina MiSeq 2x250 stranded paired-end reads without UMI barcodes. Each 250 base-pair read was sequenced from one end of the target cDNA, so that the two reads together cover the entire variable region of the Ig heavy chain. The V(D)J reading frame proceeds from the start of read 2 to the start of read 1. Read 1 is in the opposite orientation (reverse complement), and contains the C-region primer sequence. Both reads begin with a random sequences of 4 nucleotides.¶

Example Data¶

We have hosted a small subset of the data (Accession: ERR346600) on the pRESTO website in FASTQ format, with accompanying primer files and an example workflow script. The sample data set and workflow script may be downloaded from here:

Overview of the Workflow¶

In the following sections, we demonstrate each step of the workflow to move from raw sequence reads to a fully annotated repertoire of complete V(D)J sequences. The workflow is divided into four high-level tasks:

Paired-end assembly

Quality control and primer annotation

Deduplication and filtering

A graphical representation of the workflow along with the corresponding sequence of pRESTO commands is shown below.

Flowchart¶

Flowchart of processing steps. Each pRESTO tool is shown as a colored box. The workflow is divided into three primary tasks: (orange) paired-end assembly, (green) quality control and primer annotation, and deduplication and filtering (blue). The intermediate files output by each tool are not shown for the sake of brevity.¶

Commands¶

#!/usr/bin/env bash
AssemblePairs.py align -1 ERR346600_2.fastq -2 ERR346600_1.fastq \
    --coord sra --rc tail --outname M1 --log AP.log
FilterSeq.py quality -s M1_assemble-pass.fastq -q 20 --outname M1 --log FS.log
MaskPrimers.py score -s M1_quality-pass.fastq -p Greiff2014_VPrimers.fasta \
    --start 4 --mode mask --pf VPRIMER --outname M1-FWD --log MPV.log
MaskPrimers.py score -s M1-FWD_primers-pass.fastq -p Greiff2014_CPrimers.fasta \
    --start 4 --mode cut --revpr --pf CPRIMER --outname M1-REV --log MPC.log
CollapseSeq.py -s M1-REV_primers-pass.fastq -n 20 --inner --uf CPRIMER \
    --cf VPRIMER --act set --outname M1
SplitSeq.py group -s M1_collapse-unique.fastq -f DUPCOUNT --num 2 --outname M1
ParseHeaders.py table -s M1_atleast-2.fastq -f ID DUPCOUNT CPRIMER VPRIMER
ParseLog.py -l AP.log -f ID LENGTH OVERLAP ERROR PVALUE
ParseLog.py -l FS.log -f ID QUALITY
ParseLog.py -l MPV.log MPC.log -f ID PRIMER ERROR

Stern, Yaari and Vander Heiden et al, 2014 Example Files

Paired-end assembly¶

Each set of paired-ends mate-pairs is first assembled into a full length Ig sequence using the align subcommand of the AssemblePairs.py tool:

AssemblePairs.py align -1 ERR346600_2.fastq -2 ERR346600_1.fastq \
    --coord sra --rc tail --outname M1 --log AP.log

During assembly we have defined read 2 (V-region) as the head of the sequence (-1) and read 1 as the tail of the sequence (-2). The --coord argument defines the format of the sequence header so that AssemblePairs.py can properly identify mate-pairs; in this case, we use --coord sra as our headers are in the SRA/ENA format.

Note

For both the AssemblePairs.py and PairSeq.py commands using the correct --coord argument is critical for matching mate-pairs. If this was raw data from Illumina, rather than data downloaded from SRA/ENA, then the appropriate argument would be --coord illumina.

The ParseLog.py tool is then used to build a tab-delimited file of results from the AssemblePairs.py log file:

13	ParseLog.py -l AP.log -f ID LENGTH OVERLAP ERROR PVALUE

Which will containing the following columns:

Field	Description
ID	Sequence name
LENGTH	Length of the assembled sequence
OVERLAP	Length of the overlap between mate-pairs
ERROR	Mismatch rate of the overlapping region
PVALUE	P-value for the assembly

See also

Depending on the amplicon length in your data, not all mate-pairs may overlap. For the sake of simplicity, we have excluded a demonstration of assembly in such cases. pRESTO provides a couple approaches to deal with such reads. The reference subcommand of AssemblePairs.py can use the ungapped V-segment reference sequences to properly space non-overlapping reads. Or, if all else fails, the join subcommand can be used to simply stick mate-pairs together end-to-end with some intervening gap.

Quality control and primer annotation¶

Removal of low quality reads¶

Quality control begins with the identification and removal of low-quality reads using the quality subcommand of the FilterSeq.py tool. In this example, reads with mean Phred quality scores less than 20 (-q 20) are removed:

FilterSeq.py quality -s M1_assemble-pass.fastq -q 20 --outname M1 --log FS.log

The ParseLog.py tool is then used to build tab-delimited file from the FilterSeq.py log:

14	ParseLog.py -l FS.log -f ID QUALITY

Capturing the following annotations:

Field	Description
ID	Sequence name
QUALITY	Quality score

Read annotation and masking of primer regions¶

When dealing with Ig sequences, it is important to cut or mask the primers, as B cell receptors are subject to somatic hypermutation (the accumulation of point mutations in the DNA) and degenerate primer matches can look like mutations in downstream applications. The score subcommand of MaskPrimers.py is used to identify and remove the V-segment and C-region PCR primers for both reads:

MaskPrimers.py score -s M1_quality-pass.fastq -p Greiff2014_VPrimers.fasta \
    --start 4 --mode mask --pf VPRIMER --outname M1-FWD --log MPV.log
MaskPrimers.py score -s M1-FWD_primers-pass.fastq -p Greiff2014_CPrimers.fasta \
    --start 4 --mode cut --revpr --pf CPRIMER --outname M1-REV --log MPC.log

In this data set the authors have added a random sequence of 4 bp to the start of each read before the primer sequence to increase sequence diversity and the reliability of cluster calling on the Illumina platform. As such, both primers begin at position 4 (--start 4), but the C-region primer begins 4 bases from the end of the assembled read. The addition of the --revpr argument to the second MaskPrimers.py step instructs the tool to reverse complement the primer sequences and check the tail of the read. The two primer regions have also been treated differently. The V-segment primer has been masked (replaced by Ns) using the --mode mask argument to preserve the V(D)J length, while the C-region primer has been removed from the sequence using the --mode cut argument.

During each iteration of the MaskPrimers.py tool an annotation field was added with the name of the primer, with the name of the field specified by the --pf argument, such that after both iterations each sequence contains an annotation of the form:

VPRIMER=V-segment primer|CPRIMER=C-region primer

To summarize these steps, the ParseLog.py tool is used to build tab-delimited files from the two MaskPrimers.py logs:

15	ParseLog.py -l MPV.log MPC.log -f ID PRIMER ERROR

Capturing the following annotations:

Field	Description
ID	Sequence name
PRIMER	Primer name
ERROR	Primer match error rate

Note

This library was prepared in a stranded manner. Meaning, the read orientation is constant for all reads; read 1 is always the C-region end of the amplicon and read 2 is always the V-segment end. If your data is unstranded (50% of the reads are forward, 50% are reversed), then you must modify the first MaskPrimers.py step to account for this by using the align subcommand instead:

MaskPrimers.py align -s M1*quality-pass.fastq -p Greiff2014_VPrimers.fasta \
    --maxlen 30 --mode mask --pf VPRIMER --log MP1.log

This will perform a slower process of locally aligning the primers, checking the reverse compliment of each read for matches, and correcting the the output sequences to the forward orientation (V to J).

Deduplication and filtering¶

Removal of duplicate sequences¶

The last stage of the workflow involves two filtering steps to yield the final repertoire. First, the set of unique sequences is identified using the CollapseSeq.py tool, allowing for up to 20 interior N-valued positions (-n 20 and --inner), and requiring that all reads considered duplicates share the same C-region primer annotation (--uf CPRIMER). Additionally, the V-segment primer annotations of the set of duplicate reads are propagated into the annotation of each retained unique sequence (--cf VPRIMER and --act set):

CollapseSeq.py -s M1-REV_primers-pass.fastq -n 20 --inner --uf CPRIMER \
    --cf VPRIMER --act set --outname M1

Filtering to repeated sequences¶

CollapseSeq stores the count of duplicate reads for each sequence in the DUPCOUNT annotation. Following duplicate removal, the data is subset to only those unique sequence with at least two representative reads by using the group subcommand of SplitSeq.py on the count field (-f DUPCOUNT) and specifying a numeric threshold (--num 2):

SplitSeq.py group -s M1_collapse-unique.fastq -f DUPCOUNT --num 2 --outname M1

Creating an annotation table¶

Finally, the annotations, including duplicate read count (DUPCOUNT), isotype (CPRIMER) and V-segment primer (VPRIMER), of the final repertoire are then extracted from the SplitSeq.py output into a tab-delimited file using the table subcommand of ParseHeaders.py:

12	ParseHeaders.py table -s M1_atleast-2.fastq -f ID DUPCOUNT CPRIMER VPRIMER

Output files¶

The final set of sequences, which serve as input to a V(D)J reference aligner (Eg, IMGT/HighV-QUEST or IgBLAST), and tables that can be plotted for quality control are:

File	Description
M1_collapse-unique.fastq	Total unique sequences
M1_atleast-2.fastq	Unique sequences represented by at least 2 reads
M1_atleast-2_headers.tab	Annotation table of the atleast-2 file
AP_table.tab	Table of the AssemblePairs log
FS_table.tab	Table of the FilterSeq log
MPV_table.tab	Table of the V-segment MaskPrimers log
MPC_table.tab	Table of the C-region MaskPrimers log

A number of other intermediate and log files are generated during the workflow, which allows easy tracking/reversion of processing steps. These files are not listed in the table above.

UMI Barcoded Illumina MiSeq 2x250 BCR mRNA¶

Overview of Experimental Data¶

This example uses publicly available data from:

B cells populating the multiple sclerosis brain mature in the draining cervical lymph nodes.

Stern JNH, Yaari G, and Vander Heiden JA, et al.

Sci Transl Med. 2014. 6(248):248ra107. doi:10.1126/scitranslmed.3008879.

Which may be downloaded from the NCBI Sequence Read Archive under BioProject accession ID: PRJNA248475. For this example, we will use the first 25,000 sequences of sample M12 (accession: SRR1383456), which may downloaded downloaded using fastq-dump from the SRA Toolkit:

fastq-dump --split-files -X 25000 SRR1383456

Primers sequences are available online at the supplemental website for the publication.

Read Configuration¶

Schematic of the Illumina MiSeq 2x250 paired-end reads with UMI barcodes. Each 250 base-pair read was sequenced from one end of the target cDNA, so that the two reads together cover the entire variable region of the Ig heavy chain. The V(D)J reading frame proceeds from the start of read 2 to the start of read 1. Read 1 is in the opposite orientation (reverse complement), and contains a 15 nucleotide UMI barcode preceding the C-region primer sequence.¶

Example Data¶

We have hosted a small subset of the data (Accession: SRR1383456) on the pRESTO website in FASTQ format with accompanying primer files. The sample data set and workflow script may be downloaded from here:

Overview of the Workflow¶

In the following sections, we demonstrate each step of the workflow to move from raw sequence reads to a fully annotated repertoire of complete V(D)J sequences. The workflow is divided into four high-level tasks:

Quality control, UMI annotation and primer masking

Generation of UMI consensus sequences

Paired-end assembly of UMI consensus sequences

Deduplication and filtering

A graphical representation of the workflow along with the corresponding sequence of pRESTO commands is shown below.

Flowchart¶

Flowchart of processing steps. Each pRESTO tool is shown as a colored box. The workflow is divided into four primary tasks: (red) quality control, UMI annotation and primer masking; (orange) generation of UMI consensus sequences; (green) paired-end assembly of UMI consensus sequences; and (blue) deduplication and filtering to obtain the high-fidelity repertoire. Grey boxes indicate the initial and final data files. The intermediate files output by each tool are not shown for the sake of brevity.¶

Commands¶

#!/usr/bin/env bash
FilterSeq.py quality -s SRR1383456_1.fastq -q 20 --outname MS12_R1 --log FS1.log
FilterSeq.py quality -s SRR1383456_2.fastq -q 20 --outname MS12_R2 --log FS2.log
MaskPrimers.py score -s MS12_R1_quality-pass.fastq -p Stern2014_CPrimers.fasta \
    --start 15 --mode cut --barcode --outname MS12_R1 --log MP1.log
MaskPrimers.py score -s MS12_R2_quality-pass.fastq -p Stern2014_VPrimers.fasta \
    --start 0 --mode mask --outname MS12_R2 --log MP2.log
PairSeq.py -1 MS12_R1_primers-pass.fastq -2 MS12_R2_primers-pass.fastq \
    --1f BARCODE --coord sra
BuildConsensus.py -s MS12_R1_primers-pass_pair-pass.fastq --bf BARCODE --pf PRIMER \
    --prcons 0.6 --maxerror 0.1 --maxgap 0.5 --outname MS12_R1 --log BC1.log
BuildConsensus.py -s MS12_R2_primers-pass_pair-pass.fastq --bf BARCODE --pf PRIMER \
    --maxerror 0.1 --maxgap 0.5 --outname MS12_R2 --log BC2.log
PairSeq.py -1 MS12_R1_consensus-pass.fastq -2 MS12_R2_consensus-pass.fastq \
    --coord presto
AssemblePairs.py align -1 MS12_R2_consensus-pass_pair-pass.fastq \
    -2 MS12_R1_consensus-pass_pair-pass.fastq --coord presto --rc tail \
    --1f CONSCOUNT --2f CONSCOUNT PRCONS --outname MS12 --log AP.log
ParseHeaders.py collapse -s MS12_assemble-pass.fastq -f CONSCOUNT --act min
CollapseSeq.py -s MS12*reheader.fastq -n 20 --inner --uf PRCONS \
    --cf CONSCOUNT --act sum --outname MS12
SplitSeq.py group -s MS12_collapse-unique.fastq -f CONSCOUNT --num 2 --outname MS12
ParseHeaders.py table -s MS12_atleast-2.fastq -f ID PRCONS CONSCOUNT DUPCOUNT
ParseLog.py -l FS1.log FS2.log -f ID QUALITY
ParseLog.py -l MP1.log MP2.log -f ID PRIMER BARCODE ERROR
ParseLog.py -l BC1.log BC2.log -f BARCODE SEQCOUNT CONSCOUNT PRIMER PRCONS PRCOUNT \
    PRFREQ ERROR
ParseLog.py -l AP.log -f ID LENGTH OVERLAP ERROR PVALUE FIELDS1 FIELDS2

Vander Heiden et al, 2017 Example Files

Quality control, UMI annotation and primer masking¶

Removal of low quality reads¶

Quality control begins with the identification and removal of low-quality reads using the quality subcommand of the FilterSeq.py tool. In this example, reads with mean Phred quality scores less than 20 (-q 20) are removed:

FilterSeq.py quality -s SRR1383456_1.fastq -q 20 --outname MS12_R1 --log FS1.log
FilterSeq.py quality -s SRR1383456_2.fastq -q 20 --outname MS12_R2 --log FS2.log

The ParseLog.py tool is then used to extract results from the FilterSeq.py logs into tab-delimited files:

24	ParseLog.py -l FS1.log FS2.log -f ID QUALITY

Extracting the following information from the log:

Field	Description
ID	Sequence name
QUALITY	Quality score

UMI annotation and masking of primer regions¶

Next, the score subcommand of MaskPrimers.py is used to identify and remove the PCR primers for both reads. When dealing with Ig sequences, it is important to cut or mask the primers, as B cell receptors are subject to somatic hypermutation (the accumulation of point mutations in the DNA) and degenerate primer matches can look like mutations in downstream applications. The MaskPrimers tool is also used to annotate each read 1 sequence with the 15 nucleotide UMI that precedes the C-region primer (MaskPrimers score --barcode):

MaskPrimers.py score -s MS12_R1_quality-pass.fastq -p Stern2014_CPrimers.fasta \
    --start 15 --mode cut --barcode --outname MS12_R1 --log MP1.log
MaskPrimers.py score -s MS12_R2_quality-pass.fastq -p Stern2014_VPrimers.fasta \
    --start 0 --mode mask --outname MS12_R2 --log MP2.log

To summarize these steps, the ParseLog.py tool is used to build a tab-delimited file from the MaskPrimers.py log:

25	ParseLog.py -l MP1.log MP2.log -f ID PRIMER BARCODE ERROR

Containing the following information:

Field	Description
ID	Sequence name
PRIMER	Primer name
BARCODE	UMI sequence
ERROR	Primer match error rate

Note

For this data set the UMI is immediately upstream of the C-region primer. Another common approach for UMI barcoding involves placing the UMI immediately upstream of a 5’RACE template switch site. Modifying the workflow is simple for this case. You just need to replace the V-segment primers with a fasta file containing the TS sequences and move the --barcode argument to the appropriate read:

MaskPrimers.py score -s R1_quality-pass.fastq -p CPrimers.fasta \
    --start 0 --mode cut --outname R1 --log MP1.log
MaskPrimers.py score -s R2_quality-pass.fastq -p TSSites.fasta \
    --start 17 --barcode --mode cut --maxerror 0.5 \
    --outname R2 --log MP2.log

In the above we have moved the UMI annotation to read 2, increased the allowable error rate for matching the TS site (--maxerror 0.5), cut the TS site (--mode cut), and increased the size of the UMI from 15 to 17 nucleotides (--start 17).

Generation of UMI consensus sequences¶

Copying the UMI annotation across paired-end files¶

In this task, a single consensus sequence is constructed for each set of reads annotated with the same UMI barcode. As the UMI barcode is part of read 1, the BARCODE annotation identified by MaskPrimers.py must first be copied to the read 2 mate-pair of each read 1 sequence. Propogation of annotations between mate pairs is performed using PairSeq.py which also removes unpaired reads and ensures that paired reads are sorted in the same order across files:

PairSeq.py -1 MS12_R1_primers-pass.fastq -2 MS12_R2_primers-pass.fastq \
    --1f BARCODE --coord sra

Note

For both the PairSeq.py and AssemblePairs.py commands using the correct --coord argument is critical for matching mate-pairs. If this was raw data from Illumina, rather than data downloaded from SRA/ENA, then the appropriate argument would be --coord illumina.

Note

If you have followed the 5’RACE modification above, then you must also modify the first PairSeq.py step to copy the UMI from read 2 to read 1, instead of vice versa (--2f BARCODE):

PairSeq.py -1 R1_primers-pass.fastq -2 R2_primers-pass.fastq \
    --2f BARCODE --coord sra

Multiple alignment of UMI read groups¶

Before generating a consensus for a set of reads sharing a UMI barcode, the sequences must be properly aligned. Sequences may not be aligned if more than one PCR primer is identified in a UMI read group - leading to variations in the the start positions of the reads. Ideally, each set of reads originating from a single mRNA molecule should be amplified with the same primer. However, different primers in the multiplex pool may be incorporated into the same UMI read group during amplification if the primers are sufficiently similar. This type of primer misalignment can be corrected using the AlignSets.py tool. In the example data used here, this step was not necessary due to the aligned primer design for the 45 V-segment primers, though this does require that the V-segment primers be masked, rather than cut, during the MaskPrimers.py step (--mode mask).

See also

If your data requires alignment, then you can create multiple aligned UMI read groups as follows:

AlignSets.py muscle -s R1_primers-pass_pair-pass.fastq --bf BARCODE \
    --exec ~/bin/muscle --outname R1 --log AS1.log
AlignSets.py muscle -s R2_primers-pass_pair-pass.fastq --bf BARCODE \
    --exec ~/bin/muscle --outname R2 --log AS2.log

Where the --bf BARCODE defines the field containing the UMI and --exec ~/bin/muscle is the location of the MUSCLE executable.

For additional details see the section on fixing UMI alignments.

Generating UMI consensus reads¶

After alignment, a single consensus sequence is generated for each UMI barcode using BuildConsensus.py:

BuildConsensus.py -s MS12_R1_primers-pass_pair-pass.fastq --bf BARCODE --pf PRIMER \
    --prcons 0.6 --maxerror 0.1 --maxgap 0.5 --outname MS12_R1 --log BC1.log
BuildConsensus.py -s MS12_R2_primers-pass_pair-pass.fastq --bf BARCODE --pf PRIMER \
    --maxerror 0.1 --maxgap 0.5 --outname MS12_R2 --log BC2.log

To correct for UMI chemistry and sequencing errors, UMI read groups having high error statistics (mismatch rate from consensus) are removed by specifiying the --maxerror 0.1 threshold. As the accuracy of the primer assignment in read 1 is critical for correct isotype identification, additional filtering of read 1 is carried out during this step. Specifying the --prcons 0.6 threshold: (a) removes individual sequences that do not share a common primer annotation with the majority of the set, (b) removes entire read groups which have ambiguous primer assignments, and (c) constructs a consensus primer assignment for each UMI.

Note

The --maxgap 0.5 argument tells BuildConsensus.py to use a majority rule to delete any gap positions which occur in more than 50% of the reads. The --maxgap argument is not really necessary for this example data set as we did not perform a multiple alignment of the UMI read groups. However, if you have performed an alignment, then use of --maxgap during consensus generation is highly recommended.

The ParseLog.py tool is then used to build a tab-delimited file contain the consensus results:

ParseLog.py -l BC1.log BC2.log -f BARCODE SEQCOUNT CONSCOUNT PRIMER PRCONS PRCOUNT \

With the following annotations:

Field	Description
BARCODE	UMI sequence
SEQCOUNT	Number of total reads in the UMI group
CONSCOUNT	Number of reads used for the UMI consensus
PRIMER	Set of primer names in the UMI group
PRCONS	Consensus primer name
PRCOUNT	Count of primers in the UMI group
PRFREQ	Frequency of primers in the UMI group
ERROR	Average mismatch rate from consensus

Paired-end assembly of UMI consensus sequences¶

Syncronizing paired-end files¶

Following UMI consensus generation, the read 1 and read 2 files may again be out of sync due to differences in UMI read group filtering by BuildConsensus.py. To synchronize the reads another instance of PairSeq.py must be run, but without any annotation manipulation:

PairSeq.py -1 MS12_R1_consensus-pass.fastq -2 MS12_R2_consensus-pass.fastq \
    --coord presto

Assembling UMI consensus mate-pairs¶

Once the files have been synchronized, each paired-end UMI consensus sequence is assembled into a full length Ig sequence using the align subcommand of AssemblePairs.py:

AssemblePairs.py align -1 MS12_R2_consensus-pass_pair-pass.fastq \
    -2 MS12_R1_consensus-pass_pair-pass.fastq --coord presto --rc tail \
    --1f CONSCOUNT --2f CONSCOUNT PRCONS --outname MS12 --log AP.log

During assembly, the consensus isotype annotation (PRCONS) from read 1 and the number of reads used to define the consensus sequence (CONSCOUNT) for both reads are propagated into the annotations of the full length Ig sequence (--1f CONSCOUNT --2f CONSCOUNT PRCONS.

ParseLog.py is then uses to extract the results from the AssemblePairs.py log into a tab-delimited file:

27	PRFREQ ERROR

Containing the following information:

Field	Description
ID	Sequence name (UMI)
LENGTH	Length of the assembled sequence
OVERLAP	Length of the overlap between mate-pairs
ERROR	Mismatch rate of the overlapping region
PVALUE	P-value for the assembly
FIELDS1	Annotations copied from read 2 into the assembled sequence
FIELDS2	Annotations copied from read 1 into the assembled sequence

See also

Depending on the amplicon length in your data, not all mate-pairs may overlap. For the sake of simplicity, we have excluded a demonstration of assembly in such cases. pRESTO provides a couple approaches to deal with such reads. The reference subcommand of AssemblePairs.py can use the ungapped V-segment reference sequences to properly space non-overlapping reads. Or, if all else fails, the join subcommand can be used to simply stick mate-pairs together end-to-end with some intervening gap.

Deduplication and filtering¶

Combining UMI read group size annotations¶

In the final stage of the workflow, the high-fidelity Ig repertoire is obtained by a series of filtering steps. First, the annotation specifying the number of raw reads used to build each sequence (-f CONSCOUNT) is updated to be the minimum (--act min) of the forward and reverse reads using the collapse subcommand of ParseHeaders.py:

19	ParseHeaders.py collapse -s MS12_assemble-pass.fastq -f CONSCOUNT --act min

Removal of duplicate sequences¶

Second, duplicate nucleotide sequences are removed using the CollapseSeq.py tool with the requirement that duplicate sequences share the same isotype primer (--uf PRCONS). The duplicate removal step also removes sequences with a high number of interior N-valued nucleotides (-n 20 and --inner) and combines the read counts for each UMI read group (--cf CONSCOUNT and --act sum).

CollapseSeq.py -s MS12*reheader.fastq -n 20 --inner --uf PRCONS \
    --cf CONSCOUNT --act sum --outname MS12

Filtering to sequences with at least two representative reads¶

Finally, unique sequences are filtered to those with at least 2 contributing sequences using the group subcommand of SplitSeq.py, by splitting the file on the CONSCOUNT annotation with a numeric threshold (-f CONSCOUNT and --num 2):

SplitSeq.py group -s MS12_collapse-unique.fastq -f CONSCOUNT --num 2 --outname MS12

Creating an annotation table¶

For further analysis, the annotations of the final repertoire are then converted to into a table using the table subcommand of ParseHeaders.py:

ParseHeaders.py table -s MS12_atleast-2.fastq -f ID PRCONS CONSCOUNT DUPCOUNT

Output files¶

The final set of sequences, which serve as input to a V(D)J reference aligner (Eg, IMGT/HighV-QUEST or IgBLAST), and tables that can be plotted for quality control are:

File	Description
M12_collapse-unique.fastq	Total unique sequences
M12_atleast-2.fastq	Unique sequences represented by at least 2 reads
M12_atleast-2_headers.tab	Annotation table of the atleast-2 file
FS1_table.tab	Table of the read 1 FilterSeq log
FS2_table.tab	Table of the read 2 FilterSeq log
MP1_table.tab	Table of the C-region MaskPrimers log
MP2_table.tab	Table of the V-segment MaskPrimers log
BC1_table.tab	Table of the read 1 BuildConsensus log
BC2_table.tab	Table of the read 2 BuildConsensus log
AP_table.tab	Table of the AssemblePairs log

A number of other intermediate and log files are generated during the workflow, which allows easy tracking/reversion of processing steps. These files are not listed in the table above.

Performance¶

Example performance statistics for a comparable, but larger, MiSeq workflow are presented below. Performance was measured on a 64-core system with 2.3GHz AMD Opteron(TM) 6276 processors and 512GB of RAM, with memory usage measured at peak utilization. The data set contained 1,723,558 x 2 raw reads, and required matching of 1 C-region primer, 45 V-segment primers, and averaged 24.3 reads per UMI.

Line	Tool	Reads	Cores	MB	Minutes
01	R1: FilterSeq.py quality	1,723,558	10	1,219	13.0
02	R2: FilterSeq.py quality	1,723,558	10	1,219	12.8
03	R1: MaskPrimers.py score	1,722,116	10	1,221	18.9
04	R2: MaskPrimers.py score	1,684,050	10	1,221	46.7
06	PairSeq.py	1,665,584	1	4,734	25.5
07	R1: BuildConsensus.py	1,565,017	10	1,228	50.1
08	R2: BuildConsensus.py	1,565,017	10	1,229	58.6
11	AssemblePairs.py align	66,285	10	358	3.0
13	ParseHeaders.py collapse	56,104	1	88	0.4
14	CollapseSeq.py	55,480	1	822	0.7
15	SplitSeq.py group	51,047	1	88	0.2

UMI Barcoded Illumina MiSeq 325+275 paired-end 5’RACE BCR mRNA¶

Overview of Experimental Data¶

This example uses publicly available data from:

Dysregulation of B Cell Repertoire Formation in Myasthenia Gravis Patients Revealed through Deep Sequencing.

Vander Heiden JA, et al.

J Immunol. 2017 198(4):1460-1473. doi:10.4049/jimmunol.1601415.

which may be downloaded from the NCBI Sequence Read Archive under BioProject accession ID: PRJNA248475. For this example, we will use the first 25,000 sequences of sample HD09_N_AB8KB (accession: SRR4026043), which may downloaded using fastq-dump from the SRA Toolkit:

fastq-dump --split-files -X 25000 SRR4026043

Primer sequences are available online from the protocols/AbSeq directory in the Immcantation repository.

Read Configuration¶

Schematic of the Illumina MiSeq 325+275 paired-end reads with UMI barcodes. Each read was sequenced from one end of the target cDNA so that the two reads together cover the entire variable region of the Ig heavy chain. Sequencing was performed using an uneven number of cycles for read 1 and read 2 using a 2x300 kit. The V(D)J reading frame proceeds from the start of read 2 to the start of read 1. Read 1 is in the opposite orientation (reverse complement), contains a partial C-region, and is 325 nucleotides in length. Read 2 contains the 5’RACE template switch site with a 17 nucleotide UMI barcode preceding it, and is 275 nucleotides in length.¶

Example Data¶

We have hosted a small subset of the data (Accession: SRR1383456) on the pRESTO website in FASTQ format with accompanying primer files. The sample data set and workflow script may be downloaded from here:

Overview of the Workflow¶

In the following sections, we demonstrate each step of the workflow to move from raw sequence reads to a fully annotated repertoire of complete V(D)J sequences. The workflow is divided into four high-level tasks:

Quality control, UMI annotation and primer masking

Generation of UMI consensus sequences

Paired-end assembly of UMI consensus sequences

Deduplication and filtering

A graphical representation of the workflow along with the corresponding sequence of pRESTO commands is shown below.

Flowchart¶

Flowchart of processing steps. Each pRESTO tool is shown as a colored box. The workflow is divided into four primary tasks: (red) quality control, UMI annotation and primer masking; (orange) generation of UMI consensus sequences; (green) paired-end assembly of UMI consensus sequences; and (blue) deduplication and filtering to obtain the high-fidelity repertoire. Grey boxes indicate the initial and final data files. The intermediate files output by each tool are not shown for the sake of brevity.¶

Commands¶

#!/usr/bin/env bash
FilterSeq.py quality -s SRR4026043_1.fastq -q 20 --outname HD09N-R1 --log FS1.log
FilterSeq.py quality -s SRR4026043_2.fastq -q 20 --outname HD09N-R2 --log FS2.log
MaskPrimers.py score -s HD09N-R1_quality-pass.fastq -p AbSeq_R1_Human_IG_Primers.fasta \
    --start 0 --mode cut --outname HD09N-R1 --log MP1.log
MaskPrimers.py score -s HD09N-R2_quality-pass.fastq -p AbSeq_R2_TS.fasta \
    --start 17 --barcode --mode cut --maxerror 0.5 --outname HD09N-R2 --log MP2.log
PairSeq.py -1 HD09N-R1_primers-pass.fastq -2 HD09N-R2_primers-pass.fastq \
    --2f BARCODE --coord sra
BuildConsensus.py -s HD09N-R1_primers-pass_pair-pass.fastq --bf BARCODE --pf PRIMER \
    --prcons 0.6 --maxerror 0.1 --maxgap 0.5 --outname HD09N-R1 --log BC1.log
BuildConsensus.py -s HD09N-R2_primers-pass_pair-pass.fastq --bf BARCODE \
    --maxerror 0.1 --maxgap 0.5 --outname HD09N-R2 --log BC2.log
PairSeq.py -1 HD09N-R1_consensus-pass.fastq -2 HD09N-R2_consensus-pass.fastq \
    --coord presto
AssemblePairs.py sequential -1 HD09N-R2_consensus-pass_pair-pass.fastq \
    -2 HD09N-R1_consensus-pass_pair-pass.fastq -r IMGT_Human_IG_V.fasta \
    --coord presto --rc tail --scanrev --1f CONSCOUNT --2f CONSCOUNT PRCONS \
    --aligner blastn --outname HD09N-C --log AP.log
MaskPrimers.py align -s HD09N-C_assemble-pass.fastq \
    -p AbSeq_Human_IG_InternalCRegion.fasta --maxlen 100 --maxerror 0.3 \
    --mode tag --revpr --skiprc --pf CREGION --outname HD09N-C --log MP3.log 
ParseHeaders.py collapse -s HD09N-C_primers-pass.fastq -f CONSCOUNT --act min
CollapseSeq.py -s HD09N-C_primers-pass_reheader.fastq -n 20 --inner \
    --uf CREGION --cf CONSCOUNT --act sum --outname HD09N-C
SplitSeq.py group -s HD09N-C_collapse-unique.fastq \
    -f CONSCOUNT --num 2 --outname HD09N-C
ParseHeaders.py table -s HD09N-C_atleast-2.fastq -f ID CREGION CONSCOUNT DUPCOUNT
ParseLog.py -l FS1.log FS2.log -f ID QUALITY
ParseLog.py -l MP1.log MP2.log MP3.log -f ID PRIMER BARCODE ERROR
ParseLog.py -l BC1.log BC2.log -f BARCODE SEQCOUNT CONSCOUNT PRCONS PRFREQ ERROR
ParseLog.py -l AP.log -f ID REFID LENGTH OVERLAP GAP ERROR IDENTITY