SNPcalling例1バイオインフォ

#4【SNPcalling例1】Picardでマーキング

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

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

ねこ研究員
ねこ研究員

今回の記事ではPicardによるマーキングを行います!

Picardによるマーキングとはどういうステップですか?

ねこ研究員
ねこ研究員

後に続くGATKというソフトウェアの入力を作成するステップです。

具体的には、下記の2点をbamファイルのタグ領域に加えることを行なっています。

  • RGタグを加える
  • 重複したリードに「重複してますよ」タグを加える

実際のスクリプト

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

TARGET=${ARG1}
MEM=${ARG2}
DIR=analysis
PICA=~/src/picard/build/libs/picard.jar
REF=genome/genome.fa

# Note that the latest Oracle java may not be compatible with Picard.
# If that's the case, recommend to use older one.

# Add RG TAG
mkdir -p ${DIR}/pica1
java -Dpicard.useLegacyParser=false -Xmx${MEM} -jar ${PICA} \
    AddOrReplaceReadGroups \
    -I ${DIR}/bam/${TARGET}.refsort.bam \
    -O ${DIR}/pica1/${TARGET}.rg_added_sorted.bam \
    -SO coordinate \
    -RGID 1 \
    -RGLB Project \
    -RGPL illumina \
    -RGPU run \
    -RGSM ${TARGET}_all
# Remove(Mark) dupulicates
mkdir -p ${DIR}/pica2
java -Dpicard.useLegacyParser=false -Xmx${MEM} -jar ${PICA} \
    MarkDuplicates \
    -I ${DIR}/pica1/${TARGET}.rg_added_sorted.bam \
    -O ${DIR}/pica2/${TARGET}.dedupped.bam \
    -CREATE_INDEX true \
    -VALIDATION_STRINGENCY SILENT \
    -M ${DIR}/pica2/${TARGET}.output.metrics
#   -REMOVE_DUPLICATES true
#rm ${DIR}/pica1/${TARGET}.rg_added_sorted.bam

Javaで動作しています

ねこ研究員
ねこ研究員

下記1行の通り、Picardの動作にはJavaが必要です。ない方はインストールしましょう。

java -Dpicard.useLegacyParser=false -Xmx${MEM} -jar ${PICA} \

ただ、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してダウンロードをしてパスを通しておきましょう。パスを通したら念のために再度バージョン確認を。もし変わっていなかったら、デフォのパスより先に参照するようにしましょう。

最新版のJavaとPicardは相性が悪いかも ?

ねこ研究員
ねこ研究員

最新版のOracleでPicardを走らせたところ、バグりました。

バグに遭遇した場合は、デフォルトでインストールされているJavaや古いバージョン、または、Open Javaなど別のものを検討してみてください。

Javaのパラメタ部分の説明

-Dpicard?!#$”a………って何??

-Dpicard.useLegacyParser=false -Xmx${MEM}

まずは、-Dpicard.useLegacyParser=false について。

2019年6月現在、Picardは大きなバージョン変更中です。このバージョン変更は、

  • 旧バージョン:I=xxxx
  • 新バージョン:-I xxxx

のようなsyntaxの変更を伴うものです。現在、どちらのバージョンで記述してもプログラムは走りますが、旧バージョンの文法で書くと警告が出力されます。一方、新バージョンで書いても、新バージョンで書いている、という旨の余計な出力が出ます。そこで、その余計な出力をさせないために -Dpicard.useLegacyParser=false を加えています。

ねこ研究員
ねこ研究員

-Xmxは??

これは使用するメモリの指定です。例えば、-Xmx32gなどと書けばOKです。従って、変数MEMには32gなどを格納します。

16gしかない環境で32gなどとするとバグ(って停止してくれ)るから注意してください。

RGタグとは?

ねこ研究員
ねこ研究員

下記部分は、RGタグを加えています。このタグがないbamファイルは、次のステップで使うGATKの入力にならずエラーとなります。

