RNAseq基因表現量常用指標:RPKM, FPKM,TPM

在RNAseq的分析中,如何將比對到特定基因範圍內的reads數量轉換成“這段基因的基因表現量”,一直是RNAseq分析的起頭,但也是還沒有定論的部分,這邊紀錄最一開始用來轉換比對到特定基因範圍內的read數量到“此基因表現量”的方式,那就是RPKM,緊接者則是FPKM,和2012年開始提出的TPM,這三種指標某種程度來說,觀念類似,主要都有考慮到基因長度和總reads數量。

文獻回顧
第一篇提出用RPKM來代表某個基因的表現量的論文:
Mortazavi, A., Williams, B.A., McCue, K., Schaeffer, L., Wold, B.(2008).Mapping and quantifying mammalian transcriptomes by RNA-seq. Nature Methods. 5, 621-628

這一篇在2012年提出來TPM代表轉錄體表現量的論文,則有詳細對於RPKM, FPKM, TPM的公式解說
Gunter P. Wagner, Koryu Kin, Vincent J. Lynch.(2012). Measurement of mRNA abundance using RNA-seq data: RPKM measure is inconsistent among samples. Theory in Biosciences. 131(4), 281-285

基本上,在RNAseq中,實際上每一個sample的mRNA總量可以想成:
mRNA_total = \sum{mRNA_g}
而以下三種指標其實都是建立在這看法下,考慮到mRNA長度、總reads數、transcript數等來調整成代表基因表現量的數值

RPKM(Reads Per Kilobase per Million)
RPKM = \frac{10^6 * n_r}{L * N}
n_r 代表單一個gene其上的reads數量,L代表是這段基因的總長度除1000(轉換成kilobase的單位),N則是總共的reads數,通常用於singel reads

FPKM(Fragments Per Kilobase per Million)
FPKM = \frac{10^6 * n_f}{ L * N}
同上,這邊的fragments其實是代表一組paired-end reads所捕捉到的transcripts,所以其實適合用來在paired read的實驗中使用

TPM(Transcripts Per Million)
TPM = \frac{n_r * rl * 10^6}{fl_g * T}
T = \sum{\frac{r_g * rl}{fl_g}}
這邊的rl代表read lengths,然後fl_g 則代表是這基因的exon總長度

看完這三者的公式便能理解為何會說RPKM, FPKM , TPM都是一種樣本內的標準化方式(within normalization),它代表的是單一個gene或是transcript相對於所有樣本的表現總量。而當RNAseq的實驗設計較為複雜,要用來比較樣本間(between samples)各個基因表現量差異時,就會多少造成一些低表現向的基因之偏差,因為單一基因在不同組中比較時,用其mRNA的長度來做調整變顯多於,也引入偏差。

閱讀參考
What the FPKM? A review of RNA-Seq expression units
Somnath Datta, Dan Nettleton.(2016).Statistical Analysis of the next generation sequencing data. Frontiers in Probability and the Statistical Sciences, Springer
基因组学技术专题(二)—— 为什么说FPKM/RPKM是错的

R的網絡視覺化包整理: igraph, diagrammR , ggraph

整理一下R包中做網絡視覺化(network visualization)的工具包

igraph

簡介
有十年開發歷史的network graph包,其發展成熟,很多後來開發的工具包
開發者:Gábor Csárdi, Tamás Nepusz

GGally

簡介
為ggplot2的延伸套件,其中的函數ggnet2專門可以用來畫network
開發者:Barret Schloerke

network3d

簡介
是D3 javascript的network graph R套件,也支持igraph的資料結構,他可以將圖形輸出至rmarkdown, shiny或是網頁的形式。
開發者:Christopher Gandrud, JJ Allaire, Kent Russell, & CJ Yetman

Rgraphviz

簡介
主要由bioconductor項目的負責人所開發的工具包,但主要語法承襲base繪圖語法,且視覺化的效果跟plot系列的調性一樣
開發者:Jeff Gentry, Robert Gentleman, Wolfgang Huber

diagrammR

簡介
這個工具的設計很有野心,整合多種常見的network創建語法,包括graphviz等,其整套語法建構得很完整,加上文檔講解清晰,算是目前network visualization系列的工具包中文檔最詳盡用心的。
開發者: Richard Iannone, Kent Russell , JJ Allaire, Michaël Benesty

network

簡介
使用network object來將關聯性資料轉換成graph型態
開發者:Carter Butts

ggnetwork

簡介
可以將network object直接用ggplot來繪圖
開發者:François Briatte

ggraph

簡介
是目前很完整地實現ggplot2語法來設計graph, network視覺化的工具包
開發者:Thomas Lin Pedersen

