【直播】我的基因组 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来画一个呗~
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)'),counts=c(180515+181567,698322+697568,184176+185144,148387+148580,177215+177415,676816+675821))library(ggplot2)p=ggplot(dat,aes( x=type,y=counts))+geom_bar(stat="identity")print(p)
当然，mutation spectrum这个画图代码只能出一个最简单的条形图，如果你想达到下面的效果，需要学习ggplot啦！
画条形图请参考：
如果要区分染色体，可以重新考虑第1行，拿去可视化！
文：Jimmy
图文编辑：吃瓜群众