Merge pull request #171 from joergi-w/snp_indel_call

Calculate active regions based on Cigar information

Author: Irallia

include/variant_detection/snp_indel_detection.hpp

@@ -2,8 +2,30 @@
 
 #include <seqan3/std/filesystem>
 
+#include <seqan3/alphabet/cigar/cigar.hpp>
+
 /*!
  * \brief Detect Single Nucleotide Polymorphisms (SNPs) and short deletions and insertions.
  * \param[in] reads_filename - The file path where to find the sequenced reads.
+ * \param[in] min_var_length - The length above which an indel/SNP is considered a variant.
+ */
+void detect_snp_and_indel(std::filesystem::path const & reads_filename, uint64_t min_var_length);
+
+/*!
+ * \brief Extract activity from SAM records by counting indels and soft clips.
+ * \param[in,out] activity       - The activity values for one reference genome.
+ * \param[in]     cigar_sequence - The cigar string of the SAM record.
+ * \param[in]     min_var_length - The length above which an indel/SNP is considered a variant.
+ * \param[in]     ref_pos        - The start position of the alignment in the genome.
+ */
+void update_activity_for_record(std::vector<unsigned> & activity,
+                                std::vector<seqan3::cigar> const & cigar_sequence,
+                                uint64_t min_var_length,
+                                int32_t ref_pos);
+
+/*!
+ * \brief Extract active regions from activity profile.
+ * \param[in] activity - The activity values for one reference genome.
+ * \return a list of position intervals, where the activity is high.
  */
-void detect_snp_and_indel(std::filesystem::path const & reads_filename);
+std::vector<std::pair<size_t, size_t>> active_regions(std::vector<unsigned> const & activity);

src/iGenVar.cpp

@@ -150,7 +150,7 @@ void detect_variants_in_alignment_file(cmd_arguments const & args)
     if (!args.alignment_short_reads_file_path.empty() && !args.genome_file_path.empty())
     {
         seqan3::debug_stream << "Detect SNPs and indels in short reads...\n";
-        detect_snp_and_indel(args.alignment_short_reads_file_path);
+        detect_snp_and_indel(args.alignment_short_reads_file_path, args.min_var_length);
     }
 
     std::sort(junctions.begin(), junctions.end());

src/variant_detection/snp_indel_detection.cpp

@@ -1,14 +1,106 @@
 #include <iostream>
 
-#include <seqan3/alphabet/cigar/cigar.hpp>
 #include <seqan3/core/debug_stream.hpp>
 #include <seqan3/io/sam_file/input.hpp>
 
 #include "iGenVar.hpp"
 #include "variant_detection/bam_functions.hpp"          // check the sam flags
 #include "variant_detection/snp_indel_detection.hpp"    // detect_snp_and_indel
 
