Mac
こちらをご覧下さい．Mac の場合は，こちらから最新版をダウンロードして，.dmg ファイルをダブルクリックするだけす．
　アプリケーションフォルダーに「R」と「R64」という二つのアプリケーションができます．これは使っている PC のプロセッサが i386 か x86_64 によるようです．Mac であればターミナルに，

uname -p

と入力すれば，プロセッサの情報が得られます．
注意:
・バージョンが古いために関数が正常に動かないことがよくあります．そのような場合は，最新のバージョンをインストールすると問題が解決されることがよくあります．
・MacOS Big Sur で R や Rscript が Terminal から起動できないです。こちらを参考に解決しました (2021 年 4 月)。

Linux

wget https://cran.r-project.org/src/base/R-3/R-3.6.3.tar.gz

tar xzvf R-3.6.3.tar.gz

cd R-3.6.3

sudo ./configure --with-readline=no --with-x=no

sudo make

sudo make install

こちらを参照しました。

Mac の場合，R に入っているエディタを使うとワンタッチで動作確認ができます．

x <- c(1,2,3,3,3,4,4,4,4,5,5,5,5,5,5,6,6,7,7,8,9)
hist(x)

ape は分子系統解析のパッケージです．通常は，R のコンソールに以下を打ち込めばインストールされます．

install.packages("ape")

詳しくはこちらを参照してください．

R の最新版に上記の手順で ape がインストールできないことがありました．そこで root で以下の手順に従いソースコードからインストールしました．

1. こちらから ape_3.3.tar.gz をダウンロード．

wget http://cran.r-project.org/src/contrib/ape_3.3.tar.gz

2. コンパイル

R CMD INSTALL ape_3.3.tar.gz

こちらの Linux の箇所を参照しました (2015 年 6 月)．

ミラーサイトの変更

CentOS7 で ape がインストールできない場合がありました。まず、サーバーで、sudo を使って R を立ち上げました：

sudo /usr/local/bin/R

その後、ape のインストールを試しました：

install.packages("ape")

すると、以下のエラーメッセージが R から出ました。

URL 'https://cran.ism.ac.jp/src/contrib/PACKAGES' を開けません

ミラーサイトを変更したら無事インストールできました。インストールの途中でミラーサイトを JAPAN にすると、統計数理研究所が選ばれます。この時は、 統数研のサーバーに何らかの問題があったようです。

R に以下を入力して、再度、install.packages("ape") を行います。

options(repos="https://ftp.yz.yamagata-u.ac.jp/pub/cran/")

こちらを参照しました (2021 年 3 月)。

ライブラリを指定してパッケージをインストール

web tool などの場合、Rscript などから R を操作するので、ユーザーが apache になります。この場合、パッケージのインストールする場合、場所を指定する必要があります。そうしないと、terminal からの操作では ape が使えますが、web tool (apache) からの操作では ape が使えないことがあります。 apache が R を立ち上げるので、該当するライブラリに ape が入っている必要があるのです。

まず、sudoコマンドを用いて、terminal から R を立ち上げます。

sudo R

その後、R のコンソールから以下の操作を行います。

library()

以下のような表示が出ます。ライブラリ ‘/usr/lib64/R/library’ に ape が入っていませんね。

パッケージ （ライブラリ ‘/home/jinoue/R/x86_64-redhat-linux-gnu-library/3.6’ 中):

ape Analyses of Phylogenetics and Evolution
Rcpp Seamless R and C++ Integration

パッケージ （ライブラリ ‘/usr/lib64/R/library’ 中):

base The R Base Package
boot Bootstrap Functions (Originally by Angelo Canty
for S)
....

以下のコマンドによって、ライブラリ ‘/usr/lib64/R/library’ に ape をインストールします。

install.packages("ape",lib="/usr/lib64/R/library")

インストールできたか確認してみます。

library()
パッケージ （ライブラリ ‘/usr/lib64/R/library’ 中):

ape Analyses of Phylogenetics and Evolution
base The R Base Package
boot Bootstrap Functions (Originally by Angelo Canty
for S)
....

こちらを参照しました。ありがとうございました (2020 年 5 月)。

R コンソールをクリアにする

option + command + l

ターミナルの「コマンド+K」にあたります．

スクリプトを途中で止める

