2. シークエンスデータのQC・トリミング

2021.06.10 初稿

シークエンスデータにはアダプター配列が付加されており、また低クオリティの塩基も含まれているため、解析を始める前にまずデータのクオリティチェック（QC）とフィルター処理を行う必要があります。シーケンスデータのQC、フィルタリング、及びトリミングを行うツールとしてFastQCやTrimmomaticが有名ですが、これらは各ステップで組み合わせて使う必要があるため、時間がかかり非効率になるという問題点がありました。その問題点を解消し、一つのツールで全てをこなしてくれるfastpを私は使用しました。

特徴 (Githubより引用 (https://github.com/OpenGene/fastp))

0. comprehensive quality profiling for both before and after filtering data (quality curves, base contents, KMER, Q20/Q30, GC Ratio, duplication, adapter contents...)

filter out bad reads (too low quality, too short, or too many N...)
cut low quality bases for per read in its 5' and 3' by evaluating the mean quality from a sliding window (like Trimmomatic but faster).
trim all reads in front and tail
cut adapters. Adapter sequences can be automatically detected,which means you don't have to input the adapter sequences to trim them.
correct mismatched base pairs in overlapped regions of paired end reads, if one base is with high quality while the other is with ultra low quality
trim polyG in 3' ends, which is commonly seen in NovaSeq/NextSeq data. Trim polyX in 3' ends to remove unwanted polyX tailing (i.e. polyA tailing for mRNA-Seq data)
preprocess unique molecular identifer (UMI) enabled data, shift UMI to sequence name.
report JSON format result for further interpreting.
visualize quality control and filtering results on a single HTML page (like FASTQC but faster and more informative).
split the output to multiple files (0001.R1.gz, 0002.R1.gz...) to support parallel processing. Two modes can be used, limiting the total split file number, or limitting the lines of each split file.
support long reads (data from PacBio / Nanopore devices).

以下がコマンドになります。

#-q 30: クオリティスコア30, -u 40: 条件を満さない塩基が40%以上ならば除去

#paired endのデータを使用するので、-Iオプションで2つ目のファイルを追加

$ fastp -i sample_1_1.fastq -I sample_1_2.fastq -o qc_sample_1_1.fq -O qc_sample_1_2.fq -h sample_1_report.html -q 30 -u 40

解凍したfastqファイルだけでなくfastq.gzファイルを使うことができるので、ファイル要領を増やさないためにもfastq.gzで解析を進めるのがおすすめです。

以下、paired endでのQC、トリミングをfor loopで回すコマンド。

$ ls ./fastq/*_1*|sed -e ’s/fastq\/$.*$_1.*/\1/g' > FastQ_list.txt

$　for file in `cat FastQ_list.txt`

fastp --detect_adapter_for_pe -i ./fastq/${file}_1.fastq.gz -I ./fasatq/${file}_2.fastq.gz -o ./qc/${file}_1.fastq.gz -O ./qc/${file}_2.fastq.gz -h ./qc/${file}_report.html -j ./qc/${file}_report.json

done

これをシェルスクリプトにして、qsubでスパコンにジョブとして投げる。

以下がスクリプトにした時の書き方

-------------------

"qc_filename.sh"

#!/bin/bash

#$ -S /bin/sh

#$ -cwd

#$ -o ~/results_sh_eando

#$ -e ~/results_sh_eando

echo starting at

date

for file in `cat fq_gz_list.txt`

singularity exec -e /usr/local/biotools/f/fastp:0.22.0--h2e03b76_0 \

fastp --detect_adapter_for_pe \

-i ./${file}_1.fastq.gz \

-I ./${file}_2.fastq.gz \

-o ./QC_out/qc_${file}_1.fastq.gz \

-O ./QC_out/qc_${file}_2.fastq.gz \

-h ./QC_out/${file}_report.html \

-j ./QC_out/${file}_report.json \

-w 16

done

echo ending at

date

-------------------

QC_outに造られたレポートファイルで、シークエンスのクオリティとトリミングの結果を確認する。

htmlファイルは、ウェブブラウザーで開くことができる。

参考にしたウェブページ

・fastpについて：https://github.com/OpenGene/fastp#fastp

・使い方やQC結果の見方：https://kazumaxneo.hatenablog.com/entry/2018/05/21/111947