FileFormats.VCFtools.PLINK

(Compiled for the WPSG2016 by Julia M.I. Barth)

Table of contents

• Background and Aim
• Getting Started
• Exercise 1 – The VCF format
• Exercise 2 – VCFtools
• • Filtering for missing data
  • Exclude the mitochondrial scaffold and very close sites
  • Export the VCF to PLINK format
• Exercise 3 – PLINK
• • Input formats
  • Filtering for Minor allele frequency (MAF)
  • Filtering for Linkage disequilibrium
  • Export formats
• Survey

Background and Aim

Please find the intro-slides here.

With the advances of sequencing technology, reading whole genomes and re-sequencing population samples has become incredibly cheap over the last years. In the field of population genetics, this progress has led to a shift from single-marker analyses to whole-genome studies, which at the same time caused a change from analysis software using graphical user interfaces 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”.

During this workshop we will focus on the analysis of genomic data sets, rather than on their generation. However, in order to prepare your data for any biological interpretation, some initial adjustments are almost always necessary. This can be a simple change of file format, or it can be that you want to remove that one individual that has a very low sequencing depth. In addition, you might also want to include different filtering steps for different analyses, to specify varying thresholds for the proportion of missing data, or for the minor-allele frequency.

During this tutorial we will primarily use the command-line based software VCFtools and PLINK, as well as some of the basic UNIX commands you learned yesterday, and we will filter and convert population data from a raw VCF file so we can use it for analyses of e.g. population structure.

Instead of working with a data set that is mapped against a well-assembled genome (like human data), we will work with data from a non-model organism where a high-quality genome is not yet available and which is thus arranged in thousands of smaller scaffolds instead of a few chromosomes. This adds a bit of complexity to handling of the data, but is likely similar to what most of you will be doing in your own work. The example data has already gone through quality filtering (e.g. it was filtered to contain only variants that have a read depth (DP) between 3 and 30, and a minimum genotype quality (GQ) of 20) and filtering for bi-allelic markers (SNPs) only. These SNPs have been called using the variant calling software GATK. If you’re interested in how to actually get to such a data set, you might have a look at the exercises of the Genomics workshop from last week.

Getting Started

All computer programs needed for this activity are installed on the Amazon server. These are:
bcftools v1.3
VCFtools v0.1.14
PLINK v1.9p – special build of PLINK v1.9 for high contig size (request from authors, or ask Julia)
R / Rstudio / gnuplot (for plotting only, whichever you prefer)

The example data file can be found here: ~/wpsg_2016/activities/vcftools_plink
Or can be downloaded here.

How to do the exercise:
The main purpose of this activity is to get you familiar with the software listed above, to work on the command line, and to do so independently.
You could easily go through the exercise and just copy and paste the code, but that way you would most likely only wear out your Ctrl, C, and V keys without getting the most out of it. Instead, use the help options and the manual if you don’t know how to do the task.
With every question or task a link to the software manual (and very often the required section of the manual) is provided. Use those links and try to find the solution BEFORE looking into the answer box!
If you get stuck, check the answer and try to understand it. If your command was not correct, you can (and should) of course copy and paste long commands rather then type them letter-by-letter. All text underlined with gray color indicates commands that you can type or copy to the terminal.
There are always alternative options and solutions for UNIX code. E.g., if you prefer cat to less, go ahead and use cat. If you’ve found another solution that works as well: Great! Just go ahead and use your solution!

If you have questions, bother the instructors! They are here to help you! 😉

And now, let’s start!

Exercise 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, we will learn how the information in a VCF is stored, and how we can inspect it.

Open the terminal and use change directory cd to navigate to the analysis folder including the example data file.
What is the file size of the example data file data.vcf.gz? How many bytes, or MB does it have?

What is in the file?

Use head data.vcf.gz to see the first ten lines of the file.

