用R语言对vcf文件进行数据挖掘.7 测序深度覆盖度

Python08

用R语言对vcf文件进行数据挖掘.7 测序深度覆盖度,第1张

目录

vcf数据里除了位点的ATGC的对比,进行纯合/杂合判断的以外。还有一个重要的项目就是 DP ,测序深度。测序深度不仅是看测序质量的重要参考,也是对染色体倍数体以及基因拷贝数进行评估的重要指标。

一般的VCF文件都很大,用手动提取里面的信息肯定不大现实。用 vcfR 就可以轻松实现。

查看一下R读取的数据。

选取我们需要的部分也就是Genotype Section里的 DP 区域。

众所周知箱状图的特点就是(boxplot)包含了所有的信息,包括异常值outlier。正因为这个原因,这张图很大程度上受到了这些异常值的影响,变得非常难懂。自己看看还可以,用来发表文章的话肯定不行。

经过log2转换,我们可以得到理想的效果。

又或者不需要转换,而是通过过滤数据来改善箱图效果。举个例子,提取90%的信赖区间的数据来可视化。

这样也可以获得类似的结果。

目录

一般的VCF文件都很大,用手动提取里面的信息肯定不大现实。用 vcfR 就可以轻松实现。

vcfR 自带测试文件 vcfR_test 。就用这个文件来操作一下吧。

在分区 Genotype 里,通过观察 FORMAT 列可以看到一共有四种类型的数据 GT:GQ:DP:HQ ,至于这四种类型的数据个各自代表什么意思大家可以查阅知乎百度谷歌。我们可以提取出我们想要的数据类型。比方说最重要的 GT (genotype)。

同样,我们也可以提取例如 DP (测序深度Read Depth)的数字矩阵。

值的注意的是这里用到了参数 as.numeric = TRUE 使得数据自动转换成了数字。但是并不是对所有类型的数据都有效,比方说我们重复一下提取 gt 。

在没有任何报错的情况下 gt 变成了一堆毫无意义的数字,很明显不合理,不要用这些经过错误转换的数据进行下一步分析,比方说喜闻乐见的主成分分析。

在一些类型的数据里可能会出现一个以上的结果,比方说上面的 HQ 数据。

一般情况下我们只需要每一列的第一个数字

不需要samtools之类的软件我们也可以实现vcf数据读取自由,关键是可以直接写入内存进行下一步的统计分析和数据可视化,个人感觉是很有效的提高了生产力。值得花时间学习一下这个工具。