#
# A script that fully automates the RNA-seq analysis of UHR/BHR samples.
#

# Stop on any error.
set -ueo pipefail

# This is the reference genome.
REF=annot/ERCC92.fa

# These are the coordinates of the genes.
GTF=annot/ERCC92.gtf

# This the name of the index
IDX=annot/ERCC92.fa

# This is the URL for the data
URL=http://data.biostarhandbook.com/rnaseq/data/brain.tar.gz

# Download and unpack the data
echo "Downloading the data"
curl -s $URL | tar zxv

# Switch to the new directory.
cd brain

# This keeps track of the messages printed during exectuion.
RUNLOG=runlog.txt
echo "Run by `whoami` on `date`" > $RUNLOG

# Create output folder
mkdir -p bam

# Exit this script on any error.
set -euo pipefail


# Build the annotation index
hisat2-build $REF $IDX 1>> $RUNLOG 2>> $RUNLOG

for SAMPLE in HBR UHR;
do
    for REPLICATE in 1 2 3;
    do
        # Build the name of the files.
        R1=reads/${SAMPLE}_${REPLICATE}_R1.fq
        R2=reads/${SAMPLE}_${REPLICATE}_R2.fq
        BAM=bam/${SAMPLE}_${REPLICATE}.bam

        # Run the aligner.
        echo "Aligning: $BAM"
        hisat2 $IDX -1 $R1 -2 $R2 2>> $RUNLOG | samtools sort > $BAM 2>> $RUNLOG
        samtools index $BAM
    done
done

echo "Counting features with: $GTF"
featureCounts -a $GTF -g gene_name -o counts.txt  bam/HBR*.bam  bam/UHR*.bam 2>> $RUNLOG

echo "Generating simple counts."
cat counts.txt | cut -f 1,7-12 > simple_counts.txt

echo "Downloading DESeq1 script"
curl -s -O http://data.biostarhandbook.com/rnaseq/code/deseq1.r

echo "Generating results into: results_deseq1.txt "
cat simple_counts.txt | Rscript deseq1.r 3x3 > results_deseq1.txt  2>> $RUNLOG

echo "Downloading the DESeq2 script"
curl -s -O http://data.biostarhandbook.com/rnaseq/code/deseq2.r

echo "Generating results into: results_deseq2.txt "
cat simple_counts.txt | Rscript deseq2.r 3x3 > results_deseq2.txt  2>> $RUNLOG

echo "Your results are in the folder called: brain"