Post date: Oct 18, 2016 9:52:6 PM
My previous alignments with mummer are nice in that they are all against all, but they end up with pretty sparse alignments. I am considering more population genomic based tests of introgression (i.e., sliding window Fst HMM stuff, or something like that) and it would be nice to have all of the genomes aligned to the annotated L. melissa genome for this. I am using bwa mem for this. The alignments will be in /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/bwa_alignments/:
1. Alignment with bwa mem, submitted via perl ../scripts/wrap_qsub_slurm_bwa-wgs.pl *180*1.fastq.gz
cd /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/bwa_alignments/
bwa mem -t 8 -k 20 -w 100 -r 1.3 -T 30 -R '@RG\tID:anna\tPL:ILLUMINA\tLB:anna\tSM:anna' /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/melissa_genome/final.assembly.fasta /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/sequences/anna_180bp_1.fastq.gz /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/sequences/anna_180bp_2.fastq.gz > /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/bwa_alignments/aln_anna.sam 2> /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/bwa_alignments/error_anna.log
cd /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/bwa_alignments/
bwa mem -t 8 -k 20 -w 100 -r 1.3 -T 30 -R '@RG\tID:idas\tPL:ILLUMINA\tLB:idas\tSM:idas' /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/melissa_genome/final.assembly.fasta /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/sequences/idas_180bp_1.fastq.gz /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/sequences/idas_180bp_2.fastq.gz > /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/bwa_alignments/aln_idas.sam 2> /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/bwa_alignments/error_idas.log
cd /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/bwa_alignments/
bwa mem -t 8 -k 20 -w 100 -r 1.3 -T 30 -R '@RG\tID:melissa\tPL:ILLUMINA\tLB:melissa\tSM:melissa' /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/melissa_genome/final.assembly.fasta /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/sequences/melissa_180bp_1.fastq.gz /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/sequences/melissa_180bp_2.fastq.gz > /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/bwa_alignments/aln_melissa.sam 2> /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/bwa_alignments/error_melissa.log
cd /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/bwa_alignments/
bwa mem -t 8 -k 20 -w 100 -r 1.3 -T 30 -R '@RG\tID:sierra\tPL:ILLUMINA\tLB:sierra\tSM:sierra' /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/melissa_genome/final.assembly.fasta /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/sequences/sierra_180bp_1.fastq.gz /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/sequences/sierra_180bp_2.fastq.gz > /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/bwa_alignments/aln_sierra.sam 2> /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/bwa_alignments/error_sierra.log
cd /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/bwa_alignments/
bwa mem -t 8 -k 20 -w 100 -r 1.3 -T 30 -R '@RG\tID:warner\tPL:ILLUMINA\tLB:warner\tSM:warner' /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/melissa_genome/final.assembly.fasta /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/sequences/warner_180bp_1.fastq.gz /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/sequences/warner_180bp_2.fastq.gz > /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/bwa_alignments/aln_warner.sam 2> /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/bwa_alignments/error_warner.log
The alignments look pretty good at a first pass, for L. anna (presumably most divergent) 77.7% of reads aligned.
2. Compress, sort and index alignments with samtools.
perl ../scripts/wrap_qsub_slurm_sam2bam.pl aln_*sam
cd /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/bwa_alignments/
samtools view -b -S -o aln_anna.bam aln_anna.sam
samtools sort aln_anna.bam aln_anna.sorted
samtools index aln_anna.sorted.bam
3. Use picard tools to remove PCR duplicates.
perl markDups.pl *sorted.bam
An example.
java -Xmx100g -jar /uufs/chpc.utah.edu/sys/installdir/picard/2.1.1/picard.jar MarkDuplicates INPUT=aln_anna.sorted.bam OUTPUT=aln_anna.dedup.bam METRICS_FILE=aln_anna.dedup.log.txt
Re-indexed with samtools.
4. Create a dict file for the reference genome (required for GATK) a
module load picard
java -jar /uufs/chpc.utah.edu/sys/installdir/picard/2.1.1/picard.jar CreateSequenceDictionary R=/uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/melissa_genome/final.assembly.fasta O=/uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/melissa_genome/final.assembly.dict
5. Initial component of variant calling where the HaplotypeCaller is run on each sample independently to produce an intermediate *.g.vcf file (this includes genotype likelihoods). This will then be combined for true variant calling afterwards. I submitted the job to my node, using the multi-prog command.
Specifically, I ran this: sbatch subGatk.sh from here /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/bwa_alignments/
#!/bin/sh
#SBATCH --time=240:00:00
#SBATCH --nodes=1
#SBATCH --ntasks=5
#SBATCH --account=gompert-kp
#SBATCH --partition=gompert-kp
#SBATCH --job-name=gatk-hc
#SBATCH --mail-type=FAIL
#SBATCH --mail-user=zach.gompert@usu.edu
srun --multi-prog /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/bwa_alignments/gatk.conf
Which executes gatk.conf:
0 gatk0.sh
1 gatk1.sh
2 gatk2.sh
3 gatk3.sh
4 gatk4.sh
And runs these five commands (each from a gatk*.sh file):
==> gatk0.sh <==
#!/bin/sh
java -Xmx96g -jar ~/bin/GenomeAnalysisTK.jar -T HaplotypeCaller -R /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/melissa_genome/final.assembly.fasta -I /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/bwa_alignments/aln_anna.dedup.bam -o /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/variants/aln_anna.dedup.g.vcf -gt_mode DISCOVERY -hets 0.001 -mbq 30 -out_mode EMIT_VARIANTS_ONLY -ploidy 2 -stand_call_conf 50 -pcrModel AGGRESSIVE --emitRefConfidence GVCF --variant_index_type LINEAR --variant_index_parameter 128000
==> gatk1.sh <==
#!/bin/sh
java -Xmx96g -jar ~/bin/GenomeAnalysisTK.jar -T HaplotypeCaller -R /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/melissa_genome/final.assembly.fasta -I /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/bwa_alignments/aln_idas.dedup.bam -o /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/variants/aln_idas.dedup.g.vcf -gt_mode DISCOVERY -hets 0.001 -mbq 30 -out_mode EMIT_VARIANTS_ONLY -ploidy 2 -stand_call_conf 50 -pcrModel AGGRESSIVE --emitRefConfidence GVCF --variant_index_type LINEAR --variant_index_parameter 128000
==> gatk2.sh <==
#!/bin/sh
java -Xmx96g -jar ~/bin/GenomeAnalysisTK.jar -T HaplotypeCaller -R /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/melissa_genome/final.assembly.fasta -I /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/bwa_alignments/aln_melissa.dedup.bam -o /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/variants/aln_melissa.dedup.g.vcf -gt_mode DISCOVERY -hets 0.001 -mbq 30 -out_mode EMIT_VARIANTS_ONLY -ploidy 2 -stand_call_conf 50 -pcrModel AGGRESSIVE --emitRefConfidence GVCF --variant_index_type LINEAR --variant_index_parameter 128000
==> gatk3.sh <==
#!/bin/sh
java -Xmx96g -jar ~/bin/GenomeAnalysisTK.jar -T HaplotypeCaller -R /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/melissa_genome/final.assembly.fasta -I /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/bwa_alignments/aln_sierra.dedup.bam -o /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/variants/aln_sierra.dedup.g.vcf -gt_mode DISCOVERY -hets 0.001 -mbq 30 -out_mode EMIT_VARIANTS_ONLY -ploidy 2 -stand_call_conf 50 -pcrModel AGGRESSIVE --emitRefConfidence GVCF --variant_index_type LINEAR --variant_index_parameter 128000
==> gatk4.sh <==
#!/bin/sh
java -Xmx96g -jar ~/bin/GenomeAnalysisTK.jar -T HaplotypeCaller -R /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/melissa_genome/final.assembly.fasta -I /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/bwa_alignments/aln_warner.dedup.bam -o /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/variants/aln_warner.dedup.g.vcf -gt_mode DISCOVERY -hets 0.001 -mbq 30 -out_mode EMIT_VARIANTS_ONLY -ploidy 2 -stand_call_conf 50 -pcrModel AGGRESSIVE --emitRefConfidence GVCF --variant_index_type LINEAR --variant_index_parameter 128000
6. Joint variant calling with GATK.
sbatch genomeVarCall.sh
#!/bin/sh
#SBATCH --time=240:00:00
#SBATCH --nodes=1
#SBATCH --ntasks=24
#SBATCH --account=gompert-kp
#SBATCH --partition=gompert-kp
#SBATCH --job-name=gatk
#SBATCH --mail-type=FAIL
#SBATCH --mail-user=zach.gompert@usu.edu
cd /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/whole_genomes/variants
java -Xmx460g -jar ~/bin/GenomeAnalysisTK.jar -T GenotypeGVCFs -R /uufs/chpc.utah.edu/common/home/u6000989/data/lycaeides/melissa_genome/final.assembly.fasta --variant aln_anna.dedup.g.vcf --variant aln_melissa.dedup.g.vcf --variant aln_warner.dedup.g.vcf --variant aln_idas.dedup.g.vcf --variant aln_sierra.dedup.g.vcf -ploidy 2 -o LycWgVars.vcf