Merge pull request #145 from Irallia/MISC/small_corrections_2

[MISC, DOC] Style changes and some documentation

Author: Irallia

Diff for: LICENSE.md

@@ -1,4 +1,4 @@
-BSD 3-Clause License
+# BSD 3-Clause License
 
 Copyright (c) 2020, Jörg Winkler
 All rights reserved.

Diff for: README.md

@@ -4,11 +4,29 @@
 
 The official repository for the iGenVar project.
 
+iGenVar is intended to be a caller for all types of genetic variation: SNPs, indels and larger structural variations
+(insertions, deletions, inversions, translocations, CNVs, nested SVs).
+It uses both Illumina short reads and PacBio long reads for this purpose.
+
+David Heller in the Vingron lab of the MPI-MG and Tim White in the Kehr lab at BIH have both developed an SV caller for
+long read sequencing data. Instead of competing with each other, we want to join forces and combine the two tools, SVIM
+and SVIRL, into one better and more versatile tool.
+On the other hand, there were some tool developments in the Reinert Lab (FU Berlin): Vaquita, a short read SV caller,
+SViper, a refinement tool and Vaquita-LR a further development of Vaquita for long reads.
+We want to combine these approaches and use the SeqAn3 library as a basis for this new tool.
+
+## Current status:
+
+We can call insertions and deletions from long read data (SVIM methods implemented).
+For more information, see the release plan at the bottom of the page.
+
 ## Installation
 
 Instructions:
+
 1. clone this repository: `git clone --recurse-submodules https://github.com/seqan/iGenVar.git`
-    or `git clone https://github.com/seqan/iGenVar.git` and fetch the seqan3 submodule after cloning: `git submodule update --recursive --init`
+    or `git clone https://github.com/seqan/iGenVar.git` and fetch the seqan3 submodule after cloning:
+    `git submodule update --recursive --init`
 2. create a build directory and visit it: `mkdir build && cd build`
 3. run cmake: `cmake ../iGenVar`
 4. build the application: `make`
