First steps in genomic data analysis

First steps in genomic data analysis

Compiled by Julia M.I. Barth, 21 January 2020

---

Table of contents

---

• Introduction

1. The Variant Call Format (VCF)
2. Hard quality filtering of variants
   • Inspection of quality measurements
   • Application of hard filtering thresholds
3. Further quality filtering
   • Filter on minimum read depth (DP) and genotype quality (GQ)
   • Remove multiallelic SNPs and indels, monomorphic SNPs, and SNPs in the close proximity of indels
   • Remove individuals with a high amount of missing data
   • Remove variants with a high amount of missing genotypes and filter on minor allele frequency
4. Performing a principal component analysis (PCA) to inspect the data structure

---

Introduction

---

Background and Aim
This activity will guide you though some of the the initial steps in population genomic data analyses – that is, the first steps after the sequencing reads were mapped against a reference genome assembly and variants were called (if you’re interested in mapping and variant calling, check out the “Workshop on Genomics” exercises).
The aim is to get some hands-on experience working with VCF files, practice some of the considerations in producing a subset of high quality SNPs, and to explore the data structure using a principal component analysis (PCA). The aim is also to apply and expand on the UNIX and R skills that you learned in the previous activities: UNIX and R.


Learning goals

After this activity, you should be able to perform first quality checks and filtering steps, preparing a population genomic dataset file to explore e.g., population structure.

• Apply the UNIX and R skills you learned in the previous exercises.
• Learn about the VCF format and how to handle and manipulate VCF files.
• Plot different parameters to check the overall quality of the data set.
• Perform hard filtering, removing low quality variants.
• Perform further filtering to prepare the data file for analyses.
• Explore the data structure using a PCA.

Requirements
BCFtools v1.9
PLINK v1.9
EIGENSOFT smartpca v7.2.1
R / Rstudio


The data set

We will work with a data set containing genomic data of 204 individuals of Atlantic cod (Gadus morhua) sampled at seven sites in the Atlantic, North Sea, and Baltic Sea (Barth JMI et al. (2019) Disentangling structural genomic and behavioural barriers in a sea of connectivity. Mol Ecol 28, 1394–1411).

The individuals in this data set were sequenced using Illumina HiSeq paired-end sequencing to ~10x coverage, reads were mapped using BWA-mem to the Atlantic cod reference assembly, and variants were called jointly for a total number of about 860 individuals using the GATK Haplotypecaller pipeline. All filtering steps were performed per chromosome (i.e., per linkage group) for parallelized processing, followed by merging of high quality variants of all linkage groups for follow-up analyses.

In this activity, we will perform similar filtering steps, starting with the unfiltered raw VCF, but we will work with a reduced data file that contains only the 204 individuals in the study mentioned above, and only 10% randomly subsampled variants of one linkage group (LG 5) to reduce time and computer resources (the original compressed LG 5 VCF file contained 3,281,615 variants and had a file size for 11 GB). By randomly subsampling the variants, physically very close variant sites that may share the same information (they are in linkage disequilibrium) have already been excluded, for which reason this tutorial does not contain a step on filtering for linked sites. If you are interested in those steps, have a look at the tutorials from 2016 and 2018.

How to do this activity

• Connect to your Amazon instance via the terminal using SSH.
• Questions or activities are indicated with Q.
• Switch to Rstudio in the browser (http://ec2-XX-XXX-XXX-XXX.compute-1.amazonaws.com:8787) if you see this: R task .
• All text in red on gray background indicates commands or file names that you can type or copy to the terminal. Text in black on gray background refers to terminal outputs.
• Try to work independently using the provided manuals and information, discuss questions with your neighbor. If you get stuck, check the answer-box:
  

  Show me the answer!

  Oh no, only if you get really stuck!! First try to find the answer yourself!
  

There are often alternative options and other solutions to the tasks. E.g., whether you prefer cat over less doesn’t matter in most cases – go ahead and use your solution! Don’t hesitate to ask the workshop team if you have questions!

---

1) The Variant Call Format (VCF)

---

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

