Visualizing an ADMIXTURE run as a polygonal diagram

[vbnetkhio]

the calculation on (most of) evolbio's database has finished:

fst *1000 version:
https://pastebin.com/raw/UzeqH7Dr

s.dev * 1000000 version:
https://pastebin.com/raw/3HDjLwGU

phylip version
https://pastebin.com/raw/NxZUj10x

sadly, i messed up the pop numbers, so now there's a Corsico-Croat and a Germano-Greek population

the rest should be fine

[Komintasavalta]

would there be any problems with calculating fst this way:

run supervised admixture, assign each sample to it's population, and run with as many K as there is populations, and fst gets written to the output.

http://dalexander.github.io/admixtur...ure-manual.pdf

smartpca's version takes forever, this could actually be faster? because there are no unassigned samples, just the allele frequencies and fst will be calculated

Did you do LD pruning? The paper by Kerminen et al. said that they only used about 60,000 SNPs to calculate FST, even though their paper was about Finnish subpopulations.

The ADMIXTURE manual says this:

2.4 How many markers do I need to supply to ADMIXTURE?

This depends on how genetically differentiated your populations are, and on what you plan to do with the estimates. It has been noted elsewhere [4] that the number of markers needed to resolve populations in this kind of analysis is inversely proportional to the genetic distance (FST) betweeen the populations.

It is also noted in that paper that more markers are needed to perform adequate GWAS correction than are needed to simply observe the population structure.

As a rule of thumb, we have found that 10,000 markers suffice to perform GWAS correction for continentally separated populations (for example, African, Asian, and European populations FST > .05) while more like 100,000 markers are necessary when the populations are within a continent (Europe, for instance, FST < 0.01).

Using supervised ADMIXTURE to calculate FST actually works, but it's slower than SmartPCA because the five priming steps take a long time:

$ cut -d' ' -f2 uralic.pick|awk '!a[$0]++'>uralic.pop
$ admixture -j4 --supervised uralic.bed 13
****                   ADMIXTURE Version 1.3.0                  ****
****                    Copyright 2008-2015                     ****
****           David Alexander, Suyash Shringarpure,            ****
****                John  Novembre, Ken Lange                   ****
****                                                            ****
****                 Please cite our paper!                     ****
****   Information at www.genetics.ucla.edu/software/admixture  ****

Parallel execution requested.  Will use 4 threads.
Random seed: 43
Point estimation method: Block relaxation algorithm
Convergence acceleration algorithm: QuasiNewton, 3 secant conditions
Point estimation will terminate when objective function delta < 0.1
Estimation of standard errors disabled; will compute point estimates only.
Supervised analysis mode.  Examining .pop file...
Size of G: 181x597573
Performing five EM steps to prime main algorithm
1 (EM) 	Elapsed: 18.302	Loglikelihood: -7.38566e+07	(delta): 1.40226e+08
2 (EM) 	Elapsed: 18.246	Loglikelihood: -7.38526e+07	(delta): 3976.64
3 (EM) 	Elapsed: 20.037	Loglikelihood: -7.38526e+07	(delta): 0.176208
4 (EM) 	Elapsed: 20.448	Loglikelihood: -7.38526e+07	(delta): 0.0233383
5 (EM) 	Elapsed: 20.357	Loglikelihood: -7.38526e+07	(delta): 0.015767
Initial loglikelihood: -7.38526e+07
Starting main algorithm
1 (QN/Block) 	Elapsed: 12.214	Loglikelihood: -7.38526e+07	(delta): 0
Summary:
Converged in 1 iterations (113.629 sec)
Loglikelihood: -73852623.464629
Fst divergences between estimated populations:
	Pop0	Pop1	Pop2	Pop3	Pop4	Pop5	Pop6	Pop7	Pop8	Pop9	Pop10	Pop11
Pop0
Pop1	0.043
Pop2	0.041	0.015
Pop3	0.043	0.023	0.021
Pop4	0.062	0.040	0.039	0.042
Pop5	0.056	0.046	0.042	0.041	0.061
Pop6	0.047	0.017	0.018	0.026	0.043	0.050
Pop7	0.117	0.083	0.086	0.096	0.111	0.121	0.086
Pop8	0.053	0.024	0.025	0.033	0.050	0.055	0.026	0.091
Pop9	0.075	0.049	0.049	0.056	0.073	0.076	0.050	0.113	0.055
Pop10	0.066	0.033	0.035	0.044	0.061	0.069	0.035	0.096	0.040	0.064
Pop11	0.069	0.034	0.037	0.047	0.064	0.074	0.036	0.097	0.041	0.065	0.048
Pop12	0.175	0.140	0.142	0.154	0.168	0.179	0.139	0.192	0.142	0.160	0.148	0.144
Writing output files.