# Add RG TAG
mkdir -p ${DIR}/pica1
java -Dpicard.useLegacyParser=false -Xmx${MEM} -jar ${PICA} \
    AddOrReplaceReadGroups \
    -I ${DIR}/bam/${TARGET}.refsort.bam \
    -O ${DIR}/pica1/${TARGET}.rg_added_sorted.bam \
    -SO coordinate \
    -RGID 1 \
    -RGLB MEDAKA \
    -RGPL illumina \
    -RGPU run \
    -RGSM ${TARGET}_all

まず、上記の操作前後のファイルを比較してみましょう。

入力はBWAの出力のbamでした。最初の1行を覗いて見ると…

samtools view xxx.refsort.bam | less
E00441:179:HHNNFCCXY:3:1103:30086:10433 0       1       56      0       151M    *       0       0       TGTTTTGGTTTTATTTAGGGTCAAACTTGTGCATTCTTAGAGGCGTGTCTTCTTGATTGACAGGTGGGCGGACCCATGGGTCAACTTAATTCGCTTTGACCCACGCCATTATTACTATTCTTAGAACGCAGAATGTTCTCGCCTAGTCACT AAA<FFJJFFJJJJJFJJJJJJJJJFJJJJJJJJJJJJJFJJJJJJJJJJJJJJJJJJJJJJJJJFJFJJJJJJJJJJJJFJJJJJJJJJJJJJJJJJFJJ<AJJJAAFFJAJJJJJAJFJJJJJJJJJJJJJFJJJJ<JFJJJAFFJJJJ NM:i:1  MD:Z:114T36     AS:i:146        XS:i:146        XA:Z:14,+2947330,151M,1;14,+2954038,151M,1;

一方、上記コードの出力である、xxx.rg_added_sorted.bam の最初の行を覗いて見ると…

E00441:179:HHNNFCCXY:3:1103:30086:10433 0       1       56      0       151M    *       0       0       TGTTTTGGTTTTATTTAGGGTCAAACTTGTGCATTCTTAGAGGCGTGTCTTCTTGATTGACAGGTGGGCGGACCCATGGGTCAACTTAATTCGCTTTGACCCACGCCATTATTACTATTCTTAGAACGCAGAATGTTCTCGCCTAGTCACT AAA<FFJJFFJJJJJFJJJJJJJJJFJJJJJJJJJJJJJFJJJJJJJJJJJJJJJJJJJJJJJJJFJFJJJJJJJJJJJJFJJJJJJJJJJJJJJJJJFJJ<AJJJAAFFJAJJJJJAJFJJJJJJJJJJJJJFJJJJ<JFJJJAFFJJJJ XA:Z:14,+2947330,151M,1;14,+2954038,151M,1;     MD:Z:114T36     RG:Z:1  NM:i:1  AS:i:146        XS:i:146

となっています。

最後の方に RG:Z:1 というタグが付いています。RGタグをつける操作とは、たったこれだけのことを行うものです。

RGタグやその他のオプションは何を意味するの?

ねこ研究員
ねこ研究員

基本的にサンプルに関する情報を与えるだけです。

詳しくは、Picardのマニュアルを参照した方が正確でしょう。

duplicateとは?

ねこ研究員
ねこ研究員

いきなりduplicatesと言われても…

Duplicatesは例えば「PCRのバイアスで生じた重複したリード」などのことです。

これにより「本来の数より多くカウントすることを防ぐ」をします。

というのも、SNPを判断するときには「何本のリードで変異が生じていたか?」で判断するからです。

下記が、Duplicatesをマークするコードです。

# Remove(Mark) dupulicates
mkdir -p ${DIR}/pica2
java -Dpicard.useLegacyParser=false -Xmx${MEM} -jar ${PICA} \
    MarkDuplicates \
    -I ${DIR}/pica1/${TARGET}.rg_added_sorted.bam \
    -O ${DIR}/pica2/${TARGET}.dedupped.bam \
    -CREATE_INDEX true \
    -VALIDATION_STRINGENCY SILENT \
    -M ${DIR}/pica2/${TARGET}.output.metrics
