This section will focus on analyzing SNP data from North American gray wolves (Canis lupus). The original data consist of a few hundred wolf samples that were enriched for a custom capture array consisting of several thousands of regions across the wolf genome.

Our data consist of genic SNPs genotyped for three wolf ecotypes: (1) the British Columbia coastal, (2) the High Arctic, and (3) the Boreal Forest. These wolves have pretty striking ecological and behavioral differences that are in part driven by differences amongst their environments. We used these data to look at population genetic variation and signals of local adaptation amongst these ecotypes (Schweizer et al. 2015, Molecular Ecology).

Figure 5.1: Ranges and ecotypes of sampled gray wolf populations.

This data are also included in the project repo you previously cloned. Let's navigate to the wolf-snps folder now to minimize the paths we have to type when using the data. You should already be in the root project directory (~/congen-bioinformatics/), so you should just have to type the command:

cd data/wolf-snps/

If you are somewhere else other than the root project directory, first navigate back to that folder and then change to the wolf-snps folder. Or, if you're familiar enough with the file tree, navigate straight to the data folder!

I've downsampled the full data set to a subset of 35 individuals and a subset of genic SNPs. We will be working with these SNP data in a VCF file. VCF stands for variant call format and is a tab-delimited plain text file that stores genetic variants, such as SNPs, insertions, and deletions, found by comparing sequenced samples to a reference genome. Each site is described by its location, alleles, and quality scores, plus each individual's genotype, depth of coverage, and more, at that site. Note that sites can be both variable and non-variable.

The full VCF file format specifications can be found here: samtools VCF specs

VCF files also contain a header, denoted with the # character at the beginning of the line. The header can be a single line or thousands of lines and provides information on how the VCF file was generated, to which reference genome the sequence data was aligned, and definitions of what the various fields in each column of the data mean. Given all this information, VCF files can potentially be enormous in size and cumbersome to work with.

Luckily, there is software out there to help! PLINK is a widely used toolset for working with genotype data. It was originally developed for human genome-wide association studies, but it works with data from any diploid organism and can read VCF files directly. PLINK can be used to filter SNPs and individuals, assess data quality such as missingness, and calculate a variety of statistics like allele frequencies, Hardy-Weinberg equilibrium, and per-individual genotype counts. We'll be using PLINK 2, the current version of the program, which is run with the command plink2.

Unlike many of the Unix commands we've seen so far, PLINK doesn't print its results to the terminal to be piped (|) or redirected (>). Instead, it writes its results to files whose names start with a prefix you choose with the --out option, and the file extension is determined by the analysis. Every run also writes a .log file that records the command you ran and a summary of what PLINK did, which is handy for checking your work later. In general, the program is called as follows:

plink2 <input option> <input file> <analysis options> --out <output prefix>

Here are some more explicit examples of PLINK commands:

# Calculate allele frequencies for every SNP (writes input_freqs.afreq)
plink2 --vcf input.vcf --dog --freq --out input_freqs

# Report the amount of missing data per sample and per SNP (writes input_missing.smiss and input_missing.vmiss)
plink2 --vcf input.vcf --dog --missing --out input_missing

# Keep only the samples listed in a file and write them to a new VCF (writes input_subset.vcf)
plink2 --vcf input.vcf --dog --keep samples.txt --export vcf --out input_subset

grep -v "#" Filtered_NAwolf_n35_variableSites_GenicRegions.recode.vcf | cut -f1 | sort | uniq -c

For this section of the workshop, we're going to start with a simple task, calculating heterozygosity, for a single wolf ecotype. Then we're going to learn how to automate it over multiple populations through using loops and writing a script.

Before we do that, let's get some basic info about our called variants in the VCF file.

TASK 1: How many SNPs are there?

