---
title: "[第21回](https://www.iu.a.u-tokyo.ac.jp/kadota/r_seq2.html#book_JSLAB_21)のウェブ資料"
date: '最終更新: `r format(Sys.Date(), "%Y/%m/%d")`'
output: 
  html_document:
    toc: TRUE
    toc_float: TRUE
---

-   [JSLAB21.Rmd](https://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB21.Rmd)\
-   [JSLAB21.html](https://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB21.html)

# はじめに（W1）

-   元となるRNA-seqカウントデータ([JSLAB18.csv](https://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB18.csv))\
    第20回はxlsxファイル([JSLAB18.xlsx](https://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB18.xlsx))でしたが、今回はcsvファイルを読み込むやり方の紹介を兼ねて[JSLAB18.csv](https://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB18.csv)を利用します。
-   [乳酸菌NGS連載第18回のPDF](https://www.iu.a.u-tokyo.ac.jp/kadota/JSLAB_18_kadota.pdf)\
-   [<i>Lactobacillus rhamnosus</i> GG (LGG)](https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?mode=Info&id=568703)\
-   [Bang et al., J Microbiol Biotechnol., 2018](http://www.ncbi.nlm.nih.gov/pubmed/30196592)\
-   乳酸菌NGS連載第18回の図5\
    ![](JSLAB18_fig5.png)

# パラメータの取得（W2）

入力ファイル([JSLAB18.csv](https://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB18.csv))の読込から、群ごとのカウントの平均($M$)とdispersion ($\varPhi$)を得る必要最小限のスクリプトです。

```{r, message=FALSE, error=TRUE}
source("https://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB21_func.R")
in_f <- "https://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB18.csv"
nk <- c(3, 3, 3)
A <- read.csv(in_f, row.names=1, stringsAsFactors=T)
T <- colSums(A)
Anorm <- sweep(x=A, MARGIN=2, STATS=1000000/T, FUN="*")
K <- length(nk)
L <- as.vector(mapply(rep, 1:K, nk))
M <- V <- Phi <- NULL
for(i in 1:K){
  M <- cbind(M, apply(X=Anorm[, L==i], MARGIN=1, FUN=mean))
  V <- cbind(V, apply(X=Anorm[, L==i], MARGIN=1, FUN=var))
  Phi <- cbind(Phi, apply(X=Anorm[, L==i], MARGIN=1, FUN=func_phi))
}
```

# NBモデルに従う乱数（W3）

## W3.1 $M$, $V$, $\varPhi$の概要

`head`関数を用いて、最初の4行分を表示させています。この結果を眺めて、（$\phi$ \> 0と$\phi \leq$ 0の両方をもつので$i$ = 2でもよかったのですがなんとなく）$i$ = 4に相当するEBG00001128509のパラメータを用いて説明しようと決めました。

```{r, message=FALSE, error=TRUE}
head(M, n=4)
head(V, n=4)
head(Phi, n=4)
```

## W3.2 `rnbinom`関数が正常に動作する場合（$\varPhi$ \> 0）

$m_{4,3}$ = 190.1133と$\phi_{4,3}$ = 0.056272をパラメータとして与え、8個の乱数を生成と分散の表示、そして真の値（$v_{4,3}$ = 2223.958）の表示まで行っています。

```{r, message=FALSE, error=TRUE}
i <- 4
k <- 3
M[i,k]
Phi[i,k]
set.seed(2023)
output <- rnbinom(n=8, mu=M[i,k], size=1/Phi[i,k])
output
var(output)
V[i,k]
```

次に、同条件で生成させる乱数のみ8個から1000000個に変更して分散を算出し、乱数が増えれば真の値（$v_{4,3}$ = 2223.958）に近づくことを確認します。

```{r, message=FALSE, error=TRUE}
set.seed(2023)
output <- rnbinom(n=1000000, mu=M[i,k], size=1/Phi[i,k])
var(output)
```

## W3.3 `rnbinom`関数のマニュアル

```{r, message=FALSE, error=TRUE}
?rnbinom
```

## W3.4 `?rnbinom`実行直後のスクショ

①のチャンク実行後は、以下に示すような関数マニュアルが②Helpタブ上に見られます。計4つの関数についての説明がなされており、目的の関数は③で見えています。③rnbinomは、`n`, `size`, `prob`, `mu`という計4つの引数（Arguments）を受けつけることがわかりますが、それぞれの説明は④に書かれています。\
![](https://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB21_W3.4.png)

## W3.5 `rnbinom`関数マニュアルのArguments部分

W3.4の続きとして、④Argumentsの⑤`size`に関する説明が見られるようにした状態です。`size`オプションで与える値は、⑥に明記されているように整数（integer）である必要はありませんが、正の値をもたねばなりません。\
![](https://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB21_W3.5.png)

## W3.6 `rnbinom`関数がエラーとなる場合（$\phi \leq$ 0）

$m_{4,2}$ = 166.7145と$\phi_{4,2}$ = -0.0005365914をパラメータとして与え、8個の乱数を生成と分散の表示まで行っています。

```{r, message=FALSE, error=TRUE}
i <- 4
k <- 2
M[i,k]
Phi[i,k]
set.seed(2023)
output <- rnbinom(n=8, mu=M[i,k], size=1/Phi[i,k])
output
var(output)
```

# シミュレーションデータ生成の基本形（W4）

## W4.1 第20回の図4

[第20回](https://www.iu.a.u-tokyo.ac.jp/kadota/r_seq2.html#book_JSLAB_20)の原稿はモノクロですので、こちらのほうが特に(d)はわかりやすいかと思います。\
![](https://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB20_fig4.png)

## W4.2 G1群のみの$\varPhi$ \> 0を満たす経験分布を得る

W2のチャンクを実行しておけば、このチャンクをエラーなく実行できます。`k <- 1`としてG1群に限定し、計2949遺伝子のベクトルに対して、`Phi[,k] > 0`という条件判定を行い、条件を満たす要素をTRUE、そうでなければFALSEという情報からなる`obj`を作成しています。`sum(obj)`は、`obj`がTRUEの要素数を返すので、出力結果が2511というのは妥当です。その後の`head(obj)`では、`obj`の最初の6個の要素を出力しています。5番目の要素がFALSEになっていることがわかります。`head(Phi[, k])`の結果と合わせると、確かに5番目の要素が条件を満たしていないことがわかります。`sub_M <- M[obj, k]`は、2949行からなる`M[, k]`の中から、`obj`がTRUEとなる行のみを`sub_M`に格納せよというコマンドです。結果として得られる`sub_M`は、2511個の要素からなる数値ベクトルになります。`sub_Phi <- Phi[obj, k]`についても同じノリで、平均の行列`M`の部分がdispersionの行列`Phi`になっているだけです。**ここのスクリプトは、わかりやすさ（説明のしやすさ）重視で汎用性は低い**のでご注意ください。

```{r, message=FALSE, error=TRUE}
k <- 1
obj <- (Phi[,k] > 0)
sum(obj)
head(obj)
head(Phi[, k])
sub_M <- M[obj, k]
sub_Phi <- Phi[obj, k]
```

## W4.3 経験分布の記載例

[Kadota et al., Algorithms Mol Biol, 2012](https://pubmed.ncbi.nlm.nih.gov/22475125/)では、以下のような記載をしています。\
To mimic this mean-dispersion relationship in the simulation, we used an empirical distribution of these values ($\mu$ and $\phi$) calculated from Arabidopsis data available in the [NBPSeq](https://CRAN.R-project.org/package=NBPSeq) package.\
2012年の論文は、[TCC](https://bioconductor.org/packages/TCC/)に実装されているDEGES正規化法の原著論文です。TCCで生成しているシミュレーションデータは、モデル植物であるシロイヌナズナのリアルカウントデータから経験分布を得ているということですね。

## W4.4 生成予定の遺伝子数分のパラメータ情報を取得

遺伝子数が$G_{sim}$ = 4の例を示します。`s_num`は、`sample`関数を用いて、`1:length(sub_M)`の整数ベクトルの中から、`S_sim`で指定した個数を、復元抽出(`replace=TRUE`)した結果であり、計4個の整数の要素からなる数値ベクトルです。`sub_M[s_num]`は、W4.2で作成した`sub_M`の要素の中から、`s_num`で指定した要素の情報のみを抽出していることに相当します。`sub_Phi[s_num]`についても同様です。`match`関数は、この場合は`table=rownames(M))`で指定している行列`M`の行名情報(つまり`rownames(M)`)の文字列ベクトルの中から、`x=names(sub_M[s_num])`で指定した文字列ベクトル中の要素と一致するインデックス(要素番号)情報を返します。最終行の`c("A", "B", "C", "D", "E")`の中から`c("D", "B")`中の要素と一致するインデックス情報として4と2が返されているのをみると全体像がよくわかると思います。

```{r, message=FALSE, error=TRUE}
G_sim <- 4
set.seed(2023)
s_num <- sample(1:length(sub_M), G_sim, replace=TRUE)
s_num
sub_M[s_num]
sub_Phi[s_num]
match(x=names(sub_M[s_num]), table=rownames(M))
match(x=c("D", "B"), table=c("A", "B", "C", "D", "E"))
```

## W4.5 乱数を生成(べた書き版)

遺伝子数が$G_{sim}$ = 4、反復数が$n_{sim}$ = 5の例を示します。`s_num[1]` = 1909、`s_num[2]` = 1488、`s_num[3]` = 2479、`s_num[4]` = 1385が代入されて、このインデックスに対応する平均とdispersionを入力として乱数が`n_sim`個分だけ生成されます。ここはあくまでもW4.6が難解であったヒト用です。

```{r, message=FALSE, error=TRUE}
n_sim <- 5
set.seed(2023)
i <- 1
rnbinom(n=n_sim, mu=sub_M[s_num[i]], size=1/sub_Phi[s_num[i]])
i <- 2
rnbinom(n=n_sim, mu=sub_M[s_num[i]], size=1/sub_Phi[s_num[i]])
i <- 3
rnbinom(n=n_sim, mu=sub_M[s_num[i]], size=1/sub_Phi[s_num[i]])
i <- 4
rnbinom(n=n_sim, mu=sub_M[s_num[i]], size=1/sub_Phi[s_num[i]])
```

## W4.6 乱数を生成(forループ版)

遺伝子数が$G_{sim}$ = 4、反復数が$n_{sim}$ = 5の例を示します。`simdata <- NULL`は、ただの"これから`simdata`という名前のオブジェクトに具体的な数値情報を入れていくので、まずは場所の確保だけしときますよ(プレースホルダの作成)"という宣言です。`for(i in 1:G_sim)`は、`i`という変数を使って、`1`から`G_sim`まで1づつカウントアップし、`{}`で囲った範囲のコマンドを実行せよという意味です。`hoge`の中身は、ただの`rnbinom`関数実行結果であり、`n_sim`個の乱数からなる数値ベクトルです。その次の行の`rbind`関数は、既に存在する`simdata`オブジェクトの下の行に`hoge`を追加せよ、という指令です。たとえば、`i` = 1のときは、`simdata <- rbind(simdata, hoge)`実行結果として得られる`simdata`の中身は1行5列の行列となり、その中身は`s_num[1]` = 1909で乱数を発生させた結果の`hoge`と同じです。`i` = 2のときは、既に存在する1行5列の`simdata`の下の行に`s_num[2]` = 1488で乱数を発生させた結果のhogeが追加され2行5列の`simdata`になる、といった具合です。

```{r, message=FALSE, error=TRUE}
G_sim <- 4
n_sim <- 5
set.seed(2023)
simdata <- NULL
for(i in 1:G_sim){
  hoge <- rnbinom(n=n_sim, mu=sub_M[s_num[i]], size=1/sub_Phi[s_num[i]])
  simdata <- rbind(simdata, hoge)
}
simdata
```

## W4.7 行ごとの分散

まず、行列`simdata`に対して、行名と列名を付加しています。その後、`apply`関数を用いて、行ごと（`MARGIN=1`）に、分散(`FUN=var`)を計算しています。他には、たとえば列ごとに総和を計算したい場合は、`MARGIN=2`と`FUN=sum`とします。`posi`は、W4.4で得た$\varPhi$ \> 0を満たすG1群の2511個からランダム抽出した`s_num`から、遺伝子名情報の文字列の一致を頼りに同定した、大元の2949遺伝子の中でのインデックス（要素番号）情報です。最後にその`posi`の行かつG1群の列（つまり`1`列目）の分散情報を表示させています。

```{r, message=FALSE, error=TRUE}
rownames(simdata) <- paste("gene", 1:G_sim, sep="")
colnames(simdata) <- paste("rep", 1:n_sim, sep="_")
simdata
apply(X=simdata, MARGIN=1, FUN=var)
posi <- match(x=names(sub_M[s_num]), table=rownames(M))
posi
V[posi, 1]
```

## W4.8 反復数を増やして検証

反復数が$n_{sim}$ = 500000の例です。これくらい反復数を増やすと、正解にかなり近い値が得られることが実感できます。

```{r, message=FALSE, error=TRUE}
G_sim <- 4
n_sim <- 500000
set.seed(2023)
simdata <- NULL
for(i in 1:G_sim){
  hoge <- rnbinom(n=n_sim, mu=sub_M[s_num[i]], size=1/sub_Phi[s_num[i]])
  simdata <- rbind(simdata, hoge)
}
apply(X=simdata, MARGIN=1, FUN=var)
posi <- match(x=names(sub_M[s_num]), table=rownames(M))
V[posi, 1]
```

## W4.9 復元抽出の意味

ここのチャンクは、前半がW4.2、そして後半がW4.4を合わせたものになっています。乱数取得のための元情報`length(sub_M)`が`G_sim`未満であれば、sample関数実行時に非復元抽出（`replace=FALSE`）で行うとエラーが出ることを示しています。

```{r, message=FALSE, error=TRUE}
k <- 1
obj <- (Phi[,k] > 0)
sub_M <- M[obj, k]
sub_Phi <- Phi[obj, k]
length(sub_M)

G_sim <- 3000
set.seed(2023)
s_num <- sample(1:length(sub_M), G_sim, replace=FALSE)

```

## W4.10 遺伝子数10000でシミュレーションデータを生成

`G_sim` = 10000、`n_sim` = 5でシミュレーションデータを生成し、行数と列数を確認するところまでです。W2実行後であれば、ここのチャンクが正常に動作するはずです。

```{r, message=FALSE, error=TRUE}
k <- 1
G_sim <- 10000
n_sim <- 5
set.seed(2023)
obj <- (Phi[,k] > 0)
sub_M <- M[obj, k]
sub_Phi <- Phi[obj, k]
s_num <- sample(1:length(sub_M), G_sim, replace=TRUE)
simdata <- NULL
for(i in 1:G_sim){
  hoge <- rnbinom(n=n_sim, mu=sub_M[s_num[i]], size=1/sub_Phi[s_num[i]])
  simdata <- rbind(simdata, hoge)
}
rownames(simdata) <- paste("gene", 1:G_sim, sep="")
colnames(simdata) <- paste("rep", 1:n_sim, sep="_")
dim(simdata)
```

# 平均-分散プロットの基本形（W5）

[第20回](https://www.iu.a.u-tokyo.ac.jp/kadota/r_seq2.html#book_JSLAB_20)のW9.7に対応するプロットの基本形です。

## W5.1 G1群のシミュレーションデータ

リアルデータと数を揃えるべく、`G_sim` = 2949、`n_sim` = 3でシミュレーションデータを生成し、平均-分散プロットを描画します。2949行×3列のシミュレーションデータのオブジェクト`simdata`を作成する基本形はW4.10と同じです。その後は、[第20回](https://www.iu.a.u-tokyo.ac.jp/kadota/r_seq2.html#book_JSLAB_20)のW9.7と似たようなノリでプロット作業を行っています。具体的には、`simdata`の遺伝子（行）ごとの平均`M_sim`と分散`V_sim`を作成し、それらを`data.frame`関数を用いて列方向でマージした結果を`d`オブジェクトに格納しています。そして、[第20回](https://www.iu.a.u-tokyo.ac.jp/kadota/r_seq2.html#book_JSLAB_20)のW9.7.2をテンプレートとしてプロットしています。

```{r, message=FALSE, error=TRUE}
k <- 1
G_sim <- 2949
n_sim <- 3
set.seed(2023)
obj <- (Phi[,k] > 0)
sub_M <- M[obj, k]
sub_Phi <- Phi[obj, k]
s_num <- sample(1:length(sub_M), G_sim, replace=TRUE)
simdata <- NULL
for(i in 1:G_sim){
  hoge <- rnbinom(n=n_sim, mu=sub_M[s_num[i]], size=1/sub_Phi[s_num[i]])
  simdata <- rbind(simdata, hoge)
}
rownames(simdata) <- paste("gene", 1:G_sim, sep="")
colnames(simdata) <- paste("rep", 1:n_sim, sep="_")
dim(simdata)

simdata_norm <- sweep(x=simdata, MARGIN=2, STATS=1000000/colSums(simdata), FUN="*")
M_sim <- apply(X=simdata_norm, MARGIN=1, FUN=mean)
V_sim <- apply(X=simdata_norm, MARGIN=1, FUN=var)
Phi_sim <- apply(X=simdata_norm, MARGIN=1, FUN=func_phi)
d <- data.frame(MEAN = M_sim,
                VARIANCE = V_sim,
                DISPERSION = Phi_sim)

library(ggplot2)
coloor <- "magenta"
axis_range <- c(1e-01, 1e+7)
g <- ggplot(d, aes(x=MEAN, y=VARIANCE)) +
     geom_point(show.legend = FALSE, alpha=0.45, shape=19, size=0.4, col=coloor) +
     scale_x_log10() +
     scale_y_log10() +
     theme_bw() +
     geom_abline(intercept=0, slope=1, color="blue") +
     coord_cartesian(xlim=axis_range, ylim=axis_range) +
     labs(x="Mean",y="Variance")
plot(g)
d1_sim <- d
g1_sim <- g
```

## W5.2 G1群のリアルデータ

実質的に[第20回](https://www.iu.a.u-tokyo.ac.jp/kadota/r_seq2.html#book_JSLAB_20)のW9.7.2と同じですが、`data.frame`関数を用いてマージする部分を最初からG1群のみのデータにするなど若干変更しています。今回のW2を実行していれば、`M`と`V`のオブジェクトが利用可能なはずですので、ここのチャンクを実行可能です。なお、`k <- 1`は、G1群のみに限定するという意味で、抽出したい列のインデックス情報に相当します。

```{r, message=FALSE, error=TRUE}
k <- 1
d <- data.frame(MEAN = M[, k],
                VARIANCE = V[, k])

library(ggplot2)
coloor <- "magenta"
axis_range <- c(1e-01, 1e+7)
g <- ggplot(d, aes(x=MEAN, y=VARIANCE)) +
     geom_point(show.legend = FALSE, alpha=0.45, shape=19, size=0.4, col=coloor) +
     scale_x_log10() +
     scale_y_log10() +
     theme_bw() +
     geom_abline(intercept=0, slope=1, color="blue") +
     coord_cartesian(xlim=axis_range, ylim=axis_range) +
     labs(x="Mean",y="Variance")
plot(g)
d1_real <- d
g1_real <- g
```

## W5.3 G2群のシミュレーションデータ

W5.1と基本的に同じで、G2群について行っています。

```{r, message=FALSE, error=TRUE}
k <- 2
G_sim <- 2949
n_sim <- 3
set.seed(2023)
obj <- (Phi[,k] > 0)
sub_M <- M[obj, k]
sub_Phi <- Phi[obj, k]
s_num <- sample(1:length(sub_M), G_sim, replace=TRUE)
simdata <- NULL
for(i in 1:G_sim){
  hoge <- rnbinom(n=n_sim, mu=sub_M[s_num[i]], size=1/sub_Phi[s_num[i]])
  simdata <- rbind(simdata, hoge)
}
rownames(simdata) <- paste("gene", 1:G_sim, sep="")
colnames(simdata) <- paste("rep", 1:n_sim, sep="_")
dim(simdata)

simdata_norm <- sweep(x=simdata, MARGIN=2, STATS=1000000/colSums(simdata), FUN="*")
M_sim <- apply(X=simdata_norm, MARGIN=1, FUN=mean)
V_sim <- apply(X=simdata_norm, MARGIN=1, FUN=var)
Phi_sim <- apply(X=simdata_norm, MARGIN=1, FUN=func_phi)
d <- data.frame(MEAN = M_sim,
                VARIANCE = V_sim,
                DISPERSION = Phi_sim)

library(ggplot2)
coloor <- "navy"
axis_range <- c(1e-01, 1e+7)
g <- ggplot(d, aes(x=MEAN, y=VARIANCE)) +
     geom_point(show.legend = FALSE, alpha=0.45, shape=19, size=0.4, col=coloor) +
     scale_x_log10() +
     scale_y_log10() +
     theme_bw() +
     geom_abline(intercept=0, slope=1, color="blue") +
     coord_cartesian(xlim=axis_range, ylim=axis_range) +
     labs(x="Mean",y="Variance")
plot(g)
d2_sim <- d
g2_sim <- g
```

## W5.4 G2群のリアルデータ

W5.2と基本的に同じで、G2群について行っています。

```{r, message=FALSE, error=TRUE}
k <- 2
d <- data.frame(MEAN = M[, k],
                VARIANCE = V[, k])

library(ggplot2)
coloor <- "navy"
axis_range <- c(1e-01, 1e+7)
g <- ggplot(d, aes(x=MEAN, y=VARIANCE)) +
     geom_point(show.legend = FALSE, alpha=0.45, shape=19, size=0.4, col=coloor) +
     scale_x_log10() +
     scale_y_log10() +
     theme_bw() +
     geom_abline(intercept=0, slope=1, color="blue") +
     coord_cartesian(xlim=axis_range, ylim=axis_range) +
     labs(x="Mean",y="Variance")
plot(g)
d2_real <- d
g2_real <- g
```

## W5.5 複数プロットのレイアウト

作成済みの描画用オブジェクトである`g1_real`、`g1_sim`、`g2_real`、`g2_sim`が4つ存在するという前提で行います。基本的には表示させたい順に記載していくだけです。最後に`plot_layout(ncol=2)`とやっていますが、これは列数を2つに限定せよという指令です。

```{r, message=FALSE, error=TRUE, warning=FALSE}
library(patchwork)
g1_real + g2_real + g1_sim + g2_sim + plot_layout(ncol=2)
```

# 平均-分散プロットの発展形（W6）

[第20回](https://www.iu.a.u-tokyo.ac.jp/kadota/r_seq2.html#book_JSLAB_20)のW9.8に対応するプロットの発展形です。ggplot関数実行時に、`geom_abline(intercept=0.4, slope=1.6, color="gray")`として追加された灰色の実線は、計4種類のプロットを眺めて、分散が全体的に小さいG2群のプロットがほぼこの灰色の実線の右下となるようなパラメータ（つまり`intercept=0.4`と`slope=1.6`）を目視で定めています。

## W6.1 G1群のシミュレーションデータ

W5.1と[第20回](https://www.iu.a.u-tokyo.ac.jp/kadota/r_seq2.html#book_JSLAB_20)のW9.8.2をテンプレートとしてプロットしています。

```{r, message=FALSE, error=TRUE}
k <- 1
G_sim <- 2949
n_sim <- 3
set.seed(2023)
obj <- (Phi[,k] > 0)
sub_M <- M[obj, k]
sub_Phi <- Phi[obj, k]
s_num <- sample(1:length(sub_M), G_sim, replace=TRUE)
simdata <- NULL
for(i in 1:G_sim){
  hoge <- rnbinom(n=n_sim, mu=sub_M[s_num[i]], size=1/sub_Phi[s_num[i]])
  simdata <- rbind(simdata, hoge)
}
rownames(simdata) <- paste("gene", 1:G_sim, sep="")
colnames(simdata) <- paste("rep", 1:n_sim, sep="_")
dim(simdata)

simdata_norm <- sweep(x=simdata, MARGIN=2, STATS=1000000/colSums(simdata), FUN="*")
M_sim <- apply(X=simdata_norm, MARGIN=1, FUN=mean)
V_sim <- apply(X=simdata_norm, MARGIN=1, FUN=var)
Phi_sim <- apply(X=simdata_norm, MARGIN=1, FUN=func_phi)
d <- data.frame(MEAN = M_sim,
                VARIANCE = V_sim,
                DISPERSION = Phi_sim)
Vlow <- sum(d[,1] > d[,2])
Vhigh <- sum(d[,1] < d[,2])

library(ggplot2)
coloor <- "magenta"
axis_range <- c(1e-01, 1e+7)
g <- ggplot(d, aes(x=MEAN, y=VARIANCE)) +
     geom_point(show.legend = FALSE, alpha=0.45, shape=19, size=0.4, col=coloor) +
     scale_x_log10() +
     scale_y_log10() +
     theme_bw() +
     geom_abline(intercept=0, slope=1, color="blue") +
     geom_abline(intercept=0.4, slope=1.6, color="gray") +
     coord_cartesian(xlim=axis_range, ylim=axis_range) +
     labs(x="Mean",y="Variance") +
     annotate("text", x=1e+06, y=1e+00, label=Vlow) +
     annotate("text", x=1e+00, y=1e+06, label=Vhigh)
plot(g)
d1_sim <- d
g1_sim <- g
```

## W6.2 G1群のリアルデータ

W5.2と[第20回](https://www.iu.a.u-tokyo.ac.jp/kadota/r_seq2.html#book_JSLAB_20)のW9.8.2をテンプレートとしてプロットしています。

```{r, message=FALSE, error=TRUE}
k <- 1
d <- data.frame(MEAN = M[, k],
                VARIANCE = V[, k],
                DISPERSION = Phi[, k])
Vlow <- sum(d[,1] > d[,2])
Vhigh <- sum(d[,1] < d[,2])

library(ggplot2)
coloor <- "magenta"
axis_range <- c(1e-01, 1e+7)
g <- ggplot(d, aes(x=MEAN, y=VARIANCE)) +
     geom_point(show.legend = FALSE, alpha=0.45, shape=19, size=0.4, col=coloor) +
     scale_x_log10() +
     scale_y_log10() +
     theme_bw() +
     geom_abline(intercept=0, slope=1, color="blue") +
     geom_abline(intercept=0.4, slope=1.6, color="gray") +
     coord_cartesian(xlim=axis_range, ylim=axis_range) +
     labs(x="Mean",y="Variance") +
     annotate("text", x=1e+06, y=1e+00, label=Vlow) +
     annotate("text", x=1e+00, y=1e+06, label=Vhigh)
plot(g)
d1_real <- d
g1_real <- g
```

## W6.3 G2群のシミュレーションデータ

W5.3と[第20回](https://www.iu.a.u-tokyo.ac.jp/kadota/r_seq2.html#book_JSLAB_20)のW9.8.2をテンプレートとしてプロットしています。

```{r, message=FALSE, error=TRUE}
k <- 2
G_sim <- 2949
n_sim <- 3
set.seed(2023)
obj <- (Phi[,k] > 0)
sub_M <- M[obj, k]
sub_Phi <- Phi[obj, k]
s_num <- sample(1:length(sub_M), G_sim, replace=TRUE)
simdata <- NULL
for(i in 1:G_sim){
  hoge <- rnbinom(n=n_sim, mu=sub_M[s_num[i]], size=1/sub_Phi[s_num[i]])
  simdata <- rbind(simdata, hoge)
}
rownames(simdata) <- paste("gene", 1:G_sim, sep="")
colnames(simdata) <- paste("rep", 1:n_sim, sep="_")
dim(simdata)

simdata_norm <- sweep(x=simdata, MARGIN=2, STATS=1000000/colSums(simdata), FUN="*")
M_sim <- apply(X=simdata_norm, MARGIN=1, FUN=mean)
V_sim <- apply(X=simdata_norm, MARGIN=1, FUN=var)
Phi_sim <- apply(X=simdata_norm, MARGIN=1, FUN=func_phi)
d <- data.frame(MEAN = M_sim,
                VARIANCE = V_sim,
                DISPERSION = Phi_sim)
Vlow <- sum(d[,1] > d[,2])
Vhigh <- sum(d[,1] < d[,2])

library(ggplot2)
coloor <- "navy"
axis_range <- c(1e-01, 1e+7)
g <- ggplot(d, aes(x=MEAN, y=VARIANCE)) +
     geom_point(show.legend = FALSE, alpha=0.45, shape=19, size=0.4, col=coloor) +
     scale_x_log10() +
     scale_y_log10() +
     theme_bw() +
     geom_abline(intercept=0, slope=1, color="blue") +
     geom_abline(intercept=0.4, slope=1.6, color="gray") +
     coord_cartesian(xlim=axis_range, ylim=axis_range) +
     labs(x="Mean",y="Variance") +
     annotate("text", x=1e+06, y=1e+00, label=Vlow) +
     annotate("text", x=1e+00, y=1e+06, label=Vhigh)
plot(g)
d2_sim <- d
g2_sim <- g
```

## W6.4 G2群のリアルデータ

W5.4と[第20回](https://www.iu.a.u-tokyo.ac.jp/kadota/r_seq2.html#book_JSLAB_20)のW9.8.2をテンプレートにしています。

```{r, message=FALSE, error=TRUE}
k <- 2
d <- data.frame(MEAN = M[, k],
                VARIANCE = V[, k],
                DISPERSION = Phi[, k])
Vlow <- sum(d[,1] > d[,2])
Vhigh <- sum(d[,1] < d[,2])

library(ggplot2)
coloor <- "navy"
axis_range <- c(1e-01, 1e+7)
g <- ggplot(d, aes(x=MEAN, y=VARIANCE)) +
     geom_point(show.legend = FALSE, alpha=0.45, shape=19, size=0.4, col=coloor) +
     scale_x_log10() +
     scale_y_log10() +
     theme_bw() +
     geom_abline(intercept=0, slope=1, color="blue") +
     geom_abline(intercept=0.4, slope=1.6, color="gray") +
     coord_cartesian(xlim=axis_range, ylim=axis_range) +
     labs(x="Mean",y="Variance") +
     annotate("text", x=1e+06, y=1e+00, label=Vlow) +
     annotate("text", x=1e+00, y=1e+06, label=Vhigh)
plot(g)
d2_real <- d
g2_real <- g
```

## W6.5 複数プロットのレイアウトの基本形

[第20回](https://www.iu.a.u-tokyo.ac.jp/kadota/r_seq2.html#book_JSLAB_20)のW9.7.3をテンプレートにしています。作成済みの描画用オブジェクトである`g1_real`、`g1_sim`、`g2_real`、`g2_sim`が4つ存在するという前提で行います。基本的には表示させたい順に記載していくだけです。最後に`plot_layout(ncol=2)`とやっていますが、これは列数を2つに限定せよという指令です。

```{r, message=FALSE, error=TRUE, warning=FALSE}
g1_real + g2_real + g1_sim + g2_sim + plot_layout(ncol=2)
```

## W6.6 複数プロットのレイアウトの本番

[第20回](https://www.iu.a.u-tokyo.ac.jp/kadota/r_seq2.html#book_JSLAB_20)のW9.8.3をテンプレートにしています。これが図4の原形です。

```{r, message=FALSE, error=TRUE, warning=FALSE}
g <- g1_real + g2_real + g1_sim + g2_sim + plot_layout(ncol=2)
ggsave("JSLAB21.w6.6.png", g, units="px", width=1900, height=1760, device="png")
```

得られた[JSLAB21.w6.6.png](JSLAB21.w6.6.png)は、以下に見えているものと同じになります。 ![](JSLAB21.w6.6.png)

## W6.7 dispersionパラメータ$\varPhi$の値の分布

$\varPhi$ \> 0という条件を満たす$\varPhi$の分布をG1群(`k <- 1`)とG2群(`k <- 2`)それぞれについて示しています。

```{r, message=FALSE, error=TRUE, warning=FALSE}
k <- 1
obj <- (Phi[,k] > 0)
summary(Phi[obj, k])
k <- 2
obj <- (Phi[,k] > 0)
summary(Phi[obj, k])
```

G1群のリアルデータの情報は`d1_real`に、G2群のリアルデータの情報は`d2_real`に、G1群のシミュレーションデータの情報は`d1_sim`に、G2群のシミュレーションデータの情報は`d2_sim`に、それぞれ格納しているので、以下のような書き方でもかまいません。一旦`hoge`というオブジェクトにコピーして使用することで、その下の`obj`を得るコマンドや`summary`関数部分のコマンドを同じにして使いまわせるという利点があります。

```{r, message=FALSE, error=TRUE, warning=FALSE}
hoge <- d1_real
obj <- (hoge$DISPERSION > 0)
summary(hoge$DISPERSION[obj])

hoge <- d2_real
obj <- (hoge$DISPERSION > 0)
summary(hoge$DISPERSION[obj])

hoge <- d1_sim
obj <- (hoge$DISPERSION > 0)
summary(hoge$DISPERSION[obj])

hoge <- d2_sim
obj <- (hoge$DISPERSION > 0)
summary(hoge$DISPERSION[obj])
```

# シミュレーションデータ生成の発展形（W7）

## W7.1 $\varPhi$ \> 0を満たす経験分布を得る

W4.2の発展形として、計3群全ての平均$M$とdispersion$\varPhi$の情報をシミュレーションデータ作成時の経験分布として利用します。理解のしやすさを重視して、W4.2とR4.6を組み合わせたようなスクリプトにしています。W2のチャンクを実行しておけば、このチャンクをエラーなく実行できます。$K$は群数です。`append`は、既存のベクトルの末尾に要素を追加する関数であり、`c`という関数と基本的に同じ機能です。`rbind`が既存の行列の行末に新たな行を追加する機能をもっていたことを思い出せば、そのベクトル版だと納得できると思います。`for`で、`k`番目の群ごとに、最大`K`までループを回し、dispersion情報を格納した行列`Phi`の`k`番目の列について、値が0より大きい行をTRUE、それ以外をFALSEとした論理値ベクトル`obj`を得たのち、その情報を用いてフィルタリングをしています。最終的に得られる`sub_M`と`sub_Phi`の要素数は6026であり、W4.1または[第20回](https://www.iu.a.u-tokyo.ac.jp/kadota/r_seq2.html#book_JSLAB_20)の図4(d)で見えている数値と一致していることがわかります。

```{r, message=FALSE, error=TRUE, warning=FALSE}
K <- 3
set.seed(2023)
sub_M <- sub_Phi <- NULL
for(k in 1:K){
  obj <- (Phi[,k] > 0)
  sub_M <- append(sub_M, M[obj, k])
  sub_Phi <- append(sub_Phi, Phi[obj, k])
}
length(sub_M)
length(sub_Phi)
```

## W7.2 全てがnonDEGの数値行列を生成

W4.10をテンプレートとしています。遺伝子数`G_sim`と群ごとの反復数`nk`を入力として与えて、数値行列`simdata`を出力として得ています。最終行で、念のため`write.csv`関数を用いて作成したシミュレーションデータを[JSLAB21_w7.2.csv](JSLAB21_w7.2.csv)という名前で保存しています。

```{r, message=FALSE, error=TRUE, warning=FALSE}
G_sim <- 10000
nk <- c(3, 3)
n_sim <- sum(nk)
set.seed(2023)
s_num <- simdata <- NULL
for(i in 1:G_sim){
  hoge1 <- sample(1:length(sub_M), 1, replace=TRUE)
  hoge2 <- rnbinom(n=n_sim, mu=sub_M[hoge1], size=1/sub_Phi[hoge1])
  s_num <- append(s_num, hoge1)
  simdata <- rbind(simdata, hoge2)
}
rownames(simdata) <- paste("gene", 1:G_sim, sep="")
colnames(simdata) <- c(paste("G1_rep", 1:nk[1], sep=""), paste("G2_rep", 1:nk[2], sep=""))
dim(simdata)
head(simdata)
write.csv(x=simdata, file="JSLAB21_w7.2.csv")
```

## W7.3 log2(G2/G1)の算出

W7.2で作成した`simdata`オブジェクトがあるという前提で記載しています。群ごとの反復数の情報`nk`を入力として与えてクラスラベル情報`L`を作成し、それに基づいて群ごとのカウントの平均を`M_G1`および`M_G2`として得ています。その後、底が2のlog比を`log2(M_G2/M_G1)`として算出した結果を`logratio`というオブジェクトに格納しています。`summary(logratio)`実行結果には-Inf、NaN、そしてInfが表示されていますが、これらはそれぞれG1群が0、両群とも0、そしてG2群が0の場合に相当します。それゆえ、`logratio`のベクトルに対して、有限かどうかを`is.finite`関数で判定し、有限である（InfやNaNやNAではない）要素のみで`summary`関数を実行しなおしているのが`summary(logratio[obj])`になります。`logratio[!obj]`は、有限ではない（InfやNaNやNAである）要素を表示していることになります。`!obj`とすることでTRUE or FALSEが入れ替わることを利用しています。こんな感じで実際のところどの遺伝子が悪さをしているのかを知ることができます。`length(logratio[!obj])`は、有限ではない（InfやNaNやNAである）要素を表示していることになります。

```{r, message=FALSE, error=TRUE, warning=FALSE}
nk <- c(3, 3)
K <- length(nk)
L <- as.vector(mapply(rep, 1:K, nk))
hoge <- simdata
M_G1 <- apply(X=hoge[, L==1], MARGIN=1, FUN=mean)
M_G2 <- apply(X=hoge[, L==2], MARGIN=1, FUN=mean)
logratio <- log2(M_G2/M_G1)
summary(logratio)

obj <- is.finite(logratio)
summary(logratio[obj])
logratio[!obj]
length(logratio[!obj])
```

## W7.4 ヒストグラム（DEGなし）

W7.3で作成した、log2(M_G2/M_G1)の計算結果に相当する`logratio`オブジェクトが存在するという前提で記載しています。最初に、わざわざ`data.frame`形式にしているのは、[ggplot2](https://ja.wikipedia.org/wiki/Ggplot2)でのを意識しているためり、log2(M_G2/M_G1)の計算の実体である`logratio`オブジェクトの中身を`log2ratio`という列名で取り扱えるようにして、`d`というオブジェクトに格納しています。そうしているからこそ、`ggplot`関数の部分で、`aes(x=log2ratio)`のように書いて指定できるのです。このあたりのノリは、[第20回](https://www.iu.a.u-tokyo.ac.jp/kadota/r_seq2.html#book_JSLAB_20)のW8と同じです。

```{r, message=FALSE, error=TRUE, warning=FALSE}
d <- data.frame(log2ratio = logratio)
library(ggplot2)
g <- ggplot(d, aes(x=log2ratio)) +
     geom_histogram()
plot(g)
```

<br> プロットの横軸のラベルが原点の0よりも右側にずれている印象を受けたかもしれませんが、これは`summary(logratio[obj])`実行結果を見返していただければ納得できると思います。つまり、有限である（InfやNaNやNAではない）要素のみでsummary関数を実行した結果の最小値が-5.028408、そして最大値が6.203428だからです。たとえば、以下のように、軸の取りうる範囲を`c(-4, 4)`として与えれば、プロットの横軸のラベルがピタリと原点の0付近に位置するようになります。

```{r, message=FALSE, error=TRUE, warning=FALSE}
d <- data.frame(log2ratio = logratio)
library(ggplot2)
axis_range <- c(-4, 4)
g <- ggplot(d, aes(x=log2ratio)) +
     geom_histogram() + 
     coord_cartesian(xlim=axis_range)
plot(g)
```

<br> 発現変動が2倍L未満に収まっている遺伝子数を調べています。2倍未満というのは、G2群がG1群に比べて1/2倍より大きく2倍より小さい範囲ということを意味します。これは、log2スケールでは、log2(1/2)より大きくlog2(2/1)より小さい範囲ということになり、結局-1より大きく+1より小さい範囲内にある遺伝子数をカウントすればよいということになります。つまり、以下のような条件判定でよいことになります。

```{r, message=FALSE, error=TRUE, warning=FALSE}
obj <- ((logratio[obj] > -1) & (logratio[obj] < 1))
sum(obj)
```

## W7.5 DEGの導入

W7.2で作成した`simdata`オブジェクトがあるという前提で記載しています。群ごとの反復数の情報`nk`を入力として与えてクラスラベル情報`L`を作成し、それに基づいてどの群を`FC`倍するかを与えています。ここでは、`simdata`と同じ行数と列数の、どの要素を何倍するかという情報からなる`fc_mat`を作成したのち、`simdata`と`fc_mat`の行列同士を掛けることで目的を達成しています。ここでは**PDEG \> 0かつP1 \< 1.0であるという前提でシンプルに記載**しています。なお、`(G_sim*PDEG*P1)` = 2250、`(G_sim*PDEG*P1+1)` = 2251、`(G_sim*PDEG)` = 2500です。`head(simdata)`と`head(simdata_DEG)`の比較からも、ちゃんとうまく目的の箇所を`FC`倍できていることがわかります。`fc_mat[((G_sim*PDEG*P1-2):(G_sim*PDEG*P1+3)), ]`は、DEG_G1パターンとDEG_G2パターンの切り替わり部分が意図通りに作成できているかを確認しているだけです。最終行で、念のため`write.csv`関数を用いて作成したシミュレーションデータを[JSLAB21_w7.5.csv](JSLAB21_w7.5.csv)という名前で保存しています。

```{r, message=FALSE, error=TRUE, warning=FALSE}
G_sim <- 10000
nk <- c(3, 3)
PDEG <- 0.25
P1 <- 0.9
FC <- 4
K <- length(nk)
L <- as.vector(mapply(rep, 1:K, nk))
head(simdata)
fc_mat <- matrix(1, nrow=G_sim, ncol=sum(nk))
fc_mat[1:(G_sim*PDEG*P1), L==1] <- FC
fc_mat[((G_sim*PDEG*P1+1):(G_sim*PDEG)), L==2] <- FC
fc_mat[((G_sim*PDEG*P1-2):(G_sim*PDEG*P1+3)), ]
simdata_DEG <- simdata*fc_mat
head(simdata_DEG)
write.csv(x=simdata_DEG, file="JSLAB21_w7.5.csv")
```

## W7.6 ヒストグラム1（DEGあり）

W7.5で作成した`simdata_DEG`とラベル情報`L`が存在するという前提です。プロットの-2付近のやや大きめのピークがG1群で高発現のDEG_G1由来のもの、そして+2付近の小さめのピークがG2群で高発現のDEG_G2由来のものです。最後の2行は、作成したシミュレーションデータの右側に`obj`と`logratio`を追加したもの[JSLAB21_w7.6.csv](JSLAB21_w7.6.csv)という名前で保存しています。

```{r, message=FALSE, error=TRUE, warning=FALSE}
hoge <- simdata_DEG
M_G1 <- apply(X=hoge[, L==1], MARGIN=1, FUN=mean)
M_G2 <- apply(X=hoge[, L==2], MARGIN=1, FUN=mean)
logratio <- log2(M_G2/M_G1)
summary(logratio)

obj <- is.finite(logratio)
summary(logratio[obj])
logratio[!obj]
length(logratio[!obj])

d <- data.frame(log2ratio = logratio)
library(ggplot2)
axis_range <- c(-4, 4)
g <- ggplot(d, aes(x=log2ratio)) +
     geom_histogram() + 
     coord_cartesian(xlim=axis_range)
plot(g)
tmp <- cbind(simdata_DEG, obj, logratio)
write.csv(x=tmp, file="JSLAB21_w7.6.csv")
```

## W7.7 ヒストグラム2（DEGあり）

W7.6と基本的に同じですが、1～2250行目のG1群で高発現のDEG_G1、2251～2500行目のG2群で高発現のDEG_G2、そして残りのnonDEGパターンの3種類に分けてプロットするやり方です。`facet_grid`は、[第20回](https://www.iu.a.u-tokyo.ac.jp/kadota/r_seq2.html#book_JSLAB_20)のW9.2.3でも利用しています。

```{r, message=FALSE, error=TRUE, warning=FALSE}
G_sim <- 10000
PDEG <- 0.25
P1 <- 0.9
exp_pat <- rep("nonDEG", G_sim)
exp_pat[1:(G_sim*PDEG*P1)] <- "DEG_G1"
exp_pat[(G_sim*PDEG*P1+1):(G_sim*PDEG)] <- "DEG_G2"

d <- data.frame(log2ratio = logratio,
                Pattern = exp_pat)
library(ggplot2)
axis_range <- c(-4, 4)
g <- ggplot(d, aes(x=log2ratio, color=Pattern)) +
     geom_histogram() + 
     coord_cartesian(xlim=axis_range) +
     facet_grid(rows=vars(Pattern))
plot(g)
```

<br> ちなみに、上のスクリプトは`ggplot(d, aes(x=log2ratio, color=Pattern))`としましたが、これを以下で示すように`ggplot(d, aes(x=log2ratio, fill=Pattern))`のようにすると、外枠だけでなく中も塗りつぶすことができます。

```{r, message=FALSE, error=TRUE, warning=FALSE}
G_sim <- 10000
PDEG <- 0.25
P1 <- 0.9
exp_pat <- rep("nonDEG", G_sim)
exp_pat[1:(G_sim*PDEG*P1)] <- "DEG_G1"
exp_pat[(G_sim*PDEG*P1+1):(G_sim*PDEG)] <- "DEG_G2"

d <- data.frame(log2ratio = logratio,
                Pattern = exp_pat)
library(ggplot2)
axis_range <- c(-4, 4)
g <- ggplot(d, aes(x=log2ratio, fill=Pattern)) +
     geom_histogram() + 
     coord_cartesian(xlim=axis_range) +
     facet_grid(rows=vars(Pattern))
plot(g)
```

## W7.8 パターンごとのlog2(G2/G1)分布

"DEG_G1", "DEG_G2", "nonDEG"の3パターンを分けてlog比の分布を`summary`関数で得ます。`summary(logratio[obj])`までは、W7.3をテンプレートとして作成し、全遺伝子の結果を示しています。残りは、W7.7で作成した"どの行番号に相当する要素がどのパターンに属するかを定めた`exp_pat`"を利用して、特定のパターンに限定して`summary`関数を実行しています。"DEG_G1"のMEANが-2.0、"DEG_G2"のMEANが2.0、そして"nonDEG"のMEANが0.0と酷似した値になっており妥当だと判断できます。

```{r, message=FALSE, error=TRUE, warning=FALSE}
hoge <- simdata_DEG
M_G1 <- apply(X=hoge[, L==1], MARGIN=1, FUN=mean)
M_G2 <- apply(X=hoge[, L==2], MARGIN=1, FUN=mean)
logratio <- log2(M_G2/M_G1)
obj <- is.finite(logratio)
summary(logratio[obj])

summary(logratio[(exp_pat == "DEG_G1") & obj])
summary(logratio[(exp_pat == "DEG_G2") & obj])
summary(logratio[(exp_pat == "nonDEG") & obj])
```

## W7.9 検証

ここでは、一見変な結果のものを見つけ出して確認しています。まず、"DEG_G1"のlog比中で最大値(ただし実数のみ)を示す遺伝子が"gene769"であることを見出し、そのカウント情報を表示する作業を行っています。そして目視で各群のカウントの総和を算出し（本当は平均ですが、数学的には等価です）、手計算でlog比の結果を得て、`summary(logratio[(exp_pat == "DEG_G1") & obj])`の結果と同じ値になっているところまでを確認して安心しています。次に、"DEG_G2"のlog比中で最小値(ただし実数のみ)を示す遺伝子が"gene2498"であり、そのカウント情報を表示、そしてlog比の検証を行っています。話がややこしくなるので気になったヒトのみでよいですが、`which.max(logratio[(exp_pat == "DEG_G1") & obj])`の主な結果である766と"gene769"の間に数値上の不一致があると思われたのではないでしょうか。そしてその疑念は、`which.min(logratio[(exp_pat == "DEG_G2") & obj])`の主な結果である246と"gene2498"の大きな不一致によって確信に変わったと思います。これは、例えば`which.max`は最大値をとるインデックス情報が返されるのですが、入力の段階でlogratio[(exp_pat == "DEG_G1")に限定されていて、かつlog比が実数のもの(`obj`がTRUEの要素のみ)に限定された状態でのインデックス番号なのだと理解すれば納得できると思います。つまり、exp_pat == "DEG_G2")を入力として与えた段階ですでに入力要素数は250個に限定されています。そしてさらにlog比が実数のもの(`obj`がTRUEの要素のみ)に限定されているので、遺伝子名が"gene2498"だったとしても、出力結果は246番目の要素だという風になるのです。

```{r, message=FALSE, error=TRUE, warning=FALSE}
which.max(logratio[(exp_pat == "DEG_G1") & obj])
simdata_DEG["gene769",]
log2(11/4)

which.min(logratio[(exp_pat == "DEG_G2") & obj])
simdata_DEG["gene2498",]
log2(4/14)
```

# tidyverseの概要（W8）

-   [Wickham H, J Stat Software, 2014](10.18637/jss.v059.i10)\
-   [tidyverse](https://www.tidyverse.org/)
    -   [readr](https://CRAN.R-project.org/package=readr)\
    -   [tidyr](https://cran.r-project.org/package=tidyr)\
    -   [dplyr](https://cran.r-project.org/package=dplyr)\
    -   [ggplot2](https://cran.r-project.org/package=ggplot2)\
    -   [tidymodels](https://cran.r-project.org/package=tidymodels)
-   [乳酸菌NGS連載第19回のPDF](https://www.iu.a.u-tokyo.ac.jp/kadota/JSLAB_19_kadota.pdf)\
-   [乳酸菌NGS連載第18回のPDF](https://www.iu.a.u-tokyo.ac.jp/kadota/JSLAB_18_kadota.pdf)\
-   [Python](https://www.python.org/)の[pandas](https://pandas.pydata.org/)\
-   [Qiita](https://qiita.com/)の記事数
    -   [pandasがタグ付けされた記事](https://qiita.com/tags/pandas)\
    -   [tidyverseがタグ付けされた記事](https://qiita.com/tags/tidyverse)

# tidyverseによるパラメータの取得 (W9)

## W9.1 パッケージの読み込み

tidyverseは複数のパッケージのコレクションであり、`library(tidyverse)`の一発で構成するパッケージの提供するすべての関数を利用できるようになります。

```{r message=FALSE, warning=FALSE}
library(tidyverse)
```

## W9.2 csvファイルの読み込み

tidyverseではデータの読み込みや処理に際して、基本的にtibbleと呼ばれる形式を使います。一般的にRでデータを扱う際にはdata.frame形式が使われますが、tibbleはその拡張版と考えて問題ありません。さらにtidyverseの関数はdata.frame型もサポートしているため、tidyverseによる解析において必ずしもtibble形式を使わなければならない訳ではありません（data.frame形式のデータを可視化する際にggplot2を使用する場合が多いですが、ggplot2はもともとtidyverseのパッケージの1つです）。

このtibbleは"tidy data"というデータ構造の考え方に基づいて設計されています。tidy dataは解析を始めるにあたり最適な、"整然とした"データのあり方を定義しており、乱雑なデータから最終的にtidy dataへと変換することがtidyverseにおける前処理の目的と言えます。tidy dataの詳細については[tibbleのリファレンス](https://tidyr.tidyverse.org/articles/tidy-data.html)を参照してください。なお、tidy dataの定義からしてtibbleでは基本的に行名を持つことはできませんが、どうしても行名を指定したい場合`column_to_rownames`関数を使って"ふつうの列"から"行名を表す列"に変換する必要があります。一方、列名は最初の1行目が自動的に指定されます。

```{r message=FALSE, warning=FALSE}
library(tidyverse)

source_url <- "https://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB21_func.R"
input_file_url <- "https://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB18.csv"

source(source_url)
A <- read_csv(input_file_url) %>% 
  column_to_rownames(., var = names(.)[1])
```

## W9.3 パイプ演算子 "%\>%"

`%>%`はパイプ演算子と呼ばれ、データに対して複数の処理を連鎖的に行っていることを視覚的にわかりやすく記述するために導入されています。例えば、次の2つのコードは等価です。

```{r, eval=FALSE}
result <- data %>%
  function1() %>%
  function2() %>%
  function3()
```

```{r, eval=FALSE}
temp1 <- function1(data)
temp2 <- function2(temp1)
result <- function3(temp2)
```

上記の例において%\>%は次のような役割を担っています。

1.  %\>%の左側のオブジェクト(例えば、`data`)を右側の関数(例えば、`function1`)に渡す

2.  関数の戻り値(ここでは`temp1`に相当するもの)を右側の関数(例えば、`function2`)に渡す

3.  %\>%が無い時点で処理の流れを止め、結果をresultに格納する

具体的には、`%>%`は左側のオブジェクトや関数の戻り値を右側の関数の第1引数に渡しています。このとき、受け渡されるオブジェクトは省略することができますが、`.`として明示的に書くことも可能です。したがって、上のコードは次とも等価です。`%>%`に慣れない内はこのように書くのも手でしょう。

```{r, eval=FALSE}
result <- data %>%
  function1(.) %>%
  function2(.) %>%
  function3(.)
```

## W9.4 複数の列に関数を適用する

複数の列に関数を適用するには、`mutate(across())`という関数を利用すると良いでしょう。`mutate`関数は本来、新しい列を作成したり既存の列を変更したりする際に使われる関数であり、単独でも頻繁に使用されますが`across`関数と組み合わせることで、同じ処理を複数の列に対して一括に適用することもできます。次のコードはすべての列に対して正規化を施し、`Anorm`というオブジェクトに格納するものです。

```{r, eval=FALSE}
Anorm <- A %>%
  mutate(across(everything(), ~ 1000000 * . / sum(.))) 
```

ここでは処置対象の列を指定する際に、`everything`関数を使用して"Aのすべての列を対象とする"と指示しています。その他に列を選択する際便利な補助関数として、"◯◯という文字列から始まる列"を選択する`starts_with`関数や、"◯◯という文字列を含む列"を選択する`contains`関数など多くの関数が用意されています。

処理内容を記述する第2引数の部分では、`~`という記号が使用されていますが、これは単に"後続の処理を適用してほしい"と意味しています。これは無名関数（他の言語では"ラムダ式"とも呼ばれます）を用いた書き方であり、簡単な処理に対して関数を別途定義する手間を省くためのテクニックとして割と頻繁に使用されます。少し冗長にはなりますが上記のコードは次の書き方と同じ意味です。

```{r}
norm_func <- function(x) {
  return(1000000 * x / sum(x))
}

Anorm <- A %>% 
  mutate(across(everything(), norm_func))

#' ---気になる人向け---
#' パイプ演算子内の関数の引数に関数をわたしたい場合、
#' 関数オブジェクトを渡さなければならなりません。
#' 下記のnorm_funcに()がついていないのはそのためです。
#' これは主にtidyverseのパイプ演算子をサポートする
#' パッケージ"magrittr"の仕様によるものです。
```

## W9.5 データの"縦持ち"と"横持ち"とは？

同じ情報を持つデータでも、その表現方法には"横持ち"と"縦持ち"の2つの形式があります。一般的に、横持ちのデータは人が見てわかりやすい形式であり、縦持ちのデータは計算機による処理に適した形式とされています。それでは、まずはざっと目でこれらの形式を確認してみましょう。

次に、おなじみの横持ちで表されたカウントデータを見てみます。

```{r}
head(A)
```

一方、次は縦持ちに変換したデータです。

```{r}
A %>%
  rownames_to_column(var = "gene_ID") %>%
  pivot_longer(cols = -gene_ID, names_to = "sample", values_to = "count") %>%
  separate(sample, into = c("pH", "time", "replicate"), sep = "_") %>% 
  unite("pH_time", c("pH", "time"), sep = "_") %>% 
  head(A, n = 10)
```

これらのデータを比較してみると、縦持ちのデータでは、列に含まれていた"pH4.5"や"pH7"という情報が、抽象化された"pH"という列名の中に格納されていることがわかります。同様に、列名の"1h"、"24h"や"CCG"という情報は"time"という列に、各反復に付けられた番号"1〜3"は"replicate"という列にまとめて格納されています。つまり、縦持ちのデータでは、1行にあたり1つの観測値（count）と観測条件（gene_ID、pH_time、replicate）が含まれる形式を持つようになります。

縦持ちに変換することによって処理の高速化、可視化の容易さなど、計算機目線のメリットも多く存在します。しかし殊解析においては、**列名に隠れていた観測条件をデータの値として抽出し、その値に基づいて計算機が処理を変えることができるようになる**という点が大きなメリットと言えるでしょう。

今回の場合、"酸ストレス短期暴露群（pH4.5_1h）"、"酸ストレス長期暴露群（pH4.5_24h）"、"対照群（pH7_CCG）の各3回の反復において平均、分散、dispersionの3つの値を算出することがゴールです。このような時には`pivot_longer`関数と`group_by`関数の組み合わせが非常に有効です。具体的なやり方を次のW9.6で見ていきましょう。

## W9.6 データを縦持ちに変換し、グループごとに処理を適用する

再度横持ちのデータから出発し、縦持ちへの変換とグループごとの処理に用いる関数を一つずつ見ていきましょう。

```{r}
head(Anorm)
```

列名を観察すると、"(pHの値)\_(暴露した時間)\_(反復数)"という形式で統一されていることがわかります。例えば最初のpHについて注目すると、"pH4.5\_・・・"と"pH7\_・・・"の2種類存在します。よって、これらをまとめて"pH"という名前の列を作って、その中の値として"pH4.5"と"pH7"のいずれかが格納されているようにするのが変換のゴールです。同様に、2番目の暴露した時間は"time"という列に、3番目の反復については"replicate"という列にまとめることにしましょう。また、肝心の観測値は"count"という列名に格納することにします。

縦持ちデータへの変換には`pivot_longer`関数という便利なモジュールが用意されており、このように"\_"などのデリミタで分けられた列名に対して特に柔軟な操作を可能にしています。`cols = -gene_ID`は"gene_IDという列を除く全ての列"という意味であることを踏まえて、次のようにデータを縦持ちに変数することができます。実際に縦持ちに変換できているか確認しましょう。

```{r}
Anorm %>% 
  rownames_to_column(var = "gene_ID") %>%
  pivot_longer(cols = -gene_ID, names_to = c("pH", "time", "replicate"), values_to = "count", names_sep = "_") %>% 
  head()
```

無事"pH"、"time"、"replicate"という名前の列を作り、条件を値になおすことができたようです。このまま次の解析に進んでも良いですが、3つの群の条件は"pH + time"と簡潔に表されるため、今のままではやや冗長に思えます。例えばG1は"pH4.5_1h"、G2は"pH4.5_24h"のよう表した方がわかりやすく、今後の解析も楽に進められそうです。`unite`関数で"pH_time"という名前の列を作成し、"pH"と"time"の列を一つの変数としてまとめそこに格納してみよましょう。

```{r}
Anorm %>% 
  rownames_to_column(var = "gene_ID") %>%
  pivot_longer(cols = -gene_ID, names_to = c("pH", "time", "replicate"), values_to = "count", names_sep = "_") %>% 
  unite(col = "pH_time", c("pH", "time"), sep = "_") %>% 
  head()
```

以上で縦持ちのtidyなデータを用意することができました。次に各遺伝子の各群ごとで求めたい統計量を算出していきましょう。

一旦縦持ちになおしてしまえば、次のようにパイプ演算子内で処理したい部分を`group_by`関数と`ungroup`関数で挟むだけで自動的にグループ分けして処理してくれます。今回は遺伝子のグループ(`gene_ID`列)と条件のグループ(`pH_time`列)ごとに処理したいので、`group_by`関数の引数にそれらの列名を指定します。次のコード内の"\# 処理したい内容"は、`gene_ID`の値と`pH_time`の値ごとにに処理されます。つまり`replicate`の1〜3に対して各統計量が計算されるようになります。

```{r, eval=FALSE}
Anorm %>% 
  rownames_to_column(var = "gene_ID") %>%
  pivot_longer(cols = -gene_ID, names_to = c("pH", "time", "replicate"), values_to = "count", names_sep = "_") %>% 
  unite(col = "pH_time", c("pH", "time"), sep = "_") %>% 
  group_by(gene_ID, pH_time)
  # 処理したい内容
  ungroup()
  head()
```

## W9.7 列やグループごと情報を1つの値に要約する

ある列やグループにはたいてい複数の値が入っています。これらの値を処理して、平均値といった1つの値に集約したい時は`summarise`関数が有効です。"複数の値→1つの値"という処理がこの関数の大きな役目ですから、BaseRに標準装備されている関数の他に自前の関数を用いて集約させることも可能です。次のコードは各グループごとの平均`M`、分散`V`、dispersionの`Phi`を算出するコードですが、`Phi = func_phi()` のように、先に読み込んだ関数を利用して処理を行っていることに注目してください。

```{r}
Anorm %>%
  rownames_to_column(var = "gene_ID") %>%
  pivot_longer(cols = -gene_ID, names_to = c("pH", "time", "replicate"), values_to = "count", names_sep = "_") %>% 
  unite(col = "pH_time", c("pH", "time"), sep = "_") %>%
  group_by(gene_ID, pH_time) %>%
  summarize(
    M = mean(count),
    V = var(count),
    Phi = func_phi(count)
  ) %>% 
  ungroup()
```

## W9.8 データを横持ちになおす

縦持ち形式のデータが計算機による処理に適していることは先に見てきましたが、一方でヒトにとっては読みにくい形式と言えます。今後の解析をスムーズにすすめるためにも、最後にデータを横持ちになおすことにしましょう。`pivot_wider`関数を用います。`pivot_wider`関数は`names_ftom`で指定された列に格納された値を列名へと変換する機能を提供しています。`names_sep`引数で指定されたデリミタを用いて`values_from`引数で指定された列名を結合して新たな列を作成し、`values_from`引数の各列の値をそれらの列に格納します。

```{r}
Anorm %>%
  rownames_to_column(var = "gene_ID") %>%
  pivot_longer(cols = -gene_ID, names_to = c("pH", "time", "replicate"), values_to = "count", names_sep = "_") %>% 
  unite(col = "pH_time", c("pH", "time"), sep = "_") %>%
  group_by(gene_ID, pH_time) %>%
  summarize(
    M = mean(count),
    V = var(count),
    Phi = func_phi(count)
  ) %>% 
  ungroup() %>% 
  pivot_wider(names_from = pH_time, values_from = c("M", "V","Phi"), names_sep = "_")
```

## W9.9 列の名前を変更する

現在、各実験条件は"(pHの値)\_(暴露した時間)"という形式で表されていますが、これをさらに簡潔に記述できるように"G1"、"G2"、"G3"という形に列名を変更しておきましょう。この処理のメリットは、もちろんデータの視認性の向上もありますが、tidyverseでは「末尾がG1で終わる列のMとPhiの値を〇〇して・・・」のように列名を指定して解析を行う場合が多いため、コードを書く際の手間を減らしミスを少なくする、という面でも大きいと言えます。

列名を変更する最もシンプルな方法は`rename`関数を用いるものです。しかしこの方法で一つ一つ名前を変えていくのは手間ですし、ミスの温床になり得ます。今回のように特定の条件を満たす列名を一括で変更したいときは、`rename_with`関数を用いるのが有効です。以下のコードでは、末尾が"pH4.5_1h"、"\_pH4.5_24h"、"\_pH7_CCG"で終わる列をそれぞれ一括で選択し、Base-Rに含まれる`gsub`関数により正規表現を用いて変更を行っています。

```{r}
Anorm %>%
  rownames_to_column(var = "gene_ID") %>%
  pivot_longer(cols = -gene_ID, names_to = c("pH", "time", "replicate"), values_to = "count", names_sep = "_") %>% 
  unite(col = "pH_time", c("pH", "time"), sep = "_") %>%
  group_by(gene_ID, pH_time) %>%
  summarize(
    M = mean(count),
    V = var(count),
    Phi = func_phi(count)
  ) %>% 
  ungroup() %>% 
  pivot_wider(names_from = pH_time, values_from = c("M", "V","Phi"), names_sep = "_") %>% 
  rename_with(~gsub("_pH4.5_1h$", "_G1", .), ends_with("_pH4.5_1h")) %>% 
  rename_with(~gsub("_pH4.5_24h$", "_G2", .), ends_with("_pH4.5_24h")) %>% 
  rename_with(~gsub("_pH7_CCG$", "_G3", .), ends_with("_pH7_CCG")) 
```

上のコードで一旦`~`を用いて無名関数を置いて`gsub`関数を呼び出している点が気になった人がいるかもしれません。W9.4でも同様な書き方をしていますが気になる人向けに説明すると、この理由はパイプ演算子内の関数の引数に関数を渡したいという時、関数オブジェクトの形で呼び出す必要があるのですが(つまり、`gsub`のみで呼び出すのが記法としては正しい)、それでは`gsub("_pH4.5_1h$", "_G1", .)`のように**"どのように文字列の変更を行うか"を引数で指定することができません**。引数をとれる`gsub(…)`の形で関数を呼び出すために、無名関数と`.`を用いて`rename_with`関数に渡している、というわけです。

## W9.10 まとめ

以上の説明を踏まえ、W2で示したコードをtidyverse流で記述した最終的なコードは次のようになります。

```{r message=FALSE, warning=FALSE}
library(tidyverse)

source_url <- "https://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB21_func.R"
input_file_url <- "https://www.iu.a.u-tokyo.ac.jp/kadota/book/JSLAB18.csv"

source(source_url)
A <- read_csv(input_file_url) %>% 
  column_to_rownames(., var = names(.)[1])

# 正規化
Anorm <- A %>%
  mutate(across(everything(), ~ 1000000 * . / sum(.)))

# 各群の統計量を算出
result <- Anorm %>%
  rownames_to_column(var = "gene_ID") %>%
  pivot_longer(cols = -gene_ID, names_to = c("pH", "time", "replicate"), values_to = "count", names_sep = "_") %>% 
  unite(col = "pH_time", c("pH", "time"), sep = "_") %>%
  group_by(gene_ID, pH_time) %>%
  summarize(
    M = mean(count),
    V = var(count),
    Phi = func_phi(count)
  ) %>%
  ungroup() %>%
  pivot_wider(names_from = pH_time, values_from = c(M, V, Phi), names_sep = "_") %>%
  rename_with(~gsub("_pH4.5_1h$", "_G1", .), ends_with("_pH4.5_1h")) %>% 
  rename_with(~gsub("_pH4.5_24h$", "_G2", .), ends_with("_pH4.5_24h")) %>% 
  rename_with(~gsub("_pH7_CCG$", "_G3", .), ends_with("_pH7_CCG"))
```

# G1群のみの$\varPhi$ \> 0を満たす経験分布 (W10)

## W10.1 列の選択

列を選択して抽出するためには、select関数を使用します。コードの中で単に列名を指定することも可能ですが、先に登場した`starts_with`関数や`ends_with`関数などの補助関数と組み合わせて使用することもできます。また、特定の列を除外する場合は、列名や列を取得する補助関数の前に`-`を付けるだけでその列を除外した結果を得ることもできます。例えば次のコードは、`M_G1`と`Phi_G1`を除いたすべての列を表示しています。

```{r}
result %>% 
  select(-M_G1, -Phi_G1)
```

## W10.2 行の選択

特定の条件を満たす行の選択と抽出には、filter関数を使用します。複雑な条件を指定したい場合、`<`(より小さい)や、`==`(等しい)などの比較演算子や、OR演算子`|`、AND演算子`&`といった論理演算子を組み合わせて表現できます。

select関数とfilter関数は、tidyverseを用いてデータを操作する際に頻繁に使用します。例えば遺伝子名と各群における統計量が格納されているresultから末尾が"G2"で終わる列を抽出し、その中でPhiが0より大きく、平均が50以下という条件を満たす遺伝子名を抽出するには、次のように記述します。

```{r}
result %>%
  select(gene_ID, ends_with("G2")) %>% 
  filter(Phi_G2 > 0 & M_G2 <= 50) %>% 
  select(gene_ID)
```

## W10.3 まとめ

以上の説明を踏まえ、W4.2で示したコードをtidyverse流で記述した最終的なコードは次のようになります。

```{r}
sub_G1_M_Phi <- result %>% 
  select(M_G1, Phi_G1) %>% 
  filter(Phi_G1 > 0)
head(sub_G1_M_Phi)
```

# 乱数を生成 (W11)

## W11.1 データからサンプリングする

データから行をサンプリングするには、`slice_sample`関数を使用します。BaseRの`sample`関数と同様に、`replace`引数を`TRUE`に指定することで復元抽出も可能です。

## W11.2 1行ずつ処理を適用する

先に見てきた`group_by`関数は、列を指定することでその値を基にグループ分けして処理を行ってくれるものでした。それに対して、`rowwise`関数を適用すると、以降の一つ一つの行を個別のグループとして処理を行ってくれます。使い方も`group_by`関数と同様に、行ごとに行ってほしい処理を`rowwise()`と`ungroup()`で挟むだけです。

`rowwise`関数を適用するとしないとでは得られる結果がかなり異なってしまいます。`rowwise`関数の役割を理解するために、`rowwise`関数を使用しなかった場合のコードの挙動を見ていきましょう。

```{r}
G_sim <- 4
n_sim <- 5
set.seed(2023)
sub_G1_M_Phi %>% 
  slice_sample(n = G_sim, replace = TRUE) %>% 
  mutate(rep = list(rnbinom(n = n_sim, mu = .[["M_G1"]], size = 1 / .[["Phi_G1"]]))) %>%
  unnest_wider(col = rep, names_sep = "_") %>%
  mutate(gene_ID = paste0("gene", 1:nrow(.))) %>%
  column_to_rownames("gene_ID") %>% 
  head()
```

上記のように、rep_1からrep_5の各列で延々と同じ値が繰り返されるという結果になりました。これはどういうことでしょうか。実は3行目の`rnbinom`関数の仕様による問題が起こっています。

上記のコードでは、`rnbinom`関数の`mu`と`size`引数に要素数4(= `G_sim`)の"M_G1"と"Phi_G1"が渡されますが、生成される乱数の個数はあくまで引数に指定したn個のみのようです。以下の例から、引数に要素数5のベクトルを渡しても整数を一つ渡しても、生成される乱数の個数は、`n`引数で指定した同じ3つであることを確認してください。

```{r}
set.seed(2023)
rnbinom(n = 3, mu =5, size = 8)
set.seed(2023)
# muとsizeに要素数5のベクトルを渡してみる
rnbinom(n = 3, mu = c(5, 3, 4, 9, 3), size = c(8, 2, 5, 1, 4))
```

このように`rowwise`関数を適用しない場合、`rnbinom`関数の引数にベクトルが渡され、関数による処理が1回だけ行われることになります。しかしその結果生成されるのは、あくまで要素数n個のベクトルです。そのベクトルが各行の`rep`という列に一様にコピーされることが原因で同じ値が生成されたように見えるというわけです。

`rnbinom`関数の引数に渡したいものは、あくまで各行における"M_G1"や"Phi_G1"といった列の**ひとつの**値です。列をベクトルとしてまるごと引数に渡すのではなく、各行ごとに列の値を取り出して処理に渡したい場合に、`rowwise`関数が有効であることがわかるでしょう。

## W11.3 ネストされた値を行方向に展開する

data.frame内の値にリストやベクトルを格納できる様に、tibbleの値にもリストやベクトルを格納することができます。以下の例では、各行のMとPhiをrnbinom関数に渡して得た乱数のベクトルがrepという列に格納されています。

```{r}
G_sim <- 4
n_sim <- 5
set.seed(2023)
sub_G1_M_Phi %>% 
  slice_sample(n = G_sim, replace = TRUE) %>% 
  rowwise() %>%
  mutate(rep = list(rnbinom(n = n_sim, mu = .[["M_G1"]], size = 1 / .[["Phi_G1"]]))) %>%
  ungroup()
```

これらのベクトルを横方向に展開するために、`unnest_wider`関数を用いることができます。上の例ではひとつの値に格納されていた要素数5のベクトルが、下の例で行方向に一つずつ展開されていることを確認してください。

```{r}
G_sim <- 4
n_sim <- 5
set.seed(2023)
sub_G1_M_Phi %>% 
  slice_sample(n = G_sim, replace = TRUE) %>% 
  rowwise() %>%
  mutate(rep = list(rnbinom(n = n_sim, mu = .[["M_G1"]], size = 1 / .[["Phi_G1"]]))) %>%
  ungroup() %>% 
  unnest_wider(col = rep, names_sep = "_")
```

## W11.4 まとめ

以上の説明を踏まえ、W4.6で示したコードをtidyverse流で記述した最終的なコードは次のようになります。

```{r}
G_sim <- 4
n_sim <- 5
set.seed(2023)
simdata <- sub_G1_M_Phi %>% 
  slice_sample(n = G_sim, replace = TRUE) %>% 
  rowwise() %>%
  mutate(rep = list(rnbinom(n = n_sim, mu = .[["M_G1"]], size = 1 / .[["Phi_G1"]]))) %>%
  ungroup() %>%
  unnest_wider(col = rep, names_sep = "_") %>%
  mutate(gene_ID = paste0("gene", 1:nrow(.))) %>%
  column_to_rownames("gene_ID") %>% 
  select(-"M_G1", -"Phi_G1")
simdata
```

# 遺伝子数10000のデータ生成 (W12)

W4.10に示したコードをtidyverse流に書くと次のようになります。G1群の平均とdispersionが格納されたオブジェクト`sub_G1_M_Phi`から10000行復元抽出し、それらをパラメータとして乱数を5個生成し行方向に展開することですべてがnonDEGである場合のシミュレーションデータとしています。実行した結果から、W4.10と同様に10000行(=遺伝子数)、5列(=反復数)のデータが生成されていることを確認してください。

```{r}
G_sim <- 10000
n_sim <- 5

set.seed(2023)

simdata <- sub_G1_M_Phi %>% 
  slice_sample(n = G_sim, replace = TRUE) %>% 
  rowwise() %>%
  mutate(rep = list(rnbinom(n = n_sim, mu = .[["M_G1"]], size = 1 / .[["Phi_G1"]]))) %>%
  unnest_wider(col = rep, names_sep = "_") %>%
  mutate(gene_ID = paste0("gene", 1:nrow(.))) %>%
  column_to_rownames("gene_ID") %>% 
  select(-"M_G1", -"Phi_G1")
dim(simdata)
```

# DEGありのデータ生成 (W13)

## W13.1 $\varPhi$ \> 0を満たす経験分布を得る

W7.1に示したコードをtidyverse流に書くと次のようになります。G1〜G3のすべての群から、乱数生成のパラメータに使用するMとPhiをデータから抽出し、縦持ちに変換してPhiが0より大きいものを抽出、そして縦持ちのまま`sub_M_Phi`に格納しています。`sub_M_Phi`の行数がW7.1の結果と同じであることを確認してください。このような処理において、tidyverseの可読性の高さや縦持ちの便利さが際立つと言えるでしょう。

```{r}
sub_M_Phi <- result %>% 
  select(-starts_with("V")) %>% 
  pivot_longer(cols = c(starts_with("M"), starts_with("Phi")), names_to = c(".value", "group"), names_sep = "_") %>% 
  filter(Phi > 0)

nrow(sub_M_Phi)
```

## W13.2 全てがnonDEGの数値行列を生成

W7.2に示したコードをtidyverse流に書くと次のようになります。W13.1で得た`sub_G1_M_Phi`から10000行を反復抽出を許してサンプリングし、それらの値をパラメータとし5つの乱数を生成することで全てnonDEGのシミュレーションデータとしています。

```{r}
G_sim <- 10000
nk <- c(3, 3)
set.seed(2023)
simdata <- sub_G1_M_Phi %>% 
  slice_sample(n = G_sim, replace = TRUE) %>% 
  rowwise() %>%
  mutate(G1 = list(rnbinom(n = nk[1], mu = .data[["M_G1"]], size = 1 / .data[["Phi_G1"]])),
         G2 = list(rnbinom(n = nk[2], mu = .data[["M_G1"]], size = 1 / .data[["Phi_G1"]]))) %>%
  unnest_wider(col = G1, names_sep = "_") %>%
  unnest_wider(col = G2, names_sep = "_") %>%
  select(-"M_G1", -"Phi_G1")
head(simdata)
```

## W13.3 DEGの導入

W7.5に示したコードをtidyverse流に書くと次のようになります。`mutate`関数を用いて、nonDEGで生成したシミュレーションデータ`simdata`から、G1で高発現するDEGを1〜2250行目、G2で高発現するDEGを2251〜2500行目として導入しています。

```{r}
G_sim <- 10000
nk <- c(3, 3)
PDEG <- 0.25
P1 <- 0.9
FC <- 4
```

```{r}
simdata_DEG <- simdata %>% 
  mutate(across(starts_with("G1"), ~ ifelse(row_number() <= G_sim * PDEG * P1, . * FC, .))) %>%
  mutate(across(starts_with("G2"), ~ ifelse(row_number() > G_sim * PDEG * P1 & row_number() <= G_sim * PDEG, . * FC, .)))

write.csv(x=simdata_DEG, file="JSLAB21_w7.5_tidyverse.csv")
```

## W13.4 ヒストグラムによるDEGありシミュレーションデータの確認

W7.6に示したコードをtidyverse流に書くと次のようになります。tidyverseの記法によって生成したDEGありのシミュレーションデータのヒストグラムと、W7.6のヒストグラムが同様の形状を示していることを確認してください。

```{r}
M <- simdata_DEG %>%
  rowwise() %>%
  mutate(
    M_G1 = mean(c_across(starts_with("G1"))),
    M_G2 = mean(c_across(starts_with("G2")))) %>% 
  ungroup() %>% 
  select(starts_with("M")) %>% 
  mutate(logratio = log2(M_G2/M_G1))
```

```{r message=FALSE, warning=FALSE}
axis_range <- c(-4, 4)
g <- ggplot(M, aes(x=logratio)) +
     geom_histogram() + 
     coord_cartesian(xlim=axis_range)
plot(g)
```
