From 8140259974cf8f6f16d00527e8d660446e2ce233 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Sun, 10 Mar 2024 13:47:34 -0400 Subject: [PATCH] r1183: added lr:hq; fixed transition * Added the lr:hq preset suggested by Nanopore developers (#1127) * Fixed transition scoring. It did not work with presets. * Cleaned up preset documentation --- align.c | 16 ++++++++------- main.c | 10 ++++----- minimap.h | 2 +- minimap2.1 | 59 +++++++++++++++++++++++++++++++++++++++++++++--------- options.c | 15 +++++++------- 5 files changed, 73 insertions(+), 29 deletions(-) diff --git a/align.c b/align.c index f13d264..5a0af66 100644 --- a/align.c +++ b/align.c @@ -23,13 +23,14 @@ static void ksw_gen_simple_mat(int m, int8_t *mat, int8_t a, int8_t b, int8_t sc static void ksw_gen_ts_mat(int m, int8_t *mat, int8_t a, int8_t b, int8_t transition, int8_t sc_ambi) { - assert(m==5); - ksw_gen_simple_mat(m,mat,a,b,sc_ambi); + assert(m == 5); + ksw_gen_simple_mat(m, mat, a, b, sc_ambi); + if (transition == 0 || transition == b) return; transition = transition > 0? -transition : transition; - mat[0*m+2]=transition; // A->G - mat[1*m+3]=transition; // C->T - mat[2*m+0]=transition; // G->A - mat[3*m+1]=transition; // T->C + mat[0 * m + 2] = transition; // A->G + mat[1 * m + 3] = transition; // C->T + mat[2 * m + 0] = transition; // G->A + mat[3 * m + 1] = transition; // T->C } static inline void mm_seq_rev(uint32_t len, uint8_t *seq) @@ -334,7 +335,8 @@ static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint for (i = 0; i < qlen; ++i) fputc("ACGTN"[qseq[i]], stderr); fputc('\n', stderr); } - if (opt->b != opt->transition) flag |= KSW_EZ_GENERIC_SC; + if (opt->transition != 0 && opt->b != opt->transition) + flag |= KSW_EZ_GENERIC_SC; if (opt->max_sw_mat > 0 && (int64_t)tlen * qlen > opt->max_sw_mat) { ksw_reset_extz(ez); ez->zdropped = 1; diff --git a/main.c b/main.c index d50b565..09c551e 100644 --- a/main.c +++ b/main.c @@ -368,12 +368,12 @@ int main(int argc, char *argv[]) fprintf(fp_help, " --version show version number\n"); fprintf(fp_help, " Preset:\n"); fprintf(fp_help, " -x STR preset (always applied before other options; see minimap2.1 for details) []\n"); - fprintf(fp_help, " - map-pb/map-ont/map-iclr-prerender/map-iclr - PacBio/Nanopore/ICLR vs reference mapping\n"); - fprintf(fp_help, " - map-hifi - PacBio HiFi reads vs reference mapping\n"); - fprintf(fp_help, " - ava-pb/ava-ont - PacBio/Nanopore read overlap\n"); + fprintf(fp_help, " - lr:hq - accurate long reads (error rate <1%%) against a reference genome\n"); + fprintf(fp_help, " - splice/splice:hq - spliced alignment for long reads/accurate long reads\n"); fprintf(fp_help, " - asm5/asm10/asm20 - asm-to-ref mapping, for ~0.1/1/5%% sequence divergence\n"); - fprintf(fp_help, " - splice/splice:hq - long-read/Pacbio-CCS spliced alignment\n"); - fprintf(fp_help, " - sr - genomic short-read mapping\n"); + fprintf(fp_help, " - sr - short reads against a reference\n"); + fprintf(fp_help, " - map-pb/map-hifi/map-ont/map-iclr - CLR/HiFi/Nanopore/ICLR vs reference mapping\n"); + fprintf(fp_help, " - ava-pb/ava-ont - PacBio CLR/Nanopore read overlap\n"); fprintf(fp_help, "\nSee `man ./minimap2.1' for detailed description of these and other advanced command-line options.\n"); return fp_help == stdout? 0 : 1; } diff --git a/minimap.h b/minimap.h index 59f6fa0..a8929c3 100644 --- a/minimap.h +++ b/minimap.h @@ -5,7 +5,7 @@ #include #include -#define MM_VERSION "2.26-r1175" +#define MM_VERSION "2.26-r1183-dirty" #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 diff --git a/minimap2.1 b/minimap2.1 index a88dc6c..c2e68d9 100644 --- a/minimap2.1 +++ b/minimap2.1 @@ -343,6 +343,10 @@ Matching score [2] .BI -B \ INT Mismatching penalty [4] .TP +.BI -b \ INT +Mismatching penalty for transitions [same as +.BR -B ]. +.TP .BI -O \ INT1[,INT2] Gap open penalty [4,24]. If .I INT2 @@ -356,10 +360,19 @@ costs .RI min{ O1 + k * E1 , O2 + k * E2 }. In the splice mode, the second gap penalties are not used. .TP +.BI -J \ INT +Splice model [1]. 0 for the original minimap2 splice model that always penalizes non-GT-AG splicing; +1 for the miniprot model that considers non-GT-AG. Option +.B -C +has no effect with the default +.BR -J1 . +.BR -J0 . +.TP .BI -C \ INT Cost for a non-canonical GT-AG splicing (effective with -.BR --splice ) -[0] +.B --splice +.BR -J0 ) +[0]. .TP .BI -z \ INT1[,INT2] Truncate an alignment if the running alignment score drops too quickly along @@ -569,15 +582,43 @@ are: Align noisy long reads of ~10% error rate to a reference genome. This is the default mode. .TP +.B lr:hq +Align accurate long reads (error rate <1%) to a reference genome +.RB ( -k19 +.B -w19 -U50,500 +.BR -g10k ). +This was recommended by ONT developers for recent Nanopore reads +produced with chemistry v14 that can reach ~99% in accuracy. +It was shown to work better for accurate Nanopore reads +than +.BR map-hifi . +.TP .B map-hifi Align PacBio high-fidelity (HiFi) reads to a reference genome -.RB ( -k19 -.B -w19 -U50,500 -g10k -A1 -B4 -O6,26 -E2,1 +.RB ( -xlr:hq +.B -A1 -B4 -O6,26 -E2,1 .BR -s200 ). +It differs from +.B lr:hq +only in scoring. It has not been tested whether +.B lr:hq +would work better for PacBio HiFi reads. .TP .B map-pb Align older PacBio continuous long (CLR) reads to a reference genome .RB ( -Hk19 ). +Note that this data type is effectively deprecated by HiFi. +Unless you work on very old data, you probably want to use +.B map-hifi +or +.BR lr:hq . +.TP +.B map-iclr +Align Illumina Complete Long Reads (ICLR) to a reference genome +.RB ( -k19 +.B -B6 -b4 +.BR -O10,50 ). +This was recommended by Illumina developers. .TP .B asm5 Long assembly to reference mapping @@ -585,21 +626,21 @@ Long assembly to reference mapping .B -w19 -U50,500 --rmq -r1k,100k -g10k -A1 -B19 -O39,81 -E3,1 -s200 -z200 .BR -N50 ). Typically, the alignment will not extend to regions with 5% or higher sequence -divergence. Only use this preset if the average divergence is far below 5%. +divergence. Use this preset if the average divergence is not much higher than 0.1%. .TP .B asm10 Long assembly to reference mapping .RB ( -k19 .B -w19 -U50,500 --rmq -r1k,100k -g10k -A1 -B9 -O16,41 -E2,1 -s200 -z200 .BR -N50 ). -Up to 10% sequence divergence. +Use this if the average divergence is around 1%. .TP .B asm20 Long assembly to reference mapping .RB ( -k19 .B -w10 -U50,500 --rmq -r1k,100k -g10k -A1 -B4 -O6,26 -E2,1 -s200 -z200 .BR -N50 ). -Up to 20% sequence divergence. +Use this if the average divergence is around several percent. .TP .B splice Long-read spliced alignment @@ -615,13 +656,13 @@ costs are different during chaining; 4) the computation of the tag ignores introns to demote hits to pseudogenes. .TP .B splice:hq -Long-read splice alignment for PacBio CCS reads +Spliced alignment for accurate long RNA-seq reads such as PacBio iso-seq .RB ( -xsplice .B -C5 -O6,24 .BR -B4 ). .TP .B sr -Short single-end reads without splicing +Short-read alignment without splicing .RB ( -k21 .B -w11 --sr --frag=yes -A2 -B8 -O12,32 -E2,1 -b0 -r100 -p.5 -N20 -f1000,5000 -n2 -m25 .B -s40 -g100 -2K50m --heap-sort=yes diff --git a/options.c b/options.c index cd63111..4ed6d34 100644 --- a/options.c +++ b/options.c @@ -45,7 +45,7 @@ void mm_mapopt_init(mm_mapopt_t *opt) opt->alt_drop = 0.15f; opt->a = 2, opt->b = 4, opt->q = 4, opt->e = 2, opt->q2 = 24, opt->e2 = 1; - opt->transition = opt->b; + opt->transition = 0; opt->sc_ambi = 1; opt->zdrop = 400, opt->zdrop_inv = 200; opt->end_bonus = -1; @@ -91,7 +91,7 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo) if (preset == 0) { mm_idxopt_init(io); mm_mapopt_init(mo); - } else if (strcmp(preset, "map-ont") == 0) { // this is the same as the default + } else if (strcmp(preset, "lr") == 0 || strcmp(preset, "map-ont") == 0) { // this is the same as the default } else if (strcmp(preset, "ava-ont") == 0) { io->flag = 0, io->k = 15, io->w = 5; mo->flag |= MM_F_ALL_CHAINS | MM_F_NO_DIAG | MM_F_NO_DUAL | MM_F_NO_LJOIN; @@ -106,13 +106,14 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo) mo->min_chain_score = 100, mo->pri_ratio = 0.0f, mo->max_chain_skip = 25; mo->bw_long = mo->bw; mo->occ_dist = 0; - } else if (strcmp(preset, "map-hifi") == 0 || strcmp(preset, "map-ccs") == 0) { + } else if (strcmp(preset, "lr:hq") == 0 || strcmp(preset, "map-hifi") == 0 || strcmp(preset, "map-ccs") == 0) { io->flag = 0, io->k = 19, io->w = 19; mo->max_gap = 10000; - mo->a = 1, mo->b = 4, mo->q = 6, mo->q2 = 26, mo->e = 2, mo->e2 = 1; - mo->occ_dist = 500; mo->min_mid_occ = 50, mo->max_mid_occ = 500; - mo->min_dp_max = 200; + if (strcmp(preset, "map-hifi") == 0 || strcmp(preset, "map-ccs") == 0) { + mo->a = 1, mo->b = 4, mo->q = 6, mo->q2 = 26, mo->e = 2, mo->e2 = 1; + mo->min_dp_max = 200; + } } else if (strcmp(preset, "map-iclr-prerender") == 0) { io->flag = 0, io->k = 15; mo->b = 6, mo->transition = 1; @@ -165,7 +166,7 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo) mo->junc_bonus = 9; mo->zdrop = 200, mo->zdrop_inv = 100; // because mo->a is halved if (strcmp(preset, "splice:hq") == 0) - mo->junc_bonus = 5, mo->b = 4, mo->q = 6, mo->q2 = 24; + mo->noncan = 5, mo->b = 4, mo->q = 6, mo->q2 = 24; } else return -1; return 0; }