################################
#This pipeline is for single sample variants detection including
#Single Nucleotide Variants and Short Indels from 
#Whole Genome Sequencing (WGS) or Whole Exome Sequencing (WES)
#with a minimum 10x of average depth of coverage.
#
#Sequencing platform: illumina HiSeq/MiSeq
#Pair-end reads with length from 50bp to 150bp.
################################

#env settings, Java Runtime required. (openjdk-7 is used here) 
CUROVERSE_HOME=/mnt/curoverse
APP=$CUROVERSE_HOME/app
DATA=$CUROVERSE_HOME/data

################################
#prepare applications
################################

#1. BWA
wget http://downloads.sourceforge.net/project/bio-bwa/bwa-0.7.9a.tar.bz2
tar -jxvf bwa-0.7.9a.tar.bz2
cd bwa-0.7.9a && make && cd $CUROVERSE_HOME
mkdir -p $APP/bwa/0.7.9a/
find bwa-0.7.9a -executable -type f -print0 | xargs -0 -I {} cp {} $APP/bwa/0.7.9a/
rm -fr bwa* 

#2. samtools
wget http://downloads.sourceforge.net/project/samtools/samtools/0.1.19/samtools-0.1.19.tar.bz2
tar -jxvf samtools-0.1.19.tar.bz2
cd samtools-0.1.19 && make && cd $CUROVERSE_HOME
mkdir -p $APP/samtools/0.1.19/
find  samtools-0.1.19 -executable -type f -print0 | xargs -0 -I {} mv {} $APP/samtools/0.1.19/
rm -fr samtools*

