EvalAdmix: Difference between revisions

Revision as of 13:24, 21 July 2019

evalAdmix allows to evaluate the model fit of an admixture analysis. It takes as input the genotype data (either called genotypes in plink files or genotype likelihoods beagle files) used in the admixture analysis and the estimated ancestral frequency and admixture propotions files (P and Q files).

The output is a pairwise correlation of residuals matrix between individuals. The correlation will be 0 in case of a good fit of the data to the admixture model. When something is wrong, individuals from the same population will be positively correlated; and individuals from different populationts but that share one or more ancestral populations as admixture sources will have a negative correlation. Positive correlation between a pair of individuals might also be due to relatedness.

Error creating thumbnail: File missing

Download and Installation

evalAdmix can be installed from github

git clone https://github.com/GenisGE/evalAdmix.git
cd evalAdmix
make

Quick start

./evalAdmix -beagle inputBeagleFile.gz -fname myoutfiles.fopt.gz -qname myoutfiles.qopt -P 10

./evalAdmix -plink inputPlinkPrefix -fname myoutfiles.fopt.gz -qname myoutfiles.qopt -P 10

-beagle beagle file of genotype likelihoods
-plink binary plink file prefix with genotype data
-fname file with ancestral frequencies
-qname file with admixture proportions
-o prefix of output file names
-P Number of threads used

Parameters

./evalAdmix

Arguments:
	Required:
		-plink path to binary plink file (excluding the .bed)
		or
	      	-beagle path to beagle file containing genotype likelihoods (alternative to -plink)
		
		-fname path to ancestral population frequencies file
	       	-qname path to admixture proportions file
		
	Optional:       
	
	       -o name of the output file
	       
	 Setup (optional):
	 
	       -P 1 number of threads
	       -autosomeMax 23	 autosome ends with this chromsome
	       -nIts 5	 number of iterations to do for frequency correction; if set to 0 calculates correlation without correction (fast but biased)
	       -useSites 1.0	 proportion of sites to use to calculate correlation of residuals
	       -useInds filename     path to tab delimited file with first column containing all individuals ID and second column with 1 to include individual in analysis and 0 otherwise (default all individuals are included)
	       -misTol 0.05 	 tolerance for considering site as missing when using genotype likelihoods. Use same value as used in NGSadmix to keep compatibility when using genotype likelihoods (-beagle)
	       -minMaf 0.05 	 minimum minor allele frequency to keep site. Use same value as used in NGSadmix to keep compatibility when using genotype likelihoods (-beagle)

Input Files

Plink

Genotype data files in binary PLINK format (.bed .fam .bim).

Beagle genotype likelhoods

The input file contains genotype likelihoods in a .beagle file format [1]. and can be compressed with gzip.

BAM files

If you have BAM files you can use ANGSD to produce genotype likelihoods in .beagle format. Please see Creation of Beagle files with ANGSD

VCF files

If you already have made a VCF file that contains genotype likehood information then it should be possible to convert .vcf files with genotype likelihoods to .beagle file via vcftools [2]

vcftools --vcf input.vcf --out test --BEAGLE-GL --chr 1,2

Chromosome has to be specified.

You can also use bcftools' [3] 'query' option for generating a .beagle file from a .vcf file.

Output Files

The analysis performed by evalAdmix produces 1 file, containing a tab delimited N times N symmetric correlation matrix, where column i in line j contains the correlation of residuals between individual i and j, and the diagonal values (self-correlation) are set to NA:

NA 0.008609 -0.006919 0.002731 0.020224
0.008609 NA 0.000033 0.004968 -0.008470
-0.006919 0.000033 NA 0.006982 0.005664
0.002731 0.004968 0.006982 NA 0.000521
0.020224 -0.008470 0.005664 0.000521 NA

Run command example

Genotype data

Download the input file containing genotype in binary plink format

wget http://pontus.popgen.dk/albrecht/open/admixTjeck/plink/admixTjeck2.bed
wget http://pontus.popgen.dk/albrecht/open/admixTjeck/plink/admixTjeck2.bim
wget http://pontus.popgen.dk/albrecht/open/admixTjeck/plink/admixTjeck2.fam

Run ADMIXTURE [4] to obtain admixture proprotions

admixture admixTjeck2.bed 3

Run evalAdmix

./evalAdmix -plink admixTjeck2 -fname admixTjeck2.3.P -qname admixTjeck2.3.Q -P 20

Plot results in R

# read population labels and estimated admixture proportions
pop<-read.table("admixTjeck2.fam")
q<-read.table("admixTjeck2.3.Q")


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


r<-as.matrix(read.table("output.corres.txt"))

# Plot correlation of residuals
source("NicePlotCorRes.R")
plotCorRes(cor_mat = r, pop = as.vector(pop[,2]), title="Evaluation of 1000G admixture proportions with K=3", max_z=0.1, min_z=-0.1)

Low depth sequencing data

Download the input file containing genotype likelihoods in beagle format

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

Execute NGSadmix to obtain admixture proportions

./NGSadmix -likes Demo2input.gz -K 3 -P 20 -o myoutfiles -minMaf 0.05

Input file = Demo2input.gz

Ancestral Populations K=3

Computer cores = 20 (-P 20).

Output prefix = myoutfiles (-o myoutfiles)

SNPs with MAF > 5% (-minMaf 0.05)