RNAseq: Normalization

Han Y, Gao S, Muegge K, Zhang W, Zhou B. (2015) Advanced Applications of RNA Sequencing and Challenges. Bioinform Biol Insights 9(Suppl 1):29-46. [article]

計算好基因的表達量,之後的data normalization其實是很重要的一環,這一步驟要考慮的細節很多,像是transcripts size、GC content、sequencing depth、sequencing error rate、insert size。

在評估各種標準化的方式,可以使用measurement error model來評估其標準化的效力。

目前有很多標準化的方式,像是quantile normalization,這方法在處理microarray資料上很有效,在RNAseq中能提高data quality,即使是low amount的RNA。

在R語言中的EDASeq,其先使用within-lane normalization再來between-lane normalization,可以有效降低GC-content所造成的問題

Lowess normalization的方法在microRNA 的資料處理蠻有效的。

目前在處理RNAseq資料上的normalization依舊需要有更好的方法被發展出來。

進一步閱讀:基因的長度是否會影響RNAseq中的表現差異分析呢?

RNAseq: Reads counting

Han Y, Gao S, Muegge K, Zhang W, Zhou B. (2015) Advanced Applications of RNA Sequencing and Challenges. Bioinform Biol Insights 9(Suppl 1):29-46. [article]

在RNAseq資料裡面,常使用對應到某個基因序列的reads量來代表其基因表達量,當然要直接把reads數量對應成基因表現量,這邊一定要仔細考慮在定序前檢體處理的步驟都有相關,比如是否是strand-specific,這些都要列入考慮的因素!

這邊有一些工具可以計算reads counting的量,如Bedtools,給定其GFF和bam檔,其可以在特定的區段計算出counts,其內建是會計算兩個strands上有興趣區域的reads一起計算,但這能調整。

另一款工具HTseq則是有把strands-specific概念在預設計算中,其只會把reads對應到特定的strands上。


在R語言裡,easyRNAseq使用上簡易上手、summarizeOverlap為GenomicRanges裡頭的韓式、featureCount (Rsubread)裡頭則採用了非常有效率的chromosome hashing和feature blocking 的方式。

除了使用工具不同,計算出來的count值也會不同,所使用的gene model 同樣影響到後續的計算,尤其是在跨過junction的reads上,有論文發現在使用RefGene、Ensembl、UCSC三種版本對照序列,發現其中21958個常見基因只有16.3%的基因其算出的count相同,其中9.3%基因有大於50%的差異。

RNAseq: Reads mapping

Han Y, Gao S, Muegge K, Zhang W, Zhou B. (2015) Advanced Applications of RNA Sequencing and Challenges. Bioinform Biol Insights 9(Suppl 1):29-46. [article]

在將RNAseq 定序出來的reads經過preprocessing後,接下來的步驟便是要將其組成的contigs對應到參考序列上,如何將短片段的序列對應到正確的位置一直以來都是生物資訊領域重要的問題!目前已經有很多不錯的軟體能處理這類問題,像是ELAND、SOAP、SOAP2、MAQ、Bowtie、BWA、ZOOM、STAR。這裡有一篇比較這些工具的論文可以參考看看。

Bao S, Jiang R, Kwan W, Wang B, Ma X, Song YQ. Evaluation of next-generation sequencing software in mapping and assembly. J Hum Genet. 2011;56(6):406–14

這些工具在處理一般序列區域的功能都相當不錯,但遇到像是poly(A) tail或是exon-intron splicing junction就必須要針對splicing處理的工具像是BLAT、TopHat、GEM、MapSplice。

在對應回參考序列另一個問題就是序列的polymorphism,這會造成multiple-aligned reads的問題,在較短的reads於序列重複性較高的區域這問題會更嚴重,而在一些重複性較小的reads上,可以利用附近的reads來解決這類型的問題。

Cloonan N, Forrest AR, Kolle G, et al. Stem cell transcriptome profiling via massive-scale mRNA sequencing. Nat Methods. 2008;5(7):613–9

Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008;5(7):621–8

另一種解決方式就是將paired-end reads 延升成200-500bp的長度,在進行mapping。

