TrimmomaticによるBGISEQデータのアダプター配列除去方法について

こんにちは、BGIシーケンサーからの出力結果で問い合わせが増えていますので、ブログで紹介します。 RNAseq解析ではシーケンス時の読み込まれたアダプター配列の除去を行います。 MGI社が提供するBGISEQ/DNBSEQ/MGISEQは、特別にアダプター配列を指定する必要があります。

以下のマニュアルにadapter配列が示されていますので、FASTAファイルを作成してForward_filterとReverse_filterを表記します。「/1」or「/2」を表記することで、ForwardとReverseそれぞれのトリミング処理を指定できます。複数指定することも可能です。 https://en.mgitech.cn/Download/download_file/id/71

atapters.fa

>Forward_filter/1 
AAGTCGGAGGCCAAGCGGTCTTAGGAAGACAA 
>Reverse_filter/2 
AAGTCGGATCGTAGCCATGTCGTTCTGTGAGCCAAGGAGTTG

次にTrimmomaticを実行して、作成したadapters.faのパスを指定します。 Trimmomaticをインストール前の場合は、インストールします。
USADELLAB.org - Trimmomatic: A flexible read trimming tool for Illumina NGS data

trim.sh

#!/bin/bash

java -jar /home/centos/amelieff/soft/Trimmomatic-0.33/trimmomatic-0.33.jar \
    PE                    \  # ペアエンドかシングルエンドか指定
    -threads 4            \  # 実行時のthreads数を指定
    -phred64              \  # FASTQのクオリティスコアの指定。BGIデータでは必須で指定し、基本的に64を選択
    -trimlog log.txt      \  # logファイルの指定
    path/to/input_R1.fastq.gz  path/to/input_R2.fastq.gz  \   # 入力ファイル１と２を指定
    paired_input_R1.fastq.gz  unpaired_input_R1.fastq.gz  \   # 入力ファイル１のペアエンドと片方のみ残った場合の出力ファイル
    paired_input_R2.fastq.gz   unpaired_input_R2.fastq.gz \   # 入力ファイル２のペアエンドと片方のみ残った場合の出力ファイル
    ILLUMINACLIP:adapters.fa:2:30:10      \  # ここで作成したアダプター配列のパスを指定
    LEADING:20             \  # リード先頭からクオリティスコアが指定数未満切り捨て
    TRAILING:20            \  # リード末端からクオリティスコアが指定数未満切り捨て
    SLIDINGWINDOW:4:15     \  # ウィンドウサイズ:平均クオリティスコアを指定して、値未満切り捨て
    MINLEN:36              \  # 値を満たさないリードを除去

上記ファイルを実行した時にクリッピングシーケンスが表示されます。

Using Long Clipping Sequence: 'AAGTCGGAGGCCAAGCGGTCTTAGGAAGACAA'
Using Long Clipping Sequence: 'AAGTCGGATCGTAGCCATGTCGTTCTGTGAGCCAAGGAGTTG'
ILLUMINACLIP: Using 0 prefix pairs, 0 forward/reverse sequences, 1 forward only sequences, 1 reverse only sequences

以上でアダプター配列の除去が可能となります。