macでインフォマティクス

macでインフォマティクス

NGS関連のインフォマティクス情報についてまとめています。

RNA seq 其の4 おかしなデータを除いて再挑戦

30hのデータがおかしいので、それ以外の時間で分析を行う。tophatでのマッピングは其の2で行ったので、cuffdiffの定量だけやり直せばよい。これだけでかなり時間を短縮できる。

 

cuffdiff -o cuffdiff_except-_for_30h -p 20 Chlamydomonas_reinhardtii.v3.1.34.gtf -L 0h-1,0h-2,2h,4h,8h,12h,18h,24h,45h,60h 0h-1_out/accepted_hits.bam 0h-2_out/accepted_hits.bam 2h_out/accepted_hits.bam 4h_out/accepted_hits.bam 8h_out/accepted_hits.bam 12h_out/accepted_hits.bam 18h_out/accepted_hits.bam 24h_out/accepted_hits.bam 45h_out/accepted_hits.bam 60h_out/accepted_hits.bam

 

 Rを起動

library(cummeRbund) 

 cuff <- readCufflinks("cuffdiff_out")

genediff <- diffData(genes(cuff))

genediff2 <- genediff[(genediff$significant == 'yes'),]

write.table(genediff2, file="result3.csv", sep=",")

genediffdataID2 <- genediff2[,1]

genediffdataID2 <- t(genediffdataID2)

uesakaGenes2 <-getGenes(cuff, genediffdataID2)

k.means <-csCluster(uesakaGenes2, k=9)

k.means.plot <- csClusterPlot(k.means)

k.means.plot

 

k=9

f:id:kazumaxneo:20170329202356j:plain

 

 

 

k=6

f:id:kazumaxneo:20170329202520j:plain

 

 

k=4

f:id:kazumaxneo:20170329202602j:plain

 

k=9の 結果をよくよく見ると、30h以外にも24で発現が大きく落ちる遺伝子グループがある(1のウィンドウ)。もしかしたら30hの応答も生物的な応答かもしれない。 とりあえずk=9では十分なグループ分けができていない可能性もある。一度k=12で解析してみた。

k=12 f:id:kazumaxneo:20170330114409j:plain

 time 4h、12h、24h、45h、60hでもそれぞれ大きく落ち込む遺伝子グループがある。0hはそれ以前の変動がわからないが0h特異的に発現が少ない遺伝子群かもしれない。

 

k=20まで上げてみよう。

K=20

f:id:kazumaxneo:20170330114818j:plain

 まだまだパターンがある気がするし、人間の目で見てかぶったパターンも増えたように思える。nitrate responseのサンプルなので,0hから発現がグンと上がっている遺伝子が気になる。k=20のウィンドウ12の遺伝子をまとめた。

f:id:kazumaxneo:20170330133452j:plain

Nar1;2という硝酸同化系遺伝子が見事に応答している。非常に妥当な結果である。

 

 

 

 

ここで最初に立ち返ってt=30のデータはどうなのかを考える。例えばRNAの分解が進んだサンプルのrRNA/mRNA比率が代わり、特に低発現の遺伝子の発現パターンにばらつきが出るようになる。まあそれはそうだろう。分解スピードが同じと仮定し細胞内に10コピー存在するmRNAと100コピー存在するmRNAを考えると、先に0になるのは10コピーのmRNAの方だろう。このようなサンプルからRNAをとってシーケンスすると、分解してないサンプルと比べた時の比率が無限大になる。たとえFPKM 0.1と仮定して数値を加えたとしても、分解していないサンプルとの比率を正確に求めることはできない。また、分解しきっていなくても、低発現遺伝子の定量がばらつく要因はいくらでも考えられる。というのも、RNAの分解もシーケンス時の逆転写、PCR、シーケンスも(世の中のこと全て)確率的にしか進行しないため、初発の数値が低いまま進むとブレが大きくなるからだと解釈できる。例えると、サイコロを6回しか振らないと、6の目が3回出て、1は1回も出ない(オレすごい!)ようなことが起こるが、サイコロを6万回ふるとどの目も限りなく1/6に近い値がでる(オレ普通だった....)。それでも1/6に近づかないなら、サイコロが正6面体になっていないか重心がずれている可能性がある。

 

まずは色々な時間で散布図を書いてみる。

 

csScatter(genes(cuff),"X2h","X4h",smooth=T)

f:id:kazumaxneo:20170330121125j:plain

 

csScatter(genes(cuff),"X2h","X8h",smooth=T) 

f:id:kazumaxneo:20170330121152j:plain

 

 csScatter(genes(cuff),"X12h","X18h",smooth=T)

f:id:kazumaxneo:20170330121314j:plain

 

 csScatter(genes(cuff),"X18h","X45h",smooth=T) 

f:id:kazumaxneo:20170330121330j:plain

 

csScatter(genes(cuff),"X45h","X60h",smooth=T)

f:id:kazumaxneo:20170330121518j:plain

 

csScatter(genes(cuff),"X24h","X45h",smooth=T)

f:id:kazumaxneo:20170330121529j:plain

 

 

この中では45hが他と一番違うように感じる。さて、では30との比較だとどうなるか?

csScatter(genes(cuff),"X24h","X30h",smooth=T)

f:id:kazumaxneo:20170330121840j:plain

 

csScatter(genes(cuff),"X18h","X30h",smooth=T)

f:id:kazumaxneo:20170330121910j:plain

 

csScatter(genes(cuff),"X2h","X30h",smooth=T)

f:id:kazumaxneo:20170330122053j:plain

思ったほどばらつきはでかくなっていない。まあ対数表記だからか。しかし、どの時間と比較しても左下の方に弧を描くように描画されているプロット出る。説明できないけどアーティファクトっぽいなあ。

 

 

 

 

k=9の結果をヒートマップ化する。

f:id:kazumaxneo:20170329202356j:plain

k.means

出てきた数値をMiにコピー。エクセルに読み込んでフィルターをかけ、欲しいウィンドウだけ取り出してMiで保存。Rに再読み込み。

 

ウィンドウ1、8、3の遺伝子でヒートマップを出してみる。

csv <- read.table("1_3_8.txt" , header=F, sep="\t")

csv2 <- t(csv)

csv2data <-getGenes(cuff, csv2)

csHeatmap(csv2data)

f:id:kazumaxneo:20170330124213j:plain

 

ウィンドウ6、7の遺伝子でヒートマップを出してみる。

 

csv <- read.table("1_3_8.txt" , header=F, sep="\t")

csv2 <- t(csv)

csv2data <-getGenes(cuff, csv2)

csHeatmap(csv2data)

f:id:kazumaxneo:20170330124352j:plain

 

 思ったほどまとまらない。k=20の19と20のウィンドウだとどうだろうか?

f:id:kazumaxneo:20170330124751j:plain

 数が少ないのもあってさすがにこれは綺麗にまとまっている。

 

 

 

 

 

RNA seq 其の3 ヒートマップ作成

前回k-means法でクラスタリングするところまでやった。

 

k=9の結果が以下の通り。

 

f:id:kazumaxneo:20170329161538j:plain

k=9の結果についてヒートマップを書いて見る。

全部を1つのヒートマップに書くとしんどいので、分けて書く。まずは1、2、3のデータパターンが似たものをヒートマップ化する。

 

k.meansと打つと分析結果が表示されるので、

$silinfo

$silinfo$widths

                  cluster neighbor     sil_width

CHLREDRAFT_138936       1        2  0.4642057561

CHLREDRAFT_13257        1        2  0.4273940431

CHLREDRAFT_206075       1        2  0.4219634073

CHLREDRAFT_195343       1        2  0.4197930334

CHLREDRAFT_167270       1        2  0.4180063266

CHLREDRAFT_178353       1        2  0.4158867327

CHLREDRAFT_120099       1        2  0.4109362252

と書いてあるところをテキストエディタにコピーし、エクセルに読み込む。区分けはカンマかタブ。cluster neighbor のカラムでフィルターをかけ、1,2,3のclusterだけ取り出してcsv形式で保存 (1-3.csv)。

#読み込み

csv <- read.table("1-3.csv" , header=F, sep="\t")

 

前回と同じ

csv2 <- t(csv)

csv2data <-getGenes(cuff, csv2)

csHeatmap(csv2data)

f:id:kazumaxneo:20170329180206j:plain

まだデータが多い。

 

1のクラスタだけで描いてみる。

f:id:kazumaxneo:20170329180405j:plain

ようやく見える解像度になってきた。しかしデータに問題がある可能性がある30hのデータの落ち込みで分類されたぽいクラスタ1なので、生物的な応答パターンを見えているか怪しい。30hのデータを除いてcufflinksの定量からやり直してみる。

 

 

 次回へ

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

RNA seq 其の2 マッピングから統計処理まで

基本的な流れは以下のリンクが参考になる。

github.com

 

今回はSRAに登録されているクラミドモナスのtime courseデータを使う。クラミドモナスはEnsemblにドラフトだがfastaもgtfも登録されているので非モデル生物より解析は楽になる。

クラミドモナスゲノムのダウンロードは

Chlamydomonas reinhardtii - Ensembl Genomes 35

のページの中ほどにある

Download DNA sequence (FASTA)

をクリックする。FTPで繋がり、任意のクロモソーム配列をダウンロードできる。

 

 

2、rRNAリードの除去

最初にrRNAにマッピングして、rRNAリードを除く。

 

rRNAのインデックス作成

bowtie2-build -f rRNA.fa rRNA

-f fastaファイルを指定

インデックス名をrRNA。この場合、rRNA~という6つのファイルが出力される。

 

 

rRNAへのマッピング

bowtie2 -x rRNA -1 R1.fastq -2 2.fastq -S rRNA_mapped.sam --un-conc rRNA_unmapped.fastq

-x 先ほどつけたインデックス名

-1 ペアードエンドデータのR1

-2 ペアードエンドデータのR2

-S マッピングされたデータ出力名(sam形式)

--un-conc マッピングされなかったfastqのペア。シングルなら--unを使う。上の場合はrRNA_unmapped 1.fastqとrRNA_unmapped 2.fastqが出力される。

 

ファイルサイズは2.03GBから2.02GBまでしか減らなかった。polyAから逆転写しているデータと思われる。ということで、このデータではrRNAの除去ステップは行わない。

 

 

 

3、クロモソームへのマッピング

faファイルとgtfファイルはEnsembl genomeの

http://plants.ensembl.org/info/website/ftp/index.html

Chlamydomonas reinhardtii からダウンロード。

 

インデックス作成

bowtie2-build -f Chlamydomonas_reinhardtii.v3.1.dna.toplevel.fa Chlamydomonas_reinhardtii.v3.1.dna.toplevel

 

ペアードエンドッピング

tophat -G Chlamydomonas_reinhardtii.v3.1.34.gtf -p 8 -o output --no-novel-juncs Chlamydomonas_reinhardtii.v3.1.dna.toplevel R1.fastq R2.fastq

