测试项目:不加限制条件的情况下,最多的overlap数目
在30为:
3031
3132 33
32 3330 33 像一个数据。。。
结果少,不正确。检查方法。
改变了过滤和max的方法,识别出的数量有所增加。
1 1到32都很少,需要和别的再比。
将目前最好的30和别的数据库比对,结果也不算特别多。
130 的sam8G 能出40个,但是33的fa更大,应该是33有问题,而最小的30都能出4000个结果,比他大的31 32 33应该 5000 6000 8000,有2w个结果。
因为可能是坐标系的问题,pysam读的坐标和sam的差1,因此根据之前的测试,end得到的本身是能取到的最后一位,因此在v5更改了坐标,更改后在最大的文件中能看到多了很多,但是在31依旧不够。
132:32依旧很少,
3231:31多一点但是作用不大
31尝试全跑:
30:
30依旧不够
31:
3133 2w才能够,但是差的有点多了
测试变化坐标前后取出的值有无重复
1.变化前的:v4
1uniq一去少一半
2.变化后:v5
1掉到5000,
3.看重复
1 1能看得出,bedtools的比对是严格的,且谁是-a都无影响
1重新bwa以后33的数量确实变多了
1从fa产出sam的数量看:
13G产生 8G sam 能有5000 uniq的bp
25 11 5000
32 12 5000
50 31 2W
这样的比例才合理
但是30 31 32 都是一样的脚本跑的,3132就很少。而且30 31 是一样的细胞类型
1现在尝试看1.和所有三个数据库比2.跑untrite
跑完以后是:
1nohup python full_human_1_6.py -i /mnt/T30/wus/brantch_point_human/Mercer_data/SRR1049833_mapped_clear.sam -f /mnt/x110/guosy/Database/hg19/samtools-index/hg19.fa -o /mnt/T30/wus/brantch_point_human/Mercer_data/SRR1049833_1.6 &
重跑如上
bwa会输出多余的东西在sam中,需要nohup来限制,从中清理可以用删除的方式,sed '//d'
1 1查看指定行的办法:
1. 从第1000行开始,显示2000行。即显示1000~2999行
cat input_file | tail -n +1000 | head -n 2000
2. 显示 1000行到3000行
cat input_file | head -n 3000 | tail -n +1001
*注意两种方法的顺序
分解:
tail -n 1000:显示最后1000行
tail -n +1000:从1000行开始显示,显示1000行以后的
head -n 1000:显示前面1000行
同时测试在各个数据库中的表现:
30首先是30,能看到在最大数据中的表现和在给出文章大的数据库中表现几乎一致,但是这些数据库对特殊类型的支持可能都比较一般。
31 32
1虽然在不同数据库中表现出差异,但差异没有太大,依旧不够数量。
15.9w的数据库,
1.7更新之后过滤出的大小和之前一样。本来在igv看的one_but_multiple全是一个方向感觉对的,没想到还是不对。
检查v5并测试别的跑1.7的是否改变了
需要测试:31-33取出的序列有多少在intron上,2.比对正确的那些在igv看是否就是在intron
1对31看,确实是在intron
再看也是正确的。
1而且与方向无关
2可见取到的是正确的,过滤出的条数也很多,如果没错可能只能是全是外显子环形。
问题是31的在内含子的更多
1如果取出来并看一下,能解决这个问题。
1这个产生了更多,不该是这个结果,
1能看到这个是没被选中的,
回到初步过滤文件中看到是对的上的,末尾也没错 grep 'SRR1049831.45901126' one_pair_SRR1049831_mapped.sam
1检查原始sam,一种是没取全,一种是mercer的错误,排序检查?但是其他数据库结果证明他应该没错,
1原始sam查看发现也没错,只能在这个位点对mercer进行grep
1grep不到
1大规模找一下发现没11929的
找第一个点发现没序列
1chr10的这几个点全找不到,是没序列的那种找不到。
chr1又找到一个点,这个点显然也对
1这个点依旧找不到:
1猜测:他的取的是250nt的,这个比较长超出去了,但是我的应该兼容他的才对,而不是现在的比他的少。
grep 'SRR1049831.30557431' ../one_pair_SRR1049831_mapped.sam
grep 'SRR1049831.30557431' ../../SRR1049831_mapped.sam
1这么看过滤的也没问题,这个位点应该就是bp
找最接近的:
1问题是没序列,那么需要把最原始的sam加进来,看这个位置到底是没有还是没过滤出。
1bwa原本的没有这个点。
1再换个位点验证
1 1没有重叠区
并且原始sam中也没这个位点
1只剩‘为什么在多个数据库中表现一致’这件事存疑,其他的办法可尝试换数据试试
33的结果也不多
1发现个有问题的:
1但是经查看发现确实两个都是read1
1 1 1不知为何反向 。
用过滤的方式结果过于依赖bwa的比对,一旦比对和别人的不同,最终结果也会不同。
1能力看到和mercer比得上的确实没问题
1确实没问题
网友评论