アメリエフのブログ

バイオインフォマティクスの紹介と社員の日々
年末年始休業のお知らせ
平素よりアメリエフ株式会社のサービスをご利用頂きまして、誠にありがとうございます。年末年始の休業および営業日についてお知らせいたします。

【休業期間】
2013年12月28日(土)〜2014年1月5日(日)まで
【営業開始】
2013年は1月6日(月)より通常営業

年末年始のお問い合わせにつきましては、6日以降に順次対応させていただきます。
ご不便をおかけ致しますが、何卒ご了承くださいますよう宜しくお願いいたします。

また、この一年のご愛顧に心より感謝とお礼を申し上げます。
来年も皆様にとって、すばらしい年となることを社員一同お祈り申し上げます。

アメリエフ株式会社 社員一同
| akb | 会社のこと | 17:42 | comments(0) | - |
CircosでSelfChainを描く(1)
Circosを使ってみようでご紹介したCircosを使って、
以下のようなヒトのChainSelfの図を描いてみたいと思います。



手順
(1)データファイル用意
(2)設定ファイル作成
(3)Circos実行

今回は「(1)データファイル用意」を行います。

ChainSelfは、ヒトゲノムを自分自身にアライメントして、
類似領域を探したデータです。

データはUCSCから入手可能です。
UCSC hgdownloadから、chainSelf.txt.gzをダウンロードして
解凍します。

ChainSelfのように、ゲノム上の2領域間をつなぐデータは
CircosではLinkトラックという形式で表します。

Linkトラックとして読み込むために、ダウンロードしたファイルから
染色体1、スタート位置1、エンド位置1、
染色体2、スタート位置2、エンド位置2
だけ抜き出したファイルを作成します。

また、Circosの表記に合わせ、染色体名中のchrをhsに置換します。
さらに、全データだと数が多いため、スコアが10Mより大きい箇所にだけ
絞り込みました。

以下のコマンドを実行すると、上記の変換を行うことができます。
$ awk '$2>10000000{print $3, $5, $6, $7, $10, $11}' chainSelf.txt | sed -e "s/chr/hs/g" > data.txt

これでデータファイルの作成は完了です。

次回は、Circosの設定ファイルを作成したいと思います。
| hat | バイオインフォマティクス | 15:09 | comments(0) | - |
HGVDについて (1)
先月の11月12日、日本人のゲノム情報データベースHuman Genetic Variation Databaseが公開されました。このことは一般紙でも報道されました。
以前の連載でも、ゲノム情報のデータベース、特に人種や民族別のデータベースが、高スループットの遺伝子解析機器を用いた解析に有用だと少し触れました。ぜひ、今後このデータベースが日本人の遺伝子研究に活用され、様々な病気の治療に役立って欲しいと思います。

このデータベース上のSNPの半分以上は、dbSNPなどこれまでの既知変異データベースに含まれていない日本人特異的なものだということで、日本人の遺伝性疾患の研究に、日本人特有のデータベースがいかに重要かということがわかります。
一方で、アレル頻度が非常に低いSNPも大量に含まれているとも言われています。


次回では、せっかくなので、実際に中身を覗いてみます。
ホームページによると、このデータベースは日本人1208名のエクソームシーケンシングデータから発見されたアリルやジェノタイプの頻度が載っています。
他のゲノム情報データベースと比べてみますと、HapMapではサンプル270人、1000人ゲノムプロジェクトではサンプル2577人ということですので、日本人集団だけで1208人というのは、大規模に感じられます。
そのほか、データベースには信頼性の指標となりそうな、SNPをカバーできたサンプル人数や平均カバレージなどの情報が含まれていると説明されています。
これらの数字を見てみたいと思います。


HGVDについて (2) サンプル数編
HGVDについて (3) 平均depth編


データベースの名前を「Human Genome Variation Database」と記載しておりましたが、これは別のデータベースの名前でした。
正しくは「Human Genetic Variation Database」でした。
ブログをご覧くださった皆様、並びに関係者様各位にご迷惑をおかけしてしまい申し訳ございません。
| kubo | バイオインフォマティクス | 14:22 | 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) | - |
コーヒーマシン
ネスレのコーヒーマシン、バリスタ(左)とドルチェグスト(右)が
我が社にきました。



バリスタはインスタントコーヒーの粉を使うタイプ、
ドルチェグストはカプセルを入れ替えるタイプで、
おいしいコーヒーを気軽に飲むことができるので
カフェイン中毒に拍車がかかっています。

