ncbi fastq dump error in loop

Question

I have a very basic objective which I want to give list of id to a ncbi-srafastq tool kit which will prefetch the id which is in .sra file extension and then it will run fastq-dump on the file to split the files into fastq files, single file if there is single end file or paired end file if its paired end or two set of files.

This is my code

#!/bin/bash

input=$1

a=$(basename "$input")
id=${a/.id/}

# Download SRA data
prefetch "$id"

# Extract SRA ID from the downloaded data
#sraID=$(ls -1 | grep -v run.sh)
sraID=$(ls -1 | grep -v run.sh | sed 's/\.sra$//')
# Validate the downloaded SRA data
vdb-validate "$sraID"

# Dump fastq data from the SRA file
fastq-dump "$sraID"  --gzip

# Concatenate fastq files for read1, read2, and unpaired reads (if available)
for s in "$sraID"
do
  echo "$s"
  if [ "$s" != "$id" ]
  then
    [ -f "${s}_1.fastq" ] && cat "${s}_1.fastq" >> "${id}_R1.fastq"
    [ -f "${s}_2.fastq" ] && cat "${s}_2.fastq" >> "${id}_R2.fastq"
    [ -f "${s}.fastq" ] && cat "${s}.fastq" >> "${id}.fastq"
  fi
done

# Compress fastq files using pigz
#for i in "${id}"*.fastq
#do
  #echo "$i"
 # pigz "$i" & sleep 1
#done
#wait

# Remove the original uncompressed fastq file if both read1 and read2 files are present
[ -f "${id}_R1.fastq.gz" ] && [ -f "${id}_R2.fastq.gz" ] && [ -f "${id}.fastq.gz" ] && rm "${id}.fastq.gz"

# Test deleting intermediate folder
echo "Space before tmp folder delete"
echo "Removing tmp folder $sraID "
ls *
rm -rf "$sraID"*
ls *
echo 'done'

I run using as small example file as such

./run.sh SRR12332899

The error I get is this

2023-08-05T12:37:32 prefetch.3.0.6: Current preference is set to retrieve SRA Normalized Format files with full base quality scores.
2023-08-05T12:37:33 prefetch.3.0.6: 1) Downloading 'SRR12332899'...
2023-08-05T12:37:33 prefetch.3.0.6: SRA Normalized Format file is being retrieved, if this is different from your preference, it may be due to current file availability.
2023-08-05T12:37:33 prefetch.3.0.6:  Downloading via HTTPS...
2023-08-05T12:40:48 prefetch.3.0.6:  HTTPS download succeed
2023-08-05T12:40:48 prefetch.3.0.6:  'SRR12332899' is valid
2023-08-05T12:40:48 prefetch.3.0.6: 1) 'SRR12332899' was downloaded successfully
2023-08-05T12:40:48 prefetch.3.0.6: 'SRR12332899' has 0 unresolved dependencies
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Database 'SRR12332899.sra' metadata: md5 ok
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Table 'SEQUENCE' metadata: md5 ok
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Column 'ALTREAD': checksums ok
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Column 'QUALITY': checksums ok
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Column 'READ': checksums ok
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Column 'SPOT_GROUP': checksums ok
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Database 'sra/SRR12332899.sra' contains only unaligned reads
2023-08-05T12:40:48 vdb-validate.3.0.6 info: Database 'SRR12332899.sra' is consistent
2023-08-05T12:40:49 fastq-dump.3.0.6 err: item not found while constructing within virtual database module - the path 'sra' cannot be opened as database or table
fastq-dump quit with error code 3
sra
Space before tmp folder delete
Removing tmp folder sra
run.sh

sra:
SRR12332899.sra
run.sh
done

If I do just

prefetch SRR12332899 and then fastq-dump SRR12332899

I get a SRR12332899.fastq SRR12332899.sra

This is for one file since I want to download list of files using the script which is not working.

Any suggestion would be really helpful

Answer 1

The main problem is that your script is inadvertently treating a set of file and directory names as a single SRA accession number, all in one variable. You should update it to work with each accession separately and make sure they really are accessions.

---

On this line you're taking all the names of things in your working directory (minus the script), each on a separate line, and trimming .sra off the end of each line:

sraID=$(ls -1 | grep -v run.sh | sed 's/\.sra$//')

In your case that list was apparently just a directory named "sra". (Even though I don't have that directory on my end, you can see it in your output, from the first of two ls * at the end, where you have sra/SRR12332899.sra.) So later on it tries to fastq-dump sra, which gives you that "cannot be opened as database or table" error. (vdb-validate can take a directory as input, so it doesn't complain.) If I run your script from scratch in an empty directory with SRR12332899, it actually seems to run OK, aside from removing everything it creates at the end. But if you run it with a second accession you see the problem, since it then tries things like vdb-validate "several accessions and filenames separated by newlines".

I think a proper fix would involve a bunch of changes, but here are some tips, at least:

• use subdirectories to keep things organized so you don't have to juggle multiple types of file in one place
• set -x at the start of a script will show you each command as it happens, which can be very helpful for troubleshooting
• be careful with quotes and whitespace in bash. If you want to loop over each word in a variable you can do for thing in $var; do ...; done (unquoted) but if you run into filenames with whitespace, it breaks them up. Another pattern is command | while read thing; do ...; done which will go line-by-line. (In that case if you were working with a variable storing the command output, you might do something like echo "$output" | while..., quoted that time, since you'd want to preserve the newlines.)
• ls * will expand each file/directory name as separate arguments (so bash really ends up calling ls "sra" ...whatever else...), and when ls sees directory names is lists their contents too. If you just want to show all the files/directories you can leave off the asterisk altogether, or if you wanted to give directory names but not have it list their contents, you can use ls -d.
• you don't necessarily have to call ls to store the names as a block of text. You can loop over each filename, including those with spaces, directly in bash (e.g. for thing in *; do echo "path: '$thing'"; done)
• be careful with the rm command, especially combined with wildcards, and especially with -rf. If sraID happened to be empty somehow, think of what rm -rf "$sraID"* would actually mean. (And yes, scripts can remove themselves!)