[Seqan-dev] AlignedReadStore

[Ruolin Liu <rliu0606@gmail.com>]

Dear seqan developers, 

I realized that they might be a bug in the read counting example of the RNA-Seq tutorial. The beginPos and endPos of the AlignedReadStoreElement are in the gapped space of the reference genome while the Gff Records beginPos and endPos are supposed to be in the ungapped space of the reference genome. This relates to my observations that the AlignedReadStoreElement does not actually contain any coordinates from the ungapped space. This makes the beginPos and endPos of them are different from the beginPos and endPos of the BamAlignmentRecord.  And it creates a lot of confusion. I would actually suggest adding two extra fields for AlignedReadStoreElement, such as beginSourcePos and endSourcePos. 

Below is a copy of the code from the example. I commented on some lines. 

void countReadsPerGene(String<unsigned> & readsPerGene, String<TIntervalTree> const & intervalTrees, TStore const & store)
{
    resize(readsPerGene, length(store.annotationStore), 0);
    String<TId> result;
    int numAlignments = length(store.alignedReadStore);

    // iterate aligned reads and get search their begin and end positions
    SEQAN_OMP_PRAGMA(parallel for private (result))
    for (int i = 0; i < numAlignments; ++i)
    {
        TAlignedRead const & ar = store.alignedReadStore[i];
        TPos queryBegin = _min(ar.beginPos, ar.endPos); // In gapped space
        TPos queryEnd = _max(ar.beginPos, ar.endPos); // In gapped space

        // search read-overlapping genes
        findIntervals(result, intervalTrees[ar.contigId] /*Ungapped space*/, queryBegin, queryEnd);

        // increase read counter for each overlapping annotation given the id in the interval tree
        for (unsigned j = 0; j < length(result); ++j)
        {
            SEQAN_OMP_PRAGMA(atomic)
            readsPerGene[result[j]] += 1;
        }
    }
}

Best,

Ruolin