#   -REMOVE_DUPLICATES true
#rm ${DIR}/pica1/${TARGET}.rg_added_sorted.bam

上記コードにより出力された xxxx.dedupped.bamの例は下記です。

E00441:179:HHNNFCCXY:3:1114:25611:31547	1024	1	238	0	151M	ACAAGTGAACATTCTGTAAATTTATTGAAGATAGATTTTCTGAAAAGAAGTCATTGTGGTTCTGTGAGTTGTGTTCAGTTTTTTCTTGTATGAAGTAAGGTTTGAGTCAGGTTTGTGGCAGAGTCTGAATCTCCTCCAATTTCTATTTCTT	AAFFFJJJJJJJAJJJJJJJJJJJJJJJFFJJJJJJFJJJJJJJJJFFJJ7FJF<JJFJJJJJJJJFJAJJJJFJFFJJJJJJJJJJJJJJJJJJJJJJFFFJJFJJJJJJJJJJJJJJFFJJ7FJJFFFFJJF7F<F-AJAJJ<FFAJJJ	XA:Z:14,+2947512,151M,0;14,+2954220,84M2D27M1D10M1D30M,5;	MD:Z:151	PG:Z:MarkDuplicates	RG:Z:1	NM:i:0	AS:i:151	XS:i:151
E00441:179:HHNNFCCXY:3:1103:30086:10433 0 1 56 0 151M * 0 0 TGTTTTGGTTTTATTTAGGGTCAAACTTGTGCATTCTTAGAGGCGTGTCTTCTTGATTGACAGGTGGGCGGACCCATGGGTCAACTTAATTCGCTTTGACCCACGCCATTATTACTATTCTTAGAACGCAGAATGTTCTCGCCTAGTCACT AAA<FFJJFFJJJJJFJJJJJJJJJFJJJJJJJJJJJJJFJJJJJJJJJJJJJJJJJJJJJJJJJFJFJJJJJJJJJJJJFJJJJJJJJJJJJJJJJJFJJ<AJJJAAFFJAJJJJJAJFJJJJJJJJJJJJJFJJJJ<JFJJJAFFJJJJ XA:Z:14,+2947330,151M,1;14,+2954038,151M,1; MD:Z:114T36 PG:Z:MarkDuplicates RG:Z:1 NM:i:1 AS:i:146 XS:i:146

どちらの行でも、PG:Z:MarkDuplicates のタグがついています。これは、Duplicates処理をしましたという意味です。上記2行の違いは、2カラム目のFLAG部分です。Duplicatesと判断された場合、ここに1024が入ります。一方、そうでない場合は、1024ではない数値が入ります。この解説についてはマニュアルに下記のように解説されています。

The tool’s main output is a new SAM or BAM file, in which duplicates have been identified in the SAM flags field for each read. Duplicates are marked with the hexadecimal value of 0x0400, which corresponds to a decimal value of 1024. If you are not familiar with this type of annotation, please see the following blog post for additional information.

ねこ研究員
ねこ研究員

Duplicatesは除かなくていいの?!

後のプロセスで除くことにして、ここではマーキングだけすることが多いです。ただ、上記コードのコメントアウトのようにして、この時点で除去することもできます。

ねこ研究員
ねこ研究員

その他のオプションについては?

こちらもPicardのマニュアルを参照してみてください。

Picardのインストール方法は?

ねこ研究員
ねこ研究員

基本的にJavaで動くスクリプトはインストールが簡単です。

PicardのWebサイトの右上からダウンロードして、展開します。もしくはPicard開発の大元 Broad InstituteのGitHubのPicardのリポジトリからダウンロードしましょう。

その中にあるjarファイルは既に使える状態です。ですから、picard.jar の path を上記のコード部分に差し替えてあげればすぐに実施が可能です。

 

ねこ研究員
ねこ研究員

今回はここまでです!次回はGATKその1です💡