SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
blasr read length filter? ewilbanks Bioinformatics 1 07-17-2014 09:28 AM
Is there an easy way to filter SAM/BAM files with NM>0 dontkme Bioinformatics 3 04-24-2014 05:44 AM
How to filter a SAM/BAM file by bp alisrpp Bioinformatics 5 01-17-2014 12:11 PM
how to filter CCS by number of passes (not by long read length)? metheuse Pacific Biosciences 1 08-29-2013 12:27 PM
Filter paired end BAM file based on iSize Leif Bergsagel Bioinformatics 2 12-16-2010 11:50 AM

Reply
 
Thread Tools
Old 05-08-2015, 05:58 AM   #1
krapulaxdoctor
Member
 
Location: Netherlands

Join Date: May 2015
Posts: 20
Default Filter BAM/SAM based on read length

Dear All,

I have a BAM file that contains reads with different length.
Could you please tell me how can I filter only those reads that are longer/shorter than a given value and write it in a new BAM file?

I have checked existing forum topics, but I did not find any that could help me.
I also tried NGSutils, but it just does not work ( installation failed, the commands/programs i.e. bamutils are "not found").

Unfortunately I have no programming/bioinfo background. I'm just a biologist who have to do everything , So if it is possible please be a bit more detailed.

Thank you in advance,
krapulaxdoctor is offline   Reply With Quote
Old 05-08-2015, 08:23 AM   #2
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

You can do that with Reformat:

reformat.sh in=x.sam out=y.sam minlength=50 maxlength=200

It works with bam as well if samtools is installed and in the path. Reformat is in the BBMap package. To install that, you just unzip it and untar it.
Brian Bushnell is offline   Reply With Quote
Old 05-08-2015, 09:38 AM   #3
krapulaxdoctor
Member
 
Location: Netherlands

Join Date: May 2015
Posts: 20
Default

Hi,

Thank you for the quick response.
I downloaded the BBMap tar from here: http://sourceforge.net/projects/bbma...e=typ_redirect

I used tar -xzf BBMap.tar.gz command.

Do I have to install it somehow? Because when i try to run it the way you wrote is not working. (I have samtools and Java installed )

I tried this command ( in the directory where the program was extracted):
sh /home/user/Downloads/bbmap/reformat.sh in=/home/user/Desktop/allbams/file.bam out=/home/user/Desktop/fileout.bam minlength=1 maxlength=41
then it gave me the following error:
/home/user/Downloads/bbmap/reformat.sh: 4: /home/user/Downloads/bbmap/reformat.sh: Syntax error: "(" unexpected

Do I have to put these programs to PATH? Or how exacty do I have to run them?

Sorry for the trouble,
krapulaxdoctor is offline   Reply With Quote
Old 05-08-2015, 09:43 AM   #4
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,795
Default

No install is required for BBMap tools. Brian uses bash shell by default.

Can you try:

Code:
$ /home/user/Downloads/bbmap/reformat.sh in=/home/user/Desktop/allbams/file.bam out=/home/user/Desktop/fileout.bam minlength=1 maxlength=41
If you use "sh" (this is ubuntu?) then it will use dash.
GenoMax is offline   Reply With Quote
Old 05-08-2015, 10:25 AM   #5
krapulaxdoctor
Member
 
Location: Netherlands

Join Date: May 2015
Posts: 20
Default

Thank you for the response GenoMax.
The code helped. However, I have samtools installed but the program gave an error message with my BAM files. Now i try to convert them to SAM to see if it works.
krapulaxdoctor is offline   Reply With Quote
Old 05-08-2015, 10:39 AM   #6
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,795
Default

Ah yes. You have to convert the BAM file to SAM. I only focused on the error you had seen.
GenoMax is offline   Reply With Quote
Old 05-08-2015, 10:48 AM   #7
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Quote:
Originally Posted by krapulaxdoctor View Post
Thank you for the response GenoMax.
The code helped. However, I have samtools installed but the program gave an error message with my BAM files. Now i try to convert them to SAM to see if it works.
Can you post the error message? If samtools is installed and in the path, it should work fine.
Brian Bushnell is offline   Reply With Quote
Old 05-08-2015, 11:20 AM   #8
krapulaxdoctor
Member
 
Location: Netherlands

Join Date: May 2015
Posts: 20
Default

Yes, Brian, here it is:

java -ea -Xmx200m -cp /home/user/Downloads/bbmap/current/ jgi.ReformatReads in=/home/user/Desktop/allbams/RNA.bam out=/home/user/Desktop/RNAout.bam minlength=1 maxlength=41
Executing jgi.ReformatReads [in=/home/user/Desktop/allbams/RNA.bam, out=/home/user/Desktop/RNA.bam, minlength=1, maxlength=41]

