Post date: Oct 22, 2019 2:40:32 AM
As part of a Trans. Phil. Soc. special issue paper we plan to revisit the "genomic islands" on LG 8 for T. curi ("cuesta ridge") from the Nature E&E paper. I was able to track down my alignments and variant calling to the melanic (non-flipped, version 3) genome of T. cristinae. They are here:
/uufs/chpc.utah.edu/common/home/gompert-group2/data/timema_sp_wgs/
We have 40 samples from two hosts for curi, given in CuriSamples.txt.
The hosts are C = Ceanothus and CY = Cyprus?
I want to edit/modify the reference to be more curi like. Thus, I am recalling variants with only T. curi and then will create a new consensus sequence, as described in Sam's post.
Here is what I have done.
1. Combine g.vcf tiles. I ran
sbatch runComb.sh
which executes the jobs in cjobs.txt:
java -Xmx80g -jar ~/bin/GenomeAnalysisTK.jar -T CombineGVCFs -R /uufs/chpc.utah.edu/common/home/u6000989/data/timema/tcrDovetail/version3/mod_map_timema_06Jun2016_RvNkF702.fasta --variant uniqe_merge_moe04A1.g.vcf --variant uniqe_merge_moe04A2.g.vcf --variant uniqe_merge_moe04A3.g.vcf --variant uniqe_merge_moe04B1.g.vcf --variant uniqe_merge_moe04B2.g.vcf --variant uniqe_merge_moe04B3.g.vcf --variant uniqe_merge_moe04C1.g.vcf --variant uniqe_merge_moe04C2.g.vcf --variant uniqe_merge_moe04C3.g.vcf --variant uniqe_merge_moe04D1.g.vcf -o tcuri_combinded_0.g.vcf
java -Xmx80g -jar ~/bin/GenomeAnalysisTK.jar -T CombineGVCFs -R /uufs/chpc.utah.edu/common/home/u6000989/data/timema/tcrDovetail/version3/mod_map_timema_06Jun2016_RvNkF702.fasta --variant uniqe_merge_moe04D2.g.vcf --variant uniqe_merge_moe04D3.g.vcf --variant uniqe_merge_moe04E1.g.vcf --variant uniqe_merge_moe04E2.g.vcf --variant uniqe_merge_moe04F1.g.vcf --variant uniqe_merge_moe04F2.g.vcf --variant uniqe_merge_moe04G1.g.vcf --variant uniqe_merge_moe04G2.g.vcf --variant uniqe_merge_moe04H1.g.vcf --variant uniqe_merge_moe04H2.g.vcf -o tcuri_combinded_1.g.vcf
java -Xmx80g -jar ~/bin/GenomeAnalysisTK.jar -T CombineGVCFs -R /uufs/chpc.utah.edu/common/home/u6000989/data/timema/tcrDovetail/version3/mod_map_timema_06Jun2016_RvNkF702.fasta --variant uniqe_merge_moe07A10.g.vcf --variant uniqe_merge_moe07A11.g.vcf --variant uniqe_merge_moe07A9.g.vcf --variant uniqe_merge_moe07B10.g.vcf --variant uniqe_merge_moe07B11.g.vcf --variant uniqe_merge_moe07B9.g.vcf --variant uniqe_merge_moe07C10.g.vcf --variant uniqe_merge_moe07C11.g.vcf --variant uniqe_merge_moe07C9.g.vcf --variant uniqe_merge_moe07D10.g.vcf -o tcuri_combinded_2.g.vcf
java -Xmx80g -jar ~/bin/GenomeAnalysisTK.jar -T CombineGVCFs -R /uufs/chpc.utah.edu/common/home/u6000989/data/timema/tcrDovetail/version3/mod_map_timema_06Jun2016_RvNkF702.fasta --variant uniqe_merge_moe07D9.g.vcf --variant uniqe_merge_moe07E10.g.vcf --variant uniqe_merge_moe07E9.g.vcf --variant uniqe_merge_moe07F10.g.vcf --variant uniqe_merge_moe07F9.g.vcf --variant uniqe_merge_moe07G10.g.vcf --variant uniqe_merge_moe07G9.g.vcf --variant uniqe_merge_moe07H10.g.vcf --variant uniqe_merge_moe07H8.g.vcf --variant uniqe_merge_moe07H9.g.vcf -o tcuri_combinded_3.g.vcf
2. Joint variant calling with GATK The Genome Analysis Toolkit (GATK) v3.5-0-g36282e4. Ran,
sbatch jointVarCallTcur1.sh
Which runs,
#!/bin/sh
#SBATCH --time=210:00:00
#SBATCH --nodes=1
#SBATCH --ntasks=28
#SBATCH --account=gompert-kp
#SBATCH --partition=gompert-kp
#SBATCH --job-name=gatk-geno
#SBATCH --mail-type=FAIL
#SBATCH --mail-user=zach.gompert@usu.edu
module load gatk
cd /uufs/chpc.utah.edu/common/home/gompert-group2/data/timema_sp_wgs/variants
java -Xmx460g -jar ~/bin/GenomeAnalysisTK.jar -T GenotypeGVCFs -R /uufs/chpc.utah.edu/common/home/u6000989/data/timema/tcrDovetail/version3/mod_map_timema_06Jun2016_RvNkF702.fasta -hets 0.001 -nt 32 --variant tcuri_combinded_0.g.vcf --variant tcuri_combinded_1.g.vcf --variant tcuri_combinded_2.g.vcf --variant tcuri_combinded_3.g.vcf -o tcuri_wgs_variants1.vcf
3. Normalize indels and filter. Note that normalization didn't do anything, this must have been accomplished well enough by GATK.
bcftools norm -f /uufs/chpc.utah.edu/common/home/u6000989/data/timema/tcrDovetail/version3/mod_map_timema_06Jun2016_RvNkF702.fasta tcuri_wgs_variants1.vcf -o tcuri_wgs_variants1.norm.vcf
bcftools filter -g 3 -G 10 -e 'DP<5 || QUAL<30 || AF<0.5 || MQ<20' -o filtered_tcuri_wgs_variants1.norm.vcf tcuri_wgs_variants1.norm.vcf
grep -c ^ScR tcuri_wgs_variants1.vcf
39532366
grep -c ^ScR filtered_tcuri_wgs_variants1.vcf
28228629
grep -c ^ScR filtered_tcuri_wgs_variants1.norm.vcf
28228629
3. Create the consensus genome. This is in /uufs/chpc.utah.edu/common/home/gompert-group2/data/timema_sp_wgs/refs/.
bgzip filtered_tcuri_wgs_variants1.vcf
bcftools index filtered_tcuri_wgs_variants1.vcf.gz
bcftools consensus -f ./mod_map_timema_06Jun2016_RvNkF702.fasta filtered_tcuri_wgs_variants1.vcf.gz > tcuri_refgenome.fasta
bwa index -a bwtsw tcuri_refgenome.fasta