basic_t2019_1.md

簡単な計算の実行

3 + 8
30 / 23
3005 %% 3
sin(100)
log(20)

オブジェクト

x <- 3.5
x
x / 2

関数の調べ方（ヘルプ）

?log

• Description（概要）： 関数が何をしてくれるかを短くまとめた説明
• Usage（構文）：関数の呼び出しの書き方の例。
• Arguments（引数）：関数がとる引数。R が引数に渡すべき値として想定している情報のタイプ、その情報の使い方をまとめたリスト
• Details（詳細）：関数とその動作についての詳しい説明
• Value（戻り値）：関数を実行したときに返される情報についての説明
• See Also（参照）：関連する R 関数の短いリスト
• Examples（使用例）：関数の使用例のコード

オブジェクト名の大文字と小文字は区別される

x <- 3.5
x
X <- -1.2
X
x

作業ディレクトリの設定

windowsの場合

dir.create("C:/introR2019")
setwd("C:/introR2019")

mac の場合

dir.create("introR2019")
setwd("introR2019")

ワークスペース

ls()

データの読み書き

library(readr)
mat <- matrix(1:9, nrow=3)
colnames(mat) <- c("A","B","C")
df <- as.data.frame(mat)
write_csv(x=df, path="sample01.csv", col_names=TRUE)
write_tsv(x=df, path="sample01.tsv", col_names=TRUE)
mat1 <- read_csv("sample01.csv",col_names=TRUE)
mat2 <- read_tsv("sample01.tsv",col_names=TRUE)
mat1
mat2

演習

student <- c("georg", "niiyan", "paul", "sawada", "shima", "sito", "yshira")
math <- c(63, 61, 46, 55, 35, 91, 51)
sci <- c(32, 42, 51, 45, 11, 46, 33)
eng <- c(59, 58, 54, 63, 67, 55, 78)
group <- c("A", "A", "A", "B", "B", "B", "B")

mean(math)
student[max(sci)==sci]
mean(eng[group=="A"])
mean(eng[group=="B"])

Bflag <- group=="B"
B_stu <- student[Bflag]
B_stu[sort.list(math[Bflag],decreasing = TRUE)[2]]

いろいろなプロット

head(iris)
library(ggplot2)
ggplot(iris,aes(x=Sepal.Length,y=Sepal.Width,colour=Species))+
  geom_point()

ggplot(iris,aes(x=Sepal.Length, y=Sepal.Width,colour=Species))+
  geom_point()+
  scale_colour_manual(values = c("forestgreen","royalblue","orange2"))

統計グラフの色

ここから主成分分析

iris_pc <-prcomp(iris[,-5],center = FALSE)
iris_pc_df <-data.frame(iris_pc$x,Species=iris$Species)
ggplot(iris_pc_df,aes(x=PC1,y=PC2,colour=Species))+
  geom_point()

#検算
iris_pc <-prcomp(iris[,-5],center = FALSE)
z1 <-iris_pc$x
z2 <- as.matrix(iris[,-5]) %*% iris_pc$rotation
head(z1)
head(z2)

ここからクラスター分析

x <- iris[,-5]
d_iris <- dist(x)

dmat <- as.matrix(d_iris)
dmat[1:5,1:5]
sqrt(sum((x[1,]-x[2,])^2)) #サンプル1と2の距離
sqrt(sum((x[1,]-x[3,])^2)) #サンプル1と3の距離

hc_iris <-hclust(d_iris,method = "ward.D")
plot(hc_iris, labels  = FALSE, ann = FALSE)

iris_pc_df$cluster <- factor(cutree(hc_iris,k=3))
ggplot(iris_pc_df,aes(x=PC1,y=PC2,colour=cluster))+
  geom_point()

dplyr の基本

library(dplyr)
iris %>%
  group_by(Species) %>%
  summarise(mean(Sepal.Length))
iris %>%
  group_by(Species) %>%
  summarise(mean(Sepal.Length),sd(Sepal.Length))

#混同行列
iris_pc_df %>%
  group_by(Species,cluster) %>%
  summarise(length(PC1))

