One of the most common types of bioinformatics analysis if differential gene expression using RNA-Seq (ribonucleic acid sequence) data. In this type of analysis, RNA molecules are extracted from tissue samples under two different conditions, such as two time points during development, or wild-type (normal) and mutant individuals. The amount of RNA transcribed from some genes will differ across the two conditions. From this we can infer that expression of those genes is somehow related to the conditions being compared.

The extracted RNA is sequenced, producing files containing one sequence for each RNA molecule (in theory), such as "ACUGGCAAUCGGAAAUA...". The differential analysis pipeline has the following high-level stages:

If the counts for a given gene differ significantly between one condition and another, then we conclude that this gene is regulated (turned on or off) in response to the condition. A complete pipeline to perform this analysis on an HPC cluster is available at https://github.com/auerlab/cnc-emdiff. A similar pipeline that will run fairly quickly on a single computer is available at https://github.com/auerlab/fasda/tree/main/Test. This latter is a set of test scripts for FASDA, a program that performs the comparison of counts for each gene.

Below is a simple shell script implementation of the AmrPlusPlus pipeline, which, according to their website, is used to "characterize the content and relative abundance of sequences of interest from the DNA of a given sample or set of samples".

People can use this pipeline by uploading their data to the developer's website, or by installing the pipeline to their own Docker container or Galaxy server.

In reality, the analysis is performed by the following command-line tools, which are developed by other parties and freely available:

The role of AmrPlusPlus is to coordinate the operation of these tools. AmrPlusPlus is itself a script.

If you don't want to be dependent on their web service, a Galaxy server, or their Docker containers, or if you would like greater control over and understanding of the analysis pipeline, or if you want to use the newer versions of tools such as samtools, you can easily write your own script to run the above commands.

Also, when developing our own pipeline, we can substitute other tools that perform the same function, such as Cutadapt in place of Trimmomatic, or Bowtie in place of BWA for alignment.

All of these tools are designed for "short-read" DNA sequences (on the order of 100 base pair per fragment). When we take control of the process rather than rely on someone else's pipeline, we open the possibility of developing an analogous pipeline using a different set of tools for "long-read" sequences (on the order of 1000 base pair per fragment).

For our purposes, we install the above commands via FreeBSD ports and/or pkgsrc (on CentOS and Mac OS X). Then we just write a Unix shell script to implement the pipeline for our data.

Note that this is a real pipeline used for research at the UWM School of Freshwater Science.

It is not important whether you understand genomics analysis for this example. Simply look at how the script uses loops and other scripting constructs to see how the material you just learned can be used in actual research. I.e., don't worry about what cutadapt and bwa are doing with the data. Just see how they are run within the pipeline script, using redirection, command line arguments, etc. Also read the comments within the script for a deeper understanding of what the conditionals and loops are doing.

#!/bin/sh -e

# Get gene fraction threshold from user
printf "Resistome threshold? "
read threshold

############################################################################
# 1. Enumerate input files

raw_files="SRR*.fastq"

############################################################################
# 2.  Quality control: Remove adapter sequences from raw data

for file in $raw_files; do
    output_file=trimmed-$file
    # If the output file already exists, assume cutadapt was already run
    # successfully.  Remove trimmed-* before running this script to force
    # cutadapt to run again.
    if [ ! -e $output_file ]; then
        cutadapt $file > $output_file
    else
        printf "$raw already processed by cutadapt.\n"
    fi
done

############################################################################
# 3. If sequences are from a host organism, remove host dna

# Index resistance gene database
if [ ! -e megares_database_v1.01.fasta.ann ]; then
    bwa index megares_database_v1.01.fasta
fi

############################################################################
# 4. Align to target database with bwa mem.

for file in $raw_files; do
    # Output is an aligned sam file.  Replace trimmed- prefix with aligned-
    # and replace .fastq suffix with .sam
    output_file=aligned-${file%.fastq}.sam
    if [ ! -e $output_file ]; then
        printf "\nRunning bwa-mem on $file...\n"
        bwa mem megares_database_v1.01.fasta trimmed-$file > $output_file
    else
        printf "$file already processed by bwa mem\n"
    fi
done

############################################################################
# 5. Resistome analysis.

aligned_files=aligned-*.sam
for file in $aligned_files; do
    if [ ! -e ${file%.sam}group.tsv ]; then
        printf "\nRunning resistome analysis on $file...\n"
        resistome -ref_fp megares_database_v1.01.fasta -sam_fp $file \
            -annot_fp megares_annotations_v1.01.csv \
            -gene_fp ${file%.sam}gene.tsv \
            -group_fp ${file%.sam}group.tsv \
            -class_fp ${file%.sam}class.tsv \
            -mech_fp ${file%.sam}mech.tsv \
            -t $threshold
    else
        printf "$file already processed by resistome.\n"
    fi
done

############################################################################
# 6. Rarefaction analysis?

I generally write a companion to every analysis script to remove output files and allow a fresh start for the next attempt. Many programs will detect when output already exists and either not attempt to rewrite it, or behave differently in other ways.

#!/bin/sh -e

rm -f trimmed-* aligned-* aligned-*.tsv megares*.fasta.*