Aide-memoire

Bio-infomaticsで困った事の備忘録的なサムシング

自分が困った事をまとめた備忘録的なサムシング

kallisto/sleuth使いたい人生だった v2

昔Kallistoの使い方の記事を書いたんですけど，だいぶバージョンが古い奴だったので

## インストール

```
% wget https://github.com/pachterlab/kallisto/releases/download/v0.43.0/kallisto_linux-v0.43.0.tar.gz
% tar zxfv kallisto_linux-v0.43.0.tar.gz
% cd kallisto_linux-v0.43.0
% ./kallisto
kallisto 0.43.0

Usage: kallisto <CMD> [arguments] ..

Where <CMD> can be one of:

index Builds a kallisto index
quant Runs the quantification algorithm
pseudo Runs the pseudoalignment step
h5dump Converts HDF5-formatted results to plaintext
version Prints version information
cite Prints citation information

Running kallisto <CMD> without arguments prints usage information for <CMD>
```

僕は~/opt/binにすべてのソフトを入れて，パスを通しているので，

```
% cp ./kallisto ~/opt/bin
% cd ~
% kallisto
kallisto 0.43.0

Usage: kallisto <CMD> [arguments] ..

Where <CMD> can be one of:

index Builds a kallisto index
quant Runs the quantification algorithm
pseudo Runs the pseudoalignment step
h5dump Converts HDF5-formatted results to plaintext
version Prints version information
cite Prints citation information

Running kallisto <CMD> without arguments prints usage information for <CMD>
```

となります．

## 実際に動かす
例えば，とある生物全CDS配列 (transcripts.fa)と，RNA-Seqデータ (sample_R1.fq, sample_R2.fq)があるとします．

まずはkallistoのindexを作成します．

```
% kallisto index
kallisto 0.43.0
Builds a kallisto index

Usage: kallisto index [arguments] FASTA-files

Required argument:
-i, --index=STRING Filename for the kallisto index to be constructed

Optional argument:
-k, --kmer-size=INT k-mer (odd) length (default: 31, max value: 31)
--make-unique Replace repeated target names with unique names

% kallisto index -i transcripts.fa.kallisto transcripts.fa
```

こちらにpaired endなデータをpseudoalignmentしてTPMを計算します．

```
% kallisto quant
kallisto 0.43.0
Computes equivalence classes for reads and quantifies abundances

Usage: kallisto quant [arguments] FASTQ-files

Required arguments:
-i, --index=STRING Filename for the kallisto index to be used for
quantification
-o, --output-dir=STRING Directory to write output to

Optional arguments:
--bias Perform sequence based bias correction
-b, --bootstrap-samples=INT Number of bootstrap samples (default: 0)
--seed=INT Seed for the bootstrap sampling (default: 42)
--plaintext Output plaintext instead of HDF5
--single Quantify single-end reads
--fr-stranded Strand specific reads, first read forward
--rf-stranded Strand specific reads, first read reverse
-l, --fragment-length=DOUBLE Estimated average fragment length
-s, --sd=DOUBLE Estimated standard deviation of fragment length
(default: value is estimated from the input data)
-t, --threads=INT Number of threads to use (default: 1)
--pseudobam Output pseudoalignments in SAM format to stdout

%kallisto quant -i transcripts.fa.kallisto -o sample.kallisto --bias -b 100 -t 32 sample_R1.fq sample_R2.fq
```

-t でコア数を指定しますが、基本的に1分程度で終ります．

 

 

 

=========

とろあえずここまで。。。

後で続き書きます．

 

 

 

2016年

今年も残る事数時間，ちょっとだけ今年の振り返り．

 

 

しょっぱなは学部の卒論を書くために奔走したりでいきなりクライマックス．気づいたらボスからちょっと手伝わないかと言われ解析したら論文の共著に (

Genome sequencing of a single tardigrade Hypsibius dujardini individual : Scientific Data)．

中盤はちょっと中だるみしてしまって解析がうまく行かなかったり実験も失敗したりで，色んな事が重なりなかなか闇の深い時代．

終盤は今書いている論文のマラソン40キロ地点ぐらいで，学部時代からやってきた解析が無駄にならなかったと思っていたり．

 

なんか特に良い事はなかったけど，悪い事はたっぷりとまぁ平均ぐらいの年だった気がする．来年は論文のリバイスがしょっぱなで待っていたり，中盤は博士進学とかがあったり，修論が最後に控えていたりとなかなか山ばっかりで，その中でもう一本ぐらい論文掛けたらいいなぁって．

 

きっと来年は今年よりよい年でありますように．

 

You're beautiful, You're one of the kind, Smile More.

 

-abs

 

 

PS. 

今年は一個も記事を書かなかったので来年は最低5個ぐらい掛けるようになりたいです．(願い)

 

