・お題：他にも遺伝子発現解析の例があったので、そちらもなぞってみたい。
・次の例は以下。とても詳しく記載してあるので、きちんとしたことは元サイトをご確認いただきたい。
・元論文は以下。
・メスマウスから単離したCD24+CD29hi乳腺幹細胞-濃縮基底細胞（MaSC/basal）、CD24+CD29loCD61+管腔前成細胞-濃縮（LP）及びCD24+CD29loCD61-熟腔-濃縮（ML）細胞集団の遺伝子発現レベルを決定し、さらに、これらの初代細胞タイプと培養MaSC/Basal由来マンモスフィア細胞及びCommaD-βGeo（CommaDβ）細胞株との比較を行ったらしい。
・まず、ライブラリをアクチベートする。
> library(limma)
> library(Glimma)
> library(edgeR)
> library(Mus.musculus)
Error in library(Mus.musculus) :
‘Mus.musculus’ という名前のパッケージはありません
・Mus.musculusというパッケージがないらしい。Bioconductorでマウスの遺伝子発現を解析するときに使うパッケージっぽい。
・とりあえずインストールする。
> if (!require("BiocManager", quietly = TRUE))
+ install.packages("BiocManager")
> BiocManager::install("Mus.musculus")
・これでしばらく待つとインストールできる。もう一度ライブラリをアクチベートする。
> library(Mus.musculus)
・今度は問題なく行けた。今回は、公共データをダウンロードするところから始まる。
・まず、urlというオブジェクトにデータのアドレスを格納。
> url <- "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file"
・このアドレスが何なのかというと、プログラム経由でncbiからデータをとって来るのに使えるアドレスらしい。acc=の後のところを弄ると、他のデータセットもダウンロードできるみたい（引数？の意味は以下のサイトの"Construct a URL"をご確認いただきたい）。
・次に、データをダウンロードする。utilsパッケージのdownload.file関数を使う。
> utils::download.file(url, destfile="GSE63310_RAW.tar", mode="wb")
・すると、GSE63310_RAW.tarというアーカイブがカレントディレクトリにできる。このままだとデータを使えないので、ファイルを取り出す。すると、カレントディレクトリにたくさんのファイルができる。
> utils::untar("GSE63310_RAW.tar", exdir = ".")
・gzファイルは圧縮されたファイル形式なので、解凍する必要がある。Rのutlisにあるgunzipで解凍する。すると、カレントディレクトリに解凍済ファイルがたくさんできる。
> files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt", "GSM1545538_purep53.txt",
+ "GSM1545539_JMS8-2.txt", "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
+ "GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt", "GSM1545545_JMS9-P8c.txt")
> for(i in paste(files, ".gz", sep=""))
+ R.utils::gunzip(i, overwrite=TRUE)
・これらのテキストファイルには、生の遺伝子レベルのカウントが入っている。試しに1つ読んでみる。
> read.delim(files[1], nrow=5)
・1列目に遺伝子のIDが、3列目にカウントデータが入っているらしい。次に、readDGE関数でカウントデータを読み込む。引数の意味に関しては、以下の解説を参照していただきたい。読み込むファイル名、遺伝子が含まれている列番号、カウントデータが含まれている列番号を指定すると、まとめてDGEListクラスのファイルを作ってくれる。
> x <- readDGE(files, columns=c(1,3))
・サンプル名がGSM1545545_****というのがクドいので、簡単にする。colnamesで元のサンプル名を抽出し、頭12文字（GSM1545545_）を削り、samplenamesと名付ける。
> samplenames <- substring(colnames(x), 12, nchar(colnames(x)))
> samplenames
[1] "10_6_5_11" "9_6_5_11" "purep53" "JMS8-2" "JMS8-3" "JMS8-4" "JMS8-5" "JMS9-P7c" "JMS9-P8c"
> colnames(x) <- samplenames#colnamesにsamplenamesを付けなおす。
・このままでは、どのサンプルがなんのサンプルなのか分からない。Nasal、LP、MLのいずれのサンプルなのか、factorで与える。
> group <- as.factor(c("LP", "ML", "Basal", "Basal", "ML", "LP", "Basal", "ML", "LP"))
> x$samples$group <- group
・更に、バッチ情報として、laneを作成し、こちらもfactorとして付け加える。
> lane <- as.factor(rep(c("L004","L006","L008"), c(3,4,2)))
> x$samples$lane <- lane
> x$samples
・次に、どの遺伝子IDがなんの遺伝子なのか、情報を加える。
> geneid <- rownames(x)#もともとrownamesがgeneIDになっていた
> genes <- select(Mus.musculus, keys=geneid, columns=c("SYMBOL", "TXCHROM"), keytype="ENTREZID")#情報はMus.musculusに入っている。geneidをキーとして、ENTREZIDを探し、SYMBOLとTXCHROMをとってきたい。
> head(genes)
・複数の遺伝子IDが同じ遺伝子を指している（重複している）ケースがあるので、処理しておきたい。今回は、染色体の一つを選択して重複するアノテーションを持つ遺伝子を代表させる。
#duplicated(genes$ENTREZID)で重複の有無のBoolを取得し、!で反転して選抜している感じ。
> genes <- genes[!duplicated(genes$ENTREZID),]
> x$genes <- genes
・ここまでやると、データセットが以下のようになる。
> x
・次に前処理に進むわけだけれど、そのままのデータで処理するのではなく、edgeRのlog cpmを使って処理する。
> lcpm <- cpm(x, log=TRUE)
> L <- mean(x$samples$lib.size) * 1e-6#ライブラリサイズ/mil.の平均値
> M <- median(x$samples$lib.size) * 1e-6#ライブラリサイズ/mil.の中央値
> c(L, M)
[1] 45.48855 51.36862
・lcpmがどんな分布になっているか見てみる。せっかくなので、tidyverseを使って作図してみる。
> library(tidyverse)
> lcpm.cutoff <- log2(10/M + 2/L)
> df_lcpm<- lcpm %>% data.frame() %>% gather()
> ggplot(df_lcpm, aes(x=value, color = key))+
+ geom_density()+
+ geom_vline( aes(xintercept = lcpm.cutoff), linetype = "dashed")+
+ xlab("Log-cpm") +
+ ggtitle("Raw Data") +
+ scale_x_continuous(limits = c(-10, 15))+
+ theme_classic()
・かなり小さい値に偏っていることが分かる。それも、すべてのサンプルで発現量が低いものが目立つ。発現量の小さすぎるデータは解析の邪魔になるので、除去する。edgeRパッケージのfilterByExpr関数を使う。
> keep.exprs <- filterByExpr(x, group = group)
> keep.exprs %>% head()
・keep.exprsの中身には、どの遺伝子を解析のために残すのか、という情報が含まれている。データセットxをkeep.exprsでフィルタリングする。ちなみに、keep.lib.sizes = FALSEはフィルタリングしたあとにライブラリサイズを再計算する引数らしい。
> x %>% dim()#フィルタリング前の遺伝子数を見てみる。
[1] 27179 9
> x <- x[keep.exprs,, keep.lib.sizes=FALSE]
> x %>% dim()#フィルタリング後の遺伝子数を見てみる。
[1] 16624 9
・もう一度分布をみてみる（グラフを描くコードは割愛）。明らかに分布が変わった。
lcpm <- cpm(x, log = TRUE)#log-cpmを計算しなおす。
・次に、遺伝子発現の分布をnormalizeする。その前に、箱ひげ図でnormalize前の分布をみておく。
> ggplot(df_lcpm, aes(y=value, fill=key))+
+ geom_boxplot()+
+ ggtitle("Un-normalized Data") +
+ guides(x = "none")+#x軸は邪魔なので消す。
+ theme_classic()
・何もしなくてもとても揃っている。次に、normalizeする。TMM normalizeに関しては、以下を参照。
> x <- calcNormFactors(x, method = "TMM")#TMMでnormalizeする。
> x$samples$norm.factors#補正するための係数を確認。
[1] 0.8943956 1.0250186 1.0459005 1.0458455 1.0162707 0.9217132 0.9961959 1.0861026 0.9839203
・次に、normalized版の箱ひげ図も描く。グラフを描くコードは割愛。気持ち揃った…？
> lcpm <- cpm(x, log=TRUE)
つづく。