Found samtools.
Input is being processed as unpaired
[E::sam_parse1] missing SAM header
[W::sam_read1] parse error at line 3
[main_samview] truncated file.
Exception in thread "Thread-4" java.lang.RuntimeException: java.io.IOException: Broken pipe
at stream.ReadStreamByteWriter.run(ReadStreamByteWriter.java:31)
Caused by: java.io.IOException: Broken pipe
at java.io.FileOutputStream.writeBytes(Native Method)
at java.io.FileOutputStream.write(FileOutputStream.java:345)
at java.io.BufferedOutputStream.write(BufferedOutputStream.java:122)
at stream.ReadStreamByteWriter.writeSam(ReadStreamByteWriter.java:533)
at stream.ReadStreamByteWriter.processJobs(ReadStreamByteWriter.java:86)
at stream.ReadStreamByteWriter.run2(ReadStreamByteWriter.java:41)
at stream.ReadStreamByteWriter.run(ReadStreamByteWriter.java:27)
Exception in thread "main" java.lang.RuntimeException: Writing to a terminated thread.
at stream.ConcurrentGenericReadOutputStream.write(ConcurrentGenericReadOutputStream.java:198)
at stream.ConcurrentGenericReadOutputStream.addDisordered(ConcurrentGenericReadOutputStream.java:193)
at stream.ConcurrentGenericReadOutputStream.add(ConcurrentGenericReadOutputStream.java:97)
at jgi.ReformatReads.process(ReformatReads.java:892)
at jgi.ReformatReads.main(ReformatReads.java:46)

Hope it helps.
krapulaxdoctor is offline   Reply With Quote
Old 05-08-2015, 12:45 PM   #9
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Those errors indicate that the bam file is broken. Most of that text comes from Reformat, but these 3 lines:
Code:
[E::sam_parse1] missing SAM header
[W::sam_read1] parse error at line 3
[main_samview] truncated file.
...come from samtools. Are you able to convert it to a sam successfully?
Brian Bushnell is offline   Reply With Quote
Old 05-09-2015, 04:12 AM   #10
krapulaxdoctor
Member
 
Location: Netherlands

Join Date: May 2015
Posts: 20
Default

Yes,
samtools can convert them into SAM with the:
Quote:
$ samtools view -h -o out.sam in.bam
BBmap can deal with these SAM files.
krapulaxdoctor is offline   Reply With Quote
Old 05-11-2015, 11:30 AM   #11
krapulaxdoctor
Member
 
Location: Netherlands

Join Date: May 2015
Posts: 20
Default

Actually, you were right. There must be a problem with the header. It can be converted to SAM but I cant convert it back to BAM. Also, I cant convert it into BED. Do you know any way to fix the headers? ( unfortunately I have only these files with the crappy header ) Or to generate new header for the SAM files from scratch?
Sorry if I ask too many questions.
krapulaxdoctor is offline   Reply With Quote
Old 05-11-2015, 11:37 AM   #12
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,795
Default

Can you post the first few lines of your sam file?

Code:
$ samtools view -h your.bam | more
GenoMax is offline   Reply With Quote
Old 05-11-2015, 12:04 PM   #13
krapulaxdoctor
Member
 
Location: Netherlands

Join Date: May 2015
Posts: 20
Default

@SQ SN:chr1 LN:195471971
@SQ SN:chr2 LN:182113224
@SQ SN:chr3 LN:160039680
@SQ SN:chr4 LN:156508116
@SQ SN:chr5 LN:151834684
@SQ SN:chr6 LN:149736546
@SQ SN:chr7 LN:145441459
@SQ SN:chr8 LN:129401213
@SQ SN:chr9 LN:124595110
@SQ SN:chr10 LN:130694993
@SQ SN:chr11 LN:122082543
@SQ SN:chr12 LN:120129022
@SQ SN:chr13 LN:120421639
@SQ SN:chr14 LN:124902244
@SQ SN:chr15 LN:104043685
@SQ SN:chr16 LN:98207768
@SQ SN:chr17 LN:94987271
@SQ SN:chr18 LN:90702639
@SQ SN:chr19 LN:61431566
@SQ SN:chrM LN:16299
@SQ SN:chrX LN:171031299
@SQ SN:chrY LN:91744698
HISEQLN:122:HCW3JADXX:2:2109:6425:52017 272 chr1 3001923 0
43M * 0 0 TCTGATTATTATGTGTCAGGAGGAATTTCTTTTCTGGTCCAAT
JJJJJJJJIIIJJJJJJJJJJJJJJJJIJJHHHHHFFFFFCCC AS:i:0 XN:i:0 XM:i:0 XO:i:0
XG:i:0 NM:i:0 MD:Z:43 YT:Z:UU NH:i:20 CC:Z:= CP:i:48081495 XS:A:- HI:i:0
krapulaxdoctor is offline   Reply With Quote
Old 05-11-2015, 12:12 PM   #14
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Well, it's missing the HD line. Try adding this to the beginning of one of the sam files:
"@HD VN:1.3 SO:unsorted"

