Rを使ってSNPデータ入りのvcfファイルからDAPCをする方法を教わったのでメモ。
解析にはRstudioを使っています。
デモデータは202座位33個体、参照URLの閲覧日は本記事を作成した2020年12月19日現在のものとします。
ちなみに、DAPCとPCAのインプットは各遺伝子座あたり1SNPにしないといけないので、事前に各遺伝子座から１SNPをランダム抽出しておきましょう。
最初にvcfをDAPCに使えるデータ形式に変換する
> install.packages("vcfR")
> library(vcfR)
> vcf_file <- "hoge.vcf"
> vcf <- read.vcfR(vcf_file, verbose = FALSE) # read vcf
> gt <- extract.gt(vcf) # extract genotypes
> gt <- as.data.frame(t(gt), stringsAsFactors = FALSE) # transpose and data.frame
> gt[1:6,1:6] #縦軸は個体ID、横軸は座位 #1/1はアリルの組み合わせ
186:220 262:219 569:90 637:222 845:35 1039:161
101_pair 0/0 0/0 0/0 0/0 0/0 0/0
102_pair 0/0 0/0 0/1 0/0 0/0 0/0
104_pair 0/0 0/0 1/1 0/0 0/0 0/0
105_pair 0/0 0/0 0/0 0/0 0/0 0/0
106_pair 0/0 0/0 0/0 0/0 0/0 0/0
107_pair 0/1 0/0 0/0 0/0 0/0 0/0
> Hoge <- vcfR2genind(vcf, sep = "[|/]")
#”Hoge”の箇所にgenindデータのファイル名として自分で使いたいものを入力します
> Hoge #中身を見てみると
/// GENIND OBJECT /////////
// 33 individuals; 202 loci; 404 alleles; size: 167.6 Kb
// Basic content
@tab: 33 x 404 matrix of allele counts
@loc.n.all: number of alleles per locus (range: 2-2)
@loc.fac: locus factor for the 404 columns of @tab
@all.names: list of allele names for each locus
@ploidy: ploidy of each individual (range: 2-2)
@type: codom
@call: adegenet::df2genind(X = t(x), sep = sep)
// Optional content
- empty -
DAPCスタート
> install.packages("adegenet")
> library(adegenet)
> Hoge@tab[1:6,1:6] #行列データの中身を一部表示（ゲノムデータは膨大なので一部だけ）
#メモ：「data名@属性」でデータの中身のデータを属性ごとに表示できる
#大事そうなので後で調べる
#以下のデータは本当は行列形式で出てきます
186:220.0 186:220.1 262:219.0 262:219.1 569:90.0
101_pair 2 0 2 0 2
102_pair 2 0 2 0 1
104_pair 2 0 2 0 0
105_pair 2 0 2 0 2
106_pair 2 0 2 0 2
107_pair 1 1 2 0 2
569:90.1
101_pair 0
102_pair 1
104_pair 2
105_pair 0
106_pair 0
107_pair 0
> Hoge@loc.n.all[1:6]
186:220 262:219 569:90 637:222 845:35 1039:161
2 2 2 2 2 2
> Hoge.group <- find.clusters(Hoge)
Choose the number PCs to retain (>= 1): #最大推定クラスター数（今回は最大個体数にしてみた）
33↓
#BIC値グラフが表示される #BIC値が小さいものが最適なので今回は3が最適値
Choose the number of clusters (>=2): #最適クラスター数を入力
3
> Hoge.group
$Kstat
K=1 K=2 K=3
121.35494 100.11358 50.47074
$stat
K=3
50.47074
$grp#以下のデータは本当は行列形式で出てきます
101_pair 102_pair 104_pair 105_pair 106_pair 107_pair
2 2 2 2 2 2
108_pair 109_pair 110_pair 111_pair 112_pair 113_pair
2 2 2 2 2 2
114_pair 115_pair 116_pair 117_pair 118_pair 119_pair
2 2 2 2 2 2
120_pair 122_pair 123_pair 124_pair 125_pair 126_pair
2 2 1 1 1 1
127_pair 128_pair 129_pair 130_pair 131_pair 132_pair
3 3 3 3 3 3
2512_pair 2817_pair 4413_pair
2 2 2
#K＝3で分けるのが最適っていうのと、各個体がどのグループに入るのかが表示される
> dapc.Hoge.groups <- dapc(Hoge,Hoge.group$grp) #DAPC関数にSNPデータと上記で行ったグループ分けの情報を代入
Choose the number PCs to retain (>=1):
33
Choose the number discriminant functions to retain (>=1):
1 #いくつの関数で説明できるか、表示された図を参考にして見てみる #今回は一つ
> dapc.Hoge.groups
#################################################
# Discriminant Analysis of Principal Components #
#################################################
class: dapc
$call: dapc.genind(x = Hoge, pop = Hoge.group$grp)
$n.pca: 30 first PCs of PCA used
$n.da: 1 discriminant functions saved
$var (proportion of conserved variance): 1
$eig (eigenvalues): 1.52e+33 2.258e+30 0.1183 vector length content
1 $eig 3 eigenvalues
2 $grp 33 prior group assignment
3 $prior 3 prior group probabilities
4 $assign 33 posterior group assignment
5 $pca.cent 404 centring vector of PCA
6 $pca.norm 404 scaling vector of PCA
7 $pca.eig 30 eigenvalues of PCA
data.frame nrow ncol
1 $tab 33 30
2 $means 3 30
3 $loadings 30 1
4 $ind.coord 33 1
5 $grp.coord 3 1
6 $posterior 33 3
7 $pca.loadings 404 30
8 $var.contr 404 1
content
1 retained PCs of PCA
2 group means
3 loadings of variables
4 coordinates of individuals (principal components)
5 coordinates of groups
6 posterior membership probabilities
7 PCA loadings of original variables
8 contribution of original variables
#結果の図を描く
> compoplot(dapc.Hoge.groups,main="Hoge groups",legend=TRUE, lab="",border=NA,
xlab="",las=2,col.lab="white",cex.axis=0.7,cex.main=1)
> Hoge.group$grp #個体データと、グループ分けの対応関係を見る
#以下のデータは本当は行列形式で出てくるのでもう少し見やすいです #ごめんなさい
101_pair 102_pair 104_pair 105_pair 106_pair 107_pair
2 2 2 2 2 2
108_pair 109_pair 110_pair 111_pair 112_pair 113_pair
2 2 2 2 2 2
114_pair 115_pair 116_pair 117_pair 118_pair 119_pair
2 2 2 2 2 2
120_pair 122_pair 123_pair 124_pair 125_pair 126_pair
2 2 1 1 1 1
127_pair 128_pair 129_pair 130_pair 131_pair 132_pair
3 3 3 3 3 3
2512_pair 2817_pair 4413_pair
2 2 2
Levels: 1 2 3
これでおしまい！
