試驗資料
chr14_1.fastq chr14_2.fastq (1.47G each one .gz)
chr14.fasta (28M .gz)
chr14.fastq檔案可以在GAGE下載
chr14.fasta檔案可以在UCSC下載

軟體的版本：
bwa-0.7.17 gatk-4.0.2.1 samtools-1.5

試驗過程

1.生成索引檔案

bwa index chr14.fa

構建索引時需要注意的問題：bwa構建索引有兩種演算法，兩種演算法都是基於BWT的，這兩種演算法通過引數-a is 和-a bwtsw進行選擇。其中-a bwtsw對於短的參考序列是不工作的，必須要大於等於10Mb；-a is是預設引數，這個引數不適用於大的參考序列，必須要小於等於2G。
所以這裡使用預設引數即可

2.尋找輸入reads檔案的SA座標

bwa aln chr14_fa/chr14.fa chr14_1.fq -t 8 > chr14_1.sai
bwa aln chr14_fa/chr14.fa chr14_2.fq -t 8 > chr14_2.sai

3.生成sam,並新增標頭檔案@RG

bwa sampe chr14_fa/chr14.fa -r "@RG\tID:<chr14>\tLB:<human_chr14>\tSM:<human_no.1_chr14>\tPL:ILLUMINA" chr14_1.sai chr14_2.sai
 
 chr14_1.fq chr14_2.fq > chr14.sam

如果多樣品call SNP , 更改成相應樣品名稱，否則將會被當成一個樣品處理

4.對sam檔案進行重排序
4.1 成fa檔案的字典檔案，在fa檔案所在目錄執行，否則會提示’ReorderSam No reference sequence dictionary found’

samtools faidx chr14.fa
gatk CreateSequenceDictionary -R chr14_fa/chr14.fa -O chr14.dict

4.2排序

gatk --java-options "-Xmx8G"
 
 ReorderSam -I chr14.sam -O chr14.sorted.sam -R chr14_fa/chr14.fa

5.sam轉bamw

samtools view [email protected] 8 -bS  chr14.sam -o chr14.bam

bam是二進位制檔案,運算速度快,這裡使用8執行緒跑

6.對bam檔案進行排序

gatk --java-options "-Xmx8G" SortSam -I chr14.bam  --SORT_ORDER coordinate -O chr14.sorted.bam

7.Duplicates Marking
測序原理是隨機打斷，那麼理論上出現兩條完全相同的read的概率是非常低的，而且建庫時PCR擴增存在偏向性，因此標出完全相同的read
(METRICS_FILE是輸出檔案)

gatk --java-options "-Xmx8G" MarkDuplicates -I chr14.sorted.bam --REMOVE_DUPLICATES false --MAX_FILE_HANDLES_FOR_READ_ENDS_MAP 8000 \
-O chr14.sorted.dup.bam --METRICS_FILE chr14.sorted.metrics

8.bam檔案生成索引
samtools index chr14.sorted.dup.bam

9.GATK變異檢測
9.1生成raw vcf 檔案

gatk --java-options "-Xmx8G" HaplotypeCaller -R chr14_fa/chr14.fa -I chr14.sorted.dup.bam -O chr14.raw.vcf 

生成的vcf檔案張這樣：

#source=HaplotypeCaller
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  <human_no.1_chr14>
chr14   16016893        .       GAC     G       18.76   .       AC=2;AF=1.00;AN=2;DP=2;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=
60.00;QD=9.38;SOR=0.693  GT:AD:DP:GQ:PL  1/1:0,2:2:6:55,6,0
chr14   16032753        .       A       G       68.77   .       AC=1;AF=0.500;AN=2;BaseQRankSum=-0.508;ClippingRankSum=0.000;DP=9;Exces
sHet=3.0103;FS=3.522;MLEAC=1;MLEAF=0.500;MQ=55.11;MQRankSum=-1.383;QD=7.64;ReadPosRankSum=-0.765;SOR=0.105       GT:AD:DP:GQ:PL  0/1:5,
4:9:97:97,0,185
chr14   16032848        .       C       T       133.77  .       AC=1;AF=0.500;AN=2;BaseQRankSum=0.153;ClippingRankSum=0.000;DP=13;Exces
sHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=53.32;MQRankSum=-2.308;QD=10.29;ReadPosRankSum=-1.949;SOR=0.527      GT:AD:DP:GQ:PL  0/1:7,
6:13:99:162,0,231
chr14   16032862        .       A       G       14.91   .       AC=1;AF=0.500;AN=2;BaseQRankSum=0.792;ClippingRankSum=0.000;DP=8;Excess
Het=3.0103;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=54.11;MQRankSum=-2.100;QD=1.86;ReadPosRankSum=-0.464;SOR=0.307        GT:AD:DP:GQ:PL  0/1:6,
2:8:43:43,0,223
chr14   16032864        .       G       A       186.77  .       AC=1;AF=0.500;AN=2;BaseQRankSum=-0.060;ClippingRankSum=0.000;DP=7;Exces
sHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=57.19;MQRankSum=-0.366;QD=26.68;ReadPosRankSum=0.876;SOR=1.022       GT:AD:DP:GQ:PL  0/1:2,
5:7:58:215,0,58
chr14   16032877        .       A       T       158.82  .       AC=1;AF=0.500;AN=2;BaseQRankSum=-0.366;ClippingRankSum=0.000;DP=7;Exces
sHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=57.19;MQRankSum=-0.180;QD=22.69;ReadPosRankSum=0.000;SOR=0.527       GT:AD:DP:GQ:PL  0/1:1,
6:7:19:187,0,19

9.2select SNP

gatk --java-options "-Xmx8G"  SelectVariants -R chr14_fa/chr14.fa -select-type SNP --variant chr14.raw.vcf -O chr14_snp.vcf 

9.3.select indel

gatk --java-options "-Xmx8G"  SelectVariants -R chr14_fa/chr14.fa -select-type INDEL --variant chr14.raw.vcf -O chr14_indel.vcf

9.2和9.3按需執行

reference:
https://www.jianshu.com/p/938d362fc48d
https://www.plob.org/article/7009.html