#!/usr/bin/env bash


#https://www.biostars.org/p/91020/

usage='
Use star to align fastq files against a genome
Usage: star_align.sh <igenome> <fastq_files>...

Options:
'
#-c        Cache results

eval $(echo "$usage" | ~/bin/docopts/docopts -h - -A dopts : "$@")
#eval $(echo "$usage" | ~/bin/docopts/docopts -h -  : /projects/bioinfo/igenomes/Mus_musculus/Ensembl/GRCm38 /projects/bioinfo/holger/projects/florio_11b_2nd_batch/lanereps_pooled/arhgap11b_1.fastq.gz)

#for fastqFile in ${fastq_files[@]} ; do
#    echo processing $fastqFile
#done
#
#echo $igenome


## build index if not present



export star_index="$igenome/Sequence/StarIndex/genome"
export gtfFile="$igenome/Annotation/Genes/genes.gtf"
head $gtfFile

if [ ! -f $gtfFile ]; then
    >&2 echo "gtf '$gtfFile' does not exis"; exit 1;
fi

## basic 
#http://www.homolog.us/blogs/blog/2012/11/02/star-really-kick-ass-rna-seq-aligner/

if ! -d "${igenome}/Sequence/StarIndex" ]; then

mailme "${project}: creating STAR index for $igenome"
mkdir ${igenome}/Sequence/StarIndex

cmd="STAR --runMode genomeGenerate --genomeDir ${igenome}/Sequence/StarIndex --genomeFastaFiles ${igenome}/Sequence/WholeGenomeFasta/genome.fa --runThreadN 10"
#eval $cmd
#STAR --runMode genomeGenerate --genomeDir ${igenome}/Sequence/StarIndex --genomeFastaFiles ${igenome}/Sequence/Chromosomes/*.fa --runThreadN 10
mysub "${project}_star_index" "$cmd" -n 5 -R span[hosts=1] -q medium  | blockScript

## prevent modification
chmod -R -w ${igenome}/Sequence/StarIndex

fi



if [ ! -f $gtfFile ]; then
    >&2 echo "gtf '$gtfFile' does not exis"; exit 1;
fi

if [ -z "$(which STAR)" ]; then
    >&2 echo "no tomcat binary in PATH"; exit 1;
fi


echo "running STAR using igenome '${igenome}' for the following files"
ll $fastqFiles

#fastqFiles=$(ls $baseDir/treps_pooled/*fastq.gz)

for fastqFile in $fastqFiles ; do
    echo "submitting tophat job for $fastqFile"

    # DEBUG fastqFile=/projects/bioinfo/holger/projects/helin/mouse/trimmed/mouse_big_cyst_rep4_ca.fastq.gz
    fastqBaseName=$(basename ${fastqFile%%.fastq.gz})
    outputdir=$fastqBaseName

    ## uniquely mapping reads only:   http:/seqanswers.com/forums/showthread.php?s=93da11109ae12a773a128a637d69defe&t=3464
    ###     tophat -p6  -G $gtfFile -g1 -o test $bowtie_gindex $fastqFile

    ## note --outSAMstrandField intronMotif is required for cuffdiff compatiblity (xs flag)
    mysub "${project}__star__${fastqBaseName}" "
#    tophat -p6  -G $gtfFile -g1 -o $outputdir $bowtie_gindex $fastqFile
    STAR --genomeDir $star_index --readFilesIn /path/to/read1 $fastqFile --runThreadN 6 --outFileNamePrefix $fastqBaseName --outSAMtype BAM SortedByCoordinate --outSAMstrandField intronMotif --sjdbGTFfile
    " -n 5 -R span[hosts=1] -q long | joblist .tophatjobs
done

wait4jobs .tophatjobs



dge_bam_correlate .

## create tophat mapping report
## todo adjust report to STAR
#spin.R ${NGS_TOOLS}/dge_workflow/tophat_qc.R .

mailme "$project: STAR done in $(pwd)"