megahit/compact_sequence.cpp at master · shaman-narayanasamy/megahit · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
/**
 * @file compact_sequence.cpp
 * @brief
 * @author Yu Peng (ypeng@cs.hku.hk)
 * @version 1.0.0
 * @date 2011-08-02
 * @modified by Dinghua Li
 * @date 2014-10-02
 */

#include "compact_sequence.h"

#include <string>
#include <algorithm>
#include "bit_operation.h"

#include <bitset>
#include <iostream>

using namespace std;

const CompactSequence &CompactSequence::Append(const CompactSequence &compact_seq, int offset, size_t length)
{
    if (length == std::string::npos || length > compact_seq.size() - offset)
        length = compact_seq.size() - offset;

    uint32_t old_size = size();
    resize(old_size + length);

    for (unsigned i = 0; i < length; ++i)
        set_base(i + old_size, compact_seq[i + offset]);
    return *this;
}

const CompactSequence &CompactSequence::Append(const std::string &seq, int offset, size_t length)
{
    if (length == std::string::npos || length > seq.size() - offset)
        length = seq.size() - offset;

    uint32_t old_size = size();
    resize(old_size + length);

    for (unsigned i = 0; i < length; ++i)
        set_base(i + old_size, seq[i + offset]);
    return *this;
}

const CompactSequence &CompactSequence::Append(uint8_t ch)
{
    if (data_[data_.size()-1] > 0) {
        data_[data_.size()-2] &= (1 << ((4 - data_[data_.size()-1]) << 1)) - 1;
        data_[data_.size()-2] |= (ch&3) << ((4 - data_[data_.size()-1]) << 1);
        --data_[data_.size()-1];
    } else {
        data_[data_.size()-1] = ch & 3;
        data_.push_back(3);
    }
    return *this;
}

const CompactSequence &CompactSequence::ReverseComplement()
{
    uint8_t remain_bases = data_[data_.size()-1];
    reverse(data_.begin(), data_.begin() + data_.size() - 1);
    for (unsigned i = 0; i+1 < data_.size(); ++i)
        data_[i] = bit_operation::ReverseComplement(data_[i]);
    for (unsigned i = 0; i+2 < data_.size(); ++i)
        data_[i] = (uint8_t(data_[i]) >> (remain_bases<<1)) | (uint8_t(data_[i+1]) << ((4-remain_bases)<<1));
    data_[data_.size()-2] >>= (remain_bases << 1);
    return *this;
}

const CompactSequence &CompactSequence::Reverse()
{
    uint8_t remain_bases = data_[data_.size()-1];
    reverse(data_.begin(), data_.begin() + data_.size() - 1);
    for (unsigned i = 0; i+1 < data_.size(); ++i)
        data_[i] = bit_operation::Reverse(data_[i]);
    for (unsigned i = 0; i+2 < data_.size(); ++i)
        data_[i] = (uint8_t(data_[i]) >> (remain_bases<<1)) | (uint8_t(data_[i+1]) << ((4-remain_bases)<<1));
    data_[data_.size()-2] >>= (remain_bases << 1);
    return *this;
}

std::string CompactSequence::ToDNAString()
{
    static char acgt[] = "ACGT";
    std::string dna;
    dna.resize(this->size());
    for (unsigned i = 0; i < this->size(); ++i) {
        dna[i] = acgt[this->get_base(i)];
    }
    return dna;
}