#构建ref genome index file
bwa index genome.fasta genome
vim sam_file.sh
#!/usr/bin/env bash
for i in {s1,d1,fuben,muben}
do
#将重测序的clean data 回帖到参考基因组上
bwa mem -t 28 -M ./data/genome.fasta ./data/${i}_1.fq ./data/${i}_2.fq >./${i}/${i}.sam
#转换为bam文件
samtools view -@ 28 -S -b ./${i}/${i}.sam >./${i}/${i}.bam
#统计比对结果
samtools flagstat ./${i}/${i}.bam >./${i}/${i}.log
#进行sort
samtools sort -@ 28 ./${i}/${i}.bam ./${i}/${i}.sort
done
nohup bash sam_file.sh &
#使用GATK4 call variants
vim GATK.sh
#!/usr/bin/bash
for i in {s1,d1,fuben,muben}
do
#将sort后的bam文件加上head信息
java -jar /biosoftware/picard-tools-1.134/picard.jar AddOrReplaceReadGroups I=${i}/${i}.sort.bam O=${i}/${i}_picard.bam SO=coordinate ID=${i} LB=${i} SM=${i} PL=illumina PU=run
#建立 bam索引和参考基因组的库文件
samtools index ${i}/${i}_picard.bam
java -jar /biosoftware/picard-tools-1.134/picard.jar CreateSequenceDictionary R=data/assembly.fasta O=data/assembly.dict
samtools faidx data/assembly.fasta
#过量扩增的reads并不是基因组自身固有序列，不能作为变异检测的证据，因此，要尽量去除这些由PCR扩增所形成的duplicates
java -jar /biosoftware/picard-tools-1.134/picard.jar MarkDuplicates I=./${i}/${i}_picard.bam O=./${i}/${i}.dedup.bam CREATE_INDEX=true VALIDATION_STRINGENCY=SILENT M=./${i}/${i}.metrics
#把所有的样本的bam文件都用HaplotypeCaller的GVCF模式分别分析变异
/data2/yanxu2016/gatk-4.0.0.0/gatk --java-options -Xmx60g HaplotypeCaller --emit-ref-confidence GVCF -R ./data/assembly.fasta -I ./${i}/${i}.dedup.bam -O ./gatk_out/${i}.HC.g.vcf
done
#合并gvcf
/data2/yanxu2016/gatk-4.0.0.0/gatk --java-options -Xmx60g CombineGVCFs -V ./gatk_out/fuben.HC.g.vcf -V ./gatk_out/muben.HC.g.vcf -V ./gatk_out/d1.HC.g.vcf -V ./gatk_out/s1.HC.g.vcf -R ./data/assembly.fasta -O ./gatk_out/Reseq_all_fmds.HC.vcf
#进行分型
/data2/yanxu2016/gatk-4.0.0.0/gatk --java-options -Xmx60g GenotypeGVCFs -V Reseq_all_dfsm.HC.vcf -R ../data/assembly.fasta -O all_Reseq_genetype.HC.vcf
#使用samtools call variants 综合GATK和samtools的结果获得可靠变异
#之前安装的是1.1版本的samtools，1.1版本缺少对基因型信息统计的部分信息的筛选功能，所以重新装1.9版本，总是报错，查到github的project页面发现有些库不全。
./samtools-1.9/samtools mpileup -t AD,ADF,ADR,DP,SP -ugf ../data/assembly.fasta ../fuben/fuben.dedup.bam ../muben/muben.dedup.bam ../d1/d1.dedup.bam ../s1/s1.dedup.bam | bcftools call -vm > variants.samtools.raw1.vcf
未完待续
