NGS で読んだ配列にはアダプターが付いています。解析する前に生データからこのアダプターを取り除いたり、低クオリティのリードは除去したりする必要があります。
今回はその工程を紹介します。
アダプターってなに、という人はこのサイトを見てみてね。
今回使うソフトウェアはこちら。
公式マニュアルはこちら
http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/TrimmomaticManual_V0.32.pdf
Trimmomatic について日本語のブログで紹介してくれている人がほかにもいるので、データの数が少ない場合には以下のブログを見てもらうので十分かもしれません。
このブログではRAD-seq やMIG-seq で大量の個体データを扱う場合に便利な方法を紹介していきます。
この方法を使うときはアナコンダさん経由でTrrimomatic を入れてあげる必要があるので、まだ入れていないという人はまずアナコンダを入れてください。
【Linux/ubuntsu】Windows10でバイオインフォマティクス系ソフトウェアを使うための準備 - 研究メモ (hatenablog.com)
この記事の３段階目に入れ方は書いてあります。
#install
conda install -c bioconda trimmomatic
#main script
trimmomatic PE -threads スレッド数 \ #PE はペアエンド構造が壊れたファイルを捨てる設定
-phred33 \ #raw read のクオリティを記述している言語。シーケンサーによって違うので、データが来たら調べて書き換える。基本は33（これの調べ方はわかりません）。
-trimlog ログ名.txt
生データ名(R1).fastq.gz 生データ(R2).fastq.gz \ #インプット生データ
アウトプット名_pair_R1.fastq アウトプット名_unpair_R1.fastq \ #アウトプットの指定
アウトプット名_pair_R2.fastq アウトプット名_unpair_R2.fastq \ #アウトプットの指定
ILLUMINACLIP:adapter.fasta:2:30:10 \
#ILLUMINACLIP:<fastaWithAdaptersEtc>:<seed mismatches>:<palindrome clip threshold>:<simple clip threshold>
#アダプターの切り方を決める。詳しくはマニュアルの６ページ目へ。
SLIDINGWINDOW:4:15 \
#ウインドウサイズと平均クオリティの閾値を設定
#4:15 であれば4bp ごとに平均クオリティ値を出し、クオリティ値が15を切ったら以降の配列を切る
CROP:指定数値 \ #クオリティ値に関係なく、リードの先頭からこの数だけ塩基を保持する
HEADCROP:指定数値 \ #クオリティ値に関係なく、リードの先頭からここで指定した数だけ塩基を除去
MINLEN:指定数値 \ #この長さ以下の配列は除去
わたしが使っていたオプションは以上ですが、Trimmomatic には他にもたくさんのオプションがあります。
LEADING: リードの先頭が指定したクオリティ値未満の場合は除去
TRAILING: リードの末端が指定したクオリティ値未満の場合は除去AVGQUAL: リードの平均クオリティ値が指定したクオリティ値未満の場合はリードを除去
これらは論文を読んでいてもよく見かけるので要検討です。
で、ここからが大量データを扱う際のコツなのですが、メインスクリプトを見てもらうとわかるように、このソフトウェアは解析サンプル名を一つ一つ指定する必要があります。RAD-seq などで数十～数百個体扱う場合、いちいち打ってられません。
そこでfor 文を使って、解析時に解析サンプル名を自動で書き換えていきます。
生データは大体名前に規則性があるので、それを利用します。
例えば今回
Q1-01C_S3_L001_R1_001.fastq.gz, Q1-02C_S4_L001_R1_001.fastq.gz,
Q1-03C_S5_L001_R1_001.fastq.gz... および
Q1-01C_S3_L001_R2_001.fastq.gz, Q1-02C_S4_L001_R2_001.fastq.gz,
Q1-03C_S5_L001_R2_001.fastq.gz...
というペアエンド生データの入ったsamples という名前のフォルダがあるとします。
for file in ./samples/*_R1_001.fastq.gz
#for file in ＝ファイルに対して処理をしてください
#* (ワイルドカード) は「なんでも入ります」の意味。今回は.1.fastq.gz という拡張子のすべてのサンプルについて処理してください、の意味。
dosmpl=${file%_S*_L001_R1_001.fastq.gz}
#サンプル番地をsmpl と定義
#「%」は「～より前」を意味するので、今回はfastq.gz ファイル名の中で_S より前（例：Q1-01C）を指定
smpl=${smpl#./sample/}
#「＃」は「~より後ろ」を意味するので、今回はファイル名の中で_Sより後（例：3_L001_R1_001.fastq.gz）は切り捨てられる
snum=${file%_L001_R1_001.fastq.gz}
#サンプル番号をsnum と定義し、snum は_L001_R1_001.fastq.gz の前を指定。
snum=${snum#*_S} #残った部分の内、「_S より後ろ」がsnum（例：Q1-03C_S5）
echo ${smpl} #確認のため、smpl の書き出し
echo ${snum} #確認のため、snum の書き出し
#ここまでで、サンプル名を自動的に取得する設定ができる
実際にやるときは自分のデータに合わせて書き換えてね
あとはこれを先のメインスクリプトに投入
trimmomatic PE -threads スレッド数 \
-phred33 \
-trimlog ログ名.txt \
./samples/ ${smpl}_S${snum}_L001_R1_001.fastq.gz \
./samples/ ${smpl}_S${snum}_L001_R2_001.fastq.gz \
./output_folder/${smpl}_pair_R1.fastq \
./output_folder/${smpl}_unpair_R1.fastq \
./output_folder/ ${smpl}_pair_R2.fastq \
./output_folder/${smpl}_unpair_R2.fastq \
ILLUMINACLIP:adapter.fasta:2:30:10 \
#【注意】フォルダ名の前の「./」はない方がいいかもしれないです！自分のPC でやった時は必要でしたが、サーバ上でやった時はこれがあると動きませんでした。
SLIDINGWINDOW:4:15 \
CROP:指定数値 \
HEADCROP:指定数値 \
MINLEN:指定数値 \
done
pigz -p output_folder ./output_folder/*.fastq
#最後にできたファイルを高速でガンジップ化
今日のところはこんな感じで。
そういえば、なぜこの方法だとアナコンダさんが必要なのかというと、Trimmomatic の元のjava ファイルだとfor 文による${}の置き換えが上手くいかないから。
スパコンでの解析の時にアナコンダさん経由せずにやろうとしたらいつものスクリプトで動かなくて泣きそうになったことがあったので。
ただ、スパコンは色々とイレギュラーなトラブルがあるので、もしかするとjava のせいではなくスパコンのせいかもしれないのですが…。
まあ、アナコンダさん入れといて損したことは今のところあんまりないのでよかったら使ってみてください。
おしまい。