@@ -19,5 +37,11 @@ Instructions:
 (Built using the [SeqAn3 App Template](https://github.com/seqan/app-template))
 
 We created small examples, which you can use to test our app:
-`./bin/iGenVar -i ./test/data/paired_end_short_read_mini_example.sam -j ./test/data/single_end_mini_example.sam `
-`-o ./test/data/output.vcf --method cigar_string --method split_read --min_var_length 5`
+```bash
+./bin/iGenVar -i ./test/data/paired_end_short_read_mini_example.sam -j ./test/data/single_end_mini_example.sam \
+-o ./test/data/output.vcf --method cigar_string --method split_read --min_var_length 5
+```
+
+## Release plan:
+
+<p align="center"><img height="500" src="https://github.com/seqan/iGenVar/tree/master/doc/ReleasePlan.png"></p>

Diff for: doc/ReleasePlan.png

@@ -43,5 +43,4 @@ int junction_distance(Junction const & lhs, Junction const & rhs);
  * \details For the algorithms we use the library hclust.
  * \see https://lionel.kr.hs-niederrhein.de/~dalitz/data/hclust/ (last access 01.06.2021).
  */
-std::vector<Cluster> hierarchical_clustering_method(std::vector<Junction> const & junctions,
-                                                    double clustering_cutoff);
+std::vector<Cluster> hierarchical_clustering_method(std::vector<Junction> const & junctions, double clustering_cutoff);

Diff for: include/modules/clustering/hierarchical_clustering_method.hpp

@@ -5,7 +5,7 @@ add_library ("${PROJECT_NAME}_lib" STATIC modules/clustering/hierarchical_cluste
                                           modules/clustering/simple_clustering_method.cpp
                                           modules/sv_detection_methods/analyze_cigar_method.cpp
                                           modules/sv_detection_methods/analyze_read_pair_method.cpp
-                                          modules/sv_detection_methods/analyze_sa_tag_method.cpp
+                                          modules/sv_detection_methods/analyze_split_read_method.cpp
                                           structures/aligned_segment.cpp
                                           structures/breakend.cpp
                                           structures/cluster.cpp

Diff for: include/modules/sv_detection_methods/analyze_sa_tag_method.hpp renamed to include/modules/sv_detection_methods/analyze_split_read_method.hpp

@@ -116,8 +116,7 @@ inline std::vector<Junction> subsample_partition(std::vector<Junction> const & p
     return subsample;
 }
 
-std::vector<Cluster> hierarchical_clustering_method(std::vector<Junction> const & junctions,
-                                                    double clustering_cutoff)
+std::vector<Cluster> hierarchical_clustering_method(std::vector<Junction> const & junctions, double clustering_cutoff)
 {
     auto partitions = partition_junctions(junctions);
     std::vector<Cluster> clusters{};

Diff for: src/CMakeLists.txt

@@ -1,4 +1,4 @@
-#include "modules/sv_detection_methods/analyze_sa_tag_method.hpp"
+#include "modules/sv_detection_methods/analyze_split_read_method.hpp"
 
 #include <seqan3/core/debug_stream.hpp>
 
@@ -97,7 +97,7 @@ void analyze_aligned_segments(std::vector<AlignedSegment> const & aligned_segmen
             // map to different reference sequences (e.g. translocation, interspersed duplication),
             // have a large distance on the reference (e.g. deletion, inversion, tandem duplication), or
             // have a large distance on the read (e.g. insertion)
-            if (current.ref_name != next.ref_name ||
+            if (current.ref_name != next.ref_name || //TODO / QUESTION (irallia 23.06.2021): What about translocation on the same ref?
                 std::abs(distance_on_ref) >= min_length ||
                 distance_on_read >= min_length)
             {
@@ -114,7 +114,8 @@ void analyze_aligned_segments(std::vector<AlignedSegment> const & aligned_segmen
                 }
                 else
                 {
-                    auto inserted_bases = query_sequence | seqan3::views::slice(current.get_query_end(), next.get_query_start());
+                    auto inserted_bases = query_sequence | seqan3::views::slice(current.get_query_end(),
+                                                                                next.get_query_start());
                     junctions.emplace_back(mate1, mate2, inserted_bases, read_name);
                 }
                 seqan3::debug_stream << "BND: " << junctions.back() << "\n";
@@ -139,6 +140,7 @@ void analyze_sa_tag(std::string const & query_name,
     strand strand = (hasFlagReverseComplement(flag) ? strand::reverse : strand::forward);
     aligned_segments.push_back(AlignedSegment{strand, ref_name, pos, mapq, cigar});
     retrieve_aligned_segments(sa_tag, aligned_segments);
+    // sort by query start, query end, mapping quality (in this order):
     std::sort(aligned_segments.begin(), aligned_segments.end());
     analyze_aligned_segments(aligned_segments,
                              junctions,

Diff for: src/modules/clustering/hierarchical_clustering_method.cpp

@@ -2,6 +2,7 @@
 
 std::unordered_map<std::string, detection_methods> enumeration_names(detection_methods)
 {
+// LCOV_EXCL_START
 return std::unordered_map<std::string, detection_methods>{{"0", detection_methods::cigar_string},
                                                           {"cigar_string", detection_methods::cigar_string},
                                                           {"1", detection_methods::split_read},
@@ -10,27 +11,31 @@ return std::unordered_map<std::string, detection_methods>{{"0", detection_method
                                                           {"read_pairs", detection_methods::read_pairs},
                                                           {"3", detection_methods::read_depth},
                                                           {"read_depth", detection_methods::read_depth}};
+// LCOV_EXCL_STOP
 };
 
 std::unordered_map<std::string, clustering_methods> enumeration_names(clustering_methods)
 {
+// LCOV_EXCL_START
 return std::unordered_map<std::string,
-                            clustering_methods>{{"0", clustering_methods::simple_clustering},
-                                                {"simple_clustering",
-                                                clustering_methods::simple_clustering},
-                                                {"1", clustering_methods::hierarchical_clustering},
-                                                {"hierarchical_clustering",
-                                                clustering_methods::hierarchical_clustering},
-                                                {"2", clustering_methods::self_balancing_binary_tree},
-                                                {"self_balancing_binary_tree",
-                                                clustering_methods::self_balancing_binary_tree},
-                                                {"3", clustering_methods::candidate_selection_based_on_voting},
-                                                {"candidate_selection_based_on_voting",
-                                                clustering_methods::candidate_selection_based_on_voting}};
+                          clustering_methods>{{"0", clustering_methods::simple_clustering},
+                                              {"simple_clustering",
+                                              clustering_methods::simple_clustering},
+                                              {"1", clustering_methods::hierarchical_clustering},
+                                              {"hierarchical_clustering",
+                                              clustering_methods::hierarchical_clustering},
+                                              {"2", clustering_methods::self_balancing_binary_tree},
+                                              {"self_balancing_binary_tree",
+                                              clustering_methods::self_balancing_binary_tree},
+                                              {"3", clustering_methods::candidate_selection_based_on_voting},
+                                              {"candidate_selection_based_on_voting",
+                                              clustering_methods::candidate_selection_based_on_voting}};
+// LCOV_EXCL_STOP
 };
 
 std::unordered_map<std::string, refinement_methods> enumeration_names(refinement_methods)
 {
+// LCOV_EXCL_START
 return std::unordered_map<std::string, refinement_methods>{{"0", refinement_methods::no_refinement},
                                                            {"no_refinement",
                                                            refinement_methods::no_refinement},
@@ -40,4 +45,5 @@ return std::unordered_map<std::string, refinement_methods>{{"0", refinement_meth
                                                            {"2", refinement_methods::sVirl_refinement_method},
                                                            {"sVirl_refinement_method",
                                                            refinement_methods::sVirl_refinement_method}};
+// LCOV_EXCL_STOP
 };

Diff for: src/modules/sv_detection_methods/analyze_sa_tag_method.cpp renamed to src/modules/sv_detection_methods/analyze_split_read_method.cpp

@@ -3,10 +3,10 @@
 #include <seqan3/core/debug_stream.hpp>
 #include <seqan3/io/sam_file/input.hpp>         // SAM/BAM support (seqan3::sam_file_input)
 
-#include "modules/sv_detection_methods/analyze_cigar_method.hpp"    // for the split read method
-#include "modules/sv_detection_methods/analyze_read_pair_method.hpp"// for the read pair method
-#include "modules/sv_detection_methods/analyze_sa_tag_method.hpp"   // for the cigar string method
-#include "variant_detection/bam_functions.hpp"                      // for hasFlag* functions
+#include "modules/sv_detection_methods/analyze_cigar_method.hpp"        // for the split read method
+#include "modules/sv_detection_methods/analyze_read_pair_method.hpp"    // for the read pair method
+#include "modules/sv_detection_methods/analyze_split_read_method.hpp"   // for the cigar string method
+#include "variant_detection/bam_functions.hpp"                          // for hasFlag* functions
 
 using seqan3::operator""_tag;
 
@@ -167,8 +167,7 @@ void detect_junctions_in_long_reads_sam_file(std::vector<Junction> & junctions,
                         }
                     }
                     break;
-                case detection_methods::read_pairs: // There are no read pairs in long reads.
-                    break;
+                // There are no read pairs in long reads.
                 case detection_methods::read_depth: // Detect junctions from read depth evidence
                     seqan3::debug_stream << "The read depth method for long reads is not yet implemented.\n";
                     break;

Diff for: src/variant_detection/method_enums.cpp

@@ -8,3 +8,5 @@ add_api_test (detection_test.cpp)
 add_api_test (clustering_test.cpp)
 
 # add_api_test (refinement_test.cpp)
+
+add_api_test (structures_test.cpp)

Diff for: src/variant_detection/variant_detection.cpp

@@ -3,8 +3,8 @@
 #include <seqan3/alphabet/cigar/cigar.hpp>
 #include <seqan3/io/sam_file/sam_flag.hpp>
 
-#include "modules/sv_detection_methods/analyze_cigar_method.hpp"    // for the split read method
-#include "modules/sv_detection_methods/analyze_sa_tag_method.hpp"   // for the cigar string method
+#include "modules/sv_detection_methods/analyze_cigar_method.hpp"        // for the split read method
+#include "modules/sv_detection_methods/analyze_split_read_method.hpp"   // for the cigar string method
 
 using seqan3::operator""_cigar_operation;
 using seqan3::operator""_dna5;

Diff for: test/api/CMakeLists.txt

@@ -0,0 +1,90 @@
+#include <gtest/gtest.h>
+
+#include "structures/breakend.hpp"
+#include "variant_detection/method_enums.hpp"
+
+/* tests for method_enums */
+
+TEST(structures, method_enums_detection_methods)
+{
+    std::unordered_map<std::string,
+                       detection_methods> cigar_string_mapping = enumeration_names(detection_methods::cigar_string);
+    EXPECT_EQ(detection_methods::cigar_string, cigar_string_mapping["cigar_string"]);
+    EXPECT_EQ(detection_methods::cigar_string, cigar_string_mapping["0"]);
+
+    std::unordered_map<std::string,
+                       detection_methods> split_read_mapping = enumeration_names(detection_methods::split_read);
+    EXPECT_EQ(detection_methods::split_read, split_read_mapping["split_read"]);
+    EXPECT_EQ(detection_methods::split_read, split_read_mapping["1"]);
+
+    std::unordered_map<std::string,
+                       detection_methods> read_pairs_mapping = enumeration_names(detection_methods::read_pairs);
+    EXPECT_EQ(detection_methods::read_pairs, read_pairs_mapping["read_pairs"]);
+    EXPECT_EQ(detection_methods::read_pairs, read_pairs_mapping["2"]);
+
+    std::unordered_map<std::string,
+                       detection_methods> read_depth_mapping = enumeration_names(detection_methods::read_depth);
+    EXPECT_EQ(detection_methods::read_depth, read_depth_mapping["read_depth"]);
+    EXPECT_EQ(detection_methods::read_depth, read_depth_mapping["3"]);
+}
+
+TEST(structures, method_enums_clustering_methods)
+{
+    std::unordered_map<std::string, clustering_methods> simple_clustering_mapping
+                                        = enumeration_names(clustering_methods::simple_clustering);
+    EXPECT_EQ(clustering_methods::simple_clustering, simple_clustering_mapping["simple_clustering"]);
+    EXPECT_EQ(clustering_methods::simple_clustering, simple_clustering_mapping["0"]);
+
+    std::unordered_map<std::string, clustering_methods> hierarchical_clustering_mapping
+                                        = enumeration_names(clustering_methods::hierarchical_clustering);
+    EXPECT_EQ(clustering_methods::hierarchical_clustering, hierarchical_clustering_mapping["hierarchical_clustering"]);
+    EXPECT_EQ(clustering_methods::hierarchical_clustering, hierarchical_clustering_mapping["1"]);
+
+    std::unordered_map<std::string, clustering_methods> self_balancing_binary_tree_mapping
+                                        = enumeration_names(clustering_methods::self_balancing_binary_tree);
+    EXPECT_EQ(clustering_methods::self_balancing_binary_tree,
+              self_balancing_binary_tree_mapping["self_balancing_binary_tree"]);
+    EXPECT_EQ(clustering_methods::self_balancing_binary_tree, self_balancing_binary_tree_mapping["2"]);
+
+    std::unordered_map<std::string, clustering_methods> candidate_selection_based_on_voting_mapping
+                                        = enumeration_names(clustering_methods::candidate_selection_based_on_voting);
+    EXPECT_EQ(clustering_methods::candidate_selection_based_on_voting,
+              candidate_selection_based_on_voting_mapping["candidate_selection_based_on_voting"]);
+    EXPECT_EQ(clustering_methods::candidate_selection_based_on_voting,
+              candidate_selection_based_on_voting_mapping["3"]);
+}
+
+TEST(structures, method_enums_refinement_methods)
+{
+    std::unordered_map<std::string, refinement_methods> no_refinement_mapping
+                                        = enumeration_names(refinement_methods::no_refinement);
+    EXPECT_EQ(refinement_methods::no_refinement, no_refinement_mapping["no_refinement"]);
+    EXPECT_EQ(refinement_methods::no_refinement, no_refinement_mapping["0"]);
+
+    std::unordered_map<std::string, refinement_methods> sViper_refinement_method_mapping
+                                        = enumeration_names(refinement_methods::sViper_refinement_method);
+    EXPECT_EQ(refinement_methods::sViper_refinement_method,
+              sViper_refinement_method_mapping["sViper_refinement_method"]);
+    EXPECT_EQ(refinement_methods::sViper_refinement_method, sViper_refinement_method_mapping["1"]);
+
+    std::unordered_map<std::string, refinement_methods> sVirl_refinement_method_mapping
+                                        = enumeration_names(refinement_methods::sVirl_refinement_method);
+    EXPECT_EQ(refinement_methods::sVirl_refinement_method, sVirl_refinement_method_mapping["sVirl_refinement_method"]);
+    EXPECT_EQ(refinement_methods::sVirl_refinement_method, sVirl_refinement_method_mapping["2"]);
+}
+
+/* tests for junctions */
+
+TEST(structures, breakend_flip_orientation)
+{
+    Breakend forward_breakend{"chr1", 42, strand::forward};
+    Breakend reverse_breakend{"chr1", 42, strand::reverse};
+
+    EXPECT_NE(forward_breakend, reverse_breakend);
+    forward_breakend.flip_orientation();
+    EXPECT_EQ(forward_breakend, reverse_breakend); // both are reverse now
+    reverse_breakend.flip_orientation();
+    EXPECT_NE(forward_breakend, reverse_breakend);
+    forward_breakend.flip_orientation();
+    EXPECT_EQ(forward_breakend, reverse_breakend); // both are forward now
+}