-void detect_snp_and_indel(std::filesystem::path const & reads_filename)
+void update_activity_for_record(std::vector<unsigned> & activity,
+                                std::vector<seqan3::cigar> const & cigar_sequence,
+                                uint64_t const min_var_length,
+                                int32_t ref_pos)
+{
+    // Track the current position within the read.
+    unsigned read_pos = 0;
+
+    // Step through the CIGAR string.
+    for (auto && [length, operation] : cigar_sequence)
+    {
+        // Skip long variants.
+        if (length >= min_var_length)
+        {
+            read_pos += length;
+            ref_pos += length;
+            continue;
+        }
+
+        // Case distinction for cigar elements.
+        switch (operation.to_char())
+        {
+            case 'I':
+            case 'S': // insertion or soft clip
+            {
+                activity[ref_pos] += 1u;
+                read_pos += length;
+                break;
+            }
+            case 'D': // deletion
+            {
+                for (unsigned idx = 0; idx < length; ++idx)
+                    ++activity[ref_pos + idx];
+                ref_pos += length;
+                break;
+            }
+            default: // ignore match, mismatch, hard clipping, padding, and skipped region
+            { // TODO (joergi-w 30.09.2021) mismatches are ignored but should be verified with the reference
+                ref_pos += length;
+                read_pos += length;
+            }
+        }
+    }
+}
+
+std::vector<std::pair<size_t, size_t>> active_regions(std::vector<unsigned> const & activity)
+{
+    // Skip reference sequences that are not contained in the reads.
+    if (activity.empty())
+        return {};
+
+    // Some constants for the sliding window approach.
+    size_t const window_width = 5u; // another idea: (min_var_length + 1) / 2;
+    size_t const activity_threshold = 2; // TODO (joergi-w 25.10.2021) depend on read coverage
+
+    // Set up the result vector.
+    std::vector<std::pair<size_t, size_t>> regions{};
+
+    size_t window_score = 0; // sum of the activities inside the window
+    size_t pos;              // the current position in the genome sequence
+    bool active = false;     // whether `pos` is inside an active region
+
+    // Initialisation of a window of length min_var_length.
+    for (pos = 0; pos < window_width && pos < activity.size(); ++pos)
+        window_score += activity[pos];
+
+    while (pos < activity.size())
+    {
+        // Check for start or end of an active region.
+        if (!active && window_score >= activity_threshold)
+        {
+            regions.emplace_back(pos - window_width, 0);
+            active = true;
+        }
+        else if (active && window_score < activity_threshold)
+        {
+            std::get<1>(regions.back()) = pos - 1u;
+            active = false;
+        }
+
+        // Advance the sliding window.
+        window_score -= activity[pos - window_width];
+        window_score += activity[pos];
+        ++pos;
+    }
+
+    // Handle the case if the genome ends in an active region.
+    if (active)
+        std::get<1>(regions.back()) = pos - 1u;
+
+    return regions;
+}
+
+void detect_snp_and_indel(std::filesystem::path const & reads_filename, uint64_t min_var_length)
 {
     // Get the header information and set the necessary fields.
     using sam_fields = seqan3::fields<seqan3::field::flag,       // 2: FLAG
@@ -18,17 +110,14 @@ void detect_snp_and_indel(std::filesystem::path const & reads_filename)
                                       seqan3::field::cigar>;     // 6: CIGAR
 
     seqan3::sam_file_input reads_file{reads_filename, sam_fields{}};
-    auto const & header = reads_file.header();
+    auto & header = reads_file.header();
 
     // Check that the file is sorted.
     if (header.sorting != "coordinate")
         throw seqan3::format_error{"ERROR: Input file must be sorted by coordinate (e.g. samtools sort)"};
 
     // Prepare one activity vector per reference sequence.
-    std::vector<std::vector<unsigned>> activity{};
-    activity.resize(header.ref_id_info.size()); // allocate the number of reference sequences
-    for (size_t idx = 0; idx < activity.size(); ++idx)
-        activity[idx].resize(std::get<0>(header.ref_id_info[idx])); // allocate the length of the reference sequence
+    std::vector<std::vector<unsigned>> activity_vec(header.ref_id_info.size()); // number of ref sequences
 
     for (auto && record : reads_file)
     {
@@ -47,48 +136,18 @@ void detect_snp_and_indel(std::filesystem::path const & reads_filename)
             continue;
         }
 
-        // Track the current position within the read.
-        unsigned read_pos = 0;
-
-        // Step through the CIGAR string.
-        for (auto && [length, operation] : record.cigar_sequence())
-        {
-            // Skip long variants.
-            if (length >= 30) // TODO (joergi-w 30.09.2021) Use the min_length parameter here.
-            {
-                read_pos += length;
-                ref_pos += length;
-                continue;
-            }
+        // Check if we have already allocated the length of the reference sequence.
+        if (activity_vec[ref_id].empty())
+            activity_vec[ref_id].resize(std::get<0>(header.ref_id_info[ref_id])); // length of reference sequence
 
-            // Case distinction for cigar elements.
-            seqan3::cigar::operation const op = operation;
-            switch (op.to_char())
-            {
-                case 'I':
-                case 'S': // insertion or soft clip
-                {
-                    activity[ref_id][ref_pos] += length;
-                    read_pos += length;
-                    break;
-                }
-                case 'D': // deletion
-                {
-                    for (unsigned idx = 0; idx < length; ++idx)
-                        ++activity[ref_id][ref_pos + idx];
-                    ref_pos += length;
-                    break;
-                }
-                default: // ignore match, mismatch, hard clipping, padding, and skipped region
-                { // TODO (joergi-w 30.09.2021) mismatches are ignored but should be verified with the reference
-                    ref_pos += length;
-                    read_pos += length;
-                }
-            }
-        }
+        update_activity_for_record(activity_vec[ref_id], record.cigar_sequence(), min_var_length, ref_pos);
     }
 