kallisto/sleuth使いたい人生だった

[更新版を書いています]

abs-0110.hatenablog.com

 

 

自分のメモ用と，あとは要望があったので

 

最近インフォ界隈で発現量解析のソフトで盛り上がりを見せているソフトの[kallisto]とそのパイプライン下流の発現変動の有意差検定を行う[sleuth]．早さと正確さが売りになっているみたいです．

patcherラボのブログ記事を見るとそれぞれのアルゴリズムがわかりますkallisto sleuth

 

僕なりの理解ですと，kallistoではトランスリプトームをk-mer (デフォだと31)にわけてde Bruijn graphを作って，それに対してfastqを当ててる感じです．有名なtophat/bowtieなどは全bpをアライメントしていくので非常に遅くなってしまうのですが，kallistoはこれらのアルゴリズムと比べて爆発的に早くなっています．そのため，なんとbootstrapを使えるようになったのです．

最近うちのラボでもこのkallistoが流行っているんですけど，このパイプライン下流にedgeR/DESeqを使うとあんまりよくないという記事を見たので，せっかくだからpatcherラボが作ったkallisto用の解析ソフトのsleuthも使ってみたいという感じで個々数日いじくりまわしていました．このsleuthもcuffdiffと比べて爆発的に早く，あくまでも個人的な感覚ですけど約10倍近く早いです．アルゴリズム的には，RNA-Seqのfastqをkallistoでトランスクリプトームにマップして得られるexpected_countsから(bootstrapの結果を使って)統計的にtechnical noiseとbiological noise を除いた各transcriptの真のcount値を算出して，(正規分布に乗っ取っているはずなので)そこから有意差を検定しています．

 

そんな二つのソフトのインストールから使い方をとりあえず載せておきます．

 

Kallisto

まだ生まれたてのソフトなので (最新がver0.42.3)かなり更新が早いです．

Downloadより最新版のbundle，もしくはsource codeをDLしてください．なおlinux機だと数verが無いので，野良buildしなければ行けません．ですが，僕みたいにrequirement (gcc ver > 4.8; zlib; cmake; HDF5; 僕はgccが4.4でした)が入らない場合は古いバージョンで我慢してください．

  Bundle (linux; version 0.42.1)

% wget 

https://github.com/pachterlab/kallisto/releases/download/v0.42.1/kallisto_linux-v0.42.1.tar.gz

% tar zxfv https://github.com/pachterlab/kallisto/releases/download/v0.42.1/kallisto_linux-v0.42.1.tar.gz

% cd kallisto_linux-v0.42.1

## ソフトは統一して~/opt/binに入れています

% cp kallisto ~/opt/bin

 

## Transcriptome index作成

% kallisto index -i <index_name> <transcript_fasta>

# Paired end

% kallisto quant -i <index_name> -o <output_dir> -b <bootstrap_num> R1.fq R2.fq

# Single read

% kallisto quant -i <index_name> -o <output_dir> -b <bootstrap_num> --single -l <fragment_length>

 

single endですと--singleと-l を入れないと行けません．fragment lengthはシークエンスライブラリーを作成した人に聞いてください．以下使う可能性が高いオプションです．なおこの先のsleuthは-bを使います (制作者は大体100あれば大丈夫という風に言っています)．

-t      　   thread数の指定 (v0.42.2以降)
--pseudobam  alignmentのSAMでの出力
--bias       bias補正を行う
-b           bootstrap数

 

これを前サンプルに対して行い(僕はshell scriptで全部bgで一気に流しています)，一つのdirectoryにまとめておきます．このまとめておくのがsleuthに必要です．ここでは~/k_resultsとします．

 

sleuth

kallistoより遅れて出されているんですけど同様に更新が早いのでupdateに気をつけてください．現状v0.27.3です．

とりあえずインストールからです．

source("http://bioconductor.org/biocLite.R")
biocLite("rhdf5")
install.packages("devtools")
devtools::install_github("pachterlab/sleuth")

僕はこのdevtoolsをインストールするときにgit2rが入らなくて結構苦労しました，結局libiconvを再びインストールしたら通りました．

 

ここから解析に入りますが，sample-to-conditionのファイルが必要です (一例です)．

sample condition
A1 Control
A2 Control
A3 Control
B1 Sample
B2 Sample
B3 Sample

 

そして解析に入ります．

## ライブラリの読み込み
library(sleuth)

## 結果ファイルのdirectory
base_dir <- "~/k_results"

# サンプル名の取得
sample_id <- dir(file.path(base_dir))
# サンプル名とそれぞれのkallisto結果directoryの連結
kal_dirs <- sapply(sample_id, function(id) file.path(base_dir, id))
# sample-to-condition table
s2c<-read.delim("sample_table.txt", header=T)

