読者です 読者をやめる 読者になる 読者になる

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も登録されているので非モデル生物より解析は楽になる。

 

 

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解析ができる。というのも、de novoのトランスクリプトーム解析ならアノテーションがないため、アノテーション、ID変換等を駆使していかなかれば、ならないがchlamydomonasの情報はEMBLの運営するEnsembl Genomesからダウンロードしているもので、すでにEnsembl のaccession IDが付いているため。

 

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 トレーニング

類似ソフトにDNA plotter、genome diagram、snap geneなどがあるが、機能の豊富さではcircosが抜きん出てる。ユーザー数も多く、フィードバックも得やすい。難点として、機能が豊富な分使いこなすにはかなり勉強が必要。

 

幸い、素晴らしく丁寧に解説されたトレーニングページが用意されている。http://circos.ca/documentation/tutorials/ を見て、1歩ずつ理解していけばよい。

 

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

brew install circos

使い方

circos.confファイルを作成して

circos -conf circos.conf

 

トレーニング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に続く。

 

 

 

 

 

 

 

 

 

 

indelの検出ツール

バクテリアでパフォーマンス比較したペーパーが出ている。

 

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

 

 

Break-dancer-max Chen et al. (2009)

PEM法

 

macOSX環境でのインストールも可能。ペアードエンドリードの挙動から構造変化を予測する。インストールはめんどくさかった、cpanでたくさんのperlモジュールをインストールする必要があった。 

 

マニュアルで以下のように説明されている。

There are two sequential steps to run BreakDancer. First, the “bam2cfg.pl” script analyzes the first few thousand reads in each input BAM file and generates a configuration file containing statistics for each read group (RG; see Unit 11.10). Second, the “breakdancer-max” program uses these statistics to identify and cluster DRPs, ultimately yielding a list of putative SVs

上の文章にあるようにperlスクリプトで.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がpprなデータだと解析できない (ランした時にpoor quality dataと出る)。 どうしてもランしたければ、オプション"-v"をdefaultの1から2-10まであげると一応ランはできる。

 

#構造変化のコール

 breakdancer-max output.cfg > output.txt

 

 

 

LUMPY layer et al. (2014)

split-read approaches  + read-pair technologies + CNV

 

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

brew install homebrew/science/lumpy-sv

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

f:id:kazumaxneo:20170327095430j:plain

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

 

 下のような流れでコールする。

bwa index -a bwtsw 6803.fasta
bwa mem -R "@RG\tID:X\tLB:Y\tSM:Z\tPL:ILLUMINA" 6803.fasta T1second_R1.fastq T1second_R2.fastq |samtools view -S -b - > sample.bam
samtools view -b -F 1294 sample.bam > sample.discordants.unsorted.bam
samtools view -h sample.bam | ./extractSplitReads_BwaMem -i stdin | samtools view -Sb - > sample.splitters.unsorted.bam
samtools sort sample.discordants.unsorted.bam sample.discordants
samtools sort sample.splitters.unsorted.bam sample.splitters
lumpyexpress -B sample.bam -S sample.splitters.bam -D sample.discordants.bam -o sample.vcf

こんな感じで使う。

途中のextractSplitReads_BwaMemが見つからなければ、LUMPYのソースコードに入っている。それにパスを通すか解析ディレクトリにコピーして(上のように)、使う。

 

 

Hydra-sv

PEM法

マニュアル: 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

 

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

 

 

Pindel Ye et al. (2009)

split-read approaches.

 

macOSX環境ではエラがー大量に出てインストールできなかったのでcentOSサーバーにインストールした。ショートリードのペアードエンドデータからindelを予測するツール。下は公式ページから。

f:id:kazumaxneo:20170327095850j:plain

黄色はdouble strand breakの修復過程で新規に追加されたリファレンスに存在しない配列。pinldeの新バージョンではこれも含めて予測できると記載されている。

http://gmt.genome.wustl.edu/packages/pindel/install.html にインストールの仕方が書いてある。

linuxでroot環境で

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

 

パスが通ってないければbashrcかbash_profileに記載してsourceしておく。

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

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

ランに必要なものはリファンレス配列(.fasta)とそのindexファイル(.fasta.fai)。またそのリファレンスにmapして作ったbamファイルとそのindexファイル(.bam.bai)。

bamファイルについては-fの引数で指定するsimulated_config.txtファイルの中に記載する。 テキスト内には以下のことが書かれている。

simulated_sample_1.bam  250     SAMPLE1

simulated_sample_2.bam  250     SAMPLE2

