はじめに（W1）

パラメータの取得（W2）

入力ファイル(JSLAB18.csv)の読込から、群ごとのカウントの平均(\(M\))とdispersion (\(\varPhi\))を得る必要最小限のスクリプトです。

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のパラメータを用いて説明しようと決めました。

head(M, n=4)
##                       [,1]       [,2]        [,3]
## EBG00001128470   219.35372   139.9754   163.20042
## EBG00001128476    88.02994   116.9911    45.61295
## EBG00001128500 26661.18736 23549.0411 35571.20466
## EBG00001128509   259.62436   166.7145   190.11327
head(V, n=4)
##                       [,1]         [,2]         [,3]
## EBG00001128470   13304.976     246.7638 4.272069e+03
## EBG00001128476    2165.165     178.4106 1.521660e+01
## EBG00001128500 2736641.821 2342808.8423 5.120825e+07
## EBG00001128509   29046.424     151.8007 2.223958e+03
head(Phi, n=4)
##                       [,1]          [,2]        [,3]
## EBG00001128470 0.271959587  0.0054503107  0.15426937
## EBG00001128476 0.268042701  0.0044874537 -0.01460983
## EBG00001128500 0.003812477  0.0041821822  0.04044278
## EBG00001128509 0.427073383 -0.0005365914  0.05627206

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

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

i <- 4
k <- 3
M[i,k]
## EBG00001128509 
##       190.1133
Phi[i,k]
## EBG00001128509 
##     0.05627206
set.seed(2023)
output <- rnbinom(n=8, mu=M[i,k], size=1/Phi[i,k])
output
## [1] 177 175 153 174 185 206 202 216
var(output)
## [1] 427.4286
V[i,k]
## EBG00001128509 
##       2223.958

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

set.seed(2023)
output <- rnbinom(n=1000000, mu=M[i,k], size=1/Phi[i,k])
var(output)
## [1] 2224.638

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

?rnbinom

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

①のチャンク実行後は、以下に示すような関数マニュアルが②Helpタブ上に見られます。計4つの関数についての説明がなされており、目的の関数は③で見えています。③rnbinomは、n, size, prob, muという計4つの引数（Arguments）を受けつけることがわかりますが、それぞれの説明は④に書かれています。

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

W3.4の続きとして、④Argumentsの⑤sizeに関する説明が見られるようにした状態です。sizeオプションで与える値は、⑥に明記されているように整数（integer）である必要はありませんが、正の値をもたねばなりません。

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

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

i <- 4
k <- 2
M[i,k]
## EBG00001128509 
##       166.7145
Phi[i,k]
## EBG00001128509 
##  -0.0005365914
set.seed(2023)
output <- rnbinom(n=8, mu=M[i,k], size=1/Phi[i,k])
## Warning in rnbinom(n = 8, mu = M[i, k], size = 1/Phi[i, k]): NAs produced
output
## [1] NaN NaN NaN NaN NaN NaN NaN NaN
var(output)
## [1] NA

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

W4.1 第20回の図4

第20回の原稿はモノクロですので、こちらのほうが特に(d)はわかりやすいかと思います。

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になっているだけです。ここのスクリプトは、わかりやすさ（説明のしやすさ）重視で汎用性は低いのでご注意ください。

k <- 1
obj <- (Phi[,k] > 0)
sum(obj)
## [1] 2511
head(obj)
## EBG00001128470 EBG00001128476 EBG00001128500 EBG00001128509 EBG00001128529 
##           TRUE           TRUE           TRUE           TRUE          FALSE 
##      LGG_00001 
##           TRUE
head(Phi[, k])
## EBG00001128470 EBG00001128476 EBG00001128500 EBG00001128509 EBG00001128529 
##    0.271959587    0.268042701    0.003812477    0.427073383    0.000000000 
##      LGG_00001 
##    0.002775528
sub_M <- M[obj, k]
sub_Phi <- Phi[obj, k]

W4.3 経験分布の記載例

Kadota et al., Algorithms Mol Biol, 2012では、以下のような記載をしています。
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 package.
2012年の論文は、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が返されているのをみると全体像がよくわかると思います。

