Answers.Body: <p>The helpdesk for my institution's supercomputers showed me the following flags:</p>
  <pre><code>g++ seqan_test.cpp -DSEQAN_HAS_ZLIB -lz -lpthread
  </code></pre>
  <p>This fixed the error. I'm now running into different errors and am suspicious that I am getting them because I am using SeqAn2 rather than SeqAn3. I'm going to try HTSlib.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-06-29T18:24:05.813Z
Answers.LastActivityDate: 2021-06-29T18:24:05.813Z
Answers.UserDisplayName: annabelperry
Answers.Body: <p><a href="https://github.com/samtools/htslib" rel="nofollow noreferrer">Htslib</a> is the official library to access alignment files. It supports more formats. The typical way to read through a SAM/BAM/CRAM:</p>
  <pre class="lang-cpp prettyprint-override"><code>#include &quot;htslib/sam.h&quot;
  int main(int argc, char *argv[])
  {
    samFile *fp = sam_open(argv[1], &quot;r&quot;);
    bam_hdr_t *h = sam_hdr_read(fp);
    bam1_t *b = bam_init1();
    while (sam_read1(fp, h, b) &gt;= 0) {
      if (b-&gt;core.tid &lt; 0) continue;
      int end_pos = bam_endpos(b);
      printf(&quot;%s\t%d\t%d\n&quot;, h-&gt;target_name[b-&gt;core.tid], b-&gt;core.pos, end_pos);
    }
    bam_destroy1(b);
    bam_hdr_destroy(h);
    sam_close(fp);
  }
  </code></pre>
  <p>PS: the above started as a comment for &quot;relevant but minor or transient information to&quot; the original question. I didn't consider htslib a direct answer to the question, but because this had been changed to an answer, I added some details anyway. By the way, seqan2 has largely been discontinued. Use seqan3.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-06-29T16:22:51.407Z
Answers.Text: Thank you for this! I tested this script and it worked. I will try to edit it to suit my needs. Also, I am using SeqAn2 because SeqAn3 is not available on my institution's supercomputers. If HTSlib doesn't pan out or if no one can figure this error out, I will email them to request that it be installed.
Answers.UserDisplayName: annabelperry
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-06-28T09:16:04.920Z
Answers.LastActivityDate: 2021-06-30T10:31:31.487Z
Answers.LastEditDate: 2021-06-30T10:31:31.487Z
Answers.UserDisplayName: gringer
Answers.Body: <p>I'm not quite sure why assigning a string to a variable would lead to a complaint about file extensions. From the <a href="https://seqan.readthedocs.io/en/master/Tutorial/InputOutput/SamAndBamIO.html#tutorial-io-sam-bam-io" rel="nofollow noreferrer">SeqAn BAM tutorial</a>, it looks like you need to use the <em>BamFileIn</em> object for reading BAM/SAM files:</p>
  <pre><code>#include &lt;seqan/bam_io.h&gt;
  using namespace seqan;
  
  int main(){
      CharString bamFileName = &quot;copyPolly_SRR6511930.bam&quot;;
      // Open input file, BamFileIn can read SAM and BAM files.
      BamFileIn bamFileIn;
      if (!open(bamFileIn, toCString(bamFileName))) {
          std::cerr &lt;&lt; &quot;ERROR: Could not open &quot; &lt;&lt; bamFileName &lt;&lt; std::endl;
          return 1;
      }
      return 0;
  }
  </code></pre>
  <p>Edit: after seeing the updated code... I've got nothing, sorry. It all seems like reasonable code that shouldn't be generating errors.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-06-29T16:05:15.550Z
