【原】【直播】我的基因组 45：SNV突变(6种)频谱的制作

受热心读者的委托，特意为他讲解一下SNV突变(6种)频谱的制作，同时欢迎大家留言其它需求！

突变频谱呢，就是对含有SNV的VCF格式的文件进行一个统计。

全基因组SNP突变可以分成6类（C>A, C>G, C>T, A>C, A>G, A>T）。肯定会有人问为什么是六类？

以A:T>C:G为例，此种类型SNP突变包括A>C和T>G。由于测序数据既可比对到参考基因组的正链，也可比对到参考基因组的负链，当T>C类型突变出现在参考基因组正链上，A>G类型突变即在参考基因组负链的相同位置，所以将T>C和A>G划分成一类,换句话说我们只考虑正链的突变形式，参考碱基只允许有C或者T，因为它们等价于G或者A。所以全基因组SNP突变可以分成这6类。

很明显，我们只需要考虑VCF文件的第4,5行即可！

 cat realign.vcf |grep -v INDEL |grep -v "^#" |cut -f 1-5 |head




cat realign.vcf |grep -v INDEL |grep -v "^#" |cut -f 4,5|sort |uniq -c |grep -v ","

我们过滤掉了多种变异形式的SNV，比如T，突变成G或者C！最后的结果如下：

一般来说，是要可视化一下的，我用R语言的ggplot来画一个呗~

1. dat <- data.frame(type=c('C>A(G>T)','C>T(G>A)','C>G(G>C)','T>A(A>T)','T>G(A>C)','T>C(A>G)'),

2. counts=c(180515+181567,698322+697568,184176+185144,148387+148580,177215+177415,676816+675821)

3. )

4. library(ggplot2)

5. p=ggplot(dat,aes( x=type,y=counts))+geom_bar(stat="identity")

6. print(p)

当然，mutation spectrum这个画图代码只能出一个最简单的条形图，如果你想达到下面的效果，需要学习ggplot啦！

画条形图请参考：http://docs./0.9.3.1/geom_bar.html 

如果要区分染色体，可以重新考虑第1行，拿去可视化！

http://www./1619.html

文：Jimmy

图文编辑：吃瓜群众