fixing boost anycast error

CMakeLists.txt

@@ -1,5 +1,5 @@
 cmake_minimum_required(VERSION 3.22.1)
-project(RNAnue VERSION 0.2.0)
+project(RNAnue VERSION 0.2.1)
 set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_CXX_STANDARD_REQUIRED True)
 set(CMAKE_CXX_FLAGS -fopenmp)

Dockerfile

@@ -15,7 +15,7 @@ LABEL authors="Richard A. Schaefer"
 RUN apt-get -y update && apt-get -y upgrade
 RUN apt-get install -y curl build-essential cmake git pkg-config
 RUN apt-get install -y libbz2-dev zlib1g-dev libncurses5-dev liblzma-dev
-RUN apt-get install -y libboost-all-dev gdb
+RUN apt-get install -y libboost-all-dev gdb libmpfr-dev liblapacke-dev
 
 # install htslib
 WORKDIR /
@@ -25,24 +25,30 @@ WORKDIR /htslib-1.20
 RUN ./configure
 RUN make
 RUN make install
+RUN export PKG_CONFIG_PATH=/usr/local/lib/pkgconfig
 
 # install segemehl
 WORKDIR /
 RUN curl -L http://legacy.bioinf.uni-leipzig.de/Software/segemehl/downloads/segemehl-0.3.4.tar.gz -o segemehl-0.3.4.tar.gz
 RUN tar -xvf segemehl-0.3.4.tar.gz && rm segemehl-0.3.4.tar.gz
 WORKDIR /segemehl-0.3.4
-RUN export PKG_CONFIG_PATH=/usr/local/lib/pkgconfig
 RUN make all
 RUN cp segemehl.x /usr/local/bin
 
 # install ViennaRNA
+WORKDIR /
+RUN curl -L https://www.tbi.univie.ac.at/RNA/download/sourcecode/2_6_x/ViennaRNA-2.6.4.tar.gz -o ViennaRNA-2.6.4.tar.gz
+RUN tar -xvf ViennaRNA-2.6.4.tar.gz && rm ViennaRNA-2.6.4.tar.gz
+WORKDIR /ViennaRNA-2.6.4
+RUN ./configure
+RUN make && make install
 
 # retrieve RNAnue
 WORKDIR /
 RUN git clone https://github.com/Ibvt/RNAnue.git
 WORKDIR /RNAnue
 
-# retrieve SeqAn
+# retrieve SeqAn3
 RUN curl -L https://github.com/seqan/seqan3/releases/download/3.3.0/seqan3-3.3.0-Source.tar.xz -o seqan3-3.3.0-Source.tar.xz
 RUN tar -xvf seqan3-3.3.0-Source.tar.xz && rm seqan3-3.3.0-Source.tar.xz
 RUN mv seqan3-3.3.0-Source seqan3

README.md

