Level1_3_R.markdown

次世代シークエンサーDRY解析教本（細胞工学別冊）の「Level1 [3] Rの使い方」で利用したソースコード

はじめに

GitHubのページ（https://github.com/kokitsuyuzaki/Level1_3_R）に書かれた内容が、plotlyの仕様変更により動かなくなったため、改訂版を作成した（2016/2/23）

DRY解析本の他の箇所の動かない箇所は、学研メディカル秀潤社のページで正誤表が出ているため、確認されたし

DRY解析本の執筆段階では、ploltyはブラウザ越しで使うWebアプリケーションで、手元のデータをplotlyで可視化する場合、向こうのサーバに毎回転送する必要があった

そのため、plotlyのアカウント作成、APIの設定などめんどくさいステップが多かったと思う（https://github.com/kokitsuyuzaki/Level1_3_Rのplolty.Rに相当）

ploltyは現在、javascriptライブラリ「plotly.js」としてオープンソースソフトウェアとして利用できるようになったため、めんどくさい作業は必要なくなった

そのため、以下はplotly.jsで作り直した

以下のコードは何も考えずコピー&ペーストするだけで動くはずですが、もし、動かなくなった場合、koki.tsuyuzaki [at] gmail.comまでご一報ください

演習① : 準備

演習①では、R/R-Studioのインストール、拡張パッケージのインストールの仕方を説明した

# R-Studioのデモ
# 一様乱数を10個生成
d <- runif(10)
# プロット
plot(d)

getwd() # 作業ディレクトリ
setwd("~/Desktop") # ディレクトリ変更

# CRANパッケージのインストールの仕方（Update all/some/none?と聞かれたらaと入力）
install.packages("kernlab", type="source")

# Bioconductorパッケージのインストールの仕方
source("http://bioconductor.org/biocLite.R")
biocLite("meshr", ask=FALSE)

# GitHubパッケージのインストールの仕方
install.packages("devtools")
library("devtools")
devtools::install_github("rikenbit/CCIPCA", force=TRUE)

演習② : Rの基本操作

演習②では、Rの基本的な操作方法について説明した

# 四則演算（数字）
1 + 1 # 足し算
3 - 1 # 引き算
6 * 2 # 掛け算
4 / 2 # 割り算
2^5 # 2の5乗
log10(5) # 5の常用対数値

# オブジェクトに格納
A <- 1 + 1 # オブジェクトに格納
A # オブジェクト名を打つと、値が画面に出力される（print(A)も同様）
A * 2 # (1 + 1) × 2を実行した事と同じ

# 値の型
try(1 + "1") # 数値と文字列を足すことはできない
paste0("ABCDE", "1") # 文字列同士を連結することはできる
paste0("ABCDE", 1) # paste0関数の仕様で、文字列と数値を連結することはできる

# 真偽
1 == 1 # これはTRUE（==は同じという意味）
1 == 2 # これはFALSE
"A" == "A" # これはTRUE
"A" == "B" # これはFALSE
"A" != "B" # これはTRUE（!=は違うという意味）
# TRUE/FALSEはif文と一緒に使うことが多い
X <- 1
Y <- 2
if(X == Y){
  print("同じ")
}else{
  print("違うよ！")
}

####################### ベクトル ########################
# c()でベクトルを作成
A <- c(2, 3, 1)
B <- c(-2, 1, -2)
# 最大値、最小値、平均値、中央値
max(A)
min(A)
mean(B)
median(B)
# 要素ごとの掛け算（他の言語と違うところ）
A * B
# ベクトル同士の内積
A %*% B
# 文字でもベクトル
C <- c("A", "B", "C", "A", "AA")
# Aがどこの場所にあるか
which(C == "A")
# lengthと組み合わせて使うことが多い
length(which(C == "A"))

######################### 因子 ##########################
# 同じものが複数
D <- factor(c("A", "A", "B", "B", "B"))
nlevels(D) # 2つある
levels(D) # A, Bがある

######################### 行列 ##########################
# matrix()で行列を作成
E <- matrix(runif(6), nrow=2, ncol=3)
E[1,]
E[,2]
# 行数、列数
nrow(E)
ncol(E)
dim(E)
# cbind、rbind
cbind(A, B)
rbind(E, A, B)
# 行レベルでの和、列レベルでの和・平均値
rowSums(E)
colSums(E)
rowMeans(E)
colMeans(E)

##################### データフレーム #####################
# 行列は一つの型しか扱えない
# 複数の型（数値、文字、TRUE/FALSE、因子など）を扱いたいときはデータフレーム
data(iris)
head(iris)
iris$Petal.Length # $で各列にアクセスできる

######################## リスト #########################
# もはや行列で表現できないようなデータもリストで表現できる
G <- list(X = c(1,2,3), Y = matrix(runif(9), nrow=3), Z = TRUE)
G$X # $で各要素にアクセスできる
G$Y
G$Z

############### Rオブジェクトに対してのあれこれ #############
is(iris)
str(iris)
str(G)
x <- runif(10)
y <- runif(10)
lr <- lm(y ~ x)
str(lr)
lr$coefficients

