Add faidx_seq_len64(), fai_adjust_region() interfaces and faidx tests

.gitattributes

@@ -20,3 +20,7 @@ README.md       export-ignore
 # Remove the text attribute from index_dos.sam, so that the line separators
 # for the test file don't get converted into Unix format.
 test/index_dos.sam -text
+
+# Remove the text attribute from various faidx test files
+test/faidx/faidx*.fa* -text
+test/faidx/fastqs*.fq* -text

.gitignore

@@ -44,6 +44,8 @@ shlib-exports-*.txt
 /bgzip
 /htsfile
 /tabix
+/test/faidx/*.tmp*
+/test/faidx/FAIL*
 /test/fieldarith
 /test/hfile
 /test/hts_endian
@@ -59,6 +61,7 @@ shlib-exports-*.txt
 /test/test-bcf_set_variant_type
 /test/test_bgzf
 /test/test_expr
+/test/test_faidx
 /test/test_index
 /test/test_introspection
 /test/test_kfunc

Makefile

@@ -80,6 +80,7 @@ BUILT_TEST_PROGRAMS = \
 	test/sam \
 	test/test_bgzf \
 	test/test_expr \
+	test/test_faidx \
 	test/test_kfunc \
 	test/test_kstring \
 	test/test_mod \
@@ -583,12 +584,13 @@ check test: all $(HTSCODECS_TEST_TARGETS)
 	fi
 	test/test_bgzf test/bgziptest.txt
 	test/test-parse-reg -t test/colons.bam
+	cd test/faidx && ./test-faidx.sh faidx.tst
 	cd test/sam_filter && ./filter.sh filter.tst
 	cd test/tabix && ./test-tabix.sh tabix.tst
 	cd test/mpileup && ./test-pileup.sh mpileup.tst
 	cd test/fastq && ./test-fastq.sh
 	cd test/base_mods && ./base-mods.sh base-mods.tst
-	REF_PATH=: test/sam test/ce.fa test/faidx.fa test/fastqs.fq
+	REF_PATH=: test/sam test/ce.fa test/faidx/faidx.fa test/faidx/fastqs.fq
 	test/test-regidx
 	cd test && REF_PATH=: ./test.pl $${TEST_OPTS:-}
 
@@ -622,6 +624,9 @@ test/test_bgzf: test/test_bgzf.o libhts.a
 test/test_expr: test/test_expr.o libhts.a
 	$(CC) $(LDFLAGS) -o $@ test/test_expr.o libhts.a -lz $(LIBS) -lpthread
 
+test/test_faidx: test/test_faidx.o libhts.a
+	$(CC) $(LDFLAGS) -o $@ test/test_faidx.o libhts.a -lz $(LIBS) -lpthread
+
 test/test_kfunc: test/test_kfunc.o libhts.a
 	$(CC) $(LDFLAGS) -o $@ test/test_kfunc.o libhts.a -lz $(LIBS) -lpthread
 
@@ -739,6 +744,7 @@ test/test-regidx.o: test/test-regidx.c config.h $(htslib_kstring_h) $(htslib_reg
 test/test_str2int.o: test/test_str2int.c config.h $(textutils_internal_h)
 test/test_time_funcs.o: test/test_time_funcs.c config.h $(hts_time_funcs_h)
 test/test_view.o: test/test_view.c config.h $(cram_h) $(htslib_sam_h) $(htslib_vcf_h) $(htslib_hts_log_h)
+test/test_faidx.o: test/test_faidx.c config.h $(htslib_faidx_h)
 test/test_index.o: test/test_index.c config.h $(htslib_sam_h) $(htslib_vcf_h)
 test/test-vcf-api.o: test/test-vcf-api.c config.h $(htslib_hts_h) $(htslib_vcf_h) $(htslib_kstring_h) $(htslib_kseq_h)
 test/test-vcf-sweep.o: test/test-vcf-sweep.c config.h $(htslib_vcf_sweep_h)
@@ -845,7 +851,9 @@ htslib-uninstalled.pc: htslib.pc.tmp
 
 
 testclean:
-	-rm -f test/*.tmp test/*.tmp.* test/longrefs/*.tmp.* test/tabix/*.tmp.* test/tabix/FAIL* header-exports.txt shlib-exports-$(SHLIB_FLAVOUR).txt
+	-rm -f test/*.tmp test/*.tmp.* test/faidx/*.tmp* test/faidx/FAIL* \
+               test/longrefs/*.tmp.* test/tabix/*.tmp.* test/tabix/FAIL* \
+               header-exports.txt shlib-exports-$(SHLIB_FLAVOUR).txt
 	-rm -rf htscodecs/tests/test.out
 
 # Only remove this in git checkouts

faidx.c

@@ -844,49 +844,83 @@ const char *faidx_iseq(const faidx_t *fai, int i)
     return fai->name[i];
 }
 
-int faidx_seq_len(const faidx_t *fai, const char *seq)
+hts_pos_t faidx_seq_len64(const faidx_t *fai, const char *seq)
 {
     khint_t k = kh_get(s, fai->hash, seq);
     if ( k == kh_end(fai->hash) ) return -1;
     return kh_val(fai->hash, k).len;
 }
 
-static int faidx_adjust_position(const faidx_t *fai, faidx1_t *val, const char *c_name, hts_pos_t *p_beg_i, hts_pos_t *p_end_i, hts_pos_t *len) {
+int faidx_seq_len(const faidx_t *fai, const char *seq)
+{
+    hts_pos_t len = faidx_seq_len64(fai, seq);
+    return len < INT_MAX ? len : INT_MAX;
+}
+
+static int faidx_adjust_position(const faidx_t *fai, int end_adjust,
+                                 faidx1_t *val_out, const char *c_name,
+                                 hts_pos_t *p_beg_i, hts_pos_t *p_end_i,
+                                 hts_pos_t *len) {
     khiter_t iter;
+    faidx1_t *val;
 
     // Adjust position
     iter = kh_get(s, fai->hash, c_name);
 
     if (iter == kh_end(fai->hash)) {
-        *len = -2;
+        if (len)
+            *len = -2;
         hts_log_error("The sequence \"%s\" was not found", c_name);
         return 1;
     }
 
-    *val = kh_value(fai->hash, iter);
+    val = &kh_value(fai->hash, iter);
+
+    if (val_out)
+        *val_out = *val;
 
     if(*p_end_i < *p_beg_i)
         *p_beg_i = *p_end_i;
 
     if(*p_beg_i < 0)
         *p_beg_i = 0;
     else if(val->len <= *p_beg_i)
-        *p_beg_i = val->len - 1;
+        *p_beg_i = val->len;
 
     if(*p_end_i < 0)
         *p_end_i = 0;
     else if(val->len <= *p_end_i)
-        *p_end_i = val->len - 1;
+        *p_end_i = val->len - end_adjust;
 
     return 0;
 }
 
+int fai_adjust_region(const faidx_t *fai, int tid,
+                      hts_pos_t *beg, hts_pos_t *end)
+{
+    hts_pos_t orig_beg, orig_end;
+
+    if (!fai || !beg || !end || tid < 0 || tid >= fai->n)
+        return -1;
+
+    orig_beg = *beg;
+    orig_end = *end;
+    if (faidx_adjust_position(fai, 0, NULL, fai->name[tid], beg, end, NULL) != 0) {
+        hts_log_error("Inconsistent faidx internal state - couldn't find \"%s\"",
+                      fai->name[tid]);
+        return -1;
+    }
+
+    return ((orig_beg != *beg ? 1 : 0) |
+            (orig_end != *end && orig_end < HTS_POS_MAX ? 2 : 0));
+}
+
 char *faidx_fetch_seq64(const faidx_t *fai, const char *c_name, hts_pos_t p_beg_i, hts_pos_t p_end_i, hts_pos_t *len)
 {
     faidx1_t val;
 
     // Adjust position
-    if (faidx_adjust_position(fai, &val, c_name, &p_beg_i, &p_end_i, len)) {
+    if (faidx_adjust_position(fai, 1, &val, c_name, &p_beg_i, &p_end_i, len)) {
         return NULL;
     }
 
@@ -907,7 +941,7 @@ char *faidx_fetch_qual64(const faidx_t *fai, const char *c_name, hts_pos_t p_beg
     faidx1_t val;
 
     // Adjust position
-    if (faidx_adjust_position(fai, &val, c_name, &p_beg_i, &p_end_i, len)) {
+    if (faidx_adjust_position(fai, 1, &val, c_name, &p_beg_i, &p_end_i, len)) {
         return NULL;
     }
 

htslib/faidx.h

@@ -292,7 +292,22 @@ int faidx_nseq(const faidx_t *fai);
 HTSLIB_EXPORT
 const char *faidx_iseq(const faidx_t *fai, int i);
 
-/// Return sequence length, -1 if not present
+/// Return sequence length
+/** @param  fai  Pointer to the faidx_t struct
+    @param  seq  Name of the sequence
+    @return Sequence length, or -1 if not present
+*/
+HTSLIB_EXPORT
+hts_pos_t faidx_seq_len64(const faidx_t *fai, const char *seq);
+
+/// Return sequence length
+/** @param  fai  Pointer to the faidx_t struct
+    @param  seq  Name of the sequence
+    @return Sequence length, or -1 if not present
+
+    @deprecated This funtion cannot handle very long sequences.
+                Use faidx_seq_len64() instead.
+*/
 HTSLIB_EXPORT
 int faidx_seq_len(const faidx_t *fai, const char *seq);
 
