バイオインフォマティクスってかっこいい

完全wetな研究室でdryに興味が出てしまったB4 のBlogです。最近は仮想通貨、ブックメーカーとかにも興味がある悪い子です。

single-end read, paired-end readsのおはなし

きっかけ

paired-endのサンプルのoverlapを取り除く方法を検討してたらいろいろ整理されてる記事に巡り合えたので和訳しつつ理解を深める

single-end, paired-end

fastqファイルには、シーケンシングのされ方で2種類ある。UCSCなんかからfastqをダウンロードする際にも2種類ある。

  • single-end
  • paired-end

の2つだ。

片側からのみ配列が読まれた場合、それはsingle‐end readとなる。一方、両側から配列を読むことのできる技術も存在する。 そうした技術を用いて読まれると、2つの末端から読まれたpairのreadができる。こうしたものをpaired-end readsという。 GAIIx, HiSeq and MiSeqといったillumina社の機器では、paired-endが標準的なpracticeとなっている。

single-endの場合、length (L) がfragment length (F)より短いか確認する必要がある。そのような場合、その配列はDNAを使い果たしていることになる。標準的なilluminaのフラグメントライブラリーはF ~ 450bpを用いているが、これは一定の値ではなく、フラグメントごとに異なる。

paired-endの場合、Fが二つのリード長を合わせた値より十分に大きくなければならない。リードに挟まれた部分の大きさは多少の幅がある。

今年に入っての変化

しかし今年は、illuminaの体系は少し変わった。 まず、リードの長さはHiSeqで>150bp (GAIIxでも同様) MiSeqでは>250bpへと変わった。 これは、標準のライブラリサイズであるF~450がとても小さいものとなり、paired end readsがoverlapするようになることを意味する。 2つ目に、酵素を用いたNextera のライブラリ調整システムによって、以前のTruSeqシステムに比べ、広い範囲のFサイズのライブラリが作れるようになった。Nexteraを用いると、同じライブラリ中で100bp~900bpのF値を作ることができる。ゆえに、overlapするリード、overlapしないリードが同じライブラリ中で存在しうる。

paired-end readsで肝心なことは、実際には不可能な長さのシーケンスが可能にならないとしても、longer readの恩恵を受けることができる点だ。長さFのフラグメントから得られたpaired-end read(長さLの2本のリード)は一連の中間の配列がわからず、中間の長さも大まかにしかわからないことを除けば、長さFのsingle-end readのようなものだ。

しかし、多くのソフトウェアツールはoverlapするpairを与えると、誤った結果を出してしまう。そのため、overlapしうる場合は、長いsingle-end readに変えれば、多くのツールはより良い結果をより早く出すことができる。

Tools

PEAR (Paired-End Read Merger)

COPE (Connecting Overlapping Paired End reads)

SeqPrep

FLASH (Fast Length Adjustment of Short Reads to Improve Genome Assemblies)

fastq-join (part of ea-utils)

PANDAseq

stitch (now defunct, merged into PANDAseq)

mergePairs.py

bamファイルがpaired end かsingle end なのかを調べる

Rsamtools を使った判定法

> packageVersion("Rsamtools")
[1] ‘1.20.5’
library(Rsamtools)

# ファイルのPATHを指定
fl <- system.file("extdata", "ex1.bam", package="Rsamtools")

> quickBamFlagSummary(fl)
                                group |    nb of |    nb of | mean / max
                                   of |  records |   unique | records per
                              records | in group |   QNAMEs | unique QNAME
All records........................ A |     3307 |     1699 | 1.95 / 2
  o template has single segment.... S |        0 |        0 |   NA / NA
  o template has multiple segments. M |     3307 |     1699 | 1.95 / 2
      - first segment.............. F |     1654 |     1654 |    1 / 1
      - last segment............... L |     1653 |     1653 |    1 / 1
      - other segment.............. O |        0 |        0 |   NA / NA

Note that (S, M) is a partitioning of A, and (F, L, O) is a partitioning of M.
Indentation reflects this.

Details for group M:
  o record is mapped.............. M1 |     3271 |     1699 | 1.93 / 2
      - primary alignment......... M2 |     3271 |     1699 | 1.93 / 2
      - secondary alignment....... M3 |        0 |        0 |   NA / NA
  o record is unmapped............ M4 |       36 |       36 |    1 / 1

Details for group F:
  o record is mapped.............. F1 |     1641 |     1641 |    1 / 1
      - primary alignment......... F2 |     1641 |     1641 |    1 / 1
      - secondary alignment....... F3 |        0 |        0 |   NA / NA
  o record is unmapped............ F4 |       13 |       13 |    1 / 1

