fork/minimap2
1
0
Fork 0
mirror of https://github.com/lh3/minimap2.git synced 2026-09-19 02:28:01 +08:00
Code Issues Packages Projects Releases Wiki Activity

r1183: added lr:hq; fixed transition

Browse Source 
* Added the lr:hq preset suggested by Nanopore developers (#1127)
 * Fixed transition scoring. It did not work with presets.
 * Cleaned up preset documentation
This commit is contained in:
Heng Li
2024-03-10 13:47:34 -04:00
parent bc588c0eeb
commit 8140259974
5 changed files with 73 additions and 29 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
align.c
View File

@@ -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;

10
main.c
View File

@@ -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;
}

2
minimap.h
View File

@@ -5,7 +5,7 @@
#include <stdio.h>
#include <sys/types.h>
#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

59
minimap2.1
View File

@@ -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

15
options.c
View File

@@ -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;
}