Q1 Open the terminal and navigate to the analysis folder including the data file for this activity: /home/popgen/workshop_materials/21_first_steps_genomics.
What is the file size of cod204.lg05.1.vcf.gz? How many bytes and how many MB does it have?

Q2 What is written in the first ten lines of the VCF file cod204.lg05.1.vcf.gz?

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

These first lines, starting with ##, are part of the “VCF header”, which stores 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=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed">: The number of of often the alternative allele is occuring.
• ##FORMAT<ID=GT,Number=1,Type=String,Description="Genotype">: The genotypes coded as allele values separated by / (or | if the data is phased).

The allele values are 0 for the reference allele and 1 for the alternative allele. If there are more alternative alleles the second alternative allele will have the value 2, the third 3, etc. Missing data is denoted by dots ..

• ##contig<ID=LG01,length=28303952> : Part of a list of chromosomes, contigs, or scaffolds of the reference assembly.

For more information about the different tags, have a look at the VCF file format description.

The very last line of the header specifies the column names for the genotype data. The first eight columns are mandatory, follow by “FORMAT” and the sample IDs:
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT AVE1409002 AVE1409003 AVE1409004 [...] TVE130512 TVE130513 TVE130514

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.

Q3 Find a command in BCFtools that allows you to separate the header from the large body of genotype data. Save the header in a separate text file named cod204.lg05.1.header.txt.
To find the command, have a look at the BCFtools manual, or just type bcftools on the command line to see the available options.

Q4 What is the file size of cod204.lg05.1.header.txt, compared to the whole VCF file?

Q5 Open the header on the command line in a text editor. Add a new line below the line ##fileformat=VCFv4.2 stating the date of today: ##fileDate=20200121 .
Save the new header as: cod204.lg05.1.header.ed.txt. Use BCFtools to replace (or “rehead”) the VCF file cod204.lg05.1.vcf.gz with the new edited header and save the re-headed VCF as cod204.lg05.1.ed.vcf.gz.
Make sure your changes took effect and your file is still valid by visually inspecting the file and the file size of cod204.lg05.1.ed.vcf.gz.

Tip: Type bcftools to see general options for the program; then use bcftools reheader to see the more specific options for the reheader function.

In the genotype part of the VCF file, one line corresponds to one variant site. 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(s)). 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 for each variant site is given (e.g., AC is the allele count in genotypes for the “ALT” allele(s), AF the allele frequency of the “ALT” allele, and AN is the total number of alleles, etc. We will use some of the “INFO” field specifications later on for hard quality filtering).

The “FORMAT” column specifies the fields and order of the data that is given with each individual for each genotype and the remaining columns list each individual with the values corresponding to the “FORMAT” field. The first “FORMAT” field is always the genotype (GT), the remaining fields are variable, but in our data file, GT is followed by AD (allelic depths for the ref and alt alleles), DP (read depth), GQ (genotype quality), and PL (phred-scaled genotype likelihood):

This is how the first variant of the genotype part of the VCF file looks like:
LG05 338 . A G 93.01 . AC=1;AF=0.0005862;AN=408;BaseQRankSum=0.674;ClippingRankSum=0;DP=8975;ExcessHet=3.0357;FS=0;InbreedingCoeff=-0.0104;MLEAC=1;MLEAF=0.0005862;MQ=56.57;MQRankSum=0;QD=11.63;ReadPosRankSum=1.24;SOR=0.148;AC_Het=1;AC_Hom=0;AC_Hemi=0 GT:AD:DP:GQ:PL 0/0:20,0:20:57:0,57,855 0/0:9,0:9:27:0,27,345 0/0:10,0:10:30:0,30,381 [...] 0/0:7,0:7:21:0,21,248

Q6 How many variant sites are stored in the VCF?

Q7 Examine the first two variant sites in the VCF file and answer the following questions:

A) What kind of variants do you see (e.g., SNP, insertion/deletion)?

B) What is the genotype (GT) of the first individual (AVE1409002) at the first variant site (LG05 338)?

C) What is the genotype read depth (DP, specified as “FMT/DP” to distinguish it from the mean DP “INFO/DP” per variant site) for the first three individuals (AVE1409002, AVE1409003, AVE1409004) at the first variant site (LG05 338)?

