added motifTree

Makefile

@@ -51,6 +51,7 @@ SOURCES = align\
 		merge\
 		merge_candidate_variants\
 		motif_suffix_tree\
+		motif_tree\
 		ordered_bcf_overlap_matcher\
 		ordered_region_overlap_matcher\
 		partition\

filter.h

@@ -69,6 +69,8 @@
 #define VT_VARIANT_LEN   (39|VT_INT|VT_BCF_OP)
 #define VT_N_FILTER      (40|VT_INT|VT_BCF_OP)
 
+//problems will arise once you pass 63.
+
 #define VT_UNKNOWN -1
 
 /**

hfilter.cpp

@@ -157,7 +157,7 @@ class Igor : Program
     {
         std::clog << "\n";
         std::cerr << "stats: no. of variants filtered     " << no_variants_filtered << "\n";
-        std::cerr << "       total no. of variants         " << no_variants << "\n";std::clog << "\n";
+        std::cerr << "       total no. of variants        " << no_variants << "\n";std::clog << "\n";
     }
 
     void hfilter()

motif_tree.cpp

@@ -0,0 +1,175 @@
+/* The MIT License
+
+   Copyright (c) 2014 Adrian Tan <atks@umich.edu>
+
+   Permission is hereby granted, free of charge, to any person obtaining a copy
+   of this software and associated documentation files (the "Software"), to deal
+   in the Software without restriction, including without limitation the rights
+   to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+   copies of the Software, and to permit persons to whom the Software is
+   furnished to do so, subject to the following conditions:
+
+   The above copyright notice and this permission notice shall be included in
+   all copies or substantial portions of the Software.
+
+   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+   AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+   OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+   THE SOFTWARE.
+*/
+
+#include "motif_tree.h"
+
+#define A 1
+#define C 2
+#define G 4
+#define T 8
+#define N 15
+
+/**
+ * Constructor.
+ */
+MotifTree::MotifTree()
+{
+    uint32_t size = 1<<16;
+
+    std::cerr << "size : " << size << "\n";
+
+
+    tree = (uint64_t *) malloc(sizeof(uint64_t)*size);
+
+    uint64_t value = 0;
+    //you only want to loop through multiples of 2
+    for (uint32_t i=0; i<size; ++i)
+    {
+        uint32_t c = canonical(i);
+        tree[i] = ((uint64_t) c)<<32;
+        std::cerr << i << ":" << c << ":" << tree[i] << "\n";
+        if (i==1) exit(1);
+    }
+
+
+};
+
+/**
+ * Destructor.
+ */
+MotifTree::~MotifTree()
+{
+    if (tree) delete tree;
+};
+
+/**
+ * Construct suffix tree based on sequence.
+ */
+void MotifTree::set_sequence(char* sequence)
+{
+    //translate sequence to binary form
+    uint32_t len = strlen(sequence);
+
+
+
+};
+
+/**
+ * Construct suffix tree based on sequence up to max_motif_len.
+ */
+void MotifTree::set_sequence(char* sequence, int32_t max_motif_len)
+{
+
+};
+
+/**
+ * Gets candidate motifs up to max_motif_len.
+ */
+void MotifTree::get_candidate_motifs(std::vector<CandidateMotif>& candidate_motifs)
+{
+
+};
+
+/**
+ * Get canonical representation.
+ */
+uint32_t MotifTree::canonical(uint32_t motif)
+{
+    uint32_t cmotif = motif;
+    uint32_t smotif = motif;
+    std::cerr << "\t" << 0 << ") " << smotif << " - ";
+    print(smotif);
+    std::cerr << "\n";
+    for (uint32_t i=1; i<8; ++i)
+    {
+        smotif = shift1(smotif);
+        std::cerr << "\t" << i << ") " << smotif << " - ";
+        print(smotif);
+        std::cerr << "\n";
+        cmotif = smotif<cmotif ? smotif : cmotif;
+    }
+
+    return cmotif;
+}
+
+/**
+ * Adds a suffix of sequence from start to end.
+ */
+void MotifTree::add_suffix(char* sequence, int32_t start, int32_t end)
+{
+
+};
+
+/**
+ * Converts base to index.
+ */
+int32_t MotifTree::base2index(char base)
+{
+    switch (base)
+    {
+        case 'A':
+            return A;
+            break;
+        case 'C':
+            return C;
+            break;
+        case 'G':
+            return G;
+            break;
+        case 'T':
+            return T;
+            break;
+        default:
+            return N;
+    }
+};
+
+/**
+ * Print sequence.
+ */
+void MotifTree::print(uint32_t seq)
+{
+    uint8_t *seq_ptr = (uint8_t*) &seq;
+
+    std::cerr << index2base(seq_ptr[0] & 0xF);
+    std::cerr << index2base(seq_ptr[0] >> 4 & 0xF);
+    std::cerr << index2base(seq_ptr[1] & 0xF);
+    std::cerr << index2base(seq_ptr[1] >> 4 & 0xF);
+    std::cerr << index2base(seq_ptr[2] & 0xF);
+    std::cerr << index2base(seq_ptr[2] >> 4 & 0xF);
+    std::cerr << index2base(seq_ptr[3] & 0xF);
+    std::cerr << index2base(seq_ptr[3] >> 4 & 0xF);
+
+//    for (uint32_t i=0; i<8; ++i)
+//    {
+//        //seqi(s, i) ((s)[(i)>>1] >> ((~(i)&1)<<2) & 0xf)
+//        std::cerr << "(" << seqi(seq_ptr, i) << ")" ;
+//        std::cerr << index2base(seqi(seq_ptr, i));
+//    }
+}
+
+#undef A
+#undef C
+#undef G
+#undef T
+#undef N

