#!/bin/sh

#Data descriptions can be found here:
#https://www.ncbi.nlm.nih.gov/Traces/study/?acc=SRP091330

#copy data dir to desktop and extract
cd ~/Desktop
cp ~/Data/Slch04_demo.tar.gz .
tar -xvf Slch04_demo.tar.gz
rm Slch04_demo.tar.gz
cd Slch04_demo

### 2. Process data files ###
#reads should be checked for quality and contaminants with FastQC (previous class on QC)

#input_file-directory=$1
CPU=1   #this can be changed on multi-core machines

#process reference genome sequence file for building index for hisat2
gffread ITAG3.10_gene_models.gff -o ITAG3.10_gene_models.gtf -T
python ~/Programs/hisat2-2.1.0/extract_splice_sites.py ./ITAG3.10_gene_models.gtf > splicesites.txt
hisat2-build --ss splicesites.txt S_lycopersicum_chromosomes.3.00_ch04.fa S_lycopersicum_chromosomes.3.00_ch04

#map reads with hisat2
for file in `dir -d *_ch4.fastq` ; do

    #create output file name
    samfile=`echo "$file" | sed 's/.fastq/.sam/'`
    metrics=`echo "$file" | sed 's/.fastq/.metrics/'`

    #run mapping with hisat2
    hisat2 --max-intronlen 20000 --dta --met-file $metrics -p $CPU -x S_lycopersicum_chromosomes.3.00_ch04 -U $file -S $samfile
done

#convert sam files to bam files to save space and sort
ls *.sam |parallel -j $CPU samtools view -Sb -o {.}.bam {}
rm *.sam
ls *.bam |parallel -j $CPU samtools sort -o {.}.sort.bam {}
rm *4.bam
ls *.sort.bam |parallel -j $CPU samtools flagstat {} ">" {.}.flagstat
cat *flagstat |grep "mapped (" |sed 's/.*(\(.*\)%.*/\1/g' |awk '{sum+=$1} END { print "Average = ",sum/NR}' > average_mapping.txt

#run stringtie and produce counts table for DE analysis with edgeR or DESeq
for file in `dir -d *.sort.bam` ; do
    
    outfile=`echo "$file" | sed 's/.bam/.gtf/'`
    outdir=`echo "$file" |sed 's/.bam//'`
    stringtie  -e -B -p $CPU -G  ITAG3.10_gene_models.gtf -o ballgown/$outdir/$outfile $file

done

python3 ~/Scripts/prepDE.py -i ballgown -g gene_count_matrix.csv -t transcript_count_matrix.csv 

#import to R with DESeqDataSetFromMatrix (DEseq) and DGEList (EdgeR) in later class