---

2) Hard quality filtering of variants

---

To remove low-quality variants (e.g., variants with a low read depth or variants only supported by poorly aligned reads) from our data set, we will apply hard quality filtering. In the terminology of the authors of GATK, “hard filtering” refers to using specific filtering cutoffs for specific annotations. This stands opposed to using the so-called Variant Quality Score Recalibration (VQSR) method, which does not use any specific “hard cutoffs” but combines information from multiple annotations from a known, highly validated variant training set in a machine learning algorithm. Variants that are not passing the chosen thresholds or the VQSR filtering can either be completely removed or they can be kept in the VCF, but a flag is then added to them noting that they did not pass the filter. Since often such a training set is not available, we will apply “hard quality filtering” here.
The following part of the activity is adjusted from this “GATK tutorial“.

---

2.1) Inspection of quality measurements

Before we choose the measurements and thresholds for hard quality filtering, let’s start by having a look at some of these measurements and the distributions of their values.

We will consider the following seven measurements, which are given for each variant site in the “INFO” column of the VCF (more information about each of these measurements is given further below):
Fisher Strand (FS), Strand Odds Ratio (SOR), RMS Mapping Quality (MQ), Mapping Quality Rank Sum Test (MQRankSum), Quality by depth (QD), Read Position Rank Sum Test (ReadPosRankSum), and the combined depth across samples (INFO/DP).

Below, the measurements of the first variant site are highlighted:
LG05 338 . A G 93.01 . AC=1;AF=0.0005862;AN=408;BaseQRankSum=0.674;ClippingRankSum=0;DP=8975;ExcessHet=3.0357;FS=0;InbreedingCoeff=-0.0104;MLEAC=1;MLEAF=0.0005862;MQ=56.57;MQRankSum=0;QD=11.63;ReadPosRankSum=1.24;SOR=0.148;AC_Het=1;AC_Hom=0;AC_Hemi=0 GT:AD:DP:GQ:PL 0/0:20,0:20:57:0,57,855 0/0:9,0:9:27:0,27,345 0/0:10,0:10:30:0,30,381 [...] 0/0:7,0:7:21:0,21,248

Q8 Use the BCFtools “query” option to extract the Fisher Strand (FS) values of all variants and save the output in a separate text file named cod204.lg05.1_FS.txt.

Tip: Look at the examples given in the manual. The FS values should be saved as one column, with one line for the value of each variant.

Q9 Inspect the file. How many columns are in the file? How many lines? Does the number correspond to the number of variant sites in the VCF? Do the first values correspond to the FS values of the first variants in the VCF file? Apply your UNIX skills!

R task1 Now lets switch to R studio to plot the FS values. Open Rstudio in your internet browser:
http://ec2-XXX-XXX-XXX-XXX.compute-1.amazonaws.com:8787

Open the provided script “gda_plotting.R” in R and follow the points 2) and 2.1) in that R script.

Fisher Strand (FS), as well as Strand Odds Ratio (SOR), are measurements for strand bias – that is a sequencing bias, in which the reads of one DNA strand (e.g., the forward strand) support one genotype (e.g., a heterozygous call), while the reads for the reverse strand support e.g., a homozygous genotype. In cases of extreme strand bias, SNPs are not correctly called.

• Fisher Strand (FS) is the Phred-scaled probability that there is strand bias at that site. It tells us whether the alternate allele was seen more or less often on the forward or reverse strand than the reference allele. A value close to 0 indicates little or no strand bias. GATK recommends hard filtering of variants with FS greater than 60, but this is a very conservative value that will only filter out cases of extreme bias and leave many false positives in the data set.

Q10 Based on the above information and the plotted distribution of FS, choose a FS filtering threshold to remove low quality variants.

Q11 Use the BCFtools “query” option to jointly extract the Fisher Strand (FS), Strand Odds Ratio (SOR), Mapping Quality Rank Sum Test (MQRankSum), Read Position Rank Sum Test (ReadPosRankSum), Quality by depth (QD), RMS Mapping Quality (MQ),and the combined depth across samples (INFO/DP) and save the values separated by the tab character in a separate text file named cod204.lg05.1_FS.SOR.MQRS.RPRS.QD.MQ.DP.txt.