Holt RA, Jones SJ. The new paradigm of flow cell sequencing. Genome Res. 2008; 18(6):839–46.
Hillier LW, Marth GT, Quinlan AR, et al. Whole-genome sequencing and variant discovery in C. elegans. Nat Methods. 2008;5(2):183–8.
Campbell PJ, Stephens PJ, Pleasance ED, et al. Identification of somatically acquired rearrangements in cancer using genome-wide massively parallel paired- end sequencing. Nat Genet. 2008;40(6):722–9

RNAseq’s 資料前處理:Quality Control

RNAseq資料處理的步驟非常繁複,而第一步就是要看定序出來的品質如何,品質不好的話是否要去掉差的定序reads,特意設計疊加的reads是否要先處理完後再繼續下一步的aligment、是否有被污染、過多的adaptors(比較短的library中容易碰到的問題)。

最常見用來評估定序品質的工具為FastQC

FASTQC 是一款用來分析RNA-seq FASTQ 檔案的軟體,可以在command line環境下執行.

分析後會取得的資料:

  • 檔案名稱:最原始分析前之檔案名稱
  • 檔案類型
  • 解碼(encoding):其品質數值所表示的ASCII
  • 全部定序數量(Total Sequence):總共有多少定序片段
  • 過濾定序(filtered sequence):在Casava模式下,在分析過程中片段序列被標誌且移除
  • 定序長度(Sequence Length):最短和最長的定序序列
  • GC比例(GC content)

FastQC 執行完,會給出11張圖:

Per Base Sequence Quality 每個base平均品質
screenshot.png

    1. 分析內容
      1. 可以用來看是否有特別低品質的定序序列,通常都是因為照向品質不好所造成的,像是位在定序晶片邊緣的點
      2. 達不同品質指標,圖標會顯示之圖例
        • 警告(Warming):當平均品質小於27/0.2% error rate
        • 失敗(Failure):當平均品質小於20/1% error rate
      3. 備註:經過quality trimming後可能會有偏差,或是看不出error batch?

Per Sequence Quality Scores 每條序列中不同位點的品質
screenshot.png

    1. 分析內容
      1. 顯示序列裡不同位置平均品質範圍(只要是fastq中之資料)
      2. BoxWhisker type Plot圖中標示的解釋:
        • 中間紅色線:中間值(median)
        • 黃色長條:四分位距內的值(25-75%)
        • 上下whiskers: 各代表10%和90%的點
        • 藍色線:代表平均值(mean)
        • 背景色
          1. 綠色very good quality
          2. 橘色reasonable quality
          3. 紅色poor quality
        • 標題上是FastQC猜測此定序所使用的方式,在虛列品質好時比較容易猜測錯誤.
      3. 達不同品質指標,圖標會顯示之圖例
        • 警告(Warming):當25%之任何序列位點<10或是median值<25
        • 失敗(Failure):當25%之任何序列位點<5或是median值<20
      4. 問題導因:
        • 原因一: 跑定序的時間過長,造成後面定序出來的品質下降
          1. 處理方式: Quality Trimming(依據平均品質來決定裁掉的長度)
          2. 通常在adapter read-through也會有問題
        • 原因二:在一開始的時候就產生問題,可能有泡泡跑過定序的夾層,這時候就不適合trim,會因此拋棄過多好品質的序列
        • 原因三:太低的coverage可能會造成false negative

Per Base Sequence Content序列中不同位置其ATCG所佔的比例

screenshot.png

    1. 分析內容:
      1. 在一個隨機的library中,其序列中不同base應該不會有太大的差別.
      2. 通常在開頭的地方容易有不平均的ATCG比例,通常在Library製備過程中priming使用到random hexamers(幾乎所有RNA-seq的製備都會使用),或是在fragmentation過程中使用transposases,通常會在reads開頭有偏差出現,使得在有enrichment在5’端的different k-mers處,這種技術性造成的偏差無法藉由trimming消掉,但在downstream analysis中也不會造成太大的問題
      3. 達不同品質指標,圖標會顯示之圖例
        • 警告(Warming):當任何一處的ATCG其差異超過10%
        • 失敗(Failure):當任何一處的ATCG其差異超過20%
      4. 問題導因:
        • 原因一:Overrepresented sequence
          1. 可能是做定序所使用的adapter dimers或是rRNA
        • 原因二:Biases Fragmentation
          1. 因製備Library中所使用的random hexamers在前12bases會有selected biase,但這在後續的process可以處理掉,不太會影響後續的分析
        • 原因三:Biases Composition Libraries
          1. 大多數的library本身就會有些composition的偏差,尤其是使用sodium bisulphite製備的,其會造成cytosine變成thymines,所以會使得C的比值偏高
        • 原因四:Extreme Trimming
          1. 在尾端會出現,因為很積極地將跟adapter match的序列去除掉

