Tambetsin paperin näytteiden datan saa ladatua täältä:
https://evolbio.ut.ee/Tambets2018. Tiedostot saa yhdistettyä PLINKin kanssa:
wget evolbio.ut.ee/Tambets2018/Tambets_et_al_2018_{1M,660k,650k,610k}.{bed,bim,fam}
printf %s\\n Tambets_et_al_2018_{660k,610k}>mergelist
plink --allow-no-sex --bfile Tambets_et_al_2018_1M --merge-list mergelist --make-bed --out yhdistetty
Tambets_et_al_2018_650k ei suostunut yhdistymään muiden kanssa:
$ plink --allow-no-sex --bfile Tambets_et_al_2018_1M --bmerge Tambets_et_al_2018_650k --make-bed --out testi
[...]
2 more multiple-position warnings: see log file.
Error: 167034 variants with 3+ alleles present.
* If you believe this is due to strand inconsistency, try --flip with
testi.missnp.
Sen sai yhdistettyä näin:
a=yhdistetty;b=Tambets_et_al_2018_650k
plink --allow-no-sex --bfile $a --bmerge $b --make-bed --out merge
plink --bfile $b --flip merge-merge.missnp --make-bed --out $b-flip
plink --allow-no-sex --bfile $a --bmerge $b-flip --make-bed --out merge
plink --bfile $b-flip --exclude merge-merge.missnp --make-bed --out $b-filtered
plink --bfile $a --bmerge $b-filtered --make-bed --out yhdistetty2 --allow-no-sex
Latasin PLINK 1.9:n tältä sivulta:
https://www.cog-genomics.org/plink/. Alkuperäisellä Harvardin sivulla olevan PLINK 1.9:n Mac-binaarit eivät toimineet koneellani. PLINK 2:ssa piti korvata `--bmerge` asetuksella `--pmerge`, mutta senkin jälkeen sain seuraavan virheen, koska PLINK 2 käyttää uutta
pgen-pvar-psam -formaattia: `Error: Failed to open Tambets_et_al_2018_650k.psam`.
Kokeilin sitten tehdä PCA:n kummallakin PLINKillä ja SmartPCA:lla:
x=yhdistetty2;plink --bfile $x --indep-pairwise 50 10 .5 --make-bed --out $x.2;plink --bfile $x.2 --extract $x.2.prune.in --geno .1 --make-bed --out $x.3;plink --bfile $x.3 --pca
curl -Ls pastebin.com/raw/md6FrMM4|tr -d \\r>popnames # taulukosta S1B
awk 'NR==FNR{a[$1]=$2;next}{$1=a[$1]}1' popnames plink.eigenvec>plink.eigenvec.2
f=yhdistetty2.3;smartpca -p <(printf %s\\n genotypename:\ $f.bed snpname:\ $f.bim indivname:\ $f.fam evecoutname:\ $f.evec evaloutname:\ $f.eval numoutlieriter:\ 0 autoshrink:\ YES)
awk 'NR==FNR{a[$1]=$2;next}{$1=a[$1]":"$1};NF--' popnames <(sed 1d $f.evec|sed 's/[^:]*://')>$f.evec.2
En osaa vielä tehdä filtteröintiä tai imputointia kunnolla, ja tein tässä varmaan muutenkin jotain väärin. Ainakin se on epäilyttävää, että yksi burjatti on suunnilleen samassa kohdassa PC1:llä kuin yksi tataari. Mansien sijoittuminen vaikuttaa kuitenkin vastaavan Tambetsin ADMIXTURE-analyysiä, jossa mansien siperialainen perimä vaihteli suunnilleen hantien tasosta suomalaisten tasoon.
Yllä SmartPCA:n kuvaajassa on korkeampi PC1:n ja PC2:n selittämän varianssin prosentti, koska SmartPCA palautti 10 ominaisvektoria mutta `plink --pca` palautti 20.
Tein kummankin SmartPCA:n ja `plink --pca`:n kohdalla klusteroinnin niin, että kerroin ensin sarakkeet ominaisarvojen neliöjuurilla, ja otin sitten klusteroinnissa huomioon vain ensimmäiset kuusi ulottuvuutta. G25:ssäkin klusterointi menee usein sekaisin skaalaamattomilla koordinaateilla.
- Koodi: Valitse kaikki
library(tidyverse)
t=read.table("plink.eigenvec.2")
eig=as.double(readLines("plink.eigenval"))
pct=paste0("PC",seq(length(eig))," (",sprintf("%.1f",100*eig/sum(eig)),"%)")
t=cbind(t[,c(1,2)],t(t(t[,-c(1,2)])*sqrt(eig)))
t$k=cutree(hclust(dist(t[,c(3:8)])),k=16)
t$lab=paste0(t[,1],":",t[,2])
ggplot(t,aes(x=V3,y=V4))+
geom_polygon(data=t%>%group_by(k)%>%slice(chull(V3,V4)),alpha=.2,aes(color=as.factor(k),fill=as.factor(k)),size=.3)+
geom_point(aes(color=as.factor(k)),size=.5)+
geom_text(aes(label=lab,color=as.factor(k)),size=2,vjust=-.7)+
ggtitle("plink --pca")+
labs(x=pct[1],y=pct[2])+
scale_x_continuous(breaks=seq(-10,10,.1),expand=expansion(mult=.12))+
scale_y_continuous(breaks=seq(-10,10,.1),expand=expansion(mult=.03))+
scale_color_manual(values=hcl(head(seq(15,375,length=length(unique(t$k))+1),-1),90,50))+
theme(
aspect.ratio=1/sqrt(2),
axis.text=element_text(color="black",size=6),
axis.ticks.length=unit(0,"pt"),
axis.ticks.x=element_blank(),
axis.ticks.y=element_blank(),
axis.title=element_text(color="black",size=8),
legend.position="none",
panel.background=element_rect(fill="white"),
panel.border=element_rect(color="gray75",fill=NA,size=.4),
panel.grid.major=element_line(color="gray70",size=.15),
plot.title=element_text(size=10)
)
ggsave("a.png")
system("mogrify -trim -bordercolor white -border 16x16 a.png")