Tip: Each measurement should be saved in one column separated by tabs, with each variant on a single line.

R task2 Now lets switch again to R studio to plot the distributions of all measurements. Use the provided script “gda_plotting.R” and follow the point 2.1.1).

Q12 As you did for the FS values above, use the distribution plots and the information about each quality measurement given below to come up with a meaningful filtering threshold.

• Strand Odds Ratio (SOR) is also a measurement for strand bias, but it takes into account the ratios of reads that cover both alleles and is therefore a better measurement for sites that have more reads in one direction than the other. SOR values range from 0 to greater than 9, a value close to 0 indicates little or no strand bias. GATK recommends hard filtering of variants with SOR greater than 3.
• RMS Mapping Quality (MQ) and Mapping Quality Rank Sum Test (MQRankSum) are both measures for mapping quality. The Root Mean Square (RMS) mapping quality (MQ) of reads across all samples provides an estimation of the overall mapping quality to support a variant call. It includes the standard deviation (SD) of the mapping qualities, where a low SD indicates values close to the mean. Good mapping qualities are around MQ 60. GATK recommends hard filtering of variants with MQ less than 40.
• The Mapping Quality Rank Sum Test (MQRankSum) compares the mapping qualities of reads supporting the reference and the alternate allele. If all reads map correctly, the expected MQRankSum is 0. Negative values indicate a lower MQ of the reads supporting the alternate allele, compared to reads supporting the reference allele. MQRankSum values usually range from about -10.5 to 6.5. GATK recommends hard filtering of extreme outlier variants with MQRankSum less than -12.5.
• Quality by depth (QD) is the “confidence in a call” (i.e., the QUAL score), normalized by the read depth at that site. GATK recommends hard filtering of variants with QD below 2.
• The Read Position Rank Sum Test (ReadPosRankSum) measures the relative position of the reference versus the alternative allele within reads. E.g., SNPs towards the ends of reads may be more likely a sequencing error. A value around 0 means there is little to no difference in where the alleles are found relative to the ends of reads. GATK recommends hard filtering of variants with ReadPosRankSum less than -8.0.
• The combined depth across samples (INFO/DP) is simply the sum of the read depth (i.e., the genotype depth (FORMAT/DP) fields) over all samples. Variants that have a high mean depth across all samples are indicative of either repetitive regions or paralogs. A good rule of thumb is to remove variants that show twice the mean read depth. You can calculate the mean read depth in R: mean(cod204.lg05.1_hf.values$DP). Now, don’t be shocked about the high read depth. This field has not been updated after reducing the vcf from 860 individuals to 204 individuals, so the information is given for all 860 individauls. We can use this field nevertheless for excluding variants with read depth larger than twice the mean.

---

2.2) Application of hard filtering thresholds

We will now apply the filtering thresholds and remove variants that don’t match the chosen thresholds.

Q13 Using the BCFtools “filter” option, come up with a command that will produce a high-quality data file. You may use the thresholds you chose, or adapt them according to the suggestions above.
Don’t forget to specify the input file cod204.lg05.1.vcf.gz, the output format (we’ll use -O z since we would like to save the data file again as compressed VCF), and the output filename (-o cod204.lg05.1.hf.vcf.gz). Type bcftools filter in the command line to see the options.

Tip: It is possible to either use -e and the logical operator “or” || to exclude low-quality variants from the file, or use -i and the logical operator “and” && to include high-quality variants.

Running the filtering will take about five minutes, so time for a break and some coffee!

Q14 Check the filtered file:

• What is the file size?
• How many variant sites have been removed? Tip: see Q6
• Has the FS filter been applied correctly? What is the largest FS value?
  Tip: use bcftools query as you did above (Q8) but instead of saving the output, pipe it directly to sort. Use the correct flags with “sort” (see the options with sort --help and only check the first values by piping the output again to head.

R task3 Extract again all quality measurements (FS, SOR, MQ, MQRS, QD, RPRS, DP) as you did above (see Q11), save the output in a new text file named cod204.lg05.1.hf_FS.SOR.MQRS.RPRS.QD.MQ.DP.txt, and plot the distributions using Rstudio and the script “gda_plotting.R” (section 2.1.2).

---

3) Further quality filtering

