#!/usr/bin/env bash

# Run featureCounts on the results of the alignments. We assume
# That you ran zika-getdata.sh then zika-align.sh
set -ueo pipefail

# Collect the output of commands here.
RUNLOG=runlog.txt
echo "# See the $RUNLOG file for run time messages"

# These are the MOCK numbers
MOCK=(SRR3191542 SRR3191543 SRR3194428 SRR3194429)

# These are the ZIKV numbers.
ZIKV=(SRR3191544 SRR3191545 SRR3194430 SRR3194431)

for SRR in ${MOCK[@]}; do
    MOCK_BAMS+="bam/${SRR}.bam "
done

for SRR in ${ZIKV[@]}; do
    ZIKV_BAMS+="bam/${SRR}.bam "
done

# Run the featurecounts program.
GTF=refs/GRCh38.gtf
featureCounts 2>> $RUNLOG -p -g gene_name -a $GTF -o counts.txt $MOCK_BAMS $ZIKV_BAMS

# Simplify the counts.
cat counts.txt | cut -f 1,7-14 > sample_counts.txt

#
# We have commented out the scripts here but you can re-enable them
# if you want to get a DE results right after counting.
#
# The experimental design is pairwise comparisons 4x4.
#
# echo "*** Running DESeq1."
# curl -sO http://data.biostarhandbook.com/rnaseq/code/deseq1.r
# cat sample_counts.txt | Rscript deseq1.r 4x4 2>> $RUNLOG > results_deseq1.txt

# echo "*** Running DESeq2."
# curl -sO http://data.biostarhandbook.com/rnaseq/code/deseq2.r
# cat sample_counts.txt | Rscript deseq2.r 4x4 2>> $RUNLOG > results_deseq2.txt

# echo "*** Running EdgeR."
# curl -sO http://data.biostarhandbook.com/rnaseq/code/edger.r
# cat sample_counts.txt | Rscript edger.r 4x4 2>> $RUNLOG > results_edger.txt