[med-svn] [hilive] 01/03: New upstream version 1.0
Andreas Tille
tille at debian.org
Tue Aug 22 14:00:51 UTC 2017
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository hilive.
commit 9eef8b7b361c1f0bcbfb649fb1b6a1e58b077857
Author: Andreas Tille <tille at debian.org>
Date: Tue Aug 22 09:54:20 2017 +0200
New upstream version 1.0
---
CMakeLists.txt | 67 ++--
CONTRIBUTORS | 12 +-
LICENSE | 2 +-
README.md | 80 +++--
lib/alignmentSettings.h | 773 ++++++++++++++++++++++++++++++++++++++++++
lib/alnread.cpp | 874 +++++++++++++++++++++++++++++++++++-------------
lib/alnread.h | 91 +++--
lib/alnstream.cpp | 675 ++++++++++++++++++++++---------------
lib/alnstream.h | 39 ++-
lib/argument_parser.cpp | 778 +++++++++++++++++++++++++++++++-----------
lib/argument_parser.h | 481 +++++++++++++++++++++++++-
lib/config.h.in | 7 -
lib/definitions.h | 378 +++++++++++++++------
lib/global_variables.h | 7 +
lib/headers.h | 9 +-
lib/kindex.cpp | 263 ++++++++-------
lib/kindex.h | 36 +-
lib/parallel.cpp | 183 +++++-----
lib/parallel.h | 80 +++--
lib/tools.cpp | 296 ++++++++--------
lib/tools.h | 161 +++++++--
lib/tools_static.cpp | 181 ++++++++++
lib/tools_static.h | 171 ++++++++++
tools/build_index.cpp | 147 +++-----
tools/hilive.cpp | 265 +++++++++------
tools/hilive_out.cpp | 66 ++++
26 files changed, 4518 insertions(+), 1604 deletions(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 4c6b203..1c3eced 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -5,44 +5,16 @@ cmake_minimum_required (VERSION 2.8)
project (HiLive)
# Set the version number
-set (HiLive_VERSION_MAJOR 0)
-set (HiLive_VERSION_MINOR 3)
-
-# Set the k-mer length
-set (HiLive_K 15 CACHE INT "Set the k-mer length for index and mapper")
+add_definitions(-DHiLive_VERSION_MAJOR=1)
+add_definitions(-DHiLive_VERSION_MINOR=0)
# Set flags for compilation
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fPIC -g -pthread -W -Wall -std=gnu++11 -O0")
-
-
-# configure a header file to pass some of the CMake settings to the source code
-configure_file (
- "${PROJECT_SOURCE_DIR}/lib/config.h.in"
- "${PROJECT_BINARY_DIR}/config.h" )
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fPIC -g -pthread -W -Wall -std=gnu++14 -O0")
# add the binary tree to the search path for include files
include_directories("${PROJECT_BINARY_DIR}")
-
-#############################
-### setup seqan libraries ###
-
-set (CMAKE_MODULE_PATH "/home/schulzeja/masterthesis/tools/seqan-src/util/cmake") # adjust this to [pathToSeqanCode]/util/cmake
-set (SEQAN_INCLUDE_PATH "/home/schulzeja/masterthesis/tools/seqan-src/include/") # adjust this to [pathToSeqanCode]/include
-
-# Configure SeqAn, enabling features for libbz2 and zlib.
-#set (SEQAN_FIND_DEPENDENCIES ZLIB BZip2) # original version from seqan tutorial
-set (SEQAN_FIND_DEPENDENCIES BZip2)
-find_package (SeqAn REQUIRED)
-
-# Add include directories, defines, and flags for SeqAn (and its dependencies).
-include_directories (${SEQAN_INCLUDE_DIRS})
-add_definitions (${SEQAN_DEFINITIONS})
-set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${SEQAN_CXX_FLAGS}")
-
-
-
#############################
### setup Boost libraries ###
@@ -53,32 +25,38 @@ find_package( Boost COMPONENTS system filesystem program_options REQUIRED )
include_directories( ${Boost_INCLUDE_DIR} )
-
############################
### setup Zlib and lz4 libraries ###
find_package( ZLIB REQUIRED )
-set (LZ4_PATH /usr/local/lib) # possibly adjust this to [pathToLz4]/lib if using downloaded lz4 source code
+#set (LZ4_PATH /usr/local/lib) # possibly adjust this to [pathToLz4]/lib if using downloaded lz4 source code
include_directories(${LZ4_PATH})
link_directories(${LZ4_PATH})
set(CompressionLibs "${ZLIB_LIBRARIES};lz4")
-#############################################
-### download & build the Samtools library ###
-file(DOWNLOAD "https://sourceforge.net/projects/samtools/files/samtools/1.2/samtools-1.2.tar.bz2/download" ${CMAKE_CURRENT_BINARY_DIR}/samtools-1.2.tar.bz2)
-execute_process(COMMAND tar xjf ${CMAKE_CURRENT_BINARY_DIR}/samtools-1.2.tar.bz2 )
-execute_process(COMMAND make -s -C samtools-1.2 )
-include_directories( ${CMAKE_CURRENT_BINARY_DIR}/samtools-1.2 )
-include_directories( ${CMAKE_CURRENT_BINARY_DIR}/samtools-1.2/htslib-1.2.1 )
-link_directories( ${CMAKE_CURRENT_BINARY_DIR}/samtools-1.2 ${CMAKE_CURRENT_BINARY_DIR}/samtools-1.2/htslib-1.2.1 )
+#############################
+### setup seqan libraries ### needs to be done AFTER searching for Zlib
+
+set (CMAKE_MODULE_PATH "/usr/local/lib/seqan/util/cmake") # adjust this to [pathToSeqanCode]/util/cmake
+set (SEQAN_INCLUDE_PATH "/usr/local/lib/seqan/include/") # adjust this to [pathToSeqanCode]/include
+
+# Configure SeqAn, enabling features for libbz2 and zlib.
+#set (SEQAN_FIND_DEPENDENCIES ZLIB BZip2) # original version from seqan tutorial
+set (SEQAN_FIND_DEPENDENCIES BZip2)
+find_package (SeqAn REQUIRED)
+
+# Add include directories, defines, and flags for SeqAn (and its dependencies).
+include_directories (${SEQAN_INCLUDE_DIRS})
+add_definitions (${SEQAN_DEFINITIONS})
+set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${SEQAN_CXX_FLAGS}")
##############################
### setup HiLive libraries ###
include_directories("${PROJECT_SOURCE_DIR}/lib")
# make a list of HiLive libraries
-set (LIB_NAMES tools alnread alnstream illumina_parsers kindex parallel argument_parser)
+set (LIB_NAMES tools_static tools alnread alnstream illumina_parsers kindex parallel argument_parser)
set(LIB_LIST "")
foreach (x ${LIB_NAMES})
list(APPEND LIB_LIST "lib/${x}.cpp")
@@ -90,11 +68,14 @@ add_library(HiLiveLibs ${LIB_LIST})
### Build the executables ###
add_executable (hilive tools/hilive.cpp)
-target_link_libraries (hilive HiLiveLibs ${CompressionLibs} ${Boost_LIBRARIES} ${SEQAN_LIBRARIES} hts)
+target_link_libraries (hilive HiLiveLibs ${CompressionLibs} ${Boost_LIBRARIES} ${SEQAN_LIBRARIES})
add_executable(hilive-build tools/build_index.cpp )
target_link_libraries(hilive-build HiLiveLibs ${CompressionLibs} ${Boost_LIBRARIES} ${SEQAN_LIBRARIES})
+add_executable(hilive-out tools/hilive_out.cpp )
+target_link_libraries(hilive-out HiLiveLibs ${CompressionLibs} ${Boost_LIBRARIES} ${SEQAN_LIBRARIES})
+
#####################
### for debugging ###
diff --git a/CONTRIBUTORS b/CONTRIBUTORS
index a9e3fc2..7f43b48 100644
--- a/CONTRIBUTORS
+++ b/CONTRIBUTORS
@@ -5,4 +5,14 @@ Martin S. Lindner <martin (at) mail-linder.de>
* Project founder, implemented versions 0.1 and 0.2
Jakob M. Schulze <jakobschulze (at) arcor.de>
- * developed version 0.2 to 0.3
+ * development v0.2 -> v0.3
+ * development support v0.3 -> v1.0
+
+ Tobias P. Loka <lokat (at) rki.de>
+ * development v0.3 -> v1.0
+
+ Simon H. Tausch <tauschs (at) rki.de>
+ * continuous development support
+
+ Kristina Kirsten
+ * development support v0.3 -> v1.0 (real-time SAM/BAM output)
\ No newline at end of file
diff --git a/LICENSE b/LICENSE
index dfe02f6..ce60a52 100644
--- a/LICENSE
+++ b/LICENSE
@@ -1,4 +1,4 @@
-Copyright (c) 2015-2016, Martin S. Lindner and the HiLive contributors. See CONTRIBUTORS for more info.
+Copyright (c) 2015-2017, Martin S. Lindner and the HiLive contributors. See CONTRIBUTORS for more info.
All rights reserved.
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
diff --git a/README.md b/README.md
index 6b47601..3d1e925 100644
--- a/README.md
+++ b/README.md
@@ -33,20 +33,21 @@ Make sure that the following dependencies are installed:
* zlib
* lz4
-If using a local version of lz4 then adjust path in CMakeLists.txt line 61.
+If using a local version of lz4 then adjust path in CMakeLists.txt line 32.
---
-You also need to download a specific seqan build. Navigate to a folder where
-you want to store the seqan library and run:
+You also need to download download seqan.
+Cloning the repository makes it possible to switch between different versions:
git clone https://github.com/seqan/seqan.git
- cd seqan && git checkout 9119aa6
+
+HiLive 1.0 was tested with SeqAn version 2.3.2.
---
Check out the HiLive source code from the project website and adjust the paths of the
-seqan module in the file CMakeLists.txt in the hilive folder (line 31 and 32).
+seqan module in the file CMakeLists.txt in the hilive folder (line 41 and 42).
Then, compile HiLive with:
cd [hilive-code]
@@ -54,36 +55,36 @@ Then, compile HiLive with:
cmake ..
make
----
-
-To compile HiLive with a different k-mer size than k=15 make the following
-adjustment (here: k=10):
-
- cmake -DHiLive\_K 10 ..
- make
-
-
Usage
-----
-HiLive has two components:
+HiLive has three components:
* ``hilive-build`` builds the k-mer index of the reference genome
* ``hilive`` the read mapper itself
+ * ``hilive-out`` executable to produce output files
---
#### Using hilive-build:
-Building a k-mer index from FASTA file input.fa to output file input.fa.kix:
+Building a k-mer index from FASTA file input.fa to output file input.fa.kix with k-mer weight 15:
- hilive-build input.fa
+ hilive-build input.fa 15
Building an index from a large reference genome. Here is makes sense to use trimming,
i.e. removing k-mers from the index that occurr more than 1000 times (for example) in
the index. The index is written into the file trimmed.kix
- hilive-build -t 1000 -o trimmed.kix input.fa
+ hilive-build -t 1000 -o trimmed.kix input.fa 15
+
+For gapped k-mers, use the -p parameter to specify the gap positions.
+For example, for a gap pattern of 1101110011, type:
+
+ hilive-build -p 3 7 8 input.fa 7
+
+With the current index structure, we strongly recommend to use a maximum k-mer weight of 15 because of huge disk space and memory requirements for large
+k-mers.
---
@@ -93,15 +94,33 @@ To map reads in a 100bp run using default settings:
hilive /path/to/BaseCalls /path/to/index.kix 100 /path/to/outputFolder
-Hint: to concatenate the resulting sam files of a lane to a single file using samtools you could use the following commands:
+For an overview of additional parameters, type
+
+ hilive --help
- for file in /path/to/outputFolder/L001/s_1_*.sam; do samtools view -bS $file > ${file%.sam}.bam; done
- samtools cat /path/to/outputFolder/L001/s_1_*.bam > /path/to/outputFolder/L001/allAlignments.bam
- samtools view -SH /path/to/outputFolder/L001/s_1_1101.sam > /path/to/outputFolder/L001/allAlignments.sam
- samtools view /path/to/outputFolder/L001/allAlignments.bam >> /path/to/outputFolder/L001/allAlignments.sam
+To prevent errors during argument parsing we recommend to set optional parameters AFTER the positional options:
+
+ hilive BC_DIR INDEX CYCLES OUTDIR [options]
+
+However, if unexpected parsing errors occur, please try to specify all parameters with the "--"-syntax (e.g. --BCDIR /path/to/BaseCalls) instead of using positional arguments. This is also necessary when loading (some of the) positional arguments from a settings file instead of using the command line.
---
+#### Using hilive-out:
+
+To create a SAM or BAM alignment output from existing temporary files in HiLive, type:
+
+ hilive-out --settings /path/to/temp/dir/hilive_settings.xml
+
+This will output the alignment results of the last cycle based on the settings that were specified for the related HiLive run.
+To produce output files for other cycles, e.g. 50, 70, 90, type:
+
+ hilive-out --settings /path/to/temp/dir/hilive_settings.xml --output-cycles 50 70 90
+
+Please note, that the temporary files for the respective cycles must be present in the temp folder.
+This is only the case if the --keep-files parameter and/or the --output-cycles parameter for the respective cycles was activated in the HiLive run.
+
+
#### Demultiplexing:
To map reads from multiplexed sequencing runs, you can provide HiLive with the barcode sequences from your Sample Sheet.
@@ -110,17 +129,11 @@ If you use double indexing, please concatenate both indices in the correct order
hilive /path/to/BaseCalls /path/to/index.kix 107 /path/to/outputFolder -b barcode1 -b barcode2 ...
-Reads containing one of the given barcodes will be written to the resulting samfiles. The corresponding barcode sequence is stored in the BC-field of the samfile.
-
-Hint: You can write the results of each sample into a separate files using the following command:
-
- for barcode in barcode1 barcode2 ... ; do grep "^@" allAlignments.sam > allAlignments_$barcode.sam && grep "BC:Z:$barcode" allAlignments.sam >> allAlignments_$barcode.sam; done
+One output file will be produced for each barcode. To get alignments with undetermined barcodes, activate the --keep-all-barcodes parameter.
+Dual barcodes must be delimited with "-" (e.g., -b ATCGTGAT-TAGTTAGC for a 2x8bp barcode).
---
-We recommend to adjust the numbers of threads used by HiLive with -n. If possible,
-make use of all threads on the machine!
-
Please consult the project website for more details on the parameters!
@@ -136,8 +149,7 @@ Contact
Please consult the HiLive project website for questions!
If this does not help, please feel free to consult:
- * [Martin S. Lindner](mailto:martin at mail-lindner.de "Mail Martin S. Lindner") (technical contact)
- * [Bernhard Y. Renard](mailto:renardb at rki.de "Mail Bernhard Y. Renard") (project head)
-
+ * Technical support <hilive.team (at) gmail.com> (technical contact)
+ * Bernhard Y. Renard <renardb (at) rki.de> (project head)
also see CONTRIBUTORS for a complete list of contributors and their contact information
diff --git a/lib/alignmentSettings.h b/lib/alignmentSettings.h
new file mode 100644
index 0000000..3b934ea
--- /dev/null
+++ b/lib/alignmentSettings.h
@@ -0,0 +1,773 @@
+#ifndef ALIGNMENTSETTINGS_H
+#define ALIGNMENTSETTINGS_H
+
+#include "headers.h"
+#include "definitions.h"
+#include "tools_static.h"
+
+// Data structure to store the alignment settings
+class AlignmentSettings {
+
+private:
+
+ // kmer gap positions
+ Unmodifiable<std::vector<unsigned>> kmer_gaps;
+
+ // reverse gap positions
+ Unmodifiable<std::vector<unsigned>> rev_kmer_gaps;
+
+ // PARAMETER: kmer span (automatically computed from kmer_weight and kmer_gaps)
+ Unmodifiable<uint8_t> kmer_span;
+
+ // PARAMETER: Weight of the k-mers
+ Unmodifiable<uint8_t> kmer_weight;
+
+ // VARIABLE: maximum number of consecutive gaps in the gap pattern (will be computed at runtime)
+ Unmodifiable<CountType> max_consecutive_gaps;
+
+ // PARAMETER: Base Call quality cutoff, treat BC with quality < bc_cutoff as miscall
+ Unmodifiable<CountType> min_qual;
+
+ // PARAMETER: max. insert/deletion size
+ Unmodifiable<DiffType> window;
+
+ // PARAMETER: minimum number of errors allowed in alignment
+ Unmodifiable<CountType> min_errors;
+
+ // SWITCH: discard One-hit-wonders
+ Unmodifiable<bool> discard_ohw;
+
+ // PARAMETER: first cycle to discard one-hit-wonders
+ Unmodifiable<CountType> start_ohw;
+
+ // PARAMETER: All-Best-N-Scores-Mode::N
+ Unmodifiable<CountType> best_n;
+
+ // PARAMETER: temporary directory for the streamed alignment
+ Unmodifiable<std::string> temp_dir;
+
+ // SWITCH: write sam/bam output or not
+ Unmodifiable<bool> write_bam;
+
+ // PARAMETER: Cycles for intermediate SAM/BAM output
+ Unmodifiable<std::vector<uint16_t>> output_cycles;
+
+ // SWITCH: Keep the old alignment files of previous cycles
+ Unmodifiable<bool> keep_aln_files;
+
+ // PARAMETER: Memory block size for the input and output buffer in the streamed alignment
+ Unmodifiable<uint64_t> block_size;
+
+ // PARAMETER: Compression format for alignment files
+ Unmodifiable<uint8_t> compression_format;
+
+ // PARAMETER: list of lanes to process
+ Unmodifiable<std::vector<uint16_t>> lanes;
+
+ // PARAMETER: list of tiles to process
+ Unmodifiable<std::vector<uint16_t>> tiles;
+
+ // PARAMETER: root directory of hilive run
+ Unmodifiable<std::string> root;
+
+ // PARAMETER: path to the index file
+ Unmodifiable<std::string> index_fname;
+
+ // PARAMETER: read length of all reads (including barcodes)
+ Unmodifiable<CountType> cycles;
+
+ //PARAMETER: Stores the barcodes defined by the user. The inner vector contains the single fragments of multi-barcodes.
+ Unmodifiable<std::vector<std::vector<std::string>>> barcodeVector;
+
+ // PARAMETER: directory in which to create the output directory structure
+ Unmodifiable<std::string> out_dir;
+
+ // PARAMETER: number of threads to use
+ Unmodifiable<CountType> num_threads;
+
+ // SWITCH: activate extended CIGAR annotation
+ Unmodifiable<bool> extended_cigar;
+
+ /**
+ * Contains the read information of the sequencing machine (as SequenceElement objects). Includes sequence reads and barcodes.
+ * Arbitrary numbers and orders of reads are supported. The summed length of all elements must equal the number of sequencing cycles.
+ * @author Tobias Loka
+ */
+ Unmodifiable<std::vector<SequenceElement>> seqs;
+
+ // Number of mates (information taken from the seqLengths parameter), (Hint: corresponding indeces are 1-based)
+ Unmodifiable<uint16_t> mates;
+
+ // PARAMETER: number of allowed errors for the single barcodes
+ Unmodifiable<std::vector<uint16_t>> barcode_errors;
+
+ // SWITCH: if true, keep all barcodes (disables barcode filtering).
+ Unmodifiable<bool> keep_all_barcodes;
+
+ Unmodifiable<AlignmentMode> mode;
+
+ Unmodifiable<float> min_as_ratio;
+
+ template<typename T>
+ bool set_unmodifiable(Unmodifiable<T> & unmodifiable, T value, std::string variable_name) {
+ try {
+ unmodifiable.set(value);
+ }
+ catch (unmodifiable_error& e) {
+// std::cerr << e.what() << " (" << variable_name << ")." << std::endl;
+ variable_name.length(); // TODO: just to remove warnings. Remove variable_name string when finished.
+ return false ;
+ }
+ return true;
+ }
+
+ template<typename T>
+ T get_unmodifiable(Unmodifiable<T> unmodifiable, std::string variable_name, bool allow_unset = false) {
+ try {
+ return unmodifiable.get(allow_unset);
+ }
+ catch (unmodifiable_error& e) {
+ std::cerr << e.what() << " (" << variable_name << ")." << std::endl;
+ return T();
+ }
+ }
+
+ std::vector<std::string> xmlParse_barcodeVector() {
+ std::vector<std::string> bc_strings;
+ for ( CountType i = 0; i < get_barcodeVector().size(); i++ ) {
+ bc_strings.push_back( get_barcodeString(i) );
+ }
+ return bc_strings;
+ }
+
+ std::vector<std::string> xmlParse_seqs() {
+ std::vector<std::string> seq_vector;
+ for ( auto el : get_seqs() ) {
+ std::string seq_string;
+ seq_string += std::to_string(el.length);
+ seq_string += el.mate == 0 ? "B" : "R";
+ seq_vector.push_back(seq_string);
+ }
+ return seq_vector;
+ }
+
+ void set_barcodeVector(std::vector<std::vector<std::string> > value) {
+ set_unmodifiable(barcodeVector, value, "barcodeVector");
+ }
+
+ void set_seqs(std::vector<SequenceElement> value) {
+ set_unmodifiable(seqs, value, "seqs");
+ }
+
+ void set_block_size(uint64_t value) {
+ set_unmodifiable(block_size, value, "block_size");
+ }
+
+ void set_mates(uint16_t value) {
+ set_unmodifiable(mates, value, "mates");
+ }
+
+ AlignmentMode get_mode() {
+ return get_unmodifiable(mode, "mode");
+ }
+
+ void set_mode(AlignmentMode value, CountType bestn = 0) {
+ set_unmodifiable(mode, value, "mode");
+ set_unmodifiable(best_n, bestn, "best_n");
+ }
+
+public:
+
+ /**
+ * Create a property tree that is filled with all (relevant) settings.
+ * @return Property tree containing all settings.
+ * @author Tobias Loka
+ */
+ boost::property_tree::ptree to_ptree() {
+
+ boost::property_tree::ptree xml_out;
+
+ // General settings
+ xml_out.add_child("settings.lanes", getXMLnode_vector ( get_lanes() ));
+ xml_out.add_child("settings.tiles", getXMLnode_vector ( get_tiles() ));
+ xml_out.add_child("settings.min_errors", getXMLnode (get_min_errors() ));
+ xml_out.add_child("settings.cycles", getXMLnode ( get_cycles() ));
+ xml_out.add_child("settings.sequences", getXMLnode_vector ( xmlParse_seqs() ));
+
+ // Barcode settings
+ xml_out.add_child("settings.barcodes.sequences", getXMLnode_vector( xmlParse_barcodeVector() ));
+ xml_out.add_child("settings.barcodes.errors", getXMLnode_vector ( get_barcode_errors() ));
+ xml_out.add_child("settings.barcodes.keep_all", getXMLnode ( get_keep_all_barcodes() ));
+
+ // Alignment mode
+ std::string mode = std::string(1, char(get_mode()));
+ if ( get_mode() == AlignmentMode::BESTN )
+ mode += std::to_string(get_best_n());
+ xml_out.add_child("settings.mode", getXMLnode ( mode ));
+
+ // Paths
+ xml_out.add_child("settings.paths.temp_dir", getXMLnode ( get_temp_dir() ));
+ xml_out.add_child("settings.paths.out_dir", getXMLnode ( get_out_dir() ));
+ xml_out.add_child("settings.paths.root", getXMLnode ( get_root() ));
+ xml_out.add_child("settings.paths.index", getXMLnode ( get_index_fname() ));
+
+ // Output settings
+ xml_out.add_child("settings.out.bam", getXMLnode ( get_write_bam() ));
+ xml_out.add_child("settings.out.cycles", getXMLnode_vector ( get_output_cycles() ));
+ xml_out.add_child("settings.out.extended_cigar", getXMLnode ( get_extended_cigar() ));
+ xml_out.add_child("settings.out.min_as_ratio", getXMLnode ( get_min_as_ratio()) );
+
+ // Technical settings
+ xml_out.add_child("settings.technical.num_threads", getXMLnode ( get_num_threads() ));
+ xml_out.add_child("settings.technical.keep_aln_files", getXMLnode ( get_keep_aln_files() ));
+ xml_out.add_child("settings.technical.block_size", getXMLnode ( get_block_size() ));
+ xml_out.add_child("settings.technical.compression_format", getXMLnode ( get_compression_format() ));
+
+ // Alignment algorithm settings
+ xml_out.add_child("settings.align.min_qual", getXMLnode (get_min_qual() ));
+ xml_out.add_child("settings.align.window", getXMLnode (get_window() ));
+ xml_out.add_child("settings.align.discard_ohw", getXMLnode ( get_discard_ohw() ));
+ xml_out.add_child("settings.align.start_ohw", getXMLnode ( get_start_ohw() ));
+
+ return xml_out;
+ }
+
+ void set_barcodes ( std::vector< std::string > barcodeArg ) {
+
+ // Get the barcode length(s) from the seqs vector
+ std::vector<uint16_t> barcode_lengths;
+ for ( uint16_t seq_num = 0; seq_num < get_seqs().size(); seq_num++ ) {
+ if ( getSeqById(seq_num).isBarcode() )
+ barcode_lengths.push_back( getSeqById(seq_num).length );
+ }
+
+ // Fill 2D-vector for internal storage
+ std::vector<std::vector<std::string> > barcodeVector;
+ for ( auto barcode = barcodeArg.begin(); barcode != barcodeArg.end(); ++barcode) {
+
+ // Check if all characters in the current barcode are valid
+ std::string valid_chars = seq_chars + "-";
+ for(CountType i = 0; i != (*barcode).length(); i++){
+ char c = (*barcode)[i];
+ if ( valid_chars.find(c) == std::string::npos )
+ throw std::runtime_error("Invalid character '" + std::string(1,c) + "' in barcode sequence " + *barcode + ".");
+ }
+
+ // Split barcode string into fragments
+ std::vector<std::string> fragments;
+ split(*barcode, '-', fragments);
+
+ // Check correct length of all fragments
+ if ( barcode_lengths.size() != fragments.size())
+ throw std::runtime_error("Wrong fragment length in barcode " + *barcode);
+
+ // Check correct number of fragments
+ for ( uint16_t num = 0; num != fragments.size(); num++ ) {
+ if ( fragments[num].length() != barcode_lengths[num] ) {
+ throw std::runtime_error("Wrong number of fragments for barcode " + *barcode);
+ }
+ }
+
+ // Push barcode to the final 2D vector
+ barcodeVector.push_back(fragments);
+ }
+
+ set_barcodeVector(barcodeVector);
+
+ }
+
+ std::vector<std::vector<std::string> > get_barcodeVector() {
+ return get_unmodifiable(barcodeVector, "barcodeVector", true);
+ }
+
+ std::string format_barcode(std::string unformatted_barcode) {
+ CountType pos = 0;
+ for ( auto el : get_seqs() ) {
+ if ( el.mate == 0 ) {
+ pos+=el.length;
+ if ( unformatted_barcode.length() >= pos )
+ unformatted_barcode.insert(pos++, "-");
+ else
+ break;
+ }
+ }
+
+ return unformatted_barcode.substr(0,pos-1);
+
+ }
+
+ std::string get_barcodeString(CountType index) {
+
+ // invalid index
+ if ( index >= get_barcodeVector().size() ) {
+ return "";
+ }
+
+ else {
+
+ std::vector<std::string> bc_vec = get_barcodeVector()[index];
+
+ std::stringstream ss;
+ for ( auto fragment : bc_vec ) {
+ ss << fragment;
+ }
+ std::string barcode_string = ss.str();
+ return format_barcode(barcode_string);
+ }
+ }
+
+ void set_read_structure ( std::vector<std::string> read_argument ) {
+
+ // Init variables
+ CountType lenSum = 0;
+ CountType length = 0;
+ std::string length_string = "";
+ char type;
+ unsigned mates = 0;
+ std::vector<SequenceElement> temp;
+
+ // Iterate through input vector
+ for ( auto read = read_argument.begin(); read != read_argument.end(); ++read ) {
+
+ // Split string into fragment length and fragment type (B or R)
+ length_string = (*read).substr(0,(*read).length()-1);
+ type = (*(*read).rbegin());
+
+ if ( length_string.find_first_not_of("0123456789")!=std::string::npos ) {
+ throw std::runtime_error("Invalid length for read fragment " + *read + ". Please only use unsigned integer values.");
+ }
+
+ try{
+ length = CountType(std::atol(length_string.c_str()));
+ } catch( std::bad_cast & ex ){
+ std::cerr << "Error while casting length " << length_string << " to type uint16_t." << std::endl;
+ throw ex;
+ }
+
+ if ( type!='B' && type!='R' ) {
+ std::stringstream ss;
+ ss << "\'" << type << "\' is no valid read type. Please use \'R\' for sequencing reads or \'B\' for barcode reads.";
+ throw std::runtime_error(ss.str());
+
+ }
+
+ temp.push_back(SequenceElement(temp.size(), (type == 'R') ? ++mates : 0, length));
+ lenSum += length;
+
+ }
+ set_seqs(temp);
+ set_mates(mates);
+
+ if ( lenSum!=get_cycles() ) {
+ throw std::runtime_error("Sum of defined reads does not equal the given number of cycles.");
+ }
+
+ }
+
+ /**
+ * Get a SequenceElement object from the seqs vector by using the id
+ * @param id The id of the SequenceElement.
+ * @return The respective SequenceElement object for the given id.
+ * @author Tobias Loka
+ */
+ SequenceElement getSeqById(CountType id) {return seqs.get()[id];}
+
+ /**
+ * Get a SequenceElement object from the seqs vector by using the mate number
+ * @param id The mate number of the SequenceElement.
+ * @return The respective SequenceElement object for the given mate number. NULLSEQ if mate==0 (barcodes).
+ * @author Tobias Loka
+ */
+ SequenceElement getSeqByMate(CountType mate) {
+ if ( mate == 0 ) return NULLSEQ;
+ auto the_seq = seqs.get();
+ for (uint16_t i = 0; i != the_seq.size(); i++) {
+ if(the_seq[i].mate == mate) return the_seq[i];
+ }
+ return NULLSEQ;
+ }
+
+ std::vector<unsigned> get_kmer_gaps() {
+ return get_unmodifiable(kmer_gaps, "kmer_gaps", true);
+ }
+
+ std::vector<unsigned> get_rev_kmer_gaps() {
+ return get_unmodifiable(rev_kmer_gaps, "rev_kmer_gaps", true);
+ }
+
+ bool set_kmer( uint8_t kmer_weight, std::vector<unsigned> gaps ) {
+
+ if ( gaps.size() > 0 ) {
+ // Prepare user-defined list of gap positions (sort and erase duplicates)
+ std::sort(gaps.begin(), gaps.end());
+ gaps.erase( std::unique(gaps.begin(), gaps.end()), gaps.end());
+
+ // Weight and gap positions not consistent
+ if ( kmer_weight + gaps.size() <= *(std::max_element(gaps.begin(), gaps.end())) || *(std::min_element(gaps.begin(), gaps.end())) <= 1 ) {
+ std::cerr << "Warning: k-mer weight and gap pattern not consistent. Ensure that the first gap positions is >1 and" <<
+ "the maximal gap positions is lower than the total length of the k-mer pattern." << std::endl;
+ return false;
+ }
+ }
+
+ // Set k-mer variables
+ set_unmodifiable(this->kmer_weight, kmer_weight, "kmer_weight");
+ set_unmodifiable(this->kmer_gaps, gaps, "kmer_gaps");
+ set_unmodifiable(this->kmer_span, uint8_t(kmer_weight + gaps.size()), "kmer_span");
+
+ std::vector<unsigned> rev_kmer_gaps;
+ for ( auto gap:gaps ) {
+ rev_kmer_gaps.push_back(this->kmer_span - gap + 1);
+ }
+ std::reverse(rev_kmer_gaps.begin(), rev_kmer_gaps.end());
+ set_unmodifiable(this->rev_kmer_gaps, rev_kmer_gaps, "rev_kmer_gaps");
+
+ // Compute maximal consecutive gaps in gap pattern
+ CountType current_consecutive_gaps = 0;
+ CountType last_gap = 0;
+ CountType temp_max_consecutive_gaps = 0;
+
+ for ( unsigned el : this->get_kmer_gaps() ) {
+
+ // init first gap
+ if ( last_gap == 0 ) {
+ current_consecutive_gaps = 1;
+ last_gap = el;
+ continue;
+ }
+
+ // handle consecutive gaps
+ else if ( el == unsigned( last_gap + 1 ) ){
+ current_consecutive_gaps += 1;
+ last_gap = el;
+ }
+
+ // handle end of gap region
+ else {
+ temp_max_consecutive_gaps = std::max ( temp_max_consecutive_gaps, current_consecutive_gaps );
+ current_consecutive_gaps = 1;
+ last_gap = el;
+ }
+
+ }
+ set_unmodifiable(this->max_consecutive_gaps, std::max ( temp_max_consecutive_gaps, current_consecutive_gaps ), "max_consecutive_gaps");
+
+ return true;
+ }
+
+ uint8_t get_kmer_span() {
+ return get_unmodifiable(kmer_span, "kmer_span");
+ }
+
+ uint8_t get_kmer_weight() {
+ return get_unmodifiable(kmer_weight, "kmer_weight");
+ }
+
+ void set_min_qual(CountType value) {
+ set_unmodifiable(min_qual, value, "min_qual");
+ }
+
+ CountType get_min_qual() {
+ return get_unmodifiable(min_qual, "min_qual");
+ }
+
+ void set_window(DiffType value) {
+ set_unmodifiable(window, value, "window");
+ }
+
+ DiffType get_window() {
+ return get_unmodifiable(window, "window");
+ }
+
+ void set_min_errors(CountType value) {
+ set_unmodifiable(min_errors, value, "min_errors");
+ }
+
+ CountType get_min_errors() {
+ return get_unmodifiable(min_errors, "min_errors");
+ }
+
+ void disable_ohw(bool value) {
+ set_unmodifiable(discard_ohw, !value, "discard_ohw");
+ }
+
+ bool get_discard_ohw() {
+ return get_unmodifiable(discard_ohw, "discard_ohw");
+ }
+
+ void set_start_ohw(CountType value) {
+ set_unmodifiable(start_ohw, value, "start_ohw");
+ }
+
+ CountType get_start_ohw() {
+ return get_unmodifiable(start_ohw, "start_ohw");
+ }
+
+ bool get_any_best_hit_mode() {
+ return (get_unmodifiable(mode, "mode")==AlignmentMode::ANYBEST);
+ }
+
+ bool get_all_hit_mode() {
+ return (get_mode()==AlignmentMode::ALL);
+ }
+
+ bool get_all_best_hit_mode() {
+ return (get_mode()==AlignmentMode::ALLBEST);
+ }
+
+ bool get_all_best_n_scores_mode() {
+ return (get_mode()==AlignmentMode::BESTN);
+ }
+
+ CountType get_best_n() {
+ return get_unmodifiable(best_n, "best_n", true);
+ }
+
+ void set_temp_dir(std::string value) {
+ set_unmodifiable(temp_dir, value, "temp_dir");
+ }
+
+ std::string get_temp_dir() {
+ return get_unmodifiable(temp_dir, "temp_dir");
+ }
+
+ void set_write_bam(bool value) {
+ set_unmodifiable(write_bam, value, "write_bam");
+ }
+
+ bool get_write_bam() {
+ return get_unmodifiable(write_bam, "write_bam");
+ }
+
+ void set_output_cycles(std::vector<uint16_t> cycles) {
+ std::vector<uint16_t> the_cycles;
+ for ( auto it = cycles.begin(); it != cycles.end(); ++it ) {
+ if ( *it > get_cycles() )
+ the_cycles.push_back(get_cycles());
+ else
+ the_cycles.push_back(*it);
+ }
+ std::sort(the_cycles.begin(), the_cycles.end());
+ the_cycles.erase( std::unique(the_cycles.begin(), the_cycles.end()), the_cycles.end());
+ set_unmodifiable(output_cycles, the_cycles, "output_cycles");
+ }
+
+ std::vector<uint16_t> get_output_cycles() {
+ return get_unmodifiable(output_cycles, "output_cycles", true);
+ }
+
+ bool is_output_cycle(CountType cycle) {
+ auto out_cycles = get_output_cycles();
+ if ( std::find(out_cycles.begin(), out_cycles.end(), cycle) == out_cycles.end() )
+ return false;
+ return true;
+ }
+
+ void set_keep_aln_files(bool value) {
+ set_unmodifiable(keep_aln_files, value, "keep_aln_files");
+ }
+
+ bool get_keep_aln_files() {
+ return get_unmodifiable(keep_aln_files, "keep_aln_files");
+ }
+
+ void set_block_size(std::string value) {
+
+ uint64_t size;
+ char type = 'B';
+
+ // Split value to size and type
+ if ( value.find_first_of("BKM") != std::string::npos ) {
+ type = *value.rbegin();
+ value = value.substr(0,value.length()-1);
+ }
+
+ if ( value.find_first_not_of("0123456789")!=std::string::npos ) {
+ throw std::runtime_error("Invalid block size " + value + ". Please only use unsigned integer values.");
+ }
+
+ try{
+ size = uint64_t(std::atol(value.c_str()));
+ } catch( std::bad_cast & ex ){
+ std::cerr << "Error while casting length " << value << " to type uint16_t." << std::endl;
+ throw ex;
+ }
+
+ if ( type == 'B' )
+ set_block_size(size);
+ else if ( type == 'K' )
+ set_block_size(size*1024);
+ else if ( type == 'M' )
+ set_block_size(size*1024*1024);
+ else
+ throw std::runtime_error("Invalid block size type. Only 'B' (Bytes), 'K' (Kilobytes) or 'M' (Megabytes) are permitted.");
+
+ }
+
+ uint64_t get_block_size() {
+ return get_unmodifiable(block_size, "block_size");
+ }
+
+ void set_compression_format(uint16_t value) {
+ if ( value > 2 )
+ value = 2;
+ uint8_t one_byte_value = value;
+ set_unmodifiable(compression_format, one_byte_value, "compression_format");
+ }
+
+ uint8_t get_compression_format() {
+ return get_unmodifiable(compression_format, "compression_format");
+ }
+
+ void set_lanes(std::vector<uint16_t> value) {
+ std::sort( value.begin(), value.end() );
+ value.erase( std::unique( value.begin(), value.end() ), value.end() );
+ set_unmodifiable(lanes, value, "lanes");
+ }
+
+ std::vector<uint16_t> get_lanes() {
+ return get_unmodifiable(lanes, "lanes", true);
+ }
+
+ void set_mode(std::string value) {
+
+ // All hit mode
+ if ( value == "ALL" || value == "A" ) {
+ set_mode(AlignmentMode::ALL);
+ }
+
+ // Best N scores mode
+ else if ( value.substr(0,5) == "BESTN" || value.substr(0,1) == "N" ) {
+
+ std::string bestn = value.substr(0,5) == "BESTN" ? value.substr(5) : value.substr(1);
+
+ if ( bestn.find_first_not_of("0123456789")!=std::string::npos ) {
+ throw std::runtime_error("Invalid alignment mode: " + value + ".");
+ }
+ try{
+ set_mode(AlignmentMode::BESTN, CountType(std::atol(bestn.c_str())));
+ } catch( std::bad_cast & ex ){
+ std::cerr << "Error while casting length " << bestn << " to type uint16_t." << std::endl;
+ throw ex;
+ }
+ }
+
+ // All best mode
+ else if ( value == "ALLBEST" || value == "H" ) {
+ set_mode(AlignmentMode::ALLBEST);
+ }
+
+ // All hit mode
+ else if ( value == "ANYBEST" || value == "B" ) {
+ set_mode(AlignmentMode::ANYBEST);
+ }
+
+ // Unknown mode
+ else {
+ throw std::runtime_error("Invalid alignment mode: " + value + ".");
+ }
+ }
+
+ void set_tiles(std::vector<uint16_t> value) {
+ std::sort( value.begin(), value.end() );
+ value.erase( std::unique( value.begin(), value.end() ), value.end() );
+ set_unmodifiable(tiles, value, "tiles");
+ }
+
+ std::vector<uint16_t> get_tiles() {
+ return get_unmodifiable(tiles, "tiles", true);
+ }
+
+ void set_root(std::string value) {
+ set_unmodifiable(root, value, "root");
+ }
+
+ std::string get_root() {
+ return get_unmodifiable(root, "root");
+ }
+
+ void set_index_fname(std::string value) {
+ set_unmodifiable(index_fname, value, "index_fname");
+ }
+
+ std::string get_index_fname() {
+ return get_unmodifiable(index_fname, "index_fname");
+ }
+
+ void set_cycles(CountType value) {
+ set_unmodifiable(cycles, value, "cycles");
+ }
+
+ CountType get_cycles() {
+ return get_unmodifiable(cycles, "cycles");
+ }
+
+ void set_out_dir(std::string value) {
+ set_unmodifiable(out_dir, value, "out_dir");
+ }
+
+ std::string get_out_dir() {
+ return get_unmodifiable(out_dir, "out_dir");
+ }
+
+ void set_num_threads(CountType value) {
+ set_unmodifiable(num_threads, value, "num_threads");
+ }
+
+ CountType get_num_threads() {
+ return get_unmodifiable(num_threads, "num_threads");
+ }
+
+ std::vector<SequenceElement> get_seqs() {
+ return get_unmodifiable(seqs, "seqs", true);
+ }
+
+ uint16_t get_mates() {
+ return get_unmodifiable(mates, "mates", true);
+ }
+
+ void set_barcode_errors(std::vector<uint16_t> value) {
+ set_unmodifiable(barcode_errors, value, "barcode_errors");
+ }
+
+ std::vector<uint16_t> get_barcode_errors() {
+ return get_unmodifiable(barcode_errors, "barcode_errors", true);
+ }
+
+ void set_keep_all_barcodes(bool value) {
+ set_unmodifiable(keep_all_barcodes, value, "keep_all_barcodes");
+ }
+
+ bool get_keep_all_barcodes() {
+ if ( get_barcodeVector().size() == 0 )
+ return true;
+ return get_unmodifiable(keep_all_barcodes, "keep_all_barcodes");
+ }
+
+ void set_extended_cigar(bool value) {
+ set_unmodifiable(extended_cigar, value, "extended_cigar");
+ }
+
+ bool get_extended_cigar() {
+ return get_unmodifiable(extended_cigar, "extended_cigar");
+ }
+
+ CountType get_max_consecutive_gaps() {
+ return get_unmodifiable(max_consecutive_gaps, "max_consecutive_gaps");
+ }
+
+ float get_min_as_ratio() {
+ return get_unmodifiable(min_as_ratio, "min_as_ratio");
+ }
+
+ void set_min_as_ratio(float value) {
+ if ( value > 1.0f )
+ value = 1.0f;
+ if ( value < 0.0f )
+ value = 0.0f;
+ set_unmodifiable(min_as_ratio, value, "min_as_ratio");
+ }
+
+};
+
+#endif
diff --git a/lib/alnread.cpp b/lib/alnread.cpp
index 06ebea8..72f096b 100644
--- a/lib/alnread.cpp
+++ b/lib/alnread.cpp
@@ -1,53 +1,78 @@
#include "alnread.h"
-seqan::String<seqan::CigarElement<> > Seed::returnSeqanCigarString() {
+seqan::String<seqan::CigarElement<> > Seed::returnSeqanCigarString(unsigned* nm_i) {
typedef seqan::String<seqan::CigarElement<> > TSeqanCigarString;
TSeqanCigarString seqanCigarString;
seqan::CigarElement<> cigarElem;
int last_offset = 0;
- for (CigarVector::const_iterator it = cigar_data.begin(); it != cigar_data.end(); ++it) { // for every element in cigar_data
- if (it == cigar_data.begin() && (*it).offset==NO_MATCH) { // Alignment begins with NO_MATCH region => Softclipped start
+ for (CigarVector::const_iterator it = cigar_data.begin(); it != cigar_data.end(); ++it) {
+
+ // Alignment begins with NO_MATCH region => Softclipped start
+ if (it == cigar_data.begin() && (*it).offset==NO_MATCH ) {
cigarElem.operation='S';
cigarElem.count=(*it).length;
seqan::appendValue(seqanCigarString, cigarElem);
continue;
}
- if (it == --cigar_data.end() && (*it).offset==NO_MATCH) { // Alignment ends with NO_MATCH region => Softclipped end
+
+ // Alignment ends with NO_MATCH region => Softclipped end
+ if (it == --cigar_data.end() && (*it).offset==NO_MATCH ) {
cigarElem.operation='S';
cigarElem.count=(*it).length;
seqan::appendValue(seqanCigarString, cigarElem);
continue;
}
- if ((*it).offset==NO_MATCH) { // ??? -> Mismatch
- cigarElem.operation='M';
+
+ // Alignment ends with NO_MATCH + TRIMMED_MATCH region => Softclipped end
+ if ( it == --(--cigar_data.end()) && (*it).offset==NO_MATCH && (--cigar_data.end())->offset==TRIMMED_MATCH ) {
+ cigarElem.operation='S';
+ cigarElem.count=( (*it).length + (*(++it)).length );
+ seqan::appendValue(seqanCigarString, cigarElem);
+ continue;
+ }
+
+ // Mismatch region
+ if ((*it).offset==NO_MATCH) {
+ cigarElem.operation = globalAlignmentSettings.get_extended_cigar() ? 'X' : 'M';
cigarElem.count=(*it).length;
+ (*nm_i) += (*it).length;
seqan::appendValue(seqanCigarString, cigarElem);
continue;
}
- if ((*it).offset!=NO_MATCH) { // ??? -> Match
- if (last_offset == (*it).offset) { // no offset change
- cigarElem.operation='M';
+
+ // Match region
+ if ((*it).offset!=NO_MATCH) {
+
+ // no offset change
+ if (last_offset == (*it).offset) {
+ cigarElem.operation = globalAlignmentSettings.get_extended_cigar() ? '=' : 'M';
cigarElem.count=(*it).length;
seqan::appendValue(seqanCigarString, cigarElem);
last_offset = (*it).offset;
continue;
}
- if (last_offset < (*it).offset) { // offset gets bigger => reference in alignment is longer than read => deletion in read (and thereby cigar string)
+
+ // offset gets bigger => reference in alignment is longer than read => deletion in read (and thereby cigar string)
+ if (last_offset < (*it).offset) {
cigarElem.operation='D';
cigarElem.count=(*it).offset - last_offset;
+ (*nm_i) += (*it).offset - last_offset;
seqan::appendValue(seqanCigarString, cigarElem);
- cigarElem.operation='M';
+ cigarElem.operation = globalAlignmentSettings.get_extended_cigar() ? '=' : 'M';
cigarElem.count=(*it).length;
seqan::appendValue(seqanCigarString, cigarElem);
last_offset = (*it).offset;
continue;
}
- if (last_offset > (*it).offset) { // offset gets smaller => reference in alignment is smaller than read => insertion in read (and thereby cigar string)
+
+ // offset gets smaller => reference in alignment is smaller than read => insertion in read (and thereby cigar string)
+ if (last_offset > (*it).offset) {
cigarElem.operation='I';
cigarElem.count=last_offset - (*it).offset;
+ (*nm_i) += last_offset - (*it).offset;
seqan::appendValue(seqanCigarString, cigarElem);
- cigarElem.operation='M';
+ cigarElem.operation = globalAlignmentSettings.get_extended_cigar() ? '=' : 'M';
cigarElem.count=(*it).length;
seqan::appendValue(seqanCigarString, cigarElem);
last_offset = (*it).offset;
@@ -56,9 +81,9 @@ seqan::String<seqan::CigarElement<> > Seed::returnSeqanCigarString() {
}
}
- // collapse Neighboring match regions
+ // collapse Neighboring identical regions
for (unsigned k = 1; k<length(seqanCigarString); k++)
- if ((seqanCigarString[k-1].operation == 'M') && (seqanCigarString[k].operation == 'M')) {
+ if ( seqanCigarString[k-1].operation == seqanCigarString[k].operation ) {
unsigned temp = seqanCigarString[k].count;
erase(seqanCigarString, k);
seqanCigarString[k-1].count += temp;
@@ -71,6 +96,25 @@ seqan::String<seqan::CigarElement<> > Seed::returnSeqanCigarString() {
return seqanCigarString;
}
+void Seed::cout(){
+ std::cout << "----- SEED START -----" << std::endl;
+ std::cout << "gid: " << this->gid << std::endl;
+ std::cout << "start_pos: " << this->start_pos << std::endl;
+ std::cout << "num_matches: " << this->num_matches << std::endl;
+ std::cout << "CIGAR: ";
+ for ( auto el : this->cigar_data ) {
+ std::cout << el.length;
+ if ( el.offset == NO_MATCH ) {
+ std::cout << "X ";
+ } else if (el.offset == TRIMMED_MATCH ) {
+ std::cout << "T";
+ } else {
+ std::cout << "M(" << el.offset << ") ";
+ }
+ }
+ std::cout << std::endl << "------ SEED END ------" << std::endl;
+}
+
uint16_t Seed::serialize_size() {
// calculate total size
uint16_t total_size = 0;
@@ -165,8 +209,8 @@ uint16_t Seed::deserialize(char* d) {
}
-void ReadAlignment::set_rlen(CountType r) {
- rlen = r;
+void ReadAlignment::set_total_cycles(CountType c) {
+ total_cycles = c;
}
@@ -179,6 +223,8 @@ uint64_t ReadAlignment::serialize_size() {
total_size += sizeof(CountType); // the sequence length
total_size += sequenceStoreVector.size()*(sizeof(uint8_t)); // the sequence information itself
+ total_size += sizeof(CountType); // The barcode length
+ total_size += barcodeStoreVector.size()*(sizeof(uint8_t)); // the barcode sequence information
// total number of seeds
total_size += sizeof(uint32_t);
@@ -223,6 +269,16 @@ std::vector<char> ReadAlignment::serialize() {
d += sizeof(uint8_t);
}
+ // write the barcode length
+ memcpy(d,&barcodeLen,sizeof(CountType));
+ d += sizeof(CountType);
+
+ // write the barcodeStoreVector
+ for (auto it = barcodeStoreVector.begin(); it != barcodeStoreVector.end(); ++it) {
+ memcpy(d,&(*it),sizeof(uint8_t));
+ d += sizeof(uint8_t);
+ }
+
// write the number of seeds
memcpy(d,&num_seeds,sizeof(uint32_t));
d += sizeof(uint32_t);
@@ -266,7 +322,7 @@ uint64_t ReadAlignment::deserialize(char* d) {
bytes += sizeof(CountType);
// read the sequence
- unsigned seqVec_size = (unsigned) std::ceil((float) sequenceLen / 4.0);
+ unsigned seqVec_size = (unsigned) std::ceil((float) sequenceLen / 2.0);
sequenceStoreVector.clear();
sequenceStoreVector.reserve(seqVec_size);
for (unsigned i = 0; i <seqVec_size; ++i) {
@@ -276,6 +332,22 @@ uint64_t ReadAlignment::deserialize(char* d) {
sequenceStoreVector.push_back(elem);
}
+ // read the barcode length
+ barcodeLen = 0;
+ memcpy(&barcodeLen,d+bytes,sizeof(CountType));
+ bytes += sizeof(CountType);
+
+ // read the barcode
+ unsigned barVec_size = (unsigned) std::ceil((float) barcodeLen / 2.0);
+ barcodeStoreVector.clear();
+ barcodeStoreVector.reserve(barVec_size);
+ for (unsigned i = 0; i <barVec_size; ++i) {
+ uint8_t elem;
+ memcpy(&(elem),d+bytes,sizeof(uint8_t));
+ bytes += sizeof(uint8_t);
+ barcodeStoreVector.push_back(elem);
+ }
+
// read the number of seeds
uint32_t num_seeds = 0;
memcpy(&num_seeds,d+bytes,sizeof(uint32_t));
@@ -305,72 +377,116 @@ uint64_t ReadAlignment::deserialize(char* d) {
// convert and return sequence of the seed as string
-std::string ReadAlignment::getSequenceString(AlignmentSettings & settings) {
- std::string seq = "";
- // append one 4 base block at a time
- for (unsigned i = 0; i<sequenceStoreVector.size(); i++)
- seq.append(unhash(sequenceStoreVector[i], 4)); // 4 because 4 nucleotides fit in one element of sequenceStoreVector, namely a uint8_t
-
- // delete last few bases, because sequence is only len long
- for (unsigned i = sequenceLen; i<4*sequenceStoreVector.size(); ++i)
- seq.pop_back();
-
- // return sequence without barcode
- return seq.substr(0,settings.seqlen);
+std::string ReadAlignment::getSequenceString() {
+
+ std::string seq = "";
+ uint8_t two_bit_mask = 3;
+ uint8_t four_bit_mask = 15;
+
+ // iterate through all sequence bytes
+ for (unsigned i = 0; i<sequenceLen; i++ ) {
+
+ // Four-bit representation of the next base call (2-bit qual + 2-bit nucl)
+ uint8_t next = ( sequenceStoreVector[i/2] >> ( (i%2 == 0) * 4) ) & four_bit_mask;
+
+ if ( next < 4 ) { // two-bit qual == 0 --> N-call
+ seq.append("N");
+ } else { // two-bit qual > 0 --> write nucleotide
+ seq += revtwobit_repr(next & two_bit_mask);
+ }
+ }
+
+ // return barcode sequence
+ return seq;
}
-// convert and return barcode
-std::string ReadAlignment::getBarcodeString(AlignmentSettings & settings) {
- std::string seq = "";
- // append one 4 base block at a time
- for (unsigned i = 0; i<sequenceStoreVector.size(); i++)
- seq.append(unhash(sequenceStoreVector[i], 4)); // 4 because 4 nucleotides fit in one element of sequenceStoreVector, namely a uint8_t
+std::string ReadAlignment::getBarcodeString() {
+
+ std::string seq = "";
+ uint8_t two_bit_mask = 3;
+ uint8_t four_bit_mask = 15;
+
+ // iterate through all sequence bytes
+ for (unsigned i = 0; i<barcodeLen; i++ ) {
+
+ // Four-bit representation of the next base call (2-bit qual + 2-bit nucl)
+ uint8_t next = ( barcodeStoreVector[i/2] >> ( (i%2 == 0) * 4) ) & four_bit_mask;
- // delete last few bases, because sequence is only len long
- for (unsigned i = sequenceLen; i<4*sequenceStoreVector.size(); ++i)
- seq.pop_back();
+ if ( next < 4 ) { // two-bit qual == 0 --> N-call
+ seq.append("N");
+ } else { // two-bit qual > 0 --> write nucleotide
+ seq += revtwobit_repr(next & two_bit_mask);
+ }
+ }
- // return only barcode of sequence
- return seq.substr(settings.seqlen);
+ // return barcode sequence
+ return seq;
}
// append one nucleotide to sequenceStoreVector
-void ReadAlignment::appendNucleotideToSequenceStoreVector(char nuc) {
+void ReadAlignment::appendNucleotideToSequenceStoreVector(char bc, bool appendToBarCode) {
+
+ uint8_t nucl = bc & 3;
+ uint8_t qual = bc >> 2;
+
+ // Convert nucl and qual to four-bit value (first two bits describing the quality (0=N; 1=invalid; 2=valid); second two bits describing the nucleotide)
+ uint8_t four_bit_repr = ( ((qual != 0) + (qual > globalAlignmentSettings.get_min_qual())) << 2 ) | nucl;
+
+ CountType & len = appendToBarCode ? barcodeLen : sequenceLen;
+ std::vector<uint8_t> & seqVector = appendToBarCode ? barcodeStoreVector : sequenceStoreVector;
+
// check if all bits from sequenceStoreVector are used
- if (sequenceLen % 4 == 0) { // yes, all used => new 8 Bit block needs to be created
- uint8_t newBlock = twobit_repr(nuc);
- newBlock = newBlock << 6; // 'empty' bits are on the right side
- sequenceStoreVector.push_back(newBlock);
- ++sequenceLen;
+ if (len % 2 == 0) { // yes, all used => new 8 Bit block needs to be created
+ uint8_t newBlock = four_bit_repr << 4; // 'empty' bits are on the right side
+ seqVector.push_back(newBlock);
+ ++len;
}
- else { // not all bits are used. At least the last two are unoccupied
- // append new 2 bits to the right of the old bits
- if (sequenceLen % 4 == 1) // 6 bits empty
- sequenceStoreVector.back() = sequenceStoreVector.back() >> 4;
- if (sequenceLen % 4 == 2) // 4 bits empty
- sequenceStoreVector.back() = sequenceStoreVector.back() >> 2;
- sequenceStoreVector.back() = sequenceStoreVector.back() | twobit_repr(nuc);
- ++sequenceLen;
- // shift used bits until they are left aligned
- if (sequenceLen % 4 == 2) // 4 bits empty
- sequenceStoreVector.back() = sequenceStoreVector.back() << 4;
- if (sequenceLen % 4 == 3) // 2 bits empty
- sequenceStoreVector.back() = sequenceStoreVector.back() << 2;
+
+ else { // not all bits are used. There is enough space for the new basecall
+ seqVector.back() = seqVector.back() | four_bit_repr;
+ ++len;
}
}
// helper function for add_new_seeds
bool seed_compare_pos (const USeed & i, const USeed & j) {
- return (i->start_pos < j->start_pos);
+ if ( i->start_pos == j->start_pos )
+ return i->gid < j->gid;
+ return (i->start_pos < j->start_pos);
}
// Create new seeds from a list of kmer positions and add to current seeds
-void ReadAlignment::add_new_seeds(GenomePosListType& pos, std::vector<bool> & posWasUsedForExtension, AlignmentSettings & settings) {
- SeedVecIt sit = seeds.begin();
+void ReadAlignment::add_new_seeds(GenomePosListType& pos, std::vector<bool> & posWasUsedForExtension) {
+
+ SeedVecIt sit = seeds.begin();
+
+ CigarVector front;
+ CountType num_matches_placeholder = 0;
+
+ // If PLACEHOLDER exist, start with its CIGAR vector
+ if ( seeds.size() > 0 && (*sit)->gid == TRIMMED ) {
+
+ front = (*sit)->cigar_data;
+ num_matches_placeholder = (*sit)->num_matches;
+
+ seeds.pop_front();
+ sit = seeds.begin();
+
+ }
+
+ // If no PLACEHOLDER exist, create initial CIGAR vector
+ else {
+
+ if ( cycle > globalAlignmentSettings.get_kmer_span() ) {
+ front.emplace_back(cycle-globalAlignmentSettings.get_kmer_span(), NO_MATCH);
+ }
+ front.emplace_back(0,0);
+
+ }
for(GenomePosListIt it = pos.begin(); it != pos.end(); ++it) {
if (posWasUsedForExtension[it - pos.begin()]) // if current reference hit was used at least once for seed extension
@@ -378,19 +494,22 @@ void ReadAlignment::add_new_seeds(GenomePosListType& pos, std::vector<bool> & po
USeed s (new Seed);
s->gid = it->gid;
- s->start_pos = it->pos - (cycle-settings.kmer_span);
- s->num_matches = K_HiLive;
- s->cigar_data.clear();
- if (cycle-settings.kmer_span > 0)
- s->cigar_data.emplace_back(cycle-settings.kmer_span,NO_MATCH);
+ s->start_pos = it->pos - (cycle-globalAlignmentSettings.get_kmer_span());
+ s->num_matches = globalAlignmentSettings.get_kmer_weight() + num_matches_placeholder;
+ s->cigar_data = front;
// set correct matches and mismatches depending on kmer mask
- std::vector<unsigned> gapVec = settings.kmer_gaps;
- gapVec.push_back(settings.kmer_span+1);
+ std::vector<unsigned> gapVec = s->start_pos > 0 ? globalAlignmentSettings.get_kmer_gaps() : globalAlignmentSettings.get_rev_kmer_gaps();
+
+ gapVec.push_back(globalAlignmentSettings.get_kmer_span()+1);
unsigned lastProcessedGapPosition = 0;
for (unsigned gapIndex = 0, nextGap; gapIndex < gapVec.size(); ++gapIndex) {
nextGap = gapVec[gapIndex];
- if (nextGap - lastProcessedGapPosition - 1 > 0)
+
+ // Join first kmer match region with the existing match region at the end of the CIGAR string
+ if ( lastProcessedGapPosition == 0 )
+ s->cigar_data.back().length += (nextGap - lastProcessedGapPosition - 1);
+ else if (nextGap - lastProcessedGapPosition - 1 > 0)
s->cigar_data.emplace_back(nextGap - lastProcessedGapPosition - 1,0);
lastProcessedGapPosition = nextGap;
s->cigar_data.emplace_back(1,NO_MATCH);
@@ -398,6 +517,7 @@ void ReadAlignment::add_new_seeds(GenomePosListType& pos, std::vector<bool> & po
s->cigar_data.pop_back(); // remove tailing NO_MATCH from the for-loop
// insert seed into sorted list of seeds
+ // PLACEHOLDER does not have to be considered here because it was converted during seed creation
while (sit != seeds.end() && seed_compare_pos(*sit, s)) // if seed exists and elem has larger starting position than (*sit)
++sit;
sit = seeds.insert(sit, std::move(s));
@@ -406,113 +526,223 @@ void ReadAlignment::add_new_seeds(GenomePosListType& pos, std::vector<bool> & po
// Extend or create a placeholder seed for read with only trimmed matches
-void ReadAlignment::create_placeholder_seed(AlignmentSettings & settings) {
- // if no seed exist, create one (only in the correct cycle and if parameter is true).
- if(seeds.size()==0){
- USeed s (new Seed);
- s->gid = TRIMMED;
- s->num_matches = 1;
- s->cigar_data.clear();
- s->cigar_data.emplace_back(settings.kmer_span,TRIMMED_MATCH);
-
- // this is always sorted, because there is only one seed now
- seeds.push_back(std::move(s));
+void ReadAlignment::create_placeholder_seed() {
+
+ // Don't create PLACEHOLDER if already exist
+ if ( minErrors_in_region( cycle - globalAlignmentSettings.get_kmer_span(), 1) > globalAlignmentSettings.get_min_errors() ) {
+ return;
+ }
+
+ // Don't create PLACEHOLDER if already existing
+ if ( seeds.size() > 0 && (*seeds.begin())->gid == TRIMMED ) {
+ return;
}
+
+ USeed s (new Seed);
+ s->gid = TRIMMED;
+ s->num_matches = 1;
+ s->cigar_data.clear();
+ if ( cycle > globalAlignmentSettings.get_kmer_span() )
+ s->cigar_data.emplace_back(cycle - globalAlignmentSettings.get_kmer_span(), NO_MATCH);
+ s->cigar_data.emplace_back(1,0);
+
+ // Put PLACEHOLDER to the first position of the vector
+ seeds.push_front(std::move(s));
+
}
+CountType ReadAlignment::minErrors_in_region(CountType region_length, CountType border, CountType offset_change) {
-// convert a placeholder seed to a set of normal seeds
-void ReadAlignment::convertPlaceholder(GenomePosListType& pos, AlignmentSettings & settings){
- if(seeds.size()!=1)
- return;
- auto s = seeds.begin();
- auto matches = (*s)->num_matches;
- if(cycle - settings.kmer_span > 0 && (*s)->gid==TRIMMED){
- seeds.clear();
- for(auto p = pos.begin(); p != pos.end(); ++p){
- USeed newSeed (new Seed());
- newSeed->gid = p->gid;
- newSeed->start_pos = p->pos - (cycle - settings.kmer_span);
- newSeed->num_matches = matches;
- newSeed->cigar_data.clear();
- newSeed->cigar_data.emplace_back(cycle - settings.kmer_span,NO_MATCH);
-
- // set correct matches and mismatches depending on kmer mask
- std::vector<unsigned> gapVec = settings.kmer_gaps;
- gapVec.push_back(settings.kmer_span); // without +1 because the last match will be inserted in extend_alignment
- unsigned lastProcessedGapPosition = 0;
- for (unsigned gapIndex = 0, nextGap; gapIndex < gapVec.size(); ++gapIndex) {
- nextGap = gapVec[gapIndex];
- if (nextGap - lastProcessedGapPosition - 1 > 0)
- newSeed->cigar_data.emplace_back(nextGap - lastProcessedGapPosition - 1,0);
- lastProcessedGapPosition = nextGap;
- newSeed->cigar_data.emplace_back(1,NO_MATCH);
- }
- newSeed->cigar_data.pop_back(); // remove tailing NO_MATCH from the for-loop
+ // Border must not be larger than 2
+ if ( border > 2 )
+ return 0;
- // insert into sorted list. The list is empty and pos vector is already sorted.
- // therefore I only push back
- seeds.push_back(std::move(newSeed));
- }
+ // Estimate no errors for regions of length 0
+ if ( region_length == 0 ) {
+ return offset_change;
}
+
+ // Magic formula to estimate the minimal number of matches (supports gapped/spaced kmers)
+ int minErr = ( 1 - border ) + ( ( region_length + border + globalAlignmentSettings.get_kmer_span() + globalAlignmentSettings.get_max_consecutive_gaps() - 2 ) / (globalAlignmentSettings.get_kmer_span() + globalAlignmentSettings.get_max_consecutive_gaps()) );
+ minErr = std::max( minErr, int(offset_change) );
+
+ // Catch negative values
+ if ( minErr < 0 )
+ return 0;
+
+ return CountType(minErr);
}
-// filter seeds based on filtering mode and q gram lemma. Also calls add_new_seeds.
-void ReadAlignment::filterAndCreateNewSeeds(AlignmentSettings & settings, GenomePosListType & pos, std::vector<bool> & posWasUsedForExtension) {
- // TODO this function misbehaves when using demultiplexing. Reads might be preferred, kicking others only to get discarded due to barcode in the end.
- // compute possible remaining matches
- int possibleRemainingMatches = settings.seqlen - cycle + settings.kmer_span - 1;
- if (cycle == settings.seqlen)
- possibleRemainingMatches = 0;
-
- // compute threshold for number of matching bases a seed must already have to have a chance of getting printed in the end
- // adjust threshold via adapted q-gram-lemma
- int num_matches_threshold = settings.seqlen - settings.min_errors*(settings.kmer_span) - possibleRemainingMatches;
-
- if ((settings.all_best_hit_mode) || (settings.any_best_hit_mode)) {
- int max_num_matches = 0;
- for(SeedVecIt sd = seeds.begin() ; sd !=seeds.end(); ++sd)
- max_num_matches = std::max(max_num_matches, (int) (*sd)->num_matches);
- num_matches_threshold = std::max(num_matches_threshold, (int) max_num_matches - possibleRemainingMatches);
- }
- if (settings.all_best_n_scores_mode && settings.best_n > 0) {
- std::vector<CountType> num_matches_vector;
- for (SeedVecIt it=seeds.begin(); it!=seeds.end();++it)
- num_matches_vector.push_back((*it)->num_matches);
- std::sort(num_matches_vector.begin(), num_matches_vector.end(), std::greater<CountType>()); // sort in decreasing order
- unsigned numberOfUniques = std::unique(num_matches_vector.begin(), num_matches_vector.end()) - num_matches_vector.begin(); // uniques are in front of vector now
- if (numberOfUniques > settings.best_n) {
- CountType nth_best_num_matches = num_matches_vector[settings.best_n - 1];
- num_matches_threshold = std::max(num_matches_threshold, (int) nth_best_num_matches - possibleRemainingMatches);
+CountType ReadAlignment::min_errors(USeed & s) {
+
+ CigarVector* c = &(s->cigar_data);
+
+ // Catch elements with length 1 beforehand to save runtime. There can't be any errors in this case.
+ if ( (*c).size() <= 1 )
+ return 0;
+
+ CountType minErr = 0;
+ CountType border = 0;
+ CountType region_length = 0;
+
+ DiffType last_offset = 0;
+
+ // Iterate through all CIGAR elements
+ for ( auto cig_el = (*c).begin(); cig_el != (*c).end(); ++cig_el ) {
+
+ if ( cig_el == (--(*c).end()) && cig_el->offset == TRIMMED_MATCH ) {
+ border += 1;
+ continue;
+ }
+
+ // Finish error region if CIGAR MATCH element spans a complete k-mer
+ if ( cig_el->offset != NO_MATCH && cig_el->length >= ( globalAlignmentSettings.get_kmer_span() -1 ) ) {
+ CountType offset_change = ( cig_el->offset > last_offset ) ? ( cig_el->offset - last_offset ) : ( last_offset - cig_el->offset );
+ minErr += minErrors_in_region(region_length, border, offset_change);
+ region_length = 0;
+ border = 0;
+ last_offset = cig_el->offset;
+ continue;
+ }
+
+ // Init or continue error region for NO_MATCH and too short MATCH elements
+ else {
+ region_length += cig_el->length;
+ border += ( cig_el == (*c).begin() );
+ border += ( cig_el == ( --( (*c).end() ) ) );
+ }
}
- }
- // else nothing has to be done for all-hit-mode
-
+
+ DiffType final_offset = NO_MATCH;
+
+ // Compute offset of the last match CIGAR element
+ for ( auto cig_el = (*c).rbegin(); final_offset == NO_MATCH || final_offset == TRIMMED_MATCH; ++cig_el ) {
+ final_offset = (*cig_el).offset;
+ }
+ CountType offset_change = ( final_offset > last_offset ) ? ( final_offset - last_offset ) : ( last_offset - final_offset );
+
+ // Finish last region
+ minErr += minErrors_in_region(region_length, border, offset_change);
+
+ return minErr;
+}
+
+
+// filter seeds based on filtering mode and q gram lemma. Also calls add_new_seeds.
+void ReadAlignment::filterAndCreateNewSeeds(GenomePosListType & pos, std::vector<bool> & posWasUsedForExtension) {
+
+ // Compute the number of maximum estimated number of errors for the remaining cycles
+ CountType possible_remaining_errors = minErrors_in_region( total_cycles - cycle, 1);
+
+ CountType min_num_errors = globalAlignmentSettings.get_min_errors();
+ CountType max_num_matches = 0; // only required for any best mode in last cycle
+
+ // Compute the number of errors to remove a seed in any_best and all_best mode
+ if ( globalAlignmentSettings.get_all_best_hit_mode() || globalAlignmentSettings.get_any_best_hit_mode() ) {
+ for(SeedVecIt sd = seeds.begin() ; sd !=seeds.end(); ++sd) {
+
+ // Ignore PLACEHOLDER seed
+ if ( (*sd)->gid == TRIMMED ) {
+ continue;
+ }
+
+ // If seed has the lowest maximal number of errors set the values
+ CountType max_seed_errors = min_errors(*sd) + possible_remaining_errors;
+ if ( max_seed_errors < min_num_errors ) {
+ min_num_errors = max_seed_errors;
+ max_num_matches = (*sd)->num_matches;
+ }
+ else if ( max_seed_errors == min_num_errors ) {
+ max_num_matches = std::max( max_num_matches, (*sd)->num_matches );
+ }
+
+ }
+ }
+
+ // Fill the vector containing the best n scores for seed filtering decisions in all_best_n mode
+ else if ( globalAlignmentSettings.get_all_best_n_scores_mode() && globalAlignmentSettings.get_best_n() > 0 ) {
+
+ std::set<CountType> all_min_errors;
+ for(SeedVecIt sd = seeds.begin() ; sd !=seeds.end(); ++sd) {
+ if ( (*sd)->gid == TRIMMED ) {
+ continue;
+ }
+ all_min_errors.insert( min_errors(*sd) );
+ }
+
+ auto it = all_min_errors.begin();
+ if ( all_min_errors.size() > 0 ) {
+ std::advance(it, std::min( int(globalAlignmentSettings.get_best_n() - 1 ) , int (all_min_errors.size() - 1 ) ) );
+ min_num_errors = (*it) + possible_remaining_errors;
+ }
+ }
+
+ // All hit mode: Only consider the min_errors parameter
+ else {
+ min_num_errors = globalAlignmentSettings.get_min_errors();
+ }
// delete all seeds which do not reach threshold
SeedVecIt it=seeds.begin();
bool foundHit = false;
+
while ( it!=seeds.end()) {
- if ((*it)->num_matches < num_matches_threshold)
+
+ // Handle PLACEHOLDER seeds separately
+ if ( (*it)->gid == TRIMMED ) {
+
+ // Filter if last cycle
+ if ( cycle == total_cycles ) {
+ it = seeds.erase(it);
+ continue;
+ }
+
+ // Keep it otherwise
+ ++it;
+ continue;
+ }
+
+ CountType seed_errors = min_errors(*it);
+
+ // Filter all seeds that have more errors than the threshold
+ if ( seed_errors > min_num_errors ) {
it = seeds.erase(it);
- else if (settings.discard_ohw && (cycle>settings.start_ohw) && ((*it)->num_matches <= K_HiLive)) // remove one-hit-wonders
+ continue;
+ }
+
+ // Filter One-hit-Wonders
+ else if ( globalAlignmentSettings.get_discard_ohw() && (cycle>globalAlignmentSettings.get_start_ohw()) && ((*it)->num_matches <= globalAlignmentSettings.get_kmer_weight()) && ( (*it)->cigar_data.back().length > globalAlignmentSettings.get_max_consecutive_gaps() ) ) {
it = seeds.erase(it);
- else if (cycle == settings.seqlen && settings.any_best_hit_mode && foundHit)
+ continue;
+ }
+
+ // Filter suboptimal alignments in the last cycle for any_best mode
+ else if (cycle == total_cycles && globalAlignmentSettings.get_any_best_hit_mode() && ( (*it)->num_matches < max_num_matches || foundHit ) ) {
it = seeds.erase(it);
+ continue;
+ }
+
else
++it;
+
+ // If not filtered, a hit was found
foundHit = true;
}
- // if a new seed would have a chance then create it
- if (num_matches_threshold <= 1) // new seed would have 1 more match than expected by possibleRemainingMatches
- add_new_seeds(pos, posWasUsedForExtension, settings);
+ // Create new seeds if they have a chance to stay below the error threshold (Consider the number of matches given by a PLACEHOLDER seed)
+ CountType placeholder_matches = 0;
+ if ( seeds.size() > 0 && (*seeds.begin())->gid == TRIMMED ) {
+ placeholder_matches = (*seeds.begin())->num_matches;
+ }
+ if ( pos.size() != 0 && cycle < total_cycles && minErrors_in_region( cycle - placeholder_matches - globalAlignmentSettings.get_kmer_span(), 1) <= min_num_errors ) {
+ add_new_seeds(pos, posWasUsedForExtension);
+ }
}
// updates cigar_data accordingly to a new matching kmer
-void ReadAlignment::addMatchingKmer(USeed & s, DiffType offset, AlignmentSettings & settings) {
+void ReadAlignment::addMatchingKmer(USeed & s, DiffType offset) {
s->cigar_data.emplace_back(1,offset);
s->num_matches += 1;
@@ -530,7 +760,7 @@ void ReadAlignment::addMatchingKmer(USeed & s, DiffType offset, AlignmentSetting
//// split last kmer-span bases in single CigarElements ////
CigarVector::iterator it = --(s->cigar_data.end());
unsigned summedLength = 1;
- while (summedLength < settings.kmer_span) {
+ while (summedLength < globalAlignmentSettings.get_kmer_span()) {
++summedLength;
--it;
if ((*it).length > 1) {
@@ -565,7 +795,8 @@ void ReadAlignment::addMatchingKmer(USeed & s, DiffType offset, AlignmentSetting
unsigned positionInKmer = 1;
while (it != s->cigar_data.end()) {
// if positionInKmer is not in kmer_gaps and cigar element was match
- if ((*it).offset == NO_MATCH && std::find(settings.kmer_gaps.begin(), settings.kmer_gaps.end(), positionInKmer) == settings.kmer_gaps.end()) {
+ std::vector<unsigned> kmer_gaps = s->start_pos > 0 ? globalAlignmentSettings.get_kmer_gaps() : globalAlignmentSettings.get_rev_kmer_gaps();
+ if ( (*it).offset == NO_MATCH && std::find(kmer_gaps.begin(), kmer_gaps.end(), positionInKmer) == kmer_gaps.end()) {
(*it).offset = offset;
s->num_matches += 1;
}
@@ -590,144 +821,249 @@ void ReadAlignment::addMatchingKmer(USeed & s, DiffType offset, AlignmentSetting
// Extend an existing seed (to be precise, extend the CIGAR vector / data).
-bool ReadAlignment::extendSeed(USeed & s, DiffType offset, AlignmentSettings & settings){
- // Extend CIGAR for TRIMMED k-mers
- if ((offset == TRIMMED_MATCH) || ((offset == NO_MATCH) && (s->cigar_data.back().offset == TRIMMED_MATCH))){
- s->cigar_data.back().length += 1;
- if (s->cigar_data.back().offset != NO_MATCH){
+bool ReadAlignment::extendSeed(USeed & s, DiffType offset){
+
+ // TODO: Do we need to handle this?
+ // Extend placeholder seed only with trimmed k-mer
+ if ( s->gid == TRIMMED ) {
+
+ // Extend PLACEHOLDER when last k-mer is trimmed
+ if ( offset == TRIMMED_MATCH ) {
+ s->cigar_data.back().length += 1;
s->num_matches += 1;
return true;
}
+
+ // Everthing else can not happen, but if so, just return false.
return false;
}
-
- // if last region was match region
- if(s->cigar_data.back().offset != NO_MATCH){
- // MATCH -> NO_MATCH
- if(offset == NO_MATCH){
- // new mismatch region
- s->cigar_data.emplace_back(1,NO_MATCH);
- return false;
- }
- // MATCH -> MATCH Extend match region with new match.
- else {
- int offset_change = offset - s->cigar_data.rbegin()->offset;
- // without any mismatch in between there cannot be a valid offset_change other than 0
- if (offset_change!=0) {
- // new mismatch region
- s->cigar_data.emplace_back(1,NO_MATCH);
- return false;
- }
- // else: extend current match region
- addMatchingKmer(s, offset, settings);
- return true;
- }
- }
- // so last region was mismatch region
- else {
- // NO_MATCH -> NO_MATCH Extend mismatch region with new mismatch.
- if(offset == NO_MATCH){
- s->cigar_data.back().length += 1;
- return false;
- }
- // NO_MATCH -> MATCH
- else {
+
+ // Extend CIGAR for TRIMMED k-mers
+ if ( offset == TRIMMED_MATCH ) {
+
+ if ( s->cigar_data.back().offset == TRIMMED_MATCH ) {
+ s->cigar_data.back().length += 1;
+ }
+ else if ( s->cigar_data.back().offset == NO_MATCH ){
+ s->cigar_data.emplace_back(1,TRIMMED_MATCH);
+ }
+ else {
+ addMatchingKmer(s, s->cigar_data.back().offset);
+ }
+ return true;
+
+
+ }
+
+
+ // Extend CIGAR for NO_MATCH k-mer
+ else if ( offset == NO_MATCH ) {
+
+ // NO_MATCH --> NO_MATCH
+ if ( s->cigar_data.back().offset == NO_MATCH ) {
+ s->cigar_data.back().length += 1;
+ return false;
+ }
+
+ // TRIMMED_MATCH --> NO_MATCH
+ else if ( s->cigar_data.back().offset == TRIMMED_MATCH ) {
+ CountType trimmed_length = s->cigar_data.back().length;
+ s->cigar_data.erase( (--(s->cigar_data.end())) );
+ s->cigar_data.back().length += (trimmed_length + 1);
+ }
+
+ // MATCH --> NO_MATCH
+ else {
+ s->cigar_data.emplace_back(1, NO_MATCH);
+ return false;
+ }
+
+ }
+
+
+ // Extend CIGAR for MATCH k-mer
+ else {
+
+ // NO_MATCH --> MATCH
+ if ( s->cigar_data.back().offset == NO_MATCH ) {
+
assert((++(s->cigar_data.rbegin()))->offset != TRIMMED_MATCH && (++(s->cigar_data.rbegin()))->offset != NO_MATCH);
int offset_change = offset - (++(s->cigar_data.rbegin()))->offset;
+
// If there is an offset change, I need to have seen the appropriate mismatches before.
if ( offset_change == 0
- || ((offset_change < 0) && (s->cigar_data.back().length >= -offset_change + settings.kmer_span - 1)) // Insertion in read
- || ((offset_change > 0) && (s->cigar_data.back().length >= settings.kmer_span - 1 )) ) { // Deletion in read
- addMatchingKmer(s, offset, settings);
+ || ((offset_change < 0) && (s->cigar_data.back().length >= -offset_change + globalAlignmentSettings.get_kmer_span() - 1)) // Insertion in read
+ || ((offset_change > 0) && (s->cigar_data.back().length >= globalAlignmentSettings.get_kmer_span() - 1 )) ) { // Deletion in read
+ addMatchingKmer(s, offset);
return true;
+
+ // Appropriate mismatches not existing: Extend mismatch area
} else {
- // criteria not fulfilled: extend existing mismatch area.
s->cigar_data.back().length += 1;
return false;
}
- }
- }
+ }
+
+ // TRIMMED_MATCH --> MATCH
+ else if ( s->cigar_data.back().offset == TRIMMED_MATCH ) {
+
+ CountType trimmed_length = s->cigar_data.back().length;
+ s->cigar_data.erase( (--(s->cigar_data.end())) );
+
+ int offset_change = offset - (++(s->cigar_data.rbegin()))->offset;
+
+ // Insertion or deletion in read
+ if ( offset_change != 0 ) {
+
+ // Offset change is only considered for Insertions (negative offset change)
+ int considered_offsetChange = (offset_change < 0) ? offset_change : 0;
+
+ // Move TRIMMED_MATCHES from TRIMMED_MATCH region to previous NO_MATCH region such that the offset criteria are fulfilled.
+ // If this is not possible, count all trimmed MATCHes and current MATCH as NO_MATCH.
+
+ if ( (s->cigar_data.back().length < globalAlignmentSettings.get_kmer_span() - considered_offsetChange - 1) ) {
+
+ CountType required_nomatches = std::max(0, int(globalAlignmentSettings.get_kmer_span()) - considered_offsetChange - 1 - s->cigar_data.back().length);
+
+ if ( trimmed_length >= required_nomatches ) {
+ trimmed_length -= required_nomatches;
+ s->cigar_data.back().length += required_nomatches;
+ }
+
+ else {
+ s->cigar_data.back().length += (trimmed_length + 1);
+ return false;
+ }
+ }
+ }
+
+ // Add all previous TRIMMED k-mers as MATCH k-mers.
+ for ( CountType i = 0; i < trimmed_length + 1; i++ ) {
+ addMatchingKmer(s, offset);
+ }
+ return true;
+
+ }
+
+ // MATCH --> MATCH
+ else {
+
+ int offset_change = offset - s->cigar_data.rbegin()->offset;
+
+ // without any mismatch in between there cannot be a valid offset_change other than 0
+ if (offset_change!=0) {
+
+ // new mismatch region
+ s->cigar_data.emplace_back(1,NO_MATCH);
+ return false;
+
+ }
+
+ // else: extend current match region
+ addMatchingKmer(s, offset);
+ return true;
+ }
+ }
+
+ // Default: Should not be reached.
+ return false;
+
}
-void ReadAlignment::extend_alignment(char bc, KixRun* index, AlignmentSettings* settings) {
+void ReadAlignment::extend_alignment(char bc, KixRun* index, bool testRead) {
// move to the next cycle
cycle += 1;
+ // cycle is not allowed to be > total_cycles
+ assert( total_cycles >= cycle );
+
// update the last k-mer
uint8_t qual = ((bc >> 2) & 63); // get bits 3-8
- if ( (bc == 0) || (qual < settings->min_qual) ){ // no call if all 0 bits or quality below threshold
- last_invalid = cycle; // TODO append an N as basecall? Could be a bad idea
+ if ( (bc == 0) || (qual < globalAlignmentSettings.get_min_qual()) ){ // no call if all 0 bits or quality below threshold
+ last_invalid = last_invalid > cycle ? last_invalid : cycle; // TODO append an N as basecall? Could be a bad idea
}
- unsigned mask = 3;
- if (flags != 0) // if read is valid
- appendNucleotideToSequenceStoreVector(revtwobit_repr(bc & mask)); // get the nucleotide as an actual character, disregarding the quality
- // do not update the alignments when reading barcode
- if (cycle > settings->seqlen) {
- if (cycle == rlen)
- if (std::find(settings->barcodeVector.begin(), settings->barcodeVector.end(), getBarcodeString(*settings)) == settings->barcodeVector.end())
- this->disable(*settings);
- return;
- }
+ if (flags != 0) // if read is valid
+ appendNucleotideToSequenceStoreVector(bc); // get the nucleotide as an actual character, disregarding the quality
// do not update the alignments when reading the first kmer_span-1 cycles
- if (cycle < settings->kmer_span)
+ if (cycle < globalAlignmentSettings.get_kmer_span())
return;
// update the alignments
GenomePosListType pos;
std::vector<bool> posWasUsedForExtension;
// if last kmer of read is not valid
- if (!( last_invalid+settings->kmer_span-1 < cycle )) {
+ if (!( last_invalid+globalAlignmentSettings.get_kmer_span()-1 < cycle )) {
// write a NO_MATCH
for (auto sit = seeds.begin(); sit != seeds.end(); ++sit)
- extendSeed(*sit, NO_MATCH, *settings);
+ extendSeed(*sit, NO_MATCH);
+
+ // Remove placeholder if exist
+ if ( seeds.size() > 0 && (*seeds.begin())->gid == TRIMMED )
+ seeds.pop_front();
}
else {
// get all occurrences of last_kmer (fwd & rc) from index
- const std::string sequence = getSequenceString(*settings);
- std::string::const_iterator it_lastKmer = sequence.end() - settings->kmer_span;
+ const std::string sequence = getSequenceString();
+ std::string::const_iterator it_lastKmer = sequence.end() - globalAlignmentSettings.get_kmer_span();
HashIntoType last_kmer = 0;
- hash_fw(it_lastKmer, sequence.end(), last_kmer, *settings);
- pos = index->retrieve_positions(sequence.substr(sequence.length()-settings->kmer_span), *settings);
+ hash_fw(it_lastKmer, sequence.end(), last_kmer);
+ pos = index->retrieve_positions(sequence.substr(sequence.length()-globalAlignmentSettings.get_kmer_span()));
posWasUsedForExtension.resize(pos.size(), false);
// check if the current k-mer was trimmed in the index
if ( (pos.size() == 1) && ((*pos.begin()).gid == TRIMMED) ) {
- if(seeds.size()==0 && cycle==settings->kmer_span ){ // Create placeholder seed if first k-mer is trimmed
- create_placeholder_seed(*settings);
- } else
- // pretend that all existing seeds could be extended
- for(auto sd = seeds.begin() ; sd !=seeds.end(); ++sd)
- extendSeed(*sd, TRIMMED_MATCH, *settings);
+
+ // pretend that all existing seeds could be extended
+ for(auto sd = seeds.begin() ; sd !=seeds.end(); ++sd)
+ extendSeed(*sd, TRIMMED_MATCH);
+
+ if ( seeds.size() == 0 || (*seeds.begin())->gid != TRIMMED )
+ create_placeholder_seed();
// clear the pos list so nothing bad happens in the next steps
pos.clear();
}
// not trimmed in the index --> try to extend existing seeds
else {
- // Convert placeholder seeds (if any) to "real" seeds
- if(seeds.size()==1 && (*(seeds.begin()))->gid==TRIMMED){
- convertPlaceholder(pos, *settings);
- }
+
// find support for each candidate: iterate over seed candidates and positions simultaneously
auto cPos1 = pos.begin(), cPos2 = pos.begin(); // sliding window [cPos1, cPos2)
+// if ( pos.size() > 0 ) {
+// std::cout << cPos1->pos << std::endl;
+// }
for (auto cSeed = seeds.begin(); cSeed!=seeds.end(); ++cSeed ) {
- PositionType seed_pos = (*cSeed)->start_pos + cycle -settings->kmer_span;
+
+ // Don't handle PLACEHOLDER seed
+ if ( (*cSeed)->gid == TRIMMED ) {
+ continue;
+ }
+
+ // Compute the last offset of the current seed
+ PositionType last_offset = prev((*cSeed)->cigar_data.end())->offset;
+ if(last_offset == NO_MATCH) {
+ last_offset = prev(prev((*cSeed)->cigar_data.end()))->offset;
+ } else if ( last_offset == TRIMMED_MATCH ) {
+ last_offset = prev(prev(prev((*cSeed)->cigar_data.end())))->offset;
+ }
+
+ // Compute the optimal match position for the next k-mer
+ PositionType seed_pos = (*cSeed)->start_pos + cycle - globalAlignmentSettings.get_kmer_span() + last_offset;
// adjust the window in the position list
- while( (cPos1!=pos.end()) && (cPos1->pos < seed_pos - settings->window) )
+ while( (cPos1!=pos.end()) && (cPos1->pos < seed_pos - globalAlignmentSettings.get_window()) )
++cPos1;
- while( (cPos2!=pos.end()) && (cPos2->pos <= seed_pos + settings->window) )
+ while( (cPos2!=pos.end()) && (cPos2->pos <= seed_pos + globalAlignmentSettings.get_window()) )
++cPos2;
// search all positions in the window for the best matching extension of the seed
- DiffType best_offset = settings->window+1; // set larger than search window
+ DiffType best_offset = globalAlignmentSettings.get_window()+1; // set larger than search window
GenomePosListIt best_match = cPos2; // set behind the last element of the window
for(GenomePosListIt kmerHitIt = cPos1; kmerHitIt!=cPos2; ++kmerHitIt)
if (kmerHitIt->gid == (*cSeed)->gid){
@@ -741,29 +1077,91 @@ void ReadAlignment::extend_alignment(char bc, KixRun* index, AlignmentSettings*
// check if a best match was found for this seed
if (best_match != cPos2) {
- if(extendSeed(*cSeed, best_offset, *settings))
+ if(extendSeed(*cSeed, best_offset + last_offset))
// if pos was used as match, mark it so that later it does not get converted into a new seed
posWasUsedForExtension[best_match-pos.begin()] = true;
}
else{
// no position found to extend the current seed
- extendSeed(*cSeed, NO_MATCH, *settings);
+ extendSeed(*cSeed, NO_MATCH);
}
+
} // END: for(seeds...)
} // END: not trimmed
} // END: if last kmer is valid
- filterAndCreateNewSeeds(*settings, pos, posWasUsedForExtension);
+ filterAndCreateNewSeeds(pos, posWasUsedForExtension);
+
+ if ( testRead ) {
+ for ( auto seed = seeds.begin(); seed != seeds.end(); ++seed ) {
+ (*seed)->cout();
+ std::cout << "Seed's min errors: " << min_errors((*seed)) << std::endl;;
+ }
+ }
return;
}
+CountType ReadAlignment::getBarcodeIndex() {
+
+ // Get the barcodes of the read
+ std::string read_bc = getBarcodeString();
+
+ uint16_t fragment_errors = 0;
+ uint16_t fragment_pos = 0;
+ uint16_t fragment_num = 0;
+ uint16_t matching_bc = NO_MATCH;
+
+ // Iterate through all user-defined (multi-)barcodes
+ // That's quite complicated since the read barcodes are consecutive and the user barcodes are divided in vectors. // TODO: change that?
+ for ( uint16_t barcodeIndex = 0; barcodeIndex < globalAlignmentSettings.get_barcodeVector().size(); barcodeIndex++ ) {
+
+ // reset values for the barcode
+ fragment_errors = 0;
+ fragment_pos = 0;
+ fragment_num = 0;
+ matching_bc = barcodeIndex;
+
+ // for each base of the read barcode
+ for ( uint16_t nucl = 0; nucl < read_bc.length(); nucl++ ) {
+
+ // reset values for each barcode fragment
+ if ( fragment_pos >= (globalAlignmentSettings.get_barcodeVector()[barcodeIndex])[fragment_num].length() ) {
+ fragment_pos = 0;
+ fragment_num += 1;
+ fragment_errors = 0;
+ assert( fragment_num < (globalAlignmentSettings.get_barcodeVector()[barcodeIndex]).size() );
+ }
+
+ // compare nucleotides and increase the number of fragment errors if not equal
+ if ( read_bc.at(nucl) != (globalAlignmentSettings.get_barcodeVector()[barcodeIndex])[fragment_num].at(fragment_pos) ) {
+ fragment_errors++;
+ }
+
+ // if too many errors in a fragment, break the loop for the barcode
+ if ( fragment_errors > globalAlignmentSettings.get_barcode_errors()[fragment_num] ) {
+ matching_bc = NO_MATCH;
+ break;
+ }
+
+ fragment_pos += 1; // increment the fragment position
+
+ }
+
+ // if one barcode fulfilled the criteria, we can stop.
+ if ( matching_bc != NO_MATCH )
+ break;
+ }
+
+ return matching_bc;
+}
+
// disable this alignment, i.e. delete all seeds and set the last_invalid indicator to the
// end of the read. --> This read will not be aligned and consumes almost no space.
-void ReadAlignment::disable(AlignmentSettings & settings) {
- last_invalid = settings.seqlen;
+void ReadAlignment::disable() {
+ last_invalid = total_cycles;
seeds.clear();
flags = 0;
sequenceLen=0;
@@ -772,13 +1170,13 @@ void ReadAlignment::disable(AlignmentSettings & settings) {
// obtain start position of a seed according to SAM (leftmost)
-PositionType ReadAlignment::get_SAM_start_pos(USeed & sd, AlignmentSettings & settings) {
+PositionType ReadAlignment::get_SAM_start_pos(USeed & sd) {
PositionType pos = sd->start_pos;
if (pos < 0) {
if (sd->cigar_data.back().offset == NO_MATCH)
- pos = -pos - settings.seqlen + settings.kmer_span - (++sd->cigar_data.rbegin())->offset;
+ pos = -pos - total_cycles + globalAlignmentSettings.get_kmer_span() - (++sd->cigar_data.rbegin())->offset;
else
- pos = -pos - settings.seqlen + settings.kmer_span - (sd->cigar_data.rbegin())->offset;
+ pos = -pos - total_cycles + globalAlignmentSettings.get_kmer_span() - (sd->cigar_data.rbegin())->offset;
}
return pos;
}
diff --git a/lib/alnread.h b/lib/alnread.h
index 040a213..3b42a3d 100644
--- a/lib/alnread.h
+++ b/lib/alnread.h
@@ -30,7 +30,7 @@ struct Seed {
CigarVector cigar_data;
// return Seqans String of CigarElement
- seqan::String<seqan::CigarElement<> > returnSeqanCigarString();
+ seqan::String<seqan::CigarElement<> > returnSeqanCigarString(unsigned* nm_i);
// get the size of the serialized object
uint16_t serialize_size();
@@ -40,6 +40,8 @@ struct Seed {
// deserialize (read) data from a char vector
uint16_t deserialize(char* d);
+
+ void cout();
};
@@ -61,30 +63,46 @@ typedef SeedVec::iterator SeedVecIt;
class ReadAlignment {
private:
+
// read length
- CountType rlen;
+ CountType total_cycles;
- // sequence of the read so far, saved as vector<uint8_t> so interpretation is not that trivial. Contains barcode
+ // sequence of the read so far, saved as vector<uint8_t> so interpretation is not that trivial.
CountType sequenceLen=0;
std::vector<uint8_t> sequenceStoreVector;
- // Extend or create a placeholder seed for read with only trimmed matches
- void create_placeholder_seed(AlignmentSettings & settings);
+ // sequence of the barcode so far, saved as vector<uint8_t> so interpretation is not that trivial
+ CountType barcodeLen=0;
+ std::vector<uint8_t> barcodeStoreVector;
- // convert a placeholder seed to a set of normal seeds
- void convertPlaceholder(GenomePosListType& pos, AlignmentSettings & settings);
+ // Extend or create a placeholder seed for read with only trimmed matches
+ void create_placeholder_seed();
// Create new seeds from a list of kmer positions and add to current seeds
- void add_new_seeds(GenomePosListType& pos, std::vector<bool> & posWasUsedForExtension, AlignmentSettings & settings);
+ void add_new_seeds(GenomePosListType& pos, std::vector<bool> & posWasUsedForExtension);
+
+ /**
+ * This function is the modified pigeonhole principle holding for both spaced and unspaced kmers.
+ * It computes the minimum number of errors in an error region of a given CIGAR vector.
+ * An error region is a region that is surrounded by MATCH elements of length >= ( kmer_span - 1 ).
+ * The error region cannot contain MATCH elements of length >= ( kmer_span - 1 ).
+ *
+ * @param region_length Sum of all (!) elements within the error region, including involved MATCH elements.
+ * @param border Number of included borders of the CIGAR vector (begin and/or end). Must be in [0,2].
+ * @param Absolute number (positive) of the offset change during a region
+ * @return Minimum number of errors that caused a region of the given length.
+ * @author Tobias Loka, Jakob Schulze
+ */
+ CountType minErrors_in_region(CountType region_length, CountType border, CountType offset_change=0 );
// filter seeds based on filtering mode and q gram lemma. Also calls add_new_seeds.
- void filterAndCreateNewSeeds(AlignmentSettings & settings, GenomePosListType & pos, std::vector<bool> & posWasUsedForExtension);
+ void filterAndCreateNewSeeds(GenomePosListType & pos, std::vector<bool> & posWasUsedForExtension);
// updates cigar_data accordingly to a new matching kmer
- void addMatchingKmer(USeed & s, DiffType offset, AlignmentSettings & settings);
+ void addMatchingKmer(USeed & s, DiffType offset);
// Extend an existing CIGAR string for a seed based on a new basecall. return false if last CIGAR element after extension is mismatch area (NO_MATCH), true otherwise.
- bool extendSeed(USeed & s, DiffType offset, AlignmentSettings & settings);
+ bool extendSeed(USeed & s, DiffType offset);
public: // have everything public until the apropriate access functions are available
@@ -104,7 +122,7 @@ class ReadAlignment {
CountType max_num_matches;
// set the read_length
- void set_rlen(CountType r);
+ void set_total_cycles(CountType c);
// get the size of the serialized object
uint64_t serialize_size();
@@ -116,22 +134,51 @@ class ReadAlignment {
uint64_t deserialize(char* d);
// convert and return sequence of the read as string (without barcode)
- std::string getSequenceString(AlignmentSettings & settings);
-
- // convert and return sequence of the barcode
- std::string getBarcodeString(AlignmentSettings & settings);
-
- // append one nucleotide to sequenceStoreVector
- void appendNucleotideToSequenceStoreVector(char nuc);
+ std::string getSequenceString();
+
+ /**
+ * Convert and return sequence of the barcode. Multiple barcodes are concatenated (without delimiter).
+ * @return The Barcode as string
+ * @author Tobias Loka
+ */
+ std::string getBarcodeString();
+
+ /**
+ * Check whether the barcode of this read fulfills the criteria of at least one user-defined barcode.
+ * The nucleotides are only compared pairwise, not allowing for Indels.
+ * @return The index of the matching barcode in globalAlignmentSettings.multiBarcodeVector. NO_MATCH, if none.
+ * Also return NO_MATCH, if demultiplexing is not activated.
+ * @author Tobias Loka
+ */
+ CountType getBarcodeIndex() ;
+
+
+ /**
+ * Append one nucleotide to sequenceStoreVector
+ * @param nucl The nucleotide. Must be 2-bit-formatted.
+ * @param appendToBarcode If true, the nucleotide is appended to the barcode instead of the read sequence (default: false).
+ * @return
+ * @author Jakob Schulze
+ */
+ void appendNucleotideToSequenceStoreVector(char bc, bool appendToBarcode=false);
// extend the alignment by one basecall using reference database index
- void extend_alignment(char bc, KixRun* index, AlignmentSettings* settings);
+ void extend_alignment(char bc, KixRun* index, bool testRead=false);
// disable this alignment
- void disable(AlignmentSettings & settings);
+ void disable();
// obtain start position of a seed according to SAM (leftmost)
- PositionType get_SAM_start_pos(USeed & sd, AlignmentSettings & settings);
+ PositionType get_SAM_start_pos(USeed & sd);
+
+ /**
+ * Compute the minimum number of errors for a seed by using the modified pigeonhole principle implemented in ReadAlignment::minErrors_in_region.
+ *
+ * @param s The seed.
+ * @return The minimum number of errors for the given seed.
+ * @author Tobias Loka, Jakob Schulze
+ */
+ CountType min_errors(USeed & s);
}; // END class ReadAlignment
diff --git a/lib/alnstream.cpp b/lib/alnstream.cpp
index 88dc524..be5ad35 100644
--- a/lib/alnstream.cpp
+++ b/lib/alnstream.cpp
@@ -6,8 +6,8 @@
//-------------------------------------------------------------------//
// new output Alignment Stream class
-oAlnStream::oAlnStream(uint16_t ln, uint16_t tl, uint16_t cl, std::string rt, CountType rl, uint32_t nr, uint64_t bs, uint8_t fmt):
- lane(ln), tile(tl), cycle(cl), root(rt), rlen(rl), num_reads(nr), num_written(0), buffer(bs,0), buf_size(bs), buf_pos(0), format(fmt), ofile(NULL), ozfile(Z_NULL) {}
+oAlnStream::oAlnStream(uint16_t ln, uint16_t tl, uint16_t cl, CountType rl, uint32_t nr, uint64_t bs, uint8_t fmt):
+ lane(ln), tile(tl), cycle(cl), rlen(rl), num_reads(nr), num_written(0), buffer(bs,0), buf_size(bs), buf_pos(0), format(fmt), ofile(NULL), ozfile(Z_NULL) {}
// write function for lz4 compression
@@ -62,11 +62,6 @@ uint64_t oAlnStream::open(std::string fname) {
total_size += sizeof(uint16_t); // tile
total_size += sizeof(CountType); // cycle
- // root directory name
- uint16_t root_size = root.size();
- total_size += sizeof(uint16_t);
- total_size += root_size;
-
// read length
total_size += sizeof(CountType);
@@ -89,13 +84,6 @@ uint64_t oAlnStream::open(std::string fname) {
memcpy(d,&cycle,sizeof(CountType));
d += sizeof(CountType);
- // root directory name
- memcpy(d,&root_size,sizeof(uint16_t));
- d += sizeof(uint16_t);
-
- memcpy(d,root.c_str(),root_size);
- d += root_size;
-
// write the read length
memcpy(d,&rlen,sizeof(CountType));
d += sizeof(CountType);
@@ -149,7 +137,8 @@ uint64_t oAlnStream::write_alignment(ReadAlignment * al) {
case 1: written = gzwrite(ozfile, buffer.data(), buffer.size()); break;
case 2: written = lz4write(buffer.data(), buffer.size()); break;
}
- assert(written == buf_size);
+ if(written != buf_size)
+ throw std::runtime_error("Could not write out buffer 1 in oAlnStream::write_alignment.");
// copy remaining data
memcpy(buffer.data(),temp.data()+first_part,sizeof(uint32_t)-first_part);
@@ -172,7 +161,8 @@ uint64_t oAlnStream::write_alignment(ReadAlignment * al) {
case 1: written = gzwrite(ozfile, buffer.data(), buffer.size()); break;
case 2: written = lz4write(buffer.data(), buffer.size()); break;
}
- assert(written == buf_size);
+ if(written != buf_size)
+ throw std::runtime_error("Could not write out buffer 2 in oAlnStream::write_alignment.");
buf_pos = 0;
}
}
@@ -192,7 +182,8 @@ bool oAlnStream::close() {
case 1: written = gzwrite(ozfile, buffer.data(), buf_pos); break;
case 2: written = lz4write(buffer.data(), buf_pos); break;
}
- assert(written == buf_pos);
+ if(written != buf_pos)
+ throw std::runtime_error("Could not write out buffer in oAlnStream::close.");
buf_pos = 0;
if (num_written == num_reads) {
switch (format) {
@@ -220,7 +211,7 @@ bool oAlnStream::close() {
// new Alignment Stream class
iAlnStream::iAlnStream(uint64_t bs, uint8_t fmt):
- lane(0), tile(0), cycle(0), root(""), rlen(0), num_reads(0), num_loaded(0), buffer(bs,0), buf_size(bs), buf_pos(bs), format(fmt), ifile(NULL), izfile(Z_NULL) {}
+ lane(0), tile(0), cycle(0), rlen(0), num_reads(0), num_loaded(0), buffer(bs,0), buf_size(bs), buf_pos(bs), format(fmt), ifile(NULL), izfile(Z_NULL) {}
// read function for lz4 decompression, reads one block of data
@@ -273,7 +264,6 @@ uint64_t iAlnStream::open(std::string fname) {
// load the header:
uint64_t bytes = 0;
- uint16_t root_size;
switch (format) {
case 0: case 2:
{
@@ -283,14 +273,6 @@ uint64_t iAlnStream::open(std::string fname) {
bytes += fread(&tile,sizeof(uint16_t),1,ifile);
// read the cycle
bytes += fread(&cycle,sizeof(CountType),1,ifile);
- // root directory name size
- bytes += fread(&root_size,sizeof(uint16_t),1,ifile);
- // root name
- char * tmp = new char[root_size+1];
- bytes += fread(tmp,1,root_size,ifile);
- tmp[root_size] = 0; // make the string null-terminated
- root = tmp;
- delete tmp;
// read the read length
bytes += fread(&rlen,sizeof(CountType),1,ifile);
// read the number of reads
@@ -305,14 +287,6 @@ uint64_t iAlnStream::open(std::string fname) {
bytes += gzread(izfile,&tile,sizeof(uint16_t));
// read the cycle
bytes += gzread(izfile,&cycle,sizeof(CountType));
- // root directory name size
- bytes += gzread(izfile,&root_size,sizeof(uint16_t));
- // root name
- char * tmp = new char[root_size+1];
- bytes += gzread(izfile,tmp,root_size);
- tmp[root_size] = 0; // make the string null-terminated
- root = tmp;
- delete tmp;
// read the read length
bytes += gzread(izfile,&rlen,sizeof(CountType));
// read the number of reads
@@ -349,19 +323,15 @@ ReadAlignment* iAlnStream::get_alignment() {
memcpy(temp.data(),buffer.data()+buf_pos,first_part);
// load new buffer
- uint64_t loaded;
switch (format) {
case 0:
- loaded = fread(buffer.data(),1,buf_size,ifile);
- assert( (loaded == buf_size) || feof(ifile) );
+ fread(buffer.data(),1,buf_size,ifile);
break;
case 1:
- loaded = gzread(izfile,buffer.data(),buf_size);
- assert( (loaded == buf_size) || gzeof(izfile) );
+ gzread(izfile,buffer.data(),buf_size);
break;
case 2:
- loaded = lz4read_block();
- assert( loaded>0 || feof(ifile) );
+ lz4read_block();
break;
}
@@ -382,20 +352,16 @@ ReadAlignment* iAlnStream::get_alignment() {
// read new buffer from disk if necessary
if(buf_pos >= buf_size){
- uint64_t loaded;
switch (format) {
case 0:
- loaded = fread(buffer.data(),1,buf_size,ifile);
- assert( (loaded == buf_size) || feof(ifile) );
- break;
+ fread(buffer.data(),1,buf_size,ifile);
+ break;
case 1:
- loaded = gzread(izfile,buffer.data(),buf_size);
- assert( (loaded == buf_size) || gzeof(izfile) );
- break;
+ gzread(izfile,buffer.data(),buf_size);
+ break;
case 2:
- loaded = lz4read_block();
- assert( loaded>0 || feof(ifile) );
- break;
+ lz4read_block();
+ break;
}
buf_pos = 0;
}
@@ -403,7 +369,7 @@ ReadAlignment* iAlnStream::get_alignment() {
// finally, deserialize the alignment
ReadAlignment* ra = new ReadAlignment();
- ra->set_rlen(rlen);
+ ra->set_total_cycles(rlen);
ra->deserialize(data.data());
num_loaded++;
@@ -445,79 +411,72 @@ StreamedAlignment& StreamedAlignment::operator=(const StreamedAlignment& other)
lane = other.lane;
tile = other.tile;
- root = other.root;
rlen = other.rlen;
return *this;
}
-std::string StreamedAlignment::get_bcl_file(uint16_t cycle) {
+std::string StreamedAlignment::get_bcl_file(uint16_t cycle, uint16_t read_number) {
std::ostringstream path_stream;
-
- path_stream << root << "/L00" << lane << "/C" << cycle << ".1/s_"<< lane <<"_" << tile << ".bcl";
+ path_stream << globalAlignmentSettings.get_root() << "/L00" << lane << "/C" << getSeqCycle(cycle, read_number) << ".1/s_"<< lane <<"_" << tile << ".bcl";
return path_stream.str();
}
// get the path to the alignment file. The alignment file is located in
// <base>/L00<lane>/s_<lane>_<tile>.<cycle>.align
-// if base == "": base = root
-std::string StreamedAlignment::get_alignment_file(uint16_t cycle, std::string base){
+// if base == "": base = globalAlignmentSettings.get_root()
+std::string StreamedAlignment::get_alignment_file(uint16_t cycle, uint16_t mate, std::string base){
if (base == "") {
- base = root;
+ base = globalAlignmentSettings.get_root();
}
std::ostringstream path_stream;
- path_stream << base << "/L00" << lane << "/s_"<< lane << "_" << tile << "." << cycle << ".align";
+ path_stream << base << "/L00" << lane << "/s_"<< lane << "_" << tile << "." << mate << "."<< cycle << ".align";
return path_stream.str();
}
std::string StreamedAlignment::get_filter_file() {
std::ostringstream path_stream;
- path_stream << root << "/L00" << lane << "/s_"<< lane << "_" << tile << ".filter";
+ path_stream << globalAlignmentSettings.get_root() << "/L00" << lane << "/s_"<< lane << "_" << tile << ".filter";
return path_stream.str();
}
-// create directories required to store the alignment files (only if not stored in root)
-void StreamedAlignment::create_directories(AlignmentSettings* settings) {
+// create directories required to store the alignment files (only if not stored in globalAlignmentSettings.get_root())
+void StreamedAlignment::create_directories() {
std::ostringstream path_stream;
- if (settings->temp_dir == "") {
- path_stream << root;
+ if (globalAlignmentSettings.get_temp_dir() == "") {
+ path_stream << globalAlignmentSettings.get_root();
}
else {
- path_stream << settings->temp_dir;
+ path_stream << globalAlignmentSettings.get_temp_dir();
}
path_stream << "/L00" << lane;
boost::filesystem::create_directories(path_stream.str());
-
- std::ostringstream sam_stream;
- sam_stream << settings->out_dir;
- sam_stream << "/L00" << lane;
-
- boost::filesystem::create_directories(sam_stream.str());
+ boost::filesystem::create_directories(globalAlignmentSettings.get_out_dir());
}
// initialize empty alignment. Creates alignment files for a virtual Cycle 0
-void StreamedAlignment::init_alignment(AlignmentSettings* settings) {
- std::string out_fname = get_alignment_file(0, settings->temp_dir);
+void StreamedAlignment::init_alignment(uint16_t mate) {
+ std::string out_fname = get_alignment_file(0, mate, globalAlignmentSettings.get_temp_dir());
// get the number of reads in this tile by looking in the first bcl file
- std::string first_cycle = get_bcl_file(1);
+ std::string first_cycle = get_bcl_file(1, 0);
// extract the number of reads
uint32_t num_reads = num_reads_from_bcl(first_cycle);
// open output alignment stream
- oAlnStream output (lane, tile, 0, root, rlen, num_reads, settings->block_size, settings->compression_format);
+ oAlnStream output (lane, tile, 0, rlen, num_reads, globalAlignmentSettings.get_block_size(), globalAlignmentSettings.get_compression_format());
output.open(out_fname);
// write empty read alignments for each read
for (uint32_t i = 0; i < num_reads; ++i) {
ReadAlignment * ra = new ReadAlignment();
- ra->set_rlen(rlen);
+ ra->set_total_cycles(rlen);
output.write_alignment(ra);
delete ra;
}
@@ -530,20 +489,19 @@ void StreamedAlignment::init_alignment(AlignmentSettings* settings) {
// extend an existing alignment from cycle <cycle-1> to <cycle>. returns the number of seeds
-uint64_t StreamedAlignment::extend_alignment(uint16_t cycle, KixRun* index, AlignmentSettings* settings) {
+uint64_t StreamedAlignment::extend_alignment(uint16_t cycle, uint16_t read_no, uint16_t mate, KixRun* index) {
// 1. Open the input file
//-----------------------
- std::string in_fname = get_alignment_file(cycle-1, settings->temp_dir);
- std::string bcl_fname = get_bcl_file(cycle);
+ std::string in_fname = get_alignment_file(cycle-1, mate, globalAlignmentSettings.get_temp_dir());
+ std::string bcl_fname = get_bcl_file(cycle, read_no); // TODO: correct cycle
std::string filter_fname = get_filter_file();
- iAlnStream input ( settings->block_size, settings->compression_format );
+ iAlnStream input ( globalAlignmentSettings.get_block_size(), globalAlignmentSettings.get_compression_format() );
input.open(in_fname);
assert(input.get_cycle() == cycle-1);
assert(input.get_lane() == lane);
assert(input.get_tile() == tile);
- assert(input.get_root() == root);
assert(input.get_rlen() == rlen);
uint32_t num_reads = input.get_num_reads();
@@ -551,20 +509,15 @@ uint64_t StreamedAlignment::extend_alignment(uint16_t cycle, KixRun* index, Alig
// 2. Open output stream
//----------------------------------------------------------
- std::string out_fname = get_alignment_file(cycle, settings->temp_dir);
- oAlnStream output (lane, tile, cycle, root, rlen, num_reads, settings->block_size, settings->compression_format);
+ std::string out_fname = get_alignment_file(cycle, mate, globalAlignmentSettings.get_temp_dir());
+ oAlnStream output (lane, tile, cycle, rlen, num_reads, globalAlignmentSettings.get_block_size(), globalAlignmentSettings.get_compression_format());
output.open(out_fname);
-
// 3. Read the full BCL file (this is not too much)
//-------------------------------------------------
BclParser basecalls;
basecalls.open(bcl_fname);
-
- // extract the number of reads from the BCL file
- uint32_t num_base_calls = basecalls.size();
- assert(num_base_calls == num_reads);
// 4. Load the filter flags if filter file is available
@@ -585,21 +538,24 @@ uint64_t StreamedAlignment::extend_alignment(uint16_t cycle, KixRun* index, Alig
//-------------------------------------------------
uint64_t num_seeds = 0;
for (uint64_t i = 0; i < num_reads; ++i) {
+
+ bool testRead = false;
+
ReadAlignment* ra = input.get_alignment();
if (filters.size() > 0 && filters.has_next()) {
// filter file was found -> apply filter
if(filters.next()) {
- ra->extend_alignment(basecalls.next(), index, settings);
+ ra->extend_alignment(basecalls.next(), index, testRead);
num_seeds += ra->seeds.size();
}
else {
basecalls.next();
- ra->disable(*settings);
+ ra->disable();
}
}
// filter file was not found -> treat every alignment as valid
else {
- ra->extend_alignment(basecalls.next(), index, settings);
+ ra->extend_alignment(basecalls.next(), index, testRead);
num_seeds += ra->seeds.size();
}
@@ -607,6 +563,7 @@ uint64_t StreamedAlignment::extend_alignment(uint16_t cycle, KixRun* index, Alig
delete ra;
}
+
// 6. Close files
//-------------------------------------------------
if (!(input.close() && output.close())) {
@@ -615,199 +572,371 @@ uint64_t StreamedAlignment::extend_alignment(uint16_t cycle, KixRun* index, Alig
// 7. Delete old alignment file, if requested
//-------------------------------------------
- if ( ! settings->keep_aln_files ) {
+ if ( ! ( globalAlignmentSettings.get_keep_aln_files() || globalAlignmentSettings.is_output_cycle(cycle-1)) ) {
std::remove(in_fname.c_str());
}
return num_seeds;
}
-//-------------------------------------------------------------------//
-//------ Streamed SAM generation -----------------------------------//
-//-------------------------------------------------------------------//
+void StreamedAlignment::extend_barcode(uint16_t bc_cycle, uint16_t read_cycle, uint16_t read_no, uint16_t mate) {
-uint64_t alignments_to_sam(uint16_t ln, uint16_t tl, std::string rt, CountType rl, KixRun* index, AlignmentSettings* settings) {
- // set the file names
- std::string temp;
- if (settings->temp_dir == "")
- temp = rt;
- else
- temp = settings->temp_dir;
-
- std::string sam_dir = settings->out_dir;
- std::string filter_fname = filter_name(rt, ln, tl);
- std::string alignment_fname = alignment_name(temp, ln, tl, rl);
- std::string sam_fname = sam_tile_name(sam_dir, ln, tl, settings->write_bam);
-
- // check if files exist
- if ( !file_exists(alignment_fname) )
- throw std::runtime_error(std::string("Could not create SAM file. Alignment file not found: ")+ alignment_fname);
- if ( !file_exists(filter_fname) )
- std::cerr << "Could not find .filter file: " << filter_fname << std::endl;
-
- // open the alignment file
- iAlnStream input ( settings->block_size, settings->compression_format );
- input.open(alignment_fname);
- uint64_t num_reads = input.get_num_reads();
-
- // open the filter file, if applicable
- FilterParser filters;
- if (file_exists(filter_fname)) {
- filters.open(filter_fname);
- // extract the number of reads from the filter file
- uint32_t num_reads_filter = filters.size();
-
- if (num_reads != num_reads_filter){
- std::string msg = std::string("Number of reads in filter file (") + std::to_string(num_reads_filter) + ") does not match the number of reads in the BCL file (" + std::to_string(num_reads) + ").";
- throw std::length_error(msg.c_str());
- }
- }
+ // 1. Open the input file
+ //-----------------------
+ std::string in_fname = get_alignment_file(read_cycle, mate, globalAlignmentSettings.get_temp_dir());
+ std::string bcl_fname = get_bcl_file(bc_cycle, read_no);
+ std::string filter_fname = get_filter_file();
- /////////////////////////////////////////////
- // generate sam file
+ iAlnStream input ( globalAlignmentSettings.get_block_size(), globalAlignmentSettings.get_compression_format() );
+ input.open(in_fname);
+ assert(input.get_cycle() == read_cycle);
+ assert(input.get_lane() == lane);
+ assert(input.get_tile() == tile);
- // set BamFileOut object
- seqan::CharString samFileName = sam_fname;
- seqan::BamFileOut samFileOut(seqan::toCString(samFileName));
- seqan::StringSet<seqan::CharString> referenceNames;
- for (unsigned i=0; i<seqan::length(index->seq_names); i++) {
- seqan::appendValue(referenceNames, index->seq_names[i]);
- }
- seqan::NameStoreCache<seqan::StringSet<seqan::CharString> > referenceNamesCache(referenceNames);
- seqan::BamIOContext<seqan::StringSet<seqan::CharString> > bamIOContext(referenceNames, referenceNamesCache);
- samFileOut.context = bamIOContext;
+ uint32_t num_reads = input.get_num_reads();
- /////////////////
- // set SAM header
- seqan::BamHeaderRecord headerRecord;
-
- // @HD header.
- seqan::clear(headerRecord);
- headerRecord.type = seqan::BAM_HEADER_FIRST;
- seqan::resize(headerRecord.tags, 2);
- headerRecord.tags[0].i1 = "VN";
- headerRecord.tags[0].i2 = "1.5";
- headerRecord.tags[1].i1 = "GO";
- headerRecord.tags[1].i2 = "query";
- seqan::writeHeader(samFileOut, headerRecord);
-
- // @SQ header.
- std::stringstream ss;
- for ( uint64_t i = 0; i < index->seq_names.size(); i++ ) {
- ss.str(std::string()); // clear string stream
- seqan::clear(headerRecord);
- headerRecord.type = seqan::BAM_HEADER_REFERENCE;
- seqan::resize(headerRecord.tags, 2);
- headerRecord.tags[0].i1 = "SN";
- headerRecord.tags[0].i2 = index->seq_names[i];
-
- headerRecord.tags[1].i1 = "LN";
- ss << index->seq_lengths[i];
- headerRecord.tags[1].i2 = ss.str();
-
- seqan::writeHeader(samFileOut, headerRecord);
- }
-
- // @PG header.
- seqan::clear(headerRecord);
- headerRecord.type = seqan::BAM_HEADER_PROGRAM;
- seqan::resize(headerRecord.tags, 3);
- headerRecord.tags[0].i1 = "ID";
- headerRecord.tags[0].i2 = "hilive";
- headerRecord.tags[1].i1 = "PN";
- headerRecord.tags[1].i2 = "HiLive";
- headerRecord.tags[2].i1 = "VN";
- ss.str(std::string());
- ss << HiLive_VERSION_MAJOR << "." << HiLive_VERSION_MINOR;
- headerRecord.tags[2].i2 = ss.str();
- seqan::writeHeader(samFileOut, headerRecord);
-
+ // 2. Open output stream
+ //----------------------------------------------------------
+ std::string out_fname = in_fname + ".temp";
+ oAlnStream output (lane, tile, read_cycle, input.get_rlen(), num_reads, globalAlignmentSettings.get_block_size(), globalAlignmentSettings.get_compression_format());
+ output.open(out_fname);
- /////////////////
- // prepare sam entries
- uint64_t num_alignments = 0;
- for (uint64_t i = 0; i < num_reads; ++i) {
- ReadAlignment * ra = input.get_alignment();
-
- // if either: filter file is open and all filter flags were loaded and the filter flag is > 0
- // or: the filter file is not available
- // then proceed else skip
- if (!((filters.size()>0 && filters.next()) || filters.size() == 0))
- continue;
-
- num_alignments += ra->seeds.size();
-
- // Read name format <instrument‐name>:<run ID>:<flowcell ID>:<lane‐number>:<tile‐number>:<x‐pos>:<y‐pos>
- //readname << "<instrument>:<run-ID>:<flowcell-ID>:" << ln << ":" << tl << ":<xpos>:<ypos>:" << i;
- std::stringstream readname;
- readname << "lane." << ln << "|tile." << tl << "|read." << i;
-
- /////////////////
- // set sam entries
- seqan::BamAlignmentRecord record;
- for (SeedVecIt it = ra->seeds.begin(); it != ra->seeds.end(); ++it) {
- seqan::clear(record);
-
- record.qName = readname.str();
-
- record.rID = (*it)->gid;
-
- record.beginPos = ra->get_SAM_start_pos(*it, *settings)-1; // seqan expects 0-based positions, but in HiLive we use 1-based
- if (record.beginPos < 0)
- continue;
-
- record.cigar = (*it)->returnSeqanCigarString();
-
- // flag and seq
- record.flag = 0;
- record.seq = ra->getSequenceString(*settings);
- if ((*it)->start_pos < 0) { // if read matched reverse complementary
- seqan::reverseComplement(record.seq);
- record.flag |= 16;
- }
- if (it != ra->seeds.begin()) { // if current seed is secondary alignment
- record.flag |= 256;
- seqan::clear(record.seq);
- record.qual = "*";
- }
-
-
- // check if cigar string sums up to read length
- unsigned cigarElemSum = 0;
- unsigned deletionSum = 0;
- for (seqan::Iterator<seqan::String<seqan::CigarElement<> > >::Type elem = seqan::begin(record.cigar); elem != end(record.cigar); ++elem) {
- if ((elem->operation == 'M') || (elem->operation == 'I') || (elem->operation == 'S') || (elem->operation == '=') || (elem->operation == 'X'))
- cigarElemSum += elem->count;
- if (elem->operation == 'D')
- deletionSum += elem->count;
- }
- if (cigarElemSum != settings->seqlen) {
- std::cerr << "WARNING: Excluded an alignment of read " << record.qName << " at position " << ra->get_SAM_start_pos(*it, *settings) << " because its cigar vector had length " << cigarElemSum << std::endl;
- continue;
- }
-
- // tags
- seqan::BamTagsDict dict;
- seqan::appendTagValue(dict, "AS", (*it)->num_matches);
- if (settings->seqlen < settings->rlen) // if demultiplexing is on
- seqan::appendTagValue(dict, "BC", ra->getBarcodeString(*settings));
- seqan::appendTagValue(dict, "NM", deletionSum + settings->seqlen - (*it)->num_matches);
- record.tags = seqan::host(dict);
-
-
- // write record to disk
- seqan::writeRecord(samFileOut, record);
- }
- }
-
- std::ofstream statsfile;
- statsfile.open(sam_fname+".stats");
- statsfile << "Number of reads\t" << num_reads << std::endl;
- statsfile << "Number of alignments\t" << num_alignments << std::endl;
- statsfile.close();
+ // 3. Read the full BCL file (this is not too much)
+ //-------------------------------------------------
+ BclParser basecalls;
+ basecalls.open(bcl_fname);
+
+
+ // 4. Extend barcode sequence
+ //-------------------------------------------------
+ for (uint64_t i = 0; i < num_reads; ++i) {
+ char bc = basecalls.next();
+ ReadAlignment* ra = input.get_alignment();
+ ra->appendNucleotideToSequenceStoreVector(bc, true);
+
+ // filter invalid barcodes if new barcode fragment is completed
+ // TODO: Is done for each mate. Check if it's worth to change it (runtime should not be too high?)
+ if ( !globalAlignmentSettings.get_keep_all_barcodes() && bc_cycle == globalAlignmentSettings.get_seqs()[read_no].length && ra->getBarcodeIndex() == NO_MATCH ) {
+ ra->disable();
+ }
+
+ output.write_alignment(ra);
+ delete ra;
+ }
+
+ // 5. Close files
+ //-------------------------------------------------
+ if (!(input.close() && output.close())) {
+ std::cerr << "Could not finish alignment!" << std::endl;
+ }
+
+ // 6. Move temp out file to the original file.
+ //-------------------------------------------
+ std::rename(out_fname.c_str(), in_fname.c_str());
+
+}
+
+
+//-------------------------------------------------------------------//
+//------ Streamed SAM generation -----------------------------------//
+//-------------------------------------------------------------------//
+
+uint64_t alignments_to_sam(std::vector<uint16_t> lns, std::vector<uint16_t> tls, KixRun* index, CountType cycle) {
+
+ std::ofstream logfile;
+ logfile.open( get_out_log_name(), std::ofstream::app );
+
+ logfile << "Start to write output for cycle " << std::to_string(cycle) << "." << std::endl;
+ logfile << "Only alignments with a relative alignment score of >= " << std::to_string(globalAlignmentSettings.get_min_as_ratio()) << " are considered." << std::endl;
+
+ // Fill list of specified barcodes
+ std::vector<std::string> barcodes;
+ for ( unsigned i = 0; i < globalAlignmentSettings.get_barcodeVector().size(); i++ ) {
+ barcodes.push_back(globalAlignmentSettings.get_barcodeString(i));
+ }
+
+ // Get the finished cycles for each mate
+ std::vector<CountType> mateCycles;
+ for ( CountType mate = 1; mate <= globalAlignmentSettings.get_mates(); mate++ ) {
+ mateCycles.push_back( getMateCycle( mate, cycle ) );
+ }
+
+ // Init the bamIOContext (the same object can be used for all output streams)
+ seqan::StringSet<seqan::CharString> referenceNames;
+ seqan::NameStoreCache<seqan::StringSet<seqan::CharString> > referenceNamesCache(referenceNames);
+ seqan::BamIOContext<seqan::StringSet<seqan::CharString> > bamIOContext(referenceNames, referenceNamesCache);
+
+ seqan::contigNames(bamIOContext) = index->seq_names;
+ seqan::contigLengths(bamIOContext) = index->seq_lengths;
+
+ // Init the header (the same object can be used for all output streams)
+ seqan::BamHeader header = getBamHeader();
+
+ // Vector that contains the output streams (USE POINTERS !)
+ std::vector<std::unique_ptr<seqan::BamFileOut>> bfos;
+
+ // Init output stream for each barcode (plus undetermined if keep_all_barcodes is set)
+ for ( unsigned barcode=0; barcode < (barcodes.size() + 1); barcode ++) {
+ if ( barcode < barcodes.size() || globalAlignmentSettings.get_keep_all_barcodes() ) {
+
+ std::string barcode_string = ( barcode == barcodes.size() ) ? "undetermined" : barcodes[barcode];
+
+ // Open file in Bam output stream and write the header
+ std::unique_ptr<seqan::BamFileOut> bfo( new seqan::BamFileOut(getBamTempFileName(barcode_string, cycle).c_str()));
+ bfos.push_back( std::move(bfo) );
+ bfos[barcode]->context = bamIOContext;
+ seqan::writeHeader(*bfos[barcode], header);
+ }
+ }
+
+
+ ////////////////////////////////////////////////////
+ // Main loop //////////////////////////////////////
+ ////////////////////////////////////////////////////
+
+ uint64_t num_alignments = 0;
+ unsigned totalNumberOfReads = 0;
+
+ // for all lanes
+ /////////////////////////////////////////////////////////////////////////////
+ for (auto ln:lns) {
+
+ // for all tiles
+ /////////////////////////////////////////////////////////////////////////////
+ for (auto tl:tls) {
+
+ // set the filter file
+ std::string filter_fname = filter_name(ln, tl);
+ FilterParser filters;
+ if (file_exists(filter_fname)) {
+ filters.open(filter_fname);
+ }
+ else
+ logfile << "Could not find .filter file: " << filter_fname << ". Treated all reads as valid."<< std::endl;
+
+
+ // set the alignment files
+ std::vector<iAlnStream*> alignmentFiles;
+ unsigned numberOfAlignments = 0;
+ for (unsigned mateIndex = 1; mateIndex <= mateCycles.size(); mateIndex++) {
+
+ if ( globalAlignmentSettings.getSeqByMate(mateIndex) == NULLSEQ )
+ return 0;
+
+ CountType mateCycle = mateCycles[mateIndex-1];
+
+ // Open alignment file
+ std::string alignment_fname = alignment_name(ln, tl, mateCycle, mateIndex);
+ if ( !file_exists(alignment_fname) ) {
+ logfile << "Alignment file not found: " << alignment_fname << ". Ignored all related alignments." << std::endl;
+ continue;
+ }
+ iAlnStream* input = new iAlnStream( globalAlignmentSettings.get_block_size(), globalAlignmentSettings.get_compression_format() );
+ input->open(alignment_fname);
+
+ // compare number of reads in alignment file with number of reads in filter file, if filter file exists
+ if (file_exists(filter_fname) && input->get_num_reads() != filters.size()) {
+ std::string msg = std::string("Number of reads in filter file (") + std::to_string(filters.size()) + ") does not match the number of reads in the alignment file (" + std::to_string(input->get_num_reads()) + ").";
+ logfile << "Writing output for cycle " << std::to_string(cycle) << "not successful: " << msg.c_str() << std::endl;
+ throw std::length_error(msg.c_str());
+ }
+
+ // compare number of reads in alignment file with number of reads in previous alignment file
+ if (mateIndex != 1 && input->get_num_reads() != numberOfAlignments) {
+ std::string msg = std::string("Number of reads in alignment file (") + std::to_string(input->get_num_reads()) + ") does not match the number of reads in previous alignment file (" + std::to_string(input->get_num_reads()) + ").";
+ logfile << "Writing output for cycle " << std::to_string(cycle) << "not successful: " << msg.c_str() << std::endl;
+ throw std::length_error(msg.c_str());
+ }
+
+ numberOfAlignments = input->get_num_reads(); // set this after last if-then construct
+ alignmentFiles.push_back(input);
+ }
+
+ totalNumberOfReads += numberOfAlignments;
+
+ // for all reads in a tile
+ /////////////////////////////////////////////////////////////////////////////
+ for (uint64_t i = 0; i < numberOfAlignments; ++i) {
+
+ std::vector<ReadAlignment*> mateAlignments;
+ for (auto e:alignmentFiles) {
+ mateAlignments.push_back(e->get_alignment());
+ }
+
+ // if the filter file is available and the filter flag is 0 then skip
+ if (filters.size() != 0 && filters.next() == false)
+ continue;
+
+ // compute barcode sequence as it should be written to BC tag
+ std::string barcode = globalAlignmentSettings.format_barcode(mateAlignments[0]->getBarcodeString());
+
+ // Barcode index for the read
+ CountType barcodeIndex = mateAlignments[0]->getBarcodeIndex();
+
+ // If read has undetermined barcode and keep_all_barcodes is not set, skip this read
+ if ( barcodeIndex == NO_MATCH && !globalAlignmentSettings.get_keep_all_barcodes() )
+ continue;
+ else if ( barcodeIndex == NO_MATCH )
+ barcodeIndex = barcodes.size(); // this is the index for the "undetermined" output stream
+
+ // setup QNAME
+ // Read name format <instrument‐name>:<run ID>:<flowcell ID>:<lane‐number>:<tile‐number>:<x‐pos>:<y‐pos>
+ // readname << "<instrument>:<run-ID>:<flowcell-ID>:" << ln << ":" << tl << ":<xpos>:<ypos>:" << i;
+ //TODO: where do we get the Illumina read coordinate from?
+ std::stringstream readname;
+ readname << "lane." << ln << "|tile." << tl << "|read." << i;
+
+
+ // for all mates
+ /////////////////////////////////////////////////////////////////////////////
+ seqan::BamAlignmentRecord record;
+ unsigned printedMates = 0;
+ for (unsigned mateAlignmentIndex=0; mateAlignmentIndex < mateAlignments.size(); ++mateAlignmentIndex) {
+
+ // for all seeds
+ /////////////////////////////////////////////////////////////////////////////
+ for (SeedVecIt it = mateAlignments[mateAlignmentIndex]->seeds.begin(); it != mateAlignments[mateAlignmentIndex]->seeds.end(); ) {
+ if ( (*it)->gid == TRIMMED ) {
+ //TODO: count trimmed reads for output stats
+ ++it;
+ continue;
+ }
+
+ seqan::clear(record);
+
+ record.qName = readname.str();
+
+ record.rID = (*it)->gid;
+
+ record.beginPos = mateAlignments[mateAlignmentIndex]->get_SAM_start_pos(*it)-1; // seqan expects 0-based positions, but in HiLive we use 1-based
+ if (record.beginPos < 0) {
+ it = mateAlignments[mateAlignmentIndex]->seeds.erase(it);
+ continue;
+ }
+
+ unsigned nm_i = 0;
+ record.cigar = (*it)->returnSeqanCigarString(&nm_i);
+
+ // flag and seq
+ record.flag = 0;
+ record.seq = mateAlignments[mateAlignmentIndex]->getSequenceString();
+ if ((*it)->start_pos < 0) { // if read matched reverse complementary
+ seqan::reverseComplement(record.seq);
+ record.flag |= 16;
+ }
+ if (printedMates >= mateAlignments.size()) { // if current seed is secondary alignment
+ record.flag |= 256;
+ seqan::clear(record.seq);
+ seqan::clear(record.qual);
+ }
+ if (globalAlignmentSettings.get_mates() > 1) { // if there are more than two mates
+ record.flag |= 1;
+ if (mateAlignmentIndex == 0) {
+ record.flag |= 64;
+ } else if (mateAlignmentIndex == mateAlignments.size()-1) {
+ record.flag |= 128;
+ } else {
+ record.flag |= 192; // 64 + 128
+ }
+
+ bool eachMateAligned = true;
+ for (auto e:mateAlignments)
+ eachMateAligned = eachMateAligned && e->seeds.size() > 0;
+ if (eachMateAligned)
+ record.flag |= 2;
+ }
+
+
+ // check if cigar string sums up to read length
+ // TODO Potentially conflicts with the 'eachMateAligned' flag if done here.
+ unsigned cigarElemSum = 0;
+ unsigned deletionSum = 0;
+ unsigned supposed_cigar_length = mateCycles[mateAlignmentIndex];
+ CountType min_as_score = mateCycles[mateAlignmentIndex] * globalAlignmentSettings.get_min_as_ratio();
+
+ unsigned as_score = (*it)->num_matches;
+ for (seqan::Iterator<seqan::String<seqan::CigarElement<> > >::Type elem = seqan::begin(record.cigar); elem != end(record.cigar); ++elem) {
+ if ((elem->operation == 'M') || (elem->operation == 'I') || (elem->operation == 'S') || (elem->operation == '=') || (elem->operation == 'X'))
+ cigarElemSum += elem->count;
+ if ( (elem->operation == 'I') )
+ as_score += elem->count - 1;
+ if (elem->operation == 'D') {
+ deletionSum += elem->count;
+ as_score -= 1;
+ }
+ }
+ if (cigarElemSum != supposed_cigar_length) {
+ logfile << "WARNING: Excluded an alignment of read " << record.qName << " at position " << mateAlignments[mateAlignmentIndex]->get_SAM_start_pos(*it) << " because its cigar vector had length " << cigarElemSum << std::endl;
+ it = mateAlignments[mateAlignmentIndex]->seeds.erase(it);
+ continue;
+ }
+ if (deletionSum >= supposed_cigar_length) {
+ logfile << "WARNING: Excluded an alignment of read " << record.qName << " at position " << mateAlignments[mateAlignmentIndex]->get_SAM_start_pos(*it) << " because its cigar vector had " << deletionSum << " deletions" << std::endl;
+ it = mateAlignments[mateAlignmentIndex]->seeds.erase(it);
+ continue;
+ }
+
+ if (as_score < min_as_score) {
+// logfile << "WARNING: Excluded an alignment of read " << record.qName << " at position " << mateAlignments[mateAlignmentIndex]->get_SAM_start_pos(*it) << " because its AS:i score was below the threshold." << std::endl;
+ it = mateAlignments[mateAlignmentIndex]->seeds.erase(it);
+ continue;
+ }
+
+ // tags
+ seqan::BamTagsDict dict;
+ seqan::appendTagValue(dict, "AS", ( as_score ) );
+ if (barcode!="") { // if demultiplexing is on
+ seqan::appendTagValue(dict, "BC", barcode);
+ }
+ seqan::appendTagValue(dict, "NM", nm_i);
+ record.tags = seqan::host(dict);
+
+
+ // write record to disk
+ seqan::writeRecord(*bfos[barcodeIndex], record);
+
+ ++num_alignments;
+ ++printedMates;
+ ++it;
+ }
+ }
+
+ for (auto e:mateAlignments)
+ delete e;
+ }
+ for (auto e:alignmentFiles)
+ delete e;
+ }
+
+ }
+
+ // Init output stream for each barcode (plus undetermined if keep_all_barcodes is set)
+ for ( unsigned barcode=0; barcode < barcodes.size() + 1; barcode ++) {
+ if ( barcode < barcodes.size() || globalAlignmentSettings.get_keep_all_barcodes() ) {
+
+ std::string barcode_string = ( barcode == barcodes.size() ) ? "undetermined" : barcodes[barcode];
+
+ std::rename(getBamTempFileName(barcode_string, cycle).c_str(), getBamFileName(barcode_string, cycle).c_str());
+ }
+ }
+
+ // TODO maybe find a way to generate statsfiles when generating multiple output files.
+ std::ofstream statsfile;
+ std::string statsfile_fname = globalAlignmentSettings.get_out_dir() + "/hilive_out" + "_cycle" + std::to_string(cycle) + ".stats";
+ statsfile.open( statsfile_fname );
+ statsfile << "Number of reads\t" << totalNumberOfReads << std::endl;
+ statsfile << "Number of alignments\t" << num_alignments << std::endl;
+ statsfile.close();
+
+ // TODO What todo with completely trimmed reads?
+ // --> Amount could for example be written in the stats file(s)
+
+ logfile << "Writing output for cycle " << std::to_string(cycle) << " finished." << std::endl;
+ logfile.close();
+
+ return 1;
- return 0;
}
diff --git a/lib/alnstream.h b/lib/alnstream.h
index 613feff..e55ac38 100644
--- a/lib/alnstream.h
+++ b/lib/alnstream.h
@@ -3,13 +3,11 @@
#include "headers.h"
#include "definitions.h"
+#include "global_variables.h"
#include "kindex.h"
#include "tools.h"
#include "alnread.h"
#include "illumina_parsers.h"
-#include <seqan/basic.h>
-#include <seqan/sequence.h>
-#include <seqan/bam_io.h>
// Output alignment stream: write alignments to file one by one
class oAlnStream {
@@ -17,7 +15,6 @@ class oAlnStream {
uint16_t lane;
uint16_t tile;
uint16_t cycle;
- std::string root;
CountType rlen;
uint32_t num_reads;
@@ -48,7 +45,7 @@ class oAlnStream {
public:
// constructor initializes all member variables
- oAlnStream(uint16_t ln, uint16_t tl, uint16_t cl, std::string rt, CountType rl, uint32_t nr, uint64_t bs, uint8_t fmt);
+ oAlnStream(uint16_t ln, uint16_t tl, uint16_t cl, CountType rl, uint32_t nr, uint64_t bs, uint8_t fmt);
// open Alignment stream file and write header
uint64_t open(std::string fname);
@@ -69,7 +66,6 @@ class iAlnStream {
uint16_t lane;
uint16_t tile;
uint16_t cycle;
- std::string root;
CountType rlen;
uint32_t num_reads;
@@ -116,7 +112,6 @@ class iAlnStream {
inline uint16_t get_lane() {return lane;};
inline uint16_t get_tile() {return tile;};
inline uint16_t get_cycle() {return cycle;};
- inline std::string get_root() {return root;};
inline CountType get_rlen() {return rlen;};
inline uint32_t get_num_reads() {return num_reads;};
inline uint32_t get_num_loaded() {return num_loaded;};
@@ -132,7 +127,6 @@ class StreamedAlignment {
// dataset information
uint16_t lane;
uint16_t tile;
- std::string root; // the BaseCalls directory
CountType rlen;
// fetch the next read from the input stream
@@ -142,30 +136,42 @@ class StreamedAlignment {
uint64_t write_alignment(ReadAlignment& ral);
// get the path to the bcl file of a given cycle
- std::string get_bcl_file(uint16_t cycle);
+ std::string get_bcl_file(uint16_t cycle, uint16_t read_number);
// get the path to the alignment file. The alignment file is located in
// <base>/L00<lane>/s_<lane>_<tile>.<cycle>.align
- // if base == "": base = root
- std::string get_alignment_file(uint16_t cycle, std::string base = "");
+ // if base == "": base = globalAlignmentSettings.get_root()
+ std::string get_alignment_file(uint16_t cycle, uint16_t mate, std::string base = "");
// get the path to the filter file. The illumina filter information is located in
// <root>/L00<lane>/s_<lane>_<tile>.filter
std::string get_filter_file();
public:
- StreamedAlignment(uint16_t ln, uint16_t tl, std::string rt, CountType rl): lane(ln), tile(tl), root(rt), rlen(rl) {};
+ StreamedAlignment(uint16_t ln, uint16_t tl, CountType rl): lane(ln), tile(tl), rlen(rl) {};
StreamedAlignment& operator=(const StreamedAlignment& other);
// create directories required to store the alignment files (only if not stored in root)
- void create_directories(AlignmentSettings* settings);
+ void create_directories();
// initialize empty alignment. Creates files for a virtual Cycle 0
- void init_alignment(AlignmentSettings* settings);
+ void init_alignment(uint16_t mate);
// extend an existing alignment from cycle <cycle-1> to <cycle>
- uint64_t extend_alignment(uint16_t cycle, KixRun* index, AlignmentSettings* settings);
+ uint64_t extend_alignment(uint16_t cycle, uint16_t read_no, uint16_t mate, KixRun* index);
+
+ /**
+ * Extend the barcode for all reads with the information of the current sequencing cycle.
+ * @param bc_cycle The cycle of the barcode read.
+ * @param read_cycle The last handled cycle for the respective mate (should always be 0 or the full length)
+ * @param read_no The number of the sequence read for which the barcode will be extended (:= index in globalAlignmentSettings.seqs).
+ * @param mate The read mate to extend the barcode.
+ * @return
+ * @author Tobias Loka
+ */
+ void extend_barcode(uint16_t bc_cycle, uint16_t read_cycle, uint16_t read_no, uint16_t mate);
+
}; /* END class StreamedAlignment */
@@ -175,7 +181,6 @@ class StreamedAlignment {
//------ Streamed SAM generation -----------------------------------//
//-------------------------------------------------------------------//
-uint64_t alignments_to_sam(uint16_t ln, uint16_t tl, std::string rt, CountType rl, KixRun* index, AlignmentSettings* settings);
-
+uint64_t alignments_to_sam(std::vector<uint16_t> lns, std::vector<uint16_t> tls, KixRun* index, CountType cycle);
#endif /* ALNSTREAM_H */
diff --git a/lib/argument_parser.cpp b/lib/argument_parser.cpp
index 499e631..075858f 100644
--- a/lib/argument_parser.cpp
+++ b/lib/argument_parser.cpp
@@ -1,231 +1,625 @@
#include "argument_parser.h"
+
namespace po = boost::program_options;
+//|||||||||||||||||||||||||||||||||||||||||||||||||||||||
+//-----ArgumentParser------------------------------------
+//|||||||||||||||||||||||||||||||||||||||||||||||||||||||
+
+ArgumentParser::ArgumentParser(int argC, char const ** argV):argc(argC),argv(argV){}
+
+//|||||||||||||||||||||||||||||||||||||||||||||||||||||||
+//-----BuildIndexArgumentParser--------------------------
+//|||||||||||||||||||||||||||||||||||||||||||||||||||||||
+
+po::options_description BuildIndexArgumentParser::general_options() {
+ po::options_description general("General");
+ general.add_options()
+ ("help,h", "Print this help message and exit")
+ ("license", "Print licensing information and exit");
+ return general;
+}
+
+po::options_description BuildIndexArgumentParser::positional_options() {
+ po::options_description parameters("Required");
+ parameters.add_options()
+ ("INPUT", po::value<std::string>()->required(), "Reference genomes in (multi-) FASTA format.")
+ ("KMER_WEIGHT", po::value<uint16_t>()->required(), "Number of non-gap positions in a k-mer (For ungapped k-mers this is the k-mer size).");
+ return parameters;
+}
+
+po::options_description BuildIndexArgumentParser::build_options() {
+ po::options_description options("Options");
+ options.add_options()
+ ("outfile,o", po::value<std::string>(), "Set output file name [Default: INPUT.kix]")
+ ("trim,t", po::value<unsigned>(&trim)->default_value(0), "Ignore k-mers with more than t occurrences. [Default: no limit]")
+ ("gap-positions,p", po::value< std::vector<unsigned> >()->multitoken()->composing(), "Gap positions in the k-mer pattern (example: -p 3 6 7 for 1101100111 with k=7). [Default: ungapped]")
+ ("do-not-convert-spaces", po::bool_switch(&do_not_convert_spaces)->default_value(false), "Do not convert all spaces in reference ids to underscores [Default: converting is on]")
+ ("trim-after-space", po::bool_switch(&trim_ids)->default_value(false), "Trim all reference ids after first space [Default: false]");
+ return options;
+}
+
+void BuildIndexArgumentParser::init_help(po::options_description visible_options) {
+
+ std::stringstream help_message;
+
+ help_message << "Copyright (c) 2015-2017, Martin S. Lindner and the HiLive contributors. See CONTRIBUTORS for more info." << std::endl;
+ help_message << "All rights reserved" << std::endl << std::endl;
+ help_message << "HiLive is open-source software. Check with --license for details." << std::endl << std::endl;
+ help_message << "Usage: " << std::endl << " hilive-build INPUT KMER_WEIGHT [options]" << std::endl << std::endl;
+ help_message << "Required:" << std::endl;
+ help_message << " INPUT Reference genomes in (multi-) FASTA format." << std::endl;
+ help_message << " KMER_WEIGHT Number of non-gap positions in a k-mer (For ungapped k-mers this is the k-mer size)." << std::endl;
+
+ help_message << visible_options;
+
+ help = help_message.str();
+}
+
+bool BuildIndexArgumentParser::set_positional_variables(po::variables_map vm) {
+
+ // Name of the input fasta file
+ fasta_name = vm["INPUT"].as<std::string>();
+
+
+ // User-defined k-mer weight
+ kmer_weight = vm["KMER_WEIGHT"].as<uint16_t>();
+
+ // Check if input file exists
+ if ( !file_exists(fasta_name) ){
+ std::cerr << "Input error: Could not find input file " << fasta_name << std::endl;
+ return false;
+ }
+
+ // Check for maximal k-mer size
+ CountType maxKmerWeight = sizeof(HashIntoType)*4;
+ if ( kmer_weight > maxKmerWeight ) {
+ std::cerr << "K-mer weight is too high. Maximal k-mer weight is " << maxKmerWeight << "." << std::endl;
+ return false;
+ }
+
+ return true;
+}
+
+bool BuildIndexArgumentParser::set_build_variables(po::variables_map vm) {
+
+ if ( vm.count("gap-positions") ) {
+ gap_positions = vm["gap-positions"].as< std::vector <unsigned> >();
+ }
+
+ // Init the k-mer structure
+ if ( ! globalAlignmentSettings.set_kmer(kmer_weight, gap_positions) ) {
+ return false;
+ }
+
+ // Name of the ouput .kix file
+ if (vm.count("outfile")) {
+ index_name = vm["outfile"].as<std::string>();
+ } else {
+ index_name = fasta_name + std::string(".kix");
+ }
+
+ return true;
+}
+
+void BuildIndexArgumentParser::report() {
+
+ std::cout << "K-mer weight: " << (uint16_t) globalAlignmentSettings.get_kmer_weight() << std::endl;
+ std::cout << "K-mer span: " << (uint16_t) globalAlignmentSettings.get_kmer_span() << std::endl;
+ std::cout << "K-mer gap positions: ";
+
+ if ( globalAlignmentSettings.get_kmer_gaps().size() > 0 ) {
+ for ( auto pos : globalAlignmentSettings.get_kmer_gaps() ) {
+ if ( pos != *(globalAlignmentSettings.get_kmer_gaps().begin()) )
+ std::cout << ",";
+ std::cout << (uint16_t) pos;
+ }
+ std::cout << std::endl;
+ } else {
+ std::cout << "-" << std::endl;
+ }
+ std::cout << std::endl;
+
+}
+
+int BuildIndexArgumentParser::parseCommandLineArguments() {
+
+ po::options_description gen_opt = general_options();
+ po::options_description pos_opt = positional_options();
+ po::options_description build_opt = build_options();
+
+ po::options_description cmdline_options;
+ cmdline_options.add(pos_opt).add(gen_opt).add(build_opt);
+
+ po::options_description visible_options;
+ visible_options.add(gen_opt).add(build_opt);
+
+ init_help(visible_options);
+
+ po::positional_options_description p;
+ p.add("INPUT", 1);
+ p.add("KMER_WEIGHT", 1);
+
+ po::variables_map vm;
+ try {
+ // parse arguments
+ po::store(po::command_line_parser(argc, argv).
+ options(cmdline_options).positional(p).run(), vm);
+ // first check if -h or --help was called
+ if (vm.count("help")) {
+ printHelp();
+ return 1;
+ }
+ // first check if --license was called
+ if (vm.count("license")) {
+ printLicense();
+ return 1;
+ }
+
+ // then check arguments
+ po::notify(vm);
+ }
+ catch ( po::required_option& e ) {
+ std::cerr << "Missing Parameter: " << e.what() << std::endl;
+ return -1;
+ }
+ catch( po::error& e) {
+ std::cerr << "Error while parsing command line options: " << e.what() << std::endl;
+ return -1;
+ }
+
+ if ( !set_positional_variables(vm) ) {
+ return -1;
+ }
-int parseCommandLineArguments(AlignmentSettings & settings, std::string license, int argc, char const ** argv)
-{
- po::options_description general("General");
- general.add_options()
- ("help,h", "Print this help message and exit")
- ("license", "Print licensing information and exit");
-
- po::options_description parameters("Parameters");
- parameters.add_options()
- ("BC_DIR", po::value<std::string>(&settings.root)->required(), "Illumina BaseCalls directory")
- ("INDEX", po::value<std::string>(&settings.index_fname)->required(), "Path to k-mer index")
- ("CYCLES", po::value<CountType>(&settings.rlen)->required(), "Number of cycles")
- ("OUTDIR", po::value<std::string>(&settings.out_dir), "Directory to store sam files in [Default: temporary or BaseCalls directory");
-
- po::options_description io_settings("IO settings");
- io_settings.add_options()
- ("temp", po::value<std::string>(&settings.temp_dir)->default_value(""), "Temporary directory for the alignment files [Default: use BaseCalls directory]")
- //("bam,B", po::bool_switch(&settings.write_bam)->default_value(false), "Create BAM files instead of SAM files [Default: false]")
- ("keep-files,k", po::bool_switch(&settings.keep_aln_files)->default_value(false), "Keep intermediate alignment files [Default: false]")
- ("lanes,l", po::value< std::vector<uint16_t> >()->multitoken()->composing(), "Select lane [Default: all lanes]")
- ("tiles,t", po::value< std::vector<uint16_t> >()->multitoken()->composing(), "Select tile numbers [Default: all tiles]")
- ("barcodes,b", po::value< std::vector<std::string> >()->multitoken()->composing(), "Enumerate barcodes (must have same length) for demultiplexing, i.e. -b AGGATC -b CCCTTT [Default: no demultiplexing]");
-
- po::options_description alignment("Alignment settings");
- alignment.add_options()
- ("min-errors,e", po::value<CountType>(&settings.min_errors)->default_value(2), "Number of errors tolerated in read alignment [Default: 2]")
- ("all-best-hit,H", po::bool_switch()->default_value(false), "Report all of the best alignments for each read")
- ("any-best-hit", po::bool_switch(), "Report one of the best alignments for each read (default)")
- ("all-best-n-scores,N", po::value<CountType>(&settings.best_n), "Report all alignments of the N best alignment scores for each read")
- ("all-hits,A", po::bool_switch()->default_value(false), "Report all valid alignments for each read")
- ("disable-ohw-filter", po::bool_switch(&settings.discard_ohw)->default_value(true), "disable the One-Hit Wonder filter [Default: false]")
- ("start-ohw", po::value<CountType>(&settings.start_ohw)->default_value(K_HiLive+5), "First cycle to apply One-Hit Wonder filter [Default: K+5]")
- ("window,w", po::value<DiffType>(&settings.window)->default_value(5), "Set the window size to search for alignment extension, i.e. maximum total insertion/deletion size [Default: 5]")
- ("min-quality", po::value<CountType>(&settings.min_qual)->default_value(1), "Minimum allowed basecall quality [Default: 1]");
-
- po::options_description technical("Technical settings");
- technical.add_options()
- ("block-size", po::value<uint64_t>()->default_value(64*1024*1024), "Block size for the alignment input/output stream in Bytes. Use -K or -M to specify in Kilobytes or Megabytes")
- (",K", po::bool_switch()->default_value(false), "Interpret the block-size argument as Kilobytes instead of Bytes")
- (",M", po::bool_switch()->default_value(false), "Interpret the block-size argument as Megabytes instead of Bytes")
- ("compression,c", po::value<uint8_t>(&settings.compression_format)->default_value(2), "Compress alignment files. 0: no compression (default) 1: Deflate (smaller) 2: LZ4 (faster)")
- ("num-threads,n", po::value<int>(), "Number of threads to spawn [Default: all available]");
+ if ( !set_build_variables(vm) ) {
+ return -1;
+ }
+ report();
+
+ return 0;
+}
+
+//|||||||||||||||||||||||||||||||||||||||||||||||||||||||
+//-----HiLiveArgumentParser------------------------------
+//|||||||||||||||||||||||||||||||||||||||||||||||||||||||
+
+po::options_description HiLiveArgumentParser::general_options() {
+ po::options_description general("General");
+ general.add_options()
+ ("help,h", "Print this help message and exit")
+ ("license", "Print licensing information and exit")
+ ("settings,s", po::value<std::string>(), "Load settings from file. If command line arguments are given additionally, they are prefered.")
+ ("runinfo", po::value<std::string>(), "Path to runInfo.xml for parsing read and index lengths [Default (if activated): BC_DIR/../../RunInfo.xml]");
+
+ return general;
+}
+
+po::options_description HiLiveArgumentParser::positional_options() {
+ po::options_description parameters("Parameters");
+ parameters.add_options()
+ ("BC_DIR", po::value<std::string>(), "Illumina BaseCalls directory")
+ ("INDEX", po::value<std::string>(), "Path to k-mer index")
+ ("CYCLES", po::value<CountType>(), "Number of cycles")
+ ("OUTDIR", po::value<std::string>(), "Directory to store sam files in [Default: temporary or BaseCalls directory");
+
+ return parameters;
+}
+
+po::options_description HiLiveArgumentParser::io_options() {
+ po::options_description io_settings("IO settings");
+ io_settings.add_options()
+ ("temp", po::value<std::string>(), "Temporary directory for the alignment files [Default: use BaseCalls directory]")
+ ("bam,B", po::bool_switch(), "Create BAM files instead of SAM files [Default: false]")
+ ("output-cycles,O", po::value<std::vector<CountType>>()->multitoken()->composing(), "Cycles for alignment output. The respective temporary files are kept. [Default: last cycle]")
+ ("extended-cigar", po::bool_switch(), "Activate extended CIGAR format (= and X instead of only M) in output files [Default: false]")
+ ("keep-files,k", po::bool_switch(), "Keep intermediate alignment files [Default: false]")
+ ("min-as-ratio", po::value<float>(), "Minimum alignment score (relative to the current read length) for alignments to be reported (0-1) [Default: 0 - Report all alignments]")
+ ("lanes,l", po::value< std::vector<uint16_t> >()->multitoken()->composing(), "Select lane [Default: all lanes]")
+ ("tiles,t", po::value< std::vector<uint16_t> >()->multitoken()->composing(), "Select tile numbers [Default: all tiles]")
+ ("reads,r", po::value< std::vector<std::string> >()->multitoken()->composing(), "Enumerate read lengths and type. Example: -r 101R 8B 8B 101R equals paired-end sequencing with 2x101bp reads and 2x8bp barcodes. Overwrites information of runInfo.xml. [Default: single end reads without barcodes]");
+ return io_settings;
+}
+
+po::options_description HiLiveArgumentParser::alignment_options() {
+ po::options_description alignment("Alignment settings");
+ alignment.add_options()
+ ("min-errors,e", po::value<CountType>(), "Number of errors tolerated in read alignment [Default: 2]")
+ ("mode,m", po::value<std::string>(), "Alignment mode. [ALL|A]: Report all alignments; [BESTN#|N#]: Report alignments of the best # scores; "
+ "[ALLBEST|H]: Report all alignments with the best score (similar to N1); [ANYBEST|B]: Report one best alignment (default)")
+ ("disable-ohw-filter", po::bool_switch(), "Disable the One-Hit Wonder filter [Default: false]")
+ ("start-ohw", po::value<CountType>(), "First cycle to apply One-Hit Wonder filter [Default: 20]")
+ ("window,w", po::value<DiffType>(), "Set the window size to search for alignment extension, i.e. maximum total insertion/deletion size [Default: 5]")
+ ("min-quality", po::value<CountType>(), "Minimum allowed basecall quality [Default: 1]")
+ ("barcodes,b", po::value< std::vector<std::string> >()->multitoken()->composing(), "Enumerate barcodes (must have same length) for demultiplexing, e.g. -b AGGATC -b CCCTTT [Default: no demultiplexing]")
+ ("barcode-errors,E", po::value< std::vector<uint16_t> >()->multitoken()->composing(), "Enumerate the number of tolerated errors (only SNPs) for each barcode fragment, e.g. -E 2 2 [Default: 1 per fragment]")
+ ("keep-all-barcodes", po::bool_switch()->default_value(false), "Align and output all barcodes [Default: false]");
+ return alignment;
+}
+
+po::options_description HiLiveArgumentParser::technical_options() {
+ po::options_description technical("Technical settings");
+ technical.add_options()
+ ("block-size", po::value<std::string>(), "Block size for the alignment input/output stream in Bytes. Append 'K' or 'M' to specify in Kilobytes or Megabytes, respectively (e.g. '--block-size 64M' for 64 Megabytes)")
+ ("compression,c", po::value<uint16_t>(), "Compress alignment files. 0: no compression 1: Deflate (smaller) 2: LZ4 (faster; default)")
+ ("num-threads,n", po::value<CountType>(), "Number of threads to spawn [Default: all available]");
+ return technical;
+}
+
+void HiLiveArgumentParser::init_help(po::options_description visible_options) {
+
+ std::stringstream help_message;
+
+ help_message << "Copyright (c) 2015-2017, Martin S. Lindner and the HiLive contributors. See CONTRIBUTORS for more info." << std::endl;
+ help_message << "All rights reserved" << std::endl << std::endl;
+ help_message << "HiLive is open-source software. Check with --license for details." << std::endl << std::endl;
+
+ help_message << "Usage: " << std::endl << " hilive BC_DIR INDEX CYCLES OUTDIR [options]" << std::endl << std::endl;
+ help_message << "Required:" << std::endl;
+ help_message << " BC_DIR Illumina BaseCalls directory of the sequencing run to analyze" << std::endl;
+ help_message << " INDEX Path to k-mer index file (*.kix)" << std::endl;
+ help_message << " CYCLES Total number of sequencing cycles" << std::endl;
+ help_message << " OUTDIR Output directory" << std::endl;
+
+ help_message << visible_options;
+
+ help = help_message.str();
+}
+
+bool HiLiveArgumentParser::checkPaths() {
+ if (!file_exists(globalAlignmentSettings.get_index_fname())){
+ std::cerr << "Input error: Could not find k-mer index file " << globalAlignmentSettings.get_index_fname() << std::endl;
+ return false;
+ }
+
+ std::size_t found = globalAlignmentSettings.get_root().find("BaseCalls");
+ if (!(found != std::string::npos && found >= globalAlignmentSettings.get_root().size()-10)) {
+ std::cerr << "Warning: BaseCalls directory seems to be invalid: " << globalAlignmentSettings.get_root() << std::endl;
+ }
+
+ if (!is_directory(globalAlignmentSettings.get_root())){
+ std::cerr << "Input error: Could not find BaseCalls directory " << globalAlignmentSettings.get_root() << std::endl;
+ return false;
+ }
+
+ for ( uint16_t ln : globalAlignmentSettings.get_lanes() ) {
+ std::string ln_dir = globalAlignmentSettings.get_root();
+ if ( ln < 10 )
+ ln_dir += "/L00";
+ else if ( ln < 100 )
+ ln_dir += "/L0";
+ else
+ ln_dir += "/L";
+ ln_dir += std::to_string(ln);
+ if (!is_directory(ln_dir)){
+ std::cerr << "Input error: Could not find location of Lane " << ln << ": " << ln_dir << std::endl;
+ return false;
+ }
+ }
+
+ if ( !is_directory(globalAlignmentSettings.get_temp_dir())) {
+ boost::filesystem::create_directories(globalAlignmentSettings.get_temp_dir());
+ }
+
+ if ( !is_directory(globalAlignmentSettings.get_out_dir())) {
+ boost::filesystem::create_directories(globalAlignmentSettings.get_out_dir());
+ }
+
+ return true;
+}
+
+void HiLiveArgumentParser::report() {
+ std::cout << "Running HiLive with " << globalAlignmentSettings.get_num_threads() << " thread(s)." << std::endl;
+ std::cout << "BaseCalls directory: " << globalAlignmentSettings.get_root() << std::endl;
+ if (globalAlignmentSettings.get_temp_dir() != "") {
+ std::cout << "Temporary directory: " << globalAlignmentSettings.get_temp_dir() << std::endl;
+ }
+ if (!globalAlignmentSettings.get_write_bam())
+ std::cout << "SAM output directory: " << globalAlignmentSettings.get_out_dir() << std::endl;
+ else
+ std::cout << "BAM output directory: " << globalAlignmentSettings.get_out_dir() << std::endl;
+ std::cout << "Lanes: ";
+ for ( uint16_t ln : globalAlignmentSettings.get_lanes() )
+ std::cout << ln << " ";
+ std::cout << std::endl;
+ std::cout << "K-mer index: " << globalAlignmentSettings.get_index_fname() << std::endl;
+ std::cout << "Read lengths: ";
+ std::string barcode_suffix;
+ for ( uint16_t read = 0; read != globalAlignmentSettings.get_seqs().size(); read ++) {
+ std::cout << globalAlignmentSettings.getSeqById(read).length;
+ barcode_suffix = globalAlignmentSettings.getSeqById(read).isBarcode() ? "B" : "R";
+ std::cout << barcode_suffix << " ";
+ }
+ std::cout << std::endl;
+ std::cout << "Mapping error: " << globalAlignmentSettings.get_min_errors() << std::endl;
+ if (globalAlignmentSettings.get_any_best_hit_mode())
+ std::cout << "Mapping mode: Any-Best-Hit-Mode" << std::endl;
+ else if (globalAlignmentSettings.get_all_best_hit_mode())
+ std::cout << "Mapping mode: All-Best-Hit-Mode" << std::endl;
+ else if (globalAlignmentSettings.get_all_best_n_scores_mode())
+ std::cout << "Mapping mode: All-Best-N-Scores-Mode with N=" << globalAlignmentSettings.get_best_n() << std::endl;
+ else
+ std::cout << "Mapping mode: All-Hits-Mode" << std::endl;
+ std::cout << std::endl;
+}
+
+bool HiLiveArgumentParser::parseRunInfo(po::variables_map vm) {
+
+ if ( ! vm.count("runinfo"))
+ return false;
+
+ boost::property_tree::ptree tree;
+ read_xml(tree, vm["runinfo"].as<std::string>());
+
+ using boost::property_tree::ptree;
+
+ if (!tree.empty() && tree.count("RunInfo")!=0) {
+ ptree ptree_RunInfo = tree.get_child("RunInfo");
+
+ if (ptree_RunInfo.count("Run")!=0) {
+ ptree ptree_Run = ptree_RunInfo.get_child("Run");
+
+ if (ptree_Run.count("Reads")!=0) {
+ ptree ptree_Reads = ptree_Run.get_child("Reads");
+ unsigned lenSum = 0;
+ std::string sequences = "";
+ std::vector<SequenceElement> tempSeqs;
+ CountType num_cycles = 0;
+ unsigned mates = 0;
+ for (const auto &read : ptree_Reads) {
+ sequences += read.second.get<std::string>("<xmlattr>.NumCycles");
+ num_cycles += read.second.get<unsigned>("<xmlattr>.NumCycles");
+ sequences += read.second.get<std::string>("<xmlattr>.IsIndexedRead") == "N" ? "R " : "B ";
+ tempSeqs.push_back(SequenceElement(tempSeqs.size(), (read.second.get<std::string>("<xmlattr>.IsIndexedRead") == "N") ? ++mates : 0, read.second.get<unsigned>("<xmlattr>.NumCycles")));
+ lenSum += read.second.get<unsigned>("<xmlattr>.NumCycles");
+ }
+ if ( sequences != "" )
+ runInfo_settings.put("settings.sequences", sequences);
+ runInfo_settings.put("settings.cycles", num_cycles);
+
+ if (ptree_Run.count("FlowcellLayout")!=0) {
+ ptree ptree_FlowcellLayout = ptree_Run.get_child("FlowcellLayout");
+
+ std::vector<uint16_t> lanes_vec(ptree_FlowcellLayout.get<unsigned>("<xmlattr>.LaneCount"));
+ std::iota(lanes_vec.begin(), lanes_vec.end(), 1);
+ std::string lanes = "";
+
+ for ( auto l : lanes_vec )
+ lanes += std::to_string(l) + " ";
+ runInfo_settings.put("settings.lanes", lanes);
+
+ std::vector<uint16_t> tiles_vec;
+ std::string tiles;
+
+ for (uint16_t l = 1; l <= ptree_FlowcellLayout.get<unsigned>("<xmlattr>.SurfaceCount"); l++)
+ for (uint16_t s = 1; s <= ptree_FlowcellLayout.get<unsigned>("<xmlattr>.SwathCount"); s++)
+ for (uint16_t t = 1; t <= ptree_FlowcellLayout.get<unsigned>("<xmlattr>.TileCount"); t++)
+ tiles += std::to_string(l*1000 + s*100 + t) + " ";
+ }
+ }
+ }
+ }
+ return true;
+}
+
+int HiLiveArgumentParser::parseCommandLineArguments() {
+
+ this->set_required_parameters();
+
+ // Init general options
+ po::options_description gen_opt = general_options();
+ po::variables_map vm;
+
+ // Init all other options
+ po::options_description pos_opt = positional_options();
+ po::options_description io_opt = io_options();
+ po::options_description align_opt = alignment_options();
+ po::options_description tech_opt = technical_options();
+
+ // All command line options
po::options_description cmdline_options;
- cmdline_options.add(general).add(parameters).add(io_settings).add(alignment).add(technical);
+ cmdline_options.add(gen_opt).add(pos_opt).add(io_opt).add(align_opt).add(tech_opt);
+ // Options visible in the help
po::options_description visible_options;
- visible_options.add(general).add(io_settings).add(alignment).add(technical);
-
- std::stringstream help_message;
- help_message << "HiLive v"<< HiLive_VERSION_MAJOR << "." << HiLive_VERSION_MINOR << " - Realtime Alignment of Illumina Reads" << std::endl;
- help_message << "Copyright (c) 2015, Martin S. Lindner" << std::endl;
- help_message << "HiLive is open-source software. Check with --license for details." << std::endl << std::endl;
- help_message << "Fixed k-mer size: " << K_HiLive << std::endl << std::endl;
- help_message << "Usage: " << std::string(argv[0]) << " [options] BC_DIR INDEX CYCLES OUTDIR" << std::endl;
- help_message << " BC_DIR Illumina BaseCalls directory of the sequencing run to analyze" << std::endl;
- help_message << " INDEX Path to k-mer index file (*.kix)" << std::endl;
- help_message << " CYCLES Total number of cycles for read 1" << std::endl;
- help_message << " OUTDIR Directory to store sam files in" << std::endl;
- help_message << visible_options << std::endl;
+ visible_options.add(gen_opt).add(io_opt).add(align_opt).add(tech_opt);
- po::positional_options_description p;
- p.add("BC_DIR", 1);
- p.add("INDEX", 1);
- p.add("CYCLES", 1);
- p.add("OUTDIR", 1);
-
-
- // parse the arguments
+ init_help(visible_options);
- po::variables_map vm;
+ // First parameter iteration for general options (includes help, license and input settings file.
try {
- // parse arguments
- po::store(po::command_line_parser(argc, argv).options(cmdline_options).positional(p).run(), vm);
+ po::store(po::command_line_parser(argc, argv).options(gen_opt).allow_unregistered().run(), vm);
+
// first check if -h or --help was called
if (vm.count("help")) {
- std::cout << help_message.str();
- return 1;
+ printHelp();
+ return 1;
}
+
// first check if --license was called
if (vm.count("license")) {
- std::cout << license << std::endl;
- return 1;
+ printLicense();
+ return 1;
}
+ vm.notify();
+
+ // Load input settings if exist
+ if ( vm.count("settings")) {
+ if ( ! read_xml(input_settings, vm["settings"].as<std::string>()) ) {
+ std::cerr << "Input settings file not found: " << vm["settings"].as<std::string>() << std::endl;
+ return -1;
+ }
+ } else if ( isRequired("settings") ) {
+ throw po::required_option("settings");
+ }
+
+ if ( vm.count("runinfo") ) {
+ if ( ! parseRunInfo(vm) ) {
+ std::cerr << "Error while parsing Run Info file: " << vm["runinfo"].as<std::string>() << std::endl;
+ return -1;
+
+ }
+ } else if ( isRequired("runinfo") ) {
+ throw po::required_option("runinfo");
+ }
+ }
+ catch( po::error& e) {
+ std::cerr << "Error while parsing command line options: " << std::endl << e.what() << std::endl;
+ return -1;
+ }
+
+ po::positional_options_description p;
+ p.add("BC_DIR", 1);
+ p.add("INDEX", 1);
+ p.add("CYCLES", 1);
+ p.add("OUTDIR", 1);
+
+ // Parse all command line arguments to cmd_settings
+ try {
+ // parse arguments
+ po::store(po::command_line_parser(argc, argv).options(cmdline_options).positional(p).run(), cmd_settings);
+
// then check arguments
- po::notify(vm);
+ po::notify(cmd_settings);
}
- catch ( boost::program_options::required_option& e ) {
- std::cerr << "Missing Parameter: " << e.what() << std::endl << std::endl;
- std::cout << help_message.str();
+ catch ( po::required_option& e ) {
+ std::cerr << "Missing Parameter: " << e.what() << std::endl;
return -1;
}
- catch( boost::program_options::error& e) {
- std::cerr << "Error while parsing command line options: " << e.what() << std::endl << std::endl;
- std::cout << help_message.str();
+ catch( po::error& e) {
+ std::cerr << "Error while parsing command line options: " << e.what() << std::endl;
return -1;
}
- if (vm.count("OUTDIR"))
- settings.out_dir = vm["OUTDIR"].as<std::string>();
- else {
- if (settings.temp_dir == "")
- settings.out_dir = settings.root;
- else
- settings.out_dir = settings.temp_dir;
+ // Set all options from command line and input files
+ if ( !set_options() ) {
+ return -1;
}
- if (vm.count("lanes"))
- settings.lanes = vm["lanes"].as< std::vector<uint16_t> >();
- else
- settings.lanes = all_lanes();
+ if ( !checkPaths() ) {
+ return -1;
+ }
- if (vm.count("tiles"))
- settings.tiles = vm["tiles"].as< std::vector<uint16_t> >();
- else
- settings.tiles = all_tiles();
+ // Report the basic settings
+ report();
- settings.barcodeVector.clear();
- settings.seqlen = settings.rlen;
- if (vm.count("barcodes")) {
- settings.barcodeVector = vm["barcodes"].as< std::vector<std::string> >();
- settings.seqlen = settings.rlen - settings.barcodeVector[0].size();
- }
+ return 0;
+}
- if (vm["all-hits"].as<bool>()) {
- // all hits: disable other modes
- settings.any_best_hit_mode = false;
- settings.all_best_hit_mode = false;
- settings.all_best_n_scores_mode = false;
- } else if (vm["all-best-hit"].as<bool>()) {
- // enable all-best-hit mode and disable others
- settings.any_best_hit_mode = false;
- settings.all_best_hit_mode = true;
- settings.all_best_n_scores_mode = false;
- } else if (vm.count("all-best-n-scores")) {
- // enable any-best-n mode and get parameter
- settings.any_best_hit_mode = false;
- settings.all_best_hit_mode = false;
- settings.all_best_n_scores_mode = true;
- settings.best_n = vm["all-best-n-scores"].as<CountType>();
- } else { // the default behaviour
- // enable any-best-hit mode and disable others
- settings.any_best_hit_mode = true;
- settings.all_best_hit_mode = false;
- settings.all_best_n_scores_mode = false;
- }
-
- if (vm["-M"].as<bool>())
- settings.block_size = vm["block-size"].as<uint64_t>()*1024*1024;
- else if (vm["-K"].as<bool>())
- settings.block_size = vm["block-size"].as<uint64_t>()*1024;
- else
- settings.block_size = vm["block-size"].as<uint64_t>();
-
- if (vm.count("num-threads"))
- settings.num_threads = vm["num-threads"].as<int>();
- else { // the default case, meaning as much as physically useful
- uint32_t n_cpu = std::thread::hardware_concurrency();
- if (n_cpu > 1)
- settings.num_threads = n_cpu;
- else
- settings.num_threads = 1;
- }
+bool HiLiveArgumentParser::set_options() {
+ try {
- // check paths and file names
- if (!file_exists(settings.index_fname)){
- std::cerr << "Input error: Could not find k-mer index file " << settings.index_fname << std::endl;
- return -1;
- }
+ // Set positional arguments
+ set_option<std::string>("BC_DIR", "settings.paths.root", "", &AlignmentSettings::set_root);
+ set_option<std::string>("INDEX", "settings.paths.index", "", &AlignmentSettings::set_index_fname);
+ set_option<CountType>("CYCLES", "settings.cycles", 0, &AlignmentSettings::set_cycles);
+ set_option<std::string>("OUTDIR", "settings.paths.out_dir", "", &AlignmentSettings::set_out_dir);
- std::size_t found = settings.root.find("BaseCalls");
- if (!(found != std::string::npos && found >= settings.root.size()-10)) {
- std::cerr << "Warning: BaseCalls directory seems to be invalid: " << settings.root << std::endl;
- }
+ // Set I/O options
+ set_option<std::string>("temp", "settings.paths.temp_dir", "", &AlignmentSettings::set_temp_dir);
+ set_option<bool>("bam", "settings.out.bam", false, &AlignmentSettings::set_write_bam);
- if (!is_directory(settings.root)){
- std::cerr << "Input error: Could not find BaseCalls directory " << settings.root << std::endl;
- return -1;
- }
+ std::vector<CountType> output_cycles = {globalAlignmentSettings.get_cycles()};
+ set_option<std::vector<CountType>>("output-cycles", "settings.out.cycles", output_cycles, &AlignmentSettings::set_output_cycles);
+ set_option<bool>("extended-cigar", "settings.out.extended_cigar", false, &AlignmentSettings::set_extended_cigar);
+ set_option<bool>("keep-files", "settings.technical.keep_aln_files", false, &AlignmentSettings::set_keep_aln_files);
+ set_option<float>("min-as-ratio", "settings.out.min_as_ratio", 0.0f, &AlignmentSettings::set_min_as_ratio);
+ set_option<std::vector<uint16_t>>("lanes", "settings.lanes", all_lanes(), &AlignmentSettings::set_lanes);
+ set_option<std::vector<uint16_t>>("tiles", "settings.tiles", all_tiles(), &AlignmentSettings::set_tiles);
- for ( uint16_t ln : settings.lanes ) {
- std::string ln_dir = settings.root;
- if ( ln < 10 )
- ln_dir += "/L00";
- else if ( ln < 100 )
- ln_dir += "/L0";
- else
- ln_dir += "/L";
- ln_dir += std::to_string(ln);
- if (!is_directory(ln_dir)){
- std::cerr << "Input error: Could not find location of Lane " << ln << ": " << ln_dir << std::endl;
- return -1;
- }
- }
+ // Set alignment options
+ std::vector<std::string> default_read_structure;
+ default_read_structure.push_back(std::to_string(globalAlignmentSettings.get_cycles()) + "R");
+ set_option<std::vector<std::string>>("reads", "settings.sequences", default_read_structure, &AlignmentSettings::set_read_structure);
+ set_option<CountType>("min-errors", "settings.min_errors", 2, &AlignmentSettings::set_min_errors);
+ set_option<std::string>("mode", "settings.mode", "ANYBEST", &AlignmentSettings::set_mode);
+ set_option<bool>("disable-ohw-filter", "settings.align.discard_ohw", false, &AlignmentSettings::disable_ohw);
+ set_option<CountType>("start-ohw", "settings.align.start_ohw", 20, &AlignmentSettings::set_start_ohw);
+ set_option<DiffType>("window", "settings.align.window", 5, &AlignmentSettings::set_window);
+ set_option<CountType>("min-quality", "settings.align.min_qual", 1, &AlignmentSettings::set_min_qual);
- // Report the basic settings
- std::cout << "Running HiLive with " << settings.num_threads << " thread(s)." << std::endl;
- std::cout << "BaseCalls directory: " << settings.root << std::endl;
- if (settings.temp_dir != "") {
- std::cout << "Temporary directory: " << settings.temp_dir << std::endl;
- }
- //if (!settings.write_bam)
- std::cout << "SAM output directory: " << settings.out_dir << std::endl;
- //else
- //std::cout << "BAM output directory: " << settings.out_dir << std::endl;
- std::cout << "Lanes: ";
- for ( uint16_t ln : settings.lanes )
- std::cout << ln << " ";
- std::cout << std::endl;
- std::cout << "K-mer index: " << settings.index_fname << std::endl;
- std::cout << "Read length: " << settings.rlen << std::endl;
- std::cout << "Mapping error: " << settings.min_errors << std::endl;
- if (settings.any_best_hit_mode)
- std::cout << "Mapping mode: Any-Best-Hit-Mode" << std::endl;
- else if (settings.all_best_hit_mode)
- std::cout << "Mapping mode: All-Best-Hit-Mode" << std::endl;
- else if (settings.all_best_n_scores_mode)
- std::cout << "Mapping mode: All-Best-N-Scores-Mode with N=" << settings.best_n << std::endl;
- else
- std::cout << "Mapping mode: All-Hits-Mode" << std::endl;
- std::cout << std::endl;
+ std::vector<std::string> barcode_sequences_default;
+ set_option<std::vector<std::string>>("barcodes", "settings.barcodes.sequences", barcode_sequences_default, &AlignmentSettings::set_barcodes);
- return 0;
+ std::vector<CountType> barcode_errors_default = {2};
+ set_option<std::vector<CountType>>("barcode-errors", "settings.barcodes.errors", barcode_errors_default, &AlignmentSettings::set_barcode_errors);
+
+ set_option<bool>("keep-all-barcodes", "settings.barcodes.keep_all", false, &AlignmentSettings::set_keep_all_barcodes);
+
+ // Set technical options
+ set_option<std::string>("block-size", "settings.technical.block_size", "64M", &AlignmentSettings::set_block_size);
+ set_option<uint16_t>("compression", "settings.technical.compression_format", 2, &AlignmentSettings::set_compression_format);
+
+ CountType n_cpu = std::thread::hardware_concurrency();
+ CountType n_threads_default = 1;
+ if (n_cpu > 1)
+ n_threads_default = std::min( n_cpu, CountType( globalAlignmentSettings.get_lanes().size() * globalAlignmentSettings.get_tiles().size() ) ) ;
+ set_option<CountType>("num-threads", "settings.technical.num_threads", n_threads_default, &AlignmentSettings::set_num_threads);
+
+ } catch ( std::exception & ex ) {
+ std::cerr << "Error while parsing options: " << std::endl << ex.what() << std::endl;
+ return false;
+ }
+ return true;
+}
+
+//|||||||||||||||||||||||||||||||||||||||||||||||||||||||
+//-----HiLiveOutArgumentParser--------------------------
+//|||||||||||||||||||||||||||||||||||||||||||||||||||||||
+
+void HiLiveOutArgumentParser::init_help(po::options_description visible_options) {
+
+ std::stringstream help_message;
+
+ help_message << "Copyright (c) 2015-2017, Martin S. Lindner and the HiLive contributors. See CONTRIBUTORS for more info." << std::endl;
+ help_message << "All rights reserved" << std::endl << std::endl;
+ help_message << "HiLive is open-source software. Check with --license for details." << std::endl << std::endl;
+
+ help_message << "Usage: " << std::endl << " hilive-out --settings /path/to/settings/file [options]" << std::endl << std::endl;
+ help_message << "Required:" << std::endl;
+ help_message << " settings Path to a HiLive settings file (by default, the file is in the temp directory of the respective run)" << std::endl;
+ help_message << std::endl << "All parameters can be set as for the HiLive main program." << std::endl;
+ help_message << "By default, only output files for the last cycle are produced." << std::endl;
+ help_message << "Use the --output-cycles parameter to declare different cycle numbers (will only work if --keep-files was activated for the respective HiLive run)" << std::endl;
+
+ help_message << visible_options;
+
+ help = help_message.str();
+}
+
+void HiLiveOutArgumentParser::report() {
+ if (globalAlignmentSettings.get_temp_dir() != "") {
+ std::cout << "Temporary directory: " << globalAlignmentSettings.get_temp_dir() << std::endl;
+ }
+ if (!globalAlignmentSettings.get_write_bam())
+ std::cout << "SAM output directory: " << globalAlignmentSettings.get_out_dir() << std::endl;
+ else
+ std::cout << "BAM output directory: " << globalAlignmentSettings.get_out_dir() << std::endl;
+ std::cout << "Lanes: ";
+ for ( uint16_t ln : globalAlignmentSettings.get_lanes() )
+ std::cout << ln << " ";
+ std::cout << std::endl;
+ std::cout << "K-mer index: " << globalAlignmentSettings.get_index_fname() << std::endl;
+ std::cout << "Total Read lengths: ";
+ std::string barcode_suffix;
+ for ( uint16_t read = 0; read != globalAlignmentSettings.get_seqs().size(); read ++) {
+ std::cout << globalAlignmentSettings.getSeqById(read).length;
+ barcode_suffix = globalAlignmentSettings.getSeqById(read).isBarcode() ? "B" : "R";
+ std::cout << barcode_suffix << " ";
+ }
+ std::cout << std::endl;
+ std::cout << "Mapping error: " << globalAlignmentSettings.get_min_errors() << std::endl;
+ if (globalAlignmentSettings.get_any_best_hit_mode())
+ std::cout << "Mapping mode: Any-Best-Hit-Mode" << std::endl;
+ else if (globalAlignmentSettings.get_all_best_hit_mode())
+ std::cout << "Mapping mode: All-Best-Hit-Mode" << std::endl;
+ else if (globalAlignmentSettings.get_all_best_n_scores_mode())
+ std::cout << "Mapping mode: All-Best-N-Scores-Mode with N=" << globalAlignmentSettings.get_best_n() << std::endl;
+ else
+ std::cout << "Mapping mode: All-Hits-Mode" << std::endl;
+ std::cout << std::endl;
}
diff --git a/lib/argument_parser.h b/lib/argument_parser.h
index bbeca10..ff3212d 100644
--- a/lib/argument_parser.h
+++ b/lib/argument_parser.h
@@ -1,9 +1,482 @@
-#include <iostream>
#include <boost/program_options.hpp>
#include "headers.h"
#include "definitions.h"
-#include "kindex.h"
-#include "alnstream.h"
+#include "global_variables.h"
#include "parallel.h"
-int parseCommandLineArguments(AlignmentSettings & settings, std::string license, int argc, char const ** argv);
+namespace po = boost::program_options;
+
+
+/**
+ * Interface for the argument parsers of the different executables.
+ * @author Tobias Loka
+ */
+class ArgumentParser {
+
+protected:
+
+ /** Number of command line arguments. */
+ int argc;
+
+ /** List of command line arguments. */
+ char const ** argv;
+
+ /** Vector of required options. */
+ std::vector<std::string> required_options;
+
+ /** Map of command line arguments. */
+ po::variables_map cmd_settings;
+
+ /** Settings property tree obtained from the runInfo. */
+ boost::property_tree::ptree runInfo_settings;
+
+ /** Settings property tree obtained from an settings input file. */
+ boost::property_tree::ptree input_settings;
+
+ /** Bool describing whether an input settings file was specified. */
+ bool has_input_settings = false;
+
+ /** Program license. */
+ std::string license = "Copyright (c) 2015-2017, Martin S. Lindner and the HiLive contributors. See CONTRIBUTORS for more info.\n"
+ "All rights reserved.\n"
+ "\n"
+ "Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:\n"
+ "\n"
+ "1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.\n"
+ "\n"
+ "2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.\n"
+ "\n"
+ "3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.\n"
+ "\n"
+ "THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 'AS IS' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, [...]
+
+ /** Help output (differs for each executable). */
+ std::string help;
+
+ /**
+ * (virtual) function that is called from the executable's main() to parse the command line arguments.
+ * @return Error code ( 0: continue program; 1: exit with success; -1: exit with error)
+ * @author Tobias Loka
+ */
+ virtual int parseCommandLineArguments() = 0;
+
+ /**
+ * Print the license.
+ * @author Tobias Loka
+ */
+ void printLicense(){ std::cout << license << std::endl; };
+
+ /**
+ * Init the help text by using the options that should be visible to the user.
+ * @param visible_options The options that shall be visible for the user.
+ * @return help message as string
+ * @author Tobias Loka
+ */
+ virtual void init_help(po::options_description visible_options) = 0;
+
+ /**
+ * Print the license.
+ * @author Tobias Loka
+ */
+ void printHelp(){ std::cout << help << std::endl; };
+
+ /**
+ * Report the most important settings to the console.
+ * @author Martin Lindner
+ */
+ virtual void report() = 0;
+
+ virtual void set_required_parameters() {};
+
+ /**
+ * Check if an option is required.
+ * @param vm_key Key for the virtual map for the parameter to check (equals cmd line option name).
+ * @return true, if the option is required.
+ */
+ bool isRequired( std::string vm_key) {
+ if ( std::find(required_options.begin(), required_options.end(), vm_key) == required_options.end() )
+ return false;
+ return true;
+ }
+
+ /**
+ * Set an option in the globalAlignmentSettings.
+ * Thereby, the different input sources have the following priority:
+ * 1. Command line argument
+ * 2. RunInfo file
+ * 3. Input settings file
+ * @param vm_key Key for the command line's variables map.
+ * @param settings_key Key for the RunInfo and Input settings property tree (must be similar for both files)
+ * @param default_value A default value that is set if the variable is set in none of the input sources (this is not considered if "required" is true)
+ * @param function Function that is called to set the variable in globalAlignmentSettings. Must be of type void &AlignmentSettings::*(T).
+ * @param required true, if the option must be set by the user (from one of the input sources)
+ * @author Tobias Loka
+ */
+ template<class T> void set_option(std::string vm_key, std::string settings_key, T default_value, void (AlignmentSettings::*function)(T)) {
+ set_option_impl(vm_key, settings_key, default_value, function, isRequired(vm_key), static_cast<T*>(0));
+ }
+
+ /**
+ * General implementation of set_option(...).
+ */
+ template<class T> void set_option_impl(std::string vm_key, std::string settings_key, T default_value, void (AlignmentSettings::*function)(T), bool required, T*) {
+
+ T value = default_value;
+ bool was_set = false;
+
+ boost::optional<T> isv = input_settings.get_optional<T>(settings_key);
+ boost::optional<T> rsv = runInfo_settings.get_optional<T>(settings_key);
+
+
+ // User parameter -> first priority
+ if ( cmd_settings.count(vm_key) ) {
+ value = cmd_settings[vm_key].as<T>();
+ was_set = true;
+ }
+
+ // Settings file -> second priority
+ else if ( isv ) {
+ value = isv.get();
+ was_set = true;
+ }
+
+ // RunInfo file -> third priority
+ else if ( rsv ) {
+ value = rsv.get();
+ was_set = true;
+ }
+
+ // Throw exception if unset
+ if ( required && !was_set )
+ throw po::required_option(vm_key);
+
+ // Otherwise set value
+ auto binded_function = std::bind(function, &globalAlignmentSettings, std::placeholders::_1);
+ binded_function(value);
+
+ }
+
+ /**
+ * Overload of set_option_impl for bool data type.
+ */
+ void set_option_impl(std::string vm_key, std::string settings_key, bool default_value, void (AlignmentSettings::*function)(bool), bool required, bool *) {
+
+ bool value = default_value;
+ bool was_set = false;
+
+ boost::optional<bool> isv = input_settings.get_optional<bool>(settings_key);
+ boost::optional<bool> rsv = runInfo_settings.get_optional<bool>(settings_key);
+
+ if ( cmd_settings[vm_key].as<bool>() != default_value ) {
+ value = !default_value;
+ was_set = true;
+ }
+
+ else if ( isv && isv.get() != default_value ) {
+ value = !default_value;
+ was_set = true;
+ }
+
+ else if ( rsv && rsv.get() != default_value ) {
+ value = !default_value;
+ was_set = true;
+ }
+
+
+ if ( required && !was_set )
+ throw po::required_option(vm_key);
+
+ // Otherwise set value
+ auto binded_function = std::bind(function, &globalAlignmentSettings, std::placeholders::_1);
+ binded_function(value);
+
+ }
+
+ /** Overload of set_option_impl for std::vector data types. */
+ template<class T> void set_option_impl(std::string vm_key, std::string settings_key, std::vector<T> default_value, void (AlignmentSettings::*function)(std::vector<T>), bool required, std::vector<T> *) {
+
+ std::vector<T> value;
+ bool was_set = false;
+
+ auto sub_isv = input_settings.get_child_optional(settings_key);
+ auto sub_rsv = input_settings.get_child_optional(settings_key);
+
+
+ // User parameter -> first priority
+ if ( cmd_settings.count(vm_key) ) {
+ value = cmd_settings[vm_key].as<std::vector<T>>();
+ was_set = true;
+ }
+
+ // Settings file -> third priority
+ else if ( sub_isv && sub_isv.get().count("el") ) {
+ for ( auto& v : sub_isv.get() ) {
+ if ( v.first == "el" )
+ value.push_back(v.second.get_value<T>());
+ }
+ was_set = true;
+ }
+
+ // RunInfo file -> second priority
+ else if ( sub_rsv && sub_rsv.get().count("el") ) {
+ for ( auto& v : sub_rsv.get() ) {
+ if ( v.first == "el" )
+ value.push_back(v.second.get_value<T>());
+ }
+ was_set = true;
+ }
+
+ // Throw exception if unset
+ if ( required && !was_set )
+ throw po::required_option(vm_key);
+
+ if ( !was_set) {
+ value = default_value;
+ }
+
+ // Otherwise set value
+ auto binded_function = std::bind(function, &globalAlignmentSettings, std::placeholders::_1);
+ binded_function(value);
+ }
+
+public:
+
+ /**
+ * Default constructor.
+ * @argC Number of command line arguments.
+ * @argV List of command line arguments.
+ * @author Tobias Loka
+ */
+ explicit ArgumentParser(int argC, char const ** argV);
+
+ /**
+ * Virtual destructor.
+ * @author Tobias Loka
+ */
+ virtual ~ArgumentParser(){};
+
+};
+
+/**
+ * Class to parse arguments for HiLive build.
+ */
+class BuildIndexArgumentParser : public ArgumentParser {
+
+ uint16_t kmer_weight;
+
+ std::vector<unsigned> gap_positions = {};
+
+ /**
+ * Use the constructor of the inherited ArgumentParser class.
+ */
+ using ArgumentParser::ArgumentParser;
+
+ /**
+ * General options of HiLive build.
+ * @return Option descriptor containing all general options that can be set by the user.
+ * @author Martin Lindner
+ */
+ po::options_description general_options();
+
+ /**
+ * Positional options of HiLive build.
+ * @return Option descriptor containing all positional options that must be set by the user.
+ * @author Martin Lindner
+ */
+ po::options_description positional_options();
+
+ /**
+ * Build options of HiLive build.
+ * @return Option descriptor containing all positional options that must be set by the user.
+ * @author Martin Lindner
+ */
+ po::options_description build_options();
+
+ /**
+ * Set all variables for the positional command arguments.
+ * @param vm The variables map containing the user parameters.
+ * @return true on success, false otherwise
+ * @author Tobias Loka
+ */
+ bool set_positional_variables(po::variables_map vm);
+
+ /**
+ * Set all variables for the build arguments.
+ * @param vm The variables map containing the user parameters.
+ * @return true on success, false otherwise
+ * @author Tobias Loka
+ */
+ bool set_build_variables(po::variables_map vm);
+
+ void report() override;
+
+ void init_help(po::options_description visible_options) override;
+
+public:
+
+ // name of the index file
+ std::string index_name;
+
+ // name of the input fasta file
+ std::string fasta_name;
+
+ // trimming parameter
+ unsigned trim;
+
+ // do_not_convert_spaces_switch
+ bool do_not_convert_spaces;
+
+ // trim_ids switch
+ bool trim_ids;
+
+ int parseCommandLineArguments() override;
+
+};
+
+
+/**
+ * Class to parse arguments for HiLive.
+ */
+class HiLiveArgumentParser : public ArgumentParser {
+
+protected:
+
+ /**
+ * Use the constructor of the inherited ArgumentParser class.
+ */
+ using ArgumentParser::ArgumentParser;
+
+ /**
+ * General options of HiLive.
+ * @return Option descriptor containing all general options that can be set by the user.
+ * @author Martin Lindner
+ */
+ po::options_description general_options();
+
+ /**
+ * Positional options of HiLive.
+ * @return Option descriptor containing all positional options that must be set by the user.
+ * @author Martin Lindner
+ */
+ po::options_description positional_options();
+
+ /**
+ * I/O options of HiLive.
+ * @return Option descriptor containing all I/O options that can be set by the user.
+ * @author Martin Lindner
+ */
+ po::options_description io_options();
+
+ /**
+ * Alignment options of HiLive.
+ * @return Option descriptor containing all alignment options that can be set by the user.
+ * @author Martin Lindner
+ */
+ po::options_description alignment_options();
+
+ /**
+ * Technical options of HiLive.
+ * @return Option descriptor containing all technical options that can be set by the user.
+ * @author Martin Lindner
+ */
+ po::options_description technical_options();
+
+ /**
+ * Check all paths that are relevant for the functionality of HiLive.
+ * @return true if all paths and files are accessible
+ * @author Jakob Schulze
+ */
+ bool checkPaths();
+
+ /**
+ * Parse Lanes, Tiles and read fragments from a RunInfo.xml file.
+ * @param vm The variables map containing the user parameters.
+ * @return true on success, false otherwise
+ * @author Jakob Schulze
+ */
+ bool parseRunInfo(po::variables_map vm);
+
+ virtual void report() override;
+
+ void init_help(po::options_description visible_options) override;
+
+ bool set_options();
+
+ virtual void set_required_parameters() override { required_options = {"BC_DIR", "INDEX", "CYCLES"}; }
+
+public:
+
+ int parseCommandLineArguments() override;
+
+};
+
+/**
+ * Class to parse arguments for HiLive out.
+ */
+
+class HiLiveOutArgumentParser : public HiLiveArgumentParser {
+
+ using HiLiveArgumentParser::HiLiveArgumentParser;
+
+ void init_help(po::options_description visible_options) override;
+
+ void report() override;
+
+ void set_required_parameters() override { required_options = {"settings", "INDEX"}; };
+};
+
+//class HiLiveOutArgumentParser : public ArgumentParser {
+//
+// /**
+// * Use the constructor of the inherited ArgumentParser class.
+// */
+// using ArgumentParser::ArgumentParser;
+//
+// /**
+// * General options of HiLive build.
+// * @return Option descriptor containing all general options that can be set by the user.
+// * @author Martin Lindner
+// */
+// po::options_description general_options();
+//
+// /**
+// * Positional options of HiLive build.
+// * @return Option descriptor containing all positional options that must be set by the user.
+// * @author Martin Lindner
+// */
+// po::options_description positional_options();
+//
+// /**
+// * Build options of HiLive build.
+// * @return Option descriptor containing all positional options that must be set by the user.
+// * @author Martin Lindner
+// */
+// po::options_description output_options();
+//
+// /**
+// * Set all variables for the positional command arguments.
+// * @param vm The variables map containing the user parameters.
+// * @return true on success, false otherwise
+// * @author Tobias Loka
+// */
+// bool set_positional_variables(po::variables_map vm);
+//
+// /**
+// * Set all variables for the build arguments.
+// * @param vm The variables map containing the user parameters.
+// * @return true on success, false otherwise
+// * @author Tobias Loka
+// */
+// bool set_output_variables(po::variables_map vm);
+//
+// void report() override;
+//
+// void init_help(po::options_description visible_options) override;
+//
+// bool set_options() override;
+//
+//public:
+//
+// int parseCommandLineArguments() override;
+//
+//};
diff --git a/lib/config.h.in b/lib/config.h.in
deleted file mode 100644
index 16fee5b..0000000
--- a/lib/config.h.in
+++ /dev/null
@@ -1,7 +0,0 @@
-// the k-mer length
-#define K_HiLive @HiLive_K@
-
-// the HiLive Version Number
-#define HiLive_VERSION_MAJOR @HiLive_VERSION_MAJOR@
-#define HiLive_VERSION_MINOR @HiLive_VERSION_MINOR@
-
diff --git a/lib/definitions.h b/lib/definitions.h
index 5956458..1a72e8f 100644
--- a/lib/definitions.h
+++ b/lib/definitions.h
@@ -4,182 +4,340 @@
#include "headers.h"
-
-// bit representation of A/C/G/T.
+/////////////////////////////////////////////
+////////// Sequences / Nucleotides //////////
+/////////////////////////////////////////////
+
+/**
+ * Two-bit representation of a nucleotide.
+ * @param ch Nucleotide as char
+ * @return 2-bit representation
+ */
#define twobit_repr(ch) ((toupper(ch)) == 'A' ? 0LL : \
(toupper(ch)) == 'C' ? 1LL : \
(toupper(ch)) == 'G' ? 2LL : 3LL)
-// complement bit representation of A/C/G/T.
+/**
+ * Complementary two-bit representation of a nucleotide.
+ * @param ch Nucleotide as char
+ * @return Complementary 2-bit representation
+ */
#define twobit_comp(ch) ((toupper(ch)) == 'A' ? 3LL : \
(toupper(ch)) == 'C' ? 2LL : \
(toupper(ch)) == 'G' ? 1LL : 0LL)
-// bit representation to character
+/**
+ * Nucleotide of a 2-bit representation.
+ * @param n 2-bit representation of a nucleotide.
+ * @return Nucleotide as char
+ */
#define revtwobit_repr(n) ((n) == 0 ? 'A' : \
(n) == 1 ? 'C' : \
(n) == 2 ? 'G' : 'T')
-
-
-
-// Allowed characters in sequence
+/**
+ * Supported nucleotides.
+ */
const std::string seq_chars = "ACGTacgt";
-// largest number we're going to hash into. (8 bytes/64 bits/32 nt)
-// probably 32 bit/16 nt are enough here
-typedef uint64_t HashIntoType;
-// construct a mask to truncate a binary representation of a k-mer to length K
-const HashIntoType MASK = HashIntoType(pow(4,K_HiLive))-1;
+////////////////////////////////////////
+////////// Genome Identifiers //////////
+////////////////////////////////////////
-// identifiers for genome sequences
+/**
+ * Type for the identifier of genomes (gid).
+ */
typedef uint32_t GenomeIdType;
+
+/**
+ * Constant variable to tag a k-mer as "trimmed".
+ */
const GenomeIdType TRIMMED = std::numeric_limits<GenomeIdType>::max();
-// list of genome identifiers
+/**
+ * A list of Genome Ids
+ */
typedef std::vector<GenomeIdType> GenomeIdListType;
-// list of strings
-typedef std::vector<std::string> StringListType;
-// position in a genome
+//////////////////////////////////////
+////////// Genome Positions //////////
+//////////////////////////////////////
+
+/**
+ * Type for positions in a genome.
+ */
typedef int32_t PositionType;
-// pair of genome ID and position
+/**
+ * A pair of genome ID and position.
+ */
struct GenomePosType {
+
GenomeIdType gid;
PositionType pos;
GenomePosType()=default;
GenomePosType(GenomeIdType g, PositionType p): gid(g), pos(p) {};
};
-//typedef std::tuple<GenomeIdType,PositionType> GenomePosType;
-// size of a pair of genome ID and position
+/**
+ * Size of a pair of genome ID and position (in bytes)
+ */
const uint64_t GenomePos_size = sizeof(GenomeIdType) + sizeof(PositionType);
-// vector of ID:position pairs
+/**
+ * A vector of GenomePosTypes.
+ */
typedef std::vector<GenomePosType> GenomePosListType;
-// iterator on GenomePosList
+/**
+ * Iterator on GenomePosListType.
+ */
typedef GenomePosListType::iterator GenomePosListIt;
-// the k-mer index array
-const HashIntoType n_kmer = pow(4,K_HiLive);
-typedef std::array<GenomePosListType,n_kmer> KmerIndexType;
-
-// small counters
-typedef uint16_t CountType;
-
-// difference between k-mer position in the read and matching position in the reference
-typedef int16_t DiffType;
-
-// define a mismatch as max(DiffType)
-const DiffType NO_MATCH = std::numeric_limits<DiffType>::max();
-
-// define a trimmed position as max(DiffType)-1
-const DiffType TRIMMED_MATCH = std::numeric_limits<DiffType>::max()-1;
-
-// one element in Cigar vector containing match/mismatch information about consecutive k-mers
-struct CigarElement {
- CountType length;
- DiffType offset;
- CigarElement (CountType l, DiffType o): length(l), offset(o) {};
- CigarElement (): length(0), offset(NO_MATCH) {};
-};
-
-// CigarVector containing CIGAR string like information about the alignments
-typedef std::list<CigarElement> CigarVector;
-
+/////////////////////////////////
+////////// K-mer index //////////
+/////////////////////////////////
+/**
+ * Type to hash k-mers into.
+ * This type also limits the k-mer weight (currently to 32).
+ */
+typedef uint64_t HashIntoType;
-// all user parameters are stored in the alignment settings
-struct AlignmentSettings {
- // HARD CODED: kmer gap structure (this is not used anywhere)
- //std::string kmer_structure = "11111110111110111";
- std::string kmer_structure = "111111111111111";
+/**
+ * K-mer index type.
+ */
+typedef std::vector<GenomePosListType> KmerIndexType;
- // HARD CODED: kmer gap positions (one-based)
- //std::vector<unsigned> kmer_gaps = {8, 14};
- std::vector<unsigned> kmer_gaps;
+/**
+ * A lightweight type for storing the index.
+ */
+typedef std::vector<char*> KixRunDB;
- // HARD CODED: kmer span (kmer weight is K_HiLive)
- //unsigned kmer_span = K_HiLive+2;
- unsigned kmer_span = K_HiLive;
- // PARAMETER: Base Call quality cutoff, treat BC with quality < bc_cutoff as miscall
- CountType min_qual;
+////////////////////////////////////////
+////////// Integer data types //////////
+////////////////////////////////////////
- // PARAMETER: max. insert/deletion size
- DiffType window;
+/**
+ * Type for small counters.
+ */
+typedef uint16_t CountType;
- // PARAMETER: minimum number of errors allowed in alignment
- CountType min_errors;
+/**
+* Difference between k-mer position in the read and matching position in the reference.
+*/
+typedef int16_t DiffType;
- // SWITCH: discard One-hit-wonders
- bool discard_ohw;
- // PARAMETER: first cycle to discard one-hit-wonders
- CountType start_ohw;
+////////////////////////////////////////
+////////// Offset definitions //////////
+////////////////////////////////////////
- // SWITCH: Best-Hit-Mode
- bool any_best_hit_mode;
+/**
+ * Define a mismatch as maximum value of DiffType.
+ */
+const DiffType NO_MATCH = std::numeric_limits<DiffType>::max();
- // SWITCH: Best-Hit-Mode
- bool all_best_hit_mode;
+/**
+ * Define a trimmed match maximum value of DiffType -1.
+ */
+const DiffType TRIMMED_MATCH = std::numeric_limits<DiffType>::max()-1;
- // SWITCH: Best-N-Mode
- bool all_best_n_scores_mode;
- // PARAMETER: Best-N-Mode::N
- CountType best_n;
+////////////////////////////////////
+////////// CIGAR elements //////////
+////////////////////////////////////
- // PARAMETER: temporary directory for the streamed alignment
- std::string temp_dir;
+/**
+ * One (internal) CIGAR element.
+ */
+struct CigarElement {
- // SWITCH: write sam/bam output or not
- bool write_bam=false;
+ /** Length of the region. */
+ CountType length;
- // SWITCH: Keep the old alignment files of previous cycles
- bool keep_aln_files;
+ /** Offset of the region to the original start pos (created by InDels). */
+ DiffType offset;
- // PARAMETER: Memory block size for the input and output buffer in the streamed alignment
- uint64_t block_size;
+ CigarElement (CountType l, DiffType o): length(l), offset(o) {};
+ CigarElement (): length(0), offset(NO_MATCH) {};
+};
- // PARAMETER: Compression format for alignment files
- uint8_t compression_format;
+/**
+ * Vector of CIGAR elements, representing the alignment information for one seed.
+ */
+typedef std::list<CigarElement> CigarVector;
- // PARAMETER: list of lanes to process
- std::vector<uint16_t> lanes;
-
- // PARAMETER: list of tiles to process
- std::vector<uint16_t> tiles;
- // PARAMETER: root directory of hilive run
- std::string root;
+///////////////////////////////////////
+////////// Sequence Elements //////////
+///////////////////////////////////////
+
+/**
+ * Information about the sequences.
+ * One element can be a read or a barcode.
+ * @author Tobias Loka
+ */
+struct SequenceElement {
+
+ /** The id of the read. Equals the position in the argument list and in the AlignmentSettings::seqs vector (0-based). */
+ CountType id;
+
+ /** The mate number. 0 for barcodes, increasing for sequence reads in the given order (1-based). */
+ CountType mate;
+
+ /** The length of the respective read. */
+ CountType length;
+
+ /**
+ * Constructor of a SequenceElement NULL object.
+ * @author Tobias Loka
+ */
+ SequenceElement () : id(0), mate(0), length(0) {};
+
+ /**
+ * Constructor of a valid SequenceElement object.
+ * @param id The id of the read.
+ * @param m The mate number of the read (0 for barcodes, incrementing for sequence reads)
+ * @param l The length of the read
+ * @author Tobias Loka
+ */
+ SequenceElement (CountType id, CountType m, CountType l): id(id), mate(m), length(l) {};
+
+ /**
+ * Check whether the SequenceElement object is a barcode or not.
+ * @return true, if SequenceElement is a barcode. False if not.
+ * @author Tobias Loka
+ */
+ bool isBarcode() { return (mate==0);}
+};
- // PARAMETER: path to the index file
- std::string index_fname;
+/**
+ * Check if two Sequence elements are equal.
+ */
+inline bool operator==(const SequenceElement l, const SequenceElement r) {return (l.length==r.length) && (l.mate==r.mate) && (l.id==r.id);}
+
+/**
+ * Checks if two sequence elements are not equal.
+ */
+inline bool operator!=(const SequenceElement l, const SequenceElement r) {return !(l==r);}
+
+/**
+ * An undefined sequence element (NULL element).
+ */
+const SequenceElement NULLSEQ = SequenceElement();
+
+
+////////////////////////////////////////////
+////////// Unmodifiable variables //////////
+////////////////////////////////////////////
+
+/**
+ * Exception specialization for Unmodifiable data types.
+ * @author Tobias Loka
+ */
+class unmodifiable_error : public std::logic_error
+{
+public:
+ using std::logic_error::logic_error;
+};
- // PARAMETER: read length of all reads (including barcodes)
- CountType rlen;
+/**
+ * Template to define data types that can only be set once.
+ * @type T Data type of the unmodifiable object.
+ * @author Tobias Loka
+ */
+template <typename T>
+class Unmodifiable {
+
+private:
+
+ /** The unmodifiable object. */
+ T unmodifiable_object;
+
+ /** Flag to check if the object was already set once. */
+ bool setFlag = false;
+
+public:
+
+ /** Constructor without setting the object (to only declare the object).*/
+ Unmodifiable(){ }
+
+ /** Constructor with setting the object (to init the object).*/
+ Unmodifiable(T object) {
+ unmodifiable_object = object;
+ }
+
+ /** Automatic cast to of the unmodifiable to the object type. */
+ operator T() { return unmodifiable_object; }
+
+ /**
+ * Set the unmodifiable object (will only work once!).
+ * @param object The object to be copied to this unmodifiable data type.
+ * @return true if setting was successful
+ * @author Tobias Loka
+ */
+ void set(T object) {
+ if ( isSet() ) {
+ throw unmodifiable_error("Tried to modify unmodifiable object");
+ }
+
+ unmodifiable_object = object;
+ setFlag = true;
+ }
+
+ /**
+ * Check if the object was already set.
+ * @return true if the object was already set.
+ * @author Tobias Loka
+ */
+ bool isSet() {
+ return setFlag;
+ }
+
+ /**
+ * Return a copy of the unmodifiable object.
+ * @param allow_unset if false, an exception is thrown when the object was not set before. Should only be true for
+ * objects that require access to certain properties before their initialization (e.g. to check a container's size
+ * without knowing if the container was already set).
+ * @return (copy/value of) the unmodifiable object
+ * @author Tobias Loka
+ */
+ T get(bool allow_unset = false ) {
+ if ( ! isSet() && ! allow_unset) {
+ throw unmodifiable_error("Tried to access uninitialized object");
+ }
+
+ return unmodifiable_object;
+ }
- // PARAMETER: length of the sequence of all reads (excluding barcodes)
- CountType seqlen;
+};
- // PARAMETER: vector containing all barcodes of the reads which should be outputted
- std::vector<std::string> barcodeVector;
+///////////////////////////////////////
+////////// Other definitions //////////
+///////////////////////////////////////
- // PARAMETER: directory in which to create the output directory structure
- std::string out_dir;
+/**
+ * A list of strings
+ */
+typedef std::vector<std::string> StringListType;
- // PARAMETER: number of threads to use
- CountType num_threads;
+/**
+ * The different alignment modes.
+ * @author Tobias Loka
+ */
+enum AlignmentMode:char {
+ ALL='A',
+ ALLBEST='H',
+ ANYBEST='B',
+ BESTN='N',
+ UNKNOWN='U'
};
-
-
#endif /* DEFINITIONS_H */
diff --git a/lib/global_variables.h b/lib/global_variables.h
new file mode 100644
index 0000000..28ac3cc
--- /dev/null
+++ b/lib/global_variables.h
@@ -0,0 +1,7 @@
+#ifndef GLOBAL_VARIABLES_H
+#define GLOBAL_VARIABLES_H
+
+class AlignmentSettings;
+
+extern AlignmentSettings globalAlignmentSettings;
+#endif
diff --git a/lib/headers.h b/lib/headers.h
index ca54a5c..6f05f62 100644
--- a/lib/headers.h
+++ b/lib/headers.h
@@ -1,4 +1,3 @@
-#include "config.h"
#include <iostream>
#include <sstream>
#include <fstream>
@@ -10,6 +9,7 @@
#include <string>
#include <cstring>
#include <math.h>
+#include <numeric>
#include <assert.h>
#include <tuple>
#include <algorithm>
@@ -24,3 +24,10 @@
#include <stdexcept>
#include <utility>
#include <boost/filesystem.hpp>
+#include <boost/property_tree/xml_parser.hpp>
+#include <boost/property_tree/ptree.hpp>
+#include <boost/optional/optional.hpp>
+#include <boost/foreach.hpp>
+#include <seqan/basic.h>
+#include <seqan/sequence.h>
+#include <seqan/bam_io.h>
diff --git a/lib/kindex.cpp b/lib/kindex.cpp
index 77d90b7..167b449 100644
--- a/lib/kindex.cpp
+++ b/lib/kindex.cpp
@@ -1,12 +1,17 @@
#include "kindex.h"
-int KixBuild::add_fasta(const std::string &fname, AlignmentSettings & settings, bool convert_spaces, bool trim_ids) {
+KixBuild::KixBuild() { // default constructor
+ this->db.resize(pow(4,globalAlignmentSettings.get_kmer_weight()));
+}
+
+
+int KixBuild::add_fasta(const std::string &fname, bool convert_spaces, bool trim_ids) {
GenomeIdListType added_ids;
- return add_fasta(fname, added_ids, settings, convert_spaces, trim_ids);
+ return add_fasta(fname, added_ids, convert_spaces, trim_ids);
}
-int KixBuild::add_fasta(const std::string &fname, GenomeIdListType &ids, AlignmentSettings & settings, bool convert_spaces, bool trim_ids) {
+int KixBuild::add_fasta(const std::string &fname, GenomeIdListType &ids, bool convert_spaces, bool trim_ids) {
std::ios::sync_with_stdio(false);
std::ifstream::sync_with_stdio(false);
@@ -43,19 +48,19 @@ int KixBuild::add_fasta(const std::string &fname, GenomeIdListType &ids, Alignme
added_ids.push_back(seq_id);
seq_names.push_back(seq_name);
- seq_lengths.push_back(0); // gets later corrected to settings.kmer_span - 1
+ seq_lengths.push_back(0); // gets later corrected to globalAlignmentSettings.get_kmer_span() - 1
assert(seq_names.size() == num_seq);
startNewSequence = true;
}
else { // sequence line
if (startNewSequence) {
- if (line.length() < settings.kmer_span)
+ if (line.length() < globalAlignmentSettings.get_kmer_span())
continue; // ignore sequences shorter than K
- start_sequence(line, tailingKmer, sequencePosition, settings);
+ start_sequence(line, tailingKmer, sequencePosition);
startNewSequence = false;
}
else
- continue_sequence(line, tailingKmer, sequencePosition, settings);
+ continue_sequence(line, tailingKmer, sequencePosition);
}
}
infile.close();
@@ -68,9 +73,9 @@ int KixBuild::add_fasta(const std::string &fname, GenomeIdListType &ids, Alignme
A new ID is created for this sequence.
Return: sequence ID
*/
-GenomeIdType KixBuild::start_sequence(const std::string &s, std::string& tailingKmer, PositionType& sequencePosition, AlignmentSettings & settings) {
+GenomeIdType KixBuild::start_sequence(const std::string &s, std::string& tailingKmer, PositionType& sequencePosition) {
assert(seq_names.size() == num_seq);
- assert(s.length() >= settings.kmer_span);
+ assert(s.length() >= globalAlignmentSettings.get_kmer_span());
// add sequence kmers to index
sequencePosition = 0; // use 1-based positions (to allow for negative positions)
@@ -78,23 +83,23 @@ GenomeIdType KixBuild::start_sequence(const std::string &s, std::string& tailing
std::string::const_iterator last_invalid;
HashIntoType fw; // forward k-mer
- seq_lengths.back()=settings.kmer_span-1;
- for (; it_s < s.end()-settings.kmer_span+1; ++it_s) {
+ seq_lengths.back()=globalAlignmentSettings.get_kmer_span()-1;
+ for (; it_s < s.end()-globalAlignmentSettings.get_kmer_span()+1; ++it_s) {
++sequencePosition; // use 1-based positions (to allow for negative positions)
seq_lengths.back()+=1;
- last_invalid = hash_fw(it_s, s.end(), fw, settings);
+ last_invalid = hash_fw(it_s, s.end(), fw);
// add k-mer to database
if (last_invalid < it_s)
add_kmer(fw,num_seq-1,sequencePosition);
else {
- unsigned jumplength = std::min(last_invalid - it_s, s.end() - settings.kmer_span - it_s);
+ unsigned jumplength = std::min(last_invalid - it_s, s.end() - globalAlignmentSettings.get_kmer_span() - it_s);
sequencePosition += jumplength;
seq_lengths.back() += jumplength;
it_s = last_invalid;
}
}
- tailingKmer = s.substr(s.length()-settings.kmer_span);
+ tailingKmer = s.substr(s.length()-globalAlignmentSettings.get_kmer_span());
return num_seq;
}
@@ -102,7 +107,7 @@ GenomeIdType KixBuild::start_sequence(const std::string &s, std::string& tailing
/* Continue adding all k-mers in a sequence string s to the database.
Return: sequence ID
*/
-GenomeIdType KixBuild::continue_sequence(const std::string &s, std::string& tailingKmer, PositionType& sequencePosition, AlignmentSettings & settings) {
+GenomeIdType KixBuild::continue_sequence(const std::string &s, std::string& tailingKmer, PositionType& sequencePosition) {
assert(seq_names.size() == num_seq);
std::string concatString = tailingKmer + s;
@@ -111,22 +116,22 @@ GenomeIdType KixBuild::continue_sequence(const std::string &s, std::string& tail
std::string::const_iterator last_invalid;
HashIntoType fw; // forward k-mer
- for (; it_s < concatString.end()-settings.kmer_span+1; ++it_s) {
+ for (; it_s < concatString.end()-globalAlignmentSettings.get_kmer_span()+1; ++it_s) {
++sequencePosition; // use 1-based positions (to allow for negative positions)
seq_lengths.back()+=1;
- last_invalid = hash_fw(it_s, concatString.end(), fw, settings);
+ last_invalid = hash_fw(it_s, concatString.end(), fw);
// add k-mer to database
if (last_invalid < it_s)
add_kmer(fw,num_seq-1,sequencePosition);
else {
- unsigned jumplength = std::min(last_invalid - it_s, concatString.end()- settings.kmer_span - it_s);
+ unsigned jumplength = std::min(last_invalid - it_s, concatString.end()- globalAlignmentSettings.get_kmer_span() - it_s);
sequencePosition += jumplength;
seq_lengths.back() += jumplength;
it_s = last_invalid;
}
}
- tailingKmer = concatString.substr(concatString.length()-settings.kmer_span);
+ tailingKmer = concatString.substr(concatString.length()-globalAlignmentSettings.get_kmer_span());
return num_seq;
}
@@ -161,6 +166,7 @@ uint64_t KixBuild::trim(uint64_t max_count) {
std::vector<char> KixBuild::serialize() {
+
// first of all, sort the database entries by position
for (auto it = db.begin(); it != db.end(); ++it)
std::sort(it->begin(), it->end(), gp_compare);
@@ -168,9 +174,15 @@ std::vector<char> KixBuild::serialize() {
// calculate total size
unsigned long int total_size = 0;
- // K_HiLive itself
+ // The k-mer weight itself
+ total_size += 1;
+
+ // The number of gaps
total_size += 1;
+ // The gaps themselves
+ total_size += globalAlignmentSettings.get_kmer_gaps().size();
+
// total number of sequences in database
total_size += sizeof(GenomeIdType);
@@ -200,10 +212,22 @@ std::vector<char> KixBuild::serialize() {
char* d = data.data();
// write K
- uint8_t kk = K_HiLive;
+ uint8_t kk = globalAlignmentSettings.get_kmer_weight();
memcpy(d,&kk,1);
d++;
+ // number of gaps
+ std::vector<unsigned> kmer_gaps = globalAlignmentSettings.get_kmer_gaps();
+ uint8_t gap_num = kmer_gaps.size();
+ memcpy(d,&gap_num,1);
+ d++;
+
+ // The gaps themselves
+ for ( uint8_t gap : kmer_gaps) {
+ memcpy(d,&gap,1);
+ d++;
+ }
+
// total number of sequences in database
memcpy(d,&num_seq,sizeof(GenomeIdType));
d += sizeof(GenomeIdType);
@@ -276,121 +300,34 @@ uint64_t KixBuild::serialize_file(std::string f) {
}
-uint64_t KixBuild::deserialize(char* d) {
+uint64_t KixRun::deserialize(char* d) {
+
// the total number of bytes read
uint64_t bytes = 0;
- // read and check K
+ // read k-mer weight
uint8_t kk;
memcpy(&kk,d+bytes,1);
bytes++;
- assert(K_HiLive == kk);
-
- // read total number of sequences in database
- memcpy(&num_seq,d+bytes,sizeof(GenomeIdType));
- bytes += sizeof(GenomeIdType);
-
-
- // sequence names
- seq_names.clear();
- seq_lengths.clear();
- for (uint32_t i = 0; i < num_seq; i++) {
- uint16_t nm_length;
- memcpy(&nm_length,d+bytes,sizeof(uint16_t));
- bytes += sizeof(uint16_t);
-
- char * tmp = new char[nm_length+1];
- memcpy(tmp,d+bytes,nm_length);
- tmp[nm_length] = 0; // make the string null-terminated
- seq_names.push_back(tmp);
- delete tmp;
- bytes += nm_length;
- }
- for (uint32_t i = 0; i < num_seq; i++) {
- uint64_t seq_len;
- memcpy(&seq_len,d+bytes,sizeof(uint64_t));
- bytes += sizeof(uint64_t);
- seq_lengths.push_back(seq_len);
- }
+ this->kmer_weight = kk;
-
- // database entries
- for (auto it = db.begin(); it != db.end(); ++it) {
- // number of positions
- uint32_t num_positions;
- memcpy(&num_positions,d+bytes,sizeof(uint32_t));
- bytes += sizeof(uint32_t);
- (*it).clear();
- (*it).reserve(num_positions);
-
- // genome ID and position
- for( uint32_t i = 0; i < num_positions; i++) {
- GenomeIdType gid;
- PositionType pos;
- GenomePosType gp;
-
- memcpy(&gid,d+bytes,sizeof(GenomeIdType));
- bytes += sizeof(GenomeIdType);
-
- memcpy(&pos,d+bytes,sizeof(PositionType));
- bytes += sizeof(PositionType);
-
- gp.gid = gid;
- gp.pos = pos;
-
- (*it).push_back(gp);
- }
- }
-
- return bytes;
-}
-
-
-uint64_t KixBuild::deserialize_file(std::string f) {
- std::string fname = f;
-
- // obtain file size
- unsigned long int size = get_filesize(fname);
-
- // open binary file
- FILE* ifile;
- ifile = fopen(fname.c_str(), "rb");
-
- if (!ifile) {
- std::cerr << "Error reading from file " << fname << ": Could not open file." << std::endl;
- return 0;
- }
-
- // allocate memory
- std::vector<char> sdata (size);
-
- // read all data
- unsigned long int read = fread(sdata.data(), 1, size, ifile);
-
- // close file
- fclose (ifile);
+ // read number of gaps in k-mer pattern
+ uint8_t gap_num;
+ memcpy(&gap_num,d+bytes,1);
+ bytes++;
- if (read != size){
- std::cerr << "Error reading from file " << fname << ": File size: " << size << " bytes. Read: " << read << " bytes." << std::endl;
- return 0;
+ // read k-mer pattern
+ for ( uint8_t i = 0; i < gap_num; i++ ) {
+ uint8_t gap;
+ memcpy(&gap, d+bytes, 1);
+ kmer_gaps.push_back(gap);
+ bytes++;
}
- // deserialize data
- deserialize(sdata.data());
- return read;
-}
-
-
-uint64_t KixRun::deserialize(char* d) {
- // the total number of bytes read
- uint64_t bytes = 0;
-
- // read and check K
- uint8_t kk;
- memcpy(&kk,d+bytes,1);
- bytes++;
- assert(kk == K_HiLive);
+// globalAlignmentSettings.set_kmer_weight(this->kmer_weight);
+ store_kmer();
+ this->db.resize(pow(4,globalAlignmentSettings.get_kmer_weight()));
// read total number of sequences in database
memcpy(&num_seq,d+bytes,sizeof(GenomeIdType));
@@ -450,14 +387,22 @@ uint64_t KixRun::deserialize_file(std::string f) {
}
+// return k-mer weight read from index
+uint8_t KixRun::get_kmer_weight() {
+ return(this->kmer_weight);
+}
+
/* Retrieve all occurrences (fwd & rc) of kmer in the reference from the index */
-GenomePosListType KixRun::retrieve_positions(std::string kmerSpan, AlignmentSettings & settings) {
+GenomePosListType KixRun::retrieve_positions(std::string kmerSpan) {
// get the reverse complement of the kmer
HashIntoType fwHashValue;
HashIntoType rcHashValue;
- hash(kmerSpan.c_str(), fwHashValue, rcHashValue, settings);
+ hash(kmerSpan.c_str(), fwHashValue, rcHashValue);
+// std::cout << "FW: " << fwHashValue << std::endl;
+// std::cout << "REV: " << rcHashValue << std::endl;
+
// obtain the list of positions for each k-mer
char* fwd_begin = db[fwHashValue];
uint32_t fwd_len;
@@ -499,3 +444,67 @@ GenomePosListType KixRun::retrieve_positions(std::string kmerSpan, AlignmentSett
return pos;
}
+
+/**
+ * get kix-header information only to not load the complete index
+ * structure:
+ * 8 bit (k-mer size)
+ * 32 bit (# sequences)
+ * for each sequence: 16 bit (sequence name length) + sequence name
+ * for each sequence: 32 bit (length of sequence)
+ */
+uint64_t KixRun::get_header_information(std::string f) {
+ std::string fname = f;
+
+ // open binary file
+ FILE* fi;
+ fi = fopen(fname.c_str(), "rb");
+
+ if (!fi) {
+ std::cerr << "Error reading binary file " << fname << ": Could not open file." << std::endl;
+ return 0;
+ }
+
+ uint16_t seq_name_len;
+ uint64_t seq_amount, seq_len;
+ int c;
+
+ // sequence names & lengths
+ seq_names.clear();
+ seq_lengths.clear();
+ fread(&kmer_weight,1,1,fi);
+
+
+ // Ignore k-mer gaps
+ uint8_t gap_num;
+ fread(&gap_num,1,sizeof(uint8_t),fi);
+
+ // read k-mer pattern
+ for ( uint8_t i = 0; i < gap_num; i++ ) {
+ uint8_t gap;
+ fread(&gap, 1, sizeof(uint8_t), fi);
+ }
+
+ fread(&seq_amount,4,1,fi);
+
+ //get the names for every sequence
+ for( unsigned a = 0; a < seq_amount; a = a + 1 ){
+ fread(&seq_name_len,2,1,fi);
+ std::string seq_name;
+ for( unsigned b = 0; b < seq_name_len; b = b + 1 ){
+ c = fgetc(fi);
+ seq_name.push_back(c);
+ }
+ seq_names.push_back(seq_name);
+ }
+
+ for( unsigned a = 0; a < seq_amount; a = a + 1 ) {
+ fread(&seq_len,4,1,fi);
+ }
+ seq_lengths.push_back(seq_len);
+
+
+ fclose(fi);
+ return seq_amount;
+}
+
diff --git a/lib/kindex.h b/lib/kindex.h
index 5b238fa..be76bc6 100644
--- a/lib/kindex.h
+++ b/lib/kindex.h
@@ -3,6 +3,7 @@
#include "headers.h"
#include "definitions.h"
+#include "global_variables.h"
#include "tools.h"
@@ -18,13 +19,16 @@ class KixBuild {
public:
+ // constructor resizing db (see below) to match the number of possible k-mers
+ KixBuild();
+
// add k-mers of all sequences in FASTA file
- int add_fasta(const std::string &fname, GenomeIdListType &ids, AlignmentSettings & settings, bool convert_spaces, bool trim_ids);
- int add_fasta(const std::string &fname, AlignmentSettings & settings, bool convert_spaces, bool trim_ids);
+ int add_fasta(const std::string &fname, GenomeIdListType &ids, bool convert_spaces, bool trim_ids);
+ int add_fasta(const std::string &fname, bool convert_spaces, bool trim_ids);
// add all k-mers in a string sequence to the database
- GenomeIdType start_sequence(const std::string &s, std::string& tailingKmer, PositionType& sequencePosition, AlignmentSettings & settings);
- GenomeIdType continue_sequence(const std::string &s, std::string& tailingKmer, PositionType& sequencePosition, AlignmentSettings & settings);
+ GenomeIdType start_sequence(const std::string &s, std::string& tailingKmer, PositionType& sequencePosition);
+ GenomeIdType continue_sequence(const std::string &s, std::string& tailingKmer, PositionType& sequencePosition);
// trim the database: remove kmers with more than max_count occurrences
uint64_t trim(uint64_t max_count);
@@ -35,14 +39,8 @@ class KixBuild {
// serialize and store the KixBuild to a file
uint64_t serialize_file(std::string f);
- // deserialize KixBuild
- uint64_t deserialize(char* d);
-
- // load and deserialize KixBuild from file
- uint64_t deserialize_file(std::string f);
-
- GenomeIdType num_seq; // total number of sequences in the database
+ GenomeIdType num_seq=0; // total number of sequences in the database
KmerIndexType db; // the database structure itself
StringListType seq_names; // names of the sequences in the database
std::vector<uint32_t> seq_lengths; // lengths of the sequences in the database
@@ -56,22 +54,30 @@ class KixBuild {
//-------------------------------------------------------------------//
-typedef std::array<char*,n_kmer> KixRunDB;
-
class KixRun {
+ private:
+ uint8_t kmer_weight; // k-mer weight read from file
+ std::vector<unsigned> kmer_gaps;
public:
// pointer to the matching positions for a k-mer
char* kmer(HashIntoType kmer);
// retrieve all fwd and rc occurrences of kmer in the index
- GenomePosListType retrieve_positions(std::string kmerSpan, AlignmentSettings & settings);
+ GenomePosListType retrieve_positions(std::string kmerSpan);
- // deserialize Kix
+ // deserialize Kix, also sets kmer_weight and globalAlignmentSettings.kmer_weight
uint64_t deserialize(char* d);
// load and deserialize Kix from file
uint64_t deserialize_file(std::string f);
+ // return k-mer weight of the k-mers in the index
+ uint8_t get_kmer_weight();
+
+ uint64_t get_header_information(std::string f);
+
+ void store_kmer() { globalAlignmentSettings.set_kmer(kmer_weight, kmer_gaps); };
+
// Database content
GenomeIdType num_seq; // total number of sequences in the database
StringListType seq_names; // names of the sequences in the database
diff --git a/lib/parallel.cpp b/lib/parallel.cpp
index 6d54171..1897195 100644
--- a/lib/parallel.cpp
+++ b/lib/parallel.cpp
@@ -1,14 +1,12 @@
#include "parallel.h"
-
-
std::ostream& operator<<(std::ostream& os, const Task& t)
{
- os << "Lane " << t.lane << " Tile " << t.tile << " Cycle " << t.cycle;
+ std::string mate = t.seqEl.mate == 0 ? "b" : std::to_string(t.seqEl.mate);
+ os << "Lane " << t.lane << " Tile " << t.tile << " Cycle " << mate << "." << t.cycle;
return os;
}
-
// Add element to the task list
void TaskQueue::push(Task t) {
std::lock_guard<std::mutex> lk(m);
@@ -34,8 +32,6 @@ uint64_t TaskQueue::size() {
return tasks.size();
}
-
-
// create a vector with all lane numbers
std::vector<uint16_t> all_lanes() {
std::vector<uint16_t> ln;
@@ -49,7 +45,6 @@ std::vector<uint16_t> one_lane(uint16_t l) {
return std::vector<uint16_t> (1,l);
}
-
// create a vector with all tile numbers
std::vector<uint16_t> all_tiles() {
std::vector<uint16_t> tl;
@@ -69,33 +64,31 @@ std::vector<uint16_t> one_tile(uint16_t t) {
return std::vector<uint16_t> (1,t);
}
-
-// initialize agenda with root directory and read length only (all lanes, all tiles)
-Agenda::Agenda (std::string rt, uint16_t rl) {
+// initialize agenda with read length only (all lanes, all tiles)
+Agenda::Agenda (uint16_t rl) {
// add lanes 1-8 to the list
std::vector<uint16_t> ln = all_lanes();
// call the tiles constructor
- Agenda(rt, rl, ln);
+ Agenda(rl, ln);
}
-// initialize agenda with root directory, read length, and lanes (all tiles)
-Agenda::Agenda (std::string rt, uint16_t rl, std::vector<uint16_t> ln) {
+// initialize agenda with read length and lanes (all tiles)
+Agenda::Agenda (uint16_t rl, std::vector<uint16_t> ln) {
// add all tiles to the list
std::vector<uint16_t> tl = all_tiles();
// call the full constructor
- Agenda (rt, rl, ln, tl);
+ Agenda (rl, ln, tl);
}
-// initialize agenda with root directory, read length, lanes, and tiles
-Agenda::Agenda (std::string rt, uint16_t rl, std::vector<uint16_t> ln, std::vector<uint16_t> tl) {
+// initialize agenda with read length, lanes, and tiles
+Agenda::Agenda (uint16_t rl, std::vector<uint16_t> ln, std::vector<uint16_t> tl) {
- root = rt;
rlen = rl;
lanes = ln;
tiles = tl;
@@ -116,63 +109,62 @@ Agenda::Agenda (std::string rt, uint16_t rl, std::vector<uint16_t> ln, std::vect
// check for BCL files and update item status
void Agenda::update_status () {
- // iterate over lanes
- for (uint16_t ln_id = 0; ln_id < items.size(); ++ln_id) {
-
- // iterate over all tiles
- for (uint16_t tl_id = 0; tl_id < items[ln_id].size(); ++tl_id) {
-
- // get the first cycle that is not in the FINISHED status
- uint16_t first_unfinished = 0;
- while ( (first_unfinished < items[ln_id][tl_id].size()) && (items[ln_id][tl_id][first_unfinished] == FINISHED)) {
- first_unfinished++;
- }
-
- // if there is one, check if there is a BCL file available
- if ((first_unfinished != items[ln_id][tl_id].size()) && (items[ln_id][tl_id][first_unfinished] == WAITING)) {
- std::string this_fname = bcl_name(root, lanes[ln_id], tiles[tl_id], first_unfinished+1);
- // only change the status if the file exists
- if ( file_exists(this_fname) ) {
- // TODO: probably find a way to check if the machine currently writes to that file
- items[ln_id][tl_id][first_unfinished] = BCL_AVAILABLE;
+ // iterate over lanes
+ for (uint16_t ln_id = 0; ln_id < items.size(); ++ln_id) {
+
+ // iterate over all tiles
+ for (uint16_t tl_id = 0; tl_id < items[ln_id].size(); ++tl_id) {
+
+ // get the first cycle that is not in the FINISHED status
+ uint16_t first_unfinished = 0;
+ while ( (first_unfinished < items[ln_id][tl_id].size()) && (items[ln_id][tl_id][first_unfinished] == FINISHED)) {
+ first_unfinished++;
+ }
+
+ // if there is one, check if there is a BCL file available
+ if ((first_unfinished != items[ln_id][tl_id].size()) && (items[ln_id][tl_id][first_unfinished] == WAITING)) {
+ std::string this_fname = bcl_name(lanes[ln_id], tiles[tl_id], first_unfinished+1);
+ // only change the status if the file exists
+ if ( file_exists(this_fname) ) {
+ // TODO: probably find a way to check if the machine currently writes to that file
+ items[ln_id][tl_id][first_unfinished] = BCL_AVAILABLE;
+ }
+ }
+ }
}
- }
-
- }
-
- }
-
}
-
// generate a new task from the agenda
Task Agenda::get_task(){
- // iterate over lanes
- for (uint16_t ln_id = 0; ln_id < items.size(); ++ln_id) {
-
- // iterate over all tiles
- for (uint16_t tl_id = 0; tl_id < items[ln_id].size(); ++tl_id) {
-
- // check if there is a cycle with an unprocessed BCL file
- uint16_t unprocessed = 0;
- while ( (unprocessed < items[ln_id][tl_id].size()) && (items[ln_id][tl_id][unprocessed] != BCL_AVAILABLE)) {
- unprocessed++;
- }
-
- // generate a new task if there is an unprocessed BCL file
- if ( unprocessed != items[ln_id][tl_id].size() ) {
- Task t (lanes[ln_id], tiles[tl_id], unprocessed+1, rlen, root);
- return t;
- }
-
- }
-
- }
- // return indicator that no new task could be created
- return NO_TASK;
+ // iterate over lanes
+ for (uint16_t ln_id = 0; ln_id < items.size(); ++ln_id) {
+
+ // iterate over all tiles
+ for (uint16_t tl_id = 0; tl_id < items[ln_id].size(); ++tl_id) {
+
+ // check if there is a cycle with an unprocessed BCL file
+ uint16_t unprocessed = 0;
+ while ( (unprocessed < items[ln_id][tl_id].size()) && (items[ln_id][tl_id][unprocessed] != BCL_AVAILABLE)) {
+ unprocessed++;
+ }
+
+ // generate a new task if there is an unprocessed BCL file
+ if ( unprocessed != items[ln_id][tl_id].size() ) {
+ uint16_t cycle = unprocessed + 1;
+ uint16_t read_no = 0;
+ while ( cycle > globalAlignmentSettings.getSeqById(read_no).length) {
+ cycle -= globalAlignmentSettings.getSeqById(read_no).length;
+ read_no += 1;
+ }
+ Task t (lanes[ln_id], tiles[tl_id], globalAlignmentSettings.getSeqById(read_no), cycle);
+ return t;
+ }
+ }
+ }
+ // return indicator that no new task could be created
+ return NO_TASK;
}
-
// set a status
void Agenda::set_status(Task t, ItemStatus status) {
// get the lane index
@@ -193,12 +185,11 @@ void Agenda::set_status(Task t, ItemStatus status) {
if ( (t.cycle > rlen) || (t.cycle == 0) ) {
throw std::out_of_range("Cycle out of range.");
}
- uint16_t cl_id = t.cycle -1;
+ uint16_t cl_id = getSeqCycle(t.cycle,t.seqEl.id) -1;
items[ln_id][tl_id][cl_id] = status;
}
-
// get the status of a task
ItemStatus Agenda::get_status(Task t) {
// get the lane index
@@ -224,7 +215,6 @@ ItemStatus Agenda::get_status(Task t) {
return items[ln_id][tl_id][cl_id];
}
-
// check if all items of the agenda were processed, if possible
bool Agenda::finished() {
// check for each tile if either all cycles are finished OR there is a failed status item
@@ -246,13 +236,47 @@ bool Agenda::finished() {
return true;
}
+// check if all items of the agenda were processed, if possible
+bool Agenda::finished( CountType cycle ) {
+ // check for each tile if either all cycles are finished OR there is a failed status item
+ for (uint16_t ln_id = 0; ln_id < items.size(); ++ln_id) {
+ for (uint16_t tl_id = 0; tl_id < items[ln_id].size(); ++tl_id) {
+ for (uint16_t cl_id = 0; cl_id < cycle; ++cl_id) {
+ ItemStatus s = items[ln_id][tl_id][cl_id];
+ if ( s == FAILED ) {
+ // the rest of the tile is "allowed" to be unprocessed --> skip
+ continue;
+ }
+ else if (s != FINISHED) {
+ // otherwise any other status means that the agenda is not finished
+ return false;
+ }
+ }
+ }
+ }
+ return true;
+}
+
+bool Agenda::cycle_available( CountType cycle ) {
+
+ if ( cycle == 0 || cycle > rlen )
+ return false;
+
+ for (uint16_t ln_id = 0; ln_id < items.size(); ++ln_id) {
+ for (uint16_t tl_id = 0; tl_id < items[ln_id].size(); ++tl_id) {
+ if ( items[ln_id][tl_id][cycle-1] == WAITING )
+ return false;
+ }
+ }
+
+ return true;
+}
// the total number of tasks on the agenda
uint32_t Agenda::task_count() {
return lanes.size() * tiles.size() * rlen;
}
-
// the total number of finished tasks on the agenda
uint32_t Agenda::tasks_finished() {
uint32_t num_finished = 0;
@@ -268,20 +292,3 @@ uint32_t Agenda::tasks_finished() {
}
return num_finished;
}
-
-
-// generate a complete TaskQueue with tasks to generate SAM files
-// SAM files can only be generated for tiles where all cycles are completed
-std::vector<Task> Agenda::get_SAM_tasks() {
- std::vector<Task> tv;
- // find tiles that are completely mapped
- for (uint16_t ln_id = 0; ln_id < items.size(); ++ln_id) {
- for (uint16_t tl_id = 0; tl_id < items[ln_id].size(); ++tl_id) {
- if ( items[ln_id][tl_id][rlen-1] == FINISHED ) {
- tv.push_back(Task(lanes[ln_id],tiles[tl_id],rlen,rlen,root));
- }
- }
- }
-
- return tv;
-}
diff --git a/lib/parallel.h b/lib/parallel.h
index c404592..912578d 100644
--- a/lib/parallel.h
+++ b/lib/parallel.h
@@ -8,27 +8,63 @@
//------ Threading tools --------------------------------------------//
-// Task data structure. Contains all information for a thread to
-// process a BCL file.
+/**
+ * Task data structure. Contains all information for a thread to process a BCL file.
+ * @author Martin Lindner
+ */
struct Task {
- // dataset information
+ /** The lane of the task. */
uint16_t lane;
+ /** The tile of the task. */
uint16_t tile;
+ /** Struct containing the read properties (Barcode vs. sequence; length; mate). */
+ SequenceElement seqEl;
+ /** Current cycle of the particular read (in general, this does NOT equal the sequencing cycle!). Must be <=seqEl.length. */
uint16_t cycle;
- uint16_t rlen;
- std::string root;
-
- // constructor initializes all member variables
- Task(uint16_t ln, uint16_t tl, uint16_t cl, uint16_t rl, std::string rt): lane(ln), tile(tl), cycle(cl), rlen(rl), root(rt) {};
- // overload << operator for printing
+ /**
+ * Constructor for a NULL task.
+ * @author Tobias Loka
+ */
+ Task() : lane(255), tile(0), seqEl(NULLSEQ), cycle(0) {};
+
+ /**
+ * Constructor for a valid task.
+ * @param ln The lane number.
+ * @param tl The tile number.
+ * @param seq The respective seqEl element for the current read containing information about length, type (barcode vs. sequence), mate number ...
+ * @param cl The cycle of the current read (in general, this does NOT equal the sequencing cycle!). Must be <=seqEl.length.
+ * @author Martin Lindner
+ */
+ Task(uint16_t ln, uint16_t tl, SequenceElement seq, uint16_t cl):
+ lane(ln), tile(tl), seqEl(seq), cycle(cl) {};
+
+ /**
+ * Overload of the << operator. Defines the cout form of a task.
+ * @author Martin Lindner
+ */
friend std::ostream& operator<<(std::ostream& os, const Task& t);
};
-inline bool operator==(const Task& l, const Task& r){ return (r.lane==l.lane)&&(r.tile==l.tile)&&(r.cycle==l.cycle)&&(r.rlen==l.rlen)&&(r.root==l.root); }
+/**
+ * Overload of the == operator.
+ * @return true, if all fields/variables of the compared tasks equal.
+ * @author Martin Lindner
+ */
+inline bool operator==(const Task& l, const Task& r){ return (r.lane==l.lane)&&(r.tile==l.tile)&&(r.cycle==l.cycle)&&(r.seqEl==l.seqEl); }
+
+/**
+ * Overload of the != operator.
+ * @return true, if at least one field/variable of the compared tasks is different.
+ * @author Martin Lindner
+ */
inline bool operator!=(const Task& l, const Task& r){ return !(l==r); }
-const Task NO_TASK (255,0,0,0,"");
+/**
+ * Definition of a NULL task.
+ * @author Martin Lindner
+ */
+const Task NO_TASK (255,0,NULLSEQ,0);
// Task queue data structure. Manages a list of task objects in a thread safe way.
class TaskQueue {
@@ -49,7 +85,6 @@ class TaskQueue {
uint64_t size();
};
-
// Agenda item status
typedef uint8_t ItemStatus;
const ItemStatus WAITING = 0;
@@ -70,20 +105,19 @@ class Agenda {
std::vector< std::vector< std::vector<ItemStatus> > > items;
// dataset information
- std::string root;
uint16_t rlen;
std::vector<uint16_t> lanes;
std::vector<uint16_t> tiles;
public:
- // initialize agenda with root directory and read length only (all lanes, all tiles)
- Agenda (std::string rt, uint16_t rl);
+ // initialize agenda with read length only (all lanes, all tiles)
+ Agenda (uint16_t rl);
- // initialize agenda with root directory, read length, and lanes (all tiles)
- Agenda (std::string rt, uint16_t rl, std::vector<uint16_t> ln);
+ // initialize agenda with read length and lanes (all tiles)
+ Agenda (uint16_t rl, std::vector<uint16_t> ln);
- // initialize agenda with root directory, read length, lanes, and tiles
- Agenda (std::string rt, uint16_t rl, std::vector<uint16_t> ln, std::vector<uint16_t> tl);
+ // initialize agenda with read length, lanes and tiles
+ Agenda (uint16_t rl, std::vector<uint16_t> ln, std::vector<uint16_t> tl);
// check for BCL files and update item status
void update_status();
@@ -100,6 +134,10 @@ class Agenda {
// check if all items of the agenda were processed, if possible
bool finished();
+ bool finished( CountType cycle );
+
+ bool cycle_available( CountType cycle );
+
// the total number of tasks on the agenda
uint32_t task_count();
@@ -111,8 +149,6 @@ class Agenda {
std::vector<Task> get_SAM_tasks();
};
-
-
// create a vector with all lane numbers
std::vector<uint16_t> all_lanes();
@@ -125,6 +161,4 @@ std::vector<uint16_t> all_tiles();
// create a vector with one tile number
std::vector<uint16_t> one_tile(uint16_t t);
-
-
#endif /* PARALLEL_H */
diff --git a/lib/tools.cpp b/lib/tools.cpp
index 34a1e72..1d6403c 100644
--- a/lib/tools.cpp
+++ b/lib/tools.cpp
@@ -1,144 +1,29 @@
#include "tools.h"
-// compares two genome positions by position (not by genome id)
-bool gp_compare (GenomePosType i,GenomePosType j) {
- return (i.pos < j.pos);
-}
-
-// reads a binary file from hdd and stores its raw content in a char vector
-std::vector<char> read_binary_file(const std::string &fname) {
-
- // get file size
- uint64_t size = get_filesize(fname);
-
- // open binary file
- FILE* f;
- f = fopen(fname.c_str(), "rb");
-
- if (!f) {
- std::cerr << "Error reading binary file " << fname << ": Could not open file." << std::endl;
- return std::vector<char>();
- }
-
- // allocate memory
- std::vector<char> data (size);
-
- // read all data at once
- uint64_t read = fread(data.data(), 1, size, f);
-
- if (read != size){
- std::cerr << "Error reading binary file " << fname << ": File size: " << size << " bytes. Read: " << read << " bytes." << std::endl;
- return std::vector<char>();
- }
-
- fclose(f);
-
- return data;
-}
-
-
-// checks if a directory with the given path exists
-bool is_directory(const std::string &path) {
- if ( boost::filesystem::exists(path) ) {
- if ( boost::filesystem::is_directory(path) ) {
- return true;
- }
- else {
- return false;
- }
- }
- else {
- return false;
- }
-}
-
-// checks if a file exists
-bool file_exists(const std::string &fname) {
- return boost::filesystem::exists(fname);
- /*if (FILE *f = fopen(fname.c_str(), "r")) {
- fclose(f);
- return true;
- }
- else {
- return false;
- }*/
-}
-
-
-// writes a char vector into a binary file
-uint64_t write_binary_file(const std::string &fname, const std::vector<char> & data) {
-
- // open binary file
- FILE* ofile;
- ofile = fopen(fname.c_str(), "wb");
-
- if (!ofile) {
- std::cerr << "Error serializing object to file " << fname << ": Could not open file for writing." << std::endl;
- return 0;
- }
-
- // write all data
- uint64_t written = fwrite(data.data(), 1, data.size(), ofile);
-
- // close file
- fclose(ofile);
-
- if (written != data.size()){
- std::cerr << "Error serializing object to file " << fname << ": Total size: " << data.size() << " bytes. Written: " << written << " bytes." << std::endl;
- }
-
- return written;
-}
-
-
-
-// extract the number of reads from a BCL file
-uint32_t num_reads_from_bcl(std::string bcl) {
- // open BCL file of first cycle
- FILE* ifile;
- ifile = fopen(bcl.c_str(), "rb");
-
- if (!ifile) {
- std::cerr << "Error reading BCL file " << bcl << ": Could not open file." << std::endl;
- return 0;
- }
-
- // extract the number of reads
- uint32_t num_reads;
- assert(fread(&num_reads, 1, sizeof(uint32_t), ifile));
-
- // close file
- fclose (ifile);
-
- return num_reads;
-}
-
-/* get the size (in bytes) of a file */
-std::ifstream::pos_type get_filesize(const std::string &fname)
-{
- std::ifstream in(fname, std::ios::binary | std::ios::ate);
- return in.tellg();
-}
+///////////////////////////////////
+////////// K-mer Hashing //////////
+///////////////////////////////////
-/* calculates the first forward and reverse complement k-mer in the
- string <kmer> and returns the canonical representation. */
-HashIntoType hash(const char * kmer, HashIntoType& _h, HashIntoType& _r, AlignmentSettings & settings)
+HashIntoType hash(const char * kmer, HashIntoType& _h, HashIntoType& _r)
{
- assert(strlen(kmer) >= settings.kmer_span);
+ assert(strlen(kmer) >= globalAlignmentSettings.get_kmer_span());
HashIntoType h = 0, r = 0;
h |= twobit_repr(kmer[0]);
- r |= twobit_comp(kmer[settings.kmer_span-1]);
+ r |= twobit_comp(kmer[globalAlignmentSettings.get_kmer_span()-1]);
- for (unsigned int i = 1, j = settings.kmer_span-2; i < settings.kmer_span; i++, j--) {
- // if i not gap position
- if (std::find(settings.kmer_gaps.begin(), settings.kmer_gaps.end(), i+1) == settings.kmer_gaps.end()) {
+ for (unsigned int i = 1, j = globalAlignmentSettings.get_kmer_span()-2; i < globalAlignmentSettings.get_kmer_span(); i++, j--) {
+ // if i not gap position
+
+ auto gaps_vec = globalAlignmentSettings.get_kmer_gaps();
+ if (std::find(gaps_vec.begin(), gaps_vec.end(), i+1) == gaps_vec.end()) {
h = h << 2;
h |= twobit_repr(kmer[i]);
r = r << 2;
r |= twobit_comp(kmer[j]);
+
}
}
@@ -148,25 +33,27 @@ HashIntoType hash(const char * kmer, HashIntoType& _h, HashIntoType& _r, Alignme
return (h)<(r)?h:r;
}
-/* calculates the first forward k-mer in the string <kmer> */
-std::string::const_iterator hash_fw(std::string::const_iterator it, std::string::const_iterator end, HashIntoType& _h, AlignmentSettings & settings)
+std::string::const_iterator hash_fw(std::string::const_iterator it, std::string::const_iterator end, HashIntoType& _h)
{
- assert(it+settings.kmer_span-1 < end);
+ if (!(it+globalAlignmentSettings.get_kmer_span()-1 < end)) {
+ std::cerr << "Error: hash_fw was called using an begin position which had not at least kmer_span bases behind it." << std::endl;
+ }
HashIntoType h = 0;
std::string::const_iterator last_invalid = it-1;
h |= twobit_repr(*it);
- std::string::const_iterator kmerEnd = it+settings.kmer_span;
+ std::string::const_iterator kmerEnd = it+globalAlignmentSettings.get_kmer_span();
++it;
int positionInKmer = 2;
+ auto kmer_gaps = globalAlignmentSettings.get_kmer_gaps();
for (; it != kmerEnd; ++it, ++positionInKmer) {
- if (std::find(settings.kmer_gaps.begin(), settings.kmer_gaps.end(), positionInKmer) != settings.kmer_gaps.end())
+ if (std::find(kmer_gaps.begin(), kmer_gaps.end(), positionInKmer) != kmer_gaps.end())
continue;
h = h << 2;
h |= twobit_repr(*it);
if ( seq_chars.find(*it) == std::string::npos ) {
- last_invalid = it+settings.kmer_span-1;
+ last_invalid = it+globalAlignmentSettings.get_kmer_span()-1;
}
}
@@ -174,8 +61,6 @@ std::string::const_iterator hash_fw(std::string::const_iterator it, std::string:
return last_invalid;
}
-
-/* returns the sequence of a k-mer */
std::string unhash(HashIntoType myHash, unsigned hashLen)
{
std::string kmer = "";
@@ -191,53 +76,134 @@ std::string unhash(HashIntoType myHash, unsigned hashLen)
-// file name construction functions
+////////////////////////////////////////////
+////////// File name construction //////////
+////////////////////////////////////////////
-// construct BCL file name from: root, lane, tile, cycle
-std::string bcl_name(std::string rt, uint16_t ln, uint16_t tl, uint16_t cl) {
+std::string bcl_name(uint16_t ln, uint16_t tl, uint16_t cl) {
std::ostringstream path_stream;
- path_stream << rt << "/L00" << ln << "/C" << cl << ".1/s_"<< ln <<"_" << tl << ".bcl";
+ path_stream << globalAlignmentSettings.get_root() << "/L00" << ln << "/C" << cl << ".1/s_"<< ln <<"_" << tl << ".bcl";
return path_stream.str();
}
-
-// construct alignment file name from: root, lane, tile, cycle
-std::string alignment_name(std::string rt, uint16_t ln, uint16_t tl, uint16_t cl) {
+std::string alignment_name(uint16_t ln, uint16_t tl, uint16_t cl, uint16_t mt){
std::ostringstream path_stream;
- path_stream << rt << "/L00" << ln << "/s_"<< ln << "_" << tl << "." << cl << ".align";
+ std::string base = globalAlignmentSettings.get_temp_dir() != "" ? globalAlignmentSettings.get_temp_dir() : globalAlignmentSettings.get_root();
+ path_stream << base << "/L00" << ln << "/s_"<< ln << "_" << tl << "." << mt << "."<< cl << ".align";
return path_stream.str();
}
-// construct tile-wise SAM file name from: root, lane, tile
-std::string sam_tile_name(std::string rt, uint16_t ln, uint16_t tl, bool write_bam) {
- std::ostringstream path_stream;
- if (write_bam)
- path_stream << rt << "/L00" << ln << "/s_"<< ln << "_" << tl << ".bam";
- else
- path_stream << rt << "/L00" << ln << "/s_"<< ln << "_" << tl << ".sam";
- return path_stream.str();
+uint16_t getSeqCycle(uint16_t cycle, uint16_t seq_id) {
+ uint16_t seq_cycle = cycle;
+ for ( int i = 0; i < seq_id; i++ )
+ seq_cycle += globalAlignmentSettings.getSeqById(i).length;
+ return seq_cycle;
}
-// construct lane-wise SAM file name from: root, lane
-std::string sam_lane_name(std::string rt, uint16_t ln, bool write_bam) {
- std::ostringstream path_stream;
- if (write_bam)
- path_stream << rt << "/L00" << ln << "/s_"<< ln << ".bam";
- else
- path_stream << rt << "/L00" << ln << "/s_"<< ln << ".sam";
- return path_stream.str();
+uint16_t getMateCycle( uint16_t mate_number, uint16_t seq_cycle ) {
+
+ // Invalid mate
+ if ( mate_number == 0 || mate_number > globalAlignmentSettings.get_mates() )
+ return 0;
+
+ // Iterate through all sequence elements (including barcodes)
+ for ( CountType id = 0; id < globalAlignmentSettings.get_seqs().size(); id++ ) {
+
+ // Current sequence element
+ SequenceElement seq = globalAlignmentSettings.getSeqById(id);
+
+ // Seq is mate of interest
+ if ( seq.mate == mate_number )
+ return ( seq.length > seq_cycle ? seq_cycle : seq.length );
+
+ // Not enough cycles left to reach mate of interest
+ else if ( seq.length >= seq_cycle )
+ return 0;
+
+ // Reduce number of cycles by the Seq length
+ else
+ seq_cycle -= seq.length;
+
+ }
+
+ // Should not be reached
+ return 0;
}
-// construct filter file name from: root, lane, tile
-std::string filter_name(std::string rt, uint16_t ln, uint16_t tl) {
+
+std::string filter_name(uint16_t ln, uint16_t tl) {
std::ostringstream path_stream;
- path_stream << rt << "/L00" << ln << "/s_"<< ln << "_" << tl << ".filter";
+ path_stream << globalAlignmentSettings.get_root() << "/L00" << ln << "/s_"<< ln << "_" << tl << ".filter";
return path_stream.str();
}
-// construct position file name from: root, lane, tile
-std::string position_name(std::string rt, uint16_t ln, uint16_t tl) {
+
+std::string position_name(uint16_t ln, uint16_t tl) {
std::ostringstream path_stream;
- path_stream << rt << "../L00" << ln << "/s_"<< ln << "_" << tl << ".clocs";
+ path_stream << globalAlignmentSettings.get_root() << "../L00" << ln << "/s_"<< ln << "_" << tl << ".clocs";
return path_stream.str();
}
+
+
+std::string get_settings_name() {
+ std::ostringstream path_stream;
+ std::string base = globalAlignmentSettings.get_temp_dir() != "" ? globalAlignmentSettings.get_temp_dir() : globalAlignmentSettings.get_root();
+ path_stream << base << "/hilive_settings.xml";
+ return path_stream.str();
+}
+
+
+std::string get_out_log_name() {
+ return ( globalAlignmentSettings.get_out_dir() + "/hilive_out.log" );
+}
+
+
+
+////////////////////////////////////
+////////// SAM/BAM output //////////
+////////////////////////////////////
+
+seqan::BamHeader getBamHeader() {
+ std::stringstream ss;
+ ss.str(std::string());
+ ss << HiLive_VERSION_MAJOR << "." << HiLive_VERSION_MINOR;
+
+ seqan::BamHeader header;
+ resize(header, 2);
+
+ // @HD header.
+ seqan::resize(header[0].tags, 2);
+ header[0].type = seqan::BAM_HEADER_FIRST;
+ header[0].tags[0].i1 = "VN";
+ header[0].tags[0].i2 = "1.5";
+ header[0].tags[1].i1 = "GO";
+ header[0].tags[1].i2 = "query";
+
+ // @PG header.
+ seqan::resize(header[1].tags, 3);
+ header[1].type = seqan::BAM_HEADER_PROGRAM;
+ header[1].tags[0].i1 = "ID";
+ header[1].tags[0].i2 = "hilive";
+ header[1].tags[1].i1 = "PN";
+ header[1].tags[1].i2 = "HiLive";
+ header[1].tags[2].i1 = "VN";
+ header[1].tags[2].i2 = ss.str();
+
+ return header;
+}
+
+
+std::string getBamTempFileName(std::string barcode, CountType cycle) {
+ std::ostringstream fname;
+ std::string file_suffix = globalAlignmentSettings.get_write_bam() ? ".bam" : ".sam";
+ fname << globalAlignmentSettings.get_out_dir() << "/hilive_out_" << "cycle" << std::to_string(cycle) << "_" << barcode << ".temp" << file_suffix;
+ return fname.str();
+}
+
+
+std::string getBamFileName(std::string barcode, CountType cycle) {
+ std::ostringstream fname;
+ std::string file_suffix = globalAlignmentSettings.get_write_bam() ? ".bam" : ".sam";
+ fname << globalAlignmentSettings.get_out_dir() << "/hilive_out_" << "cycle" << std::to_string(cycle) << "_" << barcode << file_suffix;
+ return fname.str();
+}
diff --git a/lib/tools.h b/lib/tools.h
index b8060c8..06a81b0 100644
--- a/lib/tools.h
+++ b/lib/tools.h
@@ -1,46 +1,147 @@
+/**
+ * This class provides functions that are dependent of the alignmentSettings!
+ * For functions that are not dependent of any other HiLive class, please use the tools_static class!
+ * Please do NOT add further includes to this file since this will lead to unwanted dependencies!
+ */
+
#ifndef TOOLS_H
#define TOOLS_H
-#include "headers.h"
-#include "definitions.h"
-#include "kindex.h"
+/* DONT ADD ANY INCLUDES */
+#include "tools_static.h"
+#include "alignmentSettings.h"
+#include "global_variables.h"
+/* DONT ADD ANY INCLUDES */
+
+
+///////////////////////////////////
+////////// K-mer Hashing //////////
+///////////////////////////////////
+
+/**
+ * Calculate the first forward and reverse complement k-mer in the string <kmer>.
+ * @param kmer Input sequence.
+ * @param _h Reference to forward hash variable.
+ * @param _r Reference to reverse hash variable.
+ * @return The larger hash value (TODO: why?)
+ */
+HashIntoType hash(const char * kmer, HashIntoType& _h, HashIntoType& _r);
+
+/**
+ * Calculates the first forward k-mer in the string <kmer>.
+ * @param it Iterator of the input sequence.
+ * @param end End of the iterator of the input sequence.
+ * @param _h Reference to the forward hash variable.
+ * @return Iterator pointing at the last invalid base.
+ */
+std::string::const_iterator hash_fw(std::string::const_iterator it, std::string::const_iterator end, HashIntoType& _h);
+
+/**
+ * Calculate the sequence from a hash value.
+ * @param myHash The input hash value.
+ * @param hashLen Length (weight) of the hashed sequence.
+ * @return The unhashed sequence.
+ */
+std::string unhash(HashIntoType myHash, unsigned hashLen=globalAlignmentSettings.get_kmer_weight());
+
+
+////////////////////////////////////////////
+////////// File name construction //////////
+////////////////////////////////////////////
+
+/**
+ * Get the name of a bcl file.
+ * @param ln The lane number.
+ * @param tl The tile number.
+ * @param cl The sequencing cycle.
+ * @return Path to the bcl file.
+ */
+std::string bcl_name(uint16_t ln, uint16_t tl, uint16_t cl);
+
+/**
+ * Get the name of an alignment file.
+ * @param ln The lane number.
+ * @param tl The tile number.
+ * @param cl The cycle for the respective mate.
+ * @param mt The mate number.
+ * @return Path to the alignment file.
+ */
+std::string alignment_name(uint16_t ln, uint16_t tl, uint16_t cl, uint16_t mt);
+/**
+ * Get the name of a filter file.
+ * @param ln The lane number.
+ * @param tl The tile number.
+ * @return Path to the filter file.
+ */
+std::string filter_name(uint16_t ln, uint16_t tl);
-// compare function to sort GenomePosType objects by position
-bool gp_compare (GenomePosType i,GenomePosType j);
+/**
+ * Get the name of a clocs file.
+ * @param ln The lane number.
+ * @param tl The tile number.
+ * @return Path to the clocs file.
+ */
+std::string position_name(uint16_t ln, uint16_t tl);
-// calculates the total size of a file in bytes
-std::ifstream::pos_type get_filesize(const std::string &fname);
+/**
+ * Get the name of the settings file.
+ * @return Path to the settings file.
+ */
+std::string get_settings_name();
-// checks if a directory with the given name exists
-bool is_directory(const std::string &path);
+/**
+ * Get the name of the output log file.
+ * @return Path to the output log file.
+ */
+std::string get_out_log_name();
-// checks if a file exists
-bool file_exists(const std::string &fname);
+/** Get the current sequencing cycle using the current alignment cycle and read number.
+ * @param cycle The read cycle.
+ * @param seq_id The sequence id (:= id of the respective element in globalAlignmentSettings::seqs)
+ * @return The sequencing cycle.
+ * @author Tobias Loka
+ */
+uint16_t getSeqCycle(uint16_t cycle, uint16_t seq_id=1);
-// reads a binary file from hdd and stores its raw content in a char vector
-std::vector<char> read_binary_file(const std::string &fname);
+/**
+ * Get the cycle of a mate for a given sequencing cycle.
+ * When the mate is completely finished in the given cycle, return its total sequence length.
+ * @param mate_number Mate of interest.
+ * @param seq_cycle The sequencing cycle.
+ * @return Cycle of the mate in the given sequencing cycle.
+ * @author Tobias Loka
+ */
+uint16_t getMateCycle( uint16_t mate_number, uint16_t seq_cycle );
-// extract the number of reads from a BCL file
-uint32_t num_reads_from_bcl(std::string bcl);
-// writes a char vector into a binary file
-uint64_t write_binary_file(const std::string &fname, const std::vector<char> & data);
+////////////////////////////////////
+////////// SAM/BAM output //////////
+////////////////////////////////////
-//------ Hashing helper functions ---------------------------------//
-HashIntoType hash(const char * kmer, HashIntoType& _h, HashIntoType& _r, AlignmentSettings & settings);
-std::string::const_iterator hash_fw(std::string::const_iterator it, std::string::const_iterator end, HashIntoType& _h, AlignmentSettings & settings);
-//HashIntoType rc(HashIntoType fw);
-/* returns the sequence of a k-mer */
-std::string unhash(HashIntoType myHash, unsigned hashLen=K_HiLive);
+/**
+ * Get the header for a SAM/BAM output file.
+ * @return The BAM header.
+ * @author Tobias Loka
+ */
+seqan::BamHeader getBamHeader();
-// file name construction functions
-std::string bcl_name(std::string rt, uint16_t ln, uint16_t tl, uint16_t cl);
-std::string alignment_name(std::string rt, uint16_t ln, uint16_t tl, uint16_t cl);
-std::string filter_name(std::string rt, uint16_t ln, uint16_t tl);
-std::string position_name(std::string rt, uint16_t ln, uint16_t tl);
-std::string sam_tile_name(std::string rt, uint16_t ln, uint16_t tl, bool write_bam);
-std::string sam_lane_name(std::string rt, uint16_t ln, bool write_bam);
+/**
+ * Name of a temporary SAM/BAM file (for the time it is written).
+ * @param barcode Barcode of the output file (or "undetermined" for undetermined reads)
+ * @param cycle The output cycle.
+ * @return Name of the temporary output file for writing.
+ * @author Tobias Loka
+ */
+std::string getBamTempFileName(std::string barcode, CountType cycle);
+/**
+ * Final name of a SAM/BAM file.
+ * @param barcode Barcode of the output file (or "undetermined" for undetermined reads)
+ * @param cycle The output cycle.
+ * @return Name of the final output file.
+ * @author Tobias Loka
+ */
+std::string getBamFileName(std::string barcode, CountType cycle);
#endif /* TOOLS_H */
diff --git a/lib/tools_static.cpp b/lib/tools_static.cpp
new file mode 100644
index 0000000..71d5f1a
--- /dev/null
+++ b/lib/tools_static.cpp
@@ -0,0 +1,181 @@
+#include "tools_static.h"
+
+
+/////////////////////////////////
+////////// Comparators //////////
+/////////////////////////////////
+
+bool gp_compare (GenomePosType i,GenomePosType j) {
+ if ( i.pos == j.pos )
+ return i.gid < j.gid;
+ return (i.pos < j.pos);
+}
+
+
+/////////////////////////////////////
+////////// Type convertion //////////
+/////////////////////////////////////
+
+void split(const std::string &s, char delim, std::vector<std::string> &elems) {
+ std::stringstream ss;
+ ss.str(s);
+ std::string item;
+ while (std::getline(ss, item, delim)) {
+ elems.push_back(item);
+ }
+}
+
+
+///////////////////////////////////
+////////// File handling //////////
+///////////////////////////////////
+
+std::ifstream::pos_type get_filesize(const std::string &fname)
+{
+ std::ifstream in(fname, std::ios::binary | std::ios::ate);
+ return in.tellg();
+}
+
+bool is_directory(const std::string &path) {
+ if ( boost::filesystem::exists(path) ) {
+ if ( boost::filesystem::is_directory(path) ) {
+ return true;
+ }
+ else {
+ return false;
+ }
+ }
+ else {
+ return false;
+ }
+}
+
+bool file_exists(const std::string &fname) {
+ return boost::filesystem::exists(fname);
+
+}
+
+std::string absolute_path(std::string fname) {
+ boost::filesystem::path input_path(fname);
+ return boost::filesystem::canonical(fname).string();
+}
+
+std::vector<char> read_binary_file(const std::string &fname) {
+
+ // get file size
+ uint64_t size = get_filesize(fname);
+
+ // open binary file
+ FILE* f;
+ f = fopen(fname.c_str(), "rb");
+
+ if (!f) {
+ std::cerr << "Error reading binary file " << fname << ": Could not open file." << std::endl;
+ return std::vector<char>();
+ }
+
+ // allocate memory
+ std::vector<char> data (size);
+
+ // read all data at once
+ uint64_t read = fread(data.data(), 1, size, f);
+
+ if (read != size){
+ std::cerr << "Error reading binary file " << fname << ": File size: " << size << " bytes. Read: " << read << " bytes." << std::endl;
+ return std::vector<char>();
+ }
+
+ fclose(f);
+
+ return data;
+}
+
+uint64_t write_binary_file(const std::string &fname, const std::vector<char> & data) {
+
+ // open binary file
+ FILE* ofile;
+ ofile = fopen(fname.c_str(), "wb");
+
+ if (!ofile) {
+ std::cerr << "Error serializing object to file " << fname << ": Could not open file for writing." << std::endl;
+ return 1;
+ }
+
+ // write all data
+ uint64_t written = fwrite(data.data(), 1, data.size(), ofile);
+
+ // close file
+ fclose(ofile);
+
+ if (written != data.size()){
+ std::cerr << "Error serializing object to file " << fname << ": Total size: " << data.size() << " bytes. Written: " << written << " bytes." << std::endl;
+ }
+
+ return written;
+}
+
+
+////////////////////////////////////////////////
+////////// Property trees / XML files //////////
+////////////////////////////////////////////////
+
+bool read_xml(boost::property_tree::ptree & xml_in, std::string xml_fname) {
+
+ if ( !file_exists(xml_fname) ) {
+ std::cout << "XML file not found: " << xml_fname << std::endl;
+ return false;
+ }
+
+ try {
+ boost::property_tree::read_xml (xml_fname, xml_in);
+ } catch ( const std::exception &ex) {
+ std::cerr << "Error loading xml file " << xml_fname << ": " << std::endl << ex.what() << std::endl;
+ return false;
+ }
+
+ return true;
+
+}
+
+bool write_xml(boost::property_tree::ptree & xml_out, std::string xml_fname) {
+
+ try {
+ boost::property_tree::write_xml( xml_fname, xml_out );
+ } catch ( const std::exception &ex ) {
+ std::cerr << "Error writing xml file " << xml_fname << ": " << std::endl << ex.what() << std::endl;
+ return false;
+ }
+
+ return true;
+
+}
+
+
+/////////////////////////////////
+////////// Other stuff //////////
+/////////////////////////////////
+
+uint32_t num_reads_from_bcl(std::string bcl) {
+ // open BCL file of first cycle
+ FILE* ifile;
+ ifile = fopen(bcl.c_str(), "rb");
+
+ if (!ifile) {
+ std::cerr << "Error reading BCL file " << bcl << ": Could not open file." << std::endl;
+ return 0;
+ }
+
+ // extract the number of reads
+ uint32_t num_reads;
+ bool res = fread(&num_reads, 1, sizeof(uint32_t), ifile);
+ if (!res) {
+ std::cerr << "Error extracting number of reads from BCL file " << bcl << std::endl;
+ return 0;
+ }
+
+ // close file
+ fclose (ifile);
+
+ return num_reads;
+}
+
diff --git a/lib/tools_static.h b/lib/tools_static.h
new file mode 100644
index 0000000..acab5dd
--- /dev/null
+++ b/lib/tools_static.h
@@ -0,0 +1,171 @@
+/**
+ * This class provides functions that are independent of any other HiLive class.
+ * Please do NOT add further includes to this file since this will lead to unwanted dependencies!
+ */
+
+#ifndef TOOLS_STATIC_H
+#define TOOLS_STATIC_H
+
+/* DONT ADD ANY INCLUDES */
+#include "headers.h"
+#include "definitions.h"
+/* DONT ADD ANY INCLUDES */
+
+
+/////////////////////////////////
+////////// Comparators //////////
+/////////////////////////////////
+
+/**
+ * Compare function to sort GenomePosType objects by position.
+ * If position if equal, compare by gid.
+ * @param i First position to compare
+ * @param j Second position to compare
+ * @return true, if first position is "smaller" than second position.
+ */
+bool gp_compare (GenomePosType i,GenomePosType j);
+
+
+/////////////////////////////////////
+////////// Type convertion //////////
+/////////////////////////////////////
+
+/**
+ * Split a std::string to a std::vector<std::string>.
+ * @param s Reference to the input string.
+ * @param delim A split delimiter.
+ * @param elems The target vector.
+ * @author Tobias Loka
+ */
+void split(const std::string &s, char delim, std::vector<std::string> &elems);
+
+
+///////////////////////////////////
+////////// File handling //////////
+///////////////////////////////////
+
+
+/**
+ * Get total size of a file (in bytes)
+ * @param fname Name of the file.
+ * @return Size of the file.
+ */
+std::ifstream::pos_type get_filesize(const std::string &fname);
+
+/**
+ * Check if a given path is a directory.
+ * @param Path of interest.
+ * @return true, if the given path is a directory.
+ */
+bool is_directory(const std::string &path);
+
+/**
+ * Check if a given path is a file.
+ * @param Path of interest.
+ * @return true, if the given path is a file.
+ */
+bool file_exists(const std::string &fname);
+
+/**
+ * Convert a relative to an absolute path.
+ * @param fname Input path.
+ * @return Absolute path to fname.
+ * @author Tobias Loka
+ * TODO: Not tested and used yet.
+ */
+std::string absolute_path(std::string fname);
+
+
+/**
+ * Read a binary file and stores its content in a char vector.
+ * @param fname Path to the file.
+ * @return All data from the file as char vector.
+ */
+std::vector<char> read_binary_file(const std::string &fname);
+
+/**
+ * Write data from a char vector into a binary file.
+ * @param fname Path to the file.
+ * @param data Data to be saved in the file.
+ * @return Number of written bytes.
+ */
+uint64_t write_binary_file(const std::string &fname, const std::vector<char> & data);
+
+
+////////////////////////////////////////////////
+////////// Property trees / XML files //////////
+////////////////////////////////////////////////
+
+/**
+ * Read a file in XML format. Results are stored as property tree.
+ * @param xml_in Reference to the property tree to store the XML data.
+ * @param xml_fname Name of the input file.
+ * @return true on success
+ * @author Tobias Loka
+ */
+bool read_xml(boost::property_tree::ptree & xml_in, std::string xml_fname);
+
+/**
+ * Write a property tree to an XML file.
+ * @param xml_out Property tree that contains the data.
+ * @param xml_fname Name of the output file.
+ * @return true on success
+ * @author Tobias Loka
+ */
+bool write_xml(boost::property_tree::ptree & xml_out, std::string xml_fname);
+
+/**
+ * Convert a variable of a non-vector type to a property tree.
+ * @param variable The variable to convert.
+ * @return The property tree for the input variable
+ * @author Tobias Loka
+ * TODO: check if the exception handling makes sense.
+ */
+/** Convert a variable to an XML node. T must be a data type that can be cast to a string-like output format. */
+template<typename T> boost::property_tree::ptree getXMLnode (T variable) {
+
+ boost::property_tree::ptree node;
+
+ try {
+ node.put("", variable);
+ } catch ( const std::exception &ex ) {
+ std::cerr << "Failed to convert variable to XML output format." << std::endl;
+ }
+
+ return node;
+
+}
+
+/**
+ * Convert a variable of a vector type to a property tree.
+ * The subnodes have key "el".
+ * @param vector The vector to convert.
+ * @return The property tree for the input variable
+ * @author Tobias Loka
+ */
+/** Convert a vector to an XML node. T must be a data type that can be cast to a string-like output format. */
+template<typename T> boost::property_tree::ptree getXMLnode_vector (std::vector<T> vector) {
+
+ boost::property_tree::ptree node;
+
+ for ( auto el = vector.begin(); el != vector.end(); ++el ) {
+ node.add_child("el", getXMLnode ( *el ));
+ }
+
+ return node;
+
+}
+
+
+/////////////////////////////////
+////////// Other stuff //////////
+/////////////////////////////////
+
+/**
+ * Extract the number of reads from a BCL file.
+ * @param bcl Path to the bcl file.
+ * @return Number of reads in the bcl file.
+ */
+uint32_t num_reads_from_bcl(std::string bcl);
+
+#endif /* TOOLS_STATIC_H */
diff --git a/tools/build_index.cpp b/tools/build_index.cpp
index 090d886..7f4fbe3 100644
--- a/tools/build_index.cpp
+++ b/tools/build_index.cpp
@@ -3,124 +3,55 @@
#include "../lib/headers.h"
#include "../lib/definitions.h"
#include "../lib/kindex.h"
+#include "../lib/argument_parser.h"
-namespace po = boost::program_options;
+AlignmentSettings globalAlignmentSettings;
-std::string license =
-"Copyright (c) 2015-2016, Martin S. Lindner and the HiLive contributors. See CONTRIBUTORS for more info.\n"
-"All rights reserved.\n"
-"\n"
-"Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:\n"
-"\n"
-"1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.\n"
-"\n"
-"2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.\n"
-"\n"
-"3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.\n"
-"\n"
-"THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 'AS IS' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PRO [...]
+/**
+ * Main function that organizes the overall structure of the program.
+ * @param argc Number of arguments
+ * @param argv Argument array
+ * @return 0 on success, other numbers on error
+ */
+int main(int argc, const char* argv[]) {
+ // Program start output.
+ std::cout << std::endl << "------" << std::endl << "HiLive Index Builder v" << HiLive_VERSION_MAJOR << "." << HiLive_VERSION_MINOR <<
+ " - Build Index for Realtime Alignment of Illumina Reads" << std::endl << "------" << std::endl << std::endl;
+ // Init argument parser
+ BuildIndexArgumentParser argumentParser(argc, argv);
+ // Parse command line arguments
+ int parser_returnStatus = argumentParser.parseCommandLineArguments();
-int main(int argc, char* argv[]) {
+ // Successful execution of "help" or "license"
+ if ( parser_returnStatus == 1 ) {
+ exit(EXIT_SUCCESS);
+ }
- unsigned trim;
- bool do_not_convert_spaces;
- bool trim_ids;
+ // Parsing error
+ else if ( parser_returnStatus == -1 ) {
+ std::cout << "Parsing of command line options failed. For help, type 'hilive-build --help'." << std::endl;
+ exit(EXIT_FAILURE);
+ }
- // setting up the command line interface
- po::options_description general("General");
- general.add_options()
- ("help,h", "Print this help message and exit")
- ("license", "Print licensing information and exit");
+ // Build index
+ std::cout << "Create index from file " << argumentParser.fasta_name << " ..." << std::endl;
+ KixBuild* index = new KixBuild();
+ index->add_fasta(argumentParser.fasta_name, !argumentParser.do_not_convert_spaces, argumentParser.trim_ids);
- po::options_description parameters("Parameters");
- parameters.add_options()
- ("INPUT", po::value<std::string>()->required(), "Input reference genome (fasta file)");
+ // Trim index
+ if ( argumentParser.trim > 0) {
+ uint64_t trimmed = index->trim(argumentParser.trim);
+ std::cout << "Removed " << trimmed << " k-mer positions from the database." << std::endl;
+ }
- po::options_description options("Options");
- options.add_options()
- ("outfile,o", po::value<std::string>(), "Set output file name [Default: INPUT.kix]")
- ("trim,t", po::value<unsigned>(&trim)->default_value(0), "Ignore k-mers with more than t occurrences. [Default: no limit]")
- ("do-not-convert-spaces", po::bool_switch(&do_not_convert_spaces)->default_value(false), "Do not convert all spaces in reference ids to underscores [Default: converting is on]")
- ("trim-after-space", po::bool_switch(&trim_ids)->default_value(false), "Trim all reference ids after first space [Default: false]");
+ // Write index to file
+ std::cout << "Writing index to file " << argumentParser.index_name << std::endl;
+ index->serialize_file(argumentParser.index_name);
- po::options_description cmdline_options;
- cmdline_options.add(general).add(parameters).add(options);
+ delete index;
- po::options_description visible_options;
- visible_options.add(general).add(options);
-
- std::stringstream help_message;
- help_message << "HiLive index builder v"<< HiLive_VERSION_MAJOR << "." << HiLive_VERSION_MINOR << std::endl;
- help_message << "Index creation tool for HiLive - Realtime Alignment of Illumina Reads" << std::endl;
- help_message << "Copyright (c) 2015, Martin S. Lindner" << std::endl;
- help_message << "HiLive is open-source software. Check with --license for details." << std::endl << std::endl;
- help_message << "Fixed k-mer size: " << K_HiLive << std::endl << std::endl;
- help_message << "Usage: " << std::string(argv[0]) << " [options] INPUT" << std::endl;
- help_message << " INPUT Reference genomes in (multi-) FASTA format" << std::endl;
-
- help_message << visible_options << std::endl;
-
- po::positional_options_description p;
- p.add("INPUT", 1);
-
- po::variables_map vm;
- try {
- // parse arguments
- po::store(po::command_line_parser(argc, argv).
- options(cmdline_options).positional(p).run(), vm);
- // first check if -h or --help was called
- if (vm.count("help")) {
- std::cout << help_message.str();
- return 1;
- }
- // first check if --license was called
- if (vm.count("license")) {
- std::cout << license << std::endl;
- return 1;
- }
-
- // then check arguments
- po::notify(vm);
- }
- catch ( boost::program_options::required_option& e ) {
- std::cerr << "Missing Parameter: " << e.what() << std::endl << std::endl;
- std::cout << help_message.str();
- return -1;
- }
- catch( boost::program_options::error& e) {
- std::cerr << "Error while parsing command line options: " << e.what() << std::endl << std::endl;
- std::cout << help_message.str();
- return -1;
- }
-
-
- std::string fasta_name = vm["INPUT"].as<std::string>();
-
- std::string index_name;
- if (vm.count("outfile")) {
- // use the output file name if provided as argument
- index_name = vm["outfile"].as<std::string>();
- } else {
- // construct it from the input file name otherwise
- index_name = fasta_name + std::string(".kix");
- }
-
- std::cout << "Creating index with K_HiLive=" << K_HiLive << " from file " << fasta_name << std::endl;
- KixBuild* index = new KixBuild();
- AlignmentSettings settings; // for hard coded gapped kmer structure
- index->add_fasta(fasta_name, settings, !do_not_convert_spaces, trim_ids);
-
- if (trim > 0) {
- uint64_t trimmed = index->trim(trim);
- std::cout << "Removed " << trimmed << " k-mer positions from the database." << std::endl;
- }
-
-
- std::cout << "Writing index to file " << index_name << std::endl;
- index->serialize_file(index_name);
-
- delete index;
+ return EXIT_SUCCESS;
}
diff --git a/tools/hilive.cpp b/tools/hilive.cpp
index feb7349..4239082 100644
--- a/tools/hilive.cpp
+++ b/tools/hilive.cpp
@@ -1,47 +1,63 @@
#include "../lib/headers.h"
#include "../lib/definitions.h"
+#include "../lib/global_variables.h"
#include "../lib/kindex.h"
#include "../lib/alnstream.h"
#include "../lib/parallel.h"
#include "../lib/argument_parser.h"
-std::string license =
-"Copyright (c) 2015-2016, Martin S. Lindner and the HiLive contributors. See CONTRIBUTORS for more info.\n"
-"All rights reserved.\n"
-"\n"
-"Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:\n"
-"\n"
-"1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.\n"
-"\n"
-"2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.\n"
-"\n"
-"3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.\n"
-"\n"
-"THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 'AS IS' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PRO [...]
+AlignmentSettings globalAlignmentSettings;
+/**
+ * Worker function for the alignment threads.
+ * @param tasks Reference to the "to do" task queue
+ * @param finished Reference to the "finished" task queue
+ * @param failed Reference to the "failed" task queue
+ * @param idx Pointer to the index object
+ * @param surrender Control flag (threads stop if true)
+ */
+void worker (TaskQueue & tasks, TaskQueue & finished, TaskQueue & failed, KixRun* idx, bool & surrender ) {
-
-
-// the worker function for the threads
-void worker (TaskQueue & tasks, TaskQueue & finished, TaskQueue & failed, AlignmentSettings* settings, KixRun* idx, bool & surrender ) {
-
- // loop that keeps on running until the surrender flag is set
+ // Continue until surrender flag is set
while ( !surrender ) {
- // try to obtain a new task
+
+ // Try to obtain a new task
Task t = tasks.pop();
+
+ // If "to do" task was found
if ( t != NO_TASK ) {
+
// Execute the task
bool success = true;
try {
- StreamedAlignment s (t.lane, t.tile, t.root, t.rlen);
+
+ StreamedAlignment s (t.lane, t.tile, t.seqEl.length);
uint64_t num_seeds;
- num_seeds = s.extend_alignment(t.cycle,idx,settings);
- std::cout << "Task [" << t << "]: Found " << num_seeds << " seeds." << std::endl;
+
+ // Seed extension if current read is sequence fragment.
+ if ( !t.seqEl.isBarcode() ) {
+ num_seeds = s.extend_alignment(t.cycle,t.seqEl.id,t.seqEl.mate,idx);
+ std::cout << "Task [" << t << "]: Found " << num_seeds << " seeds." << std::endl;
+
+ }
+
+ // Barcode extension if current read is barcode fragment
+ else {
+ CountType mate = 1;
+ for ( ; mate <= globalAlignmentSettings.get_mates(); mate++ ) {
+ SequenceElement seqEl = globalAlignmentSettings.getSeqByMate(mate);
+ CountType current_mate_cycle = t.seqEl.id < seqEl.id ? 0 : seqEl.length;
+ s.extend_barcode(t.cycle, current_mate_cycle, t.seqEl.id, mate);
+ }
+ std::cout << "Task [" << t << "]: Extended barcode of " << --mate << " mates." << std::endl;
+ }
}
catch (const std::exception &e) {
std::cerr << "Failed to finish task [" << t << "]: " << e.what() << std::endl;
success = false;
}
+
+ // Push the task in the correct Task Queue (Finished or Failed)
if (success) {
finished.push(t);
}
@@ -55,106 +71,161 @@ void worker (TaskQueue & tasks, TaskQueue & finished, TaskQueue & failed, Alignm
std::this_thread::sleep_for (std::chrono::milliseconds(100));
}
}
-
}
+/**
+ * Worker function for the output thread.
+ * Should only be called by one single thread!
+ * @param agenda Reference to the agenda that organizes the tasks
+ * @param idx Pointer to the index object
+ * @param surrender Control flag (threads stop if true)
+ * @author Tobias Loka
+ */
+void output_worker( Agenda & agenda, KixRun* idx, bool & surrender ) {
+
+ // Get the output cycles (must be sorted!)
+ std::vector<CountType> output_cycles = globalAlignmentSettings.get_output_cycles();
+
+ // Delayed surrender
+ bool alignments_finished = false;
+
+ // Continue as long as there are output cycles left
+ while ( output_cycles.size() > 0 ) {
+
+ alignments_finished = surrender;
+
+ // Create output for the next cycle if all related tasks are finished.
+ CountType next_cycle = *(output_cycles.begin());
+ if ( agenda.finished(next_cycle) ) {
+
+ try {
+ if ( !alignments_to_sam(globalAlignmentSettings.get_lanes(), globalAlignmentSettings.get_tiles(), idx, next_cycle) )
+ std::cerr << "Writing output for cycle " << std::to_string(next_cycle) << " failed." << std::endl;
+ else {
+ std::cout << "Wrote output for cycle " << std::to_string(next_cycle) << "." << std::endl;
+ }
+ }
+ catch ( std::exception & e ) {
+ std::cerr << "Writing output for cycle " << std::to_string(next_cycle) << " failed: " << e.what() << std::endl;
+ }
+
+ // Remove current cycle from the list of output cycles.
+ output_cycles.erase(output_cycles.begin());
+ }
+
+ // If the the surrender flag set but no cycle was handled, remove the cycle from the list
+ if ( alignments_finished && output_cycles.size() > 0 && next_cycle == *(output_cycles.begin()) ) {
+ std::cerr << "Writing output for cycle " << std::to_string(next_cycle) << " failed: Not all tasks finished." << std::endl;
+ output_cycles.erase(output_cycles.begin());
+ }
+
+ // Sleep a second
+ std::this_thread::sleep_for (std::chrono::milliseconds(1000));
+
+ }
-// create a special SAM worker, that writes out a SAM file for a tile
-void sam_worker (TaskQueue & tasks, AlignmentSettings* settings, KixRun* idx) {
+}
- // loop that keeps on running until the surrender flag is set
- while ( true ) {
- // try to obtain a new task
- Task t = tasks.pop();
- if ( t != NO_TASK ) {
- // Execute the task
- alignments_to_sam(t.lane,t.tile,t.root,t.rlen, idx,settings);
- }
- else {
- return;
- }
- }
+/**
+ * Main function that organizes the overall structure of the program.
+ * @param argc Number of arguments
+ * @param argv Argument array
+ * @return 0 on success, other numbers on error
+ */
+int main(int argc, const char* argv[]) {
-}
+ // Variable for runtime measurement
+ time_t t_start = time(NULL);
+ // Program start output
+ std::cout << std::endl << "------" << std::endl << "HiLive v"<< HiLive_VERSION_MAJOR << "." << HiLive_VERSION_MINOR <<
+ " - Realtime Alignment of Illumina Reads" << std::endl << "------" << std::endl<< std::endl;
+ // Parse command line arguments
+ HiLiveArgumentParser argumentParser(argc, argv);
+ int parser_returnStatus = argumentParser.parseCommandLineArguments();
-int main(int argc, const char* argv[]) {
- time_t t_start = time(NULL);
+ // Successful execution of "help" or "license"
+ if ( parser_returnStatus == 1 ) {
+ exit(EXIT_SUCCESS);
+ }
- // parse the command line arguments, store results in settings
- AlignmentSettings settings;
- if (parseCommandLineArguments(settings, license, argc, argv)) // returns true if error occured
- return 1;
+ // Parsing error
+ else if ( parser_returnStatus == -1 ) {
+ std::cout << "Parsing of command line options failed. For help, type 'hilive --help'." << std::endl;
+ exit(EXIT_FAILURE);
+ }
- // load the index
- std::cout << "Loading Index" << std::endl;
+ // Load the index
+ std::cout << "Loading Index ... " << std::endl;
KixRun* index = new KixRun();
- index->deserialize_file(settings.index_fname);
+ index->deserialize_file(globalAlignmentSettings.get_index_fname());
+
+ // Report loaded k-mer properties
+ std::cout << std::endl;
+ std::cout << "kmer span: " << std::to_string(globalAlignmentSettings.get_kmer_span()) << std::endl;
+ std::cout << "kmer weight: " << std::to_string(globalAlignmentSettings.get_kmer_weight()) << std::endl;
+ std::cout << "kmer gaps: ";
+ for ( auto gap : globalAlignmentSettings.get_kmer_gaps() ) {
+ std::cout << gap << " ";
+ }
+ std::cout << std::endl << std::endl;
+ // Write the alignment settings to an XML file
+ boost::property_tree::ptree xml_out = globalAlignmentSettings.to_ptree();
+ if ( ! write_xml(xml_out, get_settings_name()) )
+ exit(EXIT_FAILURE);
// Create the overall agenda
- Agenda agenda (settings.root, settings.rlen, settings.lanes, settings.tiles);
-
+ Agenda agenda (globalAlignmentSettings.get_cycles(), globalAlignmentSettings.get_lanes(), globalAlignmentSettings.get_tiles());
- // prepare the alignment
- std::cout << "Initializing Alignment files. Waiting for the first cycle to finish." << std::endl;
- bool first_cycle_available = false;
-
- // wait for the first cycle to be written. Attention - this loop will wait infinitely long if no first cycle is found
- while ( !first_cycle_available ) {
- // check for new BCL files and update the agenda status
+ // Wait for the first cycle to be written
+ std::cout << "Waiting for the first cycle to finish..." << std::endl;
+ while ( ! agenda.cycle_available(1) ) {
agenda.update_status();
+ std::this_thread::sleep_for(std::chrono::milliseconds(1000));
+ }
- // check if the first cycle is available for all tiles
- first_cycle_available = true;
- for ( auto ln : settings.lanes ) {
- for ( auto tl : settings.tiles ) {
- if ( agenda.get_status(Task(ln,tl,1,settings.rlen,"")) != BCL_AVAILABLE) {
- first_cycle_available = false;
- }
+ // Write empty alignment file for each tile and for each sequence read
+ std::cout << "Initializing Alignment files..." << std::endl;
+ for (uint16_t ln : globalAlignmentSettings.get_lanes()) {
+ for (uint16_t tl : globalAlignmentSettings.get_tiles()) {
+ CountType mate = 1;
+ for ( ; mate <= globalAlignmentSettings.get_mates(); mate++ ) {
+ StreamedAlignment s (ln, tl, globalAlignmentSettings.getSeqByMate(mate).length);
+ s.create_directories();
+ s.init_alignment(mate);
}
}
-
- // take a small break
- std::this_thread::sleep_for(std::chrono::milliseconds(1000));
}
std::cout << "First cycle complete. Starting alignment." << std::endl;
- // write empty alignment file for each tile
- for (uint16_t ln : settings.lanes) {
- for (uint16_t tl : settings.tiles) {
- StreamedAlignment s (ln, tl, settings.root, settings.rlen);
- s.create_directories(&settings);
- s.init_alignment(&settings);
- }
- }
- if (settings.temp_dir != "" && !is_directory(settings.temp_dir)){
- std::cerr << "Error: Could not find temporary directory " << settings.temp_dir << std::endl;
- return -1;
- }
// Set up the queues
TaskQueue toDoQ;
TaskQueue finishedQ;
TaskQueue failedQ;
- // Create the threads
- std::cout << "Creating " << settings.num_threads << " threads." << std::endl;
+ // Create the alignment threads
+ std::cout << "Creating " << globalAlignmentSettings.get_num_threads() << " threads." << std::endl;
bool surrender = false;
std::vector<std::thread> workers;
- for (int i = 0; i < settings.num_threads; i++) {
- workers.push_back(std::thread(worker, std::ref(toDoQ), std::ref(finishedQ), std::ref(failedQ), &settings, index, std::ref(surrender)));
+ for (int i = 0; i < globalAlignmentSettings.get_num_threads(); i++) {
+ workers.push_back(std::thread(worker, std::ref(toDoQ), std::ref(finishedQ), std::ref(failedQ), index, std::ref(surrender)));
}
+ // Create the output thread
+ std::thread output_thread = std::thread( output_worker, std::ref(agenda), index, std::ref(surrender) );
+
// Process all tasks on the agenda
while ( !agenda.finished() ) {
+
// check for new BCL files and update the agenda status
agenda.update_status();
- // fill the ToDo queue with tasks from the agenda
+ // fill the To Do queue with tasks from the agenda
while(true) {
Task t = agenda.get_task();
if (t == NO_TASK)
@@ -196,28 +267,12 @@ int main(int argc, const char* argv[]) {
for (auto& w : workers) {
w.join();
}
- std::cout << "All threads joined." << std::endl;
-
-
- std::cout << "Writing SAM files." << std::endl;
- // Create individual SAM files for every tile
- TaskQueue sam_tasks;
- std::vector<Task> tv = agenda.get_SAM_tasks();
- for ( auto t: tv ) {
- sam_tasks.push(t);
- }
-
- workers.clear();
- for (int i = 0; i < settings.num_threads; i++) {
- workers.push_back(std::thread(sam_worker, std::ref(sam_tasks), &settings, index));
- }
-
- for (auto& w : workers) {
- w.join();
- }
+ output_thread.join();
+ std::cout << "All threads joined." << std::endl;
+ std::cout << "Total mapping time: " << time(NULL) - t_start << " s" << std::endl << std::endl;
+ delete index;
std::cout << "Total run time: " << time(NULL) - t_start << " s" << std::endl;
- delete index;
- return 1;
+ exit(EXIT_SUCCESS);
}
diff --git a/tools/hilive_out.cpp b/tools/hilive_out.cpp
new file mode 100644
index 0000000..5e24b3c
--- /dev/null
+++ b/tools/hilive_out.cpp
@@ -0,0 +1,66 @@
+#include <boost/program_options.hpp>
+
+#include "../lib/headers.h"
+#include "../lib/definitions.h"
+#include "../lib/kindex.h"
+#include "../lib/alnstream.h"
+#include "../lib/parallel.h"
+#include "../lib/argument_parser.h"
+
+namespace po = boost::program_options;
+
+AlignmentSettings globalAlignmentSettings;
+
+/**
+ * Main function that organizes the overall structure of the program.
+ * @param argc Number of arguments
+ * @param argv Argument array
+ * @return 0 on success, other numbers on error
+ */
+int main(int argc, const char* argv[]) {
+
+ // Program start output
+ std::cout << std::endl << "------" << std::endl << "HiLive Output Tool v"<< HiLive_VERSION_MAJOR << "." << HiLive_VERSION_MINOR <<
+ " - Output of Realtime Alignments of Illumina Reads" << std::endl << "------" << std::endl<< std::endl;
+
+ // Parse the command line arguments
+ HiLiveOutArgumentParser argumentParser(argc, argv);
+ int parser_returnStatus = argumentParser.parseCommandLineArguments();
+
+ // Successful execution of "help" or "license"
+ if ( parser_returnStatus == 1 ) {
+ exit(EXIT_SUCCESS);
+ }
+
+ // Parsing error
+ else if ( parser_returnStatus == -1 ) {
+ std::cout << "Parsing of command line options failed. For help, type 'hilive-out --help'." << std::endl;
+ exit(EXIT_FAILURE);
+ }
+
+ // load the index
+ std::cout << "Loading Index Header..." << std::endl;
+ KixRun* index = new KixRun();
+
+ index->get_header_information(globalAlignmentSettings.get_index_fname());
+ index->store_kmer();
+
+ std::cout << "Start writing ouput." << std::endl;
+
+ for ( CountType cycle : globalAlignmentSettings.get_output_cycles() ) {
+ try {
+ if ( alignments_to_sam(globalAlignmentSettings.get_lanes(), globalAlignmentSettings.get_tiles(), index, cycle) )
+ std::cout << "Cycle " << std::to_string(cycle) << " ... " << "success." << std::endl;
+ else
+ std::cout << "Cycle " << std::to_string(cycle) << " ... " << "failed." << std::endl;
+ } catch ( std::exception & ex ) {
+ std::cout << "Cycle " << std::to_string(cycle) << " ... " << "failed: " << std::endl << ex.what() << std::endl << std::endl;
+ }
+ }
+
+ std::cout << "Finished." << std::endl;
+
+ delete index;
+
+ return EXIT_SUCCESS;
+}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/hilive.git
More information about the debian-med-commit
mailing list