diff --git a/FAQ.md b/FAQ.md new file mode 100644 index 0000000..94a1bbc --- /dev/null +++ b/FAQ.md @@ -0,0 +1,46 @@ +#### 1. Alignment different with option `-a` or `-c`? + +Without `-a`, `-c` or `--cs`, minimap2 only finds *approximate* mapping +locations without detailed base alignment. In particular, the start and end +positions of the alignment are impricise. With one of those options, minimap2 +will perform base alignment, which is generally more accurate but is much +slower. + +#### 2. How to map Illumina short reads to noisy long reads? + +No good solutions. The better approach is to assemble short reads into contigs +and then map noisy reads to contigs. + +#### 3. The output SAM doesn't have a header. + +By default, minimap2 indexes 4 billion reference bases (4Gb) in a batch and map +all reads against each reference batch. Given a reference longer than 4Gb, +minimap2 is unable to see all the sequences and thus can't produce a correct +SAM header. In this case, minimap2 doesn't output any SAM header. There are two +solutions to this issue. First, you may increase option `-I` to, for example, +`-I8g` to index more reference bases in a batch. This is preferred if your +machine has enough memory. Second, if your machines doesn't have enough memory +to hold the reference index, you can use the `--split-prefix` option in a +command line like: +```sh +minimap2 -ax map-ont --split-prefix=tmp ref.fa reads.fq +``` +This second approach uses less memory, but it is slower and requires temporary +disk space. + +#### 4. The output SAM is malformatted. + +This typically happens when you use nohup to wrap a minimap2 command line. +Nohup is discouraged as it breaks piping. If you have to use nohup, please +specify an output file with option `-o`. + +#### 5. How to output one alignment per read? + +You can use `--secondary=no` to suppress secondary alignments (aka multiple +mappings), but you can't suppress supplementary alignment (aka split or +chimeric alignment) this way. You can use samtools to filter out these +alignments: +```sh +minimap2 -ax map-out ref.fa reads.fq | samtools view -F0x900 +``` +However, this is discouraged as supplementary alignment is informative. diff --git a/Makefile b/Makefile index bcf32dd..93ca59b 100644 --- a/Makefile +++ b/Makefile @@ -34,6 +34,16 @@ else #if aarch64 is defined endif endif +ifneq ($(asan),) + CFLAGS+=-fsanitize=address + LIBS+=-fsanitize=address +endif + +ifneq ($(tsan),) + CFLAGS+=-fsanitize=thread + LIBS+=-fsanitize=thread +endif + .PHONY:all extra clean depend .SUFFIXES:.c .o diff --git a/README.md b/README.md index 5e39daf..addeb2c 100644 --- a/README.md +++ b/README.md @@ -315,9 +315,10 @@ highlighted in bold. The description may help to tune minimap2 parameters. ### Getting help Manpage [minimap2.1][manpage] provides detailed description of minimap2 -command line options and optional tags. If you encounter bugs or have further -questions or requests, you can raise an issue at the [issue page][issue]. -There is not a specific mailing list for the time being. +command line options and optional tags. The [FAQ](FAQ.md) page answers several +frequently asked questions. If you encounter bugs or have further questions or +requests, you can raise an issue at the [issue page][issue]. There is not a +specific mailing list for the time being. ### Citing minimap2 diff --git a/chain.c b/chain.c index be9d51b..95abf5e 100644 --- a/chain.c +++ b/chain.c @@ -155,8 +155,8 @@ mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int m memcpy(&a[k], &b[w[i].y>>32], n * sizeof(mm128_t)); k += n; } - memcpy(u, u2, n_u * 8); - memcpy(b, a, k * sizeof(mm128_t)); // write _a_ to _b_ and deallocate _a_ because _a_ is oversized, sometimes a lot + if (n_u) memcpy(u, u2, n_u * 8); + if (k) memcpy(b, a, k * sizeof(mm128_t)); // write _a_ to _b_ and deallocate _a_ because _a_ is oversized, sometimes a lot kfree(km, a); kfree(km, w); kfree(km, u2); return b; } diff --git a/example.c b/example.c index 80e493c..ca1fdbe 100644 --- a/example.c +++ b/example.c @@ -35,6 +35,8 @@ int main(int argc, char *argv[]) while ((mi = mm_idx_reader_read(r, n_threads)) != 0) { // traverse each part of the index mm_mapopt_update(&mopt, mi); // this sets the maximum minimizer occurrence; TODO: set a better default in mm_mapopt_init()! mm_tbuf_t *tbuf = mm_tbuf_init(); // thread buffer; for multi-threading, allocate one tbuf for each thread + gzrewind(f); + kseq_rewind(ks); while (kseq_read(ks) >= 0) { // each kseq_read() call reads one query sequence mm_reg1_t *reg; int j, i, n_reg; diff --git a/format.c b/format.c index 6d9e7ab..0dd943a 100644 --- a/format.c +++ b/format.c @@ -79,11 +79,11 @@ static char *mm_escape(char *s) return s; } -static void sam_write_rg_line(kstring_t *str, const char *s) +static int sam_write_rg_line(kstring_t *str, const char *s) { char *p, *q, *r, *rg_line = 0; memset(mm_rg_id, 0, 256); - if (s == 0) return; + if (s == 0) return 0; if (strstr(s, "@RG") != s) { if (mm_verbose >= 1) fprintf(stderr, "[ERROR] the read group line is not started with @RG\n"); goto err_set_rg; @@ -108,20 +108,23 @@ static void sam_write_rg_line(kstring_t *str, const char *s) for (q = p, r = mm_rg_id; *q && *q != '\t' && *q != '\n'; ++q) *r++ = *q; mm_sprintf_lite(str, "%s\n", rg_line); + return 0; err_set_rg: free(rg_line); + return -1; } -void mm_write_sam_hdr(const mm_idx_t *idx, const char *rg, const char *ver, int argc, char *argv[]) +int mm_write_sam_hdr(const mm_idx_t *idx, const char *rg, const char *ver, int argc, char *argv[]) { kstring_t str = {0,0,0}; + int ret = 0; if (idx) { uint32_t i; for (i = 0; i < idx->n_seq; ++i) mm_sprintf_lite(&str, "@SQ\tSN:%s\tLN:%d\n", idx->seq[i].name, idx->seq[i].len); } - if (rg) sam_write_rg_line(&str, rg); + if (rg) ret = sam_write_rg_line(&str, rg); mm_sprintf_lite(&str, "@PG\tID:minimap2\tPN:minimap2"); if (ver) mm_sprintf_lite(&str, "\tVN:%s", ver); if (argc > 1) { @@ -132,6 +135,7 @@ void mm_write_sam_hdr(const mm_idx_t *idx, const char *rg, const char *ver, int } mm_err_puts(str.s); free(str.s); + return ret; } static void write_cs_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq, const mm_reg1_t *r, char *tmp, int no_iden, int write_tag) diff --git a/index.c b/index.c index 164e8e5..6a49bb1 100644 --- a/index.c +++ b/index.c @@ -117,8 +117,8 @@ void mm_idx_stat(const mm_idx_t *mi) if (kh_key(h, k)&1) ++n1; } } - fprintf(stderr, "[M::%s::%.3f*%.2f] distinct minimizers: %d (%.2f%% are singletons); average occurrences: %.3lf; average spacing: %.3lf\n", - __func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), n, 100.0*n1/n, (double)sum / n, (double)len / sum); + fprintf(stderr, "[M::%s::%.3f*%.2f] distinct minimizers: %d (%.2f%% are singletons); average occurrences: %.3lf; average spacing: %.3lf; total length: %ld\n", + __func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), n, 100.0*n1/n, (double)sum / n, (double)len / sum, (long)len); } int mm_idx_index_name(mm_idx_t *mi) diff --git a/kalloc.c b/kalloc.c index b36c333..8499552 100644 --- a/kalloc.c +++ b/kalloc.c @@ -18,15 +18,14 @@ * | | | | * p=p->ptr->ptr->ptr->ptr p->ptr p->ptr->ptr p->ptr->ptr->ptr */ - -#define MIN_CORE_SIZE 0x80000 - typedef struct header_t { size_t size; struct header_t *ptr; } header_t; typedef struct { + void *par; + size_t min_core_size; header_t base, *loop_head, *core_head; /* base is a zero-sized block always kept in the loop */ } kmem_t; @@ -36,31 +35,39 @@ static void panic(const char *s) abort(); } -void *km_init(void) +void *km_init2(void *km_par, size_t min_core_size) { - return calloc(1, sizeof(kmem_t)); + kmem_t *km; + km = (kmem_t*)kcalloc(km_par, 1, sizeof(kmem_t)); + km->par = km_par; + km->min_core_size = min_core_size > 0? min_core_size : 0x80000; + return (void*)km; } +void *km_init(void) { return km_init2(0, 0); } + void km_destroy(void *_km) { kmem_t *km = (kmem_t*)_km; + void *km_par; header_t *p, *q; if (km == NULL) return; + km_par = km->par; for (p = km->core_head; p != NULL;) { q = p->ptr; - free(p); + kfree(km_par, p); p = q; } - free(km); + kfree(km_par, km); } static header_t *morecore(kmem_t *km, size_t nu) { header_t *q; size_t bytes, *p; - nu = (nu + 1 + (MIN_CORE_SIZE - 1)) / MIN_CORE_SIZE * MIN_CORE_SIZE; /* the first +1 for core header */ + nu = (nu + 1 + (km->min_core_size - 1)) / km->min_core_size * km->min_core_size; /* the first +1 for core header */ bytes = nu * sizeof(header_t); - q = (header_t*)malloc(bytes); + q = (header_t*)kmalloc(km->par, bytes); if (!q) panic("[morecore] insufficient memory"); q->ptr = km->core_head, q->size = nu, km->core_head = q; p = (size_t*)(q + 1); @@ -125,7 +132,7 @@ void *kmalloc(void *_km, size_t n_bytes) if (n_bytes == 0) return 0; if (km == NULL) return malloc(n_bytes); - n_units = (n_bytes + sizeof(size_t) + sizeof(header_t) - 1) / sizeof(header_t) + 1; + n_units = (n_bytes + sizeof(size_t) + sizeof(header_t) - 1) / sizeof(header_t); /* header+n_bytes requires at least this number of units */ if (!(q = km->loop_head)) /* the first time when kmalloc() is called, intialize it */ q = km->loop_head = km->base.ptr = &km->base; @@ -160,18 +167,18 @@ void *kcalloc(void *_km, size_t count, size_t size) void *krealloc(void *_km, void *ap, size_t n_bytes) // TODO: this can be made more efficient in principle { kmem_t *km = (kmem_t*)_km; - size_t n_units, *p, *q; + size_t cap, *p, *q; if (n_bytes == 0) { kfree(km, ap); return 0; } if (km == NULL) return realloc(ap, n_bytes); if (ap == NULL) return kmalloc(km, n_bytes); - n_units = (n_bytes + sizeof(size_t) + sizeof(header_t) - 1) / sizeof(header_t); p = (size_t*)ap - 1; - if (*p >= n_units) return ap; /* TODO: this prevents shrinking */ + cap = (*p) * sizeof(header_t) - sizeof(size_t); + if (cap >= n_bytes) return ap; /* TODO: this prevents shrinking */ q = (size_t*)kmalloc(km, n_bytes); - memcpy(q, ap, (*p - 1) * sizeof(header_t)); + memcpy(q, ap, cap); kfree(km, ap); return q; } diff --git a/kalloc.h b/kalloc.h index e891892..6d72e4e 100644 --- a/kalloc.h +++ b/kalloc.h @@ -17,6 +17,7 @@ void *kcalloc(void *km, size_t count, size_t size); void kfree(void *km, void *ptr); void *km_init(void); +void *km_init2(void *km_par, size_t min_core_size); void km_destroy(void *km); void km_stat(const void *_km, km_stat_t *s); @@ -24,4 +25,13 @@ void km_stat(const void *_km, km_stat_t *s); } #endif +#define KMALLOC(km, ptr, len) ((ptr) = (__typeof__(ptr))kmalloc((km), (len) * sizeof(*(ptr)))) +#define KCALLOC(km, ptr, len) ((ptr) = (__typeof__(ptr))kcalloc((km), (len), sizeof(*(ptr)))) +#define KREALLOC(km, ptr, len) ((ptr) = (__typeof__(ptr))krealloc((km), (ptr), (len) * sizeof(*(ptr)))) + +#define KEXPAND(km, a, m) do { \ + (m) = (m) >= 4? (m) + ((m)>>1) : 16; \ + KREALLOC((km), (a), (m)); \ + } while (0) + #endif diff --git a/main.c b/main.c index 29257e4..63730c1 100644 --- a/main.c +++ b/main.c @@ -1,12 +1,13 @@ #include #include #include +#include #include "bseq.h" #include "minimap.h" #include "mmpriv.h" #include "ketopt.h" -#define MM_VERSION "2.17-r943-dirty" +#define MM_VERSION "2.17-r963-dirty" #ifdef __linux__ #include @@ -172,7 +173,7 @@ int main(int argc, char *argv[]) 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\n", o.arg); + fprintf(stderr, "[ERROR]\033[1;31m failed to write the output to file '%s'\033[0m: %s\n", o.arg, strerror(errno)); exit(1); } } @@ -208,6 +209,7 @@ int main(int argc, char *argv[]) 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 == 314) { // --frag yes_or_no(&opt, MM_F_FRAG_MODE, o.longidx, o.arg, 1); @@ -322,11 +324,11 @@ 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: PacBio/Nanopore vs reference mapping\n"); - fprintf(fp_help, " - ava-pb/ava-ont: PacBio/Nanopore read overlap\n"); - fprintf(fp_help, " - asm5/asm10/asm20: asm-to-ref mapping, for ~0.1/1/5%% sequence divergence\n"); - fprintf(fp_help, " - splice: long-read spliced alignment\n"); - fprintf(fp_help, " - sr: genomic short-read mapping\n"); + fprintf(fp_help, " - map-pb/map-ont - PacBio/Nanopore vs reference mapping\n"); + fprintf(fp_help, " - ava-pb/ava-ont - PacBio/Nanopore read overlap\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, "\nSee `man ./minimap2.1' for detailed description of these and other advanced command-line options.\n"); return fp_help == stdout? 0 : 1; } @@ -337,7 +339,7 @@ int main(int argc, char *argv[]) } idx_rdr = mm_idx_reader_open(argv[o.ind], &ipt, fnw); if (idx_rdr == 0) { - fprintf(stderr, "[ERROR] failed to open file '%s'\n", argv[o.ind]); + 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) { @@ -355,13 +357,19 @@ int main(int argc, char *argv[]) return 1; } if ((opt.flag & MM_F_OUT_SAM) && idx_rdr->n_parts == 1) { + int ret; if (mm_idx_reader_eof(idx_rdr)) { - mm_write_sam_hdr(mi, rg, MM_VERSION, argc, argv); + ret = mm_write_sam_hdr(mi, rg, MM_VERSION, argc, argv); } else { - mm_write_sam_hdr(0, rg, MM_VERSION, argc, argv); + 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) { + mm_idx_destroy(mi); + mm_idx_reader_close(idx_rdr); + return 1; + } } if (mm_verbose >= 3) fprintf(stderr, "[M::%s::%.3f*%.2f] loaded/built the index for %d target sequence(s)\n", @@ -384,7 +392,7 @@ int main(int argc, char *argv[]) mm_split_merge(argc - (o.ind + 1), (const char**)&argv[o.ind + 1], &opt, n_parts); if (fflush(stdout) == EOF) { - fprintf(stderr, "[ERROR] failed to write the results\n"); + perror("[ERROR] failed to write the results"); exit(EXIT_FAILURE); } diff --git a/map.c b/map.c index d33405c..58e74fc 100644 --- a/map.c +++ b/map.c @@ -1,6 +1,7 @@ #include #include #include +#include #include "kthread.h" #include "kvec.h" #include "kalloc.h" @@ -622,7 +623,7 @@ static mm_bseq_file_t **open_bseqs(int n, const char **fn) 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'\n", fn[i]); + fprintf(stderr, "ERROR: failed to open file '%s': %s\n", fn[i], strerror(errno)); for (j = 0; j < i; ++j) mm_bseq_close(fp[j]); free(fp); diff --git a/minimap2.1 b/minimap2.1 index 21b84bd..ff0b922 100644 --- a/minimap2.1 +++ b/minimap2.1 @@ -638,6 +638,7 @@ s2 i Chaining score of the best secondary chain NM i Total number of mismatches and gaps in the alignment MD Z To generate the ref sequence in the alignment AS i DP alignment score +SA Z List of other supplementary alignments ms i DP score of the max scoring segment in the alignment nn i Number of ambiguous bases in the alignment ts A Transcript strand (splice mode only) diff --git a/misc.c b/misc.c index f6b8c5f..ddfa9de 100644 --- a/misc.c +++ b/misc.c @@ -125,7 +125,7 @@ void mm_err_puts(const char *str) int ret; ret = puts(str); if (ret == EOF) { - fprintf(stderr, "[ERROR] failed to write the results\n"); + perror("[ERROR] failed to write the results"); exit(EXIT_FAILURE); } } @@ -135,7 +135,7 @@ void mm_err_fwrite(const void *p, size_t size, size_t nitems, FILE *fp) int ret; ret = fwrite(p, size, nitems, fp); if (ret == EOF) { - fprintf(stderr, "[ERROR] failed to write data\n"); + perror("[ERROR] failed to write data"); exit(EXIT_FAILURE); } } @@ -145,7 +145,7 @@ void mm_err_fread(void *p, size_t size, size_t nitems, FILE *fp) int ret; ret = fread(p, size, nitems, fp); if (ret == EOF) { - fprintf(stderr, "[ERROR] failed to read data\n"); + perror("[ERROR] failed to read data"); exit(EXIT_FAILURE); } } diff --git a/misc/paftools.js b/misc/paftools.js index a4a8c0e..d7bea8c 100755 --- a/misc/paftools.js +++ b/misc/paftools.js @@ -1,6 +1,6 @@ #!/usr/bin/env k8 -var paftools_version = '2.17-r941'; +var paftools_version = '2.17-r949-dirty'; /***************************** ***** Library functions ***** @@ -1509,6 +1509,7 @@ function paf_gff2bed(args) var colors = { 'protein_coding':'0,128,255', + 'mRNA':'0,128,255', 'lincRNA':'0,192,0', 'snRNA':'0,192,0', 'miRNA':'0,192,0', @@ -1541,8 +1542,8 @@ function paf_gff2bed(args) print(a[0][0], st, en, name, 1000, a[0][3], cds_st, cds_en, color, a.length, sizes.join(",") + ",", starts.join(",") + ","); } - var re_gtf = /(transcript_id|transcript_type|transcript_biotype|gene_name|transcript_name) "([^"]+)";/g; - var re_gff3 = /(transcript_id|transcript_type|transcript_biotype|gene_name|transcript_name)=([^;]+)/g; + var re_gtf = /\b(transcript_id|transcript_type|transcript_biotype|gene_name|gene_id|gbkey|transcript_name) "([^"]+)";/g; + var re_gff3 = /\b(transcript_id|transcript_type|transcript_biotype|gene_name|gene_id|gbkey|transcript_name)=([^;]+)/g; var buf = new Bytes(); var file = args[getopt.ind] == '-'? new File() : new File(args[getopt.ind]); @@ -1559,19 +1560,19 @@ function paf_gff2bed(args) if (t[2] != "CDS" && t[2] != "exon") continue; t[3] = parseInt(t[3]) - 1; t[4] = parseInt(t[4]); - var id = null, type = "", gname = "N/A", biotype = "", m, tname = "N/A"; + var id = null, type = "", name = "N/A", biotype = "", m, tname = "N/A"; while ((m = re_gtf.exec(t[8])) != null) { if (m[1] == "transcript_id") id = m[2]; else if (m[1] == "transcript_type") type = m[2]; - else if (m[1] == "transcript_biotype") biotype = m[2]; - else if (m[1] == "gene_name") name = m[2]; + else if (m[1] == "transcript_biotype" || m[1] == "gbkey") biotype = m[2]; + else if (m[1] == "gene_name" || m[1] == "gene_id") name = m[2]; else if (m[1] == "transcript_name") tname = m[2]; } while ((m = re_gff3.exec(t[8])) != null) { if (m[1] == "transcript_id") id = m[2]; else if (m[1] == "transcript_type") type = m[2]; - else if (m[1] == "transcript_biotype") biotype = m[2]; - else if (m[1] == "gene_name") name = m[2]; + else if (m[1] == "transcript_biotype" || m[1] == "gbkey") biotype = m[2]; + else if (m[1] == "gene_name" || m[1] == "gene_id") name = m[2]; else if (m[1] == "transcript_name") tname = m[2]; } if (type == "" && biotype != "") type = biotype; diff --git a/mmpriv.h b/mmpriv.h index c9b91cd..c1164e2 100644 --- a/mmpriv.h +++ b/mmpriv.h @@ -59,7 +59,7 @@ uint32_t ks_ksmall_uint32_t(size_t n, uint32_t arr[], size_t kk); void mm_sketch(void *km, const char *str, int len, int w, int k, uint32_t rid, int is_hpc, mm128_v *p); -void mm_write_sam_hdr(const mm_idx_t *mi, const char *rg, const char *ver, int argc, char *argv[]); +int mm_write_sam_hdr(const mm_idx_t *mi, const char *rg, const char *ver, int argc, char *argv[]); void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int opt_flag); void mm_write_paf3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int opt_flag, int rep_len); void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, int n_regs, const mm_reg1_t *regs); diff --git a/python/mappy.pyx b/python/mappy.pyx index f036450..3273514 100644 --- a/python/mappy.pyx +++ b/python/mappy.pyx @@ -113,6 +113,7 @@ cdef class Aligner: cdef cmappy.mm_mapopt_t map_opt def __cinit__(self, fn_idx_in=None, preset=None, k=None, w=None, min_cnt=None, min_chain_score=None, min_dp_score=None, bw=None, best_n=None, n_threads=3, fn_idx_out=None, max_frag_len=None, extra_flags=None, seq=None, scoring=None): + self._idx = NULL cmappy.mm_set_opt(NULL, &self.idx_opt, &self.map_opt) # set the default options if preset is not None: cmappy.mm_set_opt(str.encode(preset), &self.idx_opt, &self.map_opt) # apply preset @@ -170,6 +171,7 @@ cdef class Aligner: cdef void *km cdef cmappy.mm_mapopt_t map_opt + if self._idx == NULL: return map_opt = self.map_opt if max_frag_len is not None: map_opt.max_frag_len = max_frag_len if extra_flags is not None: map_opt.flag |= extra_flags @@ -186,27 +188,36 @@ cdef class Aligner: _seq2 = seq2 if isinstance(seq2, bytes) else seq2.encode() regs = cmappy.mm_map_aux(self._idx, _seq, _seq2, &n_regs, b._b, &map_opt) - for i in range(n_regs): - cmappy.mm_reg2hitpy(self._idx, ®s[i], &h) - cigar, _cs, _MD = [], '', '' - for k in range(h.n_cigar32): # convert the 32-bit CIGAR encoding to Python array - c = h.cigar32[k] - cigar.append([c>>4, c&0xf]) - if cs or MD: # generate the cs and/or the MD tag, if requested - if cs: - l_cs_str = cmappy.mm_gen_cs(km, &cs_str, &m_cs_str, self._idx, ®s[i], _seq, 1) - _cs = cs_str[:l_cs_str] if isinstance(cs_str, str) else cs_str[:l_cs_str].decode() - if MD: - l_cs_str = cmappy.mm_gen_MD(km, &cs_str, &m_cs_str, self._idx, ®s[i], _seq) - _MD = cs_str[:l_cs_str] if isinstance(cs_str, str) else cs_str[:l_cs_str].decode() - yield Alignment(h.ctg, h.ctg_len, h.ctg_start, h.ctg_end, h.strand, h.qry_start, h.qry_end, h.mapq, cigar, h.is_primary, h.mlen, h.blen, h.NM, h.trans_strand, h.seg_id, _cs, _MD) - cmappy.mm_free_reg1(®s[i]) - free(regs) - free(cs_str) + try: + i = 0 + while i < n_regs: + cmappy.mm_reg2hitpy(self._idx, ®s[i], &h) + cigar, _cs, _MD = [], '', '' + for k in range(h.n_cigar32): # convert the 32-bit CIGAR encoding to Python array + c = h.cigar32[k] + cigar.append([c>>4, c&0xf]) + if cs or MD: # generate the cs and/or the MD tag, if requested + if cs: + l_cs_str = cmappy.mm_gen_cs(km, &cs_str, &m_cs_str, self._idx, ®s[i], _seq, 1) + _cs = cs_str[:l_cs_str] if isinstance(cs_str, str) else cs_str[:l_cs_str].decode() + if MD: + l_cs_str = cmappy.mm_gen_MD(km, &cs_str, &m_cs_str, self._idx, ®s[i], _seq) + _MD = cs_str[:l_cs_str] if isinstance(cs_str, str) else cs_str[:l_cs_str].decode() + yield Alignment(h.ctg, h.ctg_len, h.ctg_start, h.ctg_end, h.strand, h.qry_start, h.qry_end, h.mapq, cigar, h.is_primary, h.mlen, h.blen, h.NM, h.trans_strand, h.seg_id, _cs, _MD) + cmappy.mm_free_reg1(®s[i]) + i += 1 + finally: + while i < n_regs: + cmappy.mm_free_reg1(®s[i]) + i += 1 + free(regs) + free(cs_str) def seq(self, str name, int start=0, int end=0x7fffffff): cdef int l - cdef char *s = cmappy.mappy_fetch_seq(self._idx, name.encode(), start, end, &l) + cdef char *s + if self._idx == NULL: return + s = cmappy.mappy_fetch_seq(self._idx, name.encode(), start, end, &l) if l == 0: return None r = s[:l] if isinstance(s, str) else s[:l].decode() free(s) @@ -224,6 +235,7 @@ cdef class Aligner: @property def seq_names(self): cdef char *p + if self._idx == NULL: return sn = [] for i in range(self._idx.n_seq): p = self._idx.seq[i].name diff --git a/splitidx.c b/splitidx.c index 69f4864..507e498 100644 --- a/splitidx.c +++ b/splitidx.c @@ -2,6 +2,7 @@ #include #include #include +#include #include "mmpriv.h" FILE *mm_split_init(const char *prefix, const mm_idx_t *mi) @@ -13,7 +14,7 @@ FILE *mm_split_init(const char *prefix, const mm_idx_t *mi) sprintf(fn, "%s.%.4d.tmp", prefix, mi->index); if ((fp = fopen(fn, "wb")) == NULL) { if (mm_verbose >= 1) - fprintf(stderr, "[ERROR]\033[1;31m failed to write to temporary file '%s'\033[0m\n", fn); + fprintf(stderr, "[ERROR]\033[1;31m failed to write to temporary file '%s'\033[0m: %s\n", fn, strerror(errno)); exit(1); } mm_err_fwrite(&k, 4, 1, fp); @@ -41,7 +42,7 @@ mm_idx_t *mm_split_merge_prep(const char *prefix, int n_splits, FILE **fp, uint3 sprintf(fn, "%s.%.4d.tmp", prefix, i); if ((fp[i] = fopen(fn, "rb")) == 0) { if (mm_verbose >= 1) - fprintf(stderr, "ERROR: failed to open temporary file '%s'\n", fn); + fprintf(stderr, "ERROR: failed to open temporary file '%s': %s\n", fn, strerror(errno)); for (j = 0; j < i; ++j) fclose(fp[j]); free(fn);