Asymmetric trimmomatic output with paired-end sequencing reads

2016-05-12
Default Category

运行trimmomatic的默认参数

1

java -jar trimmomatic.jar PE -threads 16 -phred33 sample_1.fastq sample_2.fastq sample_trimmed_paired_1.fastq.gz sample_trimmed_unpaired_1.fastq.gz sample_trimmed_paired_2.fastq.gz sample_trimmed_unpaired_2.fastq.gz ILLUMINACLIP:adapters/TruSeq3-PE-2.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

输出文件中

sample_trimmed_unpaired_1.fastq.gz是sample_trimmed_unpaired_2.fastq.gz的十多倍，异常的大。1文件（forward reads）中unpaired reads非常多，显著多余2文件（reverse reads）中的unpaired reads。

阅读更多

当awk碰到百分号%数值时

2016-05-05
Default Category

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

阅读更多

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

2016-04-03
Default Category

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

2016-04-02
Default Category

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
阅读更多