前期准备,生成vcf格式

../../bowtie2 -x ../index/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

1.把突变记录的vcf文件区分成 INDEL和SNP条目

grep 'INDEL' tmp.vcf | grep -v '#'   # INDEL条目
grep -v -e 'INDEL' -v -e '#' tmp.vcf   # SNP条目

2.统计INDEL和SNP条目的各自的平均测序深度

cat tmp.vcf | grep -v "^#" | grep 'INDEL' | cut -f 8 | cut -d ";" -f 4 | awk '{print substr($0, 4)}' | awk '{sum += $0} END {print sum/NR}'   # INDEL条目
cat tmp.vcf | grep -v "^#" | grep -v 'INDEL' | cut -f 8 | cut -d ";" -f 1 | awk '{print substr($0, 4)}' | awk '{sum += $0} END {print sum/NR}'   # SNP条目

3.把INDEL条目再区分成insertion和deletion情况

grep -v '^#' tmp.vcf | grep 'INDEL' | awk '{if ($4 > $5) print $0}' | less -SN   
grep -v '^#' tmp.vcf | grep 'INDEL' | awk '{if ($4 < $5) print $0}' | less -SN

4.统计SNP条目的突变组合分布频率

grep -v '^#' tmp.vcf | grep -v 'INDEL' | cut -f 4-5 | tr '\t' '-' | sort | uniq -c

5.找到基因型不是 1/1 的条目，个数

grep -v '^#' tmp.vcf | awk '{if ($10 !~ "1/1") print}'
grep -v '^#' tmp.vcf | awk '{if ($10 !~ "1/1") print}' | wc -l

6.筛选测序深度大于20的条目

egrep 'DP=2[1-9]|DP=[3-9][0-9]' tmp.vcf | less -SN

7.筛选变异位点质量值大于30的条目

grep -v '^#' tmp.vcf | awk '{if ($6 > 30) print}' | less -SN

8.组合筛选变异位点质量值大于30并且深度大于20的条目

grep -v '^#' tmp.vcf | awk '{if ($6 > 30) print}' | egrep 'DP=2[1-9]|DP=[3-9][0-9]' | less -SN

9.理解DP4=4,7,11,18 这样的字段，就是 Number of high-quality ref-forward , ref-reverse, alt-forward and alt-reverse bases 计算每个变异位点的 AF

grep -v "^#" tmp.vcf | cut -f 8 | egrep -o 'DP4=[0-9]+,[0-9]+,[0-9]+,[0-9]+' | egrep -o '[0-9]+,[0-9]+,[0-9]+,[0-9]+' | awk -v FS="," '{print($3+$4)/($1+$2+$3+$4)}'

10.在前面步骤的bam文件里面找到这个vcf文件的某一个突变位点的测序深度表明的那些reads，并且在IGV里面可视化bam和vcf定位到该变异位点。

less -SN tmp.vcf
samtools view tmp.bam | awk '{if ($6 != "*" && $4 <= 1104 && substr($10, 1104-$4+1, 1) == "A") print}' | less -SN
samtools view tmp.bam | awk '{if ($6 != "*" && $4 <= 1344 && substr($10, 1344-$4+1, 1) == "T") print}' | less -SN

其他思考题

1.vcf的全称是什么？是用来记录什么信息的标准格式的文本？
Variant Call Format
记录序列的变异信息
2.一般选用哪个指令查看vcf文件，为什么不用vim
less -SN,文件太大了,不用vim处理
3.vcf文件以’##’开头的是什么信息？请认真查看这些信息。’#’开头的是什么信息？
'##'------注释信息
'#'----标题
4.Vcf文件除头信息，每行有多少列，请详细叙述每行的含义！请准确记忆
每行有8个固定列。第9列开始是基因型信息。
①.#CHROM
②.POS
③.ID
④.REF
⑤.ALT
⑥.QUAL
⑦；FILTER
⑧.INFO
5.理解format列和样本列的对应关系以及GT AD DP的含义。
format列和样本列是一一对应的,使用冒号分隔。
GT: 样本基因型。
AD: reads覆盖度，
DP: sample中该位点的覆盖度
6.Vcf文件第三列如果不是’.’，出现的rs号的id是什么？
dbSNP数据库中该SNP的ID号。
7.Vcf文件的ref，alt列和样本列的0/1 1/1 或者1/2的联系？
0 表示样品中有ref的allele；1 表示样品中variant的allele；2表示有第二个variant的allele
8.Vcf文件一般用什么软件生成？请至少说出两个软件。请注意不同软件生成的vcf格式的稍有不同的地方。
GATK
9.Vcf文件一般都比较大，压缩后的.gz文件用什么指令直接查看而不用解压后查看？
zcat
10.了解gvcf是什么格式，gvcf全称是什么？他与vcf有什么前后联系
Genome Variant Call Format。GVCF文件包含所有位点的记录，既包括存在变异的位点，也包括不存在变异的位点,是为了方便后续的群体分析
11.把alt列出现’,’的行提取出来
12.请将chrid、postion、ref、alt、format、样本列切割出来生成一个文本

grep -v "^##" tmp.vcf | cut -f 1,2,4,5,9- > tmp_cut.txt

13.将一个含snp,indel信息的vcf拆成一个只含snp,一个只含indel信息的2个vcf文件。可借鉴软件
14.用指令操作indel的vcf文件，提取indel长度>4的变异行数，存成一个文本.
15.用Vcftools过滤vcf文件，如maf 设置成0.05， depth设置成5-20，统计过滤前后的变异位点的总个数
16.利用vcftools提取每个样本每一个位点的变异信息和深度信息，生成一个矩阵的文件，至少含义以下信息

image.png

1/提取出变异位点上样本有纯和突变的行数
2/统计一下1号染色体上的变异总个数。
3/提取一下BRCA1基因上发生变异的行数，如果是人类wes变异结果文件
4/统计vcf文件各样本的缺失率，如果是多个样本的群体call结果。
进阶思考题：可视化vcf中变异位点质量值，横坐标是质量值，纵坐标是百分比或者该质量值时的变异位点的总个数，可以化成线图或者点图（可能需要学一些代码，还有R画图）

VCF格式文件的shell小练习