README.md

@@ -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

+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

@@ -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

@@ -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

@@ -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

@@ -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'