isoseqAlgn.hpp

Go to the documentation of this file.

/*
 * Copyright (c) 2024-2025 Anthony J. Greenberg and Rebekah Rogers
 *
 * Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
 *
 * 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
 * THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS
 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
 * IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
 * THE POSSIBILITY OF SUCH DAMAGE.
 */


 
#pragma once
 
#include <memory>
#include <utility>
#include <vector>
#include <array>
#include <set>
#include <unordered_map>
#include <string>
#include <filesystem>
 
#include "sam.h"
#include "bgzf.h"
 
namespace isaSpace {
    struct BamAndGffFiles;
    struct StatsAndFastqFiles;
    struct TokenAttibuteListPair;
    struct MappedReadInterval;
    struct MappedReadMatchStatus;
    struct BinomialWindowParameters;
    struct ReadExonCoverage;
    struct ReadPortion;
    struct BAMsecondary;
    struct BAMheaderDeleter;
    struct BAMrecordDeleter;
    struct BGZFhandleDeleter;
    class  BAMsafeReader;
    class  ExonGroup;
    class  ReadMatchWindowBIC;
    class  BAMrecord;
    class  BAMtoGenome;
    class  BAMfile;

    struct BamAndGffFiles {
        std::string bamFileName;
        std::string gffFileName;
    };

    struct StatsAndFastqFiles {
        std::string statsFileName;
        std::string fastqFileName;
    };

    struct TokenAttibuteListPair {
        std::string tokenName;
        std::vector<std::string> attributeList;
    };

    struct MappedReadInterval {
        hts_pos_t readStart{-1};
        hts_pos_t readEnd{-1};
        hts_pos_t referenceStart{-1};
        hts_pos_t referenceEnd{-1};
    };

    struct MappedReadMatchStatus {
        hts_pos_t mapStart{-1};
        std::vector<float> matchStatus;
    };

    struct BinomialWindowParameters {
        float currentProbability{0.0F};
        float alternativeProbability{0.0F};
        float bicDifferenceThreshold{0.0F};
        std::vector< std::pair<float, hts_pos_t> >::difference_type windowSize{0};
    };

    struct ReadExonCoverage {
        std::string chromosomeName;
        std::string readName;
        hts_pos_t alignmentStart{0};
        hts_pos_t alignmentEnd{0};
        hts_pos_t bestAlignmentStart{0};
        hts_pos_t bestAlignmentEnd{0};
        uint32_t firstSoftClipLength{0};
        uint16_t nSecondaryAlignments{0};
        uint16_t nGoodSecondaryAlignments{0};
        uint16_t nLocalReversedAlignments{0};
        uint16_t nExons{0};
        std::string geneName;
        char strand{'+'};
        hts_pos_t firstExonLength{-1};
        hts_pos_t firstExonStart{-1};
        hts_pos_t lastExonEnd{-1};
        std::vector<float> exonCoverageScores;
        std::vector<float> bestExonCoverageScores;
    };

    struct ReadPortion {
        size_t start{0};
        size_t end{0};
        std::string originalName;
    };

    struct BAMsecondary {
        hts_pos_t mapStart{-1};
        hts_pos_t mapEnd{-1};
        bool sameStrandAsPrimary{false};
        std::vector<uint32_t> cigar;
    };

    struct BAMheaderDeleter {
        void operator()(sam_hdr_t *header) const { sam_hdr_destroy(header); }
    };

    struct BAMrecordDeleter {
        void operator()(bam1_t *bamRecord) const { bam_destroy1(bamRecord); }
    };

    struct BGZFhandleDeleter {
        void operator()(BGZF *bgzfHandle) const { bgzf_close(bgzfHandle); }
    };

    class BAMsafeReader {
    public:
        BAMsafeReader() = default;
        BAMsafeReader(const std::string &bamFileName);
        BAMsafeReader(const BAMsafeReader &toCopy) = delete;
        BAMsafeReader& operator=(const BAMsafeReader &toCopy) = delete;
        BAMsafeReader(BAMsafeReader &&toMove) noexcept {
            *this = std::move(toMove);
        };