Details for group L:
  o record is mapped.............. L1 |     1630 |     1630 |    1 / 1
      - primary alignment......... L2 |     1630 |     1630 |    1 / 1
      - secondary alignment....... L3 |        0 |        0 |   NA / NA
  o record is unmapped............ L4 |       23 |       23 |    1 / 1

single segmentだけ:シングルエンド multiple segmentだけ:ペアエンド

とりあえず知りたい

testPairedEndBam(fl)
[1] TRUE

Matplotlib subplot の仕方

グラフを並べて表示したいとき、 2つの方法がある。ほかにもあるかもだけど2つの方法を知っている。

matplotlib.pyplot.subplotsを使う

fig, axes = plt.subplots(figsize=(10, 10), nrows=2, ncols=4, subplot_kw={'adjustable': 'box-forced'})

axes には nrows , ncols で指定したmatrixが入る。この場合、

array([[<matplotlib.axes._subplots.AxesSubplot object at 0x0000023A5A4D6F28>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x0000023A59495EF0>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x0000023A582FFFD0>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x0000023A59DBFE10>],
       [<matplotlib.axes._subplots.AxesSubplot object at 0x0000023A58305208>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x0000023A57F77AC8>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x0000023A59BABC18>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x0000023A58763908>]], dtype=object)

matplotlib.axes._subplots.AxesSubplot object なるものが作成される。

ここにどしどしグラフを埋め込むことができる。

axes[1,2].bar(left, height, width=0.02, yerr=0.1, ecolor="black", tick_label=label,capsize=10, align ='center')

f:id:kimoppy126:20180207191522p:plain

matplotlib.pyplot.subplot を使う

微妙に違う。複数形じゃなくなっている。

# 2×2の1枚目
plt.subplot(221)
# 2×2の2枚目
plt.subplot(222)
# 2×2の3枚目
plt.subplot(223)
# 2×2の4枚目
plt.subplot(224)
plt.ylim(ymax = 1.0, ymin = 0)
plt.bar(left, height, width=0.01, yerr=0.1, ecolor="black", tick_label=label,capsize=10, align ='center')

f:id:kimoppy126:20180207191736p:plain

以上。久しぶりの投稿でした。

Vagrantでのssh接続

ssh 接続にはパスワード認証方式と、公開鍵認証の二通りの接続方法がある。

公開鍵認証

  • vagrantでは、ゲストOSの初回起動時にホストOS側の鍵情報を自動で変更し、ゲストOSと暗号の再調整を行っている。
vagrant ssh

vagrant ssh-config

で置換後に使う秘密鍵ファイルの場所を取得し、

ssh -p [ポート番号] -l vagrant -i /path/to/private_key [IPアドレス]

で同じことができる。

もしくは、起動時に置き換えないように設定。Vagrantfileに以下の一行を追加する。

config.ssh.insert_key = false

複数のマシンで共通のsshキーを使用したい場合なんかはこうしないといけない。

ssh コマンドで秘密鍵ファイルの場所を~/.ssh/config 中に追記しておけば、明示的に指定する必要はなくなる。

vagrant ssh-config >> ~/.ssh/config
  • 鍵認証による ssh を行おうとすると ssh する元のサーバで公開鍵を生成し、ssh 先のサーバに持って行って.ssh/authorized_keysに追記する必要がある

公開鍵の生成

ssh-keygen -t rsa

公開鍵の.ssh/authorized_keys へのコピー

ssh-copy-id root@192.168.100.20

パスワード認証

  • パスワード認証の場合、Vagrantfile中でパスワードとユーザー名を設定できる。
  config.ssh.username = 'vagrant'  
  config.ssh.password = 'vagrant'  

Selenium備忘録

今友達からbed365 というサイのスクレイピングを任されている。ブックメーカーで有名なサイトだ。ここのデータを集めて解析したいらしい。

スクレイピングでは基本ChromeDriverを使って行っているのだが、久しぶりに動かしたら、以下のようなエラーが度々表示された。

[2084:15336:0122/010511.788:ERROR:process_metrics.cc(105)] NOT IMPLEMENTED

スクレイピング自体はできていたからそんなに気にしていなかったのだが、そろそろ気になってきたので調べてみた。

github.com

これだ。。

Seleniumの問題じゃなくChromeの問題らしい。v63でのバグだとか。。

このissueでは、v62に変更することが推奨されている。