テイクアウト可のカフェが近所に無数にあるのですが、
わざわざ買いに行くこともなくなりました。
たぶん全然違うのでしょうが、シーケンサーを買われた先生の
気持ちが少しだけわかったような気がしました。
| hat | よもやま話 | 14:30 | comments(0) | - |
Phased vs. Unphased
VCFファイルの一番右側には、以下のようにジェノタイプ情報が記されています。

----------------------------------------
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA00001 NA00002 NA00003
20 14370 rs6054257 G A 29 PASS NS=3;DP=14;AF=0.5;DB;H2 GT:GQ:DP:HQ 0|0:48:1:51,51 1|0:48:8:51,51 1/1:43:5:.,.
----------------------------------------

このジェノタイプ情報が、パイプ(|)で区切られている場合とスラッシュ(/)で区切られている場合があります。例えば、1|0とか0/1のように。
今までこの違いを意識したことはなかったのですが、どうやらちゃんとしたルールがあるようです。

"|"で区切られているデータは、フェーズされている(Phasedな)データ
"/"で区切られているデータは、フェーズされていない(Unphasedな)データ
といいます。

詳細は1000人ゲノムプロジェクトのWebページで解説されています。
http://www.1000genomes.org/node/101
| heshi | 次世代シーケンサー解析 | 17:33 | comments(0) | - |
DeNovoGear
tokunagaです。

有言実行です。
第2回使ってみたシリーズです。

今回ご紹介するのは

DeNovoGear

サイト
http://denovogear.weebly.com/index.html
論文
http://www.nature.com/nmeth/journal/v10/n10/full/nmeth.2611.html


trioデータからDeNovoのmutationやIndelを検出するツールです。
likelihood-based errorというモデルを使用し、
false positiveを減らしているようです。

Inputファイルはpedファイルとbcfファイルです。

pedファイルについては過去のブログ記事をご参照ください。

bcfファイルは簡単に言うとVCFファイルのバイナリ版で、
検出したvariantの情報が記載されています。
詳しい説明はこちらに。
http://samtools.sourceforge.net/mpileup.shtml

このDeNovoGearでは

1.trioデータから新規変異を検出する(X染色体については別解析)
2.pairデータから新規変異を検出する
3.Phasing denovo mutations

といったパターンの解析が出来るようです。

テストデータもありましたので、1のtrioデータで実行してみました。
実行コマンドは以下です。
$ ./denovogear dnm auto --bcf sample_CEU.bcf --ped ¥
 sample_CEU.ped

上記を実行しますとCHILD_IDのサンプルに特異的な変異を検出してくれます。
outputの指定なしだと以下のように結果が出力されました。


Denovogear0.5.3
Created SNP lookup table
First mrate: 1 last: 1
First code: 6 last: 6
First target string: AA/AA/AA last: TT/TT/TT
First tref: 0.0002388 last: 0.99301

Created indel lookup table First code: 6 last: 6
First target string: RR/RR/RR last: DD/DD/DD
First prior: 0.05 last: 0.114

Created paired lookup table
First target string: AA/AA last: TT/TT
First prior 1 last: 1
DENOVO-SNP CHILD_ID: NA12878_vald-sorted.bam.bam chr: 2 pos: 214668360 ref: G alt: A maxlike_null: 3.95324e-12 pp_null: 0.000399465 tgt_null(child/mom/dad): GG/GG/GG snpcode: 1 code: 6 maxlike_dnm: 9.9301e-09 pp_dnm: 0.999601 tgt_dnm(child/mom/dad): AG/GG/GG lookup: 4 flag: 0 READ_DEPTH child: 48 dad: 76 mom: 34 MAPPING_QUALITY child: 59 dad: 59 mom: 59

Total number of SNP sites interrogated: 36
Total number of SNP sites passing read-depth filters: 36
Total number of INDEL sites interrogated: 1
Total number of INDEL sites passing read-depth filters: 1
Total number of Paired sample sites interrogated: 0
Total number of Paired sample sites passing read-depth filters: 0
Done !


--out_vcfオプションで
vcfファイルを出力することも可能です。

ご興味のある方は試されてみてはいかがでしょうか?
| tokunaga | バイオインフォマティクス | 17:46 | comments(0) | - |
第28回バイオインフォマティクス勉強会「NGSデータ解析サーバ体験セミナー」満員御礼
第28回バイオインフォマティクス勉強会「NGSデータ解析サーバ体験セミナー」はお申込みが定員に達しましたのでキャンセル待ちとさせていただいております。

たくさんのお申込み、ありがとうございました!
| akb | 勉強会 | 17:10 | 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) | - |
Pythonでグラフを書こう その2
Pythonでグラフを書こう、の第2回です。
前回の記事では、グラフの描画パッケージであるmatplotlibのインストールを行いました。今回は、簡単なグラフを書いてみたいと思います。