tidygraph

簡介
承襲ggraph,但著眼與如何用tidyverse的語法來處理graph/network data structure
開發者:Thomas Lin Pedersen

總結

各種包的差異上,主要是在graphics、graph data manipulation、layout三大方面的差異,graphics效果上可以大致分成三大類,承襲base plot的效果、ggplot grammer和javascript功能三類,graph data manipulation上來說,diagrammR和ggraph在圖形資料結構的處理上都支持關聯性資料、matrix資料等模式,目前只對igraph, diagrammR這兩個包有比較深入的使用,在視覺化的效果上,後者的確優秀許多。

R視覺化套件:grid簡介

grid套件是由Paul Murrell所寫的,而目前知名的ggplot2便是奠基於grid的工作上,相對於ggplot2的設計,grid可以讓使用者從比較底層的方式去控制視覺化物件,且最新版的R已經屬於base package,所以一安裝R就會安裝好grid。

想要了解更深入grid背後的邏輯和設計可以看Paul Murrell所寫的書R Graphics,有提供免費的網頁版本可以閱讀。

grid架構視覺化圖形的邏輯主要有兩個:
1. 直接底層控制畫線、正方形、文字、圓形、點,使用x-y座標控制
2. 藉由控制視野變化來調整視覺化效果(Viewports)
screenshot.png

PBS(Portable Batch System)系統效調細節筆記:使用叢集係統來做百萬級抗原變異點的預測

之前有測試幾種加速運算的方法,用python裡的線程搭配外掛程式來躲避GIL的效能限制等,但不同種類的架構能應付的問題(高I/O, 高運算…),這邊嘗試使用手邊另一個叢集系統來做計算,利用PBS中分發任務的能力來調用底下的電腦群們來做運算,過程也是一堆的大坑,主要的框架是用shell script寫成的,因為手邊這個叢集電腦python是舊版的,然後還得確定下面的電腦們每個環境都要乾淨,所以就先試者用最簡單的方式來加速計算,最後效調的結果可以在一周內分析完3425328個蛋白序列,雖然還是不夠快。(不確定是否是最佳做法,可能偏“土炮”)

整個過程有幾個重點:
1. 使用shell指令來處理大檔案的技巧
2. PBS系統的理解

第一部分在shell script的撰寫上,相對於高階語言如python,就少掉很多如線城池、共享變數等等用來暫時存放運算資料的方法,這時候就會需要比較複雜的代碼來實現想要的邏輯。技術上可以先把script用最小單位來撰寫,接者在測試成功後,慢慢合成一個大的,因為當牽涉越多層的呼叫使用,需多未知的東西狀態出現,尤其是想要用&的方式來加速(因為目前這叢集系統無法安裝linux parell指令)。

這邊可以堆間使用awk和sed的功能來做不雜一點的計算和資料搬移(他們很輕且通常很快),尤其搭配awk裡面的條件語法和函數,計算速度可以因此稍微提生一點,且能用“一行”來解決在高階語法可能要很多行代碼才搞定的字串處理問題。

第二部分關於PBS語法和功能的理解,雖然每個叢集系統多少會不同(就跟每個linux管理員管理的想法不同一樣,會對自己的系統加入不同的觀念),但指令的說明是差不多的。

三大關於PBS的指令分別為:pbsnodes, qstat, qsub

pbsnodes: 用來看目前叢集系統中,各個節點的資源和任務條波狀況

screenshot.pngqstat:用來看目前任務的運行狀況和單個任務所使用的資源

screenshot.png

qsub:提交任務到PBS系統上

screenshot.png

 

 

還有意外發現使用parent-child的shell process中,當程式的架構互叫的方式可能會造成無法預期的意外,比如簡單的file descriptor的指定,在parent和child process中分開執行,且一對多的時候運行的邏輯很怪,這部份可能要更深入底層才會理解這部分的坑。

R空間視覺化的學習資源(GIS or Spatial Visualization)

R在空間視覺化上在最近有很多的提升,尤其是整合如Google map api, leaflet等含有多種美麗地圖風格的package,這邊整理一下相關的學習資源。

 

R語言官方針對空間資料的文檔

CRAN Spatial Data Task View
裡面也有各種處理圖資的package可以使用


Roger Bivand. Applied Spatial Data Analysis with R. (2008). Springer
Chris Brunsdon, Lex Comber. An Introduction to R for Spatial Analysis and Mapping. (2015). SAGE

文章
R Spatial Cheatsheet
Zev Ross. Manipulating and mapping US Census data in R using the acs, tigris and leaflet packages.

Gene Transcription Regulation Database(GTRD):轉錄因子調控資料庫

