shell script template for RNA-seq processing on Odyssey cluster

compbio_hw2_rnaseq_template.sh

#!/bin/bash
#SBATCH -n 1
#SBATCH -t 400
#SBATCH -p serial_requeue
#SBATCH --mem=35000
#SBATCH --mail-type=END,FAIL
#SBATCH --mail-user=my_username@mail.harvard.edu

#take filename from command line argument
FILENAME=$1
#FILENAME = SRR1592185

#load required modules
module load centos6/STAR-2.3.1
module load centos6/samtools-0.1.19
module load centos6/HTSeq-0.6.1_python-2.7.3
module load bio/pysam-0.7.4_python-2.7.3 

#store files in /scratch directory for more disk space and faster I/O
mkdir -p /scratch/${FILENAME}
pushd /scratch/${FILENAME}

#Align reads from FASTQ files to produce SAM files
STAR --genomeDir /n/stat115/db/Mus_musculus/UCSC/mm9/Sequence/STARIndex/ --readFilesIn  /n/stat115/hws/2/${FILENAME}.fastq --runThreadN 8 --outSAMstrandField intronMotif --outFilterScoreMinOverLread 0.5 --outFilterMatchNminOverLread 0.5 --outSAMtype BAM Unsorted
#mv Aligned.out.sam ${FILENAME}.sam
#convert SAM to BAM
#samtools view -bS ${FILENAME}.sam > ${FILENAME}.bam
#rm ${FILENAME}.sam
#no need to sort BAM file according to STAR manual.
#samtools sort ${FILENAME}.bam ${FILENAME}.sort
#count reads from each gene using HTSeq
htseq-count -f bam -r pos -s no Aligned.out.bam /n/stat115/db/Mus_musculus/UCSC/mm9/Annotation/Genes/genes.gtf > ${FILENAME}.count
#move result files back to original directory
popd
mv /scratch/${FILENAME}/${FILENAME}.count ${FILENAME}.count
rm /scratch/${FILENAME}/Aligned.out.bam
mv /scratch/${FILENAME} scratch
#now download the count files to local machine for further analysis with DESeq (bioconductor)