アメリエフのブログ

バイオインフォマティクスの紹介と社員の日々
重複値の削除や抽出
Rで、ベクトルなどに含まれる重複した値をユニークにしたいとき、unique()を使用します。
> a
[1] "AAA" "BBB" "CCC" "AAA"
> unique(a)
[1] "AAA" "BBB" "CCC"

linuxコマンドの、隣り合っている重複行をユニークにするuniqコマンドとよく混同して、uniqと書いてしまうのは私だけでしょうか。
$ cat file.txt
cat
cat
dog
mouse
$ uniq file.txt
cat
dog
mouse

linuxのuniqコマンドでは「-d」オプションで重複している値を取得することができます。
$ uniq -d file.txt
cat

これが頭にあった私は、Rのunique()でも同じようにオプションで指定するのだろうと思い込み、unique()のヘルプを舐めるように読みましたが、そんなオプションはありませんでした。
しかしRのヘルプは親切なので、一番最後の方で、「重複した値の取得にはduplicated()を使えるよ!」と教えてくれていました。助かります。
duplicated()は、真偽値を返します。要素が1つしかないときはFALSE、同じ要素が2回以上出現するとき、2回目以降の要素についてTRUEを返します。

> a
[1] "AAA" "BBB" "CCC" "AAA"
> duplicated(a)
[1] FALSE FALSE FALSE TRUE
> a[duplicated(a)]
[1] "AAA"
| kubo | 統計解析ソフトR | 15:29 | comments(0) | - |
seq2pathwayでアノテーション(後)
Bioconductorのseq2pathwayパッケージを使ってパスウェイのアノテーションを行います。
※seq2pathwayのインストールはこちらの記事をご覧ください。

seq2pathwayパッケージに含まれているテストデータを使ってテストしてみましょう。
library(seq2pathway)
data(dat_chip)


こんなデータです。遺伝子名とスコアのセットです。
head(dat_chip)



これらの遺伝子に対して、パスウェイをアノテートします。
遺伝子情報は、seq2pathwayに入っているMsigDB_C5データを使います。
data(MsigDB_C5,package="seq2pathway.data")
result_FAIME <- gene2pathway_test(dat=dat_chip, DataBase=MsigDB_C5, FisherTest=T, EmpiricalTest=T, method="FAIME", alpha=5, logCheck=F, na.rm=F)


結果は、2つのデータフレームから成るリストです。
names(result_FAIME)



データフレームgene2pathway_result.FETの1列目、3列目、9列目にパスウェイ名、Fisher検定のp値、遺伝子のリストが入っています。
head(result_FAIME$gene2pathway_result.FET[,c(1,3,9)])


※クリックで拡大します

gene2pathwayではGeneOntologyのアノテーションも行えます。

遺伝子・パスウェイ・GeneOntologyのアノテーションを行うBioconductorのパッケージには、以前紹介したBiomaRtReactomePAなどさまざまなパッケージがありますが、gene2pathwayも使いやすいと思いました。
| hat | 統計解析ソフトR | 15:14 | comments(0) | - |
seq2pathwayでアノテーション(前)
ChIP-seqやRIP-seqでピークが得られた後は、遺伝子をアノテーションしたり、それらの遺伝子がどんなパスウェイに関連するか調べたりしたいですね。

Bioconductorのパッケージseq2pathwayを使って、ゲノム領域に遺伝子をアノテーションしたり、遺伝子にパスウェイをアノテートしたりすることができます。

seq2pathwayのインストールはおなじみbiocLiteで行います。
source("http://bioconductor.org/biocLite.R")
biocLite("seq2pathway")


seq2pathwayパッケージに含まれているテストデータを使ってテストしてみましょう。
library(seq2pathway)
data(Chipseq_Peak_demo)


こんなデータです。
head(Chipseq_Peak_demo)



これらの領域に対して、遺伝子をアノテートします。
Chipseq_anno <- runseq2gene(inputfile=Chipseq_Peak_demo, genome="hg19", adjacent=F, SNP=F, search_radius=150000, PromoterStop=F, NearestTwoDirection=F)


※実行時にPython2.7が必要というメッセージが出た場合は、https://www.python.org/getit/からお使いの環境に合わせたPython2.7をインストールしてください。

このような結果が得られます。領域に対し、遺伝子のアノテーションがつきました。
head(Chipseq_anno[[1]])


※クリックで拡大します

次回はこのseq2pathwayを使ってパスウェイのアノテーションを行ってみたいと思います。
| hat | 統計解析ソフトR | 13:30 | comments(0) | - |
ヘッダに注意
バイオインフォには関係ないのですが、少しつまづいたことです。

ReactomePAを使っています。
hatさんがかかれた紹介記事では、enrichPathway()でエンリッチされているパスウェイを取得した結果を棒グラフに描画していますが、summary()することでデータフレームにも変換できます。
> library(ReactomePA)
> library(DOSE)
> data(geneList)
> de <- names(geneList)[abs(geneList)>1]
> require(ReactomePA)
> x <- enrichPathway(gene=de, pvalueCutoff=0.05, readable=T)
> pathway.df <- summary(x)


データフレームはパスウェイ上にエンリッチされている遺伝子やp値などがまとめられている一覧表になっています。
> names(pathway.df)
[1] "ID" "Description" "GeneRatio" "BgRatio" "pvalue"
[6] "p.adjust" "qvalue" "geneID" "Count"


エクセルで開いてみたらわかりやすそうですので、データフレームをタブ区切りのテキストファイルに書き出しました。
> write.table(pathway.df, file = "reactomepa.txt", quote = FALSE, sep="¥t", col.names=TRUE)

このファイルをエクセルで開こうとするとエラーが起きてしまいます。
どうやら、ヘッダーが「ID」で始まるファイルを、エクセルは別のフォーマットのファイルだと判断することが理由のようです。

エクセルで開く前にヘッダを書き換えるか、最初からxlsxパッケージなどを使ってエクセル形式で出力したほうが良いようです。
| kubo | 統計解析ソフトR | 15:13 | comments(0) | - |
paste関数を使う
短い間ですが、弊社のトレーニングのTAを担当させていただいていた時期があります。
RのTAをやるなかで、とても便利なのに、初めてRに触れる方が戸惑われることが多いと感じたポイントのひとつに paste 関数があります。
今回は paste 関数をあれこれいじくりまわしてみます。

