RNA-seq基本流程

下圖是一個(gè)大概的RNA-seq基本流程

[圖片上傳失敗...(image-8ca9c-1586161495252)]

把RNA破碎成小片段，然后將RNA轉(zhuǎn)變成一條cDNA郊酒，這一步需要用到反轉(zhuǎn)錄酶 reverse transcriptase (RT) 才能用RNA作為模板合成DNA。

不論是轉(zhuǎn)錄還是反轉(zhuǎn)錄都需要引物需五。通常如果我們要mRNA，那就可以用oligo-dT作為RT的引物轧坎，但是用它有兩個(gè)問題宏邮，第一個(gè)是只能反轉(zhuǎn)錄那些有A尾巴的RNA，第二個(gè)問題是RT不是一個(gè)高度持續(xù)性的聚合酶缸血，可能讓轉(zhuǎn)錄提前發(fā)生終止蜜氨，造成的結(jié)果就是3'端要比5'端reads富集，這樣就會(huì)使得后續(xù)定量分析帶來bias捎泻。

另一種常用的引物稱為隨機(jī)引物飒炎，隨機(jī)引物的好處是沒有A尾巴的諸如ncRNA也被留下了，而且不會(huì)存在明顯的3'端偏差笆豁。但是很多研究也發(fā)現(xiàn)郎汪，所謂的隨機(jī)引物根本就不隨機(jī)赤赊，這也是測序結(jié)果中，通常前6個(gè)堿基的GC含量分布特別不均勻的原因煞赢。這幾個(gè)堿基GC含量均勻很可能不是接頭或者barcode那些東西抛计，其實(shí)是Illumina 測序RT這一步的random hexamer priming 造成的bias，很多人在處理數(shù)據(jù)的時(shí)候會(huì)把這幾個(gè)堿基去掉照筑，其實(shí)很多時(shí)候真多RNA-seq數(shù)據(jù)去不去掉基本什么影響吹截，不過開頭如果有低質(zhì)量的堿基倒是應(yīng)該去掉。

隨后是第二條鏈合成凝危，這一步用是DNA聚合酶波俄，以剛才合成的第一條鏈作為模板。

接下來就是在序列兩端加上接頭媒抠，加接頭一方面是為了讓機(jī)器可以識(shí)別這些序列弟断，把這些序列固定咏花；二是為了讓多個(gè)樣品可以同時(shí)上機(jī)趴生，平攤每個(gè)樣品的測序價(jià)格。雙端測序?yàn)榱俗宺ead從兩邊開始延伸昏翰，也需要在兩端有所需的引物苍匆。下表是常見接頭種類。

[圖片上傳失敗...(image-67ffaa-1586161495252)]

所謂雙端測序棚菊，因?yàn)楹芏鄷r(shí)候read的長度要短于insert浸踩，為了增加覆蓋度于是就想出了從insert兩端同時(shí)測序的辦法。使得測序深度增加的同時(shí)也能夠用來判斷isoform方向统求。

對(duì)于illumina數(shù)據(jù)检碗，有一條5-3的universal adaptor；還有一條是3-5的indexed adatpor码邻，這條引物含有特異常的barcode折剃。需要說明的是，在雙端測序中像屋，如果insert 不是足夠長怕犁，那么R1可能就會(huì)測到R2的引物，同時(shí)R2可能會(huì)測到R1引物的反向互補(bǔ)序列己莺。

大概的意思就是下面兩張圖奏甫。

[圖片上傳失敗...(image-5ff8aa-1586161495252)]

加了接頭以后進(jìn)行PCR的擴(kuò)增。擴(kuò)增后就開始測序凌受，測序的過程如下圖所示阵子。

[圖片上傳失敗...(image-3d82d-1586161495252)]

[圖片上傳失敗...(image-3fb685-1586161495252)]

測序的基本思想是機(jī)器識(shí)別四種堿基發(fā)出的不同顏色的熒光，可以理解為一個(gè)flow cell 立著非常多序列胜蛉，機(jī)器一層一層掃過去挠进，通過識(shí)別熒光而判斷這一層每個(gè)序列的堿基是什么智蝠。

因?yàn)橐粋€(gè)cell密密麻麻的全是熒光信號(hào)，機(jī)器并不是總能把每一個(gè)判斷的非常準(zhǔn)確奈梳，如果某一個(gè)熒光信號(hào)沒有那么清晰杈湾，這個(gè)堿基的測序質(zhì)量就比較低，如下圖攘须。

