loop for samtools to deal with multiple files

Question

[hfang4]

Hello,

I have a script to run samtools pileup for multiple files. The script worked if I ran one file at a time without loop, but it gave me error saying "Illegal variable name" with the loop. So, I think it is something wrong with the loop. I have to use the old version of samtools with the pileup for the downstream analysis. Please anyone help me out with the issue.

 #! /bin/csh
cd /share/gwascotton/GBS_data/pileup_trial/
source /usr/local/apps/samtools/samtools017a.csh
for i in $(ls *.bam);do
samtools pileup -c -l /share/gwascotton/GBS_data/merged_files.sorted.bam.pileup -f /share/gwascotton/GBS_data/Ghirsutum_458_v1.0.fa /share/gwascotton/GBS_data/pileup_trial/"$i" > /share/gwascotton/GBS_data/new/"$i".pileup
done ;

Thank you.

fhzh

Answer 1

[Kevin Blighe]

First get a listing of the BAM files and save in BAM.list, then:

#!/bin/bash
paste BAM.list | while read BAM ;
do
    samtools pileup -c -l /share/gwascotton/GBS_data/merged_files.sorted.bam.pileup -f /share/gwascotton/GBS_data/Ghirsutum_458_v1.0.fa /share/gwascotton/GBS_data/pileup_trial/"${BAM}" > /share/gwascotton/GBS_data/new/"${BAM}".pileup
done ;

[hfang4]

Thank you for the code. I think I have to use csh because I have to submit my job to high performance computing (HPC) server and using the old version of samtools which the source code give the path to it. I modified your code a little bit as attached below, but I got error "while: Expression Syntax". Any ideas? Thanks

 #!/usr/bin/csh

cd /share/gwascotton/GBS_data/pileup_trial/
source /usr/local/apps/samtools/samtools017a.csh

ls *.bam > /share/gwascotton/GBS_data/pileup_trial/BAM.list ;
paste BAM.list | while read BAM ; 
do
    samtools pileup -c -l /share/gwascotton/GBS_data/merged_files.sorted.bam.pileup -f /share/gwascotton/GBS_data/Ghirsutum_458_v1.0. "${BAM}" > /share/gwascotton/GBS_data/new/"${BAM}".pileup
done;

[GenoMax]

Every modern computing environment will support bash. So you should be able to use it. You will need to figure out what /usr/local/apps/samtools/samtools017a.csh is doing and replace that with equivalent bash commands/script.

[Kevin Blighe]

Oh dear, I noticed the old version of SAMtools as judged by the fact that you are using pileup, which has been replaced with mpileup.

If you really need to use the c-shell, then try:

#!/bin/csh
foreach BAM (`cat BAM.list`) ;
    samtools pileup -c -l /share/gwascotton/GBS_data/merged_files.sorted.bam.pileup -f /share/gwascotton/GBS_data/Ghirsutum_458_v1.0.fa /share/gwascotton/GBS_data/pileup_trial/"${BAM}" > /share/gwascotton/GBS_data/new/"${BAM}".pileup
end ;

NB - cat BAM.list is surrounded by backticks (not single quotes).

Answer 2

[dariober]

I think the csh shell doesn't support that for loop syntax (see for example https://www.cyberciti.biz/faq/csh-shell-scripting-loop-example/ ). I would use Bash unless you have to use csh.

[hfang4]

Thanks for the comments, but the link you provided can't be opened.

I tried bath script. It worked somehow because it produced .pileup files, but I also got error messages “/usr/local/apps/samtools/samtools017a.csh: line2: setenv: command not found” and “/usr/local/apps/samtools/samtools017a.csh: line8: syntax error: unexpected end of file”.

I think that because the bath script did not apply the samtools017a.csh which I wanted.

Anyway, thanks for your insight.

[GenoMax]

Link in @Dario's post fixed.

[Kevin Blighe]

That's correct, but the syntax for looping in c-shell is:

foreach VAR (`cat file.list`) ;
    echo "do something" ;
