如何用vcftools从VCF文件中提取某条染色体信息

Posted 2023-04-12

参考技术A vcftools --gzvcf input.vcf --chr n --recode – recode-INFO-all --stdout | gzip -c > output.vcf.gz
说明：
–gzvcf：处理压缩格式的vcf文件（可替换为–vcf）
–chr n：选择染色体n，例：–chr 1
–recode：重新编码为vcf文件,有过滤操作都要加上--recode
–recode-INFO-all：将输出的文件保存所有INFO信息
–stdout：标准输出，后接管道命令
–gzip -c：压缩

--max-missing
--max-missing的取值是0-1，为1时表示某个位点上所有的样本必须都有基因型，一个样本的基因型都不能缺。所以这个选项可以理解为：能分型的样本占总样本的比例至少为多少。

基本的思想就是利用数据流重定向，把原来输出到屏幕上的数据定向">"到文件里

用R语言对vcf文件进行数据挖掘.11 CNV分析

参考技术A 目录

在之前的文章里介绍了如何通过直方图来可视化等位杂合碱基的比例来判断物种的染色体倍数性。在本文里会继续向下挖掘，介绍如何可视化染色体上的拷贝数变化(CNVs)。

和前文一样的操作，使用包自带的数据。

我们需要去除过高和过低深度的数据。和前文的操作一样，提取vcf文件里的深度数据"AD"。

然后过滤出10%~90%的数据，当然此处可以根据实际情况进行微调。然后对第一种出现频率最高的碱基进行可视化。(一般情况下一个位点上会有两种碱基，具体参考前文。)

同样也可以对出现频率第二高的碱基进行同样的操作，这里节约篇幅就省略了。

为了避免复杂的基于AD比例的模型假设，程序里设计了非参数估计法来计算峰值。计算完了以后可以直接对染色体进行拆分以后可视化进行校验。

根据尺寸把染色体分割成合适的大小

然后用 freq_peak 函数计算峰值。并对数据进行处理，去掉负数和Na值。

计算到此为止，可以可视化实际数据来验证计算的正确性。

仔细想一下，峰值计算的结果其实就是CNV的结果。这里根据窗口大小把染色体分成了若干段。(那么是不是可以给每一段 CDS进行细分然后计算出每一个CDS的具体数字呢????)

当然也可以把所有样本组合到一起。