Oups, only gibberish!
It’s a compressed file (indicated by the ending .gz). Compressing VCF files with gzip (or bgzip and indexing it with tabix) is the standard way in which VCF files are stored.
We could decompress the file, which would take some time and inflate the file size to several GB, so we’re not going to do that. Instead, we can use another tool to display the file content: “bcftools”, a software for variant calling and manipulating of files in VCF format.
Have a look at the manual here: bcftools manual, or just type bcftools on the command line to see available options.

Do you find a command that would allow us to have a look into (“view”) the compressed VCF?

Ok, lets test that one:
bcftools view data.vcf.gz

Ah…. what’s happening?? It’s printing the whole file on the screen!

Press the “panic-button” ctrl-c (control c) to abort the command.

Can you come up with something that we can combine with bcftools view data.vcf.gz to have a sneak-peek into the file without printing everything on the screen? Try something you learned yesterday!

With these commands we can already see that we indeed are dealing with a VCF file:
##fileformat=VCFv4.1

We also see some more information about the VCF file, e.g. when it was created, which reference file was used, etc.
##fileDate=20140121
##reference=file:///work/jobs/11187331/reference.scf.fasta

These lines are part of the “VCF header”, which stores useful meta-information about the data and always starts with ##. 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.
Can you find a command in “bcftools” on how to separate the header from the large genotype data in order to open and edit it with a text editor?
Type bcftools view to see the options.

Found something? Try it!

Uh – it’s only the header, but still a lot of information! Save the header to a file called data.header.txt instead of printing it on the screen so you can have a look at it with a text editor.

What’s the file size of data.header.txt that we extracted from the VCF, compared to the whole VCF file?

If you like you can download the header to look at it with a text editor, or you can examine it right here on the command line. How would you do that?

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; some reads may have been filtered”> – The read depth.
##FORMAT<ID=GT,Number=1,Type=String,Description=”Genotype”> The genotypes, which, as we will see below are encoded 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=scf7180003948298,length=1427> A list of 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 IND01 IND02 IND03 IND04 IND05 AFR01 AFR02 AFR03 AFR04 AFR05

After this line the genotypes are specified, but this is not part of the header anymore.
Any idea how we could have a sneak-peek into the genotype data without printing everything on the screen?

Tip: We know that the last line of the header starts with #CHROM, so we could grab CHROM and the following line by specifying grep -A1.

This is the last header line, and the first genotype line:

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT IND01 IND02 IND03 IND04 IND05 AFR01 AFR02 AFR03 AFR04 AFR05
scf7180003948298 203 . T C 2899.83 PASS AC=3;AF=0.793;AN=4;BaseQRankSum=0.731;ClippingRankSum=0.731;DP=151;FS=6.8;GQ_MEAN=16;GQ_STDDEV=21.53;InbreedingCoeff=0.2421;MLEAC=47;MLEAF=0.81;MQ=60;MQ0=0;MQRankSum=0.731;NCC=11;QD=34.24;ReadPosRankSum=0.804 GT:AD:DP:GQ:PL 1/1:0,8:8:24:327,24,0 ./.:0,4:4:12:167,12,0 ./.:0,0:0:.:. ./.:0,0:0:.:. ./.:0,5:5:15:199,15,0 ./.:0,4:4:12:167,12,0 0/1:4,1:5:27:27,0,103 ./.:0,2:2:6:84,6,0 ./.:0,0:0:.:. ./.:5,0:5:15:0,15,186

In the genotype part of the VCF file, one line corresponds to one variant. It starts with the chromosome ID, the position, and an ID (missing here (.)) of the variant, and is followed by the bases of the two alleles (REF for the reference allele, and ALT for the alternative allele). The 6th column 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).

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

B) What genotype (GT) does the first individual IND01 have?

C) I mentioned the field “read depth” (DP) above. Can you tell me the read depth for the first three individuals IND01-IND03?

______________________________________________________________________________________________

Exercise 2 – VCFtools

