数据集上的试验"/>
BWA0.7+Samtools1.5+GATK4.0在小数据集上的试验
试验数据
chr14_1.fastq chr14_2.fastq (1.47G per .gz)
chr14.fasta (28M per .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 ,<SAMPLE_NAME> 更改成相应样品名称,否则将会被当成一个样品处理
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转bam
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:
.html
更多推荐
BWA0.7+Samtools1.5+GATK4.0在小数据集上的试验
发布评论