アメリエフのブログ

バイオインフォマティクスの紹介と社員の日々
GRCh38
ヒトの新しいリファレンスゲノムGRCh38が公開されました。
GRCh37(≒hg19)と比較してどれくらい変わったのでしょうか?

まずは、各染色体の長さや、N(不明塩基)の数がGRCh37→GRCh38で
どのように変わったか調べてみました。



染色体長の変化
・各染色体につき最大30Mbpくらい増えたり減ったりしている
・ゲノム全体では7.4Mbp減少

N数の変化
・ゲノム全体で83Mb減少
・ゲノムに占めるNの割合は7.6%から4.9%に減少

GenomeRef:Announcing GRCh38によると、GRCh37で
セントロメア領域にあった大量のNが、GRCh38では減らせている
そうです。
また、1000ゲノムプロジェクトと一致しない領域やIndelが
修正されているそうです(1q21、10q11、chr9など)。

新しいリファレンスゲノムで新しい知見がたくさん見つかることを
期待しています。
| hat | バイオインフォマティクス | 16:47 | comments(0) | - |
HGVDについて (2) サンプル数編
前の記事
次の記事

Human Genetic Variation Databaseで公開している日本人ゲノムデータの、データの確からしさの判断に使えそうな項目を見てみます。

配布されている圧縮ファイルをダウンロードし、解凍しますと、READMEとタブ区切りテキストの2つのファイルができます。

READMEによると、日本人から検出された478,228個のSNPについて、1行につき1SNPの以下の情報が記載されています。
hgvdのカラム説明
■#Sample(サンプル数):データベースに含まれる1,208人のうち、N人でそのSNPのポジションをシーケンスできたという意味です。Altのアリルが検出されたサンプルの数ではありません。

このうち、#Sample、Mean_depthが、確からしさを判断するのに使えそうです。

まず#Sampleの分布を見てみます。
HGVDのサンプル数ヒストグラム
(※クリックで大きくなります)

まず左の図をご覧ください。
横軸が各SNPをシーケンスできたサンプル数、縦軸がその登場回数です。
200〜500サンプル、あるいはそれ以上でシーケンスできているSNPが多いようですが、案外サンプル数が少ないSNPも登録されているようです。
試しに、サンプル数が100より少ないSNPでヒストグラムを書いたのが右図です。40サンプル以下、その中でも5サンプルより少ないサンプルでしか読まれていないSNPが多数あるようです。

ちなみに、データベースのうち、100サンプル未満でしかシーケンスされていないSNPをすべて除くと、約95.0%のSNPが残ります。

次はMean_depthの分布を見てみます。



データベースの名前を「Human Genome Variation Database」と記載しておりましたが、これは別のデータベースの名前でした。
正しくは「Human Genetic Variation Database」でした。
ブログをご覧くださった皆様、並びに関係者様各位にご迷惑をおかけしてしまい申し訳ございません。
| kubo | バイオインフォマティクス | 15:58 | comments(2) | - |
レゴ
先日、社内研修がありました。

その中で「アメリエフが提供する未来をレゴブロックで作る」
というグループワークがあり、大変に盛り上がりました。

「社内」を作りこむ人、「社外との関わり」が気になる人、
なぜか一心に車を作る人など、各自の個性が出たように思います。

作品はこちらです。



この後、意味付けや鑑賞をして、皆で会社のビジョンを考えました。

なんか可愛いのもいますね。
| hat | 会社のこと | 15:56 | comments(0) | - |
社員研修に伴う休業のお知らせ
本日1月17日(金)12:00より、全社員研修のため営業をお休みをさせていただきます。

ご不便をお掛け致しますが、何卒ご了承くださいますよう宜しくお願い申し上げます。


| akb | 会社のこと | 09:59 | comments(0) | - |
bwa memの-Mオプション
bwaバージョン0.7が去年の2月にリリースされて、
ほぼ1年が経とうとしています。
もう0.7に切り替えた方も多いのではないでしょうか。

bwa 0.7では、memというコマンドが使えるようになりました。

memは従来のbwaswの後継で、ロングリード(〜1Mbp)に使える
コマンドですが、bwaswよりも速度・精度とも向上しているとのことです。

また、ショートリード(70-100bp)に対する精度もよいそうです。
ショートリードはこれまで「fastqごとにaln」→「samse/sampe」で
実行するのが主流で、中間ファイルができたりして面倒くさかったので、
mem一発でやれると便利ですね。

このように素晴らしいmemですが、1点注意があります。
memには-Mというオプションがありますが、bwaの結果を
PicardやGATKで処理する場合は、この-Mオプションを必ずつけて
実行してください。

$ bwa mem -M GENOME.fa 1.FASTQ 2.FASTQ > BAM

何故かと言うと、memはローカルアライメントを行い、
キメラリードを検出することができるようになりました。
キメラリードは、BAMのflagに0x800(supplementary alignment)が
立ちます。
ですが、PicardのMarkDuplicateはこの仕様に対応していないため、
このBAMを入力するとエラーになってしまうのです。

