2

I was working on an array job for a small pipeline, and I happened to need a way to execute a specific command based on file size. I found this post and similar which describe how to do it. At the moment I'm using the following:

find $d/*.fasta -size +100M -exec sh -c '
     chromap -i -r $1 -o $1.index

     chromap --preset hic -x $1.index -r $1 -1 $d/hi-c/${ID}_1.fq.gz -2 $d/hi-c/${ID}_2.fq.gz --SAM -o /dev/stdout -t 48 | \
         samtools view -bS -@ 48 | samtools sort -n -@ 48 | samtools view -h | sed -e "s/\/.//" | samtools view -bS -o ${ID}.bam -@ 48
     ' sh {} \;

which, aside from the bioinfo and tools used that all work, it appears to execute only the first command line — chromap -i -r $1 -o $1.index. Then, for some reason, once it gets into the second set of instructions the script returns the following:

Cannot find sequence file /hi-c/_1.fq.gz

indicating it does not have knowledge of the environment variables I have successfully used so far, or it cannot compute more than two operations at the time? I have no clues... I tried also something more simple e.g.

mkdir $d/scaffolding
find $d/*.fasta -size +100M -exec sh -c '
     chromap -i -r $1 -o $1.index && mv $1 $1.index $d/scaffolding
' sh {} \;

but Bash complains about: mv: the destination '/scaffolding' is not a directory.

What should I do to get either one (or both) to work? Am I missing something, please if someone has some insights on this issue let me know! Thanks in advance.

2 Answers 2

8
export d ID
find -H "$d"/*.fasta -prune -size +100M -exec bash -o pipefail -c '
  ret=0
  for file do
    chromap -i -r "$file" -o "$file.index" &&
      chromap --preset hic \
              -x "$file.index" \
              -r "$file" \
              -1 "$d/hi-c/${ID}_1.fq.gz" \
              -2 "$d/hi-c/${ID}_2.fq.gz" \
              --SAM -o /dev/stdout -t 48 |
        samtools view -bS -@ 48 |
        samtools sort -n -@ 48 |
        samtools view -h |
        sed "s:/.::" |
        samtools view -bS -o "${ID}.bam" -@ 48 ||
      ret=$?
    done
    exit "$ret"
  ' bash {} +

Addressing several issues:

  • unless you export $d and $ID the sh started by find won't be able to see them
  • Without -H, for those fasta files that are symlinks, find would check the size of the symlink rather than that of the fasta file
  • Without -prune, for those that are directories, find would descend into them. You'd probably want to exclude them altogether with ! -type d or only include regular files with -type f (which includes symlinks to regular files with -H). Adding those checks does not remove the need for -prune (or -maxdepth 0 if using GNU find or compatible).
  • You forgot to quote the parameter expansions. You can get away without in zsh, but not sh nor bash.
  • You're not checking the status of commands before running the next one. Here I've added the pipefail option so that pipelines return failure if any of the commands failed. That's now a standard sh option, but there are some sh implementations such as dash that still don't support it, so I've replaced sh with bash though you could use any other sh implementation that supports pipefail.
  • You had single quotes within single quotes, which meant the s/\/.// was outside quotes, which meant it was the same as s//.// which is not valid sed code. Here we're using double quotes for the inner quotes, and switch from / to : for the separator to avoid having to escape the /. Note that s:/.:: removes the first occurrence of / followed by any single character with nothing. To remove /. literally, that's s:/\.:: and s:/\.::g to remove all occurrences.
  • with the -exec sh ... {} ';', you're running one sh per file and also the failure of any of those is not reported. Replacing with -exec sh ... {} + addresses both.

If you could switch to zsh, then most of those issues could be easily avoided.

set -o pipefail
errors=0
for file (*.fasta(N-.LM+100)) {
  chromap -i -r $file -o $file.index &&
      chromap --preset hic \
              -x $file.index \
              -r $file \
              -1 $d/hi-c/${ID}_1.fq.gz \
              -2 $d/hi-c/${ID}_2.fq.gz \
              --SAM -o /dev/stdout -t 48 |
        samtools view -bS -@ 48 |
        samtools sort -n -@ 48 |
        samtools view -h |
        sed 's:/.::' |
        samtools view -bS -o ${ID}.bam -@ 48 ||
      (( error++ ))
}
  • zsh doesn't have that misfeature of the Bourne shell whereby unquoted expansions are subject to split+glob so you can omit the quotes.
  • zsh globs have most of find's feature builtins so you don't have to combine shell globs (your *.fasta) with find.
  • that means, no need for nested quotes, that means no need to run a separate shell and the need to export variables to it, that means we can more easily return an error count.
  • zsh does support pipefail (contrary to some sh implementations).
3
  • Does the x && y || z not still suffer from “run z if x succeeds but y fails”? I think I would expect if x; then y; else z; fi instead, but I probably misunderstood. Commented Aug 10 at 19:28
  • @D.BenKnoble well no. Here we want to run z unless x and y succeed. It does read like "x succeeds and y succeeds or else run z (and here we don't care about z's success). With if, that would have to be if ! { x && y; }; then z; fi. Or if ! x || ! y; then z; fi. Commented Aug 10 at 21:31
  • Got it, thanks! Commented Aug 11 at 13:47
3

You assume a couple of shell variables in your sh -c inline script are inherited via the environment. These are d and ID. If the shell does not export these into the environment before find is called, the sh -c script will not see them and will substitute empty strings in their places.

So, before calling find, ensure that you have

export d ID

... in your script.

Alternatively, set them for the find command when you call it, like so:

d=$d ID=$ID find "$d"/*.fasta -size +100M -exec ...

Also remember to double quote all variable expansions in your scripts:

export d ID

find "$d"/*.fasta -size +100M -exec sh -c '
    chromap -i -r "$1" -o "$1.index"

    chromap --preset hic -x "$1.index" -r "$1" -1 "$d/hi-c/${ID}_1.fq.gz" -2 "$d/hi-c/${ID}_2.fq.gz" --SAM -o /dev/stdout -t 48 |
    samtools view -bS -@ 48 |
    samtools sort -n -@ 48 |
    samtools view -h | sed -e 's/\/.//' |
    samtools view -bS -o "${ID}.bam" -@ 48' sh {} \;
0

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .