fastpを用いたアダプタートリミングとクオリティコントロール
fastp
fastpは、アダプタートリミングとクオリティコントロールを同時に行う、高速なNGSデータ前処理パイプラインである。 fastqファイルに対応しており、fq.gzでもランが可能。 詳しくはfastpのgithubを参照
インストール
condaによるインストールがシンプルかつ簡単。 git cloneによるダウンロードとコンパイルを用いたバイナリのインストールも可能(著者は試していない)。
#condaによるインストール。著者はpython3.9で走らせてるが、割とどのバージョンでも動く。 conda install -c bioconda -y fastp conda update -c bioconda -y fastp #githubよりソースコードをダウンロードしてコンパイル(pathを通す必要がある) git clone https://github.com/OpenGene/fastp.git cd fastp make sudo make install #インストールされたか確認。インストールされているとバージョンが出力される。 fastp -v #fastp 0.22.0
fastpの実行
ここでは、基本的なオプションを説明する。
fastpはアダプター配列を自動検出する。自動検出はオプション-A
をつけることで無効化できる。
※自動検出は内部プログラムでなにが実行されているかわかりずらいので、嫌う人もいる。
著者はアダプター配列がわかれば-A
をつけて配列指定、わからなければ自動検出で実行している。
#githubよりソースコードをダウンロード。\はコードの改行(見やすくなる)。 fastp \ -i DNBSEQ_2021-A1_Read1.fq.gz \ #インプット1 -I DNBSEQ_2021-A1_Read2.fq.gz \ #インプット2 -o A1_Read1_tremmed.fastq.gz \ #アウトプット1 -O A1_Read2_tremmed.fastq.gz \ #アウトプット2 -h report.html \ #トリミング後のクオリティ値。htmlで開ける。 -j report.json \ #トリミング後のクオリティ値。json方式。 -q 20 \ #phred quality scoreの平均が20以下(エラー率1%以上)の領域を除去(低クオリティを除去) -t 1 \ #3'末端の1塩基を強制トリミング(フォワード) -T 1 \ #3'末端の1塩基を強制トリミング(リバース) -l 20 \ #20bp以下の配列を強制除去 -w 8 \ #スレッド数 --adapter_sequence AAGTCGGAGGCCAAGCGGTCTTAGGAAGACAA \ #トリミングするアダプターの指定 --adapter_sequence_r2 AAGTCGGATCGTAGCCATGTCGTTCTGTGAGCCAAGGAGTTG \ #トリミングするアダプターの指定 >& fastp.log & #実行後にログを作成
以下はfor文による自動化。フォルダ内のサンプルをすべて処理する。 使用するシーケンサーによってraw dataの出力形式が異なるので、その部分については適宜改変が必要。 まずサンプルリストを作成し、そのリスト内のサンプルをすべて処理するスクリプト。
#サンプルリストの作成。_R1_001.fastq.gzの部分は使用するサンプル名の後ろ側。自身のサンプルに合わせて改変。 #以下はMiseq出力ファイル形式。Miseq出力のraw dataであれば、コピペでいける。 ls *_L001_R1_001.fastq.gz | sed -e 's/_L001_R1_001.fastq.gz//g' > list.txt #listの中身を確認 less list.txt #トリミング&フィルタリング for file in `cat list.txt`; do fastp -i ${file}_L001_R1_001.fastq.gz -I ${file}_L001_R2_001.fastq.gz -o ${file}_R1_tremmed.fastq.gz -O ${file}_R2_tremmed.fastq.gz -h ${file}.html -j ${file}.json -q 20 -t 1 -T 1 -l 20 -w 8 done
参考(2023/04/05)
cutadaptを用いたアダプター配列のトリミング
cutadapt
cutadapt はpythonで書かれたハイスループットシーケンシングリードからアダプター配列、プライマー、ポリA テールおよびその他のタイプの不要なシーケンスを検出して削除するプログラムである。 詳しくはこちらを参照
インストール
condaによるインストールが簡単。 ホームページでは仮想環境の立ち上げを推奨しているが、他のソフトとはあまり干渉しない(著者はまだ経験していない)。
#condaによるインストール。著者はpython3.9で走らせてるが、割とどのバージョンでも動く。 conda install -c bioconda cutadapt #pipによるインストール pip install cutadapt #インストールされたか確認。インストールされているとバージョンが出力される。 cutadapt --version
cutadaptの実行
cutadaptは指定された配列を除去する。 以下で基本的な方法を記載するが、プライマー除去の場合、3'末端側は配列を相補佐に反転させなければならない。 ※アダプター配列がわからない場合は、自動検出して除去できるfastpを推奨する。 fastpを用いたアダプタートリミングとクオリティコントロール
#3' 末端からアダプター配列を除去する。 cutadapt -g AGGATTAGATACCCTGGTA input.fastq > output.trimmed.fastq #5' 末端からアダプター配列を除去する。 cutadapt -g AGGATTAGATACCCTGGTA input.fastq > output.trimmed.fastq #3' 末端と 5' 末端からアダプター配列を除去す cutadapt -b AGGATTAGATACCCTGGTA input.fastq > output.trimmed.fastq #リードとアダプター配列の許容するオーバーラップ塩基数を指定する。デフォルトは 3 。オーバーラップを変更したい場合は -O オプションを利用する。 #-e 0.2をつけることでミスマッチ20%を許容する。 cutadapt -O 5 -e 0.2 -b AGGATTAGATACCCTGGTA input.fastq > output.trimmed.fastq
以下はフォワード、リバースの両方のプライマーを除去するプログラム。 まず、変数にプライマーを代入し、相補佐側の配列を作成する。
FWD='AGGATTAGATACCCTGGTA' FWDRC=$(echo "${FWD}" | tr ACGTMRYKVHDBacgtmrykvhdb TGCAKYRMBDHVtgcakyrmbdhv | rev) REV='CRRCACGAGCTGACGAC' REVRC=$(echo "${REV}" | tr ACGTMRYKVHDBacgtmrykvhdb TGCAKYRMBDHVtgcakyrmbdhv | rev
cutadaptで代入したプライマーを除去
cutadapt \ -j 4 \ #スレッド数 -a ${FWDRC} \ #R1の3’末端から除く配列 -A ${REVRC} \ #R3の3’末端から除く配列 -g ${FWD} \ #R1の5’末端から除く配列 -G ${REV} \ #R2の5’末端から除く配列 -n 10 \ #繰り返し除去する回数 --match-read-wildcards \ #プライマー配列の縮重コードを翻訳 --minimum-length 20 \ #トリム後20bp以下の配列を除去 --discard-untrimmed \ #トリムが行われなかった配列を除去 -o out.1.fastq \ #アウトプット1 -p out.2.fastq \ #アウトプット2 read1.fastq.gz \ #インプット1 read2.fastq.gz #インプット2