整个流程需要软件：
bwa+picard+gatk+annovar
软件的安装与使用大家可以自行百度。。。。
代码如下：
##质控
fastp -i /Data2/SNAKEMAKE/Data/NA12878-I33_case_R1.fastq.gz -o NA12878-I33_case_R1_clean.fastq.gz -I NA12878-I33_case_R2.fastq.gz -O NA12878-I33_case_R2_clean.fastq.gz
#比对
bwa mem -M ref.fasta -t 2 Z1_R1.fq Z1_R2.fq > aligned_reads_Z1.sam
#排序
java -Xmx5g -jar picard.jar SortSam I=aligned_reads_Z1.sam O=sorted_reads_Z1.bam SORT_ORDER=coordinate
#加头
java -Xmx5g -jar picard.jar AddOrReplaceReadGroups I=sorted_reads_Z1.bam O=sorted_reads_head_Z1.bam ID=group1 SM=Z1 PL=illumina LB=lib1 PU=unit1
#建立索引
java -Xmx5g -jar picard.jar BuildBamIndex INPUT=sorted_reads_head_Z1.bam
#去重复
java -Xmx5g -jar picard.jar MarkDuplicates I=sorted_reads_head_Z1.bam O=Z1.dedup_reads.bam M=Z1.marked_dup_metrics.txt VALIDATION_STRINGENCY=LENIENT MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=1000
#建立新bam索引
java -Xmx5g -jar picard.jar BuildBamIndex INPUT=Z1.dedup_reads.bam
#显示比对结果
samtools flagstat Z1.dedup_reads.bam > map_Z1.stat
#生成gvcf文件
gatk HaplotypeCaller --native-pair-hmm-threads 2 -R REF.fasta -I Z1.dedup_reads.bam -pairHMM FASTEST_AVAILABLE --emit-ref-confidence GVCF --create-output-variant-index True -stand-call-conf 30 --genotyping-mode DISCOVERY -O Z1.raw.gvcf
##生成vcf文件
gatk HaplotypeCaller --native-pair-hmm-threads 2 -R REF.fasta -I Z1.dedup_reads.bam --create-output-variant-index True -stand-call-conf 30 --genotyping-mode DISCOVERY -O Z1.raw.vcf
#结合gvcf文件变为vcf文件
/gatk CombineGVCFs -R REF.fasta --variant Z1.raw.gvcf --variant Z2.raw.gvcf -O combine_variants.vcf
#操作
gatk GenotypeGVCFs -R ../../ref/P.fulvidraco.female.genome.fasta --variant combine_variants.vcf -O raw_variants.vcf
#输出snp
gatk SelectVariants -R REF.fasta -V raw_variants.vcf -select-type SNP -O raw_snp.vcf
#输出indel
gatk SelectVariants -R REF.fasta -V raw_variants.vcf -select-type INDEL -O raw_indel.vcf
#过滤snp和indel
gatk VariantFiltration -R REF.fasta -V raw_snp.vcf --filter-expression "QD < 2.0 || FS > 60.0 || MQ < 40.0 || MQRankSum<-12.5 || ReadPosRankSum < -8.0" --filter-name "filter" -O filtered_snp.vcf
gatk VariantFiltration -R REF.fasta -V raw_indel.vcf --filter-expression "QD < 2.0 || FS > 200.0 || ReadPosRankSum <-20.0" --filter-name "filter" -O filtered_indel.vcf
##生成annovar输入文件
convert2annovar.pl -format vcf4 filtered.vcf > filtered.avinput
##进行单项注释
annotate_variation.pl --geneanno --dbtype refGene --buildver hg19 control_filtered_snp.avinput /Data2/SNAKEMAKE/annovar/humandb/
##进行联合注释
table_annovar.pl case_filtered_snp.avinput /Data2/SNAKEMAKE/annovar/humandb/ -buildver hg19 -out case_snp -remove -protocol refGene,cosmic70,snp142,avsift,esp5400_ea,esp5400_all -operation g,f,f,f,f,f -nastring NA &
代码仅供参考
请铭记马克思主义活的灵魂：具体问题具体分析！！！
