zzh
/
minimap2
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

添加了一些统计运行时间的变量，测试性能

This commit is contained in:
zzh 2023-08-03 03:12:38 +08:00
parent 1641a66d15
commit 023c18ba32
10 changed files with 988 additions and 456 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
.vscode/extensions.json vendored 100644
View File

@ -0,0 +1,5 @@
{
"recommendations": [
"mhutchie.git-graph"
]
}

8
.vscode/launch.json vendored
View File

@ -10,12 +10,16 @@
"type": "cppdbg",
"request": "launch",
"program": "${workspaceRoot}/minimap2",
"args": ["-x",
"args": [
"-x",
"ava-ont",
"-t",
"1",
"/public/home/zzh/work/3gseq/TGM-2021YFF/Acinetobacter_pittii.fastq",
"/public/home/zzh/work/3gseq/TGM-2021YFF/Acinetobacter_pittii.fastq"],
"/public/home/zzh/work/3gseq/TGM-2021YFF/Acinetobacter_pittii.fastq",
"-o",
"reads.paf"
],
"cwd": "${workspaceFolder}", // 当前工作路径：当前文件所在的工作空间
}
]

6
.vscode/settings.json vendored 100644
View File

@ -0,0 +1,6 @@
{
"files.associations": {
"minimap.h": "c",
"time.h": "c"
}
}

3
Makefile
View File

@ -1,4 +1,5 @@
CFLAGS= -g -Wall -O2 -Wc++-compat #-Wextra
CFLAGS= -g -Wall -Wc++-compat #-Wextra
#CFLAGS= -g -Wall -O2 -Wc++-compat #-Wextra
CPPFLAGS= -DHAVE_KALLOC
INCLUDES=
OBJS= kthread.o kalloc.o misc.o bseq.o sketch.o sdust.o options.o index.o \

11
index.c
View File

@ -15,6 +15,11 @@
#include "kvec.h"
#include "khash.h"
#ifdef ANALYSIS_PERF
extern int64_t get_mseconds();
extern int64_t time_mm_idx_reader_read;
#endif
#define idx_hash(a) ((a)>>1)
#define idx_eq(a, b) ((a)>>1 == (b)>>1)
KHASH_INIT(idx, uint64_t, uint64_t, 1, idx_hash, idx_eq)
@ -608,6 +613,9 @@ void mm_idx_reader_close(mm_idx_reader_t *r)
mm_idx_t *mm_idx_reader_read(mm_idx_reader_t *r, int n_threads)
{
mm_idx_t *mi;
#ifdef ANALYSIS_PERF
int64_t tmp_cur_time = get_mseconds();
#endif
if (r->is_idx) {
mi = mm_idx_load(r->fp.idx);
if (mi && mm_verbose >= 2 && (mi->k != r->opt.k || mi->w != r->opt.w || (mi->flag&MM_I_HPC) != (r->opt.flag&MM_I_HPC)))
@ -618,6 +626,9 @@ mm_idx_t *mm_idx_reader_read(mm_idx_reader_t *r, int n_threads)
if (r->fp_out) mm_idx_dump(r->fp_out, mi);
mi->index = r->n_parts++;
}
#ifdef ANALYSIS_PERF
time_mm_idx_reader_read += get_mseconds() - tmp_cur_time;
#endif
return mi;
}

22
lchain.c
View File

