SNPデータ入りのvcfファイルからPCAをする方法を教わったのでメモ。
デモデータは202座位33個体、参照URLの閲覧日は本記事を作成した2020年12月19日現在のものとします。
PCA とは
→A Tutorial on Principal Component Analysis pca.dvi (cmu.edu)
まずvcfファイルをインプットに用いるgenind形式に変換
> 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、横軸は座位
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
> install.packages("adegenet")
> library(adegenet)
> 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と一緒なので、DAPCとPCAは一緒にやってしまうと楽です
PCA解析を行う
> Hoge.tab <- tab(Hoge, freq=T, NA.method="mean")
#freqは論理値。もしTRUEならヒストグラムは結果のcounts成分である度数の表示。もしFALSEならdensity成分である確率密度がプロットされ，総面積は1になる。breaksが等間隔である(そしてprobabilityが指定されない)ときだけ既定値はTRUE。
#参照：R基本統計関数マニュアル （https://cran.r-project.org/doc/contrib/manuals-jp/Mase-Rstatman.pdf）
#何を言っているのかよくわかりません
#NA.method はある遺伝子座位において欠損値があった（あるSNPが欠失している個体がいた）場合に置換する関数
#NA.methodの置換方法は以下→ asis: NAをそのままにしておく; mean: 平均対立遺伝子頻度で置換する; zero: ゼロで置換する
#参照URL（https://rdrr.io/cran/adegenet/man/tab.html）
>pca.Hoge <- dudi.pca(Hoge.tab, scale=F,scannf=F)
> pca.Hoge #結果の書き出し #結果の中身の説明は以下のURLにて#https://www.rdocumentation.org/packages/ade4/versions/1.7-16/topics/dudi.pca
Duality diagramm
class: pca dudi
$call: dudi.pca(df = Hoge.tab, scale = F, scannf = F)
$nf: 2 axis-components saved
$rank: 30
eigen values: 19.53 12.71 0.7152 0.4582 0.2895 ...
vector length mode content
1 $cw 404 numeric column weights
2 $lw 33 numeric row weights
3 $eig 30 numeric eigen values
data.frame nrow ncol content
1 $tab 33 404 modified array
2 $li 33 2 row coordinates #これが主成分
3 $l1 33 2 row normed scores
4 $co 404 2 column coordinates
5 $c1 404 2 column normed scores
other elements: cent norm
> plot(pca.Hoge$li[,1],pca.Hoge$li[,2]) #第一主成分軸と第二種成分軸を図にプロット
>colorplot(pca.Hoge$li, pca.Hoge$li, transp=TRUE, cex=2, xlab="PC 1", ylab="PC 2")
#遺伝的な構成に従ってカラープロット
#$li: the principal components of the analysis; these are the synthetic variables summarizing the genetic diversity, usually visualized using scatterplots.
#参考: tutorial-basics.pdf (r-project.org)
#シンボルサイズの変更：cex=(size)（例：cex=2, ならシンボルサイズはデフォの二倍）
#軸ラベルのフォントサイズ：cex.lab = (size)
#軸の数字等（ラベル）の大きさを設定する：cex.axis = 1.8
主成分軸がどの程度変数を説明できているかの度合いを見る方法
> pca.Hoge$eig
[1] 1.953455e+01 1.270779e+01 7.152222e-01 4.581928e-01
[5] 2.895304e-01 2.853223e-01 2.421248e-01 1.962387e-01
[9] 1.928476e-01 1.587351e-01 1.426831e-01 9.987635e-02
[13] 9.782591e-02 7.733217e-02 7.604947e-02 5.470648e-02
[17] 5.238435e-02 4.004565e-02 3.448543e-02 2.720587e-02
[21] 2.296187e-02 1.752990e-02 1.408798e-02 1.146402e-02
[25] 1.040276e-02 5.194359e-03 2.667595e-03 7.002338e-04
[29] 5.573597e-04 3.094007e-05
#eig→Eigenvalue（固有値）：データの全情報量（データの全てのばらつき）のうち，それぞれの主成分がどのくらいの分量の情報を表現しているかを示す
#参照URL（https://www.ai.u-hyogo.ac.jp/~arima/lectures/JT-15.pdf）
> pca.Hoge$eig/sum(pca.Hoge$eig)
#固有値全量に対してその固有値が何割を説明しているかを書き出し
#この値は浮動小数点表記法 （floating point expression）で書かれている
#故に第一主成分は5.492055e-01=5.49..×(10の-1乗)=0.549=54.9% 程度を説明していることになる
[1] 5.492055e-01 3.572741e-01 2.010816e-02 1.288189e-02
[5] 8.140022e-03 8.021711e-03 6.807233e-03 5.517167e-03
[9] 5.421825e-03 4.462769e-03 4.011474e-03 2.807980e-03
[13] 2.750333e-03 2.174160e-03 2.138098e-03 1.538049e-03
[17] 1.472763e-03 1.125866e-03 9.695427e-04 7.648812e-04
[21] 6.455631e-04 4.928456e-04 3.960774e-04 3.223059e-04
[25] 2.924691e-04 1.460371e-04 7.499827e-05 1.968677e-05
[29] 1.566992e-05 8.698666e-07
これでPCAはおしまい
【追記1】
ずっと勘違いしていたのですが、
>Hoge.group <- find.clusters(Hoge)
Choose the number PCs to retain (>= 1):
最大予想クラスター数を入力
自分は解析個体数を入れました
Choose the number of clusters (>=2):
DAPC と同様、表示された図を見てBIC が一番大きくなるクラスター数を選択
の工程はDAPC ですね。申し訳ありませんでした。
【追記2】
どれが何のサンプルなのかわかるように表示したいなと思って方法を探してみました
データ解析・マイニングとR言語 (doshisha.ac.jp)
ここに載っている方法を使うと、プロットを点ではなくてサンプル番号で表示できます。
>plot(pca.Hoge$li[,1],pca.Hoge$li[,2],type="n")
>text(pca.Hoge$li[,1],pca.Hoge$li[,2])
