PCAngsdTutorial

We will go through a some examples on how to use PCAngsd with visualization of the data.

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

Demo 1: Allele Frequencies

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

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

Admixture based on 1st two PC

Let's try to use the PCA to infer admixture proportions based on the first 2 principal components. For the optimization we will use a small penalty on the admixture proportions (alpha).

$PCANGSD -beagle $IN_DIR/Demo1input.gz -o $OUT_DIR/Demo3PCANGSD -admix -admix_alpha 50

Plot the results in R

#open R
pop<-read.table("pop.info",as.is=T)

###########CHECK NAME OF FILE###################
q<-read.table("Demo3PCANGSD.K3.a50.0.qopt")

## order according to population
ord<-order(pop[,1])
barplot(t(q)[,ord],col=2:10,space=0,border=NA,xlab="Individuals",ylab="Admixture proportions")
text(tapply(1:nrow(pop),pop[ord,1],mean),-0.05,unique(pop[ord,1]),xpd=T)
abline(v=cumsum(sapply(unique(pop[ord,1]),function(x){sum(pop[ord,1]==x)})),col=1,lwd=1.2)

## close R

Inbreeding in the admixed individuals

Inbreeding in admixed samples is usually not possible to estimate using standard approaches. Let's try to estimate the inbreeding coefficient of the samples using the average allele frequency.

$PCANGSD -beagle $IN_DIR/Demo1input.gz -o $OUT_DIR/Demo3PCANGSD -inbreed 2 -iter 0

Join names and results, sort the values and look at the results

paste pop.info Demo3PCANGSD.inbreed | LC_ALL=C sort -k3g

The third column is an estimate of the inbreeding coefficient (allowing for negative)

Inbreeding with individual allele frequencies

Now let's try to estimate the inbreeding coefficient of the samples by using the individual allele frequencies predicted by the PCA

$PCANGSD -beagle $IN_DIR/Demo1input.gz -o $OUT_DIR/Demo4PCANGSD -inbreed 2

Join names and results, sort the values and look at the results paste pop.info Demo4PCANGSD.inbreed | LC_ALL=C sort -k3g

Demo 2: Selection

For very resent selection we can look within closely related individuals for example within Europeans.

The objective is to use PCangsd to estimate the covariance matrix while jointly estimating the individuals allele frequencies.

Data file:

• Genotype likelihoods in Beagle format
• ~150k random SNPs with maf > 5%
• Four EU populations with ~100 individuals in each

CEU	European (mostly British)
GBR	Great Britain
IBS	Iberian/Spain
TSI	Italien

Download the beagle genotype likelihood input file

There is a file with the positions and sample information. (cp $ThePath/PCangsd/data/eu1000g.sample.Info)

And a beagle file (EU1000=$ThePath/PCangsd/data/eu1000g.small.beagle.gz)

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

*ANDERS*

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

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

mv Demo2input.gz $IN_DIR

mv Demo2sample.info $IN_DIR

*ANDERS*

wc Demo2sample.info
N=424 #one line for header

#run PCANGSD
$PCANGSD -beagle $EU1000 -o EUsmall -threads 4
#-n $N -threads 20 # from previous year...

Plot the results in R

## R
 cov <- as.matrix(read.table("EUsmall.cov"))

 e<-eigen(cov)
 ID<-read.table("eu1000g.sample.Info",head=T)
 plot(e$vectors[,1:2],col=ID$POP)

 legend("topleft",fill=1:4,levels(ID$POP))

Since the European individuals in 1000G are not simple homogeneous disjoint populations it is hard to use PBS/FST or similar statistics to infer selection based on populating differences. However, PCA offers a good description of the differences between individuals which out having the define disjoint groups.

Now let's try to use the PC to infer selection along the genome based on the PCA

$PCANGSD -beagle $EU1000 -o EUsmall -selection -sites_save #-n $N

plot the results

## function for QQplot
qqchi<-function(x,...){
lambda<-round(median(x)/qchisq(0.5,1),2)
  qqplot(qchisq((1:length(x)-0.5)/(length(x)),1),x,ylab="Observed",xlab="Expected",...);abline(0,1,col=2,lwd=2)
legend("topleft",paste("lambda=",lambda))
}

### read in seleciton statistics (chi2 distributed)
s<-scan("EUsmall.selection.gz")
## make QQ plot to QC the test statistics
qqchi(s)

# convert test statistic to p-value
pval<-1-pchisq(s,1)

## read positions (hg38)
p<-read.table("EUsmall.sites",colC=c("factor","integer"),sep="_")

names(p)<-c("chr","pos")

## make manhatten plot
plot(-log10(pval),col=p$chr,xlab="Chromosomes",main="Manhatten plot")

## zoom into region
 w<-range(which(pval<1e-7)) + c(-100,100)
 keep<-w[1]:w[2]
 plot(p$pos[keep],-log10(pval[keep]),col=p$chr[keep],xlab="HG38 Position chr2")

## see the position of the most significant SNP
 p$pos[which.max(s)]

Now let try to use the PC to infer selection along the genome based on the PCA

python $PCAngsd -beagle $EU1000 -o EUsmall -n $N -selection -sites_save

plot the results

## function for QQplot
qqchi<-function(x,...){
lambda<-round(median(x)/qchisq(0.5,1),2)
  qqplot(qchisq((1:length(x)-0.5)/(length(x)),1),x,ylab="Observed",xlab="Expected",...);abline(0,1,col=2,lwd=2)
legend("topleft",paste("lambda=",lambda))
}

### read in seleciton statistics (chi2 distributed)
s<-scan("EUsmall.selection.gz")
## make QQ plot to QC the test statistics
qqchi(s)

# convert test statistic to p-value
pval<-1-pchisq(s,1)

## read positions (hg38)
p<-read.table("EUsmall.sites",colC=c("factor","integer"),sep="_")

names(p)<-c("chr","pos")

## make manhatten plot
plot(-log10(pval),col=p$chr,xlab="Chromosomes",main="Manhatten plot")

## zoom into region
 w<-range(which(pval<1e-7)) + c(-100,100)
 keep<-w[1]:w[2]
 plot(p$pos[keep],-log10(pval[keep]),col=p$chr[keep],xlab="HG38 Position chr2")

## see the position of the most significant SNP
 p$pos[which.max(s)]