G_sim <- 4
set.seed(2023)
s_num <- sample(1:length(sub_M), G_sim, replace=TRUE)
s_num
## [1] 1909 1488 2479 1385
sub_M[s_num]
## LGG_02238 LGG_01737 LGG_02907 LGG_01623 
## 875.46580  60.89932 274.83665  89.63483
sub_Phi[s_num]
##  LGG_02238  LGG_01737  LGG_02907  LGG_01623 
## 0.22292266 0.16477002 0.11164579 0.04881805
match(x=names(sub_M[s_num]), table=rownames(M))
## [1] 2243 1742 2912 1628
match(x=c("D", "B"), table=c("A", "B", "C", "D", "E"))
## [1] 4 2

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が難解であったヒト用です。

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]])
## [1] 738 705 445 643 733
i <- 2
rnbinom(n=n_sim, mu=sub_M[s_num[i]], size=1/sub_Phi[s_num[i]])
## [1] 67 65 73 92 41
i <- 3
rnbinom(n=n_sim, mu=sub_M[s_num[i]], size=1/sub_Phi[s_num[i]])
## [1] 219 198 343 171 204
i <- 4
rnbinom(n=n_sim, mu=sub_M[s_num[i]], size=1/sub_Phi[s_num[i]])
## [1] 106 100  94 108 113

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になる、といった具合です。

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
##      [,1] [,2] [,3] [,4] [,5]
## hoge  738  705  445  643  733
## hoge   67   65   73   92   41
## hoge  219  198  343  171  204
## hoge  106  100   94  108  113

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列目）の分散情報を表示させています。

rownames(simdata) <- paste("gene", 1:G_sim, sep="")
colnames(simdata) <- paste("rep", 1:n_sim, sep="_")
simdata
##       rep_1 rep_2 rep_3 rep_4 rep_5
## gene1   738   705   445   643   733
## gene2    67    65    73    92    41
## gene3   219   198   343   171   204
## gene4   106   100    94   108   113
apply(X=simdata, MARGIN=1, FUN=var)
##   gene1   gene2   gene3   gene4 
## 14923.2   334.8  4506.5    54.2
posi <- match(x=names(sub_M[s_num]), table=rownames(M))
posi
## [1] 2243 1742 2912 1628
V[posi, 1]
##   LGG_02238   LGG_01737   LGG_02907   LGG_01623 
## 171732.3881    671.9864   8708.0225    481.8588

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

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

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)
##        hoge        hoge        hoge        hoge 
## 171990.2270    670.3646   8688.4658    482.0801
posi <- match(x=names(sub_M[s_num]), table=rownames(M))
V[posi, 1]
##   LGG_02238   LGG_01737   LGG_02907   LGG_01623 
## 171732.3881    671.9864   8708.0225    481.8588

W4.9 復元抽出の意味

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

k <- 1
obj <- (Phi[,k] > 0)
sub_M <- M[obj, k]
sub_Phi <- Phi[obj, k]
length(sub_M)
## [1] 2511
G_sim <- 3000
set.seed(2023)
s_num <- sample(1:length(sub_M), G_sim, replace=FALSE)
## Error in sample.int(length(x), size, replace, prob): cannot take a sample larger than the population when 'replace = FALSE'

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

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

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)
## [1] 10000     5

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

第20回のW9.7に対応するプロットの基本形です。

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

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

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)
## [1] 2949    3
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)
## Warning: package 'ggplot2' was built under R version 4.1.2
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)
## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: Transformation introduced infinite values in continuous y-axis

d1_sim <- d
g1_sim <- g

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

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

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)
## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: Transformation introduced infinite values in continuous y-axis

d1_real <- d
g1_real <- g

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

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

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)
## [1] 2949    3
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)
## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: Transformation introduced infinite values in continuous y-axis

d2_sim <- d
g2_sim <- g

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

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

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)
## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: Transformation introduced infinite values in continuous y-axis

d2_real <- d
g2_real <- g

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

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

library(patchwork)
g1_real + g2_real + g1_sim + g2_sim + plot_layout(ncol=2)

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

第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回のW9.8.2をテンプレートとしてプロットしています。

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)
## [1] 2949    3
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)
## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: Transformation introduced infinite values in continuous y-axis

d1_sim <- d
g1_sim <- g

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

W5.2と第20回のW9.8.2をテンプレートとしてプロットしています。

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)
## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: Transformation introduced infinite values in continuous y-axis

d1_real <- d
g1_real <- g

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

W5.3と第20回のW9.8.2をテンプレートとしてプロットしています。

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)
## [1] 2949    3
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)
## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: Transformation introduced infinite values in continuous y-axis

d2_sim <- d
g2_sim <- g

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

W5.4と第20回のW9.8.2をテンプレートにしています。

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)
## Warning: Transformation introduced infinite values in continuous x-axis
## Warning: Transformation introduced infinite values in continuous y-axis

d2_real <- d
g2_real <- g

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

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

g1_real + g2_real + g1_sim + g2_sim + plot_layout(ncol=2)