q()

です．perl の「exit;」 に相当します．

作業ディレクトリの内容を見る

dir()

作業ディレクトリのアドレスを表示する

getwd()

作業ディレクトリを変更する

setwd('/Users/JunInoue/Downloads/histR')

ファイルを読み込んでデータフレームを作成する．
read.table('xxx.txt', header = T, sep = '\t')　とすると，データが欠けている場合　(例: エクセルで空欄) は NA と認識．詳しくは竹中先生のページ「補足: 結束値などの扱い」を参照してください．

read.table()

read.table でデータを df と入力して読み込んだ後に，因子になってしまった変数を文字列に変換する．こちらを参照しました．

df$ID<- as.character(df$ID)

データの変数名をリストアップする

names()

要素の種類をリストアップする

unique()

データ内にある変数の種類を確認する

str()

変数の確認

ls() または objects()

解析結果を .dat ファイルとして保存する

save(aValues, bValues, nValues, file="tancha.dat")

.dat ファイルを読み込む

load("tancha.dat")

グラフをPDF ファイルとして保存する

dev.copy2pdf(out.type="pdf",file="tancha.pdf")

データ内にある変数の種類を確認するパターンマッチ grep regexpr．こちらを参考にしました．ここでは files ベクトル内部で ".gb"で終わる要素を抽出しています．

grep("\\.gb$", files)

日本語の文字化け

plot した図で文字化けする場合は，こちらを参考に処置してください．

対数減数曲線
例題として，WGD 由来重複遺伝子の対数減数モデルの曲線を描きます．

Force et al. (1999) は Nei and Roychoudhury (1973) が使っている対数曲線の式に基づいて，重複遺伝子の機能が失われる確率を「1-e^(-2ut)」と提唱しています．t は世代時間，u は null mutation rate です．t を時間と考えると，ここで用いられた e^(-2ut) は重複遺伝子の欠失率をモデル化した式と見なすことができます．この式，

y= e^(-2ut)

を R を使ってグラフ化してみましょう．突然変異率　u は 1E-6 とします．
1E-6 は指数標記で 1*10^(-6) と同じです．

## グラフをプロットします．xlim = c(0, 5e+6) は表示される x 軸の範囲です．

lineExponential = function(x) exp(-2*(1E-6)*x)
plot (lineExponential, xlim = c(0, 5E+6), 0,col = "blue")

## WGD が生じた時点 x=0 の残存率は 1 となっていることを確認します．

x=0
exp(-2*(1E-6)*x)

## 100 万年後に残っている割合を見てみましょう．0.135... になっています．

x=1E+6
exp(-2*1E-6*x)

x <- seq(0, 3, 0.05)
plot(x,dlnorm(x, mean=0.068, sd=sqrt(0.24)), type="n")
curve(dlnorm(x, mean=0.068, sd=sqrt(0.24)), type="l",add=T)

以下のようなカーブが描けるはずです．
こちらを参考にしました．Log-noarmal distribution についてはこちらを参考にしてください．

こちらをご覧下さい．

freq <- c(70,30)
names(freq) <- c("好き","嫌い")
barplot(freq, col="blue", xlab="好き嫌い", ylab="度数")
abline(h=0)

をテキストファイルとして保存します．メニューバーのファイル/ソースを読み込むを選んで，上記のテキストファイルを読み込みます．

a <- c(1, 3, 5)
b <- c(2, 4, 6)
#x <- matrix(c(1, 2, 3, 4, 5, 6), 2, 3)
x <- rbind(a, b)
x
barplot(x, horiz=T, beside=T, names=1:3) # Horizontal
#barplot(x, horiz=T, names=1:3) # Horizontal bar

例題 1
ヒストグラムを作成します．?hist と入力して，最後の例題を見てみましょう．R にはすでに例題が入っていて，

hist(islands)

とすることで，ヒストグラムを作成してくれます．例えば，x というベクトルを作成して，

x <- c(1,1,2,3,4,4,5)

hist(x)

としても，ヒストグラムが作成されます．こちらに詳しい説明があります．

最頻値
hist 関数を使いますが，plot させません．

xHist <- hist(x, plot=FALSE)
which(xHist$counts==max(xHist$counts))

こちらを参照しました (2013 年 4 月)．