# kallistoの出力をcondition依存でsleuth object (so)として読み込む
so <- sleuth_prep(kal_dirs, s2c, ~ condition)

# <標準出力>
reading in kallisto results
............
normalizing est_counts
xxxxxx targets passed the filter
normalizing tpm
merging in metadata
normalizing bootstrap samples
# </標準出力>

# soをモデルにfitさせる
so <- sleuth_fit(so)

#<標準出力>
summarizing bootstraps
fitting measurement error models
shrinkage estimation
computing variance of betas
# </標準出力>

# fitされたモデルの確認．
models(so)

# <標準出力>
> models(so)
[ full ]
formula: ~condition
coefficients:
 (Intercept)
 conditionSample
no tests found.
# </標準出力>

# test
so <- sleuth_test(so, which_beta = 'conditionSample')
models(so)

# <標準出力>
> models(so)
[ full ]
formula: ~condition
coefficients:
 (Intercept)
 conditionSample
tests:
 conditionSample
# </標準出力>

# 結果の出力
results_table <- sleuth_results(so, 'conditionSample')

ここで色々とわかりにくいのがmodels(so)とsleuth_test，sleuth_resultsです．

models(so)はsleuth_fitでfitされたモデルの一覧をくれます．そのときの結果 (上記のsample_to_conditionを一例に)はこのようになります．

> models(so)
[ full ]
formula: ~condition
coefficients:
 (Intercept)
 conditionSample
no tests found

これはControlとSampleの二群間比較なのでこうなりますが，複数群間でやろうとするとこうなります． 

> models(so)
[ full ]
formula: ~condition
coefficients:
 (Intercept)
 conditionSample1
 conditionSample2
 conditionSample3
　　・
　　・
　　・
no tests found

# </標準出力>ですが，一つ注意しないといけないのは，これらはすべて"Control"との比較でのモデルです．現在のversionではこうですが，今後多群間の前組み合わせの比較が実装されるかも知れません．どのconditionXXXがどのサンプル間を比較しているのかを確認するにはdesign_matirx()を使えばわかります．  

 またsleuth_testでは，which_betaを指定する事でどの(design_matrixに基づいた)組み合わせで検定を行うかを指定します．*****一回のsleuth_testでは一つのwhich_betaしか指定できませんが，which_betaを変えながら複数回sleuth_testを行う事は可能です*****．そして気づくのに時間かかりましたが，このsleuth_testしなければsleuth_resultsで書き出せません．

 

 

ここまでいけば，q-value (BH補正されたp-value)がでるのでそこから機能解析やエンリッチメント解析などに進めます．また，sleuthにはデフォルトで様々なplot関数があるので，これを使ってもよし，shinyというさむしんぐ (sleuth_live)を使ってGUIで解析する事ができます．

 

 

この記事で一番言いたかったのは，

Kallistoは早い．

Sleuthは早い．

どっちもラップトップでも早い．

やべぇ．

 

です．

 

 

 

 

mclblastlineを使いたい

一つのFastaファイルの中で，配列類似性に基づいてクラスタリングしたい

→mcl (http://micans.org/mcl/)を使う

 

いろんなやりかたあるけど，mclblastlineを使ったら楽かもしれない？

## protein配列

formatdb -i input.fa -p T

 

blastall -p blastp -i input.fa -d input.fa -o input.fa.blast -m 8 -e 1e-10

mclblastline --mcl-I=1.0 --blast-m9  input.fa.blast

  # --mcl-te=20 でコア数指定

 

cytoscape用のsifファイルを作りたい

mcxdump -imx out.input.fa.blast.I10 -dump-sif A -tabr input.fa.blast.ta

 

という感じ

メモ：無名ハッシュを使ったときの配列のアクセス

これの後半のところで使おうとして、毎回迷ってるなぁって思ったから、メモ代わりに。

https://gist.github.com/t12968yy/8621f304c76c790caa46

my $tree;

for {

  push @{$tree}, $hoge;

}

 

foreach my $factor ( @{  $tree } ) {

 print $factor."¥n";

}

 

 

 

 

Perl :: Storable

perlで構造体を保存するための裏技っていうかあれ。

大抵はハッシュとかを保存したかったりでよく使います．

% cpanm Storable;

% emacs -nw foo.pl 

 

$!/usr/env/bin perl

use strict;

use warnings;

 

use Data::Dumper;

use Storable qw(store, retrieve);

 

##### 保存したいハッシュを %hashとします。

 

#構造体を出力

store \%hash, "<output file name>";

 

#別のコードで読み込むとき

my %hash = %{retrieve "<output file name>"};

 

 これで結構構造体とかを使い回したりしてます、でっかいデータを読み込むのを毎回やってたら時間食いますしね。