iris_pc_df %>%
  group_by(Species,cluster) %>%
  tally()

library(tidyr)
iris_pc_df %>%
  group_by(Species,cluster) %>%
  tally() %>%
  spread(cluster,n)

iris_pc_df %>%
  group_by(Species,cluster) %>%
  tally() %>%
  spread(cluster,n,fill=0)

エンテロタイピング

library(tidyverse)
dat <- read_tsv("MetaHIT_SangerSamples.genus.txt")

Y <- t(dat[-c(1:2),-1])
prout <- prcomp(Y,rank. = 2)

dY <- dist(Y)
clout <- hclust(dY,method = "ward.D")
clindex <- cutree(clout,k=3)

prdf <- data.frame(prout$x) %>%
  mutate(cluster=factor(clindex))

ggplot(prdf,aes(PC1,PC2,colour=cluster))+
  geom_point()+
  theme_bw()

colnames(Y) <- dat$X1[-c(1:2)]

dfY <- as.data.frame(Y) %>%
  mutate(cluster=factor(clindex))

ggplot(dfY,aes(cluster,Bacteroides))+
  geom_boxplot()

ggplot(dfY,aes(cluster,Ruminococcus))+
  geom_boxplot()

ggplot(dfY,aes(cluster,Prevotella))+
  geom_boxplot()

centerdf <-
  gather(dfY, genus, value,-cluster) %>%
  group_by(genus,cluster) %>%
  summarise(value=mean(value)) %>%
  group_by(cluster) %>%
  arrange(desc(value)) %>%
  mutate(order=row_number()) %>%
  group_by(genus) %>%
  filter(any(order<10))

ggplot(centerdf,aes(x=cluster,y=value))+
  geom_col()+
  facet_wrap(~genus,scales = "free")

生存時間分析

library(survival)
head(aml)
sf_aml <-survfit(Surv(time,status)~x,data = aml)
plot(sf_aml)

#シミュレーション
n <- 100
x<-rnorm(n,2,3)
d<-rep(1,n)
sfnorm <-survfit(Surv(x,d)~1) #カプランマイヤープロット
plot(sfnorm,conf.int = FALSE)
curve(1-pnorm(x,2,3),add=TRUE,col="red") #生存関数と対応
hist(x,freq = FALSE) #ヒストグラム
curve(dnorm(x,2,3),add=TRUE,col="red") #密度関数と対応
cdfnorm <- ecdf(x) #経験分布関数
curve(cdfnorm(x),-4,4)
curve(pnorm(x,2,3),add=TRUE,col="red") #分布関数と対応

gliomaSet の分析

データの読み込み

gliomaSet_express <-read_csv("gliomaSet_express.csv")
gliomaSet_surv <-read_csv("gliomaSet_surv.csv")
library(dplyr)
gliomaSet <- left_join(gliomaSet_surv,gliomaSet_express,by="X") #データを結合
#gliomaSet <- merge(gliomaSet_surv,gliomaSet_express,by="X") #でもよい
gliomaSet[1:7,1:7] #1から7列目の1から7行目を表示

いろいろなプロット

#いろいろなプロット
ggplot(gliomaSet,aes(x=g1,y=g2))+
  geom_point() #散布図
#plot(g2~g1,data=gliomaSet) #と同じ
ggplot(gliomaSet,aes(x=Age,fill=Gender))+
  geom_histogram() #ヒストグラム
ggplot(gliomaSet,aes(x=Age,colour=Gender))+
  geom_freqpoly() #ヒストグラムを折れ線で表示したもの
ggplot(gliomaSet,aes(x=Age,colour=Gender))+
  geom_freqpoly()+
  geom_rug() #データのある位置をバーコード状に表示したもの
ggplot(gliomaSet,aes(x=Gender,y=g3))+
  geom_boxplot() #箱ひげ図
#boxplot(g3~Gender,data = gliomaSet) #と同じ
ggplot(gliomaSet,aes(x=Gender,y=g3))+
  geom_violin() #バイオリンプロット
