JunctionCompare.hpp

Go to the documentation of this file.

// Copyright 2010 Complete Genomics, Inc.
//
// Licensed under the Apache License, Version 2.0 (the "License"); you
// may not use this file except in compliance with the License. You
// may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
// implied. See the License for the specific language governing
// permissions and limitations under the License.

#ifndef CGA_TOOLS_JUNCTION_COMPARE_HPP_
#define CGA_TOOLS_JUNCTION_COMPARE_HPP_



#include "cgatools/core.hpp"
#include "cgatools/junctions/Junction.hpp"
#include "cgatools/util/GenericHistogram.hpp"

// Standard library.
#include <string>
#include <vector>
#include <map>


namespace cgatools { namespace junctions {

    class JunctionComparatorWithinContig
    {
    public:
        static const char SEPARATOR = '\t';

        typedef std::vector<size_t> JunctionIndices;
        typedef std::pair<int,int> JunctionGroupId;
        typedef std::map< JunctionGroupId, JunctionIndices > JunctionGroups;

        JunctionComparatorWithinContig(
            const ReferenceGenome& reference,
            JunctionFile& file1,
            JunctionFile& file2);

        void compareWithinContig(
            const JunctionFile& file1,
            const JunctionFile& file2,
            const JunctionIndices& indexes1,
            const JunctionIndices& indexes2,
            JunctionFile& result1and2,
            JunctionFile& result1not2,
            std::ostream &report
            ) const;

        void compare();

        // Compare two junction files 1 and 2.
        // 1and2 = junctions in 1 that are also in 2.
        // 2and1 = junctions in 2 that are also in 1.
        // 1not2 = junctions in 1 that are not in 2.
        // 2not1 = junctions in 2 that are not in 1.
        // Semi-junctions are ignored and don't appear in the output.
        // This function also canonicalizes the inputs.

        const JunctionFile &getResult1and2() const { return result1and2_; }
        const JunctionFile &getResult1not2() const { return result1not2_; }
        const JunctionFile &getResult2and1() const { return result2and1_; }
        const JunctionFile &getResult2not1() const { return result2not1_; }

    protected:

        void collectJunctionGroups(JunctionFile &junctFile, JunctionGroups &result);

        size_t distanceTolerance_;

        const ReferenceGenome& reference_;

        JunctionFile& file1_;
        JunctionFile& file2_;

        JunctionGroups junctionGroups1_;
        JunctionGroups junctionGroups2_;

        JunctionFile result1and2_;
        JunctionFile result1not2_;
        JunctionFile result2and1_;
        JunctionFile result2not1_;
    };

    #define ADD_FNAME(names, N) names[N] = #N

    class JunctionCounterFields {
    public:
        enum Fields {
            TOTAL_JUNCTIONS = 0,
            LEFT_SIDES_EQUAL,
            RIGHT_SIDES_EQUAL,
            BOTH_SIDES_EQUAL,
            LEFT_UNKNOWN,
            RIGHT_UNKNOWN,
            STRAND_BOTH_POSITIVE,
            STRAND_LEFT_NEGATIVE,
            STRAND_RIGHT_NEGATIVE,
            STRAND_BOTH_NEGATIVE,

            NUMBER_OF_FIELDS
        };
        typedef boost::array<std::string,NUMBER_OF_FIELDS> Names;

        static Names getFieldNames()
        {
            Names names;
            ADD_FNAME(names, TOTAL_JUNCTIONS);
            ADD_FNAME(names, LEFT_SIDES_EQUAL);
            ADD_FNAME(names, RIGHT_SIDES_EQUAL);
            ADD_FNAME(names, BOTH_SIDES_EQUAL);
            ADD_FNAME(names, LEFT_UNKNOWN);
            ADD_FNAME(names, RIGHT_UNKNOWN);
            ADD_FNAME(names, STRAND_BOTH_POSITIVE);
            ADD_FNAME(names, STRAND_LEFT_NEGATIVE);
            ADD_FNAME(names, STRAND_RIGHT_NEGATIVE);
            ADD_FNAME(names, STRAND_BOTH_NEGATIVE);
            return names;
        }
    };