Answers.Text: Sorry about that, I did not include all of my code. I have updated my post with the rest of the code which led to my error.
Answers.UserDisplayName: annabelperry
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-06-28T23:12:37.507Z
Answers.LastActivityDate: 2021-06-30T10:32:27.130Z
Answers.LastEditDate: 2021-06-30T10:32:27.130Z
Answers.UserDisplayName: gringer
Answers.Body: <p>I think <code>SEQAN_HAS_ZLIB</code> flag not works correctly as it supposed to be.
  So you can check the seqan library to see if <code>SEQAN_HAS_ZLIB</code> works correctly, and seek the reason of extension error.</p>
  <p>Seqan library is a mature library which using template to define general-purpose classes and structs.</p>
  <p>1.The reason for extension error.</p>
  <p>Step into the seqan library:</p>
  <p>When you call BamfileIn: <code>BamFileIn bamFileIn(toCString(bamFileName));</code></p>
  <p>An object of BamFileIn struct is created, and <code>bamFileName</code> file is opened by the struct construction function.</p>
  <p>And when this object is created, some extra information is wrapped with the object to check the type of file.</p>
  <p>This information include suffix categories. In bam_file.h header file.</p>
  <pre><code>template &lt;typename TDirection, typename TSpec&gt;
  struct FileFormat&lt;FormattedFile&lt;Bam, TDirection, TSpec&gt; &gt;
  {
  #if SEQAN_HAS_ZLIB
      typedef TagSelector&lt;
                  TagList&lt;Bam,
                  TagList&lt;Sam
                  &gt; &gt;
              &gt; Type;
  #else
      typedef Sam Type;
  #endif
  };
  </code></pre>
  <p>As you can see, if <code>SEQAN_HAS_ZLIB</code> is setted, &quot;Bam&quot; file will be added to TagSelector which will be wrapped into the BamFileIn struct and redefined by <code>typedef</code> as <code>format</code>.</p>
  <p>Then _open() function is called to open Bam file, and extension is checked in guess_format.h</p>
  <pre><code>template &lt;typename TFilename, typename TTagList&gt;
  inline bool guessFormatFromFilename(TFilename const &amp;fname, TagSelector&lt;TTagList&gt; &amp;format)
  {
      typedef typename TTagList::Type TFormat;
  
      if (value(format) == -1 || isEqual(format, TFormat()))
      {
          // if tagId is set to -1 (auto-detect) or the current format (TFormatTag) then test for TFormatTag format
          if (guessFormatFromFilename(fname, TFormat()))
          {
              assign(format, TFormat());
              return true;
          }
      }
      return guessFormatFromFilename(fname, static_cast&lt;typename TagSelector&lt;TTagList&gt;::Base &amp;&gt;(format));
  }
  
  template &lt;typename TFilename, typename TFormat_&gt;
  inline bool guessFormatFromFilename(TFilename const &amp; fileName, Tag&lt;TFormat_&gt; /*formatTag*/)
  {
      typedef typename Value&lt;TFilename&gt;::Type                                     TValue;
      typedef ModifiedString&lt;TFilename const, ModView&lt;FunctorLowcase&lt;TValue&gt; &gt; &gt;    TLowcase;
      typedef Tag&lt;TFormat_&gt;                                                       TFormat;
  
      TLowcase lowcaseFileName(fileName);
      for (unsigned i = 0; i &lt; sizeof(FileExtensions&lt;TFormat&gt;::VALUE) / sizeof(char*); ++i)
          if (endsWith(lowcaseFileName, lowerString(FileExtensions&lt;TFormat&gt;::VALUE[i])))
              return true;
  
      return false;
  }
  </code></pre>
  <p>The above code is two instance of  <code>guessFormatFromFilename</code> function.
  It's pretty complex here, but if <code>SEQAN_HAS_ZLIB</code> is set correctly, <code>TagSelector</code> will include bam and sam, and the first <code>guessFormatFromFilename</code> above will return true in if condition, and the error will not happend. But if <code>SEQAN_HAS_ZLIB</code> is not working, second formt <code>guessFormatFromFilename</code> function will be called and only &quot;sam&quot; suffix is exist in <code>TFormat</code>, and it will be exhibited by <code>FileExtensions&lt;TFormat&gt;::VALUE[i]</code>. You can using std::cerr to print the suffix on your screen to see if there is only &quot;sam&quot; suffix is exist in <code>TFormat</code>.
  Althrough you add -DSEAN_HAS_ZLIB, I think you could check the code in library to see if it really works fine.</p>
  <p>If <code>SEQAN_HAS_ZLIB</code> is working, you could check other variables or functiones related with open bam process.These function are in <code>bam_file.h</code>, <code>formatted_file.h</code>, <code>fundamental_tag.h</code> and so on.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-12-31T09:11:26.517Z
