アメリエフのブログ

バイオインフォマティクスの紹介と社員の日々
重複値の削除や抽出
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) | - |
Excelの1セルにコンマ区切り文字列を入れる
解析結果を出力する際、複数の数値をコンマで結合して羅列したいことがあります。

例えばBEDフォーマットでは、エクソンの長さや開始位置を示すblockSizesやblockStartsには「468,69,147,159」や「0,608,1434,2245」のような値を入れることになっています。
※BEDフォーマットって何?という方は、過去記事「BEDフォーマット完全解説」をご覧ください。

このようなファイルをExcelで開くと、コンマが桁区切りと認識され、1つの数値に変換されてしまうことがあります。

下の図の赤で囲ったセルには「152,159,198,136,456」という文字列が入っています。
「152」「159」「198」「136」「456」という5つの数値をコンマで結合した文字列を入れたかったのに、「152159198136456」という数値として認識されてしまいました。


対策として、Excelにファイルを読み込むときにデータ型を明示的に文字列として指定すれば良いのですが、毎回やるのは少し面倒くさいですね。

もっと簡単な対策として「152, 159, 198, 136, 456」のようにコンマの後に半角空白を入れるのがおすすめです。これで勝手に変換されなくなります。

もしくは「152,159,198,136,456,」のように、最後にもう一つコンマをつけてもいいです。UCSC GenomeBrowserのBED12のblockSizesやblockStartsには最後にコンマがついています。
| hat | バイオインフォマティクス | 14:53 | comments(0) | - |
Vimで不可視文字を表示させる方法
こんにちは、まだ都営三田線に乗ったことがない久保(kubor)です。

先日CentOS6.7がリリースされましたね。
個人的に、Vim7.4に対応したのが嬉しいです。
嬉しいので、VimのTipsを書かせてください。

「Vimで不可視文字を表示させる方法」です。

Vimに限らず、テキストエディタでコーディングしているとタブやスペースが見えず、不安になると思います。怖いですよね。

Vimでは、~/.vimrcに以下の2行を書いて、表示させることが出来ます。
set list
set listchars=tab:»-,trail:-,nbsp:%,eol:↲

set listで、不可視文字を表示させます。
set listcharsで、表示させる(置き換える)文字を設定します。
表示文字の設定はお好みでどうぞ。

listcharsについて、項目の対応は以下です。
他にも有りますが基本的にはこの4つで十分だと思います。

  • tab:タブ文字を表示
  • trail:行末のスペースを表示
  • nbsp:ノーブレークスペースを表示
  • eol:改行を表示



余分なスペースもわかるし、タブとスペースの区別も出来て便利です。

ちなみに、一時的にこれらを非表示にしたいときは、
コマンドで対応するといいと思います。
:set nolist
| kubor | バイオインフォマティクス | 14:43 | 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) | - |
「次世代シークエンサ(NGS)ハンズオン講習会」参加者追加募集のお知らせ
 昨年に引き続き、「平成27年度NGSハンズオン講習会」にて一部の講義でアメリエフ社員が講師を勤めさせていただいております。
 この度、B日程(2015年8月26日(水)〜8月28日(金)の3日間)につきまして、追加で受講の申込を受付けることとなりましたので、お知らせいたします。

次世代シークエンサ(NGS)ハンズオン講習会

○受付期間: 〜8月18日(火)12時00分
○受講定員:30名程度(応募者多数の場合は、途中で募集締切となります)
○主催:
科学技術振興機構バイオサイエンスデータベースセンター(NBDC)
東京大学大学院農学生命科学研究科アグリバイオインフォマティクス教育研究ユニット
○お申し込みはこちらから:
http://biosciencedbc.jp/human/human-resources/workshop/h27-2

 8月初旬のA日程では、「NGS解析」について約90名の方にご受講いただきました。日程の都合でA日程でのご参加が難しかった方や、ご興味のある方は是非ご参加いただければと存じます。
| ymm | 会社のこと | 10:29 | comments(0) | - |
四半期研修
アメリエフでは四半期ごとに社内研修を行います。
7月中に行事などが立て込んでいたため時期がずれ込んでしまいましたが、上半期が終わったので八月の頭に研修を行いました。
スイカを食べながら、これまでを振り返ったり、今後の見通しや予定の話しました。

今年初のスイカだという人も多く、最初は「少し多いんじゃないか」という声もありましたが、気づけば一個丸ごと食べつくしていました。
食べる際に種を飲み込むか飲み込まないかで、議論があったりもしましたが、みなさんはスイカの種は飲み込む派、それとも飲み込まない派ですか?
| kubo | 会社のこと | 15:57 | comments(0) | - |
NGSハンズオン講習会
7月22日(水) から次世代シーケンサ(NGS)ハンズオン講習会が開催されています。

本講習会は、バイオサイエンスデータベースセンター(NBDC)の主催で、NGS解析を学びたい方を対象に行われており、弊社からは講師として2名、TAとして2名が参加しています。
今週までに、スクリプト言語として以下の講義を担当しました。
7/24(金) シェルスクリプト
7/27(月) Perl入門
7/28(火) Python入門

今週はNGS解析の講義が行われています。
8/3(月) NGS解析基礎
8/4(火) Reseq、変異解析
8/5(水) RNA-seq
NGSハンズオン
NGS解析については、8/26〜8/28に2回目の講義が行われます。

弊社では受託解析だけでなくバイオインフォマティクス勉強会や出張トレーニングも行っておりますので、今回のハンズオン講習会に参加できなかった方や参加されている方でもっと勉強したいという方は、ぜひご利用ください。
| onouek | スクール | 14:42 | 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) | - |
      1
2345678
9101112131415
16171819202122
23242526272829
3031     
<< August 2015 >>

このページの先頭へ