GOALS
- To learn how to use VCFtools to filter a VCF file for missing data, genotype depth, locus quality score, minor allele frequency, and genotype call depth
- To learn how to use vcflib to filter FreeBayes VCF files generated with RAD data
- To filter a VCF file for HWE within populations
- How to decompose a VCF into SNPs and INDELs and
- How to use a haplotyping script to further filter SNPs for paralogs and genotyping errors.
Tutorial
Welcome to the SNP filtering exercise. For the first part of the exercise, the filtering steps should work on almost any VCF file.
For the second part of the exercise, we are going to assume you are working with a VCF file that was generated by
FreeBayes. Note that other SNP callers can be configured to include the same annotations.
Let’s find our way back to your original working directory and make a new filtering directory
mkdir filtering
cd filtering
Now, let’s download some data to look at.
curl -L -o data.zip https://www.dropbox.com/sh/bf9jxviaoq57s5v/AAD2Kv5SPpHlZ7LC7sBz4va8a?dl=1
unzip data.zip
ll
total 165620
-rwxr--r--. 1 jpuritz users 109633 Mar 6 14:56 BR_004-RG.bam
-rwxr--r--. 1 jpuritz users 247496 Mar 6 14:57 BR_004-RG.bam.bai
-rwxr--r--. 1 jpuritz users 120045 Mar 6 15:14 BR_006-RG.bam
-rwxr--r--. 1 jpuritz users 247712 Mar 6 15:14 BR_006-RG.bam.bai
-rw-r--r--. 1 jpuritz users 78979977 Mar 6 16:15 data.zip
drwxr-xr-x. 2 jpuritz users 21 Mar 6 16:16 __MACOSX
-rwxr--r--. 1 jpuritz users 399 Mar 6 15:08 popmap
-rwxr--r--. 1 jpuritz users 68264393 Mar 6 13:40 raw.vcf.gz
-rwxr--r--. 1 jpuritz users 6804314 Mar 6 14:49 reference.fasta
-rwxr--r--. 1 jpuritz users 1085387 Mar 6 14:49 reference.fasta.amb
-rwxr--r--. 1 jpuritz users 1068884 Mar 6 14:49 reference.fasta.ann
-rwxr--r--. 1 jpuritz users 6379720 Mar 6 14:49 reference.fasta.bwt
-rwxr--r--. 1 jpuritz users 353544 Mar 6 14:49 reference.fasta.clstr
-rwxr--r--. 1 jpuritz users 976388 Mar 6 14:49 reference.fasta.fai
-rwxr--r--. 1 jpuritz users 1594913 Mar 6 14:49 reference.fasta.pac
-rwxr--r--. 1 jpuritz users 3189872 Mar 6 14:49 reference.fasta.sa
-rwxr--r--. 1 jpuritz users 137209 Mar 6 14:30 stats.out
To start, we are going to use the program VCFtools (http://vcftools.sourceforge.net) to filter our vcf file. This program has a binary executable
and has several perl scripts as well that are useful for filtering.
I find it much more useful to use version 0.1.11, since it has more useful filtering commands (I think). Let’s load that version
This raw.vcf file is going to have a lot of erroneous variant calls and a lot of variants that are only present in one individual.
To make this file more manageable, let’s start by applying three step filter. We are going to only keep variants that have been successfully genotyped in
50% of individuals, a minimum quality score of 30, and a minor allele count of 3.
vcftools --gzvcf raw.vcf.gz --max-missing 0.5 --mac 3 --minQ 30 --recode --recode-INFO-all --out raw.g5mac3
In this code, we call vcftools, feed it a vcf file after the --vcf flag, --max-missing 0.5 tells it to filter genotypes called below 50% (across all individuals)
the --mac 3 flag tells it to filter SNPs that have a minor allele count less than 3.
This is relative to genotypes, so it has to be called in at least 1 homozygote and 1 heterozygote or 3 heterozygotes.
The --recode flag tells the program to write a new vcf file with the filters, --recode-INFO-all keeps all the INFO flags from the old vcf file in the new one.
Lastly, --out designates the name of the output
The output will scroll through a lot of lines, but should end like:
After filtering, kept 40 out of 40 Individuals
After filtering, kept 78434 out of a possible 147540 Sites
Outputting VCF file... Done
Run Time = 40.00 seconds
Those two simple filters got rid of 50% of the data and will make the next filtering steps run much faster.
We now have a filtered VCF called raw.g5mac3.recode.vcf. There is also a logfile generated called raw.g5mac3.log The next filter we will apply is a minimum depth for a genotype call and a minimum mean depth
vcftools --vcf raw.g5mac3.recode.vcf --minDP 3 --recode --recode-INFO-all --out raw.g5mac3dp3
This command will recode genotypes that have less than 3 reads. I’ll give you a second to take a deep breath. Yes, we are keeping genotypes with as few as 3 reads. We talked about this in the lecture portion of this course, but the short answer is that sophisticated multisample variant callers like FreeBayes and GATK can confidently call genotypes with few reads because variants are assessed across all samples simultaneously. So, the genotype is based on three reads AND prior information from all reads from all individuals. Relax. We will do plenty of other filtering steps!
Don’t believe me do you? I’ve made a script to help evaluate the potential errors.
curl -L -O https://github.com/jpuritz/dDocent/raw/master/scripts/ErrorCount.sh
chmod +x ErrorCount.sh
./ErrorCount.sh raw.g5mac3dp3.recode.vcf
This script counts the number of potential genotyping errors due to low read depth
It report a low range, based on a 50% binomial probability of observing the second allele in a heterozygote and a high range based on a 25% probability.
Potential genotyping errors from genotypes from only 1 read range from 0 to 0.0
Potential genotyping errors from genotypes from only 2 reads range from 0 to 0.0
Potential genotyping errors from genotypes from only 3 reads range from 15986 to 53714.22
Potential genotyping errors from genotypes from only 4 reads range from 6230 to 31502.04
Potential genotyping errors from genotypes from only 5 reads range from 2493 to 18914
40 number of individuals and 78434 equals 3137360 total genotypes
Total genotypes not counting missing data 2380094
Total potential error rate is between 0.0103815227466 and 0.0437504821238
SCORCHED EARTH SCENARIO
WHAT IF ALL LOW DEPTH HOMOZYGOTE GENOTYPES ARE ERRORS?????
The total SCORCHED EARTH error rate is 0.129149100834.
Right now, the maximum error rate for our VCF file because of genotypes less than 5 reads is less than 5%. See, nothing to worry about.
The next step is to get rid of individuals that did not sequence well. We can do this by assessing individual levels of missing data.
vcftools --vcf raw.g5mac3dp3.recode.vcf --missing-indv
This will create an output called out.imiss. Let’s examine it.
cat out.imiss
INDV N_DATA N_GENOTYPES_FILTERED N_MISS F_MISS
BR_002 78434 0 13063 0.166548
BR_004 78434 0 16084 0.205064
BR_006 78434 0 25029 0.319109
BR_009 78434 0 30481 0.38862
BR_013 78434 0 69317 0.883762
BR_015 78434 0 8861 0.112974
BR_016 78434 0 29789 0.379797
BR_021 78434 0 17422 0.222123
BR_023 78434 0 43913 0.559872
BR_024 78434 0 24220 0.308795
BR_025 78434 0 21998 0.280465
BR_028 78434 0 26786 0.34151
BR_030 78434 0 74724 0.952699
BR_031 78434 0 26488 0.337711
BR_040 78434 0 19492 0.248515
BR_041 78434 0 17107 0.218107
BR_043 78434 0 16384 0.208889
BR_046 78434 0 28770 0.366805
BR_047 78434 0 13258 0.169034
BR_048 78434 0 24505 0.312428
WL_031 78434 0 22566 0.287707
WL_032 78434 0 22604 0.288191
WL_054 78434 0 32902 0.419486
WL_056 78434 0 34106 0.434837
WL_057 78434 0 37556 0.478823
WL_058 78434 0 31448 0.400949
WL_061 78434 0 35671 0.45479
WL_064 78434 0 47816 0.609634
WL_066 78434 0 10062 0.128286
WL_067 78434 0 47940 0.611215
WL_069 78434 0 38260 0.487799
WL_070 78434 0 21188 0.270138
WL_071 78434 0 16692 0.212816
WL_072 78434 0 46347 0.590904
WL_076 78434 0 78178 0.996736
WL_077 78434 0 55193 0.703687
WL_078 78434 0 54400 0.693577
WL_079 78434 0 19457 0.248068
WL_080 78434 0 30076 0.383456
WL_081 78434 0 30334 0.386746
You can see that some individuals have as high as 99.6% missing data. We definitely want to filter those out. Let’s take a look at a histogram
mawk '!/IN/' out.imiss | cut -f5 > totalmissing
gnuplot << \EOF
set terminal dumb size 120, 30
set autoscale
unset label
set title "Histogram of % missing data per individual"
set ylabel "Number of Occurrences"
set xlabel "% of missing data"
#set yr [0:100000]
binwidth=0.01
bin(x,width)=width*floor(x/width) + binwidth/2.0
plot 'totalmissing' using (bin($1,binwidth)):(1.0) smooth freq with boxes
pause -1
EOF
Histogram of % missing data per individual
Number of Occurrences
3 ++----------+---------***---------***-----------+------------+-----------+-----------+-----------+----------++
+ + * * * * + 'totalmissing' using (bin($1,binwidth)):(1.0) ****** +
| * * * * |
| * * * * |
| * * * * |
| * * * * |
2.5 ++ * * * * ++
| * * * * |
| * * * * |
| * * * * |
| * * * * |
| * * * * |
2 ++ *************** * **** * * ++
| * * * ** * * ** * * * |
| * * * ** * * ** * * * |
| * * * ** * * ** * * * |
| * * * ** * * ** * * * |
| * * * ** * * ** * * * |
1.5 ++ * * * ** * * ** * * * ++
| * * * ** * * ** * * * |
| * * * ** * * ** * * * |
| * * * ** * * ** * * * |
| * * * ** * * ** * * * |
+ * * +* ** * * ** * * * + + + + + +
1 +*************************************************************************************************************
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
% of missing data
Most of the individuals have less than 0.5 missing data. That is probably a good cutoff to use for the moment.
Now we need to create a list of individuals with more than 50% missing data. Anyone have any ideas?
We can use mawk to do it.
mawk '$5 > 0.5' out.imiss | cut -f1 > lowDP.indv
Take a moment to think about what this code is doing.
Now that we have a list of individuals to remove, we can feed that directly into VCFtools for filtering.
vcftools --vcf raw.g5mac3dp3.recode.vcf --remove lowDP.indv --recode --recode-INFO-all --out raw.g5mac3dplm
As you can see from the output, this removed 9 individuals.
I’ve included a script called filter_missing_ind.sh that will automate this process for you in the future. Try it out.
curl -L -O https://github.com/jpuritz/dDocent/raw/master/scripts/filter_missing_ind.sh
chmod +x filter_missing_ind.sh
./filter_missing_ind.sh raw.g5mac3dp3.recode.vcf DP3g95maf05
The command always follows the structure of filter_missing_ind.sh vcf_to_filter name_prefix_for_new_vcf The script prints out a histogram like the one above and also calculates the 85% for missing data. Enter “no” Now that we have removed poor coverage individuals, we can restrict the data to variants called in a high percentage of individuals and filter by mean depth of genotypes
vcftools --vcf raw.g5mac3dplm.recode.vcf --max-missing 0.95 --maf 0.05 --recode --recode-INFO-all --out DP3g95maf05 --min-meanDP 20
This leaves us with about 12,754 loci in our filtered vcf file.
This applied a genotype call rate (95%) across all individuals. With two localities, this is sufficient, but when you have multiple localities being sampled You are also going to want to filter by a population specific call rate. VCFtools won’t calculate this directly, but it is an easy workaround. First we need a file to define localities (populations). Most programs want the file to have two tab separated columns. First with the sample name, second with population assignment. I’ve already made one for this exercise.
cat popmap
BR_002 BR
BR_004 BR
BR_006 BR
BR_009 BR
BR_013 BR
BR_015 BR
BR_016 BR
BR_021 BR
BR_023 BR
BR_024 BR
BR_025 BR
BR_028 BR
BR_030 BR
BR_031 BR
BR_040 BR
BR_041 BR
BR_043 BR
BR_046 BR
BR_047 BR
BR_048 BR
WL_031 WL
WL_032 WL
WL_054 WL
WL_056 WL
WL_057 WL
WL_058 WL
WL_061 WL
WL_064 WL
WL_066 WL
WL_067 WL
WL_069 WL
WL_070 WL
WL_071 WL
WL_072 WL
WL_076 WL
WL_077 WL
WL_078 WL
WL_079 WL
WL_080 WL
WL_081 WL
Now we need to create two lists that have just the individual names for each population
mawk '$2 == "BR"' popmap > 1.keep && mawk '$2 == "WL"' popmap > 2.keep
The above line demonstrates the use of && to simultaneous execute two tasks.
Next, we use VCFtools to estimate missing data for loci in each population
vcftools --vcf DP3g95maf05.recode.vcf --keep 1.keep --missing-site --out 1
vcftools --vcf DP3g95maf05.recode.vcf --keep 2.keep --missing-site --out 2
This will generate files named 1.lmiss and 2.lmiss
They follow this format
head -3 1.lmiss
CHR POS N_DATA N_GENOTYPE_FILTERED N_MISS F_MISS
E1_L101 9 34 0 0 0
E1_L101 15 34 0 0 0
I added extra tabs to make this easier to read, but what we are interested in is that last column with is the percentage of missing data for that locus. We can combine the two files and make a list of loci about the threshold of 10% missing data to remove. Note this is double the overall rate of missing data.
cat 1.lmiss 2.lmiss | mawk '!/CHR/' | mawk '$6 > 0.1' | cut -f1,2 >> badloci
Who can walk us through that line of code?
We then feed this file back into VCFtools to remove any of the loci
vcftools --vcf DP3g95maf05.recode.vcf --exclude-positions badloci --recode --recode-INFO-all --out DP3g95p5maf05
Again, we only had two populations so our overall filter caught all of these. However, this will not be the case in multi-locality studies I also have made a script to automate this process as well. It’s called pop_missing_filter.sh Executing it with no parameters will give you the usage.
curl -L -O https://github.com/jpuritz/dDocent/raw/master/scripts/pop_missing_filter.sh
chmod +x pop_missing_filter.sh
./pop_missing_filter.sh
Usage is pop_missing_filter vcffile popmap percent_missing_per_pop number_of_pops_for_cutoff name_for_output
From this point forward, the filtering steps assume that the vcf file was generated by FreeBayes
Note that other SNP callers can be configured to include the similar annotations.
FreeBayes outputs a lot of information about a locus in the VCF file, using this information and the properties of RADseq, we add some sophisticated filters to the data. Let’s take a look at the header of our VCF file and take a quick look at all the information.
mawk '/#/' DP3g95maf05.recode.vcf
This will output several lines of INFO tags, I have highlighted a few below:
INFO=<ID=NS,Number=1,Type=Integer,Description="Number of samples with data">
INFO=<ID=DP,Number=1,Type=Integer,Description="Total read depth at the locus">
INFO=<ID=QR,Number=1,Type=Integer,Description="Reference allele quality sum in phred">
INFO=<ID=QA,Number=A,Type=Integer,Description="Alternate allele quality sum in phred">
INFO=<ID=SRF,Number=1,Type=Integer,Description="Number of reference observations on the forward strand">
INFO=<ID=SRR,Number=1,Type=Integer,Description="Number of reference observations on the reverse strand">
INFO=<ID=SAF,Number=A,Type=Integer,Description="Number of alternate observations on the forward strand">
INFO=<ID=SAR,Number=A,Type=Integer,Description="Number of alternate observations on the reverse strand">
INFO=<ID=AB,Number=A,Type=Float,Description="Allele balance at heterozygous sites: a number between 0 and 1 representing the ratio of reads showing the reference allele to all reads, considering only reads from individuals called as heterozygous">
INFO=<ID=TYPE,Number=A,Type=String,Description="The type of allele, either snp, mnp, ins, del, or complex.">
INFO=<ID=CIGAR,Number=A,Type=String,Description="The extended CIGAR representation of each alternate allele, with the exception that '=' is replaced by 'M' to ease VCF parsing. Note that INDEL alleles do not have the first matched base (which is provided by default, per the spec) referred to by the CIGAR.">
INFO=<ID=MQM,Number=A,Type=Float,Description="Mean mapping quality of observed alternate alleles">
INFO=<ID=MQMR,Number=1,Type=Float,Description="Mean mapping quality of observed reference alleles">
INFO=<ID=PAIRED,Number=A,Type=Float,Description="Proportion of observed alternate alleles which are supported by properly paired read fragments">
INFO=<ID=PAIREDR,Number=1,Type=Float,Description="Proportion of observed reference alleles which are supported by properly paired read fragments">
FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
FORMAT=<ID=GQ,Number=1,Type=Float,Description="Genotype Quality, the Phred-scaled marginal (or unconditional) probability of the called genotype">
FORMAT=<ID=GL,Number=G,Type=Float,Description="Genotype Likelihood, log10-scaled likelihoods of the data given the called genotype for each possible genotype generated from the reference and alternate alleles given the sample ploidy">
FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
FORMAT=<ID=RO,Number=1,Type=Integer,Description="Reference allele observation count">
FORMAT=<ID=QR,Number=1,Type=Integer,Description="Sum of quality of the reference observations">
FORMAT=<ID=AO,Number=A,Type=Integer,Description="Alternate allele observation count">
FORMAT=<ID=QA,Number=A,Type=Integer,Description="Sum of quality of the alternate observations">
The first filter we will apply will be on allele balance. Allele balance is: a number between 0 and 1 representing the ratio of reads showing the reference allele to all reads, considering only reads from individuals called as heterozygous Because RADseq targets specific locations of the genome, we expect that the allele balance in our data (for real loci) should be close to 0.5 We can use the vcffilter program from vcflib. (https://github.com/ekg/vcflib) Typing it with no parameters will give you the usage.
vcffilter
usage: vcffilter [options] <vcf file>
options:
-f, --info-filter specifies a filter to apply to the info fields of records,
removes alleles which do not pass the filter
-g, --genotype-filter specifies a filter to apply to the genotype fields of records
-k, --keep-info used in conjunction with '-g', keeps variant info, but removes genotype
-s, --filter-sites filter entire records, not just alleles
-t, --tag-pass tag vcf records as positively filtered with this tag, print all records
-F, --tag-fail tag vcf records as negatively filtered with this tag, print all records
-A, --append-filter append the existing filter tag, don't just replace it
-a, --allele-tag apply -t on a per-allele basis. adds or sets the corresponding INFO field tag
-v, --invert inverts the filter, e.g. grep -v
-o, --or use logical OR instead of AND to combine filters
-r, --region specify a region on which to target the filtering, requires a BGZF
compressed file which has been indexed with tabix. any number of
regions may be specified.
Let’s use our first filter
vcffilter -s -f "AB > 0.25 & AB < 0.75 | AB < 0.01" DP3g95p5maf05.recode.vcf > DP3g95p5maf05.fil1.vcf
vcffilter works with simple conditional statements, so this filters out loci with an allele balance below 0.25 and above 0.75. However, it does include those that are close to zero.
The last condition is to catch loci that are fixed variants (all individuals are homozygous for one of the two variants).
The -s tells the filter to apply to sites, not just alleles
To see how many loci are now in the VCF file, you could feed it into VCFtools or you can just use a simple mawk statement
mawk '!/#/' DP3g95p5maf05.recode.vcf | wc -l
12754
mawk '!/#/' DP3g95p5maf05.fil1.vcf | wc -l
9678
You’ll notice that we’ve filtered a lot of loci. In my experience though, I find that most of these tend to be errors of some kind. However, this will be data dependent. I encourage you to explore your own data sets.
The next filter we will apply filters out sites that have reads from both strands. For GWAS and even RNAseq, this can be a good thing.
Unless you are using super small genomic fragment or really long reads (MiSeq). A SNP should be covered only by forward or only reverse reads.
vcffilter -f "SAF / SAR > 100 & SRF / SRR > 100 | SAR / SAF > 100 & SRR / SRF > 100" -s DP3g95p5maf05.fil1.vcf > DP3g95p5maf05.fil2.vcf
The filter is based on proportions, so that a few extraneous reads won’t remove an entire locus. In plain english, it’s keeping loci that have over 100 times more forward alternate reads than reverse alternate reads and 100 times more forward reference reads than reverse reference reads along with the reciprocal.
mawk '!/#/' DP3g95p5maf05.fil2.vcf | wc -l
9491
That only removes a small proportion of loci, but these loci are likely to be paralogs, microbe contamination, or weird PCR chimeras.
The program SAMtools is a great way to visualize alignments right from the terminal.
samtools tview BR_006-RG.bam reference.fasta -p E28188_L151
11 21 31 41 51 61 71 81 91 101 111 121 131
AATTCTCAGAGCTAGAGTGGGGACGGCAGTTGGTAGAGGGTACAGCAGTTCTAAAAACATGTAGAAATTTTCTCTTCAACTCGCTCCTACGGCCACAGCGTTCACTCCACATACACAAATTGTACACCAAAACATAGGAAAAG
...........S...........Y.K......S.........G.......K.........S............................Y........Y....W.........................M...G.........
..........................................G.......G......................................T.....
..........................................G.......T............................................
...........G...........T.T......C.........G.................C...............................
..........................................G.......T....................................G.......
..........................................G.......T............................................
...........G...........T.T......C.........G.................C..................................
..........................................G.......T............................................
...........G.................C............G.......G.................******........A...............T..
,,,,,,g,,,,,,,,,,,,,,,,,******,,,,,,,,a,,,,,,,,,,,,,,,t,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,g,,,,,,,,,
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,t,,,,,,,,,,,,,,,,,,,,,,,,,,,,a,,,,,,,,,,,g,,,,,,,,,
t,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,t,,,,t,,,,,,,,,,,,,,,,,,,,,,,,,c,,,g,,,,,,,,,
,,,,,,,,,,c,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,a,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,g,,,,,,,,,
g,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,t,,,,,,,,,,,,,t,,,,,,,,,,,,,,,,,,,,,,,,,c,,,g,,,,,,,,,
t,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,g,,,,,,,,,,t,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,g,,,,,,,,,
g,,,,,,,,,a,,,,,,,,,,,,,,,,,,,,,,,,,,,,t,,,,,,,,,,,,,t,,,,,,,,,,,,,,,,,,,,,,,,,c,,,g,,,,,,,,,
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,t,,,,,,,,,,,,,,,,,,,,,,,,,,,,a,,,,,,,,,,,g,,,,,,,,,
,,,,,,,,,,,,,,,,,,,g,,,,,,,,,,,,,,,,,,,,,,t,,,,,,,,,,,,,,,,,,,,,,,,,,,,a,,,,,,,,,,,g,,,,,,,,,
t,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,t,,,,t,,,,,,,,,,,,,,,,,,,,,,,,,c,,,g,,,,,,,,,
g,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,t,,,,,,,,,,,,,t,,,,,,,,,a,,,,,,,,,,,,,,,c,,,g,,,,,,,,,
g,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,t,,,,,,,,,,,,,t,,,,,,,,,,,,,,,,,,,,,,,,,c,,,g,,,,,,,,,
,,,,,,,,,,c,g,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,t,,,,,,,,,,,,,,,,,,,,,,,g,,,,,,,,,
g,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,t,,,,,,,,,,,,,t,,,,,,,,,,,,,,,,,,,,,,,,,c,,,g,,,,,,,,,
t,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,t,,,,t,,,,,,,,,,,,,,,,,,,,,,,,,c,,,g,,,,,,,,,
g,,,,,,t,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,t,,,,,,,,,,,,,t,,,,,,,,,,,,,,,,,,,,,,,,,c,,,g,,,,,,,,,
g,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,t,,,,,,,,,,,,,t,,,,,,,,,,,,,,,,,,,,,,,,,c,,,g,,,,,,,,,
As you can see this is a mess. There appear to be over two haplotypes here.
For more info on how to use samtools tview press the question mark while you are in the window.
The next filter looks at the ratio of mapping qualities between reference and alternate alleles
vcffilter -f "MQM / MQMR > 0.9 & MQM / MQMR < 1.05" DP3g95p5maf05.fil2.vcf > DP3g95p5maf05.fil3.vcf
The rationale here is that, again, because RADseq loci and alleles all should start from the same genomic location there should not be large discrepancy between the mapping qualities of two alleles.
mawk '!/#/' DP3g95p5maf05.fil3.vcf | wc -l
9229
This filters away less than 3% of the variants, but likely need to be filtered. Let’s take a look at one
samtools tview BR_004-RG.bam reference.fasta -p E20_L173
1 11 21 31 41 51 61 71 81 91 101 111 121 131
NAATTCATCTGTTGCAGGCAGCTCACACTTGCAGCCTCGGCTCGCACCAGCAGAGCAGCCGTAGAATACTTAGTTTAATAGAATGGCTTGGCATTTNNNNNNNNNNCATGAGGTTGTTATTCTCAGAAGACTAATCACAGACA
.......Y.........YM....WS...Y....S...R....R....................................................C ....................G................
.......T.........TC....TG...C....G...A....A.................................. ,,,,,,,,,,,,,,,,,,,,g,,,,,,,,,,,,,,,,
.......T.........TC.....G...C....G...A....A.................................. ,,,,,,,,,,,,,,,,,,,,g,,,,,,,,,,,,,,,,
.......T.........TC....TG...C....G...A....A.................................. ,,,,,,,,,,,,,,,,,,,,g,,,,,,,,,,,,,,,,
.......T.........TC.....G...C....G...A....A.................................. ,,,,,,,,,,,,,,,,,,,,g,,,,,,,,,,,,,,,,
.......T.........TC.....G...C....G...A....A.................................. ,,,,,,,,,,,,,,,,,,,,g,,,,,,,,,,,,,,,,
.......T.........TC....TG...C....G...A....A..................................
.......T.........TC.....G...C....G...A....A..................................
.......T.........TC.....G...C....G...A....A..................................
.......T.........TC....TG...C....G...A....A..................................
.......T.........TC....TG...C....G...A....A..................................
..........C....................................................................................C
.......T.........TC.....G...C....G...A....A..................................
..........C....................................................................................C
.......T.........TC.....G...C....G...A....A..................................
......GT.........TC.....G...C....G...A....A..................................
.......T.........TC....TG...C....G...A....A..................................
There is a large amount of clipping going on here for the variant alleles likely why the mapping quality is low for them. You can also see that there are three different alleles present here. Press SHIFT+L to scroll further down the alignment. You can see that some of the polymorphism is also link to a cut site variant. All things that should be avoided.
Yet another filter that can be applied is whether or not their is a discrepancy in the properly paired status of for reads supporting reference or alternate alleles.
vcffilter -f "PAIRED > 0.05 & PAIREDR > 0.05 & PAIREDR / PAIRED < 1.75 & PAIREDR / PAIRED > 0.25 | PAIRED < 0.05 & PAIREDR < 0.05" -s DP3g95p5maf05.fil3.vcf > DP3g95p5maf05.fil4.vcf
Since de novo assembly is not perfect, some loci will only have unpaired reads mapping to them. This is not a problem. The problem occurs when all the reads supporting the reference allele are paired but not supporting the alternate allele. That is indicative of a problem.
mawk '!/#/' DP3g95p5maf05.fil4.vcf | wc -l
9166
Our loci count keeps dwindling, but our signal to noise ration keeps increasing. Let’s look at an example of what we filtered.
samtools tview BR_006-RG.bam reference.fasta -p E4407_L138
This output doesn’t paste well to the terminal, but you can see the clear discrepancy between mapping status and allele status. This could be indicative of cut site polymorphism or paralogs.
The next filter we will apply is to look at the ration of locus quality score to depth
Heng Li found some interesting results about how quality score and locus depth are related to each other in real and spurious variant calls
See his preprint here (http://arxiv.org/pdf/1404.0929.pdf)
Also see this great blog post about it here (http://bcb.io/2014/05/12/wgs-trio-variant-evaluation/) I REALLY recommend following that blog. Brad Chapman’s group is really good.
In short, with whole genome samples, it was found that high coverage can lead to inflated locus quality scores. Heng proposed that for read depths greater than the mean depth plus 2-3 times
the square root of mean depth that the quality score will be twice as large as the depth in real variants and below that value for false variants.
I actually found that this is a little too conservative for RADseq data, likely because the reads aren’t randomly distributed across contigs. I implement two filters based on this idea.
the first is removing any locus that has a quality score below 1/4 of the depth.
vcffilter -f "QUAL / DP > 0.25" DP3g95p5maf05.fil4.vcf > DP3g95p5maf05.fil5.vcf
The second is more complicated. The first step is to create a list of the depth of each locus
cut -f8 DP3g95p5maf05.fil5.vcf | grep -oe "DP=[0-9]*" | sed -s 's/DP=//g' > DP3g95p5maf05.fil5.DEPTH
The second step is to create a list of quality scores.
mawk '!/#/' DP3g95p5maf05.fil5.vcf | cut -f1,2,6 > DP3g95p5maf05.fil5.vcf.loci.qual
Next step is to calculate the mean depth
mawk '{ sum += $1; n++ } END { if (n > 0) print sum / n; }' DP3g95p5maf05.fil5.DEPTH
1952.82
Now the the mean plus 3X the square of the mean
python -c "print int(1952+3*(1952**0.5))"
2084
Next we paste the depth and quality files together and find the loci above the cutoff that do not have quality scores 2 times the depth
paste DP3g95p5maf05.fil5.vcf.loci.qual DP3g95p5maf05.fil5.DEPTH | mawk -v x=2084 '$4 > x' | mawk '$3 < 2 * $4' > DP3g95p5maf05.fil5.lowQDloci
Now we can remove those sites and recalculate the depth across loci with VCFtools
vcftools --vcf DP3g95p5maf05.fil5.vcf --site-depth --exclude-positions DP3g95p5maf05.fil5.lowQDloci --out DP3g95p5maf05.fil5
Now let’s take VCFtools output and cut it to only the depth scores
cut -f3 DP3g95p5maf05.fil5.ldepth > DP3g95p5maf05.fil5.site.depth
Now let’s calculate the average depth by dividing the above file by the number of individuals 31
mawk '!/D/' DP3g95p5maf05.fil5.site.depth | mawk -v x=31 '{print $1/x}' > meandepthpersite
Let’s plot the data as a histogram
gnuplot << \EOF
set terminal dumb size 120, 30
set autoscale
set xrange [10:150]
unset label
set title "Histogram of mean depth per site"
set ylabel "Number of Occurrences"
set xlabel "Mean Depth"
binwidth=1
bin(x,width)=width*floor(x/width) + binwidth/2.0
set xtics 5
plot 'meandepthpersite' using (bin($1,binwidth)):(1.0) smooth freq with boxes
pause -1
EOF
Histogram of mean depth per site
Number of Occurrences
250 ++--+---+---+---+--+---+---+---+---+---+---+---+---+---+--+---+---+---+---+---+---+---+---+--+---+---+---+--++
+ + + + + + + + + + + + + +'meandepthpersite' using (bin($1,binwidth)):(1.0)+****** +
| ** ** |
| *** ** |
| ******** |
200 ++ ********* * ++
| ** ********** * |
| ************* ** * |
| **************** * * |
| ******************* ** |
150 ++ ******************* *** ++
| ************************* |
| *************************** |
| ******************************* |
100 ++ ********************************** ++
| ********************************** |
| ************************************** * |
| ************************************** * ** |
| ********************************************* |
50 ++ ************************************************ ++
| ************************************************ ** ** |
| *********************************************************** ** |
| **************************************************************** ** ** ** |
+ + ************************************************************************************** ******+********
0 ++--+---******************************************************************************************************
10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125130 135 140 145 150
Mean Depth
Loci that have high mean depth are indicative of either paralogs or multicopy loci. Either way we want to remove them. Here, I’d remove all loci above a mean depth of 102.5. Now we can combine both filters to produce another VCF file
vcftools --vcf DP3g95p5maf05.fil5.vcf --recode-INFO-all --out DP3g95p5maf05.FIL --max-meanDP 102.5 --exclude-positions DP3g95p5maf05.fil5.lowQDloci --recode
In the end, VCFtools kept 8417 out of a possible 9164 Sites. BTW, I’ve also written a script to automate the filterings steps described in steps 23-44. It’s called dDocent_filters. It will go through the filtering steps and recode a log file for you for each of the steps, including the depth histogram.
curl -L -O https://github.com/jpuritz/dDocent/raw/master/scripts/dDocent_filters
chmod +x dDocent_filters
./dDocent_filters
This script will automatically filter a FreeBayes generated VCF file using criteria related to site depth,
quality versus depth, strand representation, allelic balance at heterzygous individuals, and paired read representation.
The script assumes that loci and individuals with low call rates (or depth) have already been removed.
Contact Jon Puritz (jpuritz@gmail.com) for questions and see script comments for more details on particular filters
Usage is sh FB_filters.sh VCF_file Output_prefix
The next filter to apply is HWE. Heng Li also found that HWE is another excellent filter to remove erroneous variant calls. We don’t want to apply it across the board, since population structure will create departures from HWE as well. We need to apply this by population. I’ve included a perl script written by Chris Hollenbeck, one of the PhD student’s in my current lab that will do this for us.
curl -L -O https://github.com/jpuritz/dDocent/raw/master/scripts/filter_hwe_by_pop.pl
chmod +x filter_hwe_by_pop.pl
./filter_hwe_by_pop.pl
Usage:
filter_hwe_by_pop.pl -v <vcffile> -p <popmap> [options]
Options: -v <vcffile> input vcf file -p <popmap> tab-separated file of
samples and population designations -h [hwe] minimum Hardy-Weinberg
p-value cutoff for SNPs -c [cutoff] proportion of all populations that a
locus can be below HWE cutoff without being filtered -o [out] name of
outfile
Options:
-v, --vcffile
VCF input file
-p, --popmap
File with names of individuals and population designations, one
per line
-h, --hwe
Minimum cutoff for Hardy-Weinberg p-value (for test as
implemented in vcftools) [Default: 0.001]
-c, --cutoff
Proportion of all populations that a locus can be below HWE
cutoff without being filtered. For example, choosing 0.5 will
filter SNPs that are below the p-value threshold in 50% or more
of the populations. [Default: 0.25]
-o, --out
Name of outfile, by vcftools conventions (will be named
X.recode.vcf)
Let’s filter our SNPs by population specific HWE First, we need to convert our variant calls to SNPs To do this we will use another command from vcflib called vcfallelicprimatives
vcfallelicprimitives DP3g95p5maf05.FIL.recode.vcf --keep-info --keep-geno > DP3g95p5maf05.prim.vcf
This will decompose complex variant calls into phased SNP and INDEL genotypes and keep the INFO flags for loci and genotypes. Next, we can feed this VCF file into VCFtools to remove indels.
vcftools --vcf DP3g95p5maf05.prim.vcf --remove-indels --recode --recode-INFO-all --out SNP.DP3g95p5maf05
We now have 8379 SNP calls in our new VCF. Now, let’s apply the HWE filter
./filter_hwe_by_pop.pl -v SNP.DP3g95p5maf05.recode.vcf -p popmap -o SNP.DP3g95p5maf05.HWE -h 0.01
Processing population: BR (20 inds)
Processing population: WL (20 inds)
Outputting results of HWE test for filtered loci to 'filtered.hwe'
Kept 8176 of a possible 8379 loci (filtered 203 loci)
Note, I would not normally use such a high -h value. It’s purely for this example. Typically, errors would have a low p-value and would be present in many populations.
We have now created a thoroughly filtered VCF, and we should have confidence in these SNP calls.
However, our lab is currently developing one more script, called rad_haplotyper.
This tool takes a VCF file of SNPs and will parse through BAM files looking to link SNPs into haplotypes along paired reads.
curl -L -O https://raw.githubusercontent.com/chollenbeck/rad_haplotyper/master/rad_haplotyper.pl
chmod +x rad_haplotyper.pl
Note, this script requires several Perl libraries. See the README here It has a lot of options, let’s take a look
perl rad_haplotyper.pl
Usage:
perl rad_haplotyper.pl -v <vcffile> -r <reference> [options]
Options: -v <vcffile> input vcf file
-r <reference> reference genome
-s [samples] optionally specify an individual sample to be haplotyped
-u [snp_cutoff] remove loci with more than a specified number of SNPs
-h [hap_cutoff] remove loci with more than a specified number of haplotypes relative to SNPs
-m [miss_cutoff] cutoff for proportion of missing data for loci to be included in the output
-mp [max_paralog_inds] cutoff for excluding possible paralogs
-ml [max_low_cov_inds] cutoff for excluding loci with low coverage or genotyping errors
-d [depth] sampling depth used by the algorithm to build haplotypes
-g [genepop] genepop file for population output
-p [popmap] population map for organizing Genepop file
-t [tsvfile] tsv file for linkage map output
-a [imafile] IMa file output
-p1 [parent1] first parent in the mapping cross
-p2 [parent2] second parent in the mapping cross
-x [threads] number of threads to use for the analysis
-n use indels
-e debug
Options:
-v, --vcffile
VCF input file
-r, --reference
Reference genome (FASTA format)
-s, --samples
Individual samples to use in the analysis - can be used multiple
times for multiple individuals [Default: All]
-u, --cutoff
Excludes loci with more than the specified number of SNPs
[Default: No filter]
-h, --hap_count
Excludes loci with more than the specified number of haplotypes
relative to number of SNPs. Excluding forces other than mutation
(i.e. recombination) the maximum number of haplotypes should be
one more than the number of SNPs at the locus. The value
provided is the number of haplotypes allowed in excess of the
number of SNPs, which allows that mechanisms other than mutation
may have influenced the number of haplotypes in the population.
[Default: 100]
-x, --threads
Run in parallel across individuals with a specified number of
threads
-n, --indels
Includes indels that are the only polymorphism at the locus
(tag)
-d, --depth
Specify a depth of sampling reads for building haplotypes
[Default: 20]
-m, --miss_cutoff
Proportion of missing data cutoff for removing loci from the
final output. For example, to keep only loci with successful
haplotype builds in 95% of individuals, enter 0.95. [Default:
0.9]
-mp, --max_paralog_inds
Count cutoff for removing loci that are possible paralogs from
the final output. The value is the maximum allowable number of
individuals with more than the expected number of haplotypes
[Default: No filter]
-ml, --max_low_cov_inds
Count cutoff for removing loci with low coverage or genotyping
errors from the final output. The value is the maximum allowable
number of individuals with less than the expected number of
haplotypes [Default: No filter]
-g, --genepop
Writes a genepop file using haplotypes. Must provide the name of
the genepop file.
-a, --ima
Writes a IMa file using haplotypes. Must provide the name of the
IMa file.
-p, --popmap
Tab-separated file of individuals and their population
designation, one per line (required for Genepop output)
-t, --tsvfile
Writes a tsv file using haplotypes - for mapping crosses only.
Must provide the name of the tsv file.
-p1, --parent1
Parent 1 of the mapping cross (must be specified if writing a
tsv file)
-p2, --parent2
Parent 2 of the mapping cross (must be specified if writing a
tsv file)
-e, --debug
Output extra logs for debugging purposes
We don’t have enough time to go into depth with all these options and this tool is still under development.
It also take some substantial resources to run. I will simulate running this for you.
#rad_haplotyper.pl -v SNP.DP3g95p5maf05.HWE.recode.vcf -x 40 -mp 1 -u 20 -ml 4 -n -r reference.fasta
Note, this will not actually run. It needs all the BAM files to proceed.
Here’s the simulated output
Removed 0 loci (0 SNPs) with more than 20 SNPs at a locus
Building haplotypes for BR_024
Building haplotypes for BR_028
Building haplotypes for WL_054
Building haplotypes for BR_016
Building haplotypes for BR_009
Building haplotypes for BR_006
Building haplotypes for BR_041
Building haplotypes for BR_040
Building haplotypes for BR_046
Building haplotypes for BR_031
Building haplotypes for BR_025
Building haplotypes for BR_002
Building haplotypes for WL_058
Building haplotypes for WL_057
Building haplotypes for WL_061
Building haplotypes for WL_069
Building haplotypes for WL_070
Building haplotypes for BR_048
Building haplotypes for WL_031
Building haplotypes for WL_056
Building haplotypes for BR_047
Building haplotypes for WL_079
Building haplotypes for WL_080
Building haplotypes for WL_032
Building haplotypes for WL_071
Building haplotypes for WL_081
Building haplotypes for BR_004
Building haplotypes for BR_021
Building haplotypes for BR_015
Building haplotypes for BR_043
Building haplotypes for WL_066
Filtered 26 loci below missing data cutoff
Filtered 66 possible paralogs
Filtered 17 loci with low coverage or genotyping errors
Filtered 0 loci with an excess of haplotypes
The script found another 109 loci to remove from our file. Besides this output to the terminal, the script outputs a file called stats.out
head stats.out
Locus Sites Haplotypes Inds_Haplotyped Total_Inds Prop_Haplotyped Status Poss_Paralog Low_Cov/Geno_Err Miss_Geno Comment
E10001_L101 1 2 30 31 0.968 PASSED 0 0 1
E10003_L101 7 9 30 31 0.968 PASSED 1 0 0
E10004_L101 - - - - - FILTERED0 0 1 Complex
E10008_L101 1 2 30 31 0.968 PASSED 0 0 1
E10013_L142 3 6 30 31 0.968 PASSED 0 0 1
E10014_L117 2 3 31 31 1.000 PASSED 0 0 0
E10024_L101 1 2 31 31 1.000 PASSED 0 0 0
E10029_L101 1 2 31 31 1.000 PASSED 0 0 0
We can use this file to create a list of loci to filter
grep FILTERED stats.out | mawk '!/Complex/' | cut -f1 > loci.to.remove
Now that we have the list we can parse through the VCF file and remove the bad RAD loci I’ve made a simple script to do this remove.bad.hap.loci.sh
curl -L -O https://github.com/jpuritz/dDocent/raw/master/scripts/remove.bad.hap.loci.sh
chmod +x remove.bad.hap.loci.sh
./remove.bad.hap.loci.sh loci.to.remove SNP.DP3g95p5maf05.HWE.recode.vcf
This produces a FINAL FINAL FINAL filtered VCF file SNP.DP3g95p5maf05.HWE.filtered.vcf
mawk '!/#/' SNP.DP3g95p5maf05.HWE.filtered.vcf | wc -l
We’re left with 7,666 SNPs! How many possible errors?
./ErrorCount.sh SNP.DP3g95p5maf05.HWE.filtered.vcf
This script counts the number of potential genotyping errors due to low read depth
It report a low range, based on a 50% binomial probability of observing the second allele in a heterozygote and a high range based on a 25% probability.
Potential genotyping errors from genotypes from only 1 read range from 0 to 0.0
Potential genotyping errors from genotypes from only 2 reads range from 0 to 0.0
Potential genotyping errors from genotypes from only 3 reads range from 302 to 1014.72
Potential genotyping errors from genotypes from only 4 reads range from 162 to 822.232
Potential genotyping errors from genotypes from only 5 reads range from 88 to 669
31 number of individuals and 7666 equals 237646 total genotypes
Total genotypes not counting missing data 237081
Total potential error rate is between 0.00232831816974 and 0.0105700245908
SCORCHED EARTH SCENARIO
WHAT IF ALL LOW DEPTH HOMOZYGOTE GENOTYPES ARE ERRORS?????
The total SCORCHED EARTH error rate is 0.0330857386294.
Congrats! You’ve finished the Filtering Tutorial