FastNGSadmix
This page contains information about the program fastNGSadmix, a very fast tool for finding admixture proportions from NGS data of a single individual and a method for doing PCA of NGS data, using the estimated admixture proportions. It is based on genotype likelihoods. It also read plink files. The admixture estimation part is written in C++ and the PCA part is written in R.
Download
The program can be downloaded from github:
https://github.com/e-jorsboe/fastNGSadmix
git clone https://github.com/e-jorsboe/fastNGSadmix.git; cd fastNGSadmix make
So far it has only been tested on Linux systems.
A data folder that has an already made reference panel (more on this in the Making a reference panel section) for the admixture estimation and the genotypes of the reference individuals in plink format for the PCA analysis, as well as a .sites file for when generating genotype likelihoods with ANGSD. And an example folder with example data for trying out the method, can be downloaded from:
wget popgen.dk/software/download/fastNGSadmix/data.tar.gz wget popgen.dk/software/download/fastNGSadmix/example.tar.gz
Use curl if you are on a MAC
They can be unpacked thus:
tar -xzf data.tar.gz tar -xzf example.tar.gz
Quick start
#beagle file with genotype likelihoods GL=example/yriFrenchHan_depth05.beagle.gz #plink file with genotypes for a single individual PLINKFILE=example/NA20502_TSI #frequencies from reference panel REF=data/refPanel_humanOrigins_7worldPops.txt #number of individuals in each population of reference panel NIND=data/nInd_humanOrigins_7worldPops.txt #Estimate admixture proportions from the $GL genotype likelihoods file ./fastNGSadmix -likes $GL -fname $REF -Nname $NIND -out yriFrenchHan_depth05 -whichPops French,Han,Yoruba #perform PCA based on the calculated admixture proportions from the previous command and using the $GL genotype likelihoods file Rscript R/fastNGSadmixPCA.R -likes $GL -qopt yriFrenchHan_depth05.qopt -out yriFrenchHan_depth05 -ref data/humanOrigins_7worldPops #Estimate admixture proportions from the $PLINKFILE plink file ./fastNGSadmix -plink $PLINKFILE -fname $REF -Nname $NIND -out NA20502_TSI -whichPops French,Han,Yoruba #perform PCA based on the calculated admixture proportions from the previous command and using the $PLINKFILE plink file Rscript R/fastNGSadmixPCA.R -plinkFile $PLINKFILE -qopt NA20502_TSI.qopt -out NA20502_TSI -ref data/humanOrigins_7worldPops
Quick start - 1000 genomes reference panel
A reference panel was created using the 1000 genomes data, taking 20 individuls from the following populations: Utah Residents (CEPH) with Northern and Western European Ancestry (CEU), Han Chinese (CHB), Finnish (FIN), Peruvian (PEL), Punjabi (PJL) and Yoruba (YRI), in order to cover most of the worldwide genetic diversity. For the Peruvians I made sure to only select individauls with no European admixture. The sites were filtered so MAF > 0.05 and missingness < 0.01 and only SNPs were included, leaving 6 676 750 sites.
First download these files:
wget popgen.dk/software/download/fastNGSadmix/data1000genomes.tar.gz wget popgen.dk/software/download/fastNGSadmix/example.tar.gz
Use curl if you are on a MAC
They can be unpacked thus:
tar -xzf data1000genomes.tar.gz tar -xzf example.tar.gz
And then you can easily run your analysis on a .bam file - requires ANGSD:
#example bam file with genotype data BAM=example/smallNA12874.mapped.ILLUMINA.bwa.CEU.low_coverage.20130415.bam #frequencies from reference panel REF=data1000genomes/refPanel_1000genomesRefPanel.txt #sites file to be used with ANGSD SITES=data1000genomes/1000genomesRefPanel.sites #number of individuals in each population of reference panel NIND=data1000genomes/nInd_1000genomesRefPanel.txt #genotypes of the 1000 genomes data - for the PCA part GENO=data1000genomes/1000genomesRefPanel #use ANGSD for generating a beagle file with the genotypes likelihoods ./angsd -i $BAM -GL 2 -sites $SITES -doGlf 2 -doMajorMinor 3 -minMapQ 30 -minQ 20 -doDepth 1 -doCounts 1 -out CEU.low_coverage #Estimate admixture proportions from the $GL genotype likelihoods file ./fastNGSadmix -likes CEU.low_coverage.beagle.gz -fname $REF -Nname $NIND -out CEU.low_coverage -whichPops CEU,CHB,YRI,PEL #do PCA based on the calculated admixture proportions from the previous command, using the $GL genotype likelihoods file Rscript R/fastNGSadmixPCA.R -likes CEU.low_coverage.beagle.gz -qopt CEU.low_coverage.qopt -out CEU.low_coverage -ref $GENO
Please be aware if you are not a server with a lot of RAM, that the PCA part might consume a lot of RAM, especially if you are including all ref sites for the PCA (changing the settings to -overlapRef 0).
Input Files
Input files are genotype likelihoods in the genotype likelihood beagle input file format [1]. Or called genotypes in the binary plink files (*.bed) format [2] We recommend ANGSD for easy transformation of Next-generation sequencing data to the beagle format and plink2 for handling plink files. See the ANGSD wiki for installing and other documentation
The example below shows how to make a beagle file of genotype likelihoods using ANGSD.
#BAM/CRAM file BAM=example/smallNA12874.mapped.ILLUMINA.bwa.CEU.low_coverage.20130415.bam #sites in the referene panel including major and minor allele SITES=data/humanOrigins_7worldPops.sites #run ANGSD ./angsd -i $BAM -GL 2 -sites $SITES -doGlf 2 -doMajorMinor 3 -minMapQ 30 -minQ 20 -doDepth 1 -doCounts 1 -out outName
Example of a beagle genotype likelihood input file for 1 individual.
marker allele1 allele2 Ind0 Ind0 Ind0 1_14000023 1 0 0.941 0.058 0.000 1_14000072 2 3 0.709 0.177 0.112 1_14000113 0 2 0.855 0.106 0.037 1_14000202 2 0 0.835 0.104 0.060 ...
Where the marker column is recommended to be chr_pos, instead of rs ID - but custom ids can be used.
If one wants to analyse a haploid organism, one should create a regular beagle file (with 6 columns), but only fill out the two first genotype likelihood columns, for respectively specifying the genotype likelihood of allele1 and allele2.
Example of haploid data in a beagle file:
marker allele1 allele2 Ind0 Ind0 Ind0 1_752566 A G 0.50 0.50 0 1_842013 T G 0.99 0.01 0 1_891021 A G 0.99 0.01 0 ...
The reference panel looks like this:
id chr pos name A0_freq A1 French Han Chukchi Karitiana Papuan Sindhi Yoruba 1_752566 1 752566 rs3094315 G A 0.166 0.061 0.37 0.0833 0.071 0.306 0.671 1_842013 1 842013 rs7419119 G T 0.22 0.106 0.565 0.083 0.036 0.361 0.143 1_891021 1 891021 rs13302957 G A 0.060 0.197 0.326 0.75 0.893 0.056 0.264 1_903426 1 903426 rs6696609 C T 0.62 0.636 0.705 0.25 0.179 0.583 0.257 ...
Where again we have the id column (like marker column in beagle file) that is recommended to be chr_pos which is being used for detecting the overlap with the input - custom ids can be used. The frequencies have to be in direction of the A0_freq allele. Custom ids can be used for the id column in the .beagle and refPanel files, as the program matches sites on id_A0_A1 (alleles are sorted).
It should also be noted that the program quits if there are duplicate sites in the input or the reference panel.
The number of individuals in each reference population looks like this:
French Han Chukchi Karitiana Papuan Sindhi Yoruba 25 33 23 12 14 18 70
These number are allowed to be floats, and the names have to match the names in the reference panel!
A provided .sites file has been included, in the data folder.
It looks like this:
1 752566 G A 1 842013 G T 1 891021 G A 1 903426 C T 1 949654 A G ...
vcf files
Via vcftools [3] it should be possible to convert .vcf files with genotype likelihoods to .beagle files:
vcftools --vcf test.vcf --out test --BEAGLE-GL --chr 1,2
Chromosome has to be specified, try specifying multiple chromosomes.
Or use bcftools' [4] 'query' option for generating a .beagle file from a .vcf file.
Running fastNGSadmix
An example of running a command with fastNGSadmix with a beagle file:
#genotype likelihoods from sequencing file GL=example/yriFrenchHan_depth05.beagle.gz #frequencies from reference panel REF=data/refPanel_humanOrigins_7worldPops.txt #number of individuals in each population of reference panel NIND=data/nInd_humanOrigins_7worldPops.txt ./fastNGSadmix -likes $GL -fname $REF -Nname $NIND -out yriFrenchHan_depth05 -whichPops all
You can pick which populations should be analyzed via the "-whichPops" option, where you write the names of the population comma separated, or "all" if you want to include all populations.
It should also be noted that the program quits if there are duplicate sites (based on chr_pos_A0_A1 (alleles are sorted) ID) in the input or the reference panel.
It then produces two files yriFrenchHan_depth05.qopt with the admixture proportions and a log file yriFrenchHan_depth05.log
Or with a plink file, with also the -whichPops option:
PLINKFILE=example/NA20502_TSI ./fastNGSadmix -plink $PLINKFILE -fname $REF -Nname $NIND -out NA20502_TSI -whichPops French,Han,Yoruba
A whole list of options can be explored by running fastNGSadmix without any input:
./fastNGSadmix
Must be specified:
- -likes [char*] (either -likes or -plink)
Path of Beagle likelihood file of one input individual.
- -plink [char*] (either -likes or -plink)
Path of Plink file of one input individual - in the binary bed format.
- -Nname [char*]
Path of Number of individuals in each reference population file, must have names of the populations.
- -fname [char*]
Path of reference panel with population frequencies. Use "-whichPops" to specify which population to include.
- -whichPops [char*]
This option tells which populations from the reference panel to include in the analysis, must be comma separated (pop1,pop2,..) if "all", all populations in the reference will be included.
Optional:
- -haploid [FLAG] (Only write '-haploid' in command line)
Raise this flag if haploid organism being analyzed, first two cols of beagle file will be used - does not work for plink! Only write "-haploid".
- -out [char*]
Prefix for the output files .qopt and .log.
- -printFreq [int]
This option prints the admixture adjusted allele frequencies of reference panel + input individual. Disabled per default can be enabled by setting this to 1.
- -doAdjust [int]
To use the unadjusted approach set "-doAdjust 0". By default the method adjusting the frequencies is used (see more in the paper),
- -seed [int]
Set seed for random generation of initial guesses in EM and for bootstrap.
- -method [int]
Set this option to enable acceleration of the EM algorithm, enabled per default, set it to 0 for unaccelerated EM.
- -maf [float]
Set this option to filter away sites from the reference panel with lower minor allele frequency in any of the analyzed populations, default value is 0.
- -Qconv [int]
Set this option for faster inference of admixture proportions, the "-Qconv" option can be set to 1, this bases the converge criteria on change in the admixture proportions values. The threshold of this can be set with "-Qtol". This is less precise than the likelihood based convergence.
- -Qtol [float]
Set the threshold for changes in Q values. By default the "-Qtol" threshold is 0.0000001, it should not be put lower than this and generally this is only for when wanting a fast overview of the data, as it is less precise than the likelihood based convergence.
- -tol [float]
Set the tolerance for convergence based on likelihoods, per default 0.00001, can only be adjusted for non accelerated EM, the tolerance threshold is 1e-7 for the accelerated EM algorithm.
- -maxiter [int]
Set maximum number of EM iterations, per default it is set to 2000.
- -conv [int]
Specifies the number of convergence runs, with a new random starting point for each run. This is useful to test for convergence, testing that the same Q values are estimated for each run. The program will execute a maximum of 10 times.
- -boot [int]
Specifies the number of bootstrap runs, where random sites are sampled (with replacement) for each run. This is useful for generating a confidence interval of your estimate, for example doing the empirical 0.025 and 0.975 quantiles. The maximum number of bootstraps is 10000.
- -randomBoot [int]
Specify if the program takes random Q starting points for each bootstrap, instead of using converged upon estimate, default value 0.
Outputs
fastNGSadmix produces two files, a .qopt file, with the estimated admixture proportions, with names of which population they are inferred for. If n number of bootstraps are run, there will be n+2 rows in the .qopt file the first two rows with the names of the populations analyzed and the converged upon estimates, and then n rows with the bootstrapping estimates.
cat NA20502_TSI.qopt French Han Yoruba 1.0000 0.0000 0.0000
And if doing 10 bootstraps:
./fastNGSadmix -plink $PLINKFILE -fname $REF -Nname $NIND -out NA20502_TSI_boot -whichPops French,Han,Yoruba -boot 10 -seed 1
cat NA20502_TSI_boot.qopt French Han Yoruba 1.0000 0.0000 0.0000 1.0000 0.0000 0.0000 1.0000 0.0000 0.0000 1.0000 0.0000 0.0000 1.0000 0.0000 0.0000 1.0000 0.0000 0.0000 1.0000 0.0000 0.0000 1.0000 0.0000 0.0000 1.0000 0.0000 0.0000 1.0000 0.0000 0.0000 1.0000 0.0000 0.0000
Where the first row of numbers is the converged upon estimates, and the rest are the bootstrapping runs.
It also produces a .log file with all the information about the run, as well as the converged upon estimate:
cat NA20502_TSI_boot.log
Input: likes=(null) plink=example/NA20502_TSI Nname=data/nInd_humanOrigins_7worldPops.txt fname=data/refPanel_humanOrigins_7worldPops.txt outfiles=NA20502_TSI_boot
Setup: seed=1 method=1
The accelerated EM has been chosen
The adjusted method has been chosen
Convergence: maxIter=2000 tol=0.00000010
The following number of bootstraps have been chosen: 10
Input has this many sites without missing data 441695
Ref has this many sites 442769
Overlap: of 441695 sites between input and ref
nPop=3
Opening nInd file: data/refPanel_humanOrigins_7worldPops.txt with nPop=3
Chosen pop French
N = 25.000000
Chosen pop Han
N = 33.000000
Chosen pop Yoruba
N = 70.000000
This many iterations 30 for run 0
At this bootstrapping: 1 out of: 10
At this bootstrapping: 2 out of: 10
At this bootstrapping: 3 out of: 10
At this bootstrapping: 4 out of: 10
At this bootstrapping: 5 out of: 10
At this bootstrapping: 6 out of: 10
At this bootstrapping: 7 out of: 10
At this bootstrapping: 8 out of: 10
At this bootstrapping: 9 out of: 10
At this bootstrapping: 10 out of: 10
best like -304264.169822 after 0!
Q 0.999980 Q 0.000010 Q 0.000010 after 0!
Estimated Q = 0.999980 0.000010 0.000010
best like -304264.169822 after 0 runs!
FIRST row of .qopt file is BEST estimated Q, rest are nBoot bootstrapping Qs
[ALL done] cpu-time used = 14.96 sec
[ALL done] walltime used = 15.00 sec
A barplot can easily be generated using R/makeBarplot.R:
Rscript R/makeBarplot.R NA20502_TSI_boot.qopt
PCA
After having run the admixture estimation, a PCA analysis can be run using the estimated admixture proportions to account for population structure in the PCA. The PCA method is implemented in R and is run using the script fastNGSAdmixPCA.R in the R folder. It requires the R package snpStats, it requires the R package parallel, if multi-threaded faster analysis is wanted.
The method needs the estimated admixture proportions, the analyzed beagle or plink file as well as the genotypes of the reference panel used (as binary plink files).
For the provided genotypes the .fam file must have the group or population specified as FID, meaning first column and then the individual ID as the second column:
French F1 0 0 1 1 French F2 0 0 1 1 Yoruba Y1 0 0 2 1 Yoruba Y2 0 0 1 1 Yoruba Y3 0 0 2 1 ...
Arguments should be given like this '-OPTION ARGUMENT', where there is a space between option and the argument.
GL=example/yriFrenchHan_depth05.beagle.gz ./fastNGSadmix -likes $GL -fname $REF -Nname $NIND -out yriFrenchHan_depth05 -whichPops French,Yoruba,Han Rscript R/fastNGSadmixPCA.R -likes $GL -qopt yriFrenchHan_depth05.qopt -out yriFrenchHan_depth05 -ref data/humanOrigins_7worldPops
Or using analyzed plink files.
./fastNGSadmix -plink $PLINKFILE -fname $REF -Nname $NIND -out NA20502_TSI -whichPops French,Yoruba,Han Rscript R/fastNGSadmixPCA.R -plinkFile $PLINKFILE -qopt NA20502_TSI.qopt -out NA20502_TSI -ref data/humanOrigins_7worldPops
By default the populations of the analyzed *.qopt file are used. It can be specified which PCs should be plotted (default PC1 and PC2), by for example giving -PCs 2,3 as an argument (if PC2 and PC3), the PCA plot is plotted as *_PCAplot.pdf. The script also generates barplots (*_quantile_admixBarplot.png) of the admixture proportions supplied with confidence intervals (in the case of bootstraps), as well as generating files of the covariance matrix (*_covar.txt), eigenvectors (*_eigenvecs.txt) and eigenvalues (*_eigenvals.txt) used for the PCA plot. Plotting can be disabled by setting -doPlots 0.
The covariance matrix and eigenvectors have the ids of the individuals, the input individual has the id "SAMPLE". An *_indi.txt file is also created with the individual id and population/groupd id of each individual.
For faster analyses when analysing multiple individuals with the same reference panel of genotypes, the option saveCovar can be set -saveCovar 1, this stores the covariance matrix for the reference panel indivduals, this will be a lot quicker for big reference panel. Also one can specify -multiCores X, where parallel:::mcapply will then be used with X cores, this is require the package snpStats.
The PCA can also be done using the method from Fumagalli et al., 2013, that is done by setting '-ngsTools 1'. The method can also be run only using the estimated prior, that means calculating the genotype probability given only the admixture proportions and population specific frequencies and NOT the sequencing data, for more see supplementaries of Jørsboe et al., 2017. This is done by setting '-onlyPrior 1'. Setting '-overlapRef 0' means that the covariance matrix for the reference panel genotypes is calculated using all markers in the reference genotypes (per default only overlapping markers are used for calculating the reference panel covariance matrix).
To find out which options there are run wihtout any input:
Rscript R/fastNGSadmixPCA.R
The snpStats R package, is needed for running this script.
For combining these analyses plotting a joint PCA plot, procrustes analyses can be done, there is a script for this:
Rscript R/doProcastesPCA.R
Procrustes analysis is needed as the calculated covariance matrix differs, when the inputted sample differs. One should be careful doing procrustes if the PCA plots differ much, for instance if different number of sites are used for the ref genotypes.
It will take a directory of *_covar.txt files. The .fam file of the plink files used for the reference genotypes. Then a file with all the samples/files you want to include in the analyses, these should fit with the name of the *_covar.txt files (excluding _covar.txt), these should also all be run with the same reference panel and they should all be run with the same populations. And then lastly an argument specifying the name of the group being procrusted.
The vegan and parallel R packages are needed.
Here it can be done like this:
echo -e "yriFrenchHan_depth05\nNA20502_TSI" > procrustes.list Rscript R/doProcastesPCA.R . data/humanOrigins_7worldPops.fam procrustes.list Samples
Easy Example for analysing multiple beagle files
First we create a new folder with the .beagle files we want to analyse:
mkdir test cp example/NA12763.mapped.ILLUMINA.bwa.CEU.low_coverage.20130502.bam.beagle.gz test/ cp example/yriFrenchHan_depth05.beagle.gz test/
For analysing all the *.beagle files in a directory using a one-line command one can use xargs:
find test/ | egrep "beagle$|beagle.gz$" | xargs -I % bash -c './fastNGSadmix -likes % -fname data/refPanel_humanOrigins_7worldPops.txt -Nname data/nInd_humanOrigins_7worldPops.txt -out `basename %` -whichPops all;'
fastNGSadmix can also be run in parallel giving the argument -P X to xargs, where X denotes the number of cores. Be aware that this command will analyse ALL *.beagle or *.beagle.gz files in the "test/" directory, in order to see which files will be analysed run this command first:
find test/ | egrep "beagle$|beagle.gz$"
And then for doing the PCA (also analyses all *.beagle or *.beagle.gz files in the specified directory):
find test/ | egrep "beagle$|beagle.gz$" | xargs -I % bash -c 'Rscript R/fastNGSadmixPCA.R -likes % -ref data/humanOrigins_7worldPops -qopt `basename %`.qopt -out test/`basename %`'
The PCA plots can be compared using procrustes analysis, using R/doProcastesPCA.R.
ls test/ | grep "_covar.txt$" | sed -e 's/_covar.txt//g' > procrustes.list Rscript R/doProcastesPCA.R test/ data/humanOrigins_7worldPops.fam procrustes.list Samples
Making a reference panel
The program needs a reference panel of population specific frequencies, populations for which admixture proportions should be estimated for. This reference panel can be derived from 1000 Genomes or HGDP, or other data. The program also needs a file telling the size of each reference panel population. It is recommended to run a missingness filter on the ref data first, as it is implicitly assumed in the model that there is full data for most sites, for example in plink using --geno 0.05. There is an R script plinkToRef.R in the R folder, that can convert a plink file to a reference panel fine and size of reference populations file.
Example:
Rscript R/plinkToRef.R data/humanOrigins_7worldPops
This generates 3 files, a reference panel named data/refPanel_humanOrigins_7worldPops.txt, a number of individuals file called data/nInd_humanOrigins_7worldPops.txt and a data/humanOrigins_7worldPops.sites file with chr pos and minor major alleles, for being used with ANGSD. The provided .fam file must have the group or population specified as FID, meaning first column and then the individual ID as the second column:
French F1 0 0 1 1 French F2 0 0 1 1 Yoruba Y1 0 0 2 1 Yoruba Y2 0 0 1 1 Yoruba Y3 0 0 2 1 ...
A second argument can be given, telling if duplicate sites (based on chr_pos_A0_A1 ID) should be removed from the reference panel created, this argument has to be 0 (no) or 1 (yes), 1 pr. default.
Rscript R/plinkToRef.R data/humanOrigins_7worldPops 1
A third argument of a MAF threshold can be supplied, meaning all sites where the MAF is below this value is removed.
Rscript R/plinkToRef.R data/humanOrigins_7worldPops 1 0.05
The snpStats R package, is also needed for running the plinkToRef.R script.
It should be noted that reading in plink files especially big ones might require a lot of RAM, why doing it on a server might be preferable!
Making a reference panel with NGSadmix
We can use NGSadmix for generating a reference panel, which we can then use for fast analysis of individual samples.
For instance if we want to construct a reference panel from 1000 genomes beagle data of (CEU, CHB, JHB and YRI), we first do an analysis of the reference panel with NGSadmix using K=3:
The beagle file should have chr_pos as marker ID for this to work!
gunzip example/ceuChbJhbYriChr1.beagle.gz ./ngsadmix32 -likes example/ceuChbJhbYriChr1.beagle -K 3 -printInfo 1 -outfiles ceuChbJhbYriChr1.beagle
We can then create a reference panel from the .fopt.gz file, with three populations:
echo "id chr pos name A0_freq A1 K1 K2 K3" > refPanel_ceuChbJhbYriChr1.beagle.txt cut -f1,2,3 example/ceuChbJhbYriChr1.beagle > tmp.beagle gunzip ceuChbJhbYriChr1.beagle.fopt.gz ## find overlap of .beagle and .filter files and creates tmp.ref with 6 first column of refPanel_ Rscript example/makeRefNGSadmix.R tmp.beagle ceuChbJhbYriChr1.beagle.filter paste tmp.ref ceuChbJhbYriChr1.beagle.fopt >> refPanel_ceuChbJhbYriChr1.beagle.txt rm tmp.ref tmp.beagle gzip ceuChbJhbYriChr1.beagle.fopt
And then we can create a number of individuals in each reference population file, summing each column of the .qopt file:
echo "K1 K2 K3" > nInd_ceuChbJhbYriChr1.beagle.txt paste -d" " <(cut -f1 -d" " ceuChbJhbYriChr1.beagle.qopt | paste -sd+ | bc) <(cut -f2 -d" " ceuChbJhbYriChr1.beagle.qopt | paste -sd+ | bc) <(cut -f3 -d" " ceuChbJhbYriChr1.beagle.qopt | paste -sd+ | bc) >> nInd_ceuChbJhbYriChr1.beagle.txt
Then we can run fastNGSadmix using this reference panel
#genotype likelihoods from sequencing file GL=example/NA12763.mapped.ILLUMINA.bwa.CEU.low_coverage.20130502.bam.beagle.gz ./fastNGSadmix -likes $GL -fname refPanel_ceuChbJhbYriChr1.beagle.txt -Nname nInd_ceuChbJhbYriChr1.beagle.txt -out NA12763_CEU_K3 -whichPops all
The only thing that should be changed when building a reference panel, from another beagle file, is the value of 'K' and the K columns in the reference panel and nInd file should be set accordingly.
Making a reference panel with ADMIXTURE
A reference panel can be generated making use of ADMIXTURE.
This could for instance be a clustered analysis, where instead of inferring admixture proportions for populations, we instead estimate admixture proportions for clusters of populations.
For instance if we want to detect the african and non-african component on an individual, we might first do an analysis of the reference panel with ADMIXTURE using K=2:
admixture data/humanOrigins_7worldPops.bed 2
It might be necessary to run ADMIXTURE more than once in order to make sure that it has converged.
We can then create a reference panel from the .P file, with two populations:
echo "id chr pos name A0_freq A1 K1 K2" > refPanel_humanOrigins_7worldPops.2.P.txt
paste -d" " <( awk -F " " ' { print $1"_"$4, $1, $4, $2, $6, $5 } ' data/humanOrigins_7worldPops.bim ) humanOrigins_7worldPops.2.P >> refPanel_humanOrigins_7worldPops.2.P.txt
And then we can create a number of individuals in each reference population file, summing each column of the .Q file:
echo "K1 K2" > nInd_humanOrigins_7worldPops.2.Q.txt paste -d" " <(cut -f1 -d" " humanOrigins_7worldPops.2.Q | paste -sd+ | bc) <(cut -f2 -d" " humanOrigins_7worldPops.2.Q | paste -sd+ | bc) >> nInd_humanOrigins_7worldPops.2.Q.txt
Then we can run fastNGSadmix using this reference panel
#genotype likelihoods from sequencing file GL=example/NA12763.mapped.ILLUMINA.bwa.CEU.low_coverage.20130502.bam.beagle.gz ./fastNGSadmix -likes $GL -fname refPanel_humanOrigins_7worldPops.2.P.txt -Nname nInd_humanOrigins_7worldPops.2.Q.txt -out NA12763_CEU_K2 -whichPops all
If you run fastNGSadmix multiple times, and then would like to cat the .Q file from ADMIXTURE with your results you can do the following, assuming your files are named NAXXXXX_CEU_K2.qopt:
cat <(cat humanOrigins_7worldPops.2.Q ) <(find . | grep "NA....._CEU_K2.qopt$" | xargs -I % bash -c 'cat % | head -n 2 | tail -n 1;') > humanOrigins_7worldPops_fastNGSadmix.2.Q
To find out which files from fastNGSadmix have been been merged with the ADMIXTURE .Q file and in which order, run:
find . | grep "NA....._CEU_K2.qopt$"
The only thing that should be changed for constructing another reference panel, is the value of 'K' and K columns of the reference panel should be set accordingly.
Human reference panels
fastNGSadmix already comes with a premade reference panel.
refPanel_humanOrigins_7worldPops
This reference panel is made from Lazaridis et al. (2014) where the curated Human Origins dataset was selected.
The dataset was lifted to hg19 using the program liftOver, The SNPs then got rs IDs, using 1000 Genomes data.
Furthermore sites with more than 5 % missing and a MAF below 5 % were removed, and only autosomal sites were kept. 7 populations were selected French, Han, Chukchi, Karitiana, Papuan, Sindhi and Yoruba to have representation for most of the world. Furthermore it was made sure that there were only unadmixed individuals within each population.
Creating other reference panels
Other reference panels can easily be made from the Human Origins data:
wget popgen.dk/software/download/fastNGSadmix/humanOrigins_ALL.tar.gz tar -xzf humanOrigins_ALL.tar.gz
Where I have updated the ids with the group/population-id for the FID and individual-id as the IID as well as the snp-ids as chr_pos. Then I would recommend to keep only autosomal sites and do a MAF and geno filter (one can play around with these filters!):
plink --bfile humanOrigins_ALL/humanOrigins_ALL --make-bed --out humanOrigins_ALL/humanOrigins_ALLV2 --maf 0.05 --geno 0.05 --chr 1-22
And then the desired individuals/populations can be extracted from the humanOrigins_ALL/humanOrigins_ALLV2.* plink files using plink, using the --keep option. And it can then easily be turned into a reference panel using plinkToRef.R. (Please notice that this command takes up a lot of RAM!)
Rscript R/plinkToRef.R humanOrigins_ALL/humanOrigins_ALLV2
Other reference panels can be created for instance based on 1000 Genomes or own data, first convert the data into plink files and then use the script plinkToRef.R.
Citation
The program has been published in Bioinformatics:
Emil Jørsboe, Kristian Hanghøj, Anders Albrechtsen; fastNGSadmix: Admixture proportions and principal component analysis of a single NGS sample, Bioinformatics, btx474, https://doi.org/10.1093/bioinformatics/btx474
Log
version 1.0: Stable initially released version of fastNGSadmix - 07-09-2017
version 1.1: Beagle input and reference panel do not have to be sorted genomically, only find overlap based on ID - 23-09-2017
version 1.2: Can also handle and analyse haploid data, only with beagle input files - 03-02-2018