1.paste関数とは

複数の文字列同士を結合する関数です。

> paste("hello", "world")
[1] "hello world"

なぜpaste関数が必要かというと、Rではこのように文字列を単にprint関数に渡すことができないためです。

> print("hello","world")
以下にエラー print.default("hello", "world") : 'digits' 引数が不正です
追加情報: 警告メッセージ:
In print.default("hello", "world") : 強制変換により NA が生成されました


対話モードだとprint関数を省略することも多いので、こんな書き方をされるかもしれません。いずれにせよ、複数の文字列を並べただけでは結合できません。
> "hello","world"
エラー: 予想外の ',' です in ""hello","

paste関数は大事です。

2.paste関数の基本

使い方の簡単な例をお見せします。
数字も、文字列に変換されて結合されます。

> paste(108, "yen")
[1] "108 yen"

結合した文字列の間に挟まる文字は、sep オプションで指定可能です。指定がない場合は半角スペースが入ります。

> paste(108, "yen", sep="-")
[1] "108-yen"

文字を挟みたくないときは、sep=""と指定します。

> paste(108, "yen", sep="")
[1] "108yen"


引数には、複数の要素があるベクトルを渡すこともできます。この場合、それぞれの要素について結合されます。

■ 1対複数
> pre
[1] "gene"
> str
[1] "a" "b" "c" "d" "e"
> paste(pre, str)
[1] "gene a" "gene b" "gene c" "gene d" "gene e"

■ 複数対複数(要素の数が同じ場合)
> num
[1] 1 2 3 4 5
> suf
[1] "st" "nd" "rd" "th" "th"
> paste(num, suf, sep="")
[1] "1st" "2nd" "3rd" "4th" "5th"

■ 複数対複数(要素の数が一致しない場合)…要素が足りないほうのベクトルが使いまわされます。
> paste(num, suf[1:4], sep="")
[1] "1st" "2nd" "3rd" "4th" "5st"

ベクトル同士ではなく、ひとつのベクトルの中の各要素を結合する場合は、collapseオプションを使用します。
collapseオプションに、結合する要素同士の間に挟まれる文字を指定します。

> day
[1] "2009" "July" "1"
> paste(day, collapse="-")
[1] "2009-July-1"

sepオプションとcollapseオプションを両方指定したらどうなるでしょう?
> paste(num, suf, sep="", collapse="/")
[1] "1st/2nd/3rd/4th/5th"

ベクトル同士を結合した後、各要素を結合していますね。

もう少し実践的な例として、for文とグラフでの使い方をみていきます。
続きを読む >>
| kubo | 統計解析ソフトR | 15:16 | comments(0) | - |
heatmapの横に色をつけるオプション
heatmap関数のRowSideColors/ColSideColorsオプションの紹介です。
gplotsパッケージのheatmap.2やheatmap3パッケージでも同様の機能があります。
ヒートマップのデンドログラムの横に、カラフルなサイドバーを付ける機能です。

テスト用のマトリクスを作ります。
(結果は実行するたびにランダムに変わります)
mat<-matrix(sample(100,70,replace=T),ncol=7)
rownames(mat)<-c("geneA","geneB","geneC","geneD","geneE","geneF","geneG","geneH","geneI","geneJ")
colnames(mat)<-c("sample1","sample2","sample3","sample4","sample5","sample6","sample7")

このマトリクスでヒートマップを作成するとこうなります。
heatmap(mat)


以下は左にあるデンドログラムの横に、サンプルごとに色を付ける場合の例です。

RowSideColorsオプションに、色名(16進数での指定も可)のベクトルを指定します。
この色名のベクトルは、マトリクスの行数と同じ長さである必要があります。
nrow(mat)
[1] 10

また、色はマトリクスの各行に対して指定していますので、ヒートマップではクラスタリングされた順に並びかえられて描画されます。
rownames(mat)
[1] "geneA" "geneB" "geneC" "geneD" "geneE" "geneF" "geneG" "geneH" "geneI" "geneJ"   ←この順に色が指定されます
rcol<-c("red","blue","green","goldenrod","violetred","salmon","orange","springgreen","pink","black")
heatmap(mat,RowSideColors=rcol)


クラスタリングされたグループごとに色をつけたい場合は、マトリクスをあらかじめクラスタリングしたうえでcutree関数を使います。
mat.d<-dist(mat)
mat.h<-hclust(mat.d)
mat.c<-cutree(mat.h,k=3)
cutree関数のkで指定した数値と同じ長さの色名のベクトルをRowSideColorsに指定します。
clust.col<-c("red","blue","orange")
heatmap(mat,RowSideColors=clust.col[mat.c])


ColSideColorsの使い方も同様です。
なお、heatmap.2のlmatオプションを使用した場合のColSideColors / RowSideColorsの作図順については、以前記事を書きましたのでご参照ください。
| kubo | 統計解析ソフトR | 15:39 | comments(0) | - |
一部のデータだけ抜き出す
Rで解析をしていると、データフレーム(やマトリクスやベクトル)から、条件を満たすデータだけを抜き出したいことがよくあります。
そんなときはsubset関数を使います。

例として、みんな大好きirisを使います。
> head(iris)
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 5.1 3.5 1.4 0.2 setosa
2 4.9 3.0 1.4 0.2 setosa
3 4.7 3.2 1.3 0.2 setosa
4 4.6 3.1 1.5 0.2 setosa
5 5.0 3.6 1.4 0.2 setosa
6 5.4 3.9 1.7 0.4 setosa

続きを読む >>
| hat | 統計解析ソフトR | 15:12 | comments(0) | - |
heatmap.2のlmatオプション
Rにはデフォルトでもヒートマップを作成するheatmap関数がありますが、私はheatmap関数より多機能なgplotsパッケージのheatmap.2関数をよく使います。

