GATK4.1 call SNP

GATK4.0 和之前的版本相比還是有較大的不同,更加趨于流程化。

軟件安裝

1 wget https://github.com/broadinstitute/gatk/releases/download/4.1.5.0/gatk-4.1.5.0.zip 
2 unzip gatk-4.1.5.0.zip
#該安裝步驟僅供參考,需要加入環(huán)境變量才可以直接調用

下載GATK

wget https://github.com/broadinstitute/gatk/releases/download/4.1.9.0/gatk-4.1.9.0.zip
unzip gatk-4.1.9.0.zip
#進入環(huán)境變量修改PATH,很重要
sudo vi /etc/profile
sudo gedit  /etc/profile
#添加下面的PATH
# export PATH=$PATH:/home/zhiwufy/biosoft/gatk-4.1.9.0
source /etc/profile

下載后解壓

image
# 在該目錄下
conda env create -n gatk -f gatkcondaenv.yml -y
# 本次安裝 -y參數(shù) 出錯,原因未知 去掉 -y

檢查是否成功

python -c "import vqsr_cnn"
# 出現(xiàn):Using TensorFlow backend,則成功
# 出現(xiàn):ImportError: No module named vqsr_cnn,則失敗

#報錯:ModuleNotFoundError: No module named 'numpy.testing.decorators'
#原因 numpy 版本不匹配!
pip uninstall numpy
pip install numpy==1.18 # > 1.18
pip install scipy==1.1.0
python -c "import vqsr_cnn"
#通過驗證

如果以后要更新gatk

需要先卸載掉原來的小環(huán)境,再重新安裝,否則會產(chǎn)生報錯,比如:

CondaValueError: prefix already exists: /anaconda2/envs/gatk

安裝的正確方法:

source deactivate
conda env remove -n gatk
# 然后重新運行上面的安裝代碼

另外gatk運行還需要依賴Java

先進入到gatk小環(huán)境,再安裝openjdk即可

apt-cache search openjdk
conda install -c conda-forge openjdk
sudo apt-get install openjdk-8-jdk

GATK 簡單說明

## 幫助信息
gat --help

## 列出所有的工具
gatk --list

## 工具的說明,比如以VariantAnnotator 為例
gatk VariantAnnotator --help

GATK分析簡要流程

  • 所需數(shù)據(jù) : ref.fa

      • reads1.fq
      • reads2.fq
  • 建立索引(本人使用的是hisat2,結果差異不大)

bwa index ref.fa
samtools  faidx ref.fa
gatk CreateSequenceDictionary -R ref.fa -O ref.dict
gatk CreateSequenceDictionary  -R Tifrunner2.fasta  
##
-R Input reference fasta or fasta.gz  Required
-O  輸出文件
  • 比對

## bwa 比對
bwa mem -t 4 -R '@RG\tID:id1\tPL:illumina\tSM:test' ref.fa test_1.fq test_2.fq | samtools view -bS - >test.bam

##參數(shù)
-R 設置reads group,gatk必須要的信息,其中ID,PL和SM信息是必須要的

## 排序
samtools sort -@ 3 -o test.sorted.bam test.bam
rm test.bam

2. 通過管道命令直接鏈接samtools

bowtie2 -x genome_index -1 input_1.fq -2 input_2.fq | samtools view -bS | samtools sort > output.bam

這條命令把bowtie2 生成的sam文件通過管道|傳遞到samtools,將sam轉換為bam文件,省去中間sam文件的空間占用

GATK 要求read group的格式

ID = Read group identifier

每一個read group 獨有的ID,每一對reads 均有一個獨特的ID,可以自定義命名;

PL = Platform

測序平臺;ILLUMINA, SOLID, LS454, HELICOS and PACBIO,不區(qū)分大小寫;

SM = sample

reads屬于的樣品名;SM要設定正確,因為GATK產(chǎn)生的VCF文件也使用這個名字;

LB = DNA preparation library identifier

