PCAngsdTutorial

We will go through a simple and more complex example on how to use PCAngsd with visualization of the data.

PCA with admixture aware priors

This example will perform a PCA analysis on 1000 genotype likelihoods.

Every time you open a new terminal window, set the paths to the program and the input file.

Set the path to PCAngsd

PCANGSD=~/Software/pcangsd

Test the link

ls $PCAngsd

Create directories

Create the directories that will be used for working:

mkdir Demo

cd Demo

mkdir Data

mkdir Results

Set the paths to your local directories

IN_DIR=Demo/Data

OUT_DIR=Demo/Results

Test the links

ls $IN_DIR

ls $OUT_DIR

Download the beagle genotype likelihood input file

PCAngsd uses Genotype Likelihoods (GLs) in .beagle format as input. The input file has been created for you.

Download the files and move them to your input folder (for example, $IN_DIR):

*ANDERS*

wget popgen.dk/software/download/NGSadmix/data/Demo1input.gz

wget popgen.dk/software/download/NGSadmix/data/Demo1pop.info

mv Demo1input.gz $IN_DIR

mv Demo1pop.info $IN_DIR

*ANDERS*

View the genotype likelihood beagle file

• In general, the first three columns of a beagle file contain marker name and the two alleles, allele1 and allele2, present in the locus (in beagle A=0, C=1, G=2, T=3). All following columns contain genotype likelihoods (three columns for each individual: first GL for homozygote for allele1, then GL for heterozygote and then GL for homozygote for allele2). Note that the GL values sum to one per site for each individual. This is just a normalization of the genotype likelihoods in order to avoid underflow problems in the beagle software, but it does not mean that they are genotype probabilities.

• In order to see the first 10 columns and 10 lines of the input file, type:

gunzip -c $IN_DIR/Demo1input.gz | head -n 10 | cut -f 1-10 | column -t

• Use this command to count the number of lines of the input file. The number of lines, indicates the number of loci for which there are GLs plus one (as the command includes the count of the header line):

gunzip -c $IN_DIR/Demo1input.gz | wc -l

View population information file

To view a summary of the population information file, cut the first column, sort and count:

cut -f 1 -d " " $IN_DIR/Demo1pop.info | sort | uniq -c

Lets make a population label file and place it in the output directory

cut -f1 -d" " $IN_DIR/Demo1pop.info > $OUT_DIR/poplabel

Run PCAngsd

The program estimates the covariance matrix that can then be used for PCA.

Estimating Individual Allele Frequencies

$PCANGSD -beagle $IN_DIR/Demo1input.gz -o $OUT_DIR/Demo1PCANGSD

Plot the results in R

#open R
pop<-read.table("pop.info")
C <- as.matrix(read.table("Demo1PCANGSD.cov"))
e <- eigen(C)
pdf("PCAngsd1.pdf")
plot(e$vectors[,1:2],col=pop[,1],xlab="PC1",ylab="PC2", main="individual allele frequency")
legend("top",fill=1:5,levels(pop[,1]))
dev.off()
## close R

To view the plot: evince PCAngsd1.pdf

Without Estimating Individual Allele Frequencies

Try the same analysis but without estimating individual allele frequencies. This is the same as using the first iteration of the algorithm.

$PCANGSD -beagle $IN_DIR/Demo1input.gz -o $OUT_DIR/Demo2PCANGSD -iter 0

Plot the results in R

#open R
pop<-read.table("pop.info")
C <- as.matrix(read.table("Demo2PCANGSD.cov"))
e <- eigen(C)
pdf("PCAngsd2.pdf")
plot(e$vectors[,1:2],col=pop[,1],xlab="PC1",ylab="PC2",main="joint allele frequency")
legend("top",fill=1:5,levels(pop[,1]))
dev.off()
## close R

To view the plot: evince PCAngsd2.pdf