---

After hard quality filtering, we have a data set containing only variant sites that we trust with high confidence. However, so far we have only removed variants that had skewed values across all samples – we haven’t yet looked at the individual genotypes at each variant site. For example, we may have kept a variant that has a high quality across all 204 individuals, but there are three individuals where the read depth (FMT/DP) or the genotype quality (GQ) for the individual genotype is very low, so we don’t trust the calls for these individuals.

LG05 338 . A G 93.01 . AC=1;AF=0.0005862;AN=408;BaseQRankSum=0.674;ClippingRankSum=0;DP=8975;ExcessHet=3.0357;FS=0;InbreedingCoeff=-0.0104;MLEAC=1;MLEAF=0.0005862;MQ=56.57;MQRankSum=0;QD=11.63;ReadPosRankSum=1.24;SOR=0.148;AC_Het=1;AC_Hom=0;AC_Hemi=0 GT:AD:DP:GQ:PL 0/0:20,0:20:57:0,57,855 0/0:9,0:9:27:0,27,345 0/0:10,0:10:30:0,30,381 0/0:5,0:5:15:0,15,209 [...] 0/0:7,0:7:21:0,21,248

Furthermore, we will exclude SNPs that are in the close proximity of indels. That is, because the proper alignment of sequences around indels can often be problematic and produce false positives.

For example, dependent on the “gap opening penalty”, a region can be aligned to have four indels (‐) and a SNP.
CTTTGAACCAAATTT----ACTGGAGTTCTGTTGCC
CTTTGAACC---TTTGCTAACT--AGGT-TGTTGCC

Or it can be aligned to have three indels (‐) and four SNPs.
CTTTGAACCAAATTT--TCAGGAGTTCTGTTGCC
CTTTGAACC---TTTGCTAACTAGGT-TGTTGCC

In addition, we will also remove monomorphic SNPs (i.e., SNPs that were called because they differ from the reference, but that are not polymorphic in our sample).

Then, since we are only interested in biallelic SNPs for the downstream analyses in this activity, we will remove indels and multiallelic SNPs (for data sets were indels are also going to be analysed, it is recommended to subset indels form SNPs as the very first step and adjust the above hard filtering thresholds separately).

And last but not least, we will remove individuals with high amounts of missing data, as well as variants with lots of missing genotypes, and we will filter on the minor allele frequency (MAF).

---

3.1) Filter on minimum read depth (DP) and genotype quality (GQ)

We will again use bcftools filter and filter reads that have a read depth below 3 or a genotype quality below 20. The difference is, that instead of the whole variant site, we are now considering single genotypes (the information in the “FORMAT” fields”) using -e 'FMT/DP<3 | FMT/GQ<20' and we will not remove the whole variant site, but only set the respective genotypes to missing (./.) by using the bcftools filter -S . command.
Use the already “hard-filtered” data file as input cod204.lg05.1.hf.vcf.gz, specify a compressed VCF as the output format -O z, and name the new output file (-o cod204.lg05.1.hf.DP3.GQ20.vcf.gz).

Q15 Puzzle the code together with the help of the above information, bcftools filter, and the manual. It takes ~3 min to run the command.

Q16 Let’s have a brief visual check whether our filter has been applied correctly:

With the following command we extract genotypes of samples that have a read depth lower than three -i 'FMT/DP<3' and in the output we would like to see the genotype (GT), followed by the read depth (DP) and the genotype quality (GQ) for that sample -f '[GT=%GT\tDP=%DP\tGQ=%GQ\t]\n'.

Before you execute the following command, have a brief thought about the expected output regarding the genotypes.

bcftools query -i 'FMT/DP<3' -f '[GT=%GT\tDP=%DP\tGQ=%GQ\t]\n' cod204.lg05.1.hf.DP3.GQ20.vcf.gz | less -S