基本的な使い方はheatmap関数と同じです。
gplotsパッケージをインストールしていない場合はインストールから始めます。
install.packages("gplots")
library(gplots)
引数はnumericのmatrixです。
ここではRに組み込まれているデータセットmtcarsをmatrixに変換して使います。
input<-as.matrix(mtcars)
heatmap.2(input)
heatmap.2で作図すると、デフォルトでは左上にカラーキー、右上にinputのcolをクラスタリングしたデンドログラム、左(左下)にinputのrowのデンドログラム、右(右下)にヒートマップが描画されます。この配置はheatmap.2のlmatオプションで変えることができます。
lmatオプションの使い方はlayout関数と同じで、作図の順番を指定するためのmatrixを指定して使います(heatmap.2内部でlayout関数を使っています)。作図される順は下図の通りです。
heatmap2(default)lmatを変更すると、layout関数のheights、widthsオプションに相当するlheiオプションとlwidオプションが要求されるので、適当な数字を指定する必要があります。

カラーキーを右上に移動してみます。
lm<-matrix(c(0,2,3,1,4,0),ncol=3)
heatmap.2(input,lmat=lm,lwid=c(1,5,2),lhei=c(1,4))

ただし、 ColSideColors/RowSideColorsオプションを使うと、これらの作図順が“屬砲覆蝓以降のの順番がずれます。
heatmap.2(input,ColSideColors=col)
# colで任意の色を指定します
heatmap2(sidecol)
両方使用した場合、RowSideColors ColSideColorsの順になります。
私はここにひっかかって苦労したので、紹介してみました。
| kubo | 統計解析ソフトR | 14:25 | comments(0) | - |
RSQLite依存パッケージインストールのエラーについて
ある環境で使用していたとあるBioconductorパッケージを、別の環境でも使うためインストールしようとしたらエラーが出ました。
具体的にはggbioとbiovizBase、他にもhatさんご紹介のReactomePAもダメでした。

エラー文を見てみると、「RSQLite」パッケージの名前が出てきます。
調べたところ、どうやらこのRSQLiteパッケージは昨年10月にアップデートされています。その結果、RSQLiteに依存しているパッケージに影響があったようです。

最新のR-3.1以降のバージョンで使えるbiocLiteでは、既にこの問題に対応済みのパッケージをインストールできるようです。
これ以前のバージョンのRでは、既にインストール済みのパッケージは(何かのはずみでRSQLiteをバージョンアップしない限り)動作に問題がありませんが、新しくRSQLiteに依存しているパッケージをインストールすることができないようです。

今まで問題なかったものが、突然うまくいかなくなると焦りますね。
私の場合は、最新のバージョンのRを使うことにしました。
| kubo | 統計解析ソフトR | 15:18 | comments(0) | - |
続・Rでgrep
以前、Rのgrep()についての記事で、grep()は文字列のベクトルにのみ使えるとご紹介しましたが、今回はマトリクスやデータフレームではどうしたらよいか?という記事です。

例として以下のような文字列のマトリクスを用意しました。
1列目が季節、2列目が英語の月名の12行の文字列のマトリクスです。
> monthly
  season month
[1,] "winter" "January"
[2,] "winter" "February"
[3,] "spring" "March"
[4,] "spring" "April"
[5,] "spring" "May"
[6,] "summer" "June"
[7,] "summer" "July"
[8,] "summer" "August"
[9,] "autumn" "September"
[10,] "autumn" "October"
[11,] "autumn" "November"
[12,] "winter" "December"


さて突然ですが、牡蠣はお好きでしょうか?
嘘か本当か、牡蠣はrのつかない月に食べてはいけないと言われています。
そこで、用意したマトリクスから、牡蠣を食べてはいけない月を抜き出してみようと思います。
まず、月名のベクトル、つまりマトリクスの2列目のみを取得します。
> monthly[,2]
[1] "January" "February" "March" "April"
[5] "May" "June" "July" "August"
[9] "September" "October" "November" "December"


さっそく、このベクトルからrのつく月をgrep()で抜き出しましょう。
続きを読む >>
| kubo | 統計解析ソフトR | 16:00 | comments(2) | - |
Rでgrep
Rのオブジェクトから、特定のキーワードを含む要素だけ抜き出したいときはどうしたらいいでしょうか。

ベクトルの場合はgrep()関数が使用できます(grepの仲間にgrepl、regexpr、gregexpr、regexec、検索後に置換を行うsub、gsubがあります)。
このgrepは文字列からなるベクトルのみに使用できます。

たとえば、Rで図を作成するときに色名指定に使用できる色名の一覧を返してくれるcolors()関数というのがあります。
> iro <- colors()
> head(iro)

[1] "white" "alicered" "antiquewhite" "antiquewhite1" "antiquewhite2" "antiquewhite3"
> is.vector(iro)
[1] TRUE
> class(iro)
[1] "character"
データ構造は文字列から成るベクトルなので、これでgrepを試してみます。
続きを読む >>
| kubo | 統計解析ソフトR | 14:43 | comments(0) | - |
遺伝子にGeneOntologyのアノテーションをつける(2)
前回はbiomaRtパッケージを使って、「ADH1B」「ALDH2」遺伝子にGeneOntologyのアノテーションをつけました。

今回は、RamiGOパッケージを使って、それらのGeneOntologyIDがどのような関係にあるか、図示してみたいと思います

RamiGOパッケージをインストールします。
>biocLite('RamiGO')

前回取得したデータから、GeneOntologyIDだけを抜き出します。
>go <- data$go_id

これらのGeneOntologyIDについて、GeneOntology階層中の位置を描画します。
>library(RamiGO)
>getAmigoTree(goIDs=go, color='yellow', filename='gene2go')


gene2go.png
 ←クリックで拡大します

局在(Cellular Component)については細胞質基質やミトコンドリア マトリックス、分子機能(Molecular Function)については酸化還元、プロセス(Biological Process)についてはアルコール代謝に関係がありそうです。

このようにGeneOntologyでアノテーションをつけることによって、(例えばある病気で特異的に発現する)いくつかの遺伝子について、どのような傾向があるかを調べることができます。

