次世代シーケンサーで得られた fastq データの精製方法を紹介します．

A. Trimmomatic: アダプター配列の除去
B. SolexaQA 2: 質の悪い配列を取り除く．
C. condetri: PCR による余分な重複を取り除く．

2018 年 11 月の時点で，A の解析を勧められました．B と C は古い技術のようです．つまり、A. Trimmomatic だけやれば良い、とアッセンブルに詳しい同僚から助言をもらいました。

Trimmomatic は，次世代シーケンサーで得られた fastq データから，アダプター配列を取り除くプログラムです．日本語の解説として，binoinformatics を参考にしました．manual も参考になります

Trimmomatic は，Illumina 最新式のシーケンサーから得られたデータの解析に不向きなようです．Trimmomatic は 4 色式で得られたデータ用のトリミングツールなので，Novaseq や NextSeq など，最新の 2 色式で算出されたデータの解析には適していないそうです．

以下のコマンドでヘルプが出力されます．

java -jar trimmomatic-0.35.jar

Slurm で解析する場合の job script は，以下です．前に記述したコマンドが先に実行されます．クオリティの低い部分か，あるいはアダプターのどちらを先に削除するのか，適宜調節する必要があります．以下は，アダプターだけでなく，クオリティの低い部分もトリムしています．

#!/bin/bash
#SBATCH --job-name=prefetch
#SBATCH --mail-user="jun.inoue@oist.jp"
#SBATCH --partition=compute
#SBATCH --mem=10G
#SBATCH --cpus-per-task=4
#SBATCH --ntasks=1 # 1 task

java -jar /work/SatohU/0inoue/software/Trimmomatic-0.38/trimmomatic-0.38.jar \
PE \
-threads 4 \
-phred33 \
-trimlog log.txt \
SRR5760179_sub1.fq \
SRR5760179_sub2.fq \
# .fq あるいは .fq.gz など圧縮ファイルでも OK。
paired_output_1.fq \
unpaired_output_1.fq \
paired_output_2.fq \
unpaired_output_2.fq \
# 2:30:10 は一致率．adapters は以下を参照．
ILLUMINACLIP:/address/adapters.fa:2:30:10 \
# LEADING と TRAILING は先端と末端のクオリティを示す．
# クオリティ 20 以下になったら，それよりもアダプター側を削除．
LEADING:20 \
TRAILING:20 \
# 4 塩基づつで平均 score 15 以下を切断．
SLIDINGWINDOW:4:15 \
# 36 塩基以下は捨てる．
MINLEN:36

Index 1 (i7) Adapters
CAAGCAGAAGACGGCATACGAGAT[i7]GTCTCGTGGGCTCGGAGATGTGTATAAGAGACA

の場合，i7 はバーコードシーケンスであり，これを含めた全体で試薬固有の配列になっているそうです．adapters.fa の記述は以下のような感じです：

[cluster:adapters]$ cat TruSeq_LT+HT.consensus.fa
>Read 1
AGATCGGAAGAGCACACGTCTGAACTCCAGTCA
>Read 2
AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT

解析が終了したら，fastqc でアウトプットファイル，paired_output_1.fq と paired_output_2.fq をチェックしてみましょう．問題がなければ，Trinity で解析を行います．

結果の解釈
binoinformatics に丁寧な解説があります．

Adapter Content
グラフの最も右側の数 bp が 5% 以下ぐらいであれば，アダプター配列は除去されたと考えて良いようです．

Per base sequence content
ATGC の量が異なるのは，mRNA の場合，片側しか読んでいないためです．グラフ左端方がバラバラなのは，シグナルが一定していないためと考えられ，これは毎回生じるそうです．この状態が良いのか悪いのかはわかりませんが，気になるようであれば，削除して使っても良いようです．

エラー対策
slurm　というスケジューラーで解析した場合に，slurmstepd: error: slurm_receive_msg: というエラーが出ることがあります．

TrimmomaticPE: Completed successfully
slurmstepd: error: slurm_receive_msg: Socket timed out on send/recv operation

しかし，TrimmomaticPE: Completed successfully が出た後なので，結果は問題なく使えるようです．