You may also check the genotypes with a read depth of exactly 3. bcftools query -i 'FMT/DP=3' -f '[GT=%GT\tDP=%DP\tGQ=%GQ\t]\n' cod204.lg05.1.hf.DP3.GQ20.vcf.gz | less -S.
You will see that there are only few genotypes with a high genotype quality (GQ) that are not set to missing.

---

3.2) Remove multiallelic SNPs and indels, monomorphic SNPs, and SNPs in the close proximity of indels

• To remove monomorphic SNPs, we will use bcftools filter as before to exclude -e all sites at which no alternative alleles are called for any of the samples AC==0 and all sites at which only alternative alleles are called AC==AN. In addition, to remove SNPs that are within 10 bp of indels, we add the flag --SnpGap 10.
• To reduce the data file to contain only biallelic SNPs, we use the bcftools view command using the flags -m2 -M2 -v snps, which set both the minimum -m and maximum -M number of alleles to 2, and keeps only snps via -v snps.

You can execute both commands in one command line by piping | the output of the first command (bcftools filter …) to the second command (bcftools view …). This way you only have to specify the input file cod204.lg05.1.hf.DP3.GQ20.vcf.gz with the first command, and the output format -O z and output file name -o cod204.lg05.1.hf.DP3.GQ20.allele.vcf.gz with the second command.

Q17 Apply the command by using the information above, bcftools filter and bcftools view, and the manual. It will take ~2 min to run.

Q18 Verify that the filters have been applied by visually inspecting the file (e.g., using zless -S).

Q19 How many variants are left in our data file after all the above filtering steps have been applied?

---

3.3) Remove individuals with a high amount of missing data

Sometimes, sequencing doesn’t work as well for some individuals as it does for others. Such individuals will be marked by low quality genotypes and low read depth, resulting in high amounts of missing data after the above filtering steps.

We can check the amount of missing data by using the bcftools stats command. Have a look at the options by typing bcftools stats in the terminal or check the manual for what it can do. With -s - we can request stats for all samples.

Run bcftools stats -s - cod204.lg05.1.hf.DP3.GQ20.allele.vcf.gz and quickly scroll through the large output.

Find the section PSC, Per-sample counts, which displays some summary statistics per sample. The last column named [14]nMissing displays the number of missing genotypes per sample.
There is also other valuable information that allows us for example to assess the effectiveness of the filtering steps above (e.g., in clolumn [9]nIndels, also under PSC the number should be 0 since we removed all indels).

Also have a quick look at the number of [11]nSingletons, listed in column 11. We will get back to that in part 3.4) of the activity.

Q20 Use your UNIX (and regex) skills to extract only the 3rd and 14th column of the PSC output to get a list of sample names and numbers of missing genotypes. Save the file as: cod204.lg05.1.hf.DP3.GQ20.allele.imiss.
Visually check that the file is as expected.

R task4 Lets switch to R studio to plot the missingness per sample. Use the provided script “gda_plotting.R” and follow the point 3) of the R script.

Q21 Apply the bcftools view -s command to remove the three individuals with more than 70% missing data from the data file cod204.lg05.1.hf.DP3.GQ20.allele.vcf.gz .
Save the output as compressed VCF by using option -O z and specify the output file name -o cod204.lg05.1.hf.DP3.GQ20.allele.missi.vcf.gz with option -o. It takes ~1 min to run this command.

Q22 Count the number of samples in cod204.lg05.1.hf.DP3.GQ20.allele.missi.vcf.gz. Have the three individuals been successfully removed?

---

3.4) Remove variants with a high amount of missing genotypes and filter on minor allele frequency

In the missingness plot above, we saw that the data set in general has a lot of missing data. This is because we set genotypes with low quality and low read depth to missing. As the penultimate step, we will now remove all variants that have more than 20% missing data.

As the last step, we will remove the singletons that you saw above in the bcftools stats command. Singletons – alleles that are observed only once – and other alleles that have a low minor allele frequency (MAF) frequently present false calls, which can be obstructive in downstream analyses by adding noise to the data (note however, that you may also need them e.g., for demographic analyses).
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 and can therefore be removed. Thus, 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.

