代码
bowtie2 -x ../reference/lambda_virus -1 reads_1.fq -2 reads_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条目的各自的平均测序深度
cat tmp.vcf|grep -v "^#"|grep INDEL|cut -f8|cut -d";" -f4|cut -d"=" -f2|awk '{ sum += $0; } END { print sum/NR }'
3.把INDEL条目再区分成insertion和deletion情况
cat tmp.vcf|grep -v "^#"|grep INDEL|awk '{if(length($4)>length($5)) print }' |less -S
cat 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.txt
cat > snp_perc.sh
cat uniq.txt|while read line
do
total=`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`%
done
bash snp_perc.sh
5.找到基因型不是 1/1 的条目，个数
cat tmp.vcf|grep -v '1/1'|grep -v '^#'|less -S
cat 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.txt
cat > target_20.sh
cat linenum.txt|while read line
do
echo $line
cat tmp.vcf|grep -v '^#'|sed -n ${line}p >> target_20.txt
done
bash target_20.sh
7.筛选变异位点质量值大于30的条目
cat tmp.vcf|grep -v "^#"|awk '{if($6>30)print}'|less -S
8.组合筛选变异位点质量值大于30并且深度大于20的条目
cat target_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
cat tmp.vcf|grep -v '^#'|cut -f 8|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.txt
paste -d -s tmp.txt AF.txt |less -S
