ターミナルを開いて、local/work/以下にRディレクトリをさくせいして、Rにて作業する。ちなみにMacsまでのデータはlocal/workにある。
cd local/work mkdir R cd R Rとすると以下のようになる。
R version 3.1.0 (2014-04-10) -- "Spring Dance" Copyright (C) 2014 The R Foundation for Statistical Computing Platform: x86_64-unknown-linux-gnu (64-bit) R は、自由なソフトウェアであり、「完全に無保証」です。 一定の条件に従えば、自由にこれを再配布することができます。 配布条件の詳細に関しては、'license()' あるいは 'licence()' と入力してください。 R は多くの貢献者による共同プロジェクトです。 詳しくは 'contributors()' と入力してください。 また、R や R のパッケージを出版物で引用する際の形式については 'citation()' と入力してください。 'demo()' と入力すればデモをみることができます。 'help()' とすればオンラインヘルプが出ます。 'help.start()' で HTML ブラウザによるヘルプがみられます。 'q()' と入力すれば R を終了します。 >となってRが使えるようになる。Rが起動している時は>が左についている状態にある。簡単なコマンドについてはググれば色々出てくる。計算や変数への代入などにある程度慣れてから、bioconductorをインストールしてみた。とりあえず書いてある通りやってみる。
>source("http://bioconductor.org/biocLite.R")
URL 'http://www.bioconductor.org/packages/2.14/bioc/src/contrib/BiocInstaller_1.13.3.tar.gz' を試しています
Content type 'application/x-gzip' length 14185 bytes (13 Kb)
開かれた URL
==================================================
downloaded 13 Kb
* installing *source* package ‘BiocInstaller’ ...
** R
** inst
** preparing package for lazy loading
** help
*** installing help indices
** building package indices
** testing if installed package can be loaded
Bioconductor version 2.14 (BiocInstaller 1.13.3), ?biocLite for help
* DONE (BiocInstaller)
The downloaded source packages are in
‘/tmp/RtmpIpHodt/downloaded_packages’
Updating HTML index of packages in '.Library'
Making 'packages.html' ... done
Bioconductor version 2.14 (BiocInstaller 1.13.3), ?biocLite for help
>
となる。数分かかる。ChIPpeakAnnoというパッケージをダウンロードする。
>biocLite("ChIPpeakAnno")
BioC_mirror: http://bioconductor.org
Using Bioconductor version 2.14 (BiocInstaller 1.13.3), R version 3.1.0.
Installing package(s) 'ChIPpeakAnno’
以下色々続く
結構かかる５分くらい。つぎにChIPpeakAnnoを使うにはRを起動するたびに以下のコマンドを叩く必要があるらしい
>library("ChIPpeakAnno")
要求されたパッケージ biomaRt をロード中です
要求されたパッケージ IRanges をロード中です
要求されたパッケージ Biostrings をロード中です
要求されたパッケージ XVector をロード中です
10秒くらいでおわる。次に解析手順。
MACSから出てきたピークのbedファイルをread.tableで読み込む。
それをRのパッケージで使いやすい形式にBED2RangedDataで変換する。
>SoxDay0.df <- read.table("../SoxDay0_peaks.bed", header= FALSE)
> head(SoxDay0.df)
V1 V2 V3 V4 V5
1 chr1 3062864 3063136 MACS_peak_1 81.45
2 chr1 3482902 3483152 MACS_peak_2 65.14
3 chr1 4150841 4151273 MACS_peak_3 222.85
4 chr1 4660129 4660269 MACS_peak_4 53.63
5 chr1 4802472 4802926 MACS_peak_5 193.69
6 chr1 4855364 4855775 MACS_peak_6 270.83
> SoxDay0.gr <- BED2RangedData(SoxDay0.df, header= FALSE)
どんな中身かを少しだけ見るにはhead()とする。
>head(SoxDay0.gr)
RangedData with 6 rows and 2 value columns across 23 spaces
space ranges | strand score
<factor> <IRanges> | <character> <numeric>
MACS_peak_1 1 [3062864, 3063136] | + 81.45
MACS_peak_2 1 [3482902, 3483152] | + 65.14
MACS_peak_3 1 [4150841, 4151273] | + 222.85
MACS_peak_4 1 [4660129, 4660269] | + 53.63
MACS_peak_5 1 [4802472, 4802926] | + 193.69
MACS_peak_6 1 [4855364, 4855775] | + 270.83
このフォーマットにしておくと今後色んな解析を行うのに便利らしい。ここまでのデータをセーブしておく。(save())>save(list=ls(), file = "./bed2rangedData.rdat”)次にTSSのデータがパッケージ内に入っているのでdata()としてそれを読み込む。annotatePeakInBatchを使ってそのTSSデータと先ほど作成したRangedDataを結びつけてEnsEMBL Gene ID、遺伝子までの距離、上流・下流の情報をアノテーションする。
>data(TSS.mouse.GRCm38)
>SoxDay0.anno <- annotatePeakInBatch(RangedData(SoxDay0.gr), AnnotationData = TSS.mouse.GRCm38,output="both")
>head(SoxDay0.anno)
RangedData with 6 rows and 9 value columns across 23 spaces
space ranges | peak
<factor> <IRanges> | <character>
MACS_peak_1 ENSMUSG00000090025 1 [ 3062864, 3063136] | MACS_peak_1
MACS_peak_10 ENSMUSG00000033740 1 [ 6458167, 6458317] | MACS_peak_10
MACS_peak_100 ENSMUSG00000089964 1 [74824905, 74825354] | MACS_peak_100
MACS_peak_101 ENSMUSG00000082012 1 [76036406, 76036593] | MACS_peak_101
MACS_peak_106 ENSMUSG00000094946 1 [83844615, 83844958] | MACS_peak_106
MACS_peak_107 ENSMUSG00000062590 1 [86187886, 86188233] | MACS_peak_107
strand feature start_position
<character> <character> <numeric>
MACS_peak_1 ENSMUSG00000090025 + ENSMUSG00000090025 3054233
MACS_peak_10 ENSMUSG00000033740 + ENSMUSG00000033740 6487231
MACS_peak_100 ENSMUSG00000089964 + ENSMUSG00000089964 74850895
MACS_peak_101 ENSMUSG00000082012 + ENSMUSG00000082012 75636390
MACS_peak_106 ENSMUSG00000094946 + ENSMUSG00000094946 83795912
MACS_peak_107 ENSMUSG00000062590 + ENSMUSG00000062590 86154780
end_position insideFeature distancetoFeature
<numeric> <character> <numeric>
MACS_peak_1 ENSMUSG00000090025 3054733 downstream 8631
MACS_peak_10 ENSMUSG00000033740 6860940 upstream -29064
MACS_peak_100 ENSMUSG00000089964 74879953 upstream -25990
MACS_peak_101 ENSMUSG00000082012 75637179 downstream 400016
MACS_peak_106 ENSMUSG00000094946 83796032 downstream 48703
MACS_peak_107 ENSMUSG00000062590 86278284 inside 33106
shortestDistance fromOverlappingOrNearest
<numeric> <character>
MACS_peak_1 ENSMUSG00000090025 8131 NearestStart
MACS_peak_10 ENSMUSG00000033740 28914 NearestStart
MACS_peak_100 ENSMUSG00000089964 25541 NearestStart
MACS_peak_101 ENSMUSG00000082012 399227 NearestStart
MACS_peak_106 ENSMUSG00000094946 48583 NearestStart
MACS_peak_107 ENSMUSG00000062590 33106 NearestStart
上流や下流の情報でまとめたテーブルを作ってみる。
>SoxDay0.anno.table <- table(as.data.frame(SoxDay0.anno)$insideFeature)
>head(SoxDay0.anno.table)
downstream includeFeature inside overlapEnd overlapStart
1197 1 1228 14 14
upstream
1209
これを画像として出力させるには以下のようにする。
>png("SoxDay0.anno.barplot.png")
>barplot(SoxDay0.anno.table)
>dev.off()
ちょっと待つnull device
1
とでる。local/work/R内にSoxDay0.anno.barplot.pngが出来てる
これをscpでPCに転送
scp kosugi@サーバー名:/home/kosugi/local/work/R/SoxDay0.anno.barplot.png /Users/TK/Analysis/ChIPpeakAnno/などとすれば見れる。こんな感じ。
