Question: Filtering out the spliced reads from bam file.
1
gravatar for EVR
4.3 years ago by
EVR570
Earth
EVR570 wrote:

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

spliced-reads rna-seq bam • 3.0k views
ADD COMMENTlink modified 4.3 years ago by Brian Bushnell17k • written 4.3 years ago by EVR570

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

ADD REPLYlink written 4.3 years ago by Devon Ryan97k

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?

ADD REPLYlink modified 4.3 years ago • written 4.3 years ago by EVR570

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".

ADD REPLYlink written 4.3 years ago by Devon Ryan97k

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.

ADD REPLYlink written 4.3 years ago by dariober11k

Also look at this previous thread for some ideas:

ADD REPLYlink modified 4.3 years ago • written 4.3 years ago by Alastair Kerr5.2k
0
gravatar for Pierre Lindenbaum
4.3 years ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum131k wrote:

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);
ADD COMMENTlink modified 4.3 years ago • written 4.3 years ago by Pierre Lindenbaum131k
0
gravatar for Brian Bushnell
4.3 years ago by
Walnut Creek, USA
Brian Bushnell17k wrote:

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.

ADD COMMENTlink modified 4.3 years ago • written 4.3 years ago by Brian Bushnell17k

Apparently, the maxdellen parameter doesn't work anymore

ADD REPLYlink written 2.3 years ago by caggtaagtat1.3k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1403 users visited in the last hour