[圖片上傳失敗...(image-e29313-1586161495252)]

有的時(shí)候漆撞，如果一大片點(diǎn)都是同一種熒光，機(jī)器也可能犯暈于宙，不知道到底哪一個(gè)熒光屬于哪一個(gè)序列浮驳。這種情況尤其是在序列的前幾個(gè)堿基容易發(fā)生。

The sequencing machine uses the first few bases to establish where the cDNA fragments are on the flow cell. If all of the bases in one part of the flow cell are all the same, like 'C', and all show up green in the picture, then the colors will bleed together and it will not be clear where exactly all of the reads are. In contrast, if you have a lot of different colors in a region, it's easier to determine where each one is, even with a little color bleed.

[圖片上傳失敗...(image-604ea7-1586161495252)]

鏈特異性測序

和普通的RNAseq不同捞魁，鏈特異性測序可以保留最初產(chǎn)生RNA的方向至会，普通建庫方式為什么不行呢？因?yàn)閭鹘y(tǒng)建庫方式通過兩個(gè)接頭的ligation把RNA已經(jīng)變成了雙鏈DNA谱俭，最后的文庫中一部被測序的鏈對(duì)應(yīng)正義鏈（sense strand）奉件，一部分被測序的鏈測是反義鏈。

鏈特異性建庫方式有不止一種昆著，對(duì)應(yīng)到不同的軟件又有不同的叫法县貌，下面是幾種稱呼。要記住的是dUTP 測序方式的名字是fr-firstrand凑懂，也是RF煤痕。至于具體的read方向接下來通過更詳細(xì)的IGV截圖說明問題。

[圖片上傳失敗...(image-32a649-1586161495251)]

[圖片上傳失敗...(image-b9bc15-1586161495251)]

鏈特異性建庫方式（以目前最常用的dUTP為例接谨，如下圖所示）首先利用隨機(jī)引物合成RNA的一條cDNA鏈摆碉，在合成第二條鏈的時(shí)候用dUTP代替dTTP，加adaptor后用UDGase處理脓豪，將有U的第二條cDNA降解掉巷帝。

[圖片上傳失敗...(image-7bbd5-1586161495251)]

[圖片上傳失敗...(image-3eedb-1586161495251)]

這樣最后的insert DNA fragment都是來自于第一條cDNA，也就是dUTP叫fr-firststrand的原因跑揉。對(duì)于dUTP數(shù)據(jù)锅睛，tophat的參數(shù)應(yīng)該為 –library-type fr-firststrand。這里的first-strand cDNA可不是RNA strand历谍，在使用htseq-count 時(shí)现拒，真正的正義鏈應(yīng)該是使用參數(shù) -s reverse 得到的結(jié)果。

正正反反不清楚

說到鏈特異性測序望侈，實(shí)在讓人困惑的是各種鏈的概念印蔬，尤其是翻譯成中文就更說不清了。

DNA 的正鏈和負(fù)鏈脱衙，就是那兩條反向互補(bǔ)的鏈侥猬。參考基因組給出的那個(gè)鏈就是所謂的正鏈（forward）例驹，另一條鏈?zhǔn)欠存湥╮everse）。但是這正反一定不能和正義鏈（sense strand）反義鏈（antisense strand）混淆退唠，兩條互補(bǔ)的DNA鏈其中一條攜帶編碼蛋白質(zhì)信息的鏈稱為正義鏈鹃锈，另一條與之互補(bǔ)的稱為反義鏈。但是攜帶編碼信息的正義鏈不是模板瞧预，只是因?yàn)樗男蛄泻蚏NA相同屎债，正義鏈也是編碼鏈。而反義鏈雖然和RNA反向互補(bǔ)垢油，但它可是真正給RNA當(dāng)模板的鏈盆驹，因此反義鏈也是模板鏈。

總結(jié)兩點(diǎn)

1. 正義鏈（sense strand）= 編碼鏈（coding strand）= 非模板鏈

2. forward strand 上可以同時(shí)有sense strand 和 antisense strand滩愁。因?yàn)檫@完全是兩個(gè)不同的概念躯喇。

