Trimmomatic の keepBothReads 設定

はじめに

生命科学系の博士課程の学生です。
Trimmomatic を使っていて気づいた点があり、これに関して日本語での記事が見つからなかったので投稿してみました。誰かの助けになればと思います。

Trimmomatic

Trimmomatic は次世代シーケンサーのシーケンシングデータのクオリティコントロールに用いられるツールの一つです。
シーケンシングライブラリーには目的の配列に加えて、シーケンシングそのものに必要なアダプターシーケンスやインデックス配列が含まれています。これがデータに含まれている場合、下流の解析の効率が下がることがあります。Trimmomatic はこのような余分な配列をトリミングすると同時に、シーケンシングクオリティが低いリードもトリミングしてくれます。
基本的な使い方はこちらのウェブサイトなどが参考になります。

Issue

ライブラリサイズに対してリード長を長めにとったペアエンドのシーケンスデータに対し、下記のようなオプションで Trimmomatic を使うと Both Surviving のリードが著しく低くなります。Reverse Only Surviving が低く出るのも特徴的です。
筆者の手元では、300-500 bp のライブラリに対して 100-150 bp で読むと起こりました。

trimmomatic \
    PE \
    -threads 4 \
    -phred33 \
    ${dir}/${line1}_L1_1.fq.gz \
    ${dir}/${line1}_L1_2.fq.gz \
    result/trimming/paired_${line1}_L1_1_trim.fastq \
    result/trimming/unpaired_${line1}_L1_1_trim.fastq \
    result/trimming/paired_${line1}_L1_2_trim.fastq \
    result/trimming/unpaired_${line1}_L1_2_trim.fastq \
    ILLUMINACLIP:${path to trimmomatic}/trimmomatic/adapters/NexteraPE-PE.fa:2:10:10 \
    MINLEN:30 \
    LEADING:3 \
    TRAILING:3 \
    SLIDINGWINDOW:4:15

ILLUMINACLIP: Using 1 prefix pairs, 4 forward/reverse sequences, 0 forward only sequences, 0 reverse only sequences
Input Read Pairs: 50748982 Both Surviving: 22590181 (44.51%) Forward Only Surviving: 27764030 (54.71%) Reverse Only Surviving: 13682 (0.03%) Dropped: 381089 (0.75%)

結論

ILLUMINACLIP の引数の最後の keepBothReads を True にしてください。

trimmomatic \
    PE \
    -threads 4 \
    -phred33 \
    ${dir}/${line1}_L1_1.fq.gz \
    ${dir}/${line1}_L1_2.fq.gz \
    result/trimming/paired_${line1}_L1_1_trim.fastq \
    result/trimming/unpaired_${line1}_L1_1_trim.fastq \
    result/trimming/paired_${line1}_L1_2_trim.fastq \
    result/trimming/unpaired_${line1}_L1_2_trim.fastq \
    ILLUMINACLIP:${path to trimmomatic}/trimmomatic/adapters/NexteraPE-PE.fa:2:10:10:8:TRUE \
    MINLEN:30 \
    LEADING:3 \
    TRAILING:3 \
    SLIDINGWINDOW:4:15

理由

そもそもどうしてアダプター配列がコンタミするのでしょうか？ライブラリーの調整方法にも依存しますが、大概ライブラリの構成はこのように中央のインサートにアダプター、インデックスが結合する形になっていて、インサート配列の端から読むようになっています。下記の図のようにインサート長に対してリードの長さが短ければシーケンスデータにアダプターなどは含まれないでしょう。

しかし、インサート長がリード長に対して短い場合はどうでしょうか？

あ。。。リードがインサートを超えてアダプターやインデックスまで突き抜けてしまっています。このような時にリード内にアダプター配列がコンタミします。蛇足ですがリードがあまりにも長くライブラリの先端まで読み切ってしまうと、その先は polyG として出力されていました (筆者の経験です。シーケンサーの仕様によると思います)。

さて、このような時、片方のリードだけでインサート領域を丸々カバーできてしまっていますよね？この場合、Forward と Reverse のリードに含まれている情報は丸々重複しています。
Trimmomatic の話に戻りますが、このようなとき、言い換えれば "ペアエンドで Forward リードにアダプター配列が検出された場合"、 Trimmomatic は 対応する Reverse のリードを抹消します。
では、相方が抹消された Forward のリードは一体どこへ？

"unpaired" な fastq として出力されます！！！
何も気づかずに下流の解析に paired なデータだけ用いてしまった場合、かなりのリードが失われてしまうことになります。。

これを防ぐのが keepBothReads 設定です。ILLUMINACLIP の引数の最後の値です。

ILLUMINACLIP:<fastaWithAdaptersEtc>:<seed mismatches>:<palindrome clip
threshold>:<simple clip threshold>:<minAdapterLength>:<keepBothReads>

ここに TRUE と入力するとペアエンドで forward 配列にアダプターが含まれていても Reverse のリードが捨てられずに済みます。

要するに

ペアエンドで Trimmomatic を使うとき、特にライブラリサイズに対して長めのリード長をとっているときには keepBothReads 設定をオンにしましょう。
さもなくばアダプターが含まれたリードが問答無用で "unpaired" に入れられてしまいます。

最後に

一般的にライブラリを過剰な長さで読んでしまうことはあまりないと思うのですが、ライブラリそのものが短いケースや異なるライブラリを同一フローセルで読むときなどは気を付けましょう。