-o 出力ディレクトリ名

-G gtfファイル。

-p スレッド数。多すぎるとエラーになるので8くらいに留める。

--no-novel-juncs アノテーションにない新規のジャンクションを破棄

最後にindex名、R1とR2のfastqを指定している。

 

 

 

 

4、FPKM算出

Cufflinksを使って計算する。

cufflinks -o output_FPKM -G genes.gtf -N --compatible-hits-norm mapping_0h/accepted_hits.bam 

-G 既存のアイソフォームのみの計算の指示

N --compatible-hits-norm 既存のアイソフォームにマッピングされるものだけでFPKMを計算

 -o FPKMの出力ディレクトリ名

最後にtophatで出力したbamファイルを指定する。複数結果がある時は順番に実行する。

 

出力されたらFPKMのカラムをcutで抜き出す

cut -f 10 mapping_0h/accepted_hits.bam > 0h_FPKM.txt

cut -f 10 mapping_2h/accepted_hits.bam > 2h_FPKM.txt

 cut -f 10 mapping_0h/accepted_hits.bam > 4h_FPKM.txt

 

pasteコマンドでカラムを結合。

paste 0h_FPKM.txt 2h_FPKM.txt 4h_FPKM.txt > 0-4h_FPKM.txt

 

0-4h_FPKM.txtをエクセルかRに読み込んで経時変化をグラフ化する。

 

 5、Cuffdiffによる発現量に差がある遺伝子の抽出

cuffdiff -o cuffdiff_out -p 20 Chlamydomonas_reinhardtii.v3.1.34.gtf -L 0h-1,0h-2,2h,4h,8h,12h,18h,24h,30h,45h,60h 0h-1_out/accepted_hits.bam 0h-2_out/accepted_hits.bam 2h_out/accepted_hits.bam 4h_out/accepted_hits.bam 8h_out/accepted_hits.bam 12h_out/accepted_hits.bam 18h_out/accepted_hits.bam 24h_out/accepted_hits.bam 30h_out/accepted_hits.bam 45h_out/accepted_hits.bam 60h_out/accepted_hits.bam

 

ファイルが多いとかなりの時間がかかる。

下のような大量のファイルが出力される。

f:id:kazumaxneo:20170328180325j:plain

 

 

 

 

おまけ、 Go termによる分類

今回使用したChlamydomonasのgtfファイルにはEntrezのgene_id情報がある。これがあればすぐにGo解析ができる。それはChlamydomonasの情報はEMBLの運営するEnsembl Genomesからダウンロードしており、すでにEnsembl のaccession IDが付いているからである。それに対して、De novoのトランスクリプトーム解析ならアノテーションがないため、アノテーション、ID変換等を駆使していかなかればGo termをつけてカテゴリ分類することはできない。

 

AgBase GORetrieverを貼ったプレーンテキストをアップロードして、Ensembl accession IDを選択。チェックは必要に応じて外すが、全選択してjobをなげても問題ない。しばらく待つとab-deliminatedファイルが

ダウンロードできるリンクが発生する。クリックしてダウンロードする。ただし、GOの分類はそもそもまだまだ主観性の強い分類であり、分類が間違っている可能性がある。どの階層を選ぶかによって結果も変わってくるため、確度の高い分析を行うのは厳しい。あくまで大雑把にどのような応答になっているのか掴むための手法と捉えた方がよさそうである。

 

 

7、統計解析

Rで行う。まずは準備。

CummuneRbandのインストー

参照したページ。こんなことができるようです。

CummeRbund - An R package for persistent storage, analysis, and visualization of RNA-Seq from cufflinks output

 

 

#ターミナルでcufflinks解析結果ディレクトリの上に移動して、(つまり上の解析時そのままの位置で)  Rを起動(ターミナルでRと打ち、エンター)

R

 

library(cummeRbund) #ライブラリ読み込み

cuff <- readCufflinks("cuffdiff_out") #cufflinks解析結果のディレクトリ全体を読み込んでくれる

#これでCuffSet オブジェクトが生成された。

 

cuff

#正常に読み込まれていれば、cuffオブジェクトを分析しcuflinksの解析結果の概要が表示される。サンプルは0-60hのタイムコース実験の11サンプル。

CuffSet instance with:

11 samples

14487 genes

14553 isoforms 0

TSS 0

CDS 0

promoters 0

splicing 0

relCDS

 

 

上のgenesの分布

csBoxplot(genes(cuff))

f:id:kazumaxneo:20170329112140j:plain

上の図作成は

box <- csBoxplot(genes(cuff))

print(box)

のようにやってもよい。この場合はまずcsBoxplot(genes(cuff))の結果をboxに代入し、それをpirntコマンドで画面に描画している。

#さらにpngで保存するなら

png("box.png")

plot(box)

dev.off()

これを同時にうつ。

 

 

分子生物学ではあまり使わないので、box plot(箱ひげ図)の見かたについて触れておく。(参考)

1、箱の中央付近のヨコ線 → データxの中央値。平均値ではない。また必ずしもboxの中央にあるわけでもない点に注意。

 

2、箱のヨコ線 → データxの第1四分位数(下側)と第3四分位数(上側)

分かりにくいので言い換える。箱の底辺の線の位置は、データの中央値と一番小さい値との間の中央値。つまり、中央値を最大値と考え、その最大値から下半分のデータについての中央値を求め、その位置を箱の低線の位置としている。この1/4番目の区切りの中央値を第1四分位数という。同様に箱の上辺の線の位置はデータの中央値と一番大きい値との間の中央値(=第3四分位数)。

 

3、箱の上下に向かって縦方向に伸びている線(ヒゲ) → データの最大値と最小値。

 ここで言う最大値と最小値はデータの真の最大値と最小値ではなくて外れ値を除外したもの。ヒゲの下方向の長さは

箱の底辺 - (箱の高さ)x 1.5

で出たの範囲で最も小さな値となる。

上の式を難しく言い換えると

第1四分位数-1.5 x (第3四分位数-第1四分位数)

 

ヒゲの上方向の長さは

第3四分位数+1.5*(第3四分位数-第1四分位数)

となる。その範囲に入らなかった値が外れ値として黒いプロットで表現されている。

 

 

箱ひげ図は品質管理で用いられたりしている。RNA seqの品質を管理するものとして捉えればよい。箱の位置が被らないようなデータは有意な差がある可能性が高い。30hのデータを見ると、中央値は他のデータと大差ないが箱の底辺は遥か下まで伸びている。言い換えると第1四分位数が他のサンプルと大きく異なっている。

 

isoformのbox plotを書く。

csBoxplot(isoforms(cuff))

f:id:kazumaxneo:20170329121448j:plain

アイソフォームでも30hのデータはおかしい。これが生物的な応答なのかサンプルの品質なのかは後で追及する。

 

 

 

csDensity(genes(cuff)) #ヒストグラム

f:id:kazumaxneo:20170329112234j:plain

 

 

2サンプル間の散布図

samples(object)

を使う。

samples(cuff)

   sample_index sample_name sample_name parameter value

1             1       X0h_1        <NA>      <NA>  <NA>

2             2       X0h_2        <NA>      <NA>  <NA>

3             3         X2h        <NA>      <NA>  <NA>

4             4         X4h        <NA>      <NA>  <NA>

5             5         X8h        <NA>      <NA>  <NA>

6             6        X12h        <NA>      <NA>  <NA>

7             7        X18h        <NA>      <NA>  <NA>

8             8        X24h        <NA>      <NA>  <NA>

9             9        X45h        <NA>      <NA>  <NA>

10           10        X60h        <NA>      <NA>  <NA>

 

なので、サンプルX2hとX4hを使って散布図を書くなら

csScatter(genes(cuff),"X2h","X4h",smooth=T)

f:id:kazumaxneo:20170330113904j:plain

 

 

 

 

csDendro(genes(cuff)) #樹形図

f:id:kazumaxneo:20170329112344j:plain

 

 

 

 csScatterMatrix(genes(cuff))

f:id:kazumaxneo:20170329125505j:plain

 

 

dispersionPlot(genes(cuff))

f:id:kazumaxneo:20170329125957j:plain

 

 

 csVolcanoMatrix(genes(cuff))

f:id:kazumaxneo:20170329113458j:plain

 

 

一部分を拡大

f:id:kazumaxneo:20170329113602j:plain

横軸はlog2(fold change)

縦軸はlog10(p value)

 

 

 

 

 #genediffというオブジェクトに差のある遺伝子データ(genes)を読み込む。

genediff <- diffData(genes(cuff))

 

#最初の部分だけを表示

head(genediff,3)

gene_id sample_1 sample_2 status value_1 value_2 log2_fold_change test_stat p_value q_value significant

1 AAB93440 X0h_1 X0h_2 NOTEST 1.34796 4.71785 1.80735 0 1 1 1 no

2 AAB93441 X0h_1 X0h_2 NOTEST 0.00000 0.00000 0.00000 0 2 1 1 no

3 AAB93442 X0h_1 X0h_2 NOTEST 0.00000 0.00000 0.00000 0 3 1 1 no

 

ここでRの基本的なコマンドを学習しておく。

 

 

 

 

 

 

#一番前にgene_idがあることに注目して差のある遺伝子のIDのカラムだけ抽出、genediffdataIDに代入。

genediffdataID <- genediff[,1] 

 

#中身をチェック

head(genediffdataID,10)

[1] "AAB93440" "AAB93441" "AAB93442"

[4] "AAB93443" "AAB93444" "AAB93445"

[7] "AAB93446" "AAB93447" "CHLREDRAFT_100005"

[10] "CHLREDRAFT_100008"

#headはカラムが少ないとこのように1つの行に3行分表示したりする。分かりにくい。

ファイルに保存してexcelで開いて見ると

write.table(genediffdataID, file="result.csv", sep=",")

f:id:kazumaxneo:20170329154145j:plain

こんな感じ。

 

 #データ解析には、縦のカラムではだめで行ベクトルに縦横変換する必要がある。縦横変換は

genediffdataID2 <- t(genediffdataID)

 

#ファイルを保存して

write.table(genediffdataID2, file="result.csv", sep=",")

 

#excelで開くと縦横変換されている。

f:id:kazumaxneo:20170329153630j:plain

 #ちなみにこれをターミナルからカッコよく開くなら

open result.csv #csvのデフォルトアプリで開く

 

#uesakaGenesオブジェクトにcuffオブジェクトからこのgenediffdataID2を持つ遺伝子の情報だけ代入

uesakaGenes <-getGenes(cuff, genediffdataID2)

 

#uesakaGenesをもとにヒートマップを書かせる

csHeatmap(uesakaGenes)

f:id:kazumaxneo:20170329155326j:plain

 

#数が多すぎてよくわからない。発現量が多い遺伝子だけに絞り込む。

#そのため差のある遺伝子のIDのカラムだけ抽出のところからやり直す。

genediff <- diffData(genes(cuff))

 

#統計的に優位なものだけ取ってくる。(significant列がyes)