Per Sequence GC Content 測量GC在序列中的比例
screenshot.png

    1. 分析內容:
      1. 其主要是分析不同序列的GC之分布圖,應該要跟理論分布差不多
      2. 比較陡或是尖:可能是被adapter汙染
      3. 比較寬:可能是被不同物種所污染
    2. 達不同品質指標,圖標會顯示之圖例:
      1. 警告:偏離理論值15%
      2. 失敗:偏離理論值30%

Per Base N Content 平均每個位點壞品質(N)

screenshot.png

    1. 分析內容:
      1. 當某個序列即位點定序品質很差,則會用N代表
      2. 通常當某個部分其N所佔的比例上升的話,其會使後面的分析無法獲得有效的base cal.
    2. 達不同品質指標,圖標會顯示之圖例:
      1. 警告:任一個位置其N所佔比例大於5%
      2. 失敗:任一個位置其N所佔比例大於20%
    3. 備註原因:
      1. 原因一:
        • 通常是整體定序品質不佳,可以檢查一下特定bin的coverage,可能有某一個bin其定到的序列非常少
      2. 原因二:
        • 有嚴重偏差的LIBRARY,在Per Base Sequence Content也可以看出來

 

   Sequence Length Distribution 定序長度分布圖
screenshot.png

    1. 分析內容:
      1. 顯示平均定序的長度,這部分就要先了解是使用何種定序方式,有的會有單一長度,有的則是不同
      2. 達不同品質指標,圖標會顯示之圖例:
        • 警告:假如定序的長度有不一樣的話
        • 失敗:任一個定序長度是零

Duplicate Sequences 重複序列
screenshot.png

    1. 分析內容:
      1. 在一個均勻的library中序列應該只會出現一次,很低程度的重複率可能暗示高coverage,高度序列重複則可能是library中有某種程度的enrichment bias
      2. 為了使分析的計算效率提高,這個分析模組只會分析在前100000的序列,另外在超過10個重複的序列,其會被放在一起呈現
      3. 比較長的定序讀數,通常會造成一定程度的低估重複率
      4. 藍色線:代表完整序列其重複率, 紅色線:代表去掉重複後不同序列的比例
  1. Overrepresented Sequences
  2. Adapter Content
  3. Kmer Content
  4. Per Tile Sequence Quality

去掉品質差Reads的工具Trimmomatic

可以用來處理illumina機器裡的adaptor,其也可以將low-quality的base或是N base去掉,其可以處理single-end和paired-end的資料。此程式是java base的軟體,在command line下操作。

 

處理特殊設計的paired-end 實驗設計:BBMerge FLASH

有些實驗設計時,會將paired-end read設計成疊在一起的樣子,而這兩款程式可以將重疊的paired-end reads合成單一條較長的reads,且能計算出innert-size的分布圖。

偷偷查了一下FLASH的使用趨勢,似乎暴增人數使用!
screenshot.png

RNAseq可用來研究的主題

Han Y, Gao S, Muegge K, Zhang W, Zhou B. (2015) Advanced Applications of RNA Sequencing and Challenges. Bioinform Biol Insights 9(Suppl 1):29-46. [article]

RNAseq的資料為high resolution到single base,可應用範圍寬廣

RNAseq的技術可以用來區分許多不同種類的transcripts像是long non-coding RNA、miRNA、siRNA、snRNA、piRNA,以及他們如何調控蛋白質的translation或是調控chromatin states

RNAseq也可以用來研究enhancer region的transcripts表現量,這可能跟epigenetic gene regulation相關,可以稱這類型的transcripts為enhancer RNA

從RNAseq的資料可以取得transciption start sites、alternative promoter usage、mRNA isoforms (alternative splicing)、mRNA 穩定性(premature transciption termination at the 3′ )

Locating regulatory element

Allele-specific expression、disease-associated single nucleotide polymorphisms(SNP)、fusion gene for causal variants in cancer

Transcription of endogenous retrotransposons and parasitic repeat elements that influence the transcription of neighboring genes or result in somatic mosaicism in the brain

