整個(gè)流程需要軟件：
bwa+picard+gatk+annovar
軟件的安裝與使用大家可以自行百度。。。。
代碼如下：
##質(zhì)控
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
#比對(duì)
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
#去重復(fù)
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
#顯示比對(duì)結(jié)果
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
#結(jié)合gvcf文件變?yōu)関cf文件
/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
##進(jìn)行單項(xiàng)注釋
annotate_variation.pl --geneanno --dbtype refGene --buildver hg19 control_filtered_snp.avinput /Data2/SNAKEMAKE/annovar/humandb/
##進(jìn)行聯(lián)合注釋
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 &
代碼僅供參考
請(qǐng)銘記馬克思主義活的靈魂：具體問題具體分析?。?！