Answers.LastActivityDate: 2023-01-01T08:45:26.457Z
Answers.LastEditDate: 2023-01-01T08:45:26.457Z
Answers.UserDisplayName: xin king
Body: <p>I am trying to open .bam files in C++ to extract reads occurring at specific scaffolds and loci. I essentially want to call &quot;samtools view sample.bam -o sample.sam scaffold:pos-pos&quot; from C++. I have tried system(&quot;samtools view sample.bam -o sample.sam scaffold:pos-pos&quot;) and can explain the issues if you want to know, but for the sake of conciseness I will leave the details out of this post. So, I am using SeqAn. When I run my script, I get the following error:</p>
  <pre><code>terminate called after throwing an instance of 'seqan::UnknownExtensionError'
    what():  Unknown file extension of copyPolly_SRR6511930.bam: iostream error
  </code></pre>
  <p>I am running my program on my institution's supercomputer clusters (Texas A&amp;M Terra cluster), to which I do <strong>not</strong> have sudo access. I have loaded the following modules:</p>
  <pre><code>GCC/9.3.0
  SeqAn/2.4.0-GCCcore-9.3.0
  zlib/1.2.11-GCCcore-9.3.0
  </code></pre>
  <p>Here is my minimal script:</p>
  <pre><code>#include &lt;seqan/bam_io.h&gt;
  #include &lt;zlib.h&gt;
  
  using namespace seqan;
  
  int main(){
      CharString bamFileName = &quot;copyPolly_SRR6511930.bam&quot;;
  // Open input file.
      BamFileIn bamFileIn(toCString(bamFileName));
      if (!open(bamFileIn, toCString(bamFileName))){
          std::cerr &lt;&lt; &quot;ERROR: Could not open &quot; &lt;&lt; bamFileName &lt;&lt; std::endl;
          return 1;
      }
  
      unsigned numUnmappedReads = 0;
  
      try{
          // Read header.
          BamHeader header;
          readHeader(header, bamFileIn);
  
          // Read records.
          BamAlignmentRecord record;
          while (!atEnd(bamFileIn)){
              readRecord(record, bamFileIn);
              if (hasFlagUnmapped(record))
                  numUnmappedReads += 1;
          }
      }
      catch (Exception const &amp; e)
      {
          std::cout &lt;&lt; &quot;ERROR: &quot; &lt;&lt; e.what() &lt;&lt; std::endl;
          return 1;
      }
      return 0;
  }
  </code></pre>
  <p>I have tried the following solutions:</p>
  <p>-Adding <code>#define SEQAN_HAS_ZLIB (1)</code> to my script's header</p>
  <p>-Adding <code>-D SEQAN_HAS_ZLIB</code> to my compilation such that it reads <code>g++ script.cpp -D SEQAN_HAS_ZLIB</code></p>
  <p>-Adding <code>-lz</code> to my compilation such that it reads <code>g++ script.cpp -lz</code></p>
  <p>I am not dead-set on using SeqAn for this project, so if anyone knows of another C++ library which  extracts specific reads from .bam files, please share.</p>
  <p>Update: I've never used HTSlib and it seems to have a fairly steep learning curve, so if anyone has an answer for the SeqAn error, that would still be very much appreciated.</p>
  
