实验记录 | scATAC-seq数据的比对（五）

午睡醒来有一些坏脾气。最讨厌不思考，不独立的人。于是发了一些脾气。有被打扰到谢谢！查文献难道不是研究生的基本素养吗？人还是要主动学习的？因为知识不是自己跑过来的，因为老师不是自己跑过来的。而是需要你自己主动去争取的。在这里，我还是非常感谢肖老师的，她教会了我研究生的一些基本的本领和基本的原则，这对于本科生而言，足够了。
我还是挺感恩我的本科老师的。我从来不放弃思考，从来不放弃主动求知识。

---

今天要解决负链的问题。
1、就是说，如果我们的比对结果比对到负链，那么会是怎样的一种呈现方式？
参考链接：https://www.jianshu.com/p/56aa023decd1

###### forward strand
samtools view -b -f 128 -F 16 a.bam > a.fwd1.bam
samtools view -b -f 80 a.bam > a.fwd2.bam
samtools merge -f fwd.bam a.fwd1.bam a.fwd2.bam

###### reverse strand
samtools view -b -f 144 a.bam > a.rev1.bam
samtools view -b -f 64 -F 16 a.bam > a.rev2.bam
samtools merge -f rev.bam a.rev1.bam a.rev2.bam

需要理解这里的flag值的意思：
（1）flag=128，在双端测序中，该reads为双端测序的第2条。
（2）flag=16,与该序列反向互补的序列能够比对上。
（3）flag=80=64+16,在双端测序中，该reads为双旦测序的第1条，且与该read反向互补的比对到染色体上。
总结一下，上面筛选的即是先将双端测序的两条reads区分开，对于其中一条，其本身能够比对上，其反向互补不能。另一条则是，其反向互补能够比对上。
（1）flag=144=128+16，在双端测序中，其为该reads的第2条，且比对上的是其反向互补序列。
（2）flag=64,其为该reads的第1条，且比对上的是其本身，其反向互补序列无法比对上。

（我觉得上面的还要再加值，比如该reads是paired-end中的一个，paired的每个都能够正确的比对到参考基因组上。）
2、上面的结果中提到了双端测序中的第1条，第2条的reads？这个和我们基因组上的5’-3’有啥关联？我之前的理解是第2条就是负链的方向？这样对吗？

首先，双端测序中，R1和R2的关系是互为反向互补序列。
其次，无论是R1还是R2序列，在测序的时候都是按照5’-3’的方向进行测序的
然后，正负链的主要区别是，正链是我们通常的5’-3’的顺序；负链是3’-5’的顺序。所以，虽然坐标上都是从左边往右数，但是有方向上的区别。
所以，如果这条reads是比对到负链上的，它给出的坐标也是最左边比对的那个坐标。

3、这样的话，我们再理一下这个作者提供的这个代码。

作者在介绍中有说：

这里需要再看看刚刚画的示意图，双端测序中的r2测得的是5’->3’方向的reads，筛选出 比对到reverse链上的r1方向的reads，这些reads对应的就是负链，而r1测序所得的reads是3’->5’方向，与r2方向互补，这不就是reads在负链上的碱基序列么？

我觉得可能，在一般的测序中，这种解释是合理的。但是通过观察我们参考手册中的示意图，我们会发现对于我们的ATAC-seq的数据，无论是哪一端，都是从5’-3’的测序方向。

所以，我不同意这个作者在这里提到的选取比对到负链上的reads的方法，但是给我提供了一些借鉴的思路。
对于同样都是5’-3’测序的双端的reads，其连两者在于参考基因组比对的时候的情况应该是这样的？
（1）R1能够比对上，R2的反向互补序列能够比对上（我们来看看在我们的数据中，情况是否是这样的！）。==>正链的情况
（2）R1的反向序列能够比对上，R2能够比对上 ==>负链的情况。

因此如果要进行flag筛选的话，应该是怎样的？
（1）我们目前正链比对的数据有99和163这两种，且都是正链。

• 99=64+32+2+1，在paired序列中的第一条，与该reads成对的reads的反向互补能够比对上，paired reads每个都能正确的比对到参考序列上，该read是成对的reads中的一个。
• 163=128+32+2+1，在paired序列中的第二条，与该reads成对的reads的反向互补能够比对上，paired reads每个都能正确的比对到参考序列上，该read是成对的reads中的一个。

4、这里又出现了一些新的问题：如何区分flag=16和flag=32的情况？

