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 [email protected] 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 [email protected] 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