深入R語言Object-Oriented Programming(二):開發ggplot2延伸套件所使用的ggproto系統(未完)

目前ggplot2基於grid建立了個龐大的繪圖系統,因此當然會需要使用oop來管理日漸複雜的ggplot2 “world",在ggplot2 2.0.0 版提供了官方的文檔在於ggplot2底層的資訊,在Wickham Hadley一開始寫ggplot時,沒有預想到如今如此廣泛的使用和許多套件及開發者的參與,所以在底層物件撰寫上有許多可改進的地方,以幫助跨套件的移植性和讓開發者能使用ggplot2裡面的基礎往上做出更豐富的視覺化。

在ggplot2 2.0.0釋出了ggproto的詳細文檔,而ggproto是基於proto這個oop系統修改而成的“微”面向編成,具有跨套件的繼承性(cross-package inheritance),另外,因為ggplot2為基於grib的繪圖系統,所以到這邊開始需要基礎的grid知識,會比較能融會貫通兩者的關聯性,這邊先專注解釋ggproto裡面的架構,之後會在grid的部分詳加說明。

這部分可能先對R中的OOP系統多少理解一下,會比較能理解他們語法的初衷,可參考前一篇文章

建立ggproto物件

# Construct new_Geom class inheired from Geom class
new_Geom <- ggproto("new_Geom", Geom,
required_aes = c("x", "y"),
defaut_aes = aes(fill="yello"),
draw_key_key = draw_key_polygon,
draw_group = function(data, panel_scales, coord){
coords <- coord$transform(data, panel_scales)
grid::polygonGrob(
x = coords_df$x,
y = coords_df$y,
gp = grid::gpar(col = coords_df$colour,
fill = coords_df$fill,
alpha = coords_df$alpha)
)
}
view raw Geom hosted with ❤ by GitHub
# create ggproto object
# 會發現語法承襲自proto的系統,不需要像是RC system的field等
new_ggproto_object <- ggproto("Class_name",
parent_ggproto_class,
x = 0
f1_method = function(self,n){
self$x <- self$x + n
self$x
}

目前已有ggplot2::Geom, ggplot2::Stat, ggplot2::Position, ggplot2::Scale這四個內建的ggproto類別,而在做延伸ggplot2視覺化的函數時可以用上面四個類別來起手。

Geom類別
– draw_panel(self, data, panel_scales, coord) or draw_group(self, data, panel_scales, coord)
– required_aes
– default_aes

Coord類別
– aspect
– labels
– render_bg
– render_fg
– render_axis_h
– render_axis_v
– range
– train
– transform
– distance
– is_linear

Facet類別
– compute_layer
– map_data
_ draw_panel

Stat類別
– finish_layer
– setup_params
– setup_data
– required_aws
– default_area

ggpubr:提供“可發表品質”的圖表(ggplot2-Based Publication Ready Plots )

最近在看ggplot2/grid底層視覺化的架構,不小心碰到這個頗美的套件ggpubr,標榜“ggplot2-Based Publication Ready Plots ”,他把一些在圖表發表前會遇到的細節處理得頗好,使用了ggrepel來處理標號擁擠,且整體的theme將ggplot2常見會忘記設的參數提供更好的預設,此套件幾乎包含了所有常用的圖表:density ploy, histogram, box plot, violin plot, dot plot… 。推薦他在pie chart上的處理,這邊是原本ggplot2畫出來效果比較差的圖表,另外,提供如cleveland’s dot plots的繪製和scatter plot側邊density的顯示。

這邊是ggpubr套件詳細的內容列表,這整套件可以當作是advance ggplot2用法可以創作出的效果展示!

其中竟然有提供基因資料在比較表現量很重要的MA plot

 

 

探索資料庫應用(三):UALCAN快速查找TCGA基因表現和存活分析的資料庫

screenshot.png

UALCAN是一個可以快速查找TCGA裡面使用RNAseq資料為基礎,提供基因表現量和存活分析的資料庫,是基於PERL-CGI、javascript和css所搭建的(有點酷,他發在neoplasia這本雜誌中的文章,直白地寫使用何種程式語言)。其有四大優點:

1. 方便快速的取得癌症相關的轉錄資料
2. 可以用來快速驗證biomarker的gene
3. 提供表現量和存活分析以及把分析的結果作圖
4. 提供基因資料的超連結到GeneCards、Pubmed、TargetScan、The human protein atlas。

操作介面簡介

提供兩種模式來查找在TCGA的RNAseq資料,可以輸入特定想查的gene,在不同種癌症種類種的表現量

screenshot.png

輸入想查找的基因(可以一次輸入多個基因,使用,相隔即可),結果會提供表現量、存活分析、基本資訊超連結。

screenshot.png

點選進去看分析結果,且這些結果都可以直接下載

screenshot.pngscreenshot.png

論文鏈結:

UALCAN: A Portal for Facilitating Tumor Subgroup Gene Expression and Survival Analyses(2017)

次世代定序資料分析學習資源(NGS analysis learning material)

update: 20200619

整理一些自學次世代定序分析(NGS/Massively Parallel Sequencing)的學習資源:
老實說,門外漢剛入門要做次世代分析真的是很苦的一件事,因為通常這類人都是沒有電腦資訊背景的,但次世代定序的分析所需要的技能其實蠻廣的,不論是統計知識、計算機知識、基因體知識等等,而這領域進步太快,即是到了現在大部分東西還是不懂,建議入門的人可以先參加實體的工作坊一次(花錢撞墻),然後可以看一些網路課程如john hopkin在coursera上開的,同時要準備一台好的電腦盡量實際動手操作,有了自己的研究題目後再慢慢往比較細的主題鑽研

中研院蔡怡陞老師 次世代定序課程
蔡怡陞老師之前在wellcome trust sanger institute工作過,對於次世代定序的資料分析非常有經驗,下面的evolution and genomic workshop也是他部落格裡面推薦的!

哈佛大學生物資訊中心

Evolution and Genomic workshop

加拿大Ontario Institute for Cancer Research 研究中心bioinformatic.ca 研討會資料
這是ontario institute for cancer research中心裡生物資訊團隊的網站,他們每年都會辦許多的工作坊,資料都會公開在網路上,也是頗不錯的學習資源

bioconductor courses
裡面主要都是在R語言中分析次世代定序資料的各種套件使用法,且屬於偏資料後段部分的分析

美國能源部系統生物學知識部門 Kbase


2017. The Biostar Handbook. biostar
很簡單的一本書,之前有稍微介紹過, 適合當工具書查,不適合當學習主軸
2014. RNAseq data analysis: a practical approach. CRC
對於做RNAseq分析的人,是一本好入手的書,雖然部分內容已經開始落後
2014. Statisitcal analysis of next generation sequencing. springer
在統計分析細節講得很細…..假如沒有要很深入數理統計面可以避過

網路MOOC
coursera
John Hopkins: Become a next generation sequencing data scientist(難度易)
1. Introduction to Genomic Technologies
2. Genomic Data Science with Galaxy
3. Python for Genomic Data Science
4. Algorithms for DNA Sequencing
5. Command Line Tools for Genomic Data Science(推)
6. Bioconductor for Genomic Data Science(推)
7. Statistics for Genomic Data Science(推)

edx
Harvard: Biomedical series(難度中部分偏難)
1. Introduction to Bioconductor: Annotation and Analysis of Genomes and Genomic Assays(推)
2. High-performance Computing for Reproducible Genomics(推)
3. Case Studies in Functional Genomics

基因組瀏覽器整理(Genome Browser tools)

screenshot.png
做複雜的次世代定序資料分析,少不了使用基因組瀏覽器來(Genome Browser),基因組瀏覽器是圖形化的工具用來展示生物序列資訊,和不同的基因組特徵(genomic feature)。通常由長條狀的展示框(track)組成,以基因序列的位置作為線性軸,目前並沒有一套標準用來展示各式各樣的基因組特徵,所以頗需要一定的經驗才能上手。

通常基因組瀏覽器可以用來視覺化fasta、bed、gtf/gff、sam/bam檔,且大多數的基因組瀏覽器都有預先放置人類和老鼠的基因組參考序列。

最陽春的方式來瀏覽bam檔可以使用samtools tview的功能,但需要較複雜的探索性分析還是需要基因組瀏覽器的幫助,下面是目前還在active維護的基因組瀏覽器,大多都具備相當不錯的表現能力,幾乎都以java寫的,所以在安裝的時候要注意java相關的依賴,在啟動的時候大多都可以調整記憶體需求,也能利用remote的方式傳送影像。

IGV
上次更新時間:2016/11/1
screenshot.png
由Broad institue所開發,元老級的基因瀏覽器,是目前“市佔率”最大的基因組瀏覽器,在2013年的時候性能有一次大升級。

Artemis
上次更新時間:2017/1/31
由wellcome trust sanger institue所開發,可同時組裡genebank和ensemb格式,且其能手動去調整註釋,另外,也是ensembl內部在組裡資料所使用的軟體

SeqMonk
上次更新時間:2017/6/2
由開發fastqc軟體的Babraham institue所開發,對於許多chipseq, methylation等3d structure資料有支持

iobio
上次更新時間:2017/7

screenshot.png
唯一提供線上基因組資料瀏覽和分析的工具,整體軟體使用起來很舒服,為University of Utah研究團隊Marth lab所開發的軟體,在展示variant時,還會提供各種variant score來做為搜尋disease-causing variant幫助使用

IGB
上次更新時間:2016/6/7
在2001就被Affymetrix開發出來了,直到現在都還在維持中,他的使用手冊是這幾個genome browser中算是寫得最詳盡和清楚的,可以用來瀏覽Chipseq、RNAseq等等輸出資料

tablet
上次更新時間:2016/9/6

tablet1
由The James Hutton Institue所開發的瀏覽器

htmlwidgets: 將javascripts的視覺化功能帶入R中

這幾年Rstudio團隊一系列的努力,其中一個方向便是導入熱門的前端語言如javascript(D3.js, three.js都是很厲害的資料視覺化套件),將其在渲染資料的能力和互動性帶入到R中,雖然想要利用htmlwidge來導入javascript資料視覺化的力量其門檻相對較高,畢竟需要知道javascript的語法和背景知識,目前使用htmlwidge引入javascript資料視覺化到R之相關R套件有networkD3, threejs, rglwidget, DiagrammeR, rbokeh, visNetwork, Plotly。

htmlwidgets提供一套簡便的框架,讓R可以使用javascript library,達到下面四大目的:
1. 將Javascript 資料視覺化的套件能在R console中使用,就如同一般畫R plots的方式
2. 可以將視覺化的圖嵌入在R markdown之中
3. 可以在shiny應用中使用
4. 可以在Rstudio中輸出成單獨的web pages

 

screenshot.png

R包推薦tidyverse:健達出奇蛋般一次滿足所有需求,broom/dplyr/ggplot2/purrr/readr/stringr/…

不得不說,當開始依賴上hadley wickham帝國裡面的語法和封包時,幾乎在處理資料時都會一次load上好幾個“tidyverse”的包,如今hadley wickham終於推出了一個同名專輯同名封包tidyverse一次可以把常用的"tidyverse"包讀進去R裡面使用。

這一方面代表者新時代的來臨,畢竟"tidyverse"風格的語法跟原本的base R有相當大的差異,但“使用者為大”,新風格帶來的是進步和規範,當然減少不少“技術債”,同時,也會刺激大家研讀tidyverse包的開發邏輯,以便讓自己所寫的封包可以被納入tidyverse神殿之中。

screenshot.png

裡頭最核心的包有ggplot2, tibble, tidyr, readr, purrr, dplyr, stringr

個人覺得tidyverse帶給我們除了所謂家庭號分享包的便利外,主要是傳達一整個代碼和分析的“價值觀”和“原則”,擁有原則和價值觀加上努力,長期來說,會讓我們往好的方向前進

  1. 重複使用已經使用的資料結構
  2. 將簡單的函數整合到pipe中
  3. 擁抱函數編程(functional programming)的快感
  4. 為人類服務(寫code讓別人讓自己看得懂xd)

參考閱讀:

Tidyverse:R 語言學習之旅的新起點

The tidyr tools manifesto

The tidyverse style guide

關於fastq的Illumina sequence identifiers(instrutment type information)

在複雜的RNAseq甚至最近的single cell sequencing,在fastq上所顯示的資訊越來越複雜,且很多統計上的方法其實也會使用到這些資訊,所以多少要理解一下而要理解這些資訊就需要知道測序晶片上面的一些命名,雖然不同型號的有所差異,但概念差不多,可以由下圖體會一下感覺

screenshot.png

下面則可以看相關fastq上面每個reads上都會有的編碼,和相關意思的對照

screenshot.pngIllumina sequence identifiers

閱讀參考:
Question: illumina instrument type from fastq?
https://www.biostars.org/p/198143/

FASTQ format
https://en.wikipedia.org/wiki/FASTQ_format

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是错的