BAM文件插入片段长度与碱基质量均值绘图异常求助

初始实现与问题 #

我使用PySam计算BAM文件中reads的插入片段长度与平均碱基质量，代码如下：

import math
import pysam
import numpy as np
import matplotlib.pyplot as plt

bam_file = pysam.AlignmentFile("mybam.bam", "rb")

mean_base_qualities = []
insert_sizes = []
zipped = []
for read in bam_file:
    if read.is_paired and read.is_proper_pair:
        # length
        mean_bq = np.mean(read.query_qualities)
        read_length = read.query_alignment_length
        tlen = read.template_length
        insert_size = int(math.sqrt(math.pow(tlen, 2))-read_length)
        if insert_size > 0:
            zipped.append((mean_bq, insert_size))
            mean_base_qualities.append(mean_bq)
            insert_sizes.append(insert_size)


plt.plot(mean_base_qualities,insert_sizes)
plt.show()

最终得到两个大型列表，示例如下：

insert_sizes = [291, 188, 165, 212, 186, 127, 240, 166, 120, 303, 255, ...]

mean_base_qualities = [36.24324324324324, 34.73831775700935, 35.3448275862069, 30.64788732394366, 36.467289719626166, 22.673267326732674, 36.54216867469879, 31.20754716981132, 25.085714285714285, 34.65625,  ...]

绘图后得到明显异常的图表，尝试独立排序列表后绘图仍无改善，其他绘图方式也均失败。

更新后的代码与新问题 #

根据建议，我对平均碱基质量按插入长度进行聚合均值计算，代码如下：

import math
import pysam
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import make_interp_spline

bam_file = pysam.AlignmentFile("mybam.bam", "rb")

mean_base_qualities = []
insert_sizes = []
zipped = []

i = 0
# No longer needed to use read.query_alignment_length for substracting to TLEN value
# TLEN represents the insert size properly
for read in bam_file:
    # Only processing paired proper reads
    if read.is_paired and read.is_proper_pair and read.is_read1:
        if read.template_length != 0:
            i+=1
            mean_bq = np.mean(read.query_qualities)
            i_size = int(math.sqrt(math.pow(read.template_length, 2)))
            zipped.append((i_size, mean_bq))
            mean_base_qualities.append(mean_bq)
            insert_sizes.append(i_size)

res = [(key, statistics.mean(map(itemgetter(1), ele)))
       for key, ele in groupby(sorted(zipped, key = itemgetter(0)),
                                                key = itemgetter(0))]

tuples = zip(*res)
list1, list2 = [list(tuple) for tuple in  tuples]  

xnew = np.linspace(min(list1), max(list1), len(list1))
gfg = make_interp_spline(list1, list2, k=3)
y_new = gfg(xnew)
plt.plot(xnew, y_new)
plt.show()

但生成的图表仍然不合理：出现了碱基质量均值的负值，而实际我的碱基质量均值范围为18.677至37.0，插入片段长度范围为2至25381844。恳请帮忙找出问题所在！

内容的提问来源于stack exchange，提问作者P. Solar