BWA0.7+Samtools1.5+GATK4.0在小数据集上的试验

试验数据
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 -@ 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

BWA0.7+Samtools1.5+GATK4.0在小数据集上的试验

猜你喜欢