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 Copy 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 Copy dictファイルの作成方法 picard CreateSequenceDictionary R=genome.fa O=genome.dict Copy 2. STARを使ったMapping Indexの準備 STAR \ --runMode genomeGenerate \ --genomeDir hg38_index \ --genomeFastaFiles Homo_sapiens_assembly38.fasta \ --sjdbGTFfile UCSC.hg38.gtf \ Copy 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 Copy 最近の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 Copy 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 Copy 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 Copy gatk --java-options "-Xmx4G" \ ApplyBQSR \ -R Homo_sapiens_assembly38.fasta \ -I output_split.bam \ -bqsr output_recal.table \ -O output_recal.bam Copy 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 Copy また、このパイプラインではGVCFファイルを出力します。GVCFファイルは変異が検出されなかった領域の情報を含んだVCFファイルの亜種です。詳しくはこちらを参照してください。 7. GVCFからVCFへの変換 GVCFファイルのままでもいいですが、基本的に欲しいのは変異のデータだけなのでVCFフォーマットに変換します。 gatk --java-options "-Xmx4G" \ GenotypeGVCFs \ -R Homo_sapiens_assembly38 \ -V output.g.vcf \ -O output.vcf Copy 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 Copy

この記事に関するIssueをGithubで作成する

gatk4のRNA-seq bestpractice

TL;DR

全体の流れ

環境

1. 必要なファイルの準備

必要なファイルのダウンロード

faiの作成方法

dictファイルの作成方法

2. STARを使ったMapping

Indexの準備

two pass Mapping

3. PCR duplicatesの除去

4. SplitNCigarReads

5. BQSRの計算、適用

6. 変異の検出

7. GVCFからVCFへの変換

8. 変異のフィルタリング

Read Next