古いChromeは以下から入手できる。
https://www.slimjet.com/chrome/google-chrome-old-version.php

v63アンインストール。

v62 インストール。


なおった

IPython データサイエンスクックブック memo (2)

2.4 Workflow using git branch

git stash

commitしていない変更の一時退避

git stash pop

変更を戻す

2.5 High reproductivce, interactive computing

  • ファイルの命名規則ディレクトリ構造を一貫性のあるものにする。
  • 全てのソフトウェアスタックの正確なバージョンを記録する
  • dill, Joblibの利用。JoblibはNumpyにも配したメモ化パターンが存在。
  • pararell, multiprocessing, Joblibなどの並列計算ライブラリを使った並列化の検討

2.6 High quality python code

  • assert の組み込み
  • クラスはなるべく避ける
  • 関数の引数にはキーワード引数を。

2.7 Unit test using nose

%%writefile -a [file]

でファイルに追記可能。

test_xxx.py モジュールをxxx.py モジュールに付随させる

テスト開始時:setup()
終了時:teardown()
テストを始める前と終了した後で環境が変わらないようにする。

nosetests

で実行

2.8 debug with IPython

script.pyっでインポートされているextscript.pyの20行目にブレークポイントを設定し、デバッカ制御下でscript.py実行

%run -d -b extscript.py:20 script

コマンドラインPythonスクリプトとして実行されるようになっている、GUIに統合されている、等の場合は

from IPython import embed
embed()

を組み込めば、そこでIPythonデバッガが起動する。

IPython データサイエンスクックブック memo (1)

1.1 Introduction

%%writefile

IPython magicコマンド。テキストファイルの作製。

インラインの数式は$...$ を使って記述。
独立した数式は $$...$$ を使って記述。

HTML()SVG()YouTubeVideo()

nbviewerを使ってIPythonで生成したJSONテキストを公開可能。

1.2 Data Analysis with IPython

df = pd.read_csv(url, index_col='Date', parse_dates=False, dayfirst=False)
  • pase_dates
    01/01/2013 → 2013-01-01
  • dayfirst
    2013-02-01 → 2013-01-02
df.describe()

サマリー表示

days = np.array(['Monday', 'Tuesday', 'Wednesday', 
                 'Thursday', 'Friday', 'Saturday', 
                 'Sunday'])
df['Weekday'] = days[df.index.weekday]

number → names

df_week.ix[days].plot(lw=3, figsize=(6,4));

-lw
line width - figsize figure size

.ix で行っているのは並べ替えのみ。

1.3 Numpy

%precision 3

ipython のmagicコマンド。数値出力を第3位までに。

n = 10000
x = [random.random() for _ in range(n)]
y = [random.random() for _ in range(n)]
%\time it [x[i] ; y[i] for i in range(n)]
1 loop, best of 3: 250 ms per loop
  • %timeit
    実行時間の計測。複数回実行、最速の実行時間を表示。
%timeit sum(x)  # pure Python
%timeit np.sum(xa)  # NumPy
The slowest run took 25.57 times longer than the fastest. This could mean that an intermediate result is being cached.
100 loops, best of 3: 6.69 ms per loop
The slowest run took 59.18 times longer than the fastest. This could mean that an intermediate result is being cached.
100 loops, best of 3: 1.53 ms per loop

numpy使うとこのくらい早くなる。

x[:1000,None]

Noneを使えば次元を増やせる。

1.4 custom magic command

from IPython.core.magic import (register_line_magic, 
                                register_cell_magic)
  • register_line_magic
    行magicコマンドの定義
  • register_cell_magic 列magicコマンドの定義
def load_ipython_extension(ipython):
    ipython.register_magic_function(function, 'cell')

この関数を組み込んでおけば実行時に読まれてmagicコマンドとして登録される。 IPython のInteractiveShellインスタンスを引数に受け取る。

cythonmagicrmagic

cython の関数をセル内で定義、インポート可能。

1.5 Interactive computing

定義したmagicコマンドの引数の指定

  • ipython実行時に指定
ipython --profile=cookbook --RandomMagics.text='Your number is {n}.' --RandomMagics.max=10 --RandomMagics.seed=1
  • 設定ファイルで設定 ~/.ipython/profile_default_ipython_config.py を編集
c.RandomMagics.text = 'random{n}'

Configurableクラス > 設定ファイル  いずれも継承可能

1.6 IPython Kernel

IPythonでカーネルを作ることができる。用途:IPythonのユーザーインターフェースを書き換える。

IPython version5.1.0 ではカーネルの起動の仕方がわからなかった。要調査。