Yaossg/hisat-3n
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
main
hisat-3n/hisat2_repeat.cpp

997 lines
36 KiB
C++
Raw Permalink Normal View History

initial commit
2025-01-18 13:09:52 +00:00
/*
* 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 <fstream>
#include <string>
#include <cassert>
#include <getopt.h>
#include <algorithm>
#include "assert_helpers.h"
#include "endian_swap.h"
#include "formats.h"
#include "sequence_io.h"
#include "tokenize.h"
#include "timer.h"
#include "ref_read.h"
#include "filebuf.h"
#include "reference.h"
#include "ds.h"
#include "gfm.h"
#include "aligner_sw.h"
#include "aligner_result.h"
#include "search_globals.h"
#include "scoring.h"
#include "mask.h"
#include "repeat_builder.h"
#include "utility_3n.h"
/**
* \file Driver for the bowtie-build indexing tool.
*/
#include <algorithm>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <vector>
MemoryTally gMemTally;
// Build parameters
int verbose;
static int sanityCheck;
static int format;
static TIndexOffU bmax;
static TIndexOffU bmaxMultSqrt;
static uint32_t bmaxDivN;
static int dcv;
static int noDc;
static int entireSA;
static int seed;
static int showVersion;
// GFM parameters
static int32_t lineRate;
static bool lineRate_provided;
static int32_t linesPerSide;
static int32_t offRate;
static int32_t ftabChars;
static int32_t localOffRate;
static int32_t localFtabChars;
static int bigEndian;
static bool autoMem;
static int nthreads; // number of pthreads operating concurrently
static string wrapper;
static TIndexOffU seed_length;
static TIndexOffU seed_count;
static TIndexOffU repeat_count;
static TIndexOffU min_repeat_length;
static TIndexOffU max_repeat_length;
static EList<pair<TIndexOffU, TIndexOffU> > repeat_length_pair;
static TIndexOffU max_repeat_edit;
static TIndexOffU max_repeat_matchlen;
static bool symmetric_extend;
static bool repeat_indel;
static bool forward_only;
static bool CGtoTG;
static string repeat_str1;
static string repeat_str2;
TIndexOffU max_seed_mm;
TIndexOffU max_seed_repeat;
TIndexOffU max_seed_extlen;
static bool save_sa;
static bool load_sa;
bool threeN = false;
char convertedFrom;
char convertedTo;
char convertedFromComplement;
char convertedToComplement;
bool base_change_entered;
static void resetOptions() {
verbose = true; // be talkative (default)
sanityCheck = 0; // do slow sanity checks
format = FASTA; // input sequence format
bmax = OFF_MASK; // max blockwise SA bucket size
bmaxMultSqrt = OFF_MASK; // same, as multplier of sqrt(n)
bmaxDivN = 4; // same, as divisor of n
dcv = 1024; // bwise SA difference-cover sample sz
noDc = 0; // disable difference-cover sample
entireSA = 0; // 1 = disable blockwise SA
seed = 0; // srandom seed
showVersion = 0; // just print version and quit?
// GFM parameters
lineRate = GFM<TIndexOffU>::default_lineRate_gfm;
lineRate_provided = false;
linesPerSide = 1; // 1 64-byte line on a side
offRate = 4; // sample 1 out of 16 SA elts
ftabChars = 10; // 10 chars in initial lookup table
localOffRate = 3;
localFtabChars = 6;
bigEndian = 0; // little endian
autoMem = true; // automatically adjust memory usage parameters
nthreads = 1;
seed_length = 50;
seed_count = 5;
min_repeat_length = 100;
max_repeat_length = numeric_limits<uint16_t>::max();
repeat_length_pair.clear();
repeat_count = 5;
max_repeat_edit = 10;
max_repeat_matchlen = min_repeat_length / 2; // half of repeat_length
repeat_indel = false;
symmetric_extend = true;
forward_only = false;
CGtoTG = false;
max_seed_mm = 5;
max_seed_repeat = 5;
max_seed_extlen = 25;
save_sa = false;
load_sa = false;
wrapper.clear();
threeN = false;
convertedFrom = 'C';
convertedTo = 'T';
convertedFromComplement = asc2dnacomp[convertedFrom];
convertedToComplement = asc2dnacomp[convertedTo];
base_change_entered = false;
}
// Argument constants for getopts
enum {
ARG_BMAX = 256,
ARG_BMAX_MULT,
ARG_BMAX_DIV,
ARG_DCV,
ARG_SEED,
ARG_CUTOFF,
ARG_PMAP,
ARG_NTOA,
ARG_USAGE,
ARG_REVERSE_EACH,
ARG_SA,
ARG_SEED_LENGTH,
ARG_SEED_COUNT,
ARG_MIN_REPEAT_LENGTH,
ARG_MAX_REPEAT_LENGTH,
ARG_REPEAT_LENGTH,
ARG_REPEAT_COUNT,
ARG_REPEAT_EDIT,
ARG_REPEAT_MATCHLEN,
ARG_REPEAT_INDEL,
ARG_WRAPPER,
ARG_ASYMMETRIC_EXTEND,
ARG_FORWARD_ONLY,
ARG_CGTOTG,
ARG_REPEAT_STR1,
ARG_REPEAT_STR2,
ARG_MAX_SEED_MM,
ARG_MAX_SEED_REPEAT,
ARG_MAX_SEED_EXTLEN,
ARG_SAVE_SA,
ARG_LOAD_SA,
ARG_3N,
ARG_BASE_CHANGE
};
/**
* Print a detailed usage message to the provided output stream.
*/
static void printUsage(ostream& out) {
out << "HISAT2 version " << string(HISAT2_VERSION).c_str() << " by Chanhee Park <parkchanhee@gmail.com> and Daehwan Kim <infphilo@gmail.com>" << endl;
string tool_name = "hisat2-repeat";
out << "Usage: " << tool_name << " [options]* <reference_in>" << endl
<< " reference_in comma-separated list of files with ref sequences" << endl
<< "Options:" << endl
<< " -c reference sequences given on cmd line (as" << endl
<< " <reference_in>)" << endl;
if(wrapper == "basic-0") {
out << " --large-index force generated index to be 'large', even if ref" << endl
<< " has fewer than 4 billion nucleotides" << endl;
}
out << " -a/--noauto disable automatic -p/--bmax/--dcv memory-fitting" << endl
<< " -p <int> number of threads" << endl
<< " --bmax <int> max bucket sz for blockwise suffix-array builder" << endl
<< " --bmaxdivn <int> max bucket sz as divisor of ref len (default: 4)" << endl
<< " --dcv <int> diff-cover period for blockwise (default: 1024)" << endl
<< " --nodc disable diff-cover (algorithm becomes quadratic)" << endl
<< " --seed-length <int> seed length (default: 50)" << endl
<< " --seed-count <int> seed count (default: 5)" << endl
<< " --min-repeat-length <int> minimum repeat length (default: 100)" << endl
<< " --max-repeat-length <int> maximum repeat length (default: 65535)" << endl
<< " --repeat-length <int>-<int>[,<int>-<int>] minimum-maximum repeat length pairs" << endl
<< " --repeat-count <int> minimum repeat count (default: 5)" << endl
<< " --repeat-edit <int> maximum repeat edit distance (default: 10)" << endl
<< " --repeat-matchlen <int>" << endl
<< " --repeat-indel" << endl
<< " --repeat-str1" << endl
<< " --repeat-str2" << endl
<< " --asymmetric-extend extend seeds asymmetrically" << endl
<< " --forward-only use forward strand only" << endl
<< " --CGtoTG change CG to TG" << endl
<< " --max-seed-mm <int>" << endl
<< " --max-seed-repeat <int>" << endl
<< " --max-seed-extlen <int>" << endl
<< " --save-sa" << endl
<< " --load-sa" << endl
<< " --3N make 3N repeat database" << endl
<< " --base-change <chr,chr> the converted nucleotide and converted to nucleotide (default:C,T)" << endl
<< " -q/--quiet disable verbose output (for debugging)" << endl
<< " -h/--help print detailed description of tool and its options" << endl
<< " --usage print this usage message" << endl
<< " --version print version information and quit" << endl
;
}
static const char *short_options = "qrap:h?nscfl:i:o:t:h:3C";
static struct option long_options[] = {
{(char*)"quiet", no_argument, 0, 'q'},
{(char*)"sanity", no_argument, 0, 's'},
{(char*)"threads", required_argument, 0, 'p'},
{(char*)"little", no_argument, &bigEndian, 0},
{(char*)"big", no_argument, &bigEndian, 1},
{(char*)"bmax", required_argument, 0, ARG_BMAX},
{(char*)"bmaxmultsqrt", required_argument, 0, ARG_BMAX_MULT},
{(char*)"bmaxdivn", required_argument, 0, ARG_BMAX_DIV},
{(char*)"dcv", required_argument, 0, ARG_DCV},
{(char*)"nodc", no_argument, &noDc, 1},
{(char*)"seed", required_argument, 0, ARG_SEED},
{(char*)"entiresa", no_argument, &entireSA, 1},
{(char*)"version", no_argument, &showVersion, 1},
{(char*)"noauto", no_argument, 0, 'a'},
{(char*)"noblocks", required_argument, 0, 'n'},
{(char*)"linerate", required_argument, 0, 'l'},
{(char*)"linesperside", required_argument, 0, 'i'},
{(char*)"usage", no_argument, 0, ARG_USAGE},
{(char*)"seed-length", required_argument, 0, ARG_SEED_LENGTH},
{(char*)"seed-count", required_argument, 0, ARG_SEED_COUNT},
{(char*)"min-repeat-length", required_argument, 0, ARG_MIN_REPEAT_LENGTH},
{(char*)"max-repeat-length", required_argument, 0, ARG_MAX_REPEAT_LENGTH},
{(char*)"repeat-length", required_argument, 0, ARG_REPEAT_LENGTH},
{(char*)"repeat-count", required_argument, 0, ARG_REPEAT_COUNT},
{(char*)"repeat-edit", required_argument, 0, ARG_REPEAT_EDIT},
{(char*)"repeat-matchlen",required_argument, 0, ARG_REPEAT_MATCHLEN},
{(char*)"repeat-indel", no_argument, 0, ARG_REPEAT_INDEL},
{(char*)"wrapper", required_argument, 0, ARG_WRAPPER},
{(char*)"asymmetric-extend", no_argument, 0, ARG_ASYMMETRIC_EXTEND},
{(char*)"forward-only", no_argument, 0, ARG_FORWARD_ONLY},
{(char*)"CGtoTG", no_argument, 0, ARG_CGTOTG},
{(char*)"repeat-str1", required_argument, 0, ARG_REPEAT_STR1},
{(char*)"repeat-str2", required_argument, 0, ARG_REPEAT_STR2},
{(char*)"max-seed-mm", required_argument, 0, ARG_MAX_SEED_MM},
{(char*)"max-seed-repeat",required_argument, 0, ARG_MAX_SEED_REPEAT},
{(char*)"max-seed-extlen",required_argument, 0, ARG_MAX_SEED_EXTLEN},
{(char*)"save-sa", no_argument, 0, ARG_SAVE_SA},
{(char*)"load-sa", no_argument, 0, ARG_LOAD_SA},
{(char*)"3N", no_argument, 0, ARG_3N},
{(char*)"base-change", required_argument, 0, ARG_BASE_CHANGE},
{(char*)0, 0, 0, 0} // terminator
};
/**
* Parse an int out of optarg and enforce that it be at least 'lower';
* if it is less than 'lower', then output the given error message and
* exit with an error and a usage message.
*/
template<typename T>
static T parseNumber(const char *str, T lower, const char *errmsg)
{
char *endPtr= NULL;
T t = (T)strtoll(str, &endPtr, 10);
if(endPtr != NULL) {
if(t < lower) {
cerr << errmsg << endl;
printUsage(cerr);
throw 1;
}
return t;
}
cerr << errmsg << endl;
printUsage(cerr);
throw 1;
return -1;
}
template<typename T>
static T parseNumber(T lower, const char *errmsg) {
return parseNumber(optarg, lower, errmsg);
}
static void parsePair(EList<pair<TIndexOffU, TIndexOffU> >& repeat_pair)
{
string tok;
istringstream ss(optarg);
while(getline(ss, tok, ',')) {
if(tok.empty()) {
continue;
}
TIndexOffU min_len, max_len;
size_t pos = tok.find('-');
if(pos == string::npos) {
// min
min_len = parseNumber<TIndexOffU>(tok.c_str(), 1, "min-repeat-length must be at least 1");
max_len = numeric_limits<uint16_t>::max();
} else if(pos == tok.length() - 1) {
// min-
min_len = parseNumber<TIndexOffU>(tok.substr(0, pos).c_str(), 1, "min-repeat-length must be at least 1");
max_len = numeric_limits<uint16_t>::max();
} else if(pos == 0) {
// -max
// not support?
min_len = 100;
max_len = parseNumber<TIndexOffU>(tok.substr(pos + 1).c_str(), 1, "max-repeat-length must be at least 1");
} else {
min_len = parseNumber<TIndexOffU>(tok.substr(0, pos).c_str(), 1, "min-repeat-length must be at least 1");
max_len = parseNumber<TIndexOffU>(tok.substr(pos + 1).c_str(), 1, "max-repeat-length must be at least 1");
}
if(min_len > max_len) {
printUsage(cerr);
throw 1;
}
if(max_len > numeric_limits<uint16_t>::max()) {
printUsage(cerr);
throw 1;
}
repeat_pair.push_back(pair<TIndexOffU, TIndexOffU>(min_len, max_len));
}
}
/**
* Read command-line arguments
*/
static void parseOptions(int argc, const char **argv) {
int option_index = 0;
int next_option;
bool saw_max_repeat_matchlen = false;
do {
next_option = getopt_long(
argc, const_cast<char**>(argv),
short_options, long_options, &option_index);
switch (next_option) {
case ARG_WRAPPER:
wrapper = optarg;
break;
case 'f': format = FASTA; break;
case 'c': format = CMDLINE; break;
//case 'p': packed = true; break;
case 'C':
cerr << "Error: -C specified but Bowtie 2 does not support colorspace input." << endl;
throw 1;
break;
case 'l':
lineRate = parseNumber<int>(3, "-l/--lineRate arg must be at least 3");
lineRate_provided = true;
break;
case 'i':
linesPerSide = parseNumber<int>(1, "-i/--linesPerSide arg must be at least 1");
break;
case 'o':
offRate = parseNumber<int>(0, "-o/--offRate arg must be at least 0");
break;
case 'n':
// all f-s is used to mean "not set", so put 'e' on end
bmax = 0xfffffffe;
break;
case 'h':
case ARG_USAGE:
printUsage(cout);
throw 0;
break;
case ARG_BMAX:
bmax = parseNumber<TIndexOffU>(1, "--bmax arg must be at least 1");
bmaxMultSqrt = OFF_MASK; // don't use multSqrt
bmaxDivN = 0xffffffff; // don't use multSqrt
break;
case ARG_BMAX_MULT:
bmaxMultSqrt = parseNumber<TIndexOffU>(1, "--bmaxmultsqrt arg must be at least 1");
bmax = OFF_MASK; // don't use bmax
bmaxDivN = 0xffffffff; // don't use multSqrt
break;
case ARG_BMAX_DIV:
bmaxDivN = parseNumber<uint32_t>(1, "--bmaxdivn arg must be at least 1");
bmax = OFF_MASK; // don't use bmax
bmaxMultSqrt = OFF_MASK; // don't use multSqrt
break;
case ARG_DCV:
dcv = parseNumber<int>(3, "--dcv arg must be at least 3");
break;
case ARG_SEED:
seed = parseNumber<int>(0, "--seed arg must be at least 0");
break;
case ARG_SEED_LENGTH:
seed_length = parseNumber<TIndexOffU>(1, "--seed-length arg must be at least 1");
break;
case ARG_SEED_COUNT:
seed_count = parseNumber<TIndexOffU>(2, "--repeat-count arg must be at least 2");
break;
case ARG_MIN_REPEAT_LENGTH:
min_repeat_length = parseNumber<TIndexOffU>(1, "--min-repeat-length arg must be at least 1");
break;
case ARG_MAX_REPEAT_LENGTH:
max_repeat_length = parseNumber<TIndexOffU>(1, "--max-repeat-length arg must be at least 1");
break;
case ARG_REPEAT_LENGTH:
parsePair(repeat_length_pair);
break;
case ARG_REPEAT_COUNT:
repeat_count = parseNumber<TIndexOffU>(2, "--repeat-count arg must be at least 2");
break;
case ARG_REPEAT_EDIT:
max_repeat_edit = parseNumber<TIndexOffU>(0, "--repeat-edit arg must be at least 0");
break;
case ARG_REPEAT_MATCHLEN:
max_repeat_matchlen = parseNumber<TIndexOffU>(0, "--repeat-matchlen arg must be at least 0");
saw_max_repeat_matchlen = true;
break;
case ARG_REPEAT_INDEL:
repeat_indel = true;
break;
case ARG_ASYMMETRIC_EXTEND:
symmetric_extend = false;
break;
case ARG_FORWARD_ONLY:
forward_only = true;
break;
case ARG_CGTOTG:
CGtoTG = true;
break;
case ARG_REPEAT_STR1:
repeat_str1 = optarg;
break;
case ARG_REPEAT_STR2:
repeat_str2 = optarg;
break;
case ARG_MAX_SEED_MM:
max_seed_mm = parseNumber<TIndexOffU>(1, "--max_seed_mm arg must be at least 1");
break;
case ARG_MAX_SEED_REPEAT:
max_seed_repeat = parseNumber<TIndexOffU>(5, "--max_seed_repeat arg must be at least 5");
break;
case ARG_MAX_SEED_EXTLEN:
max_seed_extlen = parseNumber<TIndexOffU>(0, "--max_seed_extlen arg must be at least 0");
break;
case ARG_SAVE_SA:
save_sa = true;
break;
case ARG_LOAD_SA:
load_sa = true;
break;
case ARG_3N: {
threeN = true;
break;
}
case ARG_BASE_CHANGE: {
EList<string> args;
tokenize(optarg, ",", args);
if(args.size() != 2) {
cerr << "Error: expected 2 comma-separated "
<< "arguments to --base-change option, got " << args.size() << endl;
throw 1;
}
getConversion(args[0][0], args[1][0], convertedFrom, convertedTo);
string s = "ACGT";
if ((s.find(convertedFrom) == std::string::npos) || (s.find(convertedTo) == std::string::npos)) {
cerr << "Please enter the nucleotide in 'ACGT' for --base-change option." << endl;
throw 1;
}
if (convertedFrom == convertedTo) {
cerr << "Please enter two different base for --base-change option. If you wish to build the repeat database without nucleotide conversion, please do not use --base-change and --3N options." << endl;
throw 1;
}
base_change_entered = true;
}
case 'a': autoMem = false; break;
case 'q': verbose = false; break;
case 's': sanityCheck = true; break;
case 'p':
nthreads = parseNumber<int>(1, "-p arg must be at least 1");
break;
case -1: /* Done with options. */
break;
case 0:
if (long_options[option_index].flag != 0)
break;
default:
printUsage(cerr);
throw 1;
}
} while(next_option != -1);
if(bmax < 40) {
cerr << "Warning: specified bmax is very small (" << bmax << "). This can lead to" << endl
<< "extremely slow performance and memory exhaustion. Perhaps you meant to specify" << endl
<< "a small --bmaxdivn?" << endl;
}
if(!saw_max_repeat_matchlen) {
max_repeat_matchlen = min_repeat_length / 2;
}
}
extern void initializeCntLut();
extern void initializeCntBit();
ConvertMatrix3N baseChange;
/**
* Drive the index construction process and optionally sanity-check the
* result.
*/
template<typename TStr>
static void driver(
const string& infile,
EList<string>& infiles,
const string& outfile,
bool packed,
bool forward_only,
bool CGtoTG)
{
initializeCntLut();
initializeCntBit();
EList<FileBuf*> is(MISC_CAT);
bool bisulfite = false;
bool nsToAs = false;
RefReadInParams refparams(false, false /* reverse */, nsToAs, bisulfite);
assert_gt(infiles.size(), 0);
if(format == CMDLINE) {
// Adapt sequence strings to stringstreams open for input
stringstream *ss = new stringstream();
for(size_t i = 0; i < infiles.size(); i++) {
(*ss) << ">" << i << endl << infiles[i].c_str() << endl;
}
FileBuf *fb = new FileBuf(ss);
assert(fb != NULL);
assert(!fb->eof());
assert(fb->get() == '>');
ASSERT_ONLY(fb->reset());
assert(!fb->eof());
is.push_back(fb);
} else {
// Adapt sequence files to ifstreams
for(size_t i = 0; i < infiles.size(); i++) {
FILE *f = fopen(infiles[i].c_str(), "r");
if (f == NULL) {
cerr << "Error: could not open "<< infiles[i].c_str() << endl;
throw 1;
}
FileBuf *fb = new FileBuf(f);
assert(fb != NULL);
if(fb->peek() == -1 || fb->eof()) {
cerr << "Warning: Empty fasta file: '" << infile.c_str() << "'" << endl;
continue;
}
assert(!fb->eof());
assert(fb->get() == '>');
ASSERT_ONLY(fb->reset());
assert(!fb->eof());
is.push_back(fb);
}
}
if(is.empty()) {
cerr << "Warning: All fasta inputs were empty" << endl;
throw 1;
}
// Vector for the ordered list of "records" comprising the input
// sequences. A record represents a stretch of unambiguous
// characters in one of the input sequences.
EList<RefRecord> szs;
EList<string> ref_names;
std::pair<size_t, size_t> sztot;
{
if(verbose) cerr << "Reading reference sizes" << endl;
Timer _t(cerr, " Time reading reference sizes: ", verbose);
sztot = BitPairReference::szsFromFasta(is, "", bigEndian, refparams, szs, sanityCheck, &ref_names);
}
assert_gt(sztot.first, 0);
assert_gt(sztot.second, 0);
assert_gt(szs.size(), 0);
// Compose text strings into single string
cerr << "Calculating joined length" << endl;
TIndexOffU jlen = 0;
for(unsigned int i = 0; i < szs.size(); i++) {
jlen += (TIndexOffU)szs[i].len;
}
// assert_geq(jlen, sztot);
cerr << " Joined length: " << jlen << endl;
TStr s;
{
bool both_strand = forward_only ? false : true;
cerr << "Reserving space for joined string" << endl;
cerr << "Joining reference sequences" << endl;
Timer timer(cerr, " Time to join reference sequences: ", verbose);
GFM<TIndexOffU>::join<TStr>(
is,
szs,
(TIndexOffU) sztot.first,
refparams,
seed,
s,
both_strand, // include reverse complemented sequence
CGtoTG); //Change CG to TG
}
TStr sOriginal;
if (threeN) {
baseChange.restoreNormal();
bool both_strand = forward_only ? false : true;
for (int i = 0; i < is.size(); i++) {
is[i]->reset();
}
GFM<TIndexOffU>::join<TStr>(
is,
szs,
(TIndexOffU) sztot.first,
refparams,
seed,
sOriginal,
both_strand, // include reverse complemented sequence
CGtoTG); //Change CG to TG
baseChange.restoreConversion();
long long int guessLen = s.length() / 2;
for (TIndexOffU i = 0; i < guessLen; i++) {
int nt = sOriginal[guessLen - i - 1];
assert_range(0, 3, nt);
s[guessLen + i] = dnacomp[nt];
}
} else {
sOriginal = s;
}
// Successfully obtained joined reference string
#ifndef NDEBUG
if(forward_only) {
assert_geq(s.length(), jlen);
} else {
assert_geq(s.length(), jlen << 1);
}
#endif
BitPackedArray suffix_array;
bool sa_file_exist = false;
string sa_fname = outfile + ".rep.sa";
if(load_sa) {
ifstream fp(sa_fname, std::ifstream::binary);
sa_file_exist = fp.is_open();
}
if(load_sa && sa_file_exist) {
cerr << "Load SA from " << sa_fname << endl;
suffix_array.readFile(sa_fname);
} else {
suffix_array.init(s.length() + 1);
if(bmax != (TIndexOffU) OFF_MASK) {
// VMSG_NL("bmax according to bmax setting: " << bmax);
} else if(bmaxDivN != (uint32_t) OFF_MASK) {
bmax = max<uint32_t>(jlen / (bmaxDivN * nthreads), 1);
// VMSG_NL("bmax according to bmaxDivN setting: " << bmax);
} else {
bmax = (uint32_t) sqrt(s.length());
// VMSG_NL("bmax defaulted to: " << bmax);
}
int iter = 0;
bool first = true;
bool passMemExc = false, sanity = false;
// Look for bmax/dcv parameters that work.
while(true) {
if(!first && bmax < 40 && passMemExc) {
cerr << "Could not find appropriate bmax/dcv settings for building this index." << endl;
cerr << "Please try indexing this reference on a computer with more memory." << endl;
if (sizeof(void *) == 4) {
cerr << "If this computer has more than 4 GB of memory, try using a 64-bit executable;" << endl
<< "this executable is 32-bit." << endl;
}
throw 1;
}
if(dcv > 4096) dcv = 4096;
if((iter % 6) == 5 && dcv < 4096 && dcv != 0) {
dcv <<= 1; // double difference-cover period
} else {
bmax -= (bmax >> 2); // reduce by 25%
}
iter++;
suffix_array.reset();
try {
cerr << "Using parameters --bmax " << bmax << endl;
if(dcv == 0) {
cerr << " and *no difference cover*" << endl;
} else {
cerr << " --dcv " << dcv << endl;
}
{
cerr << " Doing ahead-of-time memory usage test" << endl;
// Make a quick-and-dirty attempt to force a bad_alloc iff
// we would have thrown one eventually as part of
// constructing the DifferenceCoverSample
dcv <<= 1;
TIndexOffU sz = (TIndexOffU) DifferenceCoverSample<TStr>::simulateAllocs(s, dcv >> 1);
if(nthreads > 1) sz *= (nthreads + 1);
AutoArray<uint8_t> tmp(sz, EBWT_CAT);
dcv >>= 1;
// Likewise with the KarkkainenBlockwiseSA
sz = (TIndexOffU) KarkkainenBlockwiseSA<TStr>::simulateAllocs(s, bmax);
AutoArray<uint8_t> tmp2(sz, EBWT_CAT);
// Grab another 20 MB out of caution
AutoArray<uint32_t> extra(20 * 1024 * 1024, EBWT_CAT);
// If we made it here without throwing bad_alloc, then we
// passed the memory-usage stress test
cerr << " Passed! Constructing with these parameters: --bmax " << bmax << " --dcv " << dcv
<< endl;
cerr << "" << endl;
}
cerr << "Constructing suffix-array element generator" << endl;
KarkkainenBlockwiseSA<TStr> bsa(s,
bmax,
nthreads,
dcv,
seed,
sanity,
passMemExc,
false /* verbose */,
outfile);
assert(bsa.suffixItrIsReset());
assert_eq(bsa.size(), s.length() + 1);
TIndexOffU count = 0;
while(count < s.length() + 1) {
TIndexOffU saElt = bsa.nextSuffix();
count++;
if(count && (count % 10000000 == 0)) {
cerr << "SA count " << count << endl;
}
if(saElt == s.length()) {
assert_eq(count, s.length() + 1);
break;
}
suffix_array.pushBack(saElt);
}
break;
} catch (bad_alloc &e) {
if(passMemExc) {
cerr << " Ran out of memory; automatically trying more memory-economical parameters." << endl;
} else {
cerr << "Out of memory while constructing suffix array. Please try using a smaller" << endl
<< "number of blocks by specifying a smaller --bmax or a larger --bmaxdivn" << endl;
throw 1;
}
}
first = false;
}
if(save_sa) {
suffix_array.writeFile(sa_fname);
}
}
cerr << "suffix_array: " << endl;
suffix_array.dump();
// Build Repeats
{
if(repeat_length_pair.empty()) {
repeat_length_pair.push_back(pair<TIndexOffU, TIndexOffU>(min_repeat_length, max_repeat_length));
}
RepeatParameter rp;
rp.seed_len = seed_length;
rp.seed_count = seed_count;
rp.seed_mm = max_seed_mm;
rp.max_edit = max_repeat_edit;
rp.symmetric_extend = symmetric_extend;
rp.extend_unit_len = max_seed_extlen;
for(size_t i = 0; i < repeat_length_pair.size(); i++) {
rp.min_repeat_len = repeat_length_pair[i].first;
rp.max_repeat_len = repeat_length_pair[i].second;
rp.repeat_count = repeat_count;
rp.append_result = (i != 0);
RepeatBuilder<TStr> repeatBuilder(s,
sOriginal,
szs,
ref_names,
forward_only,
outfile);
cerr << "RepeatBuilder: " << outfile << " " << rp.min_repeat_len << "-" << rp.max_repeat_len << endl;
{
Timer _t(cerr, " Time reading suffix array: ", verbose);
repeatBuilder.readSA(rp, suffix_array);
}
repeatBuilder.build(rp);
repeatBuilder.saveFile(rp);
}
}
}
static const char *argv0 = NULL;
extern "C" {
/**
* main function. Parses command-line arguments.
*/
int hisat2_repeat(int argc, const char **argv) {
string outfile;
try {
// Reset all global state, including getopt state
opterr = optind = 1;
resetOptions();
string infile;
EList<string> infiles(MISC_CAT);
parseOptions(argc, argv);
argv0 = argv[0];
if(showVersion) {
cout << argv0 << " version " << string(HISAT2_VERSION).c_str() << endl;
if(sizeof(void*) == 4) {
cout << "32-bit" << endl;
} else if(sizeof(void*) == 8) {
cout << "64-bit" << endl;
} else {
cout << "Neither 32- nor 64-bit: sizeof(void*) = " << sizeof(void*) << endl;
}
cout << "Built on " << BUILD_HOST << endl;
cout << BUILD_TIME << endl;
cout << "Compiler: " << COMPILER_VERSION << endl;
cout << "Options: " << COMPILER_OPTIONS << endl;
cout << "Sizeof {int, long, long long, void*, size_t, off_t}: {"
<< sizeof(int)
<< ", " << sizeof(long) << ", " << sizeof(long long)
<< ", " << sizeof(void *) << ", " << sizeof(size_t)
<< ", " << sizeof(off_t) << "}" << endl;
return 0;
}
if (!threeN && base_change_entered) {
cerr << "To build hisat-3n repeat database, please add argument --3N. To build the hisat2 repeat database, please remove the argument --base-change." << endl;
printUsage(cerr);
throw 1;
}
if (threeN) {
convertedFromComplement = asc2dnacomp[convertedFrom];
convertedToComplement = asc2dnacomp[convertedTo];
}
// Get input filename
if(optind >= argc) {
cerr << "No input sequence or sequence file specified!" << endl;
printUsage(cerr);
return 1;
}
infile = argv[optind++];
if(optind >= argc) {
cerr << "No output file specified!" << endl;
printUsage(cerr);
return 1;
}
outfile = argv[optind++];
tokenize(infile, ",", infiles);
if(infiles.size() < 1) {
cerr << "Tokenized input file list was empty!" << endl;
printUsage(cerr);
return 1;
}
// Optionally summarize
if(verbose) {
cerr << "Settings:" << endl;
// << " Output files: \"" << outfile.c_str() << ".*." << gfm_ext << "\"" << endl;
cerr << " Endianness: " << (bigEndian? "big":"little") << endl
<< " Actual local endianness: " << (currentlyBigEndian()? "big":"little") << endl
<< " Sanity checking: " << (sanityCheck? "enabled":"disabled") << endl;
#ifdef NDEBUG
cerr << " Assertions: disabled" << endl;
#else
cerr << " Assertions: enabled" << endl;
#endif
cerr << " Random seed: " << seed << endl;
cerr << " Sizeofs: void*:" << sizeof(void*) << ", int:" << sizeof(int) << ", long:" << sizeof(long) << ", size_t:" << sizeof(size_t) << endl;
cerr << "Input files DNA, " << file_format_names[format].c_str() << ":" << endl;
for(size_t i = 0; i < infiles.size(); i++) {
cerr << " " << infiles[i].c_str() << endl;
}
}
// Seed random number generator
srand(seed);
{
Timer timer(cerr, "Total time for call to driver() for forward index: ", verbose);
try {
int nloop = threeN ? 2 : 1; // if threeN == true, nloop = 2. else one loop
for (int i = 0; i < nloop; i++) {
string tag = "";
if (threeN) {
tag += ".3n.";
if (i == 0) {
tag += convertedFrom;
tag += convertedTo;
baseChange.convert(convertedFrom, convertedTo);
} else {
tag += convertedFromComplement;
tag += convertedToComplement;
baseChange.convert(convertedFromComplement, convertedToComplement);
}
string indexFilename = outfile + tag + ".rep.fa";
if (fileExist(indexFilename)) {
cerr << "*** Find repeat database for " << outfile + tag << "，skip this repeat database building process." << endl;
cerr << " To re-build your hisat-3n repeat database, please delete the old index manually before running hisat2-repeat or hisat-3n-build." << endl;
continue;
}
}
driver<SString<char> >(infile, infiles, outfile + tag, false, forward_only, CGtoTG);
}
} catch(bad_alloc& e) {
if(autoMem) {
cerr << "Switching to a packed string representation." << endl;
} else {
throw e;
}
}
}
return 0;
} catch(std::exception& e) {
cerr << "Error: Encountered exception: '" << e.what() << "'" << endl;
cerr << "Command: ";
for(int i = 0; i < argc; i++) cerr << argv[i] << " ";
cerr << endl;
return 1;
} catch(int e) {
if(e != 0) {
cerr << "Error: Encountered internal HISAT2 exception (#" << e << ")" << endl;
cerr << "Command: ";
for(int i = 0; i < argc; i++) cerr << argv[i] << " ";
cerr << endl;
}
return e;
}
}
}
Reference in New Issue Copy Permalink