File "/usr/local/opt/anaconda/lib/python2.7/site-packages/ipyrad-0.1.71-py2.7.egg/ipyrad/assemble/cluster_across.py", line 427, in singlecat
icatg[iloc] = catarr[iloc, :icatg[iloc].shape[1], :]
raise TypeError("Can't broadcast %s -> %s" % (target_shape, count))
TypeError: Can't broadcast (4, 4) -> (1, 150, 4)
Found the bug in step 3 for reference mapping. I'm cleaning it up and verifying it's fixed.
OK. There were two problems that were causing the high indels in refmapping. Problem 1) I was fsck'ing up the writing to clust.gz so the last unmapped cluster and the first mapped cluster were not separated by
//. Easy fix.
2) The other thing that's causing high indels is an artifact of the hackish nature of the simulated mitogenome we're using. I generated this by randomly grabbing reads from the
sim_rad data and randomly sticking them into the mt genome from the mouse. Good enough for government work, but randomly some of the simulated fragments overlap enough so that bedtools merges them, in which case they are joined as one cluster and then muscle tries to align them with obvious consequences.
For the sim data this affects 3 clusters total, i don't think there's anything we can do about it unless we simulate a better mitogenome. I'll push the fix for the first problem momentarily.
OK. Thanks. Sounds good. Well have to start testing on real data more to make sure it's working. It's good that the filter is excluding the bad clusters.
Yah, once i get glenn's gbs data run all the way through
denovo i'll start testing refmapping. should be interesting. you in idayho?
Ipyrad workshop is Sunday
So theres a little time left for testing
But ask the sim datasets should be working now
And I'll do preview mode on an empirical dataset
It seems to be working better in so far as the alignment bug is fixed. But why do the samples recover different numbers of refmatched reads? It seems they should all match to the same regions of the reference, but I'm getting between 50 and 80 clusters recovered among the different samples.
Is this a problem with how smalt identifies homology, or with our treatment of the mapping results?
I think the problem has to do with refmapping sticking a bunch of inserts or Ns in the reads somewhere. This is what I see in the .loci file from a finished refmap run
1C_0 TGTTTCCTAGACGGGAGATCATC-----TGAAATTCGGTT-TTATCCTTGACCTGGTGCCACCGCCTCAATCTCTTGATTTCAGACCAAAAGAA
1D_0 TGTTTCCTAGACGGGAGTTCATCNNNNNTGAAATTCGGTTNTTATCCTTGACCTGGTGCCACCGCCTCAATCTCTTGATTTCAGACCAAAAGAA
2E_0 GGTTTCCTAGACGGGAGTTCATC----NTGAAATTCGGTTNTTATCCTTGACCTGGTGCCACCGCCTCAATCTCTTGATTTCAGACCAAAAGAA
2G_0 TGTTTCCTAGACGGGAGTTCATCNNNNNTGAAATTCGGTTNTTATCCTTGACCTGGTGCCACCGCCTCAATCTCTTGATTTCAGACCAAAAGAA
3I_0 TGTTTCCAAGACGGGAGTTCATCNNNNNTGAAATTCGGTTNTTATCCTTGACCTGGTGCCACCGCCTCAATCTCTTGATTTCAGACCAAAAGAA
// - - - |
1B_0 TCGTGATGGTGCTGACTCAATTGGTACTCTGGCACGAT-----TATGCCGCATAGCACCACCCGTGCCGGATTGGCACCGNNNNNNNNNACAAGTGTTAGGC
1C_0 TCGTGATGGTGCTGACTCAATTGGTACTCTGGCACGAT-----TATGCCGCATAGCACCACCCGTGCCGGATTGGCACCGNNNNNNNNNACAAGTGTTAGGC
2G_0 TCGTGATGGTGCTGACTCAATTGGTACTCTGGCACGAT-----TATGCCGCATAGCACCACCCGTGCCGGATTGGCACCGNNNNNNNNNACAAGTGTTAGGC
3I_0 TCGTGATGGTGCTGACTCAATTGGTACTCTGGCACGATNNNNNTATGCCGCATAGCACCACCCGTGCCGGATTGGCACCG-----NNNNACAAGTGTTAGGC
// |
ugh, something is going wrong here. Reads with than many Ns should have been removed in step5 anyway...
I'm trying to get a sense of the refmapping code. Do you refmap the reads before we dereplicate them?
I think I'll have to skip doing refseq for the workshop. It seems close, but there's definitely some major bugs to work out. I should have tested with it earlier. We should plan a meetup some time soon so we can really dive into the code together.
Refmap happens before dereplication, but then it dereplicates the mapped reads itself.
What was the idea behind doing it after vs before?
Increasing clustering threshold solves some of the problem (at 0.90 you only see very few that are fsck).
Right now we use clustering threshold for both vsearch clustering and refmap sequence similarity. It's a simplifying shortcut, but not sure if it's justifiable. We might need to add a parameter for refmap sequence similarity if it turns out they need to be tuned independently. Needs more testing.
I guess the idea of doing refmap before dereplication was about potentially using quality scores during the mapping process.
It also simplifies the code a bunch. Dereplication happens inside clustall() after merging pairs, so we also want to map paired reads independently.
Note from Julian's talk: random oligos to deal with pcr duplicates.
At clustering threshold of .9 all the reads that show this bad behavior are a result of simulated reads overlapping regions. I'm investigating lower clustering thresholds (#126). In practice this shouldn't be a major issue. Also, i'm not sure if i'd expect all samples to recover the same number of refmapped sequences (because of the way the sim genome was constructed (by randomly selecting reads from simrad.gz)). The simulated mt genome isn't necessarily going to do anything very interesting, it's more a proof of concept.
It doesn't use qualscores now, but i wanted to have the option for the future (V2.0), more importantly is the unmerged pairs.
I created a new
sim_mt_genome that isn't so much of a shit show. It has 24 sequences inserted from sim_rad_test_r1_ so that none are overlapping, no high indels errors even at low clustering thresholds.
I also changed the default smalt wordlen (in hackersonly) so it works better on the sim data. Will have to experiment on real data to get a better sense of default, but my experience is 13-16 is best.
oh hell yeah, the implemented the new solver in conda (>4.0), it's really fast.
hells yeah. Sounds awesome. I'll play around with the new stuff while I watch Julian's presentation.
did you push that yet? I see smalt_wordlen=8 in my version
8 is the optimal value for the simulated data, so you've got the newest ver
I might need you to push a mac conda version before tomorrow of whatever version we end settled on by the end of today.
I made some minor changes to refmap, before you made fixes. didn't find any conflicts on pull, but maybe you should check that I didn't mess anything up.
I can calculate inner mate distance for mapped pairs, rather than hard coding it (dumb). On real data do you think it's better to use max inner mate distance or mean +2 standard deviations? I'm inclined toward the latter.
Not slower, just the max insert size could be a significant outlier