Filtering out the spliced reads from bam file.

Question

[EVR]

Hi,

I have bam file and I would like to filter out spliced reads from the bam file. How can I achieve that directly on bam file. Kindly guide me. Thanks in advance

[Devon Ryan]

Write a script (likely with pysam) to do it.

[EVR]

What is pysam? I already tried with following code

 samtools view  scaffoldxxx.bam | awk '!($6 ~ /N/)' > skipped_scaffold_xxx.bam

But somehow the output is not in the bam format. Should I export it to bed format and later convert back to BAM format?

[Devon Ryan]

You'll need to include the header and pipe the output of that to samtools, since it's SAM format. I'm sure you can google "pysam".

[dariober]

The output you get is not in bam format for two reasons: 1) samtools view doesn't output the header by default and 2) It outputs in sam, not bam format by default.

[barslmn]

You can do the filtering in samtools.

 samtools view -e 'cigar!~"N"' scaffoldxxx.bam  -bo skipped_scaffold_xxx.bam

Or if you want to use awk you need to print out with header

 samtools view -h  scaffoldxxx.bam | awk '/^@/ {print;next} !($6 ~ /N/)' | samtools view -bo  skipped_scaffold_xxx.bam

[Alastair Kerr]

Also look at this previous thread for some ideas:

Answer 1

[WouterDeCoster]

I needed the inverse, selecting only spliced reads. Using python, save this as select_spliced_reads.py and execute as python only_spliced_reads.py -b input.bam -o spliced.bam

import pysam
from argparse import ArgumentParser


def main():
    args = get_args()
    get_spliced_reads(args)


def get_args():
    parser = ArgumentParser("")
    parser.add_argument("-b", "--bam", help="input bam", required=True)
    parser.add_argument("-o", "--output", help="output bam", required=True)
    parser.add_argument("-i", "--introns", help="minimal number of introns", default=1, type=int)
    return parser.parse_args()


def get_spliced_reads(args):
    bam_file = pysam.AlignmentFile(args.bam, "rb")
    out_file = pysam.AlignmentFile(args.output, "wb", template=bam_file)
    for read in bam_file.fetch():
        if sum(operation == 3 for (operation, _) in read.cigartuples) > args.introns:
            out_file.write(read)
    bam_file.close()
    out_file.close()


if __name__ == "__main__":
    main()

Answer 2

[Pierre Lindenbaum]

Using picard FilterSamReads http://broadinstitute.github.io/picard/command-line-overview.html#Overview and a javascript filter (not tested)

java -jar picard.jar FilterSamReads \
       I=input.bam \ 
       O=output.bam \
      JAVASCRIPT_FILE=filter.js

with filter.js

function accept(r) {
if(r.getReadUnmappedFlag()) return true;
var i,c= r.getCigar();
if(c==null) return true;
for( i=0;i< c. numCigarElements() ;++i) {
  if(  c. getCigarElement(i).getOperator().name().equals("N") ) return false;
  }
return true;
}

accept(record);

Answer 3

[Brian Bushnell]

Using the BBMap package:

reformat.sh in=mapped.bam out=filtered.bam maxdellen=50

You can set "maxdellen" to whatever length deletion event you consider the minimum to signify splicing, which depends on the organism.

[caggtaagtat]

Apparently, the maxdellen parameter doesn't work anymore

[predeus]

It is now known as dellenfilter. Very convenient way to get rid of spliced reads.