genediff2 <- genediff[(genediff$significant == 'yes'),]

 

genediff2

#うてば何行残っているかわかる。2527行だけ残っていた。

#ファイルに保存してexcelで開くと

f:id:kazumaxneo:20170329160516j:plain

#右端の列がyesになっていることが確認できる。

あとは同じような流れ。

 

genediffdataID2 <- genediff2[,1]

genediffdataID2 <- t(genediffdataID2)

uesakaGenes2 <-getGenes(cuff, genediffdataID2)

 

 

 

 #Rでは k-means によるクラスタリングを行う関数が標準で用意されている。

#はじめてk-means法を聞く人が概要を掴むならiOSアプリの

アルゴリズム図鑑

アルゴリズム図鑑

  • Moriteru Ishida
  • 教育
  • 無料

 の中のk-meansの説明で素晴らしく分かりやすく説明されている。

 

 

k.means <-csCluster(uesakaGenes2, k=4)

k.means.plot <- csClusterPlot(k.means)

k.means.plot

f:id:kazumaxneo:20170329161325j:plain

 

 

#k=6にすると

k.means <-csCluster(uesakaGenes2, k=6) 

k.means.plot <- csClusterPlot(k.means) 

k.means.plot

f:id:kazumaxneo:20170329161726j:plain

 

 

 

 

 

#k=8にすると

k.means <-csCluster(uesakaGenes2, k=8) 

k.means.plot <- csClusterPlot(k.means) 

k.means.plot

f:id:kazumaxneo:20170329161455j:plain

 

 

#k=9にすると

k.means <-csCluster(uesakaGenes2, k=9) 

k.means.plot <- csClusterPlot(k.means) 

k.means.plot

f:id:kazumaxneo:20170329161538j:plain

 

 

#k=12にすると

k.means <-csCluster(uesakaGenes2, k=12) 

k.means.plot <- csClusterPlot(k.means) 

k.means.plot 

f:id:kazumaxneo:20170329161601j:plain

これはさすがにやりすぎ。直感的には9が過不足なく綺麗に分けられている気がするがどうだろう。

 

30hで沈むデータがかなりあるのがわかる。box plotで30hだけおかしかったのはこれを反映しているのだろう。ぶっちゃけて言うと、おそらく30hのデータを除かないと正しい解析はできないだろう。

 

次回に続く。

 

 

 

 

 

 

 

 

 

RNA seq 其の1 クオリティチェック

1、クオリティチェック

変なリードが混じっていると、本来シーケンスされた領域以外の遺伝子にマッピングされたりして解析結果を歪める恐れがある。そのため、シーケンスが終わりBCLからfastqデータが出力されたら、一度クオリティをチェックする。

クオリティがイマイチな場合、リード低クオリティ領域のトリミングや、リードの除去を行う必要がある。

 

トリミングはphread quality scoreに基づいて行う。クオリティ分布をグラフ化したり、指定した閾値にしたがって5'側や3'側をトリムするソフトはたくさん存在するし、有償の次世代データ解析ソフトにもほぼ間違いなく実装されている。

オススメはfastx-toolkit

使い方だが、

https://bi.biopapyrus.net/rnaseq/qc/fastx-toolkit.html

の東大の方?の運営されているHPにわかりやすくまとめられている。

(いつも勉強させてもらっています)

インストールはhomebrewを使って

brew install homebrew/science/fastx_toolkit

使い方は先ほどのページ通りです。

 

Rならqrqcが色々見れて便利。

ターミナルを起動し、シーケンスfastqがあるディレクトリに移動してRを起動(osx 10.11以上でRを標準搭載)

$R #R起動のコマンド

>install.packages("qrqc", dependencies = TRUE)

インストール参考ページ

 

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

fq <- readSeqFile("DRR016695_1_fill.fastq") #データの読み込み

makeReport(fq)

これで実行ディレクトリにフォルダができているはず。その中に下のデータが全て保存されている。

 

勉強を兼ねて個別に実行してみる。

クオリティ分布を見るには

qualPlot(fq)

 

f:id:kazumaxneo:20170327122303j:plain

 

塩基の出現頻度は

basePlot(fq)

f:id:kazumaxneo:20170327122350j:plain

 

リード長の分布

seqlenPlot(fq)

f:id:kazumaxneo:20170327122411j:plain

 

 

 

K-mer 頻度を見るには(アセンブリの参考になる)

kmerKLPlot(fq)

f:id:kazumaxneo:20170327122732j:plain

 

 

 

複数ウィンドウを残しながら進めるなら 

dev.new() #これで新しい白紙ウィンドウができる。

qualPlot(fq)

dev.new()

basePlot(fq)

dev.new()

seqlenPlot(fq)

dev.new()

kmerKLPlot(fq)

結果は以下のようになる。

f:id:kazumaxneo:20170327123215j:plain

 

 

あとから追加編集したくなった場合、例えば2番ウィンドウをアクティブにするには dev.set(2)

 

 

アクティブなデバイスを閉じる

dev.off()

すべてのデバイスを閉じる

graphics.off()  

 

実際のトリミングはastx-toolkitで行えばよい。処理スピードも早いし、柔軟なトリミングが可能。問題はどこまで貴重なデータを捨てるか?ということになる。真核生物のRNA seqでは最低2000万リードは必要と言われたりする。これは細胞内で100万分の1分子数だけ発現している(=低発現のmRNA)を20リードシーケンスでき、ある程度の統計的処理に対応できる数ということをどこかで聞いた。

時々低発現のmRNAはstochasticな発現なのかよく分からんし高発現のものだけ捉えられたらいいというような発現を聞くことがあるが、ome解析なのだから低発現のRNAも含めてシーケンスされているのが理想である。場合によっては隠れた登場人物を見過ごして結論を出してしまう恐れもあるのだから。

 

ということで、トリミングを欲張ればデータがなくなるし、トリムを全くしないと、エラーリッチなリード異なるmRNAにマッピングされる恐れがある。phread quality scoreが20の場合、経験則で1%の確率でエラーであることがわかっている(画像データから塩基を決定する時のアナログ->デジタル変換時の間違いなど?)。

1%というと低いようだが、100bpシーケンスすれば、1リードにつき1塩基は間違っている可能性が高いことになる。Q=10なら エラー率10%となる。100bpシーケンスすれば、10塩基くらい間違っている。リピートに対して十分な精度があると言えるだろうか?

ということで、Q=20くらいでトリムする人が多いのではないだろうか?

ただし個人が開発したプログラムによってはペアが同一数あるのを前提に作られているソフトもある。トリムすると、ペアの個数が異なりその後の処理が行えない場合もある。その時はhomemadeなスクリプトを書いてペアードエンドのリード数を少ない方に合わせる必要がある。私は海外の次世代フォーラムみつけたperlのスクリプトを使ってます。実行するには

perl To_Remove_Unpaired_Reads.pl trimmed_R1.fastq trimmed_R2.fastq

To_Remove_Unpaired_Reads.plは自分がコピペして保存した適当な名前。

引数としてトリムして数が異なってしまったR1とR2のfastqを指定する。

 

 

ただし、ペアの個数がより現実的な問題になるには、何らかの理由によってR2だけクオリティが悪くてR2のリードがR1より随分減ったときなどである。ペアリードとしてマッピングするなら、数が少ない方にリード数を揃える必要があり、ますますリードが減ってしまう。ぶっちゃけそんなデータを解析に使うな、となるが、高額の研究費を使って初めてRNA seqやる人ならそこで諦める人はいないだろう。

ということでトリミングなしでマッピングしてまた問題になったりする。書いているときりがないが、クオリティ以前の問題で、非モデル生物からRNAを無理やりとってRIN値の低いサンプルを使ってRNA seqすることもあるかもしれない。rRNA/mRNA値が変わると特に低発現のmRNAの分布のばらつきが大きくなることがわかっている。

RNA seqをしっかりやるには簡単ではない。

 

脱線してしまったが、クオリティ解析が終わればマッピング解析を行う。それは次回に

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

CIRCOS トレーニング2

バクテリアアセンブルデータをCIRCOSで描画してみる。

 

最初に、ゲノムサイズが小さいのでメモリのサイズを小さくする必要がある。main.confファイルの中のchromosomes_unitsを10万に変更。

chromosomes_units = 100000

 

 

main.confファイルの内容は以下の通り。

> $cat main.conf

# 1.2 IDEOGRAM LABELS, TICKS, AND MODULARIZING CONFIGURATION

karyotype = /Users/user/Documents/circos_training/plectonema/chr_orf.txt

chromosomes_units = 100000

 

<<include ideogram.conf>>

<<include ticks.conf>>

 

<image>

<<include etc/image.conf>>                

</image>

 

<<include etc/colors_fonts_patterns.conf>> 

<<include etc/housekeeping.conf>> 

 

 

上で読み込んでいるchr_orf.txtは以下のように記述して保存している。

chr - chr1 1 0 6176365 chr1

 

ゲノムが1つで環状なので、1行だけになる。6176365がゲノムサイズ。

 

 

 

 

 

 

そのほか5つのファイルを準備する。

IDEOGRAM.CONF

IDEOGRAM.POSITION.CONF

IDEOGRAM.LABEL.CONF

TICKS.CONF

BANDS.CONF

 

 

それぞれのファイルの中身は以下の通り。

> $cat IDEOGRAM.CONF

<ideogram>

#boxの情報

<spacing>

break   = 0.5r

default = 0u #これを0にすることで環状になる

</spacing>

 

<<include ideogram.position.conf>>

<<include ideogram.label.conf>>

<<include bands.conf>>

 

radius*       = 0.825r

 

</ideogram>

 

 

> $cat IDEOGRAM.POSITION.CONF

radius           = 0.775r

thickness        = 50p

fill             = yes

fill_color       = white

stroke_thickness = 2

stroke_color     = black

 

 

> $cat IDEOGRAM.LABEL.CONF

show_label       = yes

label_font       = default

label_radius     = dims(image,radius)-30p

label_size       = 24

label_parallel   = yes

label_case       = lower

label_format     = eval(sprintf("chr%s",var(label)))

 

 

> $cat TICKS.CONF

show_ticks          = yes

show_tick_labels    = yes

 

<ticks>

radius           = dims(ideogram,radius_outer)

orientation      = out

label_multiplier = 1e-4

color            = black

size             = 20p

thickness        = 3p

label_offset     = 5p

<tick>

spacing        = .1u

show_label     = no

</tick>

 

<tick>

spacing        = .5u

show_label     = yes 

label_size     = 18p

format         = %d

</tick>

 

<tick>

spacing        = 5u

show_label     = yes

label_size     = 40p

format         = %d

</tick>

 

</ticks>

 

 

> $cat BANDS.CONF

show_bands            = yes

fill_bands            = yes

band_stroke_thickness = 2

band_stroke_color     = white

band_transparency     = 0

 

 

 

 

 

 

