CLC genomics workbench (7.0)でRNA seq解析 - macでインフォマティクス

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を使って調べてみる。

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

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

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

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

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

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

・CLC7を起動

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

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

f:id:kazumaxneo:20170418140342p:plain

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

・最後にFile -> Save で保存。

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

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

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

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

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

３、モチーフ検索ツールでサンプリングした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。

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

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

・~ 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以外のコマンドツールではどうか？

１、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に変えて