bwa mem: fail to locate the index files

Question

[little_more]

Hi all,

I have a problem with running bwa in parallel. I'm using this script (as kindly suggested in this topic):

# reference
hg1kv37=/dir/to/reference/and/index/files/human_g1k_v37.fasta

# create a function 
bwatosam() {
  id="$1"

  bwa mem -t 10 -R "@RG\tID:$id\tSM:$id\tPL:ILLUMINA\tLB:$id_exome" -v 1 -M $hg1kv37 "$id"_R1_fq.gz "$id"_R2_fq.gz \
  | samtools view -Sb -o "$id".bam
}
export -f bwatosam

# --plus enables {%...} to remove a substring    
ls *_R1_fq.gz |
  parallel --plus bwatosam {%_R1_fq.gz}

And I keep getting the following error:

[E::bwa_idx_load_from_disk] fail to locate the index files

Although the contents of a /dir/to/reference/and/index/files/ are:

human_g1k_v37.fasta.bwt  
human_g1k_v37.dict
human_g1k_v37.fasta.fai  
human_g1k_v37.fasta
human_g1k_v37.fasta.amb    
human_g1k_v37.fasta.ann        
human_g1k_v37.fasta.idx.1.bt2
human_g1k_v37.fasta.idx.4.bt2
human_g1k_v37.fasta.idx.rev.1.bt2                  
human_g1k_v37.fasta.idx.rev.2.bt2      
human_g1k_v37.fasta.idx.2.bt2      
human_g1k_v37.fasta.idx.3.bt2      
human_g1k_v37.fasta.pac
human_g1k_v37.fasta.sa

The index files were created some time ago but they should be OK because if I run a loop over all files (using a list with their names):

for line in `cat ids_list.txt`
do
bwa mem \
-t 10 \
-M $hg1kv37 \
"$line"_R1_fq.gz \
"$line"_R2_fq.gz \
-v 1 -R '@RG\tID:$line\tSM:$line\tPL:ILLUMINA\tLB:$line_exome' \
-M | samtools view -Sb - > $line.bam
done

everything works. Any ideas?

[from the mountains]

try option -v 4 to get better debugging output.

Answer 1

[ATpoint]

As far as I know parallel does not see variables defined outside of the function that it runs. I therefore always add outside variables like:

INDEX=/foo/bar/index.fa

function BWA {

  Basename=$1
  Index=$2

  bwa mem ${Index} ${Basename}R1.fastq.gz ${Basename}R2.fastq.gz > (...)

}; export -f BWA

ls *R1.fastq.gz | awk -F "R1" '{print $1 | "sort -u"}' | parallel "BWA {} ${INDEX}"

[Pierre Lindenbaum]

and as specified in the previous answers, you should use a workflow manager : C: using GNU parallel for bwa mem and samtools

[little_more]

it works now, thank you very much!

Answer 2

[ole.tange]

GNU Parallel does not run its jobs in the same shell session as it is started. That means that non-exported variables and non-exported functions cannot be seen.

However, env_parallel can make a copy of the environment in the current session and pass that to GNU Parallel. By using env_parallel --session only new definitions will be given to GNU Parallel.

# start recording new environment names
env_parallel --session

# reference
hg1kv37=/dir/to/reference/and/index/files/human_g1k_v37.fasta

# create a function 
bwatosam() {
  id="$1"

  bwa mem -t 10 -R "@RG\tID:$id\tSM:$id\tPL:ILLUMINA\tLB:$id_exome" -v 1 -M $hg1kv37 "$id"_R1_fq.gz "$id"_R2_fq.gz |
    samtools view -Sb -o "$id".bam
}
# export is not needed due to --session
# export -f bwatosam

# --plus enables {%...} to remove a substring    
ls *_R1_fq.gz |
  env_parallel --plus bwatosam {%_R1_fq.gz}

# finish the session (not needed if this is the end of the script)
env_parallel --end-session