来源：http://ddrv.cn/a/185068
代码
bowtie2 -x../reference/lambda_virus -1reads_1.fq -2reads_2.fq|samtools sort -@ 5 -o tmp.bam -
bcftools mpileup -f ../reference/lambda_virus.fa tmp.bam |bcftools call -vm > tmp.vcf
VCF 4.2格式中，关于SNP和INDEL，已没有 Dels 这个标签，无INDEL即为SNP；
1.把突变记录的vcf文件区分成 INDEL和SNP条目
cat tmp.vcf|grep'INDEL'|less -S cat tmp.vcf|grep-v'INDEL'|less -S
2.统计INDEL和SNP条目的各自的平均测序深度
cattmp.vcf|grep -v"^#"|grep INDEL|cut -f8|cut -d";"-f4|cut -d"="-f2|awk'{ sum +=$0; } END { print sum/NR }'
3.把INDEL条目再区分成insertion和deletion情况
cattmp.vcf|grep -v"^#"|grep INDEL|awk'{if(length($4)>length($5)) print }'|less -Scat tmp.vcf|grep -v"^#"|grep INDEL|awk'{if(length($5)>length($4)) print }'|less -S
4.统计SNP条目的突变组合分布频率
cat tmp.vcf|grep -v'^#'|grep -v'INDEL'|cut -f 4,5|sed's/\s\+//g'|sort|uniq >uniq.txtcat > snp_perc.shcat uniq.txt|whilereadlinedototal=`cat tmp.vcf|grep -v'^#'|grep -v'INDEL'|cut -f 4,5|wc -l`sub=`cat tmp.vcf|grep -v'^#'|grep -v'INDEL'|cut -f 4,5|sed's/\s\+//g'|grep$line|wc -l`echo$line$sub`echo"scale=2;$sub/$total*100"|bc`%donebash snp_perc.sh
5.找到基因型不是 1/1 的条目，个数
cat tmp.vcf|grep-v'1/1'|grep-v'^#'|less -Scat tmp.vcf|grep-v'1/1'|grep-v'^#'|wc -l
6.筛选测序深度大于20的条目
#####注意本题在脚本中的输出方式应该是>>，如果变量的使用不熟练，注意看脚本中的颜色变化cat tmp.vcf|grep-o -E'DP='[0-9]+|grep-Eo'[0-9]+'|awk'{if($0>20) print NR}'> linenum.txtcat > target_20.shcat linenum.txt|whilereadlinedoecho $line cat tmp.vcf|grep-v'^#'|sed -n ${line}p >> target_20.txtdonebash target_20.sh
7.筛选变异位点质量值大于30的条目
cattmp.vcf|grep -v"^#"|awk'{if($6>30)print}'|less -S
8.组合筛选变异位点质量值大于30并且深度大于20的条目
cattarget_20.txt|grep -v"^#"|awk'{if($6>30)print}'|less -S
9.理解DP4=4,7,11,18 这样的字段，就是 Number of high-quality ref-forward , ref-reverse, alt-forward and alt-reverse bases 计算每个变异位点的 AF
cattmp.vcf|grep -v'^#'|cut -f8|grep -Eo'DP4=[0-9]+,[0-9]+,[0-9]+,[0-9]+'|grep -Eo'[0-9]+,[0-9]+,[0-9]+,[0-9]+'|awk'BEGIN{FS=","}{print ($3+$4)/($1+$2+$3+$4)}'> AF.txtpaste -d -s tmp.txt AF.txt |less -S
10.在前面步骤的bam文件里面找到这个vcf文件的某一个突变位点的测序深度表明的那些reads，并且在IGV里面可视化bam和vcf定位到该变异位点。
