SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Convert merged BAM back to per lane BAM or FASTQ file danielsbrewer Bioinformatics 6 10-03-2013 07:29 AM
Looking up specific position in bam Ooinp Bioinformatics 2 12-09-2011 01:58 PM
Get allele frequencies for specific coordinates from a .bam file mehc Bioinformatics 1 10-28-2011 12:05 PM
Tool which splits bam files into specific genomic intervals? Alex Coventry Bioinformatics 2 09-15-2011 02:23 AM
BAM : select only alignment of a defined length NicoBxl Bioinformatics 0 08-18-2011 03:29 AM

Reply
 
Thread Tools
Old 12-27-2011, 04:23 PM   #1
joseph
Member
 
Location: ca

Join Date: Feb 2008
Posts: 39
Default subsetting a bam file with specific alingment length

Hi
I would like to extract a bam file with alignments length = 30.
I tried a script that I found on BioStar:
http://biostar.stackexchange.com/questions/12433/answer
Code:
samtools view -h in.bam | perl -lane '$l = 0; $F[5] =~ s/(\d+)[MX=DN]/$l+=$1/eg; print if $l = 30 or /^@/' | samtools view -bS - > out.bam
but, the out.bam still has alignments with length other than 30.
my in.bam was created by bowtie:
Code:
bowtie /indexes/dmelr5.42 -v 2 -k 5 --best --strata -S -t my.fastq | samtools view -uS -o my.bam -
Please help
Thanks
Joseph



http://biostar.stackexchange.com/questions/12433/answer
joseph is offline   Reply With Quote
Old 12-27-2011, 07:41 PM   #2
jbrwn
Member
 
Location: Denver, CO

Join Date: Mar 2011
Posts: 37
Default

rather than that perl, pipe it into this:
Code:
awk 'BEGIN{OFS=FS="\t"}/^@/{print}!/^@/{if(length($10)==30) print}'
it'll print the header and any line where the sequence (column 10) is 30 characters.
jbrwn is offline   Reply With Quote
Old 12-28-2011, 12:45 AM   #3
cedance
Senior Member
 
Location: Germany

Join Date: Feb 2011
Posts: 108
Default

Both your links to biostar are broken. Are you sure its not at least 30?
cedance is offline   Reply With Quote
Old 12-28-2011, 07:32 PM   #4
joseph
Member
 
Location: ca

Join Date: Feb 2008
Posts: 39
Default

Quote:
Originally Posted by jbrwn View Post
rather than that perl, pipe it into this:
Code:
awk 'BEGIN{OFS=FS="\t"}/^@/{print}!/^@/{if(length($10)==30) print}'
it'll print the header and any line where the sequence (column 10) is 30 characters.
thanks. It worked.
joseph is offline   Reply With Quote
Old 12-28-2011, 07:39 PM   #5
joseph
Member
 
Location: ca

Join Date: Feb 2008
Posts: 39
Default

Quote:
Originally Posted by cedance View Post
Both your links to biostar are broken. Are you sure its not at least 30?
I changed (print if $l < 60) to (print if $l = 30)
here is the link again
http://biostar.stackexchange.com/que...pecific-length
joseph is offline   Reply With Quote
Reply

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:02 PM.


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