This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Samtools: How To Find Tandem Hits Of Same Query Sequence In An Indexed Bam/Sam File?

I would like to know how can I find tandem hits of the same query sequence in an indexed bam/sam file?

By this I mean hits of the same sequence in the same location, one hit next to each other, non-overlapping, given a maximum distance D between them, with no other hits from other sequences in the BAM file in between. I assume that given this would be queried in an indexed BAM files, these should be two consecutive entries. As a graphical example below, I am interested in finding the two consecutive hits of sequence 'aaaaaaaaaa' in the target genome:

----------------------------------------------------------------------------
    xxxxxxxxx  aaaaaaaaaa       aaaaaaaaaa       bbbbbbbbb    ccccccccc
  ddddddddd                                       eeeeeeeeeeeee
samtools bam
• 3,827 views • link
updated by Fred Fred • written by Ahdf-Lell-Kocks • Ahdf-Lell-Kocks

It's not clear to me. Can you give a graphical example please ?

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

It sounds like: find all non-overlapping intervals in the BAM which are not more distant than D from each other. You are correct that these should be adjacent in a sorted file, but you have to account for overlaps if these can exist.

• 0 views • link
written by Michael Dondrup Michael Dondrup

Are you really interested in finding repeated hits, or is there a biological question that you are trying to answer?

• 0 views • link
written by Sean Davis Sean Davis

how about dummying up a sam file which has some such hits in it?

• 0 views • link
written by Malcolm.Cook Malcolm.Cook

I still don't get "hits of the same sequence in the same location". What is the "sequence" ? what is the "hit" ?. Can your replace those words with "Reference Sequence" and "Short Read" please.

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Pierre, agreed more detail would be valuable. My "solution" below assumes "sequence" is that of the "Short Read" and "hit" is the "alignment" (being one row in the sam file). Hoping I'm right....

• 0 views • link
written by Malcolm.Cook Malcolm.Cook

1 answer

I'm not sure to understand exactly what you want to do, but the following script may help.

It takes a sorted input sam file and outputs the lines when the consecutive ones corresponds to identical sequences with different coordinates.

#!/usr/bin/perl -w
use strict;
my $prevChr   = -1;
my $prevStart = -1;
my $prevSequence = -1;
my $prevLine = -1;

while(<>){
    if(!/^@/){
        chomp;
        my @tab = split(/\t/);
        my ($chr,$start,$sequence) = @tab[2,3,9];
        if($prevStart ne -1 && $prevChr ne -1 && 
           $prevSequence eq $sequence && 
           $prevChr eq $chr && $start != $prevStart){
            print $prevLine."\n";
            print $_."\n";
        }

        $prevChr   = $chr;
        $prevStart = $start;
        $prevSequence = $sequence;
        $prevLine = $_;
    }
}

You just have to do something like that:

$ samtools view sorted_file.bam | script.pl
• 1 views • link
written by Fred Fred

Log in to answer this question.

More posts like this