node.txtを解析を行うディレクトリ(/Users/user/Documents/circos_training/)に保存。ファイルの中身は下記の通り。

cat node.txt

 

chr - node1 node1 0 449142 chr1

chr - node2 node2 0 398826 chr2

chr - node3 node3 0 350896 chr3

chr - node4 node4 0 246509 chr4

chr - node5 node5 0 235822 chr5

chr - node6 node6 0 214352 chr6

chr - node7 node7 0 181306 chr7

chr - node8 node8 0 141627 chr8

chr - node9 node9 0 70584 chr9

chr - node10 node10 0 28980 chr10

 

 

実行する。

circos -conf main.conf

 正常に終了するとpngファイルが出力される。

 

f:id:kazumaxneo:20170327102027j:plain

うまくできた。

 

 

 

 

 

今度はorf情報を追加してみる。orfの描画は1例としてtile機能で表現できると記載されている。

http://circos.ca/documentation/tutorials/2d_tracks/tiles/lesson

 

<plots> </plot>で囲んだ領域にorf情報を記述する。main.confは以下のようになる。

 

 

 

# 1.2 IDEOGRAM LABELS, TICKS, AND MODULARIZING CONFIGURATION

karyotype = /Users/user/Documents/circos_training/plectonema/chr_orf.txt

chromosomes_units = 100000

 

<<include ideogram.conf>>

<<include ticks.conf>>

 

<image>

<<include etc/image.conf>>                

</image>

 

<plot>

type            = tile

layers_overflow = hide

 

 <plot>

 file        = sense_orf.txt

 r1          = 0.95r

 r0          = 0.95r

 orientation = out

 layers      = 1 #レイヤー1

 margin      = 0.01u

 thickness   = 50

 padding     = 1

 stroke_thickness = .5

 stroke_color     = black

 </plots>

 

 <plot>

 file        = antisense_orf.txt

 r1          = 0.90r

 r0          = 0.90r

 orientation = out

 layers      = 1 #レイヤー2

 margin      = 0.01u

 thickness   = 50

 padding     = 1

 stroke_thickness = .5

 stroke_color     = red

 </plots>

 

</plots>

 

 

<<include etc/colors_fonts_patterns.conf>> 

<<include etc/housekeeping.conf>> 

 

 

orf情報は以下のようなファイルを準備した。

> $cat sense_orf.txt

chr1 130 1470 + color=black

chr1 1550 2029 + color=black

chr1 4431 5162 + color=black

...

chr1 5782854 5782925 + color=black

chr1 5936365 5936436 + color=black

chr1 5936365 5936436 + color=black

 

 

> $cat antisense_orf.txt

chr1 2080 2418 - color=red

chr1 2408 2698 - color=red

chr1 2761 4080 - color=red

...

chr1 5772159 5772229 - color=red

chr1 6077090 6077162 - color=red

chr1 6077090 6077162 - color=red

 

 

 

今回はゲノム配列を一旦RASTサーバーでorf解析し、生じたgtfファイルから加工して作った。加工はスクリプトを使わずともexcelですぐ修正できる。

gtfがない場合、UCSCからダウンロードする。

 

なければNCBIかかDDBJからgene bankファイルを落とし変換する。もちろん自前のperlのスクプトで処理してもいい。


 

 

IDEOGRAMのboxだが、orfのboxがあればなしでもいいかもしれない。だが視覚的には輪郭線んに相当するものがあったほうがコントラストもでて見やすくなる。そこで今回は細いIDEOGRAMをorfの外に描画する。IDEOGRAM.LABEL.CONFを以下のように修正した。

 

> $cat IDEOGRAM.LABEL.CONF

radius           = 0.775r

thickness        = 20p

fill             = yes

fill_color       = white

stroke_thickness = 1

stroke_color     = black

 

 

 

実行する。

circos -conf main.conf

 正常に終了するとpngファイルが出力される。

f:id:kazumaxneo:20170327102105j:plain

次はRNA seqデータを表示してみる。全部は表示しないが、このようなものを作成できる。

 

 

 

 

 

RNA seqデータ。

f:id:kazumaxneo:20170327102156j:plain

ゲノムの繰り返し配列(LASTを使った)

f:id:kazumaxneo:20170327102205j:plain

今の所使う予定はないが、ヒストグラムや散布図も描画することができる。

自分が必要な機能だけ図を描きながら学んでいくのが手っ取り早いと思われる。

 

 

 

CIRCOS トレーニング

circosは多様な機能を備えたゲノムのビジュアル化ツールである。類似ソフトにDNA plotter、genome diagram、snap geneなどがあるが、機能の豊富さではcircosが抜きん出てる。circos独自の機能として、遺伝子の相関や相同性を線で表現してlinkageを表す機能が充実していることが挙げられる。

例えばキーワード"circos"で画像検索して確認してください。→をクリック(search)

 

難点として、機能が豊富な分使いこなすにはかなり勉強が必要である。

 

インストール(mac OS Sierra 10.12.2で検証した)

brew install circos

 

ランするには、circos.confファイルを作成してターミナルで下のように打つだけである。

circos -conf circos.conf

 この.confファイル中に図を書くのに必要なパラメータを記載する。circosはその情報に従ってランしているわけである。

confファイルには非常に多様な情報を記載することができるが、プログラムに慣れていないといきなり理解するのは難しい。circos公式サイトのチュートリアルを真似て、1つずつ理解していくのが一番である。一日頑張れば、欲しい図は大体かけると思う。

下にチュートリアルのリンクを載せておく。http://circos.ca/documentation/tutorials/ 

 

トレーニング1

circos.confの中身は、トレーニングページに公開されている。例えば練習1の

http://circos.ca/documentation/tutorials/quick_start/ticks_and_labels/configuration のページの # 1.1 MINIMUM CIRCOS CONFIGURATION ... ... <<include etc/housekeeping.conf>>

までを全てプレーンテキストエディタ(mi とか)に 貼り付ける。circos.confとして保存して、

circos -conf circos.conf

と実行する。 正常に終わるとpngsvgファイルができる。DNA plotterならマウスオーバーで 情報を表示したり編集したりできるが、他ソフトと異なりcircosの結果は画像で出力される。後から編集することはできない。練習ページ1の出力結果が下になる。

f:id:kazumaxneo:20170327101206j:plain

 

 

上の図のパラメータは以下の通り

default = 0.005r  #box間のスペース

radius = 0.90r  #サークルのサイズ

thickness = 20p  # boxの厚み

fill = yes #box内部のカラー

stroke_color = dgrey

stroke_thickness = 2p #box輪郭線の厚み

 

 

 

上のパラメータを

default = 0.05r

fill             = no

stroke_thickness = 10p

に変更すると、下のようになる。

f:id:kazumaxneo:20170327101348j:plain

 

 

circos.confの中の

# Chromosome name, size and color definition

の次の行に

karyotype = /Users/user/karyotype.human.txt

という行がある。このkaryotype.human.txtはあくまで

例だが、中身は以下のような構造をしている。

 

chr - hs1 1 0 249250621 chr1

chr - hs2 2 0 243199373 chr2

chr - hs3 3 0 198022430 chr3

.

.中略

.

band hs1 p36.33 p36.33 0 2300000 gneg

band hs1 p36.32 p36.32 2300000 5400000 gpos25

band hs1 p36.31 p36.31 5400000 7200000 gneg

band hs1 p36.23 p36.23 7200000 9200000 gpos25

.

.中略

.

band hsY q11.223 q11.223 22100000 26200000 gpos50

band hsY q11.23 q11.23 26200000 28800000 gneg

band hsY q12 q12 28800000 59373566 gvar

 

886行目まである(karyotype.human.txtはcircosのexampleファイルの中に入っている)。

circos.confのkaryotype = /Users/user/karyotype.human.txt

の行は、karyotype.human.txtのパスを指定している。karyotype.human.txtを保存したパスに修正しておく。パスの修正はその都度行う必要がある。はじめのうちはミスしやすいところなので注意する。

karyotype.human.txtはhttps://github.com/vigsterkr/circos/tree/master/data

からダウンロードできる。文字をコピペしてkaryotype.human.txtとして保存。

 

 

トレーニング2

各クロモソームにアノテーション情報や目盛りをつける。

 

http://circos.ca/documentation/tutorials/quick_start/ticks_and_labels/configuration

 

の3つのコードを

CIRCOS.CONF

IDEOGRAM.CONF

TICKS.CONF

それぞれ上の名称で保存する。実行はcircos -conf circos.conf

下の画像が出力されるはず。

f:id:kazumaxneo:20170327101611j:plain

CIRCOS.CONF の中の

 

<<include ideogram.conf>>

<<include ticks.conf>>

 

 の2行でIDEOGRAM.CONFとTICKS.CONFを読み込んでいる。

boxの情報はideogram.confファイルの中の

radius           = 0.90r

thickness        = 20p

fill             = yes

stroke_color     = dgrey

stroke_thickness = 2p

に書かれている。

 

アノテーション情報は

show_label       = yes

# see etc/fonts.conf for list of font names

label_font       = default 

label_radius     = dims(image,radius) - 60p

label_size       = 30

label_parallel   = yes

に書かれている。

 

 

 

目盛りはticks.confファイルの中の

show_ticks          = yes

show_tick_labels    = yes

 

<ticks>

radius           = 1r

color            = black

thickness        = 2p

 

でメモリの厚みなどを規定。さらに、下の方の

<tick> #1種類目

spacing        = 5u

size           = 10p

</tick>

 

<tick>#2種類目

spacing        = 25u

size           = 15p

show_label     = yes

label_size     = 20p

label_offset   = 10p

format         = %d

</tick>

 

 

の2つの<tick> </tick>の中で、2種類の目盛りを指定している。

f:id:kazumaxneo:20170327101703j:plain

 

上の図のように、5ポイントごとの目盛りと、25ポイントごとの目盛りが描画される。1つ目の目盛りは文字なしなので、

show_label     = yes

label_size     = 20p

label_offset   = 10p

format         = %d

の行は書かれていない。

 

 

肝心のクロモソーム情報はkaryotype.human.txtから読み込んでいる。karyotype.human.txtファイル冒頭の

 

chr - hs1 1 0 249250621 chr1

chr - hs2 2 0 243199373 chr2

chr - hs3 3 0 198022430 chr3

.

.中略

.

chr - hs22 22 0 51304566 chr22

chr - hsX x 0 155270560 chrx

chr - hsY y 0 59373566 chry

 

を元にboxの名前とサイズが決められている。

 

 

 

 

 

 

トレーニング3

各クロモソームのサイズや向き、場所を変える。

 

http://circos.ca/documentation/tutorials/quick_start/selection_and_scale/configuration

から

の3つのコード

CIRCOS.CONF

IDEOGRAM.CONF

TICKS.CONF

を保存する。実行はcircos -conf circos.conf

下の画像が出力されるはず。

f:id:kazumaxneo:20170327101801p:plain

CIRCOS.CONFファイルの中の

