VCF activity

The VCF format, filtering, and first analyses

Julia M.I. Barth, 22 January 2018

Background and Aim

With the advances of sequencing technology, reading whole genomes and re-sequencing population samples has become incredibly cheap over the last years. This has led to a shift from single-marker analyses to whole-genome studies, which at the same time caused a change from graphical user interfaces (GUI) software to command-line based programs. Such programs are much faster and more suitable for handling very large data files, but it might take some time to get used to ‘not seeing your data and analyses’.

Learning goals

• Apply the UNIX skills you learned in the previous exercises.
• Learn about the VCF format and how to handle and manipulate VCF files.
• Perform a PCA with genomic data.

Requirements
bcftools v1.6
VCFtools v0.1.14
EIGENSOFT smartpca v6.0.1
R / Rstudio (for plotting)

Example data

The example data is a VCF file containing biallelic SNP data of several individuals in two populations. The data was produced using Illumina sequencing, mapped against a reference genome, and already filtered to contain only high quality SNPs (e.g. high read depth and genotype quality).

How to do this activity
Questions or tasks are indicated with Q.
All text in red underlined with gray color indicates commands or file names that you can type or copy to the terminal. Text in black underlined with gray color refers to terminal outputs.
Try to work independently, use the provided manuals, discuss with your neighbor. If you get stuck, check the answer-box:

There are always alternative options and solutions. Whether you prefer cat over less, or whether you found another solution to the task, go ahead and use your way of doing it!

Table of contents

1. The VCF format
2. VCFtools
• Filter missing data
• Filter on minor allele frequency
• Exclude chromosome 7 and linked sites
• Export the VCF to PLINK format
3. Perform a PCA

1) The VCF format

A VCF (Variant Call Format) file is a standardized text file format that is used to store genetic variation calls such as SNPs or insertions/deletions of multiple individuals. The full format specifications can be found here.
In the following first part of the exercise, you will learn how the information in a VCF is stored, and how you can inspect it.

Q Open the terminal and navigate to the analysis folder including the example data file: ~/workshop_materials/01_unix_intro/vcf_activity
What is the file size of data.vcf.gz? How many bytes and how many MB does it have?

Q What is written on the first ten lines of the VCF file ‘data.vcf.gz’ file?

The first line is already telling you that this is indeed a VCF file:
##fileformat=VCFv4.1

You can also see some more information about the VCF file, e.g., when it was created and which reference file was used:
##fileDate=20160104
##reference=file:///work/jobs/11187331/nuc_reference.fasta

These lines, starting with ##, are part of the ‘VCF header’, which stores useful meta-information about the data. If you use zless you can scroll down to see more of the header and also the genotype data, which make up the majority of the file.

There can be a lot of information in the header. Different meta-information lines, e.g. ‘INFO’, ‘FILTER’, ‘FORMAT’, ‘ALT’ describe types of information that are given for every variant in the VCF file. A few examples:
##INFO<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered"> :The read depth.
##FORMAT<ID=GT,Number=1,Type=String,Description="Genotype"> :The genotypes coded as allele values separated by /. The allele values are 0 for the reference allele and 1 for the alternative allele. Missing data is denoted as ..
##contig<ID=LG01,length=28303952> : A list of chromosomes, contigs, or scaffolds.

The very last line of the header specifies the column names for the genotype data. The first 8 columns are mandatory, follow by ‘FORMAT’ and the sample IDs:
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT BS001 BS002 BS003 BS004 BS005 BS006 BS007 BS008 BS009 BS010 NS001 NS002 NS003 NS004 NS005 NS006 NS007 NS008 NS009 NS010

You may want to edit these header-lines to add changes that you did to your VCF file, or you received a VCF file from colleagues and wish to inspect the header in more detail to see what they did.

Q Find a command in ‘bcftools’ that allows you to separate the header from the large genotype data. Save the header in a seperate text file named data.header.txt.
To find the command have a look at the bcftools manual, or just type bcftools --help on the command line to see available options.

Q What’s the file size of data.header.txt, compared to the whole VCF file?

Q Open the header in a text editor, change the name of the individual ‘BS004’ to ‘BS014’, and save the new header as: data.header.ed.txt. Replace (‘re-head’) the VCF file ‘data.vcf.gz’ with the new edited header and save the re-headed VCF as data.ed.vcf.gz. Make sure your changes took effect and your file is still valid by inspecting the individual names and the file size of ‘data.ed.vcf.gz’.