Now that we know the VCF format, we can start to adjust our file for certain analyses. As I 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 is a package that has various functions to manipulate, inspect, filter, and merge VCF files. VCFtools can also output statistics such as heterozygosity, allele frequencies, or Fst.

Go to the terminal and type man vcftools. Have a look at the “SYNOPSIS” to get to know the general commands needed to run VCFtools. If this appears cryptic, have a look at the “EXAMPLES”. Once finished, exit the man pages by typing q.

Then, let us begin to work on the data.vcf.gz file using VCFtools!

2.1) Filtering for missing data

Individuals and/or variants with a lot of missing data (often due to low ‘sequence coverage’) are not very informative and usually are not included in further analyses. In order to exclude such non-informative data and at the same time also reduce file size and speed-up analyses, the first thing we’ll do is to exclude all individuals and variants that have too much missing data.

Let’s start with checking how much missing data our individuals have.
Have a look at the vcftools manual website or the man vcftools and find a command to be used to get an estimate of the missingness per individual.
Try it!

Have a look at the text file output data.imiss! Can you spot the odd one out?

IND04 has much more missing data then all other individuals, so we’ll exclude it from further analysis.
Can you find a command that would allow us to do so? Have a look into the vcftools manual or the man vcftools again!

Keep this command in mind, we’ll use it at the next step!!

Along with filtering individuals with few data, we also want to exclude variants with too much missing data. How much missing data we allow will depend on the analyses that we’re planing to do, but for a lot of analyses (such as population structure) we don’t want variants for which half or more of the data is missing. By excluding such variants, we can fairly reduce the file size and increase the speed of follow-up analyses!

Have another look into the VCFtools manual or the man vcftools to see if there is a command that would allow us to exclude all variants with 50% or more missing data.

Now, combine the two commands to remove the one individual with a lot of missing data and to exclude sites with more than 50% missing data, and write the data to a new compressed VCF file called data.noIND04.miss0.5.vcf.gz. Make sure to also include all the INFO fields and values in the new VCF!

Tip: Look at the examples at the beginning of the VCFtools manual to find a solution on how to output a new VCF file to standard out (stdout) and then compress it with gzip.

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 new reduced data.noIND04.miss0.5.vcf.gz to the original data.vcf.gz VCF file.

Answer the following questions:
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: You can use bcftools to check the individual IDs! Use view and grep as in Exercise 1, or have a look at the bcftools manual to see whether there is a more direct command to output the sample IDs. Check the file size like you did in Exercise 1, and use word-count (wc) to count the variants.

2.2) Exclude the mitochondrial scaffold and very close 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 our downstream analysis we want to exclude the mitochondrial scaffold from the VCF. We also want to exclude sites that are very close to each other and thus likely share the same information.

Make a new VCF without the mitochondrial scaffold “scf7180003949713”, and in which sites closer then 10 bp from each other are excluded. Use the file data.noIND04.miss0.5.vcf.gz as input and name the new file data.noIND04.miss0.5.noMT.thin10.vcf.gz. Compress the output and include all INFO fields as before. Check the manual for the necessary commands!

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

2.3) Export the VCF to PLINK format

As a last step, we’ll use VCFtools to export the VCF into a format that can be read by the software PLINK.

Since we have thousands of scaffolds rather than only a few chromosomes, we should specify a map that defines every scaffold as one “chromosome”, otherwise the scaffold ID information will be lost during the conversion to the PLINK format. This is optional, but needed if you would like to perform analysis per scaffold as e.g. filtering for Linkage disequilibrium. We can use any alphanumeric name, so for simplicity we just use the scaffold IDs also as chromosome IDs. The chromosome map must contain every scaffold ID found in the file, so we just extract the scaffold IDs from the VCF file (not the header, see reason above!) with:

We use cut to “cut” the first column (-f 1) of the file (check cut –help for more info):
bcftools view -H data.noIND04.miss0.5.noMT.thin10.vcf.gz | cut -f 1