################### データ保存・読み込み ##################
# 保存
# 基本的にTSVかCSVで保存しよう（セルを後で認識できるから）
write.table(iris, "iris.tsv", sep="\t")
write.table(iris, "iris.csv", sep=",")
# 画面に出力される文字をそのままファイルに保存しておきたいとき
sink("iris.txt")
iris
sink()
# Rのオブジェクトそのものを保存しておきたいとき
save(iris, file="iris.RData")

# 読み込み
iris1 <- read.table("iris.tsv", sep="\t")
iris2 <- read.table("iris.csv", sep=",")
# 今まで宣言してきたRオブジェクトを見てみる
ls()
# 一度消してみる
rm(iris)
ls()
load("iris.RData")
# 消したirisが復活している
ls()

# read.table時に文字がfactor化されてしまう（よくはまる箇所）
write.table(data.frame(moji = c("A", "B", "C", "D", "E"), suuji = runif(5)), "factor_test.tsv", sep="\t")
H <- read.table("factor_test.tsv")
H[,1]
# strigsAsFactorsをFALSEとしておくことでfactor化を防ぐことができる
H <- read.table("factor_test.tsv", stringsAsFactors = FALSE)
H[,1]

演習③ : データの前処理

演習③では、データ解析の前に不可欠な作業である「前処理」について説明した

# データをとってくる
download.file("http://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE54006&format=file&file=GSE54006%5Fumitab%2Etxt%2Egz", destfile="Count.txt.gz") # 発現量データ
download.file("http://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE54006&format=file&file=GSE54006%5Fexperimental%5Fdesign%2Etxt%2Egz", destfile="ExpDesign.txt.gz") # 実験情報データ

# ファイルを解凍する
try(system("gunzip -f Count.txt.gz"))
try(system("gunzip -f ExpDesign.txt.gz"))

# データを読み込む
Count <- read.table("Count.txt", row.names=1, header=TRUE, sep="\t")
ExpDesign <- read.table("ExpDesign.txt", row.names=1, header=TRUE, sep="\t", skip=6, stringsAsFactors = FALSE)

# どのくらいのサイズの行列か見てみる
dim(Count) # 20190 × 4590
dim(ExpDesign) # 4598 × 14

# 上の行だけ見て、内容を推測する
head(Count)
head(ExpDesign)


################### Count ###################
which(colSums(Count) == 0) # 90個、発現量0のウェルがある
ncol(Count) # 4590ウェル
Count <- Count[, which(colSums(Count) != 0)] # ゼロ細胞を取り除く
ncol(Count) # 4500ウェルまで減る

# IDが書かれた列を追加する
Count <- t(Count)
# 列名をWellにする
Count <- data.frame(Count, Well=sub("X", "", rownames(Count)))



################# ExpDesign #################
# 列名をWellにする
colnames(ExpDesign)[14] <- "Well"

levels(ExpDesign$group_name) # Cell typeに関する情報が、group_name列に記述されている
which(ExpDesign$group_name == "B cell") # B cell (CD19+/B220+) : 48細胞
which(ExpDesign$group_name == "CD8+pDC") # pCD (CD11c+/CD8a+/PDCA-1+) : 96細胞
which(ExpDesign$group_name == "monocyte_or_neutrophil") # Monocyte (Gr-1+/CD11b+) : 48細胞
which(ExpDesign$group_name == "NK_cell") # NK (NK-1.1+) : 48細胞

# 4種類の細胞のデータだけを取り出す
ExpDesign <- rbind(
  ExpDesign[which(ExpDesign$group_name == "B cell"), ],
  ExpDesign[which(ExpDesign$group_name == "CD8+pDC"), ],
  ExpDesign[which(ExpDesign$group_name == "monocyte_or_neutrophil"), ],
  ExpDesign[which(ExpDesign$group_name == "NK_cell"), ]
  )
# 48 + 96 + 48 + 48 = 240になってる
nrow(ExpDesign)
# 1細胞のみ入っているものだけ取り出す
sort(unique(ExpDesign$number_of_cells))
ExpDesign <- ExpDesign[which(ExpDesign$number_of_cells == "1"), ]
# 1細胞だけである事を確認
sort(unique(ExpDesign$number_of_cells))
# 237に減る
nrow(ExpDesign)
# 必要なデータだけ、取り出す
ExpDesign <- ExpDesign[, c(11,14)]


############## ExpDesignとCount #############
# ExpDesignとCountをマージする
Count <- merge(ExpDesign, Count, by="Well")
# 保存する
save(Count, file="Count.Rdata")

演習④ : データ解析

演習④では、データ解析（主成分分析）について説明した

# データ読み込み
load("Count.Rdata")

# 主成分分析
pca_data <- t(Count[, 3:ncol(Count)])
pca <- prcomp(pca_data, scale=TRUE) # 正規化（平均値0, 分散1にする）

# データ保存
save(pca, file="pca.Rdata")

演習⑤ : データの可視化

演習⑤では、データをプロットして、細胞がどのように分布しているか可視化した

まず、rglパッケージのplot3dを使った