@@ -16,7 +16,7 @@ located in $PATH.
 * [Boost C++ Libraries](https://www.boost.org/) (v1.7.2)
 * [SeqAn](https://github.com/seqan/seqan3) (v3.3.0)
 * [Segemehl](http://www.bioinf.uni-leipzig.de/Software/segemehl/) (v0.3.4)
-* [Vienna Package](https://www.tbi.univie.ac.at/RNA/#binary_packages) (v2.4.17)
+* [Vienna Package](https://www.tbi.univie.ac.at/RNA/#binary_packages) (>= v2.4.17)
 
 ### CMake
 RNAnue is build using CMake. At first, clone the repository and change into the source directory.

include/Analysis.hpp

@@ -5,33 +5,73 @@
 #include <iostream>
 #include <map>
 #include <vector>
+#include <algorithm>
+#include <random>
+#include <cmath>
 
 // Boost
 #include <boost/program_options.hpp>
 #include <boost/filesystem.hpp>
 #include <boost/property_tree/ptree.hpp>
 #include <boost/accumulators/accumulators.hpp>
 #include <boost/accumulators/statistics.hpp>
+#include <boost/math/distributions/binomial.hpp>
+
+
+// SeqAn3
+#include <seqan3/io/sam_file/all.hpp>
+#include <seqan3/core/debug_stream.hpp>
 
 // Class
 #include "IBPTree.hpp"
 
+// define tags
+using seqan3::operator""_tag;
+
 namespace po = boost::program_options;
 namespace fs = boost::filesystem;
 namespace pt = boost::property_tree;
+namespace ma = boost::math;
 
 class Analysis {
     public:
         Analysis(po::variables_map params);
         ~Analysis();
 
-        void start(pt::ptree sample);
-        void process(std::string& single, std::string& splits);
+        void start(pt::ptree sample, pt::ptree condition);
+        void processSingle(std::string& single);
+        void processSplits(std::string& splits, std::string& output);
+        void normalize(); // normalize the frequencies to 1
+
+        // write output files (of the analysis)
+        void writeInteractionsHeader(std::ofstream& fout);
+        void writeAllIntsHeader(std::ofstream& fout);
+        void addToAllIntsHeader(std::ofstream& fout, std::string key);
+        void writeAllInts();
+
+        // other operations
+        void addToFreqMap(std::pair<std::string,std::string> key, double value);
+        void addToSuppReads(dtp::IntPartner p1, dtp::IntPartner p2);
+        void addToCompl(dtp::IntPartner p1, dtp::IntPartner p2, double complementarity);
+        void addToHybEnergies(dtp::IntPartner p1, dtp::IntPartner p2, double hybenergy);
+
+        // calculate filters
+        double calcGCS(std::vector<double>& complementarities);
+        double calcGHS(std::vector<double>& hybenergies);
+        double calcStat(dtp::IntKey key, int x);
 
     private:
         po::variables_map params;
         IBPTree features;
-        std::map<std::pair<std::string, std::string>, double> pdf;
+        std::map<std::pair<std::string,std::string>,double> freq; // strand, name
+        std::string condition; // buffers the current condition
+        // maps for storing filters and suppreads
+        std::map<dtp::IntKey, std::vector<double>> suppreads;
+        std::map<dtp::IntKey, std::vector<std::vector<double>>> complementarities;
+        std::map<dtp::IntKey, std::vector<std::vector<double>>> hybenergies;
+
+        int repcount; // counter for current replicate
+        int readcount; // total number of reads
 };
 
 #endif //RNANUE_ANALYSIS_HPP

include/Clustering.hpp

@@ -24,13 +24,13 @@ namespace fs = boost::filesystem;
 
 struct Segment {
     std::string refid;
-    uint32_t flag;
+    seqan3::sam_flag flag;
     uint32_t start;
     uint32_t end;
 
-    Segment() : refid(""), flag(0), start(0), end(0) {};
-    Segment(std::string _refid, uint32_t _flag, unsigned int _start, unsigned int _end) :
-            refid(_refid), flag(_flag), start(_start), end(_end) {};
+    Segment() : refid(""), flag(seqan3::sam_flag{}), start(0), end(0) {};
+    Segment(std::string refid, seqan3::sam_flag flag, unsigned int start, unsigned int end) :
+            refid(refid), flag(flag), start(start), end(end) {};
 };
 
 struct Cluster {

include/Config.hpp

@@ -1,4 +1,4 @@
 // the configured options and settings for Tutorial
 #define RNAnue_VERSION_MAJOR 0
 #define RNAnue_VERSION_MINOR 2
-#define RNAnue_VERSION_PATCH 0
+#define RNAnue_VERSION_PATCH 1

include/DataTypes.hpp

@@ -79,6 +79,26 @@ namespace dtp {
     };
     using StatsMap = std::map<std::string, StatsFields>;
     using SpliceJunctions = std::map<std::string, std::vector<std::pair<size_t,size_t>>>;
+
+    // Analysis
+    using PdfMap = std::map<std::pair<std::string, std::string>, double>;
+    struct IntPartner {
+        std::string partner;
+        std::string orientation;
+        IntPartner() : partner(""), orientation("") {}
+        IntPartner(std::string partner, std::string orientation) : partner(partner), orientation(orientation) {}
+
+        bool operator<(const IntPartner& other) const {
+            if(partner < other.partner) {
+                return true;
+            } else if(partner == other.partner) {
+                return orientation < other.orientation;
+            } else {
+                return false;
+            }
+        }
+    };
+    using IntKey = std::vector<IntPartner>;
 }
 
 #endif //RNANUE_DATATYPES_HPP

include/IBPTree.hpp

@@ -33,10 +33,12 @@ class IBPTree {
         // tree operations
         Node* getRoot(IntervalData& data);
         void insert(IntervalData& data);
+        void update(Node* node); // updates the keys in parent
         void insertIter(Node* node, IntervalData& data);
         void splitNode(Node* node, int index);
-        std::vector<IntervalData*> search(std::string chrom, dtp::Interval interval);
-        void searchIter(Node* node, const dtp::Interval& interval, std::vector<IntervalData*>& results);
+        std::vector<std::pair<Node*, IntervalData*>> search(std::string chrom, dtp::Interval interval);
+        void searchIter(Node* node, const dtp::Interval& interval,
+                        std::vector<std::pair<Node*,IntervalData*>>& results);
         bool isOverlapping(dtp::Interval intvl1, dtp::Interval intvl2);
 
 

include/Node.hpp

@@ -13,7 +13,7 @@
 class IntervalData {
     public:
         IntervalData(std::string chrom, char strand, std::string id, std::string name,
-                 std::string biotype, dtp::Interval interval);
+                 std::string biotype, dtp::Interval interval, IntervalData* split);
         ~IntervalData();
 
         // operator overloading
@@ -35,6 +35,8 @@ class IntervalData {
         void setInterval(dtp::Interval interval);
         dtp::SpliceJunctions getJunctions() const;
         void setJunctions(dtp::SpliceJunctions junctions);
+        IntervalData* getSplit() const;
+        void setSplit(IntervalData* split);
 
         // operations
         void addJunction(std::string name, std::pair<size_t,size_t> junction);
@@ -50,6 +52,7 @@ class IntervalData {
         std::string biotype;
         dtp::Interval interval;
         dtp::SpliceJunctions junctions;
+        IntervalData* split;
 };
 
 class Node {
@@ -60,7 +63,6 @@ class Node {
         std::vector<Node*> children;
         Node* next; // link to the next node
         Node* parent; // link to the parent node
-        std::vector<Node*> splits; // link to the corresponding split node
         bool isLeaf;
 
         // constructor

include/Utility.hpp

@@ -41,12 +41,13 @@ namespace helper {
     // general operations
     bool withinRange(int a, int b, int range);
     std::string removeNonPrintable(const std::string str);
-
-
-
     std::string getTime(); // reports the current time
 }
 
+namespace stats {
+    double median(std::vector<double>& values);
+}
+
 // sequence input/output
 namespace seqIO {
     bool filterReads(auto& qual, int quality, int minlen);