ggplot(gliomaSet,aes(x=g1,y=g2,colour=Gender))+
  geom_point() #散布図（色付き）

#help(package="ggplot2") #他に様々なグラフがある。ヘルプから探してみよう

ggplot(gliomaSet,aes(x=g1,y=g2,colour=Gender))+
  geom_point()+
  theme_bw() #白地のテーマ

ggplot(gliomaSet,aes(x=g1,y=g2,colour=Gender))+
  geom_point()+
  theme_bw(base_size = 20) #フォントサイズの変更

ggplot(gliomaSet,aes(x=g1,y=g2,colour=Gender))+
  geom_point()+
  scale_color_manual(values = c("royalblue","orange2")) #色を変更

ここから主成分分析　

#df_expression <- gliomaSet[,6:105]
df_expression <-select(gliomaSet,g1:g100) #g1~g100の列だけ取り出す
pc_gli <-prcomp(df_expression,center = FALSE) #主成分分析
z1 <-pc_gli$x #主成分得点
w <- pc_gli$rotation #重み
x <- as.matrix(df_expression) #行列として認識させる
z2 <- x %*% w #行列の掛け算
z1[1:7,1:5] #主成分得点
z2[1:7,1:5]　#主成分得点
#結果が一致するはず

pc_gli <- prcomp(df_expression) #中心が0になるようにして主成分分析
pc_gli_df <-data.frame(pc_gli$x) #主成分得点を取り出す

ggplot(pc_gli_df,aes(x=PC1,y=PC2))+ #第一、第二主成分の散布図
  geom_point()

ここからクラスター分析

d_expression <- dist(df_expression) #距離行列の作成
dmat <- as.matrix(d_expression) #行列として認識させる
dmat[1:5,1:5]
#距離の計算の仕方を確認
sqrt(sum((df_expression[1,]-df_expression[2,])^2)) #サンプル1と2の距離
sqrt(sum((df_expression[1,]-df_expression[3,])^2)) #サンプル1と3の距離

hc_expression <- hclust(d_expression, method = "ward.D")
plot(hc_expression) #樹形図の描画

cluster <- cutree(hc_expression,k=2) #クラスターを2つに分ける
head(cluster)
cluster <- factor(cluster) #ファクター（因子）として認識
head(cluster)

ggplot(pc_gli_df,aes(x=PC1,y=PC2,colour=cluster))+
  geom_point() #散布図

ggplot(pc_gli_df,aes(x=PC1,y=PC2,colour=cluster))+
  geom_point(size=3) + #散布図2
  theme_bw()

gliomaSet <- mutate(gliomaSet,cluster=cluster) #gliomaSetにclusterという列を追加
gliomaSet <- mutate(gliomaSet,Status2=factor(Status)) #gliomaSetにStatus2という列を追加
#gliomaSet$cluster <- cluster #でも良い
#gliomaSet$Status2 <- factor(gliomaSet$Status) #でも良い

ggplot(gliomaSet,aes(x=cluster,fill=Status2))+
  geom_bar() #棒グラフ

ggplot(gliomaSet,aes(x=cluster,fill=Status2))+
  geom_bar(position = "fill") #帯グラフ（積み上げ棒グラフを割合表示に）

ここから生存時間分析

library(survival)
#生存関数の推定はsurvfit関数
sf_gli <-survfit(Surv(Time,Status)~cluster,data = gliomaSet)
Surv(gliomaSet$Time,gliomaSet$Status)

library(broom)
sf_gli_df <- tidy(sf_gli) #データフレームとして使用
head(sf_gli_df)
ggplot(sf_gli_df,aes(x=time,y=estimate,colour=strata))+
  geom_step() #階段状の折れ線プロット

#ログランク検定
#帰無仮説：「生存時間に差はない」で検定
sd_gli <-survdiff(Surv(Time,Status)~cluster,data = gliomaSet_surv)
print(sd_gli)

以下複雑な集計

group_by(gliomaSet,cluster) %>%
  summarise(mean(g1))

out <- group_by(gliomaSet,cluster) %>%
  summarise_if(is.numeric,mean)