例題 2
．histR.tar.gz

　ダウンロードして得られた histSimple.R　を R アイコンに　ドラッグ&ドロップして R のエディターで開きます．コマンド+ A で全部を選択した後，コマンド+リターンで R を走らせます．

例題 3
1000a.dat.tar.gz
.dat として保存されたデータを読み込んでヒストグラムを書きます (2015 年 6 月)．

> cc <- load("~/.../1000a.dat")
> str(cc)
chr [1:4] "aValues" "bValues" "dValues" "nValues"
> aValues
[1] 0.0025884912 0.0023656549 ....
[1000] 0.0026137924
> median(log(aValues))
[1] -5.964303
> hist(log(aValues),breaks=50)

>Z <- stats::rnorm(10000)
>table(cut(Z, breaks = -6:6))
>as.integer(table(cut(Z, breaks = -6:6)))

こちらを参照しています．意外とこのテクニックは見つかりません．サイトがなくなったときのために，コマンドだけ移させていただいています．

Lattice (ラティス) パッケージを使うと，複数のグラフを一枚の図として並べるなど，R を使ったグラフィックスがさらに便利になります．

install.packages("lattice")

によって lattice パッケージをインストールしてください．その後，ダウンロードして得られた histLattice.R を走らせると，以下のような図が得られるはずです．
histLatticeR.tar.gz

d <- read.table(file="infile.txt",header = TRUE)

### Retrieve the recoreds with distance for less than 2
d_less2 <- d[d$Dist < 2,]

################## Make histograms
pdf("histLattice.pdf")
library(lattice)
histogram( ~ Dist | Class, data = d_less2,
layout = c(1,3),nint = 30, xlab = "Distance", strip = FALSE,
strip.left = TRUE, ylab="Frequencies")
dev.off()

################## Calcurate the averages of distances
## Check the variables
names(d_less2)
unique(d_less2$Class)

## Retrieve the lines of each class
typeR <- d_less2[d_less2$Class == "typeR",]
typeM <- d_less2[d_less2$Class == "typeM",]
typeO <- d_less2[d_less2$Class == "typeO",]

## Calculate the average of each class
mean(typeR$Dist)
mean(typeM$Dist)
mean(typeO$Dist)

いくつか存在する変数の相互関係を考察するには，総当たりで散布図を表示する対散布図が便利です．以下の行をコピーして R のコンソールに貼付けると，下に示した図が得られます．

pairs (iris[1:4], main = "Anderson's Iris Data -- 3 species", pch = 21, bg = c("red", "green3", "blue")[unclass(iris$Species)], panel=panel.smooth)

iris[1:4]
iris というデータはすでにデータフレームとして R に組み込まれています．

bg = c("red", "green3", "blue")[unclass(iris$Species)],
種ごとに色分けをしています．Species というコラムには 3 種があります．「iris」 と打ち込んで，データフレームを確認してください．

panel=panel.smooth
散布図の点の傾向を示す曲線が描かれます．

詳しくは，コンソールに「?pairs」と打ち込むと得られる説明を読んでください．

nls, glm 関数
ある分布のプロットを，ロジスティック関数 (logistic model) ，指数関数 (exponential model)，左右非対称の S 字曲線の関数 (sigmoidal model) へ回帰します．長いのでこちらのページに移動しました [2013 年 1 月]．

library(nlmrt)

time = c(-353.0, -276, -191, -159, -154, -70, 0, 0, 0, 0, 0, 0)
gene_num = c( 7019, 1040, 868, 839, 786, 610, 630, 505, 563, 577, 476, 551)
threeRxTime <- time[[1]]

plot(time, gene_num, xaxp = c(-450, 0, 9),
xlim=c(-450,0), pch=20, ylim=c(0,9000))
par(new=T)

##################################################################
# Double exponential:
#F(t)=a *e^(-b *(t+353))+c *e^(-d*(t+353))
#F(t)=5739*e^(-0.448*(t+353))+1280*e^(-0.00237*(t+353))
# <- Rough estimation by manual

dat <- data.frame(n=gene_num, t=time)
doubleExponential <- nlxb( n~a*exp(-b*(t-threeRxTime))+
c*exp(-d*(t-threeRxTime))+e,
data=dat, start=c(a=6000,b=0.01,c=1000,d=0.002,e=100), lower=0,
upper=c(a=10000,b=0.1,c=10000,d=0.1,e=10000))

