|
################################
|
|
#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
|
|
################################
|