Yaossg/hisat-3n
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
a3817bf0b7
hisat-3n/repeat_builder.cpp
Yaossg 5e601d0401 initial commit
2025-01-18 21:09:52 +08:00

4772 lines
161 KiB
C++
Raw Blame History

/*
* Copyright 2018, Chanhee Park <parkchanhee@gmail.com> and Daehwan Kim <infphilo@gmail.com>
*
* This file is part of HISAT 2.
*
* HISAT 2 is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* HISAT 2 is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with HISAT 2. If not, see <http://www.gnu.org/licenses/>.
*/
#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
#include "timer.h"
#include "aligner_sw.h"
#include "aligner_result.h"
#include "scoring.h"
#include "sstring.h"
#include "repeat_builder.h"
#include "repeat_kmer.h"
#include "bit_packed_array.h"
unsigned int levenshtein_distance(const std::string& s1, const std::string& s2)
{
const std::size_t len1 = s1.size(), len2 = s2.size();
std::vector<unsigned int> col(len2+1), prevCol(len2+1);
for (unsigned int i = 0; i < prevCol.size(); i++)
prevCol[i] = i;
for (unsigned int i = 0; i < len1; i++) {
col[0] = i+1;
for (unsigned int j = 0; j < len2; j++)
// note that std::min({arg1, arg2, arg3}) works only in C++11,
// for C++98 use std::min(std::min(arg1, arg2), arg3)
col[j+1] = std::min( std::min(prevCol[1 + j] + 1, col[j] + 1), prevCol[j] + (s1[i]==s2[j] ? 0 : 1) );
col.swap(prevCol);
}
return prevCol[len2];
}
unsigned int hamming_distance(const string& s1, const string& s2)
{
assert_eq(s1.length(), s2.length());
uint32_t cnt = 0;
for(size_t i = 0; i < s1.length(); i++) {
if(s1[i] != s2[i]) {
cnt++;
}
}
return cnt;
}
string reverse(const string& str)
{
string rev = str;
size_t str_len = str.length();
for(size_t i = 0; i < str_len; i++) {
rev[i] = str[str_len - i - 1];
}
return rev;
}
string reverse_complement(const string& str)
{
string rev = str;
size_t str_len = str.length();
for(size_t i = 0; i < str_len; i++) {
char nt = str[str_len - i - 1];
rev[i] = asc2dnacomp[nt];
}
return rev;
}
template<typename TStr>
TIndexOffU getLCP(const TStr& s,
CoordHelper& coordHelper,
TIndexOffU a,
TIndexOffU b,
TIndexOffU max_len = 1000)
{
TIndexOffU a_end = coordHelper.getEnd(a);
TIndexOffU b_end = coordHelper.getEnd(b);
assert_leq(a_end, s.length());
assert_leq(b_end, s.length());
TIndexOffU k = 0;
while((a + k) < a_end && (b + k) < b_end && k < max_len) {
if(s[a + k] != s[b + k]) {
break;
}
k++;
}
return k;
}
template<typename TStr>
bool isSameSequenceUpto(const TStr& s,
CoordHelper& coordHelper,
TIndexOffU a,
TIndexOffU b,
TIndexOffU upto)
{
TIndexOffU a_end = coordHelper.getEnd(a);
TIndexOffU b_end = coordHelper.getEnd(b);
assert_leq(a_end, s.length());
assert_leq(b_end, s.length());
if(a + upto > a_end || b + upto > b_end)
return false;
for(int k = upto - 1; k >= 0; k--) {
if(s[a + k] != s[b + k]) {
return false;
}
}
return true;
}
bool isSenseDominant(CoordHelper& coordHelper,
const EList<TIndexOffU>& positions,
size_t seed_len)
{
// count the number of seeds on the sense strand
size_t sense_mer_count = 0;
for(size_t i = 0; i < positions.size(); i++) {
if(positions[i] < coordHelper.forward_length())
sense_mer_count++;
}
// there exists two sets of seeds that are essentially the same
// due to sense and antisense strands except palindromic sequences
// choose one and skip the other one
assert_leq(sense_mer_count, positions.size());
size_t antisense_mer_count = positions.size() - sense_mer_count;
if(sense_mer_count < antisense_mer_count) {
return false;
} else if(sense_mer_count == antisense_mer_count) {
assert_geq(positions.back(), coordHelper.forward_length());
TIndexOffU sense_pos = coordHelper.length() - positions.back() - seed_len;
if(positions[0] > sense_pos) return false;
}
return true;
}
static const Range EMPTY_RANGE = Range(1, 0);
struct RepeatRange {
RepeatRange() {
forward = true;
};
RepeatRange(Range r, int id) :
range(r), rg_id(id), forward(true) {};
RepeatRange(Range r, int id, bool fw) :
range(r), rg_id(id), forward(fw) {};
Range range;
int rg_id;
bool forward;
};
Range reverseRange(const Range& r, TIndexOffU size)
{
size_t len = r.second - r.first;
Range rc;
rc.first = size - r.second;
rc.second = rc.first + len;
return rc;
}
string reverseComplement(const string& str)
{
string rc;
for(TIndexOffU si = str.length(); si > 0; si--) {
rc.push_back(asc2dnacomp[str[si - 1]]);
}
return rc;
}
string toMDZ(const EList<Edit>& edits, const string& read)
{
StackedAln stk;
BTDnaString btread;
BTString buf;
btread.install(read.c_str(), true);
stk.init(btread, edits, 0, 0, 0, 0);
stk.buildMdz();
stk.writeMdz(&buf, NULL);
return string(buf.toZBuf());
}
string applyEdits(const string& ref, size_t read_len, const EList<Edit>& edits, const Coord& coord)
{
string read;
size_t ref_pos = coord.off();
size_t read_pos = 0;
for(size_t i = 0; i < edits.size(); i++) {
for(; read_pos < edits[i].pos; read_pos++, ref_pos++) {
read.push_back(ref[ref_pos]);
}
if(edits[i].type == EDIT_TYPE_READ_GAP) {
// delete on read
ref_pos++;
} else if(edits[i].type == EDIT_TYPE_REF_GAP) {
// insert on read
read.push_back(edits[i].qchr);
read_pos++;
} else if(edits[i].type == EDIT_TYPE_MM) {
assert_eq(ref[ref_pos], edits[i].chr);
read.push_back(edits[i].qchr);
read_pos++;
ref_pos++;
} else {
assert(false);
}
}
for(; read_pos < read_len; read_pos++, ref_pos++) {
read.push_back(ref[ref_pos]);
}
return read;
}
template<typename index_t>
static bool compareRepeatCoordByJoinedOff(const RepeatCoord<index_t>& a, const RepeatCoord<index_t>& b)
{
return a.joinedOff < b.joinedOff;
}
template<typename TStr>
string getString(const TStr& ref, TIndexOffU start, size_t len)
{
ASSERT_ONLY(const size_t ref_len = ref.length());
assert_leq(start + len, ref_len);
string s;
for(size_t i = 0; i < len; i++) {
char nt = "ACGT"[ref[start + i]];
s.push_back(nt);
}
return s;
}
template<typename TStr>
void getString(const TStr& ref, TIndexOffU start, size_t len, string& s)
{
s.clear();
ASSERT_ONLY(const size_t ref_len = ref.length());
assert_leq(start + len, ref_len);
for(size_t i = 0; i < len; i++) {
char nt = "ACGT"[ref[start + i]];
s.push_back(nt);
}
}
template<typename TStr>
inline uint8_t getSequenceBase(const TStr& ref, TIndexOffU pos)
{
assert_lt(pos, ref.length());
return (uint8_t)ref[pos];
}
template<typename TStr>
void dump_tstr(const TStr& s)
{
static int print_width = 60;
size_t s_len = s.length();
for(size_t i = 0; i < s_len; i += print_width) {
string buf;
for(size_t j = 0; (j < print_width) && (i + j < s_len); j++) {
buf.push_back("ACGTN"[s[i + j]]);
}
cerr << buf << endl;
}
cerr << endl;
}
size_t getMaxMatchLen(const EList<Edit>& edits, const size_t read_len)
{
size_t max_length = 0;
uint32_t last_edit_pos = 0;
uint32_t len = 0;
if (edits.size() == 0) {
// no edits --> exact match
return read_len;
}
for(size_t i = 0; i < edits.size(); i++) {
if (last_edit_pos > edits[i].pos) {
continue;
}
len = edits[i].pos - last_edit_pos;
if(len > max_length) {
max_length = len;
}
last_edit_pos = edits[i].pos + 1;
}
if (last_edit_pos < read_len) {
len = read_len - last_edit_pos;
if(len > max_length) {
max_length = len;
}
}
return max_length;
}
int max_index(size_t array[4])
{
int max_idx = 0;
for(size_t i = 1; i < 4; i++) {
if(array[max_idx] < array[i]) {
max_idx = i;
}
}
return max_idx;
}
CoordHelper::CoordHelper(TIndexOffU length,
TIndexOffU forward_length,
const EList<RefRecord>& szs,
const EList<string>& ref_names) :
length_(length),
forward_length_(forward_length),
szs_(szs),
ref_namelines_(ref_names)
{
// build ref_names_ from ref_namelines_
buildNames();
buildJoinedFragment();
}
CoordHelper::~CoordHelper()
{
}
void CoordHelper::buildNames()
{
ref_names_.resize(ref_namelines_.size());
for(size_t i = 0; i < ref_namelines_.size(); i++) {
string& nameline = ref_namelines_[i];
for(size_t j = 0; j < nameline.length(); j++) {
char n = nameline[j];
if(n == ' ') {
break;
}
ref_names_[i].push_back(n);
}
}
}
int CoordHelper::mapJoinedOffToSeq(TIndexOffU joinedOff)
{
/* search from cached_list */
if(num_cached_ > 0) {
for(size_t i = 0; i < num_cached_; i++) {
Fragments *frag = &cached_[i];
if(frag->contain(joinedOff)) {
return frag->frag_id;
}
}
/* fall through */
}
/* search list */
int top = 0;
int bot = fraglist_.size() - 1;
int pos = 0;
Fragments *frag = &fraglist_[pos];
while((bot - top) > 1) {
pos = top + ((bot - top) >> 1);
frag = &fraglist_[pos];
if (joinedOff < frag->joinedOff) {
bot = pos;
} else {
top = pos;
}
}
frag = &fraglist_[top];
if (frag->contain(joinedOff)) {
// update cache
if (num_cached_ < CACHE_SIZE_JOINEDFRG) {
cached_[num_cached_] = *frag;
num_cached_++;
} else {
cached_[victim_] = *frag;
victim_++; // round-robin
victim_ %= CACHE_SIZE_JOINEDFRG;
}
return top;
}
return -1;
}
int CoordHelper::getGenomeCoord(TIndexOffU joinedOff,
string& chr_name,
TIndexOffU& pos_in_chr)
{
int seq_id = mapJoinedOffToSeq(joinedOff);
if (seq_id < 0) {
return -1;
}
Fragments *frag = &fraglist_[seq_id];
TIndexOffU offset = joinedOff - frag->joinedOff;
pos_in_chr = frag->seqOff + offset;
chr_name = ref_names_[frag->seq_id];
return 0;
}
void CoordHelper::buildJoinedFragment()
{
TIndexOffU acc_joinedOff = 0;
TIndexOffU acc_seqOff = 0;
TIndexOffU seq_id = 0;
TIndexOffU frag_id = 0;
for(size_t i = 0; i < szs_.size(); i++) {
const RefRecord& rec = szs_[i];
if(rec.len == 0) {
continue;
}
if (rec.first) {
acc_seqOff = 0;
seq_id++;
}
fraglist_.expand();
fraglist_.back().joinedOff = acc_joinedOff;
fraglist_.back().length = rec.len;
fraglist_.back().frag_id = frag_id++;
fraglist_.back().seq_id = seq_id - 1;
fraglist_.back().first = rec.first;
acc_seqOff += rec.off;
fraglist_.back().seqOff = acc_seqOff;
acc_joinedOff += rec.len;
acc_seqOff += rec.len;
}
// Add Last Fragment(empty)
fraglist_.expand();
fraglist_.back().joinedOff = acc_joinedOff;
fraglist_.back().length = 0;
fraglist_.back().seqOff = acc_seqOff;
fraglist_.back().first = false;
fraglist_.back().frag_id = frag_id;
fraglist_.back().seq_id = seq_id;
}
TIndexOffU CoordHelper::getEnd(TIndexOffU e) {
assert_lt(e, length_)
TIndexOffU end = 0;
if(e < forward_length_) {
int frag_id = mapJoinedOffToSeq(e);
assert_geq(frag_id, 0);
end = fraglist_[frag_id].joinedOff + fraglist_[frag_id].length;
} else {
// ReverseComplement
// a, b are positions w.r.t reverse complement string.
// fragment map is based on forward string
int frag_id = mapJoinedOffToSeq(length_ - e - 1);
assert_geq(frag_id, 0);
end = length_ - fraglist_[frag_id].joinedOff;
}
assert_leq(end, length_);
return end;
}
TIndexOffU CoordHelper::getStart(TIndexOffU e) {
assert_lt(e, length_)
TIndexOffU start = 0;
if(e < forward_length_) {
int frag_id = mapJoinedOffToSeq(e);
assert_geq(frag_id, 0);
start = fraglist_[frag_id].joinedOff;
} else {
// ReverseComplement
// a, b are positions w.r.t reverse complement string.
// fragment map is based on forward string
int frag_id = mapJoinedOffToSeq(length_ - e - 1);
assert_geq(frag_id, 0);
start = length_ - (fraglist_[frag_id].joinedOff + fraglist_[frag_id].length);
}
assert_leq(start, length_);
return start;
}
template<typename TStr>
void SeedExt::getExtendedSeedSequence(const TStr& s,
string& seq) const
{
seq.clear();
TIndexOffU prev_end = orig_pos.first;
for(size_t j = 0; j < left_gaps.size(); j++) {
TIndexOffU curr_end = orig_pos.first - left_gaps[j].first;
assert_leq(curr_end, prev_end);
if(curr_end < prev_end) {
seq = getString(s, curr_end, prev_end - curr_end) + seq;
}
int gap_len = left_gaps[j].second;
assert_neq(gap_len, 0);
if(gap_len > 0) { // deletion
string gap_str(gap_len, '-');
seq = gap_str + seq;
} else {
curr_end += gap_len;
}
prev_end = curr_end;
}
assert_leq(pos.first, prev_end);
if(pos.first < prev_end) {
seq = getString(s, pos.first, prev_end - pos.first) + seq;
}
if(orig_pos.second - orig_pos.first > 0)
seq += getString(s, orig_pos.first, orig_pos.second - orig_pos.first);
TIndexOffU prev_begin = orig_pos.second;
for(size_t j = 0; j < right_gaps.size(); j++) {
TIndexOffU curr_begin = orig_pos.second + right_gaps[j].first;
assert_leq(prev_begin, curr_begin);
if(prev_begin < curr_begin) {
seq += getString(s, prev_begin, curr_begin - prev_begin);
}
int gap_len = right_gaps[j].second;
assert_neq(gap_len, 0);
if(gap_len > 0) { // deletion
string gap_str(gap_len, '-');
seq += gap_str;
} else {
curr_begin -= gap_len;
}
prev_begin = curr_begin;
}
assert_leq(prev_begin, pos.second);
if(prev_begin < pos.second) {
seq += getString(s, prev_begin, pos.second - prev_begin);
}
}
SeedSNP *lookup_add_SNP(EList<SeedSNP *>& repeat_snps, SeedSNP& snp)
{
for(size_t i = 0; i < repeat_snps.size(); i++) {
if(*repeat_snps[i] == snp) {
return repeat_snps[i];
}
}
repeat_snps.expand();
repeat_snps.back() = new SeedSNP();
*(repeat_snps.back()) = snp;
return repeat_snps.back();
}
template<typename TStr>
void SeedExt::generateSNPs(const TStr &s, const string& consensus, EList<SeedSNP *>& repeat_snps)
{
EList<pair<TIndexOffU, int> > gaps;
// Merge Gaps
{
TIndexOffU left_ext_len = getLeftExtLength();
for(size_t g = 0; g < left_gaps.size(); g++) {
gaps.expand();
gaps.back().first = left_ext_len - left_gaps[g].first + left_gaps[g].second;
gaps.back().second = left_gaps[g].second;
}
TIndexOffU right_base = orig_pos.second - pos.first;
for(size_t g = 0; g < right_gaps.size(); g++) {
gaps.expand();
gaps.back().first = right_base + right_gaps[g].first;
gaps.back().second = right_gaps[g].second;
}
gaps.sort();
}
//
string con_ref;
string seq_read;
TIndexOffU prev_con_pos = consensus_pos.first;
TIndexOffU prev_seq_pos = pos.first;
for(size_t g = 0; g < gaps.size(); g++) {
TIndexOffU curr_seq_pos = pos.first + gaps[g].first;
TIndexOffU curr_con_pos = prev_con_pos + (curr_seq_pos - prev_seq_pos);
con_ref = consensus.substr(prev_con_pos, curr_con_pos - prev_con_pos);
seq_read = getString(s, prev_seq_pos, curr_seq_pos - prev_seq_pos);
for(size_t l = 0; l < con_ref.length(); l++) {
if(seq_read[l] != con_ref[l]) {
// Single
SeedSNP snp;
snp.init(EDIT_TYPE_MM, prev_con_pos + l, 1, seq_read[l]);
// lookup repeat_snp
snps.expand();
snps.back() = lookup_add_SNP(repeat_snps, snp);
}
}
int gap_len = gaps[g].second;
assert_neq(gap_len, 0);
if(gap_len > 0) {
// Deletion (gap on read)
string gap_str(gap_len, '-');
SeedSNP snp;
snp.init(EDIT_TYPE_READ_GAP, curr_con_pos, gap_len, consensus.substr(curr_con_pos, gap_len));
snps.expand();
snps.back() = lookup_add_SNP(repeat_snps, snp);
curr_con_pos += gap_len;
} else {
// Insertion (gap on reference)
string gap_str(abs(gap_len), '-');
SeedSNP snp;
snp.init(EDIT_TYPE_REF_GAP, curr_con_pos, abs(gap_len), getString(s, curr_seq_pos, abs(gap_len)));
snps.expand();
snps.back() = lookup_add_SNP(repeat_snps, snp);
curr_seq_pos += abs(gap_len);
}
prev_con_pos = curr_con_pos;
prev_seq_pos = curr_seq_pos;
}
assert_eq(consensus_pos.second - prev_con_pos, pos.second - prev_seq_pos);
con_ref = consensus.substr(prev_con_pos, consensus_pos.second - prev_con_pos);
seq_read = getString(s, prev_seq_pos, pos.second - prev_seq_pos);
for(size_t l = 0; l < con_ref.length(); l++) {
if(seq_read[l] != con_ref[l]) {
// Single
SeedSNP snp;
snp.init(EDIT_TYPE_MM, prev_con_pos + l, 1, seq_read[l]);
snps.expand();
snps.back() = lookup_add_SNP(repeat_snps, snp);
}
}
}
bool seedCmp(const SeedExt& s, const SeedExt& s2)
{
if(s.getLength() != s2.getLength())
return s.getLength() > s2.getLength();
if(s.pos.first != s2.pos.first)
return s.pos.first < s2.pos.first;
if(s.pos.second != s2.pos.second)
return s.pos.second < s2.pos.second;
return false;
}
template<typename TStr>
void RB_Repeat::init(const RepeatParameter& rp,
const TStr& s,
CoordHelper& coordHelper,
const RB_SubSA& subSA,
const RB_RepeatBase& repeatBase)
{
consensus_ = repeatBase.seq;
assert_geq(consensus_.length(), rp.min_repeat_len);
assert_eq(consensus_.length(), rp.min_repeat_len + repeatBase.nodes.size() - 1);
EList<TIndexOffU> positions, next_positions;
const EList<TIndexOffU>& repeat_index = subSA.getRepeatIndex();
seeds_.clear();
for(size_t n = 0; n < repeatBase.nodes.size(); n++) {
TIndexOffU idx = repeatBase.nodes[n];
TIndexOffU saBegin = repeat_index[idx];
TIndexOffU saEnd = (idx + 1 < repeat_index.size() ? repeat_index[idx+1] : subSA.size());
next_positions.clear();
for(size_t sa = saBegin; sa < saEnd; sa++) {
TIndexOffU left = subSA.get(sa);
TIndexOffU right = left + subSA.seed_len();
next_positions.push_back(left);
if(left > 0) {
size_t idx = positions.bsearchLoBound(left - 1);
if(idx < positions.size() && positions[idx] == left - 1)
continue;
}
seeds_.expand();
SeedExt& seed = seeds_.back();
seed.reset();
seed.orig_pos = pair<TIndexOffU, TIndexOffU>(left, right);
seed.pos = seed.orig_pos;
seed.consensus_pos.first = n;
seed.consensus_pos.second = n + rp.min_repeat_len;
seed.bound = pair<TIndexOffU, TIndexOffU>(coordHelper.getStart(left), coordHelper.getEnd(left));
#ifndef NDEBUG
string tmp_str = getString(s, seed.pos.first, seed.pos.second - seed.pos.first);
assert_eq(tmp_str, consensus_.substr(n, seed.pos.second - seed.pos.first));
#endif
for(size_t p = seed.consensus_pos.second; p < consensus_.length(); p++) {
size_t pos = seed.pos.second;
if(pos >= seed.bound.second)
break;
int nt = getSequenceBase(s, pos);
if("ACGT"[nt] != consensus_[p])
break;
seed.pos.second++;
seed.consensus_pos.second++;
}
}
positions = next_positions;
}
internal_update();
}
template<typename TStr>
void RB_Repeat::extendConsensus(const RepeatParameter& rp,
const TStr& s)
{
size_t remain = seeds_.size();
EList<string> left_consensuses, right_consensuses, empty_consensuses;
EList<size_t> ed_seed_nums;
const TIndexOffU default_max_ext_len = (rp.max_repeat_len - rp.seed_len) / 2;
const TIndexOffU seed_mm = 1;
string left_ext_consensus = "", right_ext_consensus = "";
empty_consensuses.resize(seed_mm + 1);
empty_consensuses.fill("");
while(remain >= rp.repeat_count) {
for(size_t i = 0; i < remain; i++) {
seeds_[i].done = false;
seeds_[i].curr_ext_len = 0;
}
// extend seeds in left or right direction
TIndexOffU max_ext_len = min(default_max_ext_len, (TIndexOffU)(rp.max_repeat_len - consensus_.length()));
get_consensus_seq(s,
seeds_,
0, // seed begin
remain, // seed end
max_ext_len,
max_ext_len,
seed_mm,
rp,
ed_seed_nums,
&left_consensuses,
&right_consensuses);
size_t allowed_seed_mm = 0;
left_ext_consensus.clear();
right_ext_consensus.clear();
for(int i = 0 /* (int)seed_mm */; i >= 0; i--) {
size_t extlen = (ed_seed_nums[i] < rp.repeat_count ? 0 : (left_consensuses[i].length() + right_consensuses[i].length()));
// if(extlen <= 0 || extlen < max_ext_len * i / seed_mm)
//if(extlen < max_ext_len * 2)
// continue;
if(extlen <= 0)
continue;
left_ext_consensus = left_consensuses[i];
right_ext_consensus = right_consensuses[i];
allowed_seed_mm = (size_t)i;
assert_gt(left_ext_consensus.length(), 0);
assert_gt(right_ext_consensus.length(), 0);
break;
}
size_t num_passed_seeds = 0;
if(left_ext_consensus.length() > 0 && right_ext_consensus.length()) {
consensus_ = reverse(left_ext_consensus) + consensus_;
consensus_ += right_ext_consensus;
#if 0
calc_edit_dist(s,
seeds_,
0,
remain,
left ? consensuses : empty_consensuses,
left ? empty_consensuses : consensuses,
allowed_seed_mm);
#endif
// update seeds
for(size_t i = 0; i < seeds_.size(); i++) {
SeedExt& seed = seeds_[i];
if(i < remain) {
if(seed.ed <= allowed_seed_mm) {
num_passed_seeds++;
seed.done = true;
seed.total_ed += seed.ed;
seed.pos.first -= left_ext_consensus.length();
seed.pos.second += right_ext_consensus.length();
seed.consensus_pos.first = 0;
seed.consensus_pos.second = consensus_.length();
if(seed.left_gaps.size() > 0 &&
seed.left_gaps.back().first >= seed.orig_pos.first - seed.pos.first) {
int gap_len = seed.left_gaps.back().second;
seed.pos.first += gap_len;
}
if(seed.right_gaps.size() > 0 &&
seed.right_gaps.back().first >= seed.pos.second - seed.orig_pos.second) {
int gap_len = seed.right_gaps.back().second;
seed.pos.second -= gap_len;
}
}
}
}
// move up "done" seeds
size_t j = 0;
for(size_t i = 0; i < remain; i++) {
if(!seeds_[i].done) continue;
assert_geq(i, j);
if(i > j) {
SeedExt temp = seeds_[j];
seeds_[j] = seeds_[i];
seeds_[i] = temp;
// Find next "undone" seed
j++;
while(j < i && seeds_[j].done) {
j++;
}
assert(j < remain && !seeds_[j].done);
} else {
j = i + 1;
}
}
}
remain = num_passed_seeds;
break;
} // while(remain >= rp.repeat_count)
#ifndef NDEBUG
// make sure seed positions are unique
EList<pair<TIndexOffU, TIndexOffU> > seed_poses;
for(size_t i = 0; i < seeds_.size(); i++) {
seed_poses.expand();
seed_poses.back() = seeds_[i].orig_pos;
}
seed_poses.sort();
for(size_t i = 0; i + 1 < seed_poses.size(); i++) {
if(seed_poses[i].first == seed_poses[i+1].first) {
assert_lt(seed_poses[i].second, seed_poses[i+1].second);
} else {
assert_lt(seed_poses[i].first, seed_poses[i+1].first);
}
}
for(size_t i = 0; i < seeds_.size(); i++) {
assert(seeds_[i].valid());
}
#endif
RB_Repeat::internal_update();
// check repeats within a repeat sequence
self_repeat_ = false;
for(size_t i = 0; i + 1 < seed_ranges_.size(); i++) {
const RB_AlleleCoord& range = seed_ranges_[i];
for(size_t j = i + 1; j < seed_ranges_.size(); j++) {
RB_AlleleCoord& range2 = seed_ranges_[j];
if(range.right <= range2.left)
break;
self_repeat_ = true;
}
}
}
template<typename TStr>
void RB_Repeat::getNextRepeat(const RepeatParameter& rp,
const TStr& s,
RB_Repeat& o)
{
o.reset();
size_t i = 0;
for(i = 0; i < seeds_.size() && seeds_[i].done; i++);
if(i >= seeds_.size())
return;
for(size_t j = i; j < seeds_.size(); j++) {
assert(!seeds_[j].done);
o.seeds_.expand();
o.seeds_.back() = seeds_[j];
}
seeds_.resizeExact(i);
internal_update();
assert_gt(o.seeds_.size(), 0);
o.consensus_ = getString(s, o.seeds_[0].orig_pos.first, o.seeds_[0].orig_pos.second - o.seeds_[0].orig_pos.first);
for(size_t j = 0; j < o.seeds_.size(); j++) {
o.seeds_[j].pos = o.seeds_[j].orig_pos;
o.seeds_[j].left_gaps.clear();
o.seeds_[j].right_gaps.clear();
o.seeds_[j].ed = o.seeds_[j].total_ed = 0;
}
}
void RB_Repeat::internal_update()
{
sort(seeds_.begin(), seeds_.end(), seedCmp);
seed_ranges_.resizeExact(seeds_.size());
for(size_t i = 0; i < seeds_.size(); i++) {
seed_ranges_[i].left = seeds_[i].pos.first;
seed_ranges_[i].right = seeds_[i].pos.second;
seed_ranges_[i].idx = i;
}
seed_ranges_.sort();
// check repeats within a repeat sequence
size_t remove_count = 0;
for(size_t i = 0; i + 1 < seed_ranges_.size(); i++) {
const RB_AlleleCoord& range = seed_ranges_[i];
if(range.left == numeric_limits<size_t>::max())
continue;
for(size_t j = i + 1; j < seed_ranges_.size(); j++) {
RB_AlleleCoord& range2 = seed_ranges_[j];
if(range2.left == numeric_limits<size_t>::max())
continue;
if(range.right <= range2.left)
break;
if(range.right >= range2.right) {
range2.left = numeric_limits<size_t>::max();
remove_count++;
}
}
}
if(remove_count <= 0)
return;
for(size_t i = 0; i < seed_ranges_.size(); i++) {
if(seed_ranges_[i].left == numeric_limits<size_t>::max())
seeds_[seed_ranges_[i].idx].reset();
}
sort(seeds_.begin(), seeds_.end(), seedCmp);
#ifndef NDEBUG
for(size_t i = 0; i < seeds_.size(); i++) {
if(i < seeds_.size() - remove_count) {
assert_lt(seeds_[i].pos.first, seeds_[i].pos.second);
} else {
assert_eq(seeds_[i].pos.first, seeds_[i].pos.second);
}
}
#endif
seeds_.resize(seeds_.size() - remove_count);
seed_ranges_.resize(seeds_.size());
for(size_t i = 0; i < seeds_.size(); i++) {
seed_ranges_[i].left = seeds_[i].pos.first;
seed_ranges_[i].right = seeds_[i].pos.second;
seed_ranges_[i].idx = i;
}
seed_ranges_.sort();
}
bool RB_Repeat::contain(TIndexOffU left, TIndexOffU right) const
{
size_t l = 0, r = seed_ranges_.size();
while(l < r) {
size_t m = (l + r) / 2;
const RB_AlleleCoord& coord = seed_ranges_[m];
if(right <= coord.left) {
r = m;
} else if(left >= coord.right) {
l = m + 1;
} else {
return coord.left <= left && right <= coord.right;
}
}
return false;
}
template<typename TStr>
void RB_Repeat::saveSeedExtension(const RepeatParameter& rp,
const TStr& s,
CoordHelper& coordHelper,
TIndexOffU grp_id,
ostream& fp,
size_t& total_repeat_seq_len,
size_t& total_allele_seq_len) const
{
// apply color, which is compatible with linux commands such as cat and less -r
#if 1
const string red = "\033[31m", reset = "\033[0m", resetline = "\033[4m", redline = "\033[4;31m";
#else
const string red = "", reset = "", resetline = "", redline = "";
#endif
bool show_exterior_seq = true;
const size_t max_show_seq_len = 700;
size_t total_count = 0;
for(size_t i = 0; i < seeds_.size(); i++) {
const SeedExt& seed = seeds_[i];
size_t ext_len = seed.getLength();
if(ext_len < rp.min_repeat_len) continue;
total_allele_seq_len += ext_len;
total_count++;
bool sense_strand = seed.pos.first < coordHelper.forward_length();
fp << setw(6) << repeat_id_ << " " << setw(5) << seeds_.size();
fp << " " << setw(4) << i;
fp << " " << setw(4) << ext_len;
fp << " " << (sense_strand ? '+' : '-');
fp << " " << setw(10) << seed.pos.first << " " << setw(10) << seed.pos.second;
fp << " " << setw(10) << seed.orig_pos.first << " " << setw(10) << seed.orig_pos.second;
string chr_name;
TIndexOffU pos_in_chr;
if(sense_strand) {
coordHelper.getGenomeCoord(seed.pos.first, chr_name, pos_in_chr);
} else {
coordHelper.getGenomeCoord(s.length() - seed.pos.first - (seed.pos.second - seed.pos.first), chr_name, pos_in_chr);
}
fp << " " << setw(5) << chr_name << ":" << setw(10) << std::left << pos_in_chr << std::right;
if(sense_strand) {
coordHelper.getGenomeCoord(seed.pos.second, chr_name, pos_in_chr);
} else {
coordHelper.getGenomeCoord(s.length() - seed.pos.second - (seed.pos.second - seed.pos.first), chr_name, pos_in_chr);
}
fp << " " << setw(5) << chr_name << ":" << setw(10) << std::left << pos_in_chr << std::right;
if(!seed.aligned) {
fp << endl;
continue;
}
string deststr = "";
seed.getExtendedSeedSequence(s, deststr);
// add exterior sequences
if(seed.consensus_pos.first > 0) {
TIndexOffU seq_pos, seq_len;
if(seed.pos.first >= seed.consensus_pos.first) {
seq_pos = seed.pos.first - seed.consensus_pos.first;
seq_len = seed.consensus_pos.first;
} else {
seq_pos = 0;
seq_len = seed.pos.first;
}
deststr = getString(s, seq_pos, seq_len) + deststr;
if(seq_len < seed.consensus_pos.first) {
deststr = string(seed.consensus_pos.first - seq_len, 'N') + deststr;
}
}
if(seed.consensus_pos.second < consensus_.length()) {
deststr += getString(s, seed.pos.second, consensus_.length() - seed.consensus_pos.second);
}
assert_eq(consensus_.length(), deststr.length());
fp << " ";
// print sequence w.r.t. the current group
for(size_t j = 0; j < min(consensus_.length(), max_show_seq_len); j++) {
bool outside = j < seed.consensus_pos.first || j >= seed.consensus_pos.second;
bool different = (consensus_[j] != deststr[j]);
if(outside) {
if(show_exterior_seq) {
if(different) {
fp << redline;
fp << deststr[j];
fp << reset;
} else {
fp << resetline;
fp << deststr[j];
fp << reset;
}
} else {
if(j < seed.consensus_pos.first)
fp << " ";
}
} else {
if(different) fp << red;
fp << deststr[j];
if(different) fp << reset;
}
}
#if 0
fp << "\t";
for(size_t ei = 0; ei < seed.edits.size(); ei++) {
const Edit& edit = seed.edits[ei];
if(ei > 0) fp << ",";
fp << edit;
if (edit.snpID != std::numeric_limits<uint32_t>::max()) {
fp << "@" << edit.snpID;
}
}
#endif
fp << endl;
}
if(total_count > 0) fp << setw(5) << total_count << endl << endl;
total_repeat_seq_len += consensus_.length();
}
template<typename TStr>
size_t calc_edit_dist(const TStr& s,
const SeedExt& seed,
const SeedExt& seed2,
size_t left_ext,
size_t right_ext,
size_t max_ed)
{
if(seed.bound.first + left_ext > seed.pos.first ||
seed.pos.second + right_ext > seed.bound.second ||
seed2.bound.first + left_ext > seed2.pos.first ||
seed2.pos.second + right_ext > seed2.bound.second)
return max_ed + 1;
size_t ed = 0;
for(size_t i = 0; i < left_ext; i++) {
int ch = getSequenceBase(s, seed.pos.first - i - 1);
assert_range(0, 3, ch);
int ch2 = getSequenceBase(s, seed2.pos.first - i - 1);
assert_range(0, 3, ch2);
if(ch != ch2) ed++;
if(ed > max_ed)
return ed;
}
for(size_t i = 0; i < right_ext; i++) {
int ch = getSequenceBase(s, seed.pos.second + i);
assert_range(0, 3, ch);
int ch2 = getSequenceBase(s, seed2.pos.second + i);
assert_range(0, 3, ch2);
if(ch != ch2) ed++;
if(ed > max_ed)
return ed;
}
return ed;
}
size_t extract_kmer(const string& seq, size_t offset, size_t k = 5)
{
assert_leq(offset + k, seq.length());
size_t kmer = 0;
for(size_t i = 0; i < k; i++) {
kmer = (kmer << 2) | asc2dna[seq[offset + i]];
}
return kmer;
}
size_t next_kmer(size_t kmer, char base, size_t k = 5)
{
kmer &= ((1 << ((k-1)*2))) - 1;
kmer = (kmer << 2) | asc2dna[base];
return kmer;
}
void build_kmer_table(const string& consensus, EList<pair<size_t, size_t> >& kmer_table, size_t k = 5)
{
kmer_table.clear();
if(consensus.length() < k)
return;
size_t kmer = 0;
for(size_t i = 0; i + k <= consensus.length(); i++) {
if(i == 0) {
kmer = extract_kmer(consensus, i, k);
} else {
kmer = next_kmer(kmer, consensus[i+k-1], k);
}
kmer_table.expand();
kmer_table.back().first = kmer;
kmer_table.back().second = i;
}
kmer_table.sort();
}
void find_gap_pos(const string& s,
const string& s2,
EList<size_t>& ed,
EList<size_t>& ed2,
bool del,
size_t gap_len,
size_t max_mm,
size_t& gap_pos,
size_t& mm,
bool debug = false)
{
assert_eq(s.length(), s2.length());
size_t seq_len = s.length();
ed.resizeExact(seq_len); ed.fill(max_mm + 1);
ed2.resizeExact(seq_len); ed2.fill(max_mm + 1);
// from left to right
for(size_t i = 0; i < seq_len; i++) {
size_t add = (s[i] == s2[i] ? 0 : 1);
if(i == 0) ed[i] = add;
else ed[i] = ed[i-1] + add;
if(ed[i] >= max_mm + 1) break;
}
// from right to left
size_t s_sub = del ? 0 : gap_len;
size_t s2_sub = del ? gap_len : 0;
for(int i = seq_len - 1; i >= gap_len; i--) {
size_t add = (s[i - s_sub] == s2[i - s2_sub] ? 0 : 1);
if(i == seq_len - 1) ed2[i] = add;
else ed2[i] = ed2[i+1] + add;
if(ed2[i] > max_mm) break;
}
if(debug) {
cout << s << endl << s2 << endl;
for(size_t i = 0; i < ed.size(); i++) {
cout << (ed[i] % 10);
}
cout << endl;
for(size_t i = 0; i < ed2.size(); i++) {
cout << (ed2[i] % 10);
}
cout << endl;
}
size_t min_mm = ed2[gap_len];
int min_mm_i = -1;
assert_eq(ed.size(), ed2.size());
for(size_t i = 0; i + gap_len + 1 < ed.size(); i++) {
if(ed[i] > max_mm) break;
size_t cur_mm = ed[i] + ed2[i + gap_len + 1];
if(cur_mm < min_mm) {
min_mm = cur_mm;
min_mm_i = i;
}
}
mm = min_mm;
if(mm > max_mm)
return;
gap_pos = (size_t)(min_mm_i + 1);
}
void align_with_one_gap(const string& s,
const EList<pair<size_t, size_t> >& s_kmer_table,
const string& s2,
EList<size_t>& counts,
size_t max_gap,
size_t max_mm,
size_t& mm,
size_t& gap_pos,
int& gap_len,
size_t k = 5)
{
mm = max_mm + 1;
assert_eq(s.length(), s2.length());
counts.resizeExact(max_gap * 2 + 1);
counts.fillZero();
size_t max_count = 0, max_count_i = 0;
for(size_t i = 0; i + k <= s2.length(); i += 1) {
pair<size_t, size_t> kmer(0, 0);
kmer.first = extract_kmer(s2, i, k);
size_t lb = s_kmer_table.bsearchLoBound(kmer);
while(lb < s_kmer_table.size() && kmer.first == s_kmer_table[lb].first) {
int gap = (int)s_kmer_table[lb].second - (int)i;
if(gap != 0 && abs(gap) < max_gap) {
size_t gap_i = (size_t)(gap + max_gap);
counts[gap_i]++;
if(counts[gap_i] > max_count) {
max_count = counts[gap_i];
max_count_i = gap_i;
}
}
lb++;
}
}
if(max_count <= 0)
return;
assert_lt(max_count_i, counts.size());
int gap = (int)max_count_i - (int)max_gap;
assert_neq(gap, 0);
size_t abs_gap = abs(gap);
bool del = gap > 0;
EList<size_t> ed, ed2;
find_gap_pos(s,
s2,
ed,
ed2,
del,
abs_gap,
max_mm,
gap_pos,
mm);
if(mm > max_mm)
return;
gap_len = del ? abs_gap : -abs_gap;
#ifndef NDEBUG
assert_leq(mm, max_mm);
string ds = s, ds2 = s2;
size_t debug_mm = 0;
if(del) ds.erase(gap_pos, abs_gap);
else ds2.erase(gap_pos, abs_gap);
for(size_t i = 0; i < min(ds.length(), ds2.length()); i++) {
if(ds[i] != ds2[i])
debug_mm++;
}
assert_eq(debug_mm, mm);
#endif
}
template<typename TStr>
void calc_edit_dist(const TStr& s,
EList<SeedExt>& seeds,
size_t sb,
size_t se,
const EList<string>& left_consensuses,
const EList<string>& right_consensuses,
uint32_t max_ed)
{
string left_consensus = left_consensuses[max_ed];
string right_consensus = right_consensuses[max_ed];
EList<pair<size_t, size_t> > left_kmer_table, right_kmer_table;
build_kmer_table(left_consensus, left_kmer_table);
build_kmer_table(right_consensus, right_kmer_table);
string left_seq, right_seq;
const size_t max_gap = 10;
EList<size_t> counts;
size_t left_ext = left_consensus.length();
size_t right_ext = right_consensus.length();
for(size_t i = sb; i < se; i++) {
SeedExt& seed = seeds[i];
if(seed.bound.first + left_ext > seed.pos.first ||
seed.pos.second + right_ext > seed.bound.second) {
seed.ed = max_ed + 1;
continue;
}
size_t left_ed = 0;
if(left_ext > 0) {
getString(s, seed.pos.first - left_ext, left_ext, left_seq);
reverse(left_seq.begin(), left_seq.end());
for(size_t j = 0; j < left_ext; j++) {
if(left_seq[j] != left_consensus[j]) left_ed++;
if(left_ed <= max_ed && j < left_consensuses[left_ed].length()) {
seed.curr_ext_len = j + 1;
}
}
if(left_ed > max_ed) {
size_t gap_pos = 0;
int gap_len = 0;
align_with_one_gap(left_consensus,
left_kmer_table,
left_seq,
counts,
min(left_ed - max_ed, max_gap),
max_ed,
left_ed,
gap_pos,
gap_len);
if(left_ed <= max_ed) {
seed.left_gaps.expand();
seed.left_gaps.back().first = seed.getLeftExtLength() + gap_pos;
seed.left_gaps.back().second = gap_len;
}
}
} else {
left_seq.clear();
}
size_t right_ed = 0;
if(right_ext > 0) {
getString(s, seed.pos.second, right_ext, right_seq);
for(size_t j = 0; j < right_ext; j++) {
if(right_seq[j] != right_consensus[j]) right_ed++;
if(right_ed <= max_ed && j < right_consensuses[right_ed].length()) {
seed.curr_ext_len = j + 1;
}
}
if(right_ed > max_ed) {
size_t gap_pos = 0;
int gap_len = 0;
align_with_one_gap(right_consensus,
right_kmer_table,
right_seq,
counts,
min(right_ed - max_ed, max_gap),
max_ed,
right_ed,
gap_pos,
gap_len);
if(right_ed <= max_ed) {
seed.right_gaps.expand();
seed.right_gaps.back().first = seed.getRightExtLength() + gap_pos;
seed.right_gaps.back().second = gap_len;
}
}
} else {
right_seq.clear();
}
seed.ed = left_ed + right_ed;
}
}
template<typename TStr>
void RB_Repeat::get_consensus_seq(const TStr& s,
EList<SeedExt>& seeds,
size_t sb,
size_t se,
size_t min_left_ext,
size_t min_right_ext,
size_t max_ed,
const RepeatParameter& rp,
EList<size_t>& ed_seed_nums,
EList<string>* left_consensuses,
EList<string>* right_consensuses) const
{
assert_lt(sb, se);
assert_geq(se - sb, rp.seed_count);
assert_leq(se, seeds.size());
if(left_consensuses != NULL) {
left_consensuses->clear(); left_consensuses->resize(max_ed + 1);
for(size_t i = 0; i < max_ed + 1; i++) {
(*left_consensuses)[i].clear();
}
}
if(right_consensuses != NULL) {
right_consensuses->clear(); right_consensuses->resize(max_ed + 1);
for(size_t i = 0; i < max_ed + 1; i++) {
(*right_consensuses)[i].clear();
}
}
assert_eq(min_left_ext, min_right_ext);
EList<string> seqs;
seqs.reserveExact(seeds.size());
size_t max_i = 0, max_count = 0;
while(min_left_ext > 0) {
seqs.clear();
max_i = 0;
max_count = 0;
for(size_t i = sb; i < se; i++) {
const SeedExt& seed = seeds[i];
if(seeds[i].bound.first + min_left_ext > seed.pos.first ||
seed.pos.second + min_right_ext > seed.bound.second)
continue;
seqs.expand();
if(min_left_ext > 0) {
getString(s, seed.pos.first - min_left_ext, min_left_ext, seqs.back());
}
if(min_right_ext > 0) {
seqs.back() += getString(s, seed.pos.second, min_right_ext);
}
}
seqs.sort();
for(size_t i = 0; i + max_count < seqs.size();) {
size_t count = 1;
for(size_t j = i + 1; j < seqs.size(); j++) {
if(seqs[i] == seqs[j]) count++;
else break;
}
if(count >= max_count) {
max_count = count;
max_i = i;
}
i = i + count;
}
if(max_count >= rp.seed_count)
break;
min_left_ext--;
min_right_ext--;
}
// update ed
// extend consensus string
ed_seed_nums.resize(max_ed + 1); ed_seed_nums.fillZero();
EList<size_t> next_ed_seed_nums; next_ed_seed_nums.resize(max_ed + 1);
if(max_count < rp.seed_count)
return;
for(size_t i = sb; i < se; i++) seeds[i].ed = 0;
size_t seed_ext_len = 0;
while(seed_ext_len < max(min_left_ext, min_right_ext)) {
// select base to be used for extension
uint8_t left_ext_base = 0;
if(seed_ext_len < min_left_ext) left_ext_base = seqs[max_i][min_left_ext - seed_ext_len - 1];
uint8_t right_ext_base = 0;
if(seed_ext_len < min_right_ext) right_ext_base = seqs[max_i][min_left_ext + seed_ext_len];
// estimate extended ed
next_ed_seed_nums.fillZero();
for(size_t i = sb; i < se; i++) {
uint32_t est_ed = seeds[i].ed;
if(seed_ext_len < min_left_ext) {
if(seeds[i].bound.first + seed_ext_len + 1 <= seeds[i].pos.first) {
TIndexOffU left_pos = seeds[i].pos.first - seed_ext_len - 1;
int ch = getSequenceBase(s, left_pos);
assert_range(0, 3, ch);
if ("ACGT"[ch] != left_ext_base) {
est_ed++;
}
} else {
est_ed = max_ed + 1;
}
}
if(seed_ext_len < min_right_ext) {
TIndexOffU right_pos = seeds[i].pos.second + seed_ext_len;
if(right_pos >= seeds[i].bound.second) {
est_ed = max_ed + 1;
} else {
int ch = getSequenceBase(s, right_pos);
assert_range(0, 3, ch);
if ("ACGT"[ch] != right_ext_base) {
est_ed++;
}
}
}
if(est_ed <= max_ed) {
next_ed_seed_nums[est_ed]++;
}
}
for(size_t i = 1; i < next_ed_seed_nums.size(); i++) next_ed_seed_nums[i] += next_ed_seed_nums[i-1];
if(next_ed_seed_nums[max_ed] < rp.repeat_count) {
break;
}
for(size_t i = sb; i < se; i++) {
if(seed_ext_len < min_left_ext) {
if(seeds[i].bound.first + seed_ext_len + 1 <= seeds[i].pos.first) {
TIndexOffU left_pos = seeds[i].pos.first - seed_ext_len - 1;
int ch = getSequenceBase(s, left_pos);
assert_range(0, 3, ch);
if ("ACGT"[ch] != left_ext_base) {
seeds[i].ed++;
}
} else {
seeds[i].ed = max_ed + 1;
}
}
if(seed_ext_len < min_right_ext) {
TIndexOffU right_pos = seeds[i].pos.second + seed_ext_len;
if(right_pos >= seeds[i].bound.second) {
seeds[i].ed = max_ed + 1;
} else {
int ch = getSequenceBase(s, right_pos);
assert_range(0, 3, ch);
if ("ACGT"[ch] != right_ext_base) {
seeds[i].ed++;
}
}
}
}
for(size_t i = 0; i < next_ed_seed_nums.size(); i++) {
if(next_ed_seed_nums[i] < rp.repeat_count)
continue;
ed_seed_nums[i] = next_ed_seed_nums[i];
if(seed_ext_len < min_left_ext) {
assert(left_consensuses != NULL);
(*left_consensuses)[i] += left_ext_base;
}
if(seed_ext_len < min_right_ext) {
assert(right_consensuses != NULL);
(*right_consensuses)[i] += right_ext_base;
}
}
seed_ext_len++;
}
}
template<typename TStr>
void RB_RepeatExt::get_consensus_seq(const TStr& s,
EList<SeedExt>& seeds,
size_t sb,
size_t se,
size_t min_left_ext,
size_t min_right_ext,
size_t max_ed,
const RepeatParameter& rp,
EList<size_t>& ed_seed_nums,
EList<string>* left_consensuses,
EList<string>* right_consensuses) const
{
assert_lt(sb, se);
assert_leq(se, seeds.size());
if(left_consensuses != NULL) {
left_consensuses->clear(); left_consensuses->resize(max_ed + 1);
for(size_t i = 0; i < max_ed + 1; i++) {
(*left_consensuses)[i].clear();
}
}
if(right_consensuses != NULL) {
right_consensuses->clear(); right_consensuses->resize(max_ed + 1);
for(size_t i = 0; i < max_ed + 1; i++) {
(*right_consensuses)[i].clear();
}
}
// cluster sequences
EList<size_t> belongto;
belongto.resizeExact(se - sb);
for(size_t i = 0; i < belongto.size(); i++) belongto[i] = i;
for(size_t i = 0; i + 1 < belongto.size(); i++) {
for(size_t j = i + 1; j < belongto.size(); j++) {
if(belongto[j] != j) continue;
size_t ed = calc_edit_dist(s,
seeds[sb + i],
seeds[sb + j],
min_left_ext,
min_right_ext,
max_ed + 1);
if(ed <= max_ed + 1) {
belongto[j] = belongto[i];
}
}
}
// find the maximum group that has the most sequences
EList<size_t> vote;
vote.resizeExact(seeds.size());
vote.fillZero();
size_t max_group = 0;
for(size_t i = 0; i < belongto.size(); i++) {
size_t cur_group = belongto[i];
assert_lt(cur_group, vote.size());
vote[cur_group]++;
if(cur_group != max_group && vote[cur_group] > vote[max_group]) {
max_group = cur_group;
}
}
// reuse vote to store seeds
EList<size_t>& consensus_group = vote;
consensus_group.clear();
for(size_t i = 0; i < belongto.size(); i++) {
if(belongto[i] == max_group)
consensus_group.push_back(i);
}
// update ed
// extend consensus string
ed_seed_nums.resize(max_ed + 1); ed_seed_nums.fillZero();
EList<size_t> next_ed_seed_nums; next_ed_seed_nums.resize(max_ed + 1);
for(size_t i = sb; i < se; i++) seeds[i].ed = 0;
size_t seed_ext_len = 0;
while(seed_ext_len < max(min_left_ext, min_right_ext)) {
// count base
size_t l_count[4] = {0, };
size_t r_count[4] = {0, };
for(size_t i = 0; i < consensus_group.size(); i++) {
size_t si = sb + consensus_group[i];
int ch;
if(seed_ext_len < min_left_ext) {
if(seeds[i].bound.first + seed_ext_len + 1 <= seeds[i].pos.first) {
ch = getSequenceBase(s, seeds[si].pos.first - seed_ext_len - 1);
assert_range(0, 3, ch);
l_count[ch]++;
}
}
if(seed_ext_len < min_right_ext) {
if(seeds[i].bound.second + seed_ext_len) {
ch = getSequenceBase(s, seeds[si].pos.second + seed_ext_len);
assert_range(0, 3, ch);
r_count[ch]++;
}
}
}
// select base to be used for extension
uint8_t left_ext_base = 0;
if(seed_ext_len < min_left_ext) left_ext_base = max_index(l_count);
uint8_t right_ext_base = 0;
if(seed_ext_len < min_right_ext) right_ext_base = max_index(r_count);
// estimate extended ed
next_ed_seed_nums.fillZero();
for(size_t i = sb; i < se; i++) {
uint32_t est_ed = seeds[i].ed;
if(seed_ext_len < min_left_ext) {
if(seeds[i].bound.first + seed_ext_len + 1 <= seeds[i].pos.first) {
TIndexOffU left_pos = seeds[i].pos.first - seed_ext_len - 1;
int ch = getSequenceBase(s, left_pos);
assert_range(0, 3, ch);
if (ch != left_ext_base) {
est_ed++;
}
} else {
est_ed = max_ed + 1;
}
}
if(seed_ext_len < min_right_ext) {
TIndexOffU right_pos = seeds[i].pos.second + seed_ext_len;
if(right_pos >= seeds[i].bound.second) {
est_ed = max_ed + 1;
} else {
int ch = getSequenceBase(s, right_pos);
assert_range(0, 3, ch);
if (ch != right_ext_base) {
est_ed++;
}
}
}
if(est_ed <= max_ed) {
next_ed_seed_nums[est_ed]++;
}
}
for(size_t i = 1; i < next_ed_seed_nums.size(); i++) next_ed_seed_nums[i] += next_ed_seed_nums[i-1];
if(next_ed_seed_nums[max_ed] < rp.repeat_count) {
break;
}
for(size_t i = sb; i < se; i++) {
if(seed_ext_len < min_left_ext) {
if(seeds[i].bound.first + seed_ext_len + 1 <= seeds[i].pos.first) {
TIndexOffU left_pos = seeds[i].pos.first - seed_ext_len - 1;
int ch = getSequenceBase(s, left_pos);
assert_range(0, 3, ch);
if (ch != left_ext_base) {
seeds[i].ed++;
}
} else {
seeds[i].ed = max_ed + 1;
}
}
if(seed_ext_len < min_right_ext) {
TIndexOffU right_pos = seeds[i].pos.second + seed_ext_len;
if(right_pos >= seeds[i].bound.second) {
seeds[i].ed = max_ed + 1;
} else {
int ch = getSequenceBase(s, right_pos);
assert_range(0, 3, ch);
if (ch != right_ext_base) {
seeds[i].ed++;
}
}
}
}
for(size_t i = 0; i < next_ed_seed_nums.size(); i++) {
if(next_ed_seed_nums[i] < rp.repeat_count)
continue;
ed_seed_nums[i] = next_ed_seed_nums[i];
if(seed_ext_len < min_left_ext) {
assert(left_consensuses != NULL);
(*left_consensuses)[i] += (char)("ACGT"[left_ext_base]);
}
if(seed_ext_len < min_right_ext) {
assert(right_consensuses != NULL);
(*right_consensuses)[i] += (char)("ACGT"[right_ext_base]);
}
}
seed_ext_len++;
}
}
EList<size_t> RB_Repeat::ca_ed_;
EList<size_t> RB_Repeat::ca_ed2_;
string RB_Repeat::ca_s_;
string RB_Repeat::ca_s2_;
size_t RB_Repeat::seed_merge_tried = 0;
size_t RB_Repeat::seed_merged = 0;
template<typename TStr>
void RB_RepeatExt::extendConsensus(const RepeatParameter& rp,
const TStr& s)
{
size_t remain = seeds_.size();
EList<string> consensuses, empty_consensuses;
EList<size_t> ed_seed_nums;
bool left = true;
const TIndexOffU default_max_ext_len = 100;
const TIndexOffU seed_mm = 4;
string ext_consensus = "";
empty_consensuses.resize(seed_mm + 1);
empty_consensuses.fill("");
while(remain >= rp.repeat_count) {
for(size_t i = 0; i < remain; i++) {
seeds_[i].done = false;
seeds_[i].curr_ext_len = 0;
}
// extend seeds in left or right direction
TIndexOffU max_ext_len = min(default_max_ext_len, (TIndexOffU)(rp.max_repeat_len - consensus_.length()));
get_consensus_seq(s,
seeds_,
0, // seed begin
remain, // seed end
left ? max_ext_len : 0,
left ? 0 : max_ext_len,
seed_mm,
rp,
ed_seed_nums,
left ? &consensuses : NULL,
left ? NULL : &consensuses);
size_t allowed_seed_mm = 0;
ext_consensus.clear();
for(int i = (int)seed_mm; i >= 0; i--) {
size_t extlen = (ed_seed_nums[i] < rp.repeat_count ? 0 : consensuses[i].length());
if(extlen <= 0 || extlen < max_ext_len * i / seed_mm)
continue;
if(i > 0 && consensuses[i].length() <= consensuses[i-1].length() + 5)
continue;
ext_consensus = consensuses[i];
allowed_seed_mm = (size_t)i;
assert(ext_consensus.length() > 0);
break;
}
size_t num_passed_seeds = 0;
if(ext_consensus.length() > 0) {
if(left) consensus_ = reverse(ext_consensus) + consensus_;
else consensus_ += ext_consensus;
calc_edit_dist(s,
seeds_,
0,
remain,
left ? consensuses : empty_consensuses,
left ? empty_consensuses : consensuses,
allowed_seed_mm);
// update seeds
for(size_t i = 0; i < seeds_.size(); i++) {
SeedExt& seed = seeds_[i];
if(i < remain) {
if(seed.ed <= allowed_seed_mm) {
num_passed_seeds++;
seed.done = true;
seed.total_ed += seed.ed;
if(left) {
if(seed.left_gaps.size() > 0 &&
seed.left_gaps.back().first >= seed.orig_pos.first - seed.pos.first) {
int gap_len = seed.left_gaps.back().second;
seed.pos.first += gap_len;
}
seed.pos.first -= ext_consensus.length();
seed.consensus_pos.first = 0;
seed.consensus_pos.second = consensus_.length();
} else {
if(seed.right_gaps.size() > 0 &&
seed.right_gaps.back().first >= seed.pos.second - seed.orig_pos.second) {
int gap_len = seed.right_gaps.back().second;
seed.pos.second -= gap_len;
}
seed.pos.second += ext_consensus.length();
seed.consensus_pos.second = consensus_.length();
}
} else {
if(left) {
assert_leq(seed.curr_ext_len, ext_consensus.length());
TIndexOffU adjust = ext_consensus.length() - seed.curr_ext_len;
seed.consensus_pos.first += adjust;
seed.consensus_pos.second += ext_consensus.length();
assert_leq(seed.curr_ext_len, seed.pos.first);
seed.pos.first -= seed.curr_ext_len;
} else {
assert_leq(seed.curr_ext_len, ext_consensus.length());
seed.consensus_pos.second += seed.curr_ext_len;
seed.pos.second += seed.curr_ext_len;
}
}
} else {
if(left) {
seed.consensus_pos.first += ext_consensus.length();
seed.consensus_pos.second += ext_consensus.length();
}
}
}
// move up "done" seeds
size_t j = 0;
for(size_t i = 0; i < remain; i++) {
if(!seeds_[i].done) continue;
assert_geq(i, j);
if(i > j) {
SeedExt temp = seeds_[j];
seeds_[j] = seeds_[i];
seeds_[i] = temp;
// Find next "undone" seed
j++;
while(j < i && seeds_[j].done) {
j++;
}
assert(j < remain && !seeds_[j].done);
} else {
j = i + 1;
}
}
}
remain = num_passed_seeds;
if(remain < rp.repeat_count) {
if(left) {
left = false;
remain = seeds_.size();
}
}
} // while(remain >= rp.repeat_count)
#ifndef NDEBUG
// make sure seed positions are unique
EList<pair<TIndexOffU, TIndexOffU> > seed_poses;
for(size_t i = 0; i < seeds_.size(); i++) {
seed_poses.expand();
seed_poses.back() = seeds_[i].orig_pos;
}
seed_poses.sort();
for(size_t i = 0; i + 1 < seed_poses.size(); i++) {
if(seed_poses[i].first == seed_poses[i+1].first) {
assert_lt(seed_poses[i].second, seed_poses[i+1].second);
} else {
assert_lt(seed_poses[i].first, seed_poses[i+1].first);
}
}
for(size_t i = 0; i < seeds_.size(); i++) {
assert(seeds_[i].valid());
}
#endif
internal_update();
// check repeats within a repeat sequence
self_repeat_ = false;
for(size_t i = 0; i + 1 < seed_ranges_.size(); i++) {
const RB_AlleleCoord& range = seed_ranges_[i];
for(size_t j = i + 1; j < seed_ranges_.size(); j++) {
RB_AlleleCoord& range2 = seed_ranges_[j];
if(range.right <= range2.left)
break;
self_repeat_ = true;
}
}
}
bool RB_Repeat::overlap(const RB_Repeat& o,
bool& contain,
bool& left,
size_t& seed_i,
size_t& seed_j,
bool debug) const
{
contain = left = false;
seed_i = seed_j = 0;
size_t p = 0, p2 = 0;
while(p < seed_ranges_.size() && p2 < o.seed_ranges_.size()) {
const RB_AlleleCoord& range = seed_ranges_[p];
const SeedExt& seed = seeds_[range.idx];
Range range_ = seed.getExtendedRange(consensus_.length());
RB_AlleleCoord ex_range(range_.first, range_.second, p);
bool representative = (float)range.len() >= consensus_.length() * 0.95f;
const RB_AlleleCoord& range2 = o.seed_ranges_[p2];
const SeedExt& seed2 = o.seeds_[range2.idx];
Range range2_ = seed2.getExtendedRange(o.consensus().length());
RB_AlleleCoord ex_range2(range2_.first, range2_.second, p2);
bool representative2 = (float)range2.len() >= o.consensus_.length() * 0.95f;
seed_i = range.idx;
seed_j = range2.idx;
const size_t relax = 5;
if(representative && representative2) {
if(ex_range.overlap_len(ex_range2) > 0) {
if(ex_range.contain(ex_range2, relax)) {
contain = true;
left = true;
} else if(ex_range2.contain(ex_range, relax)) {
contain = true;
left = false;
} else {
left = ex_range.left <= ex_range2.left;
}
return true;
}
} else if(representative) {
if(range2.contain(ex_range, relax)) {
contain = true;
left = false;
return true;
}
} else if(representative2) {
if(range.contain(ex_range2, relax)) {
contain = true;
left = true;
return true;
}
}
if(range.right <= range2.right) p++;
if(range2.right <= range.right) p2++;
}
return false;
}
void RB_Repeat::showInfo(const RepeatParameter& rp,
CoordHelper& coordHelper) const
{
cerr << "\trepeat id: " << repeat_id_ << endl;
cerr << "\tnumber of alleles: " << seeds_.size() << endl;
cerr << "\tconsensus length: " << consensus_.length() << endl;
cerr << consensus_ << endl;
for(size_t i = 0; i < seeds_.size(); i++) {
const SeedExt& seed = seeds_[i];
size_t ext_len = seed.getLength();
if(ext_len < rp.min_repeat_len) continue;
bool sense_strand = seed.pos.first < coordHelper.forward_length();
cerr << "\t\t" << setw(4) << i << " " << setw(4) << ext_len;
cerr << " " << (sense_strand ? '+' : '-');
cerr << " " << setw(10) << seed.pos.first << " " << setw(10) << seed.pos.second;
cerr << " " << setw(4) << seed.consensus_pos.first << " " << setw(4) << seed.consensus_pos.second;
cerr << " " << (seed.aligned ? "aligned" : "unaligned");
cerr << endl;
}
}
template<typename TStr>
void RB_Repeat::generateSNPs(const RepeatParameter& rp, const TStr& s, TIndexOffU grp_id) {
EList<SeedExt>& seeds = this->seeds();
for(size_t i = 0; i < seeds.size(); i++) {
SeedExt& seed = seeds[i];
if(!seed.aligned) {
continue;
}
if(seed.getLength() < rp.min_repeat_len) {
continue;
}
assert_eq(seed.getLength(), seed.pos.second - seed.pos.first);
seed.generateSNPs(s, consensus(), snps());
}
if(snps().size() > 0) {
sort(snps_.begin(), snps().end(), SeedSNP::cmpSeedSNPByPos);
}
}
void RB_Repeat::saveSNPs(ofstream &fp, TIndexOffU grp_id, TIndexOffU& snp_id_base)
{
string chr_name = "rep" + to_string(repeat_id_);
for(size_t j = 0; j < snps_.size(); j++) {
SeedSNP& snp = *snps_[j];
// assign SNPid
snp.id = (snp_id_base++);
fp << "rps" << snp.id;
fp << "\t";
if(snp.type == EDIT_TYPE_MM) {
fp << "single";
} else if(snp.type == EDIT_TYPE_READ_GAP) {
fp << "deletion";
} else if(snp.type == EDIT_TYPE_REF_GAP) {
fp << "insertion";
} else {
assert(false);
}
fp << "\t";
fp << chr_name;
fp << "\t";
fp << snp.pos;
fp << "\t";
if(snp.type == EDIT_TYPE_MM || snp.type == EDIT_TYPE_REF_GAP) {
fp << snp.base;
} else if(snp.type == EDIT_TYPE_READ_GAP) {
fp << snp.len;
} else {
assert(false);
}
fp << endl;
}
}
float RB_RepeatExt::mergeable(const RB_Repeat& o) const
{
const EList<RB_AlleleCoord>& ranges = seed_ranges();
const EList<RB_AlleleCoord>& ranges2 = o.seed_ranges();
size_t num_overlap_bp = 0;
size_t p = 0, p2 = 0;
while(p < ranges.size() && p2 < ranges2.size()) {
const RB_AlleleCoord& range = ranges[p];
const RB_AlleleCoord& range2 = ranges2[p2];
TIndexOffU overlap = range.overlap_len(range2);
num_overlap_bp += overlap;
if(range.right <= range2.right) p++;
else p2++;
}
size_t num_total_bp = 0, num_total_bp2 = 0;
for(size_t r = 0; r < ranges.size(); r++) num_total_bp += (ranges[r].right - ranges[r].left);
for(size_t r = 0; r < ranges2.size(); r++) num_total_bp2 += (ranges2[r].right - ranges2[r].left);
float portion = float(num_overlap_bp) / float(min(num_total_bp, num_total_bp2));
return portion;
}
inline void get_next_range(const EList<int>& offsets, size_t i, Range& r, float& avg)
{
r.first = i;
for(; r.first < offsets.size() && offsets[r.first] < 0; r.first++);
r.second = r.first + 1;
avg = 0.0f;
if(r.first < offsets.size()) {
float diff = (float)offsets[r.first] - (float)r.first;
avg += diff;
}
for(; r.second < offsets.size() &&
offsets[r.second] >= 0 &&
(r.second == 0 || offsets[r.second] >= offsets[r.second - 1]);
r.second++) {
float diff = (float)offsets[r.second] - (float)r.second;
avg += diff;
}
avg /= (float)(r.second - r.first);
return;
}
template<typename TStr>
bool RB_RepeatExt::align(const RepeatParameter& rp,
const TStr& ref,
const string& s,
const EList<pair<size_t, size_t> >& s_kmer_table,
const string& s2,
EList<int>& offsets,
size_t k,
SeedExt& seed,
int consensus_approx_left,
int consensus_approx_right,
size_t left,
size_t right,
bool debug)
{
RB_Repeat::seed_merge_tried++;
seed.reset();
seed.pos.first = left;
seed.pos.second = right;
seed.orig_pos.first = seed.pos.first;
seed.orig_pos.second = seed.orig_pos.first;
seed.consensus_pos.first = 0;
seed.consensus_pos.second = consensus_.length();
seed.aligned = false;
{
const int query_len = right - left;
const int consensus_len = consensus_approx_right - consensus_approx_left;
const int abs_gap_len = max<int>(abs(consensus_len - query_len), 5);
offsets.resize(s2.length());
offsets.fill(-1);
pair<size_t, size_t> kmer(0, 0);
for(size_t i = 0; i + k <= s2.length(); i++) {
if(i == 0) {
kmer.first = extract_kmer(s2, i, k);
} else {
kmer.first = next_kmer(kmer.first, s2[i+k-1], k);
}
size_t lb = s_kmer_table.bsearchLoBound(kmer);
while(lb < s_kmer_table.size() && kmer.first == s_kmer_table[lb].first) {
int expected = (int)i + consensus_approx_left;
int real = (int)s_kmer_table[lb].second;
int abs_diff = abs(expected - real);
if(abs_diff <= abs_gap_len * 2 || debug) {
if(offsets[i] == -1) {
offsets[i] = (int)s_kmer_table[lb].second;
} else if(offsets[i] >= 0) {
offsets[i] = -2;
} else {
assert_lt(offsets[i], -1);
offsets[i] -= 1;
}
}
lb++;
}
if(offsets[i] > 0 && i + k == s2.length()) {
for(size_t j = i + 1; j < s2.length(); j++) {
offsets[j] = offsets[j-1] + 1;
}
}
}
}
if(debug) {
cerr << "initial offsets" << endl;
for(size_t j = 0; j < offsets.size(); j++) {
cout << j << ": " << offsets[j] << " " << s2[j] << ": " << (offsets[j] < 0 ? ' ' : s[offsets[j]]) << endl;
}
}
// remove inconsistent positions
Range range; float range_avg;
get_next_range(offsets, 0, range, range_avg);
while(range.second < offsets.size()) {
Range range2; float range_avg2;
get_next_range(offsets, range.second, range2, range_avg2);
if(range2.first >= offsets.size())
break;
assert_leq(range.second, range2.first);
float abs_diff = abs(range_avg - range_avg2);
if(offsets[range.second - 1] > offsets[range2.first] ||
(abs_diff > 10.0f && (range.second - range.first < 5 || range2.second - range2.first < 5)) ||
abs_diff > 20.0f) {
if(range.second - range.first < range2.second - range2.first) {
for(size_t i = range.first; i < range.second; i++) {
offsets[i] = -1;
}
range = range2;
range_avg = range_avg2;
} else {
for(size_t i = range2.first; i < range2.second; i++) {
offsets[i] = -1;
}
}
} else {
range = range2;
range_avg = range_avg2;
}
}
bool weighted_avg_inited = false;
float weighted_avg = -1.0f;
for(size_t i = 0; i < offsets.size(); i++) {
if(offsets[i] < 0)
continue;
float diff = (float)offsets[i] - (float)i;
if(weighted_avg_inited) {
if(abs(diff - weighted_avg) > 20.0f) {
offsets[i] = -1;
continue;
}
}
if(weighted_avg < 0.0f) {
weighted_avg = diff;
weighted_avg_inited = true;
} else {
weighted_avg = 0.8f * weighted_avg + 0.2f * diff;
}
}
if(debug) {
cerr << "after filtering" << endl;
for(size_t j = 0; j < offsets.size(); j++) {
cout << j << ": " << offsets[j] << " " << s2[j] << ": " << (offsets[j] < 0 ? ' ' : s[offsets[j]]) << endl;
}
}
int i = 0;
while(i < offsets.size()) {
// skip non-negative offsets
for(; i < offsets.size() && offsets[i] >= 0; i++);
if(i >= offsets.size()) break;
int j = i;
for(; j < offsets.size(); j++) {
if(offsets[j] >= 0)
break;
}
assert(i >= offsets.size() || offsets[i] < 0);
assert(j >= offsets.size() || offsets[j] >= 0);
if(i > 0 && j < offsets.size()) {
i -= 1;
int left = offsets[i], right = offsets[j];
assert_geq(left, 0);
if(left > right) return false;
assert_leq(left, right);
int ref_len = right - left + 1;
int query_len = j - i + 1;
if(query_len == ref_len) { // match
for(size_t i2 = i + 1; i2 < j; i2++)
offsets[i2] = offsets[i2-1] + 1;
} else { // deletion or insertion
bool del = query_len < ref_len;
size_t gap_len = del ? ref_len - query_len : query_len - ref_len;
size_t max_len = max(ref_len, query_len);
const size_t max_mm = max_len / 25 + 1;
const size_t very_max_mm = max<size_t>(max_len / 2, max_mm);
ca_s_ = s.substr(left, max_len);
ca_s2_ = s2.substr(i, max_len);
size_t gap_pos = 0, mm = very_max_mm + 1;
find_gap_pos(ca_s_,
ca_s2_,
ca_ed_,
ca_ed2_,
del,
gap_len,
very_max_mm,
gap_pos,
mm,
debug);
if(mm > very_max_mm) {
return false;
}
assert_lt(gap_pos, query_len);
if(del) {
for(size_t i2 = i + 1; i2 < j; i2++) {
if(i2 - i == gap_pos) {
offsets[i2] = offsets[i2-1] + gap_len;
} else {
offsets[i2] = offsets[i2-1] + 1;
}
}
} else {
for(size_t i2 = i + 1; i2 < j; i2++) {
if(i2 - i >= gap_pos && i2 - i < gap_pos + gap_len) {
offsets[i2] = offsets[i2-1];
} else {
offsets[i2] = offsets[i2-1] + 1;
}
}
}
}
}
i = j;
}
if(debug) {
cerr << "final offsets" << endl;
for(size_t j = 0; j < offsets.size(); j++) {
cout << j << ": " << offsets[j] << " " << s2[j] << ": " << (offsets[j] < 0 ? ' ' : s[offsets[j]]) << endl;
}
}
#ifndef NDEBUG
for(size_t i = 1; i < offsets.size(); i++) {
if(offsets[i-1] < 0 || offsets[i] < 0) continue;
assert_leq(offsets[i-1], offsets[i]);
}
#endif
size_t b, e;
for(b = 0; b < offsets.size() && offsets[b] < 0; b++);
if(b >= offsets.size()) return false;
assert_lt((size_t)b, offsets.size());
for(e = offsets.size() - 1; e > b && offsets[e] < 0; e--);
if(b == e) return false;
for(size_t p = b; p <= e; p++) {
if(offsets[p] < 0) return false;
}
// try to fill the ends
if(b > 0) {
int mm = 0, pb = (int)b;
for(int i = pb - 1; i >= 0; i--) {
assert_geq(offsets[i+1], 0);
if(offsets[i+1] == 0) break;
if(s2[i] != s[offsets[i+1] - 1])
mm++;
if(pb - i < 25 * (mm - 1))
break;
offsets[i] = offsets[i+1] - 1;
b = (size_t)i;
}
}
if(e + 1 < offsets.size()) {
int mm = 0, prev_end = (int)e;
for(int i = prev_end + 1; i < offsets.size(); i++) {
assert_geq(offsets[i-1], 0);
if(offsets[i-1] + 1 >= s.length()) break;
if(s2[i] != s[offsets[i-1] + 1])
mm++;
if(i - prev_end < 25 * (mm - 1))
break;
offsets[i] = offsets[i-1] + 1;
e = (size_t)i;
}
}
assert_geq(seed.pos.first, (TIndexOffU)b);
seed.pos.first += (TIndexOffU)b;
seed.pos.second = seed.pos.first + e - b + 1;
seed.orig_pos.first = seed.pos.first;
seed.orig_pos.second = seed.pos.first;
seed.consensus_pos.first = offsets[b];
seed.consensus_pos.second = offsets[e] + 1;
seed.aligned = true;
for(int p = b; p < e; p++) {
assert_geq(offsets[p], 0);
assert_leq(offsets[p], offsets[p+1]);
if(offsets[p] + 1 == offsets[p+1]) { // match
continue;
} else if(offsets[p] + 1 < offsets[p+1]) { // deletion
seed.right_gaps.expand();
seed.right_gaps.back().first = p + 1 - b;
seed.right_gaps.back().second = offsets[p+1] - offsets[p] - 1;
} else {
assert_eq(offsets[p], offsets[p+1]);
int p2 = p + 1;
for(; offsets[p2] == offsets[p2+1]; p2++);
seed.right_gaps.expand();
seed.right_gaps.back().first = p + 1 - b;
seed.right_gaps.back().second = (int)p - (int)p2;
p = p2;
}
}
if(debug) {
string consensus_str = consensus_.substr(seed.consensus_pos.first, seed.consensus_pos.second - seed.consensus_pos.first);
string seed_str;
seed.getExtendedSeedSequence(ref, seed_str);
assert_eq(consensus_str.length(), seed_str.length());
cerr << "consensus: " << consensus_str << endl;
cerr << "seed : " << seed_str << endl;
}
RB_Repeat::seed_merged++;
return true;
}
bool RB_RepeatExt::isSelfRepeat(const RepeatParameter& rp,
const string& s,
const EList<pair<size_t, size_t> >& s_kmer_table,
EList<int>& offsets,
size_t k,
bool debug)
{
offsets.resize(s.length());
offsets.fill(-1);
pair<size_t, size_t> kmer(0, 0);
for(size_t i = 0; i + k <= s.length(); i++) {
if(i == 0) {
kmer.first = extract_kmer(s, i, k);
} else {
kmer.first = next_kmer(kmer.first, s[i+k-1], k);
}
size_t lb = s_kmer_table.bsearchLoBound(kmer);
while(lb < s_kmer_table.size() && kmer.first == s_kmer_table[lb].first) {
if(offsets[i] == -1) {
offsets[i] = (int)s_kmer_table[lb].second;
} else if(offsets[i] >= 0) {
offsets[i] = -2;
} else {
assert_lt(offsets[i], -1);
offsets[i] -= 1;
}
lb++;
}
if(offsets[i] > 0 && i + k == s.length()) {
for(size_t j = i + 1; j < s.length(); j++) {
offsets[j] = offsets[j-1] + 1;
}
}
}
if(debug) {
cerr << "offsets" << endl;
for(size_t j = 0; j < offsets.size(); j++) {
cout << j << ": " << offsets[j] << " " << s[j] << ": " << (offsets[j] < 0 ? ' ' : s[offsets[j]]) << endl;
}
}
const size_t min_repeat_size = 100;
size_t repeat_count = 0;
for(size_t i = 0; i + min_repeat_size < offsets.size(); i++) {
if(offsets[i] >= -1)
continue;
size_t j = i + 1;
for(; j < offsets.size() && offsets[j] <= -2; j++);
if(j - i >= min_repeat_size)
repeat_count++;
i = j;
}
return repeat_count >= 2;
}
template<typename TStr>
bool RB_RepeatExt::merge(const RepeatParameter& rp,
const TStr& s,
RB_SWAligner& swalginer,
const RB_RepeatExt& o,
bool contain,
size_t seed_i,
size_t seed_j,
bool debug)
{
// construct a new consensus sequence
string prev_consensus = consensus_;
size_t consensus_add_len = 0;
{
assert_lt(seed_i, seeds_.size());
assert_lt(seed_j, o.seeds_.size());
const SeedExt& seed = seeds_[seed_i];
const SeedExt& oseed = o.seeds_[seed_j];
Range range = seed.getExtendedRange(consensus_.length());
Range orange = oseed.getExtendedRange(o.consensus_.length());
assert_leq(range.first, orange.first + 10);
consensus_add_len = (int)orange.first - (int)range.first;
if(!contain) {
if(range.second <= orange.first) {
cerr << "something wrong: " << range.first << "-" << range.second << " ";
cerr << orange.first << "-" << orange.second << " " << o.consensus_.length() << endl;
assert(false);
}
if(range.second > orange.first &&
range.second - orange.first < o.consensus_.length()) {
consensus_ += o.consensus_.substr(range.second - orange.first);
}
}
}
EList<pair<size_t, size_t> > merge_list;
merge_list.reserveExact(o.seed_ranges_.size());
size_t p = 0, p2 = 0;
const size_t relax = 5;
while(p < seed_ranges_.size() && p2 < o.seed_ranges_.size()) {
const RB_AlleleCoord& range = seed_ranges_[p];
assert_lt(range.idx, seeds_.size());
const RB_AlleleCoord& range2 = o.seed_ranges_[p2];
assert_lt(range2.idx, o.seeds_.size());
if(range.contain(range2, relax)) {
merge_list.expand();
merge_list.back().first = p;
merge_list.back().second = p2;
if(debug) {
cerr << p << ":" << range.left << "-" << range.right << " > ";
cerr << p2 << ":" << range2.left << "-" << range2.right << endl;
}
} else if(range2.contain(range, relax)) {
merge_list.expand();
merge_list.back().first = p;
merge_list.back().second = p2;
if(debug) {
cerr << p << ":" << range.left << "-" << range.right << " < ";
cerr << p2 << ":" << range2.left << "-" << range2.right << endl;
}
} else {
TIndexOffU overlap_len = range.overlap_len(range2);
bool stored = !merge_list.empty() && merge_list.back().first == p;
bool stored2 = !merge_list.empty() && merge_list.back().second == p2;
if(overlap_len > 0) {
if(!stored && !stored2) {
merge_list.expand();
merge_list.back().first = p;
merge_list.back().second = p2;
if(debug) {
cerr << p << ":" << range.left << "-" << range.right << " ol ";
cerr << p2 << ":" << range2.left << "-" << range2.right << endl;
}
}
} else {
if(range2.right <= range.left) {
if(!stored2) {
merge_list.expand();
merge_list.back().first = numeric_limits<size_t>::max();
merge_list.back().second = p2;
if(debug) {
cerr << p << ":" << range.left << "-" << range.right << " <> ";
cerr << p2 << ":" << range2.left << "-" << range2.right << endl;
}
}
} else {
assert_leq(range.right, range2.left);
if(!stored) {
merge_list.expand();
merge_list.back().first = p;
merge_list.back().second = numeric_limits<size_t>::max();
if(debug) {
cerr << p << ":" << range.left << "-" << range.right << " <> ";
cerr << p2 << ":" << range2.left << "-" << range2.right << endl;
}
}
}
}
}
if(range.right <= range2.right) p++;
if(range2.right <= range.right) p2++;
}
assert(p == seeds_.size() || p2 == o.seeds_.size());
while(p < seed_ranges_.size()) {
bool stored = !merge_list.empty() && merge_list.back().first == p;
if(!stored) {
const RB_AlleleCoord& range = seed_ranges_[p];
merge_list.expand();
merge_list.back().first = p;
merge_list.back().second = numeric_limits<size_t>::max();
if(debug) {
cerr << p << ":" << range.left << "-" << range.right << " : <> " << endl;
}
}
p++;
}
while(p2 < o.seed_ranges_.size()) {
bool stored2 = !merge_list.empty() && merge_list.back().second == p2;
if(!stored2) {
const RB_AlleleCoord& range2 = o.seed_ranges_[p2];
merge_list.expand();
merge_list.back().first = numeric_limits<size_t>::max();
merge_list.back().second = p2;
if(debug) {
cerr << ": <> " << p2 << ":" << range2.left << "-" << range2.right << endl;
}
}
p2++;
}
assert(!merge_list.empty());
if(debug) {
for(size_t i = 0; i < merge_list.size(); i++) {
cerr << "merge list:" << endl;
cerr << "\t" << (merge_list[i].first < seed_ranges_.size() ? (int)merge_list[i].first : -1);
cerr << " " << (merge_list[i].second < o.seed_ranges_.size() ? (int)merge_list[i].second : -1) << endl;
}
}
const size_t kmer_len = 12;
EList<pair<size_t, size_t> > kmer_table;
build_kmer_table(consensus_, kmer_table, kmer_len);
EList<int> offsets;
bool self_repeat = isSelfRepeat(rp,
consensus_,
kmer_table,
offsets,
kmer_len,
debug);
if(self_repeat) {
consensus_ = prev_consensus;
return false;
}
string seq;
for(size_t i = 0; i < merge_list.size(); i++) {
size_t seed_id = (merge_list[i].first < seed_ranges_.size() ? seed_ranges_[merge_list[i].first].idx : merge_list[i].first);
size_t oseed_id = (merge_list[i].second < o.seed_ranges_.size() ? o.seed_ranges_[merge_list[i].second].idx : merge_list[i].second);
if(seed_id < seeds_.size()) {
if(oseed_id >= o.seeds_.size())
continue;
else if(seed_ranges_[merge_list[i].first].contain(o.seed_ranges_[merge_list[i].second]))
continue;
}
const SeedExt* seed = (seed_id < seeds_.size() ? &seeds_[seed_id] : NULL);
const SeedExt* oseed = (oseed_id < o.seeds_.size() ? &o.seeds_[oseed_id] : NULL);
assert(seed != NULL || oseed != NULL);
size_t left = (seed != NULL ? seed->pos.first : numeric_limits<size_t>::max());
size_t right = (seed != NULL ? seed->pos.second : 0);
int consensus_approx_left = (seed != NULL ? seed->consensus_pos.first : numeric_limits<int>::max());
int consensus_approx_right = (seed != NULL ? seed->consensus_pos.second : 0);
if(oseed != NULL) {
if(oseed->pos.first < left) {
left = oseed->pos.first;
}
if(oseed->pos.second > right) {
right = oseed->pos.second;
}
int oconsensus_approx_left = oseed->consensus_pos.first + consensus_add_len;
if(oconsensus_approx_left < consensus_approx_left) {
consensus_approx_left = oconsensus_approx_left;
}
int oconsensus_approx_right = oseed->consensus_pos.second + consensus_add_len;
if(oconsensus_approx_right > consensus_approx_right) {
consensus_approx_right = oconsensus_approx_right;
}
}
getString(s, left, right - left, seq);
if(seed_id >= seeds_.size()) {
seed_id = seeds_.size();
seeds_.expand();
}
/* bool succ = */ align(rp,
s,
consensus_,
kmer_table,
seq,
offsets,
kmer_len,
seeds_[seed_id],
consensus_approx_left,
consensus_approx_right,
left,
right,
debug && right - left == 1080);
}
while(true) {
internal_update();
size_t remove_count = 0;
for(size_t i = 0; i + 1 < seed_ranges_.size(); i++) {
RB_AlleleCoord& range = seed_ranges_[i];
if(range.left == numeric_limits<size_t>::max())
break;
for(size_t j = i + 1; j < seed_ranges_.size(); j++) {
RB_AlleleCoord& range2 = seed_ranges_[j];
if(range2.left == numeric_limits<size_t>::max())
break;
if(range.right <= range2.left)
break;
getString(s, range.left, range2.right - range.left, seq);
/* bool succ = */ align(rp,
s,
consensus_,
kmer_table,
seq,
offsets,
kmer_len,
seeds_[range.idx],
seeds_[range.idx].consensus_pos.first,
seeds_[range2.idx].consensus_pos.second,
range.left,
range2.right,
debug && range.left == 692422);
range.left = seeds_[range.idx].pos.first;
range.right = seeds_[range.idx].pos.second;
seeds_[range2.idx].reset();
range2.left = numeric_limits<size_t>::max();
remove_count++;
}
}
if(remove_count <= 0) break;
sort(seeds_.begin(), seeds_.end(), seedCmp);
seeds_.resize(seeds_.size() - remove_count);
}
return true;
}
#define DMAX std::numeric_limits<double>::max()
RB_SWAligner::RB_SWAligner()
{
rnd_.init(0);
}
RB_SWAligner::~RB_SWAligner()
{
if(sc_) {
delete sc_;
}
}
void RB_SWAligner::init_dyn(const RepeatParameter& rp)
{
const int MM_PEN = 3;
// const int MM_PEN = 6;
const int GAP_PEN_LIN = 2;
// const int GAP_PEN_LIN = (((MM_PEN) * rpt_edit_ + 1) * 1.0);
const int GAP_PEN_CON = 4;
// const int GAP_PEN_CON = (((MM_PEN) * rpt_edit_ + 1) * 1.0);
const int MAX_PEN = MAX_I16;
scoreMin_.init(SIMPLE_FUNC_LINEAR, rp.max_edit * MM_PEN * -1.0, 0.0);
nCeil_.init(SIMPLE_FUNC_LINEAR, 0.0, 0.0);
penCanIntronLen_.init(SIMPLE_FUNC_LOG, -8, 1);
penNoncanIntronLen_.init(SIMPLE_FUNC_LOG, -8, 1);
sc_ = new Scoring(
DEFAULT_MATCH_BONUS, // constant reward for match
DEFAULT_MM_PENALTY_TYPE, // how to penalize mismatches
MM_PEN, // max mm penalty
MM_PEN, // min mm penalty
MAX_PEN, // max sc penalty
MAX_PEN, // min sc penalty
scoreMin_, // min score as function of read len
nCeil_, // max # Ns as function of read len
DEFAULT_N_PENALTY_TYPE, // how to penalize Ns in the read
DEFAULT_N_PENALTY, // constant if N pelanty is a constant
DEFAULT_N_CAT_PAIR, // whether to concat mates before N filtering
GAP_PEN_CON, // constant coeff for read gap cost
GAP_PEN_CON, // constant coeff for ref gap cost
GAP_PEN_LIN, // linear coeff for read gap cost
GAP_PEN_LIN, // linear coeff for ref gap cost
1 /* gGapBarrier */ // # rows at top/bot only entered diagonally
);
}
void RB_SWAligner::makePadString(const string& ref,
const string& read,
string& pad,
size_t len)
{
pad.resize(len);
for(size_t i = 0; i < len; i++) {
// shift A->C, C->G, G->T, T->A
pad[i] = "CGTA"[asc2dna[ref[i]]];
if(read[i] == pad[i]) {
// shift
pad[i] = "CGTA"[asc2dna[pad[i]]];
}
}
int head_len = len / 2;
size_t pad_start = len - head_len;
for(size_t i = 0; i < head_len; i++) {
if(read[i] == pad[pad_start + i]) {
// shift
pad[pad_start + i] = "CGTA"[asc2dna[pad[pad_start + i]]];
}
}
}
int RB_SWAligner::alignStrings(const string &ref,
const string &read,
EList<Edit>& edits,
Coord& coord)
{
// Prepare Strings
// Read -> BTDnaString
// Ref -> bit-encoded string
//SwAligner swa;
BTDnaString btread;
BTString btqual;
BTString btref;
BTString btref2;
BTDnaString btreadrc;
BTString btqualrc;
string qual = "";
for(size_t i = 0; i < read.length(); i++) {
qual.push_back('I');
}
#if 0
cerr << "REF : " << ref << endl;
cerr << "READ: " << read << endl;
cerr << "QUAL: " << qual << endl;
#endif
btread.install(read.c_str(), true);
btreadrc = btread;
btreadrc.reverseComp();
btqual.install(qual.c_str());
btqualrc = btqual;
btref.install(ref.c_str());
TAlScore min_score = sc_->scoreMin.f<TAlScore >((double)btread.length());
btref2 = btref;
size_t nceil = 0;
// size_t nrow = btread.length();
// Convert reference string to mask
for(size_t i = 0; i < btref2.length(); i++) {
if(toupper(btref2[i]) == 'N') {
btref2.set(16, i);
} else {
int num = 0;
int alts[] = {4, 4, 4, 4};
decodeNuc(toupper(btref2[i]), num, alts);
assert_leq(num, 4);
assert_gt(num, 0);
btref2.set(0, i);
for(int j = 0; j < num; j++) {
btref2.set(btref2[i] | (1 << alts[j]), i);
}
}
}
bool fw = true;
uint32_t refidx = 0;
swa.initRead(
btread, // read sequence
btreadrc,
btqual, // read qualities
btqualrc,
0, // offset of first character within 'read' to consider
btread.length(), // offset of last char (exclusive) in 'read' to consider
*sc_); // local-alignment score floor
DynProgFramer dpframe(false);
size_t readgaps = 0;
size_t refgaps = 0;
size_t maxhalf = 0;
DPRect rect;
dpframe.frameSeedExtensionRect(
0, // ref offset implied by seed hit assuming no gaps
btread.length(), // length of read sequence used in DP table
btref2.length(), // length of reference
readgaps, // max # of read gaps permitted in opp mate alignment
refgaps, // max # of ref gaps permitted in opp mate alignment
(size_t)nceil, // # Ns permitted
maxhalf, // max width in either direction
rect); // DP Rectangle
assert(rect.repOk());
size_t cminlen = 2000, cpow2 = 4;
swa.initRef(
fw, // whether to align forward or revcomp read
refidx, // reference ID
rect, // DP rectangle
btref2.wbuf(), // reference strings
0, // offset of first reference char to align to
btref2.length(), // offset of last reference char to align to
btref2.length(), // length of reference sequence
*sc_, // scoring scheme
min_score, // minimum score
true, // use 8-bit SSE if positions
cminlen, // minimum length for using checkpointing scheme
cpow2, // interval b/t checkpointed diags; 1 << this
false, // triangular mini-fills?
false // is this a seed extension?
);
TAlScore best = std::numeric_limits<TAlScore>::min();
bool found = swa.align(rnd_, best);
#ifdef DEBUGLOG
cerr << "found: " << found << "\t" << best << "\t" << "minsc: " << min_score << endl;
#endif
if (found) {
#ifdef DEBUGLOG
//cerr << "CP " << "found: " << found << "\t" << best << "\t" << "minsc: " << min_score << endl;
cerr << "REF : " << ref << endl;
cerr << "READ: " << read << endl;
#endif
SwResult res;
int max_match_len = 0;
res.reset();
res.alres.init_raw_edits(&rawEdits_);
found = swa.nextAlignment(res, best, rnd_);
if (found) {
edits = res.alres.ned();
//const TRefOff ref_off = res.alres.refoff();
//const Coord& coord = res.alres.refcoord();
coord = res.alres.refcoord();
//assert_geq(genomeHit._joinedOff + coord.off(), genomeHit.refoff());
#ifdef DEBUGLOG
cerr << "num edits: " << edits.size() << endl;
cerr << "coord: " << coord.off();
cerr << ", " << coord.ref();
cerr << ", " << coord.orient();
cerr << ", " << coord.joinedOff();
cerr << endl;
Edit::print(cerr, edits); cerr << endl;
Edit::printQAlign(cerr, btread, edits);
#endif
max_match_len = getMaxMatchLen(edits, btread.length());
#ifdef DEBUGLOG
cerr << "max match length: " << max_match_len << endl;
#endif
}
#ifdef DEBUGLOG
cerr << "nextAlignment: " << found << endl;
cerr << "-------------------------" << endl;
#endif
}
return 0;
}
void RB_SWAligner::doTest(const RepeatParameter& rp,
const string& refstr,
const string& readstr)
{
init_dyn(rp);
doTestCase1(refstr,
readstr,
rp.max_edit);
}
void RB_SWAligner::doTestCase1(const string& refstr,
const string& readstr,
TIndexOffU rpt_edit)
{
cerr << "doTestCase1----------------" << endl;
EList<Edit> edits;
Coord coord;
if (refstr.length() == 0 ||
readstr.length() == 0) {
return;
}
EList<Edit> ed;
string pad;
makePadString(refstr, readstr, pad, 5);
string ref2 = pad + refstr + pad;
string read2 = pad + readstr + pad;
alignStrings(refstr, readstr, edits, coord);
size_t left = pad.length();
size_t right = left + readstr.length();
edits.reserveExact(ed.size());
for(size_t i = 0; i < ed.size(); i++) {
if(ed[i].pos >= left && ed[i].pos <= right) {
edits.push_back(ed[i]);
edits.back().pos -= left;
}
}
#if 0
RepeatGroup rg;
rg.edits = edits;
rg.coord = coord;
rg.seq = readstr;
rg.base_offset = 0;
string chr_name = "rep";
cerr << "REF : " << refstr << endl;
cerr << "READ: " << readstr << endl;
size_t snpids = 0;
rg.buildSNPs(snpids);
rg.writeSNPs(cerr, chr_name); cerr << endl;
#endif
}
template<typename TStr>
RepeatBuilder<TStr>::RepeatBuilder(TStr& s,
TStr& sOriginal,
const EList<RefRecord>& szs,
const EList<string>& ref_names,
bool forward_only,
const string& filename) :
s_(s),
sOriginal_(sOriginal),
coordHelper_(s.length(), forward_only ? s.length() : s.length() / 2, szs, ref_names),
forward_only_(forward_only),
filename_(filename),
forward_length_(forward_only ? s.length() : s.length() / 2)
{
cerr << "RepeatBuilder: " << filename_ << endl;
}
template<typename TStr>
RepeatBuilder<TStr>::~RepeatBuilder()
{
for(map<size_t, RB_Repeat*>::iterator it = repeat_map_.begin(); it != repeat_map_.end(); it++) {
delete it->second;
}
repeat_map_.clear();
}
template<typename TStr>
void RepeatBuilder<TStr>::readSA(const RepeatParameter& rp,
BlockwiseSA<TStr>& sa)
{
TIndexOffU count = 0;
subSA_.init(s_.length() + 1, rp.min_repeat_len, rp.repeat_count);
while(count < s_.length() + 1) {
TIndexOffU saElt = sa.nextSuffix();
count++;
if(count && (count % 10000000 == 0)) {
cerr << "SA count " << count << endl;
}
if(saElt == s_.length()) {
assert_eq(count, s_.length() + 1);
break;
}
subSA_.push_back(s_, coordHelper_, saElt, count == s_.length());
}
cerr << "subSA size is " << subSA_.size() << endl;
subSA_.dump();
#if 0
for(size_t i = 0; i < subSA_.size(); i++) {
TIndexOffU joinedOff = subSA_[i];
fp << setw(10) << joinedOff << " " << getString(s_, joinedOff, rp.seed_len) << endl;
}
#endif
cerr << "subSA mem Usage: " << subSA_.getMemUsage() << endl << endl;
}
template<typename TStr>
void RepeatBuilder<TStr>::readSA(const RepeatParameter &rp,
const BitPackedArray &sa)
{
TIndexOffU count = 0;
subSA_.init(s_.length() + 1, rp.min_repeat_len, rp.repeat_count);
for(size_t i = 0; i < sa.size(); i++) {
TIndexOffU saElt = sa[i];
count++;
if(count && (count % 10000000 == 0)) {
cerr << "RB count " << count << endl;
}
if(saElt == s_.length()) {
assert_eq(count, s_.length() + 1);
break;
}
subSA_.push_back(s_, coordHelper_, saElt, count == s_.length());
}
cerr << "subSA size: " << endl;
subSA_.dump();
#if 0
for(size_t i = 0; i < subSA_.size(); i++) {
TIndexOffU joinedOff = subSA_[i];
fp << setw(10) << joinedOff << " " << getString(s_, joinedOff, rp.seed_len) << endl;
}
#endif
cerr << "subSA mem Usage: " << subSA_.getMemUsage() << endl << endl;
}
template<typename TStr>
void RepeatBuilder<TStr>::build(const RepeatParameter& rp)
{
string rpt_len_str;
rpt_len_str = to_string(rp.min_repeat_len) + "-" + to_string(rp.max_repeat_len);
string seed_filename = filename_ + ".rep." + rpt_len_str + ".seed";
ofstream fp(seed_filename.c_str());
swaligner_.init_dyn(rp);
RB_RepeatManager* repeat_manager = new RB_RepeatManager;
EList<RB_RepeatBase> repeatBases;
subSA_.buildRepeatBase(s_,
coordHelper_,
rp.max_repeat_len,
repeatBases);
size_t repeat_id = 0;
for(size_t i = 0; i < repeatBases.size(); i++) {
addRepeatGroup(rp,
repeat_id,
*repeat_manager,
repeatBases[i],
fp);
}
{
// Build and test minimizer-based k-mer table
const size_t window = RB_Minimizer<TStr>::default_w;
const size_t k = RB_Minimizer<TStr>::default_k;
RB_KmerTable kmer_table;
EList<string> seqs;
seqs.reserveExact(repeat_map_.size());
for(map<size_t, RB_Repeat*>::const_iterator it = repeat_map_.begin(); it != repeat_map_.end(); it++) {
const RB_Repeat& repeat = *(it->second);
assert(repeat.satisfy(rp));
seqs.expand();
seqs.back() = repeat.consensus();
}
kmer_table.build(seqs,
window,
k);
kmer_table.dump(cerr);
cerr << endl;
string query, rc_query;
string queryOriginal;
EList<pair<uint64_t, size_t> > minimizers;
size_t total = 0, num_repeat = 0, correct = 0, false_positive = 0, false_negative = 0;
for(size_t i = 0; i + rp.min_repeat_len <= forward_length_; i += 1000) {
if(coordHelper_.getEnd(i) != coordHelper_.getEnd(i + rp.min_repeat_len))
continue;
query = getString(s_, i, rp.min_repeat_len);
queryOriginal = getString(sOriginal_, i, rp.min_repeat_len);
rc_query = reverseComplement(queryOriginal);
TIndexOffU idx = subSA_.find_repeat_idx(s_, query);
const EList<TIndexOffU>& test_repeat_index = subSA_.getRepeatIndex();
bool repeat = (idx < test_repeat_index.size());
bool est_repeat = kmer_table.isRepeat(query,
rc_query,
minimizers);
total++;
if(repeat) num_repeat++;
if(repeat == est_repeat) {
correct++;
} else {
if(est_repeat) {
false_positive++;
} else {
false_negative++;
//assert(false);
}
}
}
cerr << "total: " << total << endl;
cerr << "repeat: " << num_repeat << endl;
cerr << "correct: " << correct << endl;
cerr << "false positive: " << false_positive << endl;
cerr << "false negative: " << false_negative << endl;
cerr << endl;
ELList<RB_Alignment> position2D; EList<RB_Alignment> alignments;
size_t repeat_total = 0, repeat_aligned = 0;
const EList<TIndexOffU>& test_repeat_index = subSA_.getRepeatIndex();
size_t interval = 1;
if(test_repeat_index.size() >= 1000) {
interval = test_repeat_index.size() / 1000;
}
size_t total_alignments = 0, max_alignments = 0;
string query2, rc_query2;
string queryOriginal2, rc_queryOriginal2;
for(size_t i = 0; i < test_repeat_index.size(); i += interval) {
TIndexOffU saElt_idx = test_repeat_index[i];
TIndexOffU saElt_idx_end = (i + 1 < test_repeat_index.size() ? test_repeat_index[i+1] : subSA_.size());
TIndexOffU saElt_size = saElt_idx_end - saElt_idx;
TIndexOffU saElt = subSA_[saElt_idx];
query = getString(s_, saElt, rp.min_repeat_len);
query2 = query;
queryOriginal = getString(sOriginal_, saElt, rp.min_repeat_len);
queryOriginal2 = queryOriginal;
// introduce three mismatches into a query when the query is at lest 100-bp
if(rp.min_repeat_len >= 100) {
const size_t mid_pos1 = (size_t)(rp.min_repeat_len * 0.1);
if(query2[mid_pos1] == 'A') {
query2[mid_pos1] = 'C';
queryOriginal2[mid_pos1] ='C';
} else {
query2[mid_pos1] = 'A';
queryOriginal2[mid_pos1] = 'A';
}
const size_t mid_pos2 = (size_t)(rp.min_repeat_len * 0.5);
if(query2[mid_pos2] == 'C') {
query2[mid_pos2] = 'G';
queryOriginal2[mid_pos2] = 'G';
} else {
query2[mid_pos2] = 'C';
queryOriginal2[mid_pos1] = 'G';
}
const size_t mid_pos3 = (size_t)(rp.min_repeat_len * 0.9);
if(query2[mid_pos3] == 'G') {
query2[mid_pos3] = 'T';
queryOriginal2[mid_pos3] = 'T';
} else {
query2[mid_pos3] = 'G';
queryOriginal2[mid_pos3] = 'G';
}
}
repeat_total += saElt_size;
size_t found = 0;
size_t cur_alignments = 0;
if(kmer_table.isRepeat(query2, minimizers)) {
kmer_table.findAlignments(query2,
minimizers,
position2D,
alignments);
total_alignments += (alignments.size() * saElt_size);
cur_alignments = alignments.size();
TIndexOffU baseoff = 0;
for(size_t s = 0; s < seqs.size(); s++) {
int spos = seqs[s].find(query);
if(spos != string::npos) {
for(size_t a = 0; a < alignments.size(); a++) {
if(alignments[a].pos == baseoff + spos) {
found++;
}
}
}
baseoff += seqs[s].length();
}
assert_leq(found, 1);
}
rc_query = reverseComplement(queryOriginal);
rc_query2 = reverseComplement(queryOriginal2);
size_t rc_found = 0;
if(kmer_table.isRepeat(rc_query2, minimizers)) {
kmer_table.findAlignments(rc_query2,
minimizers,
position2D,
alignments);
total_alignments += (alignments.size() * saElt_size);
cur_alignments += alignments.size();
if(cur_alignments > max_alignments) {
max_alignments = cur_alignments;
}
TIndexOffU baseoff = 0;
for(size_t s = 0; s < seqs.size(); s++) {
int spos = seqs[s].find(rc_query);
if(spos != string::npos) {
for(size_t a = 0; a < alignments.size(); a++) {
if(alignments[a].pos == baseoff + spos) {
rc_found++;
}
}
}
baseoff += seqs[s].length();
}
assert_leq(rc_found, 1);
}
if(found + rc_found == 1) {
repeat_aligned += saElt_size;
}
}
cerr << "num repeats: " << repeat_total << endl;
cerr << "repeat aligned using minimizers: " << repeat_aligned << endl;
cerr << "average number of alignments: " << (float)total_alignments / repeat_total << endl;
cerr << "max alignment: " << max_alignments << endl;
cerr << endl;
}
cerr << "number of repeats is " << repeat_map_.size() << endl;
size_t total_rep_seq_len = 0, total_allele_seq_len = 0;
size_t i = 0;
for(map<size_t, RB_Repeat*>::iterator it = repeat_map_.begin(); it != repeat_map_.end(); it++, i++) {
RB_Repeat& repeat = *(it->second);
if(!repeat.satisfy(rp))
continue;
repeat.saveSeedExtension(rp,
s_,
coordHelper_,
i,
fp,
total_rep_seq_len,
total_allele_seq_len);
}
size_t total_qual_seeds = 0;
for(map<size_t, RB_Repeat*>::iterator it = repeat_map_.begin(); it != repeat_map_.end(); it++, i++) {
RB_Repeat& repeat = *(it->second);
EList<SeedExt>& seeds = repeat.seeds();
for(size_t i = 0; i < seeds.size(); i++) {
if(seeds[i].getLength() < rp.min_repeat_len)
continue;
total_qual_seeds++;
}
}
cerr << "total repeat sequence length: " << total_rep_seq_len << endl;
cerr << "total allele sequence length: " << total_allele_seq_len << endl;
cerr << "total number of seeds including those that are position-wise different, but sequence-wise identical: " << total_qual_seeds << endl;
cerr << endl;
fp << "total repeat sequence length: " << total_rep_seq_len << endl;
fp << "total allele sequence length: " << total_allele_seq_len << endl;
fp << "total number of seeds including those that are position-wise different, but sequence-wise identical: " << total_qual_seeds << endl;
fp.close();
delete repeat_manager;
repeat_manager = NULL;
// remove non-qualifying repeats
// and update repeat IDs for those remaining
{
map<size_t, RB_Repeat*> temp_repeat_map;
size_t i = 0;
for(map<size_t, RB_Repeat*>::iterator it = repeat_map_.begin(); it != repeat_map_.end(); it++) {
RB_Repeat* repeat = it->second;
if(!repeat->satisfy(rp)) {
delete repeat;
continue;
}
repeat->repeat_id(i);
temp_repeat_map[repeat->repeat_id()] = repeat;
i++;
}
repeat_map_ = temp_repeat_map;
}
const bool sanity_check = true;
if(sanity_check) {
EList<pair<size_t, size_t> > kmer_table;
string query;
string queryOriginal;
size_t total = 0, match = 0;
EList<TIndexOffU> positions;
const EList<TIndexOffU>& test_repeat_index = subSA_.getRepeatIndex();
size_t interval = 1;
if(test_repeat_index.size() >= 10000) {
interval = test_repeat_index.size() / 10000;
}
for(size_t i = 0; i < test_repeat_index.size(); i += interval) {
TIndexOffU saElt_idx = test_repeat_index[i];
TIndexOffU saElt_idx_end = (i + 1 < test_repeat_index.size() ? test_repeat_index[i+1] : subSA_.size());
positions.clear();
for(size_t j = saElt_idx; j < saElt_idx_end; j++) {
positions.push_back(subSA_[j]);
#ifndef NDEBUG
if(j > saElt_idx) {
TIndexOffU lcp_len = getLCP(s_,
coordHelper_,
positions[0],
positions.back(),
rp.min_repeat_len);
assert_eq(lcp_len, rp.min_repeat_len);
}
TIndexOffU saElt = subSA_[j];
TIndexOffU start = coordHelper_.getStart(saElt);
TIndexOffU start2 = coordHelper_.getStart(saElt + rp.min_repeat_len - 1);
assert_eq(start, start2);
#endif
}
TIndexOffU saElt = subSA_[saElt_idx];
size_t true_count = saElt_idx_end - saElt_idx;
getString(s_, saElt, rp.min_repeat_len, query);
getString(sOriginal_, saElt, rp.min_repeat_len, queryOriginal);
total++;
size_t count = 0, rc_count = 0;
string rc_query = reverse_complement(queryOriginal);
for(map<size_t, RB_Repeat*>::iterator it = repeat_map_.begin(); it != repeat_map_.end(); it++) {
RB_Repeat& repeat = *(it->second);
int pos = repeat.consensus().find(query);
if(pos != string::npos) {
for(size_t s = 0; s < repeat.seeds().size(); s++) {
SeedExt& seed = repeat.seeds()[s];
string seq;
seed.getExtendedSeedSequence(s_, seq);
if(seq.find(query) != string::npos)
count++;
}
}
pos = repeat.consensus().find(rc_query);
if(pos != string::npos) {
for(size_t s = 0; s < repeat.seeds().size(); s++) {
SeedExt& seed = repeat.seeds()[s];
string seq;
seed.getExtendedSeedSequence(s_, seq);
if(seq.find(rc_query) != string::npos)
rc_count++;
}
}
}
if(count == true_count || rc_count == true_count) {
match++;
} else if(total - match <= 10) {
cerr << " query: " << query << endl;
cerr << "rc_query: " << rc_query << endl;
cerr << "true count: " << true_count << endl;
cerr << "found count: " << count << endl;
cerr << "rc found count: " << rc_count << endl;
cerr << endl;
}
}
cerr << "RepeatBuilder: sanity check: " << match << " passed (out of " << total << ")" << endl << endl;
}
}
template<typename TStr>
void RepeatBuilder<TStr>::writeHaploType(const EList<SeedHP>& haplo_lists,
const EList<SeedExt>& seeds,
TIndexOffU& hapl_id_base,
const string& seq_name,
ostream &fp)
{
if(haplo_lists.size() == 0) {
return;
}
for(size_t i = 0; i < haplo_lists.size(); i++) {
const SeedHP &haploType = haplo_lists[i];
fp << "rpht" << hapl_id_base++;
fp << "\t" << seq_name;
fp << "\t" << haploType.range.first;
fp << "\t" << haploType.range.second;
fp << "\t";
assert_gt(haploType.snpIDs.size(), 0);
for (size_t j = 0; j < haploType.snpIDs.size(); j++) {
if(j) {
fp << ",";
}
fp << haploType.snpIDs[j];
}
fp << endl;
}
}
template<typename TStr>
void RepeatBuilder<TStr>::writeAllele(TIndexOffU grp_id,
TIndexOffU allele_id,
Range range,
const string& seq_name,
TIndexOffU baseoff,
const EList<SeedExt>& seeds,
ostream &fp)
{
// >rpt_name*0\trep\trep_pos\trep_len\tpos_count\t0
// chr_name:pos:direction chr_name:pos:direction
//
// >rep1*0 rep 0 100 470 0
// 22:112123123:+ 22:1232131113:+
//
size_t snp_size = seeds[range.first].snps.size();
size_t pos_size = range.second - range.first;
fp << ">";
fp << "rpt_" << grp_id << "*" << allele_id;
fp << "\t" << seq_name;
fp << "\t" << baseoff + seeds[range.first].consensus_pos.first;
fp << "\t" << seeds[range.first].consensus_pos.second - seeds[range.first].consensus_pos.first;
fp << "\t" << pos_size;
fp << "\t" << snp_size;
fp << "\t";
for(size_t i = 0; i < snp_size; i++) {
if(i > 0) {fp << ",";}
fp << "rps" << seeds[range.first].snps[i]->id;
}
fp << endl;
// print positions
for(size_t i = 0; i < pos_size; i++) {
if(i > 0 && (i % 10 == 0)) {
fp << endl;
}
if(i % 10 != 0) {
fp << " ";
}
string chr_name;
TIndexOffU pos_in_chr;
TIndexOffU joinedOff = seeds[range.first + i].pos.first;
bool fw = true;
if(joinedOff < forward_length_) {
fw = true;
} else {
fw = false;
joinedOff = s_.length() - joinedOff - seeds[range.first].getLength();
}
coordHelper_.getGenomeCoord(joinedOff, chr_name, pos_in_chr);
char direction = fw ? '+' : '-';
fp << chr_name << ":" << pos_in_chr << ":" << direction;
}
fp << endl;
}
template<typename TStr>
void RepeatBuilder<TStr>::writeSNPs(ostream& fp,
const string& rep_chr_name,
const ESet<Edit>& editset)
{
for(size_t i = 0; i < editset.size(); i++) {
const Edit& ed = editset[i];
if(ed.isMismatch()) {
fp << "rps" << ed.snpID;
fp << "\t" << "single";
fp << "\t" << rep_chr_name;
fp << "\t" << ed.pos;
fp << "\t" << ed.qchr;
fp << endl;
} else {
assert(false);
}
}
}
template<typename TStr>
void RepeatBuilder<TStr>::saveFile(const RepeatParameter& rp)
{
saveRepeats(rp);
}
bool RB_RepeatManager::checkRedundant(const RepeatParameter& rp,
const map<size_t, RB_Repeat*>& repeat_map,
const EList<TIndexOffU>& positions,
EList<size_t>& to_remove) const
{
to_remove.clear();
bool replace = false;
for(size_t i = 0; i < positions.size(); i++) {
TIndexOffU seed_pos = positions[i];
Range seed_range(seed_pos + rp.seed_len, 0);
map<Range, EList<size_t> >::const_iterator it = range_to_repeats_.upper_bound(seed_range);
if(it == range_to_repeats_.end())
continue;
for(map<Range, EList<size_t> >::const_reverse_iterator rit(it); rit != range_to_repeats_.rend(); rit++) {
Range repeat_range = rit->first;
if(repeat_range.first + rp.max_repeat_len <= seed_pos)
break;
const EList<size_t>& repeat_ids = rit->second;
assert_gt(repeat_ids.size(), 0);
for(size_t r = 0; r < repeat_ids.size(); r++) {
size_t repeat_id = repeat_ids[r];
size_t idx = to_remove.bsearchLoBound(repeat_id);
if(idx < to_remove.size() && to_remove[idx] == repeat_id)
continue;
bool overlap = seed_pos < repeat_range.second && seed_pos + rp.seed_len > repeat_range.first;
if(!overlap)
continue;
map<size_t, RB_Repeat*>::const_iterator it2 = repeat_map.find(repeat_id);
assert(it2 != repeat_map.end());
const RB_Repeat& repeat = *(it2->second);
const EList<RB_AlleleCoord>& allele_ranges = repeat.seed_ranges();
size_t num_contain = 0, num_overlap = 0, num_close = 0, num_overlap_bp = 0;
size_t p = 0, p2 = 0;
while(p < positions.size() && p2 < allele_ranges.size()) {
RB_AlleleCoord range;
range.left = positions[p];
range.right = positions[p] + rp.seed_len;
RB_AlleleCoord range2 = allele_ranges[p2];
if(range2.contain(range)) {
num_contain++;
num_overlap_bp += rp.seed_len;
} else {
TIndexOffU overlap = range2.overlap_len(range);
if(overlap > 0) {
num_overlap++;
num_overlap_bp += (range2.right - range.left);
} else if(range.right + 10 > range2.left && range2.right + 10 > range.left) {
num_close++;
}
}
if(range.right <= range2.right) p++;
else p2++;
}
// if the number of matches is >= 90% of positions in the smaller group
if((num_contain + num_overlap) * 10 + num_close * 8 >= min(positions.size(), allele_ranges.size()) * 9) {
if(positions.size() <= allele_ranges.size()) {
return true;
} else {
replace = true;
to_remove.push_back(repeat_id);
to_remove.sort();
}
}
}
}
// DK - check this out
if(replace)
break;
}
return false;
}
void RB_RepeatManager::addRepeat(const RB_Repeat* repeat)
{
const EList<RB_AlleleCoord>& allele_ranges = repeat->seed_ranges();
for(size_t i = 0; i < allele_ranges.size(); i++) {
Range allele_range(allele_ranges[i].left, allele_ranges[i].right);
addRepeat(allele_range, repeat->repeat_id());
}
}
void RB_RepeatManager::addRepeat(Range range, size_t repeat_id)
{
if(range_to_repeats_.find(range) == range_to_repeats_.end()) {
range_to_repeats_[range] = EList<size_t>();
}
EList<size_t>& repeat_ids = range_to_repeats_[range];
size_t idx = repeat_ids.bsearchLoBound(repeat_id);
if(idx < repeat_ids.size() && repeat_ids[idx] == repeat_id)
return;
repeat_ids.push_back(repeat_id);
repeat_ids.sort();
}
void RB_RepeatManager::removeRepeat(const RB_Repeat* repeat)
{
const EList<RB_AlleleCoord>& allele_ranges = repeat->seed_ranges();
for(size_t p = 0; p < allele_ranges.size(); p++) {
Range range(allele_ranges[p].left, allele_ranges[p].right);
removeRepeat(range, repeat->repeat_id());
}
}
void RB_RepeatManager::removeRepeat(Range range, size_t repeat_id)
{
EList<size_t>& repeat_ids = range_to_repeats_[range];
TIndexOffU idx = repeat_ids.bsearchLoBound(repeat_id);
if(idx < repeat_ids.size() && repeat_id == repeat_ids[idx]) {
repeat_ids.erase(idx);
if(repeat_ids.empty())
range_to_repeats_.erase(range);
}
}
void RB_RepeatManager::showInfo(const RepeatParameter& rp,
CoordHelper& coordHelper,
const map<size_t, RB_Repeat*>& repeat_map,
size_t num_case) const
{
size_t count = 0;
for(map<Range, EList<size_t> >::const_iterator it = range_to_repeats_.begin();
it != range_to_repeats_.end(); it++) {
const Range& range = it->first;
if(range.second - range.first < rp.min_repeat_len) continue;
map<Range, EList<size_t> >::const_iterator jt = it; jt++;
for(; jt != range_to_repeats_.end(); jt++) {
const Range& range2 = jt->first;
if(range2.second - range2.first < rp.min_repeat_len) continue;
if(range.second <= range2.first)
break;
if(count < num_case) {
cerr << "range (" << range.first << ", " << range.second << ") vs. range2 (";
cerr << range2.first << ", " << range2.second << ")" << endl;
for(size_t i = 0; i < it->second.size(); i++) {
cerr << "\t1 " << it->second[i] << endl;
const RB_Repeat* repeat = repeat_map.find(it->second[i])->second;
repeat->showInfo(rp, coordHelper);
}
for(size_t i = 0; i < jt->second.size(); i++) {
cerr << "\t2 " << jt->second[i] << endl;
const RB_Repeat* repeat = repeat_map.find(jt->second[i])->second;
repeat->showInfo(rp, coordHelper);
}
cerr << endl << endl;
}
count++;
}
}
cerr << "ShowInfo - count: " << count << endl;
}
template<typename TStr>
void RepeatBuilder<TStr>::reassignSeeds(const RepeatParameter& rp,
size_t repeat_bid, // repeat begin id
size_t repeat_eid, // repeat end id
EList<SeedExt>& seeds)
{
assert_lt(repeat_bid, repeat_eid);
EList<bool> updated;
updated.resizeExact(repeat_eid - repeat_bid);
updated.fillZero();
string seq;
for(size_t s = 0; s < seeds.size(); s++) {
SeedExt& seed = seeds[s];
size_t max_repeat_id = repeat_bid;
size_t max_ext_len = 0;
for(size_t i = repeat_bid; i < repeat_eid; i++) {
map<size_t, RB_Repeat*>::iterator it = repeat_map_.find(i);
assert(it != repeat_map_.end());
const RB_Repeat& repeat = *(it->second);
const string& consensus = repeat.consensus();
const size_t ext_len = (consensus.length() - rp.seed_len) / 2;
if(seed.bound.first + ext_len > seed.pos.first)
continue;
if(seed.pos.second + ext_len > seed.bound.second)
continue;
getString(s_, seed.pos.first - ext_len, rp.seed_len + ext_len * 2, seq);
assert_eq(consensus.length(), seq.length());
size_t tmp_ext_len = 0;
for(size_t j = 0; j < ext_len; j++, tmp_ext_len++) {
if(seq[ext_len - j - 1] != consensus[ext_len - j - 1] ||
seq[ext_len + rp.seed_len + j] != consensus[ext_len + rp.seed_len + j]) {
break;
}
}
if(tmp_ext_len > max_ext_len) {
max_repeat_id = i;
max_ext_len = tmp_ext_len;
}
}
if(rp.seed_len + max_ext_len * 2 >= rp.min_repeat_len) {
seed.pos.first -= max_ext_len;
seed.pos.second += max_ext_len;
map<size_t, RB_Repeat*>::iterator it = repeat_map_.find(max_repeat_id);
assert(it != repeat_map_.end());
RB_Repeat& repeat = *(it->second);
const string& consensus = repeat.consensus();
const size_t ext_len = (consensus.length() - rp.seed_len) / 2;
assert_leq(max_ext_len, ext_len);
seed.consensus_pos.first = ext_len - max_ext_len;
seed.consensus_pos.second = consensus.length() - (ext_len - max_ext_len);
assert_leq(seed.consensus_pos.second - seed.consensus_pos.first, consensus.length());
repeat.addSeed(seed);
updated[max_repeat_id - repeat_bid] = true;
}
}
for(size_t i = repeat_bid; i < repeat_eid; i++) {
if(!updated[i - repeat_bid])
continue;
map<size_t, RB_Repeat*>::iterator it = repeat_map_.find(i);
assert(it != repeat_map_.end());
RB_Repeat& repeat = *(it->second);
repeat.update();
}
}
/**
* TODO
* @brief
*
* @param rpt_seq
* @param rpt_range
*/
template<typename TStr>
void RepeatBuilder<TStr>::addRepeatGroup(const RepeatParameter& rp,
size_t& repeat_id,
RB_RepeatManager& repeat_manager,
const RB_RepeatBase& repeatBase,
ostream& fp)
{
RB_Repeat* repeat = new RB_Repeat;
repeat->repeat_id(repeat_id);
repeat->init(rp,
s_,
coordHelper_,
subSA_,
repeatBase);
assert(repeat_map_.find(repeat->repeat_id()) == repeat_map_.end());
repeat_map_[repeat->repeat_id()] = repeat;
repeat_id++;
}
template<typename TStr>
bool RepeatBuilder<TStr>::checkSequenceMergeable(const string& ref,
const string& read,
EList<Edit>& edits,
Coord& coord,
TIndexOffU rpt_len,
TIndexOffU max_edit)
{
size_t max_matchlen = 0;
EList<Edit> ed;
string pad;
swaligner_.makePadString(ref, read, pad, 5);
string ref2 = pad + ref;
string read2 = pad + read;
swaligner_.alignStrings(ref2, read2, ed, coord);
// match should start from pad string
if(coord.off() != 0) {
return false;
}
// no edits on pad string
if(ed.size() > 0 && ed[0].pos < pad.length()) {
return false;
}
size_t left = pad.length();
size_t right = left + read.length();
edits.clear();
edits.reserveExact(ed.size());
for(size_t i = 0; i < ed.size(); i++) {
if(ed[i].pos >= left && ed[i].pos <= right) {
edits.push_back(ed[i]);
edits.back().pos -= left;
}
}
max_matchlen = getMaxMatchLen(edits, read.length());
#ifdef DEBUGLOG
{
cerr << "After pad removed" << endl;
BTDnaString btread;
btread.install(read.c_str(), true);
Edit::print(cerr, edits); cerr << endl;
Edit::printQAlign(cerr, btread, edits);
}
#endif
return (max_matchlen >= rpt_len);
}
template<typename TStr>
void RepeatBuilder<TStr>::saveRepeats(const RepeatParameter &rp)
{
ios_base::openmode mode = ios_base::out;
if(rp.append_result) {
mode |= ios_base::app;
} else {
mode |= ios_base::trunc;
}
// Generate SNPs
size_t i = 0;
for(map<size_t, RB_Repeat*>::iterator it = repeat_map_.begin(); it != repeat_map_.end(); it++, i++) {
RB_Repeat& repeat = *(it->second);
if(!repeat.satisfy(rp))
continue;
// for each repeats
repeat.generateSNPs(rp, s_, i);
}
// save snp, consensus sequenuce, info
string snp_fname = filename_ + ".rep.snp";
string info_fname = filename_ + ".rep.info";
string hapl_fname = filename_ + ".rep.haplotype";
ofstream snp_fp(snp_fname.c_str(), mode);
ofstream info_fp(info_fname.c_str(), mode);
ofstream hapl_fp(hapl_fname.c_str(), mode);
const string repName = "rep" + to_string(rp.min_repeat_len) + "-" + to_string(rp.max_repeat_len);
i = 0;
TIndexOffU consensus_baseoff = 0;
TIndexOffU snp_id_base = 0;
TIndexOffU hapl_id_base = 0;
for(map<size_t, RB_Repeat*>::iterator it = repeat_map_.begin(); it != repeat_map_.end(); it++, i++) {
RB_Repeat& repeat = *(it->second);
if(!repeat.satisfy(rp))
continue;
// for each repeats
repeat.saveSNPs(snp_fp, i, snp_id_base);
saveAlleles(rp,
repName,
repeat,
info_fp,
hapl_fp,
i,
hapl_id_base,
consensus_baseoff);
consensus_baseoff += repeat.consensus().length();
}
snp_fp.close();
info_fp.close();
hapl_fp.close();
// save all consensus sequence
saveConsensus(rp, repName);
}
template<typename TStr>
void RepeatBuilder<TStr>::saveConsensus(const RepeatParameter &rp,
const string& repName) {
ios_base::openmode mode = ios_base::out;
if(rp.append_result) {
mode |= ios_base::app;
} else {
mode |= ios_base::trunc;
}
string fa_fname = filename_ + ".rep.fa";
ofstream fa_fp(fa_fname.c_str(), mode);
fa_fp << ">" << repName << endl;
size_t oskip = 0;
for(map<size_t, RB_Repeat*>::iterator it = repeat_map_.begin(); it != repeat_map_.end(); it++) {
RB_Repeat& repeat = *(it->second);
if(!repeat.satisfy(rp))
continue;
// for each repeats
const string& constr = repeat.consensus();
size_t si = 0;
size_t constr_len = constr.length();
while(si < constr_len) {
size_t out_len = std::min((size_t)(output_width - oskip), (size_t)(constr_len - si));
fa_fp << constr.substr(si, out_len);
if((oskip + out_len) == output_width) {
fa_fp << endl;
oskip = 0;
} else {
// last line
oskip = oskip + out_len;
}
si += out_len;
}
}
if(oskip) {
fa_fp << endl;
}
fa_fp.close();
}
template<typename TStr>
void RepeatBuilder<TStr>::saveAlleles(
const RepeatParameter& rp,
const string& repName,
RB_Repeat& repeat,
ofstream& fp,
ofstream& hapl_fp,
TIndexOffU grp_id,
TIndexOffU& hapl_id_base,
TIndexOffU baseoff)
{
const EList<SeedExt>& seeds = repeat.seeds();
EList<SeedHP> haplo_lists;
Range range(0, seeds.size());
int allele_id = 0;
for(size_t sb = range.first; sb < range.second;) {
size_t se = sb + 1;
for(; se < range.second; se++) {
if(!(SeedExt::isSameConsensus(seeds[sb], seeds[se])
&& SeedExt::isSameSNPs(seeds[sb], seeds[se])
&& seeds[sb].aligned == seeds[se].aligned)) {
break;
}
}
if(!seeds[sb].aligned) {
sb = se;
continue;
}
if(seeds[sb].getLength() < rp.min_repeat_len) {
sb = se;
continue;
}
// [sb, se) are same alleles
writeAllele(grp_id, allele_id, Range(sb, se),
repName, baseoff, seeds, fp);
generateHaploType(Range(sb, se), seeds, haplo_lists);
allele_id++;
sb = se;
}
// sort HaploType List by pos
haplo_lists.sort();
writeHaploType(haplo_lists, seeds, hapl_id_base, repName, hapl_fp);
}
template<typename TStr>
void RepeatBuilder<TStr>::generateHaploType(Range range, const EList<SeedExt> &seeds, EList<SeedHP> &haplo_list)
{
const EList<SeedSNP *>& snps = seeds[range.first].snps;
if(snps.size() == 0)
return;
// right-most position
TIndexOffU max_right_pos = seeds[range.first].consensus_pos.second - 1;
#ifndef NDEBUG
for(size_t i = 0; i < snps.size(); i++) {
assert_leq(snps[i]->pos, max_right_pos);
}
#endif
// create haplotypes of at least 16 bp long to prevent combinations of SNPs
// break a list of SNPs into several haplotypes if a SNP is far from the next SNP in the list
const TIndexOffU min_ht_len = 16;
size_t eb = 0, ee = 1;
while(ee < snps.size() + 1) {
if(ee == snps.size() ||
snps[eb]->pos + (min_ht_len << 1) < snps[ee]->pos) {
TIndexOffU left_pos = snps[eb]->pos;
TIndexOffU right_pos = snps[ee-1]->pos;
if(snps[ee-1]->type == EDIT_TYPE_READ_GAP) {
right_pos += snps[ee-1]->len;
}
if(left_pos + min_ht_len - 1 > right_pos) {
right_pos = left_pos + min_ht_len - 1;
}
right_pos = min<TIndexOffU>(max_right_pos, right_pos);
assert_leq(left_pos, right_pos);
SeedHP seedHP;
seedHP.range.first = left_pos;
seedHP.range.second = right_pos;
for(size_t i = eb; i < ee; i++) {
string snp_ids = "rps" + to_string(snps[i]->id);
seedHP.snpIDs.push_back(snp_ids);
}
// Add to haplo_list
bool found = false;
for(size_t i = 0; i < haplo_list.size(); i++) {
if(haplo_list[i] == seedHP) {
found = true;
break;
}
}
if(!found) {
// Add
haplo_list.push_back(seedHP);
}
eb = ee;
}
ee++;
}
}
void RB_SubSA::init(TIndexOffU sa_size,
TIndexOffU seed_len,
TIndexOffU seed_count)
{
assert_gt(sa_size, 1);
sa_size_ = sa_size;
assert_gt(seed_len, 0);
seed_len_ = seed_len;
assert_gt(seed_count, 1);
seed_count_ = seed_count;
temp_suffixes_.clear();
repeat_list_.clear();
repeat_index_.clear();
}
template<typename TStr>
void RB_SubSA::push_back(const TStr& s,
CoordHelper& coordHelper,
TIndexOffU saElt,
bool lastInput)
{
if(saElt + seed_len() <= coordHelper.getEnd(saElt)) {
if(seed_count_ == 1) {
repeat_list_.push_back(saElt);
} else {
assert_gt(seed_count_, 1);
if(temp_suffixes_.empty()) {
temp_suffixes_.push_back(saElt);
return;
}
TIndexOffU prev_saElt = temp_suffixes_.back();
// calculate common prefix length between two text.
// text1 is started from prev_saElt and text2 is started from saElt
bool same = isSameSequenceUpto(s,
coordHelper,
prev_saElt,
saElt,
seed_len_);
if(same) {
temp_suffixes_.push_back(saElt);
}
if(!same || lastInput) {
if(temp_suffixes_.size() >= seed_count_) {
repeat_index_.push_back(repeat_list_.size());
temp_suffixes_.sort();
for(size_t pi = 0; pi < temp_suffixes_.size(); pi++) {
repeat_list_.push_back(temp_suffixes_[pi]);
}
}
temp_suffixes_.clear();
if(!lastInput) {
temp_suffixes_.push_back(saElt);
} else {
temp_suffixes_.nullify();
}
}
}
}
if(lastInput) {
size_t bit = sizeof(uint32_t) * 8;
size_t num = (repeat_list_.size() + bit - 1) / bit;
done_.resizeExact(num);
done_.fillZero();
#ifndef NDEBUG
string prev_seq = "";
for(size_t i = 0; i < repeat_index_.size(); i++) {
TIndexOffU saBegin = repeat_index_[i];
TIndexOffU saEnd = (i + 1 < repeat_index_.size() ? repeat_index_[i+1] : repeat_list_.size());
string seq = getString(s, repeat_list_[saBegin], seed_len());
for(size_t j = saBegin + 1; j < saEnd; j++) {
string tmp_seq = getString(s, repeat_list_[j], seed_len());
assert_eq(seq, tmp_seq);
}
if(prev_seq != "" ) {
assert_lt(prev_seq, seq);
}
prev_seq = seq;
}
#endif
}
}
template<typename TStr>
Range RB_SubSA::find(const TStr& s,
const string& seq) const
{
TIndexOffU i = find_repeat_idx(s, seq);
if(i >= repeat_index_.size())
return Range(0, 0);
Range range(repeat_index_[i],
i + 1 < repeat_index_.size() ? repeat_index_[i+1] : repeat_list_.size());
return range;
}
template<typename TStr>
TIndexOffU RB_SubSA::find_repeat_idx(const TStr& s,
const string& seq) const
{
assert_eq(seq.length(), seed_len_);
string temp;
size_t l = 0, r = repeat_index_.size();
while(l < r) {
size_t m = (r + l) >> 1;
TIndexOffU saElt = repeat_list_[repeat_index_[m]];
getString(s, saElt, seed_len_, temp);
if(seq == temp) {
return m;
} else if(seq < temp) {
r = m;
} else {
assert(seq > temp);
l = m + 1;
}
}
return repeat_index_.size();
}
void RB_SubSA::setDone(TIndexOffU off, TIndexOffU len)
{
assert_leq(off + len, repeat_list_.size());
const TIndexOffU bit = sizeof(uint32_t) * 8;
for(TIndexOffU i = off; i < off + len; i++) {
TIndexOffU quotient = i / bit;
TIndexOffU remainder = i % bit;
assert_lt(quotient, done_.size());
uint32_t num = done_[quotient];
num = num | (1 << remainder);
done_[quotient] = num;
}
}
bool RB_SubSA::isDone(TIndexOffU off, TIndexOffU len) const
{
assert_leq(off + len, repeat_list_.size());
const TIndexOffU bit = sizeof(uint32_t) * 8;
for(TIndexOffU i = off; i < off + len; i++) {
TIndexOffU quotient = i / bit;
TIndexOffU remainder = i % bit;
assert_lt(quotient, done_.size());
uint32_t num = done_[quotient];
num = num & (1 << remainder);
if(num == 0)
return false;
}
return true;
}
template<typename TStr>
void RB_SubSA::buildRepeatBase(const TStr& s,
CoordHelper& coordHelper,
const size_t max_len,
EList<RB_RepeatBase>& repeatBases)
{
if(repeat_index_.empty())
return;
done_.fillZero();
EList<uint8_t> senseDominant;
senseDominant.resizeExact(repeat_index_.size());
senseDominant.fillZero();
EList<size_t> repeatStack;
EList<pair<TIndexOffU, TIndexOffU> > size_table;
size_table.reserveExact(repeat_index_.size() / 2 + 1);
for(size_t i = 0; i < repeat_index_.size(); i++) {
TIndexOffU begin = repeat_index_[i];
TIndexOffU end = (i + 1 < repeat_index_.size() ? repeat_index_[i+1] : repeat_list_.size());
assert_lt(begin, end);
EList<TIndexOffU> positions; positions.reserveExact(end - begin);
for(size_t j = begin; j < end; j++) positions.push_back(repeat_list_[j]);
if(!isSenseDominant(coordHelper, positions, seed_len_) && !threeN)
continue;
senseDominant[i] = 1;
size_table.expand();
size_table.back().first = end - begin;
size_table.back().second = i;
}
size_table.sort();
size_t bundle_count = 0;
string tmp_str;
EList<Range> tmp_ranges; tmp_ranges.resizeExact(4);
EList<pair<size_t, size_t> > tmp_sort_ranges; tmp_sort_ranges.resizeExact(4);
for(int64_t i = (int64_t)size_table.size() - 1; i >= 0; i--) {
TIndexOffU idx = size_table[i].second;
ASSERT_ONLY(TIndexOffU num = size_table[i].first);
assert_lt(idx, repeat_index_.size());
#ifndef NDEBUG
if(idx + 1 < repeat_index_.size()) {
assert_eq(repeat_index_[idx] + num, repeat_index_[idx+1]);
} else {
assert_eq(repeat_index_[idx] + num, repeat_list_.size());
}
#endif
TIndexOffU saBegin = repeat_index_[idx];
if(isDone(saBegin)) {
assert(isDone(saBegin, num));
continue;
}
ASSERT_ONLY(size_t rb_done = 0);
assert_lt(saBegin, repeat_list_.size());
repeatStack.push_back(idx);
while(!repeatStack.empty()) {
TIndexOffU idx = repeatStack.back();
assert(senseDominant[idx]);
repeatStack.pop_back();
assert_lt(idx, repeat_index_.size());
TIndexOffU saBegin = repeat_index_[idx];
TIndexOffU saEnd = (idx + 1 < repeat_index_.size() ? repeat_index_[idx + 1] : repeat_list_.size());
if(isDone(saBegin)) {
assert(isDone(saBegin, saEnd - saBegin));
continue;
}
TIndexOffU saElt = repeat_list_[saBegin];
size_t ri = repeatBases.size();
repeatBases.expand();
repeatBases.back().seq = getString(s, saElt, seed_len_);
repeatBases.back().nodes.clear();
repeatBases.back().nodes.push_back(idx);
ASSERT_ONLY(rb_done++);
setDone(saBegin, saEnd - saBegin);
bool left = true;
while(repeatBases[ri].seq.length() <= max_len) {
if(left) {
tmp_str = "N";
tmp_str += repeatBases[ri].seq.substr(0, seed_len_ - 1);
} else {
tmp_str = repeatBases[ri].seq.substr(repeatBases[ri].seq.length() - seed_len_ + 1, seed_len_ - 1);
tmp_str.push_back('N');
}
assert_eq(tmp_str.length(), seed_len_);
for(size_t c = 0; c < 4; c++) {
if(left) tmp_str[0] = "ACGT"[c];
else tmp_str.back() = "ACGT"[c];
assert_eq(tmp_str.length(), seed_len_);
TIndexOffU idx = find_repeat_idx(s, tmp_str);
size_t num = 0;
if(idx < repeat_index_.size()) {
if(idx + 1 < repeat_index_.size()) {
num = repeat_index_[idx+1] - repeat_index_[idx];
} else {
num = repeat_list_.size() - repeat_index_[idx];
}
}
tmp_ranges[c].first = idx;
tmp_ranges[c].second = num;
assert(num == 0 || num >= seed_count_);
tmp_sort_ranges[c].first = num;
tmp_sort_ranges[c].second = c;
if(idx == repeat_index_.size() ||
isDone(repeat_index_[idx]) ||
!senseDominant[idx]) {
#ifndef NDEBUG
if(idx < repeat_index_.size()) {
assert(isDone(repeat_index_[idx], num) ||
!senseDominant[idx]);
}
#endif
tmp_sort_ranges[c].first = 0;
tmp_ranges[c].second = 0;
}
}
tmp_sort_ranges.sort();
if(tmp_sort_ranges[3].first < seed_count_) {
if(left) {
left = false;
continue;
} else {
break;
}
}
for(size_t cc = 0; cc < 3; cc++) {
assert_leq(tmp_sort_ranges[cc].first, tmp_sort_ranges[cc+1].first);
if(tmp_sort_ranges[cc].first < seed_count_)
continue;
size_t c = tmp_sort_ranges[cc].second;
repeatStack.push_back(tmp_ranges[c].first);
}
size_t c = tmp_sort_ranges[3].second;
if(repeatBases[ri].seq.length() >= max_len) {
assert_eq(repeatBases[ri].seq.length(), max_len);
TIndexOffU idx = tmp_ranges[c].first;
assert(!isDone(repeat_index_[idx]));
repeatStack.push_back(idx);
if(left) {
left = false;
continue;
} else {
break;
}
} else {
TIndexOffU idx = tmp_ranges[c].first;
TIndexOffU num = tmp_ranges[c].second;
setDone(repeat_index_[idx], num);
if(left) {
repeatBases[ri].seq.insert(0, 1, "ACGT"[c]);
repeatBases[ri].nodes.insert(idx, 0);
} else {
repeatBases[ri].seq.push_back("ACGT"[c]);
repeatBases[ri].nodes.push_back(idx);
}
}
}
}
assert_gt(rb_done, 0);
bundle_count++;
}
cerr << "Bundle count: " << bundle_count << endl;
#ifndef NDEBUG
{
set<TIndexOffU> idx_set;
for(size_t i = 0; i < repeatBases.size(); i++) {
const RB_RepeatBase& repeatBase = repeatBases[i];
for(size_t j = 0; j < repeatBase.nodes.size(); j++) {
assert(idx_set.find(repeatBase.nodes[j]) == idx_set.end());
idx_set.insert(repeatBase.nodes[j]);
}
}
}
#endif
#if 0
{
EList<pair<size_t, size_t> > kmer_table;
string query;
size_t total = 0, match = 0;
size_t interval = 1;
if(repeat_index_.size() >= 10000) {
interval = repeat_index_.size() / 10000;
}
EList<TIndexOffU> positions;
for(size_t i = 0; i < repeat_index_.size(); i += interval) {
TIndexOffU saElt_idx = repeat_index_[i];
TIndexOffU saElt_idx_end = (i + 1 < repeat_index_.size() ? repeat_index_[i+1] : repeat_list_.size());
positions.clear();
for(size_t j = saElt_idx; j < saElt_idx_end; j++) {
positions.push_back(repeat_list_[j]);
#ifndef NDEBUG
if(j > saElt_idx) {
TIndexOffU lcp_len = getLCP(s,
coordHelper,
positions[0],
positions.back(),
seed_len_);
assert_geq(lcp_len, seed_len_);
}
TIndexOffU saElt = repeat_list_[j];
TIndexOffU start = coordHelper.getStart(saElt);
TIndexOffU start2 = coordHelper.getStart(saElt + seed_len_ - 1);
assert_eq(start, start2);
#endif
}
TIndexOffU saElt = repeat_list_[saElt_idx];
size_t true_count = saElt_idx_end - saElt_idx;
getString(s, saElt, seed_len_, query);
total++;
size_t count = 0;
for(size_t r = 0; r < repeatBases.size(); r++) {
const RB_RepeatBase& repeat = repeatBases[r];
int pos = repeat.seq.find(query);
if(pos != string::npos) {
for(size_t j = 0; j < repeat.nodes.size(); j++) {
TIndexOffU _node = repeat.nodes[j];
TIndexOffU _saElt_idx = repeat_index_[_node];
TIndexOffU _saElt_idx_end = (_node + 1 < repeat_index_.size() ? repeat_index_[_node+1] : repeat_list_.size());
TIndexOffU _saElt = repeat_list_[_saElt_idx];
string seq = getString(s, _saElt, seed_len());
if(query == seq) {
count += (_saElt_idx_end - _saElt_idx);
}
}
}
}
size_t rc_count = 0;
string rc_query = reverse_complement(query);
for(size_t r = 0; r < repeatBases.size(); r++) {
const RB_RepeatBase& repeat = repeatBases[r];
int pos = repeat.seq.find(rc_query);
if(pos != string::npos) {
for(size_t j = 0; j < repeat.nodes.size(); j++) {
TIndexOffU _node = repeat.nodes[j];
TIndexOffU _saElt_idx = repeat_index_[_node];
TIndexOffU _saElt_idx_end = (_node + 1 < repeat_index_.size() ? repeat_index_[_node+1] : repeat_list_.size());
TIndexOffU _saElt = repeat_list_[_saElt_idx];
string seq = getString(s, _saElt, seed_len());
if(rc_query == seq) {
rc_count += (_saElt_idx_end - _saElt_idx);
}
}
}
}
if(count == true_count || rc_count == true_count) {
match++;
} else if(total - match <= 10) {
cerr << "query: " << query << endl;
cerr << "true count: " << true_count << endl;
cerr << "found count: " << count << endl;
cerr << "rc found count: " << rc_count << endl;
cerr << endl;
}
}
cerr << "RB_SubSA: sanity check: " << match << " passed (out of " << total << ")" << endl << endl;
}
#endif
}
#define write_fp(x) fp.write((const char *)&(x), sizeof((x)))
void RB_SubSA::writeFile(ofstream& fp)
{
write_fp(sa_size_);
write_fp(seed_len_);
write_fp(seed_count_);
size_t sz = temp_suffixes_.size();
write_fp(sz);
for(size_t i = 0; i < sz; i++) {
write_fp(temp_suffixes_[i]);
}
sz = repeat_index_.size();
write_fp(sz);
for(size_t i = 0; i < sz; i++) {
write_fp(repeat_index_[i]);
}
sz = done_.size();
write_fp(sz);
for(size_t i = 0; i < sz; i++) {
write_fp(done_[i]);
}
}
#define read_fp(x) fp.read((char *)&(x), sizeof((x)))
void RB_SubSA::readFile(ifstream& fp)
{
TIndexOffU val;
read_fp(val);
rt_assert_eq(val, sa_size_);
read_fp(val);
rt_assert_eq(val, seed_len_);
read_fp(val);
rt_assert_eq(val, seed_count_);
size_t sz;
read_fp(sz);
temp_suffixes_.resizeExact(sz);
for(size_t i = 0; i < sz; i++) {
read_fp(temp_suffixes_[i]);
}
size_t val_sz;
// repeat_index_
read_fp(val_sz);
repeat_index_.resizeExact(val_sz);
for(size_t i = 0; i < val_sz; i++) {
read_fp(repeat_index_[i]);
}
// done
read_fp(val_sz);
done_.resizeExact(val_sz);
for(size_t i = 0; i < val_sz; i++) {
read_fp(done_[i]);
}
}
/****************************/
template class RepeatBuilder<SString<char> >;
template void dump_tstr(const SString<char>& );
template bool compareRepeatCoordByJoinedOff(const RepeatCoord<TIndexOffU>& , const RepeatCoord<TIndexOffU>&);
Reference in New Issue View Git Blame Copy Permalink