diff --git a/align.c b/align.c index 9035b42..97c29f6 100644 --- a/align.c +++ b/align.c @@ -603,7 +603,7 @@ static inline void mm_get_junc(const mm_idx_t *mi, int32_t ctg, int32_t st, int3 static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, uint8_t *qseq0[2], mm_reg1_t *r, mm_reg1_t *r2, int n_a, mm128_t *a, ksw_extz_t *ez, int splice_flag) { - int is_sr = !!(opt->flag & MM_F_SR), is_splice = !!(opt->flag & MM_F_SPLICE); + int is_sr = !!(opt->flag & MM_F_SR), is_splice = !!(opt->flag & MM_F_SPLICE), is_sr_rna = (!!(opt->flag & MM_F_SR_RNA) && is_splice); int32_t rid = a[r->as].x<<1>>33, rev = a[r->as].x>>63, as1, cnt1; uint8_t *tseq, *qseq, *junc; int32_t i, l, bw, bw_long, dropped = 0, ksw_flag = 0, rs0, re0, qs0, qe0; @@ -773,14 +773,18 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int mm_idx_getseq(mi, rid, rs, re, tseq); } mm_get_junc(mi, rid, rs, re, !!(ksw_flag&KSW_EZ_SPLICE_REV), junc); - if (is_sr) { // perform ungapped alignment + if (is_sr || (is_sr_rna && qe - qs == re - rs)) { // perform ungapped alignment + int32_t max_gapped_score = (qe - qs - 2) * opt->a - 2 * (opt->q + opt->e); assert(qe - qs == re - rs); ksw_reset_extz(ez); for (j = 0, ez->score = 0; j < qe - qs; ++j) { - if (qseq[j] >= 4 || tseq[j] >= 4) ez->score += opt->e2; + if (qseq[j] >= 4 || tseq[j] >= 4) ez->score += opt->sc_ambi > 0? -opt->sc_ambi : opt->sc_ambi; else ez->score += qseq[j] == tseq[j]? opt->a : -opt->b; } - ez->cigar = ksw_push_cigar(km, &ez->n_cigar, &ez->m_cigar, ez->cigar, MM_CIGAR_MATCH, qe - qs); + if (ez->score > max_gapped_score) + ez->cigar = ksw_push_cigar(km, &ez->n_cigar, &ez->m_cigar, ez->cigar, MM_CIGAR_MATCH, qe - qs); + else + mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, junc, mat, bw1, -1, opt->zdrop, ksw_flag|KSW_EZ_APPROX_MAX, ez); } else { // perform normal gapped alignment mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, junc, mat, bw1, -1, opt->zdrop, ksw_flag|KSW_EZ_APPROX_MAX, ez); // first pass: with approximate Z-drop } diff --git a/main.c b/main.c index cca73c7..b61ea70 100644 --- a/main.c +++ b/main.c @@ -35,7 +35,7 @@ static ko_longopt_t long_options[] = { { "splice", ko_no_argument, 310 }, { "cost-non-gt-ag", ko_required_argument, 'C' }, { "no-long-join", ko_no_argument, 312 }, - { "sr", ko_no_argument, 313 }, + { "sr", ko_optional_argument, 313 }, { "frag", ko_required_argument, 314 }, { "secondary", ko_required_argument, 315 }, { "cs", ko_optional_argument, 316 }, @@ -218,7 +218,6 @@ int main(int argc, char *argv[]) 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 (deprecated) else if (c == 320) ipt.flag |= MM_I_NO_SEQ; // --idx-no-seq @@ -257,6 +256,17 @@ int main(int argc, char *argv[]) else if (c == 501) mm_dbg_flag |= MM_DBG_SEED_FREQ; // --dbg-seed-occ else if (c == 330) { fprintf(stderr, "[WARNING] \033[1;31m --lj-min-ratio has been deprecated.\033[0m\n"); + } else if (c == 313) { // --sr + if (o.arg == 0 || strcmp(o.arg, "dna") == 0) { + opt.flag |= MM_F_SR; + } else if (strcmp(o.arg, "rna") == 0) { + opt.flag |= MM_F_SR_RNA; + } else if (strcmp(o.arg, "no") == 0) { + opt.flag &= ~(uint64_t)(MM_F_SR|MM_F_SR_RNA); + } else if (mm_verbose >= 2) { + opt.flag |= MM_F_SR; + fprintf(stderr, "[WARNING]\033[1;31m --sr only takes 'dna' or 'rna'. Invalid values are assumed to be 'dna'.\033[0m\n"); + } } else if (c == 314) { // --frag yes_or_no(&opt, MM_F_FRAG_MODE, o.longidx, o.arg, 1); } else if (c == 315) { // --secondary diff --git a/minimap.h b/minimap.h index 6b66c8f..78ce306 100644 --- a/minimap.h +++ b/minimap.h @@ -5,7 +5,7 @@ #include #include -#define MM_VERSION "2.28-r1261-dirty" +#define MM_VERSION "2.28-r1263-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 @@ -46,6 +46,7 @@ #define MM_F_SECONDARY_SEQ (0x1000000000LL) //output SEQ field for seqondary alignments using hard clipping #define MM_F_OUT_DS (0x2000000000LL) #define MM_F_WEAK_PAIRING (0x4000000000LL) +#define MM_F_SR_RNA (0x8000000000LL) #define MM_I_HPC 0x1 #define MM_I_NO_SEQ 0x2 diff --git a/minimap2.1 b/minimap2.1 index 5daac3a..395b967 100644 --- a/minimap2.1 +++ b/minimap2.1 @@ -297,11 +297,13 @@ maximum alignment gap is mostly controlled by .B --splice Enable the splice alignment mode. .TP -.B --sr -Enable short-read alignment heuristics. In the short-read mode, minimap2 -applies a second round of chaining with a higher minimizer occurrence threshold -if no good chain is found. In addition, minimap2 attempts to patch gaps between -seeds with ungapped alignment. +.BR --sr [= no | dna | rna ] +Enable short-read alignment heuristics [no]. If this option is used with no argument, +.RB ` dna ' +is set. In the DNA short-read mode, minimap2 applies a second round of chaining +with a higher minimizer occurrence threshold if no good chain is found. In +addition, minimap2 attempts to patch gaps between seeds with ungapped +alignment. .TP .BI --split-prefix \ STR Prefix to create temporary files. Typically used for a multi-part index. @@ -520,20 +522,13 @@ Copy input FASTA/Q comments to output. .B -c Generate CIGAR. In PAF, the CIGAR is written to the `cg' custom tag. .TP -.BI --cs[= STR ] +.BR --cs [= short | long ] Output the .B cs tag. -.I STR -can be either -.I short -or -.IR long . -If no -.I STR -is given, -.I short -is assumed. [none] +If no argument is given, +.RB ` short ' +is set. [none] .TP .B --MD Output the MD tag (see the SAM spec). @@ -689,7 +684,7 @@ Spliced alignment for accurate long RNA-seq reads such as PacBio iso-seq .B splice:sr Spliced alignment for short RNA-seq reads .RB ( -xsplice:hq -.B --frag=yes -m25 -s40 -2K50m --heap-sort=yes --pairing=weak +.B --frag=yes -m25 -s40 -2K100m --heap-sort=yes --pairing=weak --sr=rna .BR --secondary=no ). .TP .B sr diff --git a/options.c b/options.c index 5e7e9b8..a099eae 100644 --- a/options.c +++ b/options.c @@ -179,13 +179,13 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo) if (strcmp(preset, "splice:hq") == 0) { mo->noncan = 5, mo->b = 4, mo->q = 6, mo->q2 = 24; } else if (strcmp(preset, "splice:sr") == 0) { - mo->flag |= MM_F_NO_PRINT_2ND | MM_F_2_IO_THREADS | MM_F_HEAP_SORT | MM_F_FRAG_MODE | MM_F_WEAK_PAIRING; + mo->flag |= MM_F_NO_PRINT_2ND | MM_F_2_IO_THREADS | MM_F_HEAP_SORT | MM_F_FRAG_MODE | MM_F_WEAK_PAIRING | MM_F_SR_RNA; mo->noncan = 5, mo->b = 4, mo->q = 6, mo->q2 = 24; mo->min_chain_score = 25; mo->min_dp_max = 40; mo->pe_ori = 0<<1|1; // FR mo->best_n = 10; - mo->mini_batch_size = 50000000; + mo->mini_batch_size = 100000000; } } else return -1; return 0;