chromosomes_scale   = hs1=0.5r,/hs[234]/=0.5rn

でクロモソームのサイズを規定している。

 

chromosomes_reverse = /hs[234]/

でクロモソーム2、3、4の向きを逆転させている。

 

chromosomes_color   = hs1=red,hs2=orange,hs3=green,hs4=blue

でクロモソームの色を指定している。

 

chromosomes_radius  = hs4:0.9r

でクロモソーム4の位置を他より内側に指定している。

 

 

 

 

 

トレーニング4

各クロモソーム間のリンクを表示する練習。

 

http://circos.ca/documentation/tutorials/quick_start/links_and_rules/configuration

から

の3つのコード

CIRCOS.CONF

IDEOGRAM.CONF

TICKS.CONF

を保存する。

https://github.com/vigsterkr/circos/blob/master/data/5/segdup.txt

の文字をsegdup.txtとして保存する。

CIRCOS.CONFの58行目を

file          = segdup.txt

にパスを変更する。

 

実行はcircos -conf circos.conf

下の画像が出力されるはず。

f:id:kazumaxneo:20170327101826j:plain

 

 

 

今日はこんなところで。

トレーニング2に続く。

 

 

 

 

 

 

 

 

 

 

large indelの検出ツール

バクテリアゲノムを使い、構造変化検出ツールのパフォーマンスを比較したペーパーが出ている。

実際に導入して、パフォーマンスを比較してみる。

 

 

はじめにbrew tap でsiience系のコマンドをオフィシャルコマンドとしてインストールできるようにしておく。

brew tap homebrew/science

 

ほとんどのスクリプトはインプットとしてソート済みのbamファイルを使う。このbamの作成については、small indelとSNV検出編を確認してください。SNVとsmall indel検出前にbamのリアライメントとBQSR補正を行い、できる限り正確なアライメントを行ったbamを出力する流れを書いています。

 

R1.fqとR2.fqがfastqファイルで、ref.faがリファレンスとする。

リファレンスのindexも必要とするものが多い。samtools faidxでfasta.faiを作っておく。

samtools faidx $f

 これらの結果を1つのフォルダにまとめておき、各スクリプト実行時にそのフォルダを参照するようにする。 

 

ショートリードから構造変化を予測する手法は大きく4つに分けることができる。

1、Alignment based 

Dindel

Stampy

SAMtools

GATK

Varscan

 

2、Split-read mapping

Pindel

SV-M

 

3、Paired-end mapping

Hydra

PEMer

Breakdancer

 

4、haplotype-based

GATK-HaplotypeCaller

Platypus

 

 

 

Breakdancer Chen et al. (2009)

discordant paired read approaches

 

 breakdancerのようなread pairの手法はインサートサイズのずれリードの向き(正常なら--> <--)からindelを予測する。つまり、ペアリードのインサートサイズが一般的な正規分布に従うと仮定し、ペア間の距離が分布の外れ値になるような部位にindelがあり、リードの向きがおかしくなる(--> -->や<-- <--)部位にインバーションがあるだろうと考える方法論である。Breakdancerはこのような異常なペアリード; discordant paired read(DPR)情報に基づき構造変化を予測する。ただし、アライメントのエラーもありうるので、エラーを考量して異常なペアが1箇所に2組以上集積している部位をコールしている。

 ただしこの手法には大きな問題点がある。まずペア間の距離が大きく変わっていることが前提なのでsmall indelは検出できない。また、インサートサイズの分布によって感度が変化する。さらに構造変化後の配列に関する情報も得られないこともデメリットとなる。

 そのため現在ではread pairの手法で構造変化に関する手がかりを掴み、それをsplit-mappingやlocal de novo assemblyを組み合わせて検出する手法がいくつか報告されている。

 

2017年現在Breakdancer-maxがダウンロードできるが、依存するモジュールが多くてインストールは面倒だった(cpanでたくさんのperlモジュールをインストールする必要があった)。 

 

解析手順は、まずbamファイルをサンプリングしてインサートサイズなどを分析したファイル(.cfgファイル)を作り、それからコールを行う。

#cfgファイルを作成する(orted.bamはbwa memで作成した)。

perl bam2cfg.pl R1R2_sorted.bam > output.cfg

cfgファイルの中身は以下のようになっていた

 >cat output.cfg
readgroup:X platform:ILLUMINA map:R1R2_sorted.bam readlen:296.80 lib:Y num:9990 lower:48.86 upper:29146.26 mean:444.35 std:3920.08 SWnormality:minus infinity exe:samtools view

サンプリング自体は最初の数百行だけ行うらしく、すぐに終わる。ただしこの時qualityがpoorなデータだと解析できない (ランした時にpoor quality dataと出る)。 どうしてもランしたければ、オプション"-v"をdefaultの1から2-10まであげると一応ランはできる。

 

#構造変化のコール

 breakdancer-max output.cfg > output.txt

 人のchr20で試すと解析は10分で終わり、500箇所ほどコールされた。

 

 

Pindel Ye et al. (2009)

split-read approaches.

 

ショートリードのペアードエンドデータからdeletion(1bp-10kb)、medium insertion (1-20bp)を予測するツール。リードを分割してアライメントし、そのアライメントのパターンから構造変化を予測している。このようなアライメントをスプリットアライメントと呼び、BWA MEMやRNA-seqで使うtophatなども使っている方式となる。

手法の詳細は論文に書いてあるが、split-read mappingの代表的な論文なのでもう少し手順を説明する。

1、ペアリードの片側だけマップされるペアを探索する。この片側のリードのアライメントが肝なので、ユニークでかつミスマッチがないパーフェクトマッチのアライメントが行われる。

2、相方リードがマップされないのは、この相方リードが構造変化部位を横断するようにマップされているからと考える。この相方リードのアライメント部位を探すため、アライメントできた方のリード(アンカーと呼んでいる)からインサートサイズ平均の2倍の範囲内でマッピングされなかったリードの方がマッチする部位を探す。ただし相方リードはそのままではアライメントできないので、リードを5'側と3'側、真ん中の3ピースに分けてそれぞれでアライメントを実行する(split alignmnet)。

3、5'側と3'側のアライメント部位にずれがあれば、即ちそのずれのサイズがdeletionのサイズになる(下図a)。5'側と3'側のアライメント部位にズレはないが、真ん中のピースが浮いた状態でマッピングされていれば、その真ん中のピースがinsert配列そのもののはずである。

4、2つ以上のペアで同じ情報が検出された部位をコールする。

f:id:kazumaxneo:20170327095850j:plain

片方のリードを正確にアライメントし、その相方の位置も限定することで偽陽性を抑え且つ解析をスピードアップしている。原理上リードの長さによって検出可能なindelの長さの上限が決まる。初期の30-40bpのリードでは20bpのinsertionが上限ということらしい。

 

以下の図は解析の流れ。

f:id:kazumaxneo:20170509104834j:plain

 

http://gmt.genome.wustl.edu/packages/pindel/install.html にインストールの仕方が書いてある。macOSX環境ではインストールできなかったのでcentOSにインストールした。

git clone git://github.com/genome/pindel.git
cd pindel
./INSTALL /path/to/samtools_FOLDER/

 

デモランは以下の通りに行う。

cd demo
./pindel -f simulated_reference.fa -i simulated_config.txt -o output

ランに必要なものはリファンレス配列(.fasta)とそのindexファイル(.fasta.fai)。またそのリファレンスにmapして作ったbamファイルとそのindexファイル(.bam.bai)。このあたりは他のツールと変わらない。

pindelは情報を記載したコンフィグファイルを作って、それを元にランする。サンプルファイルたと、コンフィグファイルに相当するのはsimulated_config.txtファイルとなる。ファイル内には以下のことが書かれている。

simulated_sample_1.bam  250     SAMPLE1

simulated_sample_2.bam  250     SAMPLE2

simulated_sample_3.bam  250     SAMPLE3

データが1つなのでこれを以下のように変更した。

R1R2_sorted.bam  600     SAMPLE1

この600はインサートサイズ。ランするには

./pindel -f input.fasta -T 12 -i config.txt -o output

-T: thread

-f: ref.fa

-i: config_file

-o: output

終わると構造変化の種類ごとにアライメントファイルが出力される。結果を見るにはpindel2vcfで変換する。pindel2vcfはビルドしてでできたpindel/に入っている。

pindel2vcf -p pindel_output.txt_D -r ref.fa -R ref -d 20170509 -v pindel_output_D.vcf -e 5
pindel2vcf -p pindel_output.txt_INT_final -r ref.fa -R ref -d 20170509 -v pindel_INT_final.vcf -e 5
pindel2vcf -p pindel_output.txt_INV -r ref.fa -R ref -d 20170509 -v pindel_output_INV.vcf -e 5
pindel2vcf -p pindel_output.txt_SI -r ref.fa -R ref -d 20170509 -v pindel_output_SI.vcf -e 5
pindel2vcf -p pindel_output.txt_TD -r ref.fa -R ref -d 20170509 -v pindel_output_TD.vcf -e 5
pindel2vcf -p pindel_output.txt_BP -r ref.fa -R ref -d 20170509 -v pindel_output_BP.vcf -e 5

-r/--reference The name of the file containing the reference genome: required parameter -R/--reference_name The name and version of the reference genome: required parameter -d/--reference_date The date of the version of the reference genome used: required parameter -p/--pindel_output The name of the pindel output file containing the SVs  

-e/--min_supporting_reads  The minimum number of supporting reads required for an event to be reported (default 1)

出力のvcfフォーマットはvcf higher versionに準拠しているらしい。

-eは必須パラメータではないが、デフォルトの1では高感度すぎて明確な変異なしのシミュレートデータ(50カバレッジ)でも1000近くの偽陽性コールが検出される。indelソフトのパフォーマンス比較の論文で-e 5になっていたのでここではそれに習った。ただし、目的やリードのカバレッジによってチューニングが必要と思われる。

 

 

 

GASVpro Sindi et al. (2012)

discordant paired read approaches

ダウンロードサイト。

 

インストールは、解凍したディレクトリの中で以下のコマンドを打つだけでよい。

./install

 

入力はbamファイルとなる。ランは2段階で行う。

java -Xms512m -Xmx2048m -jar BAMToGASV.jar R1R2_sorted.bam -MAPPING_QUALITY 30 -CUTOFF_LMINLMAX SD=3

-MAPPING_QUALITY: 

-CUTOFF_LMINLMAX: 

 

 

LUMPY layer et al. (2014)

split-read + discordant paired read + CNV based approaches. 

 解析の流れはこのURLに丁寧に書いてある。

 

macOSX環境でのインストールも可能。homebrewでインストールできる。

brew install homebrew/science/lumpy-sv

ソフトクリップのマッピングファイルから構造変化を予測する手法。マッピングソフトは1Mbpまで対応したbwa memなどのマッパーを使う。デプスのないデータなどでもsensitivityが高いらしい。また複数サンプルで差分をとって精度を上げるやり方も搭載している。

