评估文库 Average Insert Size

用SOAPdenovo对Illumina paired-end进行基因组组装时需要配置文件,其中要填写每个文库的average insert size,那么如何进行average insert size大小的评估呢?

文库类型

对于基因组文库我们一般会建小库(<1K)的paired-end reads和大库的mate-pair reads,二者最主要的区别就是reads1和reads2的方向和之间的间隔大小。



现在绝大部分的主流软件都是支持将paired-end reads进行比对的,那么 mate-pair reads如何处理呢,即 mate-pair reads如何做比对?请参考我的另一篇博文 Mate-pair Reads Alignment

Insert Size

首先,什么是Insert Size呢?
对于paired-end reads来说

对于mate-pair reads来说其reads1和reads2方向指向外面,其插入大小统计需要格外注意。

基于bwa比对log文件统计插入大小

通过观察bwa软件的输出log文件发现其对每一个pair-end reads分4次读入([M::main_mem] read 2613712 sequences (200000105 bp)…),对于每一次的读入会对reads进行统计如下:

由红框标出发现占主要比例的是RF reads,进一步往下寻找analyzing insert size distribution for orientation RF…就可得到其平均插入大小。

基于比对的sam文件统计插入大小

R计算

1
2
3
4
5
6
7
8
$ cat sample.sam | cut -f9 > initial.insertsizes.txt
R
a = read.table("initial.insertsizes.txt")
a.v = a[a[,1]>0,1] # 筛选大于0的值
mn = quantile(a.v, seq(0,1,0.05))[4] #分位数计算,[4]表示取第四个分位值,15%
mx = quantile(a.v, seq(0,1,0.05))[18] #85%
mean(a.v[a.v >= mn & a.v <= mx]) # mean
sd(a.v[a.v >= mn & a.v <= mx]) # sd

可见R计算过程选择过滤掉小于等于15%和大于等于85%的值来计算平均插入大小;

awk计算

1
awk '{ if ($9 > 0) { N+=1; S+=$9; S2+=$9*$9 }} END { M=S/N; print "n="N", mean="M", stdev="sqrt ((S2-M*M*N)/(N-1))}' sample.sam

awk选择全部大于0的值计算平均插入大小;

基于sorted.bam文件

qualimap

1
qualimap bamqc -bam sample.sorted.bam --java-mem-size=300G  -c -nw 400 -hm 3 -outdir /resulted

CollectInsertSizeMetrics

1
2
3
4
5
java -jar CollectInsertSizeMetrics.jar \
I=sample.sorted.bam \
O=insert_size_metrics.txt \
H=insert_size_histogram.pdf \
M=0.5

问题

现对同一个比对产生的sam/bam文件用上述4种方法计算得出结果如下:
R :mean insert size = 260.577232343453”,”standard deviation = 27.4153198790634”;
awk: mean=250.826, stdev=51.7005;
qualimap:mean insert size = 250.8258,std insert size = 51.7005;
CollectInsertSizeMetrics:MEAN_INSERT_SIZE=260.963523,STANDARD_DEVIATION=42.809159。
qualimap 和 CollectInsertSizeMetrics 都是java封装的软件,看不到其具体计算方法,根据以上计算结果可以看出CollectInsertSizeMetrics的计算原理应该和R的一样需要过滤掉数据,qualimap和awk中发一样,所以问题最后就归结为是否需要首先过滤数值再计算平均插入大小?

在R中计算时对数据a.v做正态性检验lillie.test(a.v)

1
2
3
4
5
6
7
> library("nortest")
> lillie.test(a.v)

Lilliefors (Kolmogorov-Smirnov) normality test

data: a.v
D = 0.1541, p-value < 2.2e-16

可以看出其插入大小分布不是呈正态分布,综合考虑后还是按照R的计算结果为准。

参考资料

Question: Estimate Insert Size In Paired-End/Mate-Pair
Question: What is the difference between a Read and a Fragment in RNA-seq?
Paired-end read confusion - library, fragment or insert size?

tiramisutes wechat
欢迎关注