    template <class Fields, class Type>
    class CounterSet
    {
    public:
        CounterSet() {
            values_.assign(Type(0));
        }

        void print(std::ostream &out) const {
            typename Fields::Names fieldNames = Fields::getFieldNames();
            for (size_t i=0; i<values_.size(); ++i)
                out << fieldNames[i] << '\t' << values_[i] << std::endl;
        }

        Type &operator[] (size_t i) {CGA_ASSERT_L(i,values_.size()); return values_[i];}
        boost::array<Type,Fields::NUMBER_OF_FIELDS> values_;
    };

    class JunctionIt 
    {
    public:
        typedef Junctions::const_iterator It;
        JunctionIt() :it_(nullValue_.begin()) {}
        JunctionIt(It it) :it_(it) {}

        const Junction& operator* () const {return *it_;}
        bool operator< (const JunctionIt& other) {
            CGA_ERROR_EX("JunctionIt::operator< is disabled"); return it_<other.it_;
        }

        It                      it_;
        std::vector<JunctionIt> matchedJunctionIts_;

        bool isNull() const {return it_->id_==nullValue_.front().id_ && matchedJunctionIts_.empty();}
    public:
        static Junctions nullValue_;

        friend std::ostream& operator<<(std::ostream &ostr, const JunctionIt::It& j);
        friend std::ostream& operator<<(std::ostream &ostr, const JunctionIt& j);
    };


    //The input junctions have to be canonicalized
    class CompareBySide {
    public:
        CompareBySide(int firstSide=JUNCTION_LEFT_SIDE) 
            :firstSide_(firstSide),
            secondSide_(JUNCTION_RIGHT_SIDE-firstSide)
        {}
        static Location getLocation(const Junction& j, int side) {
            Location l = j.sideSections_[side].position_;
            //if (side==JUNCTION_LEFT_SIDE)
            //    l.offset_ += j.sideSections_[side].length_;
            return l;
        }
        bool operator() (const Junction& j1,const Junction& j2) const {
            //if (j1.id_.id_[0]=="W" && j1.id_.id_[1]=="1695") {
            //    int k=0; ++k;
            //}
            Location l1 = getLocation(j1, firstSide_);
            Location l2 = getLocation(j2, firstSide_);
            if (l1!=l2)
                return l1<l2;
            return getLocation(j1, secondSide_) < getLocation(j2, secondSide_);
        }
        bool operator() (JunctionIt::It j1,JunctionIt::It j2) const {
            return operator() (*j1,*j2);
        }
        bool operator() (const JunctionIt &j1,const JunctionIt &j2) const {
            return operator() (*j1,*j2);
        }
    protected:
        int firstSide_;
        int secondSide_;
    };

    class JuctionStatisticsCollector 
    {
    public:

        JuctionStatisticsCollector(uint32_t distanceTolerance)
            :distanceTolerance_(distanceTolerance),out_(NULL)
        {
            size_t buckets[7][3] = 
            {
                {0,101,10},
                {200, 1001, 100},
                {2000,10001,1000},
                {20000,   100001,     10000},
                {2000000, 10000001,   1000000},
                {20000000,100000001,  10000000},
                {200000000,1000000001,100000000}
            };
            
            scoreHistogram_.addRange(0,20,1);
            scoreHistogram_.addRange(20,200,10);

            for (size_t i=0; i<6; ++i)
                oneSideNeighbourDist_.addRange(buckets[i][0],buckets[i][1],buckets[i][2]);

            searchDistanceTolerance_ = distanceTolerance*5;
            twoSideNeighbourDist_.addRange(0,searchDistanceTolerance_+10,10);

            for (size_t i=0; i<4; ++i)
                clusterSizeInBases_.addRange(buckets[i][0],buckets[i][1],buckets[i][2]);

            clusterSizeInJunctions_.addRange(0,15,1);
        }

