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

Python024

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

目录

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

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

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

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

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

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

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

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

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

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

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

目录

vcfR 可以直接读取vcf格式的数据。如果同时读取参照序列fasta格式的序列文件和gff格式文件的注释文件还可以获取更完整的信息(此步骤并非必须,可以只读取vcf数据)。在此处便于重复用到了 pinfsc50 包。这个包里是植物致病微生物的基因序列测序结果。包含了一个vcf文件,一个fasta文件和一个gff文件。

这里用到参照序列的数据。

当这些数据被读取到内存的时候就可以开始对染色体名字或者其它一些东西进行修改了。由于 vcfR 更擅长对的单独染色体进行分析,所以当你的基因过大或者有很多样本的时候,建议对数据进行拆分。

读取完数据以后就可以建立 chromR ,来对数据进行详细的分析。

首先对数据进行初步的可视化,

我们在上面的图里得到很多信息,比方说测序深度(DP)的峰在500,但是拖着尾巴,这个尾巴表示数据里包含着CNV信息。然后比对质量(MQ)的峰值在60,于是我们可以以60为中心对数据进行过滤。

使用 masker 可以对数据进行过滤标记。然后可视化过滤以后的数据。

是不是顺眼多了。当然我们也可以看一下SNP的分布情况。注意右下角的图。

用 chromoqc() 可以对数据进行更完整的可视化。包括外显子内含子的分布,GC含量的分布等等。

最后可以用函数 write.vcf() 把数据输出成新的vcf文件。

介绍使用tidytext进行文本挖掘。

整洁的数据应该是这样的

对于整洁的文本数据,储存在每行中的数据通常是单个单词,但也可以是n-gram,句子或段落。

使用unnest_tokens函数对数据进行处理

简单介绍一下unnest_tokens函数:

unnest_tokens这里使用的两个基本参数。首先,输出的列名,上面是word,然后是文本来输入列(text在本例中)。

使用之后unnest_tokens,我们将每行拆分

文本分析的流程:

就是写了傲慢与偏见的那个人,说实话这部作品的确值得一看

数据来自于Jane Austen的 janeaustenr 包

linenumber 对应的是多少行, chapter 对应的是第多少章。

要将其作为一个整洁的数据集来处理,还需要将句子转化成文更加基本的格式

此函数使用 tokenizers 包将原始数据框中的每一行文本分隔为标记。默认标记化用于单词,但其他选项包括字符,n-gram,句子,行,段落或正则表达式模式周围的分隔。

也就是修改下面这个参数:

既然数据是每行一个字的格式,我们可以使用像dplyr这样的整洁工具来操作它。通常在文本分析中,我们会想要删除停用词停用词是对分析无用的词,通常是非常常见的词,例如英语中的“the”,“of”,“to”等等。我们可以用一个删除停用词(保存在tidytext数据集中stop_words)anti_join()。

我们也可以使用 dplyr count() 来查找所有书籍中最常见的单词。

可以看见,最常见的单词是 miss

进行可视化:

因为我们一直在使用整洁的工具,所以我们的字数存储在一个整洁的数据框中。这允许我们将它直接传递给ggplot2包,例如创建最常见单词的可视化