Merge pull request #113 from eldariont/FEATURE/improve_junction_detection

[FEATURE] Improve junction detection (insertion calling 1/2)

Author: eldariont

include/iGenVar.hpp

@@ -12,9 +12,10 @@ struct cmd_arguments
     std::vector<detection_methods> methods{cigar_string, split_read, read_pairs, read_depth};   // default: all methods
     clustering_methods clustering_method{hierarchical_clustering};                              // default: hierarchical clustering method
     refinement_methods refinement_method{no_refinement};                                        // default: no refinement
-    uint64_t min_var_length = 30;
-    uint64_t max_var_length = 1000000;
-    uint64_t max_tol_inserted_length = 5;
+    int32_t min_var_length = 30;
+    int32_t max_var_length = 1000000;
+    int32_t max_tol_inserted_length = 5;
+    int32_t max_overlap = 10;
 };
 
 void initialize_argument_parser(seqan3::argument_parser & parser, cmd_arguments & args);
@@ -37,9 +38,10 @@ void initialize_argument_parser(seqan3::argument_parser & parser, cmd_arguments
  *                      (0: no_refinement,
  *                       1: sViper_refinement_method,
  *                       2: sVirl_refinement_method) - *default: no refinement*\n
- *                   **args.min_var_length** - minimum length of variants to detect - *default: 30 bp*\n
- *                   **args.max_var_length** - maximum length of variants to detect - *default: 1,000,000 bp*\n
- *                   **args.max_tol_inserted_length** - longest tolerated inserted sequence at non-INS SV types - *default: 5 bp*
+ *                   **args.min_var_length** - minimum length of variants to detect (expected to be non-negative) - *default: 30 bp*\n
+ *                   **args.max_var_length** - maximum length of variants to detect (expected to be non-negative) - *default: 1,000,000 bp*\n
+ *                   **args.max_tol_inserted_length** - longest tolerated inserted sequence at non-INS SV types (expected to be non-negative) - *default: 5 bp*\n
+ *                   **args.max_overlap** - maximum overlap between alignment segments (expected to be non-negative) - *default: 10 bp*
  *
  *
  * \details Detects novel junctions from read alignment records using different detection methods.

include/modules/sv_detection_methods/analyze_cigar_method.hpp

@@ -10,7 +10,7 @@
  * \param[in]       cigar_string    - CIGAR field of the SAM/BAM file
  * \param[in]       query_sequence  - SEQ field of the SAM/BAM file
  * \param[in, out]  junctions       - vector for storing junctions
- * \param[in]       min_length      - minimum length of variants to detect (default 30 bp)
+ * \param[in]       min_length      - minimum length of variants to detect (default 30 bp, expected to be non-negative)
  *
  * \details This function steps through the CIGAR string and stores junctions with their position in reference and read.
  *          We distinguish 4 cases of CIGAR operation characters:
@@ -36,4 +36,4 @@ void analyze_cigar(std::string const & read_name,
                    std::vector<seqan3::cigar> & cigar_string,
                    seqan3::dna5_vector const & query_sequence,
                    std::vector<Junction> & junctions,
-                   uint64_t const min_length);
+                   int32_t const min_length);

include/modules/sv_detection_methods/analyze_sa_tag_method.hpp

@@ -1,6 +1,9 @@
 #pragma once
 
