genotyping with mpileup

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