motif_tree.h

@@ -0,0 +1,122 @@
+/* The MIT License
+
+   Copyright (c) 2014 Adrian Tan <atks@umich.edu>
+
+   Permission is hereby granted, free of charge, to any person obtaining a copy
+   of this software and associated documentation files (the "Software"), to deal
+   in the Software without restriction, including without limitation the rights
+   to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+   copies of the Software, and to permit persons to whom the Software is
+   furnished to do so, subject to the following conditions:
+
+   The above copyright notice and this permission notice shall be included in
+   all copies or substantial portions of the Software.
+
+   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+   AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+   OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+   THE SOFTWARE.
+*/
+
+#ifndef MOTIF_TREE_H
+#define MOTIF_TREE_H
+
+#include <cstdint>
+#include <cstring>
+#include <vector>
+#include <list>
+#include <iostream>
+#include <queue>
+#include "candidate_motif.h"
+//#include <arpa/net.h>
+
+#define A 1
+#define C 2
+#define G 4
+#define T 8
+#define N 15
+
+#define shift1(m) (((0x0FFFFFFF&(m))<<4) | ((0xF0000000&(m))>>28))
+#define shift2(m) ((0x00FFFFFF&(m)<<8) | (0xFF000000&(m)>>24))
+#define shift3(m) ((0x000FFFFF&(m)<<12) | (0xFFF00000&(m)>>20))
+#define shift4(m) ((0x0000FFFF&(m)<<16) | (0xFFFF0000&(m)>>16))
+#define shift5(m) ((0x00000FFF&(m)<<20) | (0xFFFFF000&(m)>>12))
+#define shift6(m) ((0x000000FF&(m)<<24) | (0xFFFFFF00&(m)>>8))
+#define shift7(m) ((0x0000000F&(m)<<28) | (0xFFFFFFF0&(m)>>4))
+
+#define seqi(s, i) ((s)[(i)>>1] >> ((~(i)&1)<<2) & 0xf)
+
+#define index2base(i) ("OACXGXXXTXXXXXXN"[(i)])
+
+/**
+ * Motif Suffix Tree for selecting candidate motifs.
+ */
+class MotifTree
+{
+    public:
+    uint64_t* tree;
+
+    /**
+     * Constructor.
+     */
+    MotifTree();
+
+    /**
+     * Destructor.
+     */
+    ~MotifTree();
+
+    /**
+     * Clear the suffix tree.
+     */
+    void clear();
+
+    /**
+     * Construct suffix tree based on sequence.
+     */
+    void set_sequence(char* sequence);
+
+    /**
+     * Construct suffix tree based on sequence up to max_motif_len.
+     */
+    void set_sequence(char* sequence, int32_t max_motif_len);
+
+    /**
+     * Gets candidate motifs up to max_motif_len.
+     */
+    void get_candidate_motifs(std::vector<CandidateMotif>& candidate_motifs);
+
+    /**
+     * Get canonical representation.
+     */
+    uint32_t canonical(uint32_t motif);
+
+    private:
+
+    /**
+     * Adds a suffix of sequence from start to end.
+     */
+    void add_suffix(char* sequence, int32_t start, int32_t end);
+
+    /**
+     * Converts base to index.
+     */
+    int32_t base2index(char base);
+
+    /**
+     * Print sequence.
+     */
+    void print(uint32_t seq);
+
+};
+
+#undef A
+#undef C
+#undef G
+#undef T
+#undef N
+
+#endif

view.cpp

@@ -100,7 +100,6 @@ class Igor : Program
             input_vcf_file = arg_input_vcf_file.getValue();
             output_vcf_file = arg_output_vcf_file.getValue();
             parse_intervals(intervals, arg_interval_list.getValue(), arg_intervals.getValue());
-            //read_sample_list(samples, arg_sample_list.getValue());
             fexp = arg_fexp.getValue();
             print_header = arg_print_header.getValue();
             print_header_only = arg_print_header_only.getValue();