cgeb2001's DokuWiki!

This is an old revision of the document!

You might find that in your daily computer work, there are some commands that are long and you execute quite often. Sometimes it may also be a series of commands that you execute in a row that you do quite often. It may be worth it to write a custom function that wraps these lines of code into a short, so that you do not have to type all of these commands anymore.

We have already discussed one such function in a previous section:

# easier viewing of tab separated files
function coltab {
    column -t -s $'\t' $1 | less -S
}

I will discuss here some more custom functions that I found are very useful in my daily workflow. To add these functions to your system, simply add them to your .bashrc

Selecting or removing sequences from a FASTA file

# fish out a sequence from a fasta file
function grabseq {
    fasta=$3
    to_grab=$2
    case "$1" in
        -s) seqtk subseq $fasta <(grep    "$to_grab" $fasta | sed -e 's/>//' -e 's/ .*//');;
        -l) seqtk subseq $fasta <(grep -f "$to_grab" $fasta | sed -e 's/>//' -e 's/ .*//');;
        *) echo -e "grabseq -s <seqname> <fasta> to grab a single sequence\ngrab_seq -l <seqlist> <fasta> to grab a list of sequences"
    esac
}

This is essentially a wrapper for the seqtk tool.

# select any sequence that has <pattern> in the header
# useful if you want to find a single sequence
$ grabseq -s <pattern> <FASTA> > <NEWFASTA>

# select a particular set of sequences that have <pattern> in their header and are in the <pattern_list> file
# useful if you want to multiple sequences
$ grabseq -l <pattern_list> <FASTA> > <NEWFASTA>

The next function does the opposite of grabseq. It will remove particular sequences from a FASTA file.

# remove a particular entry from a fasta file
function rmseq {
    fasta=$3
    to_rmv=$2
    case "$1" in
        -s) seqtk subseq $fasta <( grep ">" $fasta | grep -E -v "$to_rmv" | sed "s/>//" );;
        -l) seqtk subseq $fasta <( grep ">" $fasta | grep -v -f "$to_rmv" | sed "s/>//" );;
        *) echo -e "rmseq -s <seqname> <fasta> to remove a single sequence\nremove_seq -l <seqlist> <fasta> to remove a list of sequences"
    esac
}

# remove a particular sequence that has <pattern> in the header
$ rmseq -s <pattern> <FASTA> > <NEWFASTA>

# remove a particular set of sequences that have <pattern> in the header and are in the <pattern_list> file
$ rmseq -l <pattern_list> <FASTA> > <NEWFASTA>