生命科学の研究に携わっているR初心者の方にRの基本的な操作方法を学んでいただく、アメリエフ・バイオインフォマティクス・スクール「R基礎」の開催日程はこちらです。
| hat | 統計解析ソフトR | 14:26 | comments(0) | - |
遺伝子にGeneOntologyのアノテーションをつける(1)
_
_
▲▲▲ R Advent Calendar 2013に参加させていただいております。
_

Bioconductorは、バイオ解析関連のRパッケージを集めたサイトです。
750個以上(2013年12月時点)の有用なパッケージが登録されています。

複数の遺伝子の特徴の傾向を見たい時に、GeneOntologyで機能や局在を調べることがよくあります。
BioconductorのbiomaRtパッケージを使って、遺伝子にGeneOntologyのアノテーションをつけることができます。

今回は、お酒の強さに関係があるといわれる「ADH1B」「ALDH2」遺伝子に対応するGeneOntologyIDを調べてみましょう。

biomaRtパッケージをインストールします。
>source('http://bioconductor.org/biocLite.R')
>biocLite('biomaRt')


「ADH1B」「ALDH2」に対応するGeneOntologyIDを取得します。
>library(biomaRt)
>ens <- useMart('ensembl', dataset='hsapiens_gene_ensembl')
>data <- getBM(attributes=c('hgnc_symbol','go_id','name_1006'),
+filter='hgnc_symbol', values=c('ADH1B', 'ALDH2'), mart=ens)


取得したデータを見てみましょう。
>edit(data)
 ←クリックで拡大します

予想通り、アルコール代謝に関連ありそうなGeneOntologyIDがヒットしています。

次回、これらのGeneOntologyIDがどのような関係にあるか、BioconductorのRamiGOパッケージを用いて描画してみたいと思います。

生命科学の研究に携わっているR初心者の方にRの基本的な操作方法を学んでいただく、アメリエフ・バイオインフォマティクス・スクール「R基礎」の開催日程はこちらです。
| hat | 統計解析ソフトR | 14:00 | comments(0) | - |
特定の列をキーに行列を並び替える
ベクトルやリストを並び替えるときはsort()関数を使いますが
行列(マトリクス)やデータフレームを、特定の列をキーにして
並び替える場合はorder()という関数を使います。

同じジュースが、250ml・100円、350ml・120円、180ml・80円で
売られています。
1円あたりの容量が多いのはどれか、調べてみたいと思います。

まず、データを行列に入れます。

> juice <- matrix(c(250, 350, 180, 100, 120, 80), nrow=3)
> colnames(juice) <- c("volume", "price")
> juice

_____volume_price
[1,]____250___100
[2,]____350___120
[3,]____180____80

容量を価格で割った値を3列目に追加し、その列でソートします。

> juice <- cbind(juice, rate=juice[,1]/juice[,2])
> juice[order(juice[,3]),]

_____volume_price_____rate
[1,]____180____80 2.250000
[2,]____250___100 2.500000
[3,]____350___120 2.916667

僅差ですが、350mlが一番お得なようです。

なお、降順にソートしたい場合はorderのオプションに
「decreasing=TRUE」をつければOKです。

| hat | 統計解析ソフトR | 14:54 | comments(0) | - |
Rで家系図描画
家系解析を行っている方は、家系図を描く機会が多いのでは
ないでしょうか。
kinship2というパッケージを使うと、Rで家系図を描くことができます。

(1)家系図情報ファイルの作成
以下の内容のテキストファイルを作成します。
id father mother sex name
1 0 0 1 Namihei
2 0 0 2 Fune
3 1 2 2 Sazae
4 1 2 1 Katsuo
5 1 2 2 Wakame
6 0 0 1 Masuo
7 6 3 1 Tarao
項目区切りは空白でもタブでもOKです。
idが一意なid、fatherが父のid、motherが母のid、sexが性別
(男性は1、女性は2)、nameが名前です。

(2)「kinship2」パッケージのインストール
> install.packages('kinship2')

(3)ファイル読み込み、家系図描画
>library(kinship2)
> df <- read.table('c:/Users/hat/Desktop/ped.txt', header=T)
# (1)のファイルを指定します
> ped <- pedigree(id=df$id, dadid=df$father, momid=df$mother, sex=df$sex)
> plot(ped, id=df$name)




さらに、疾患の有無や生死の情報も付加することができます。
詳細はkinship2のvignettesをご覧ください。
| hat | 統計解析ソフトR | 15:09 | comments(0) | - |
あの花はどんな花
Rの非常に有名なテストデータの一つ「iris」は、
アヤメの花弁の長さ・幅、ガクの長さ・幅のデータです。
※アヤメのガクは、黄色い筋が入っていて外に垂れ下がっている部分です。

setosa(セトサ)種、versicolor(バージカラー)種、
virginica(バージニカ)種という3種について
各50サンプルの情報が入っているので、
種間で花の大きさを比較することができます。

Celebrating Wildflowersというサイトにこれら3種の
生息地や特徴や写真が載っていました。
セトサ種
アラスカなどに分布していて、葉は30センチくらい、
ガクのわりに花弁がとても小さいそうです。

バージカラー種
アメリカの北のほうに分布しているそうです。

バージニカ種
アメリカの南のほうに分布していて、葉は1m近くに
なるそうです。

北の花(ラベンダーとか)は小さく、南の花(ラフレシアとか)は
大きいイメージがありますが、やはり南にいくほどアヤメも
大きくなるのでしょうか。

irisの花弁(Petal)の長さと、ガク(Sepal)の長さの箱ひげ図を
Rで描いてみると、確かにセトサ種は「ガクのわりに花弁がとても
小さい」ようです。

> par(mfrow=c(1,2))
> boxplot(iris$Petal.Length~iris$Species, main="Petal Length")
> boxplot(iris$Sepal.Length~iris$Species, main="Sepal Length")




ぜひ日本のアヤメとも比較してみたいと思うので、
どなたか、50本のアヤメの花束を私に贈ってください!
| hat | 統計解析ソフトR | 14:57 | comments(0) | - |
SRAをRで検索
毎日寒いですね。
家ではこたつから出られなくなって、こたつの中で生活している方もいらっしゃるのではないでしょうか。

今日は、RでNCBI SRAのデータを検索するRパッケージ「SRAdb」をご紹介します。