論文閱讀:Differential gene and transcript expression analysis of RNA-seq experiments with Tophat and Cufflinks

原始論文:

Cole Trapnell,Adam Roberts,Loyal Goff,Geo Pertea,Daehwan Kim,David R Kelley,Harold Pimentel,Steven L Salzberg,John L Rinn& Lior Pachter. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks.Nature Protocols.7,562–578(2012)Doi:10.1038/nprot.2012.016

 

RNAseq的資料分析主要有兩個大目的:

  • 找尋新的基因或是transcripts
  • 比較不同狀態下的基因表現亮

而在分析步驟上,可以把它拆成三個分析流程:

  • read alignment
  • transcript assembly 或是 genome annotation
  • transcript and gene quantification

本篇論文主要介紹用期開發的程式作為workflow,並且簡單介紹每個tool之間的異同。

所使用的工具和工作流程:

Read alignment with Tophat
Reads拼回去reference的步驟非常重要,他可以找出所有可能的insertion/deletion/indel,這些資訊可以讓我們了解辨識出相對於對照序列的polymorphism,另外,對不回去的reads,可能是新的protein-coding genes或是noncoding RNA。這邊拼回去的結果,也會影響之後如何去計算transcript abundance的數值。

螢幕快照 2016-01-14 下午3.50.21

 

Transcripts assembly with Cufflinks
為了要知道各種基因的表現量,我們必須知道各種reads所對應的isoforms transcripts,這邊就必須把Reads拼組成各種特定的transcipts,這步驟牽涉到如何辨別不同種transcripts variant ,所以最後產生transfrags來表達所有可能。

cufflinks會將一些low abundance的tansfrag當作是來自於immature transcripts而丟棄,這部分是要了解的。(或許這步驟可以思考一下!),有時候我們會將所有的bam pool在一起讓他找transcripts,但這種作法其實一方面增加電腦計算的loading,另一方面,則是會讓程式在計算isoform時增加很多變數,所以cuffmerge應映兒生,他可以把individual產出的transfrag merge在一起看!

螢幕快照 2016-01-14 下午4.14.47

Differeial Expression Analysis with Cuffdiff

cuffdiff可以輸出好幾個檔案描述不同基因和transcripts的表現量差別和p value以及其名字和在基因體上的位置,另外,cuffdiff也可以把同一組基因裡面不同TSS的transcripts分組去計算他們之間比例的變化!

 

視覺化資料CummeRbound

cuffdiff提供基因和transcripts的表現量分析,且這些資料是以tab-delimited的方式輸出很方便後續的分析,而CummeRbound則是提供一個更容易操作的視覺化分析,可以將cuffdiff的計算結果做更細緻的分析。

 

替代工具:

  • Read-alignment 
    • GSNAP
    • START
    • Map-Splice
  • Transcriptome reconstruction
    • De novo transcripts assembly
  • Quantification
  • Differential expression

有趣的相關論文:

Read-alignment

Transcript reconstruction

Quantification

Differential Expression

 

論文閱讀:Optimizing and benchmarking de novo transcriptome sequencing: from library preparation to assembly evaluation

RNAseq的分析目前已經非常火紅的工具,相對於非模式生物的研究,要從頭把reference genome組起來,不如只接從transcriptome做起,在經費有限下是不錯的選擇,而如何最佳化這樣的分析是此論文的重點。

此篇論文主要有使用的分析軟體有:eggNOC、Trinity、SOAPdenovo-trans,其主要是分析蜥蜴胚胎(Madagascar ground gecko)三個發育時期的轉錄體,為了提高其在de novo assembled transcript sequences的完整性,其使用vertebrate one-to-one orthologs來作為reference。

其方法是基於調整RNA library、read lengths和insert sizes,加上使用有233個同源基因的對照組(one-to-one orthologs 來自29種species),最後使用CEGMA和BUSCO來執行completeness assessment,展現了此分析方式有效地提升的精準度。

螢幕快照 2016-01-10 上午11.54.56螢幕快照 2016-01-10 上午11.55.21

Hara Y, Tatsumi K, Yoshida M, Kajikawa E, Kiyonari H, Kuraku S. (2015) Optimizing and benchmarking de novo transcriptome sequencing: from library preparation to assembly evaluation. BMC Genomics 16(1):977. [article]