#!/bin/bash
#PBS -j oe
#PBS -S /bin/bash
#PBS -M lvclark@illinois.edu
#PBS -m abe

##### change info in this header #####
# name of compressed file
CF=sacks.201294.tgz
# name of uncompressed file
FILE=DOE_GBS_lib1_NoIndex_L002_R1_001.fastq
# name of library
LIBNAME=doe_gbs_lib1
######################################

# This script takes a compressed file containing Illumina sequence data and 
# performs several initial steps in the RADseq analysis pipeline.  First, the 
# file is uncompressed.  Second, it is filtered by quality and split by
# barcodes using the process_radtags program in the package Stacks.  Read 
# counts are extracted to a text file, and a pdf is produced containing a 
# histogram of read counts by barcode.  Third, the individual FASTQ files
# produced by process_radtags are aligned to Sorghum using bowtie2, and the
# resulting files are processed by pstacks to detect SNP variants.

# To function, several files need to be in the same directory as the script
# file: the compressed sequence file and four text files containing barcodes.
# Each barcode file contains all barcodes of a particular length, 4 nt, 5 nt,
# 6 nt, or 10 nt.  Each barcode is on a new line.  There additionally needs to 
# be a directory called Sbicolor_index, containing the files produced when 
# bowtie2-build is run on Sbicolor_79_hardmasked.fa.

# Script functions using Stacks version 0.9996.

# --Lindsay V. Clark, 13 Nov 2012.
######################################

# unzip sequences and make directories
tar -zxvf $CF
mkdir $LIBNAME
mkdir bowtie2_$LIBNAME
mkdir pstacks_$LIBNAME

# filter by quality and split by barcode, recording read counts
/home/apps/stacks/bin/process_radtags -f $FILE -o ./$LIBNAME -b barcodes4.txt -e pstI -E phred33 -c -q

sed -n 14,38p $LIBNAME/process_radtags.log > $LIBNAME/readcounts.txt

/home/apps/stacks/bin/process_radtags -f $FILE -o ./$LIBNAME -b barcodes5.txt -e pstI -E phred33 -c -q

sed -n 15,38p $LIBNAME/process_radtags.log >> $LIBNAME/readcounts.txt

/home/apps/stacks/bin/process_radtags -f $FILE -o ./$LIBNAME -b barcodes6.txt -e pstI -E phred33 -c -q

sed -n 15,38p $LIBNAME/process_radtags.log >> $LIBNAME/readcounts.txt

/home/apps/stacks/bin/process_radtags -f $FILE -o ./$LIBNAME -b barcodes10.txt -e pstI -E phred33 -c -q

sed -n 15,38p $LIBNAME/process_radtags.log >> $LIBNAME/readcounts.txt

# make a histogram of the read count distribution
cd $LIBNAME
/home/apps/R/bin/R --vanilla <<RSCRIPT
mycounts <- read.table("readcounts.txt", header=TRUE, stringsAsFactors=FALSE,
                       sep="\t")
pdf("hist.pdf",5,5)
hist(mycounts[,"Retained"], breaks=20, main="Read count distribution")
dev.off()
q()
RSCRIPT

# copy the read counts and histogram over and give library-specific names
cp readcounts.txt ../readcounts_$LIBNAME.txt
cp hist.pdf ../hist_$LIBNAME.pdf

# loop through files to do alignment and pstacks
for i in `ls -1 *.fq`
do
/home/apps/bowtie2/bowtie2 -x ../Sbicolor_index/Sbicolor_79_hardmasked -U $i -S ../bowtie2_$LIBNAME/$i.sam
/home/apps/stacks/bin/pstacks -t sam -f ../bowtie2_$LIBNAME/$i.sam -o ../pstacks_$LIBNAME
done