【生信技能树】VCF格式文件的shell小练习

首先友情宣传生信技能树

• 全国巡讲:R基础，Linux基础和RNA-seq实战演练 : 预告：12月28-30长沙站
• 广州珠江新城GEO数据挖掘滚动开班

---

题目来源：【生信技能树】VCF格式文件的shell小练习

首先使用bowtie2软件自带的测试数据生成sam/bam文件，还有vcf文件。代码如下：

mkdir -p ~/biosoft
cd ~/biosoft
wget https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh 
bash Miniconda3-latest-Linux-x86_64.sh 
source ~/.bashrc 
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes
conda create -y -n test
conda activate test
conda install -y samtools bcftools

wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.3.4.3/bowtie2-2.3.4.3-linux-x86_64.zip 
unzip bowtie2-2.3.4.3-linux-x86_64.zip 
cd ~/biosoft/bowtie2-2.3.4.3-linux-x86_64/example/reads

../../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

LINUX练习题

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

$grep '^##' tmp.vcf | grep 'INDEL'
##INFO=<ID=INDEL,Number=0,Type=Flag,Description="Indicates that the variant is an INDEL.">
$grep -v '^#' tmp.vcf | grep 'INDEL' > tmp_INDEL.vcf
$cat tmp_INDEL.vcf | wc -l
69
$grep -v '^#' tmp.vcf | grep -v 'INDEL' > tmp_SNP.vcf
$cat tmp_SNP.vcf | wc -l
36

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

$grep '^##' tmp.vcf | grep 'depth'
##INFO=<ID=DP,Number=1,Type=Integer,Description="Raw read depth">
$egrep -o 'DP=[0-9]+;' tmp_INDEL.vcf | awk -v FS='=' 'BEGIN{num=0;}{num=num+$2}END{print(num/NR)}'
40.6087
$egrep -o 'DP=[0-9]+;' tmp_SNP.vcf | awk -v FS='=' 'BEGIN{num=0;}{num=num+$2}END{print(num/NR)}'
37

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

$awk '$4>$5{print}' tmp_INDEL.vcf > tmp_deletion.vcf
$cat tmp_deletion.vcf |wc -l
68
$awk '$4<$5{print}' tmp_INDEL.vcf > tmp_insertion.vcf
$cat tmp_insertion.vcf | wc -l
1

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

$awk '{print $4"->"$5}' tmp_SNP.vcf | sort | uniq -c
      7 A->C
      1 A->G
      4 A->T
      2 C->A
      4 C->G
      3 G->A
      2 G->C
      4 G->T
      6 T->A
      1 T->C
      2 T->G

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

$grep -v '^#' tmp.vcf | grep -v '1/1' |wc -l
26
$grep -v '^#' tmp.vcf | grep -v '1/1' |less -SN

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

$egrep -o "DP=[0-9]+;" tmp.vcf |awk '{print(length($0)-4)}'| sort | uniq -c#查看测序深度位数
      2 1
    103 2
#测序深度最多只有2位数
$egrep "DP=2[1-9];|DP=[3-9][0-9];" tmp.vcf |wc -l
100
$egrep "DP=2[1-9];|DP=[3-9][0-9];" tmp.vcf |less -SN

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

$grep -v '^#' tmp.vcf |awk '$6>30{print}' > tmp_QUALmorethan30.vcf
$cat tmp_QUALmorethan30.vcf | wc -l
99
$less -SN tmp_QUALmorethan30.vcf

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

$egrep "DP=2[1-9]|DP=[3-9][0-9]" tmp_QUALmorethan30.vcf | wc -l
97
$egrep "DP=2[1-9]|DP=[3-9][0-9]" tmp_QUALmorethan30.vcf | less -SN

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

1. DP4，高质量测序碱基，4个数字分别对应ref-前，ref-后，alt-前，alt-后

• 【理解DP4=4,7,11,18】
  • 参考序列变异位点之前有4个高质量测序碱基，之后有7个高质量测序碱基;
  • 变异序列变异位点之前有11个高质量测序碱基，之后有18个高质量测序碱基

2. AF：Allel Frequency 等位基因频率
   参考：简书：刘小泽 - VCF格式

AC、AF、AN【和等位基因有关】： AC：Allele Count该位点变异的等位基因数目； AF：Allel Frequency 等位基因频率； AN：Allel Number 等位基因的总数目
【单看这个不好理解，举一个二倍体diploid例子：基因型0/1表示为杂合子，该位点只有一个等位基因发生突变，AF=0.5（在该位点只有50%的等位基因发生突变），总的等位基因数目为2；基因型1/1表示为纯合子，AC=2，AF=1，AN=2】

所以，
0/0的AF=0；
0/1的AF=0.5；
1/1的AF=1；
另外，AF标签属于INFO列

$grep -v "^#" tmp.vcf | awk '{print$10}' |awk -F":" '{print$1}' | sort | uniq -c #查询样本基因型
     26 0/1
     79 1/1
$grep -v "^#" tmp.vcf | awk '{if ($0~"1/1"){$8=$8";AF=1";print$0} else if ($0~"0/1"){$8=$8";AF=0.5";print$0}}' > tmp_addAF.vcf

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

$grep -v "^#" tmp.vcf |head -10 |less -SN

10.1 如图，选择1104突变位点，DP=43，参考序列为C，突变成A

10.2 检索bam文件
参考：徐洲更hoptop - SAMtools: SAM格式的处理利器
bam文件：
第4列：比对到的位置
第10列：segment序列

$samtools-1.9 view tmp.bam | awk '{if ($3!="*" && $4<=1104 && substr($10,1104-$4+1,1)=="A") print}' | wc -l
45#与vcf文件中DP=43不符合

10.3 tmp.vcf在IGV中可视化

IGV-vcf

10.4 tmp.bam在IGV中可视化

IGV-bam