37、鏈特異建庫

參考文獻：http://www.jianshu.com/p/a63595a41bed

http://kaopubear.top/categories/                                      ####生物信息和編程

RNA-seq基本流程

下圖是一個大概的RNA-seq基本流程

 

 

把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從兩邊開始延伸，也需要在兩端有所需的引物。下表是常見接頭種類。

所謂雙端測序，因為很多時候read的長度要短於insert，為了增加覆蓋度於是就想出了從insert兩端同時測序的辦法。使得測序深度增加的同時也能夠用來判斷isoform方向。

對於illumina數據，有一條5-3的universal adaptor；還有一條是3-5的indexed adatpor，這條引物含有特意的barcode。需要說明的是，在雙端測序中，如果insert 不是足夠長，那么R1可能就會測到R2的引物，同時R2可能會測到R1引物的反向互補序列。

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

 

 

加了接頭以后進行PCR的擴增。擴增后就開始測序，測序的過程如下圖所示。

 

 

測序的基本思想是機器識別四種鹼基發出的不同顏色的熒光，可以理解為一個flow cell 立着非常多序列，機器一層一層掃過去，通過識別熒光而判斷這一層每個序列的鹼基是什么。

因為一個cell密密麻麻的全是熒光信號，機器並不是總能把每一個判斷的非常准確，如果某一個熒光信號沒有那么清晰，這個鹼基的測序質量就比較低，如下圖。

 

 

有的時候，如果一大片點都是同一種熒光，機器也可能犯暈，不知道到底哪一個熒光屬於哪一個序列。這種情況尤其是在序列的前幾個鹼基容易發生。

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.

 

 

鏈特異性測序

和普通的RNAseq不同，鏈特異性測序可以保留最初產生RNA的方向，普通建庫方式為什么不行呢？因為傳統建庫方式通過兩個接頭的ligation把RNA已經變成了雙鏈DNA，最后的文庫中一部被測序的鏈對應正義鏈（sense strand），一部分被測序的鏈測是反義鏈。

鏈特異性建庫方式有不止一種，對應到不同的軟件又有不同的叫法，下面是幾種稱呼。要記住的是dUTP 測序方式的名字是fr-firstrand，也是RF。至於具體的read方向接下來通過更詳細的IGV截圖說明問題。

 

 

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

 

 

這樣最后的insert DNA fragment都是來自於第一條cDNA，也就是dUTP叫fr-firststrand的原因。對於dUTP數據，tophat的參數應該為 –library-type fr-firststrand。這里的first-strand cDNA可不是RNA strand，在使用htseq-count 時，真正的正義鏈應該是使用參數 -s reverse 得到的結果。

正正反反不清楚

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

DNA 的正鏈和負鏈，就是那兩條反向互補的鏈。參考基因組給出的那個鏈就是所謂的正鏈（forword），另一條鏈是反鏈（reverse）。但是這正反一定不能和正義鏈（sense strand）反義鏈（antisense strand）混淆，兩條互補的DNA鏈其中一條攜帶編碼蛋白質信息的鏈稱為正義鏈，另一條與之互補的稱為反義鏈。但是攜帶編碼信息的正義鏈不是模板，只是因為它的序列和RNA相同，正義鏈也是編碼鏈。而反義鏈雖然和RNA反向互補，但它可是真正給RNA當模板的鏈，因此反義鏈也是模板鏈。

總結兩點

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

2. forword strand 上可以同時有sense strand 和 antisense strand。因為這完全是兩個不同的概念。

寫這篇文章的原因，就是因為有人問我，鏈特異性測序數據 htseq-count 的結果是不是應該把正負鏈的基因分別在-s yes 和-s reverse 兩個參數結果中統計出來再做下游分析。這里犯的錯誤就是我們混淆了基因組正反鏈和基因正義反義鏈的概念。

dUTP到底是怎么回事

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

可以參考下面的兩個igv文件bam截圖。

首先解釋一下igv 兩個顏色參數的意義

• 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.

下面這個圖示按照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）和基因方向相同

 

 

再看下面這張圖

 

 

這張圖展示了兩個基因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來區分，會出現下圖的變化。

 

 

geng1的read全部變成了紫色，而gene2的read全部變成了粉色。

如果是非鏈特異性測序，在 first of pair of strand模式下，同一個gene相關的read顏色也是明顯混雜的。如下圖

 

 

再一次總結：

• dUTP 鏈特異性測序中，RNA 方向（gff文件中基因的方向）與read1相反，與read2相同。如果read1比對到基因組正鏈上，則對應的gene在基因組負鏈；如果read2比對到基因組正鏈則對應的gene在基因組正鏈。

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

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

• 如果想看反義鏈是否有轉錄本（比如NAT）應該用 the 1st read strand aligned with RNA ( htseq-count option -s yes，STAR ReadsPerGene.out.tab column 4)

幾個常用軟件的設置

STAR mpping 時無需特別設置，但如果不是鏈特異性數據且下游分析要用到cufflinks 則需要增加參數 --outSAMstrandField intronMotif。為的是增加一個XS標簽。

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

目的也是給比對序列添加一個XS標簽以區分方向，方面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

具體解釋參見下表

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.

參數錯了又怎樣？

到這里，會想問兩個問題。有時候我們不知道數據的建庫方式是不是鏈特異性的，如果弄錯了結果會怎么樣呢？

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

下面是兩個真是數據的count 統計情況。

對於僅僅進行基因表達定量來說，把鏈特異性數據當作普通建庫數據來處理，可以觀察第2列數據和第4列數據。具體某一個基因而言，影響不會太大，因為絕大多數反義鏈本身表達量就非常低。 不過可以注意 noFeature 和 ambiguous 這兩個值，因為基因組中存在兩個基因分別在正鏈和負鏈且又重疊的情況，不區分方向會比區分方向的ambiguous數目多一些。因為如果不能通過方向來區分到底屬於哪個基因，這樣的read就會被認為是ambiguous。

但是因為區分了方向，又會使得noFeature的數目更多一些。不過兩者總體影響不會差別非常大。如果不能判斷建庫方式，在htseq中使用-s no 參數是一個比較保險（雖然不是非常精確）的做法。

相反，如果把普通建庫方式的數據當作鏈特異性數據來處理。

比如在htseq-count中使用了-s reverse 參數，這個時候只有R2方向和基因方向相同的pair才用來算作一個count，所有R2和基因方向不同的pair就被當作no feature了。這樣的結果影響可以通過下面的表格觀察。

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

所以，如果把普通建庫的數據誤當作鏈特異性數據來處理極有可能會損失大量的數據，如果弄錯了鏈特異性建庫的方式，那坑你就沒幾個read剩下了。另外，計算出來的結果自然也會有非常大的差異，是不准確的。

 

作者：生信技能樹
鏈接：http://www.jianshu.com/p/a63595a41bed
來源：簡書
著作權歸作者所有。商業轉載請聯系作者獲得授權，非商業轉載請注明出處。