processing alignment reads in a specific region

[Lionward]

heey, I am very new to seqAn, so excuse me if the question is simple..

I had a python script for iterating over a bed file and fetching reads in a specific region and process their cigar string.

I would like to move to c++ to do that, I was wondering if it's possible to do it using seqAn, since what I found so far is itererating over the whole bam file, do you have a function similiar to pysam.Alignement.fetch(chr, start, end).
since I am only interested in the regions that I have in my bed file.

Thank you in advance!

[smehringer]

Hi @Lionward,

thanks for reaching out. Your question is very valid.

Unfortunately, the current seqan3 does not feature indexed BAM files, which is needed to extract a region. I will discuss with the team but I fear that there are currently no capacities to add this feature.

We have a former version of SeqAn, SeqAn2 that does have this feature.
You can check it out here:

• user guide to SeqAn2 https://seqan.readthedocs.io/en/main/Infrastructure/Use/index.html
• the SeqAn2 BAM IO tutorial: https://seqan.readthedocs.io/en/main/Tutorial/InputOutput/SamAndBamIO.html
• The function extracting a region viewRecords on a BAM file: Check the API documentation http://docs.seqan.de/seqan2/main/ and type viewRecords in the search field.

Best,
Svenja

[Lionward]

Hi @smehringer, thank you for your fast detailed reply, It would be really helpful to add this feature!
I am trying to get the seqan2 running... was folollwing the the documentation that you mentioned, but I encountered some issues.

So when I was using seqan3 it was working for the bamfile that I am using.

this is what I was using and it was working for the given bamfile:

#include <seqan3/core/debug_stream.hpp>
#include <seqan3/io/sam_file/all.hpp>
int main()
{
   std::string filename = "mybamfile.bam";
   seqan3::sam_file_input fin{filename};
   for (auto && record : fin)  // this will pass the record by reference (no copy)
   {
       seqan3::debug_stream << record.id() << '\n';  //  just for testing
       seqan3::debug_stream << record.sequence() << '\n';
       seqan3::debug_stream << record.flag() << '\n';
   }
}

now in the documentation of seqan2 there is the following guide code to read a bamfile.

#include <seqan/bam_io.h>

using namespace seqan2;

int main()
{
    CharString bamFileName = "mybamfile.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;
    }
}

first when I compile the current script I get the following error:
error: 'seqan2' is not a namespace-name

I downloaded the seqan library from conda using this link : https://anaconda.org/bioconda/seqan
I alos checked the version : seqan 2.4.0 hdfd78af_2 bioconda
and it seems to be correct..

When I use:
using namespace seqan; instead of using namespace seqan2;

the code compile correctly with some warning:

but the bam file is not read correctly so that I am getting the error of the If statement : ""ERROR: Could not open mybamfile.bam""

this is how I am compiling the code:

g++ -std=c++17 -g -O3 -fopenmp -DSEQAN_ENABLE_DEBUG=1 -DSEQAN_DISABLE_VERSION_CHECK=ON -lz -lpthread -o scripts/c++/cut_reads_in_TR_regions.io -I <path_to_seqan_in_conda> -I <path_to_zlib> scripts/c++/cut_reads_in_TR_regions.cpp

I've tried multiple c++ version 20,17,14 with the same behaviour.

Am I missing something there?

Would really appreciate the help.

Best
Lion

[eseiler]

Hey there,

the latest SeqAn2 release (2.4.0) is quite old. The namespace is still seqan instead of seqan2 for this version. Here are links to tutorials v2.4.0 and the api documentation v2.4.0.

SeqAn 2.4.0 will likely also have problems with newer compilers, hence it's probably best to use the seqan2 main branch - which isn't packaged.

Using CMake, there are two main ways that come to mind:

fetch_content

CMakeLists.txt

cmake_minimum_required (VERSION 3.18)
project (example_fetch_content CXX)

# Directory tree (just for reference)
# tree -L 2 .
# .
# ├── CMakeLists.txt
# └── hello_world.cpp

include (FetchContent)
FetchContent_Declare (
    seqan_fetch_content
    GIT_REPOSITORY "https://github.com/seqan/seqan"
    GIT_TAG "main"
    GIT_SHALLOW TRUE)
FetchContent_Populate (seqan_fetch_content)

find_package (SeqAn REQUIRED PATHS "${seqan_fetch_content_SOURCE_DIR}/util/cmake")

add_executable (hello_world hello_world.cpp)
target_link_libraries (hello_world PUBLIC seqan::seqan2)

as submodule

(See git submodule add part in the CMake comments for setup of a submodule).

CMakeLists.txt

cmake_minimum_required (VERSION 3.18)
project (example_submodule CXX)

# Setup submodule in bash console:
# mkdir submodules
# git submodule add --depth 1 https://github.com/seqan/seqan submodules/seqan
#
# Directory tree (just for reference)
# tree -L 2 .
# .
# ├── CMakeLists.txt
# ├── hello_world.cpp
# └── submodules
#     └── seqan

find_package (SeqAn REQUIRED PATHS "${CMAKE_CURRENT_LIST_DIR}/submodules/seqan/util/cmake")

add_executable (hello_world hello_world.cpp)
target_link_libraries (hello_world PUBLIC seqan::seqan2)

For both examples, I used
hello_world.cpp

#include <iostream>

#include <seqan/sequence.h>

int main()
{
    seqan2::String<char> text = "Hello World!\n";
    std::cout << seqan2::toCString(text);
}

to check if it works.
The documentation links for seqan2's main branch are tutorials main and api documentation main.

These are just minimal examples and will likely need to be incorporated into your existing CMake.
Don't hesitate to ask if you encounter difficulties.

Best,
Enrico