Post date: Aug 20, 2018 11:11:27 PM
I ran the alignment and variant calling from these for the field. I kept running into issues in part from using different versions of the genome. To be safe, I plan to run it all again.
Here is the original description.
I did not regenerate the IDs. This was step 1.
2. Run the alignments with bwa mem (version 0.7.17-r1188). Note that there are four runs for each sample (with two files each).
From /uufs/chpc.utah.edu/common/home/u6000989/data/timema/combind_wgs_dovetailV3/alignments_tcr/
perl ../scripts/wrap_qsub_slurm_bwa-wgs-wild.pl
Here is an example command (there are 640 jobs).
cd /uufs/chpc.utah.edu/common/home/u6000989/data/timema/combind_wgs_dovetailV3/alignments_tcr/
bwa mem -t 8 -k 20 -w 100 -r 1.3 -T 30 -R '@RG\tID:timemaR12C_10221021-69603\tPL:ILLUMINA\tLB:timemaR12C_10221021\tSM:timemaR12C_10221021' /uufs/chpc.utah.edu/common/home/u6000989/data/timema/tcrDovetail/version3/mod_map_timema_06Jun2016_RvNkF702.fasta /uufs/chpc.utah.edu/common/home/u6000989/data/timema/timema_wgrs/plate6/WTCHG_69603_272_1.fastq /uufs/chpc.utah.edu/common/home/u6000989/data/timema/timema_wgrs/plate6/WTCHG_69603_272_2.fastq > /uufs/chpc.utah.edu/common/home/u6000989/data/timema/combind_wgs_dovetailV3/alignments_tcr/aln_6_69603_timemaR12C_10221021.sam 2> /uufs/chpc.utah.edu/common/home/u6000989/data/timema/combind_wgs_dovetailV3/alignments_tcr/error_6_69603_timemaR12C_10221021.log
3. Process alignments.
First, I sorted, compressed and indexed the 640 alignment *sam files with samtools version 1.5 (htslib 1.5).
sbatch SubSam2Bam.sh
which runs Sam2BamFork.pl,
#!/usr/bin/perl
#
# conver sam to bam, then sort and index
#
use Parallel::ForkManager;
my $max = 24;
my $pm = Parallel::ForkManager->new($max);
FILES:
foreach $sam (@ARGV){
$pm->start and next FILES; ## fork
$sam =~ m/^([A-Za-z0-9_\-]+\.)sam/ or die "failed to match $fq\n";
$base = $1;
system "samtools view -b -O BAM -o $base"."bam $sam\n";
system "samtools sort -O BAM -o $base"."sorted.bam $base"."bam\n";
system "samtools index -b $base"."sorted.bam\n";
$pm->finish;
}
$pm->wait_all_children;
Next, I merged the four alignments for each individual (this made the sorted_merge*bam files).
sbatch SubMerge.sh
Which runs MergeFork.pl, shown here:
#!/usr/bin/perl
#
# merge bam files for each individual
#
use Parallel::ForkManager;
my $max = 12;
my $pm = Parallel::ForkManager->new($max);
open(IN, "ids.txt") or die "failed to read\n";
while(<IN>){
$pm->start and next; ## fork
chomp;
$id = $_;
$out = "merge_$id".".bam";
system "samtools merge -nr $out *$id*bam\n";
system "samtools sort -O BAM -o sorted_$out $out\n";
system "samtools index sorted_$out\n";
$pm->finish;
}
$pm->wait_all_children;
Finally, I removed PCR duplicates for the merged bam files (this made the uniqe_merge*bam files, which are the ones actually used for variant calling).
sbatch SubRmdup.sh
Which runs RmdupFork.pl, shown here:
#!/usr/bin/perl
#
# removes pcr duplicates
#
use Parallel::ForkManager;
my $max = 24;
my $pm = Parallel::ForkManager->new($max);
FILES:
foreach $bam (@ARGV){
$pm->start and next FILES; ## fork
$ubam = $bam;
$ubam =~ s/sorted/uniqe/ or die "failed sub $ubam\n";
system "samtools rmdup -S $bam $ubam\n";
system "samtools index -b $ubam\n";
$pm->finish;
}
$pm->wait_all_children;
4. Step 1 of variant calling with GATK HaplotypeCaller (v3.5-0-g36282e4) to generate g.vcf format files.
The job was submitted with,
sbatch RunGatk.sh from /uufs/chpc.utah.edu/common/home/u6000989/data/timema/combind_wgs_dovetailV3/alignments_tcr
It runs GatkFork.pl on un*bam, this looks like,
#!/usr/bin/perl
#
# make g.vcf files
#
use Parallel::ForkManager;
my $max = 48;
my $pm = Parallel::ForkManager->new($max);
my $idir = '/uufs/chpc.utah.edu/common/home/u6000989/data/timema/combind_wgs_dovetailV3/alignments_tcr';
my $odir = '/uufs/chpc.utah.edu/common/home/u6000989/data/timema/combind_wgs_dovetailV3/variants_tcr';
my $tdir = '/scratch/general/lustre/tcrGatk';
my $genome = "/uufs/chpc.utah.edu/common/home/u6000989/data/timema/tcrDovetail/version3/mod_map_timema_06Jun2016_RvNkF702.fasta";
FILES:
foreach $bam (@ARGV){
$pm->start and next FILES; ## fork
$out = $bam;
$out =~ s/bam/g.vcf/ or die "failed here: $out\n";
$out = "re_$out";
system "java -Xmx48g -jar ~/bin/GenomeAnalysisTK.jar -T HaplotypeCaller -R $genome -I $idir/$bam -o $tdir/$out -XL na.intervals -gt_mode DISCOVERY -hets 0.001 -mbq 20 -out_mode EMIT_VARIANTS_ONLY -ploidy 2 -stand_call_conf 50 -pcrModel AGGRESSIVE --emitRefConfidence GVCF --variant_index_type LINEAR --variant_index_parameter 128000\n";
system "cp $tdir/$out $odir/$out\n";
$pm->finish;
}
$pm->wait_all_children;
The results are the re*g.vcf files in variants_tcr.
5. Combine the g.vcf files for joint variant calling. Without doing this, joint variant calling can lag or fail. I ran,
sbatch runComb.sh
Which submits four jobs via a fork program. The four jobs are in jobs.txt and are as follows. This generates the four combinded*g.vcf files.
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 re_uniqe_merge_timemaHVA_8021022.g.vcf --variant re_uniqe_merge_timemaHVA_8021024.g.vcf --variant re_uniqe_merge_timemaHVA_8021026.g.vcf --variant re_uniqe_merge_timemaHVA_8021027.g.vcf --variant re_uniqe_merge_timemaHVA_8021028.g.vcf --variant re_uniqe_merge_timemaHVA_8021029.g.vcf --variant re_uniqe_merge_timemaHVA_8021030.g.vcf --variant re_uniqe_merge_timemaHVA_8021034.g.vcf --variant re_uniqe_merge_timemaHVA_8021035.g.vcf --variant re_uniqe_merge_timemaHVA_8021037.g.vcf --variant re_uniqe_merge_timemaHVA_8021038.g.vcf --variant re_uniqe_merge_timemaHVA_8021039.g.vcf --variant re_uniqe_merge_timemaHVA_8021040.g.vcf --variant re_uniqe_merge_timemaHVA_8091001.g.vcf --variant re_uniqe_merge_timemaHVA_8091002.g.vcf --variant re_uniqe_merge_timemaHVA_8091003.g.vcf --variant re_uniqe_merge_timemaHVA_8091004.g.vcf --variant re_uniqe_merge_timemaHVA_8091005.g.vcf --variant re_uniqe_merge_timemaHVA_8091009.g.vcf --variant re_uniqe_merge_timemaHVA_8091010.g.vcf --variant re_uniqe_merge_timemaHVC_8021001.g.vcf --variant re_uniqe_merge_timemaHVC_8021002.g.vcf --variant re_uniqe_merge_timemaHVC_8021003.g.vcf --variant re_uniqe_merge_timemaHVC_8021004.g.vcf --variant re_uniqe_merge_timemaHVC_8021005.g.vcf --variant re_uniqe_merge_timemaHVC_8021009.g.vcf --variant re_uniqe_merge_timemaHVC_8021010.g.vcf --variant re_uniqe_merge_timemaHVC_8021011.g.vcf --variant re_uniqe_merge_timemaHVC_8021012.g.vcf --variant re_uniqe_merge_timemaHVC_8021013.g.vcf --variant re_uniqe_merge_timemaHVC_8021014.g.vcf --variant re_uniqe_merge_timemaHVC_8021015.g.vcf --variant re_uniqe_merge_timemaHVC_8021016.g.vcf --variant re_uniqe_merge_timemaHVC_8021019.g.vcf --variant re_uniqe_merge_timemaHVC_8021020.g.vcf --variant re_uniqe_merge_timemaHVC_8101011.g.vcf --variant re_uniqe_merge_timemaHVC_8101012.g.vcf --variant re_uniqe_merge_timemaHVC_8101013.g.vcf --variant re_uniqe_merge_timemaLA_10160904.g.vcf --variant re_uniqe_merge_timemaLA_10160906.g.vcf --variant re_uniqe_merge_timemaLA_12070901.g.vcf --variant re_uniqe_merge_timemaLA_12070902.g.vcf --variant re_uniqe_merge_timemaLA_12070903.g.vcf --variant re_uniqe_merge_timemaLA_12070904.g.vcf --variant re_uniqe_merge_timemaLA_12070905.g.vcf --variant re_uniqe_merge_timemaLA_12070906.g.vcf --variant re_uniqe_merge_timemaLA_12070907.g.vcf --variant re_uniqe_merge_timemaLA_12070908.g.vcf -o 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 re_uniqe_merge_timemaLA_12070910.g.vcf --variant re_uniqe_merge_timemaLA_12070911.g.vcf --variant re_uniqe_merge_timemaLA_4091005.g.vcf --variant re_uniqe_merge_timemaLA_4091006.g.vcf --variant re_uniqe_merge_timemaLA_4091007.g.vcf --variant re_uniqe_merge_timemaLA_4091008.g.vcf --variant re_uniqe_merge_timemaLA_4091009.g.vcf --variant re_uniqe_merge_timemaLA_4091010.g.vcf --variant re_uniqe_merge_timemaLA_4091012.g.vcf --variant re_uniqe_merge_timemaMR1A_8051001.g.vcf --variant re_uniqe_merge_timemaMR1A_8051002.g.vcf --variant re_uniqe_merge_timemaMR1A_8051003.g.vcf --variant re_uniqe_merge_timemaMR1A_8051004.g.vcf --variant re_uniqe_merge_timemaMR1A_8051005.g.vcf --variant re_uniqe_merge_timemaMR1A_8051006.g.vcf --variant re_uniqe_merge_timemaMR1A_8051007.g.vcf --variant re_uniqe_merge_timemaMR1A_8051008.g.vcf --variant re_uniqe_merge_timemaMR1A_8051009.g.vcf --variant re_uniqe_merge_timemaMR1A_8051010.g.vcf --variant re_uniqe_merge_timemaMR1A_8051011.g.vcf --variant re_uniqe_merge_timemaMR1A_8051012.g.vcf --variant re_uniqe_merge_timemaMR1A_8051014.g.vcf --variant re_uniqe_merge_timemaMR1A_8051017.g.vcf --variant re_uniqe_merge_timemaMR1A_8051018.g.vcf --variant re_uniqe_merge_timemaMR1A_8051019.g.vcf --variant re_uniqe_merge_timemaMR1A_8051020.g.vcf --variant re_uniqe_merge_timemaMR1A_8051021.g.vcf --variant re_uniqe_merge_timemaMR1A_8051022.g.vcf --variant re_uniqe_merge_timemaMR1A_8051023.g.vcf --variant re_uniqe_merge_timemaMR1C_8111003.g.vcf --variant re_uniqe_merge_timemaMR1C_8111004.g.vcf --variant re_uniqe_merge_timemaMR1C_8111006.g.vcf --variant re_uniqe_merge_timemaMR1C_8111007.g.vcf --variant re_uniqe_merge_timemaMR1C_8111008.g.vcf --variant re_uniqe_merge_timemaMR1C_8111009.g.vcf --variant re_uniqe_merge_timemaMR1C_8111010.g.vcf --variant re_uniqe_merge_timemaMR1C_8111011.g.vcf --variant re_uniqe_merge_timemaMR1C_8111012.g.vcf --variant re_uniqe_merge_timemaMR1C_8111013.g.vcf --variant re_uniqe_merge_timemaMR1C_8111014.g.vcf --variant re_uniqe_merge_timemaMR1C_8111015.g.vcf --variant re_uniqe_merge_timemaMR1C_8111017.g.vcf --variant re_uniqe_merge_timemaMR1C_8111018.g.vcf --variant re_uniqe_merge_timemaMR1C_8111019.g.vcf --variant re_uniqe_merge_timemaMR1C_8111020.g.vcf --variant re_uniqe_merge_timemaMR1C_8111021.g.vcf --variant re_uniqe_merge_timemaMR1C_8111022.g.vcf --variant re_uniqe_merge_timemaMR1C_8111023.g.vcf -o 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 re_uniqe_merge_timemaMR1C_8111024.g.vcf --variant re_uniqe_merge_timemaPRC_12010901.g.vcf --variant re_uniqe_merge_timemaPRC_12010902.g.vcf --variant re_uniqe_merge_timemaPRC_12010903.g.vcf --variant re_uniqe_merge_timemaPRC_12010904.g.vcf --variant re_uniqe_merge_timemaPRC_12010905.g.vcf --variant re_uniqe_merge_timemaPRC_12010906.g.vcf --variant re_uniqe_merge_timemaPRC_12010907.g.vcf --variant re_uniqe_merge_timemaPRC_12010908.g.vcf --variant re_uniqe_merge_timemaPRC_12010909.g.vcf --variant re_uniqe_merge_timemaPRC_12010910.g.vcf --variant re_uniqe_merge_timemaPRC_12010911.g.vcf --variant re_uniqe_merge_timemaPRC_12010912.g.vcf --variant re_uniqe_merge_timemaPRC_12070912.g.vcf --variant re_uniqe_merge_timemaPRC_12070913.g.vcf --variant re_uniqe_merge_timemaPRC_12070914.g.vcf --variant re_uniqe_merge_timemaPRC_12070915.g.vcf --variant re_uniqe_merge_timemaPRC_12070916.g.vcf --variant re_uniqe_merge_timemaPRC_4091002.g.vcf --variant re_uniqe_merge_timemaPRC_4091003.g.vcf --variant re_uniqe_merge_timemaR12A_10201001.g.vcf --variant re_uniqe_merge_timemaR12A_10201002.g.vcf --variant re_uniqe_merge_timemaR12A_10201003.g.vcf --variant re_uniqe_merge_timemaR12A_10201004.g.vcf --variant re_uniqe_merge_timemaR12A_10201005.g.vcf --variant re_uniqe_merge_timemaR12A_10201006.g.vcf --variant re_uniqe_merge_timemaR12A_10201007.g.vcf --variant re_uniqe_merge_timemaR12A_10201008.g.vcf --variant re_uniqe_merge_timemaR12A_10201009.g.vcf --variant re_uniqe_merge_timemaR12A_10201010.g.vcf --variant re_uniqe_merge_timemaR12A_10201011.g.vcf --variant re_uniqe_merge_timemaR12A_10201012.g.vcf --variant re_uniqe_merge_timemaR12A_10201013.g.vcf --variant re_uniqe_merge_timemaR12A_10201014.g.vcf --variant re_uniqe_merge_timemaR12A_10201015.g.vcf --variant re_uniqe_merge_timemaR12A_10201016.g.vcf --variant re_uniqe_merge_timemaR12A_10201017.g.vcf --variant re_uniqe_merge_timemaR12A_10201018.g.vcf --variant re_uniqe_merge_timemaR12A_10201019.g.vcf --variant re_uniqe_merge_timemaR12A_10201020.g.vcf --variant re_uniqe_merge_timemaR12A_10201021.g.vcf --variant re_uniqe_merge_timemaR12C_10221001.g.vcf --variant re_uniqe_merge_timemaR12C_10221002.g.vcf --variant re_uniqe_merge_timemaR12C_10221003.g.vcf --variant re_uniqe_merge_timemaR12C_10221004.g.vcf --variant re_uniqe_merge_timemaR12C_10221005.g.vcf --variant re_uniqe_merge_timemaR12C_10221006.g.vcf --variant re_uniqe_merge_timemaR12C_10221007.g.vcf -o combinded_3.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 re_uniqe_merge_timemaR12C_10221008.g.vcf --variant re_uniqe_merge_timemaR12C_10221009.g.vcf --variant re_uniqe_merge_timemaR12C_10221010.g.vcf --variant re_uniqe_merge_timemaR12C_10221011.g.vcf --variant re_uniqe_merge_timemaR12C_10221012.g.vcf --variant re_uniqe_merge_timemaR12C_10221013.g.vcf --variant re_uniqe_merge_timemaR12C_10221014.g.vcf --variant re_uniqe_merge_timemaR12C_10221015.g.vcf --variant re_uniqe_merge_timemaR12C_10221016.g.vcf --variant re_uniqe_merge_timemaR12C_10221017.g.vcf --variant re_uniqe_merge_timemaR12C_10221018.g.vcf --variant re_uniqe_merge_timemaR12C_10221019.g.vcf --variant re_uniqe_merge_timemaR12C_10221020.g.vcf --variant re_uniqe_merge_timemaR12C_10221021.g.vcf -o combinded_4.g.vcf
6. Joint variant calling with GATK. From the variants_tcr sub-directory I ran
sbatch jointVarCall.sh
which runs
#!/bin/sh
#SBATCH --time=140: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/u6000989/data/timema/combind_wgs_dovetailV3/variants_tcr
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 combinded_1.g.vcf --variant combinded_2.g.vcf --variant combinded_3.g.vcf --variant combinded_4.g.vcf -o tcr_wgs_variants.vcf
This generates the key variant file, tcr_wgs_variants.vcf.