Photo by
daikuma
R:fasta形式のアミノ酸配列のアラインメントとそれに伴う系統樹の作成
library(msa)
library(seqinr)
#fastaファイルの読み込み
fasta <- readAAStringSet("〇〇.fasta",format="fasta")
#アラインメントの実行。"ClustalOmega" "Muscle"も使える
alignment <- msa(fasta,"ClustalW")
#クラスの変換。系統樹作成のために距離行列として扱えるようにしている。
alignment2 <- msaConvert(alignment, type="seqinr::alignment")
#距離行列の計算
d <- dist.alignment(alignment2, "identity")
pdf("系統樹_ward.D2-1")
hc <- hclust(d, "ward.D2")
plot(hc)
dev.off()また、途中で行った距離行列の計算は平方根で換算しているので、それが嫌な人は以下のようなコードを使うと良いかも
library(msa)
library(seqinr)
library(ape)
#fastaファイルの読み込み。ファイル名は適宜変えてもらって。
fasta <- readAAStringSet("〇〇.fasta",format="fasta")
#アラインメントの実行。"ClustalOmega" "Muscle"も使える
alignment <- msa(fasta,"ClustalW")
#クラスの変換。系統樹作成のために距離行列として扱えるようにしている。
alignment3 <- msaConvert(alignment, type="ape::AAbin")
#距離行列の計算
d2<- ape::dist.gene(alignment3, "pairwise")
pdf("系統樹_ward.D2-2")
hc <- hclust(d2, "ward.D2")
plot(hc)
dev.off()アラインメント結果を表示するのに、途中の「alignment」を入力すれば見れるが、全部見たい場合は以下のように打つと良い。二行目はpdfで結果を見るものになるが、これは普通にアラインメントとして見栄えもいいし優秀。
print(alignment,show="complete")
msaPrettyPrint(alignment, output="pdf",showLogo="none", askForOverwrite=FALSE, verbose=FALSE)masPrettyPrintは場合によっては"pdflatex not found"のようなエラーが吐かれることがあります。その場合はhttp://miktex.org/2.9/setupからMiktexというものをインストールし、以下のコマンドを打ち込むと多分使えるようになります。詳しい説明は参考文献②を参照してください。
Sys.getenv("PATH")
Sys.setenv(PATH=paste(Sys.getenv("PATH"),"C:/Users/MY_USER_NAME/AppData/Local/Programs/MiKTeX/miktex/bin/x64",sep=";"))これでFASTAファイルさえあれば一瞬でアラインメントから系統樹作成まで持っていけるはずです。
参考文献⇩
https://bioconductor.org/packages/release/bioc/vignettes/msa/inst/doc/msa.pdf