va <- doubleExponential$coefficients[["a"]]
vb <- doubleExponential$coefficients[["b"]]
vc <- doubleExponential$coefficients[["c"]]
vd <- doubleExponential$coefficients[["d"]]
ve <- doubleExponential$coefficients[["e"]]

linearDoubleExp <- function(x) va*exp(-vb*(x-time[[1]]))
+ vc*exp(-vd*(x-time[[1]]))+ve
plot (linearDoubleExp, -450, 0, add=TRUE, col ="blue")

doubleExponentialNls <- nls(gene_num ~ va*exp(-vb*(time-time[[1]]))
+ vc*exp(-d*(time-time[[1]]))+ve, start=list(d=0.001))
AIC(doubleExponentialNls)

参考文献
回帰分析入門
とてもわかりやすいです [2013 年 12 月]．

最適化問題
optim()/optimize() 関数について [2014 年 10 月]．

関数の値を最小にするパラメータの推定

バッチモードです．ターミナルから R を動かすと，perl などと組み合わせて使うことができます．以下をダウンロードして得られたフォルダに入って，

R --vanilla --quiet --args in.tre < script.R > log

と入力してください．Tree が書かれた 2 つのファイル (in.tre.png, in.tre.nwk) が得られるはずです．script.R がバッチファイルです．log はスクリーンアウトです．
batchR.tar.gz
こちらを参照しました．

# Load package ape.
library(ape)

# Read files
args <- commandArgs()
treeFile <- args[5] # in.tre

# Read tree
tr <- read.tree(treeFile)

# Open the device
png.file <- paste(treeFile, ".png", sep = "")

# Plot tree
png(png.file)
plot.phylo(tr)

# Save the tree as the newick format
nwk.file <- paste(treeFile, ".nwk", sep = "")
write.tree(tr, file=nwk.file)

dev.off()

スミルノフ・グラブス検定を用いて，root to tip length (ここでは Node10 から末端まで) が平均値から極端に異なる枝 (以下の図で * ) を検出するプログラムです．外群 (SeqSquirt) は検定から外されます．もちろん可能であれば相対速度検定を行うべきですが，計算が思ったほど簡単ではないので，ここではスミルノフ・グラブス検定を用いています．しかし，この検定方法が長枝の検出に向いてているかどうかは，検討が必要です．ご存知の方がいたら教えてください [2012 年 2 月]．
SGforBL.tar.gz

両側検定なので、p 値を計算する際の自由度は n-2 にしているようです。こちらを参照しました [2020 年 1 月]。

reroot すると，BS 値の位置がおかしくなると聞きました．以下は例題です．左の tree を d という leaf で rooting しています．しかし，私には問題ないように思えます．
rerootBS_fol.tar.gz

library(ape)
layout(matrix(1:2, 1, 2))
t=read.tree("test.tre")
plot(t)
nodelabels(t$node.label)
rerooted <- root(t,"d",r=T)
plot(rerooted)
nodelabels(rerooted$node.label)

RAxML の BS 付き tree は内群根幹 (外群直下) の分岐に BS 値が付けられていません．Notung の Rearrangement 解析などでは，この分岐を解消してしまうことがあるので，以下のスクリプトで 100% の BS 値を付けると便利です (2018 年 8 月)．
rootBS100_fol.tar.gz

$ ls
RAxMLtree.tre rootBS100.R
junINOUEpro@inouejun-no-MacBook-Pro|~/Downloads/rootBS100_fol
$ R
R version 3.3.2 (2016-10-31) -- "Sincere Pumpkin Patch"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin13.4.0 (64-bit)
．．．．．
> source("rootBS100.R")

This script draws a tree into a png file and gives color to specific leaves and labels.
treePlot.tar.gz
Move into the working directory and execute the treePlot.R file.

setwd("/Users/inouejun/Desktop/treePlot")
source("treePlot.R")

Get into the uncompress nhxReader.tar.gz file and type

perl nhxReader.pl > log.txt

rscript will read a R script, nhxReader.R. To convert nhx to newick format, I am using perl script. Colon within the node label [&&NHX:S=Percomorpha:D=N:B=100.0] should not be used for ape, read.tree function.

