Now that we have our reads aligned to a reference we can call genotypes for our samples
mpileup is an easy to use method for genotyping included in the samtools package
Run mpileup on our read alignments (sorted bam files):
run mpileup
#start and idev session if you haven't idev #copy over the exercise directory cds cd rad_intro/ cp -r /work/02260/grovesd/lonestar/intro_to_rad_2017/genotyping/ddRAD_mpileup . cd ddRAD_mpileup/ #index the reference for samtools module load samtools samtools faidx stickleback_chrom3.fasta #make a list of the bam files ls *.bam > my_bamfiles.txt #run mpileup samtools mpileup -f stickleback_chrom3.fasta -t DP,AD,ADF,ADR,SP -u -b my_bamfiles.txt > mpileup_results.bcf #now call genotypes from the mpileup results bcftools call -vmO v -o raw_calls.vcf mpileup_results.bcf
Quality filter the variant calls:
quality filter variants
#use vcftools to get information about our variant set vcftools --vcf raw_calls.vcf #returns this: VCFtools - 0.1.15 (C) Adam Auton and Anthony Marcketta 2009 Parameters as interpreted: --vcf raw_calls.vcf After filtering, kept 3 out of 3 Individuals After filtering, kept 14117 out of a possible 14117 Sites #now quality filter the raw calls bcftools filter --exclude 'QUAL < 30' raw_calls.vcf | bcftools view > qced_calls.vcf #check the new quality checked vcf vcftools --vcf qced_calls.vcf