-Mをつけて実行したBAMでは、flagに0x800の代わりに従来の0x100
(secondary alignment)が立つため、Picardで動かすことができます。

試しに、同じデータを-Mありと-Mなしで実行した結果のBAMファイルを
比較したところ、本当にflagだけが1792(0x800=2048、0x100=256、2048-256=1792)違っていました。

キメラリードを検出できるのは良いので、早くPicardが0x800に
対応してくれるといいなあと思っています。
| hat | バイオインフォマティクス | 14:32 | comments(0) | - |
CircosでSelfChainを描く(3)
CircosでSelfChainを描く(1)
CircosでSelfChainを描く(2)
のつづきです。

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


データファイルも設定ファイルも用意できましたので
いよいよ「(3)Circos実行」を行います。

次のコマンドを実行すると、circos.pngというファイルができます。
$ circos -conf test.conf



Circosはこれ以外にもいろいろな形状のトラックを描くことが
できますので、またご紹介したいと思います。
| hat | バイオインフォマティクス | 14:36 | comments(0) | - |
insert sizeを求める
先日、構造多型を解析するツールについてご紹介しました。
(Paired-end/Split-read/Complex)
これらのツールの一部は、解析の際、シーケンスデータのinsert sizeの入力が必要です。
自ら実験して得たデータならinsert sizeはわかると思いますが、公共のデータベースの公開データはリード長は記載されていてもinsert sizeは不明な場合が多いです。

そこでinsert sizeを計算してくれるツールを調べていたところ、PicardsのCollectInsertSizeMetricsがありました。
BAM/SAMファイルからinsert sizeのヒストグラムを描画してくれるツールです。
様々なパラメータがありますが、必須のパラメータは3つだけです。

1. 入力するBAM/SAM
2. insert sizeを出力するファイル
3. ヒストグラムを出力するファイル

コマンドは以下のようになります(ファイル名・バージョンは任意です)。
$ java -jar [path]/picard-tools-1.75/CollectInsertSizeMetrics.jar INPUT=sample.bam OUTPUT=output.txt HISTOGRAM_FILE=hist.pdf

OUTPUTパラメータで指定したファイル(ここではoutput.txt)には、insert sizeの平均や中央値、最大値や標準偏差など、insert sizeに関する統計量が書かれています。

HISTOGRAM_FILEで指定したファイル(ここではhist.pdf)には、下図のようなヒストグラムが生成します。
insert size distribution

また、このヒストグラムを生成したデータでは、統計量はmedian insert sizeが200bp、mean insert sizeが213.9bpと計算されました。
ヒストグラムと統計量を見比べ、このデータのinsert sizeは200 bpだと判断しました。

詳細については公式のマニュアル(こちら)をご参照ください。
| kubo | 次世代シーケンサー解析 | 15:01 | comments(0) | - |
CircosでSelfChainを描く(2)
CircosでSelfChainを描く(1)のつづきです。

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


今回は「(2)設定ファイル作成」を行います。
以下のような内容のtest.confという名前のファイルを作成します。


<<include colors_fonts_patterns.conf>>

<image>
<<include etc/image.conf>>
</image>

karyotype = karyotype.human.txt

<ideogram>
<spacing>
default = 30u
<pairwise /hs/ /hs/ >
spacing = 0.5r
</pairwise>
</spacing>
thickness = 25p
stroke_thickness = 1
stroke_color = black
fill = yes
fill_color = black
radius = 0.3r
show_label = yes
label_font = default
label_radius = dims(ideogram,radius_outer) + 10p
label_size = 28p
band_stroke_thickness = 0
band_stroke_color= black
band_transparency= 4
</ideogram>

<links>
radius = 0.7r
crest = 0.5
thickness = 1
color = red
bezier_radius = 0r
bezier_radius_purity = 0.3
<link>
file = data.txt
<rules>
<rule>
condition = var(intrachr)
color = blue
radius = 0.7r
bezier_radius = 1.7r
</rule>
</rules>
</link>
</links>

<<include etc/housekeeping.conf>>


以下のような設定を行っています。

1. 外円にヒトの染色体領域を描画する。
2. その内側にSelfChainのデータファイルを読み込んで描画する。
_a. 同一染色体上の関係は青線で外側に描画する。
_b. 異なる染色体間の関係は赤線で内側に描画する。


詳しい説明はCircosのチュートリアルをご覧ください。

次回は描画を行います。
| hat | バイオインフォマティクス | 14:44 | comments(0) | - |
ヒト Exome データ解析受託キャンペーン
ヒト Exomeデータ解析の受託キャンペーンを開始いたします!

ご注文受付期間:2014.1.14〜3.31ご発注分まで

期間限定の特別価格でのご提供となりますので、Exomeデータ解析の外注をご検討されている方は、ぜひ一度弊社ホームページをご覧ください!

ご不明な点がございましたらお気軽にお問い合わせください。
たくさんのお問い合わせお待ちしております!
| akb | 次世代シーケンサー解析 | 14:46 | comments(0) | - |
   1234
567891011
12131415161718
19202122232425
262728293031 
<< January 2014 >>

このページの先頭へ