-    // Print the raw activity profiles.
-    for (auto const & av : activity)
-        seqan3::debug_stream << av << std::endl;
+    // Extract active regions from activity profile.
+    for (size_t ref_id = 0; ref_id < activity_vec.size(); ++ref_id)
+    {
+        auto regions = active_regions(activity_vec[ref_id]);
+        if (!regions.empty())
+            seqan3::debug_stream << "Active regions of " << header.ref_ids()[ref_id] << ": " << regions << std::endl;
+    }
 }

test/api/CMakeLists.txt

@@ -9,4 +9,5 @@ add_api_test (clustering_test.cpp)
 
 # add_api_test (refinement_test.cpp)
 
+add_api_test (snp_indel_test.cpp)
 add_api_test (structures_test.cpp)

test/api/snp_indel_test.cpp

@@ -0,0 +1,40 @@
+#include <gtest/gtest.h>
+
+#include <vector>
+
+#include "variant_detection/snp_indel_detection.hpp"
+
+TEST(activity_analysis, input_empty)
+{
+    std::vector<unsigned> activity{};
+    auto regions = active_regions(activity);
+    EXPECT_TRUE(regions.empty());
+}
+
+TEST(activity_analysis, input_zero)
+{
+    std::vector<unsigned> activity(150, 0u);
+    auto regions = active_regions(activity);
+    EXPECT_TRUE(regions.empty());
+}
+
+TEST(activity_analysis, input_active_everywhere)
+{
+    std::vector<unsigned> activity(150, 9999u);
+    auto regions = active_regions(activity);
+    EXPECT_EQ(regions.size(), 1u);
+    EXPECT_EQ(regions[0], (std::pair<size_t, size_t>{0, 149}));
+}
+
+TEST(activity_analysis, input_interval)
+{
+    std::vector<unsigned> activity(150, 0u);
+    for (size_t pos = 34u; pos <= 44u; ++pos)
+        activity[pos] = 6;
+    for (size_t pos = 134u; pos <= 144u; ++pos)
+        activity[pos] = 1;
+    auto regions = active_regions(activity);
+    EXPECT_EQ(regions.size(), 2u);
+    EXPECT_EQ(regions[0], (std::pair<size_t, size_t>{30, 49}));
+    EXPECT_EQ(regions[1], (std::pair<size_t, size_t>{131, 148}));
+}

test/cli/iGenVar_cli_test.cpp

@@ -302,19 +302,8 @@ TEST_F(iGenVar_cli_test, test_genome_input)
                                          "-g", data(default_genome_file_path),
                                          "-i", data("single_end_mini_example.sam"));
     std::string const expected_err = "Detect SNPs and indels in short reads...\n"
-                                     "[0,0,0,0,0,0,0,0,0,0,48,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,0,0,0,0,"
-                                     "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,10,9,9,9,9,9,9,9,9,9,9,9,9,1,0,0,0,0,0,0,0,"
-                                     "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
-                                     "0,0,0,0,0,0,0,83,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
-                                     "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,33,0,0,0,0,0,0,0,15,0,0,0,0,0,0,0,"
-                                     "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
-                                     "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,13,2,2,2,2,2,2,2,2,"
-                                     "2,2,2,2,2,2,2,67,3,3,3,67,1,1,1,1,1,1,1,1,1,1,1,1,43,0,0,0,0,0,0,0,0,0,0,0,0,0,"
-                                     "0,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,4,4,4,4,4,4,4,4,4,4,4,4,4,0,0,0,"
-                                     "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,55,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,"
-                                     "1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
-                                     "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,52,"
-                                     "0,0,0,0,0,0,0,0,0,0]\n"
+                                     "Active regions of chr1: [(6,15),(53,74),(121,130),(176,185),(184,193),(262,304),"
+                                     "(311,319),(332,354),(364,373),(381,398),(467,476)]\n"
                                      "Start clustering...\n"
                                      "Done with clustering. Found 0 junction clusters.\n"
                                      "No refinement was selected.\n";