f:id:kazumaxneo:20170327095430j:plain

テストランはtestディレクトリにあるtest.shを実行する。

./test.sh

問題なければNA12878、honeybee、riceのindel予測結果がvcf形式で出力される。

 

 解析時は下のような流れでコールする。

bwa index -a bwtsw ref.fa
bwa mem -R "@RG ID:X LB:Y SM:Z PL:ILLUMINA" ref.fa R1.fastq R2.fastq |samtools view -S -b - > sample.bam
samtools view -@ 20 -b -F 1294 sample.bam > sample.discordants.unsorted.bam
samtools view -@ 20 -h sample.bam | ./scripts/extractSplitReads_BwaMem -i stdin | samtools view -Sb - > sample.splitters.unsorted.bam
samtools sort -@ 20 sample.discordants.unsorted.bam > sample.discordants.bam
samtools sort -@ 20 sample.splitters.unsorted.bam > sample.splitters.bam
lumpyexpress -B sample.bam -S sample.splitters.bam -D sample.discordants.bam -o sample.vcf

解析できるまで手間取った。最後の肝心のlumpyexpressがモジュールを読み込めずエラーになる。最終的に強引にrootで進めることでなんとかランできた。

samtoolsのソートでエラーになる人は、>(リダイレクト)とリダイレクトの後ろ側の.bamを決してランする。これはsamtoolsのバージョンによっては>と.bamが入らないバージョンがあるため。

 

 

 

Hydra-sv

discordant paired read approaches. 

 

マニュアル: Google Code Archive - Long-term storage for Google Code Project Hosting.

 

下のような流れでコールする。(https://code.google.com/archive/p/hydra-sv/wikis/TypicalWorkflow.wiki

bwa index -p $f -a is $f
bwa aln -t 12 $f $p > R1.sai
bwa aln -t 12 $f $q > R2.sai
bwa sampe $f R1.sai R2.sai $p $q |samtools view -Sb - > sample.tier1.queryorder.bam
samtools view -uF 2 sample.tier1.queryorder.bam |bamToFastq -bam stdin -fq1 R1_disc1.fq -fq2 R2_disc1.fq
novoindex ssuis.nix $f #Novoalignのindex
novoalign -d ssuis.nix -f R1_disc1.fq R2_disc1.fq -i 500 50 -r Random -o SAM |samtools view -Sb - > sample.tier2.queryorder.bam
samtools view -uF 2 sample.tier2.queryorder.bam |bamToFastq -bam stdin -fq1 R1_disc2.fq -fq2 R2_disc2.fq
novoalign -d ssuis.nix -f R1_disc2.fq R2_disc2.fq -i 500 50 -r Ex 1100 -t 300 -o SAM |samtools view -Sb - > sample.tier3.queryorder.bam 
bamToBed -i sample.tier3.queryorder.bam -tag NM |pairDiscordants.py -i stdin -m hydra -z 800 > sample.disc.bedpe
dedupDiscordants.py -i sample.disc.bedpe -s 3 > sample.disc.deduped.bedpe
hydra -in sample.disc.deduped.bedpe -out sample.breaks -mld 500 -mno 1500

 

ただし結果はひどくしょぼい。進め方かパラメータを間違ってる可能性がある。

 

 

 

 

 

DELLY2 Rausch et al. (2012)

split-read +discordant paired read approaches. 

 

Delly2 can discover,deletions, tandem duplications, inversions and translocations at single-nucleotide resolution in short-read massively parallel sequencing data.

deletionなどを予測する。コピーナンバーが変わるリピートのdeletionや、depthが少ないデータにも強いらしい。 コントロールのデータ(somatic sampleなど)と解析データ(tumor sampleなど)両方を使った差分解析を特徴とするが、単体のデータでも解析は可能。insertionはターゲットではないので注意する。こちらからダウンロードできる。

 

macOSX環境ではインストールできなかったのでcentOSサーバーにインストールした。

基本的にはgitからダウンロードしてmakeするだけ。

git clone --recursive https://github.com/dellytools/delly.git 
cd delly/
make all

 

ランは

./delly call -t DEL -x human.hg19.excl.tsv -o delly.bcf -g ref.fa R1R2_sorted.bam

-xはなくても問題ないが、コールを除外する領域を設定できる(例えばヒトのセントロメアなど)。

コントロールと比較するなら

./delly call -t DEL -x human.hg19.excl.tsv -o delly.bcf -g ref.fa R1R2_sorted.bam Control_R1R2_sorted.bam

Contorl~bamは変異がないと期待される体細胞シーケンスデータ。tumor~bamが変異の可能性がある腫瘍組織のシーケンスデータ。

 

ランが終わるとbinaryのbcfファイルが出力される。結果を見るにはbedtoolsで変換する必要がある。

./bcftools/bcftools view delly.bcf > delly.vcf

この時フィルターコマンドは実行することも可能。

 

5つのコールタイプがあるので逐一実行する。

./delly call -t DEL -o del.bcf -g ref.fa R1R2_sorted.bam &
./delly call -t DUP -o dup.bcf -g ref.fa R1R2_sorted.bam &
./delly call -t INV -o inv.bcf -g ref.fa R1R2_sorted.bam &
./delly call -t TRA -o tra.bcf -g ref.fa R1R2_sorted.bam &
./delly call -t INS -o ins.bcf -g ref.fa R1R2_sorted.bam &
./bcftools/bcftools view del.bcf > deletions.vcf 
./bcftools/bcftools view dup.bcf > duplications.vcf
./bcftools/bcftools view inv.bcf > inversions.vcf
./bcftools/bcftools view tra.bcf > translocations.vcf
./bcftools/bcftools view ins.bcf > insertions.vcf

 

 

Platypus Rimmer et al. (2014)

local assembly based approaches. 

nature geneticsに載った手法。原理の詳細はこちら http://nextgenseek.com/2014/07/platypus-a-new-variant-caller-that-integrates-mapping-assembly-and-haplotype-based-approaches/

 

候補部位をローカルリアライメントしてプロタイプを調べ、さらにlocal de novo assemblyも行いSNV、indel、complex polumorphismsを検出する手法。パフォーマンス比較の論文でも精度が高いと評価されている。

ふつうショートリードのアライメントは各々の塩基が一致するかに焦点を当てて行われる。この方式は計算リソースの面から考えると大量のデータを処理するには向いているが、indel周辺はミスマッチが多発してしまうことになる。そこでindelを考量して計算リソースを潤沢に使った丁寧なギャップアライメントを行う機能がいくつかの手法に実装されている(GATKなど)。ただしギャップアライメントでうまく補正できるのはせいぜい数塩基で、long indelが起きていた場合修正しきれない。Platypusはそのような補正が出来ない領域だけを対象にde novoでアセンブルを行い、long indelの検出を試みる手法である。対象領域を限定することで、解析スピードを大きく縮小していることも特徴の1つで、例えば数メガのバクテリアの解析だと早ければ1分以内でランは終わる。

 

macOSでも動く。インストールにはgccとcythonが必要。

brew install gcc #ない人だけ 
pip install cython #python2環境なので、pip3を使ってはダメ

samtoolsとhtslibも必要。

brew install homebrew/science/htslib #ない人だけ
brew install homebrew/science/samtools #ない人だけ
pip install cython #python2環境なのでpip3を使ってはダメ

GitHub - andyrimmer/Platypus: Platypus Variant Caller

からPlatypusをzipをダウンロードして解凍。中に入って

make

これで完了。

ランに必要なリファレンスのfaiファイルを作成。

samtools faidx ref.fa #faiファイル作成

 

最後にラン。ソートしたbamファイルを使う。

python bin/Platypus.py callVariants --refFile=ref.fa --bamFiles=R1R2_sorted.bam

--refFile: index付きのfastaファイル

--bamFiles: bamの指定。bamは複数入力可能。その時はbamをコンマで区切る。

パスを通しておけばPlatypus.pyだけで動く。

 

デフォルトではサポートするリードが2で検出するため、ディープにシーケンスしたデータでは感度を下げたい場合がある。サポートするリードを4に変えるなら

python bin/Platypus.py callVariants --refFile=ref.fa --bamFiles=R1R2_sorted.bam --minReads 4

他にもCPU数のオプションや(--nCPU)、リージョンを限定するオプション(--regions)、またVCFファイルを読み込ませて変異をアノテートするらしきオプション(--source)もある。--sourceは試していない。詳細はHP参照。

 

 解析速度は非常に早く、3Mゲノムだと10秒以内にrunは終わった。ヒトゲノム30xでも20分で終わるらしい。

 追記;人のデータでも試してみた。hg19のchr20をリファレンスにしてテストしたところ(gz圧縮シーケンスデータは2GBx2あったが)、10分以内に解析は終了した。 

 

fermikit Li (2015)

assembly based approaches. 

 

一度アセンブルを行なって、それからアライメントをとるプログラム。そう聞くと重そうだが、ヒトゲノムデータでもピークメモリー使用量は90GBくらいで、アセンブルも24時間くらいで終わると書いてある。またその先のアライメントは1時間半くらいで終わるとのこと。

macOSでも動く。インストールは

git clone --recursive https://github.com/lh3/fermikit.git
cd fermikit
make

ビルドしてできたfermi.kitフォルダの中に実行に必要なものが入っている。アライメントにbwaを使っている。

実行にはリファレンスのbwa indexが必要なので最初に作っておく。 

bwa index ref.fa

アセンブル条件ファイルの作成

fermi2.pl unitig -s 3m  -t 12 -l 301 -p prefix R1.fastq > prefix.mak

-s: 推定ゲノムサイズ(3mは3 Mbp)

-l: リード長 (Miseq 301bp)

-p 出力されるprefix名

ペアーリードには対応していないみたいで、ペアリード指定の方法が見つからない。

 アセンブル

make -f prefix.mak

structural variationsのコール

fermi.kit/run-calling -t 12 ref.fa prefix.mag.gz | sh

prefix.mag.gzとprefix.flt.vcf.gzというファイルができる。READMEには

`prefix.flt.vcf.gz` for filtered SNPs and short INDELs and `prefix.sv.vcf.gz` for long deletions,

に書いてある通り、prefix.flt.vcf.gzの中にSV予測結果、prefix.sv.vcf.gzの中にlong del予測結果が保存されている。playtus同様、解析時間は非常に短い。playtusやfermikitはPEM法などで予測が難しいsmall indelもコールしてくれるのが1つの利点と言える。

 

 

 

FreeBayes Garrison and Marth (2012)

gapped alignment approaches.  

 

macOSでも動く。小さなindel専用の検出ツール。ベイズ的アプローチでlndelを検出する。 倍数体ゲノムの変異検出にも対応している。ショートリードデータからhaplotypeの変異を検出するために開発された。

brewでインストール可能。

brew install FreeBayes
freebayes -f ref.fa -b R1R2_sorted.bam > out.vcf 

感度は-Cで調整できる。デフォルトの"-C 2"では感度が高すぎて、例えば50リードシーケンスリードがある部位で2リードだけッミスマッチがあると検出される。腫瘍の変異などわずかな変異を検出するには都合良いが、この条件では偽陽性も大量に出てくる。クローナルなサンプルなどで偽陽性を減らすには-Cを5ぐらいまであげるといい。 ただしはプロタイプの解析など感度が必要な解析では注意して進める。

freebayes -C 5 -f ref.fa -b R1R2_sorted.bam > out.vcf 

追記;クラミドモナス(120Mb)で解析速度を調べてみた。 データは150bp x 2 (6.8GB x 2) ) 

 だが、CPU=20にすると2分で終わり、他のツールと比べて圧倒的に高速なのがわかった。

 

