SNPcalling例1

#2【SNPcalling例1】TrimmomaticでQC

SNP callingはゲノムの変異を調べる解析です。

この記事はコースになっており、初回はこちらです。

ねこ研究員
ねこ研究員

今回の記事ではTrimmomaticを用いたリードのQCを行います!

リードのQCとはどういうステップですか?

ねこ研究員
ねこ研究員

SNP同定の前に、クオリティの低いリードを除去するステップです。

次世代シーケンサにより決められたリードの配列は必ずしも正しくありません。機械である以上、多少の間違いはあります。

Trimmomaticのステップでは、そういった間違いが多そうなリードや、また、短すぎて以降の解析の役に立たないリードを除去するなどの、fastqファイルのQuality Filteringをします。

実際のスクリプト

ねこ研究員
ねこ研究員

内容だけ知りたい方に、先に答えをお届けです。

まずは、除去前のリードの本数を数えて記録しておきましょう。

mkdir -p rep
a1=data/*R1*.gz
a2=rep/NumOfReads.txt
for i in `ls $a1`;do
    echo =================== >> $a2
    echo `basename $i` >> $a2
    expr `gzip -dc $i | wc -l` / 4 >> $a2
done
less rep/$2

続いて、Trimmomaticを走らせましょう。

少し時間がかかるのでqsubにしておきましょう。

#! /bin/bash
#$ -cwd
#$ -V
#$ -o out.txt
#$ -e error.txt
#$ -m abe

TARGET=${ARG1}
CPU=20
MINLEN=80

READ=data
DIR=analysis
Trimmomatic=${HOME}/bin/Trimmomatic-0.39/trimmomatic-0.39.jar
Illuminaclip=${HOME}/Trimmomatic-0.39/adapters/all.SE.fa
REF=genome/genome.fa

# Quality Filtering (Trimmomatic)
mkdir -p ${DIR}/fastq_clean
fastq_path=`ls ${READ}/*/${TARGET}*R1*.gz`
out_name=${DIR}/fastq_clean/`basename ${fastq_path} | cut -d_ -f1`
java -jar ${Trimmomatic} \
    PE -phred33 -threads ${CPU} \
    ${fastq_path} `echo ${fastq_path} | awk -F '_R1' '{print $1}'`_R2.fastq.gz \
    ${out_name}_R1.clean.fastq.gz ${out_name}_R1.unpaired.fastq.gz \
    ${out_name}_R2.clean.fastq.gz ${out_name}_R2.unpaired.fastq.gz \
    ILLUMINACLIP:${Illuminaclip}:2:30:10 \
    TRAILING:3 LEADING:3 SLIDINGWINDOW:4:15  MINLEN:${MINLEN}

 javaで動かします

ねこ研究員
ねこ研究員

java -jar [javaへのpath]でOK.
怖がらないで大丈夫!

javaはデフォルトで入っていることが多いと思います。

コマンドラインで java -version としてみましょう。

mac :~ user$ java -version
openjdk version "1.8.0_161"
OpenJDK Runtime Environment (build 1.8.0_161-b14)
OpenJDK 64-Bit Server VM (build 25.161-b14, mixed mode)
ねこ研究員
ねこ研究員

なお、OpenJavaよりOracleの方が早いかもです。

OracleのサイトからLinuxの*.tar.gzをDLしてダウンロードをしてパスを通しておきましょう。パスを通したら念のために再度バージョン確認を。もし変わっていなかったら、デフォのパスより先に参照するようにしましょう。

除去配列の指定

ねこ研究員
ねこ研究員

本コースの設定ではイルミナHiSeqのデータですから、HiSeqのアダプター除去をします。

 ILLUMINACLIP:${Illuminaclip}:2:30:10

上記部分が除去配列を指定する部分です。

イルミナのアダプター配列が記載されたMultifastaファイルを用意し、ファイルへのパスを指定しています。

なお、MultifastaファイルはTrimmomaticをダウンロードした際に一緒に付いてきます。もしなかったら「探す」or「作成」です。

その他の数値についてはTrimmomaticのマニュアルを参照してください。

ペアエンドの場合のTrimmomatic

ねこ研究員
ねこ研究員

シングルエンドとペアエンドで記述方法が異なるので注意です。

java -jar ${Trimmomatic} \ 
 PE -phred33 -threads ${CPU} \
 ${fastq_path} `echo ${fastq_path} | awk -F '_R1' '{print $1}'`_R2.fastq.gz \
 ${out_name}_R1.clean.fastq.gz ${out_name}_R1.unpaired.fastq.gz \
 ${out_name}_R2.clean.fastq.gz ${out_name}_R2.unpaired.fastq.gz \

まず、”PE”の記載が必要です。シングルエンドの場合は”SE”と書きます。

次に、インプットとするfastqのペアを隣り合わせに、上記コードのように指定しましょう。上記では、ペアのファイルパスを取得するために、片方のファイルパスをawkで一部削って工夫しています。

出力は4つあります。ペアエンドなので、出力は2つあります。これに加えて、除去された側のリードも出力として指定するために4つあります。

なお、ペアエンドの除去では、片方が除去対象となった場合、もう片方が除去対象でなくとも除去対象となります。

予めテストしよう

ねこ研究員
ねこ研究員

マニュアルでは36となっているけど実際どうすれば…?

 TRAILING:3 LEADING:3 SLIDINGWINDOW:4:15 MINLEN:${MINLEN}

上記コードが、除去の条件を定めている部分です。

基本的にはマニュアル内に記載されている例と同じで良いと思いますが、設定した値によってどれくらいの配列が除去されてしまうか、fastqファイルの最初の4万行ほどを抽出したファイルでテストをすると良いかもしれません。95%以上残っていれば良いでしょう。

ただし、今回はMINLENを80と厳しめにします。これは、信頼できるSNPを精度よく抽出したいからです。ちなみに、私がMINLEN:80でテストしたところ30%ほどのリードが落ちたことがあります。

結果をまとめよう

ねこ研究員
ねこ研究員

例えば下記のようなコードでQC後の残存リード数を記録しておきましょう!

mkdir -p rep
a1=analysis/fastq_clean/*R1*clean*.gz
a2=rep/NumOfReadsClean.txt
for i in `ls $a1`;do
    echo `basename $i` >> $a2
    expr `gzip -dc $i | wc -l` / 4 >> $a2
done

Trimmomaticのインストール方法

ねこ研究員
ねこ研究員

かなり簡単な部類です!  

配布サイトからバイナリファイル(つまり、*.jarのファイル)をダウンロードして、自分のローカルにおき、パスを通すだけです。

なお、パスを通さなくても、上記のコードのように直接指定しても別にOKです💡

 

ねこ研究員
ねこ研究員

今回は以上です。

次回はBWAによるマッピングを解説します。