Post date: Nov 10, 2016 11:15:55 PM
BWA Alignment
Alignment job details:
$aln ="bwa aln -n 4 -l 20 -k 2 -t 8 -q 10 -f aln"."$ind".".sai $genome $file\n";
Call Var
samtools mpileup -b bams -C 50 -d 250 -f /uufs/chpc.utah.edu/common/home/u6000989/data/callosobruchus/gbs_L14/alignments/cmacgenome.fasta -q 20 -Q 30 -I -u -g -t DP,DPR -o calloVars.bcf
~/bin/bcftools call -v -m -p 0.01 -P 0.001 -O v -o calloVars.vcf calloVars.bcf
Filter Variants based on depth, mapping quality, etc; using version with data for >= 70% of individuals, keeping 28,249 SNPs. Creates filtered2x_calloVars.vcf
perl ../scripts/vcfFilter.pl calloVars.vcf
#stringency variables
my $minCoverage = 1344; # minimum number of sequences; DP, 1344 = 2X
my $minAltRds = 10; # minimum number of sequences with the alternative allele; AC
my $notFixed = 1.0; # removes loci fixed for alt; AF
my $mq = 30; # minimum mapping quality; MQ
my $miss = 201; # maximum number of individuals with no data
my $minMaf = 6; # minimum MAF, as allele count ~ 0.005
my $ratDifDir = 0.01; ## maximum proportion of reads in reverse (relative to avg.) orientation
Filter some more based on high SNP density. SNPs must be at least 4bp apart, plus gets rid of high coverage snps. Creates morefilter_filtered2x_calloVars.vcf
#stringency variables
my $maxCoverage = 21556; # maximum depth to avoid repeats, mean + 2sd
my $mind = 4; # don't call if SNPs are closer together than this
Convert to genotype likelihood format. Creates cmacL14.gl
perl ../scripts/vcf2gl.pl 0 morefilter_filtered2x_calloVars.vcf
Split genotype likelihoods by population
perl projects/cmac_AR/scripts/splitsPops.pl projects/cmac_AR/likelihood_data/cmacL14.gl
Create allele freq by population:
sbatch ../scripts/wrap_qsub_slurm_popmod.pl ../likelihood_data/cmacseries*.gl
paremeters:
$job ="~/bin/popmod -i $file -o outp_bayes_$pop".".hdf5 -b 1000 -n 10000 -t 5 -a 0.5\n";
Create txt files of allele frequency (p) and genotypes (g)
module load hdf5 gsl gcc
estpost_popmod -s 0 -w 0 -p p ...
~/bin/estpost_popmod -p p -s 0 -w 0 -o p_L1-F100.txt outp_bayes_L1-F100.hdf5
~/bin/estpost_popmod -p p -s 0 -w 0 -o p_L2-F87.txt outp_bayes_L2-F87.hdf5
~/bin/estpost_popmod -p p -s 0 -w 0 -o p_L3-F85.txt outp_bayes_L3-F85.hdf5
~/bin/estpost_popmod -p p -s 0 -w 0 -o p_M-14.txt outp_bayes_M-14.hdf5