You can also use ADMIXTOOLS 2 to calculate FST, but it's slower than SmartPCA (the `f2m` function converts FST or f2 pairs to a square matrix):

$ R -e 'library("admixtools");f2m=function(x){t=as.data.frame(x[,1:3]);t2=rbind(t,setNames(t[,c(2,1,3)],names(t)));xtabs(t2[,3]~t2[,2]+t2[,1])};fst=fst("v44.3_HO_public",c("Besermyan","Enets","Estonian","Finnish","Hungarian","Karelian","Mansi","Mordovian","Nganasan","Saami.DG","Selkup","Udmurt","Veps"));write.csv(f2m(fst),"fst",quote=F)'

Without LD pruning, calculating FST for the 13 populations listed in my previous post took 28 seconds with SmartPCA, 52 seconds with ADMIXTOOLS 2, and 113 seconds with ADMIXTURE. After I ran `--indep-pairwise 50 10 .1`, it only took about 3 seconds with SmartPCA.

[vbnetkhio]

Did you do LD pruning? The paper by Kerminen et al. said that they only used about 60,000 SNPs to calculate FST, even though their paper was about Finnish subpopulations.

The ADMIXTURE manual says this:

2.4 How many markers do I need to supply to ADMIXTURE?

This depends on how genetically differentiated your populations are, and on what you plan to do with the estimates. It has been noted elsewhere [4] that the number of markers needed to resolve populations in this kind of analysis is inversely proportional to the genetic distance (FST) betweeen the populations.

It is also noted in that paper that more markers are needed to perform adequate GWAS correction than are needed to simply observe the population structure.

As a rule of thumb, we have found that 10,000 markers suffice to perform GWAS correction for continentally separated populations (for example, African, Asian, and European populations FST > .05) while more like 100,000 markers are necessary when the populations are within a continent (Europe, for instance, FST < 0.01).

Using supervised ADMIXTURE to calculate FST actually works, but it's slower than SmartPCA because the five priming steps take a long time:

$ cut -d' ' -f2 uralic.pick|awk '!a[$0]++'>uralic.pop
$ admixture -j4 --supervised uralic.bed 13
****                   ADMIXTURE Version 1.3.0                  ****
****                    Copyright 2008-2015                     ****
****           David Alexander, Suyash Shringarpure,            ****
****                John  Novembre, Ken Lange                   ****
****                                                            ****
****                 Please cite our paper!                     ****
****   Information at www.genetics.ucla.edu/software/admixture  ****

Parallel execution requested.  Will use 4 threads.
Random seed: 43
Point estimation method: Block relaxation algorithm
Convergence acceleration algorithm: QuasiNewton, 3 secant conditions
Point estimation will terminate when objective function delta < 0.1
Estimation of standard errors disabled; will compute point estimates only.
Supervised analysis mode.  Examining .pop file...
Size of G: 181x597573
Performing five EM steps to prime main algorithm
1 (EM) 	Elapsed: 18.302	Loglikelihood: -7.38566e+07	(delta): 1.40226e+08
2 (EM) 	Elapsed: 18.246	Loglikelihood: -7.38526e+07	(delta): 3976.64
3 (EM) 	Elapsed: 20.037	Loglikelihood: -7.38526e+07	(delta): 0.176208
4 (EM) 	Elapsed: 20.448	Loglikelihood: -7.38526e+07	(delta): 0.0233383
5 (EM) 	Elapsed: 20.357	Loglikelihood: -7.38526e+07	(delta): 0.015767
Initial loglikelihood: -7.38526e+07
Starting main algorithm
1 (QN/Block) 	Elapsed: 12.214	Loglikelihood: -7.38526e+07	(delta): 0
Summary:
Converged in 1 iterations (113.629 sec)
Loglikelihood: -73852623.464629
Fst divergences between estimated populations:
	Pop0	Pop1	Pop2	Pop3	Pop4	Pop5	Pop6	Pop7	Pop8	Pop9	Pop10	Pop11
Pop0
Pop1	0.043
Pop2	0.041	0.015
Pop3	0.043	0.023	0.021
Pop4	0.062	0.040	0.039	0.042
Pop5	0.056	0.046	0.042	0.041	0.061
Pop6	0.047	0.017	0.018	0.026	0.043	0.050
Pop7	0.117	0.083	0.086	0.096	0.111	0.121	0.086
Pop8	0.053	0.024	0.025	0.033	0.050	0.055	0.026	0.091
Pop9	0.075	0.049	0.049	0.056	0.073	0.076	0.050	0.113	0.055
Pop10	0.066	0.033	0.035	0.044	0.061	0.069	0.035	0.096	0.040	0.064
Pop11	0.069	0.034	0.037	0.047	0.064	0.074	0.036	0.097	0.041	0.065	0.048
Pop12	0.175	0.140	0.142	0.154	0.168	0.179	0.139	0.192	0.142	0.160	0.148	0.144
Writing output files.

