10.16 梳理SNP calling过程

之前跟着流程过了一遍，使用：

E.coli_K12_MG1655.fa 为参考基因组，

SRR1770413_1.fastq.gz  SRR1770419_1.fastq.gz SRR1770413_2.fastq.gz  SRR1770419_2.fastq.gz 为测序文件。

但是输出结果存在问题，今天来梳理一遍流程，看看是哪里出了问题

---

1.trimmomatic 除去接头序列和质量差的序列

2.（可做：fastqc检测序列质量）

3.索引参考基因组E.coli_K12_MG1655.fa

    bwa index E.coli_K12_MG1655.fa

4.长read使用bwa mem 将read 比对到基因组上

$bwa mem -t 80 -M -Y -R "@RG\tID:foo_lane\tPL:ILLUMINA\tSM:$sample" $reference $fq1 $fq2  > ${sample}.sam

5.sam ->bam

samtools view -@ 60 -S -b mut.sam > mut.bam

(或者4.5合并到一个语句：

$bwa mem -t 80 -M -Y -R "@RG\tID:foo_lane\tPL:ILLUMINA\tSM:$sample" $reference $fq1 $fq2 | samtools view -Sb - > ./out/bamout/${sample}.bam）

6.samtools 排序

7.GATK标记重复序列

8.samtools 为 out/bamout/${sample}.sorted.markdup.bam建立索引文件

9. gatk CreateSequenceDictionary为参考基因组建立dict

10. gatk HaplotypeCaller 进行snp calling 输出gvcf文件

11. 基因群体的joint calling：

gatk CombineGVCFs 先合并gvcf文件

gatk GenotypeGVCFs

其实还包括：

7.8步之间可以进行重比对和BQSR

11步后进行VQSR或是其他质量控制（师兄还说有imputation，就是把测序没有测出来到snp推测出来。）