写在前面的废话
经过这两周意外的断更，我明白了一个道理，人都是具有惯性的……
两周没写公众号，我似乎都习惯了懒惰，文章不想更了，笔记也懒得记了。
但是我没忍住双十一的诱惑，一不小心剁了手。于是生活告诉我，我再不努力码字就只能在接下来的日子里吃土了。
image.png
太长不看系列
链特异性建库拆分正负链：
- 单端测序
samtools view -b -f 16 test.bam > test.rev.bam
samtools view -b -F 16 test.bam > test.fwd.bam
- 双端测序
###### 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
废话超多系列
上面的东西看起来可能有点绕，我第一次看也有点懵。但是对照着samtools的flags，仔细思考一下，就会发现原来是这么个玩意啊。
双端测序的逻辑
首先需要明白双端测序的逻辑，说到这里你可能要笑了，这么简单的事情还需要讨论么？我想说的是，的确需要，我当时就是这里没弄明白，整个人蒙圈了。
所以这里如果不说清楚，跳跃幅度过大很容易影响思路。
image.png
我们知道，双端测序通常会产生两个文件，一个我们称之为read1，另一个称之为read2，测序结果会分别放到不同的fastq文件中。我们可以从fastq文件的注释行，看到这两个测序文件的区别，read1的注释行（以 @开头的行）末尾会有\1，而read2的注释行结尾则是\2
我们知道第一个文件是对应的测序片段的forward链的信息，那么第二个文件的碱基序列呢？它的碱基序列对应的是read1的反向互补链信息呢？还是单纯的反向序列呢？
我们这里画图展示一下
### 片段太短，双端测序时两个read有重合
------[f-180bp]----->
----[r2]----->
<----[r1]----
### 片段较长，双端测序时两个read没有重合
------------------[f-350bp]----------------->
----[r2]----->
<----[r1]----
从这个图里，我们知道，测序得到的第一个文件是r1序列而第二个文件是r2序列。那么，这就解答了一个疑惑，r1和r2的序列是反向互补的，不是简单的反向……
另外，我们经常说的read2对应的是我们需要知道的reads的方向，而reads1则对应的是待测reads的反向互补序列
也许你们中的许多人早就知道了，但我是最近查资料才知道 ，我也相信一定有许多人不知道这回事……
讲这个东西有什么用呢？一个是为了给下文做铺垫，另一个就是，如果你测序的reads比较长，比如双端150bp，而建库时的序列比较短，如只有30bp。这种序列是一定会被测通的，此时只需要将r2反向互补一下，并与r1汇集到一起，就从一定程度上增加了测序的深度。
fastx-toolkit其中的一个工具就可以实现这样的功能
image.png
区分正负链
通常情况下，我们是使用sam文件的flags进行区分。在这之前，首先我们需要知道几个flags的意义：
- 16：比对到reverse链上的reads
- 128：双端测序中reads
- 64：双端测序中第一个reads
-f参数：表示包含后面flags值所表示的reads
-F参数：表示不包含后面flags值所表示的reads
这里以单端链特异性测序简单举个例子：
samtools view -f 16 test.bam：只包含（-f参数）比对到reverse链（flags值：16）上的reads
接下就到了本文的重点，如何利用flags区分正负链
链特异性单端测序
比对到负链上的reads：samtools view -f 16 test.bam > rev.bam
比对到正链上的reads：samtools view -F 16 test.bam > fwd.bam
换句话说，就是排除掉那些比对到reverse链上的reads
链特异性双端测序
双端测序略显复杂，但只要把逻辑理顺了，reads区分开 ，也就是分分钟的事情了。
这里需要再看看刚刚画的示意图，双端测序中的r2测得的是5'->3'方向的reads，筛选出 比对到reverse链上的r1方向的reads，这些reads对应的就是负链，而r1测序所得的reads是3'->5'方向，与r2方向互补，这不就是reads在负链上的碱基序列么？
现在让我把刚刚说的话，用samtools的代码表示一下：
# 双端测序中属于r2这类reads的flag值为128，因为要包含这一类 所以使用-f
# 比对到reverse链上的reads的flag值为16，因为要去除这类reads，所以使用-F
# 我们需要得到属于r2且没有比对到负链上的reads
samtools view -b -f 128 -F 16 a.bam > a.fwd1.bam
# 需要属于r1这类的reads，flag=64
# 这类reads与我们想要知道的序列反向互补，因此为了得到比对到正链上的reads，我们需要对筛选出r1这类reads中比对到其对应的reverse strand的reads（flags=16）
# 16+64=80
samtools view -b -f 80 a.bam > a.fwd2.bam
# 将刚刚筛选出来的两部分属于forward链上的reads合并到一一起
$ samtools merge -f fwd.bam a.fwd1.bam a.fwd2.bam
同理，大家可以自己思考一下，如果我想获得比对到reverse链上的reads，我该如何利用samtools筛选得到呢？
如果你真的想不出来，请回到太长不看系列，那里有代码可供参考。。。
一点题外话
这两天查阅资料，验证自己的想法时，看到了思考问题的熊曾经说过的一句话，我觉得很适合警醒自己，在这里送给大家：
很多东西就是这样，你以为的明白并不是真的明白，一年前的明白和一年后的明白也不是同一个明白。我这么说，不知道你能明白还是不明白
网友评论