You can also use ADMIXTOOLS 2 to calculate FST, but it's slower than SmartPCA (the `f2m` function converts FST or f2 pairs to a square matrix):

$ R -e 'library("admixtools");f2m=function(x){t=as.data.frame(x[,1:3]);t2=rbind(t,setNames(t[,c(2,1,3)],names(t)));xtabs(t2[,3]~t2[,2]+t2[,1])};fst=fst("v44.3_HO_public",c("Besermyan","Enets","Estonian","Finnish","Hungarian","Karelian","Mansi","Mordovian","Nganasan","Saami.DG","Selkup","Udmurt","Veps"));write.csv(f2m(fst),"fst",quote=F)'

Without LD pruning, calculating FST for the 13 populations listed in my previous post took 28 seconds with SmartPCA, 52 seconds with ADMIXTOOLS 2, and 113 seconds with ADMIXTURE. After I ran `--indep-pairwise 50 10 .1`, it only took about 3 seconds with SmartPCA.

you can use the multithreaded mode, the priming steps ar much faster. e.g. you can add the flag -j8 for 8 threads.

I decided to try without LD first.
in studies pruning is usually very light, or none. a higher number of SNPs is considered much more valuable.

[vbnetkhio]

pca from the phylip fst file:


non-metric MDS, reflects geography much better:


i used the "user-supplied distance" option, it tells the program that the data is already a distance matrix, and not raw data.
the Samaritans are outliers for some reason. "and" are Amerindians from the andes.

UPGMA clustering, also with user supplied distance

[Komintasavalta]

sadly, i messed up the pop numbers, so now there's a Corsico-Croat and a Germano-Greek population

Sorry, my code printed the population names in the wrong order. I edited my post on the previous page.

pca from the phylip fst file:

You can do classical MDS instead because the input is already a distance matrix: `cmdscale(as.dist(csv.read("input.fst",row.names=1)))`.

`cmdscale(dist(df))` produces identical coordinates to `prcomp(df)$x` (except the signs of some dimensions may arbitrarily be flipped).

[vbnetkhio]

Sorry, my code printed the population names in the wrong order. I edited my post on the previous page.



You can do classical MDS instead because the input is already a distance matrix: `cmdscale(as.dist(csv.read("input.fst",row.names=1)))`.

`cmdscale(dist(df))` produces identical coordinates to `prcomp(df)$x` (except the signs of some dimensions may arbitrarily be flipped).

it's not your fault, i did it with vlookup in libreoffice.

[Lucas]

Its the new way to see genetic, very interesting!

Yes, Komin brins us everyday new beatiful visualization scripts. I wonder what he will show for a year.

[Komintasavalta]

Here's ADMIXTURE runs of Ural-Altaic populations. The clustering and nearest neighbors are still based on a matrix made by concatenating the columns of the admixture weights at different K values. I also tried calculating the clustering and nearest neighbors based on FST distances returned by SmartPCA, but it made less sense.







