smartPCA
pupulation geneticsではPCAをよく使う。外れサンプルの検出や集団構造化の補正など。使い方メモ。
まず、Plink形式のファイルからEIGENSTRAT形式に変換しないといけない。それぞれのファイルの対応は以下の通り。
hoge.bed = fuga.geno
hoge.fam = fuga.ind
hoge.bim = fuga.snp
これを変換するソフトconvertfというソフトが同梱されていて、以下のような変換ファイルを食べさせて変換させることができる。
bed2eigen.par
genotypename: hoge.bed
snpname: hoge.bim
indivname: hoge.fam
outputformat: EIGENSTRAT
genotypeoutname: fuga.geno
snpoutname: fuga.snp
indivoutname: fuga.ind
familynames: NO
convertf -p bed2eigen.par
こうするとfuga.geno/fuga.snp/fuga.indができる。次にPCAするときは同様に.parファイルを作る。
eigen2pca.par
genotypename: fuga.geno
snpname: fuga.snp
indivname: fuga.ind
evecoutname: piyo.evec
evaloutname: piyo.eval
smartpca -p eigen2pca.par
こうすると各サンプルの固有値の記載されたpiyo.evalと固有ベクトルの記載されたpiyo.evecが得られる。
さらにプロットするperlのスクリプトまでついており、
perl ploteig \
-i piyo.evec \
-c 1:2 \
-s ../out/myplot \
-x
で../out/myplot.pdf などが出力される。オプションはbed2eigen.parやeigen2pca.parファイルに記載すれば良い。fuga.indファイルに人種情報を入れておくと自動で色分けしてくれたりして楽しいが通常のfamファイルは人種情報を保持できないため後から自分でindファイルを作る必要がある。順番を変えるとまずいのでshellで一行ずつ処理するか、Rのmerge(~~~~, sort=F)とかを使うと良いと思います。
コメント
コメントを投稿