挑战生信菜鸟团VCF格式文件的shell小练习

参考生信菜鸟团VCF格式文件的shell小练习

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

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

解答

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

grep -v '#' tmp.vcf |cut -f 8|cut -d ';' -f 1|awk  -F '='  '{print $1}'|sort|uniq -c

第一题

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

cd ~/project/test/bowtie2-2.3.4.3-linux-x86_64/example/reads
grep -v '#' tmp.vcf |cut -f 8|head -10|grep 'DP='

image.png

DP有的是第一列的，有的是第四列的，解决起来比较困难。需要使用vcftools，bcftools。那就要求深度学习了。学习中。。。

第一步（废话）：安装软件

source activate rna
conda install -y vcftools bcftools

第二步：收集学习资料

https://www.jianshu.com/p/badd24cbc538
http://vcftools.sourceforge.net/man_latest.html
https://www.jianshu.com/p/e05ff3cace56
https://www.jianshu.com/p/52b2dcb601d2

image.png

一个失败的方案

vcftools --vcf tmp.vcf --site-mean-depth

image.png	image.png

顺利运行，生成了一个文件，但是文件内的结果不对

第二个失败方案

vcftools --vcf tmp.vcf --depth -c > depth_summary.txt

成功方法：用vcftools

vcftools --vcf ABC.raw.vcf --remove-indels --recode --recode-INFO-all --out SNPs_only
vcftools --vcf SNPs_only.recode.vcf --depth  -c 

vcftools --vcf ABC.raw.vcf --keep-only-indels --recode --recode-INFO-all --out indel_only
vcftools --vcf indel_only.recode.vcf --depth  -c 

image.png