CLCで行う前提でワークフローを書いてみる。
アダプター配列の確認
アダプター配列はシーケンス->画像からのシグナル取得 -> ベースコールファイル(.bcl) -> FASTQ変換の過程で自動除去されるようだが、インサートが短いペアリードファイルなどでは3'側にアダプターが存在している可能性がある。例えば以下の
https://www.researchgate.net/post/PCR_primers_trimming_for_MiSeq
でそのことが論じられている。アダプター残存による影響は解析によって変わってくる。例えばCLCの資料の
には、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。パラメータの詳細はこちらを参考
・アダプター配列分これを繰り返す。
・最後に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は先ほど1で作成したリストを選択。
・Include negative strandにチェック
・次のページでCreate TablesとAdd anotations to sequencesにチェック
・resultはsaveを選択しfinish。
・解析が終わると、モチーフが見つかったリードを記載したリストができる。
4、アダプター配列が見つかるか確認
・~ sampledファイルの中身をダブルクリックで開き、右側のリストからAnnotation typeを選択する。分かりにくいが写真の青色のMotifというアイコンの右上にAnnotation typeというタブがある。
・青色のMotifボックスにチェックをつけるとモチーフが見つかった部位が矢印で表示される。今回は人工的に下の画像の矢印部位が100%マッチするようにモチーフ配列を設計した。
同じことをRでやるなら
library(qrqc)
fastq <- readSeqFile("input.fastq")
basePlot(fastq)
先頭だけATGC頻度が極端なのでアダプターと考えられる。
kmerKLPlot(fastq) #6-mer時の塩基出現頻度
高頻出の配列がアダプターということになる。
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というのができるので、開いてみる。
先ほどの青色矢印部分の配列が除去されている。
アダプターの上流側に配列がある場合、それもまとめて除去されている。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塩基だけクオリティスコアが悪くても必ずしもトリミングはされないアルゴリズムを用いている。
#アダプター配列の除去
toolbox → NGS 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
TOP10の長さ
CLC genomics workbench 7(default condition)
TOP10の長さ
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のペアデータである。他のデータでも同様の傾向が出るのかはわからない。