最近出たばかりのSRAdbの論文はこちらです。
Zhu Y, Stephens RM, Meltzer PS, Davis SR.
SRAdb: query and use public next-generation sequencing data from within R.
BMC Bioinformatics. 2013 Jan 17;14(1):19.


・インストール
source("http://bioconductor.org/biocLite.R")
biocLite("SRAdb")

・ライブラリ読み込み
library(SRAdb)

・検索
sqlfile = getSRAdbFile()  # 10分くらいかかります
sra_con = dbConnect(SQLite(), sqlfile)
rs <- getSRA(search_terms = "Human RNA-Seq", out_types = c("run", "study"), sra_con = sra_con)
# ↑"Human RNA-Seq"で検索する場合の例
head(rs[, c(1:5, 7)])

・結果
run_alias run run_date updated_date spots run_center
1 poly_1 DRR000008 2012-09-06 19153889
2 poly_3 DRR000010 2012-09-06 18079621
3 poly_2 DRR000009 2012-09-06 17668118
4 cyto_2 DRR000012 2012-09-06 15075547
5 cyto_1 DRR000011 2009-06-12 2012-09-06 15788579 UT-MGS
6 cyto_3 DRR000013 2009-06-15 2012-09-06 15490013 UT-MGS

結果が返ってきました!


面白いと思ったのは、SRAdbには、起動中のIGVの表示領域を切り替える機能もあることです。

※予めIGVを起動しておき、「View」→「Preferences」→「Advanced」の「Enable port」にチェックが入っていて、「60151」になっていることを確認してください。

・IGVの表示領域切り替え
myBam = file.path("path/to/bam")  
# bamと同じ場所にbam.baiも置いておきます
sock <- IGVsocket()
IGVgenome(sock, "hg19")
IGVload(sock, myBam)
IGVgoto(sock, "chr1:1-1000")

今はDBCLS SRAなど便利なサイトがたくさんあるので、SRAデータの検索をあえてRで行いたい場面はあまり思いつかないのですが、面白いパッケージだと思いました。

「こたつ」ならぬ「Rに入ったままなんでもやりたい」方にはよさそうです。
| hat | 統計解析ソフトR | 14:44 | comments(0) | trackbacks(0) |
Rユーザ会に参加してきました
hatです。

先週土曜日、2012年度統計数理研究所共同研究集会「データ解析環境Rの整備と利用」の「国内ユーザーによる報告」に参加してきました。

普段はバイオの話ばかりなので、疫学やTwitter解析や感染症シミュレーションなど、他分野のお話がたくさん伺えておもしろかったです。

どこでも大容量データの取り扱いがネックになっていて、読み込み方法の工夫や並列化などで対処されているようです。