@ -5,7 +5,11 @@
#include "mmpriv.h"
#include "kalloc.h"
#include "krmq.h"
#ifdef ANALYSIS_PERF
extern int64_t get_mseconds();
extern int64_t time_mg_lchain_dp,
time_mg_chain_backtrack;
#endif
static int64_t mg_chain_bk_end(int32_t max_drop, const mm128_t *z, const int32_t *f, const int64_t *p, int32_t *t, int64_t k)
{
int64_t i = z[k].y, end_i = -1, max_i = i;
@ -151,7 +155,9 @@ mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int
int32_t *f, *t, *v, n_u, n_v, mmax_f = 0, max_drop = bw;
int64_t *p, i, j, max_ii, st = 0, n_iter = 0;
uint64_t *u;
#ifdef ANALYSIS_PERF
int64_t tmp_cur_time = get_mseconds(), tmp_diff = 0;
#endif
if (_u) *_u = 0, *n_u_ = 0;
if (n == 0 || a == 0) {
kfree(km, a);
@ -205,10 +211,20 @@ mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int
max_ii = i;
if (mmax_f < max_f) mmax_f = max_f;
}
#ifdef ANALYSIS_PERF
int64_t tmp_inner_time = get_mseconds();
#endif
u = mg_chain_backtrack(km, n, f, p, v, t, min_cnt, min_sc, max_drop, &n_u, &n_v);
#ifdef ANALYSIS_PERF
tmp_diff = get_mseconds() - tmp_inner_time;
__sync_fetch_and_add(&time_mg_chain_backtrack, tmp_diff);
#endif
*n_u_ = n_u, *_u = u; // NB: note that u[] may not be sorted by score here
kfree(km, p); kfree(km, f); kfree(km, t);
#ifdef ANALYSIS_PERF
tmp_diff = get_mseconds() - tmp_cur_time;
__sync_fetch_and_add(&time_mg_lchain_dp, tmp_diff);
#endif
if (n_u == 0) {
kfree(km, a); kfree(km, v);
return 0;

587
main.c
View File

@ -7,6 +7,38 @@
#include "mmpriv.h"
#include "ketopt.h"
// 用来调试，计算感兴趣部分的运行时间
#include "sys/time.h"
// 获取当前毫秒数
int64_t get_mseconds()
{
struct timeval tv;
gettimeofday(&tv, NULL);
return (int64_t)1000 * (tv.tv_sec + ((1e-6) * tv.tv_usec));
}
// 记录运行时间的变量
#ifdef ANALYSIS_PERF
int64_t time_mm_idx_reader_read,
time_mm_map_file_frag,
time_map_work_for_block_1,
time_map_work_for_block_2,
time_map_work_for_block_3,
time_mm_map_frag_b1,
time_mm_map_frag_b2,
time_mm_map_frag_b3,
time_mm_map_frag_b4,
time_mm_map_frag_b5,
time_mm_map_frag_b6,
time_mg_lchain_dp = 0,
time_collect_seed_hits_heap = 0,
time_collect_seed_hits = 0,
time_mg_chain_backtrack = 0;
#endif
//////////////////////////////////
#ifdef __linux__
#include <sys/resource.h>
#include <sys/time.h>
@ -18,7 +50,9 @@ void liftrlimit()
setrlimit(RLIMIT_AS, &r);
}
#else
void liftrlimit() {}
void liftrlimit()
{
}
#endif
static ko_longopt_t long_options[] = {
@ -85,18 +119,21 @@ static ko_longopt_t long_options[] = {
{"mask-level", ko_required_argument, 'M'},
{"min-dp-score", ko_required_argument, 's'},
{"sam", ko_no_argument, 'a'},
{ 0, 0, 0 }
};
{0, 0, 0}};
static inline int64_t mm_parse_num2(const char *str, char **q)
{
double x;
char *p;
x = strtod(str, &p);
if (*p == 'G' || *p == 'g') x *= 1e9, ++p;
else if (*p == 'M' || *p == 'm') x *= 1e6, ++p;
else if (*p == 'K' || *p == 'k') x *= 1e3, ++p;
if (q) *q = p;
if (*p == 'G' || *p == 'g')
x *= 1e9, ++p;
else if (*p == 'M' || *p == 'm')
x *= 1e6, ++p;
else if (*p == 'K' || *p == 'k')
x *= 1e3, ++p;
if (q)
*q = p;
return (int64_t)(x + .499);
}
@ -107,14 +144,23 @@ static inline int64_t mm_parse_num(const char *str)
static inline void yes_or_no(mm_mapopt_t *opt, int64_t flag, int long_idx, const char *arg, int yes_to_set)
{
if (yes_to_set) {
if (strcmp(arg, "yes") == 0 || strcmp(arg, "y") == 0) opt->flag |= flag;
else if (strcmp(arg, "no") == 0 || strcmp(arg, "n") == 0) opt->flag &= ~flag;
else fprintf(stderr, "[WARNING]\033[1;31m option '--%s' only accepts 'yes' or 'no'.\033[0m\n", long_options[long_idx].name);
} else {
if (strcmp(arg, "yes") == 0 || strcmp(arg, "y") == 0) opt->flag &= ~flag;
else if (strcmp(arg, "no") == 0 || strcmp(arg, "n") == 0) opt->flag |= flag;
else fprintf(stderr, "[WARNING]\033[1;31m option '--%s' only accepts 'yes' or 'no'.\033[0m\n", long_options[long_idx].name);
if (yes_to_set)
{
if (strcmp(arg, "yes") == 0 || strcmp(arg, "y") == 0)
opt->flag |= flag;
else if (strcmp(arg, "no") == 0 || strcmp(arg, "n") == 0)
opt->flag &= ~flag;
else
fprintf(stderr, "[WARNING]\033[1;31m option '--%s' only accepts 'yes' or 'no'.\033[0m\n", long_options[long_idx].name);
}
else
{
if (strcmp(arg, "yes") == 0 || strcmp(arg, "y") == 0)
opt->flag &= ~flag;
else if (strcmp(arg, "no") == 0 || strcmp(arg, "n") == 0)
opt->flag |= flag;
else
fprintf(stderr, "[WARNING]\033[1;31m option '--%s' only accepts 'yes' or 'no'.\033[0m\n", long_options[long_idx].name);
}
}
@ -135,181 +181,345 @@ int main(int argc, char *argv[])
mm_realtime0 = realtime();
mm_set_opt(0, &ipt, &opt);
while ((c = ketopt(&o, argc, argv, 1, opt_str, long_options)) >= 0) { // test command line options and apply option -x/preset first
if (c == 'x') {
if (mm_set_opt(o.arg, &ipt, &opt) < 0) {
#ifdef ANALYSIS_PERF
time_mm_idx_reader_read = 0;
time_mm_map_file_frag = 0;
time_map_work_for_block_1 = 0;
time_map_work_for_block_2 = 0;
time_map_work_for_block_3 = 0;
time_mm_map_frag_b1 = 0;
time_mm_map_frag_b2 = 0;
time_mm_map_frag_b3 = 0;
time_mm_map_frag_b4 = 0;
time_mm_map_frag_b5 = 0;
time_mm_map_frag_b6 = 0;
#endif
while ((c = ketopt(&o, argc, argv, 1, opt_str, long_options)) >= 0)
{ // test command line options and apply option -x/preset first
if (c == 'x')
{
if (mm_set_opt(o.arg, &ipt, &opt) < 0)
{
fprintf(stderr, "[ERROR] unknown preset '%s'\n", o.arg);
return 1;
}
} else if (c == ':') {
}
else if (c == ':')
{
fprintf(stderr, "[ERROR] missing option argument\n");
return 1;
} else if (c == '?') {
}
else if (c == '?')
{
fprintf(stderr, "[ERROR] unknown option in \"%s\"\n", argv[o.i - 1]);
return 1;
}
}
o = KETOPT_INIT;
while ((c = ketopt(&o, argc, argv, 1, opt_str, long_options)) >= 0) {
if (c == 'w') ipt.w = atoi(o.arg);
else if (c == 'k') ipt.k = atoi(o.arg);
else if (c == 'H') ipt.flag |= MM_I_HPC;
else if (c == 'd') fnw = o.arg; // the above are indexing related options, except -I
else if (c == 't') n_threads = atoi(o.arg);
else if (c == 'v') mm_verbose = atoi(o.arg);
else if (c == 'g') opt.max_gap = (int)mm_parse_num(o.arg);
else if (c == 'G') mm_mapopt_max_intron_len(&opt, (int)mm_parse_num(o.arg));
else if (c == 'F') opt.max_frag_len = (int)mm_parse_num(o.arg);
else if (c == 'N') old_best_n = opt.best_n, opt.best_n = atoi(o.arg);
else if (c == 'p') opt.pri_ratio = atof(o.arg);
else if (c == 'M') opt.mask_level = atof(o.arg);
else if (c == 'c') opt.flag |= MM_F_OUT_CG | MM_F_CIGAR;
else if (c == 'D') opt.flag |= MM_F_NO_DIAG;
else if (c == 'P') opt.flag |= MM_F_ALL_CHAINS;
else if (c == 'X') opt.flag |= MM_F_ALL_CHAINS | MM_F_NO_DIAG | MM_F_NO_DUAL | MM_F_NO_LJOIN; // -D -P --no-long-join --dual=no
else if (c == 'a') opt.flag |= MM_F_OUT_SAM | MM_F_CIGAR;
else if (c == 'Q') opt.flag |= MM_F_NO_QUAL;
else if (c == 'Y') opt.flag |= MM_F_SOFTCLIP;
else if (c == 'L') opt.flag |= MM_F_LONG_CIGAR;
else if (c == 'y') opt.flag |= MM_F_COPY_COMMENT;
else if (c == 'T') opt.sdust_thres = atoi(o.arg);
else if (c == 'n') opt.min_cnt = atoi(o.arg);
else if (c == 'm') opt.min_chain_score = atoi(o.arg);
else if (c == 'A') opt.a = atoi(o.arg);
else if (c == 'B') opt.b = atoi(o.arg);
else if (c == 'b') opt.transition = atoi(o.arg);
else if (c == 's') opt.min_dp_max = atoi(o.arg);
else if (c == 'C') opt.noncan = atoi(o.arg);
else if (c == 'I') ipt.batch_size = mm_parse_num(o.arg);
else if (c == 'K') opt.mini_batch_size = mm_parse_num(o.arg);
else if (c == 'e') opt.occ_dist = mm_parse_num(o.arg);
else if (c == 'R') rg = o.arg;
else if (c == 'h') fp_help = stdout;
else if (c == '2') opt.flag |= MM_F_2_IO_THREADS;
else if (c == 'J') {
while ((c = ketopt(&o, argc, argv, 1, opt_str, long_options)) >= 0)
{
if (c == 'w')
ipt.w = atoi(o.arg);
else if (c == 'k')
ipt.k = atoi(o.arg);
else if (c == 'H')
ipt.flag |= MM_I_HPC;
else if (c == 'd')
fnw = o.arg; // the above are indexing related options, except -I
else if (c == 't')
n_threads = atoi(o.arg);
else if (c == 'v')
mm_verbose = atoi(o.arg);
else if (c == 'g')
opt.max_gap = (int)mm_parse_num(o.arg);
else if (c == 'G')
mm_mapopt_max_intron_len(&opt, (int)mm_parse_num(o.arg));
else if (c == 'F')
opt.max_frag_len = (int)mm_parse_num(o.arg);
else if (c == 'N')
old_best_n = opt.best_n, opt.best_n = atoi(o.arg);
else if (c == 'p')
opt.pri_ratio = atof(o.arg);
else if (c == 'M')
opt.mask_level = atof(o.arg);
else if (c == 'c')
opt.flag |= MM_F_OUT_CG | MM_F_CIGAR;
else if (c == 'D')
opt.flag |= MM_F_NO_DIAG;
else if (c == 'P')
opt.flag |= MM_F_ALL_CHAINS;
else if (c == 'X')
opt.flag |= MM_F_ALL_CHAINS | MM_F_NO_DIAG | MM_F_NO_DUAL | MM_F_NO_LJOIN; // -D -P --no-long-join --dual=no
else if (c == 'a')
opt.flag |= MM_F_OUT_SAM | MM_F_CIGAR;
else if (c == 'Q')
opt.flag |= MM_F_NO_QUAL;
else if (c == 'Y')
opt.flag |= MM_F_SOFTCLIP;
else if (c == 'L')
opt.flag |= MM_F_LONG_CIGAR;
else if (c == 'y')
opt.flag |= MM_F_COPY_COMMENT;
else if (c == 'T')
opt.sdust_thres = atoi(o.arg);
else if (c == 'n')
opt.min_cnt = atoi(o.arg);
else if (c == 'm')
opt.min_chain_score = atoi(o.arg);
else if (c == 'A')
opt.a = atoi(o.arg);
else if (c == 'B')
opt.b = atoi(o.arg);
else if (c == 'b')
opt.transition = atoi(o.arg);
else if (c == 's')
opt.min_dp_max = atoi(o.arg);
else if (c == 'C')
opt.noncan = atoi(o.arg);
else if (c == 'I')
ipt.batch_size = mm_parse_num(o.arg);
else if (c == 'K')
opt.mini_batch_size = mm_parse_num(o.arg);
else if (c == 'e')
opt.occ_dist = mm_parse_num(o.arg);
else if (c == 'R')
rg = o.arg;
else if (c == 'h')
fp_help = stdout;
else if (c == '2')
opt.flag |= MM_F_2_IO_THREADS;
else if (c == 'J')
{
int t;
t = atoi(o.arg);
if (t == 0) opt.flag |= MM_F_SPLICE_OLD;
else if (t == 1) opt.flag &= ~MM_F_SPLICE_OLD;
} else if (c == 'o') {
if (strcmp(o.arg, "-") != 0) {
if (freopen(o.arg, "wb", stdout) == NULL) {
if (t == 0)
opt.flag |= MM_F_SPLICE_OLD;
else if (t == 1)
opt.flag &= ~MM_F_SPLICE_OLD;
}
else if (c == 'o')
{
if (strcmp(o.arg, "-") != 0)
{
if (freopen(o.arg, "wb", stdout) == NULL)
{
fprintf(stderr, "[ERROR]\033[1;31m failed to write the output to file '%s'\033[0m: %s\n", o.arg, strerror(errno));
exit(1);
}
}
}
else if (c == 300) ipt.bucket_bits = atoi(o.arg); // --bucket-bits
else if (c == 302) opt.seed = atoi(o.arg); // --seed
else if (c == 303) mm_dbg_flag |= MM_DBG_NO_KALLOC; // --no-kalloc
else if (c == 304) mm_dbg_flag |= MM_DBG_PRINT_QNAME; // --print-qname
else if (c == 306) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_SEED, n_threads = 1; // --print-seed
else if (c == 307) opt.max_chain_skip = atoi(o.arg); // --max-chain-skip
else if (c == 339) opt.max_chain_iter = atoi(o.arg); // --max-chain-iter
else if (c == 308) opt.min_ksw_len = atoi(o.arg); // --min-dp-len
else if (c == 309) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_ALN_SEQ, n_threads = 1; // --print-aln-seq
else if (c == 310) opt.flag |= MM_F_SPLICE; // --splice
else if (c == 312) opt.flag |= MM_F_NO_LJOIN; // --no-long-join
else if (c == 313) opt.flag |= MM_F_SR; // --sr
else if (c == 317) opt.end_bonus = atoi(o.arg); // --end-bonus
else if (c == 318) opt.flag |= MM_F_INDEPEND_SEG; // --no-pairing
else if (c == 320) ipt.flag |= MM_I_NO_SEQ; // --idx-no-seq
else if (c == 321) opt.anchor_ext_shift = atoi(o.arg); // --end-seed-pen
else if (c == 322) opt.flag |= MM_F_FOR_ONLY; // --for-only
else if (c == 323) opt.flag |= MM_F_REV_ONLY; // --rev-only
else if (c == 327) opt.max_clip_ratio = atof(o.arg); // --max-clip-ratio
else if (c == 328) opt.min_mid_occ = atoi(o.arg); // --min-occ-floor
else if (c == 329) opt.flag |= MM_F_OUT_MD; // --MD
else if (c == 331) opt.sc_ambi = atoi(o.arg); // --score-N
else if (c == 332) opt.flag |= MM_F_EQX; // --eqx
else if (c == 333) opt.flag |= MM_F_PAF_NO_HIT; // --paf-no-hit
else if (c == 334) opt.split_prefix = o.arg; // --split-prefix
else if (c == 335) opt.flag |= MM_F_NO_END_FLT; // --no-end-flt
else if (c == 336) opt.flag |= MM_F_HARD_MLEVEL; // --hard-mask-level
else if (c == 337) opt.max_sw_mat = mm_parse_num(o.arg); // --cap-sw-mat
else if (c == 338) opt.max_qlen = mm_parse_num(o.arg); // --max-qlen
else if (c == 340) junc_bed = o.arg; // --junc-bed
else if (c == 341) opt.junc_bonus = atoi(o.arg); // --junc-bonus
else if (c == 342) opt.flag |= MM_F_SAM_HIT_ONLY; // --sam-hit-only
else if (c == 343) opt.chain_gap_scale = atof(o.arg); // --chain-gap-scale
else if (c == 351) opt.chain_skip_scale = atof(o.arg); // --chain-skip-scale
else if (c == 344) alt_list = o.arg; // --alt
else if (c == 345) opt.alt_drop = atof(o.arg); // --alt-drop
else if (c == 346) opt.mask_len = mm_parse_num(o.arg); // --mask-len
else if (c == 348) opt.flag |= MM_F_QSTRAND | MM_F_NO_INV; // --qstrand
else if (c == 349) opt.cap_kalloc = mm_parse_num(o.arg); // --cap-kalloc
else if (c == 350) opt.q_occ_frac = atof(o.arg); // --q-occ-frac
else if (c == 352) mm_dbg_flag |= MM_DBG_PRINT_CHAIN; // --print-chains
else if (c == 353) opt.flag |= MM_F_NO_HASH_NAME; // --no-hash-name
else if (c == 354) opt.flag |= MM_F_SECONDARY_SEQ; // --secondary-seq
else if (c == 330) {
else if (c == 300)
ipt.bucket_bits = atoi(o.arg); // --bucket-bits
else if (c == 302)
opt.seed = atoi(o.arg); // --seed
else if (c == 303)
mm_dbg_flag |= MM_DBG_NO_KALLOC; // --no-kalloc
else if (c == 304)
mm_dbg_flag |= MM_DBG_PRINT_QNAME; // --print-qname
else if (c == 306)
mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_SEED, n_threads = 1; // --print-seed
else if (c == 307)
opt.max_chain_skip = atoi(o.arg); // --max-chain-skip
else if (c == 339)
opt.max_chain_iter = atoi(o.arg); // --max-chain-iter
else if (c == 308)
opt.min_ksw_len = atoi(o.arg); // --min-dp-len
else if (c == 309)
mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_ALN_SEQ, n_threads = 1; // --print-aln-seq
else if (c == 310)
opt.flag |= MM_F_SPLICE; // --splice
else if (c == 312)
opt.flag |= MM_F_NO_LJOIN; // --no-long-join
else if (c == 313)
opt.flag |= MM_F_SR; // --sr
else if (c == 317)
opt.end_bonus = atoi(o.arg); // --end-bonus
else if (c == 318)
opt.flag |= MM_F_INDEPEND_SEG; // --no-pairing
else if (c == 320)
ipt.flag |= MM_I_NO_SEQ; // --idx-no-seq
else if (c == 321)
opt.anchor_ext_shift = atoi(o.arg); // --end-seed-pen
else if (c == 322)
opt.flag |= MM_F_FOR_ONLY; // --for-only
else if (c == 323)
opt.flag |= MM_F_REV_ONLY; // --rev-only
else if (c == 327)
opt.max_clip_ratio = atof(o.arg); // --max-clip-ratio
else if (c == 328)
opt.min_mid_occ = atoi(o.arg); // --min-occ-floor
else if (c == 329)
opt.flag |= MM_F_OUT_MD; // --MD
else if (c == 331)
opt.sc_ambi = atoi(o.arg); // --score-N
else if (c == 332)
opt.flag |= MM_F_EQX; // --eqx
else if (c == 333)
opt.flag |= MM_F_PAF_NO_HIT; // --paf-no-hit
else if (c == 334)
opt.split_prefix = o.arg; // --split-prefix
else if (c == 335)
opt.flag |= MM_F_NO_END_FLT; // --no-end-flt
else if (c == 336)
opt.flag |= MM_F_HARD_MLEVEL; // --hard-mask-level
else if (c == 337)
opt.max_sw_mat = mm_parse_num(o.arg); // --cap-sw-mat
else if (c == 338)
opt.max_qlen = mm_parse_num(o.arg); // --max-qlen
else if (c == 340)
junc_bed = o.arg; // --junc-bed
else if (c == 341)
opt.junc_bonus = atoi(o.arg); // --junc-bonus
else if (c == 342)
opt.flag |= MM_F_SAM_HIT_ONLY; // --sam-hit-only
else if (c == 343)
opt.chain_gap_scale = atof(o.arg); // --chain-gap-scale
else if (c == 351)
opt.chain_skip_scale = atof(o.arg); // --chain-skip-scale
else if (c == 344)
alt_list = o.arg; // --alt
else if (c == 345)
opt.alt_drop = atof(o.arg); // --alt-drop
else if (c == 346)
opt.mask_len = mm_parse_num(o.arg); // --mask-len
else if (c == 348)
opt.flag |= MM_F_QSTRAND | MM_F_NO_INV; // --qstrand
else if (c == 349)
opt.cap_kalloc = mm_parse_num(o.arg); // --cap-kalloc
else if (c == 350)
opt.q_occ_frac = atof(o.arg); // --q-occ-frac
else if (c == 352)
mm_dbg_flag |= MM_DBG_PRINT_CHAIN; // --print-chains
else if (c == 353)
opt.flag |= MM_F_NO_HASH_NAME; // --no-hash-name
else if (c == 354)
opt.flag |= MM_F_SECONDARY_SEQ; // --secondary-seq
else if (c == 330)
{
fprintf(stderr, "[WARNING] \033[1;31m --lj-min-ratio has been deprecated.\033[0m\n");
} else if (c == 314) { // --frag
}
else if (c == 314)
{ // --frag
yes_or_no(&opt, MM_F_FRAG_MODE, o.longidx, o.arg, 1);
} else if (c == 315) { // --secondary
}
else if (c == 315)
{ // --secondary
yes_or_no(&opt, MM_F_NO_PRINT_2ND, o.longidx, o.arg, 0);
} else if (c == 316) { // --cs
}
else if (c == 316)
{ // --cs
opt.flag |= MM_F_OUT_CS | MM_F_CIGAR;
if (o.arg == 0 || strcmp(o.arg, "short") == 0) {
if (o.arg == 0 || strcmp(o.arg, "short") == 0)
{
opt.flag &= ~MM_F_OUT_CS_LONG;
} else if (strcmp(o.arg, "long") == 0) {
}
else if (strcmp(o.arg, "long") == 0)
{
opt.flag |= MM_F_OUT_CS_LONG;
} else if (strcmp(o.arg, "none") == 0) {
}
else if (strcmp(o.arg, "none") == 0)
{
opt.flag &= ~MM_F_OUT_CS;
} else if (mm_verbose >= 2) {
}
else if (mm_verbose >= 2)
{
fprintf(stderr, "[WARNING]\033[1;31m --cs only takes 'short' or 'long'. Invalid values are assumed to be 'short'.\033[0m\n");
}
} else if (c == 319) { // --splice-flank
}
else if (c == 319)
{ // --splice-flank
yes_or_no(&opt, MM_F_SPLICE_FLANK, o.longidx, o.arg, 1);
} else if (c == 324) { // --heap-sort
}
else if (c == 324)
{ // --heap-sort
yes_or_no(&opt, MM_F_HEAP_SORT, o.longidx, o.arg, 1);
} else if (c == 326) { // --dual
}
else if (c == 326)
{ // --dual
yes_or_no(&opt, MM_F_NO_DUAL, o.longidx, o.arg, 0);
} else if (c == 347) { // --rmq
if (o.arg) yes_or_no(&opt, MM_F_RMQ, o.longidx, o.arg, 1);
else opt.flag |= MM_F_RMQ;
} else if (c == 'S') {
}
else if (c == 347)
{ // --rmq
if (o.arg)
yes_or_no(&opt, MM_F_RMQ, o.longidx, o.arg, 1);
else
opt.flag |= MM_F_RMQ;
}
else if (c == 'S')
{
opt.flag |= MM_F_OUT_CS | MM_F_CIGAR | MM_F_OUT_CS_LONG;
if (mm_verbose >= 2)
fprintf(stderr, "[WARNING]\033[1;31m option -S is deprecated and may be removed in future. Please use --cs=long instead.\033[0m\n");
} else if (c == 'V') {
}
else if (c == 'V')
{
puts(MM_VERSION);
return 0;
} else if (c == 'r') {
}
else if (c == 'r')
{
opt.bw = (int)mm_parse_num2(o.arg, &s);
if (*s == ',') opt.bw_long = (int)mm_parse_num2(s + 1, &s);
} else if (c == 'U') {
if (*s == ',')
opt.bw_long = (int)mm_parse_num2(s + 1, &s);
}
else if (c == 'U')
{
opt.min_mid_occ = strtol(o.arg, &s, 10);
if (*s == ',') opt.max_mid_occ = strtol(s + 1, &s, 10);
} else if (c == 'f') {
if (*s == ',')
opt.max_mid_occ = strtol(s + 1, &s, 10);
}
else if (c == 'f')
{
double x;
char *p;
x = strtod(o.arg, &p);
if (x < 1.0) opt.mid_occ_frac = x, opt.mid_occ = 0;
else opt.mid_occ = (int)(x + .499);
if (*p == ',') opt.max_occ = (int)(strtod(p+1, &p) + .499);
} else if (c == 'u') {
if (*o.arg == 'b') opt.flag |= MM_F_SPLICE_FOR|MM_F_SPLICE_REV; // both strands
else if (*o.arg == 'f') opt.flag |= MM_F_SPLICE_FOR, opt.flag &= ~MM_F_SPLICE_REV; // match GT-AG
else if (*o.arg == 'r') opt.flag |= MM_F_SPLICE_REV, opt.flag &= ~MM_F_SPLICE_FOR; // match CT-AC (reverse complement of GT-AG)
else if (*o.arg == 'n') opt.flag &= ~(MM_F_SPLICE_FOR|MM_F_SPLICE_REV); // don't try to match the GT-AG signal
else {
if (x < 1.0)
opt.mid_occ_frac = x, opt.mid_occ = 0;
else
opt.mid_occ = (int)(x + .499);
if (*p == ',')
opt.max_occ = (int)(strtod(p + 1, &p) + .499);
}
else if (c == 'u')
{
if (*o.arg == 'b')
opt.flag |= MM_F_SPLICE_FOR | MM_F_SPLICE_REV; // both strands
else if (*o.arg == 'f')
opt.flag |= MM_F_SPLICE_FOR, opt.flag &= ~MM_F_SPLICE_REV; // match GT-AG
else if (*o.arg == 'r')
opt.flag |= MM_F_SPLICE_REV, opt.flag &= ~MM_F_SPLICE_FOR; // match CT-AC (reverse complement of GT-AG)
else if (*o.arg == 'n')
opt.flag &= ~(MM_F_SPLICE_FOR | MM_F_SPLICE_REV); // don't try to match the GT-AG signal
else
{
fprintf(stderr, "[ERROR]\033[1;31m unrecognized cDNA direction\033[0m\n");
return 1;
}
} else if (c == 'z') {
}
else if (c == 'z')
{
opt.zdrop = opt.zdrop_inv = strtol(o.arg, &s, 10);
if (*s == ',') opt.zdrop_inv = strtol(s + 1, &s, 10);
} else if (c == 'O') {
if (*s == ',')
opt.zdrop_inv = strtol(s + 1, &s, 10);
}
else if (c == 'O')
{
opt.q = opt.q2 = strtol(o.arg, &s, 10);
if (*s == ',') opt.q2 = strtol(s + 1, &s, 10);
} else if (c == 'E') {
if (*s == ',')
opt.q2 = strtol(s + 1, &s, 10);
}
else if (c == 'E')
{
opt.e = opt.e2 = strtol(o.arg, &s, 10);
if (*s == ',') opt.e2 = strtol(s + 1, &s, 10);
if (*s == ',')
opt.e2 = strtol(s + 1, &s, 10);
}
}
if ((opt.flag & MM_F_SPLICE) && (opt.flag & MM_F_FRAG_MODE)) {
if ((opt.flag & MM_F_SPLICE) && (opt.flag & MM_F_FRAG_MODE))
{
fprintf(stderr, "[ERROR]\033[1;31m --splice and --frag should not be specified at the same time.\033[0m\n");
return 1;
}
@ -317,12 +527,14 @@ int main(int argc, char *argv[])
ipt.flag |= MM_I_NO_SEQ;
if (mm_check_opt(&ipt, &opt) < 0)
return 1;
if (opt.best_n == 0) {
if (opt.best_n == 0)
{
fprintf(stderr, "[WARNING]\033[1;31m changed '-N 0' to '-N %d --secondary=no'.\033[0m\n", old_best_n);
opt.best_n = old_best_n, opt.flag |= MM_F_NO_PRINT_2ND;
}
if (argc == o.ind || fp_help == stdout) {
if (argc == o.ind || fp_help == stdout)
{
fprintf(fp_help, "Usage: minimap2 [options] <target.fa>|<target.idx> [query.fa] [...]\n");
fprintf(fp_help, "Options:\n");
fprintf(fp_help, " Indexing:\n");
@ -378,42 +590,52 @@ int main(int argc, char *argv[])
return fp_help == stdout ? 0 : 1;
}
if ((opt.flag & MM_F_SR) && argc - o.ind > 3) {
if ((opt.flag & MM_F_SR) && argc - o.ind > 3)
{
fprintf(stderr, "[ERROR] incorrect input: in the sr mode, please specify no more than two query files.\n");
return 1;
}
idx_rdr = mm_idx_reader_open(argv[o.ind], &ipt, fnw);
if (idx_rdr == 0) {
if (idx_rdr == 0)
{
fprintf(stderr, "[ERROR] failed to open file '%s': %s\n", argv[o.ind], strerror(errno));
return 1;
}
if (!idx_rdr->is_idx && fnw == 0 && argc - o.ind < 2) {
if (!idx_rdr->is_idx && fnw == 0 && argc - o.ind < 2)
{
fprintf(stderr, "[ERROR] missing input: please specify a query file to map or option -d to keep the index\n");
mm_idx_reader_close(idx_rdr);
return 1;
}
if (opt.best_n == 0 && (opt.flag & MM_F_CIGAR) && mm_verbose >= 2)
fprintf(stderr, "[WARNING]\033[1;31m `-N 0' reduces alignment accuracy. Please use --secondary=no to suppress secondary alignments.\033[0m\n");
while ((mi = mm_idx_reader_read(idx_rdr, n_threads)) != 0) {
while ((mi = mm_idx_reader_read(idx_rdr, n_threads)) != 0)
{
int ret;
if ((opt.flag & MM_F_CIGAR) && (mi->flag & MM_I_NO_SEQ)) {
if ((opt.flag & MM_F_CIGAR) && (mi->flag & MM_I_NO_SEQ))
{
fprintf(stderr, "[ERROR] the prebuilt index doesn't contain sequences.\n");
mm_idx_destroy(mi);
mm_idx_reader_close(idx_rdr);
return 1;
}
if ((opt.flag & MM_F_OUT_SAM) && idx_rdr->n_parts == 1) {
if (mm_idx_reader_eof(idx_rdr)) {
if ((opt.flag & MM_F_OUT_SAM) && idx_rdr->n_parts == 1)
{
if (mm_idx_reader_eof(idx_rdr))
{
if (opt.split_prefix == 0)
ret = mm_write_sam_hdr(mi, rg, MM_VERSION, argc, argv);
else
ret = mm_write_sam_hdr(0, rg, MM_VERSION, argc, argv);
} else {
}
else
{
ret = mm_write_sam_hdr(0, rg, MM_VERSION, argc, argv);
if (opt.split_prefix == 0 && mm_verbose >= 2)
fprintf(stderr, "[WARNING]\033[1;31m For a multi-part index, no @SQ lines will be outputted. Please use --split-prefix.\033[0m\n");
}
if (ret != 0) {
if (ret != 0)
{
mm_idx_destroy(mi);
mm_idx_reader_close(idx_rdr);
return 1;
@ -422,25 +644,37 @@ int main(int argc, char *argv[])
if (mm_verbose >= 3)
fprintf(stderr, "[M::%s::%.3f*%.2f] loaded/built the index for %d target sequence(s)\n",
__func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), mi->n_seq);
if (argc != o.ind + 1) mm_mapopt_update(&opt, mi);
if (mm_verbose >= 3) mm_idx_stat(mi);
if (junc_bed) mm_idx_bed_read(mi, junc_bed, 1);
if (alt_list) mm_idx_alt_read(mi, alt_list);
if (argc - (o.ind + 1) == 0) {
if (argc != o.ind + 1)
mm_mapopt_update(&opt, mi);
if (mm_verbose >= 3)
mm_idx_stat(mi);
if (junc_bed)
mm_idx_bed_read(mi, junc_bed, 1);
if (alt_list)
mm_idx_alt_read(mi, alt_list);
if (argc - (o.ind + 1) == 0)
{
mm_idx_destroy(mi);
continue; // no query files
}
ret = 0;
if (!(opt.flag & MM_F_FRAG_MODE)) {
for (i = o.ind + 1; i < argc; ++i) {
if (!(opt.flag & MM_F_FRAG_MODE))
{
// 如果是read / read overlap操作，那么把第一个输入的read.fq当做索引，后续的read比对到该索引
for (i = o.ind + 1; i < argc; ++i)
{
ret = mm_map_file(mi, argv[i], &opt, n_threads);
if (ret < 0) break;
if (ret < 0)
break;
}
} else {
}
else
{
ret = mm_map_file_frag(mi, argc - (o.ind + 1), (const char **)&argv[o.ind + 1], &opt, n_threads);
}
mm_idx_destroy(mi);
if (ret < 0) {
if (ret < 0)
{
fprintf(stderr, "ERROR: failed to map the query file\n");
exit(EXIT_FAILURE);
}
@ -451,17 +685,40 @@ int main(int argc, char *argv[])
if (opt.split_prefix)
mm_split_merge(argc - (o.ind + 1), (const char **)&argv[o.ind + 1], &opt, n_parts);
if (fflush(stdout) == EOF) {
if (fflush(stdout) == EOF)
{
perror("[ERROR] failed to write the results");
exit(EXIT_FAILURE);
}
if (mm_verbose >= 3) {
if (mm_verbose >= 3)
{
fprintf(stderr, "[M::%s] Version: %s\n", __func__, MM_VERSION);
fprintf(stderr, "[M::%s] CMD:", __func__);
for (i = 0; i < argc; ++i)
fprintf(stderr, " %s", argv[i]);
fprintf(stderr, "\n[M::%s] Real time: %.3f sec; CPU: %.3f sec; Peak RSS: %.3f GB\n", __func__, realtime() - mm_realtime0, cputime(), peakrss() / 1024.0 / 1024.0 / 1024.0);
}
#ifdef ANALYSIS_PERF
fprintf(stderr, "\n");
fprintf(stderr, "time_mm_idx_reader_read: %f s\n", time_mm_idx_reader_read / 1000.0);
fprintf(stderr, "time_mm_map_file_frag: %f s\n", time_mm_map_file_frag / 1000.0);
fprintf(stderr, "time_map_work_for_block_1: %f s\n", time_map_work_for_block_1 / 1000.0 / n_threads);
fprintf(stderr, "time_map_work_for_block_2: %f s\n", time_map_work_for_block_2 / 1000.0 / n_threads);
fprintf(stderr, "time_map_work_for_block_3: %f s\n", time_map_work_for_block_3 / 1000.0 / n_threads);
fprintf(stderr, "time_mm_map_frag_b1: %f s\n", time_mm_map_frag_b1 / 1000.0 / n_threads);
fprintf(stderr, "time_mm_map_frag_b2: %f s\n", time_mm_map_frag_b2 / 1000.0 / n_threads);
fprintf(stderr, "time_mm_map_frag_b3: %f s\n", time_mm_map_frag_b3 / 1000.0 / n_threads);
fprintf(stderr, "time_mm_map_frag_b4: %f s\n", time_mm_map_frag_b4 / 1000.0 / n_threads);
fprintf(stderr, "time_mm_map_frag_b5: %f s\n", time_mm_map_frag_b5 / 1000.0 / n_threads);
fprintf(stderr, "time_mm_map_frag_b6: %f s\n", time_mm_map_frag_b6 / 1000.0 / n_threads);
fprintf(stderr, "time_collect_seed_hits_heap: %f s\n", time_collect_seed_hits_heap / 1000.0 / n_threads);
fprintf(stderr, "time_collect_seed_hits: %f s\n", time_collect_seed_hits / 1000.0 / n_threads);
fprintf(stderr, "time_mg_lchain_dp: %f s\n", time_mg_lchain_dp / 1000.0 / n_threads);
fprintf(stderr, "time_mg_chain_backtrack: %f s\n", time_mg_chain_backtrack / 1000.0 / n_threads);
fprintf(stderr, "\n");
#endif
return 0;
}

501
map.c
View File

@ -10,17 +10,37 @@
#include "bseq.h"
#include "khash.h"
#ifdef ANALYSIS_PERF
extern int64_t get_mseconds();
extern int64_t time_mm_map_file_frag,
time_map_work_for_block_1,
time_map_work_for_block_2,
time_map_work_for_block_3,
time_mm_map_frag_b1,
time_mm_map_frag_b2,
time_mm_map_frag_b3,
time_mm_map_frag_b4,
time_mm_map_frag_b5,
time_mm_map_frag_b6,
time_mg_lchain_dp,
time_collect_seed_hits_heap,
time_collect_seed_hits,
time_mg_chain_backtrack;
#endif
mm_tbuf_t *mm_tbuf_init(void)
{
mm_tbuf_t *b;
b = (mm_tbuf_t *)calloc(1, sizeof(mm_tbuf_t));
if (!(mm_dbg_flag & 1)) b->km = km_init();
if (!(mm_dbg_flag & 1))
b->km = km_init();
return b;
}
void mm_tbuf_destroy(mm_tbuf_t *b)
{
if (b == 0) return;
if (b == 0)
return;
km_destroy(b->km);
free(b);
}
@ -35,22 +55,30 @@ static int mm_dust_minier(void *km, int n, mm128_t *a, int l_seq, const char *se
int n_dreg, j, k, u = 0;
const uint64_t *dreg;
sdust_buf_t *sdb;
if (sdust_thres <= 0) return n;
if (sdust_thres <= 0)
return n;
sdb = sdust_buf_init(km);
dreg = sdust_core((const uint8_t *)seq, l_seq, sdust_thres, 64, &n_dreg, sdb);
for (j = k = 0; j < n; ++j) { // squeeze out minimizers that significantly overlap with LCRs
for (j = k = 0; j < n; ++j)
{ // squeeze out minimizers that significantly overlap with LCRs
int32_t qpos = (uint32_t)a[j].y >> 1, span = a[j].x & 0xff;
int32_t s = qpos - (span - 1), e = s + span;
while (u < n_dreg && (int32_t)dreg[u] <= s) ++u;
if (u < n_dreg && (int32_t)(dreg[u]>>32) < e) {
while (u < n_dreg && (int32_t)dreg[u] <= s)
++u;
if (u < n_dreg && (int32_t)(dreg[u] >> 32) < e)
{
int v, l = 0;
for (v = u; v < n_dreg && (int32_t)(dreg[v]>>32) < e; ++v) { // iterate over LCRs overlapping this minimizer
for (v = u; v < n_dreg && (int32_t)(dreg[v] >> 32) < e; ++v)
{ // iterate over LCRs overlapping this minimizer
int ss = s > (int32_t)(dreg[v] >> 32) ? s : dreg[v] >> 32;
int ee = e < (int32_t)dreg[v] ? e : (uint32_t)dreg[v];
l += ee - ss;
}
if (l <= span>>1) a[k++] = a[j]; // keep the minimizer if less than half of it falls in masked region
} else a[k++] = a[j];
if (l <= span >> 1)
a[k++] = a[j]; // keep the minimizer if less than half of it falls in masked region
}
else
a[k++] = a[j];
}
sdust_buf_destroy(sdb);
return k; // the new size
@ -60,7 +88,8 @@ static void collect_minimizers(void *km, const mm_mapopt_t *opt, const mm_idx_t
{
int i, n, sum = 0;
mv->n = 0;
for (i = n = 0; i < n_segs; ++i) {
for (i = n = 0; i < n_segs; ++i)
{
size_t j;
mm_sketch(km, seqs[i], qlens[i], mi->w, mi->k, i, mi->flag & MM_I_HPC, mv);
for (j = n; j < mv->n; ++j)
@ -78,22 +107,32 @@ KSORT_INIT(heap, mm128_t, heap_lt)
static inline int skip_seed(int flag, uint64_t r, const mm_seed_t *q, const char *qname, int qlen, const mm_idx_t *mi, int *is_self)
{
*is_self = 0;
if (qname && (flag & (MM_F_NO_DIAG|MM_F_NO_DUAL))) {
if (qname && (flag & (MM_F_NO_DIAG | MM_F_NO_DUAL)))
{
const mm_idx_seq_t *s = &mi->seq[r >> 32];
int cmp;
cmp = strcmp(qname, s->name);
if ((flag&MM_F_NO_DIAG) && cmp == 0 && (int)s->len == qlen) {
if ((uint32_t)r>>1 == (q->q_pos>>1)) return 1; // avoid the diagnonal anchors
if ((r&1) == (q->q_pos&1)) *is_self = 1; // this flag is used to avoid spurious extension on self chain
if ((flag & MM_F_NO_DIAG) && cmp == 0 && (int)s->len == qlen)
{
if ((uint32_t)r >> 1 == (q->q_pos >> 1))
return 1; // avoid the diagnonal anchors
if ((r & 1) == (q->q_pos & 1))
*is_self = 1; // this flag is used to avoid spurious extension on self chain
}
if ((flag & MM_F_NO_DUAL) && cmp > 0) // all-vs-all mode: map once
return 1;
}
if (flag & (MM_F_FOR_ONLY|MM_F_REV_ONLY)) {
if ((r&1) == (q->q_pos&1)) { // forward strand
if (flag & MM_F_REV_ONLY) return 1;
} else {
if (flag & MM_F_FOR_ONLY) return 1;
if (flag & (MM_F_FOR_ONLY | MM_F_REV_ONLY))
{
if ((r & 1) == (q->q_pos & 1))
{ // forward strand
if (flag & MM_F_REV_ONLY)
return 1;
}
else
{
if (flag & MM_F_FOR_ONLY)
return 1;
}
}
return 0;
@ -106,44 +145,58 @@ static mm128_t *collect_seed_hits_heap(void *km, const mm_mapopt_t *opt, int max
int64_t j, n_for = 0, n_rev = 0;
mm_seed_t *m;
mm128_t *a, *heap;
#ifdef ANALYSIS_PERF
int64_t tmp_cur_time = get_mseconds(), tmp_diff = 0;
#endif
m = mm_collect_matches(km, &n_m, qlen, max_occ, opt->max_max_occ, opt->occ_dist, mi, mv, n_a, rep_len, n_mini_pos, mini_pos);
heap = (mm128_t *)kmalloc(km, n_m * sizeof(mm128_t));
a = (mm128_t *)kmalloc(km, *n_a * sizeof(mm128_t));
for (i = 0, heap_size = 0; i < n_m; ++i) {
if (m[i].n > 0) {
for (i = 0, heap_size = 0; i < n_m; ++i)
{
if (m[i].n > 0)
{
heap[heap_size].x = m[i].cr[0];
heap[heap_size].y = (uint64_t)i << 32;
++heap_size;
}
}
ks_heapmake_heap(heap_size, heap);
while (heap_size > 0) {
while (heap_size > 0)
{
mm_seed_t *q = &m[heap->y >> 32];
mm128_t *p;
uint64_t r = heap->x;
int32_t is_self, rpos = (uint32_t)r >> 1;
if (!skip_seed(opt->flag, r, q, qname, qlen, mi, &is_self)) {
if ((r&1) == (q->q_pos&1)) { // forward strand
if (!skip_seed(opt->flag, r, q, qname, qlen, mi, &is_self))
{
if ((r & 1) == (q->q_pos & 1))
{ // forward strand
p = &a[n_for++];
p->x = (r & 0xffffffff00000000ULL) | rpos;
p->y = (uint64_t)q->q_span << 32 | q->q_pos >> 1;
} else { // reverse strand
}
else
{ // reverse strand
p = &a[(*n_a) - (++n_rev)];
p->x = 1ULL << 63 | (r & 0xffffffff00000000ULL) | rpos;
p->y = (uint64_t)q->q_span << 32 | (qlen - ((q->q_pos >> 1) + 1 - q->q_span) - 1);
}
p->y |= (uint64_t)q->seg_id << MM_SEED_SEG_SHIFT;
if (q->is_tandem) p->y |= MM_SEED_TANDEM;
if (is_self) p->y |= MM_SEED_SELF;
if (q->is_tandem)
p->y |= MM_SEED_TANDEM;
if (is_self)
p->y |= MM_SEED_SELF;
}
// update the heap
if ((uint32_t)heap->y < q->n - 1) {
if ((uint32_t)heap->y < q->n - 1)
{
++heap[0].y;
heap[0].x = m[heap[0].y >> 32].cr[(uint32_t)heap[0].y];
} else {
}
else
{
heap[0] = heap[heap_size - 1];
--heap_size;
}
@ -153,15 +206,21 @@ static mm128_t *collect_seed_hits_heap(void *km, const mm_mapopt_t *opt, int max
kfree(km, heap);
// reverse anchors on the reverse strand, as they are in the descending order
for (j = 0; j < n_rev>>1; ++j) {
for (j = 0; j < n_rev >> 1; ++j)
{
mm128_t t = a[(*n_a) - 1 - j];
a[(*n_a) - 1 - j] = a[(*n_a) - (n_rev - j)];
a[(*n_a) - (n_rev - j)] = t;
}
if (*n_a > n_for + n_rev) {
if (*n_a > n_for + n_rev)
{
memmove(a + n_for, a + (*n_a) - n_rev, n_rev * sizeof(mm128_t));
*n_a = n_for + n_rev;
}
#ifdef ANALYSIS_PERF
tmp_diff = get_mseconds() - tmp_cur_time;
__sync_fetch_and_add(&time_collect_seed_hits_heap, tmp_diff);
#endif
return a;
}
@ -171,31 +230,48 @@ static mm128_t *collect_seed_hits(void *km, const mm_mapopt_t *opt, int max_occ,
int i, n_m;
mm_seed_t *m;
mm128_t *a;
#ifdef ANALYSIS_PERF
int64_t tmp_cur_time = get_mseconds(), tmp_diff = 0;
#endif
m = mm_collect_matches(km, &n_m, qlen, max_occ, opt->max_max_occ, opt->occ_dist, mi, mv, n_a, rep_len, n_mini_pos, mini_pos);
#ifdef ANALYSIS_PERF
tmp_diff = get_mseconds() - tmp_cur_time;
__sync_fetch_and_add(&time_collect_seed_hits, tmp_diff);
#endif
a = (mm128_t *)kmalloc(km, *n_a * sizeof(mm128_t));
for (i = 0, *n_a = 0; i < n_m; ++i) {
for (i = 0, *n_a = 0; i < n_m; ++i)
{
mm_seed_t *q = &m[i];
const uint64_t *r = q->cr;
uint32_t k;
for (k = 0; k < q->n; ++k) {
for (k = 0; k < q->n; ++k)
{
int32_t is_self, rpos = (uint32_t)r[k] >> 1;
mm128_t *p;
if (skip_seed(opt->flag, r[k], q, qname, qlen, mi, &is_self)) continue;
if (skip_seed(opt->flag, r[k], q, qname, qlen, mi, &is_self))
continue;
p = &a[(*n_a)++];
if ((r[k]&1) == (q->q_pos&1)) { // forward strand
if ((r[k] & 1) == (q->q_pos & 1))
{ // forward strand
p->x = (r[k] & 0xffffffff00000000ULL) | rpos;
p->y = (uint64_t)q->q_span << 32 | q->q_pos >> 1;
} else if (!(opt->flag & MM_F_QSTRAND)) { // reverse strand and not in the query-strand mode
}
else if (!(opt->flag & MM_F_QSTRAND))
{ // reverse strand and not in the query-strand mode
p->x = 1ULL << 63 | (r[k] & 0xffffffff00000000ULL) | rpos;
p->y = (uint64_t)q->q_span << 32 | (qlen - ((q->q_pos >> 1) + 1 - q->q_span) - 1);
} else { // reverse strand; query-strand
}
else
{ // reverse strand; query-strand
int32_t len = mi->seq[r[k] >> 32].len;
p->x = 1ULL << 63 | (r[k] & 0xffffffff00000000ULL) | (len - (rpos + 1 - q->q_span) - 1); // coordinate only accurate for non-HPC seeds
p->y = (uint64_t)q->q_span << 32 | q->q_pos >> 1;
}
p->y |= (uint64_t)q->seg_id << MM_SEED_SEG_SHIFT;
if (q->is_tandem) p->y |= MM_SEED_TANDEM;
if (is_self) p->y |= MM_SEED_SELF;
if (q->is_tandem)
p->y |= MM_SEED_TANDEM;
if (is_self)
p->y |= MM_SEED_SELF;
}
}
kfree(km, m);
@ -205,18 +281,23 @@ static mm128_t *collect_seed_hits(void *km, const mm_mapopt_t *opt, int max_occ,
static void chain_post(const mm_mapopt_t *opt, int max_chain_gap_ref, const mm_idx_t *mi, void *km, int qlen, int n_segs, const int *qlens, int *n_regs, mm_reg1_t *regs, mm128_t *a)
{
if (!(opt->flag & MM_F_ALL_CHAINS)) { // don't choose primary mapping(s)
if (!(opt->flag & MM_F_ALL_CHAINS))
{ // don't choose primary mapping(s)
mm_set_parent(km, opt->mask_level, opt->mask_len, *n_regs, regs, opt->a * 2 + opt->b, opt->flag & MM_F_HARD_MLEVEL, opt->alt_drop);
if (n_segs <= 1) mm_select_sub(km, opt->pri_ratio, mi->k*2, opt->best_n, 1, opt->max_gap * 0.8, n_regs, regs);
else mm_select_sub_multi(km, opt->pri_ratio, 0.2f, 0.7f, max_chain_gap_ref, mi->k*2, opt->best_n, n_segs, qlens, n_regs, regs);
if (n_segs <= 1)
mm_select_sub(km, opt->pri_ratio, mi->k * 2, opt->best_n, 1, opt->max_gap * 0.8, n_regs, regs);
else
mm_select_sub_multi(km, opt->pri_ratio, 0.2f, 0.7f, max_chain_gap_ref, mi->k * 2, opt->best_n, n_segs, qlens, n_regs, regs);
}
}
static mm_reg1_t *align_regs(const mm_mapopt_t *opt, const mm_idx_t *mi, void *km, int qlen, const char *seq, int *n_regs, mm_reg1_t *regs, mm128_t *a)
{
if (!(opt->flag & MM_F_CIGAR)) return regs;
if (!(opt->flag & MM_F_CIGAR))
return regs;
regs = mm_align_skeleton(km, opt, mi, qlen, seq, n_regs, regs, a); // this calls mm_filter_regs()
if (!(opt->flag & MM_F_ALL_CHAINS)) { // don't choose primary mapping(s)
if (!(opt->flag & MM_F_ALL_CHAINS))
{ // don't choose primary mapping(s)
mm_set_parent(km, opt->mask_level, opt->mask_len, *n_regs, regs, opt->a * 2 + opt->b, opt->flag & MM_F_HARD_MLEVEL, opt->alt_drop);
mm_select_sub(km, opt->pri_ratio, mi->k * 2, opt->best_n, 0, opt->max_gap * 0.8, n_regs, regs);
mm_set_sam_pri(*n_regs, regs);
@ -236,66 +317,105 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
mm_reg1_t *regs0;
km_stat_t kmst;
float chn_pen_gap, chn_pen_skip;
#ifdef ANALYSIS_PERF
int64_t tmp_cur_time = get_mseconds(), tmp_diff = 0;
#endif
for (i = 0, qlen_sum = 0; i < n_segs; ++i)
qlen_sum += qlens[i], n_regs[i] = 0, regs[i] = 0;
if (qlen_sum == 0 || n_segs <= 0 || n_segs > MM_MAX_SEG) return;
if (opt->max_qlen > 0 && qlen_sum > opt->max_qlen) return;
if (qlen_sum == 0 || n_segs <= 0 || n_segs > MM_MAX_SEG)
return;
if (opt->max_qlen > 0 && qlen_sum > opt->max_qlen)
return;
hash = qname && !(opt->flag & MM_F_NO_HASH_NAME) ? __ac_X31_hash_string(qname) : 0;
hash ^= __ac_Wang_hash(qlen_sum) + __ac_Wang_hash(opt->seed);
hash = __ac_Wang_hash(hash);
collect_minimizers(b->km, opt, mi, n_segs, qlens, seqs, &mv);
if (opt->q_occ_frac > 0.0f) mm_seed_mz_flt(b->km, &mv, opt->mid_occ, opt->q_occ_frac);
if (opt->flag & MM_F_HEAP_SORT) a = collect_seed_hits_heap(b->km, opt, opt->mid_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
else a = collect_seed_hits(b->km, opt, opt->mid_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
#ifdef ANALYSIS_PERF
tmp_diff = get_mseconds() - tmp_cur_time;
__sync_fetch_and_add(&time_mm_map_frag_b1, tmp_diff);
tmp_cur_time = get_mseconds();
#endif
if (opt->q_occ_frac > 0.0f)
mm_seed_mz_flt(b->km, &mv, opt->mid_occ, opt->q_occ_frac);
if (opt->flag & MM_F_HEAP_SORT)
a = collect_seed_hits_heap(b->km, opt, opt->mid_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
else
a = collect_seed_hits(b->km, opt, opt->mid_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
if (mm_dbg_flag & MM_DBG_PRINT_SEED) {
if (mm_dbg_flag & MM_DBG_PRINT_SEED)
{
fprintf(stderr, "RS\t%d\n", rep_len);
for (i = 0; i < n_a; ++i)
fprintf(stderr, "SD\t%s\t%d\t%c\t%d\t%d\t%d\n", mi->seq[a[i].x << 1 >> 33].name, (int32_t)a[i].x, "+-"[a[i].x >> 63], (int32_t)a[i].y, (int32_t)(a[i].y >> 32 & 0xff),
i == 0 ? 0 : ((int32_t)a[i].y - (int32_t)a[i - 1].y) - ((int32_t)a[i].x - (int32_t)a[i - 1].x));
}
#ifdef ANALYSIS_PERF
tmp_diff = get_mseconds() - tmp_cur_time;
__sync_fetch_and_add(&time_mm_map_frag_b2, tmp_diff);
tmp_cur_time = get_mseconds();
#endif
// set max chaining gap on the query and the reference sequence
if (is_sr)
max_chain_gap_qry = qlen_sum > opt->max_gap ? qlen_sum : opt->max_gap;
else max_chain_gap_qry = opt->max_gap;
if (opt->max_gap_ref > 0) {
else
max_chain_gap_qry = opt->max_gap;
if (opt->max_gap_ref > 0)
{
max_chain_gap_ref = opt->max_gap_ref; // always honor mm_mapopt_t::max_gap_ref if set
} else if (opt->max_frag_len > 0) {
}
else if (opt->max_frag_len > 0)
{
max_chain_gap_ref = opt->max_frag_len - qlen_sum;
if (max_chain_gap_ref < opt->max_gap) max_chain_gap_ref = opt->max_gap;
} else max_chain_gap_ref = opt->max_gap;
if (max_chain_gap_ref < opt->max_gap)
max_chain_gap_ref = opt->max_gap;
}
else
max_chain_gap_ref = opt->max_gap;
chn_pen_gap = opt->chain_gap_scale * 0.01 * mi->k;
chn_pen_skip = opt->chain_skip_scale * 0.01 * mi->k;
if (opt->flag & MM_F_RMQ) {
if (opt->flag & MM_F_RMQ)
{
a = mg_lchain_rmq(opt->max_gap, opt->rmq_inner_dist, opt->bw, opt->max_chain_skip, opt->rmq_size_cap, opt->min_cnt, opt->min_chain_score,
chn_pen_gap, chn_pen_skip, n_a, a, &n_regs0, &u, b->km);
} else {
}
else
{
a = mg_lchain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->max_chain_iter, opt->min_cnt, opt->min_chain_score,
chn_pen_gap, chn_pen_skip, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
}
if (opt->bw_long > opt->bw && (opt->flag & (MM_F_SPLICE|MM_F_SR|MM_F_NO_LJOIN)) == 0 && n_segs == 1 && n_regs0 > 1) { // re-chain/long-join for long sequences
#ifdef ANALYSIS_PERF
tmp_diff = get_mseconds() - tmp_cur_time;
__sync_fetch_and_add(&time_mm_map_frag_b3, tmp_diff);
tmp_cur_time = get_mseconds();
#endif
if (opt->bw_long > opt->bw && (opt->flag & (MM_F_SPLICE | MM_F_SR | MM_F_NO_LJOIN)) == 0 && n_segs == 1 && n_regs0 > 1)
{ // re-chain/long-join for long sequences
int32_t st = (int32_t)a[0].y, en = (int32_t)a[(int32_t)u[0] - 1].y;
if (qlen_sum - (en - st) > opt->rmq_rescue_size || en - st > qlen_sum * opt->rmq_rescue_ratio) {
if (qlen_sum - (en - st) > opt->rmq_rescue_size || en - st > qlen_sum * opt->rmq_rescue_ratio)
{
int32_t i;
for (i = 0, n_a = 0; i < n_regs0; ++i) n_a += (int32_t)u[i];
for (i = 0, n_a = 0; i < n_regs0; ++i)
n_a += (int32_t)u[i];
kfree(b->km, u);
radix_sort_128x(a, a + n_a);
a = mg_lchain_rmq(opt->max_gap, opt->rmq_inner_dist, opt->bw_long, opt->max_chain_skip, opt->rmq_size_cap, opt->min_cnt, opt->min_chain_score,
chn_pen_gap, chn_pen_skip, n_a, a, &n_regs0, &u, b->km);
}
} else if (opt->max_occ > opt->mid_occ && rep_len > 0 && !(opt->flag & MM_F_RMQ)) { // re-chain, mostly for short reads
}
else if (opt->max_occ > opt->mid_occ && rep_len > 0 && !(opt->flag & MM_F_RMQ))
{ // re-chain, mostly for short reads
int rechain = 0;
if (n_regs0 > 0) { // test if the best chain has all the segments
if (n_regs0 > 0)
{ // test if the best chain has all the segments
int n_chained_segs = 1, max = 0, max_i = -1, max_off = -1, off = 0;
for (i = 0; i < n_regs0; ++i) { // find the best chain
if (max < (int)(u[i]>>32)) max = u[i]>>32, max_i = i, max_off = off;
for (i = 0; i < n_regs0; ++i)
{ // find the best chain
if (max < (int)(u[i] >> 32))
max = u[i] >> 32, max_i = i, max_off = off;
off += (uint32_t)u[i];
}
for (i = 1; i < (int32_t)u[max_i]; ++i) // count the number of segments in the best chain
@ -303,22 +423,32 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
++n_chained_segs;
if (n_chained_segs < n_segs)
rechain = 1;
} else rechain = 1;
if (rechain) { // redo chaining with a higher max_occ threshold
}
else
rechain = 1;
if (rechain)
{ // redo chaining with a higher max_occ threshold
kfree(b->km, a);
kfree(b->km, u);
kfree(b->km, mini_pos);
if (opt->flag & MM_F_HEAP_SORT) a = collect_seed_hits_heap(b->km, opt, opt->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
else a = collect_seed_hits(b->km, opt, opt->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
if (opt->flag & MM_F_HEAP_SORT)
a = collect_seed_hits_heap(b->km, opt, opt->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
else
a = collect_seed_hits(b->km, opt, opt->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
a = mg_lchain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->max_chain_iter, opt->min_cnt, opt->min_chain_score,
chn_pen_gap, chn_pen_skip, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
}
}
b->frag_gap = max_chain_gap_ref;
b->rep_len = rep_len;
#ifdef ANALYSIS_PERF
tmp_diff = get_mseconds() - tmp_cur_time;
__sync_fetch_and_add(&time_mm_map_frag_b4, tmp_diff);
tmp_cur_time = get_mseconds();
#endif
regs0 = mm_gen_regs(b->km, hash, qlen_sum, n_regs0, u, a, !!(opt->flag & MM_F_QSTRAND));
if (mi->n_alt) {
if (mi->n_alt)
{
mm_mark_alt(mi, n_regs0, regs0);
mm_hit_sort(b->km, &n_regs0, regs0, opt->alt_drop); // this step can be merged into mm_gen_regs(); will do if this shows up in profile
}
@ -330,21 +460,30 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
i == regs0[j].as ? 0 : ((int32_t)a[i].y - (int32_t)a[i - 1].y) - ((int32_t)a[i].x - (int32_t)a[i - 1].x));
chain_post(opt, max_chain_gap_ref, mi, b->km, qlen_sum, n_segs, qlens, &n_regs0, regs0, a);
if (!is_sr && !(opt->flag&MM_F_QSTRAND)) {
if (!is_sr && !(opt->flag & MM_F_QSTRAND))
{
mm_est_err(mi, qlen_sum, n_regs0, regs0, a, n_mini_pos, mini_pos);
n_regs0 = mm_filter_strand_retained(n_regs0, regs0);
}
if (n_segs == 1) { // uni-segment
#ifdef ANALYSIS_PERF
tmp_diff = get_mseconds() - tmp_cur_time;
__sync_fetch_and_add(&time_mm_map_frag_b5, tmp_diff);
tmp_cur_time = get_mseconds();
#endif
if (n_segs == 1)
{ // uni-segment
regs0 = align_regs(opt, mi, b->km, qlens[0], seqs[0], &n_regs0, regs0, a);
regs0 = (mm_reg1_t *)realloc(regs0, sizeof(*regs0) * n_regs0);
mm_set_mapq(b->km, n_regs0, regs0, opt->min_chain_score, opt->a, rep_len, is_sr);
n_regs[0] = n_regs0, regs[0] = regs0;
} else { // multi-segment
}
else
{ // multi-segment
mm_seg_t *seg;
seg = mm_seg_gen(b->km, hash, n_segs, qlens, n_regs0, regs0, n_regs, regs, a); // split fragment chain to separate segment chains
free(regs0);
for (i = 0; i < n_segs; ++i) {
for (i = 0; i < n_segs; ++i)
{
mm_set_parent(b->km, opt->mask_level, opt->mask_len, n_regs[i], regs[i], opt->a * 2 + opt->b, opt->flag & MM_F_HARD_MLEVEL, opt->alt_drop); // update mm_reg1_t::parent
regs[i] = align_regs(opt, mi, b->km, qlens[i], seqs[i], &n_regs[i], regs[i], seg[i].a);
mm_set_mapq(b->km, n_regs[i], regs[i], opt->min_chain_score, opt->a, rep_len, is_sr);
@ -353,18 +492,24 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
if (n_segs == 2 && opt->pe_ori >= 0 && (opt->flag & MM_F_CIGAR))
mm_pair(b->km, max_chain_gap_ref, opt->pe_bonus, opt->a * 2 + opt->b, opt->a, qlens, n_regs, regs); // pairing
}
#ifdef ANALYSIS_PERF
tmp_diff = get_mseconds() - tmp_cur_time;
__sync_fetch_and_add(&time_mm_map_frag_b6, tmp_diff);
tmp_cur_time = get_mseconds();
#endif
kfree(b->km, mv.a);
kfree(b->km, a);
kfree(b->km, u);
kfree(b->km, mini_pos);
if (b->km) {
if (b->km)
{
km_stat(b->km, &kmst);
if (mm_dbg_flag & MM_DBG_PRINT_QNAME)
fprintf(stderr, "QM\t%s\t%d\tcap=%ld,nCore=%ld,largest=%ld\n", qname, qlen_sum, kmst.capacity, kmst.n_cores, kmst.largest);
assert(kmst.n_blocks == kmst.n_cores); // otherwise, there is a memory leak
if (kmst.largest > 1U<<28 || (opt->cap_kalloc > 0 && kmst.capacity > opt->cap_kalloc)) {
if (kmst.largest > 1U << 28 || (opt->cap_kalloc > 0 && kmst.capacity > opt->cap_kalloc))
{
if (mm_dbg_flag & MM_DBG_PRINT_QNAME)
fprintf(stderr, "[W::%s] reset thread-local memory after read %s\n", __func__, qname);
km_destroy(b->km);
@ -384,7 +529,8 @@ mm_reg1_t *mm_map(const mm_idx_t *mi, int qlen, const char *seq, int *n_regs, mm
* Multi-threaded mapping *
**************************/
typedef struct {
typedef struct
{
int n_processed, n_threads, n_fp;
int64_t mini_batch_size;
const mm_mapopt_t *opt;
@ -397,7 +543,8 @@ typedef struct {
FILE *fp_split, **fp_parts;
} pipeline_t;
typedef struct {
typedef struct
{
const pipeline_t *p;
int n_seq, n_frag;
mm_bseq1_t *seq;
@ -414,34 +561,56 @@ static void worker_for(void *_data, long i, int tid) // kt_for() callback
double t = 0.0;
mm_tbuf_t *b = s->buf[tid];
assert(s->n_seg[i] <= MM_MAX_SEG);
if (mm_dbg_flag & MM_DBG_PRINT_QNAME) {
if (mm_dbg_flag & MM_DBG_PRINT_QNAME)
{
fprintf(stderr, "QR\t%s\t%d\t%d\n", s->seq[off].name, tid, s->seq[off].l_seq);
t = realtime();
}
for (j = 0; j < s->n_seg[i]; ++j) {
#ifdef ANALYSIS_PERF
int64_t tmp_cur_time = get_mseconds(), tmp_diff = 0;
#endif
for (j = 0; j < s->n_seg[i]; ++j)
{
if (s->n_seg[i] == 2 && ((j == 0 && (pe_ori >> 1 & 1)) || (j == 1 && (pe_ori & 1))))
mm_revcomp_bseq(&s->seq[off + j]);
qlens[j] = s->seq[off + j].l_seq;
qseqs[j] = s->seq[off + j].seq;
}
if (s->p->opt->flag & MM_F_INDEPEND_SEG) {
for (j = 0; j < s->n_seg[i]; ++j) {
#ifdef ANALYSIS_PERF
tmp_diff = get_mseconds() - tmp_cur_time;
__sync_fetch_and_add(&time_map_work_for_block_1, tmp_diff);
tmp_cur_time = get_mseconds();
#endif
if (s->p->opt->flag & MM_F_INDEPEND_SEG)
{
for (j = 0; j < s->n_seg[i]; ++j)
{
mm_map_frag(s->p->mi, 1, &qlens[j], &qseqs[j], &s->n_reg[off + j], &s->reg[off + j], b, s->p->opt, s->seq[off + j].name);
s->rep_len[off + j] = b->rep_len;
s->frag_gap[off + j] = b->frag_gap;
}
} else {
}
else
{
mm_map_frag(s->p->mi, s->n_seg[i], qlens, qseqs, &s->n_reg[off], &s->reg[off], b, s->p->opt, s->seq[off].name);
for (j = 0; j < s->n_seg[i]; ++j) {
for (j = 0; j < s->n_seg[i]; ++j)
{
s->rep_len[off + j] = b->rep_len;
s->frag_gap[off + j] = b->frag_gap;
}
}
#ifdef ANALYSIS_PERF
tmp_diff = get_mseconds() - tmp_cur_time;
__sync_fetch_and_add(&time_map_work_for_block_2, tmp_diff);
tmp_cur_time = get_mseconds();
#endif
for (j = 0; j < s->n_seg[i]; ++j) // flip the query strand and coordinate to the original read strand
if (s->n_seg[i] == 2 && ((j == 0 && (pe_ori>>1&1)) || (j == 1 && (pe_ori&1)))) {
if (s->n_seg[i] == 2 && ((j == 0 && (pe_ori >> 1 & 1)) || (j == 1 && (pe_ori & 1))))
{
int k, t;
mm_revcomp_bseq(&s->seq[off + j]);
for (k = 0; k < s->n_reg[off + j]; ++k) {
for (k = 0; k < s->n_reg[off + j]; ++k)
{
mm_reg1_t *r = &s->reg[off + j][k];
t = r->qs;
r->qs = qlens[j] - r->qe;
@ -449,6 +618,10 @@ static void worker_for(void *_data, long i, int tid) // kt_for() callback
r->rev = !r->rev;
}
}
#ifdef ANALYSIS_PERF
tmp_diff = get_mseconds() - tmp_cur_time;
__sync_fetch_and_add(&time_map_work_for_block_3, tmp_diff);
#endif
if (mm_dbg_flag & MM_DBG_PRINT_QNAME)
fprintf(stderr, "QT\t%s\t%d\t%.6f\n", s->seq[off].name, tid, realtime() - t);
}
@ -467,12 +640,15 @@ static void merge_hits(step_t *s)
n_reg_part = qlens + max_seg;
rep_len_part = n_reg_part + s->p->n_parts;
frag_gap_part = rep_len_part + s->p->n_parts;
for (f = 0, k = k0 = 0; f < s->n_frag; ++f) {
for (f = 0, k = k0 = 0; f < s->n_frag; ++f)
{
k0 = k;
for (i = 0; i < s->n_seg[f]; ++i, ++k) {
for (i = 0; i < s->n_seg[f]; ++i, ++k)
{
int j, l, t, rep_len = 0;
qlens[i] = s->seq[k].l_seq;
for (j = 0, s->n_reg[k] = 0; j < s->p->n_parts; ++j) {
for (j = 0, s->n_reg[k] = 0; j < s->p->n_parts; ++j)
{
mm_err_fread(&n_reg_part[j], sizeof(int), 1, fp[j]);
mm_err_fread(&rep_len_part[j], sizeof(int), 1, fp[j]);
mm_err_fread(&frag_gap_part[j], sizeof(int), 1, fp[j]);
@ -481,13 +657,16 @@ static void merge_hits(step_t *s)
rep_len = rep_len_part[j];
}
s->reg[k] = CALLOC(mm_reg1_t, s->n_reg[k]);
for (j = 0, l = 0; j < s->p->n_parts; ++j) {
for (t = 0; t < n_reg_part[j]; ++t, ++l) {
for (j = 0, l = 0; j < s->p->n_parts; ++j)
{
for (t = 0; t < n_reg_part[j]; ++t, ++l)
{
mm_reg1_t *r = &s->reg[k][l];
uint32_t capacity;
mm_err_fread(r, sizeof(mm_reg1_t), 1, fp[j]);
r->rid += s->p->rid_shift[j];
if (opt->flag & MM_F_CIGAR) {
if (opt->flag & MM_F_CIGAR)
{
mm_err_fread(&capacity, 4, 1, fp[j]);
r->p = (mm_extra_t *)calloc(capacity, 4);
r->p->capacity = capacity;
@ -497,15 +676,18 @@ static void merge_hits(step_t *s)
}
if (!(opt->flag & MM_F_SR) && s->seq[k].l_seq >= opt->rank_min_len)
mm_update_dp_max(s->seq[k].l_seq, s->n_reg[k], s->reg[k], opt->rank_frac, opt->a, opt->b);
for (j = 0; j < s->n_reg[k]; ++j) {
for (j = 0; j < s->n_reg[k]; ++j)
{
mm_reg1_t *r = &s->reg[k][j];
if (r->p) r->p->dp_max2 = 0; // reset ->dp_max2 as mm_set_parent() doesn't clear it; necessary with mm_update_dp_max()
if (r->p)
r->p->dp_max2 = 0; // reset ->dp_max2 as mm_set_parent() doesn't clear it; necessary with mm_update_dp_max()
r->subsc = 0; // this may not be necessary
r->n_sub = 0; // n_sub will be an underestimate as we don't see all the chains now, but it can't be accurate anyway
}
mm_hit_sort(km, &s->n_reg[k], s->reg[k], opt->alt_drop);
mm_set_parent(km, opt->mask_level, opt->mask_len, s->n_reg[k], s->reg[k], opt->a * 2 + opt->b, opt->flag & MM_F_HARD_MLEVEL, opt->alt_drop);
if (!(opt->flag & MM_F_ALL_CHAINS)) {
if (!(opt->flag & MM_F_ALL_CHAINS))
{
mm_select_sub(km, opt->pri_ratio, s->p->mi->k * 2, opt->best_n, 0, opt->max_gap * 0.8, &s->n_reg[k], s->reg[k]);
mm_set_sam_pri(s->n_reg[k], s->reg[k]);
}
@ -522,15 +704,19 @@ static void *worker_pipeline(void *shared, int step, void *in)
{
int i, j, k;
pipeline_t *p = (pipeline_t *)shared;
if (step == 0) { // step 0: read sequences
if (step == 0)
{ // step 0: read sequences
int with_qual = (!!(p->opt->flag & MM_F_OUT_SAM) && !(p->opt->flag & MM_F_NO_QUAL));
int with_comment = !!(p->opt->flag & MM_F_COPY_COMMENT);
int frag_mode = (p->n_fp > 1 || !!(p->opt->flag & MM_F_FRAG_MODE));
step_t *s;
s = (step_t *)calloc(1, sizeof(step_t));
if (p->n_fp > 1) s->seq = mm_bseq_read_frag2(p->n_fp, p->fp, p->mini_batch_size, with_qual, with_comment, &s->n_seq);
else s->seq = mm_bseq_read3(p->fp[0], p->mini_batch_size, with_qual, with_comment, frag_mode, &s->n_seq);
if (s->seq) {
if (p->n_fp > 1)
s->seq = mm_bseq_read_frag2(p->n_fp, p->fp, p->mini_batch_size, with_qual, with_comment, &s->n_seq);
else
s->seq = mm_bseq_read3(p->fp[0], p->mini_batch_size, with_qual, with_comment, frag_mode, &s->n_seq);
if (s->seq)
{
s->p = p;
for (i = 0; i < s->n_seq; ++i)
s->seq[i].rid = p->n_processed++;
@ -544,42 +730,61 @@ static void *worker_pipeline(void *shared, int step, void *in)
s->frag_gap = s->rep_len + s->n_seq;
s->reg = (mm_reg1_t **)calloc(s->n_seq, sizeof(mm_reg1_t *));
for (i = 1, j = 0; i <= s->n_seq; ++i)
if (i == s->n_seq || !frag_mode || !mm_qname_same(s->seq[i-1].name, s->seq[i].name)) {
if (i == s->n_seq || !frag_mode || !mm_qname_same(s->seq[i - 1].name, s->seq[i].name))
{
s->n_seg[s->n_frag] = i - j;
s->seg_off[s->n_frag++] = j;
j = i;
}
return s;
} else free(s);
} else if (step == 1) { // step 1: map
if (p->n_parts > 0) merge_hits((step_t*)in);
else kt_for(p->n_threads, worker_for, in, ((step_t*)in)->n_frag);
}
else
free(s);
}
else if (step == 1)
{ // step 1: map
if (p->n_parts > 0)
merge_hits((step_t *)in);
else
kt_for(p->n_threads, worker_for, in, ((step_t *)in)->n_frag);
return in;
} else if (step == 2) { // step 2: output
}
else if (step == 2)
{ // step 2: output
void *km = 0;
step_t *s = (step_t *)in;
const mm_idx_t *mi = p->mi;
for (i = 0; i < p->n_threads; ++i) mm_tbuf_destroy(s->buf[i]);
for (i = 0; i < p->n_threads; ++i)
mm_tbuf_destroy(s->buf[i]);
free(s->buf);
if ((p->opt->flag & MM_F_OUT_CS) && !(mm_dbg_flag & MM_DBG_NO_KALLOC)) km = km_init();
for (k = 0; k < s->n_frag; ++k) {
if ((p->opt->flag & MM_F_OUT_CS) && !(mm_dbg_flag & MM_DBG_NO_KALLOC))
km = km_init();
for (k = 0; k < s->n_frag; ++k)
{
int seg_st = s->seg_off[k], seg_en = s->seg_off[k] + s->n_seg[k];
for (i = seg_st; i < seg_en; ++i) {
for (i = seg_st; i < seg_en; ++i)
{
mm_bseq1_t *t = &s->seq[i];
if (p->opt->split_prefix && p->n_parts == 0) { // then write to temporary files
if (p->opt->split_prefix && p->n_parts == 0)
{ // then write to temporary files
mm_err_fwrite(&s->n_reg[i], sizeof(int), 1, p->fp_split);
mm_err_fwrite(&s->rep_len[i], sizeof(int), 1, p->fp_split);
mm_err_fwrite(&s->frag_gap[i], sizeof(int), 1, p->fp_split);
for (j = 0; j < s->n_reg[i]; ++j) {
for (j = 0; j < s->n_reg[i]; ++j)
{
mm_reg1_t *r = &s->reg[i][j];
mm_err_fwrite(r, sizeof(mm_reg1_t), 1, p->fp_split);
if (p->opt->flag & MM_F_CIGAR) {
if (p->opt->flag & MM_F_CIGAR)
{
mm_err_fwrite(&r->p->capacity, 4, 1, p->fp_split);
mm_err_fwrite(r->p, r->p->capacity, 4, p->fp_split);
}
}
} else if (s->n_reg[i] > 0) { // the query has at least one hit
for (j = 0; j < s->n_reg[i]; ++j) {
}
else if (s->n_reg[i] > 0)
{ // the query has at least one hit
for (j = 0; j < s->n_reg[i]; ++j)
{
mm_reg1_t *r = &s->reg[i][j];
assert(!r->sam_pri || r->id == r->parent);
if ((p->opt->flag & MM_F_NO_PRINT_2ND) && r->id != r->parent)
@ -590,7 +795,9 @@ static void *worker_pipeline(void *shared, int step, void *in)
mm_write_paf3(&p->str, mi, t, r, km, p->opt->flag, s->rep_len[i]);
mm_err_puts(p->str.s);
}
} else if ((p->opt->flag & MM_F_PAF_NO_HIT) || ((p->opt->flag & MM_F_OUT_SAM) && !(p->opt->flag & MM_F_SAM_HIT_ONLY))) { // output an empty hit, if requested
}
else if ((p->opt->flag & MM_F_PAF_NO_HIT) || ((p->opt->flag & MM_F_OUT_SAM) && !(p->opt->flag & MM_F_SAM_HIT_ONLY)))
{ // output an empty hit, if requested
if (p->opt->flag & MM_F_OUT_SAM)
mm_write_sam3(&p->str, mi, t, i - seg_st, -1, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t *const *)&s->reg[seg_st], km, p->opt->flag, s->rep_len[i]);
else
@ -598,15 +805,22 @@ static void *worker_pipeline(void *shared, int step, void *in)
mm_err_puts(p->str.s);
}
}
for (i = seg_st; i < seg_en; ++i) {
for (j = 0; j < s->n_reg[i]; ++j) free(s->reg[i][j].p);
for (i = seg_st; i < seg_en; ++i)
{
for (j = 0; j < s->n_reg[i]; ++j)
free(s->reg[i][j].p);
free(s->reg[i]);
free(s->seq[i].seq); free(s->seq[i].name);
if (s->seq[i].qual) free(s->seq[i].qual);
if (s->seq[i].comment) free(s->seq[i].comment);
free(s->seq[i].seq);
free(s->seq[i].name);
if (s->seq[i].qual)
free(s->seq[i].qual);
if (s->seq[i].comment)
free(s->seq[i].comment);
}
}
free(s->reg); free(s->n_reg); free(s->seq); // seg_off, n_seg, rep_len and frag_gap were allocated with reg; no memory leak here
free(s->reg);
free(s->n_reg);
free(s->seq); // seg_off, n_seg, rep_len and frag_gap were allocated with reg; no memory leak here
km_destroy(km);
if (mm_verbose >= 3)
fprintf(stderr, "[M::%s::%.3f*%.2f] mapped %d sequences\n", __func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), s->n_seq);
@ -620,8 +834,10 @@ static mm_bseq_file_t **open_bseqs(int n, const char **fn)
mm_bseq_file_t **fp;
int i, j;
fp = (mm_bseq_file_t **)calloc(n, sizeof(mm_bseq_file_t *));
for (i = 0; i < n; ++i) {
if ((fp[i] = mm_bseq_open(fn[i])) == 0) {
for (i = 0; i < n; ++i)
{
if ((fp[i] = mm_bseq_open(fn[i])) == 0)
{
if (mm_verbose >= 1)
fprintf(stderr, "ERROR: failed to open file '%s': %s\n", fn[i], strerror(errno));
for (j = 0; j < i; ++j)
@ -635,26 +851,36 @@ static mm_bseq_file_t **open_bseqs(int n, const char **fn)
int mm_map_file_frag(const mm_idx_t *idx, int n_segs, const char **fn, const mm_mapopt_t *opt, int n_threads)
{
#ifdef ANALYSIS_PERF
int64_t tmp_cur_time = get_mseconds();
#endif
int i, pl_threads;
pipeline_t pl;
if (n_segs < 1) return -1;
if (n_segs < 1)
return -1;
memset(&pl, 0, sizeof(pipeline_t));
pl.n_fp = n_segs;
pl.fp = open_bseqs(pl.n_fp, fn);
if (pl.fp == 0) return -1;
if (pl.fp == 0)
return -1;
pl.opt = opt, pl.mi = idx;
pl.n_threads = n_threads > 1 ? n_threads : 1;
pl.mini_batch_size = opt->mini_batch_size;
if (opt->split_prefix)
pl.fp_split = mm_split_init(opt->split_prefix, idx);
pl_threads = n_threads == 1? 1 : (opt->flag&MM_F_2_IO_THREADS)? 3 : 2;
pl_threads = n_threads == 1 ? 1 : (opt->flag & MM_F_2_IO_THREADS) ? 3
: 2;
kt_pipeline(pl_threads, worker_pipeline, &pl, 3);
free(pl.str.s);
if (pl.fp_split) fclose(pl.fp_split);
if (pl.fp_split)
fclose(pl.fp_split);
for (i = 0; i < pl.n_fp; ++i)
mm_bseq_close(pl.fp[i]);
free(pl.fp);
#ifdef ANALYSIS_PERF
time_mm_map_file_frag += get_mseconds() - tmp_cur_time;
#endif
return 0;
}
@ -668,11 +894,13 @@ int mm_split_merge(int n_segs, const char **fn, const mm_mapopt_t *opt, int n_sp
int i;
pipeline_t pl;
mm_idx_t *mi;
if (n_segs < 1 || n_split_idx < 1) return -1;
if (n_segs < 1 || n_split_idx < 1)
return -1;
memset(&pl, 0, sizeof(pipeline_t));
pl.n_fp = n_segs;
pl.fp = open_bseqs(pl.n_fp, fn);
if (pl.fp == 0) return -1;
if (pl.fp == 0)
return -1;
pl.opt = opt;
pl.mini_batch_size = opt->mini_batch_size;
@ -680,7 +908,8 @@ int mm_split_merge(int n_segs, const char **fn, const mm_mapopt_t *opt, int n_sp
pl.fp_parts = CALLOC(FILE *, pl.n_parts);
pl.rid_shift = CALLOC(uint32_t, pl.n_parts);
pl.mi = mi = mm_split_merge_prep(opt->split_prefix, n_split_idx, pl.fp_parts, pl.rid_shift);
if (pl.mi == 0) {
if (pl.mi == 0)
{
free(pl.fp_parts);
free(pl.rid_shift);
return -1;

3
minimap.h
View File

@ -7,6 +7,9 @@
#define MM_VERSION "2.26-r1175"
// 用来开关调试性能分析，运行时间等信息
#define ANALYSIS_PERF 1
#define MM_F_NO_DIAG (0x001LL) // no exact diagonal hit
#define MM_F_NO_DUAL (0x002LL) // skip pairs where query name is lexicographically larger than target name
#define MM_F_CIGAR (0x004LL)

2
sketch.c
View File

@ -106,7 +106,7 @@ void mm_sketch(void *km, const char *str, int len, int w, int k, uint32_t rid, i
kmer[0] = (kmer[0] << 2 | c) & mask; // forward k-mer
kmer[1] = (kmer[1] >> 2) | (3ULL^c) << shift1; // reverse k-mer
if (kmer[0] == kmer[1]) continue; // skip "symmetric k-mers" as we don't know it strand
z = kmer[0] < kmer[1]? 0 : 1; // strand
z = kmer[0] < kmer[1]? 0 : 1; // strand // kmer的strand到底是什么意思？为什么通过比较就能确定正反？
++l;
if (l >= k && kmer_span < 256) {
info.x = hash64(kmer[z], mask) << 8 | kmer_span;