1 ： 代表这个序列采用的是PE双端测序
2： 代表这个序列和参考序列完全匹配，没有插入缺失
4： 代表这个序列没有mapping到参考序列上
8： 代表这个序列的另一端序列没有比对到参考序列上，比如这条序列是R1,它对应的R2端序列没有比对到参考序列上
16：代表这个序列比对到参考序列的负链上
32 ：代表这个序列对应的另一端序列比对到参考序列的负链上
64 ： 代表这个序列是R1端序列， read1;
128 : 代表这个序列是R2端序列，read2；
256： 代表这个序列不是主要的比对，一条序列可能比对到参考序列的多个位置，只有一个是首要的比对位置，其他都是次要的

回到前面的问题。所以163和99的flag其实是合理的，对于正链的比对情况而言。
那么对于负链的reads，我们的flag可以设置的信息是：

• 64+16+2+1=83 在paired序列中的第一条，该reads比对到负链上，paired reads每个都能正确的比对到参考序列上，该read是成对的reads中的一个。
• 128+16+2+1=147在paired序列中的第二条，该reads比对到负链上，paired reads每个都能正确的比对到参考序列上，该read是成对的reads中的一个。

5、另外，虽然咱们两个reads都是5’-3’的顺序，但是对于基因而言，其方向是自由的。其可以位于R1reads上，但是方向是3’-5’的。这个也是我刚刚突然间想到的。也就是说位于负链上的基因的意义仅仅是说，其方向是3’-5’，但是可能位于两条链。似乎这样就可以解释了。

6、我们接下来验证，在我们的数据中，配对的R1和R2的reads是反向互补的关系吗？如果R1的尾巴刚好也是R2的尾巴的互补，那么两者应该是同方向（都是5’）的测序。如果R1的尾巴刚好是R2的头部互补重合，那么两者反向测序（一个5’，一个3’）。

A00928:207:HYLCHDSXY:2:1101:5981:1078   83      chr5    178154515       60      107S44M =       178154515       -44     CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCTTTCCACCCAGAAGACGGCATACGCGATACTGCATCGTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGGTGTGAAAGATACAGAGCAAGGGATGACCCCACTGACCTGCNT      FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF,:,,,,,FF,,F,,:,F,FF,,F,,FF,FFFFF:FFFFFFFFFFFFFFFFFFFFFFFFF:FFFFF:FFFFFFFFFFFFFFFFFFF,FFFFFFF:FFFFFFFFFFFFFFFFFFFFF#F NM:i:1  MD:Z:42G1       MC:Z:47M104S    AS:i:42      XS:i:0  SA:Z:chr2,32916486,+,113S38M,0,0;
A00928:207:HYLCHDSXY:2:1101:5981:1078   163     chr5    178154515       60      47M104S =       178154515       44      GGTGTGAAAGATACAGAGCAAGGGATGACCCCACTGACCTGCGTCTGTCTCTTATACACATCTGACGCTGCCGACGACAGACGCGTCCGTAATGCTGTGATGTGTAGATCTCGGTGGTCGCCGTATCATTAAAAAGGGGGGGGGGGGGGGG      FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:,FFFFFFFFFFFFFFF NM:i:0  MD:Z:47 MC:Z:107S44M    AS:i:47 XS:i:0
A00928:207:HYLCHDSXY:2:1101:8043:1078   99      chr5    95885157        60      151M    =       95885202        196     GNTGAGCTCTGTAACAACTGAAAAGCCCCTGTGACATTTTACCTTTGAGAGTCCTAACACGGTTTGAGTGGAACAGCTGAGAAACAGCATATATATATTTTAACACCTCAAAATAGTTTGAAATGAGCCTCACAGCCTTGTTCAATCTTCA      F#FFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFF,FF:FFFFFF:FFFFFFFFFFFFFF: NM:i:1  MD:Z:1C149      MC:Z:151M       AS:i:149     XS:i:19
A00928:207:HYLCHDSXY:2:1101:8043:1078   147     chr5    95885202        60      151M    =       95885157        -196    TGAGAGTCCTAACACGGTTTGAGTGGAACAGCTGAGAAACAGCATATATATATTTTAACACCTCAAAATAGTTTGAAATGAGCCTCACAGCCTTGTTCAATCTTCAGATTACAAATAACATTGATAGCATCTCCTGTGGCCTTCAGTTAGT      FFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFF NM:i:0  MD:Z:151        MC:Z:151M       AS:i:151     XS:i:0

7、有没有发现我们刚刚推测的flag值刚好是这两者的配对。
所以接下来，究竟是怎样的关系呢？
（这点思考起来挺难的，目前的想法就是看比对到负链是啥样的情况）
（现在先把容易的那几条数据挂上去）

===>问问题可以让我们思路更加的清晰。
===>分点答题可以让思路更加的清晰。

---

之前的思路，如果被打断就很难接续下来了。现在重新开始开一篇新的博客来记录笔记。