当awk碰到百分号%数值时

awk可以通过$1,$2等,对数据进行逻辑判断或处理,比如

1
2
> echo -e "1t2n3t4" ' awk '{if($1==1){print $0}}'
1       2

但如果碰到带有百分号%的数据时,则不能直接通过加减乘除进行计算或者判断

1
2
> echo -e "1%t2%n3%t4%" ' awk '{if($1==0.01){print $0}}'
>

因为系统没有将1%转成数值0.01,无法进行判断$1是否等于0.01

样本突变频谱分析和基于突变频谱的热图聚类

对突变频谱(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()

多线程gzip压缩神器---pigz

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表示基因组上的突变密度

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

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

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