4
\$\begingroup\$

I have written a bash script that consists of multiple commands and Python scripts. The goal is to make a pipeline for detecting long non coding RNA from a certain input. Ultimately I would like to turn this into a app and host it on some bioinformatics website. One problem I am facing is using getopt tools in bash. I couldn't find a good tutorial that I understand clearly. In addition any other comments related to improving on the code is much appreciated.

#!/bin/bash
if [ "$1" == "-h" ]
then
    echo "Usage: sh $0 cuffcompare_output reference_genome blast_file"
    exit
else
wget https://github.com/TransDecoder/TransDecoder/archive/2.0.1.tar.gz && tar xvf 2.0.1 && rm -r 2.0.1
makeblastdb -in $3 -dbtype nucl -out $3.blast.out
grep '"u"' $1 | \
gffread -w transcripts_u.fa -g $2 - && \
python2.7 get_gene_length_filter.py transcripts_u.fa transcripts_u_filter.fa && \
TransDecoder-2.0.1/TransDecoder.LongOrfs -t transcripts_u_filter.fa
sed 's/ .*//' transcripts_u_filter.fa | grep ">" | sed 's/>//' > transcripts_u_filter.fa.genes
cd transcripts_u_filter.fa.transdecoder_dir
sed 's/|.*//' longest_orfs.cds | grep ">" | sed 's/>//' | uniq > longest_orfs.cds.genes
grep -v -f longest_orfs.cds.genes ../transcripts_u_filter.fa.genes > longest_orfs.cds.genes.not.genes
sed 's/^/>/' longest_orfs.cds.genes.not.genes > temp && mv temp longest_orfs.cds.genes.not.genes
python ../extract_sequences.py longest_orfs.cds.genes.not.genes ../transcripts_u_filter.fa longest_orfs.cds.genes.not.genes.fa
blastn -query longest_orfs.cds.genes.not.genes.fa -db ../$3.blast.out -out longest_orfs.cds.genes.not.genes.fa.blast.out -outfmt 6
python ../filter_sequences.py longest_orfs.cds.genes.not.genes.fa.blast.out longest_orfs.cds.genes.not.genes.fa.blast.out.filtered
grep -v -f longest_orfs.cds.genes.not.genes.fa.blast.out.filtered longest_orfs.cds.genes.not.genes.fa > lincRNA_final.fa
fi

Here is how I run it:

sh test.sh cuffcompare_out_annot_no_annot.combined.gtf /mydata/db/Brapa_sequence_v1.2.fa TE_RNA_transcripts.fa
\$\endgroup\$
  • 1
    \$\begingroup\$ By tar xvf 2.0.1 && rm -r 2.0.1 do you mean tar xvf 2.0.1.tar.gz && rm -r 2.0.1.tar.gz? \$\endgroup\$ – 200_success Feb 27 '15 at 21:44
  • \$\begingroup\$ For some reason when i wget the software it is downloaded as 2.0.1 and not as 2.0.1.tar.gz \$\endgroup\$ – upendra Feb 27 '15 at 21:54
4
\$\begingroup\$

This script looks like a wall of text and barely readable. Let's break it down from top to bottom.

Simplify the handling of -h

Since the if at the beginning to check if the first parameter is -h will exit if true, it would be simpler to drop the else part entirely:

if [ "$1" == -h ]
then
    echo "Usage: sh $0 cuffcompare_output reference_genome blast_file"
    exit
fi

I also dropped the quotes around the -h, they are not necessary.

Validate the command line arguments

The script expects 3 parameters, but doesn't validate if it really received them. I suggest to rewrite the top part like this:

help() {
    exitcode=$1
    echo "Usage: sh $0 cuffcompare_output reference_genome blast_file"
    exit $exitcode
}