        BAMsafeReader& operator=(BAMsafeReader &&toMove) noexcept;
        ~BAMsafeReader();

        [[nodiscard]] std::unique_ptr<sam_hdr_t, BAMheaderDeleter> getHeaderCopy() const;
        [[nodiscard]] std::pair<std::unique_ptr<bam1_t, BAMrecordDeleter>, int32_t> getNextRecord();
    private:
        static const uint16_t nRetries_;
        BGZF *fileHandle_{nullptr};
        std::unique_ptr<sam_hdr_t, BAMheaderDeleter> headerUPointer_{nullptr};
        std::string fileName_;
        std::filesystem::perms initialPermissions_{std::filesystem::perms::none};
    };

    class ExonGroup {
    public:
        ExonGroup() = default;
        ExonGroup(std::string geneName, const char strand, std::set< std::pair<hts_pos_t, hts_pos_t> > &exonSet);
        ExonGroup(std::string geneName, const char strand, const std::vector<std::string> &exonLinesFomGFF);
        ExonGroup(const ExonGroup &toCopy) = default;
        ExonGroup& operator=(const ExonGroup &toCopy) = default;
        ExonGroup(ExonGroup &&toMove) noexcept = default;
        ExonGroup& operator=(ExonGroup &&toMove) noexcept = default;
        ~ExonGroup() = default;

        [[nodiscard]] bool empty() const noexcept { return exonRanges_.empty(); };
        [[nodiscard]] std::string geneName() const { return geneName_; };
        [[nodiscard]] size_t nExons() const noexcept { return exonRanges_.size(); };
        [[nodiscard]] char strand() const noexcept { return isNegativeStrand_ ? '-' : '+'; };
        [[nodiscard]] std::pair<hts_pos_t, hts_pos_t> at(const size_t &idx) const {return exonRanges_.at(idx);};
        [[nodiscard]] std::pair<hts_pos_t, hts_pos_t> geneSpan() const noexcept;
        [[nodiscard]] std::pair<hts_pos_t, hts_pos_t> firstExonSpan() const noexcept;
        [[nodiscard]] hts_pos_t firstExonLength() const noexcept;
        [[nodiscard]] uint32_t firstExonAfter(const hts_pos_t &position) const noexcept;
        [[nodiscard]] uint32_t firstOverlappingExon(const hts_pos_t &position) const noexcept;
        [[nodiscard]] uint32_t lastExonBefore(const hts_pos_t &position) const noexcept;
        [[nodiscard]] uint32_t lastOverlappingExon(const hts_pos_t &position) const noexcept;
        [[nodiscard]] std::pair<hts_pos_t, hts_pos_t> getFirstIntronSpan() const;
        [[nodiscard]] std::vector<float> getExonCoverageQuality(const BAMrecord &alignment) const;
        [[nodiscard]] std::vector<float> getBestExonCoverageQuality(const BAMrecord &alignment) const;
    private:
        static const size_t strandIDidx_;
        static const char   gffDelimiter_;
        static const size_t spanStart_;
        static const size_t spanEnd_;
        std::string geneName_;
        bool isNegativeStrand_{false};
        std::vector< std::pair<hts_pos_t, hts_pos_t> > exonRanges_;
        // TODO: to implement tracking transcripts, additional vector of iterators pointing to the next exon
        // Keep the number the same and equal to the number of transcripts listed in the GFF
    };

    class ReadMatchWindowBIC {
    public:
        ReadMatchWindowBIC() = default;
        ReadMatchWindowBIC(const std::vector< std::pair<float, hts_pos_t> >::const_iterator &windowBegin, const BinomialWindowParameters &windowParameters) :
            leftProbability_{windowParameters.currentProbability},
            rightProbability_{windowParameters.alternativeProbability},
            kSuccesses_{calculateKsuccesses_(windowBegin, windowParameters)},
            nTrials_{static_cast<float>(windowParameters.windowSize)} {};

        ReadMatchWindowBIC(const ReadMatchWindowBIC &toCopy) = default;
        ReadMatchWindowBIC& operator=(const ReadMatchWindowBIC &toCopy) = default;
        ReadMatchWindowBIC(ReadMatchWindowBIC &&toMove) noexcept = default;
        ReadMatchWindowBIC& operator=(ReadMatchWindowBIC &&toMove) noexcept = default;
        ~ReadMatchWindowBIC() = default;

        [[nodiscard]] float getBICdifference() const noexcept;
    private:
        float leftProbability_{0.0F};
        float rightProbability_{0.0F};
        float kSuccesses_{0.0F};
        float nTrials_{0.0F};

        [[nodiscard]] static float calculateKsuccesses_(const std::vector< std::pair<float, hts_pos_t> >::const_iterator &windowBegin, const BinomialWindowParameters &windowParameters);
    };

    class BAMrecord {
    public:
        BAMrecord() = default;
        BAMrecord(const bam1_t *alignmentRecord, const sam_hdr_t *samHeader);
        BAMrecord(const BAMrecord &toCopy) = default;
        BAMrecord& operator=(const BAMrecord &toCopy) = default;
        BAMrecord(BAMrecord &&toMove) noexcept = default;
        BAMrecord& operator=(BAMrecord &&toMove) noexcept = default;
        ~BAMrecord() = default;

        void addSecondaryAlignment(const bam1_t *alignmentRecord, const sam_hdr_t *samHeader, const hts_pos_t localWindow = 500'000);
        [[nodiscard]] std::string getReadName() const {return readName_; };
        [[nodiscard]] hts_pos_t getMapStart() const noexcept { return mapStart_; };
        [[nodiscard]] hts_pos_t getMapEnd() const noexcept { return mapEnd_; };
        [[nodiscard]] hts_pos_t getmRNAstart() const noexcept {
            return isRev_ ? mapEnd_ : mapStart_;
        };

        [[nodiscard]] bool isRevComp() const noexcept { return isRev_; };
        [[nodiscard]] bool isMapped() const noexcept { return isMapped_; };
        [[nodiscard]] bool hasSecondaryAlignments() const noexcept { return secondaryAlignmentCount_ > 0; };
        [[nodiscard]] bool hasLocalSecondaryAlignments() const noexcept { return !localSecondaryAlignments_.empty(); };
        [[nodiscard]] uint16_t secondaryAlignmentCount() const noexcept { return secondaryAlignmentCount_; };
        [[nodiscard]] uint16_t localSecondaryAlignmentCount() const noexcept { return static_cast<uint16_t>( localSecondaryAlignments_.size() ); };
        [[nodiscard]] uint16_t localReversedSecondaryAlignmentCount() const noexcept;
        [[nodiscard]] hts_pos_t getReadLength() const noexcept { return static_cast<hts_pos_t>( sequenceAndQuality_.size() ); };
        [[nodiscard]] std::vector<uint32_t> getCIGARvector() const { return cigar_; };
        [[nodiscard]] std::string getCIGARstring() const;
        [[nodiscard]] uint32_t getFirstCIGAR() const noexcept; 
        [[nodiscard]] std::string getReferenceName() const {return referenceName_; };
        [[nodiscard]] MappedReadMatchStatus getBestReferenceMatchStatus() const;
        [[nodiscard]] std::vector< std::pair<float, hts_pos_t> > getReadCentricMatchStatus() const;
        [[nodiscard]] std::vector<MappedReadInterval> getPoorlyMappedRegions(const BinomialWindowParameters &windowParameters) const;
        [[nodiscard]] std::string getSequenceAndQuality(const MappedReadInterval &segmentBoundaries) const;
    private:
        static const uint16_t sequenceMask_;
        static const uint16_t qualityShift_;
        static const uint16_t suppSecondaryAlgn_;
        static const std::array<hts_pos_t, 10> queryConsumption_;
        static const std::array<hts_pos_t, 10> referenceConsumption_;
        static const std::array<char, 16> seqNT16str_;
        static const std::array<char, 16> complSeqNT16str_;
        static const std::array<float, 10> sequenceMatch_;
        bool isMapped_{false};
        bool isRev_{false};
        uint16_t secondaryAlignmentCount_{0};
        hts_pos_t mapStart_{0};
        hts_pos_t mapEnd_{0};
        std::string readName_;
        std::string referenceName_;
        std::vector<uint32_t> cigar_;
        std::vector<uint16_t> sequenceAndQuality_;
        std::vector<BAMsecondary> localSecondaryAlignments_;
    };