In the genotype part of the VCF file, one line corresponds to one variant. Each line starts with the chromosome ID ‘CHROM’, the position ‘POS’, and the variant ‘ID’, and is followed by the bases of the two alleles (‘REF’ for the reference allele, and ‘ALT’ for the alternative allele). The 6th column (‘QUAL’) is the phred-scaled quality score, a measure for how likely it is that this site is indeed a variant. In the next column, ‘PASS’ tells us that the variant passed all filters. In the ‘INFO’ column, information related to the header info tags is given (e.g., AC is the allele count in genotypes for the ‘ALT’ allele, AF the allele frequency of the ‘ALT’ allele, and AN is the total number of alleles, etc.). The ‘FORMAT’ column specifies the fields and order of the data that is given with each individual, and the remaining columns list each individual with the values corresponding to the ‘FORMAT’ field. The first ‘FORMAT’ field is always the genotype (GT), which in our example is followed by AD (allelic depths for the ref and alt alleles), DP (read depth), GQ (genotype quality), and PL (phred-scaled genotype likelihood):

LG03 4837 . C A 559.22 PASS AC=11;AF=0.009016;AN=1220;BaseQRankSum=0.431;ClippingRankSum=0;DP=7189;ExcessHet=4.6637;FS=43.708;InbreedingCoeff=-0.0285;MLEAC=11;MLEAF=0.009016;MQ=55.93;MQRankSum=-2.093;QD=3.37;ReadPosRankSum=0.889;SOR=3.89 GT:AD:DP:GQ:PGT:PID:PL 0/1:12,3:15:90:0|1:4829_C_T:90,0,495 0/0:18,0:18:54:.:.:0,54,702 0/0:19,0:19:54:.:.:0,54,810 0/0:13,0:13:39:.:.:0,39,533 ./.:15,0:15:3:.:.:0,3,523 ./.:17,0:17:9:.:.:0,9,603 0/0:17,0:17:51:.:.:0,51,649 0/0:13,0:13:39:.:.:0,39,501 0/0:9,0:9:27:.:.:0,27,343 0/0:10,0:10:30:.:.:0,30,395 0/0:7,0:7:21:.:.:0,21,288 ./.:8,0:8:0:.:.:0,0,255 0/0:12,0:12:33:.:.:0,33,495 0/0:14,0:14:42:.:.:0,42,479 0/0:14,0:14:39:.:.:0,39,585 ./.:6,0:6:18:.:.:0,18,247 ./.:6,0:6:15:.:.:0,15,225 0/0:21,0:21:60:.:.:0,60,900 0/0:15,0:15:42:.:.:0,42,630 0/0:11,0:11:30:.:.:0,30,450

Q Examine the first variant in the VCF file and answer the following questions:
A) What kind of variant is it (e.g. SNP, insertion/deletion (indel))? What would you expect if it was an indel?

B) What genotype (GT) has the first individual BS001 at the second SNP (LG03 8159)?

C) What is the read depth (DP) for the first three individuals BS001, BS002, and BS003 at the second SNP (LG03 8159)?

2) VCFtools

Now that we know the VCF format, we can start to adjust our file for certain analyses. As mentioned before, our VCF file is already filtered for quality, read depth, etc. (find more on variant calling and quality filtering on the Genomics workshop page).
However, you may want to use different thresholds for minor allele frequency or the proportion of missing data for specific analyses. Or you might want to filter out certain variants or chromosomes. The software VCFtools has various functions to manipulate, inspect, filter, and merge VCF files. VCFtools can also output statistics such as heterozygosity, allele frequencies, or Fst.

2.1) Filter missing data

Individuals and/or variants with a lot of missing data (often due to low sequence coverage) are not very informative and may interfere with certain analyses. Excluding individuals and/or variants with too much missing data as a first step also reduces file size and speeds up analyses.

Q How much missing data has each individual (missingness per-individual) in ‘data.ed.vcf.gz’?
Tip: Have a look at the vcftools manual website or man vcftools.

Q Using the command line, sort the individuals according to their amount of missing data and make a note which individuals have 25% or more missing genotypes.

Along with excluding individuals with few data, we also want to exclude sites for which 20% or more of the genotypes are missing.

Q Browse the vcftools manual and make a note of the command used to exclude sites on the basis of the proportion of missing data and the command used to remove individuals.

Q Using a one-line command in VCFtools, remove the two individuals with 25% missing data from the ‘data_edited.vcf.gz’ file and at the same time exclude sites with more than 20% missing genotypes. Make sure to save the output file (named data.ed.imiss25.miss20.vcf.gz) directly as a compressed VCF file using bgzip (output to standard out then pipe to bgzip). Include all VCF header info (--recode-INFO-all).
Tip: Get inspired by the examples at the beginning of the vcftools manual.
(takes about 2 minutes)

Verify that your command worked by checking the individual IDs, having a look at the new file size and the log-file entry, and by comparing the number of variants in the two VCF files:

Q How much smaller is the new file?
How many variants were excluded due to the applied missingness filters?
Do the numbers given in the log file and your calculations on the original file and the new file match?
Tip: Use a cobination of bcftools (to suppress the header) and word-count (wc) to count the variants.

2.2) Filter on minor allele frequency

