默认是没有覆盖的情况,只是去除了同start和同end的覆盖,更大范围的覆盖并不去除,实际有没有并没考虑。
选择一种查找算法来查找范围。
树的问题是最上的根节点比较麻烦,哈希表不能比较范围,索引查找比较麻烦,相对而言还是最初的二分查找比较合适。
因为都是有序的,查找时间可以保证,也能避免位置冲突带来的哈希表和索引、建树问题。
二分查找的话,先找到右边界,再-1看有没有大于左边界。
找边界:
当left>right时,也就是left = 15 right = 15,(在此代码条件下)
对右-1得左:
下标应该用first,last,之前用mid会出错。
只能取到范围内非最左端,如果end是最左端,则结果是None。
改动条件使得它只能不跨或只跨一个外显子。
搜索模块结束。和sam文件组合筛选。
需要知道start和end,在pysam里选一个方法:
可以得出只有M的长度。
直接相加和用reference_end一样
cut一下
sed一下。以便将文件转为list,用:
有了位置的list,有了比对算法,可以对sam中的位置信息进行比对了。
跑出来结果不对。突然想到intron位置在最后一次sort之后保持的倒序,但程序写的last last-1是左右位置。
正序
结果:
一对的常只取到一个。
2
目前看没有乱取,但是有少取,总是只取一半。
找第一个比对到的(seq_in_intron中)
能看到只比对到右边的,而且不完全,和之前的特征一致,都是比对到右边。
但是,输入文件是排序后的bam,输出文件中,按igv看到的,本应是连续的,但实际上:
根据igv的intron地址,回intron位置文件中找,位置并无问题。
做了不大不小的改动之后,也仅仅是多了几条。改动是使用sort过后的文件作为输入文件,以及代码中错写每次循环都重算位置。且一个read的开头结尾只能在一个intron,应该删除一下代码多判断的部分。
重新写后结果依旧没变,突然想到是不是0base和1base的问题,因为左端是对齐的,所以可能差个1就比不上了。并且检查了0819和1041,确实是完整的intron,能比到右端是1041,左端一定是0819.
gff3可能是1base,即开始的第一个位置是1,而pysam取出的序列是0base,他的start要比gff的小1,应该在start+1.而判断条件是start>=array[last-1].因为start和它本身切齐,但又小个1,所以永远不可能>=,除非是比到右端。
马上重试。
马上序列多了一倍。
bug看似被成功修正,再检查一下应该能下班。
启事是:1.如果不是勉强自己仔细看了pysam的文档并对不懂的地方查了,应该想不到这个问题,而是像之前一样以为是二分法和位置序列取的有问题。认识自己使用的工具,深入的,才能找到问题可能出现的位置。所以,认识自己使用的工具,并多读文献文档。即使这次是侥幸看到,我还是有大面的文档没有阅读。很大面积。我只是半仔细的看了可能要用的部分。所以,从头至尾的读,读作者的话是非常重要的。这次只是侥幸扫到。
2.在疲惫的时候写出的代码不严谨,设计时候的不严谨再改正时会因为思维僵化而找不到问题。并且,之前我有提示自己注意01base 的问题,所以,边做边记是个好习惯。会继续发扬下去,并且,设计时候考虑的疑点,也要如实记下。
总结:认真读文档,多读文献,设计时要严谨 ,设计时的疑点要如实记下,设计时的每一步记录要做好。
一个简单的问题实现起来用了两周。。。。
结果:最后又浏览了一下,表现尚可,大体比较准确,但不是每个都能保证正确,有些能看到到未取到。还有要改进之处。
所用脚本:seq_in_intron.py