SolexaQA 2 は次世代シーケンサから得られた生の配列データから，質の悪い配列を (大量のコピー配列も含めて) 取り除くプログラムです．どの解析もそれほど時間はかかりません．
　以前 perl スクリプトに分割されていた機能が現在は SolexaQA++ というプログラムにまとめられています．詳しくはこちらをご覧ください．

FASTQ 形式のファイルを読み込みます． まず配布されている Example を解析してみると良いでしょう．一連の解析ではたくさんのアウトファイルができるので，解析ごとに不要なものは異なるディレクトリに移動した方が良いです．

ソフトウェア: こちらから src をクリックし，SolexaQA++_v3.1.7.1.zip　など最新のファイルをダウンロードしてください．解凍して得られたディレクトリ内部に，SolexaQA.pl , DynamicTrim.pl, LengthSort.pl が入っています．

examples: 上と同じサイトから examples ディレクトリに移動すればダウンロードできます．GAIIx と MiSeq について，それぞれ良い例と悪い例がダウンロードできます．

SolexaQA.pl input_files [-p|probcutoff 0.05] [-h|phredcutoff 13] [-v|variance] [-m|minmax] [-s|sample 10000] [-b|bwa] [-d|directory path] [-sanger -solexa -illumina]

現在は SolexaQA++ というプログラムに以前 perl スクリプトに分割されていた機能がまとめられています．詳しくはこちらをご覧ください．

DynamicTrim
クオリティーの悪いリードを削ります．解析は forward, reverse の片方ずつ行います．FAQ の「How does DynamicTrim deal with paired-end data?」によると，あるリードがすべて削除されることはないので，順番がずれることはないようです．
　アウトファイルでは，.trimmed という拡張子のファイルに選定されたリードが保存されています．SRR5760179_1.fastq (ペアードエンド，24G) の解析で，1 時間ぐらいで終わりました．

SolexaQA++ dynamictrim good_MiSeq_dataset.fq -h 20
(以前は，
perl DynamicTrim.pl good_MiSeq_dataset.fq -h 20
でした)

DynamicTrim.pl input_files [-p|probcutoff 0.05] [-h|phredcutoff 13] [-b|bwa] [-d|directory path] [-sanger -solexa -illumina] [-454]

-p|probcutoff
probability value (between 0 and 1) at which base-calling error is considered too high (default; p = 0.05) *or*

-h|phredcutoff
Phred score (between 0 and 40) at which base-calling error is considered too high
20% (-h 20) だと，1% のエラー確率で，ただ悪いのを削る，30% だと 0.1% クオリティーなので良いのを残す，というイメージです．

-b|bwa
use BWA trimming algorithm

-d|directory
path to directory where output files are saved

-sanger
Sanger format (bypasses automatic format detection)

-solexa
Solexa format (bypasses automatic format detection)

-illumina
Illumina format (bypasses automatic format detection)

-454
set flag if trimming Roche 454 data (experimental feature)

アウトファイル
good_MiSeq_dataset.fq.trimmed

実際に得られた結果です．選定された質の高いリードが保存されています．以後の解析に用いるファイルです．

good_MiSeq_dataset.fq.trimmed_segments

クオリティを示す数字が出力されているようです．以降は使わないファイルです．

good_MiSeq_dataset.fq.trimmed_segments.hist.pdf

上のファイルをグラフ化したファイルです．R を使って作成されます．以降は使わないファイルです．

temp.R

R のバッチファイルです．R が使えない環境だと，このファイルがが残されたままになります．以降は使わないファイルです．

LengthSort
短いリードを削ります．

####
1. Paired end ファイルの解析．
Forward と Revers (以下では leads.left.fq と leads.right.fq) がある場合は以下のコマンドです．例題として，Trinity で配布されている reads.left.fq と reads.right.fq を用いています (trinityrnaseq_r2013-02-25/sample_data/test_Trinity_Assembly)．解析はすぐに終わります．

SolexaQA++ lengthsort reads.left.fq reads.right.fq -l 50
(以前は，
perl LengthSort.pl reads.left.fq reads.right.fq -length 50
でした)