@@ -314,6 +329,27 @@ const char *fai_parse_region(const faidx_t *fai, const char *s,
                              int *tid, hts_pos_t *beg, hts_pos_t *end,
                              int flags);
 
+/// Adjust region to the actual sequence length
+/** @param  fai   Pointer to the faidx_t struct
+    @param  tid   Sequence index, as returned by fai_parse_region()
+    @param  beg[in,out]   The start of the region (0 based)
+    @param  end[in,out]   One past end of the region (0 based)
+    @return 1, 2, or 3 if @p beg, @p end, or both are adjusted,
+            0 if @p beg and @p end are unchanged
+            -1 on error
+
+    Looks up the length of @p tid, and then adjusts the values of @p beg
+    and @p end if they fall outside the boundaries of the sequence.
+
+    If @p beg > @p end, it will be set to @p end.
+
+    The return value indicates which, if any, of the inputs have been
+    adjusted.  -1 will be returned if @p tid is not a valid sequence index.
+*/
+HTSLIB_EXPORT
+int fai_adjust_region(const faidx_t *fai, int tid,
+                      hts_pos_t *beg, hts_pos_t *end);
+
 /// Sets the cache size of the underlying BGZF compressed file
 /** @param  fai         Pointer to the faidx_t struct
  *  @param  cache_size  Selected cache size in bytes

test/faidx/ce.1.expected.fa

@@ -0,0 +1,8 @@
+>CHROMOSOME_I:5001-5125 length: 125
+AACTGGTTCAAAAACAAAAATTTTTTAAACTGTACAAACTGTCCAAAAAT
+TCGTCGTAAATCGACACACCCTTCTCATTTTTTCAAAATTTTAATTGTTT
+TCGAATGTTTTTTTTGCAGAATAAT
+>CHROMOSOME_X:101-225 length: 125
+GCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGC
+CTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCTAAGCCT
+AAGCCTAAGCCTAAGCCTAAGCCTA

test/faidx/faidx.1.expected.fa

@@ -0,0 +1,6 @@
+>trailingblank2:28-33 length: 6
+GGGCCC
+>trailingblank3:4-5 length: 2
+TA
+>bar:4-5 length: 2
+TA

test/faidx/faidx.fa.expected.fai

@@ -0,0 +1,6 @@
+	4	2	4	5
+trailingblank1	33	23	12	13
+trailingblank2	72	111	24	25
+trailingblank3	5	234	4	6
+foo	8	252	6	7
+bar	8	280	8	9

test/faidx/faidx.tst

@@ -0,0 +1,74 @@
+#    Copyright (C) 2022 Genome Research Ltd.
+#
+#    Author: Robert Davies <[email protected]>
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+# DEALINGS IN THE SOFTWARE.
+
+# First field:
+#   INIT = initialisation, not counted in testing
+#   P = expected to pass (zero return; expected output matches, if present)
+#   N = expected to return non-zero
+#   F = expected to fail
+#
+# Second field (P/N/F only):
+#   Filename of expected output.  If '.', output is not checked
+#
+# Rest:
+#   Command to execute.  $bgzip and $test_faidx are replaced with the path to
+# bgzip and test_faidx.
+
+# Index fasta
+P . $test_faidx -i faidx.fa -f faidx.fa.tmp.fai -e faidx.fa.expected.fai
+
+# Test various functions on the fasta index
+P . $test_faidx -i faidx.fa -f faidx.fa.tmp.fai -t fai_line_length -e 24 trailingblank2
+P . $test_faidx -i faidx.fa -f faidx.fa.tmp.fai	-t faidx_has_seq -e 1 foo
+P . $test_faidx -i faidx.fa -f faidx.fa.tmp.fai -t faidx_has_seq -e 0 absent
+P . $test_faidx -i faidx.fa -f faidx.fa.tmp.fai -t faidx_iseq -e trailingblank3 3
+P . $test_faidx -i faidx.fa -f faidx.fa.tmp.fai -t faidx_seq_len -e 33 trailingblank1
+P . $test_faidx -i faidx.fa -f faidx.fa.tmp.fai -t faidx_seq_len64 -e 72 trailingblank2
+
+# Index fastq
+P . $test_faidx -i fastqs.fq -f fastqs.fq.tmp.fai -e fastqs.fq.expected.fai
+
+# Test various functions on the fastq index
+P . $test_faidx -i fastqs.fq -f fastqs.fq.tmp.fai -Q -t fai_line_length -e 63 FAKE0005_3
+P . $test_faidx -i fastqs.fq -f fastqs.fq.tmp.fai -Q -t fai_line_length -e 144 SRR014849.203935_3
+P . $test_faidx -i fastqs.fq -f fastqs.fq.tmp.fai -t faidx_has_seq -e 1 SRR014849.203935_3
+P . $test_faidx -i fastqs.fq -f fastqs.fq.tmp.fai -t faidx_has_seq -e 0 absent
+P . $test_faidx -i fastqs.fq -f fastqs.fq.tmp.fai -t faidx_iseq -e FAKE0005_1 0
+P . $test_faidx -i fastqs.fq -f fastqs.fq.tmp.fai -t faidx_seq_len -e 453 FSRRS4401CM938_1
+P . $test_faidx -i fastqs.fq -f fastqs.fq.tmp.fai -t faidx_seq_len64 -e 309 FSRRS4401AOV6A_4
+
+# Fasta retrieval tests
+P faidx.1.expected.fa $test_faidx -i faidx.fa -f faidx.fa.tmp.fai trailingblank2:28-33 trailingblank3:4-5 bar:4-5
+P faidx.1.expected.fa $test_faidx -i faidx.fa -f faidx.fa.tmp.fai -t fai_fetch trailingblank2:28-33 trailingblank3:4-5 bar:4-5
+P faidx.1.expected.fa $test_faidx -i faidx.fa -f faidx.fa.tmp.fai -t faidx_fetch_seq64 trailingblank2:28-33 trailingblank3:4-5 bar:4-5
+P faidx.1.expected.fa $test_faidx -i faidx.fa -f faidx.fa.tmp.fai -t fai_adjust_region trailingblank2:28-33 trailingblank3:4-5 bar:4-5
+
+# Fastq retrieval tests
+P fastqs.1.expected.fq $test_faidx -i fastqs.fq -f fastqs.fq.tmp.fai -Q FAKE0006_1:4-12 FSRRS4401BE7HA_1:81-120 FAKE0010_2 SRR014849.50939_3:71-90
+P fastqs.1.expected.fq $test_faidx -i fastqs.fq -f fastqs.fq.tmp.fai -Q -t fai_fetch FAKE0006_1:4-12 FSRRS4401BE7HA_1:81-120 FAKE0010_2 SRR014849.50939_3:71-90
+P fastqs.1.expected.fq $test_faidx -i fastqs.fq -f fastqs.fq.tmp.fai -Q -t faidx_fetch_seq64 FAKE0006_1:4-12 FSRRS4401BE7HA_1:81-120 FAKE0010_2 SRR014849.50939_3:71-90
+P fastqs.2.expected.fa $test_faidx -i fastqs.fq -f fastqs.fq.tmp.fai FAKE0006_1:4-12 FSRRS4401BE7HA_1:81-120 FAKE0010_2 SRR014849.50939_3:71-90
+
+# Indexing and retrieval on bgzip compressed fasta
+INIT $bgzip -c < ../ce.fa > ce.fa.tmp.gz
+P . $test_faidx -i ce.fa.tmp.gz -f ce.fa.tmp.gz.fai -g ce.fa.tmp.gz.gzi -e ../ce.fa.fai
+P ce.1.expected.fa $test_faidx -i ce.fa.tmp.gz -f ce.fa.tmp.gz.fai -g ce.fa.tmp.gz.gzi CHROMOSOME_I:5001-5125 CHROMOSOME_X:101-225

test/faidx/fastqs.1.expected.fq

@@ -0,0 +1,16 @@
+@FAKE0006_1:4-12 length: 9
+TGCATGCAT
++
+{zyxwvuts
+@FSRRS4401BE7HA_1:81-120 length: 40
+GCCCGTTTGTCGATATTTGtatttaaagtaatccgtcaca
++
+c^^^YRPOSNVU\YTMMMSMRKKKRUUNNNNS[`aa```\
+@FAKE0010_2 length: 30
+gatcrywsmkhbvdnGATCRYWSMKHBVDN
++
+I?5+I?5+I?5+I?5+I?5+I?5+I?5+I?
+@SRR014849.50939_3:71-90 length: 20
+CAATAAATCAATACATAAAA
++
+\aZ\d`OY[aY[[\[[e`WP

test/faidx/fastqs.2.expected.fa

@@ -0,0 +1,8 @@
+>FAKE0006_1:4-12 length: 9
+TGCATGCAT
+>FSRRS4401BE7HA_1:81-120 length: 40
+GCCCGTTTGTCGATATTTGtatttaaagtaatccgtcaca
+>FAKE0010_2 length: 30
+gatcrywsmkhbvdnGATCRYWSMKHBVDN
+>SRR014849.50939_3:71-90 length: 20
+CAATAAATCAATACATAAAA