However, that's not required, so I'm not sure what the problem is. It might help if you encapsulate what you pasted with [ code ] [ / code ] (without the spaces) to make the tabs come through:
Code:
hello	world
vs
Code:
hello world
Brian Bushnell is offline   Reply With Quote
Old 05-11-2015, 12:12 PM   #15
krapulaxdoctor
Member
 
Location: Netherlands

Join Date: May 2015
Posts: 20
Default

And also, when I try to convert these BAM files into BED files, it generates a file, but the program (bedtools bamtobed)gives this error message:

terminate called after throwing an instance of 'std :: out_of_range'
what(): vector::_M_range_check
Aborted (core dumped)
krapulaxdoctor is offline   Reply With Quote
Old 05-11-2015, 12:14 PM   #16
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Oh, wait, is says "truncated" so presumably the problem is at the end of the file. Can you run "tail" on the file and post the last two lines?
Brian Bushnell is offline   Reply With Quote
Old 05-11-2015, 12:18 PM   #17
krapulaxdoctor
Member
 
Location: Netherlands

Join Date: May 2015
Posts: 20
Default

Quote:
Originally Posted by Brian Bushnell View Post
Oh, wait, is says "truncated" so presumably the problem is at the end of the file. Can you run "tail" on the file and post the last two lines?
How do I do this " tail " ?
Sorry im a beginner...
krapulaxdoctor is offline   Reply With Quote
Old 05-11-2015, 12:19 PM   #18
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

"tail file.sam"

That will print the last 10 lines to the console.
Brian Bushnell is offline   Reply With Quote
Old 05-11-2015, 12:24 PM   #19
krapulaxdoctor
Member
 
Location: Netherlands

Join Date: May 2015
Posts: 20
Default

Quote:
HISEQHI:525:HCYWJADXX:2:2213:8924:55099 256 * 942639 0 43M * 0 0 CAAAGGGCTGAGAAGCACTTGAAAAAATGTTCAACATCCTTAA CCCFFFFFHHHHHJJJJJJJJJJJJJJJJIIJJJJJJJJJJJJ AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:43 YT:Z:UU NH:i:20 CC:Z:chrX CP:i:128687718 XS:A:+ HI:i:17
HISEQLN:122:HCW3JADXX:2:2207:7052:25724 272 * 944767 0 43M * 0 0 TACTTACATATAATAAATAAATAAATAAATATTTTTTAAAAAA [email protected]@@ AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:43 YT:Z:UU NH:i:11 CC:Z:chr6 CP:i:52981629 XS:A:- HI:i:9
HISEQLN:121:HCYV3ADXX:1:1203:18633:64996 0 * 949324 043M * 0 0 CAGAACCCCTGAAATTGGCAAGATAGACGTCAGTGTTAGCAGA CCCFFFFFHHHHHJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ AS:i:-5 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:5G37 YT:Z:UU NH:i:20 CC:Z:chr6 CP:i:6419658 XS:A:+ HI:i:12
HISEQLN:122:HCW3JADXX:1:1112:13385:80114 272 * 949722 043M * 0 0 GGTGTCCGCTAGTGTCCTGAGGCCTGAGCGAGGGGCTCCTCTC ##A7'?DFD;BD:[email protected]?<7DD::@=1 AS:i:-2 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:11T31 YT:Z:UU NH:i:20 CC:Z:chr6 CP:i:71166409 XS:A:- HI:i:15
This is the last few lines...
krapulaxdoctor is offline   Reply With Quote
Old 05-11-2015, 12:32 PM   #20
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

Assuming all of the things that look like spaces are actually tabs (sorry, tabs often get replaced by spaces on the console), I don't see anything wrong with the sam file and I don't know what the problem is. It may have something to do with a negative number being detected where a positive number is expected, but I'm just speculating.

You could try Picard rather than Samtools, and see if you have better luck. Or, try the most recent version of Samtools, or else v0.1.19. Sometimes there's a problem with a specific version.
Brian Bushnell is offline   Reply With Quote
Reply

Tags
bam read filter

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 11:12 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2018, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO