通過雙端測序的數(shù)據(jù)比對參考基因組贿堰。
第一步準(zhǔn)備工作（這里只用了擬南芥1號染色體的數(shù)據(jù)）
bwa index ./sequences.fa（bwa建立索引）
gatk --java-options "-Xmx8G -Djava.io.tmpdir=./tmp" CreateSequenceDictionary -R ath.chr1/sequences.fa -O ath.chr1/sequences.dict（建立gatk所需要的文件抖拦。-R 參考基因組所在位置担巩。-Xmx8G -Djava.io.tmpdir=./tmp" ，內(nèi)存最大8G阱当，存儲在tmp臨時文件，gatk CreateSequenceDictionary 創(chuàng)建字典文件）
samtools faidx ath.chr1/sequences.fa（建立samtools 提取序列）
第二步準(zhǔn)備比對
mkdir mapping
bwa mem -t 4 -R '@RG\tID:s1\tSM:s1\tLB:s1\tPL:ILLUMINA' ath.chr1/sequences.fa data/sample1.1.fq data/sample1.2.fq|samtools view -bS > mapping/s1.bam（mem適合100bp以上，-t 4 4個cpu笤受，-R 輸出文件 ID是樣品名稱，LB是庫的樣品敌蜂，測序平臺箩兽，參考基因組的地址，以及雙端測序的測序文件 |前面的輸出到后面的輸入紊册，sam變成bam文件）
bwa mem -t 4 -R '@RG\tID:s2\tSM:s2\tLB:s2\tPL:ILLUMINA' ath.chr1/sequences.fa data/sample2.1.fq data/sample2.2.fq|samtools view -bS > mapping/s2.bam
samtools sort -@ 4 mapping/s1.bam -o mapping/s1.sorted.bam（按照基因組順序排序）
samtools sort -@ 4 mapping/s2.bam -o mapping/s2.sorted.bam
gatk --java-options "-Xmx4G -Djava.io.tmpdir=./tmp" MarkDuplicates -I mapping/s1.sorted.bam -O mapping/s1.MarkDup.bam -M mapping/s1.markdup_metrics.txt（對pcr duplicate進(jìn)行標(biāo)記）
gatk --java-options "-Xmx4G -Djava.io.tmpdir=./tmp" MarkDuplicates -I mapping/s2.sorted.bam -O mapping/s2.MarkDup.bam -M mapping/s2.markdup_metrics.txt
samtools index mapping/s1.MarkDup.bam（建立index比肄，否則難以比對）
samtools index mapping/s2.MarkDup.bam
第三步snp calling
nohup gatk --java-options "-Xmx4G -Djava.io.tmpdir=./tmp" HaplotypeCaller -R ath.chr1/sequences.fa -I mapping/s1.MarkDup.bam -O SNP/s1.gvcf --emit-ref-confidence GVCF -stand-call-conf 30 >SNP/s1.HaplotypeCaller.log&（HaplotypeCaller對每個樣品分別進(jìn)行變異檢測，每個樣品生成一個gvcf文件囊陡，-stand-call-conf 30 質(zhì)量高于30才可以）
nohup gatk --java-options "-Xmx4G -Djava.io.tmpdir=./tmp" HaplotypeCaller -R ath.chr1/sequences.fa -I mapping/s2.MarkDup.bam -O SNP/s2.gvcf --emit-ref-confidence GVCF -stand-call-conf 30 >SNP/s2.HaplotypeCaller.log&
gatk --java-options "-Xmx4G -Djava.io.tmpdir=./tmp" GenomicsDBImport -L 1 -V SNP/s1.gvcf -V SNP/s2.gvcf --genomicsdb-workspace-path SNP/vcfdbGenomicsDBImport(將多個樣品的gvcf信息合并為一個數(shù)據(jù)庫#-L參數(shù)后面跟著染色體編號芳绩，如果有多條染色體，-L 可以輸入多次撞反，例如 -L 1 -L 2 -L 3表示1,2,3三條染色體)
gatk --java-options "-Xmx4G -Djava.io.tmpdir=./tmp" GenotypeGVCFs -R ath.chr1/sequences.fa? -V gendb://SNP/vcfdb -O SNP/raw.vcf >SNP/raw.GenotypeGvcf.log(轉(zhuǎn)為vcf文件)
gatk --java-options "-Xmx4G -Djava.io.tmpdir=./tmp" VariantFiltration -V SNP/raw.vcf -O SNP/filter.vcf --filter-expression "QD < 2.0 || FS > 60.0" --filter-name "Fail" >SNP/filter.log（按照QD與FS進(jìn)行過濾）
將InDel與SNP分開為兩個文件
grep -v "Fail" SNP/filter.vcf >SNP/final.vcf(保留通過文件)
gatk --java-options "-Xmx8G -Djava.io.tmpdir=./tmp" SelectVariants -select-type INDEL -V SNP/final.vcf -O SNP/final.indel.vcf
gatk --java-options "-Xmx8G -Djava.io.tmpdir=./tmp" SelectVariants -select-type SNP? -V SNP/final.vcf -O SNP/final.snp.vcf
snp注釋
snpEff build -c snpEff.config ath.chr1（建庫）
snpEff -c snpEff.config -o gatk ath.chr1 SNP/final.snp.vcf > demo_anno_vcf（snp 分析）?
