生信软件31 - Bcftools操作VCF/BCF文件高级用法合集

生信与基因组学

已于 2024-08-24 12:28:32 修改

阅读量1.8k

点赞数 41

分类专栏：生物信息学软件文章标签： linux 数据分析数据挖掘

于 2024-08-24 03:00:00 首次发布

本文链接：https://blog.csdn.net/littlecomputerrobot/article/details/141474654

版权

生物信息学软件专栏收录该内容

36 篇文章 7 订阅

订阅专栏

1.Bcftools 从VCF/BCF文件中提取信息

bcftools query命令可用于提取任何VCF字段。

# 查看vcf文件包含样本名称
bcftools query -l sample.vcf

# 查看vcf文件包含样本数量
bcftools query -l sample.vcf| wc -l

# 打印POS列信息, head显示前10列
bcftools query -f '%POS\n' sample.vcf|head

# 打印CHROM POS  REF  ALT 4列信息
bcftools query -f '%CHROM %POS %REF %ALT\n' sample.vcf|head
# chr1 69270 A G

2. Bcftools 提取每个位点的等位基因频率

# 提取每个位点的等位基因频率
bcftools query -f '%CHROM %POS %AF\n' sample.vcf |head 
# chr1 69270 1


# 如果AF注释不存在，但AN和AC存在，可计算频率
bcftools query -f '%CHROM %POS %AN %AC{0}\n' sample.vcf | \
		awk '{printf "%s %s %f\n",$1,$2,$4/$3}' |head 
# 1 13380 0.000077

3. Bcftools 提取每个样本标签

# 提取CHROM  POS  基因型和PL列信息
bcftools query -f '%CHROM %POS[\t%GT\t%PL]\n' sample.vcf |head

提取每个样本标签

4. Bcftools 按固定列过滤

固定列如QUAL、FILTER、INFO可以直接过滤。

# 使用-e 'FILTER="."'表达式排除不是FILTER的位点
bcftools query -e 'FILTER="."' \
			-f '%CHROM %POS %FILTER\n' sample.vcf |head
# 1 3000150 PASS
# 1 3000151 LowQual

4. Bcftools 按质量和深度过滤

bcftools query -i 'QUAL>20 && DP>10' \
				-f '%CHROM %POS %QUAL %DP\n' sample.vcf | head
# 1 14930 31.2757 13
# 1 17538 37.9458 12

5. Bcftools mpileup生成可能基因型文件

BAM文件作为输出，输出可能基因型的BCF文件。

bcftools mpileup \
		--max-depth 10000 \
		--threads 4 \
		-f reference.fasta \
		-o genotype.bcf \
		sample.bam

参数说明：

--max-depth或-d：为比对中的每个位置设置每个输入文件的读数。在本例中，设置为10000;
--threads：设置要使用的处理器/线程的数量（n）;
--fasta-ref或-f：用于选择用于比对的faidx索引的FASTA核苷酸参考文件（reference.fasta）;
--output 或-o：输出文件（genotype.bcf）;
最后一个参数：输入BAM比对文件（sample.bam）, 可以提供多个输入文件。

6. Bcftools call生成SNP/indel变异文件

bcftools call可用于从BCF文件call SNP/indel变异。

bcftools call \
		-O v \
		--threads 4 \
		-vc --ploidy 2 \
		-p 0.05 \
		-o variants.unfiltered.vcf \
		genotype.bcf

参数说明：

--output-type或-O：用于选择输出格式。在本例中，v表示VCF；b表示生成BCF。
--threads：设置要使用的处理器/线程的数量（n）。
-vc：指定输出仅包含变体。
--ploidy：指定程序的倍性（人类为二倍体）。
--pval-threshold或-p：用于设置变异位点的p值阈值（0.05）。
--output 或-o：输出文件（variants.unfiltered.bcf）。
最后一个参数：输入BCF文件（genotype_likelihoods.bcf）。

7. Bcftools filter变异文件

bcftools filter可用于从BCF文件中过滤变异。

 bcftools filter \
		 --threads n \
		 -i '%QUAL>=20' \
		 -O v \
		 -o variants.filtered.vcf  \
		 variants.unfiltered.vcf

参数说明：

--threads：设置要使用的处理器/线程的数量（n）。
--include或-i：过滤表达式， %QUAL>=20表示筛选出质量大于或等于20位点。
--output-type或-O：输出格式，，v代表VCF。
--output 或-o：输出文件（variants.filtered.vcf）。
最后一个参数：输入BCF文件（variants.unfiltered.vcf）。

8. 检测非整倍性和污染

检测影响整个染色体的SV，如非整倍性或污染，可以直接从BAF值的总体分布中识别。

bcftools polysomy -v -o outdir/ sample.vcf

# 对于二倍体，获取第三列值缺失1.0和重复3.0的行
cat outdir/dist.dat | awk '$1=="CN" && $3!=2.0'