GATK

split-read approaches. 

多分一番有名な有名なツール。 macOSでも動く。https://software.broadinstitute.org/gatk/documentation/topic?name=methods に詳細が記載されている。

 

UnifiedGenotyperのラン

gatk -T UnifiedGenotyper -R ref.fa -I  R1R2_sorted.bam -o UnifiedGenotyper.txt 

 HaplotypeCallerのラン

gatk -T HaplotypeCaller -R ref.fa -I R1R2_sorted.bam -o HaplotypeCaller.txt

 

 

SvABA Wala et al. (2017)

local assembly approaches. 

 

gitのリンク。

macでもインストール可能とあったが、ビルド中にエラーが出たのでcentOSにインストールした。

git clone --recursive https://github.com/walaj/svaba 
cd svaba
./configure
make
make install

ビルドが終わると、src/svaba/に本体がある。

 ランは

svaba run -p 12 -t R1R2_sorted.bam -G re.fa

-t: bamファイル。

-G: リファレンス。

~indel.vcfが出力ファイル。 

-p: thread数

 

bedファイルを指定することで解析領域を限定することも可能。マニュアルには300bpまでのinsertionをコールできるとあるが、検証したところ確かに300bpまでの長さのInsertionを検出できていた。

 

 

 

VarScan Koboldt et al. (2009)

gapped alignment approaches. 

 

BLAT, Newbler, cross_match, Bowtie and Novoalignなどの多くのアライナーをサポートしている。

 

macOSでも動く。brewでインストール可能。

brew install VarScan

Varscanはpileupデータを使うため、初めにソート済みbamファイルをpileupする必要がある。

samtools mpileup -f re.fa R1R2_sorted.bam > mpileup

SNVをコール

VarScan pileup2snp mpileup > output.txt

indelをコール

VarScan pileup2indel mpileup > output.txt

オプションとして以下のようなものがある。

--min-coverage Minimum read depth at a position to make a call [8]

--min-reads2 Minimum supporting reads at a position to call variants [2]

--min-avg-qual Minimum base quality at a position to count a read [15]

--min-var-freq Minimum variant allele frequency threshold [0.01]

--min-freq-for-hom Minimum frequency to call homozygote [0.75]

--p-value Default p-value threshold for calling variants [99e-02]

--variants Report only variant (SNP/indel) positions [0]

 

pileupしたファイルは大きいので、パイプを使い1ライナースクリプトで進めるやり方がHPに紹介されている。例えばindelをコールするなら

samtools mpileup -f ref.fa R1R2_sorted.bam | VarScan pileup2indel > output.txt 

 

 

Breseq Deatherage et al. (2014)

split-read approaches.

 マニュアルリンク

 

macOSでも動く。brewでインストール可能。

brew install Breseq

indelコールは

breseq -r ref.fa -j 12 R1.fq R2.fq

-r; リファレンスを指定。fast.gbk (genebank形式)、GFF3に対応。

-j: スレッド数

-0: アウトプット

遺伝子のアノテーションがついたgff3かgenebankファイルがあるなら、検出結果を遺伝子の部位情報付きで出してくれるので非常に便利。genbankファイルは、NCBIftpサイトか、もしくは該当ページの情報を配列付きで表示させて全テキストコピーしてもいい。gff3ファイルはUCSCなどからダウンロードできる。genebankリファンレスでのランは

breseq -r ref.gbk -j 12 R1.fq R2.fq

感度調整のオプションはないが、-pをつけると多型検出ができる。helpのメッセージを載せておく。

  -p,--polymorphism-prediction     The sample is not clonal. Predict polymorphic (mixed) mutations. Setting this flag changes from CONSENSUS MODE (the default) to POLYMORPHISM MODE

 

 

解析が終わるとhtmlのまとめファイルができるのが特徴。内部でRなどの描画ライブラリを動かしているぽい。/outputの中のsummary.htmlを開くと、

f:id:kazumaxneo:20170512131212j:plain

上のようにアミノ酸置換まで含めて報告されている。 

変異部位をクリックするとアライメント結果も表示される。すごい。

f:id:kazumaxneo:20170511134951j:plain

 

coverage

f:id:kazumaxneo:20170426164310j:plain

deletion予測時は下のようにリードのデプスを見ることもできる。f:id:kazumaxneo:20170513170011j:plain

 

 

 indelの抜け漏れはこの中で一番少なかった。インフォマティクスのツールに不慣れなユーザーにも使いやすいと思われる。ただし、動作が複雑なためか他のソフトの安定性かわからないが、時折解析中にエラーを起こす。そうゆう時は大体パラレルにBreseqを動作させていた時だったりする。複数走らせないほうがいいみたい。

追記;バクテリア解析用に設計されており、真核生物ゲノムだとかなり重たい。(100Mbゲノムだと20コア使用で数日かかった)

 

 

Scalpel Narzisi et al. (2014)

assembly based approaches. 

このスクリプトは人をターゲットにしている。腫瘍サンプル特異的な変異の差分検出もできる。オプションにもdad(父)とかmam(母)とかある。人のサンプルを使うことがあれば試してみたい。一応ランの流れを書いておく。

 

解析にはbedフォーマットのポジション情報が必要(リンク)。 

bcftoolsを使ってbamファイルから変換することができる(SEQanswers)。bamファイルはあらかじめソートしてindexをつけておく。

bamToBed -i R1R2_sorted.bam > R1R2_sorted.bed 

Runにはデータベースファイルが必要。

scalpel-discovery --single --bam R1R2_sorted.bam --bed R1R2_sorted.bed --ref ref.fa

上のコマンドを実行するとoutdir/フォルダに色々なファイルができる。その中にvariants.db.dirというファイルができている。これがデータベースファイル。データベースファイルができないエラーもあるらしい

indelの解析

scalpel-export --single --db outdir/variants.db.dir --bed R1R2_sorted.bed --ref ref.fa 

人データでテストしたところ、エラーなく最後まで走ったが、コールがゼロになる。

現在のところ未解決

 

 

Breakseek Zhao et al. (2015)

breakpoint + discordant paired read approaches. 

ダウンロードリンク

 

pythonで動くため、ソースコードのビルドは必要ないがpython の2.7が必要。

またランにはsamファイルが必要。推奨される作り方がマニュアルに書いてあるのでそれに従う。bwa backtrackを使う。

bwa index -a bwtsw ref.fa
bwa aln ref.fa read1.fq > read1.sai 
bwa aln ref.fa read2.fq > read2.sai
bwa sampe ref.fa read1.sai read2.sai read1.fq read2.fq > aln_pe.sam

テスト用のsamとfaが用意されている。下記URLからダウンロードする。

https://sourceforge.net/projects/breakseek/files/toy_testset/

以下はそのファイルで進行している。 

 

bwa alnでsamを作ったら、sam_prep.pyでクロモソームごとに分割する。

python sam_prep.py -f test.sam

終わるとfasta内のクロモソームごとにディレクトリが作られ、中にsamができる。このsamファイルを使ってランする。

python breakseek.py -f testS.sam -m 600 -s 100 -c 50 -q 250 -n 12 -r test.fa

-m: インサートサイズ

-s: インサートサイズのSD。正規分布ならインサートサイズの10-20%くらいらしい。

-c: リードデプス

-q: リード長

以下のメッセージ

gathering statistics on the dataset...

が出てから、解析途中は何もlogが出てこない。初めはランに失敗したと思っていたが、強制終了せず気長に待ったところ、1時間近くかかった。3Mのバクテリアゲノムを使っただけでこれは長い...。プログラムの最適化の度合いが低いのかもしれない。

 

 samファイルはbwa alnでなくてもOK。例えばbwaの別のアルゴリズムbwa memを使うなら

bwa index -a is test.fa
bwa mem -R "@RG ID:X LB:Y SM:Z PL:ILLUMINA" -t 12 test.fa R1.fq R2.fq > R1R2.sam 

is: bwaのindexアルゴリズム。"bwasw"か"is"。

-t: スレッド数

R1.fq: fastqファイル

R2.fq: fastqファイル

 

 

 

Scanindel Yang et al. (2015)

gapped alignment + split reads and + assembly based approaches. 

performance比較の論文で高評価されているツール。

 

動作に必要なものがgitのダウンロードページに記載されている。

Softwares and Python packages:

 

上記の大半はbrewとpipコマンドですぐ導入できたが、Inchwormはmacで動作しなかった。最終的にubuntu 14.04にインストールした。まずbrewで導入できるコマンドでないものがあれば入れておく。

brew listでインストール済みリストのチェック。

brew list

 下のコマンドのうちlistででないものがあれば導入。

brew install bwa 
brew install samtools
brew install bedtools
brew install blat
brew install freebayes

pysam、PyVCF、Biopython、SciPy、NumPyもpythonのpipコマンドで導入できる。PyVCFはpythonのSNV検出スクリプト

インストールされているpythonライブラリの確認 

pip freeze 

ないものがあれば導入。

pip install pysam 
pip install PyVCF
pip install Biopython
pip install --user numpy scipy matplotlib ipython jupyter pandas sympy nose

Inchworm assemblerはTrinitiyの中のアセンブラで知ってる人も多いと思う。TrinityのソースコードHPからダウンロードしてインストールする。 上記ファイルをzipでダウンロードして解凍。中にあるInchwormディレクトリに移動。

unzip trinityrnaseq-Trinity-v2.3.2.zip
cd trinityrnaseq-Trinity-v2.3.2/Inchworm/

マニュアルに従って Inchwormをビルド。

./configure --prefix=`pwd` 
make
make install

 bin/内に実行ファイル"inchworm"ができているはず。これにパスを通すかusr/bin内にコピー。

これでようやくランができるようになる。

 

ランは

python ScanIndel.py -i sample.txt -p config.txt [options]

ScanIndelは内部でblatを動かし、自動で上のコマンドを走らせポートを開けてリファレンスとのマッチを行う。そのため、blat検索時のポート番号の指定などのオプションを備えている。

 

 

 

 

-------参考-------

scanindelはリードとリファンレスの相同部位検索にblatを用いている。blatはblastに類似した相同部位検索プログラムで、indexファイルの作成がいらないなどの特徴をもつ。詳細はこちら

 blatスタンドアローンとwebサーバー版(UCSC)がある。最初にbinaryファイルを作るところはどちらも同じ。