Filtering on minor allele frequency (MAF; the frequency of the allele that appears less often) can be done for multiple reasons. One reason is that in samples sequenced with low coverage, variants with very low MAF frequently present false calls, and can be obstructive by adding noise to the data. If the data will be used for a genome-wide association study (GWAS), usually a rather high MAF threshold is applied, as it requires very strong statistical power to make meaningful statements about very rare alleles. When interested in population structure with many samples, alleles present in only one individual do not add information for the inference of structure.
Therefore, the optimal MAF filtering threshold depends on the type of data you have, the analysis you are planning to do, the data quality, and the sample size.

Q Ignoring the quality, what would you think would be an appropriate MAF threshold to filter ‘data.ed.imiss25.miss20.vcf.gz’ for population structure analysis?

If you have no idea what an appropriate threshold could be, read hint 1:

If you’re still unsure as of which threshold to choose, read hint 2.

Q If you’ve decided on a threshold, go ahead and apply it. Make sure to save the output file (named ‘data.ed.imiss25.miss20.maf[your_threshold].noLG07.vcf.gz’) directly as a compressed VCF file using bgzip and include all VCF header info as above.

An appropriate MAF threshold is given in the answer below:

Q Go ahead and apply the proposed MAF threshold. Make sure to save the output file (named data.ed.imiss25.miss20.maf0.1.vcf.gz) directly as a compressed VCF file using bgzip and include all VCF header info as above.

Q Make sure your command worked by inspecting the new file and the log file.

2.3) Exclude chromosome 7 and linked sites

VCFtools offers various options to exclude or include certain parts of your file in the analysis. Have a look at the ‘Site filtering” section in the manual.

For population structure analysis we want to exclude sites that are very close to each other and thus likely share the same information (they are in linkage disequilibrium). In addition, we also exclude chromosome 7 as we know that a large inversion on this chromosome causes many SNPs to be linked.

Q Make a new VCF without chromosome 7, and in which sites closer than 10 bp from each other are excluded. As above, save the output file (named data.ed.imiss25.miss20.maf0.1.LD.vcf.gz) directly as a compressed VCF file using bgzip and include all VCF header info (takes about 1 min).

Have a look at the log-file; how many variants were excluded this time? If you like you can again compare the numbers to the output file using wc like above. Check whether chromosome 7 was really excluded from the file.

2.4) Export the VCF to PLINK format
As a last step, we’ll use VCFtools to export the VCF into a format that we can use to analyze population structure by performing a principal component analysis (PCA).
The software PLINK is a comprehensive genome analysis toolset with an extensive list of functions. It was originally developed for human data (hence it uses all those human terms like “family” and “father”), but the new PLINK 1.9 can also be used with genomic data of non-model organisms.
Due to time constraints PLINK is not included in this activity, but if you’re interested to learn more about PLINK and how to filter on linkage disequilibrium (LD) you can follow the WPSG2016 vcftools/PLINK activity.

Q Use VCFtools to export the final VCF (data.ed.imiss25.miss20.maf0.1.LD.vcf.gz) including the chromosome information given in the analyses folder chrom_map.txt to PLINK format. Use data.ed.imiss25.miss20.maf0.1.LD as the file prefix for the new file.

Type ls -l to see the new files. You should have two new data files (‘.ped’ and ‘.map’, and a log file).

3) Perform a PCA

The smartpca program implemented in the software EIGENSOFT uses a combination of Principal Component Analysis (PCA; a technique used to emphasize variation and bring out strong patterns in a dataset) and new statistics designed specifically for genomic data. The methods are described in detail in an excellent and accessible manuscript by the authors of the program.

The EIGENSOFT smartpca needs four input files: The ‘.ped’ (containing the genotype information) and ‘.map’ (containing the chromosome/position information) files that we exported above, a ‘.pedind’ file containing individual information, and a ‘.par’ file with the specified parameters for analysis.

Q To generate the ‘.pedind’ file, cut the first six columns of the ‘.ped’ file and save them in a new file called ‘tmp’.

Have a look at the ‘tmp’ file. It contains the family IDs and the individual IDs in the first two columns (in our case this is the same), followed by zeros that indicate ‘no information’ for the remaining columns (fatherID, motherID, sex, phenotype). In addition, we need a column to specify which individual belongs to which population.

Q Add a tab-separated 7th column specifying that BS001-BS010 belong to the population ‘BS’ and NS001-NS010 belong to the population ‘NS’ to the ‘tmp’ file and save the new file as data.ed.imiss25.miss20.maf0.1.LD.pedind. Remove the ‘tmp’ file.

A template parameter file (template.par) and a smartpca manual page (eigensoft_smartpca_readme.txt) can be found in the analysis folder.

Q Modify the template.par to include the respective file names, save it as data.ed.imiss25.miss20.maf0.1.LD.par and run the smartpca with the syntax ‘smartpca -p parfile > pca01_out’. Plot the PCA using RStudio in the browser using your own code or the provided R cheat-sheet 20180122_vcf_plot.R.

Another smartpca example on blue and yellow crater lake cichlids can be found here.