クオリティの低いリードを取り除く方法は、大きく分けて２パターンあります。1つ目は、リード全体をみてクオリティが低い場合は、リードごとすべて取り除く方法です。２つ目は、リードの3'末端のQualityをみて、クオリティが高い塩基まで削除する方法です。
ゲノムアセンブルする場合は、できるだけリードは長いほうが有利になります。そのため、前者の方法を採用した方が良いです。逆にリードの長さにあまり影響を受けない解析の場合は、後者の方法を選択しても良いでしょう。
１．リードの全体のクオリティを考慮する
FASTX-toolkit
ここでダウンロードを行う。リードの 80% 以上の塩基がクオリティ 20 未満のリードならば除去する場合
| fastq_quality_filter -q 20 -p 80 -i read_R1.fastq -o read_R1_HQ.fastq |
R
同様のことをRでやろうとするならば、以下の通りです。
in_f <- "SRR037439.fastq" out_f <- "hoge1.txt" param1 <- 20
param2 <- 0.2
library(ShortRead) fastq <- readFastq(in_f)
hoge <- as(quality(fastq), "matrix")
obj <- (rowSums(hoge < param1) <= width(fastq)*param2)fastq <- fastq[obj]
writeFastq(fastq, out_f) |
参照：http://www.iu.a.u-tokyo.ac.jp/~kadota/r_seq.htm
２．リードの末端のクオリティを考慮する
FASTX-toolkit
ここでダウンロードを行う。3' 末端からクオリティが 20 未満の塩基をトリミングし、長さが 50 bp 未満となったリードを破棄する場合
| fastq_quality_trimmer -t 20 -l 50 -i read_R1.fastq -o read_R1_HQ.fastq |
Sickele
ここでダウンロードする。Paired-end readの処理を行うには、sickeleが便利である。
クオリティが 20 未満の塩基をトリミングし、長さが 50 bp 以上のPaired-end readを出力する場合
| sickle pe -f input_file1.fastq -r input_file2.fastq -t sanger -o trimmed_output_file1.fastq
-p trimmed_output_file2.fastq -s trimmed_singles_file.fastq -q 20 -l 50 |
オプション
-f: read1 input
-r: read2 input
-o: read1 output
-p: read2 output
-s: 半端モンなread
-q: quality値指定
-l: 長さの閾値