I have found the following behavior in bxpython which sounds like a bug to me:
I have a loop where I create a new Interval tree and check how many clusters it contains for the same data. It outputs:
With max_gap: 0 #clusters is: 156
With max_gap: 1000 #clusters is: 151
With max_gap: 2000 #clusters is: 155
With max_gap: 3000 #clusters is: 155
With max_gap: 4000 #clusters is: 156
With max_gap: 5000 #clusters is: 158
Is this possible? I thought that number of clusters should decrease monotonically with the length of cluster_distance?
Just wanted to double check before I report this as a possible bug.
The relevant code is:
max_distance = [i for i in range(0,5001,1000)]
for max_gap in max_distance:
temp_tree = build_cluster_tree('chrX2.map' 10, max_gap)
print "With max_gap:", max_gap, " #clusters is:",len(temp_tree.getregions())
More infor about Bx-pythons clustertrees here: Finding and displaying short reads clustered in the genome
1 answer
It might have something to do with the minimum number of reads you defined to be considered a cluster. As you increase the cluster distance, you might be adding more reads to previously discarded clusters, allowing them to pass the minimum reads check.
I wonder if you can relate this monotonic trend to distribution of reads across your reference. A set of reads that are very evenly distributed would probably show a large number of clusters initially and a sudden decrease to a single cluster after a certain cluster distance. Whereas something that is not evenly distributed will probably show a moderate amount of clusters initially and a slow decrease to single cluster as you increase cluster distance.
Log in to answer this question.
As an aside: list comprehensions are cool and all, but
range()outputs a list anyway.