寫這篇文章的原因，就是因?yàn)橛腥藛栁蚁跬鳎溙禺愋詼y序數(shù)據(jù) htseq-count 的結(jié)果是不是應(yīng)該把正負(fù)鏈的基因分別在-s yes 和-s reverse 兩個(gè)參數(shù)結(jié)果中統(tǒng)計(jì)出來再做下游分析廉丽。這里犯的錯(cuò)誤就是我們混淆了基因組正反鏈和基因正義反義鏈的概念。

dUTP到底是怎么回事

從前文的一個(gè)圖我們可以總結(jié)出dUTP方式測序R1文件中read1 的方向和基因的方向（正義鏈）是相反的檀咙，而R2文件中的read2 方向和基因的方向是相同的雅倒。

可以參考下面的兩個(gè)igv文件bam截圖璃诀。

首先解釋一下igv 兩個(gè)顏色參數(shù)的意義

• Read strand in pastels, red for positive rightward (5' to 3') DNA strand, blue for negative leftward (reverse-complement) DNA strand, and grey for unpaired mate, mate not mapped, or otherwise unknown status.

• First-of-pair strand assignment is dependent on RNA transcript directionality and is useful for directional libraries. Displays reads or read pairs in which the forward read is first (F1 or F1R2) in red and reads or read pairs in which the reverse read is first (R1 or R1F2) in blue. Unknown status is in gray.

• For a given transcript, non-directional libraries will show a mix of red and blue reads aligning to the locus.

• Directional libraries will show reads of one color in the direction matching the transcript orientation.

下面這個(gè)圖示按照igv 顏色選項(xiàng)中的read strand 方向進(jìn)行區(qū)分弧可，可以看到所有紅色read都是在正鏈方向（注意正鏈不是正義鏈），而所有藍(lán)色的read都是負(fù)鏈方向劣欢。下面基因的方向是正鏈方向棕诵，也就是和紅色的read同向的，如果你把鼠標(biāo)放到隨意一個(gè)紅色的read上凿将，就能看到顯示的信息是second of pair校套，也就是pair中的read2（R2）；反之如果你在藍(lán)色的read上面牧抵，就會(huì)顯示信息是first of pair笛匙，也就是R1 。

總結(jié)犀变，dUTP測序中pair read 中的read1（R1）和基因方向相反妹孙，read2（R2）和基因方向相同

[圖片上傳失敗...(image-d68a15-1586161495251)]

再看下面這張圖

[圖片上傳失敗...(image-37643e-1586161495251)]

這張圖展示了兩個(gè)基因1和2，我們可以發(fā)現(xiàn)gene1的正義鏈就在正鏈上获枝，而gene2的正義鏈其實(shí)是在反鏈上蠢正。看read情況省店，a嚣崭，c兩個(gè)read雖然針對(duì)正鏈負(fù)鏈而言方向一致笨触，都是負(fù)鏈方向，但是如果把a是pair中的read1（first of pair ）雹舀，而c是pair中的read2（second of pair）芦劣。也就是說，read方向一致说榆，但一個(gè)是read1一個(gè)是read2持寄，說明這兩個(gè)read對(duì)應(yīng)的基因一定是反向的。同樣的道理娱俺，雖然b稍味，d都是兩個(gè)方向?yàn)樨?fù)鏈的read，但是b其實(shí)是所在pair的read2（second of pair）荠卷，而d是所在pair的read1（first of pair）模庐。

再次強(qiáng)調(diào)，dUTP測序中pair read 中的read1（R1）和基因方向相反油宜，read2和基因方向相同

當(dāng)使用read strand來進(jìn)行顏色區(qū)分的時(shí)候掂碱，每一個(gè)基因上兩種顏色的分布應(yīng)該相對(duì)均勻（也就是所謂的pair end）。

如果這個(gè)時(shí)候把顏色選項(xiàng)改為按照 first of pair of strand來區(qū)分慎冤，會(huì)出現(xiàn)下圖的變化疼燥。

[圖片上傳失敗...(image-4b98f7-1586161495250)]

gene1的read全部變成了藍(lán)色，而gene2的read全部變成了紅色蚁堤。

如果是非鏈特異性測序醉者，在 first of pair of strand模式下，同一個(gè)gene相關(guān)的read顏色也是明顯混雜的披诗。如下圖

[圖片上傳失敗...(image-e0fd7-1586161495250)]

再一次總結(jié)：