-#include "structures/aligned_segment.hpp"   // for struct AlignedSegment
+#include "iGenVar.hpp"                          // for struct cmd_arguments
+#include "structures/aligned_segment.hpp"       // for struct AlignedSegment
+#include "structures/junction.hpp"              // for class Junction
+#include "variant_detection/bam_functions.hpp"  // for seqan3::sam_flag and hasFlag* functions
 
 /*! \brief Splits a string by a given delimiter and stores substrings in a given container.
  *
@@ -39,11 +42,15 @@ void retrieve_aligned_segments(std::string const & sa_string, std::vector<Aligne
  * \param[in, out]  junctions           - vector for storing junctions
  * \param[in, out]  query_sequence      - SEQ field of the SAM/BAM file
  * \param[in]       read_name           - QNAME field of the SAM/BAM file
+ * \param[in]       min_length          - minimum length of variants to detect (expected to be non-negative)
+ * \param[in]       max_overlap         - maximum overlap between alignment segments (expected to be non-negative)
  */
 void analyze_aligned_segments(std::vector<AlignedSegment> const & aligned_segments,
                               std::vector<Junction> & junctions,
                               seqan3::dna5_vector const & query_sequence,
-                              std::string const & read_name);
+                              std::string const & read_name,
+                              int32_t const min_length,
+                              int32_t const max_overlap);
 
 /*! \brief Parse the SA tag from the SAM/BAM alignment of a chimeric/split-aligned read. Build
  *         [aligned_segments](\ref AlignedSegment), one for each alignment segment of the read.
@@ -57,6 +64,9 @@ void analyze_aligned_segments(std::vector<AlignedSegment> const & aligned_segmen
  * \param[in]       cigar       - CIGAR field of the SAM/BAM file
  * \param[in]       seq         - SEQ field of the SAM/BAM file
  * \param[in]       sa_tag      - SA tag, one tag from the read of the SAM/BAM file
+ * \param[in]       args        - command line arguments:\n
+ *                                **args.min_var_length** - minimum length of variants to detect (expected to be non-negative)\n
+ *                                **args.max_overlap** - maximum overlap between alignment segments (expected to be non-negative)
  * \param[in, out]  junctions   - vector for storing junctions
  */
 void analyze_sa_tag(std::string const & query_name,
@@ -67,4 +77,5 @@ void analyze_sa_tag(std::string const & query_name,
                     std::vector<seqan3::cigar> const & cigar,
                     seqan3::dna5_vector const & seq,
                     std::string const & sa_tag,
+                    cmd_arguments const & args,
                     std::vector<Junction> & junctions);

include/structures/aligned_segment.hpp

@@ -21,8 +21,10 @@ struct AlignedSegment
     int32_t mapq;
     std::vector<seqan3::cigar> cig;
 
+    //! \brief Returns the reference position of the first aligned based.
     int32_t get_reference_start() const;
 
+    //! \brief Returns the reference position of the last aligned based.
     int32_t get_reference_end() const;
 
     int32_t get_left_soft_clip() const;

include/variant_detection/variant_detection.hpp

@@ -3,47 +3,43 @@
 #include <seqan3/std/filesystem>    // for filesystem
 #include <vector>
 
-#include "method_enums.hpp"         // for enum detection_methods, clustering_methods and refinement_methods
+#include "iGenVar.hpp"              // for struct cmd_arguments
 #include "structures/junction.hpp"  // for class Junction
 
 /*! \brief Detects junctions between distant genomic positions by analyzing a short read alignment file (sam/bam). The
  *         detected junctions are stored in a vector.
  *
  * \param[in, out]  junctions - a vector of junctions
- * \param[in]       alignment_short_reads_file_path - short reads input file, path to the sam/bam file
- * \param[in]       methods - list of methods for detecting junctions (0: cigar_string,
- *                                                                     1: split_read,
- *                                                                     2: read_pairs,
- *                                                                     3: read_depth)
- * \param[in]       min_var_length - minimum length of variants to detect (default 30 bp)
+ * \param[in]       args - command line arguments:\n
+ *                         **args.alignment_short_reads_file_path** - short reads input file, path to the sam/bam file\n
+ *                         **args.methods** - list of methods for detecting junctions
+ *                            (0: cigar_string, 1: split_read, 2: read_pairs, 3: read_depth) - *default: all methods*\n
+ * 
  *
  * \details Detects junctions from the CIGAR strings and supplementary alignment tags of read alignment records.
  *          We filter unmapped alignments, secondary alignments, duplicates and alignments with low mapping quality.
  *          Then, the CIGAR string of all remaining alignments is analyzed.
  *          For primary alignments, also the split read information is analyzed.
  */
 void detect_junctions_in_short_reads_sam_file(std::vector<Junction> & junctions,
-                                              std::filesystem::path const & alignment_short_reads_file_path,
-                                              std::vector<detection_methods> const & methods,
-                                              uint64_t const min_var_length);
+                                              cmd_arguments const & args);
 
 /*! \brief Detects junctions between distant genomic positions by analyzing a long read alignment file (sam/bam). The
  *         detected junctions are stored in a vector.
  *
  * \param[in, out]  junctions - a vector of junctions
- * \param[in]       alignment_long_reads_file_path - long reads input file, path to the sam/bam file
- * \param[in]       methods - list of methods for detecting junctions (0: cigar_string,
- *                                                                     1: split_read,
- *                                                                     2: read_pairs,
- *                                                                     3: read_depth)
- * \param[in]       min_var_length - minimum length of variants to detect (default 30 bp)
+ * \param[in]       args - command line arguments:\n
+ *                         **args.alignment_long_reads_file_path** - long reads input file, path to the sam/bam file\n
+ *                         **args.methods** - list of methods for detecting junctions
+ *                            (0: cigar_string, 1: split_read, 2: read_pairs, 3: read_depth) - *default: all methods*\n
+ *                         **args.min_var_length** - minimum length of variants to detect (expected to be non-negative) - *default: 30 bp*\n
+ *                         **args.max_overlap** - maximum overlap between alignment segments (expected to be non-negative) - *default: 10 bp*
+ *
  *
  * \details Detects junctions from the CIGAR strings and supplementary alignment tags of read alignment records.
  *          We filter unmapped alignments, secondary alignments, duplicates and alignments with low mapping quality.
  *          Then, the CIGAR string of all remaining alignments is analyzed.
  *          For primary alignments, also the split read information is analyzed.
  */
 void detect_junctions_in_long_reads_sam_file(std::vector<Junction> & junctions,
-                                             std::filesystem::path const & alignment_long_reads_file_path,
-                                             std::vector<detection_methods> const & methods,
-                                             uint64_t const min_var_length);
+                                             cmd_arguments const & args);

