Changing the strand in a Bam/Sam file

[krespim]

My sequencing data has been generated using the dUTP method, which means that when viewed in IGV or other browsers the reads are shown as being in the opposite strand from which they are generated. This is usually not an issue but now I need to prepare some figures (using the fantastic bioconductor package gviz) where the direction of the reads is the key. The question is, how to reverse the strandness of the alignments in a bam/sam file?

I have been using a combination of bamToBed | awk | bedToBam but I loose the information on split reads, that is, a split read will become multiple reads. There must a more elegant way of doing it but I am missing it.

[dpryan]

There are a couple ways to do that. If the version of awk on your system is actually gawk, you can use the bitwise operators (namely "or()") to test the strand and then swap things accordingly. Some systems come with mawk, which lacks that. So, you can always use python:

#!/usr/bin/env python
import csv
import sys

f = csv.reader(sys.stdin, dialect="excel-tab")
of = csv.writer(sys.stdout, dialect="excel-tab")
for line in f :
    #Don't try to modify the header!!!
    if(line[0][0] == "@") :
        of.writerow(line)
        continue
    #Swap strand
    if((int(line[1]) & 0x10) == 0x10) :
        line[1] = int(line[1]) - 0x10
    else :
        line[1] = int(line[1]) + 0x10
    of.writerow(line)

For paired-end data, you could swap the flag that indicates the mate's orientation in the same manner. If you happen to have pysam installed, this example code could be even shorter.

[krespim]

also in awk

Thanks a lot dpryan. In the mean time some folks in lab helped to come up with a awk solution:

#!/usr/bin/awk -f

BEGIN{
	FS="\t";OFS="\t";}
{
	if($2=="16"){
		$2="0";}
	else{
		if($2=="0"){
			$2="16";}
		}
print $_;
}

[dpryan]

That's pretty similar to the awk solution I had in mind. If you end up with paired-end reads or things with secondary alignments, then you end up having to check the actual bits in the flag. Otherwise, that'll work!

[krespim]

That's pretty similar to the awk solution I had in mind. If you end up with paired-end reads or things with secondary alignments, then you end up having to check the actual bits in the flag. Otherwise, that'll work!

Thankfully this is single read and I only have the flag "0" and "16" on that field