From these runs, you can see that Tofalars have a lot of Nganasan-like ancestry. They are a people related to Tuvans who live north of the Tuvan Republic in the Eastern Sayan. Todzins are a northern subgroup of Tuvans who also have Samoyedic ancestry (https://ru.wikipedia.org/wiki/Тувинцы-тоджинцы): "As the researchers believe, the Todzhins formed as a result of the mixing of the ancient Keto-speaking and Samoyedic population with the Dubo tribes who moved to the Sayan-Altai region [7] [8]."

Zabolotniye Tatars (Swamp Tatars) are a northern subgroup of Siberian Tatars who live in the swamp regions between the Siberian Tatar and Ob-Ugric territories. In these runs, Swamp Tatars were closer to Mansi than to other Siberian Tatars. According to physical anthropology, Swamp Tatars are also similar to Khanty (https://ru.wikipedia.org/wiki/Сибирские_татары):

- the Uralic type is the main one for all groups of Siberian Tatars occupying the northern area of ​​their residence, and as a component is traced in more southern groups of Siberian Tatars.
- the South Siberian type is characteristic primarily of the Turks of the Barabinsk steppe and as an admixture is noted in almost all Siberian Tatars, with a tendency to increase in the southern, steppe groups and to decrease in the northern, forest groups.

The Central Asian type was recorded among the Barabinians. Some groups of Tobolsk and Tomsk Tatars have the Chulym type. The Zabolotnye Tatars are extremely close to the Berezovsky Khanty [4].

[Komintasavalta]

Yes, Komin brins us everyday new beatiful visualization scripts. I wonder what he will show for a year.

Yeah I just started using R 3 months ago, but I already have 230 files in my directory for R scripts.

I now also made a script for making a stacked bar chart of an ADMIXTURE run which uses the same clustering method as the scripts in this thread, where the clustering is based on a combined matrix of ADMIXTURE runs at different K values. The colors of the population labels are based on cutting the clustering tree in 12 parts.

library(tidyverse)
library(ggdendro)
library(vegan)
library(colorspace)
library(cowplot)

t=read.table("https://pastebin.com/raw/FEwYnBNb",row.names=1) # population averages from ADMIXTURE with population names in first column

t=t[,c(3,2,5,1,4)] # reorder columns (change if your input does not have five columns)
names(t)=paste0("V",1:ncol(t))

# do clustering based on a combined matrix of admixture weights at different K values
joined=sapply(3:8,function(i)read.table(paste0("uralaltaic.",i,".ave"))[,-1])%>%do.call(cbind,.)%>%set_rownames(rownames(t))
hc=hclust(dist(joined))
hc=reorder(hc,-as.matrix(t)%*%seq(ncol(t))^2)
dist=as.matrix(dist(joined))
maxdist=which(dist==max(dist))[1]
hc=reorder(hc,dist[,maxdist%%nrow(dist)]-dist[,maxdist%/%nrow(dist)+1])

# fst=read.csv("https://pastebin.com/raw/ktMkDf24",row.names=1)[rownames(t),rownames(t)]
# fst[fst<0]=0
# maxfst=which(fst==max(fst))[1] # reorder branches based on distance to the pair of populations with the highest FST distance
# hc=reorder(hclust(as.dist(fst)),fst[,maxfst%%nrow(fst)]-fst[,maxfst%/%nrow(fst)+1])
# # hc=reorder(hclust(dist(t)),-as.matrix(t)%*%exp(seq(ncol(t)))) # reorder branches based on the order of the bars

k=as.factor(cutree(hc,12))[hc$labels[hc$order]]

tree=ggdendro::dendro_data(as.dendrogram(hc),type="triangle")

p1=ggplot(ggdendro::segment(tree))+
geom_segment(aes(x=y,y=x,xend=yend,yend=xend),size=.5,lineend="round",color="gray85")+ # `lineend="round"` draws corners properly when not using `type="triangle"`
scale_x_continuous(expand=expansion(mult=c(0,.01)))+ # don't crop a few pixels from the right border of the tree
scale_y_continuous(limits=.5+c(0,nrow(t)),expand=c(0,0))+
theme(
  axis.text=element_blank(),
  axis.ticks=element_blank(),
  axis.ticks.length=unit(0,"pt"), # remove extra space normally occupied by tick marks
  axis.title=element_blank(),
  panel.background=element_rect(fill="gray30"),
  panel.grid=element_blank(),
  plot.background=element_rect(fill="gray30",color=NA), # `color=NA` removes a thin white border around the plot
  plot.margin=margin(5,5,5,0)
)

t=t[hc$labels[hc$order],]
t2=data.frame(V1=rownames(t)[row(t)],V2=colnames(t)[col(t)],V3=unname(do.call(c,t))) # an alternative to `pivot_longer` and `melt`
lab=round(100*t2$V3)
lab[lab<=1]=""

pal1=colorspace::hex(HSV(c(30,210,250,310,0),.4,.9))
pal2=colorspace::hex(HSV(c(30,210,250,310,0),.4,.2))
pal3=hex(HSV(seq(0,360,length.out=n_distinct(k)+1)%>%head(-1),.4,1))

p2=ggplot(t2,aes(x=factor(V1,level=rownames(t)),y=V3,fill=V2))+
geom_bar(stat="identity",width=1,position=position_fill(reverse=T),size=.2,color="gray10")+
geom_text(aes(label=lab),position=position_stack(vjust=.5,reverse=T),size=3.5,color="gray10")+
coord_flip()+
scale_x_discrete(expand=c(0,0))+
scale_y_discrete(expand=c(0,0))+
scale_fill_manual(values=pal1)+
theme(
  axis.text=element_text(color=pal3[k],size=11),
  axis.text.x=element_blank(),
  axis.ticks=element_blank(),
  axis.title=element_blank(),
  legend.position="none",
  plot.background=element_rect(fill="gray30",color=NA),
  panel.background=element_rect(fill="gray30"),
  plot.margin=margin(5,0,5,5)
)

cowplot::plot_grid(p2,p1,rel_widths=c(1,.4))
ggsave("a.png",height=.27*nrow(t),width=7)

[Lucas]

Yeah I just started using R 3 months ago, but I already have 230 files in my directory for R scripts.

I now also made a script for making a stacked bar chart of an ADMIXTURE run which uses the same clustering method as the scripts in this thread, where the clustering is based on a combined matrix of ADMIXTURE runs at different K values. The colors of the population labels are based on cutting the clustering tree in 12 parts.

I think your scripts should be used in academic papers. There are much better than some of their visualization of Admixture charts or PCA plots which looks like shit from 20 years ago graphically.