gatk4のRNA-seq bestpractice
TL;DR
genomic sequenceと少しRNA-seqのパイプラインは違うので、メモを。基本的には一緒なのですが、
- MappingにSTARのtwo pass Mappingを使用する
- SplitNCigarReadsを適用する
- VQSR処理をしない
の3点が異なっています。これらの処理は基本的にgatk4-rnaseq-germline-snps-indelsを参考にしています。順番に見ていきます。
全体の流れ
全体の流れは以下の通りです。
- 必要なファイルの準備
- Genomeへのmapping
- PCR duplicatesの除去
- SplitNCigarReads
- BQSRの計算、適用
- 変異の検出
- GVCFからVCFへの変換
- 変異のフィルター
環境
GATKはDockerイメージを公式が配布しているので、それを使います[参考]。
docker pull broadinstitute/gatk:4.1.3.0
docker run --rm -it broadinstitute/gatk:4.1.3.0
1. 必要なファイルの準備
とりあえずhg38を例に説明します。
必要なファイルは以下の通りです
- fasta
- fastaのインデックス(fasta.fai)
- fastaのdict(.dict)
- gtf/gffファイル(なくてもよい)
- STARのindex
- known SNPs
基本的にはGATKを提供している基本的にBroadInstituteのGoogle Cloud Platformからダウンロードすればいいと思います(参考)。STARのindexだけダウンロードできないので自力で作ります。
ということで準備しています。
必要なファイルのダウンロード
- Homo_sapiens_assembly38.fasta
- Homo_sapiens_assembly38.dict
- Homo_sapiens_assembly38.fasta.fai
- Homo_sapiens_assembly38.dbsnp138.vcf.gz
- Homo_sapiens_assembly38.known_indels.vcf.gz
- Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
それとvcf.gzに対応するidxファイルもダウンロードしておきましょう。
また、GTF/GFFファイルはUCSCのTableBrowserあたりからダウンロードしましょう。
faiの作成方法
samtools faidx genome.fa
dictファイルの作成方法
picard CreateSequenceDictionary R=genome.fa O=genome.dict
2. STARを使ったMapping
Indexの準備
STAR \
--runMode genomeGenerate \
--genomeDir hg38_index \
--genomeFastaFiles Homo_sapiens_assembly38.fasta \
--sjdbGTFfile UCSC.hg38.gtf \
indexを作るときにGFFも使うとイントロン領域とかを上手く処理してくれるらしいです。なくても大丈夫です。
two pass Mapping
STAR --twopassMode Basic \
--genomeDir hg38_index \
--readFilesCommand "gunzip -c" \
--readFilesIn input.fastq.gz \
--outSAMtype BAM SortedByCoordinate \
--outSAMattrRGline ID:input LB:lib1 PL:illumina SM:input PU:unit3 \
--outFileNamePrefix output_
samtools index output_Aligned.SortedByCoordinate.out.bam
最近のSTARはtwopassModeをBasicにすればtwopassMappingしてくれます。また、GATKを使うためにはRGタグが必要なので適当に付けてます。重要なのはIDだけなので、そこだけサンプルごとに変更しておきましょう。ついでにindexをつけておきます。
3. PCR duplicatesの除去
PCR duplicatesを除去します。picardを使うのが普通そうです。samtools markdupでもできます。
picard \
MarkDuplicates \
I=$output_Aligned.SortedByCoordinate.out.bam \
O=output_mdedup.bam \
CREATE_INDEX=true \
VALIDATION_STRINGENCY=SILENT \
M=output.metrics
4. SplitNCigarReads
CigarのNはイントロン領域を意味します。SplitNCigarReadsではその領域でリードを分割してくれます。java-optionsは環境に合わせて設定してください。
gatk \
--java-options "-Xms32G -Xmx32G" \
SplitNCigarReads \
-R Homo_sapiens_assembly38.fasta \
-I output_mdedup.bam \
-O output_split.bam
5. BQSRの計算、適用
既知の変異データをもとに、塩基スコアの再計算を行います。その後、再計算したスコアをもとにBAMファイルを修正します。
gatk --java-options "-Xmx4G" \
BaseRecalibrator \
-R Homo_sapiens_assembly38.fasta \
-I output_split.bam \
--known-sites Homo_sapiens_assembly38.dbsnp138.vcf.gz \
--known-sites Homo_sapiens_assembly38.known_indels.vcf.gz \
--known-sites Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
-O output_recal.table
gatk --java-options "-Xmx4G" \
ApplyBQSR \
-R Homo_sapiens_assembly38.fasta \
-I output_split.bam \
-bqsr output_recal.table \
-O output_recal.bam
6. 変異の検出
HaplotypeCallerを使って変異の検出を行います。HaplotypeCallerは以下の手順で変異を検出しています。
- active regionの検出
- active regionのリードをde Brujin graphでローカルアセンブリしたあと、Smith-Watermanを使用したペアワイズアラインメントで潜在的な変異の場所を検出します。
- 最後に、PairHMMでそれぞれの場所のリードにペアワイズアラインメントを行い、最終的な変異の場所を決定します。
gatk --java-options "-Xmx4G" \
HaplotypeCaller \
-R Homo_sapiens_assembly38.fasta \
-I output_recal.bam \
-ERC GVCF \
--dbsnp Homo_sapiens_assembly38.dbsnp138.vcf.gz \
-O output.g.vcf
また、このパイプラインではGVCFファイルを出力します。GVCFファイルは変異が検出されなかった領域の情報を含んだVCFファイルの亜種です。詳しくはこちらを参照してください。
7. GVCFからVCFへの変換
GVCFファイルのままでもいいですが、基本的に欲しいのは変異のデータだけなのでVCFフォーマットに変換します。
gatk --java-options "-Xmx4G" \
GenotypeGVCFs \
-R Homo_sapiens_assembly38 \
-V output.g.vcf \
-O output.vcf
8. 変異のフィルタリング
genomic seqeuenceの解析を行う際には、既知変異データをもとに変異スコアを算出し、そのデータを使ってVCFファイルのフィルターを行います。しかし、RNA-seqに関してはそのデータが存在しないのでそのステップが存在しません。単純に今回算出したVCFのデータのみに基づいて、ハードフィルタリングを行います。
gatk --java-options "-Xmx4G" \
VariantFiltration \
-R Homo_sapiens_assembly38.fasta \
-V output.vcf \
--window 35 \
--cluster 3 \
--filter-name "FS" --filter "FS > 30.0" \
--filter-name "QD" --filter "QD < 2.0" \
-O output_filtered.vcf