macでインフォマティクス

macでインフォマティクス

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

CLC genomics workbench (7.0)でRNA seq解析

CLCで行う前提でワークフローを書いてみる。

 

アダプター配列の確認

アダプター配列はシーケンス->画像からのシグナル取得 -> ベースコールファイル(.bcl) -> FASTQ変換の過程で自動除去されるようだが、インサートが短いペアリードファイルなどでは3'側にアダプターが存在している可能性がある。例えば以下の

https://www.researchgate.net/post/PCR_primers_trimming_for_MiSeq

でそのことが論じられている。アダプター残存による影響は解析によって変わってくる。例えばCLCの資料の

https://jp.illumina.com/content/dam/illumina-marketing/apac/japan/documents/pdf/2015_techsupport_session10.pdf

には、BWA MEMでマッピングした時には影響がないが(たぶんクリッピングされるから)、同じBWA でもショートリード向けのBWA Backtrackでは影響が多いと書かれている。アセンブルする場合も影響はわからない。de brujin graphを書くなら、影響は少ない可能性も十分ありえる。

 

 

アダプターが残存しているかCLCを使って調べてみる。

アダプター配列が存在しているか確認するには(オフィシャルサポートリンク)

 

 1、アダプターリストの作成

 2、fastqをダウンサンプリング

 3、モチーフ検索ツールでダウンサンプリングfastqファイル中のモチーフを検索

 4、アダプター配列が見つかることが確認できたら、アダプタートリミングツールでアダプターを除去

 

1、アダプターリストの作成

CLC7を起動

File -> New -> Trim motif list を選択。

ウィンドウが出てきたら下のAdd Rowボタンを押してアダプター配列情報を記入していく。Mismatchコストは初期値のままにしてFinish。パラメータの詳細はこちらを参考

f:id:kazumaxneo:20170418140342p:plain

 

アダプター配列分これを繰り返す。

最後にFile -> Save で保存。

 

 

2、fastqの一部データのサンプリング

・モチーフ検索ツールは1000リードまでしか対応していないみたいなので、FASTQのリードを1000だけランダムサンプリングする。手順は

・前もって読み込んだfastqファイルを右クリック -> Toolbox -> NGS core tool -> Sample reads

リード数を1000に指定してRun

10秒くらいでサンプリングは終わる。~ sampledというファイルが作られているはず。

 

 

3、モチーフ検索ツールでサンプリングしたfastqファイル中のモチーフを検索

・~ sampledファイルを右クリック -> Toolbox -> Classical Sequence Analysis -> General Sequence Analysis -> Motif search を選択。

・Motif Search TypeはMotif Listを選択。

・Motif listは先ほどで作成したリストを選択。

・Include negative strandにチェック

f:id:kazumaxneo:20170418141636p:plain

 

・次のページでCreate TablesとAdd anotations to sequencesにチェック

・resultはsaveを選択しfinish。

 ・解析が終わると、モチーフが見つかったリードを記載したリストができる。

 

 

4、アダプター配列が見つかるか確認

・~ sampledファイルの中身をダブルクリックで開き、右側のリストからAnnotation typeを選択する。分かりにくいが写真の青色のMotifというアイコンの右上にAnnotation typeというタブがある。

f:id:kazumaxneo:20170418142416p:plain

 

・青色のMotifボックスにチェックをつけるとモチーフが見つかった部位が矢印で表示される。今回は人工的に下の画像の矢印部位が100%マッチするようにモチーフ配列を設計した。

f:id:kazumaxneo:20170418144401p:plain

 

同じことをRでやるなら

library(qrqc)

fastq <- readSeqFile("input.fastq")

basePlot(fastq)

f:id:kazumaxneo:20170418154709j:plain

先頭だけATGC頻度が極端なのでアダプターと考えられる。

 

kmerKLPlot(fastq)  #6-mer時の塩基出現頻度

f:id:kazumaxneo:20170418155116j:plain

高頻出の配列がアダプターということになる。

 

k-mer長を変えるなら

fq <- readSeqFile("1.fq", kmer=TRUE, k=15) #k=15の時 

kmerKLPlot(fq) #グラフ描画

 

 

 

アダプター配列のトリミング