Note that the rscript should be the same version of your R. Otherwise, rscript cannot recognize the input file as follows:

[junINOUEpro:newickTreePlot_fol]$ rscript plotNewick.R MyTree.tre
Rscript execution error: No such file or directory

(Dec. 2017)

ape ではアミノ酸配列から NJ tree などを構築できないようなので，phangorn でやってみました．スクリプトを走らせる前に，「install.packages("phangorn")」によって，パッケージをインストールするのを忘れないでください．ただし，データによってはうまく動かないです．Ape を使って cDNA 配列 (1st, 2nd) で解析した方が良いように思えます (2017 年 1 月)．
njBS_Phangorn.tar.gz

r8s の解析を APE で行います．こちらをご覧ください．ただし，こちらによると chronogram は使われていないため，chronopl を使う必要があるみたいです．しかし以下のステップでは，枝の長さが違う tree を ultrametric tree にするだけで，まだ制約の付け方や時間軸の着いた tree 表示のやり方はわかりません．

##################

library (ape)
# ない場合は，install.packages("ape")

# get tree
data("landplants.newick") # example tree in NH format
landplants.newick
[1] "(Marchantia:0.033817,(Lycopodium:0.040281,......:0.012469):0.019372);"
tree.landplants <- read.tree(text = landplants.newick)

# plot tree
tree.landplants
plot(tree.landplants, label.offset = 0.001)

# estimate chronogram１　[lambda を 1 にしました]
chrono.plants <- chronopl(tree.landplants,1)

# plot
plot(chrono.plants, label.offset = 0.001)

# value of NPRS function for our estimated chronogram [これは動きませんでした]
NPRS.criterion(tree.landplants, chrono.plants)

http://cran.r-project.org/web/packages/TreeSim/TreeSim.pdf
にある例題を試しました．

> numbsim<-3
> age<-10
> lambda<-0.3
> mu<-0
> K<-40
> tree<- sim.bd.age(age,numbsim,lambda,mu,mrca=FALSE,complete=FALSE,K=K)
エラー: 関数 "sim.bd.age" を見つけることができませんでした
> LTT.plot(c(list(tree),list(c(age,age,age))))
エラー: 関数 "LTT.plot" を見つけることができませんでした

上のようになるはずです．以下を実行しました．

> install.packages("ape")
> install.packages("geiger")
> install.packages("TreeSim")

> library("ape")
> library("TreeSim")

するとうまく動きました．

> numbsim<-3
> age<-10
> lambda<-0.3
> mu<-0
> K<-40
> tree<- sim.bd.age(age,numbsim,lambda,mu,mrca=FALSE,complete=FALSE,K=K)
> LTT.plot(c(list(tree),list(c(age,age,age))))
>

PDF ファイルが出ました

library("ape")
library("TreeSim")

iterations <- 10

x <- 3
while (x <= iterations) {
x <- x + 1

numbsim <-3
age <-10
lambda <-0.1 * x

mu<-0
K<-40

tree<- sim.bd.age(age,numbsim,lambda,mu,mrca=FALSE,complete=FALSE,K=K)
par(new=T)
LTT.plot(c(list(tree), list(c(age,age,age)), add=TRUE))
}

expm1()
以下の二つは数学的に同じ数値を返すはずです．しかし，値が極端に少ないときは，浮動小数点に関するコンピュータの計算能力のエラーからおかしな答えが得られます．こちらによると，expm1(x) の方が exp(x)-1 に比べて正確な値を返すようです．

> expm1(1e-14)
[1] 1e-14
> exp(1e-14)-1
[1] 9.992007e-15

大きな値で計算すると，同じ答えが得られます．

> expm1(0.1)
[1] 0.1051709
> exp(0.1)-1
[1] 0.1051709

R prompt で以下のコマンドを入力してください．

source("http://bioconductor.org/biocLite.R")

有用なモジュールの追加．

biocLite("edgeR")
biocLite("DESeq")
biocLite("etc")
biocLite("Biobase")
install.packages("gplots")
install.packages("ape")

こちらをご覧下さい (2014 年 11 月)．

exp(x) = 1 + x + x^2/2 + x^3/6 + ....
exp(x) -1 = x + x^2/2 + x^3/6 + ....