對一個read group的reads進行重復序列標記時,需要使用LB來區(qū)分reads來自那條lane;有時候,同一個庫可能在不同的lane上完成測序;為了加以區(qū)分,
  同一個或不同庫只要是在不同的lane產(chǎn)生的reads都要單獨給一個ID. 一般無特殊說明,成對兒read屬于同一庫,可自定義,比如:library1

若是忘記添加read group信息還以通過 AddOrReplaceReadGroups 添加

sample name 也可以在markdup之后加
gatk AddOrReplaceReadGroups -I .bam -O .add.bam -LB library1 -PL illumina -PU pl1 -SM name

# 

##參數(shù)
-I Input file (BAM or SAM or a GA4GH url);
-O  Output file (BAM or SAM);
-LB Read-Group library;
-PL  Read-Group platform (e.g. ILLUMINA, SOLID);
-PU Read-Group platform unit (eg. run barcode);
-SM Read-Group sample name 

## 建立索引

samtools index test.sorted.markup.bam
  • 標記重復序列

gatk  MarkDuplicates -I test.sorted.bam -O test.sorted.markdup.bam -M test.sorted.markdup_metrics.txt3 ##參數(shù)
-I 排序后的一個或者多個bam或者sam文件
-M 輸出重復矩陣
-O 輸出文件
 
# 建立索引
samtools index test.sorted.markup.bam

  • 檢測變異

#在檢測變異前必須要對bam文件建立索引:samtools index 
 ##兩種方法

##(1)多樣本一起call,此次只有一個樣本,若有多個樣本,則繼續(xù)用 -I 參數(shù)添加即可
gatk --java-options -Xmx4G HaplotypeCaller -I test.sorted.markup.bam -O test.gvcf  -R ref.fa

## (2)單個樣本call,然后在合并
## 生成中間文件gvcf
gatk --java-options -Xmx4G HaplotypeCaller -I test.sorted.markup.bam -O test.gvcf -R ref.fa --emit-ref-confidence GVCF

##通過gvcf檢測變異, -V 添加上步得到的gvcf
gatk GenotypeGVCFs -R ref.fa -V test.gvcf -O test.vcf

##參數(shù)
-I BAM/SAM/CRAM file
-O  輸出文件
-R 參考基因組
--java-options: 若設置java則需要添加
-Xmx4G:內存為4G,防止內存太大
-V  A VCF file containing variants
GATK CombineGVCFs  -R GRCh38.fa -L chr10.bed --variant father_chr10.g.vcf.gz --variant mother_chr10.g.vcf.gz --variant child_chr10.g.vcf.gz -O family_chr10.g.vcf.gz

GATK CombineGVCFs  -R GRCh38.fa -L chr10.bed --variant father_chr10.g.vcf.gz --variant mother_chr10.g.vcf.gz --variant child_chr10.g.vcf.gz -O family_chr10.g.vcf.gz





  • 我用 GATK 模塊 CombineVariants 合并了 12 個樣本的 VCF 文件,之后發(fā)現(xiàn)部分位點 FORMAT 字段缺少了 AD 信息。仔細看了一下,這樣的位點都是有多個 ALT 的位點。
#正常情況下 FORMAT 字段:

GT:AD:DP:GQ:PL

#缺失 AD 的 FORMAT 字段:

GT:DP:GQ

#用 GATK 模塊 CombineVariants 合并多個樣本的 VCF 之后,DP 字段是會自動更新的,但是 AD 字段需要重新運行 VariantAnnotator , 根據(jù)合并的 GT 和 DP 重新生成新的 AD。

#所以多個樣本分別 call 變異然后進行合并需要注意這個問題。

gatk  VariantAnnotator    -R reference.fasta    -V input.vcf    -o output.vcf  

  • 提取SNP,INDEL

## 提取SNP
gatk SelectVariants -V test.vcf -O test.snp.vcf --select-type-to-include SNP