• dUTP 鏈特異性測序中撬即，RNA 方向（gff文件中基因的方向）與read1相反呈队，與read2相同剥槐。如果read1比對(duì)到基因組正鏈上，則對(duì)應(yīng)的gene在基因組負(fù)鏈宪摧；如果read2比對(duì)到基因組正鏈則對(duì)應(yīng)的gene在基因組正鏈粒竖。

• dTUP 測序方式叫做fr-firstrand（留下的是cDNA第一條鏈），也是RF几于。

• 如果dUTP鏈特異性測序蕊苗，看基因表達(dá)量應(yīng)該 counts for the 2nd read strand aligned with RNA(htseq-count option -s reverse, STAR ReadsPerGene.out.tab column 3 )

• 如果想看反義鏈?zhǔn)欠裼修D(zhuǎn)錄本（比如NAT）應(yīng)該用 the 1st read strand aligned with RNA ( htseq-count option -s yes，STAR ReadsPerGene.out.tab column 4)

幾個(gè)常用軟件的設(shè)置

STAR mpping 時(shí)無需特別設(shè)置孩革，但如果不是鏈特異性數(shù)據(jù)且下游分析要用到cufflinks 則需要增加參數(shù) --outSAMstrandField intronMotif岁歉。為的是增加一個(gè)XS標(biāo)簽。

If you have stranded RNA-seq data, you do not need to use any specific STAR options. Instead, you need to run Cufflinks with the library option --library-type options. For example, cufflinks... --library-type fr-firststrand should be used for the standard dUTP protocol, including Illumina’s stranded Tru-Seq.

hisat2 --rna-strandness RF

目的也是給比對(duì)序列添加一個(gè)XS標(biāo)簽以區(qū)分方向，方便cufflinks使用锅移。

For single-end reads, use F or R. 'F' means a read corresponds to a transcript. 'R' means a read corresponds to the reverse complemented counterpart of a transcript. For paired-end reads, use either FR or RF. With this option being used, every read alignment will have an XS attribute tag: '+' means a read belongs to a transcript on '+' strand of genome. '-' means a read belongs to a transcript on '-' strand of genome.

tophat --library-type option fr-firststrand

具體解釋參見下表

[圖片上傳失敗...(image-c810b9-1586161495249)]

htseq-count -s reverse/yes(看反義鏈)

For stranded=no, a read is considered overlapping with a feature regardless of whether it is mapped to the same or the opposite strand as the feature. For stranded=yes and single-end reads, the read has to be mapped to the same strand as the feature. For paired-end reads, the first read has to be on the same strand and the second read on the opposite strand. For stranded=reverse, these rules are reversed.

RSEM --forward-prob 0（正義鏈）1（看反義鏈）

The RNA-Seq protocol used to generate the reads is strand specific, i.e., all (upstream) reads are derived from the forward strand. This option is equivalent to --forward-prob=1.0. With this option set, if RSEM runs the Bowtie/Bowtie 2 aligner, the '--norc' Bowtie/Bowtie 2 option will be used, which disables alignment to the reverse strand of transcripts. (Default: off)

Probability of generating a read from the forward strand of a transcript. Set to 1 for a strand-specific protocol where all (upstream) reads are derived from the forward strand, 0 for a strand-specific protocol where all (upstream) read are derived from the reverse strand, or 0.5 for a non-strand-specific protocol. (Default: 0.5)

sXpress --rf-stranded / --fr-stranded(看反義鏈)

--fr eXpress only accepts alignments (single-end or paired) where the first (or only) read is aligned to the forward target sequence and the second read is aligned to the reverse-complemented target sequence. In directional sequencing, this is equivalent to second-strand only.

--rf eXpress only accepts alignments (single-end or paired) where the first (or only) read is aligned to the reverse-completemented target sequence and the second read is aligned to the forward target sequence. In directional sequencing, this is equivalent to first-strand only.

trinity --SSlibtype RF**

Trinity performs best with strand-specific data, in which case sense and antisense transcripts can be resolved.

RF: first read (/1) of fragment pair is sequenced as anti-sense (reverse(R)), and second read (/2) is in the sense strand (forward(F)); typical of the dUTP/UDG sequencing method.

參數(shù)錯(cuò)了又怎樣熔掺？

到這里，會(huì)想問兩個(gè)問題非剃。有時(shí)候我們不知道數(shù)據(jù)的建庫方式是不是鏈特異性的置逻，如果弄錯(cuò)了結(jié)果會(huì)怎么樣呢？

