FinCloud-
下图是一个大概的RNA-seq基本流程
[图片上传失败...(image-8ca9c-1586161495252)]
把RNA破碎成小片段，然后将RNA转变成一条cDNA，这一步需要用到反转录酶 reverse transcriptase (RT) 才能用RNA作为模板合成DNA。
不论是转录还是反转录都需要引物。通常如果我们要mRNA，那就可以用oligo-dT作为RT的引物，但是用它有两个问题，第一个是只能反转录那些有A尾巴的RNA，第二个问题是RT不是一个高度持续性的聚合酶，可能让转录提前发生终止，造成的结果就是3"端要比5"端reads富集，这样就会使得后续定量分析带来bias。
另一种常用的引物称为 随机引物 ，随机引物的好处是没有A尾巴的诸如ncRNA也被留下了，而且不会存在明显的3"端偏差。但是很多研究也发现，所谓的随机引物根本就不随机， 这也是测序结果中，通常前6个碱基的GC含量分布特别不均匀的原因 。这几个碱基GC含量均匀很可能不是接头或者barcode那些东西，其实是Illumina 测序RT这一步的random hexamer priming 造成的bias，很多人在处理数据的时候会把这几个碱基去掉，其实很多时候真多RNA-seq数据去不去掉基本什么影响，不过开头如果有低质量的碱基倒是应该去掉。
随后是第二条链合成，这一步用是DNA聚合酶，以刚才合成的第一条链作为模板。
接下来就是在序列两端加上接头，加接头一方面是为了让机器可以识别这些序列，把这些序列固定；二是为了让多个样品可以同时上机，平摊每个样品的测序价格。双端测序为了让read从两边开始延伸，也需要在两端有所需的引物。下表是常见接头种类。
[图片上传失败...(image-67ffaa-1586161495252)]
所谓双端测序，因为很多时候read的长度要短于insert，为了增加覆盖度于是就想出了从insert两端同时测序的办法。使得测序深度增加的同时也能够用来判断isoform方向。
对于illumina数据，有一条5-3的universal adaptor；还有一条是3-5的indexed adatpor，这条引物含有特异常的barcode。需要说明的是，在双端测序中，如果insert 不是足够长，那么R1可能就会测到R2的引物，同时R2可能会测到R1引物的反向互补序列。
大概的意思就是下面两张图。
[图片上传失败...(image-5ff8aa-1586161495252)]
加了接头以后进行PCR的扩增。扩增后就开始测序，测序的过程如下图所示。
[图片上传失败...(image-3d82d-1586161495252)]
[图片上传失败...(image-3fb685-1586161495252)]
测序的基本思想是机器识别四种碱基发出的不同颜色的荧光，可以理解为一个flow cell 立着非常多序列，机器一层一层扫过去，通过识别荧光而判断这一层每个序列的碱基是什么。
因为一个cell密密麻麻的全是荧光信号，机器并不是总能把每一个判断的非常准确，如果某一个荧光信号没有那么清晰，这个碱基的测序质量就比较低，如下图。
[图片上传失败...(image-e29313-1586161495252)]
有的时候，如果一大片点都是同一种荧光，机器也可能犯晕，不知道到底哪一个荧光属于哪一个序列。这种情况尤其是在序列的前几个碱基容易发生。
[图片上传失败...(image-604ea7-1586161495252)]
和普通的RNAseq不同，链特异性测序可以保留最初产生RNA的方向，普通建库方式为什么不行呢？因为传统建库方式通过两个接头的ligation把RNA已经变成了双链DNA，最后的文库中一部被测序的链对应正义链（sense strand），一部分被测序的链测是反义链。
链特异性建库方式有不止一种，对应到不同的软件又有不同的叫法，下面是几种称呼。 要记住的是dUTP 测序方式的名字是fr-firstrand，也是RF。 至于具体的read方向接下来通过更详细的IGV截图说明问题。
[图片上传失败...(image-32a649-1586161495251)]
[图片上传失败...(image-b9bc15-1586161495251)]
链特异性建库方式（以目前最常用的dUTP为例，如下图所示）首先利用随机引物合成RNA的一条cDNA链，在合成第二条链的时候用dUTP代替dTTP，加adaptor后用UDGase处理，将有U的第二条cDNA降解掉。
[图片上传失败...(image-7bbd5-1586161495251)]
[图片上传失败...(image-3eedb-1586161495251)]
这样 最后的insert DNA fragment都是来自于第一条cDNA，也就是dUTP叫fr-firststrand的原因 。对于dUTP数据， tophat 的参数应该为 –library-type fr-firststrand 。这里的first-strand cDNA可不是RNA strand，在使用 htseq-count 时，真正的正义链应该是使用参数 -s reverse 得到的结果。
说到链特异性测序，实在让人困惑的是 各种链的概念 ，尤其是翻译成中文就更说不清了。
DNA 的正链和负链，就是那两条反向互补的链。参考基因组给出的那个链就是所谓的正链（forward），另一条链是反链（reverse）。但是这正反一定 不能和正义链（sense strand）反义链（antisense strand）混淆 ，两条互补的DNA链其中一条携带编码蛋白质信息的链称为正义链，另一条与之互补的称为反义链。但是携带编码信息的正义链不是模板，只是因为它的序列和RNA相同，正义链也是编码链。而反义链虽然和RNA反向互补，但它可是真正给RNA当模板的链，因此反义链也是模板链。
总结两点
写这篇文章的原因，就是因为有人问我，链特异性测序数据 htseq-count 的结果是不是应该把正负链的基因分别在-s yes 和-s reverse 两个参数结果中统计出来再做下游分析。这里犯的错误就是我们混淆了基因组正反链和基因正义反义链的概念。
从前文的一个图我们可以总结出 dUTP方式测序R1文件中read1 的方向和基因的方向（正义链）是相反的，而R2文件中的read2 方向和基因的方向是相同的。
可以参考下面的两个igv文件bam截图。
首先解释一下igv 两个颜色参数的意义
下面这个图示按照igv 颜色选项中的read strand 方向进行区分，可以看到所有 红色read都是在正链方向 （注意正链不是正义链），而所有 蓝色的read都是负链方向 。下面基因的方向是正链方向，也就是和红色的read同向的，如果你把鼠标放到随意一个红色的read上，就能看到显示的信息是 second of pair ，也就是 pair中的read2（R2） ；反之如果你在蓝色的read上面，就会显示信息是first of pair，也就是R1 。
总结，dUTP测序中pair read 中的read1（R1）和基因方向相反，read2（R2）和基因方向相同
[图片上传失败...(image-d68a15-1586161495251)]
再看下面这张图
[图片上传失败...(image-37643e-1586161495251)]
这张图展示了两个基因1和2，我们可以发现 gene1的正义链就在正链上 ，而 gene2的正义链其实是在反链上 。看read情况，a，c两个read虽然针对正链负链而言方向一致，都是负链方向，但是如果把 a是pair中的read1（first of pair ） ，而 c是pair中的read2（second of pair） 。也就是说， read方向一致，但一个是read1一个是read2，说明这两个read对应的基因一定是反向的。 同样的道理，虽然 b，d都是两个方向为负链的read，但是b其实是所在pair的read2（second of pair），而d是所在pair的read1（first of pair）。
再次强调，dUTP测序中pair read 中的read1（R1）和基因方向相反，read2和基因方向相同
当使用read strand来进行颜色区分的时候，每一个基因上两种颜色的分布应该相对均匀（也就是所谓的pair end）。
如果这个时候把 颜色选项改为按照 first of pair of strand 来区分 ，会出现下图的变化。
[图片上传失败...(image-4b98f7-1586161495250)]
gene1的read全部变成了蓝色，而gene2的read全部变成了红色。
如果是非链特异性测序，在 first of pair of strand 模式下，同一个gene相关的read颜色也是明显混杂的。如下图
[图片上传失败...(image-e0fd7-1586161495250)]
再一次总结：
STAR mpping 时无需特别设置，但如果不是链特异性数据且下游分析要用到cufflinks 则需要增加参数 --outSAMstrandField intronMotif。为的是增加一个XS标签。
hisat2 --rna-strandness RF
目的也是给比对序列添加一个XS标签以区分方向，方便cufflinks使用。
tophat --library-type option fr-firststrand
具体解释参见下表
[图片上传失败...(image-c810b9-1586161495249)]
htseq-count -s reverse/yes(看反义链)
RSEM --forward-prob 0（正义链）1（看反义链）
sXpress --rf-stranded / --fr-stranded(看反义链)
trinity --SS libtype RF**
到这里，会想问两个问题。有时候我们不知道数据的建库方式是不是链特异性的，如果弄错了结果会怎么样呢？
如果你用STAR mapping 完可以用igv像上文提到的那样，看看是不是链特异性测序。
下面是两个真是数据的count 统计情况。
对于 仅仅进行基因表达定量 来说，把链特异性数据当作普通建库数据来处理，可以观察第2列数据和第4列数据。具体某一个基因而言，影响不会太大，因为绝大多数反义链本身表达量就非常低。 不过可以注意 noFeature 和 ambiguous 这两个值，因为基因组中存在两个基因分别在正链和负链且又重叠的情况，不区分方向会比区分方向的ambiguous数目多一些。因为如果不能通过方向来区分到底属于哪个基因，这样的read就会被认为是ambiguous。
但是因为区分了方向，又会使得noFeature的数目更多一些。不过两者总体影响不会差别非常大。如果不能判断建库方式，在htseq中使用-s no 参数是一个比较保险（虽然不是非常精确）的做法。
[图片上传失败...(image-471701-1586161495249)]
相反，如果把普通建库方式的数据当作链特异性数据来处理。
比如在htseq-count中使用了-s reverse 参数，这个时候 只有R2方向和基因方向相同的pair才用来算作一个count ，所有R2和基因方向不同的pair就被当作no feature了。这样的结果影响可以通过下面的表格观察。
用正常方法数出的noFeature 是6万左右，而用-s yes或者reverse数出来的noFeature 就接近46万了。将近40万的read 被丢掉。
所以，如果把普通建库的数据误当作链特异性数据来处理极有可能会损失大量的数据，如果 弄错了链特异性建库的方式 ，那坑你就没几个read剩下了。另外，计算出来的结果自然也会有非常大的差异，是不准确的。
[图片上传失败...(image-3ca037-1586161495249)]
几年前分享过一些考研经历和备考建议而被部分考研人所知，也分享过一些英语学习相关的内容而被调侃戏称为“靠谱学长”。
目前在中科院上海植物生理生态研究所硕博连读，努力在生物信息学菜鸟晋级的路上前进。
[图片上传失败...(image-c4af6d-1586161495249)]
参考资料和更多文章可以移步作者博客 http://kaopubear.top
分享好知识！