2019-02-16细菌基因组SNP calling

1. 选择一个合适的参考基因组，最好是NCBI公布的参考基因组，如Escherichia coli 选择CFT073菌株。使用bwa index建立索引

mkdir snp_calling # 创建项目工作目录
mkdir snp_calling/fqdata # 存放PE双末端clean data
bwa index CFT073.fasta # 将参考基因组放入snp_calling目录下，并在该目录下建立索引

2.使用bwa mem 进行比对。

mkdir align_bam # 创建比对结果存放目录
bwa mem -M snp_calling/fqdata/${sample}.1.fq.gz snp_calling/fqdata/${sample}.2.fq.gz \
 -o snp_calling/align_bam/${sample}.sam # 得到SAM格式的比对文件，启用-M参数，便于后续Picard处理

3. 使用samtools工具将SAM文件转换为BAM文件。

samtools view -S -b snp_calling/align_bam/${sample}.sam \
 -o snp_calling/align_bam/${sample}.bam

4.使用samtools工具，对比对结果进行排序。

samtools sort snp_calling/align_bam/${sample}.bam \ 
snp_calling/align_bam/${sample}_sorted.bam

5.使用Picard Tools检测比对质量。

mkdir snp_calling/alignment_check
java -jar picard.jar CollectGcBiasMetrics \
      -i snp_calling/align_bam/${sample}_sorted.bam \
      -o snp_calling/alignment_check/${sample}_gc_bias_metrics.txt \
      CHART=snp_calling/alignment_check/${sample}_gc_bias_metrics.pdf \
      S=snp_calling/alignment_check/${sample}_summary_metrics.txt \
      R=CFT073.fasta  # 得到GC_bias_map
java -jar picard.jar MeanQualityByCycle \
     -i snp_calling/align_bam/${sample}_sorted.bam \
     -o snp_calling/alignment_check/${sample}_mean_qual_by_cycle.txt \
      CHART=snp_calling/alignment_check/${sample}_mean_qual_by_cycle.pdf # 均值
java -jar picard.jar QualityScoreDistribution \
      -i snp_calling/align_bam/${sample}_sorted.bam \
      -o snp_calling/alignment_check/${sample}_qual_score_dist.txt \
      CHART=snp_calling/alignment_check/${sample}_qual_score_dist.pdf # 质量得分分布图

5.使用Picard Tools从比对结果中移除相同的reads。

MarkDuplicates

mkdir snp_calling/rmdup_out
java -jar picard.jar MarkDuplicates  REMOVE_DUPLICATES=true \ # 移除完全相同的reads
      -i  snp_calling/align_bam/${sample}_sorted.bam \
      -o snp_calling/rmdup_out/${sample}_sorted_rmdup.bam \ 
      -M snp_calling/rmdup_out/${sample}_marked_dup_metrics.txt # 重复reads统计

6.使用Picard Tools重新将reads比对工具，再一次建立索引。

samtools faidx CFT073.fasta # 在snp_calling目录下

6.使用samtools工具，将多个细菌基因组的sorted.bam文件合并，生成包含SNPs信息的VCF文件。

samtools mpileup -f  CFT073.fasta /w/user238/snp_calling/escherichia_coli/alined_sorted/*.bam -o raw_var.bcf 
java -jar ~/soft/bin/VarScan.v2.3.9.jar mpileup2cns /w/user238/snp_calling/escherichia_coli/raw_var.bcf --min-coverage 2 --min-var-freq 0.8 --p-value 0.005 --variants --output-vcf -o /w/user238/snp_calling/escherichia_coli/variants.vcf