开始做题前需要生成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 文件格式和内容的了解。
- 把突变记录的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
- 统计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总测序深度
- 把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
- 统计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的条目以及个数
# 先找到条目有哪些
less test2.vcf | grep -v "1/1" | grep -v "^#"
# 计数
less test2.vcf | grep -v "1/1" | grep -v "^#" | wc
- 筛选测序深度大于20的条目
# 即DP>20的条目, 反其道而行之，如下：
less test2.vcf | grep -v "DP=[1-20]" | grep -v "^#" | less -S
- 筛选变异位点质量值大于30的条目
# 质量值在第6列
less test2.vcf | grep -v "^#" | awk '{if($6>30) print}' | less -S
- 组合筛选变异位点质量值大于30并且深度大于20的条目
# 筛选第6列,且满足第7题条件
less test2.vcf | grep -v "^#" | awk '{if($6>30) print}' | grep -v "DP=[1-20]" | less -S
- 理解
DP4=4,7,11,18这样的字段，就是 Number of high-quality ref-forward , ref-reverse, alt-forward and alt-reverse bases 计算每个变异位点的 AF
# 在vcf文件中确实可以找到DP4这个tag，但是到目前为止查询的教程都没有关于这个的解释，伤脑筋。。。
- 在前面步骤的bam文件里面找到这个vcf文件的某一个突变位点的测序深度表明的那些reads，并且在IGV里面可视化bam和vcf定位到该变异位点。
其他思考题
- vcf的全称是什么？是用来记录什么信息的标准格式的文本？
全称是variant call format，是用来记录SNP,INDEL和SV变异信息的文本
- 一般选用哪个指令查看vcf文件，为什么不用vim?
一般用cat, less就行；vim可以对文本进行编辑，有误操作文件的风险
- vcf文件以’##’开头的是什么信息？请认真查看这些信息。’#’开头的是什么信息？
解释说明的信息
- 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 标签决定的。
- 理解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中该位点的覆盖度
- Vcf文件第三列如果不是’.’，出现的rs号的id是什么？
dbSNP数据库的编号
- Vcf文件的ref，alt列和样本列的0/1 1/1 或者1/2的联系？
可部分参考第5题答案
- Vcf文件一般用什么软件生成？请至少说出两个软件。请注意不同软件生成的vcf格式的稍有不同的地方。
bcftools软件、GATK软件
- Vcf文件一般都比较大，压缩后的.gz文件用什么指令直接查看而不用解压后查看？
用zcat命令
- 了解gvcf是什么格式，gvcf全称是什么？他与vcf有什么前后联系？
gVCF全称是：Genome Variant Call Format file。 gVCF和VCF的最大区别是在于gVCF文件会记录所有的点，包括哪些没有突变的点。
- 把alt列出现’,’的行提取出来
是我没理解么，ALT列不全都是碱基么，怎么会有‘，’？
- 请将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
- 将一个含snp,indel信息的vcf拆成一个只含snp,一个只含indel信息的2个vcf文件(可借鉴软件)
参考上面第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
- 用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
# 运行了一下，总是报错，我还没找到解决办法暂时。。。
- 利用vcftools提取每个样本每一个位点的变异信息和深度信息，生成一个矩阵的文件，至少含义以下信息
这部分需要学一点R语言知识，等学了再来搞它
网友评论