(don’t type this yet, see final command below)

As we have multiple variants per scaffold, we also receive every scaffold ID multiple times with the grep command. However, we only want the scaffold ID once, so we pipe the output to uniq, which reduces the list to unique IDs.
| uniq

We need the list in two columns – the original scaffold ID and the new identifier (for which we just use the scaffold ID again). We do this with the language awk to “print” the piped output twice:
| awk ‘{print $0″\t”$0}’

Using only awk ‘{print}’ would print the entire file (once). With the $ variable, we specify a certain part of the file, separated by whitespace. $0 means the entire file, $1 would be the first column, $2 the second etc. As both columns should be separated by a tab, we add “\t” in-between.

Finally, we save the output in a new text file: > data.noIND04.miss0.5.noMT.thin10.scf.chrom-map.txt

Put the parts above together and run the following command on your terminal to get the chromosome map:
bcftools view -H data.noIND04.miss0.5.noMT.thin10.vcf.gz | cut -f 1 | uniq | awk ‘{print $0″\t”$0}’ > data.noIND04.miss0.5.noMT.thin10.scf.chrom-map.txt

Have a look at the output. Does it look like the chromosome map that we need?

Use VCFtools to export our final VCF (data.noIND04.miss0.5.noMT.thin10.vcf.gz) including the chromosome information (data.noIND04.miss0.5.noMT.thin10.scf.chrom-map.txt) to PLINK format. Use data.noIND04.miss0.5.noMT.thin10 as the file prefix for the new file. Check the manual again for the commands you need.

Type ls -l to see the new files.

______________________________________________________________________________________________

Exercise 3 – PLINK

PLINK is a comprehensive genome analysis toolset with an extensive list of functions. It was originally developed for human data (hence has 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. The standard version will allow a maximum of 59 chromosomes, so here we use a special build that supports the use of up to ~64000 scaffolds or contigs.

3.1) Input formats

PLINK accepts VCF files as input, but the preferred format to work with in PLINK are the PLINK text files with endings .ped (and .map), and the smaller binary PLINK files with endings .bed (+ .bim + .fam). Do not confuse this .bed with the UCSC Genome Browser’s BED format, which is totally different. In PLINK, the input data file is split in 2 (or 3) files, each containing solely information on either genotypes and sample information (.ped), only genotypes (.bed), variant positions and information (.map / .bim), or sample information (.fam).

Have a look at the PLINK .ped and .map files you exported with VCFtools!
Tip: Turn off line-wrap when using less with the flag -S

The .ped file contains the variant information with one allele per column (meaning 2 columns are needed for the two alleles of one variant). It has one line per sample with the first six columns being:
1 Family ID
2 Within-family ID (‘IID’; cannot be ‘0’)
3 Within-family ID of father (‘0’ if missing)
4 Within-family ID of mother (‘0’ if missing)
5 Sex code (‘1’ = male, ‘2’ = female, ‘0’ = unknown)
6 Phenotype value (‘1’ = control, ‘2’ = case, ‘-9’/’0’/non-numeric = missing data)

7-n The seventh and eighth columns are the alleles for the first variant, the second variant, etc. Missing data is coded as 0 (or -9).

For our data we do not have information about the father, mother, the sex, or phenotype, so those columns are all set to missing data (0 or -9).

The .map file contains the variant positions. It has 4 columns:
1 Chromosome code. PLINK 1.9 also permits contig or scaffold names here, but most older versions do not.
2 Variant identifier
3 Position in morgans or centimorgans (optional; ‘0’ = unknown)
4 Chromosome position in bp

A detailed description of the file formats is available here.

We went through quite a few steps to obtain our file in PLINK format. Quickly verify that our data is still OK by comparing the genotypes of the first variant (scf7180003948298:263) in the new .ped/.map files with the genotypes in the original VCF file (data.vcf.gz).