如果你用STAR mapping 完可以用igv像上文提到的那樣备绽，看看是不是鏈特異性測序券坞。

下面是兩個(gè)真是數(shù)據(jù)的count 統(tǒng)計(jì)情況。

對(duì)于僅僅進(jìn)行基因表達(dá)定量來說肺素，把鏈特異性數(shù)據(jù)當(dāng)作普通建庫數(shù)據(jù)來處理恨锚，可以觀察第2列數(shù)據(jù)和第4列數(shù)據(jù)。具體某一個(gè)基因而言倍靡，影響不會(huì)太大猴伶，因?yàn)榻^大多數(shù)反義鏈本身表達(dá)量就非常低。 不過可以注意 noFeature 和 ambiguous 這兩個(gè)值塌西，因?yàn)榛蚪M中存在兩個(gè)基因分別在正鏈和負(fù)鏈且又重疊的情況他挎，不區(qū)分方向會(huì)比區(qū)分方向的ambiguous數(shù)目多一些。因?yàn)槿绻荒芡ㄟ^方向來區(qū)分到底屬于哪個(gè)基因捡需，這樣的read就會(huì)被認(rèn)為是ambiguous办桨。

但是因?yàn)閰^(qū)分了方向，又會(huì)使得noFeature的數(shù)目更多一些站辉。不過兩者總體影響不會(huì)差別非常大呢撞。如果不能判斷建庫方式，在htseq中使用-s no 參數(shù)是一個(gè)比較保險(xiǎn)（雖然不是非常精確）的做法庵寞。

[圖片上傳失敗...(image-471701-1586161495249)]

相反狸相，如果把普通建庫方式的數(shù)據(jù)當(dāng)作鏈特異性數(shù)據(jù)來處理。

比如在htseq-count中使用了-s reverse 參數(shù)捐川，這個(gè)時(shí)候只有R2方向和基因方向相同的pair才用來算作一個(gè)count，所有R2和基因方向不同的pair就被當(dāng)作no feature了逸尖。這樣的結(jié)果影響可以通過下面的表格觀察古沥。

用正常方法數(shù)出的noFeature 是6萬左右，而用-s yes或者reverse數(shù)出來的noFeature 就接近46萬了娇跟。將近40萬的read 被丟掉岩齿。

所以，如果把普通建庫的數(shù)據(jù)誤當(dāng)作鏈特異性數(shù)據(jù)來處理極有可能會(huì)損失大量的數(shù)據(jù)苞俘，如果弄錯(cuò)了鏈特異性建庫的方式盹沈，那坑你就沒幾個(gè)read剩下了。另外吃谣，計(jì)算出來的結(jié)果自然也會(huì)有非常大的差異乞封，是不準(zhǔn)確的做裙。

[圖片上傳失敗...(image-3ca037-1586161495249)]

作者：思考問題的熊

幾年前分享過一些考研經(jīng)歷和備考建議而被部分考研人所知，也分享過一些英語學(xué)習(xí)相關(guān)的內(nèi)容而被調(diào)侃戲稱為“靠譜學(xué)長”肃晚。

目前在中科院上海植物生理生態(tài)研究所碩博連讀锚贱，努力在生物信息學(xué)菜鳥晉級(jí)的路上前進(jìn)。

[圖片上傳失敗...(image-c4af6d-1586161495249)]

參考資料和更多文章可以移步作者博客 http://kaopubear.top

<label style="overflow-wrap: break-word; cursor: pointer; font-family: Tahoma, Helvetica, SimSun, sans-serif, Hei; font-style: normal; font-variant-ligatures: normal; font-variant-caps: normal; font-weight: 400; letter-spacing: normal; orphans: 2; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: 2; word-spacing: 0px; -webkit-text-stroke-width: 0px; background-color: rgb(255, 255, 255); text-decoration-style: initial; text-decoration-color: initial; font-size: 13px; color: rgb(133, 15, 15);">轉(zhuǎn)載本文請(qǐng)聯(lián)系原作者獲取授權(quán)关串，同時(shí)請(qǐng)注明本文來自王程揚(yáng)科學(xué)網(wǎng)博客拧廊。
鏈接地址：</label>http://blog.sciencenet.cn/blog-3372875-1090305.html