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条目
统计INDEL和SNP条目的各自的平均测序深度
把INDEL条目再区分成insertion和deletion情况
统计SNP条目的突变组合分布频率
找到基因型不是1/1的条目，个数
筛选测序深度大于20的条目
筛选变异位点质量值大于30的条目
组合筛选变异位点质量值大于30并且深度大于20的条目
理解DP4=4,7,11,18这样的字段，就是 Number of high-quality ref-forward , ref-reverse, alt-forward and alt-reverse bases 计算每个变异位点的 AF
在前面步骤的bam文件里面找到这个vcf文件的某一个突变位点的测序深度表明的那些reads，并且在IGV里面可视化bam和vcf定位到该变异位点。

其它思考题

vcf的全称是什么？是用来记录什么信息的标准格式的文本？

一般选用哪个指令查看vcf文件，为什么不用vim?

vcf文件以’##’开头的是什么信息？请认真查看这些信息。’#’开头的是什么信息？

Vcf文件除头信息，每行有多少列，请详细叙述每行的含义！请准确记忆。

理解format列和样本列的对应关系以及GT AD DP的含义。

Vcf文件第三列如果不是’.’，出现的rs号的id是什么？

Vcf文件的ref，alt列和样本列的0/1 1/1 或者1/2的联系？

Vcf文件一般用什么软件生成？请至少说出两个软件。请注意不同软件生成的vcf格式的稍有不同的地方。

Vcf文件一般都比较大，压缩后的.gz文件用什么指令直接查看而不用解压后查看？

了解gvcf是什么格式，gvcf全称是什么？他与vcf有什么前后联系？

把alt列出现’,’的行提取出来

请将chrid、postion、ref、alt、format、样本列切割出来生成一个文本

将一个含snp,indel信息的vcf拆成一个只含snp,一个只含indel信息的2个vcf文件。可借鉴软件

用指令操作indel的vcf文件，提取indel长度>4的变异行数，存成一个文本。

用Vcftools过滤vcf文件，如maf 设置成0.05， depth设置成5-20，统计过滤前后的变异位点的总个数

利用vcftools提取每个样本每一个位点的变异信息和深度信息，生成一个矩阵的文件，至少含义以下信息

Eg:

chrid	postion	sample_DP	sample_GT
chr1	1010	5	AA

提取出变异位点上样本有纯和突变的行数

统计一下1号染色体上的变异总个数。

提取一下BRCA1基因上发生变异的行数，如果是人类wes变异结果文件

统计vcf文件各样本的缺失率，如果是多个样本的群体call结果。

进阶思考题：可视化vcf中变异位点质量值，横坐标是质量值，纵坐标是百分比或者该质量值时的变异位点的总个数，可以化成线图或者点图（可能需要学一些代码，还有R画图）

本文作者：生信技能树团队

最后编辑于：2019.07.23 19:56:59