...
Following steps from https://github.com/zgompert/DimensionsExperiment.
Built bcftools version 1.16 and installed in
/project/evolgen/bin/.bcftoolsneeded reference genome in bzip2 format, not gzip. So I now simply have an unzipped reference genome, which I have reindexed.Completed this step with something like:
sbatch --account=evolgen --time=1-00:00 --nodes=1 --mem=8G --mail-type=ENDÂ 0_call_variants.sh(this took 12 hours and 40 minutes and 552 MB of RAM; I asked for 120GB, which likely gave me the whole node and made it a bit faster)Filtered vcf with
1_filter_variants.sh, which contains notes on the criteria that I used (could be altered to suit). This set is based on a fairly tight set of criteria (it matches what we used for a recent paper), which could be modified as needed. Note that currently there is no explicit minor allele frequency filtering. There are 5016 sites inhmax_variants_filtered.vcfCode Block ## minimum mapping quality of 30 was already enforced in bcftools mpileup ## These are written as exclusion filters # ------ INFO/DP < 952 # 2x depth overall: obtained with INFO/DP > 951 (476 x 2 = 952) ##INFO=<ID=DP,Number=1,Type=Integer,Description="Raw read depth"> # ------ INFO/AC1 < 10 ## a minimum of 10 alt reads to support a polymorphism # ------ INFO/BQBZ > 1e-5 ##INFO=<ID=BQBZ,Number=1,Type=Float,Description="Mann-Whitney U-z test of Base Quality Bias (closer to 0 is better)"> # ------ INFO/RPBZ > 1e-5 ##INFO=<ID=RPBZ,Number=1,Type=Float,Description="Mann-Whitney U-z test of Read Position Bias (closer to 0 is better)"> ## biallelic snps obtained in bcftools view ## 476 individuals (80% with data would mean 380 individuals; 380/476=0.798) ## set bcftools view to include only sites with the fraction of missingness less than 0.2
...