#!/bin/bash                                                                                                                                                   

# The goal of this script is to go from a fastq file to a count table of reads mapped to each orthologous gene. 
# Note: There are several places where users will need to put in the appropriate path for programs (e.g. tophat2, bowtie2-2.2.1, samtools, featureCounts) or files (reference sequence and annotation files). The reference sequence and annotation files for hg19 and panTro3 can be found at https://support.illumina.com/sequencing/sequencing_software/igenome.html. The orthologous exon files for hg19 and panTro3 are included as part of this dataset. 

# Step 1: Trim the read in the fastq.gz file

python fqTrim1.py  1 20 30 <(zcat ./*.fastq.gz) >(gzip -c > ./first_trim_py.gz)

echo "File has been trimmed for quality"

# Unzip the file

gunzip ./first_trim_py.gz

echo "File has been gun-zipped"

# Map reads 

export PATH="$PATH:/path_to_bowtie2-2.2.1

/path_to_tophat2/tophat2 -p 4 -G /path_to_PanTro3_or_hg19_genes/genes.gtf -m 2 -g 1\
 --library-type fr-firststrand --b2-very-sensitive -o tophat2_results --no-novel-juncs /path_to_panTro3_or_hg19_Bowtie2Index/genome ./first_trim_py > ./read_trimmed.fq.gz

echo "Finished with tophat"

# Count the number of reads that map to exons/genes in the appropriate species’ genome

/path_to_samtools/samtools view ./tophat2_results/accepted_hits.bam | /mnt/gluster/data/tools/anaconda2/bin/htseq-count -m intersection-nonempty -r pos -s\
 reverse -t exon - /path_to_panTro3_or_hg19_genes_file/genes.gtf > ./tophat2_results/htseq-count.txt
                                                       

echo "Finished counting mapped reads"

# Count the number of reads that map to orthologous exons/genes

/path_to_featureCounts/featureCounts  -p -B -t exon -g gene_id -a /path_to_panTro3_or_hg19_orthologous_exons_file/exon.ortho.correct_species.gff -o ./GeneCounts.npairs.txt ./tophat2_results/accepted_hits.bam

echo "Finished mapping to ortho exons"