    class BAMtoGenome {
    public:
        BAMtoGenome() = default;
        BAMtoGenome(const BamAndGffFiles &bamGFFfilePairNames);
        BAMtoGenome(const BAMtoGenome &toCopy) = default;
        BAMtoGenome& operator=(const BAMtoGenome &toCopy) = default;
        BAMtoGenome(BAMtoGenome &&toMove) noexcept = default;
        BAMtoGenome& operator=(BAMtoGenome &&toMove) noexcept = default;
        ~BAMtoGenome() = default;

        [[nodiscard]] size_t nChromosomes() const noexcept;
        [[nodiscard]] size_t nExonSets() const noexcept;
        void saveReadCoverageStats(const std::string &outFileName, const size_t &nThreads) const;
        void saveUnmappedRegions(const std::string &outFileName, const BinomialWindowParameters &windowParameters, const size_t &nThreads) const;
        void saveUnmappedRegions(const StatsAndFastqFiles &outFilePair, const BinomialWindowParameters &windowParameters, const size_t &nThreads) const;
    private:
        static const uint16_t suppSecondaryAlgn_;

        std::vector< std::pair<BAMrecord, ExonGroup> > readsAndExons_;

        [[nodiscard]] std::vector<ExonGroup>::const_iterator findOverlappingGene_(const std::vector<ExonGroup> &chromosomeExonGroups,
                const std::vector<ExonGroup>::const_iterator &exonGroupSearchStart, BAMrecord &alignedRead);
        void processPrimaryAlignment_(const std::string &referenceName, const BAMrecord &alignmentRecord, std::unordered_map<std::string, std::vector<ExonGroup>::const_iterator> &latestExonGroupIts);
        void processSecondaryAlignment_(const std::string &referenceName, const BAMrecord &alignmentRecord, const std::unordered_map<std::string, std::vector<ExonGroup>::const_iterator> &latestExonGroupIts);
    };

    class BAMfile {
    public:
        BAMfile() = default;
        BAMfile(const std::string &BAMfileName);
        BAMfile(const BAMfile &toCopy) = delete;
        BAMfile& operator=(const BAMfile &toCopy) = delete;
        BAMfile(BAMfile &&toMove) noexcept = default;
        BAMfile& operator=(BAMfile &&toMove) noexcept = default;
        ~BAMfile() = default;

        [[nodiscard]] size_t getPrimaryAlignmentCount() const noexcept;
        void addRemaps(const std::string &remapBAMfileName, const float &remapIdentityCutoff);
        [[nodiscard]] std::vector<std::string> saveRemappedBAM(const std::string &outputBAMfileName) const;
        [[nodiscard]] std::vector<std::string> saveSortedRemappedBAM(const std::string &outputBAMfileName) const;
    private:
        static const uint16_t secondaryOrUnmappedAlgn_;
        std::unordered_map<
            std::string,
            std::unordered_map<
                std::string,
                std::vector< std::unique_ptr<bam1_t, BAMrecordDeleter> >
            >
        > bamRecords_;
        std::unique_ptr<sam_hdr_t, BAMheaderDeleter> bamFileHeader_;
    };
}