首先依次查看每一类的序列文件
one_pair._4sam
问题
突然想到问题所在：只在变qname时候存储，反而变qname且r2也计数时未存储。
改了之后在跑就没有这个问题了。毕竟两个判断qname不等的，只有一个去存储显然不对。
改了重跑再抽样检查没有问题
one_pair_but_multiple._4sam
无误
two_pair._4sam
无误
two_pair_but_one_multiple._4sam
无误
two_pair_but_two_multiple._4sam
无误
五个分类文件内容都无误。
再改变文件格式放igv看
由于之前取的序列有问题，显然之前得到的解释igv结果也有问题，但也不是没用，经过前面的测试，现在可以在确保序列正确的情况下对每类序列进行查看。
位置在sort_and_index文件里
one_pair.sam
看来还有问题
还有空白块，看来还有问题
对其中一条：
SBSUSER:472:HBGC7ADXX:1:2211:12051:97803
文件里是两条，但还是白色，这里第一条的位置显示很奇怪
查看另一条：输入位置后能找到另一条，猜测白色可能是两段距离太远导致，并非取的少了一条。
另一条
另外还能看到扎堆的白色条带 ，猜测是因为保守区？容易被比对到？
扎堆
这几条的insert size都是九百万个碱基以上。感觉应该不合理，但是不知道是不是需要加上限制间隔长度的这种条件。
再取一个查看：
SBSUSER:472:HBGC7ADXX:1:1109:18823:9664
找没取之前的sam：
未取之前的比对文件
可见，99和147是r1 r2，所标的位置也是r1 r2相互之间的位置，而r2的supplementary则和r2相隔900w碱基。且形式也不对，不该r2开头是S supplementary开头H。
可见需要加的校验还有。【8/22认为过滤MAPQ试试，问题还存在则加限制TLEN的条件】
8/21过程记录
查看的文件：one_pair_but_multiple
结果：明白了SAM里所有比对结果的表示是统一的，都是从左到右，不必考虑原序列测序方向。
对下一步的指导：增添新的取序列限制条件
one_pair_but_multiple（一条read有超过两个的比对）：
结果好的：
感觉上最想要的结果
1
2
3
结果有问题的：
还是很多白块 insert size 900w
insert size 900w
类似结果分析：
1
2
3
4
突然发现 ，最后的是41M60H，本来我以为先M的应该在前面，但是他在最后，因为他是负链的原因？
41M
SBSUSER:472:HBGC7ADXX:1:2113:13615:43628
有可能是按链的顺序写的，因此从后到前是先41M再60S【从最后做完以后，知道这里是错的】
回头分析前面的正链，最左边的是先MS 中间SMS 最后HM 顺着的，应该是不对的。【也是错的，正确色是：与原序列是FR无关，sam里对FR是一样的的，不需要考虑原序列的方向。】
再多找几个
1
2
3
如果是按上一个的顺序，第一个的12M应该在最后，但是这个在最前面，主要是负链的顺序不清楚。
先M且在前面的也有
1
2
pos
即使看了pos的介绍还是不清楚，双端测序到底怎么看顺序。
为了得出顺序，从fq sam igv 对同一条序列进行查看
SBSUSER:472:HBGC7ADXX:1:2209:12227:36436
_1fa
_2.fa
sam
最左边，12M89H（12nt）【得出：cigar显示的是从左到右的】
_2fa显然是和sam里161的是一条，序列也没变；而_1fa则反向，反向以后从左到右，同igv结果一样。
因此，CIGAR序列显示的也是从左到右，不管是F还是R，都是一个顺序。都是从左到右
只是比对起始位置不同，像现在这个是：
起始【得出：位置和sam里一样，是从左到右的】
结尾
结论：和fa里碱基顺序无关，不管它是F还是R，sam里的比对顺序都是从左到右，cigar标示的也是从左到右，所以只要看位置和M即可。
主要是cigar一直是从左到右的，与FR无关这个重要。像sam里把R序列reverse倒是正常，也一点儿不影响比对，就相当于是两个forward序列，所以可以预测insert小的的两端都交叉比对的序列非常可靠。
由此进一步筛选程序取出的序列。
因此，这个序列并非交错比对，反而前面的是正确的
SBSUSER:472:HBGC7ADXX:1:2113:13615:43628
1
2
可见，第二条：41M60H，位置：8113 363
第一条：54S47M，位置：8113 113
先M的靠后，而后M的在前面。
这个是对的。
收获是明白了SAM里所有比对结果的表示是统一的，都是从左到右，不必考虑原序列测序方向。
结论的得出主要是看fa sam和igv显示的比对位置+参考基因组碱基顺序。
注：今天的得出的结论适宜重新整理记录（关于原fa序列顺序在SAM文件里的影响和作用），因为是记录过程，所以过程中有错误的猜测，虽然标注了，但是并不打算修改，想留下原始的过程记录。可以将重要的东西另成文章记录。记录一个过程我认为也很重要，虽然可读性十分抱歉。
可以进一步添加筛选条件修改取序列的代码
因此下一步的计划：1.根据上一章和这儿一章的问题，修改代码。
2.才看到第二个文件，把接下来的文件看了找问题
3.igv官网文档比较丰富，看一下不同颜色代表的含义，主要是空白块。【8/22已解决】
对这个文件研究的不足： 中段的HMH序列不太明白其含义。
8/22
IGV着色问题
官方文档：http://software.broadinstitute.org/software/igv/AlignmentData
1
2
3
insert size 并不是我想的insert fragement，而是插入缺失的插入。
因此空白的主要指Low mapping quality
igv上看条带也确实如此
MAPQ=0
SBSUSER:472:HBGC7ADXX:1:2110:2753:46772
MAPQ=0
已知：空白块表示MAPQ=0
待解决：1.但是不知道和insert size有无关系
2.需要过滤MAPQ分数在多少以下的reads
two_pair
1
发现：密集出现
2
不知是否是偶然
开头完全一致
一一查看过是完全不同的reads，但是开头却一致。这种特征也能被使用，以增加结果的可信程度
（不同文件，不同reads，比对起止位置一样）
挑一个序列分析：
非标准的：
1
MAPQ=0 SBSUSER:472:HBGC7ADXX:1:1209:3827:75074
3
选左边的R和右边的F，但是两条却并不是同一个fragement上测序出来的。需要看是否取错了【已解决，并未取错，而是可以提供验证的信息】
没有错
在IGV上找对应链：
75074的forward链
SBSUSER:472:HBGC7ADXX:1:1209:14484:50628reverse链
实际是一个fragement的
有趣的是看似是一对的两个却并不是，但是有同样的断点。因此猜测，是不是因为这是真的BP，才会出现这样的情况
验证：
pnext
首先，pnext这种值对我们这种情况是没什么用处的，连带的TLEN用处也不大，在这里我们主要看的位置值是下面的
位置值
后M的，位置在后，因此不是bp
再检查50628：
pos
这个下面的是FORWARD，原因未知，但是下面的显然也不是合规范的bp序列形式
提出这个问题：
F1R2
F2R1
标准情况：
正链1
正链2
负链1
负链2
观察到insert size在同一条read的两段是不同的。
重叠
非典型特征序列：
1
正链1
正链2
负链1
负链2
有必要考察TLEN的含义。
官方文档
带自己注释的翻译：标示模板长度，如果全部片段（一个read比对中，分成多个片段比对）比对到相同参考基因序列，那么TLEN绝对值=比对的end-star+1（这个信息的+-有时不是很准确，绝对值准确，因此不能拿来作为判断依据）。注意，比对上的碱基是在CIGAR上描述的，且不包含软裁剪的那部分碱基。TLEN还这样定义：最左的片段（注意是片段）为+，最右的片段为+，中间片段不定义符号（也就是不参与计算，只计算最左和最右的差值）。（但也不是万能的，对于特殊情况，还有这样的定义：）（特殊情况1：）如果片段覆盖相同的坐标（应该是说两个片段重合了），那么最左最右的选择是随机的，但必须有不同符号（这种情况下+-什么也不能代表）。（特殊情况2：）只有一条比对时TLEN记为0，或者另一条不可用时
这个信息暂时看来不好用。到此，总结需要改的代码和流程：
1.过滤MAPQ=0 的read
2.加限制条件：先M的pos>pnext
8/23
完善修改前面的问题，并提出新问题
关于加新的限制条件
改善流程
参考文件如下：
MAPQ
https://www.jianshu.com/p/6b7a442d293f
view的参数
选择命令为：
samtools view -bhS -q 1 SRR20150106_mapped > SRR20150106_mapped_withheader_mapq_nonzero.bam
有header的bam取出的特征序列文件里没header，所以后面的处理和之前一样。
结果：最直观的是文件大小的变化。
但更主要的是验证先bam 再用bam取特征序列的流程是可行的
文件结果1
文件结果2
igv结果
红线部分是去MAPQ=0前后的对比，去除后连带少了一部分read，这样比较合理，毕竟如果比对片段的一部分比对质量太差，也是不可信的。
后续也检查了其他文件，去除后结果表现很好。
空白块的小问题是解决了，大问题是M在前的pos在后，M在后的pos在前。这个需要同一个read的不同片段去比较。在取的同时处理似乎不容易，先对取出的序列文件进行二次加工。
现在的情况有：one_pair tow_pair最简单，老办法两两对比即可，
one_pair_but_multiple, 主要解决multiple问题，
multiple的问题主要是1.multiple里有很多HMH SMS，该怎么解释，并且mapping得分还是60
2.和谁比位置
3.进一步提高比对质量（MAPQ）的要求，看到很多2的
注意，对每个文件需要去除没有两位数直接两个或者多个四位数的。
multiple
需要用新的参数来帮助取
文档
解释：
结果
每一行是一条，先看后面的，写1指代有，0指代无，因此白色第一条是：MH有。在看前面的，可知M=24，H=77
那么MHM怎么表示？前后位置怎么表示？这里都看不到，因此这个参数不适用当前问题
不能知道位置只能知道有没有。因此选择用正则表达式处理。
写到一半，明天处理
8/24
正则表达式
参考1
参考2
代码
结果
对文件的二次筛选得出的结果确实是想要的。也是最标准的。
只是个demo，还将扩展。
教训：想当然的将.* 写成了* 改了一天......
8/25
扩展demo，使之能对所有文件起效
由此改变流程：先取有特征序列，再过滤取出的特征序列，做两遍。取时先不分类型，即不分之前的五种文件，只有一个中间文件temp，再过滤，过滤出5种类型。
代码思路还是和之前类似：1.同时取两行，但是只对第一行操作。2.用count变量计数来判断是第几行。不同时这次因为要与flag<200的比位置，因此需要存储flag<200的那一行，而如果第一行是flag>2000的，则所有此qname 且同read的都不要。
结果
部分主体代码
因为测试文件是multiple 的，因为会过滤所以原来比对条数>2的现在会<=2，也会出现只有一条的，这个是写的时候先让符合count=0 flag<200的输出，需优化，还未优化。
好的结果是：取出的序列符合要求，且数目大量减少，橙色圈出的是每个第一行的比对位置，以最后一个为例，第一行CIGAR:34S67M,POS:12979078 其下序列第二个： CIGAR：33M68H POS:12979202
第一行后M，078在前
第二行先M，202在后，符合需要的情况。
网友评论