生物数据格式 - vcf/bcf
2021-02-22 本文已影响0人
半夜一更
格式
VCF是Variant Call Format的简称，该格式文件是专门用于存储基因序列突变信息的一种文本文件，包括单碱基突变SNV、单核苷酸多态性SNP、InDel、拷贝数变异CNV和结构变异SV等，文件可以采取editplus或pilotedit(建议)打开查看，其二进制存储格式是BCF。vcf文件后续可以用于多种分析，包括但不限于：进化树分析、群体结构分析、PCA分析、GWAS关联分析等。vcf文件格式如下：
##fileformat=VCFv4.0
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=1000GenomesPilot-NCBI36
##phasing=partial
##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples With Data">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##INFO=<ID=AF,Number=.,Type=Float,Description="Allele Frequency">
##INFO=<ID=AA,Number=1,Type=String,Description="Ancestral Allele">
##INFO=<ID=DB,Number=0,Type=Flag,Description="dbSNP membership, build 129">
##INFO=<ID=H2,Number=0,Type=Flag,Description="HapMap2 membership">
##FILTER=<ID=q10,Description="Quality below 10">
##FILTER=<ID=s50,Description="Less than 50% of samples have data">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
##FORMAT=<ID=HQ,Number=2,Type=Integer,Description="Haplotype Quality">
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA00001 NA00002 NA00003
20 14370 rs6054257 G A 29 PASS NS=3;DP=14;AF=0.5;DB;H2 GT:GQ:DP:HQ 0|0:48:1:51,51 1|0:48:8:51,51 1/1:43:5:.,.
20 17330 . T A 3 q10 NS=3;DP=11;AF=0.017 GT:GQ:DP:HQ 0|0:49:3:58,50 0|1:3:5:65,3 0/0:41:3
20 1110696 rs6040355 A G,T 67 PASS NS=2;DP=10;AF=0.333,0.667;AA=T;DB GT:GQ:DP:HQ 1|2:21:6:23,27 2|1:2:0:18,2 2/2:35:4
20 1230237 . T . 47 PASS NS=3;DP=13;AA=T GT:GQ:DP:HQ 0|0:54:7:56,60 0|0:48:4:51,51 0/0:61:2
20 1234567 microsat1 GTCT G,GTACT 50 PASS NS=3;DP=9;AA=G GT:GQ:DP 0/1:35:4 0/2:17:2 1/1:40:3
VCF文件开头是整体的注释信息，以##作为起始，其后接FILTER、INFO、FORMAT等，其中##FILTER开头的行是vcf主体record中第七列缩写词的说明、##INFO开头的行注释vcf主体record中第8列的缩写字母、##FORMAT开头的行注释第九列的缩写字母。
在header之后，vcf主体的每一行是一条record，固定列有9列，以及在之后的第十列，它们分别是：
第一列：#CHROM，染色体号
第二列：POS，在染色体上的位置
第三列：ID，突变名称，一般只有人类才有dbSNP编号，以rs开头
第四列：REF，参考基因组碱基类型，必须大写
第五列：ALT，变异碱基类型，大写，多个以逗号分隔，‘.'表示缺失
第六列：QUAL，变异检测质量值，越高越可靠
第七列：FILTER，标记过滤结果的列：通过质控过滤标准的标记为‘PASS’，后续可用其他工具进行挑选过滤
第八列：INFO，附加信息列，附加信息的注释在header的##INFO中
第九列：FORMAT，后面信息的说明列
第十列开始为样品信息：GT=genotype、AD=碱基支持数量、DP=测序深度总和、PL=归一化后基因型的可能性、GQ=PL判读的基因型的质量值，其中当第二小的值小于99时，有必要怀疑基因型的可靠性。
获取
vcf文件基本由bam文件生成，当得到排序并建立索引的bam文件后，可以使用多种工具例如bcftools、gatk、freebayes、lumpy、delly、varscan2等处理得到。
#bcftools进行变异检测
bcftools mpileup -f ref.fna A.sorted.bam | bcftools call -c -v A.bcftools.vcf
#freebayes进行snp检测
freebayes -f ref.fna -C 5 A.sorted.bam -v A.freebayes.vcf
#GATK进行snp检测
gatk HaplotypeCaller --emit-ref-confidence GVCF -R ref.fna -I A.sorted.bam -O A.gatk.vcf.gz
gatk GenotypeGVCFs -R ref.fna -I A.gatk.vcf.gz -O A.HC.vcf.gz
#delly进行SV检测
delly call -g ref.fna -o A.delly.sv.bcf -n A.sorted.bam
delly filter -f germline -p -q 20 A.delly.sv.bcf -o A.delly.sv.filter.bcf
处理
处理vcf格式文件的软件有许多种，包括：bcftools、vcftools、gatk、python_pyvcf、plink等。
文件与内容处理
##查看
bcftools view A1.bcf.gz
##格式转换
bcftools view A1.vcf -O b -o A1.bcf.gz #vcf转换bcf
bcftools view A1.vcf -O v -o A1.bcf.gz #bcf转换vcf
##建立索引
bcftools index A1.bcf.gz #生成csi格式索引,默认
bcftools index -t A1.bcf.gz #生成tsi格式索引,-t
##查看指定区域
bcftools view A1.bcf.gz 20:1000-20000 #查看20号染色体1000-20000之间的突变位点
bcftools view A1.bcf.gz -R region.bed #利用本地文件查看
##合并
bcftools merge A1.vcf.gz B1.vcf.gz -O b -o merge.bcf.gz #合并前需要对每个样品创建索引
##统计
bcftools stats A1.bcf.gz >view.stats #统计突变种类，以及每种突变具体的数据
plot-vcfstats view.stats -p output #绘图
##筛选
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT[\t%SAMPLE=%GT]\n' A1.bcf.gz
##拆分
bcftools view -v snps A1.freebayes.bcf.gz | less -S #筛选SNPs
bcftools view -v indels A1.freebayes.bcf.gz | less -S #筛选InDels
##过滤 可以进行多种模式的过滤
bcftools filter -e "INFO/AF[0] > 0.3 & FORMAT/DP < 30" A1.bcf.gz #过滤掉等位频率>0.3同时覆盖深度小于10的突变位点
文件功能处理
##注释
# bcftools注释
bcftools annotate -a /ifs1/Database/GATK/hg38/dbsnp_146.hg38.vcf.gz -c ID,QUAL,+TAG file.vcf -o annotate.vcf
# snpeff注释
java -jar snpEff.jar databases | less #列出所有数据库
java -jar snpEff.jar databases |grep "Homo" #筛选人基因组数据库
java -jar SnpSift.jar annotate dbsnp_146.hg38.vcf.gz A1.bcf.gz >A1.anno.rs.vcf #与dbsnp进行注释得到rs号
# annovar注释
convert2annovar.pl -format vcf4old A1.snps.indel.VQSR.vcf.gz >A1.annovar.input #生成annovar格式
annotate_variation.pl --geneanno -buildver hg38 --outfile A1.geneanno.anno A1.annovar.input humandb/ #gene-based注释
# clinvar注释
convert2annovar.pl -format vcf4old A1.HC.snps.indel.VQSR.vcf.gz >A1.annovar.input
annotate_variation.pl --filter -buildver hg38 --outfile A1.clinvar.anno A1.annovar.input -dbtype clinvar_20180603 humandb/
##一致性序列 是指一条仅将突变位点进行替换其余与参考序列一致的序列，一致性序列完全根据参考序列为模板生成，是并不存在的序列，主要用于后面构建系统发育树。
bcftools consensus -f ref.fna -s Sample1 -o Sample1_consensus.fa file.vcf.gz
##VQSR
### Variant Quality Score Recalibration，是GATK的核心功能，利用机器学习算法对vcf进行过滤。利用机器学习算法对突变位点进行过滤比采用“一刀切”对所有位点处理的方式准确性更高。
———以上属个人理解与记录
