OryzaSNP Pipeline
RAP-DBから公開されている多型検出パイプラインを実行します
Analysis tools
- SamTools (v1.9)
- GATK (v4.0.11.0)
- Picard (v2.18.17)
biocondaでインストールしたsamtools/bcftoolsでlibcrypto.so.1.0.0がないと言われる という問題が起こることがあるため、上記ページを参考にしてインストールします
通常のvcfファイルは変異が入っていない場所の情報しか含まれていない
個々のサンプルだけではなく多数のサンプルを合わせて解析する場合、変異が入っていない場所の情報も必要
gvcfファイルは全ての場所の情報が含まれている
詳しくはGATKの公式のページを参照してください。
1つのSNP/INDELを1行にマージすることを意味する
A.vcf.gz + B.vcf.gz -> Merged.vcf.gz
A.vcf.gzとB.vcf.gz の多型の位置が同じとは限らない
bcftools merge A.vcf.gz B.vcf.gz -Oz -o Merged.vcf.gz
A.vcf.gz B.vcf.gz : マージするファイル
Merged.vcf.gz : 出力ファイル
-Oz : 圧縮した.gz形式で保存
例えば各染色体に分かれているvcf ファイルを1つのvcf ファイルにまとめることを意味する
A.vcf.gz + B.vcf.gz -> output.vcf.gz
bcftools concat A.vcf.gz B.vcf.gz -Oz -o output.vcf.gz
A.vcf.gz B.vcf.gz : マージするファイル
Merged.vcf.gz : 出力ファイル
-Oz : 圧縮した.gz形式で保存
複数品種の多型情報をまとめたvcfファイルを作成します
HaplotypeCallerからgvcfファイルを作成します。
gatk --java-options "-Xmx64G -Xms64G" HaplotypeCaller --input alignment.rmdup.bam --output variants.g.vcf.gz --reference genome.fa -max-alternate-alleles 2 --emit-ref-confidence GVCF
計算に時間がかかるため、染色体ごとに分けて実行します
複数の gVCF ファイルからジェノタイピング用のローカルデータベースを構築
染色体ごとに実行します。for文を使う場合を示しています。
array=(chr01 chr02 chr03 chr04 chr05 chr06 chr07 chr08 chr09 chr10 chr11 chr12)
SAMPLE_IDを指定して実行します。
例)SAMPLE_ID=SAMPLE1
for name in ${array[@]}
do
echo $name
gatk --java-options "-Xmx64G -Xms64G" GenomicsDBImport --reference genome.fa -V variants_A.g.vcf.gz -V variants_B.g.vcf.gz -V Variants_C.g.vcf.gz --genomicsdb-workspace-path "$SAMPLE_ID"_"$name" --intervals "$name"
done
計算に時間がかかるため、染色体ごとに分けて実行します
GenomicsDBImportで作成されたデータを利用してVCF ファイル(複数サンプル)を作成
array=(chr01 chr02 chr03 chr04 chr05 chr06 chr07 chr08 chr09 chr10 chr11 chr12)
SAMPLE_IDを指定して実行します。
for name in ${array[@]}
do
echo $name
gatk --java-options "-Xmx64G -Xms64G" GenotypeGVCFs --reference genome.fa -V gendb://"$SAMPLE_ID"_"$name" -G StandardAnnotation --new-qual -O "$SAMPLE_ID"_"$name".vcf.gz
done
(参考情報)1サンプルだけに対して実行する場合
gatk GenotypeGVCFs --variant variants.g.vcf.gz --output variants.genotype.vcf.gz --reference genome.fa
GBSのサンプルを含む場合、条件を変更する必要があります。
for name in ${array[@]}; do echo $name;
gatk --java-options "-Xmx64G -Xms64G" VariantFiltration \
--reference genome.fa \
--variant "$SAMPLE_ID"_"$name".vcf.gz\
--output "$SAMPLE_ID"_"$name".filter.genotype.vcf.gz \
--filter-expression "QUAL < 100 || DP < 3.0 || QD < 20.0 || FS >= 40.0 || MQ < 30.0" \
--filter-name "FILTER"
done
for name in ${array[@]}; do echo $name;
gatk --java-options "-Xmx64G -Xms64G" SelectVariants \
--reference genome.fa \
--variant "$SAMPLE_ID"_"$name".filter.genotype.vcf.gz \
--output "$SAMPLE_ID"_"$name".varonly.vcf.gz \
--exclude-filtered \
--select-type-to-include SNP \
--select-type-to-include INDEL
done
これまで染色体ごとに分けて実行していましたが、１つのvcfファイルにまとめます
bcftools concat -o "$SAMPLE_ID".varonly.vcf.gz "$SAMPLE_ID"_chr01.varonly.vcf.gz "$SAMPLE_ID"_chr02.varonly.vcf.gz "$SAMPLE_ID"_chr03.varonly.vcf.gz "$SAMPLE_ID"_chr04.varonly.vcf.gz "$SAMPLE_ID"_chr05.varonly.vcf.gz "$SAMPLE_ID"_chr06.varonly.vcf.gz "$SAMPLE_ID"_chr07.varonly.vcf.gz "$SAMPLE_ID"_chr08.varonly.vcf.gz "$SAMPLE_ID"_chr09.varonly.vcf.gz "$SAMPLE_ID"_chr10.varonly.vcf.gz "$SAMPLE_ID"_chr11.varonly.vcf.gz "$SAMPLE_ID"_chr12.varonly.vcf.gz -O z
イネは自殖のため、ヘテロの遺伝子型は除いて解析に用います
zcat "$SAMPLE_ID".varonly.vcf.gz | sed 's/1[/]0/\.\/\./g' | sed 's/0[/]1/\.\/\./g' | bgzip -c > "$SAMPLE_ID".varonly.remHet.vcf.gz
Imputationを行う前に、信頼性の高いSNPを抽出します。
vcftools --gzvcf "$SAMPLE_ID".varonly.remHet.vcf.gz --max-alleles 2 --min-alleles 2 --maxDP 1000 --minDP 3 --minQ 20 --max-missing 0.95 --recode --out "$SAMPLE_ID".varonly.remHet.bi_MQ20_DP3-1000_MS0.95
beagleを用いて欠測を補完します。
欠測の補完については、以下の資料を参照してください。
beagle gt="$SAMPLE_ID".varonly.remHet.bi_MQ20_DP3-1000_MS0.95.recode.vcf.gz out="$SAMPLE_ID".varonly.remHet.bi_MQ20_DP3-1000_MS0.95.imputed
indexファイルを作成します
tabix -p vcf "$SAMPLE_ID".varonly.remHet.bi_MQ20_DP3-1000_MS0.95.imputed.vcf.gz