include/variant_detection/variant_output.hpp

@@ -12,9 +12,9 @@
  *
  * \param[in]       clusters    - input junction clusters
  * \param[in]       args        - command line arguments:\n
- *                                **args.min_var_length** - minimum length of variants to detect - *default: 30 bp*\n
- *                                **args.max_var_length** - maximum length of variants to detect - *default: 1,000,000 bp*\n
- *                                **args.max_tol_inserted_length** - longest tolerated inserted sequence at non-INS SV types - *default: 5 bp*
+ *                                **args.min_var_length** - minimum length of variants to detect (expected to be non-negative) - *default: 30 bp* \n
+ *                                **args.max_var_length** - maximum length of variants to detect (expected to be non-negative) - *default: 1,000,000 bp*\n
+ *                                **args.max_tol_inserted_length** - longest tolerated inserted sequence at non-INS SV types (expected to be non-negative) - *default: 5 bp*
  * \param[in, out]  out_stream  - output stream
  *
  * \details Extracts genomic variants from given junction clusters.
@@ -29,9 +29,9 @@ void find_and_output_variants(std::vector<Cluster> const & clusters,
  *
  * \param[in] clusters          - input junction clusters
  * \param[in] args              - command line arguments:\n
- *                                **args.min_var_length** - minimum length of variants to detect - *default: 30 bp*\n
- *                                **args.max_var_length** - maximum length of variants to detect - *default: 1,000,000 bp*\n
- *                                **args.max_tol_inserted_length** - longest tolerated inserted sequence at non-INS SV types - *default: 5 bp*
+ *                                **args.min_var_length** - minimum length of variants to detect (expected to be non-negative) - *default: 30 bp*\n
+ *                                **args.max_var_length** - maximum length of variants to detect (expected to be non-negative) - *default: 1,000,000 bp*\n
+ *                                **args.max_tol_inserted_length** - longest tolerated inserted sequence at non-INS SV types (expected to be non-negative) - *default: 5 bp*
  * \param[in] output_file_path  - output file path
  *
  * \details Extracts genomic variants from given junction clusters.

src/iGenVar.cpp

@@ -56,14 +56,21 @@ void initialize_argument_parser(seqan3::argument_parser & parser, cmd_arguments
 
     // Options - SV specifications:
     parser.add_option(args.min_var_length, 'l', "min_var_length",
-                      "Specify what should be the minimum length of your SVs to be detected.",
+                      "Specify what should be the minimum length of your SVs to be detected. "
+                      "This value needs to be non-negative.",
                       seqan3::option_spec::advanced);
     parser.add_option(args.max_var_length, 'x', "max_var_length",
                       "Specify what should be the maximum length of your SVs to be detected. "
-                      "SVs larger than this threshold can still be output as translocations.",
+                      "SVs larger than this threshold can still be output as translocations. "
+                      "This value needs to be non-negative.",
                       seqan3::option_spec::advanced);
     parser.add_option(args.max_tol_inserted_length, 't', "max_tol_inserted_length",
-                      "Specify what should be the longest tolerated inserted sequence at sites of non-INS SVs.",
+                      "Specify what should be the longest tolerated inserted sequence at sites of non-INS SVs. "
+                      "This value needs to be non-negative.",
+                      seqan3::option_spec::advanced);
+    parser.add_option(args.max_overlap, 'p', "max_overlap",
+                      "Specify the maximum allowed overlap between two alignment segments. "
+                      "This value needs to be non-negative.",
                       seqan3::option_spec::advanced);
 }
 
@@ -77,19 +84,15 @@ void detect_variants_in_alignment_file(cmd_arguments const & args)
     {
         seqan3::debug_stream << "Detect junctions in short reads...\n";
         detect_junctions_in_short_reads_sam_file(junctions,
-                                                 args.alignment_short_reads_file_path,
-                                                 args.methods,
-                                                 args.min_var_length);
+                                                 args);
     }
 
     // long reads
     if (args.alignment_long_reads_file_path != "")
     {
         seqan3::debug_stream << "Detect junctions in long reads...\n";
         detect_junctions_in_long_reads_sam_file(junctions,
-                                                args.alignment_long_reads_file_path,
-                                                args.methods,
-                                                args.min_var_length);
+                                                args);
     }
 
     std::sort(junctions.begin(), junctions.end());