        void findOneSideNeighbours(const Junctions& junctions, int firstSide);
        void findFullNeighbours(const Junctions& junctions);
        void run(const Junctions& junctions, std::ostream &out);
    protected:
        uint32_t    distanceTolerance_;
        uint32_t    searchDistanceTolerance_;
        CounterSet<JunctionCounterFields,int64_t> counters_;
        util::GenericHistogram<size_t,size_t> scoreHistogram_;
        util::GenericHistogram<size_t,size_t> oneSideNeighbourDist_;
        util::GenericHistogram<size_t,size_t> twoSideNeighbourDist_;

        util::GenericHistogram<size_t,size_t> clusterSizeInBases_;
        util::GenericHistogram<size_t,size_t> clusterSizeInJunctions_;
        std::ostream *out_;
    };


    class JunctionIts : public std::vector<JunctionIt> 
    {
    public:
        JunctionIts() : std::vector<JunctionIt>() {}

        JunctionIts(const Junctions& junctions, const std::string &id) 
            : std::vector<JunctionIt>()
        {
            addJunctions(junctions, id);
        }

        void addJunctions(const Junctions& junctions, const std::string &id) {
            reserve(size()+junctions.size());
            for (Junctions::const_iterator it=junctions.begin(),endIt=junctions.end(); it!=endIt;++it) {
                push_back(it);
            }
        }
    };

    class FullGenomeJunctionComparator
    {
    public:
        typedef std::vector<Junctions> JunctionsArray;
        typedef std::vector<JunctionIts> JunctionsItsArray;


        FullGenomeJunctionComparator(
            uint32_t distanceTolerance,
            const std::vector<size_t> &scoreThresholds,
            std::ostream &outCommonJunctionsStream
        )
        :   distanceTolerance_(distanceTolerance),
            scoreThresholds_(scoreThresholds),
            outCommonJunctionsStream_(outCommonJunctionsStream)
        {}

        static const char SEPARATOR = '\t';

        //returns:
        // [0,1) - incompatible junctions
        // [1,2) - one side is compatible
        // [2,3) - both sides are compatible
        double getCompatibility(const Junction& j0, const Junction& j1) const {
            double result = 0;
            uint64_t totDistance = 0;
            if (size_t(j0.score_.score_+1E-5) >= scoreThresholds_[0] && 
                size_t(j1.score_.score_+1E-5) > scoreThresholds_[1]) 
            {
                for (size_t i=size_t(JUNCTION_LEFT_SIDE); i<=size_t(JUNCTION_RIGHT_SIDE); ++i) {
                    uint32_t dist = std::abs(CompareBySide::getLocation(j0,i).distanceTo(
                        CompareBySide::getLocation(j1,i)));
                    if (dist < distanceTolerance_
                        && j0.sideSections_[i].strand_ == j1.sideSections_[i].strand_)
                        ++result;
                    totDistance+=dist*dist;
                }
                result+= 10/(sqrt((double)totDistance)+20); //put the weight value into [0,1/2]
                //double totScore = j0.score_.score_*j0.score_.score_ + j1.score_.score_*j1.score_.score_;
            }
            return result;
        }

        bool isSideCompatible(const Junction& j0, const Junction& j1, int side) const {
            return  (uint32_t)std::abs(CompareBySide::getLocation(j0,side).distanceTo(
                            CompareBySide::getLocation(j1,side))) < distanceTolerance_
                    && j0.sideSections_[side].strand_ == j1.sideSections_[side].strand_
                    && size_t(j0.score_.score_+1E-5) >= scoreThresholds_[0] 
                    && size_t(j1.score_.score_+1E-5) > scoreThresholds_[1];
        }

        void compare(const JunctionsArray& junctions) const;

        void compareBySide(const JunctionsArray& junctions, int firstSide,
                           JunctionIts &oneSideCompatible,
                           JunctionIts &twoSidesCompatible,
                           JunctionIts &twoSidesIncompatible) const;
        void mergeResults(JunctionsItsArray &input, JunctionIts &output) const;
        void printResults(const JunctionIts& junctionIts, std::ostream & out) const;
    protected:
        uint32_t            distanceTolerance_;
        std::vector<size_t> scoreThresholds_;
        std::ostream &outCommonJunctionsStream_;
    };

}}


#endif //CGA_TOOLS_JUNCTION_COMPARE_HPP_

---