・お題：病変部と非病変部で何か違うのか調べたいときに、遺伝子発現の違いを調べることがあるらしい。Differential expression analysis（DEG）という解析らしいのだけれど、せっかく生物学を学んでいるのだから、その辺も少し齧ってみたい。
・私はWet専門で学んできた人間だから、この辺の解析方法にはとんと疎い。データベースにアクセスして自分で試していければ良いのだけれど、全く自信がないので、親切なHPをなぞる形で学んでいきたいと思う。もちろん1回ではとても終わりそうにない。
・ということで、以下のHPを順になぞって学んでいきたいと思う。
・この例では、RNA-seqで発現量の解析をやっているらしい。RNA-seqに関しては、以下を参考に学んだ。次世代シーケンサー大手（たぶん）のIlluminaさんの資料。要するに、RNAを逆転写して作成したcDNAの配列を読んで読んで読みまくり、リファレンスゲノムのどこがどれだけ出てきたのかを調べる方法らしい。
・今回解析の例にしたのは、以下の文献。3種類の状態のマウス（virgin, pregnant and lactating mice）に関して、それぞれ2種類の細胞（basal stem-cell enriched cells (B) and committed luminal cells (L)）をとってきて、遺伝子発現量を解析したっぽい。正しいことは元文献を読んでほしい。
・次に、データをダウンロードする。今回は、"RNA-seq analysis in R"のHPにUpしていただいているデータセットをそのまま使わせていただく。ワーキングディレクトリに"data"というフォルダを作り、以下の3つのファイルをダウンロード＆保存する。
GSE60450_Lactation-GenewiseCounts.txt
SampleInfo.txt
SampleInfo_Corrected.txt
・次に、解析に必要なRのパッケージをダウンロードする。Rでこの手の解析をする際の王道はBioconductorらしい（以下）。BiocManagerというパッケージをインストールし、個々の様々な解析パッケージを管理する感じみたい。
・以下のコードで、今回使用するパッケージをインストールした。
if (!requireNamespace("BiocManager"))
install.packages("BiocManager")
BiocManager::install(c("limma", "edgeR", "Glimma", "org.Mm.eg.db", "gplots", "RColorBrewer", "NMF", "BiasedUrn"))
・BiocManager::installは恐らくBiocManagerのinstall関数を使ってパッケージをインストールしてくれ的な意味だと思うのだけれど、前半が分からない。。調べてみると、以下のQAサイトを見つけた。気になる方は確認して欲しい。
・次に、ライブラリをロードする。一度library()を実行すると、表示がワーッとコンソールに出てくる。もう一度全く同じコマンドでライブラリをロードしてエラーが出なければ、ライブラリのロードに問題はたぶんないことを確認できる。
library(edgeR)
library(limma)
library(Glimma)
library(org.Mm.eg.db)
library(gplots)
library(RColorBrewer)
library(NMF)
・次に、データをロードする。tab区切りのデータなので、delimで読み込む。カウントデータはstringsAsFactors をFALSE、サンプル情報はTRUEにしているっぽい（factorにしてよいことが何か理解が追い付いていないので、おいおい勉強していきたい）。
#カウントデータの本体を読み込む。HPのコードのファイル名が間違っている（不要なハイフンが入っている）ので、そこは直しておく。
seqdata <- read.delim("data/GSE60450_LactationGenewiseCounts.txt", stringsAsFactors = FALSE)
#サンプル情報を読み込む。
sampleinfo <- read.delim("data/SampleInfo.txt", stringsAsFactors = TRUE)
・それぞれ中身を見てみる。まずはseqdata 。
> str(seqdata)
> View(seqdata )
・遺伝子のID（EntrezGeneID）ごとに行があり、それぞれに遺伝子の長さ（Length）が含まれる。残りの列には各実験サンプル（マウスの状態3種類×細胞2種類×2 replicates）の読み取り数が含まれているらしい。
・各実験サンプルの情報は、sampleinfoに格納されている。
> View(sampleinfo)
・この後、読み取り数を加工し、解析していくことになる。seqdataから読み取り数だけを取り込んだデータフレームcountdataを作成する。
> countdata <- seqdata[,-(1:2)]#1列目と2列目を削除したデータを作成。
> rownames(countdata) <- seqdata[,1]#インデックスにEntrezGeneIDを設定。
> View(countdata)
・列名が長すぎるので、頭の7文字だけとってきて列名を付けなおす。これでsampleinfoにあったSampleNameと列名が同じになった。
> colnames(countdata) <- substr(colnames(countdata),start=1,stop=7)
・次に、countdataをDGEListオブジェクトに変換し、yと名付ける。DGEListはedgeRという解析用パッケージでカウントデータをハンドリングするためのオブジェクト形式らしい。
> y <- DGEList(countdata)
> y
・この後、DGEListオブジェクトを弄って群間比較したりすることになる。群間比較するからには、群を分けるための情報を入れる必要がある。yでいうと$sampleのgroupという項目がそれにあたるのだけれど、すべて1になっており、群の区別がついていない。
・先に見たように、今回の実験は、マウスの状態3種類×細胞2種類の6群構成（2 replicates）になっている。したがって、マウスの状態と細胞の種類の情報を使って6つの群ラベルを作成し、groupという項目に放り込む。
・とりあえず、マウスの状態と細胞の種類の情報を使って6つの群ラベルを作成する。参考サイトに従い、「細胞の種類.マウスの状態」のgroupNameを作成する。
> groupName <- paste(sampleinfo$CellType,sampleinfo$Status,sep=".")
> groupName
[1] "luminal.virgin" "basal.virgin" "basal.pregnant"
[4] "basal.pregnant" "basal.lactate" "basal.lactate"
[7] "basal.virgin" "luminal.virgin" "luminal.pregnant"
[10] "luminal.pregnant" "luminal.lactate" "luminal.lactate"
・次に、yのsampleのgroupに、groupNameをfactorとして放り込む。
> y$samples$group <-factor(groupName)
> y$samples
・この後データをいろいろ弄っていくわけなのだけれど、どの行がなんの遺伝子なのか、全く分からない。ENTREZIDが遺伝子が何かの情報なのだけれど、数字の羅列なのでさっぱり分からない。そこで、どのENTREZIDが何の遺伝子なのか、調べておく。これにはorg.Mm.eg.dbというパッケージを使う。
・org.Mm.eg.dbには、遺伝子に関するいろいろな注釈情報がくっついている。
> columns(org.Mm.eg.db)
"ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS" "ENTREZID" "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME" "GENETYPE" "GO" "GOALL" "IPI" "MGI" "ONTOLOGY" "ONTOLOGYALL" "PATH" "PFAM" "PMID" "PROSITE" "REFSEQ"
"SYMBOL" "UNIPROT"
・この中から、ENTREZIDをキーにして、この中から"ENTREZID"、 "SYMBOL" と "GENENAME"を釣ってきたい。
> ann <- select(org.Mm.eg.db,keys=rownames(y$counts),columns=c("ENTREZID","SYMBOL","GENENAME"))
> View(ann)
・ENTREZIDをキーにして、ENTREZIDを釣ってきている理由は、本当に求めた情報が取れているかチェックするため。理屈上では、rownames(y$counts)とann$ENTREZIDが一致するはず。チェックしてみる。
> ann$ENTREZID==rownames(y$counts)
[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE ....
・TRUEがずらっと続く。あまりに分かりにくいので、TRUEが何個あるのか集計する。table関数を使う。
> table(ann$ENTREZID==rownames(y$counts))
TRUE
27179
・TRUEばかり27179個あった。問題なさそう。データセット作成の最後に、yに先ほど作った遺伝子情報をgenesとして放り込む。
> y$genes <- ann
> y
・さっきはなかったgenesという項目がyの中にできている。
つづく。
