subsetting a bam file with specific alingment length

[joseph]

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

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:

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

[jbrwn]

rather than that perl, pipe it into this:

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]

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

[joseph]

rather than that perl, pipe it into this:

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]

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