# matplotlib绘图
python outdir/dist.py

下图输出的结果图中，有三个X染色体拷贝。

9. 23andMe转换为VCF

转换只考虑SNP而忽略缺失和插入。

# 原始的23andMe结果示例
# 四列分别为：标记ID，染色体名称，位置和基因型
#rs6139074       20      63244   AA
#rs1418258       20      63799   CC
#rs6086616       20      68749   TT
#rs6039403       20      69094   AG

bcftools convert \
		--tsv2vcf sample.23andMe.gz \
		-f reference.fa \
		-s sample_name \
		-Ov \
		-o sample.vcf
		
# Rows total:     612647
# Rows skipped:   4751 （跳过基因型）
# Missing GTs:    20525 (缺失基因型)
# Hom RR:     318339
# Het RA:     165598
# Hom AA:     103420
# Het AA:     14
# Het AA数量应该比较小，因为大量的杂合等位基因（Het AA）表明输入等位基因不在正链上。

10. 合并到多个样本的VCF文件

bcftools index sample1.vcf
bcftools index sample2.vcf

# 合并为merge.vcf 
bcftools merge -Ov -o merge.vcf sample1.vcf sample2.vcf

11. VCF文件添加注释信息

输入VCF/BCF文件、fasta格式的参考基因组和可从Ensembl网站下载的GFF 3文件，输出包含注释信息的VCF文件，目前，仅支持Ensembl GFF 3文件。

bcftools csq \
		-f reference.fa \
		-g Homo_sapiens.GRCh37.82.gff3.gz \
		sample.vcf \
		-Ov \
		-o output.vcf

bcftools与三款主流变异注释软件速度和内存消耗比较。
软件比较

生信软件文章推荐

生信软件1 - 测序下机文件比对结果可视化工具 visNano

生信软件2 - 下游比对数据的统计工具 picard

生信软件3 - mapping比对bam文件质量评估工具 qualimap

生信软件4 - 拷贝数变异CNV分析软件 WisecondorX

生信软件5 - RIdeogram包绘制染色体密度图

生信软件6 - bcftools查找指定区域的变异位点信息

生信软件7 - 多线程并行运行Linux效率工具Parallel

生信软件8 - bedtools进行窗口划分、窗口GC含量、窗口测序深度和窗口SNP统计

生信软件9 - 多公共数据库数据下载软件Kingfisher

生信软件10 - DNA/RNA/蛋白多序列比对图R包ggmsa

生信软件11 - 基于ACMG的CNV注释工具ClassifyCNV

生信软件12 - 基于Symbol和ENTREZID查询基因注释的R包（easyConvert ）

生信软件13 - 基于sambamba 窗口reads计数和平均覆盖度统计

生信软件14 - bcftools提取和注释VCF文件关键信息

生信软件15 - 生信NGS数据分析强大的工具集ngs-bits

生信软件16 - 常规探针设计软件mrbait

生信软件17 - 基于fasta文件的捕获探针设计工具catch

生信软件18 - 基于docker部署Web版 Visual Studio Code

生信软件19 - vcftools高级用法技巧合辑

生信软件20 - seqkit+awk+sed+grep高级用法技巧合辑

生信软件21 - 多线程拆分NCBI-SRA文件工具pfastq-dump

生信软件22 - 测序数据5‘和3‘端reads修剪工具sickle

生信软件23 - Samtools和GATK去除PCR重复方法汇总

生信软件24 - 查询物种分类学信息和下载基因组TaxonKit和ncbi-genome-download

生信软件25 - 三代测序数据灵敏比对工具ngmlr

生信软件26 - BWA-MEM比对算法性能更好的bwa-mem2

生信软件27 - 基于python的基因注释数据查询/检索库mygene

生信软件28 - fastq与bam的reads数量计算与双端fastq配对检测工具fastq-pair

生信软件29 - 三代数据高效映射精确的长读段比对工具mapquik

生信软件30 - 快速单倍型分析工具merlin

生信软件31 - Bcftools操作VCF/BCF文件高级用法合集

1.Bcftools 从VCF/BCF文件中提取信息

2. Bcftools 提取每个位点的等位基因频率

3. Bcftools 提取每个样本标签

4. Bcftools 按固定列过滤

4. Bcftools 按质量和深度过滤

5. Bcftools mpileup生成可能基因型文件

参数说明：

6. Bcftools call生成SNP/indel变异文件

参数说明：

7. Bcftools filter变异文件

参数说明：

8. 检测非整倍性和污染

9. 23andMe转换为VCF

10. 合并到多个样本的VCF文件

11. VCF文件添加注释信息

生信软件文章推荐

更多内容请关注公众号【生信与基因组学】，定期更新生信算法和编程、基因组学、统计学、分子生物学、临床检测和深度学习等内容。