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

給生物學家的計算分析(六):shell下建立scripts

閱讀自Practical Computing for biologists by Haddock Dunn

完整建立一個scripts是幫助自己在command line完成許多事情的力氣,就不用再做很多重複的工作。

基本概念:所有在command line底下的指令本就就是一組scripts

所有在command line下的使用指令,本身就是一個scripts,只是他一開始就安裝好了,如ls、cd、rm..等等,我們來看看他們放哪邊!screenshot.png

我們使用which這個指令可以用來找尋物件在系統中的路徑,會發現這些基本指令都放在 /bin 裡面,這資料夾放滿系統各式重要的資料,當我們輸入cd按下時,系統會去尋找這幾個資料夾內的scripts,看有沒有符合的名稱,有的話就會執行!那哪邊可以知道放置這些系統會去尋找的特定資料夾呢?

screenshot.png

使用echo $PATH,就能叫出目前系統當指令輸入後會去尋找的路徑!各個路徑用:隔開。所以我們只要把建好的檔案放進去這些資料夾中,便能像是cd一般的指令來使用,但把自己弄的“土炮”scripts弄髒系統裡面核心功能資料夾,似乎不是個很好的方式,所以另一方式就是建立一個資料夾,然後把其路徑放入$PATH中,似乎是個不錯的方式!那要如何做呢?

以MacOS系統為例,每個帳號都有一個系統隱藏檔案.bash_profile,是來儲存各種個人的設定,(在linux中為.bashrc),可以在母目錄裡輸入 ls -a 看到此資料夾screenshot.png

在裡頭輸入

screenshot.png

那如何將command lines結合成scipts可以參考這篇文章如何將command lines 組合成shell script

癌症基因體閱讀心得(一)

閱讀Genetics: From Genes to Genomes, 5th edition 5th Edition by Leland H. Hartwell (Author), Michael L. Goldberg (Author), Janice A. Fischer (Author), Leroy Hood (Author), Charles F. Aquadro (Author)

 

簡單來說,癌細胞就是突破了一般細胞生長的限制,正常細胞會有以下的機制去制衡其生長:

  1. contact inhibition
    正常細胞彼此有接觸後,就會停止,在cell line中長成monolayer,而癌細胞則不會受限制,會形成所謂的transformed loci
  2. apoptosis
    大部分的細胞在受損後,假如無法修復其基因的改變,就會執行自我凋亡,但癌細胞不會,即使其基因序列已經突變,無法修復,其依舊繼續存活者
  3. autocrine stimulation
    正常細胞在沒有接受外界刺激時,其不會proliferation,但癌細胞會分泌自我刺激的激素來活化自己分化或是複製

細胞的癌化

細胞累積了一定量的突變才轉化為癌細胞,對照流行病學,年齡增長和癌症發作比率同為正相關,這樣的解釋蠻合理的,分子層面來看,這些突變可能是來自於基因中特定的strands,慢慢累積成不同的clones,而癌細胞有兩個特色:prolifereation rate上升加上癌細胞genomic instability,整個為一個惡性循環。

 

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

給生物學家的計算分析(五):shell下的文本處理

閱讀自Practical Computing for biologists by Haddock Dunn

shell環境下的文本處理效率相當高,但需要combine許多shell command line tool,且同樣的處理其實方法會有蠻多種的,可能隨者經驗增加或是使用的工具變多而調整,這邊就稍微介紹一些基本的工具,已可做非常多的應用。

shell環境裡的文本編輯器nano

在command line底下編寫文字,需要一個方便的工具,雖然之前介紹過圖形介面的TextWrangler,但直接在command line裡面編輯有時候比較方便,其實在unix核心下的系統(linux, OS X)有蠻多內建選擇的,如vi、 vim 、emac、pico,每種操作細節上都有一點不同,而這邊介紹nano,一款比較基礎的,可以直接在command line下呼叫 nano
螢幕快照 2016-01-16 上午7.35.37螢幕快照 2016-01-16 上午7.35.48

編輯器下排其實就有非常清楚的解說,關於如何離開、寫入等等,如輸入:

^X   離開
^O   寫入
^Y    向上翻頁
^V    向下翻頁

 

 

給生物學家的計算分析(四):進入shell

閱讀自Practical Computing for biologists by Haddock Dunn

command line的工作環境是一開始生物背景跨進生資分析的障礙之一,相對於已經被圖形介面(Graphical User Line)的電腦環境所嬌生慣養,要說服我們從舒適的圖形介面跳進黑壓壓的指令框框,想必要有強而有力的說詞,這邊就要說說圖形介面的壞處:
  • 生資分析常常要做複雜的步驟來處理資料,或是同樣的步驟操作多次,這些工作要是再圖形介面下不知道要用滑鼠點到什麼時候。
  • 一般在GUI的環境下,大部分的分析都不會留下log黨,讓我們去回顧可能是哪一個步驟做錯了,而在command line下則很方便能把不同分析步驟發生的事情紀錄在log中,等出問題的時候方便除錯
  • 叢集電腦系統不支持圖形化介面的操作,當我們需要使用更加計算能力時,使用圖形介面會是個問題
  • 對於軟體開發者來說,圖形介面的軟體開發複雜且不同平台使用會需要調整,所以在軟體更新上,command line操作環境下的軟體更新的較快!
以上的問題,其實很明顯地想要處理大型資料就必須要克服這個恐懼感。
基礎的unix command 知識:
  • 相對路徑和絕對路徑
  • 相對路徑的使用:
    • ~   代表home目錄
    • ..  代表前一個目錄
    • .   代表目前目錄
  • cd 可以用來在系統路徑中移動
  • 基礎指令:
指令
功能
ls
當下目錄的檔案
pwd
目前目錄之路徑
mkdir
建立資料夾
rmdir
移除資料夾
rm
移除檔案
less
看檔案
man
查看指令說明
cp
複製檔案
mv
移動檔案或是重新命名

給生物學家的計算分析(三):Regular Expression使用

思考如何寫Regular expression的thinking process
以下是資料
螢幕快照 2016-01-16 上午8.01.22
首先思考要如何整理,整理成什麼樣子
螢幕快照 2016-01-16 上午8.04.27
再來看目前資料每行的狀態
螢幕快照 2016-01-16 上午8.06.46
先看哪些是要的資訊,用(     ) 標示
(13) (Jan)uary, (1752) at (13):(53)    (-1.414)    ( 5.781)etc
先看括號裡的資料怎麼用Regular Expression表示
螢幕快照 2016-01-16 上午8.18.53
再將括號外的符號放進去
(13) (Jan)uary, (1752) at (13):(53)    (-1.414)    ( 5.781)etc
(\d+)\s+(\w{3})[\w\,\.]*\s+ (\d+)\sat\s (\d+):(\d+)\s+    ([-\d\.]+)\s+    ( [\d\.]+).*
    1               2                             3               4          5               6                         7
在調整取代裡面的括號資料順序
(13) (Jan)uary, (1752) at (13):(53)    (-1.414)    ( 5.781)etc
(\d+)\s+(\w{3})[\w\,\.]*\s+ (\d+)\sat\s (\d+):(\d+)\s+    ([-\d\.]+)\s+    ( [\d\.]+).*

    \3               \2 \.              \1               \4                 \5              \ 6                         \7

 

*:(quantifiers)至少符合0個
+:(quantifiers)至少符合1個