faToTwoBit ref.fa ref.2bit

そのあとサーバー版ではgfServerコマンドを使いポートを開ける。

gfServer -trans -mask start localhost 17775 ref.2bit & 
gfServer start localhost 17774 ref.2bit &

 localhost自分で指定する名前。

 17774は自分で指定するポート番号。

 

 

 ----エラーが出たとき-----

ただしたくさんのスクリプトを統合しているためか、エラーが出てきたときはlogファイル(gfserver.temp.log)を開きどこで止まったか確認して進める。今回はbwa memのプロセスでエラーが出た。

これはsamtools sortコマンドがバージョンにより

samtools sort input.bam input_sorted

でsorted.bamができるbwaのバージョンと

samtools sort input.bam > input_sorted.bam

にしないとsortred.bamができないバージョンがあるかららしい。この件はgithubのScanindelのところでもコメントで表示されている。logにbwa memのソートでエラーが出ていた場合、pythonスクリプトを以下のように修正。

まずScanIndel.pyをテキストエディタで開く。

94行目の

"bwa.bam"+output+".bwa.sorted"

"bwa.bam >"+output+".bwa.sorted.bam"

に変更。さらに241行目の

'temp.bam' + output + '.temp.sorted'

'temp.bam >' + output + '.temp.sorted.bam'

に変更すれば動く。

 

 

ScanIndelのランに失敗すると、gfServerが停止せずポートが開きっぱなしになる。 このままだと、次回ランでエラーになる。ポート番号を--gfServer_portオプションで指定するか、ポートを閉じておく。ポートを閉じるには、まずポート番号をlsofコマンドで確認し、

lsof -i

表示されたプロセス番号が2587だとしたら

kill 2587

 これで停止できる。その他、途中で止まるとログファイル(gfserver.temp.log)も残っている。これも次回ランでエラーを起こすぽいので、エラー内容を確認したら消しておく。 その他、シミュレーションででsmall indelを検出しようとすると、エラーログなしに解析が止まることもたくさんあった。その場合、large indelを少し加えると、ランは最後まで行くようになった。そのとき、small indelも正確に検出されていたので、何かしら小さなバグだと思われる。ランが途中で止まる人は大きなindelを擬似的に加えてやり直してみると良い(あまり懸命じゃないが)。

 

その他、シミューレーションデータで検証時、短いindelを検出しようとすると途中で止まることがあった。原因は不明だがlong indelを混ぜ込むとランは正常に終了したので、対処療法的にバグ回避はできたが原因は不明。

 

 

 

PRISM Jiang et al. (2012)

discordant paired read + split-read approaches. 

 linuxのみサポートされている。

公式ページ

 インストールは

tar xvzf PRISM.x86_64.tar.gz
cd PRISM
file bin/smapper
export PRISM_PATH=$PWD

PRISMは大きく分けて4つのステップからなる。それらのtoolはtoolkitディレクトリ中に保存されている。ただし、実際のランはrun_PRISM.shを走らせるだけでよい。

 

テストランは

./run_PRISM.sh -m500 -e30 -i /home/user/PRISM_1_1_6/sample/chr100.sam -r /home/user/PRISM_1_1_6/sample/reference/chr100.fa

samファイルとfaファイルはフルパスで指定する。

ランが成功すると、カレントディレクトリにPRISM_outputができ、中にdiscordantなリードペアのクラスターがsam形式で出力される。これが構造変化の種類ごとに分けて出力されるのでかなりたくさんのsamができる(詳細はREADMEに書いてある)。肝心のindelはsplit_all.sam_ns_rmmul_cigar_sorted_svに保存されている。

 

パスを通してからのランは下のように行う。 

run_PRISM.sh -m600 -e100 -i /home/user/input.sam -r /home/user/ref.fa -I PRISM_working_directory -O output

 -Iと-Oでそれぞれ作業ディレクトリと出力ディレクトリを明示する。インサートサイズとSDを大きく間違うと途中で止まることがあるみたいので注意する。

 

 

 

Manta Chen et al. (2016)

discordant paired read + split-read approaches. 

ダウンロードサイト

tumorとnormalを比較して、高い感受性でガン変異を検出する機能を持つ。

インストールガイドにはmacでも動作可能なはずと書かれているが、macはインストールが面倒だったので、binary版がダウンロードできるcent OSで進めた。

MantaはSun grid engineのような分散コンピュータ環境にも対応しているらしいが、ゲノムサイズがそれほど大きなければローカル環境で十分高速に動く。

はじめにbin/の中にあるconfigManta.pyにパスを通しておく。

ランは、はじめにコンフィグレーションファイルを作り、それを元に解析する2段階方式dである。

1、コンフィグディレクトリの作成。シングルのdiploidゲノムなら

configManta.py --bam NA12878_S1.bam --referenceFasta ref.fa --runDir oputput

--bam: ペアエンドのbamファイルの指定。

--referenceFasta: リファレンスfastaの指定。

--runDir: 出力ディレクト

 

2、ラン シングルdiploidデータ

python runWorkflow.py -m local -j 12

 -m: --mode=MODE  select run mode (local|sge)。ローカルなら"local"を指定。

-j: ジョブ数。SGE環境なら128まで対応。デフォルト8。

 

同じサンプルデータが複数あるときのコンフィグファイル作成 (fastqをマージしても同じ)

configManta.py --bam 1.bam --bam 2.bam --bam 3.bam --referenceFasta ref.fa --runDir oputput

 

Mantaはデフォルトではwhole genome DNA seq解析モードになっているが、オプションを指定することでRNA seq, exome、target seqの変異解析を行うことも可能になっている。

 

腫瘍と体細胞の比較時のコンフィグファイル作成

configManta.py --normalBam HCC1187BL.cram --tumorBam HCC1187C.cram --referenceFasta ref.fa --runDir oputput

 

腫瘍データのみの時のコンフィグファイル作成

configManta.py --tumorBam HCC1187C.cram --referenceFasta ref.fa --runDir oputput

 

RNA-seqデータ(まだ不完全らしい)

configManta.py --rna --Bam 1.bam --referenceFasta ref.fa --runDir oputput

RNAモードでは1データのみ現在は対応。

 

 

 

 

 

metaSV Mohiyuddin et al. (2015)

hybrid approaches

これまでに報告されてきたSV検出法を複数組み合わせ、特に検出が難しいInsertionの検出感度を上げるように工夫された手法。Pindel、Brreakdancer、gartk、lumpy、manta、cnvnator、break、Breakseq、などのSV検出結果を入力し、さらにmetaSVは独自にsoft-clipされたリード情報を使いindelを予測する。予測後に、spadesを使いアセンブルして精度を高めるプロセスも導入されている。

 

テストデータが入っており、ランは簡単にテストすることができるが、実際にmetaSVを使う場合は準備に手間がかかる。

 

現在検討中。

 

 

 ざっくり総評

 シミュレーションデータでホモの構造変化の検出テストを行うと、Fermikit、ScanIndelが最も検出率が高く、コール内容も正確だった。次点でBreseq、pindelの検出率が高かった。特にPindelは特にレアな構造変化の検出に強かった。これはリードの個数を指定して感度をカバレッジに応じて変化できる点が効いていると思われる。

small indelの検出だと、Platypusとpindelの検出感度が優れていた。ただし2つとも、デフォルトではサポートするリードが2つでバリアントを検出する設定のため、ディープにシーケンスしたデータだとエラーも拾いやすい。1st screeningではオプションをつけて感度を下げて検出する工夫も必要と思われる。

 

 

 

 

 

 

 

 

 

 

 

 

 

 

引用

-------------------------------------------------------------------------------------------------------------------------

BreakDancer 

Identification of Genomic Structural Variation from Paired-End Read Mapping

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4138716/

 

LUMPY

A probabilistic framework for structural variant discovery

https://genomebiology.biomedcentral.com/articles/10.1186/gb-2014-15-6-r84

 

Pindel

A pattern growth approach to detect break points of large deletions and medium sized insertions from paired-end short reads

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2781750/

 

GASVpro

An integrative probabilistic model for identification of structural variation in sequencing data

https://genomebiology.biomedcentral.com/articles/10.1186/gb-2012-13-3-r22

 

DELLY2

Structural variant discovery by integrated paired-end and split-read analysis

https://academic.oup.com/bioinformatics/article/28/18/i333/245403/DELLY-structural-variant-discovery-by-integrated

 

PRISM

PRISM: Pair-read informed split-read mapping for base-pair level detection of insertion, deletion and structural variants

https://academic.oup.com/bioinformatics/article/28/20/2576/202193/PRISM-Pair-read-informed-split-read-mapping-for

 

Platypus

Platypus: A New Variant Caller that Integrates Mapping, Assembly and Haplotype-based approaches

http://www.nature.com/ng/journal/v46/n8/full/ng.3036.html

 

FermiKit

FermiKit: assembly-based variant calling for Illumina resequencing data

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4757955/

 

Freebayes

Haplotype-based variant detection from short-read sequencing

https://arxiv.org/pdf/1207.3907.pdf

 

SvABA

SvABA: Genome-wide detection of structural variants and indels by local assembly

http://biorxiv.org/content/biorxiv/early/2017/02/01/105080.full.pdf

 

VarScan

VarScan: variant detection in massively parallel sequencing of individual and pooled samples

VarScan: variant detection in massively parallel sequencing of individual and pooled samples | Bioinformatics | Oxford Academic

 

Scalpel

Accurate de novo and transmitted indel detection in exome-capture data using microassembly

http://www.nature.com/nmeth/journal/v11/n10/full/nmeth.3069.html

 

Breakseek

BreakSeek: a breakpoint-based algorithm for full spectral range INDEL detection

https://academic.oup.com/nar/article/43/14/6701/2902746/BreakSeek-a-breakpoint-based-algorithm-for-full

 

INDELseek

INDELseek: detection of complex insertions and deletions from next-generation sequencing data

https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-016-3449-9

 

Breseq

Identifying structural variation in haploid microbial genomes from short-read resequencing data using breseq

https://bmcgenomics.biomedcentral.com/articles/10.1186/1471-2164-15-1039

 

ScanIndel

ScanIndel: a hybrid framework for indel detection via gapped alignment, split reads and de novo assembly

https://genomemedicine.biomedcentral.com/articles/10.1186/s13073-015-0251-2

 

PRISM

PRISM: Pair-read informed split-read mapping for base-pair level detection of insertion, deletion and structural variants

https://genomemedicine.biomedcentral.com/articles/10.1186/s13073-015-0251-2

 

Manta

Manta: rapid detection of structural variants and indels for germline and cancer sequencing applications

https://academic.oup.com/bioinformatics/article-lookup/doi/10.1093/bioinformatics/btv710

 

MetaSV

An accurate and integrative structural-variant caller for next generation sequencing

Bioinformatics first published online April 10, 2015 doi:10.1093/bioinformatics/btv204