Illuminaのシークエンスデータ（fastqファイル）の前処理についてのパイプラインです。
fastqファイルとは
fastqで必要な前処理
- +1bp（extra 1bp）の除去
- アダプター配列の除去
- 低クオリティ塩基および低クオリティリードの除去
- Iluminaのベースコールの変換（※生データのbclファイルからfastqファイルを生成すること。通常、bcl2fastqを使う）の設定によっては、上の1-3は実施済みの場合がある。ベースコール時には1-3の処理はさせずに自前で実施することをオススメする。外部委託などの場合は、生データのbclファイルもしくは非処理のfastqファイルを納品してくださいと伝える。ベースコールで特段の処理がされていないかも含めて、常にFastQCでのquality controlを実施すること。
fastqファイルのquality control（品質管理）
- FastQCが定番中の定番。基本的には、最初のfastqならびに処理過程で生成されるすべてのfastqファイルにFastQCをかけて、変なことが起きてないか確認するべきだろう。
- Javaで動く。Linux/Win/Macすべてに対応。GUIでfastqファイルをドロップしてもいいし、CUIで実施してもいい。結果はhtmlファイルで出力される。
- 基本の書式：
SAMPLE="hoge"
fastqc ${SAMPLE}.fastq
- デフォルトだと、50bp以上のデータが丸められてしまう。また、テキストベースの結果がzip圧縮されている。データを丸めず、圧縮しないで、
hoge_1.fastq.gzとhoge_2.fastq.gzというペアエンドのデータにFastQCをかける場合、次の通り。
SAMPLE="hoge"
fastqc -t 8 -o ${SAMPLE}_1 --extract --nogroup ${SAMPLE}_1.fastq.gz
fastqc -t 8 -o ${SAMPLE}_2 --extract --nogroup ${SAMPLE}_2.fastq.gz
- 入力ファイルはgz圧縮されていてもされてなくてもいい。8スレッド使用（-t）、"hoge_1"または"hoge_2"というディレクトリに出力（-o）、結果を圧縮しない（--extract）、1bpごとに出力（--nogroup)という意味。
- 次のことは表にして起こしておくといい（論文に書くときに便利）。
- 塩基長
- リード数
- GC含量
- Sequence duplication levels
- Overrepresented sequences
- htmlファイルからコピペしてもいいが、一緒に生成されるfastqc_data.txtから抜き出すこともできる。
+1bp（extra 1bp）の除去
- Illuminaのシークエンサーは、リードの最後の1bpのquality scoreを正確に計算できないため、希望のリード長+1bpでシークエンスするのが慣例。例えば150bpを読みたい場合は、151 cyclesでシークエンスする。最後の+1bp（extra 1bp）の塩基自体が周辺の塩基に比べて極端に間違っているわけではないが、正確なクオリティスコアが付与されないため、何よりもまずトリミングすることをオススメする。
- なお、bcl2fastqの時点でトリミングしている場合がある。特に外部委託データなどの場合、extra 1bpをトリミングしているかどうか不明なので、必ずFastQCを納品データにかけて、extra 1bpがあるかどうかを確認すること。
- トリミングはいろいろなソフトでできるが、Trimmomaticが便利。Javaで動く。2021年12月時点でv0.39が利用できる。
- スレッド数8、150bp PE（151 cycles x2）の場合は次の書式。
SAMPLE="hoge"
java -jar trimmomatic-0.39.jar PE -threads 8 \
${SAMPLE}_1.fastq.gz \
${SAMPLE}_2.fastq.gz \
${SAMPLE}_fw_paired_crop.fastq \
${SAMPLE}_fw_unpaired_crop.fastq \
${SAMPLE}_re_paired_crop.fastq \
${SAMPLE}_re_unpaired_crop.fastq \
CROP:150
-
hoge_fw_unpaired_crop.fastqとhoge_re_unpaired_crop.fastqは相手のいないリードを格納するファイルのため、今回の150bpにトリミングしただけでは空っぽ（0 byte）で出力される。
アダプター配列の除去
- たとえば150bp PEで読んだ場合、インサート配列（読みたいターゲット部分）が150bp未満だと、本来のターゲットではないアダプター配列までシークエンスされてしまう。3'末端に出現しているアダプター配列を除いておく必要がある。
- アダプター配列は、ライブラリの作成試薬やシークエンス機器などによって異なる。不明な場合は、Illumina Adapter Sequencesで調べておく。
- FastQCのOverrepresented sequencesでも、アダプター配列が検出されている場合があるので、ダブルチェックするといい。
- TrimmomaticのILLUMINACLIPでもアダプター配列の除去ができるが、個人的にはSkewerが直感的に使いやすいのでオススメ。2021年12月時点でv0.2.2が利用できる。Jiang et al. BMC Bioinformatics 2014。
- 150bpにCROP済のPEデータから、スレッド数8、最小リード長50bpで、3'末端側のアダプター配列を除去する場合は次の書式。
SAMPLE="hoge"
skewer-0.2.2-linux-x86_64 \
${SAMPLE}_fw_paired_crop.fastq \
${SAMPLE}_re_paired_crop.fastq \
-m pe \
-l 50 \
-o ${SAMPLE} \
-t 8
- 一緒に出力される
hoge-trimmed.logを確認して、どれくらいの割合でトリミングされたかを必ず確認すること。適切なインサートサイズのライブラリを作成できて、適切なサイクル数でのシークエンスが選択できたかの指標になる。 - skewerのデフォルトのアダプター配列は、Illuminaの長く使われてきたTruSeqなどのアダプター配列（IDT for Illumina–TruSeq DNA and RNA UD Indexes）が設定されており、forward側が
AGATCGGAAGAGCACACGTCTGAACTCCAGTCACでreverse側がAGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA。 - たとえば2021年からリリースされたllumina DNA PCR-Free Prep, Tagmentationは、タグメンテーションの過程で、
CTGTCTCTTATACACATCTとATGTGTATAAGAGACAの2種類のアダプター配列をもったリードが混在するため、マニュアルでアダプター配列を指定し、二段階でskewerをかける必要がある。
SAMPLE="hoge"
skewer-0.2.2-linux-x86_64 \
${SAMPLE}_fw_paired_crop.fastq \
${SAMPLE}_re_paired_crop.fastq \
-x CTGTCTCTTATACACATCT \
-y CTGTCTCTTATACACATCT \
-m pe \
-l 50 \
-o ${SAMPLE} \
-t 2
skewer-0.2.2-linux-x86_64 \
${SAMPLE}-trimmed-pair1.fastq \
${SAMPLE}-trimmed-pair2.fastq \
-x ATGTGTATAAGAGACA \
-y ATGTGTATAAGAGACA \
-m pe \
-l 50 \
-o ${SAMPLE}2 \
-t 8
- 先行研究などのサードパーティデータを大量に扱わなくてはならず、ひとつひとつアダプター配列を調べるのが大変な場合、EARRINGSというソフトに入れてみるのも一つの手。データ全体からアダプター配列を予測し、除去まで自動的に進行してくれる。非常に便利だけど、間違える可能性もゼロではないので、調べられる場合は調べてskewerをかけた方がいい。Wang et al. 2021 Bioinformatics。
低クオリティ塩基および低クオリティリードの除去
- クオリティスコアの悪い塩基を除去し、さらにクオリティスコアがリード全体にわたって悪い場合はそのリードそのものを除去する。
- Trimmomaticを使う。たとえば、5末端側（LEADING）ならびに3末端側（TRAILING）で、端からクオリティスコアが20未満が続く塩基をすべて削除する。さらに、連続する4塩基（4bp window）のクオリティスコアの平均が20未満になった場合（SLIDINGWINDOW）、そのウィンドウ以降の塩基をすべて除去する。そして、最小塩基長を50bpとする場合、次の書式になる。
SAMPLE="hoge"
java -jar trimmomatic-0.39.jar PE -threads 8 \
${SAMPLE}-trimmed-pair1.fastq \
${SAMPLE}-trimmed-pair2.fastq \
${SAMPLE}_fw_paired_trimmed.fastq \
${SAMPLE}_fw_unpaired_trimmed.fastq \
${SAMPLE}_re_paired_trimmed.fastq \
${SAMPLE}_re_unpaired_trimmed.fastq \
LEADING:20 \
TRAILING:20 \
SLIDINGWINDOW:4:20 \
MINLEN:50
- ここで、SLIDINGWINDOWは、たとえばこのような状況が想定される。連続する塩基のスコアが19,19,24,20,20,15,30,30,30となった場合、6番目の低品質の塩基（Q=15）はLEADINGでもTRAILINGでも除去されない。しかし、SLIDINGWINDOWを設定することで、前からの4塩基のウィンドウ（24,20,20,15）の平均値は20を切るので、このウィンドウ以降の塩基を除去できる。
- Trimmomaticが終了すると、何パーセントのリードが、forward側、reverse側、両側で除去されたかが示されるので必ず確認すること。LEADING、TRAILING、SLIDINGWINDOW、MINLENを保守的にすればするほど、除去される割合が増えていくので、その後の解析の目的に合わせて適切に設定する。ゲノムDNAなのかmRNAか、de novo assemblyなのかmappingなのか、一倍体なのか多倍体なのか、リピート配列を解析するかしないかなどで、要求される塩基のクオリティやリード長は変わってくる。
まとめ
- TruSeqアダプター配列を使用したライブラリを、150bp PE（151+151 cycles）でシークエンスして得られた
hoge_1.fastq.gz（forward側）とhoge_2.fastq.gz（reverse側）を処理するときの書式。（スレッド数8、最低リード長50bp、Qスコアの目安20）
SAMPLE="hoge"
java -jar trimmomatic-0.39.jar PE -threads 8 \
${SAMPLE}_1.fastq.gz \
${SAMPLE}_2.fastq.gz \
${SAMPLE}_fw_paired_crop.fastq \
${SAMPLE}_fw_unpaired_crop.fastq \
${SAMPLE}_re_paired_crop.fastq \
${SAMPLE}_re_unpaired_crop.fastq \
CROP:150
skewer-0.2.2-linux-x86_64 \
${SAMPLE}_fw_paired_crop.fastq \
${SAMPLE}_re_paired_crop.fastq \
-m pe \
-l 50 \
-o ${SAMPLE} \
-t 8
java -jar trimmomatic-0.39.jar PE -threads 8 \
${SAMPLE}-trimmed-pair1.fastq \
${SAMPLE}-trimmed-pair2.fastq \
${SAMPLE}_fw_paired_trimmed.fastq \
${SAMPLE}_fw_unpaired_trimmed.fastq \
${SAMPLE}_re_paired_trimmed.fastq \
${SAMPLE}_re_unpaired_trimmed.fastq \
LEADING:20 \
TRAILING:20 \
SLIDINGWINDOW:4:20 \
MINLEN:50
終わりに
今日はここまで。サンプルデータで実際にかけてみた例を後日掲載できればと思っている。（2022年1月4日、了。）
- 前処理をおこなったfastqファイルでマッピングなどをするときは、『BWAとGATKでゲノムをマッピングする 』(https://qiita.com/hayatak/items/da5d1cb096216fef87ef) を参照。（2022年1月 5日、追記。）