simulated_sample_3.bam  250     SAMPLE3

上のコマンドを実行すると、1分くらいでたくさんのファイルが出力される。indel予測結果はoutput_RPファイル内に保存されているみたい。

これを以下のように変更した。

R1R2_sorted.bam  600     SAMPLE1

この600はインサートサイズ。当たり前だがインサートサイズがリード長より短いとかだと動かない。fastaとbamのindexファイルも準備する。

samtools faidx input.fasta
samtools index R1R2_sorted.bam

 ラン。

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

 4Mのゲノムで10分くらいかかった。

 

 

DELLY Rausch et al. (2012)

split-read approaches. 

 

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

と書いてある。ショートリードのペアードエンドデータや、メイトペアのデータを使いsplit readのアプローチでdeletionなどを予測する。コピーナンバーが変わるリピートのdeletionや、depthが少ないデータにも強いらしい。

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

テストランは

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

となっているが、-xはコールを除外するセントロメアポジションを記載したファイル。今回コールを除外する領域はなかったため消す。

./delly call -t DEL -o  delly.bcf -g 6803.fasta R1R2_sorted.bam

 上のコマンドでbinaryファイルが出力されてくるので変換。フィルターコマンドは実行していない。

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

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

./delly call -t DEL -o del.bcf -g 6803.fasta R1R2_sorted.bam &
./delly call -t DUP -o dup.bcf -g 6803.fasta R1R2_sorted.bam &
./delly call -t INV -o inv.bcf -g 6803.fasta R1R2_sorted.bam &
./delly call -t TRA -o tra.bcf -g 6803.fasta R1R2_sorted.bam &
./delly call -t INS -o ins.bcf -g 6803.fasta 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)

assembly based approaches. 

 

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

解析速度は一番速いかもしれない。10秒以内に終わり、かつtrue indelのコール率も高かった。 

 

 

fermikit Li (2015)

assembly based approaches. 

 

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というファイルができる。prefix.flt.vcf.gzの中にSV予測結果が保存されている。

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

 

 

 

 

 

GATK

split-read approaches. 

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

 

 

 

VarScan Koboldt et al. (2009)

gapped alignment approaches. 

 

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

brew install VarScan

ソート済みbamファイルをpileupする。

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

SNVをコール

VarScan pileup2snp mpileup > output.txt

indelをコール

VarScan pileup2snp mpileup > output.txt

 

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

indelをコールするなら

samtools mpileup -f 6803Kazusa.fasta 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.fastq

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

-j: スレッド数

-0: アウトプット

 

解析が終わるとhtmlのまとめファイルができるのが特徴。/outputの中のsummary.htmlを開くと、

f:id:kazumaxneo:20170426164254j:plain

 

coverage

f:id:kazumaxneo:20170426164310j:plain

mutation prediction

f:id:kazumaxneo:20170426164413j:plain

 indelの抜け漏れはこの中で一番少なかった。インフォマティクスのツールに不慣れなユーザーにも使いやすいと思われる。

 

 

Scalpel Narzisi et al. (2014)

assembly based approaches. 

  

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

bcftoolsを使ってbamファイルから変換することができる(SEQanswers)。

bamToBed -i R1R2_sorted.bam > R1R2_sorted.bed 

 ただし現在のところ./scalpel-discoveryを走らせると以下のエラーがでて止まってしまう。

Command failure: extract SAM header (/home/user/scalpel-0.5.3/bamtools-2.3.0/bin/bamtools header -in /home/user/scalpel-0.5.3/outdir/bamfile.bam > /home/uesaka/scalpel-0.5.3/outdir/header.txt)...

 

色々検討したが、まだ解決していない。

 

 

Breakseek Zhao et al. (2015)

breakpoint-based algorithm

ダウンロードリンク

 

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が用意されているので動作テストにはそれを使う。

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

 

 

samができたら、以下のコマンドでsamをクロモソームごとに分割する。

python sam_prep.py -f test.sam

 

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

python breakseek.py -f testS.sam -r test.fa

実行時間はかなり長く、3Mのバクテリアで20分以上要した。また、以下のメッセージ

gathering statistics on the dataset...

が出てから、解析途中は何もlogが出てこない。

初めはランに失敗したと思っていたが、強制終了せず気長に待ったところ、2時間半かかった。100kbpのプラスミドを使っただけでこれは長いなあ...

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

引用

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

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/

 

DELLY: 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/

 

FreebayesHaplotype-based variant detection from short-read sequencing

https://arxiv.org/pdf/1207.3907.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