对突变频谱（Mutation spectrum）进行分析，可以得知样本各种类型突变（如C:G>T:A）的数量及样本是否具有某种类型突变的偏好性。

单个碱基的替换一共有六种变异类型：C>A/G>T，C>G/G>C，C>T/G>A，T>A/A>T，T>C/A>G，T>G/A>C，有可以根据发生替换的碱基类别分为两大类：颠换transversion，嘌呤与嘧啶之间的替换，转换transition，嘌呤与嘌呤或嘧啶与嘧啶之间的替换。根据各类型点突变的比例，对样本和点突变类型进行聚类分析，可以研究样本点突变类型的偏好和各样本在点突变水平上的相似程度。

1，突变频谱数据准备

从VCF或者Annovar注释之后文件中，提取突变类型，第一列为样本名，第二列为突变类型，命名为type.txt。有来自来个实验的10个样本，分别为S1–S7和CHG。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

Name    Type
S5      T>C/A>G
S5      T>C/A>G
S5      C>T/G>A
S5      T>A/A>T
S5      T>C/A>G
S1      C>T/G>A
S1      T>C/A>G
S1      C>T/G>A
S1      T>C/A>G

2，生成突变频谱图

1
2
3
4
5
6

library(ggplot2)
data1<- read.table("/path/to/type.txt",sep='t',header=T,check=F)
f1<- ggplot(data = data1, aes(x = Name, fill = Type)) + geom_bar(position = "fill") + labs(title = "Mutation Spectrum",x = "",y = "Fraction of Mutations") + theme(panel.background = element_blank(), axis.text.x  = element_text(angle=90), text = element_text(size=16) ) 
svg(file="mutation_spectrum.svg")
f1
dev.off()

Fastq文件为纯文本文件，占用的硬盘空间较大，所以一般都会将Fastq文件压缩成gz格式，很多软件也支持fastq的gz格式输入。我用过python读取gzip，非常方便。单纯的通过gzip的命令压缩fastq，效率非常非常慢，据说是没有利用整个机器的cpu。

于是我就找到了pigz这款神器，可以在压缩数据时，发挥多核多处理器的优势，简而言之就是利用多线程进行gzip任务，比单纯的gzip压缩要快很多，有人测试快了5倍多（因为gzip压缩100G的文件时间是太长了，我也就没有测试）。

pigz

官网 http://www.zlib.net/pigz/

pigz, which stands for parallel implementation of gzip, is a fully functional replacement for gzip that exploits multiple processors and multiple cores to the hilt when compressing data.

安装pigz

1
2
3
4
5
6
7

wget http://zlib.net/pigz/pigz-2.3.3.tar.gz
tar -xvzf pigz-2.3.3.tar.gz
#如果提示不是gz格式，请尝试  tar -xvf pigz-2.3.3.tar.gz
cd pigz-2.3.3.tar.gz
make
如果报错 pigz.c:(.text.startup+0xca): undefined reference to `deflateEnd' gcc
请在第八行$(CC) $(LDFLAGS) -o pigz $^ -lpthread -lm 后面添加-lz选项，表示link libz

运行pigz

pigz -h 可以看到它的command option。

运行pigz和简单，和gzip的命令差不多，比如

# -c 表示打印到标准输出std，如果没有-c选项，则会生成一个后缀为gz的压缩文件。
pigz -c file > file.gz
# -k 表示压缩后不删除源文件
pigz -k file

不光是生信，感觉整个生物领域越来越靠图吃饭了。谁的图漂亮，谁的分析就好。好吧，我也承认，图多图漂亮了，确实显得高大上哈。在我眼中图的信息量比较大，能够给人直观的表示，为了说明图的意义，我没把持住，在本文多放了几张和circos图无关的图，见谅。本文有多图，都是本人自己画的，请勿盗图，如果你想知道怎么做或者由更好的表示办法，欢迎留言讨论。

想当年，我为了表示SNP在染色体上的数目分布，用python画了24个图，每个染色体一张SNP的密度分布。那可是我第一次用matplot。

先上我年yi轻qian时做的图，那个时候我都佩服自己能想出可以用图的形式来表示突变的密度，哈哈哈,因为这样就能很直观的看出哪些地方突变频率比较高，突变频率比较高的地方，可能是突变热点区域，在肿瘤研究中常常有意义。

我合并多个文件，用cat将流重定向到一个文件，或者把一个大文件rm掉腾出空间，要进行IO，但如果这个时候服务器有进程进行IO时，同个进程同时进行IO，效率就会很慢。有时候我想把别的进程IO缓一缓，先把合并或者rm的任务有限解决掉，再继续别的进程的IO。

就google了下如何提高效率，查ionice这个命令。

ionice - 获取或设置程序的IO调度与优先级，通过设置命令或进程的IO调度优先级，加快IO效率

命令格式

跟命令时，设置命令的IO调度优先级，跟PID时，设置相应进程的IPD调度优先级。

1
2

ionice [[-c class] [-n classdata] [-t]] -p PID [PID]...
ionice [-c class] [-n classdata] [-t] COMMAND [ARG]... 

我也进行了测试，比如在有其进程进行IO工作时，我要在删除Fastq文件（R1文件20G，R2文件21G）时，用法详见下文。

1
2
3
4
5
6
7
8
9

time rm R1.fastq
real    0m37.306s
user    0m0.001s
sys     0m0.300s

time ionice -c 2 -n 0 rm R2.fastq
real    0m2.579s
user    0m0.002s
sys     0m0.682s

可以看到当提高IO的优先级后，效率还是非常快的，当然这暂时牺牲了其他进程的IO。

ionice -h用法