Tags: <sam><samtools><c++>
Title: Unable to open .bam file in C++ using SeqAn due to 'seqan::UnknownExtensionError'
ThreadText: # Unable to open .bam file in C++ using SeqAn due to 'seqan::UnknownExtensionError'
  
  (Asked by: username_0 on 2021-06-27)
  
  I am trying to open .bam files in C++ to extract reads occurring at specific scaffolds and loci. I essentially want to call "samtools view sample.bam -o sample.sam scaffold:pos-pos" from C++. I have tried system("samtools view sample.bam -o sample.sam scaffold:pos-pos") and can explain the issues if you want to know, but for the sake of conciseness I will leave the details out of this post. So, I am using SeqAn. When I run my script, I get the following error:
  
  ```
  terminate called after throwing an instance of 'seqan::UnknownExtensionError'
    what():  Unknown file extension of copyPolly_SRR6511930.bam: iostream error
  
  ```
  
  I am running my program on my institution's supercomputer clusters (Texas A&M Terra cluster), to which I do **not** have sudo access. I have loaded the following modules:
  
  ```
  GCC/9.3.0
  SeqAn/2.4.0-GCCcore-9.3.0
  zlib/1.2.11-GCCcore-9.3.0
  
  ```
  
  Here is my minimal script:
  
  ```
  #include <seqan/bam_io.h>
  #include <zlib.h>
  
  using namespace seqan;
  
  int main(){
      CharString bamFileName = "copyPolly_SRR6511930.bam";
  // Open input file.
      BamFileIn bamFileIn(toCString(bamFileName));
      if (!open(bamFileIn, toCString(bamFileName))){
          std::cerr << "ERROR: Could not open " << bamFileName << std::endl;
          return 1;
      }
  
      unsigned numUnmappedReads = 0;
  
      try{
          // Read header.
          BamHeader header;
          readHeader(header, bamFileIn);
  
          // Read records.
          BamAlignmentRecord record;
          while (!atEnd(bamFileIn)){
              readRecord(record, bamFileIn);
              if (hasFlagUnmapped(record))
                  numUnmappedReads += 1;
          }
      }
      catch (Exception const & e)
      {
          std::cout << "ERROR: " << e.what() << std::endl;
          return 1;
      }
      return 0;
  }
  
  ```
  
  I have tried the following solutions:
  
  -Adding `#define SEQAN_HAS_ZLIB (1)` to my script's header
  
  -Adding `-D SEQAN_HAS_ZLIB` to my compilation such that it reads `g++ script.cpp -D SEQAN_HAS_ZLIB`
  
  -Adding `-lz` to my compilation such that it reads `g++ script.cpp -lz`
  
  I am not dead-set on using SeqAn for this project, so if anyone knows of another C++ library which extracts specific reads from .bam files, please share.
  
  Update: I've never used HTSlib and it seems to have a fairly steep learning curve, so if anyone has an answer for the SeqAn error, that would still be very much appreciated.
  
  **Tags:** sam, samtools, c++
  
  
  ---
  
  ## 4 answers
  
  ### Answer by username_0 on 2021-06-29. Score: 3 | Accepted answer
  
  The helpdesk for my institution's supercomputers showed me the following flags:
  
  ```
  g++ seqan_test.cpp -DSEQAN_HAS_ZLIB -lz -lpthread
  
  ```
  
  This fixed the error. I'm now running into different errors and am suspicious that I am getting them because I am using SeqAn2 rather than SeqAn3. I'm going to try HTSlib.
  
  
  ---
  
  ### Answer by username_1 on 2021-06-28. Score: 3
  
  [Htslib](https://github.com/samtools/htslib) is the official library to access alignment files. It supports more formats. The typical way to read through a SAM/BAM/CRAM:
  
  ```cpp
  #include "htslib/sam.h"
  int main(int argc, char *argv[])
  {
    samFile *fp = sam_open(argv[1], "r");
    bam_hdr_t *h = sam_hdr_read(fp);
    bam1_t *b = bam_init1();
    while (sam_read1(fp, h, b) >= 0) {
      if (b->core.tid < 0) continue;
      int end_pos = bam_endpos(b);
      printf("%s\t%d\t%d\n", h->target_name[b->core.tid], b->core.pos, end_pos);
    }
    bam_destroy1(b);
    bam_hdr_destroy(h);
    sam_close(fp);
  }
  
  ```
  
  PS: the above started as a comment for "relevant but minor or transient information to" the original question. I didn't consider htslib a direct answer to the question, but because this had been changed to an answer, I added some details anyway. By the way, seqan2 has largely been discontinued. Use seqan3.
  
  > **Comments:**
  > **username_0**: Thank you for this! I tested this script and it worked. I will try to edit it to suit my needs. Also, I am using SeqAn2 because SeqAn3 is not available on my institution's supercomputers. If HTSlib doesn't pan out or if no one can figure this error out, I will email them to request that it be installed.  
  
  ---
  
  ### Answer by username_1 on 2021-06-28. Score: 1
  
  I'm not quite sure why assigning a string to a variable would lead to a complaint about file extensions. From the [SeqAn BAM tutorial](https://seqan.readthedocs.io/en/master/Tutorial/InputOutput/SamAndBamIO.html#tutorial-io-sam-bam-io), it looks like you need to use the *BamFileIn* object for reading BAM/SAM files:
  
  ```
  #include <seqan/bam_io.h>
  using namespace seqan;
  
  int main(){
      CharString bamFileName = "copyPolly_SRR6511930.bam";
      // Open input file, BamFileIn can read SAM and BAM files.
      BamFileIn bamFileIn;
      if (!open(bamFileIn, toCString(bamFileName))) {
          std::cerr << "ERROR: Could not open " << bamFileName << std::endl;
          return 1;
      }
      return 0;
  }
  
  ```
  
  Edit: after seeing the updated code... I've got nothing, sorry. It all seems like reasonable code that shouldn't be generating errors.
  
  > **Comments:**
  > **username_0**: Sorry about that, I did not include all of my code. I have updated my post with the rest of the code which led to my error.  
  
  ---
  
  ### Answer by username_2 on 2022-12-31. Score: 1
  
  I think `SEQAN_HAS_ZLIB` flag not works correctly as it supposed to be.
  So you can check the seqan library to see if `SEQAN_HAS_ZLIB` works correctly, and seek the reason of extension error.
  
  Seqan library is a mature library which using template to define general-purpose classes and structs.
  
  1.The reason for extension error.
  
  Step into the seqan library:
  
  When you call BamfileIn: `BamFileIn bamFileIn(toCString(bamFileName));`
  
  An object of BamFileIn struct is created, and `bamFileName` file is opened by the struct construction function.
  
  And when this object is created, some extra information is wrapped with the object to check the type of file.
  
  This information include suffix categories. In bam\_file.h header file.
  
  ```
  template <typename TDirection, typename TSpec>
  struct FileFormat<FormattedFile<Bam, TDirection, TSpec> >
  {
  #if SEQAN_HAS_ZLIB
      typedef TagSelector<
                  TagList<Bam,
                  TagList<Sam
                  > >
              > Type;
  #else
      typedef Sam Type;
  #endif
  };
  
  ```
  
  As you can see, if `SEQAN_HAS_ZLIB` is setted, "Bam" file will be added to TagSelector which will be wrapped into the BamFileIn struct and redefined by `typedef` as `format`.
  
  Then \_open() function is called to open Bam file, and extension is checked in guess\_format.h
  
  ```
  template <typename TFilename, typename TTagList>
  inline bool guessFormatFromFilename(TFilename const &fname, TagSelector<TTagList> &format)
  {
      typedef typename TTagList::Type TFormat;
  
      if (value(format) == -1 || isEqual(format, TFormat()))
      {
          // if tagId is set to -1 (auto-detect) or the current format (TFormatTag) then test for TFormatTag format
          if (guessFormatFromFilename(fname, TFormat()))
          {
              assign(format, TFormat());
              return true;
          }
      }
      return guessFormatFromFilename(fname, static_cast<typename TagSelector<TTagList>::Base &>(format));
  }
  
  template <typename TFilename, typename TFormat_>
  inline bool guessFormatFromFilename(TFilename const & fileName, Tag<TFormat_> /*formatTag*/)
  {
      typedef typename Value<TFilename>::Type                                     TValue;
      typedef ModifiedString<TFilename const, ModView<FunctorLowcase<TValue> > >    TLowcase;
      typedef Tag<TFormat_>                                                       TFormat;
  
      TLowcase lowcaseFileName(fileName);
      for (unsigned i = 0; i < sizeof(FileExtensions<TFormat>::VALUE) / sizeof(char*); ++i)
          if (endsWith(lowcaseFileName, lowerString(FileExtensions<TFormat>::VALUE[i])))
              return true;
  
      return false;
  }
  
  ```
  
  The above code is two instance of `guessFormatFromFilename` function.
  It's pretty complex here, but if `SEQAN_HAS_ZLIB` is set correctly, `TagSelector` will include bam and sam, and the first `guessFormatFromFilename` above will return true in if condition, and the error will not happend. But if `SEQAN_HAS_ZLIB` is not working, second formt `guessFormatFromFilename` function will be called and only "sam" suffix is exist in `TFormat`, and it will be exhibited by `FileExtensions<TFormat>::VALUE[i]`. You can using std::cerr to print the suffix on your screen to see if there is only "sam" suffix is exist in `TFormat`.
  Althrough you add -DSEAN\_HAS\_ZLIB, I think you could check the code in library to see if it really works fine.
  
  If `SEQAN_HAS_ZLIB` is working, you could check other variables or functiones related with open bam process.These function are in `bam_file.h`, `formatted_file.h`, `fundamental_tag.h` and so on.
  
  
  ---
  
  