生物信息学BioStat生信

学生信的那些事儿之十六 - 生信技能树VCF格式文件的shell

2019-07-18  本文已影响114人  李白爱吃橙子

开始做题前需要生成vcf文件先：

# 先生成BAM文件，用bowtie2
bowtie2 -x ../index/lambda_virus -1 reads_1.fq -2 reads_2.fq | samtools sort -@ 5 -o test2.bam

# 然后生成vcf文件，用bcftools，当然了前提的前提是安装了bcftools软件
bcftools mpileup -f ../reference/lambda_virus.fa test2.bam | bcftools call -vm > test2.vcf

想要顺利的把题目做完，基础的Linux命令当然需要熟练运用，更关键的是对 vcf 文件格式和内容的了解。

  1. 把突变记录的vcf文件区分成 INDEL和SNP条目
# 网上有几种答案，我比较倾向于这一种。如果是SNP,那么第四列和第五列都应该只有一个字符，否则就是InDel
less test2.vcf | grep -v '^#'| awk '{if(length($4)==1 && length($5)==1) print}' | less -S   #SNP
less test2.vcf | grep -v '^#'| awk '{if(length($4)!=1 && length($5)!=1) print}' | less -S

# 实际操作后我发现上面的代码也有问题，因为上面两行代码的值分别为36和29下面代码的行数是106,理论上两者之和应该等于106，然而明显不是。
less test2.vcf | grep -v '^#'| wc

# 经过对三种情况的对比后发现，第一行代码没问题，问题在于第二行，逻辑是有问题的，因为当第4列不等于1时，第5列可能等于1，且为InDel，所以正确的代码应该是：
less test2.vcf | grep -v '^#'| awk '{if(length($4)!=1) print}' | less -S   #INDEL
  1. 统计INDEL和SNP条目的各自的平均测序深度 
# 非常朴素的想法：把对应的部分拿出来然后加和求平均值
less test2.vcf | grep -v '^#'| awk '{if(length($4)!=1) print}' | cut -f 8 | cut -d ";" -f 4 | cut -d "=" -f 2 | awk '{sum+=$1}END{print sum}' #INDEL总测序深度

less test2.vcf | grep -v '^#'| awk '{if(length($4)==1 && length($5)==1) print}' | cut -f 8 | cut -d ";" -f 1 | cut -d "=" -f 2 | awk '{sum+=$1}END{print sum}'  #SNP总测序深度
  1. 把INDEL条目再区分成insertion和deletion情况
# 如果是insert，那么第5列比第4列字符数多；如果是deletion，那么第5列比第4列字符数少
less test2.vcf | grep -v '^#'| awk '{if(length($4)!=1) print}' | awk '{if(length($4)<length($5)) print}' | less -S   # insert

less test2.vcf | grep -v '^#'| awk '{if(length($4)!=1) print}' | awk '{if(length($4)>length($5)) print}' | less -S   # deletion
  1. 统计SNP条目的突变组合分布频率 
# 说实话一开始没看懂问题是什么意思...  意思就是看下SNP的突变类型有哪些，各种类型的比例是多少。

less test2.vcf | grep -v '^#'| awk '{if(length($4)==1 && length($5)==1) print}' | cut -f 4,5 | sort | uniq -c  #可以看都有哪些突变类型以及每种突变类型的数量
  1. 找到基因型不是 1/1 的条目以及个数
# 先找到条目有哪些
less test2.vcf | grep -v "1/1" | grep -v "^#"

# 计数
less test2.vcf | grep -v "1/1" | grep -v "^#" | wc
  1. 筛选测序深度大于20的条目
# 即DP>20的条目, 反其道而行之，如下：
less test2.vcf | grep -v "DP=[1-20]" | grep -v "^#" | less -S
  1. 筛选变异位点质量值大于30的条目
# 质量值在第6列
less test2.vcf | grep -v "^#" | awk '{if($6>30) print}' | less -S
  1. 组合筛选变异位点质量值大于30并且深度大于20的条目
# 筛选第6列,且满足第7题条件
less test2.vcf | grep -v "^#" | awk '{if($6>30) print}' | grep -v "DP=[1-20]" | less -S
  1. 理解DP4=4,7,11,18 这样的字段，就是 Number of high-quality ref-forward , ref-reverse, alt-forward and alt-reverse bases 计算每个变异位点的 AF
