subsetting a bam file with specific alingment length

joseph · 12-27-2011, 05:23 PM

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 -

jbrwn · 12-27-2011, 08:41 PM

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.

cedance · 12-28-2011, 01:45 AM

Both your links to biostar are broken. Are you sure its not at least 30?

joseph · 12-28-2011, 08:32 PM

Quote:

Originally Posted by jbrwn

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 · 12-28-2011, 08:39 PM

Quote:

Originally Posted by cedance

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

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 08:29 AM
Looking up specific position in bam	Ooinp	Bioinformatics	2	12-09-2011 02:58 PM
Get allele frequencies for specific coordinates from a .bam file	mehc	Bioinformatics	1	10-28-2011 01:05 PM
Tool which splits bam files into specific genomic intervals?	Alex Coventry	Bioinformatics	2	09-15-2011 03:23 AM
BAM : select only alignment of a defined length	NicoBxl	Bioinformatics	0	08-18-2011 04:29 AM

12-27-2011, 05:23 PM	#1
joseph Member Location: ca Join Date: Feb 2008 Posts: 39	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

12-27-2011, 08:41 PM	#2
jbrwn Member Location: Denver, CO Join Date: Mar 2011 Posts: 37	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.

12-28-2011, 01:45 AM	#3
cedance Senior Member Location: Germany Join Date: Feb 2011 Posts: 108	Both your links to biostar are broken. Are you sure its not at least 30?