- 选择一个合适的参考基因组，最好是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处理
- 使用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