上の画像の部位のトリミングを行ってトリミングツールの挙動を確かめてみる。以下の手順で進める。

File -> New -> Trim Adapter list を選択。

・アダプター配列を順番に記載してsave。リストファイルを作る。(上のモチーフ検索リスト作成と同じ要領)

・NGS core tools -> Trim sequences を選択。

・いまアダプタートリミングだけ行いたいので、クオリティトリミングや5' or 3'末端トリミングのチェックを全て外す。先ほど作ったリストファイルを読み込ませ、Finishでラン開始。

 

終わったら、FASTQファイルのtrimmedというのができるので、開いてみる。

f:id:kazumaxneo:20170418145842p:plain

先ほどの青色矢印部分の配列が除去されている。

アダプターの上流側に配列がある場合、それもまとめて除去されている。5'側がまとめて除去されるので、index配列が分かればその上流も含めて全て除去できるのは都合がいい。ただし、配列の途中にアダプターが出てくる場合、強制的に5'側がトリムされる仕様は不都合が出てくるかもしれない。例えばインサートが短くてペアリードが調整したライブラリの末端までシーケンスされていた場合(特にsmall RNAシーケンスなど)、3'側にアダプター配列が出てくる可能性がある。その場合、下流側の配列を除く必要があり、強制的に5'側が削除されるのはまずい。

 

 

CLC以外のコマンドツールではどうか?

1、cutadaptの場合

 

@DRR057320.1 MG00HS11:593

TTAGGTTTCTACAAAATGAAGATTTCGAAAGTTTATCAAAACAAAGAATCT

+

@@@DDA=DHFFBDA:EFHIII@FHHHEG?+??CFEDEDHD:FHCG;DCGIF

上のリードの太字部分がトリム対象なら、cutadapt -gコマンドで

cutadapt -g TACAAAATGAAGATTTCG 1.fq > 1_trimmed.fq #5'側トリミング

 

トリム後、アダプターの5'側も除去される。

@DRR057320.1 MG00HS11:593

AAAGTTTATCAAAACAAAGAATCT

+

G?+??CFEDEDHD:FHCG;DCGIF

こうなる。3'側をトリムするなら-gを-aに変えて

cutadapt -a TACAAAATGAAGATTTCG 1.fq > 1_trimmed.fq #3'側トリミング

リードを見ると

@DRR057320.1 MG00HS11:593

TTAGGTTTC

+

@@@DDA=DH

今度はアダプターとその下流側が除去されている。

 

 

 

 

 

 

 

 

トリミング

詳細はこちらを参照した。 

CLC genomics workbench のトリムは累積のクオリティスコアが一定以上の領域が続いた場合に取り除くらしい。つまり1塩基だけクオリティスコアが悪くても必ずしもトリミングはされないアルゴリズムを用いている。

#アダプター配列の除去

toolboxNGS core tools → Trim sequences を選択。

 

de novo アセンブリ

trinityの結果と比較する。

Trinity --seqType fq --left S1_L001_R1.fastq --right S1_L001_R2.fastq --max_memory 10G --CPU 18

 

>行をカウントしてコンティグ数を見積もる。

Trinity 1656 contigs

CLC genomics workbench 7 1956 contigs

rnaspades 1952 contigs

 

数についてはTrinityが一番少なく、rnaspadesが一番多い。

長さの分布やlongest conntigを見るとどうか?

 

trinity

f:id:kazumaxneo:20170415141857p:plain

TOP10の長さ

f:id:kazumaxneo:20170415142259p:plain

 

 

CLC genomics workbench 7(default condition)

f:id:kazumaxneo:20170415142525p:plain

TOP10の長さ

f:id:kazumaxneo:20170415142214p:plain

 

 

 

 

rnaspades

rnaspades.py -t 30 -k auto -o data -1 1_R1.fq -2 1_R2.fq )

TOP10の長さ

Name length
NODE1 19799
NODE2 18969
NODE3 17335
NODE4 15793
NODE5 15760
NODE6 13466
NODE7 13232
NODE8 12802
NODE9 12596
NODE10 12065

 

 

top10の長さはrnaspadesとtrinityが優れている感じ。今回のデータはMiseq 301bpのペアデータである。他のデータでも同様の傾向が出るのかはわからない。