Linux查看cvf格式的文件,【生信技能树】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练习题

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

$grep '^##' tmp.vcf | grep 'INDEL'

##INFO=

$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

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

$grep '^##' tmp.vcf | grep 'depth'

##INFO=

$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

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

$awk '$4>$5{print}' tmp_INDEL.vcf > tmp_deletion.vcf

$cat tmp_deletion.vcf |wc -l

68

$awk '$4 tmp_insertion.vcf

$cat tmp_insertion.vcf | wc -l

1

统计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

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

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

26

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

筛选测序深度大于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

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

$grep -v '^#' tmp.vcf |awk '$6>30{print}' > tmp_QUALmorethan30.vcf

$cat tmp_QUALmorethan30.vcf | wc -l

99

$less -SN tmp_QUALmorethan30.vcf

组合筛选变异位点质量值大于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

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

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

【理解DP4=4,7,11,18】

参考序列变异位点之前有4个高质量测序碱基，之后有7个高质量测序碱基;

变异序列变异位点之前有11个高质量测序碱基，之后有18个高质量测序碱基

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

在前面步骤的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