Skip to content
Commits on Source (3)
Show whitespace changes
Inline Side-by-side
README.md
View file @ 366d9543
......@@ -5,6 +5,13 @@ Porechop is a tool for finding and removing adapters from [Oxford Nanopore](http
Porechop also supports demultiplexing of Nanopore reads that were barcoded with the [Native Barcoding Kit](https://store.nanoporetech.com/native-barcoding-kit-1d.html), [PCR Barcoding Kit](https://store.nanoporetech.com/pcr-barcoding-kit-96.html) or [Rapid Barcoding Kit](https://store.nanoporetech.com/rapid-barcoding-sequencing-kit.html).
### Oct 2018 update: Porechop is officially unsupported
While I'm happy Porechop has so many users, it has always been a bit klugey and a pain to maintain. I don't have the time to give it the attention it deserves, so I'm going to now officially declare Porechop as abandonware (though the unanswered [issues](https://github.com/rrwick/Porechop/issues) and [pull requests](https://github.com/rrwick/Porechop/pulls) reveal that it already has been for some time). I've added a [known issues](#known-issues) section to the README to outline what I think is wrong with Porechop and how a reimplementation should look. I may someday (no promises though :stuck_out_tongue:) try to rewrite it from a blank canvas to address its faults.
# Table of contents
......@@ -24,6 +31,7 @@ Porechop also supports demultiplexing of Nanopore reads that were barcoded with
* [Verbose output](#verbose-output)
* [Known adapters](#known-adapters)
* [Full usage](#full-usage)
* [Known issues](#known-issues)
* [Acknowledgements](#acknowledgements)
* [License](#license)
......@@ -91,9 +99,6 @@ __Demultiplex barcoded reads:__<br>
__Demultiplex barcoded reads, straight from Albacore output directory:__<br>
`porechop -i albacore_dir -b output_dir`
__Demultiplex barcoded reads without trimming (appropriate if Nanopolish will be used):__<br>
`porechop -i input_reads.fastq.gz -b output_dir --untrimmed`
__Also works with FASTA:__<br>
`porechop -i input_reads.fasta -o output_reads.fasta`
......@@ -143,11 +148,9 @@ TGTTGTTGTTGTTATTGTTGTTATTGTTGTTGTATTGTTGTTATTGTTGTTGTTGTACATTGTTATTGTTGTATTGTTGT
If you run Porechop with `--discard_middle`, the reads with internal adapters will be thrown out instead of split.
This approach might make sense if you are trimming reads from a barcoded run, as chimeric reads may combine sequences from different bins. For example, consider this read:
```
BC01_rev - SEQUENCE_1 - SQK-NSK007_Y_Top - BC02_rev - SEQUENCE_2
```
SEQUENCE_1 belongs in the barcode 1 bin and SEQUENCE_2 belongs in the barcode 2 bin, so while we could split the read, we would end up with contamination from another bin. Throwing the read out with `--discard_middle` might be the better option.
If you plan on using your reads with [Nanopolish](https://github.com/jts/nanopolish), then the `--discard_middle` option is required. This is because Nanopolish first runs [`nanopolish index`](https://github.com/jts/nanopolish#data-preprocessing) to find a one-to-one association between FASTQ reads and fast5 files. If you ran Porechop _without_ `--discard_middle`, then you could end up with multiple separate FASTQ reads which are from a single fast5, and this is incompatible with Nanopolish.
This option is also recommended if you are trimming reads from a demultiplexed barcoded sequencing run. This is because chimeric reads may contain two sequences from two separate barcodes, so throwing them out is the safer option to reduce cross-barcode contamination.
### Barcode demultiplexing
......@@ -189,10 +192,10 @@ Alternately, you can run Porechop with `-b` which specifies a directory for barc
If Porechop is run without `-o` or `-b`, then it will output the trimmed reads to stdout and print its progress info to stderr. The output format of the reads will be FASTA/FASTQ based on the input reads, or else can be specified using `--format`.
The `--verbosity` option will change the output of progress info:
The `--verbosity` option will change the amount of progress info:
* `--verbosity 0` gives no progress output.
* `--verbosity 1` (the default) gives summary info about end adapter trimming and shows all instances of middle adapter splitting.
* `--verbosity 2` shows sequences and is described below.
* `--verbosity 1` (default) gives summary info about end adapter trimming and middle adapter splitting.
* `--verbosity 2` shows the actual trimmed/split sequences for each read (described more below).
* `--verbosity 3` shows tons of data (mainly for debugging).
......@@ -233,37 +236,35 @@ If you want to add your own adapter sequences to Porechop, you can do so by edit
# Full usage
```
usage: porechop [-h] -i INPUT [-o OUTPUT] [--format {auto,fasta,fastq,fasta.gz,fastq.gz}] [-v VERBOSITY]
[-t THREADS] [--version] [-b BARCODE_DIR] [--barcode_threshold BARCODE_THRESHOLD]
usage: porechop -i INPUT [-o OUTPUT] [--format {auto,fasta,fastq,fasta.gz,fastq.gz}] [-v VERBOSITY]
[-t THREADS] [-b BARCODE_DIR] [--barcode_threshold BARCODE_THRESHOLD]
[--barcode_diff BARCODE_DIFF] [--require_two_barcodes] [--untrimmed]
[--adapter_threshold ADAPTER_THRESHOLD] [--check_reads CHECK_READS]
[--scoring_scheme SCORING_SCHEME] [--end_size END_SIZE] [--min_trim_size MIN_TRIM_SIZE]
[--extra_end_trim EXTRA_END_TRIM] [--end_threshold END_THRESHOLD] [--discard_middle]
[--discard_unassigned] [--adapter_threshold ADAPTER_THRESHOLD]
[--check_reads CHECK_READS] [--scoring_scheme SCORING_SCHEME] [--end_size END_SIZE]
[--min_trim_size MIN_TRIM_SIZE] [--extra_end_trim EXTRA_END_TRIM]
[--end_threshold END_THRESHOLD] [--no_split] [--discard_middle]
[--middle_threshold MIDDLE_THRESHOLD]
[--extra_middle_trim_good_side EXTRA_MIDDLE_TRIM_GOOD_SIDE]
[--extra_middle_trim_bad_side EXTRA_MIDDLE_TRIM_BAD_SIDE]
[--min_split_read_size MIN_SPLIT_READ_SIZE]
Porechop: a tool for finding adapters in Oxford Nanopore reads, trimming them from the ends and splitting
reads with internal adapters
[--min_split_read_size MIN_SPLIT_READ_SIZE] [-h] [--version]
optional arguments:
-h, --help show this help message and exit
Porechop: a tool for finding adapters in Oxford Nanopore reads, trimming them from the ends and
splitting reads with internal adapters
Main options:
-i INPUT, --input INPUT FASTA/FASTQ of input reads or a directory which will be recursively
searched for FASTQ files (required)
-i INPUT, --input INPUT FASTA/FASTQ of input reads or a directory which will be
recursively searched for FASTQ files (required)
-o OUTPUT, --output OUTPUT Filename for FASTA or FASTQ of trimmed reads (if not set, trimmed
reads will be printed to stdout)
--format {auto,fasta,fastq,fasta.gz,fastq.gz}
Output format for the reads - if auto, the format will be chosen based
on the output filename or the input read format (default: auto)
Output format for the reads - if auto, the format will be chosen
based on the output filename or the input read format (default:
auto)
-v VERBOSITY, --verbosity VERBOSITY
Level of progress information: 0 = none, 1 = some, 2 = lots, 3 = full
- output will go to stdout if reads are saved to a file and stderr if
reads are printed to stdout (default: 1)
Level of progress information: 0 = none, 1 = some, 2 = lots, 3 =
full - output will go to stdout if reads are saved to a file and
stderr if reads are printed to stdout (default: 1)
-t THREADS, --threads THREADS Number of threads to use for adapter alignment (default: 8)
--version show program's version number and exit
Barcode binning settings:
Control the binning of reads based on barcodes (i.e. barcode demultiplexing)
......@@ -275,14 +276,15 @@ Barcode binning settings:
A read must have at least this percent identity to a barcode to be
binned (default: 75.0)
--barcode_diff BARCODE_DIFF If the difference between a read's best barcode identity and its
second-best barcode identity is less than this value, it will not be
put in a barcode bin (to exclude cases which are too close to call)
(default: 5.0)
--require_two_barcodes Reads will only be put in barcode bins if they have a strong match for
the barcode on both their start and end (default: a read can be binned
with a match at its start or end)
--untrimmed Bin reads but do not trim the ends (appropriate if reads are to be
used with Nanopolish) (default: False)
second-best barcode identity is less than this value, it will not
be put in a barcode bin (to exclude cases which are too close to
call) (default: 5.0)
--require_two_barcodes Reads will only be put in barcode bins if they have a strong match
for the barcode on both their start and end (default: a read can
be binned with a match at its start or end)
--untrimmed Bin reads but do not trim them (default: trim the reads)
--discard_unassigned Discard unassigned reads (instead of creating a "none" bin)
(default: False)
Adapter search settings:
Control how the program determines which adapter sets are present
......@@ -290,10 +292,11 @@ Adapter search settings:
--adapter_threshold ADAPTER_THRESHOLD
An adapter set has to have at least this percent identity to be
labelled as present and trimmed off (0 to 100) (default: 90.0)
--check_reads CHECK_READS This many reads will be aligned to all possible adapters to determine
which adapter sets are present (default: 10000)
--scoring_scheme SCORING_SCHEME Comma-delimited string of alignment scores: match, mismatch, gap open,
gap extend (default: 3,-6,-5,-2)
--check_reads CHECK_READS This many reads will be aligned to all possible adapters to
determine which adapter sets are present (default: 10000)
--scoring_scheme SCORING_SCHEME
Comma-delimited string of alignment scores: match, mismatch, gap
open, gap extend (default: 3,-6,-5,-2)
End adapter settings:
Control the trimming of adapters from read ends
......@@ -301,35 +304,75 @@ End adapter settings:
--end_size END_SIZE The number of base pairs at each end of the read which will be
searched for adapter sequences (default: 150)
--min_trim_size MIN_TRIM_SIZE Adapter alignments smaller than this will be ignored (default: 4)
--extra_end_trim EXTRA_END_TRIM This many additional bases will be removed next to adapters found at
the ends of reads (default: 2)
--end_threshold END_THRESHOLD Adapters at the ends of reads must have at least this percent identity
to be removed (0 to 100) (default: 75.0)
--extra_end_trim EXTRA_END_TRIM
This many additional bases will be removed next to adapters found
at the ends of reads (default: 2)
--end_threshold END_THRESHOLD Adapters at the ends of reads must have at least this percent
identity to be removed (0 to 100) (default: 75.0)
Middle adapter settings:
Control the splitting of read from middle adapters
--no_split Skip splitting reads based on middle adapters (default: split reads
when an adapter is found in the middle)
--no_split Skip splitting reads based on middle adapters (default: split
reads when an adapter is found in the middle)
--discard_middle Reads with middle adapters will be discarded (default: reads with
middle adapters are split) (this option is on by default when
outputting reads into barcode bins)
middle adapters are split) (required for reads to be used with
Nanopolish, this option is on by default when outputting reads
into barcode bins)
--middle_threshold MIDDLE_THRESHOLD
Adapters in the middle of reads must have at least this percent
identity to be found (0 to 100) (default: 85.0)
--extra_middle_trim_good_side EXTRA_MIDDLE_TRIM_GOOD_SIDE
This many additional bases will be removed next to middle adapters on
their "good" side (default: 10)
This many additional bases will be removed next to middle adapters
on their "good" side (default: 10)
--extra_middle_trim_bad_side EXTRA_MIDDLE_TRIM_BAD_SIDE
This many additional bases will be removed next to middle adapters on
their "bad" side (default: 100)
This many additional bases will be removed next to middle adapters
on their "bad" side (default: 100)
--min_split_read_size MIN_SPLIT_READ_SIZE
Post-split read pieces smaller than this many base pairs will not be
outputted (default: 1000)
Post-split read pieces smaller than this many base pairs will not
be outputted (default: 1000)
Help:
-h, --help Show this help message and exit
--version Show program's version number and exit
```
# Known issues
### Adapter search
Porechop tries to automatically determine which adapters are present by looking at the reads, but this approach has a few issues:
* As the number of kits/barcodes has grown, adapter-search part of the Porechop's pipeline has become increasingly slow.
* Porechop only does the adapter search on a subset of reads, which means there can be problems with non-randomly ordered read sets (e.g. all barcode 1 reads at the start of a file, followed by barcode 2 reads, etc).
* Many ONT adapters share common sequence with each other, making false positive adapter finds possible.
A simpler solution (and in hindsight what I should have done) would be to require the kit and/or adapters from the user. E.g. `porechop --sqk-lsk109` or `porechop --start_adapt ACGCTAGCATACGT`.
### Performance
Porechop uses [SeqAn](https://github.com/seqan/seqan) to perform its alignments in C++. This library is very flexible, but not as fast as some alternatives, such as [Edlib](https://github.com/Martinsos/edlib).
Another performance issue is that Porechop uses [ctypes](https://docs.python.org/3/library/ctypes.html) to interface with its C++ code. Function calls with ctypes can have a bit of overhead, which means that Porechop cannot use threads very efficiently (it spends too much of its time in the Python code, which is intrinsically non-parallel).
### Barcode demultiplexing
I added demultiplexing to Porechop as an afterthought – it was already looking for barcodes to trim, so why not sort reads by barcodes too? This turned out to be a very useful feature, but in hindsight I think it might be simpler (and easier to maintain) if trimming and demultiplexing functionality were in separate tools.
### Base-space problems
I've encountered a couple of issues where adapter sequences are not properly basecalled, resulting in inconsistent sequence. Porechop trims in base-space, so this is a somewhat intractable problem. See [issue #40](https://github.com/rrwick/Porechop/issues/40) for an example.
These have made me wonder if the adapter trimming should be done in signal-space instead, though that would be a more complex problem to solve. I hope that in the future ONT can integrate this kind of functionality directly into their basecallers.
# Acknowledgements
Porechop was inspired by (and largely coded during) [Porecamp Australia 2017](https://porecamp-au.github.io/). Thanks to the organisers and attendees who helped me realise that a Nanopore adapter trimmer might be a useful tool! I later met David Stoddart from Oxford Nanopore at London Calling 2017, and he helped me get many of the adapter sequences right.
......
debian/changelog
View file @ 366d9543
porechop (0.2.4+dfsg-1) unstable; urgency=medium
* Team upload.
* New upstream version.
-- Steffen Moeller <moeller@debian.org> Tue, 23 Oct 2018 13:53:48 +0200
porechop (0.2.3+dfsg-3) unstable; urgency=medium
* Team upload.
......
porechop/adapters.py
View file @ 366d9543
......@@ -43,7 +43,9 @@ class Adapter(object):
Gets the barcode name for the output files. We want a concise name, so it looks at all
options and chooses the shortest.
"""
possible_names = [self.name, self.start_sequence[0]]
possible_names = [self.name]
if self.start_sequence:
possible_names.append(self.start_sequence[0])
if self.end_sequence:
possible_names.append(self.end_sequence[0])
barcode_name = sorted(possible_names, key=lambda x: len(x))[0]
......@@ -76,28 +78,36 @@ ADAPTERS = [Adapter('SQK-NSK007',
start_sequence=('SQK-NSK007_Y_Top', 'AATGTACTTCGTTCAGTTACGTATTGCT'),
end_sequence=('SQK-NSK007_Y_Bottom', 'GCAATACGTAACTGAACGAAGT')),
Adapter('Rapid',
start_sequence=('Rapid_adapter', 'GTTTTCGCATTTATCGTGAAACGCTTTCGCGTTTTTCGTGCGCCGCTTCA')),
start_sequence=('Rapid_adapter',
'GTTTTCGCATTTATCGTGAAACGCTTTCGCGTTTTTCGTGCGCCGCTTCA')),
Adapter('RBK004_upstream',
start_sequence=('RBK004_upstream', 'AATGTACTTCGTTCAGTTACGGCTTGGGTGTTTAACC')),
Adapter('SQK-MAP006',
start_sequence=('SQK-MAP006_Y_Top_SK63', 'GGTTGTTTCTGTTGGTGCTGATATTGCT'),
end_sequence= ('SQK-MAP006_Y_Bottom_SK64', 'GCAATATCAGCACCAACAGAAA')),
Adapter('SQK-MAP006 Short',
Adapter('SQK-MAP006 short',
start_sequence=('SQK-MAP006_Short_Y_Top_LI32', 'CGGCGTCTGCTTGGGTGTTTAACCT'),
end_sequence= ('SQK-MAP006_Short_Y_Bottom_LI33', 'GGTTAAACACCCAAGCAGACGCCG')),
# The PCR adapters are used both in PCR DNA kits and some cDNA kits.
Adapter('PCR adapters 1',
start_sequence=('PCR_1_start', 'ACTTGCCTGTCGCTCTATCTTC'),
end_sequence= ('PCR_1_end', 'GAAGATAGAGCGACAGGCAAGT')),
Adapter('PCR tail 1',
start_sequence=('PCR_tail_1_start', 'TTAACCTTTCTGTTGGTGCTGATATTGC'),
end_sequence= ('PCR_tail_1_end', 'GCAATATCAGCACCAACAGAAAGGTTAA')),
Adapter('PCR adapters 2',
start_sequence=('PCR_2_start', 'TTTCTGTTGGTGCTGATATTGC'),
end_sequence= ('PCR_2_end', 'GCAATATCAGCACCAACAGAAA')),
Adapter('PCR tail 2',
start_sequence=('PCR_tail_2_start', 'TTAACCTACTTGCCTGTCGCTCTATCTTC'),
end_sequence= ('PCR_tail_2_end', 'GAAGATAGAGCGACAGGCAAGTAGGTTAA')),
Adapter('PCR adapters 3',
start_sequence=('PCR_3_start', 'TACTTGCCTGTCGCTCTATCTTC'),
end_sequence= ('PCR_3_end', 'GAAGATAGAGCGACAGGCAAGTA')),
# 1D^2 kit adapters are interesting. ONT provided the following sequences on their site:
......@@ -117,6 +127,11 @@ ADAPTERS = [Adapter('SQK-NSK007',
# adapter sequences here.
Adapter('cDNA SSP',
start_sequence=('cDNA_SSP', 'TTTCTGTTGGTGCTGATATTGCTGCCATTACGGCCGGG'),
end_sequence= ('cDNA_SSP_rev', 'CCCGGCCGTAATGGCAGCAATATCAGCACCAACAGAAA')),
# Some barcoding kits (like the native barcodes) use the rev comp barcode at the start
# of the read and the forward barcode at the end of the read.
Adapter('Barcode 1 (reverse)',
......@@ -461,11 +476,23 @@ def make_full_native_barcode_adapter(barcode_num):
end_sequence=('NB' + '%02d' % barcode_num + '_end', end_full_seq))
def make_full_rapid_barcode_adapter(barcode_num):
def make_old_full_rapid_barcode_adapter(barcode_num): # applies to SQK-RBK001
barcode = [x for x in ADAPTERS if x.name == 'Barcode ' + str(barcode_num) + ' (forward)'][0]
start_barcode_seq = barcode.start_sequence[1]
start_full_seq = 'AATGTACTTCGTTCAGTTACG' + 'TATTGCT' + start_barcode_seq + \
'GTTTTCGCATTTATCGTGAAACGCTTTCGCGTTTTTCGTGCGCCGCTTCA'
return Adapter('Rapid barcoding ' + str(barcode_num) + ' (full sequence, old)',
start_sequence=('RB' + '%02d' % barcode_num + '_full', start_full_seq))
def make_new_full_rapid_barcode_adapter(barcode_num): # applies to SQK-RBK004
barcode = [x for x in ADAPTERS if x.name == 'Barcode ' + str(barcode_num) + ' (forward)'][0]
start_barcode_seq = barcode.start_sequence[1]
start_full_seq = 'AATGTACTTCGTTCAGTTACGTATTGCT' + start_barcode_seq + 'GTTTTCGCATTTATCGTGAAACGCTTTCGCGTTTTTCGTGCGCCGCTTCA'
start_full_seq = 'AATGTACTTCGTTCAGTTACG' + 'GCTTGGGTGTTTAACC' + start_barcode_seq + \
'GTTTTCGCATTTATCGTGAAACGCTTTCGCGTTTTTCGTGCGCCGCTTCA'
return Adapter('Rapid barcoding ' + str(barcode_num) + ' (full sequence)',
return Adapter('Rapid barcoding ' + str(barcode_num) + ' (full sequence, new)',
start_sequence=('RB' + '%02d' % barcode_num + '_full', start_full_seq))
porechop/nanopore_read.py
View file @ 366d9543
......@@ -23,7 +23,13 @@ class NanoporeRead(object):
def __init__(self, name, seq, quals):
self.name = name
self.seq = seq
self.seq = seq.upper()
if self.seq.count('U') > self.seq.count('T'):
self.rna = True
self.seq = self.seq.replace('U', 'T')
else:
self.rna = False
self.quals = quals
if len(quals) < len(seq):
self.quals += '+' * (len(seq) - len(quals))
......@@ -96,6 +102,8 @@ class NanoporeRead(object):
seq = self.get_seq_with_start_end_adapters_trimmed()
if not seq: # Don't return empty sequences
return ''
if self.rna:
seq = seq.replace('T', 'U')
return ''.join(['>', self.name, '\n', add_line_breaks_to_sequence(seq, 70)])
elif discard_middle:
return ''
......@@ -106,6 +114,8 @@ class NanoporeRead(object):
if not split_read_part[0]: # Don't return empty sequences
return ''
seq = add_line_breaks_to_sequence(split_read_part[0], 70)
if self.rna:
seq = seq.replace('T', 'U')
fasta_str += ''.join(['>', read_name, '\n', seq])
return fasta_str
......@@ -119,6 +129,8 @@ class NanoporeRead(object):
quals = self.get_quals_with_start_end_adapters_trimmed()
if not seq: # Don't return empty sequences
return ''
if self.rna:
seq = seq.replace('T', 'U')
return ''.join(['@', self.name, '\n', seq, '\n+\n', quals, '\n'])
elif discard_middle:
return ''
......@@ -126,10 +138,12 @@ class NanoporeRead(object):
fastq_str = ''
for i, split_read_part in enumerate(self.get_split_read_parts(min_split_read_size)):
read_name = add_number_to_read_name(self.name, i + 1)
if not split_read_part[0]: # Don't return empty sequences
seq, qual = split_read_part[0], split_read_part[1]
if not seq: # Don't return empty sequences
return ''
fastq_str += ''.join(['@', read_name, '\n', split_read_part[0], '\n+\n',
split_read_part[1], '\n'])
if self.rna:
seq = seq.replace('T', 'U')
fastq_str += ''.join(['@', read_name, '\n', seq, '\n+\n', qual, '\n'])
return fastq_str
def align_adapter_set(self, adapter_set, end_size, scoring_scheme_vals):
......@@ -138,6 +152,7 @@ class NanoporeRead(object):
This is not to determine where to trim the reads, but rather to figure out which adapter
sets are present in the data.
"""
if adapter_set.start_sequence:
read_seq_start = self.seq[:end_size]
score, _, _, _ = align_adapter(read_seq_start, adapter_set.start_sequence[1],
scoring_scheme_vals)
......@@ -156,6 +171,8 @@ class NanoporeRead(object):
"""
read_seq_start = self.seq[:end_size]
for adapter in adapters:
if not adapter.start_sequence:
continue
full_score, partial_score, read_start, read_end = \
align_adapter(read_seq_start, adapter.start_sequence[1], scoring_scheme_vals)
if partial_score > end_threshold and read_end != end_size and \
......
porechop/porechop.py
View file @ 366d9543
......@@ -24,7 +24,8 @@ import re
from multiprocessing.dummy import Pool as ThreadPool
from collections import defaultdict
from .misc import load_fasta_or_fastq, print_table, red, bold_underline, MyHelpFormatter, int_to_str
from .adapters import ADAPTERS, make_full_native_barcode_adapter, make_full_rapid_barcode_adapter
from .adapters import ADAPTERS, make_full_native_barcode_adapter,\
make_old_full_rapid_barcode_adapter, make_new_full_rapid_barcode_adapter
from .nanopore_read import NanoporeRead
from .version import __version__
......@@ -37,16 +38,17 @@ def main():
matching_sets = find_matching_adapter_sets(check_reads, args.verbosity, args.end_size,
args.scoring_scheme_vals, args.print_dest,
args.adapter_threshold, args.threads)
matching_sets = exclude_end_adapters_for_rapid(matching_sets)
matching_sets = fix_up_1d2_sets(matching_sets)
display_adapter_set_results(matching_sets, args.verbosity, args.print_dest)
matching_sets = add_full_barcode_adapter_sets(matching_sets)
if args.barcode_dir:
forward_or_reverse_barcodes = choose_barcoding_kit(matching_sets, args.verbosity,
args.print_dest)
else:
forward_or_reverse_barcodes = None
display_adapter_set_results(matching_sets, args.verbosity, args.print_dest)
matching_sets = add_full_barcode_adapter_sets(matching_sets)
if args.verbosity > 0:
print('\n', file=args.print_dest)
......@@ -86,7 +88,7 @@ def get_arguments():
parser = argparse.ArgumentParser(description='Porechop: a tool for finding adapters in Oxford '
'Nanopore reads, trimming them from the ends and '
'splitting reads with internal adapters',
formatter_class=MyHelpFormatter)
formatter_class=MyHelpFormatter, add_help=False)
main_group = parser.add_argument_group('Main options')
main_group.add_argument('-i', '--input', required=True,
help='FASTA/FASTQ of input reads or a directory which will be '
......@@ -105,7 +107,6 @@ def get_arguments():
'a file and stderr if reads are printed to stdout')
main_group.add_argument('-t', '--threads', type=int, default=default_threads,
help='Number of threads to use for adapter alignment')
main_group.add_argument('--version', action='version', version=__version__)
barcode_group = parser.add_argument_group('Barcode binning settings',
'Control the binning of reads based on barcodes '
......@@ -127,8 +128,7 @@ def get_arguments():
'match for the barcode on both their start and end (default: '
'a read can be binned with a match at its start or end)')
barcode_group.add_argument('--untrimmed', action='store_true',
help='Bin reads but do not trim them (appropriate if reads are to '
'be used with Nanopolish) (default: trim the reads)')
help='Bin reads but do not trim them (default: trim the reads)')
barcode_group.add_argument('--discard_unassigned', action='store_true',
help='Discard unassigned reads (instead of creating a "none" bin)')
......@@ -169,9 +169,10 @@ def get_arguments():
'middle)')
middle_trim_group.add_argument('--discard_middle', action='store_true',
help='Reads with middle adapters will be discarded (default: '
'reads with middle adapters are split) (this option is '
'on by default when outputting reads into barcode bins)')
middle_trim_group.add_argument('--middle_threshold', type=float, default=85.0,
'reads with middle adapters are split) (required for '
'reads to be used with Nanopolish, this option is on by '
'default when outputting reads into barcode bins)')
middle_trim_group.add_argument('--middle_threshold', type=float, default=90.0,
help='Adapters in the middle of reads must have at least this '
'percent identity to be found (0 to 100)')
middle_trim_group.add_argument('--extra_middle_trim_good_side', type=int, default=10,
......@@ -184,6 +185,12 @@ def get_arguments():
help='Post-split read pieces smaller than this many base pairs '
'will not be outputted')
help_args = parser.add_argument_group('Help')
help_args.add_argument('-h', '--help', action='help', default=argparse.SUPPRESS,
help='Show this help message and exit')
help_args.add_argument('--version', action='version', version=__version__,
help="Show program's version number and exit")
args = parser.parse_args()
try:
......@@ -325,35 +332,43 @@ def choose_barcoding_kit(adapter_sets, verbosity, print_dest):
If the user is sorting reads by barcode bin, choose one barcode configuration (rev comp
barcodes at the start of the read or at the end of the read) and ignore the other.
"""
forward_barcodes = 0
reverse_barcodes = 0
# Tally up scores for forward and reverse barcodes.
forward_start_or_end, reverse_start_or_end = 0, 0
forward_start_and_end, reverse_start_and_end = 0, 0
for adapter_set in adapter_sets:
score = adapter_set.best_start_or_end_score()
if 'Barcode' in adapter_set.name and '(forward)' in adapter_set.name:
forward_barcodes += score
elif 'Barcode' in adapter_set.name and '(reverse)' in adapter_set.name:
reverse_barcodes += score
if forward_barcodes > reverse_barcodes:
if verbosity > 0:
print('\nBarcodes determined to be in forward orientation', file=print_dest)
return 'forward'
elif reverse_barcodes > forward_barcodes:
if verbosity > 0:
print('\nBarcodes determined to be in reverse orientation', file=print_dest)
return 'reverse'
else:
return None
if 'barcode' in adapter_set.name.lower():
if '(forward)' in adapter_set.name.lower():
forward_start_or_end += adapter_set.best_start_or_end_score()
forward_start_and_end += adapter_set.best_start_score
forward_start_and_end += adapter_set.best_end_score
elif '(reverse)' in adapter_set.name.lower():
reverse_start_or_end += adapter_set.best_start_or_end_score()
reverse_start_and_end += adapter_set.best_start_score
reverse_start_and_end += adapter_set.best_end_score
if forward_start_or_end == 0 and reverse_start_or_end == 0:
sys.exit('Error: no barcodes were found, so Porechop cannot perform barcode demultiplexing')
# If possible, make a decision using each barcode's best start OR end score.
orientation = None
if forward_start_or_end > reverse_start_or_end:
orientation = 'forward'
elif reverse_start_or_end > forward_start_or_end:
orientation = 'reverse'
# If that didn't work (i.e. it's a tie between forward and reverse), then choose based on the
# sum of both start AND end scores.
elif forward_start_and_end > reverse_start_and_end:
orientation = 'forward'
elif reverse_start_and_end > forward_start_and_end:
orientation = 'reverse'
if orientation is None:
sys.exit('Error: Porechop could not determine barcode orientation')
def exclude_end_adapters_for_rapid(matching_sets):
"""
Rapid reads shouldn't have end adapters, so we don't want to look for them if this seems to be
a rapid read set.
"""
if 'Rapid adapter' in [x.name for x in matching_sets]:
for s in matching_sets:
s.end_sequence = None
return matching_sets
if verbosity > 0:
print('\nBarcodes determined to be in ' + orientation + ' orientation', file=print_dest)
return orientation
def fix_up_1d2_sets(matching_sets):
......@@ -411,8 +426,11 @@ def add_full_barcode_adapter_sets(matching_sets):
# Rapid barcode full sequences
if all(x in matching_set_names
for x in ['SQK-NSK007', 'Rapid', 'Barcode ' + str(i) + ' (forward)']):
matching_sets.append(make_full_rapid_barcode_adapter(i))
for x in ['Rapid', 'Barcode ' + str(i) + ' (forward)']):
if 'RBK004_upstream' in matching_set_names:
matching_sets.append(make_new_full_rapid_barcode_adapter(i))
elif 'SQK-NSK007' in matching_set_names:
matching_sets.append(make_old_full_rapid_barcode_adapter(i))
return matching_sets
......@@ -424,9 +442,12 @@ def find_adapters_at_read_ends(reads, matching_sets, verbosity, end_size, extra_
if verbosity > 0:
print(bold_underline('Trimming adapters from read ends'),
file=print_dest)
name_len = max(max(len(x.start_sequence[0]) for x in matching_sets),
max(len(x.end_sequence[0]) if x.end_sequence else 0 for x in matching_sets))
name_len = max(max(len(x.start_sequence[0])
if x.start_sequence else 0 for x in matching_sets),
max(len(x.end_sequence[0])
if x.end_sequence else 0 for x in matching_sets))
for matching_set in matching_sets:
if matching_set.start_sequence:
print(' ' + matching_set.start_sequence[0].rjust(name_len) + ': ' +
red(matching_set.start_sequence[1]), file=print_dest)
if matching_set.end_sequence:
......@@ -519,14 +540,17 @@ def find_adapters_in_read_middles(reads, matching_sets, verbosity, middle_thresh
adapters = []
for matching_set in matching_sets:
if matching_set.start_sequence:
adapters.append(matching_set.start_sequence)
if matching_set.end_sequence and \
if matching_set.end_sequence:
if (not matching_set.start_sequence) or \
matching_set.end_sequence[1] != matching_set.start_sequence[1]:
adapters.append(matching_set.end_sequence)
start_sequence_names = set()
end_sequence_names = set()
for matching_set in matching_sets:
if matching_set.start_sequence:
start_sequence_names.add(matching_set.start_sequence[0])
if matching_set.end_sequence:
end_sequence_names.add(matching_set.end_sequence[0])
......
porechop/version.py
View file @ 366d9543
......@@ -2,4 +2,4 @@
The version is stored here in a separate file so it can exist in only one place.
http://stackoverflow.com/questions/458550
"""
__version__ = '0.2.3'
__version__ = '0.2.4'