#!/bin/bash

#SBATCH --array=319-384  # NEED TO CHANGE THIS!
#SBATCH --job-name=STAR_gencode # Job name
#SBATCH --nodes=1
#SBATCH --ntasks=8 # Number of cores
#SBATCH --time=3:00:00
#SBATCH --mem=48000 # Memory pool for all cores (see also --mem-per-cpu)
#SBATCH --partition=build # Partition to submit to
#SBATCH --output=arrayJob_%A_%a.out # File to which STDOUT will be written
#SBATCH --error=arrayJob_%A_%a.err # File to which STDERR will be written


start=`date +%s`

echo $HOSTNAME

echo "My SLURM_ARRAY_TASK_ID: " $SLURM_ARRAY_TASK_ID

THREADS=8

sample=`sed "${SLURM_ARRAY_TASK_ID}q;d" array_samples.txt`

echo $sample

# source /etc/profile.d/modules_sh.sh

module load star
module load samtools

which STAR >> software_versions.txt
which samtools >> software_versions.txt


# Remember, if fastq files are .gz, need to add
# --readFilesCommand zcat \

mkdir star_counts

inpath='/share/biocore/projects/Gershwin_L_UCD/Gershwin-Lebedev_Bos_taurus_Batch_TagSeq/01-HTS_Preproc/'
mappath='/share/biocore/projects/Gershwin_L_UCD/Gershwin-Lebedev_Bos_taurus_Batch_TagSeq/02-star/'
genomepath='/share/genomes/ensembl/Bos_taurus.ARS-UCD1.2/star_v2.7_index100'


mkdir ${mappath}${sample}

rm -rf ${mappath}${sample}/*

which STAR >> ${mappath}${sample}/software_versions.txt
which samtools >> ${mappath}${sample}/software_versions.txt


STAR --runThreadN 8 \
--outFilterType BySJout \
--outFilterMultimapNmax 20 \
--alignSJoverhangMin 8 \
--alignSJDBoverhangMin 1 \
--outFilterMismatchNmax 999 \
--alignIntronMin 20 \
--alignIntronMax 1000000 \
--alignMatesGapMax 1000000 \
--chimSegmentMin 20 \
--genomeDir ${genomepath} \
--outFileNamePrefix ${mappath}${sample}/${sample}. \
--outSAMtype BAM SortedByCoordinate \
--outSAMunmapped Within \
--outReadsUnmapped Fastx \
--quantMode GeneCounts \
--readFilesCommand zcat \
--readFilesIn \
${inpath}${sample}/${sample}_SE.fastq.gz \
> ${mappath}${sample}/STAR.stdout 2> ${mappath}${sample}/STAR.stderr

wait

samtools index ${mappath}${sample}/${sample}.Aligned.sortedByCoord.out.bam

echo "id ${sample}" | sed 's/ /\t/g' > ${mappath}star_counts/${sample}_counts.txt
cut -f1,3 ${mappath}${sample}/${sample}.ReadsPerGene.out.tab \
>> ${mappath}star_counts/${sample}_counts.txt

samtools flagstat ${mappath}${sample}/${sample}.Aligned.sortedByCoord.out.bam \
> ${mappath}${sample}/${sample}_flagstat.txt