matplotlibのウェブサイトには豊富なexampleが表示されていまして、任意のグラフをクリックするとそのソースが表示されます。それを、そのままコピー&ペーストするだけでグラフを書くことができます。それでは、Pie and polar charts の左側のグラフをクリックしてみましょう。するとそのグラフのページが開きます。このページの下方にソースが書かれていますので、それをそのままコピー&ペーストして、pie_chart.pyとして保存し、自分の環境(ここではCentOS)で実行してみましょう。

$ python pie_chart.py

GUIのデスクトップ上で実行した場合は、そのまま下記のような円グラフが表示されたと思います。


しかし、コマンドラインのみの環境(外部からssh等で接続している場合等)では、実行してもエラーになります。デスクトップ上に出力できないためです。そこで以下のソースを追記することで、グラフの出力先をファイルに変更することができます。

【プログラムの先頭に以下の2行追記】
import matplotlib
matplotlib.use('Agg')

【プログラムの最後に以下の1行追記】
plt.savefig("hoge.png")


この内容で実行すると、PNGフォーマットでグラフを保存できます。
このプログラムの詳細については、次回お話したいと思います。
| deda | システム | 15:03 | comments(0) | - |
もし日本がひとつの会社だったら

もし日本がひとつの会社だったら〜「日本経営」進化への提言〜

田村耕太郎氏の著書『もし日本がひとつの会社だったら〜「日本経営」進化への提言〜』では、日本が今解決すべき問題点を独自の視点で指摘しています。その中で特に印象が強かったのは、海外への独自技術のセールスを強化することです。

本書では、2009年に日本がUAEの原子力発電所建設事業を受注する機会を逃した例を挙げ、政府のセールス・マーケティング力の強化を訴えています。

今年決定した、2020年東京オリンピック招致のスピーチやPR活動を思い返してみると、政府や都の本気を見た気がします。本書が発行された2010年と比べると、政党も変わり、政府の世界に対するマーケティング力も向上したのではないでしょうか。

本書ではさらに、著者の出身地である鳥取県産のスイカをドバイ・アブダビに売り込んだ自身の経験を例に挙げ、広く知られていないローカルな良いものを見つけ出してそれをグローバルに発信する、という経営者的センスとグローバルな視点を、政治に取り入れるべきだと主張しています。

この点については、個人的に非常に共感しますが、現状では改善の余地があるように感じます。

出版年の2010年から日本も世界も大きく変わったので、若干内容が古いかもしれませんが、筆者の熱い想いが伝わる内容です。
| heshi | 書籍の紹介 | 14:45 | comments(0) | - |
新しい環境
弊社が神田に移転しまして一週間が過ぎました。

以前の環境と今の環境の違いは以下の点でしょうか。

・オフィスが広い
・オフィス周辺にリーズナブルな飲食店が多い
・ランチのお弁当が格安

…食べ物の話ばかりというツッコミはさて置き、新しいオフィスで心機一転、いままで以上に業務を加速させていきたいと考えています。

しかし、毎朝11時ぐらいに、某ビールのCM音楽を流した車が走っていくのが何とも言えません。呑みたくなってしまいますね。

それでは、detでした。
| deda | よもやま話 | 14:07 | comments(0) | - |
「第36回日本分子生物学会年会」出展のお知らせ
アメリエフ株式会社は、2013年12月3日(水)〜6日(金)に神戸国際会議場で開催される、「第36回日本分子生物学会年会」に出展いたします。皆さまのご来場をお待ち申し上げます。

■ 学会名
・第36回日本分子生物学会年会
 ホームページ:http://www.aeplan.co.jp/mbsj2013/

■ 開催期間
・2013年12月3日(水)〜6日(金)
※ブースの展示は5日までとなります

■ 会場
・神戸国際会議場(兵庫県神戸市中央区港島中町6-9-1)

■交通アクセス
・神戸新交通ポートアイランド線 市民広場駅 下車すぐ
※詳細はこちらをご参照ください。

■ 展示内容
受託解析サービスのご紹介
 超高速シーケンスデータ解析の実例
バイオインフォマティクス・スクールのご紹介
 Linux基礎
 R基礎
 Perl基礎
 NGS基礎(Resequence、RNA-seq、ChIP-seq他)
解析レポートを配布
 「トリオの変異検出解析」
バイオインフォマティクス・勉強会のご紹介
バイオインフォマティクス・スクール春季奨学生募集のご案内
| akb | 学会出展 | 15:10 | comments(0) | - |
1234567
891011121314
15161718192021
22232425262728
293031    
<< December 2013 >>

このページの先頭へ