By looking at the .ped file, you might have noticed that the “Family IDs” in column 1 (i.e. population/sampling site IDs) are the same as the “Within-family IDs” in column 2 (i.e. individual IDs). However, for analyses of population structure, these should contain information about the population from which the individual was sampled. Let us change that!

We will update the IDs using PLINK. Have a look at the general usage info in the PLINK manual or just type plink in the terminal to get information about how to use PLINK.

Once you’re familiar with the command structure, change the “Family IDs” to “IND” for the samples IND01-05 and to “AFR” for the samples AFR01-AFR05, and save the output in a new file in text format (.ped) with the file prefix data.noIND04.miss0.5.noMT.thin10.newID.
Have a look at the “Update sample information” section of the PLINK manual for how to change the IDs (see Tip below!), and at the “Data management” section for how to save as a new .ped file.

Tip: You have to make a text file with the updated IDs first.

Try again with the –aec flag!

Inspect the new files and the log file.
You may notice the following lines in the log file:

– Ambiguous sex IDs written to data.noIND04.miss0.5.noMT.thin10.newID.nosex
By default, samples with ambiguous sex have their phenotypes set to missing when analysis commands are run, and the IDs of those sampes are written to a file with the ending “.nosex”. You can add the flag –allow-no-sex to prevent this, but since we do not have phenotypes anyway, we do not bother.

– […] Performing single-pass .bed write (863943 variants, 9 people).
–file: data.noIND04.miss0.5.noMT.thin10.newID-temporary.bed +
data.noIND04.miss0.5.noMT.thin10.newID-temporary.bim +
data.noIND04.miss0.5.noMT.thin10.newID-temporary.fam written. […]

PLINK 1.9 automatically converts most non-binary formats (such as .ped, .vcf) to a temporary binary .bed format before starting analysis. With our command above, it re-codes the temporary binary file after changing the family IDs back to the text format and saves it in .ped format.
For large genomic data files it is more convenient to continue working with the smaller binary format, especially if you’re planing on multiple analyses.

Change the command above so that your file is saved in binary .bed/.bim/.fam instead of .ped/.map. Use this section of the manual for help.

Check the new files in your folder (ls -l). You’ll find the binary .bed, and the accompanying position and sample text files .bim and .fam, as well as a log file.
Have a look at the .bim and .fam files to see how they are formatted!

3.2) Filtering

3.2.1) Minor allele frequency (MAF)

We already used VCFtools to filter our file for missing data, the mitochondrial scaffold, and very close sites (within 10bp).
Now we also want to filter on minor allele frequency (MAF), the frequency of the allele that appears less often. Filtering on MAF can be done for multiple reasons. One of them 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, alleles present in only one individual do not add information to 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.
Ignoring the quality, what would you think would be an appropriate threshold to filter on for population structure analysis?

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

If you’re still lost as to which threshold to choose, read hint 2.
If you’ve decided on a threshold, go ahead and apply it. Check this section of the manual on how to filter for MAF. Save the new file in binary .bed format with a new prefix “data.noIND04.miss0.5.noMT.thin10.newID.maf.[your threshold]” and proceed to look up the answer below.

If you’ve read hint 2 and you don’t know if MAF 0.01 would filter anything or not, you may check the manual on how to filter for MAF and just filter for MAF 0.01 to see what happens! Save the output in binary .bed format and name the file test.maf0.01. Have a look at the log file. Have any variants been removed?

An appropriate MAF threshold is given in the answer below:

Go ahead and apply the proposed MAF threshold. See this section of the manual for the commands you need. Save a new .bed file with the prefix data.noIND04.miss0.5.noMT.thin10.newID.maf0.1

3.2.2) Linkage disequilibrium

PLINK has several options to estimate linkage disequilibrium (LD), the non-independent segregation of two loci, and filter variants according to certain distances and thresholds.