RはコマンドラインかWindowsのRGuiで使うものだと思っていたのですが、RStudioのような優れた開発環境や、ef-prime社のR Analytic Flowのような使い勝手のいいフロー管理ソフトがあることを知ったのも大きな収穫の一つです。
今後はぜひこのようなソフトを活用して、Rを使いこなしていきたいと思います。
| hat | 統計解析ソフトR | 09:16 | comments(0) | trackbacks(0) |
Rでつなぐ次世代オミックス情報統合解析研究会
明日、理化学研究所で開催される
“Rでつなぐ次世代オミックス情報統合解析研究会”
に参加してきます!
(詳細はコチラ

NGSの解析に特化したRを学べる事は、大変嬉しいです。
パワーアップをめざします!
| | 統計解析ソフトR | 19:04 | comments(0) | trackbacks(0) |
ご指摘ありがとうございます!(箱ひげ図)
箱ひげ図の概念に対し、ご指摘いただきました。
誠に申し訳ありません。訂正し、御詫び申し上げます。
またご指摘いただきまして、誠にありがとうございます。
下記に訂正、追加説明を行います。

箱ひげ図の概念図(訂正版)をFigure1に示します。


Figure1 箱ひげ図の概念図

御指摘がありましたように、 “ひげ”の上下末端は、“最大値”、“最小値”ではなく、それぞれ“(箱の長さ)×range(default=1.5)+第3四分位数 ”および“(箱の長さ)×range(default=1.5)-第1四分位数 ”を指します。

(ひげの定義は、書籍、参考資料によって異なります(上記の定義の他: “(箱の長さ)×range(default=1.5)+第3四分位数”内の最大数など)。そのため、上記の定義が絶対的に正しいとは申しません)

さらにRで箱ひげ図を作成する際にも何も指定することなく下記のようにコマンドを書きますと、ひげは記載されません(Figure 2)。

コマンド1
age <- c(22, 23, 24, 25, 40)
boxplot(age)


Figure2 コマンド1で作成した箱ひげ図

しかしboxplot()コマンドにはrangeというオプションがあります。こちらのrangeは、“箱の長さから何倍の場所にひげの長さを設定するか?”という意味となります。

rangeはdefaultでは1.5となっています。
そのため上記のコマンドの場合、第3四分位数=25であることから、ひげ上部の長さは、

(25-23)×1.5+25=28

となります。ageの変数内に28がないため、上のコマンドでは上のひげが表示されませんでした。

このオプションをrange=0とすると、ひげの上下は、数値データ群の“最大値”“最小値”を示します(Figure3)。

コマンド2
age <- c(22, 23, 24, 25, 40)
boxplot(age, range=0)


Figure3 コマンド2で作成した箱ひげ図


今回の御指摘、誠にありがとうございました。とても勉強になりました。
また不足の致すところがあれば、ご指摘いただけますと幸いです。
| | 統計解析ソフトR | 18:52 | comments(2) | trackbacks(0) |
中央値と箱ひげ図
昨日は、中央値について記述しました。
本日は、Rを用いた中央値の算出方法と箱ひげ図の作成について記述します。
よろしくお願いいたします。

例としまして
21歳, 23歳, 23歳, 25歳, 22歳, 60歳の方々の中央値を算出します。
中央値を算出する際は、median()というコマンドを用います。
age <- c(21, 23, 23, 25, 22, 60)
median(age)

→23(中央値です)

とても簡単に算出できました。

昨日も記述致しましたように、中央値を用いる時というのは、数値データの分布がskewed分布となる時です。せっかくだったら数値データの分布、ばらつきを見たい。
そんな時は、“箱ひげ図(box plot)”を作成します。


Figure 1 箱ひげ図概念

“箱ひげ図”の概念をFigure 1に示します。箱の中にある太線が中央値、箱の上下がそれぞれ数値データ内の3/4, 1/4位の値、そして箱から伸びている点線(これがひげ)がそれぞれ最大値、最小値を表します。
さっそくRで作成してみます。箱ひげ図は、boxplot()コマンドで作成する事ができます。
age <- c(21, 23, 23, 25, 22, 60)
age2 <- c(21, 23, 23, 25, 22)
boxplot(age, age2, names=c("with 60", "without 60"), boxwex=0.8, range=0)

今回は60歳の有無で比較してみました。上記の結果をFigure 2に示します。


Figure 2 年齢の箱ひげ図

“with 60”の箱ひげ図は、最大値が他の数値郡にたいして大きく離れています。このように箱ひげ図を用いる事で、数値データのばらつき具合を簡単に見る事ができます。

Rって便利(結局はソコ)。
| | 統計解析ソフトR | 17:59 | comments(2) | trackbacks(0) |
夏といえば!
最近鼻血がよく出るタノです。
先日Rスクリプトを作成するためのツール“Tinn-R”を紹介しました(参照)。
本日は、“Tinn-R”とRを用いてちょっと簡単なグラフ作成をしてみたいと思います。
よろしくお願いします。

突然ですが、夏といえば・・・・ダイエット!!をする方が多いですよね。私もその一人ですが、ダイエットの際に体重を管理するグラフがあればとても便利だと思います。
エクセルでも可能ではありますが、本日はRで作成します。
まず、エクセルで下図の様なcsvファイル(体重表)を作成します。今回はカロリーも入れてみました。


そして、“Tinn-R”を起動し、下記のようにスクリプトを作成します。
rm(list=ls())
graphics.off()

# read csv file
X <- read.table("weight.csv", sep=",", header=TRUE)

par(mai=c(1,1,1,1))
plot(X$Day, X$Weight, xlim=c(1,30), ylim=c(40,70), type="l", col="red", xlab="Day/day", ylab="Weight/kg")
par(new=TRUE)

plot(X$Day, X$Calorie, xlim=c(1,30), ylim=c(800,2500), type="l", col="blue", xlab="Day/day", ylab="", axes=FALSE)

axis(4)
mtext("Calories/kcal", side=4, line=2)

↓詳細はこちらをクリック♪


このスクリプトファイルを”weight.r”として適当な場所に保存し、Rで読み込みます(source()コマンドです)。
そうすると・・・・


これで、摂取カロリーと体重増加の関係がわかります!

Rは、plot()などのようにコマンド一つで、いろいろなグラフや図形を描く事ができます。
ですので、データ同士の相違や関連性を見るのに適している言語だと思います。
今回は、2つのデータを比較しましたが、もっと多くの情報も一度に比較する事ができるんですよ!!
一度、試してみてはいかがでしょうか
| | 統計解析ソフトR | 19:46 | comments(0) | trackbacks(0) |
Rでスクリプトを書くときは!!
こんにちは。
本日より、R漬けの日々を過ごします、タノです。
学生時代にやっていたとはいえ、大量のデータに四苦八苦しております。というわけで、本日は、私が愛用しているRスクリプトを書くためのR code editor、“Tinn-R”を紹介します。
Tinn-Rは、R言語用のメモ帳の様なもので、“.r”というファイルを作成する事が可能です。
例えば、Tinn-R上に下記のスクリプトを書き、“example.r”というファイルを作成します。

# Clear workspace
rm(list=ls())
# A little game
cat("Let's play a little game..¥n¥n")

# Set a number (can use a random generator later..)
k = 54
LP = 1
while (LP == 1) {
cat("¥nPick a number between 1 and 100!¥n")
cat("Selection: ")
choice = scan("",integer(0),nlines=1) # Waits for user input
if (choice == k) { # START FIRST IF
cat("¥nCORRECT!! Well done!!¥n")
LP = 0
} else {
if (choice < 1 | choice > 100) { # START SECOND IF
cat("¥nI SAID BETWEEN 1 AND 100!!!! TRY AGAIN!!¥n")
} else {
cat("¥nWRONG!!! Try again!¥n")
}
}
}


次にR(Rのインストール方法はコチラ)を起動し、ファイル→ディレクトリの変更から“example.r”の保存先にディレクトリーを設定します。
そしてsource("example.r")とコマンドを打ちます。
すると・・・・・楽しいゲームの始まりです!!!

Tinn-Rを使用することで、コマンドの確認が可能となり、スクリプトも見易くなると思います。
お試しください!

Tinn-Rは、コチラのサイトの“Tinn-R, edit code and run it in R”からダウンロードしてください。
| | 統計解析ソフトR | 20:08 | comments(1) | trackbacks(0) |
R 「関数textとlegend 5」
R 「関数textとlegend」では、文字列と凡例を書き込みました。
書きこんだ文字列が見にくいようなので、位置を調節します。
太字の部分を書き足しました。

# デバイス領域
png("110131_iris.png", height=700, width=500)
# 作図領域
par(mar=c(4,4,2,6),xpd=TRUE)
# プロット領域(散布図)
plot(iris[,1], iris[,3], xlab="iris[,1]", ylab="iris[,3]", pch=(15:17)[unclass(iris[,5])], col=rainbow(3)[unclass(iris[,5])], cex=1.5)
# 文字列の描画
names <- c(1:150)
text(iris[,1], iris[,3], names, cex=0.8, pos=4)
# 凡例
legend(8.1, 4, levels(iris[,5]), pch=15:17, col=rainbow(3), cex=1, bg='gray80')
dev.off()


※posについて
1234の順に、下左上右となります。
ここでは、右横に表示させるため「pos=4」にしました。
| きむ | 統計解析ソフトR | 10:43 | comments(0) | trackbacks(0) |
R 「関数textとlegend 4」
前回のおさらいとして、関数textで「」をプロット領域に書きます。
次に、プロット領域の外に凡例を書き込みます。
前回のプログラムに、太字の部分だけ追加しました。

# デバイス領域
png("110114_iris.png", height=700, width=500)
# 作図領域
par(mar=c(4,4,2,6),xpd=TRUE)

# プロット領域(散布図)
plot(iris[,1], iris[,3], xlab="iris[,1]", ylab="iris[,3]", pch=(15:17)[unclass(iris[,5])], col=rainbow(3)[unclass(iris[,5])], cex=1.5)
# 文字列の描画
names <- c(1:150)
text(iris[,1], iris[,3], names, cex=0.8)
text(iris[99,1], iris[99,3], "○", cex=6, col="red")
# 凡例
legend(8.1, 4, levels(iris[,5]), pch=15:17, col=rainbow(3), cex=1, bg='gray80')
dev.off()

| きむ | 統計解析ソフトR | 17:38 | comments(0) | trackbacks(0) |
R 「関数textとlegend 3」
集団から外れている個体番号を目で確認したい場合など、関数textを使うと便利です。昨日のプログラムに、太字の部分だけ追加しました。
# 散布図
png("110113_iris.png", height=600, width=400)
plot(iris[,1], iris[,3], xlab="iris[,1]", ylab="iris[,3]", pch=(15:17)[unclass(iris[,5])], col=rainbow(3)[unclass(iris[,5])], cex=1.5)
# 文字列の描画
names <- c(1:150)
text(iris[,1], iris[,3], names, cex=0.8)

# 凡例
legend(6.5, 2.5, levels(iris[,5]), pch=15:17, col=rainbow(3), cex=1.4)
dev.off()

| きむ | 統計解析ソフトR | 11:26 | comments(0) | trackbacks(0) |
R 「関数textとlegend 2」
まず、irisデータの散布図を描きます。
speciesデータ(5列目)ごとに、記号の形と色を変えます。
次に、関数legendで凡例を書き加えます。

# 読込
data(iris)
# 散布図
png("110112_iris.png", height=600, width=400)
plot(iris[,1], iris[,3], xlab="iris[,1]", ylab="iris[,3]", pch=(15:17)[unclass(iris[,5])], col=rainbow(3)[unclass(iris[,5])], cex=1.5)
# 凡例
legend(6.5, 2.5, levels(iris[,5]), pch=15:17, col=rainbow(3), cex=1.4)
dev.off()



| きむ | 統計解析ソフトR | 15:54 | comments(0) | trackbacks(0) |
R 「関数textとlegend 1」
関数textを使って、作図領域に文字列を描画する。
text(x, y, labels, ...)

関数legendを使って、凡例を追加する。
legend(x, y, legend, ...)

今週は、ふたつの関数を使用して、図の中に文字を入れていきたいと思います。
| きむ | 統計解析ソフトR | 19:11 | comments(0) | trackbacks(0) |
R「グラフィックス参考実例」
Rで作図する際、見やすい図にこだわり始めると、結構な時間がかかってしまいます。日頃からRwikiに大変お世話になっていますが、とっても素晴らしいページを発見したのでメモします。

グラフィックス参考実例集
http://www.okada.jp.org/RWiki/?%A5%B0%A5%E9%A5%D5%A5%A3%A5%C3%A5%AF%A5%B9%BB%B2%B9%CD%BC%C2%CE%E3%BD%B8

プログラムの実例が沢山載ってますワッ!
| きむ | 統計解析ソフトR | 20:21 | comments(0) | trackbacks(0) |
Rを使って図を描こう5
今日は、研究対象の集団の特性を見ていきたいと思います。
ここでの臨床データは、測定数値などの連続変数を対象にします。

臨床データのファイル名は「clinical.txt」にしました。
number Age HT BW BMI Waist HbA1c Adiponectin IL.6 Leptin
1 51 169.5 52.8 18.37784 72.5 5.2 19.6 1.9 1.3
2 54 167.0 59.2 21.22701 NA 6.3 8.4 22.6 1.2
3 55 170.5 59.4 20.43326 79.0 6.5 5.5 1.9 3.4
4 52 182.8 61.3 18.34459 70.2 5.5 6.9 0.9 2.2
:
:
単語と単語をタブで区切ったファイルです。
(上のデモでは、スペース区切りになってしまっているようです)

では、Rで作業を行います。
作業ディレクトリをファイル「clinical.txt」を置いたフォルダに変更します。
> setwd("/Users/document/study")
データフレーム「clinical」に格納します。
> clinical <- read.table("clinical.txt", header = TRUE)
前回と同じように、データフレーム「clinical」にデータファイルの中身がちゃんと入っているかどうか確認してください。

summary関数を使うと、データの統計的な要約を見るこでができます。
データフレーム名「clinical」を入力し、Enterキーを押します。
> summary(clinical)
Age HT
Min. :42.00   Min. :149.0
1st Qu.:54.00   1st Qu.:156.6
Median :60.00   Median :165.0 ‥‥
Mean :58.06   Mean :164.2
3rd Qu.:62.00   3rd Qu.:169.5
Max. :69.00   Max. :182.8

次に、臨床データのうちBMIのヒストグラムを描いてみます。
PNG形式の画像ファイル「BMI_hist.png」に出力します。
hist関数を使用します。BMIの列を指定するには、「clinical$BMI」と書いてください。
> png('BMI_hist.png')
> hist(clinical$BMI, breaks=10, main="Histogram of BMI", xlab="BMI")
> dev.off()

作業ディレクトリに「BMI_hist.png」ファイルが作成されたか確認してみましょう。


オプションの数値やラベルを変えて、お好みの図を作成してください。


お問い合わせページ
http://amelieff.jp/enquiry/index.html



| きむ | 統計解析ソフトR | 14:06 | comments(0) | trackbacks(0) |
Rを使って図を描こう4
昨日は、Rで図を描くために
1.データファイルの準備
2.データファイルの読み込み
の2ステップを行いました。

今日は、3ステップ目の「グラフの描画」です。

作図ドライバを起動して、PNG形式の画像ファイルに書き込んでいきます。
「''」の中でファイル名を「rs100.png」に指定しています。
png関数を使用します。
> png('rs100.png')

次に、boxplot関数で箱ひげ図を書いていきます。
発現量のような連続変数は縦軸、遺伝型のようなグループを表す変数は横軸になります。
> boxplot(expression~genotype, data=rs100)
最後に作図ドライバを閉じます。
> dev.off()
無事に、作業ディレクトリに「rs100.png」ファイルが作成されたか確認してみましょう。

ちょっとそっけないですね。

では、オプションを設定して見やすい図を作成してみます。
> png(filename="rs100deco.png", width = 400, height = 400)
> boxplot(expression~genotype, data=rs100, main="rs100", xlab="genotype", ylab="expression", cex.main=2, cex.axis=1.3)
>dev.off()


いかがでしょうか。
他にもいろいろ設定できますのでお試しください。

お問い合わせページ
http://amelieff.jp/enquiry/index.html


| きむ | 統計解析ソフトR | 13:10 | comments(0) | trackbacks(0) |
Rを使って図を描こう3
ここでは、Expression-QTL解析により示唆された複数のeQTL候補のうち、ひとつだけを選び出して、図を作成し保存すると想定して作業していきたいと思います。

Rで図を描くための簡単な手順は以下の通りです。
1.データファイルの準備
2.データファイルの読み込み
3.グラフの描画

1.データファイルを準備します。
ファイル名は「rs100.txt」にしました。
Name genotype expression
1 AB 41.11
2 AB 35.431565
3 AB 40.74666
4 AA 65.08612
5 AA 48.581745
6 BB 11.8881965
7 AA 56.56106
:
:

単語と単語をタブで区切ったファイルです。
(上のデモでは、スペース区切りになってしまっているようです)
私は秀丸エディタで作成しましたが、ウィンドウズのメモ帳でも作成できます。
一行目は、列固有のデータ名(変数名)を入力します。
空白には「NA」と入力してください。

2.データファイルを読み込みます。
Rでは特定の作業ディレクトリ(フォルダ)に対して、ファイルの出入力を行います。
まず、getwd関数で現在の作業ディレクトリを表示させます。
> getwd()
次に、setwd関数で作業ディレクトリを、データファイル「rs100.txt」を置いたフォルダに変更します。「""」の中でフォルダを指定しています。下では、documentフォルダ内のeQTLフォルダに変更しています。
> setwd("/Users/document/eQTL")

read.table関数でデータファイル「rs100.txt」を読み込み、データフレーム「rs100」に格納します。
> rs100 <- read.table("rs100.txt", header = TRUE)

では、データフレーム「rs100」にデータファイルの中身がちゃんと入っているかどうかチェックします。
データフレームの中身を表示させるには、データフレーム名「rs100」を入力し、Enterキーを押します。
> rs100
Name genotype expression
1 1 AB 41.11000
2 2 AB 35.43156
3 3 AB 40.74666
4 4 AA 65.08612
5 5 AA 48.58174
6 6 BB 11.88820
7 7 AA 56.56106
:
:
と表示されました。

大きいデータのチェックを行う場合はtable関数を用いると便利です。
変数ごとに、どのような値を持つかチェックすることができますが、ここでは割愛させていただきます。

次回は、3ステップ目の「グラフの描画」です。


お問い合わせページ
http://amelieff.jp/enquiry/index.html

| きむ | 統計解析ソフトR | 17:08 | comments(0) | trackbacks(0) |
Rを使って図を描こう2 -Rの準備-
RをCRAN(The Comprehensive R Archive Network)からダウンロードします。日本からですと、筑波大学のミラーサイトからダウンロードするのがよいと思います。OSによってダウンロードするファイルが違うのでご注意ください。

次に、インストールします。
Macintosh の場合は、http://aoki2.si.gunma-u.ac.jp/R/begin.html
Windows・Linux の場合は、http://stat.sm.u-tokai.ac.jp/~yama/R/install.html
に詳しく紹介されています。

Rを起動します。
Rのアイコンをダブルクリックするとウィンドウが開きます。

ウィンドウの一番下の
>

はプロンプトと呼ばれ、この後ろに指示を記述します。
さっそくRを電卓のように使ってみます。
> 1+2

と入力し、Enterキーを押すと
> 1+2
[1] 3

と計算結果を出力してくれます。
[1]は「1番目の要素」という意味です。
[ ]内には、行頭の要素が何番目の要素にあたるかが表示されます。

Rを終了するときは
> q()

と入力しEnterキーを押すと、「作業スペースを保存しますか」と質問されます。次回に引き継いで作業を行う場合は保存を選択してください。

では、明日はRにファイルを取り込む作業まで行います。


【参考】
・筑波大学のCRANミラーサイト
http://cran.md.tsukuba.ac.jp/index.html

・R による統計処理
http://aoki2.si.gunma-u.ac.jp/R/

・R - 統計解析とグラフィックスの環境
http://stat.sm.u-tokai.ac.jp/~yama/R/
| きむ | 統計解析ソフトR | 12:56 | comments(0) | trackbacks(0) |
Rを使って図を描こう1
キムです。

先週は、PLINKを用いたQTL解析やGWAS、Expression QTL解析を紹介してきました。

今日からは、ExpressionQTL解析で有意差が確認されたSNPのデータから、
Rを使って箱ひげ図を作成する作業をしていきたいと思います。

Rは統計計算とグラフィック機能を兼ね備えたフリーの統計解析ソフトウェアです。

Rによる統計解析については、私はオーム社から出版されている
「Rによる統計解析」という本を愛読させていただいております。
ryamadaの遺伝学・遺伝統計学メモというサイトも大変参考になります。

ここでは、解析結果の可視化に焦点を絞りお話していきたいと思います。

Rについてはまだまだ勉強中の身ですので、至らない点などございましたら、
ご指摘いただけますと幸いでございます。

では、次回へ続きます。


【参考】
・PLINK
http://pngu.mgh.harvard.edu/~purcell/plink/index.shtml

・ryamadaの遺伝学・遺伝統計学メモ
http://d.hatena.ne.jp/ryamada22/

| きむ | 統計解析ソフトR | 17:01 | comments(0) | trackbacks(0) |
   1234
567891011
12131415161718
19202122232425
262728293031 
<< March 2017 >>

このページの先頭へ