[ "$1" == -h ] && help 0
[ $# -lt 3 ] && help 1

I added a help function to avoid duplication, because I wasn't too reuse the same help message when exiting the program in both cases:

  1. When the first param is -h, print help and exit with 0, meaning success
  2. When there are less than 3 params, print help and exit with 1, meaning failure

Instead of if statements I used && expressions to make it shorter and simpler. The result is equivalent.

Give the parameters meaningful names

The parameters $1, $2 and $3 are scattered here and there in the wall of text, and they are not exactly easy to remember. Give them proper names right after validation:

cuffcompare_output=$1
reference_genome=$2
blast_file=$3

And use these names throughout the script.

Taming the wall of text...

Add some line breaks. Group together lines that are most tightly related. For example, for the lines that are the continuation of a pipeline, put a line break before and after would make this easier to read. And when a pipeline gets too long, add some more line breaks (escaped with \) to make the important ending part more visible, for example when redirecting output to a file.

For example, I think this is more readable:

wget https://github.com/TransDecoder/TransDecoder/archive/2.0.1.tar.gz && tar xvf 2.0.1 && rm -r 2.0.1

makeblastdb -in $3 -dbtype nucl -out $3.blast.out

grep '"u"' $1 | \
gffread -w transcripts_u.fa -g $2 - && \
python2.7 get_gene_length_filter.py transcripts_u.fa transcripts_u_filter.fa && \
TransDecoder-2.0.1/TransDecoder.LongOrfs -t transcripts_u_filter.fa

sed 's/ .*//' transcripts_u_filter.fa | grep ">" | sed 's/>//' \
> transcripts_u_filter.fa.genes

cd transcripts_u_filter.fa.transdecoder_dir

sed 's/|.*//' longest_orfs.cds | grep ">" | sed 's/>//' | uniq \
> longest_orfs.cds.genes

grep -v -f longest_orfs.cds.genes ../transcripts_u_filter.fa.genes \
> longest_orfs.cds.genes.not.genes

sed 's/^/>/' longest_orfs.cds.genes.not.genes \
> temp && mv temp longest_orfs.cds.genes.not.genes

python ../extract_sequences.py longest_orfs.cds.genes.not.genes ../transcripts_u_filter.fa longest_orfs.cds.genes.not.genes.fa

blastn -query longest_orfs.cds.genes.not.genes.fa -db ../$3.blast.out -out longest_orfs.cds.genes.not.genes.fa.blast.out -outfmt 6

python ../filter_sequences.py longest_orfs.cds.genes.not.genes.fa.blast.out longest_orfs.cds.genes.not.genes.fa.blast.out.filtered

grep -v -f longest_orfs.cds.genes.not.genes.fa.blast.out.filtered longest_orfs.cds.genes.not.genes.fa > lincRNA_final.fa

As a small tip, a shorter equivalent of grep ">" | sed 's/>//' is sed -ne 's/>//p'.

\$\endgroup\$
3
\$\begingroup\$

There isn't much sense in putting the bulk of your script inside the else clause. Just put the fi after exiting prematurely.

if [ "$1" == "-h" ]; then
    echo "Usage: sh $0 cuffcompare_output reference_genome blast_file"
    exit
fi

You use temporary files a lot. Avoid them whenever possible. For example:

wget -O- https://github.com/TransDecoder/TransDecoder/archive/2.0.1.tar.gz | tar xzvf -

That also eliminates the problem of guessing what filename wget will decide to give to the downloaded file.

You have very long commands. It makes sense to put each pipeline in its own stanza for readability.

Develop a habit to double-quote every shell variable unless you have a reason to use single quotes or no quoting. Otherwise, Very Bad Things can happen if filenames with special characters are involved.

makeblastdb -in "$3" -dbtype nucl -out "$3.blast.out"

grep '"u"' "$1" | \
    gffread -w transcripts_u.fa -g "$2" - && \
    python2.7 get_gene_length_filter.py transcripts_u.fa transcripts_u_filter.fa && \
    TransDecoder-2.0.1/TransDecoder.LongOrfs -t transcripts_u_filter.fa

sed 's/ .*//' transcripts_u_filter.fa | grep ">" | sed 's/>//' > transcripts_u_filter.fa.genes

There's a cd command buried in there! Those are tricky since they alter the context of every command that follows. If you need to do a cd at all (and I'd avoid it if possible), I suggest putting the cd in its own subshell, so that its effect is scoped to the subshell. (It's less important here, since you only do one cd, but such side-effects lead to chaos in larger scripts.)

(
    cd transcripts_u_filter.fa.transdecoder_dir

    sed 's/|.*//' longest_orfs.cds | grep ">" | sed 's/>//' | uniq > longest_orfs.cds.genes

    # The rest of the commands here…
)

This should be a grep | sed > longest_orfs.cds.genes.not.genes pipeline:

grep -v -f longest_orfs.cds.genes ../transcripts_u_filter.fa.genes > longest_orfs.cds.genes.not.genes

sed 's/^/>/' longest_orfs.cds.genes.not.genes > temp && mv temp longest_orfs.cds.genes.not.genes
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.