Post date: Feb 08, 2015 4:1:7 AM
1. Alignment with bwa to Timema draft genome v 0.3
cd /labs/evolution/data/timema/timema_experiments/within_generation/ecology_letters_gbs/bam
bwa aln -n 4 -l 20 -k 2 -t 8 -q 10 -Y -f alnCA25093-5C.sai mod_lg_timemaGenome.fasta /labs/evolution/data/timema/timema_experiments/within_generation/ecology_letters_gbs/fq/CA25093-5C.fastq
bwa samse -n 1 -r '@RG ID:tcrist-CA25093-5C' -f alnCA25093-5C.sam mod_lg_timemaGenome.fasta alnCA25093-5C.sai /labs/evolution/data/timema/timema_experiments/within_generation/ecology_letters_gbs/fq/CA25093-5C.fastq
2. convert sam to bam
3. variant calling with samtools and bcftools
cd /labs/evolution/data/timema/timema_experiments/within_generation/ecology_letters_gbs/bam
samtools mpileup -C 50 -d 50 -f mod_lg_timemaGenome.fasta -q 10 -Q 15 -D -g -I *bam > tcristGbsVars.bcf
bcftools view -c -d 0.8 -g -p 0.01 -P full -t 0.001 -v tcristGbsVars.bcf > tcristGbsVars.vcf
4. filter variants (now in /labs/evolution/projects/timema_wgexperiment/gbsgemma/variants)
perl vcfFilter tcristGbsVars.vcf
Filters: 1500 reads total, 20 reads of the alt. allele, no more then one reverse read, bi-allelic
355,265 SNPs
5. Drop low maf alleles and convert to gl format
perl subvcf2gl.pl filtered_tcristGbsVars.vcf
Filters: minimum maf 0.01
294,438 SNPs
6. Generate point estimates using the allele frequency prior (p^2, 2pq, q^2), mean or mode:
perl gl2genest.pl tcristgbs.gl af_tcristgbs.txt