if [ "$#" -lt 12 ]
then
    echo "RUN_ONE_SAMPLE_STEP1.sh [sample_name] <mate1_fastq> <mate2_fastq> [adapter_sequence_mate1] [adapter_sequence_mate2] [read_length] [minimum_fragment_size] [num_threads_per_sample] <STAR_index> <gene_model_gtf_file> <chromosome_sizes_file> [normalization_factor]"
    exit 1
fi

S=$1
shift
MATE1_FASTQ=$1
shift
MATE2_FASTQ=$1
shift
ADAPTER_SEQUENCE_MATE1=$1
shift
ADAPTER_SEQUENCE_MATE2=$1
shift
READ_LENGTH=$1
shift
MINIMUM_FRAGMENT_SIZE=$1
shift
NUM_THREADS=$1
shift
STAR_INDEX=$1
shift
GENE_MODEL=$1
shift
CHROM_SIZES=$1
shift
NORMALIZATION_FACTOR=$1
shift

mkdir -p out/$S

echo ***START*** `date`

### Filter reads ###
/ddn/gs1/home/bennettb/bin/trim_and_filter_fastq --paired -i1 $MATE1_FASTQ -i2 $MATE2_FASTQ -q 20 -o1 out/$S/filtered_reads.1.fastq -o2 out/$S/filtered_reads.2.fastq -r out/$S/filtered_reads.FilterStats.txt

### Remove adapter ###
#cutadapt -f fastq --match-read-wildcards -a $ADAPTER_SEQUENCE_MATE1 out/$S/filtered_reads.1.fastq -o out/$S/trimmed_reads.1.fastq > out/$S/adapter_trimming_report.1.txt
#cutadapt -f fastq --match-read-wildcards -a $ADAPTER_SEQUENCE_MATE2 out/$S/filtered_reads.2.fastq -o out/$S/trimmed_reads.2.fastq > out/$S/adapter_trimming_report.2.txt
#perl fix_adapter_trimming.pl out/$S/trimmed_reads.1.fastq out/$S/trimmed_reads.2.fastq $READ_LENGTH $MINIMUM_FRAGMENT_SIZE out/$S/final_reads.1.fastq out/$S/final_reads.2.fastq > out/$S/adapter_report.txt
mv out/$S/filtered_reads.1.fastq out/$S/final_reads.1.fastq
mv out/$S/filtered_reads.2.fastq out/$S/final_reads.2.fastq

### Alignment ###
mkdir -p out/$S/STAR_out
STAR --runThreadN $NUM_THREADS --genomeDir $STAR_INDEX --readFilesIn out/$S/final_reads.1.fastq out/$S/final_reads.2.fastq --outFileNamePrefix out/$S/STAR_out/
perl count_aligned_fragments.pl out/$S/STAR_out/Aligned.out.sam out/$S/number_aligned_fragments.txt
samtools view -Sb out/$S/STAR_out/Aligned.out.sam > out/$S/STAR_out/Aligned.out.bam
rm -f out/$S/STAR_out/Aligned.out.sam

### Get gene read counts ###
/ddn/gs1/home/bennettb/tools/subread-1.5.1-Linux-x86_64/bin/featureCounts -p -T $NUM_THREADS -a $GENE_MODEL -g gene_name -o out/$S/featureCounts_out.txt out/$S/STAR_out/Aligned.out.bam
perl format_counts.pl out/$S/featureCounts_out.txt out/$S/gene_read_counts.txt

### Make browser tracks ###
samtools sort -o out/$S/STAR_out/Aligned.out_sorted.bam out/$S/STAR_out/Aligned.out.bam
/ddn/gs1/home/bennettb/tools/bedtools-2.17.0/bin/genomeCoverageBed -bg -split -g $CHROM_SIZES -ibam out/$S/STAR_out/Aligned.out_sorted.bam > out/$S/coverage.bedGraph
perl normalize_coverage.pl out/$S/coverage.bedGraph $NORMALIZATION_FACTOR `cat out/$S/number_aligned_fragments.txt` out/$S/coverage_normalized.bedGraph
/ddn/gs1/home/bennettb/tools/ucsc/bedGraphToBigWig out/$S/coverage_normalized.bedGraph $CHROM_SIZES out/$S/coverage_normalized.bigWig
rm -f out/$S/STAR_out/Aligned.out_sorted.bam

### Clean up ###
rm -f out/$S/filtered_reads.1.fastq
rm -f out/$S/filtered_reads.2.fastq
rm -f out/$S/trimmed_reads.1.fastq
rm -f out/$S/trimmed_reads.2.fastq
rm -f out/$S/final_reads.1.fastq
rm -f out/$S/final_reads.2.fastq

echo ***DONE*** `date`