running samtools in shell script

[zorph]

hi
I am trying to run this shell script

for a in {seq 1 22}
samtools view -h file.bam | awk '$3=="chr$a" || /^@/' | samtools view -S - -b -o Chr$a.bam
done

when I do it on the command line it works like a charm, but when i try to throw it in a shell it doesn't give any output. Any ideas?

[dpryan]

You should specify whether you're using bash, tcsh, or something else. I recall that for bash your syntax is wrong. Also, I imagine you should escape the $3 since that could otherwise cause some issues.

[Alex Renwick]

I don't have any help for why your script didn't run, but I can tell you that
samtools can do what you want by itself.

for a in {1..22}; do 
  samtools view -b file.bam chr$a -o Chr$a.bam
done

Or, better yet, pipe the result of that "samtools view file region" command to the next step of the analysis instead of creating a bunch of new bam files.

[zorph]

You should specify whether you're using bash, tcsh, or something else. I recall that for bash your syntax is wrong. Also, I imagine you should escape the $3 since that could otherwise cause some issues.

i am using bash. The entire command is below:

for a in {1..22}
do
      samtools view -h AD3.sorted.bam | awk '$3=="chr$a" || /^@/' | samtools view -S - -b -o chr$a.bam
done

[Alex Renwick]

For me, the awk command doesn't work even from the command line. The single quotes prevent bash from expanding "chr$a" to (e.g.) "chr1". Awk is looking for the third column to literally match "chr$a".

I did this to get the awk command to work:

awk -v a=$a '$3==sprintf("chr%s",a) || /^@/ '