LengthSort.pl one single-end or two paired-end FASTQ files [-l|length 25] [-d|directory path]

-l|length
length cutoff [defaults to 25 nucleotides]
私は通常 50bp (-length 50) で解析を行っています．

-d|directory
path to directory where output files are saved

 

アウトファイル

reads.left.fq.paired1

実際の結果です．length cut off と同じかこれより長い Forward (reads.left.fq) リードが保存されています．

reads.left.fq.paired2

実際の結果です．length cut off と同じかこれより長い Reverse (reads.right.fq) リードが保存されています．ファイル名は reads.left.fq.paired2 となっていますが，実際には reads.right.fq の解析結果が保存されています．

reads.left.fq.discard

length cutoff よりも短いため取り除かれたリードです．

reads.left.fq.single

単一 fq ファイルを指定した場合のアウトファイルです．Paired end (reverse と forward の対になったファイル) でを解析したら空になりました．

reads.left.fq.summary.txt

解析結果の要約です．

その他
2 つのファイルがペアになっていない (数が合わない) と，エラーメッセージが出ます．

[i2:examples]$ perl LengthSort.pl reads.left1.fq reads.right.fq -length 50
error: files reads.left1.fq and reads.right.fq do not seem to be paired

####
2. Single ファイルの解析．
データが片方だけ (Forward と Revers ではない) の場合です．

perl LengthSort.pl good_MiSeq_dataset.fq.trimmed -length 50

アウトファイル
アウトファイルでは，.paired1 と .paired2 という拡張子のファイルに選定されたリードが保存されています．

good_MiSeq_dataset.fq.trimmed.single

実際の結果です．length cut off と同じかこれより長いリードが保存されています．

good_MiSeq_dataset.fq.trimmed.discard

length cutoff よりも短いため取り除かれたリードです．以降は使わないファイルです．

good_MiSeq_dataset.fq.trimmed.summary.txt

上二つのファイルに含まれていたリードの数です．実際の数はあってないようなので，何か他の処理も行っているようです．以降は使わないファイルです．

filterPCRdupl.pl は condetrifilter パッケージの一部です．PCR によって増幅された余分なコピー (PCR duplicate) を取り除く Perl スクリプトです． k-mer アッセンブルを採用しているそうです．ここでは filterPCRdupl_v1.01.pl を使った解析を説明します．ただし，reference 配列がないような生物の場合は，condetri による処理を行い方が良い，という意見もあります．
＊Paired end のファイルしか対応していないようです．

例題として，Trinity で配布されている reads.left.fq と reads.right.fq を用いています (trinityrnaseq_r2013-02-25/sample_data/test_Trinity_Assembly)．

perl filterPCRdupl_v1.01.pl -fastq1=reads.left.fq -fastq2=reads.right.fq -cmp 50

# filterPCRdupl.pl
# November 2010
# Author: Linnéa Smeds (linnea.smeds@ebc.uu.se)
# --------------------------------------------------------------
Description: Extract all non-redundant read pairs in a fastq file pair, by comparing the first N bases in both reads between all pairs. If there are several copies of the same pair, only the pair with the highest quality score is kept. Also saves a histogram of the copy distribution.
Usage: perl filterPCRdupl.pl -fastq1=file1 -fastq2=file2 [-prefix=s -cmp=N]

-fastq1=file -fastq2=file
Fastq files with the read pairs, (pair1 in first file and pair2
in the second). Reads must have the same order in both files.

-prefix=string
Prefix for the output files. The filtered fastq files will be
named prefix_uniq1.fastq and prefix_uniq2.fastq, and the
histogram with copy distribution will be named prefix_copy.hist.
[same prefix as file1].

-cmp=number
The number of bases from the start of the reads in each pair
that will be used for comparison [50].
ここでは -cmp 50 として，50bp の領域の類似度を見ています．

-h
Print this help message.

アウトファイル
reads.fq_uniq1.fastq

reads.left.fq から選定されたリードが保存されています．

reads.fq_uniq2.fastq

reads.right.fq から選定されたリードが保存されています．

reads.fq_copy.hist

解析結果の要約です．