在今年2017的NAR資料庫特輯中,GTRD資料庫整理人類和老鼠ChiP-seq資料並且彙整相關轉錄因子調控,使用四種負責處理ChiP-seq的peak caller程式,分別是MACS, SISSRs, GEM, PICS,並且將這些程式分析出的peak做clustering後,去掉重複的資訊。screenshot.png

目前使用起來不是很順手,每個要查看的轉錄因子都要一個個調整來看,有點不方便,似乎是使用biouml作為架構的,所以整體操作上頗有侷限性,當要查詢多個轉錄因子時會很麻煩,不過算是一個可以當作備用的轉錄因子調控資料庫,且是相對較新的,裡面用來估計peak左右可能影響的基因數量主要跟據來自於HOCOMOCO資料庫的資訊,這資料庫是同一個研究團隊在2013年的發表。screenshot.png

 

他本身可以直接下載它們整理好的interval檔案,資料欄位整理的很乾淨,並且提供那四種peak caller程式分析的資料可以挑選,或是最終其整合的版本。

screenshot.png

 

在進行變異辨認(Variant calling)前的序列品質校正(Base Quality Recalibration)流程

這邊主要是在記錄使用GATK中的RNAseq Variant Calling流程中很關鍵的一步,便是校正位點的品質資訊(Recalibrate base quality scores),這點為何重要呢?因為本質上我們就是在區分資料中單一位點是否有差異,其中最大的“障礙”,就是定序過程中的錯誤所造成的偽變異,如何去估算並且調整每一個位點下的品質數值是否是真的,就會決定我們變異辨認結果的準確程度。簡單來說,可以把Recalibrate Base Quality Score想成是Variant Calling的關鍵資料前處理步驟

關於GATK在Base Quality Recalibration的部分有詳盡的影片可以參考

實務面上重要的關鍵點:

理解在Read Group中會被用來矯正的是以lane和library為單位來思考,所以在進行Read Group的標誌的時候,會決定這一步的Recalibrate到底正不正確
(可以參考這邊GATK的文章)

在做Quality Score Recalibration主要有幾個步驟:
1. Add/Replace Read Group
把Read 資料貼上正確的Read Group(有的人會在alignment的時候就貼上,也可以在alignment後的bam檔來進行處理)

2. Markduplicates/CreatIndex(Picard)/Split’n’Trim(GATK) (調整Read alignmnet和建索引)

  1. 使用BaseRecalibrator(GATK),計算出矯正Quality Score的表

  2. 使用PrintReads(GATK),將上一步計算出的表來調整原本的Bam file,並且輸出矯正過後的Bam檔

閱讀參考:
1. GATK官方文檔和說明
https://software.broadinstitute.org/gatk/documentation/article.php?id=3891
http://gatkforums.broadinstitute.org/gatk/discussion/2801/howto-recalibrate-base-quality-scores-run-bqsr

VCF(Variant Call Format) 基因突變資料儲存格式

次世代定序技術出來後,緊接者是許多高通量資料的處理問題,從一開始單純看定序品質到後面用定序資料組成的基因體後,來看其變異資訊,此時就需要一個很好的儲存方式,這時候就有了VCF(Variant Call Format),這格式主要是由1000 Genomes Project所推動產生的,這計畫的目的希望用定序不同種族的正常人類來探討族群基因體間的異同。

一份VCF檔案大致長得像下面這樣
screenshot.png

基本上分成兩大部分:VCF header (描述所使用的工具、篩選標準、紀錄的細節名稱)和 body (變異資料本身)。screenshot.png

變異資訊基本上會有9行,分別紀錄CHR染色體號碼、POS位置、ID變異編號(是否在dbSNP中)、REF此位點在參考基因上的序列、ALT變異的序列、QUAL此序列的定序品質、FILTER是否通過所設定的篩選條件、INFO進一步的資訊、FORMAT格式資訊、SAMPLES各樣本在此位點的狀況。

會很驚訝於這樣的表示方法,其實就可以包含各種複雜的基因變異形式,如下圖所表現的樣子
screenshot.png

Reactome資料庫,可愛的視覺化小圖標(Reactome Icon Library)

之前就有稍微介紹一下Reactome,最近幾乎跟他形影不離,發現他不斷在提供一些超棒的小資源,其中就是可愛的生物小圖標,最近在練習畫複雜網絡的視覺化,他便提供了很好的素材。這邊下面的圖就是使用Cytoscape來視覺化我的項目的部分成果,整體空間的架構頗平衡的。
screenshot.png
這圖片中的背後的細胞圖標就是來自於Reactome Icon Library,頗佛心來著的。
screenshot.pngscreenshot.pngscreenshot.png