# 在vcf文件中确实可以找到DP4这个tag，但是到目前为止查询的教程都没有关于这个的解释，伤脑筋。。。
  1. 在前面步骤的bam文件里面找到这个vcf文件的某一个突变位点的测序深度表明的那些reads，并且在IGV里面可视化bam和vcf定位到该变异位点。

其他思考题

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

​ 全称是variant call format，是用来记录SNP,INDEL和SV变异信息的文本

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

​ 一般用cat, less就行；vim可以对文本进行编辑，有误操作文件的风险

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

​ 解释说明的信息

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

​ 每行固定的是8列，后面还有其他列可以添加

​ CHROM 和 POS：参考序列名和variant的位置；如果是INDEL的话，位置是INDEL的第一个碱基位置。

​ ID：variant的ID。比如在dbSNP中有该SNP的id，则会在此行给出；若没有，则用’."表示其为一个novel variant。

​ REF 和 ALT：参考序列的碱基 和 Variant的碱基。

​ QUAL：Phred格式(Phred_scaled)的质量值，表 示在该位点存在variant的可能性；该值越高，则variant的可能性越大；计算方法：

​ Phred值 = -10 * log (1-p) p为variant存在的概率; 通过计算公式可以看出值为10的表示错误概率为0.1，该位点为variant的概率为90%。

​ FILTER：使用上一个QUAL值来进行过滤的话，是不够的。GATK能使用其它的方法来进行过滤，过滤结果中通过则该值为”PASS”;

​ 若variant不可靠，则该项不为”PASS”或”.”。

​ INFO：这一行是variant的详细信息，内容很多，以下再具体详述。

​ FORMAT 和 TTG11B：这两行合起来提供了’TTG11B′这个sample的基因型的信息。’TTG11B′代表这该名称的样品，是由BAM文件中的@RG下的

​ SM 标签决定的。

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

​ GT:样品的基因型（genotype）,两个数字中间用’/"分 开，这两个数字表示双倍体的sample的基因型。0 表示样品中有ref的allele； 1 表示样品中variant的allele； 2表示有第二个variant的allele。因此： 0/0 表示sample中该位点为纯合的，和ref一致； 0/1 表示sample中该位点为杂合的，有ref和variant两个基因型； 1/1 表示sample中该位点为纯合的，和variant一致。

​ AD:AD(Allele Depth)为sample中每一种allele的reads覆盖度,在diploid中则是用逗号分割的两个值，前者对应ref基因型，后者对应variant基因型.

​ DP:Depth,为sample中该位点的覆盖度

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

​ dbSNP数据库的编号

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

​ 可部分参考第5题答案

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

​ bcftools软件、GATK软件

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

​ 用zcat命令

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

​ gVCF全称是：Genome Variant Call Format file。 gVCF和VCF的最大区别是在于gVCF文件会记录所有的点，包括哪些没有突变的点。

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

​ 是我没理解么，ALT列不全都是碱基么，怎么会有‘，’？

  1. 请将chrid、postion、ref、alt、format、样本列切割出来生成一个文本
# 相当于提取出第1，2，4，5，9，10列
less test2.vcf | grep -v "^#" | awk '{print $1,$2,$4,$5,$9,$10}' > tmp.txt
  1. 将一个含snp,indel信息的vcf拆成一个只含snp,一个只含indel信息的2个vcf文件(可借鉴软件)

​ 参考上面第1题代码，直接结果输出到一个新文件即可

  1. 用指令操作indel的vcf文件，提取indel长度>4的变异行数，存成一个文本。
# 生成indel的vcf文件先
less test2.vcf | grep -v '^#'| awk '{if(length($4)!=1) print}' > InDel.vcf

# 筛选indel长度大于4的行
cat InDel.vcf | awk '{if(length($5)>4)print}' | less -S
  1. 用Vcftools过滤vcf文件，如maf 设置成0.05， depth设置成5-20，统计过滤前后的变异位点的总个数
# 如果没有需要先下载安装vcftools软件，然后--help或者看网上教程学习相关参数：https://vcftools.github.io/man_latest.html

vcftools --vcf test2.vcf --maf 0.05 --keep-INFO DP=[5-20] --out mtest.vcf

# 运行了一下，总是报错，我还没找到解决办法暂时。。。
  1. 利用vcftools提取每个样本每一个位点的变异信息和深度信息，生成一个矩阵的文件，至少含义以下信息

​ 这部分需要学一点R语言知识，等学了再来搞它

上一篇下一篇

猜你喜欢

热点阅读