## 提取INDEL
gatk SelectVariants -V test.vcf -O test.indel.vcf --select-type-to-include INDEL

##參數(shù)
-O 輸出vcf文件
-V 輸入vcf文件
--select-type-to-include 選擇提取的變異類型{NO_VARIATION, SNP, MNP, INDEL,
                              SYMBOLIC, MIXED}
  • 對vcf文件進行過濾

gatk VariantFiltration -O test.snp.fil.vcf.temp -V test.snp.vcf --filter-expression 'QUAL < 30.0 || QD < 2.0 || FS > 60.0 ||  SOR > 4.0' \
    --filter-name lowQualFilter --cluster-window-size 10  --cluster-size 3 --missing-values-evaluate-as-failing

## 參數(shù)
-O 輸出filt.vcf文件
-V 輸入vcf文件
--filter-expression 過濾條件, VCF INFO 信息
--cluster-window-size 以10個堿基為一個窗口
--cluster-size 10個堿基為窗口,若存在3以上個則過濾
--filter-name 被過濾掉的SNP不會刪除,而是給一個標簽, 比如 Filter
--missing-values-evaluate-as-failing 當篩選標準比較多的時候,可能有一些位點沒有篩選條件當中的一條或幾條,例如下面的這個表達式;QUAL < 30.0 || QD < 2.0 || FS > 60.0 || MQ < 40.0 || HaplotypeScore > 13.0 并不一定所有位點都有這些信息,這種情況下GATK運行的時候會報很多WARNING信息,用這個參數(shù)可以把這些缺少某些FLAG的位點也給標記成沒有通過篩選的。


gatk VariantFiltration \
    -V $vcf \
    -filter "QD < 2.0" --filter-name "QD2" \
    -filter "QUAL < 30.0" --filter-name "QUAL30" \
    -filter "FS > 60.0" --filter-name "FS60" \
    -filter "MQ < 40.0" --filter-name "MQ40" \
    -filter "MQRankSum < -12.5" --filter-name "MQRankSum-12.5" \
    -filter "ReadPosRankSum < -8.0" --filter-name "ReadPosRankSum-8" \
    -O ${out_prefix}.vcf


#標出10bp范圍3個SNP的 ”SnpCluster“
gatk VariantFiltration -V 324.wgs.call.SNP.Filter-SOR.ANN.vcfFliter.vcf -cluster 3 -window 10 -O 324.wgs.call.SNP.Filter-SOR.ANN.vcfFliter.10-3filter.vcf
#去除上一步標出的SnpCluster"
gatk SelectVariants -V 324.wgs.call.SNP.Filter-SOR.ANN.vcfFliter.10-3filter.vcf -O 324.wgs.call.SNP.Filter-SOR.ANN.vcfFliter.10-3filter-2.vcf -select "FILTER == SnpCluster" --invertSelect


去除indel附近5bp范圍內的SNP

bcftools filter -g 5 -O v -o 1-SnpGap5.vcf ../324.wgs.PASS.ANN.vcf.gz
原文鏈接:https://blog.csdn.net/Gossie/article/details/109320960
  • 篩選PASS的SNP,INDEL

## 根據(jù)FILTER那列信息進行篩選
grep PASS test.snp.fil.vcf.temp >  test.snp.fil.vcf

bcftools view -f  "PASS"  zhuxing.merge.snpfiltraw.vcf >  zhuxing.snp.filter.pass.vcf

最后編輯于
?著作權歸作者所有,轉載或內容合作請聯(lián)系作者
【社區(qū)內容提示】社區(qū)部分內容疑似由AI輔助生成,瀏覽時請結合常識與多方信息審慎甄別。
平臺聲明:文章內容(如有圖片或視頻亦包括在內)由作者上傳并發(fā)布,文章內容僅代表作者本人觀點,簡書系信息發(fā)布平臺,僅提供信息存儲服務。

相關閱讀更多精彩內容

友情鏈接更多精彩內容