CMakeLists.txt

@@ -2,7 +2,7 @@ cmake_minimum_required( VERSION 2.8.2 )
 
 project( FLEXBAR )
 
-set( SEQAN_APP_VERSION "3.4.0" )
+set( SEQAN_APP_VERSION "3.5.0" )
 
 include_directories( ${FLEXBAR_SOURCE_DIR}/include )
 # link_directories( ${FLEXBAR_SOURCE_DIR}/lib )

README.md

 ## Flexbar – flexible barcode and adapter removal
 
-The program Flexbar preprocesses high-throughput sequencing data efficiently. It demultiplexes barcoded runs and removes adapter sequences. Moreover, trimming and filtering features are provided. Flexbar increases read mapping rates and improves genome as well as transcriptome assemblies. Unique molecular identifiers can be extracted in a flexible way. The program supports sequencing data in fasta and fastq format, e.g. from the Illumina platform.
+The program Flexbar preprocesses high-throughput sequencing data efficiently. It demultiplexes barcoded runs and removes adapter sequences. Several adapter removal presets for Illumina libraries are included. Flexbar computes exact overlap alignments using SIMD and multicore parallelism. Moreover, trimming and filtering features are provided, e.g. trimming of homopolymers at read ends. Flexbar increases read mapping rates and improves genome as well as transcriptome assemblies. Unique molecular identifiers can be extracted in a flexible way. The software supports data in fasta and fastq format from multiple sequencing platforms.
 
 Refer to the [manual](https://github.com/seqan/flexbar/wiki) or contact [Johannes Roehr](https://github.com/jtroehr) for support with this application.
 
@@ -22,7 +22,15 @@ See article on [PubMed](https://www.ncbi.nlm.nih.gov/pubmed/24832523)
 
 ### Download
 
-Flexbar source code as well as binaries for Linux and Mac OS can be downloaded on the [release](https://github.com/seqan/flexbar/releases) page. Please follow instructions for building or setup of binaries below. Additionally, Flexbar is available via package manager on Debian systems and in Conda. Versions before 2.4 can be found on the [old](https://sourceforge.net/projects/flexbar) page.
+Flexbar source code as well as binaries for Linux and Mac OS can be downloaded on the [release](https://github.com/seqan/flexbar/releases) page. Please follow instructions for building or setup of binaries below. Additionally, Flexbar is available via package manager on Debian systems, in Homebrew, and in Bioconda. Versions before 2.4 can be found on the [old](https://sourceforge.net/projects/flexbar) page.
+
+Installation with package managers:
+
+* Debian: `sudo apt install flexbar`
+* Homebrew: `brew install brewsci/science/flexbar`
+* Bioconda: `conda install -c bioconda flexbar`
+
+To get the latest version and best performance consider to build Flexbar from source.
 
 
 ### Building from source
@@ -30,20 +38,20 @@ Flexbar source code as well as binaries for Linux and Mac OS can be downloaded o
 Make sure that `cmake` is available, as well as development and runtime files of the TBB library 4.0 or later (Intel Threading Building Blocks). For example on Debian systems, install the packages `libtbb-dev` and `libtbb2`. Furthermore, the SeqAn library and a compiler that supports C++14 is required:
 
 * Get SeqAn library version 2.4.0 [here](https://github.com/seqan/seqan/releases/download/seqan-v2.4.0/seqan-library-2.4.0.tar.xz)
-* Download Flexbar 3.4 source code [release](https://github.com/seqan/flexbar/releases)
+* Download Flexbar 3.5.0 source code [release](https://github.com/seqan/flexbar/releases)
 
 Decompress both files:
 
-	tar xzf flexbar-3.4.0.tar.gz
+	tar xzf flexbar-3.5.0.tar.gz
 	tar xJf seqan-library-2.4.0.tar.xz
 
 Move SeqAn include folder to Flexbar:
 
-	mv seqan-library-2.4.0/include flexbar-3.4.0
+	mv seqan-library-2.4.0/include flexbar-3.5.0
 
 Use these commands for building:
 
-	cd flexbar-3.4.0
+	cd flexbar-3.5.0
 	cmake .
 	make
 
@@ -57,12 +65,12 @@ For execution of provided Flexbar binaries, the corresponding TBB library has to
 #### Linux
 Adjust lib search path to include the absolute path of the Flexbar directory containing the lib file libtbb.so.2 for the current terminal session, or permanently in shell startup scripts:
 
-	export LD_LIBRARY_PATH=/path/FlexbarDir:$LD_LIBRARY_PATH
+	export LD_LIBRARY_PATH=/YourPath/flexbar-3.5.0-linux:$LD_LIBRARY_PATH
 
 #### Mac OS
 It applies the same as for Linux. Make the file libtbb.dylib available by setting the lib search path:
 
-	export DYLD_LIBRARY_PATH=/path/FlexbarDir:$DYLD_LIBRARY_PATH
+	export DYLD_LIBRARY_PATH=/YourPath/flexbar-3.5.0-macos:$DYLD_LIBRARY_PATH
 
 
 ### Program usage
@@ -95,7 +103,7 @@ To remove adapter sequences from single-end reads, specify a file with adapters
 
 For paired-end libraries, specify both files with paired reads and a fasta file with adapters for removal. Given adapters are trimmed in right mode per default. It is recommended to activate the pair overlap detection in case of standard paired reads. This increases the sensitivity by removing very short parts of adapters if an overlap is detected for a pair.
 
-	flexbar -r r1.fq -p r2.fq -a adapters.fa -ap ON
+	flexbar -r r1.fq -p r2.fq -a a1.fa -a2 a2.fa -ap ON
 
 #### Adapter removal presets
 

debian/changelog

+flexbar (1:3.5.0-1) unstable; urgency=medium
+
+  * New upstream version
+  * debhelper-compat 12
+  * Standards-Version: 4.4.0
+  * Secure URI in copyright format
+
+ -- Andreas Tille <tille@debian.org>  Sat, 27 Jul 2019 21:27:50 +0200
+
 flexbar (1:3.4.0-2) unstable; urgency=medium
 
   * Team upload.

debian/compat

-11

debian/control

@@ -4,13 +4,13 @@ Uploaders: Andreas Tille <tille@debian.org>,
            Tony Travis <ajt@minke.ukfsn.org>
 Section: science
 Priority: optional
-Build-Depends: debhelper (>= 11~),
+Build-Depends: debhelper-compat (= 12),
                libtbb-dev,
                cmake,
                libseqan2-dev,
                zlib1g-dev,
                libbz2-dev
-Standards-Version: 4.1.5
+Standards-Version: 4.4.0
 Vcs-Browser: https://salsa.debian.org/med-team/flexbar
 Vcs-Git: https://salsa.debian.org/med-team/flexbar.git
 Homepage: https://github.com/seqan/flexbar

debian/copyright

-Format: http://www.debian.org/doc/packaging-manuals/copyright-format/1.0/
+Format: https://www.debian.org/doc/packaging-manuals/copyright-format/1.0/
 Upstream-Name: Flexbar
 Upstream-Contact: Johannes T. Roehr <johannes.roehr@fu-berlin.de>
 Source: https://sourceforge.net/projects/flexbar/files/

galaxy/flexbar.pl

+#!/usr/bin/env perl
+
+# Flexbar wrapper for Galaxy tool definition, version 3.5.0
+# Author: Johannes Roehr
+
+use warnings;
+use strict;
+
+my $format;
+my @inFiles;
+my @outFiles;
+
+my $compression = "";
+
+foreach(0..$#ARGV){
+	my $arg = $ARGV[$_];
+	
+	if($arg =~ /\.(fastq\w+)$/ || $arg =~ /\.(fastq\w+\.gz)$/ || $arg =~ /\.(fastq\w+\.bz2)$/){
+		
+		if(defined $format && $format ne $1){
+			warn "Read files should have the same format.\n";
+			exit 1;
+		}
+		$format = $1;
+		
+		my $file = $arg;
+		
+		$arg =~ s/\.fastq\w+$/\.fastq/;
+		$arg =~ s/\.fastq\w+\.gz$/\.fastq\.gz/;
+		$arg =~ s/\.fastq\w+\.bz2$/\.fastq\.bz2/;
+		
+		$compression = "GZ"  if $arg =~ /\.fastq\.gz$/;
+		$compression = "BZ2" if $arg =~ /\.fastq\.bz2$/;
+		
+		$ARGV[$_] = $arg;
+		
+		if($arg =~ /\.dat_input\w\.fastq$/ || $arg =~ /\.dat_input\w\.fastq\.gz$/ || $arg =~ /\.dat_input\w\.fastq\.bz2$/){
+			push @inFiles, $arg;
+			rename $file, $arg;
+		}
+		
+		push @outFiles, $arg if $arg =~ /\.dat\.fastq$/ || $arg =~ /\.dat\.fastq\.gz$/ || $arg =~ /\.dat\.fastq\.bz2$/;
+	}
+}
+
+my $barcoded = 0;
+
+$barcoded = 1 if $ARGV[$#ARGV] =~ /barcoded$/;
+
+my $call = join " ", @ARGV[0..($#ARGV - $barcoded)];
+
+# $call = $call ." --zip-output ". $compression if $barcoded && $compression ne "";
+
+system $call and exit 1;
+
+
+unlink $_ or warn "Could not unlink $_: $!" foreach(@inFiles);
+
+if($barcoded){
+	$format =~ s/\.gz//;
+	$format =~ s/\.bz2//;
+	
+	foreach(<$ARGV[$#ARGV]/flexbarOut*.fastq*>){
+		
+		my $file = $_;
+		
+		s/fastq$/$format/;
+		# s/fastq\.gz$/$format/;
+		# s/fastq\.bz2$/$format/;
+		
+		rename $file, $_;
+	}
+}
+else{
+	foreach(@outFiles){
+		
+		my $file = $_;
+		
+		s/\.fastq$//;
+		s/\.fastq\.gz$//;
+		s/\.fastq\.bz2$//;
+		
+		rename $file, $_;
+	}
+}
+

src/Flexbar.cpp

@@ -2,7 +2,7 @@
    
    Flexbar - flexible barcode and adapter removal
    
-   Version 3.4.0
+   Version 3.5.0
    
    BSD 3-Clause License
    
@@ -29,8 +29,8 @@ int main(int argc, const char* argv[]){
 	using namespace std;
 	using namespace seqan;
 	
-	const string version = "3.4.0";
-	const string date    = "June 2018";
+	const string version = "3.5.0";
+	const string date    = "May 2019";
 	
 	ArgumentParser parser("flexbar");
 	

src/FlexbarIO.h

@@ -295,8 +295,7 @@ void checkInputType(const std::string path, flexbar::FileFormat &format, const b
 		try{
 			if(! atEnd(seqFileIn)){
 				
-				FSeqStr seq;
-				FString id, qual;
+				FString id, seq, qual;
 				
 				readRecord(id, seq, qual, seqFileIn);
 				

src/FlexbarTypes.h

@@ -11,14 +11,15 @@ class SeqRead {
 	TSeqStr seq;
 	TString id, qual, umi;
 	
-	bool rmAdapter, rmAdapterRC, pairOverlap;
+	bool rmAdapter, rmAdapterRC, pairOverlap, poRemoval;
 	
 	SeqRead(TSeqStr& sequence, TString& seqID) :
 		seq(sequence),
 		id(seqID),
 		rmAdapter(false),
 		rmAdapterRC(false),
-		pairOverlap(false){
+		pairOverlap(false),
+		poRemoval(false){
 	}
 	
 	SeqRead(TSeqStr& sequence, TString& seqID, TString& quality) :
@@ -27,7 +28,8 @@ class SeqRead {
 		qual(quality),
 		rmAdapter(false),
 		rmAdapterRC(false),
-		pairOverlap(false){
+		pairOverlap(false),
+		poRemoval(false){
 	}
 };
 
@@ -233,6 +235,12 @@ namespace flexbar{
 		AOFF
 	};
 	
+	enum AdapterTrimmed {
+		ATON,
+		ATOFF,
+		ATONLY
+	};
+	
 	enum RunType {
 		SINGLE,
 		PAIRED,

src/Options.h

@@ -22,8 +22,8 @@ struct Options{
 	
 	bool isPaired, useAdapterFile, useNumberTag, useRemovalTag, umiTags, logStdout;
 	bool switch2Fasta, writeUnassigned, writeSingleReads, writeSingleReadsP, writeLengthDist;
-	bool useStdin, useStdout, relaxRegion, useRcTrimEnd, qtrimPostRm;
-	bool interleavedInput, htrimAdapterRm, htrimMaxFirstOnly;
+	bool useStdin, useStdout, relaxRegion, useRcTrimEnd, qtrimPostRm, addBarcodeAdapter;
+	bool interleavedInput, iupacInput, htrimAdapterRm, htrimMaxFirstOnly;
 	
 	int cutLen_begin, cutLen_end, cutLen_read, a_tail_len, b_tail_len, p_min_overlap;
 	int qtrimThresh, qtrimWinSize, a_overhang, htrimMinLength, htrimMinLength2, htrimMaxLength;
@@ -44,6 +44,7 @@ struct Options{
 	flexbar::RevCompMode     rcMode;
 	flexbar::PairOverlap     poMode;
 	flexbar::AdapterPreset   aPreset;
+	flexbar::AdapterTrimmed  aTrimmed;
 	
 	tbb::concurrent_vector<flexbar::TBar> barcodes, adapters, barcodes2, adapters2;
 	
@@ -80,10 +81,12 @@ struct Options{
 		logStdout         = false;
 		umiTags           = false;
 		interleavedInput  = false;
+		iupacInput        = false;
 		useStdin          = false;
 		useStdout         = false;
 		relaxRegion       = false;
 		useRcTrimEnd      = false;
+		addBarcodeAdapter = false;
 		qtrimPostRm       = false;
 		htrimAdapterRm    = false;
 		htrimMaxFirstOnly = false;
@@ -114,6 +117,7 @@ struct Options{
 		arc_end   = RIGHT;
 		b_end     = LTAIL;
 		aPreset   = APOFF;
+		aTrimmed  = ATON;
     }
 };
 
@@ -149,7 +153,7 @@ const std::string getFlexbarURL(){
 
 
 const std::string getFlexbarDescription(){
-	return "The program Flexbar preprocesses high-throughput sequencing data efficiently. It demultiplexes barcoded runs and removes adapter sequences. Moreover, trimming and filtering features are provided. Flexbar increases read mapping rates and improves genome as well as transcriptome assemblies. Unique molecular identifiers can be extracted in a flexible way. The program supports sequencing data in fasta and fastq format, e.g. from the Illumina platform. Refer to the manual on github.com/seqan/flexbar/wiki or contact Johannes Roehr on github.com/jtroehr for support with this application.";
+	return "The program Flexbar preprocesses high-throughput sequencing data efficiently. It demultiplexes barcoded runs and removes adapter sequences. Several adapter removal presets for Illumina libraries are included. Flexbar computes exact overlap alignments using SIMD and multicore parallelism. Moreover, trimming and filtering features are provided, e.g. trimming of homopolymers at read ends. Flexbar increases read mapping rates and improves genome as well as transcriptome assemblies. Unique molecular identifiers can be extracted in a flexible way. The software supports data in fasta and fastq format from multiple sequencing platforms. Refer to the manual on github.com/seqan/flexbar/wiki or contact Johannes Roehr on github.com/jtroehr for support with this application.";
 }
 
 
@@ -184,6 +188,7 @@ void defineOptions(seqan::ArgumentParser &parser, const std::string version, con
 	addOption(parser, ArgParseOption("r", "reads", "Fasta/q file or stdin (-) with reads that may contain barcodes.", ARG::INPUT_FILE));
 	addOption(parser, ArgParseOption("p", "reads2", "Second input file of paired reads, gz and bz2 files supported.", ARG::INPUT_FILE));
 	addOption(parser, ArgParseOption("i", "interleaved", "Interleaved format for first input set with paired reads."));
+	addOption(parser, ArgParseOption("I", "iupac", "Accept iupac symbols in reads and convert to N if not ATCG."));
 	
 	addSection(parser, "Barcode detection");
 	addOption(parser, ArgParseOption("b",  "barcodes", "Fasta file with barcodes for demultiplexing, may contain N.", ARG::INPUT_FILE));
@@ -214,7 +219,9 @@ void defineOptions(seqan::ArgumentParser &parser, const std::string version, con
 	addOption(parser, ArgParseOption("av", "adapter-min-poverlap", "Minimum overlap of paired reads for detection.", ARG::INTEGER));
 	addOption(parser, ArgParseOption("ac", "adapter-revcomp", "Include reverse complements of adapters.", ARG::STRING));
 	addOption(parser, ArgParseOption("ad", "adapter-revcomp-end", "Use different trim-end for reverse complements of adapters.", ARG::STRING));
+	addOption(parser, ArgParseOption("ab", "adapter-add-barcode", "Add reverse complement of detected barcode to adapters."));
 	addOption(parser, ArgParseOption("ar", "adapter-read-set", "Consider only single read set for adapters.", ARG::STRING));
+	addOption(parser, ArgParseOption("ak", "adapter-trimmed-out", "Modify that trimmed reads are kept.", ARG::STRING));
 	addOption(parser, ArgParseOption("ay", "adapter-cycles", "Number of adapter removal cycles.", ARG::INTEGER));
 	addOption(parser, ArgParseOption("am", "adapter-match", "Alignment match score.", ARG::INTEGER));
 	addOption(parser, ArgParseOption("ai", "adapter-mismatch", "Alignment mismatch score.", ARG::INTEGER));
@@ -257,7 +264,7 @@ void defineOptions(seqan::ArgumentParser &parser, const std::string version, con
 	addOption(parser, ArgParseOption("z", "zip-output", "Direct compression of output files.", ARG::STRING));
 	addOption(parser, ArgParseOption("1", "stdout-reads", "Write reads to stdout, tagged and interleaved if needed."));
 	addOption(parser, ArgParseOption("R", "output-reads", "Output file for reads instead of target prefix usage.", ARG::OUTPUT_FILE));
-	addOption(parser, ArgParseOption("D", "output-reads2", "Output file for reads2 instead of target prefix usage.", ARG::OUTPUT_FILE));
+	addOption(parser, ArgParseOption("P", "output-reads2", "Output file for reads2 instead of target prefix usage.", ARG::OUTPUT_FILE));
 	addOption(parser, ArgParseOption("j", "length-dist", "Generate length distribution for read output files."));
 	addOption(parser, ArgParseOption("s", "single-reads", "Write single reads for too short counterparts in pairs."));
 	addOption(parser, ArgParseOption("S", "single-reads-paired", "Write paired single reads with N for short counterparts."));
@@ -287,6 +294,8 @@ void defineOptions(seqan::ArgumentParser &parser, const std::string version, con
 	setAdvanced(parser, "adapter-min-poverlap");
 	setAdvanced(parser, "adapter-revcomp");
 	setAdvanced(parser, "adapter-revcomp-end");
+	setAdvanced(parser, "adapter-add-barcode");
+	setAdvanced(parser, "adapter-trimmed-out");
 	setAdvanced(parser, "adapter-read-set");
 	setAdvanced(parser, "adapter-cycles");
 	setAdvanced(parser, "adapter-match");
@@ -308,6 +317,7 @@ void defineOptions(seqan::ArgumentParser &parser, const std::string version, con
 	setAdvanced(parser, "bundle");
 	setAdvanced(parser, "bundles");
 	setAdvanced(parser, "interleaved");
+	setAdvanced(parser, "iupac");
 	setAdvanced(parser, "length-dist");
 	setAdvanced(parser, "single-reads");
 	setAdvanced(parser, "single-reads-paired");
@@ -359,6 +369,7 @@ void defineOptions(seqan::ArgumentParser &parser, const std::string version, con
 	
 	setValidValues(parser, "adapter-read-set", "1 2");
 	setValidValues(parser, "adapter-revcomp", "ON ONLY");
+	setValidValues(parser, "adapter-trimmed-out", "OFF ONLY");
 	setValidValues(parser, "adapter-pair-overlap", "ON SHORT ONLY");
 	setValidValues(parser, "adapter-preset", "TruSeq SmallRNA Methyl Ribo Nextera NexteraMP");
 	
@@ -403,7 +414,7 @@ void defineOptions(seqan::ArgumentParser &parser, const std::string version, con
 	addText(parser._toolDoc, "\\fBflexbar\\fP \\fB-r\\fP reads.fq \\fB-t\\fP target \\fB-q\\fP TAIL \\fB-qf\\fP i1.8", false);
 	addText(parser._toolDoc, "\\fBflexbar\\fP \\fB-r\\fP reads.fq \\fB-b\\fP barcodes.fa \\fB-bt\\fP LTAIL", false);
 	addText(parser._toolDoc, "\\fBflexbar\\fP \\fB-r\\fP reads.fq \\fB-a\\fP adapters.fa \\fB-ao\\fP 3 \\fB-ae\\fP 0.1", false);
-	addText(parser._toolDoc, "\\fBflexbar\\fP \\fB-r\\fP r1.fq \\fB-p\\fP r2.fq \\fB-a\\fP adapters.fa \\fB-ap\\fP ON", false);
+	addText(parser._toolDoc, "\\fBflexbar\\fP \\fB-r\\fP r1.fq \\fB-p\\fP r2.fq \\fB-a\\fP a1.fa \\fB-a2\\fP a2.fa \\fB-ap\\fP ON", false);
 	addText(parser._toolDoc, "\\fBflexbar\\fP \\fB-r\\fP r1.fq \\fB-p\\fP r2.fq \\fB-aa\\fP TruSeq \\fB-ap\\fP ON");
 }
 
@@ -613,6 +624,11 @@ void loadOptions(Options &o, seqan::ArgumentParser &parser){
 		}
 	}
 	
+	if(isSet(parser, "iupac")){
+		*out << "Iupac reads:           on" << endl;
+		o.iupacInput = true;
+	}
+	
 	
 	// barcode and adapter file options
 	
@@ -871,12 +887,10 @@ void loadOptions(Options &o, seqan::ArgumentParser &parser){
 		
 		if(isSet(parser, "output-reads") && (isSet(parser, "output-reads2") || o.runType == SINGLE)){
 			getOptionValue(o.outReadsFile, parser, "output-reads");
-			*out << "output-reads:          " << o.outReadsFile << endl;
 		}
 		
 		if(isSet(parser, "output-reads2") && isSet(parser, "output-reads") && o.runType == PAIRED){
 			getOptionValue(o.outReadsFile2, parser, "output-reads2");
-			*out << "output-reads2:         " << o.outReadsFile2 << endl;
 			
 			if(o.outReadsFile == o.outReadsFile2){
 				cerr << "\n" << "Output reads and reads2 file should not be the same.\n" << endl;
@@ -1043,6 +1057,22 @@ void loadOptions(Options &o, seqan::ArgumentParser &parser){
 				o.relaxRegion = true;
 			}
 			
+			if(isSet(parser, "adapter-add-barcode") && o.isPaired && o.a_end == RIGHT && o.rcMode != RCON &&
+				o.barDetect != BARCODE_READ && o.barDetect != BOFF && o.b_end == LTAIL){
+				
+				*out << "adapter-add-barcode:   on" << endl;
+				o.addBarcodeAdapter = true;
+			}
+			
+			if(isSet(parser, "adapter-trimmed-out")){
+				string a_trimmed_out;
+				getOptionValue(a_trimmed_out, parser, "adapter-trimmed-out");
+				*out << "adapter-trimmed-out:   " << a_trimmed_out << endl;
+				
+				     if(a_trimmed_out == "OFF")  o.aTrimmed = ATOFF;
+				else if(a_trimmed_out == "ONLY") o.aTrimmed = ATONLY;
+			}
+			
 			if(isSet(parser, "adapter-read-set") && o.isPaired && o.adapRm != NORMAL2){
 				string a_read_set;
 				getOptionValue(a_read_set, parser, "adapter-read-set");

src/PairedAlign.h

@@ -14,7 +14,7 @@ class PairedAlign : public tbb::filter {
 private:
 	
 	const bool m_writeUnassigned, m_twoBarcodes, m_umiTags, m_useRcTrimEnd;
-	const bool m_htrim, m_htrimAdapterRm, m_htrimMaxFirstOnly;
+	const bool m_htrim, m_htrimAdapterRm, m_htrimMaxFirstOnly, m_addBarcodeAdapter;
 	
 	const std::string m_htrimLeft, m_htrimRight;
 	
@@ -61,6 +61,7 @@ public:
 		m_umiTags(o.umiTags),
 		m_useRcTrimEnd(o.useRcTrimEnd),
 		m_writeUnassigned(o.writeUnassigned),
+		m_addBarcodeAdapter(o.addBarcodeAdapter),
 		m_htrimLeft(o.htrimLeft),
 		m_htrimRight(o.htrimRight),
 		m_htrimMinLength(o.htrimMinLength),
@@ -106,11 +107,11 @@ public:
 		using namespace flexbar;
 		
 		switch(m_barType){
-			case BARCODE_READ:         pRead->barID  = m_b1->alignSeqRead(pRead->b,  false, alBundle[0], cycle[0], idxAl[0], alMode, m_bTrimEnd); break;
-			case WITHIN_READ_REMOVAL2: pRead->barID2 = m_b2->alignSeqRead(pRead->r2, true,  alBundle[2], cycle[2], idxAl[2], alMode, m_bTrimEnd);
-			case WITHIN_READ_REMOVAL:  pRead->barID  = m_b1->alignSeqRead(pRead->r1, true,  alBundle[1], cycle[1], idxAl[1], alMode, m_bTrimEnd); break;
-			case WITHIN_READ2:         pRead->barID2 = m_b2->alignSeqRead(pRead->r2, false, alBundle[2], cycle[2], idxAl[2], alMode, m_bTrimEnd);
-			case WITHIN_READ:          pRead->barID  = m_b1->alignSeqRead(pRead->r1, false, alBundle[1], cycle[1], idxAl[1], alMode, m_bTrimEnd); break;
+			case BARCODE_READ:         pRead->barID  = m_b1->alignSeqRead(pRead->b,  false, alBundle[0], cycle[0], idxAl[0], alMode, m_bTrimEnd, ""); break;
+			case WITHIN_READ_REMOVAL2: pRead->barID2 = m_b2->alignSeqRead(pRead->r2, true,  alBundle[2], cycle[2], idxAl[2], alMode, m_bTrimEnd, "");
+			case WITHIN_READ_REMOVAL:  pRead->barID  = m_b1->alignSeqRead(pRead->r1, true,  alBundle[1], cycle[1], idxAl[1], alMode, m_bTrimEnd, ""); break;
+			case WITHIN_READ2:         pRead->barID2 = m_b2->alignSeqRead(pRead->r2, false, alBundle[2], cycle[2], idxAl[2], alMode, m_bTrimEnd, "");
+			case WITHIN_READ:          pRead->barID  = m_b1->alignSeqRead(pRead->r1, false, alBundle[1], cycle[1], idxAl[1], alMode, m_bTrimEnd, ""); break;
 			case BOFF: break;
 		}
 		
@@ -125,12 +126,49 @@ public:
 		
 		using namespace flexbar;
 		
-		if(m_adapRem != ATWO)
-			m_a1->alignSeqRead(pRead->r1, true, alBundle[0], cycle[0], idxAl[0], alMode, trimEnd);
+		if(m_adapRem != ATWO){
+			
+			TSeqStr addBarcode = "";
+			
+			if(m_addBarcodeAdapter && pRead->r2 != NULL && pRead->barID2 > 0){
+				addBarcode = m_barcodes2->at(pRead->barID2 - 1).seq;
+				
+				if(m_umiTags && pRead->r2->umi != ""){
+					unsigned int umiPos = 1;
+					
+					for(unsigned int i = 0; i < length(addBarcode); ++i){
+						if(addBarcode[i] == 'N' && length(pRead->r2->umi) > umiPos){
+							addBarcode[i] = pRead->r2->umi[umiPos++];
+						}
+					}
+				}
+				seqan::reverseComplement(addBarcode);
+			}
+			
+			m_a1->alignSeqRead(pRead->r1, true, alBundle[0], cycle[0], idxAl[0], alMode, trimEnd, addBarcode);
+		}
 		
 		if(pRead->r2 != NULL && m_adapRem != AONE){
-			if(m_adapRem != NORMAL2) m_a1->alignSeqRead(pRead->r2, true, alBundle[1], cycle[1], idxAl[1], alMode, trimEnd);
-			else                     m_a2->alignSeqRead(pRead->r2, true, alBundle[1], cycle[1], idxAl[1], alMode, trimEnd);
+			
+			TSeqStr addBarcode = "";
+			
+			if(m_addBarcodeAdapter && pRead->barID > 0){
+				addBarcode = m_barcodes->at(pRead->barID - 1).seq;
+				
+				if(m_umiTags && pRead->r1->umi != ""){
+					unsigned int umiPos = 1;
+					
+					for(unsigned int i = 0; i < length(addBarcode); ++i){
+						if(addBarcode[i] == 'N' && length(pRead->r1->umi) > umiPos){
+							addBarcode[i] = pRead->r1->umi[umiPos++];
+						}
+					}
+				}
+				seqan::reverseComplement(addBarcode);
+			}
+			
+			if(m_adapRem != NORMAL2) m_a1->alignSeqRead(pRead->r2, true, alBundle[1], cycle[1], idxAl[1], alMode, trimEnd, addBarcode);
+			else                     m_a2->alignSeqRead(pRead->r2, true, alBundle[1], cycle[1], idxAl[1], alMode, trimEnd, addBarcode);
 		}
 	}
 	

src/PairedOutput.h

@@ -26,6 +26,7 @@ private:
 	const flexbar::RunType        m_runType;
 	const flexbar::BarcodeDetect  m_barDetect;
 	const flexbar::QualTrimType   m_qtrim;
+	const flexbar::AdapterTrimmed m_aTrimmed;
 	
 	typedef SeqOutput<TSeqStr, TString>      TSeqOutput;
 	typedef SeqOutputFiles<TSeqStr, TString> TOutFiles;
@@ -50,6 +51,7 @@ public:
 		m_qtrimThresh(o.qtrimThresh),
 		m_qtrimWinSize(o.qtrimWinSize),
 		m_qtrimPostRm(o.qtrimPostRm),
+		m_aTrimmed(o.aTrimmed),
 		m_isPaired(o.isPaired),
 		m_writeUnassigned(o.writeUnassigned),
 		m_writeSingleReads(o.writeSingleReads),
@@ -239,7 +241,7 @@ public:
 		
 		using namespace flexbar;
 		
-		bool l1ok = false, l2ok = false;
+		bool r1ok = false, r2ok = false;
 		
 		switch(m_runType){
 			
@@ -253,11 +255,13 @@ public:
 							if(qualTrim(pRead->r1, m_qtrim, m_qtrimThresh, m_qtrimWinSize)) ++m_nLowPhred;
 						}
 						
-						if(length(pRead->r1->seq) >= m_minLength){
-							
-							m_outMap[pRead->barID].f1->writeRead(pRead->r1);
-						}
+						if(length(pRead->r1->seq) >= m_minLength) r1ok = true;
 						else m_outMap[pRead->barID].m_nShort_1++;
+						
+						if     (m_aTrimmed == ATOFF  &&  (pRead->r1->rmAdapter ||   pRead->r1->rmAdapterRC)) r1ok = false;
+						else if(m_aTrimmed == ATONLY && ! pRead->r1->rmAdapter && ! pRead->r1->rmAdapterRC)  r1ok = false;
+						
+						if(r1ok) m_outMap[pRead->barID].f1->writeRead(pRead->r1);
 					}
 				}
 				break;
@@ -284,14 +288,23 @@ public:
 							if(qualTrim(pRead->r2, m_qtrim, m_qtrimThresh, m_qtrimWinSize)) ++m_nLowPhred;
 						}
 						
-						if(length(pRead->r1->seq) >= m_minLength) l1ok = true;
-						if(length(pRead->r2->seq) >= m_minLength) l2ok = true;
+						if(length(pRead->r1->seq) >= m_minLength) r1ok = true;
+						if(length(pRead->r2->seq) >= m_minLength) r2ok = true;
 						
-						if(l1ok && l2ok){
+						if(! r1ok) m_outMap[outIdx].m_nShort_1++;
+						if(! r2ok) m_outMap[outIdx].m_nShort_2++;
+						
+						if     (m_aTrimmed == ATOFF  &&  (pRead->r1->rmAdapter ||   pRead->r1->rmAdapterRC ||   pRead->r1->poRemoval)) r1ok = false;
+						else if(m_aTrimmed == ATONLY && ! pRead->r1->rmAdapter && ! pRead->r1->rmAdapterRC && ! pRead->r1->poRemoval)  r1ok = false;
+						
+						if     (m_aTrimmed == ATOFF  &&  (pRead->r2->rmAdapter ||   pRead->r2->rmAdapterRC ||   pRead->r2->poRemoval)) r2ok = false;
+						else if(m_aTrimmed == ATONLY && ! pRead->r2->rmAdapter && ! pRead->r2->rmAdapterRC && ! pRead->r2->poRemoval)  r2ok = false;
+						
+						if(r1ok && r2ok){
 							m_outMap[outIdx].f1->writeRead(pRead->r1);
 							m_outMap[outIdx].f2->writeRead(pRead->r2);
 						}
-						else if(l1ok && ! l2ok){
+						else if(r1ok && ! r2ok){
 							m_nSingleReads++;
 							
 							if(m_writeSingleReads){
@@ -308,7 +321,7 @@ public:
 								m_outMap[outIdx].f2->writeRead(pRead->r2);
 							}
 						}
-						else if(! l1ok && l2ok){
+						else if(! r1ok && r2ok){
 							m_nSingleReads++;
 							
 							if(m_writeSingleReads){
@@ -325,9 +338,6 @@ public:
 								m_outMap[outIdx].f2->writeRead(pRead->r2);
 							}
 						}
-						
-						if(! l1ok) m_outMap[outIdx].m_nShort_1++;
-						if(! l2ok) m_outMap[outIdx].m_nShort_2++;
 					}
 				}
 			}

src/QualTrimming.h

@@ -44,7 +44,7 @@ unsigned qualTrimming(const TString& qual, unsigned const _cutoff, Window const
 	unsigned window = spec.size;
 	unsigned avg = 0, i = 0;
 	
-	// Work with absolute cutoff in window to avoid divisions
+	// Absolute cutoff in window to avoid divisions
 	unsigned cutoff = _cutoff * window;
 	
 	// Calculate average quality of initial window
@@ -57,9 +57,15 @@ unsigned qualTrimming(const TString& qual, unsigned const _cutoff, Window const
 		
 		// Take care only not to go over the end of the sequence. Shorten window near the end
 		avg -= getQuality(qual, i);
-		avg += i + window < length(qual) ? getQuality(qual, i + window) : 0;
+		
+		if(i + window < length(qual)){
+			avg += getQuality(qual, i + window);
+		}
+		else{
+			cutoff = _cutoff * ((length(qual) - 1) - i);
+		}
 	}
-	return i;   // i holds start of first window that turned bad
+	return i;   // holds start of first window that turned bad
 }
 
 

src/SeqAlign.h

@@ -15,7 +15,7 @@ private:
 	const flexbar::FileFormat  m_format;
 	const flexbar::PairOverlap m_poMode;
 	
-	const bool m_isBarcoding, m_writeTag, m_umiTags, m_strictRegion;
+	const bool m_isBarcoding, m_writeTag, m_umiTags, m_strictRegion, m_addBarcodeAdapter;
 	const int m_minLength, m_minOverlap, m_tailLength;
 	const float m_errorRate;
 	const unsigned int m_bundleSize;
@@ -41,6 +41,7 @@ public:
 			m_log(o.logAlign),
 			m_format(o.format),
 			m_writeTag(o.useRemovalTag),
+			m_addBarcodeAdapter(o.addBarcodeAdapter),
 			m_strictRegion(! o.relaxRegion),
 			m_bundleSize(o.bundleSize),
 			m_out(o.out),
@@ -53,7 +54,7 @@ public:
 	};
 	
 	
-	int alignSeqRead(flexbar::TSeqRead* sr, const bool performRemoval, flexbar::Alignments &alignments, flexbar::ComputeCycle &cycle, unsigned int &idxAl, const flexbar::AlignmentMode &alMode, const flexbar::TrimEnd trimEnd){
+	int alignSeqRead(flexbar::TSeqRead* sr, const bool performRemoval, flexbar::Alignments &alignments, flexbar::ComputeCycle &cycle, unsigned int &idxAl, const flexbar::AlignmentMode &alMode, const flexbar::TrimEnd trimEnd, const TSeqStr &addBarcode){
 		
 		using namespace std;
 		using namespace flexbar;
@@ -81,12 +82,18 @@ public:
 				if     (alMode == ALIGNRCOFF &&   m_queries->at(i).rcAdapter) continue;
 				else if(alMode == ALIGNRC    && ! m_queries->at(i).rcAdapter) continue;
 				
-				TSeqStr &qseq = m_queries->at(i).seq;
+				TSeqStr *qseq = &m_queries->at(i).seq;
 				TSeqStr *rseq = &seqRead.seq;
-				TSeqStr tmp;
+				TSeqStr tmp, tmpq;
+				
+				if(! m_isBarcoding && m_addBarcodeAdapter && addBarcode != ""){
+					tmpq = addBarcode;
+					append(tmpq, m_queries->at(i).seq);
+					qseq = &tmpq;
+				}
 				
 				if(trimEnd == LTAIL || trimEnd == RTAIL){
-					int tailLength  = (m_tailLength > 0) ? m_tailLength : length(qseq);
+					int tailLength  = (m_tailLength > 0) ? m_tailLength : length(*qseq);
 					
 					if(tailLength < readLength){
 						if(trimEnd == LTAIL) tmp = prefix(seqRead.seq, tailLength);
@@ -100,7 +107,7 @@ public:
 				resize(rows(alignments.aset[idxAl]), 2);
 				
 				assignSource(row(alignments.aset[idxAl], 0), *rseq);
-				assignSource(row(alignments.aset[idxAl], 1),  qseq);
+				assignSource(row(alignments.aset[idxAl], 1), *qseq);
 				
 				++idxAl;
 			}
@@ -124,6 +131,11 @@ public:
 			m_algo.alignGlobal(a, alignments, cycle, idxAl++, trimEnd);
 			
 			a.queryLength = length(m_queries->at(i).seq);
+			
+			if(! m_isBarcoding && m_addBarcodeAdapter && addBarcode != ""){
+				a.queryLength += length(addBarcode);
+			}
+			
 			a.tailLength  = (m_tailLength > 0) ? m_tailLength : a.queryLength;
 			
 			a.overlapLength = a.endPos - a.startPos;

src/SeqAlignPair.h

@@ -114,6 +114,8 @@ public:
 					
 					++m_modified;
 					
+					seqRead2.poRemoval = true;
+					
 					if(m_writeTag) append(seqRead2.id, "_Flexbar_removal_PO");
 				}
 			}
@@ -132,6 +134,8 @@ public:
 					
 					++m_modified;
 					
+					seqRead.poRemoval = true;
+					
 					if(m_writeTag) append(seqRead.id, "_Flexbar_removal_PO");
 				}
 			}

src/SeqInput.h

@@ -16,7 +16,7 @@ private:
 	const flexbar::QualTrimType m_qtrim;
 	const flexbar::FileFormat m_format;
 	
-	const bool m_preProcess, m_useStdin, m_qtrimPostRm;
+	const bool m_preProcess, m_useStdin, m_qtrimPostRm, m_iupacInput;
 	const int m_maxUncalled, m_preTrimBegin, m_preTrimEnd, m_qtrimThresh, m_qtrimWinSize;
 	tbb::atomic<unsigned long> m_nrReads, m_nrChars, m_nLowPhred;
 	
@@ -33,6 +33,7 @@ public:
 		m_qtrimThresh(o.qtrimThresh),
 		m_qtrimWinSize(o.qtrimWinSize),
 		m_qtrimPostRm(o.qtrimPostRm),
+		m_iupacInput(o.iupacInput),
 		m_format(o.format),
 		m_nrReads(0),
 		m_nrChars(0),
@@ -76,6 +77,8 @@ public:
 				reserve(seqs,     nReads);
 				reserve(uncalled, nReads);
 				
+				if(! m_iupacInput){
+					
 					if(m_format == FASTA){
 						readRecords(ids, seqs, seqFileIn, nReads);
 					}
@@ -83,6 +86,22 @@ public:
 						reserve(quals, nReads);
 						readRecords(ids, seqs, quals, seqFileIn, nReads);
 					}
+				}
+				else{
+					seqan::StringSet<seqan::IupacString> seqsIupac;
+					
+					reserve(seqsIupac, nReads);
+					
+					if(m_format == FASTA){
+						readRecords(ids, seqsIupac, seqFileIn, nReads);
+					}
+					else{
+						reserve(quals, nReads);
+						readRecords(ids, seqsIupac, quals, seqFileIn, nReads);
+					}
+					
+					seqs = seqsIupac;
+				}
 				
 				for(unsigned int i = 0; i < length(ids); ++i){
 					
@@ -133,6 +152,7 @@ public:
 						}
 					}
 				}
+				
 				m_nrReads += length(ids);
 				
 				return length(ids);