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

试验数据
fasta:hg38.fa文件可以在UCSC下载 (hg38.fa.gz 938M)
fastq非公开文件

KY18011403DNA_DHG18153-V_AHHVVHCCXY_L7_1.fq  35G
KY18011403DNA_DHG18153-V_AHHVVHCCXY_L7_2.fq  35G
KY18011403DNA_DHG18153-V_AHHVVHCCXY_L8_1.fq  36G
KY18011403DNA_DHG18153-V_AHHVVHCCXY_L8_2.fq  36G
KY18011403DNA_DHG18153-V_BHHVFHCCXY_L4_1.fq  25G
KY18011403DNA_DHG18153-V_BHHVFHCCXY_L4_2.fq  25G
KY18011403DNA_DHG18153-V_BHHVFHCCXY_L5_1.fq  25G
KY18011403DNA_DHG18153-V_BHHVFHCCXY_L5_2.fq  25G

L7_1.fq中的7表示lane的编号(PS:测试带编号)

软件的版本：
bwa-0.7.17 gatk-4.0.2.1 samtools-1.5

试验过程：
1.生成索引文件

bwa index -a bwtsw hg38.fa

构建索引时需要注意的问题：bwa构建索引有两种算法，通过参数-a is 和-a bwtsw进行选择，这里使用的是整个hg38，因此使用bwtsw算法。(bwtsw Algorithm implemented in BWT-SW. This method works with the whole human genome）

2.生成sam文件
bwa有三种比对模式
aln/samse/sampe 、 bwasw 、 mem
第1种适合illumina的100bp以下的序列，后2种适合70bp-1Mbp的序列(具体细节)，此次fastq文件的长度为154bp，因此使用mem进行比对。

为了简化命令，下面使用环境变量代替对应的文件路径：

export hg38_fa=/share/share/data/WGS-2018-2/hg38_fa
export fq_1=/share/share/data/WGS-2018-2/KY18011403DNA/KY18011403DNA_DHG18153-V_AHHVVHCCXY_L7_1.fq
export fq_2=/share/share/data/WGS-2018-2/KY18011403DNA/KY18011403DNA_DHG18153-V_AHHVVHCCXY_L7_2.fq

bwa mem $hg38_fa $fq_1 $fq_2  -t 15 -M -R "@RG\tID:<KY_L7>\tLB:<KY>\tSM:<KY_L7_h3>\tPL:ILLUMINA" -o KY_L7.sam

参数详解：
R——设置reads标头，用"\t"分割，例如：’@RG\tID:foo\tSM:bar’
t——设置线程数
M——将较短的split hits标记为secondary，与picard的markDuplicates兼容

3.对sam文件进行重排序
3.1 成fa文件的字典文件，在fa文件所在目录运行，否则会提示’ReorderSam No reference sequence dictionary found’

samtools faidx $hg38_fa
gatk CreateSequenceDictionary -R $hg38_fa -O hg38.dict

3.2排序

gatk --java-options "-Xmx8G" ReorderSam -I $result_dir/KY3_L7.sam -O $result_dir/KY3_L7.sorted.sam -R $hg38_fa

4.sam转bam

samtools view -@ 15 -bS KY3_L7.sorted.sam -o KY3_L7.bam

bam是二进制文件,运算速度快,这里使用15线程跑

5.对bam文件进行排序

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

6.合并bam文件
使用相同的方法生成不同lane的.sorted.bam文件,比如KY3_L4.sorted.bam KY3_L5.sorted.bam (这一步在之前小数据集上没有)

KY3_L7.sorted.bam  KY3_L8.sorted.bam
gatk --java-options "-Xmx8G" MergeSamFiles --USE_THREADING true -O KY3.merged.bam -I KY3_L4.sorted.bam \
 -I KY3_L5.sorted.bam  -I KY3_L7.sorted.bam  -I KY3_L8.sorted.bam

7.Duplicates Marking
测序原理是随机打断，那么理论上出现两条完全相同的read的概率是非常低的，而且建库时PCR扩增存在偏向性，因此标出完全相同的read
(METRICS_FILE是输出文件)

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

8.bam文件生成索引

samtools index -@ 8 KY3.dup.bam

9.GATK变异检测
9.1生成raw vcf 文件

gatk --java-options "-Xmx8G" HaplotypeCaller -R $hg38_fa -I KY3.dup.bam  -O KY3.raw.vcf

最终生成的vcf文件大概1G左右

-rw-r--r-- 1 javis default_os_group 1001M Mar 16 14:14 KY3.raw.vcf

内容也和之前差不多

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  <KY3>
chr1    13418   .       G       A       21.80   .       AC=1;AF=0.500;AN=2;BaseQRankSum=0.347;ClippingRankSum=0.000;DP=43;ExcessHet=3.0103;FS=28.059;MLEAC=1;MLEAF=0.500;MQ=32.58;MQRankSum=-3.321;QD=0.51;ReadPosRankSum=0.162;SOR=4.163       GT:AD:DP:GQ:PL  0/1:36,7:43:50:50,0,1103
chr1    13813   .       T       G       217.77  .       AC=1;AF=0.500;AN=2;BaseQRankSum=0.856;ClippingRankSum=0.000;DP=17;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=23.28;MQRankSum=-3.310;QD=12.81;ReadPosRankSum=1.307;SOR=1.270       GT:AD:DP:GQ:PL  0/1:10,7:17:99:246,0,447
chr1    13838   .       C       T       205.77  .       AC=1;AF=0.500;AN=2;BaseQRankSum=1.070;ClippingRankSum=0.000;DP=20;ExcessHet=3.0103;FS=5.727;MLEAC=1;MLEAF=0.500;MQ=23.41;MQRankSum=-2.966;QD=10.29;ReadPosRankSum=0.460;SOR=2.409       GT:AD:DP:GQ:PL  0/1:14,6:20:99:234,0,541

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

猜你喜欢