@@ -168,6 +171,28 @@ int main(int argc, char ** argv)
         return -1;
     }
 
+    // Check that the given parameters are non-negative.
+    if (args.min_var_length < 0)
+    {
+        seqan3::debug_stream << "[Error] You gave a negative min_var_length parameter.\n";
+        return -1;
+    }
+    if (args.max_var_length < 0)
+    {
+        seqan3::debug_stream << "[Error] You gave a negative max_var_length parameter.\n";
+        return -1;
+    }
+    if (args.max_tol_inserted_length < 0)
+    {
+        seqan3::debug_stream << "[Error] You gave a negative max_tol_inserted_length parameter.\n";
+        return -1;
+    }
+    if (args.max_overlap < 0)
+    {
+        seqan3::debug_stream << "[Error] You gave a negative max_overlap parameter.\n";
+        return -1;
+    }
+
     detect_variants_in_alignment_file(args);
 
     return 0;

src/modules/sv_detection_methods/analyze_cigar_method.cpp

@@ -15,7 +15,7 @@ void analyze_cigar(std::string const & read_name,
                    std::vector<seqan3::cigar> & cigar_string,
                    seqan3::dna5_vector const & query_sequence,
                    std::vector<Junction> & junctions,
-                   uint64_t const min_length)
+                   int32_t const min_length)
 {
     // Step through CIGAR string and store current position in reference and read
     int32_t pos_ref = query_start_pos;