read strand

[m_elena_bioinfo]

Hi NGS users,
is there a way (by using bioperl module of samtools) to extract the information about the strand of the reads using a bam alignment?
Thanx a lot!
ME

[Jose Blanca]

I don't know how to do it in perl, but I've done it in python using pysam. You can take a look at:

http://github.com/JoseBlanca/franklin/blob/master/franklin/snv/snv_annotation.py

The code that does the job is around lines 140 to 170 in the function _snvs_in_bam.

[kmcarr]

Hi NGS users,
is there a way (by using bioperl module of samtools) to extract the information about the strand of the reads using a bam alignment?
Thanx a lot!
ME

Yes, here is some sample code:

#!/usr/bin/perl

use strict;
use warnings;

use Bio::DB::Sam

# Instantiate a new Bio::DB::Sam object, passing it the BAM file and the reference FASTA.

my $sam = Bio::DB::Sam ->new(-bam=>"myPath/myBamFile.bam",
                             -fasta=>"myPath/myFastaFile.fa");


# Create a collection of alignment objects.  The arguments for get_features_by_location are seq_id (e.g. chrom name), start and end (1 based).


my @alignments = $sam->get_features_by_location(-seq_id=>'chr1',
                                                -start=>500,
                                                -end=>15000);


# Iterate over the aligments, querying each one for the strand of the read.

for my $a (@alignments) {
	
	my $strand=$a->strand;
	
	# Do something with strand info.
	
	# You can also get seq_id, start, end, cigar string and much more about each alignment.
	
}

[m_elena_bioinfo]

Thank you very much for the rapid answers! Both are very useful!

@kmcarr:in module Bio:B::Sam i found two different thing:
my $strand=$a->strand;
my $strand-query=$a->query->strand;

What is the different? The query strand (the second variable) is the strand of the read?

[kmcarr]

Thank you very much for the rapid answers! Both are very useful!

@kmcarr:in module Bio:B::Sam i found two different thing:
my $strand=$a->strand;
my $strand-query=$a->query->strand;

What is the different? The query strand (the second variable) is the strand of the read?

By my reading of the APIs they are both giving you the same information, just from slightly different points of view. In the first case you are calling the strand method on a Bio:B::BAM::Alignment object, and in the second case you are calling the strand method within a Bio:B::BAM::Query object.

In the first case strand reports whether the query sequence matches the forward (+1) or reverse (-1) strand of the reference. The second method reports whether or not the query sequence itself (+1) or the reverse complement of the query sequence (-1) matches the forward strand of the reference.

[m_elena_bioinfo]

Great Kmcarr, thanks for all the help you gave me.
I tried to compare the two codes by extracting the strand of the reads in both methods. The results were the same!

So I think that the only difference is the context and the module of bioperl that you use in your script!