#3. picard
wget http://downloads.sourceforge.net/project/picard/picard-tools/1.114/picard-tools-1.114.zip
unzip picard-tools-1.114.zip
mkdir -p $APP/picard/1.114/
mv picard-tools-1.114/* $APP/picard/1.114/
rm -fr picard*


#4. GATK
#no automatic way to install GATK.
Register and download the GATK package. unpack it to /mnt/curoverse/app/gatk/3.1-1/


################################
#prepare reference data
################################

mkdir -p $DATA/fasta/hg19/ && cd $DATA/fasta/hg19/

for i in $(seq 1 22) X Y M; do echo $i; wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr${i}.fa.gz; done
gunzip *.gz
for i in $(seq 1 22) X Y M; do cat chr${i}.fa >> hg19.fasta; done
rm -fr chr*.fasta

#create index
$APP/samtools/0.1.19/samtools faidx hg19.fasta
#create dictionary
java -Xmx8g -jar $APP/picard/0.1.19/CreateSequenceDictionary.jar R=hg19.fasta O=hg19.dict

java -Xmx8g -jar /home/hadoop/curoverse/app/picard/1.114/CreateSequenceDictionary.jar R=hg19.fasta O=hg19.dict

cd $CUROVERSE_HOME

################################
#prepare library files (dbSNP, 1000G, etc)
################################
mkdir -p $DATA/lib/hg19/ && cd $DATA/lib/hg19/

#dbSNP138
#wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/snp138Common.txt.gz

#GATK resource bundle
GATK_FTP=ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/2.8/hg19

#download library file for GATK's VariantRecalibrator
for f in hapmap_3.3.hg19.vcf.gz 1000G_omni2.5.hg19.vcf.gz 1000G_phase1.indels.hg19.vcf.gz Mills_and_1000G_gold_standard.indels.hg19.vcf.gz
do
curl --remote-name ${GATK_FTP}/${f}
done

################################
#Demo sample data (NA12878 from 1000G)
#NA12878: pair-end, WGS, 50x depth, Illumina HiSeq 2000 system  
#URL: http://www.ebi.ac.uk/ena/data/view/ERX237515
#NA12878_R1: ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR262/ERR262997/ERR262997_1.fastq.gz
#NA12878_R2: ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR262/ERR262997/ERR262997_2.fastq.gz
################################

#~100GB disk space required to download the NA12878
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR262/ERR262997/ERR262997_1.fastq.gz
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR262/ERR262997/ERR262997_2.fastq.gz

################################
#Begin of pipeline

#For testing purpose, two PE fastq files R1 and R2 were generated by pIRS (http://code.google.com/p/pirs/)
#and only chr1 of hg19 was used as reference genome

cd /mnt/

#0.1 create simulated PE Reads on 50x coverage
pirs simulate -i $DATA/hg19.chr1/hg19.chr1.fa \
-s ~/app/pIRS_111/Profiles/Base-Calling_Profiles/humNew.PE100.matrix.gz \
-d ~/app/pIRS_111/Profiles/GC-depth_Profiles/humNew.gcdep_200.dat \
-b ~/app/pIRS_111/Profiles/InDel_Profiles/phixv2.InDel.matrix \
-l 100 -x 50 -Q 33 -o A

R1=/mnt/A_100_500_1.fq.gz
R2=/mnt/A_100_500_2.fq.gz

#0.2 build genome index for BWA
$APP/app/bwa/0.7.9a/bwa index -p hg19.chr1 chr1.fa

#1. align with BWA, with 4 threads (-t 4). Reading Group is required (-R ...).
time $APP/bwa/0.7.9a/bwa mem -t 4 -R '@RG\tID:group_id\tPL:illumina\tSM:sample_id' $DATA/hg19.chr1/hg19.chr1 ${R1} ${R2} > A.chr1.sam

#2. Sort into BAM
java -Xmx4g -Djava.io.tmpdir=/tmp -jar $APP/picard/1.114/SortSam.jar \
CREATE_INDEX=True SORT_ORDER=coordinate VALIDATION_STRINGENCY=LENIENT \
INPUT=A.chr1.sam OUTPUT=A.chr1.sort.bam

#3. remove PCR duplicate
java -Xmx4g -Djava.io.tmpdir=/tmp -jar $APP/picard/1.114/MarkDuplicates.jar \
CREATE_INDEX=true REMOVE_DUPLICATES=True ASSUME_SORTED=True VALIDATION_STRINGENCY=LENIENT METRICS_FILE=/dev/null \
INPUT=A.chr1.sort.bam OUTPUT=A.chr1.sort.dedup.bam

#4. fix mate information
java -Djava.io.tmpdir=/tmp/ -jar $APP/picard/1.114/FixMateInformation.jar \
SO=coordinate VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true \
INPUT=A.chr1.sort.dedup.bam OUTPUT=A.chr1.sort.dedup.mate.bam

#5. local alignment around Indels. create de novo intervals 
java -Xmx4g -jar $APP/gatk/3.1-1/GenomeAnalysisTK.jar -R $DATA/hg19.chr1/chr1.fa \
-T RealignerTargetCreator \
-I A.chr1.sort.dedup.mate.bam -o A.intervals

java -Xmx4g -jar $APP/gatk/3.1-1/GenomeAnalysisTK.jar -R $DATA/hg19.chr1/chr1.fa \
-T IndelRealigner \
-targetIntervals A.intervals \
-I A.chr1.sort.dedup.mate.bam -o A.chr1.sort.dedup.mate.relaign.bam

#6. base recalibrator, use dbSNP139 as "truth". If we know the population of subject sample, we can 
#use specific SNP datasets extract from 1000G with same population.
java -Xmx4g -jar $APP/gatk/3.1-1/GenomeAnalysisTK.jar -R $DATA/hg19.chr1/chr1.fa \
-T BaseRecalibrator \
-knownSites $DATA/dbsnp_138.hg19.vcf \
-o A.recal \
-I A.chr1.sort.dedup.mate.relaign.bam

java -Xmx4g -jar $APP/gatk/3.1-1/GenomeAnalysisTK.jar -R $DATA/hg19.chr1/chr1.fa \
-T PrintReads \
-BQSR A.recal \
-o A.chr1.sort.dedup.mate.relaign.recal.bam \
-I A.chr1.sort.dedup.mate.relaign.bam

#7. call variants with GATK's UnifiedGenotyper (faster, lots of raw calls, must run filter after)
time java -Xmx4g -jar $APP/gatk/3.1-1/GenomeAnalysisTK.jar -R $DATA/hg19.chr1/chr1.fa \
-T UnifiedGenotyper \
-glm BOTH \
-D $DATA/dbsnp_138.hg19.vcf \
-metrics A.snps.metrics \
-stand_call_conf 50.0 -stand_emit_conf 10.0 \
-o A.ug.vcf \
-I A.chr1.sort.dedup.mate.relaign.recal.bam

#8. call variants with GATK's HaplotypeCaller (2x slower than UG, much less accurate calls)
time java -Xmx4g -jar $APP/gatk/3.1-1/GenomeAnalysisTK.jar -R $DATA/hg19.chr1/chr1.fa \
-T HaplotypeCaller \
--dbsnp $DATA/dbsnp_138.hg19.vcf \
-stand_call_conf 50.0 -stand_emit_conf 10.0 \
-o A.hc.vcf \
-I A.chr1.sort.dedup.mate.relaign.recal.bam

#End of pipeline
################################