# パッケージインストール（一度やればよい）
install.packages("rgl")

# パッケージ読み込み
library("rgl")

# データ読み込み
load("Count.Rdata")
load("pca.Rdata")

# 主成分寄与率
png(file="Contribution_PCA.png")
plot(pca$sdev/sum(pca$sdev), xlab="Principal Components", ylab="Proportion  of Variance", main="Contribution rate of the Principal Component")
dev.off()

# バイプロット
png(file="Biplot_PCA.png")
biplot(pca)
dev.off()

# 細胞をソートする際の表面マーカーで色を付ける
Label <- data.frame(group_name=Count$group_name, Well=Count$Well, Color=NA)
# B cell => オレンジ
Label$Color[which(Label$group_name == "B cell")] <- rgb(1, 0.25, 0, 0.5)
# NK T cell => 紫
Label$Color[which(Label$group_name == "NK_cell")] <- rgb(0.75, 0, 1, 0.5)
# pDC => 青
Label$Color[which(Label$group_name == "CD8+pDC")] <- rgb(0, 0, 1, 0.5)
# Monocyte => 黒
Label$Color[which(Label$group_name == "monocyte_or_neutrophil")] <- rgb(0, 0, 0, 0.5)

# 主成分分析
plot3d(pca$rotation[, 1:3], xlab="PC1", ylab="PC2", zlab="PC3", col=Label$Color)
writeWebGL(dir="PCA_3d", width=500)

外れ値細胞があったため、次にplotlyパッケージを使う（修正した箇所（2016.10.5））

#️ パッケージインストール（一度やればよい）
install.packages("devtools")
library("devtools")
devtools::install_github("ropensci/plotly", force=TRUE)

library(plotly)

# データ読み込み
load("pca.Rdata")

pca_plotly <- data.frame(
Celltype=as.factor(Label[,1]), PC1=pca$rotation[,1], PC2=pca$rotation[,2], PC3=pca$rotation[,3])

pca_plotly %>%
	plot_ly(x=~PC1, y=~PC2, z=~PC3,
    type = "scatter3d", mode = "markers",
    text = as.vector(Label[,2]),
    color =~Celltype
    )

外れ値細胞は、Well IDが14_1875のものだとわかる

14_1875を除去して再度PCAを実行し、Count2.Rdata、pca2.Rdataとして保存する

# データ読み込み
load("Count.Rdata")

# 主成分分析（14_1875を除く）
Count2 <- Count[setdiff(1:nrow(Count), which(Count$Well == "14_1875")),]
pca_data2 <- t(Count2[, 3:ncol(Count2)])
pca2 <- prcomp(pca_data2, scale=TRUE)

# データ保存
save(Count2, file="Count2.Rdata")
save(pca2, file="pca2.Rdata")

14_1875を除去した結果をrglパッケージのplot3dで可視化する

# パッケージ読み込み
library("rgl")

# データ読み込み
load("Count2.Rdata")
load("pca2.Rdata")

# 主成分寄与率（14_1875除去後）
png(file="Contribution_PCA2.png")
plot(pca2$sdev/sum(pca2$sdev), xlab="Principal Components", ylab="Proportion  of Variance", main="Contribution rate of the Principal Component")
dev.off()

# バイプロット
png(file="Biplot_PCA2.png")
biplot(pca2)
dev.off()

# 細胞をソートする際の表面マーカーで色を付ける
Label2 <- data.frame(group_name=Count2$group_name, Well=Count2$Well, Color=NA)

# B cell => オレンジ
Label2$Color[which(Label2$group_name == "B cell")] <- rgb(1, 0.25, 0, 0.5)
# NK T cell => 紫
Label2$Color[which(Label2$group_name == "NK_cell")] <- rgb(0.75, 0, 1, 0.5)
# pDC => 青
Label2$Color[which(Label2$group_name == "CD8+pDC")] <- rgb(0, 0, 1, 0.5)
# Monocyte => 黒
Label2$Color[which(Label2$group_name == "monocyte_or_neutrophil")] <- rgb(0, 0, 0, 0.5)

# 主成分分析
plot3d(pca2$rotation[, 1:3], xlab="PC1", ylab="PC2", zlab="PC3", col=Label2$Color)
writeWebGL(dir="PCA_3d_2", width=500)

14_1875を除去した結果をplotlyでも可視化する

# データ読み込み
load("pca2.Rdata")
# 名前がついたベクトルはplotlyでは受け付けない
rownames(pca2$rotation) <- NULL

# Labelからも14_1875を削除
Label2 <- Label[setdiff(1:nrow(Label), grep("14_1875", Label[,2])), ]

pca_plotly2 <- data.frame(
Celltype=as.factor(Label2[,1]), PC1=pca2$rotation[,1], PC2=pca2$rotation[,2], PC3=pca2$rotation[,3])

pca_plotly2 %>%
	plot_ly(x=~PC1, y=~PC2, z=~PC3,
    type = "scatter3d", mode = "markers",
    text = as.vector(Label2[,2]),
    color =~Celltype
    )

kokitsuyuzaki/Level1_3_R.markdown