end ;

See my solution above.

Answer 3

[arfesta]

In terminal why not just try:

  cd /share/gwascotton/GBS_data/pileup_trial/
  ls *.bam > ~/Desktop/list.of.bams
  while read file; do samtools pileup -c -l /share/gwascotton/GBS_data/merged_files.sorted.bam.pileup -f 
      /share/gwascotton/GBS_data/Ghirsutum_458_v1.0.fa /share/gwascotton/GBS_data/pileup_trial/$file >
      /share/gwascotton/GBS_data/new/$file.pileup; done < ~/Desktop/list.of.bams

[hfang4]

Thanks for the idea. I run the following codes but got another error "Ambiguous output redirect". I think we are almost there. But any idea on the new error?

 #! /bin/csh

cd /share/gwascotton/GBS_data/pileup_trial/
source /usr/local/apps/samtools/samtools017a.csh
ls *.bam > /share/gwascotton/GBS_data/pileup_trial/list.of.bams |
while read file; do
samtools pileup -c -l /share/gwascotton/GBS_data/merged_files.sorted.bam.pileup -f /share/gwascotton/GBS_data/Ghirsutum_458_v1.0.fa /share/gwascotton/GBS_data/pileup_trial/$file > /share/gwascotton/GBS_data/new/$file.pileup;
done < list.of.bams;

[GenoMax]

See if putting this line in " helps in script above. "samtools pileup -c -l /share/gwascotton/GBS_data/merged_files.sorted.bam.pileup -f /share/gwascotton/GBS_data/Ghirsutum_458_v1.0.fa /share/gwascotton/GBS_data/pileup_trial/$file > /share/gwascotton/GBS_data/new/$file.pileup"

Answer 4

[Kevin Blighe]

As per my comment from my thread above, here are the solutions for both bash and csh:

First get a listing of the BAM files and save in BAM.list, then:

BASH

#!/bin/bash
paste BAM.list | while read BAM ;
do
        samtools pileup -c -l /share/gwascotton/GBS_data/merged_files.sorted.bam.pileup -f /share/gwascotton/GBS_data/Ghirsutum_458_v1.0.fa /share/gwascotton/GBS_data/pileup_trial/"${BAM}" > /share/gwascotton/GBS_data/new/"${BAM}".pileup
done ;

c-shell

#!/bin/csh
foreach BAM (`cat BAM.list`) ;
    samtools pileup -c -l /share/gwascotton/GBS_data/merged_files.sorted.bam.pileup -f /share/gwascotton/GBS_data/Ghirsutum_458_v1.0.fa /share/gwascotton/GBS_data/pileup_trial/"${BAM}" > /share/gwascotton/GBS_data/new/"${BAM}".pileup
end ;

NB - cat BAM.list is surrounded by backticks (not single quotes).

[GenoMax]

Please don't post duplicate answers. Instead edit your original answer with this material.

[hfang4]

Thanks for the scripts. I tried several files, but neither of them works. The bash one keeps running for ever and the csh one gave me the same error message "Ambiguous output redirect".

[Kevin Blighe]

It works for me in both BASH and CSH, and even with SAMtools and a redirect. Just to be sure, cat BAM.list is enclosed by backticks, not quotes/single quotes.

This is looking like a local issue with your system architecture / program installation. What are the contents of /usr/local/apps/samtools/samtools017a.csh?

Answer 5

[hfang4]

This final version of the script worked. Thank all for your suggestions, comments, codes and insight.

#! /bin/bash

cd /share/gwascotton/GBS_data/pileup_trial/
source /usr/local/apps/samtools/samtools017a.bash for i in $(ls *.bam);
do

samtools pileup -c -l /share/gwascotton/GBS_data/merged_files.sorted.bam.pileup -f /share/gwascotton/GBS_data/Ghirsutum_458_v1.0.fa /share/gwascotton/GBS_data/pileup_trial/${i} > /share/gwascotton/GBS_data/new/${i}.pileup

done ;