Q23 Knowing that our data set contains 201 individuals with a minimum of 20 individuals per sampling site, and that we will use it to dissect population structure, what could be an appropriate MAF filtering threshold?

Q24 Use bcftools filter -e to exclude all variants that have more than 20% missing genotypes or have a minor allele frequency smaller or equal to 0.02 from the data file cod204.lg05.1.hf.DP3.GQ20.allele.missi.vcf.gz. Save the output as compressed VCF -O z and use the output file name -o cod204.lg05.1.hf.DP3.GQ20.allele.missi.miss20.maf0.02.vcf.gz.

Q25 Often the alternative allele is also the minor allel, thus we can inspect the alternative allele count (AC) to see whether our MAF filter took effect.
What would be the lowest alternative allele count to expect?

Use again bcftools query and extract the INFO field “AC”, pipe the output into sort and print only the first lines.

OPTIONAL R task: Switch to Rstudio, adjust the code of point 3) in the “gda_plotting.R” script and plot the missingness per sample after the above filtering steps. Has it changed?

Q26 How many variants (in percent) have been removed by all our filtering steps?

---

4) Performing a principal component analysis (PCA) to inspect the data structure

---

The “smartpca” program implemented in the software EIGENSOFT uses a combination of principal component analysis (PCA; a technique used to identify and illustrate the strongest components of variation in a dataset) and new statistics designed specifically for genomic data. The methods are described in Patterson et al.(2006) PLOS Genetics 2(12): e190.

The EIGENSOFT smartpca needs four input files: The “.ped” (containing the genotype information) and “.map” (containing the chromosome/position information) files that we will produce using the software PLINK1.9, a “.pedind” file containing individual information, and a ‘.par’ file specifying the analysis settings.

Q27 Convert the filterd VCF file cod204.lg05.1.hf.DP3.GQ20.allele.missi.miss20.maf0.02.vcf.gz to the PLINK “.ped/.map” format by running the following command:
plink --vcf cod204.lg05.1.hf.DP3.GQ20.allele.missi.miss20.maf0.02.vcf.gz --aec --recode --out cod204.lg05.1.hf.DP3.GQ20.allele.missi.miss20.maf0.02

If you would like to know more about this format have a look here.

Q28 The EIGENSOFT smartpca is often not happy if the chromosome IDs in the “.map” file contain anything else than numbers.
Use sed to replace the “LG” in the cod204.lg05.1.hf.DP3.GQ20.allele.missi.miss20.maf0.02.map file with “” (nothing) and save the output in a new file called tmp.
Verify that the new file is correct. If so, rename tmp again to cod204.lg05.1.hf.DP3.GQ20.allele.missi.miss20.maf0.02.map.

Q29 To generate the “.pedind” file, apply again your UNIX skills and cut the first five 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). In addition, we need a column to assign individuals to populations.

Q30 Add a tab-separated 6th column specifying that AVE1409002-AVE1409024 belong to the population “AVE”, BOR1104001-BOR1205019 belong to the population “BOR”, etc.. .
To do this, “cut” the first three characters (see cut --help) from cod204.lg05.1.hf.DP3.GQ20.allele.missi.miss20.maf0.02.ped and “paste” them to the tmp file. Save the output as cod204.lg05.1.hf.DP3.GQ20.allele.missi.miss20.maf0.02.pedind.
Remove the “tmp” file.

Q31 A template parameter file (template.par) and a smartpca manual page (eigensoft_smartpca_readme.txt) can be found in the analysis folder. Using the text editor of your choice, modify the template.par to include the respective input file names, save the par file as cod204.lg05.1.hf.DP3.GQ20.allele.missi.miss20.maf0.02.par, and run the smartpca with the syntax smartpca -p cod204.lg05.1.hf.DP3.GQ20.allele.missi.miss20.maf0.02.par > pca01.out.

R task5 Lets switch a last time to R studio to plot the PCA. Use the provided script “gda_plotting.R” and follow point 4) in the R script.

Well done!! Now it is time for socializing. My suggestion for todays pub is the “Zapa bar”.