大肠杆菌全基因组重测序变异检测小实例（侧重变异过滤）

wgsim -N 80000 -1 150 -2 150 ../Reference_genome/ecoli.fa sim_1_reads_R1.fastq sim_1_reads_R2.fastq

wgsim -N 80000 -1 150 -2 150 ../Reference_genome/ecoli.fa sim_2_reads_R1.fastq sim_2_reads_R2.fastq

wgsim -N 80000 -1 150 -2 150 ../Reference_genome/ecoli.fa sim_3_reads_R1.fastq sim_3_reads_R2.fastq

###构建索引
bwa index Reference_genome/ecoli.fa
bwa mem -t 4 -R '@RG\tID:foo\tSM:sample1' Reference_genome/ecoli.fa sim_reads/sim_1_reads_R1.fastq sim_reads/sim_1_reads_R2.fastq -o output_results/sim_1_aligned.sam
bwa mem -t 4 -R '@RG\tID:foo\tSM:sample2' Reference_genome/ecoli.fa sim_reads/sim_2_reads_R1.fastq sim_reads/sim_2_reads_R2.fastq -o output_results/sim_2_aligned.sam
bwa mem -t 4 -R '@RG\tID:foo\tSM:sample3' Reference_genome/ecoli.fa sim_reads/sim_3_reads_R1.fastq sim_reads/sim_3_reads_R2.fastq -o output_results/sim_3_aligned.sam

cd output_results
#SAM装换为BAM
samtools view -S -b -o sim_1_aligned.bam sim_1_aligned.sam 
samtools view -S -b -o sim_2_aligned.bam sim_2_aligned.sam 
samtools view -S -b -o sim_3_aligned.bam sim_3_aligned.sam 
#排序
samtools sort sim_1_aligned.bam -o sim_1_aligned.sorted.bam
samtools sort sim_2_aligned.bam -o sim_2_aligned.sorted.bam
samtools sort sim_3_aligned.bam -o sim_3_aligned.sorted.bam
#索引
samtools index sim_1_aligned.sorted.bam 
samtools index sim_2_aligned.sorted.bam 
samtools index sim_3_aligned.sorted.bam
#变异检测
time samtools mpileup -g -t DP,AD -f ../Reference_genome/ecoli.fa sim_1_aligned.sorted.bam sim_2_aligned.sorted.bam sim_3_aligned.sorted.bam > sim_variants_3sample.bcf

###其一
time bcftools call -v -c sim_variants_3sample.bcf > sim_variants_3sample.vcf
###其二
time bcftools call -f GQ,GP -vmO z sim_variants_3sample.bcf -o sim_variants_3sample_1.vcf.gz

bcftools view -H sim_variants_3sample.vcf | wc -l
6918

cd ../
mkdir vcf_handling
cd vcf_handling
vcftools --vcf ../output_results/sim_variants_3sample.vcf --freq2 --out sim_variant_AF

vcftools --vcf ../output_results/sim_variants_3sample.vcf  --depth --out sim_variant_depth

vcftools --vcf ../output_results/sim_variants_3sample.vcf --site-mean-depth --out sim_variant_sitedepth

vcftools --vcf ../output_results/sim_variants_3sample.vcf --site-quality --out sim_variant_sitequality

vcftools --vcf ../output_results/sim_variants_3sample.vcf --missing-indv --out sim_variant_missingindv

sim_variant_missingindv.imiss
sim_variant_sitequality.lqual
sim_variant_sitedepth.ldepth.mean
sim_variant_depth.idepth
sim_variant_AF.frq

setwd("../../vcf_handling/")
library(tidyverse)
var_qual<-read_delim("sim_variant_sitequality.lqual",
                     delim="\t",
                     col_names=c("chr","pos","qual"),
                     skip=1)
library(ggplot2)
a<-ggplot(var_qual,aes(qual))+
  geom_density(fill="dodgerblue1")+
  theme_light()
a1<-a+xlim(0,100)
library(ggpubr)
ggarrange(a,a1,ncol=1,nrow=2)

summary(var_depth$mean_depth)
var_depth<-read_delim("sim_variant_sitedepth.ldepth.mean",
                      delim="\t",col_names=c("chr","pos","mean_depth","var_depth"),
                      skip=1)
a<-ggplot(var_depth,aes(mean_depth))+
  geom_density(fill="dodgerblue1",color="black",alpha=0.3)+
  theme_light()
a

INDV    N_DATA  N_GENOTYPES_FILTERED    N_MISS  F_MISS
sample1 6918    0       0       0
sample2 6918    0       0       0
sample3 6918    0       0       0

var_freq<-read_delim("sim_variant_AF.frq",delim="\t",
                     col_names=c("chr","pos","nalleles","nchr","a1","a2"),skip=1)
var_freq
var_freq$maf<-var_freq%>%
  select(a1,a2)%>%
  apply(1,function(z) min(z))
a<-ggplot(var_freq,aes(maf))+
  geom_density(fill="dodgerblue1",alpha=0.3)+
  theme_light()
a

vcftools --vcf ../output_results/sim_variants_3sample.vcf --minQ 30 --min-meanDP 4 --max-meanDP 10 --minDP 4 --maxDP 10 --recode --stdout | gzip -c > sim_variants_filtered.vcf

bcftools view -H sim_variants_filtered.vcf.gz | wc -l
2212