We know that our VCF file has a number of header lines (which don't count towards the actual number of SNPs) and then a separate line for each SNP. Because of this, if we just use wc -l we will count all lines of the file, which doesn't tell us the number of SNPs.

But by using grep with the -v flag, we can find all the lines of text that don't contain a character, string, or pattern. Then we can use our handy | to get our answer:

grep -v "^#" Filtered_NAwolf_n35_variableSites_GenicRegions.recode.vcf | wc -l

This should output the number of lines in the VCF file that don't contain the "#" character:

19620

Since the only lines with "#" should be the informational header lines and all other lines represent one SNP, this means there are 19,620 SNPs represented in this VCF file.

TASK 2: Calculate heterozygosity for the coastal wolf samples.

We're going to use one of the functions of PLINK to count the homozygous and heterozygous genotypes (--sample-counts) of each individual within only the coastal wolf population. PLINK lets us specify which individuals to include in an analysis through the --keep option with a file of sample IDs.

plink2 --vcf Filtered_NAwolf_n35_variableSites_GenicRegions.recode.vcf --dog --keep pop_coastal.txt --sample-counts cols=hom,het --out Filtered_NAwolf_n35_variableSites_GenicRegions_coastal

cat Filtered_NAwolf_n35_variableSites_GenicRegions_coastal.log

wc -l pop_coastal.txt

sed '1s/.$/X/' pop_coastal.txt > pop_coastal_typo.txt

plink2 --vcf Filtered_NAwolf_n35_variableSites_GenicRegions.recode.vcf --dog --keep pop_coastal_typo.txt --sample-counts cols=hom,het --out typo_test

If we use ls to look at the output files, we can see that PLINK has provided a .log file, which contains our command and a summary of the run, as well as a .scount file that provides the per-individual genotype counts.

But does this file actually give us the heterozygosity?

head Filtered_NAwolf_n35_variableSites_GenicRegions_coastal.scount

No, the output provides the number of homozygous genotypes and the number of heterozygous genotypes for each individual, below a header line that starts with #IID. But we can use this information with a handy one-liner in awk to calculate heterozygosity, the proportion of each individual's genotyped sites that are heterozygous:

awk 'NR > 1 {print $1, $2, $3, $3/($2+$3)}' Filtered_NAwolf_n35_variableSites_GenicRegions_coastal.scount > Filtered_NAwolf_n35_variableSites_GenicRegions_coastal.scount.calc

One thing to be clear about before we go further: the denominator here is the number of SNPs in this file that the individual was genotyped at, not the number of bases in the wolf genome. Our VCF only contains sites that were variable somewhere in the full data set, so what we're calculating is the proportion of those SNPs at which an individual is heterozygous. That's a good number for comparing individuals and populations to each other here, because every sample is measured against the same set of sites. But it isn't genome-wide heterozygosity, which is reported per base pair and is a far smaller number. Calculating that would also require the invariant sites, plus a record of which positions had enough coverage to be called at all.

wc -l pop_coastal.txt Filtered_NAwolf_n35_variableSites_GenicRegions_coastal.scount Filtered_NAwolf_n35_variableSites_GenicRegions_coastal.scount.calc

sort pop_coastal.txt > pop_coastal_ids.txt

awk '{print $1}' Filtered_NAwolf_n35_variableSites_GenicRegions_coastal.scount.calc | sort > calc_coastal_ids.txt

diff pop_coastal_ids.txt calc_coastal_ids.txt

grep "^#CHROM" Filtered_NAwolf_n35_variableSites_GenicRegions.recode.vcf | tr '\t' '\n' | grep -n -x "SAMPLE_ID"

grep -v "^#" Filtered_NAwolf_n35_variableSites_GenicRegions.recode.vcf | cut -f COLUMN | cut -d: -f1 | sort | uniq -c

grep "SAMPLE_ID" Filtered_NAwolf_n35_variableSites_GenicRegions_coastal.scount.calc

Using this, we can also calculate the mean population heterozygosity using awk:

awk '{sum += $4} END {if (NR > 0) print sum / NR }' Filtered_NAwolf_n35_variableSites_GenicRegions_coastal.scount.calc

sort -k4,4n Filtered_NAwolf_n35_variableSites_GenicRegions_coastal.scount.calc

awk 'NR > 1 {print $1, $2, $3, $2/($2+$3)}' Filtered_NAwolf_n35_variableSites_GenicRegions_coastal.scount

Great! So now we have calculated the per-individual heterozygosity and the population heterozygosity for the coastal populations. Let's figure out how to repeat this set of commands for any number of additional populations.

Instead of hard-coding our command to be specific to the coastal wolves, we can use shell variables. We indicate a shell variable like so:

wolf_pop=coastal

We can check that variable assignment was successful by printing the value with the echo command:

echo ${wolf_pop}

We can use variables in any command in the terminal. For instance, we can run the same command as before to count genotypes, but instead of typing "coastal" in the command, we can use our variable:

plink2 --vcf Filtered_NAwolf_n35_variableSites_GenicRegions.recode.vcf --dog --keep pop_${wolf_pop}.txt --sample-counts cols=hom,het --out Filtered_NAwolf_n35_variableSites_GenicRegions_${wolf_pop}

This provides the same output as the command we ran earlier. So, we could assign a new value to the wolf_pop variable. Notice also how I have been consistent in how I named all of my pop files so that I can substitute an ecotype name but otherwise still have that same general file name structure.

wolf_pop=highArctic

We can check that variable assignment was successful by printing the value with the echo command:

echo ${wolf_pop}

Now that we've reassigned the wolf_pop variable to the highArctic population, we can run the same series of commands as above, but PLINK will use the samples and provide output for this population.

First, run PLINK to count the homozygous and heterozygous genotypes:

plink2 --vcf Filtered_NAwolf_n35_variableSites_GenicRegions.recode.vcf --dog --keep pop_${wolf_pop}.txt --sample-counts cols=hom,het --out Filtered_NAwolf_n35_variableSites_GenicRegions_${wolf_pop}

Then run awk to calculate heterozygosity for each individual in the highArctic sample. Again note that we use the wolf_pop variable in our input and output file names:

awk 'NR > 1 {print $1, $2, $3, $3/($2+$3)}' Filtered_NAwolf_n35_variableSites_GenicRegions_${wolf_pop}.scount > Filtered_NAwolf_n35_variableSites_GenicRegions_${wolf_pop}.scount.calc

And then awk again with the wolf_pop variable to get mean population heterozygosity:

awk '{sum += $4} END {if (NR > 0) print sum / NR }' Filtered_NAwolf_n35_variableSites_GenicRegions_${wolf_pop}.scount.calc

Now, we're going to add some of these commands with variables to a script so that we can automate the series of commands to generate per-individual and per-population heterozygosity.

A script is simply a text file where we can write and save the code and commands we want to run. Scripts give us the flexibility to perform complex operations based on conditions and ranges of values. Since a script is a text file, we're going to create a new text file and copy some commands into it. We're going to use an in-terminal text editor, vim, to do this, although there are several other options.

vim het_analysis.sh

You should now be in the vim text editor. Your screen should be cleared of commands and the cursor should be blinking on top, with tilde characters (~) on every other line. Importantly, vim has both a viewing mode and an editing mode. When you create a new file with vim, as we have here, you start out in viewing mode. We can switch to editing mode by pressing the a key. Once you do this, you'll see the text -- INSERT -- at the bottom of the screen, indicating we have switched to editing mode.

Now, we can copy and paste our commands for calculating heterozygosity into this file. Of course, there are other ways we could get the commands into the file (such as using echo and redirecting), but we'll stay simple today.

Copy the following commands:

wolf_pop=$1
plink2 --vcf Filtered_NAwolf_n35_variableSites_GenicRegions.recode.vcf --dog --keep pop_${wolf_pop}.txt --sample-counts cols=hom,het --out Filtered_NAwolf_n35_variableSites_GenicRegions_${wolf_pop}
awk 'NR > 1 {print $1, $2, $3, $3/($2+$3)}' Filtered_NAwolf_n35_variableSites_GenicRegions_${wolf_pop}.scount > Filtered_NAwolf_n35_variableSites_GenicRegions_${wolf_pop}.scount.calc
awk '{sum += $4} END {if (NR > 0) print sum / NR }' Filtered_NAwolf_n35_variableSites_GenicRegions_${wolf_pop}.scount.calc

The last three commands should be familiar to you already. They are the same commands we used above for the coastal population, but we've substituted the string "coastal" for the variable wolf_pop.

The first line, wolf_pop=$1, is our variable assignment. It tells the script to assign the variable wolf_pop the value of $1. Well what is $1? It has a dollar sign ($) preceding it, so we know it is a variable. But we haven't assigned anything to it yet.

Like other commands, your own scripts can take input options from the command line. By default, a shell script (like the one we are writing) will assign these options to internal variables: $1 for the first option, $2 for the second one, and so on. So, wolf_pop will be assigned whatever we type into the terminal when we run the script.

Let's exit out of vim to try and run our script. First, switch back to viewing mode by pressing the ESC key. You should no longer see the text -- INSERT -- at the bottom of the screen. When in viewing mode, any input we type is interpreted as internal vim commands, such as saving, exiting, etc. Let's tell vim to save the changes we've made to the file with the write command (w) and then quit (q) the text editor by typing:

:wq

TASK 3: Make the script executable and try it

We will use the chmod command to change the permissions on the script so the shell knows that its contents should be interpreted as commands. This is known as making the file executable:

chmod a+x het_analysis.sh

Now let's try to execute our file on the coastal population:

./het_analysis.sh coastal

If we use our handy ls -lht command, we can see that the three newest files are of the coastal wolf samples, and that the script we wrote printed the mean population heterozygosity to the screen.

plink2 --version

Documenting our code with comments and print statements

It is also good practice to document what each step of our code does in human-readable language. This is called commenting our code. The way to comment code differs between programming languages, but a fairly common method (used in Python, R, and the shell scripting we're doing) is with the # symbol. For languages that use this symbol, any line starting with # will be interpreted as a comment rather than a command to be run.

We can also provide some helpful context to the output of our script. Rather than just printing out a number that may be confusing to someone else using the script (or even ourselves if we re-run this same script months or years later), we can use echo within the script to print some more information, like that the number printed represents mean population heterozygosity.

Feel free to open up your script and follow along as I comment my version and add contextual print statements:

vim het_analysis.sh

TASK 4: Which population has the highest mean heterozygosity?

To find this out now, all we have to do is run our script on all three populations:

./het_analysis.sh coastal
./het_analysis.sh highArctic
./het_analysis.sh borealForest

This shows that the boreal forest population has the highest mean heterozygosity at 0.190026. Remember what that number is: the average proportion of our panel of SNPs at which an individual is heterozygous. It's the right number for ranking these three ecotypes against each other, but not comparable to heterozygosity values you'd read in a paper.

If only there was a way to automate the analysis even further... Aha! We can use for loops to call the script multiple times with different input values, in this case our 3 population strings, for the wolf_pop variable. for loops in a shell script are generated in a similar way as in other programming languages:

for <variable> in <iterable>; do <code>; done

For instance, we could do the following to loop over letters in the alphabet from A to F and print them out:

for letter in {A..F}; do echo $letter; done

Or the following to subtract 1 from each integer from 1 to 4 and print the resulting value to the screen:

for num in 1 2 3 4; do echo "The value of $num - 1 is " `expr $num - 1`; done

TASK 5: Use a for loop to calculate heterozygosity for the three wolf ecotypes

for pop in coastal highArctic borealForest; do ./het_analysis.sh $pop; done

This should produce the same output as when we ran the three commands manually above!

Hopefully you've become more familiar with how to use variables in the command-line and in scripts, how to write a script and edit it using vim, and how to use for loops to make your life easier!