This is a test version of Biostars. For the public version, visit https://www.biostars.org.
reads y-position in a bam viewer

Hi,

I m doing a BAM viewer. I wonder how to set read position in vertical axis ? I mean, given N reads , How can I display all of them to fit the space and allow smooth horizontal scrolling ?
So how to get y-coordinate given:

  • read width
  • read position
  • region viewport ( start, end )

I may need a pure function . Something like : 

int ypos(region_start, region_end, read_size, read_begin ... )

But I m not sure.. See example in IGV bellow to understand

IGV

igv bam position
• 2,952 views • link
written by sacha sacha

A genuine question. What would be 1 thing about your proposed viewer that would make a compelling use case for it when we have IGV/IGB/tablet on one end and asciigenome at the other.

• 0 views • link
written by GenoMax GenoMax

An application build with C++/Qt and seqan library under GPL license.

• 0 views • link
written by sacha sacha

Isn't that something you decide (based on what a user is doing in terms of clicking/selecting to zoom-in and out). You set the coordinates of the region user is looking at and that allows you to lookup how many reads there are in that windows and associated max depth? Not being a programmer perhaps I am missing something?

• 0 views • link
written by GenoMax GenoMax

See the picture example bellow. I would like to know where to draw the red read. ( bottom or top) enter image description here

• 0 views • link
written by sacha sacha

pack the rectangle has already answered that part. I assume you will also be providing a feature to "group reads according to read strand"? Then it may need to go at the bottom, if color represents read strand.

• 0 views • link
written by GenoMax GenoMax

You could look at the source of other genome browsers (e.g. samtools tview, IGV, etc) to find out how they approached this problem. Also, search for pack rectangles, you will find plenty of resources - I think displaying reads on a genome browser is a particular case of rectangle packing.

• 0 views • link
written by h.mon h.mon

It works ! Thanks guys for your help ! I just need some optimization now. 

       QList<QList<seqan::BamAlignmentRecord> > rows;

        seqan::BamAlignmentRecord record;
        seqan::readRecord(record, bamFileIn);

        bool newRow = true;

        for (auto& row : rows)
        {
            seqan::BamAlignmentRecord rec = row.last();
            if (int(record.beginPos  > int(rec.beginPos + seqan::length(rec.seq))))
            {
                row.append(record);
                newRow = false;
                break;
            }
        }
        if (newRow)
        {
            QList<seqan::BamAlignmentRecord> row;
            row.append(record);
            rows.append(row);
        }

enter image description here

• 0 views • link
written by sacha sacha

int(rec.beginPos + seqan::length(rec.seq)

don't use the read length, you'lll have to calculate the 'end' position using the cigar string.

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

1 answer

as said h.mon, pack the rectangle. One example of basic packing algorithm in my sources: 

    List<List<SAMRecord>> rows = new ArrayList<>();
private void add(final SAMRecord rec)
    {
    // pileup
    for(final List<SAMRecord> row:this.rows)
        {
        final SAMRecord last=row.get(row.size()-1);
        if(right2pixel.apply(last)+ this.minDistanceBetweenPairs > left2pixel.apply(rec)) continue;
        row.add(rec);
        return;
        }

    final List<SAMRecord>  row=new ArrayList<>();
    row.add(rec);
    this.rows.add(row);
               }
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Ok thanks ! I assume I must use this recursive 'add' function in the following way : 

>         for record in readBam(chromosom, start, end): 
>                add(record)
• 0 views • link
written by sacha sacha

Log in to answer this question.

More posts like this