For biallelic markers, one of the most commonly used measures for LD is the squared coefficient of correlation (r2) (Hill & Robertson, 1968). It ranges between 0 and 1, where 0 indicates no correlation (= no linkage) and 1 perfect correlation (= perfect linkage disequilibrium).

3.2.1.A) R2 statistics

Estimate the r2 values for the scaffold scf7180003983566. We want to compare only variants that are within 1000 bp of each other, but within this distance, all variants should be compared with all others.

NB! We will only use one scaffold for this exercise since calculating r2 performs and saves pairwise comparisons between loci, which can easily lead to an unmanageable huge data file!! Make sure to specify that you only want to analyze one scaffold with the –chr flag (–chr scf7180003983566).

Check the manual to see the options for analyzing individual scaffolds, and estimating r2.
Save the output with the prefix scf7180003983566.1kb.r0

Tip: Read the 9th and 10th bullet point of the LD statistics reports section to see how to specify a window size and variant distance, and how to request an output listing all r2 values.

Have a look at the output and try to understand the different columns!

Roughly, how many variants are linked according to the r2 values?
Do the r2 values between two sites change with increasing distance?

Plot a histogram of the distribution of r2 values to answer the first question and plot the distance between the compared variants against the r2 values to answer the second question.
Use the plotting options available on your Amazon instance (R, RStudio, RStudio in the Browser (http://YOUR_PUBLIC_DNS:8787)), or download the file to your computer.

In the histogram you can see that most of the SNPs are not, or only weakly correlated.
When we check the scatterplot, we see that LD (in terms of r2) decreases as a function of distance.
The main force influencing the level of LD in certain genomic regions is recombination, and the frequency of recombination between two locations depends on their distance on the chromosome. Therefore, for sites sufficiently distant from each other the frequency of crossovers between them is high enough to remove the correlation between alleles.

This method can also be applied to a whole genome to obtain a genome-wide estimate of linkage decay.

3.2.1.B) R2 pruning

We already filtered our data set for variants that are in close proximity (and thus possible LD) with VCFtools. However, with LD-based pruning we are able to not arbitrarily delete close variants, but restrict our filtering to only those SNPs that are next to each other AND genetically linked based on r2 values.

Can you come up with a command to produce a new data set where all variants within a window of 1 kb and an r2 above 0.8 are removed? Set the step-size (the number of SNPs to shift the window at each step) to 1, so all variants withing 1 kb will be compared with each other.
This is a two-command task. Use the files with prefix data.noIND04.miss0.5.noMT.thin10.newID.maf0.1 as input, and name the output with prefix data.noIND04.miss0.5.noMT.thin10.newID.maf0.1.ldkb1r0.8.

Check the Variant pruning section of the manual for a hint on how to do this.

Bonus task:
If you like, you can check again how many variants have been removed, and if the numbers are consistent with the variants in the files.
You can also produce another r2 report for the filtered file (like we did above) and plot the values again in a histogram and distance plot. In which way are the plots of the LD pruned file different from the plots we produced before?

Finally, can you find out how many variants we filtered out in total compared to the raw data file data.vcf.gz?

Congratulations! You produced a nicely filtered variant file, which can now be exported to different formats and analyzed!

3.3) Export formats

PLINK has various options to export the data into different formats.
Find a list of these formats here.

We already know how to save the data in PLINK .ped format (–recode) and in binary .bed format (–make-bed).
In the same way we can save our file in a format which is suitable to import to R and analyze for example with the R package adegenet (see the next exercise tonight!). Or you can export it in STRUCTURE format, or again in VCF format, or in several other formats.

Export the final data file data.noIND04.miss0.5.noMT.thin10.newID.maf0.1.ldkb1r0.8 in the format suitable for the R package adegenet and you have the possibility to perform, e.g., a principal component analysis (PCA) with it tonight!

Tip: We need the one “without the dominance component”.

Survey

Please let us know how you got along with the exercise in this very short opinion poll.