BioC2016: Where Software and Biology Connect 會議系列三

Bioc2016第三天的會議也相當精彩,早上的科學新知演講開場由Rob Tibshirani史丹佛著名統計學家,他的研究成果非常豐富,從早年的bootstrap到近年來不斷思考一些statistic inference selection的方法學。

 

screenshot.png

 

 

有史丹佛教授William Greenleaf 講解他在分析ATAC-seq資料的過程和背後的思維,ATAC-seq簡單來說就是用定序的技術來“間接”探討chromatin的結構,也就是epigenetics的領域,Greenleaf的背景非常數理,是應用物理學、化學出身,講話給人一種加強版的福爾摩斯,剛好其就是英國人。

 

screenshot.png

史丹佛教授Stephen Montgomery則是早期做GWIS的研究、慢慢開始有eQTL和做allele specific及gene variation間的分析

BioC2016: Where Software and Biology Connect 會議系列二

系列一

Bioconductor第二天開始的會議,就變成白天已scientific talk為主,而下午則是workshop的方式,scientific talk又會分成科學演講(主要由教授來分享)和社群的演講(偏向目前新開發之package使用分享)。

一開始由柏克萊大學Sandrine Dudoit生物統計教授分享主題為Identification of Novel Cell Types in the Brain using Single-Cell Transcriptome sequencing,其本身研究發展許多可用來標準化RNAseq的方法

再來則是史丹佛大學的Susan Holmes(其個人網站有許多之前talk的資料)分享題目為Muticomponent data inegration for the Human Microbiome,其本身研究使用無母數的分析像是SVD、bootstrap、MCMC來探討生物現象所產出的資料,最近也有研究如愛滋病病患使用藥物後其microbiome的改變。

Jenny Bryan(他個人頁面有許多學習資源鏈結)教授則是分享其開發的googlespreadsheet package 其串連goolge doc API去讓使用者可以輕鬆地將google spreadsheet輸入到R中來做進一步的處理。Michael Love給了一個有Bioconductor Workflows Following Fast, Lightweight RNA Transcript Quantifiers的talk,他本身與哈佛生統教授Rafael Irizareay在edx開了一系列biostatistic的課。發展slidify的Ramnath Vaidya則在講新的Html widget。

下午的workshop則教了許多其他的packages像是UC riverside的教授Thomas Girke來發表其所撰寫的systemPipeR,是一款可以直接在R中使用的pipeline工具(雖然部分功能有待加強),而華盛頓大學的James MacDonald則很仔細地探討各種目前常用R database其中qeury等細節的論輪(expressionSet,TxDb,OrgDb…..)

BioC2016: Where Software and Biology Connect 會議系列一

Bioconductor 算是R語言中對於生物資訊最重要的一個project,目前演進成一個龐大的社群,裡面眾多的協力者發展出很多High-throughput genomic data analysis的工具。

Bioc2016:Where Software and Biology Connect 是今年此社群最重要的聚會之一,這次會議在6/25-26號,在美國史丹佛大學舉辦,會議以Bioconductor 裡的專案為主,談論關於高通量基因資料的分析技術和分析工具,每日的上半天議程以演講為主,下午則是實作議程。
第一天為開發者會議,而第二、三天則是主要的議程。
第一天的開發者會議議程的位置在史丹佛的Jame H.clark center,也是Bio-X Program的所在地,Bio-X program裡匯集結合phyisc、biology、computer science的研究,學生也可以申請其fellowship。
symposium
會議內容在連結中有,因為飛機是中午到達San Jose airport,所以只有參加到最後一場由23andme,Robert Gentlemen的talk,關於Distributed Content Generation and Creating Pipelines for Genomic Analysis,這場talk其實在講非常觀念的東西,這問題便是有什麼很好的方式呈現關於一個基因或是關於一段序列的資訊,且能持續update且easily for maintainance。首先,如何搭建一個很好的平台有效呈現genomic data的model,且可以良好的呈現其source difference,另外,如何讓開發團隊的effort可以最小化,也是關鍵之一。
晚上的時候參加由Héctor Corrada Bravo 教授領軍的團隊,用hackthon的方式來探討如何改進epiViz(epiViz主要是一款能視覺化genomic data 分析後的資料),能讓大家用interactive的方式來交流很多圖像化的資源,活動有披薩、啤酒且備有很多紙筆可讓參加者一起腦力激盪,最後再用5min的時間簡單的報告一下,最後討論出來有針對comment authority開放的問題,是否使用location-based的方式呈現資料、如何讓comment內的資料可以被有結構的mining、如何使用將這視覺化的功能用email就能有效的溝通。
 screenshot.png

Mongodb系列一:安裝

MongoDB是開源的document資料庫,在2009年發佈第一個版本,便是NoSQL界的矚目新星,整個資料庫設計的重點便是scalable和data accessibility,以key和value的概念來儲存資料。其為JSON document database(準確講應該是BSON),且允許用ad hoc的方式來讀取巢狀的資料(如JSON資料結構的特色),沒有像是relational database需要先設計好scheme,這當然算是他的特色,相對的就要一些代價要注意。
screenshot.png

基本上,mongodb的資料庫結構如上圖所示,一個database下面會有各種collection,裡面則塞滿document,每個document裡頭的資料以json形式儲存,所以可以很“nested”。

開始摸索mongoDB的第一件事情:安裝

  可以先進入官網裡頭的說明文件非常清楚,假如是用OSX的話,可以直接用hombrew來安裝mongodb,只要直接在terminal打入

$brew install mongodb

安裝好後,要先設定好資料庫的路徑,通常預設會用/data/db,假如沒有的話,在自己下面的指令來興建

$mkdir -p /data/db

記得要把整個資料夾的讀取權限開啟,可以用chmod的指令
接者便可以運行mongod   這是mongoldb的daemon程序,負責外界與資料庫的溝通
$mongod

接者開啟另一個terminal來進入mongodb shell操作整個與資料庫的溝通,基本上只要打入

$mongodb
便可以開始rock n roll了!

screenshot

左邊是開啟的mongod daemon運行的狀態,右邊的terminal則是可以開啟mongodb shell

使用R做Network Analysis系列二:網絡資料結構和處理

network analysis就跟任何種類的分析一樣,大部分的時間會花在“clean data”,也就是把資料處理成適合分析的格式。而在network analysis中,關鍵是節點(nodes)和節點間的關係(edges),而非一般我們用observation為row,variable為column的形式。
有兩種常用的方式,來記錄一個網絡,那便是sociomatrix和edge list。
sociomatrix
  • 將所有的node之間的關係用一個巨大matrix來記錄,row代表是starting node,column代表是receiving node,但當這個network很巨大的時候,此種形式大部分的cell都會是空值。但資料小時候很一目了然。
    screenshot.png
edge list
  • 直接用一個table記錄每個paired的關係,第一個column代表是starting node,第二個column則是receiving node,這種方式比較直覺,在資料量大的時候處理速度也會比sociomatrix建立的網絡還快,因為少去了很多empty data。
    screenshot.png
 
除了單純定義好一個網絡中的node和node之間的關聯,最重要的當然是要把“意義”放進去,就是所謂的屬性(attribute),比如這個網絡是在看臉書好友網絡,而每個node屬性就要被記錄有人、身高、體重等等,而node之間的edge也會有強與弱的分別,是深交或是一般朋友,這些就是要儲存進入網絡資料之中的。
用簡圖來描述R裡面一個R network object要有的資訊
 
種類                                   描述                                                         

Nodes                               網絡中的節點                                                               
Ties                                   節點間的關聯
Node attribute                   節點的屬性
Tie attribute                       網絡中關聯的屬性
Metadata                           整個網絡的相關資訊
在R中創建一個network object

<

div>

#load the library
library(statnet)


#Creat a network object with sociometrix data format
netmat1 <- rbind(c(0,1,1,0,0),
                 c(0,0,1,1,0),
                 c(0,1,0,0,0),
                 c(0,0,0,0,0),
                 c(0,0,1,0,0))
rownames(netmat1) <- c("A","B","C","D","E")
colnames(netmat1) <- c("A","B","C","D","E")
net1 <- network(netmat1,matrix.type="adjacency")
class(net1)
summary(net1)
##Visualization
gplot(net1)
gplot(net1, vertex.col = 5, displaylabels =TRUE)
#Creat network object with edge list format
netmat2 <- rbind(c(1,2),
                 c(1,3),
                 c(2,3),
                 c(2,4),
                 c(3,2),
                 c(5,3))
net2    <- network(netmat2, matrix.type="edgelist")
network.vertex.names(net2) <- c("A","B","C","D","E")
summary(net2)


#Transform between the differnet data structure
as.sociomatrix(net1)
class(as.sociomatrix(net1))

all(as.matrix(net1) == as.sociomatrix(net1))
as.matrix(net1,matrix.type = "edgelist")

 screenshot.png

 

#-Managing Node and Tie attributes

#Node attribute
set.vertex.attribute(net1, "gender", c("F","F","M","F","M"))
net1 %v% "alldeg" <- degree(net1)
list.vertex.attributes(net1)
summary(net1)

get.vertex.attribute(net1,"gender")
net1 %v% "alldeg"

#Tie attribute
list.edge.attributes(net1)
set.edge.attribute(net1,"rndval",
                   runif(network.size(net1),0,1))
list.edge.attributes(net1)
summary(net1 %e% "rndval")
summary(get.edge.attribute(net1,"rndval"))

#example for the value netwrok and modify the edge attribute
netval1 <- rbind(c(0,2,3,0,0),
                 c(0,0,3,1,0),
                 c(0,1,0,0,0),
                 c(0,0,0,0,0),
                 c(0,0,2,0,0))
netval1 <- network(netval1,matrix.type="adjacency",ignore.eval=FALSE, names.eval="like") #use the vertex's vaule added to attribute
list.edge.attributes(netval1)
get.edge.attribute(netval1,"like")

as.sociomatrix(netval1)
as.sociomatrix(netval1,"like")

#--------------------------------------------------------------------
#igraph

#detach the statnet during the duplicate of function name
install.packages("igraph")
detach(package:statnet)
library(igraph)
#create from sociomatrix
inet1  <- graph.adjacency(netmat1)
class(inet1)
summary(inet1)
str(inet1)
#create from edge list
inet2 <- graph.edgelist(netmat2)
summary(inet2)


#modify the attribute
V(inet2)$name <- c("A","B","C","D","E")
E(inet2)$val  <- c(1:6)
summary(inet2)
str(inet2)


#--------------------------------------------------------------------
#Going back and forth between statnet and igraph
install.packages("intergraph")
library(intergraph)
class(net1)
net1igraph <- asIgraph(net1)
class(net1igraph)
str(net1igraph)


#Imporintg Netwrok data
detach("package:igraph",unload=TRUE)
library(statnet)
netmat3 <- rbind(c("A","B"),
                 c("A","C"),
                 c("B","C"),
                 c("B","D"),
                 c("C","B"),
                 c("E","C"))

net.df  <- data.frame(netmat3)
net.df
write.csv(net.df, file = "MyData.csv",
          row.names = FALSE)
net.edge <- read.csv(file = "MyData.csv")
net_import <- network(net.edge,matrix.type = "edgelist")
summary(net_import)

gden(net_import)

#Common Network Data Task

#Filtering Networks Based on Vertex or Edge Attribute Values
n1F <- get.inducedSubgraph(net1,
                           which(net1 %v% "gender" == "F"))

n1F[,]
gplot(n1F,displaylabels = TRUE)

deg <- net1%v% "alldeg"
n2  <- net1 %s% which(deg >1)
gplot(n2,displaylabels = TRUE)

#Removing Isolates
library(UserNetR)
data("ICTS_G10")
gden(ICTS_G10)
length(isolates(ICTS_G10))

n3 <- ICTS_G10
delete.vertices(n3,isolates(n3))
gden(n3)
length(isolates(n3))

#Filtering Based on Edge Values
data(DHHS)
d <- DHHS
gden(d)

op <- par(mar = rep(0,4))
gplot(d,gmode="graph",edge.lwd = d %e% 'collab',
      edge.col = "grey50", vertex.col = "lightblue",
      vertex.cex = 1.0, vertex.sides=20)

as.sociomatrix(d)[ 1:6,1:6]
list.edge.attributes(d)
as.sociomatrix(d,attrname="collab")[1:6,1:6]

table(d %e% "collab")

d.val <- as.sociomatrix(d,attrname="collab")
d.val[d.val < 3] <0
d.filt <- as.network(d.val, directed=FALSE,matrix.type="a",
                     ignore.eval=FALSE,names.eval="collab")
summary(d.filt, print.adj=FALSE)
gden(d.filt)
#Method to drawing
op <- par(mar = rep(0,4))
gplot(d.filt,gmode="graph",displaylabels = TRUE,vertex.col="lightblue",vertex.cex=1.3,
      label.cex=0.4, label.pos=5,
      displayisolates = FALSE)

#Advanced 
op <- par(mar = rep(0,4))
d.value <- as.sociomatrix(d,attrnment="collab")
gplot(d.filt,gmode="graph",thresh=2,vertex.col ="lightblue",vertex.cex=1.3,
      label.cex=0.4, label.pos=5,
      displayisolates = FALSE)
par(op)



#Transformating a directed network to non-directed network
net1mat <- symmetrize(net1, rule = "weak")
net1symm <- network(net1mat,matrix.type="adjacency")
network.vertex.names(net1symm) <- c("A","B","C","D","E")
summary(net1symm)

使用R做Network Analysis系列一:5-number summary for a network

 

使用R做Network Analysis系列二

這一系列的學習是為了之後分析複雜生物網絡所準備,也是系統生物學中很重要的研究技巧,如何做network analysis,這邊用R 來做操作

正常的descriptive statistic也會有所謂的5-number summary來描述一組資料的基本特性(最大值、最小值、中位數、25th位數、75th位數)。

同樣的道理,當我們面對一個網路,要如何了解他的基本性質也有五大特性可以來看

  • Size:
    簡單來說就是一個network裡頭的節點(可稱作nodes, vertices或是actors),這個是基本上可以讓我們知道一個網絡的大小規模。
  • Density:
    這邊的“密度”是指這網絡中點和點間聯繫數量相對於總數兩兩都連接的比值,也就是說假如有k個點,則當每個點都相互連接,則有k(k-1)的聯繫,在雙向的聯繫網絡中則算是有k(k-1)/2個聯繫,這就是網絡內點與點最大值,而實際上所有點相互間聯繫/點與點都聯繫的值就是這邊的density。
  • Components:
    這邊是指一個網絡中有幾個group,這邊的group裡的點都沒有跟彼此相互聯繫。
  • Diameter:
    這邊看的是這網路中距離最遠的兩個點,其之間經過的節點數
  • Clustering Coefficient
    這項指標就比較難理解一點,這邊看的是網絡裡面的點以特定節點為中心聚合的比例,以三個點互相連結為最小單位,簡單說就是當兩個點都有共同連接點時,這兩個點也相互連接的傾象。

screenshot.png

首先必須先安裝statnet和UserNetR兩個packages,其中UserNetR是這本User’s guide to network analysis in R 裡所有示範資料的packages,而statnet則是做network analysis的package。

library(UserNetR)
library(statnet)
data(Moreno)

先簡單把Moreno這組資料畫圖

gender &lt;- Moreno %v% &quot;gender&quot;
plot(Moreno, vertex.col = gender + 2, vertex.cex = 1.2)

用 network.size函數來看這個network object裡面的節點數

#5-number summary
 #1st: Size
network.size(Moreno)

用gden函數他會算出這network裡頭的density,一個介於0到1的數值,即目前這個網絡中所有連結數量與最大連結數量的比值

 #2nd: Density, range from 0 to 1, unidirect or direct
gden(Moreno)

直接算是這個網絡中幾個無互相聯繫的subgroup

 #3rd: Components, subgroups
components(Moreno)

算出網絡中兩者連結需要經過的最大節點數,並且可以算出整個網絡的cluster coefficient

 #4th: Diameter, how compact the network are
lgc &lt;- component.largest(Moreno,result=&quot;graph&quot;)
gd  &lt;- geodist(lgc)
max(gd$gdist)
 #5th: Cluster Coefficient, range from 0 to 1 ,measure the transitivity
gtrans(Moreno,mode=&quot;graph&quot;)

參考資料:Luke, Douglas A. A User’s Guide to Network Analysis in R

Python 3.4環境準備:pyvenv

在開始好好學習python前,當然要先瞭解怎麼管理他的套件和有無乾淨的“玩法”!

有幾個工具可以用來幫忙setup好的工作環境:
– pyenv :          可用來安裝不同python版本的隔離環境
– pyvenv:      在python 3.4之後就有的,很方便的就可以setup一個環境,也是這邊所使用的
– virtualenv:在pypi中用來建隔離環境的工具

 

雖然mac本身安裝時就內建有python 2.7,但直接大剌剌地在global environment玩耍,後果可能會墜入"dependency hell",所以不管要用初心者要用python來完成什麼有趣的專題,起頭一定要有好的習慣。

而在python裡頭,很幸運的,在python3.3版後,便內建很多很棒的套件管理pip,setuptools和venv(pyvenv)可以幫忙從一開始就維持一個好的空間,不用另外再把這些python module tool另外安裝,ㄉ可以統一用brew來處理,這邊我嘗試使用的方式來維持整個環境的整潔,整個global mac用brew來處理其內的套件,在python3所管轄的範圍內則是用python virtual environment來處理:

screenshot.png

首先,先在mac安裝homebrew套件,接者使用brew來 install python3,此種安裝方式會把pip3和setuptools及pyvene都裝好(超級佛心的),可以稍微注意一下homebrew都會把其城市安裝在哪個目錄下,基本上應該會在/usr/local/Cellar目錄下。

接者在使用pyvenv來創建虛擬環境:
pyvenv 指定要放置虛擬環境的目錄位置

$pyvenv  yout/virtual/envir/directory_name

這個創建出來的檔案夾會有這樣的結構:(基本上就是一個獨立的環境,裡面有python binary檔、所使用的module等等)
screenshot.png

那要如何進入這個虛擬環境呢?

source <venve>/bin/activate 

這樣就行了!此時promp會變成:
screenshot.png

##這時候要注意不同的shell啟動方式會不一樣,要稍微注意一下!##

screenshot.png

screenshot.png

 

dplyr使用問題:data_frames can only contain 1d atomic vectors and lists

分享在使用dplyr時遇到的狀況:screenshot.png

這是一個很平常的將檔案讀進來的步驟,但比較特別的是我使用了readr package的read_csv來讀取。

而讀進來後的nor.count長這樣

screenshot.png

接下來我做這些事情,將column重新命名,且把gene list的名字提取出來,為了等會轉換成單純的matrix,提取後在使用dplyr的select,可以直接去掉gene的這行。當nor.count為完全的numeric matrix時,我將數值做log2 transformatiopn,接下來便要開始來處理資料

screenshot.png

接者想要把資料裡的row不等於0濾掉,所以使用filter上面 (錯誤訊息發生screenshot.png

從錯誤訊息,看起來可能是nor.count的格式不符合,再來看一下現在nor.count變成什麼樣子了

screenshot.png

這時候發現一個奇怪的地方,在nor.count的structure中,看到gene下面還有一個data.frame Gene!!

原來是在tbl_df的object中,用傳統的方式subset他還是會保持tbl_df,不會因此變成atomic vector的形式

screenshot.png

像上面這種是很常見的subsetting的方法,subset完後,就照下面的方式再塞回data.frame中,但在tbl_df的object用此操作則會發生問題,這樣就會變成在data.frame裡面塞入一個data.frame,而這種形式就會造成之後的function無法處理

screenshot.png

所以面對tbl_df就必須乖乖的使用screenshot.png

把gene.list給放回去!

 

文獻閱讀:RNAseq在臨床應用的可能和挑戰(下)

原始論文:Byron, S. A., Van Keuren-Jensen, K. R., Engelthaler, D. M., Carpten, J. D., & Craig, D. W. (2016). Translating RNA sequencing into clinical diagnostics: opportunities and challenges. Nature Reviews Genetics, 17(5), 257–271. doi:10.1038/nrg.2016.10

 

論文的最後,談到這些新科技如何整合進臨床之中

第一個重點就是“新的檢測(這邊指得當然就是RNAseq)”,要能變成現行臨床工作的一部份,需要符合許多特性,可以由三個層面一步步來看,分別是Analytical validity (檢測本身的效力)、Clinical validity(臨床上的檢測能力)、Clinical utility(對於最終的臨床決策的貢獻度)。

第一步:Analytical validity

相對於實驗室裡的檢測,處理的檢體都是在control condition下,所以要能讓每一次的實驗都能reproducible的技術門檻就比較低,對比於臨床簡體的複雜性,還要能在此狀況下維持檢驗的穩定性、感測極限等等,就是第一關的重要門檻,要能在有良好的檢測標準來提供檢測sensitivity、specificity的驗證,再來是每次重複的結果應該要相似的(reproducibility),不會因為微小的變數造成數值的大波動。(robustness),如前exosome的研究便是因為其容易波動而造成無法定量。

以RNAseq來看,撇開前製定序所花的不算,光後端的分析方式、參數的差異,就會造成結果有所不同,在2013年,Genetic European Variant in Disease consortium為了解決這問題設計了一組實驗來處理這些technical reproducibility和feasible的問題,收集了465個人的lymphoblastoid cell,在七個定序中心分別執行,最後總結了一些在後端上可以用來驗證定序結果的,像是GC content、fragment size、transcript length、percentage of reads mapped to annotated exon (這部分可以參考所謂的annotation QC)。除此之外,由FDA領頭的SEQC studyABRF study也有大量研究為了解決和提出相關標準。

Su, Z. et al. A comprehensive assessment of RNA-seq accuracy, reproducibility and information content by the Sequencing Quality Control Consortium. Nat. Biotechnol. 32, 903–914 (2014).

Li, S. et al. Multi-platform assessment of transcriptome profiling using RNA-seq in the ABRF next-generation sequencing study. Nat. Biotechnol. 32, 915–925 (2014).

 

第二步: Clinical validity

當檢測技術能在臨床檢體上保持者可重複性時,接下來檢測就要達成可以實際區分biological difference的區別力。

第三步: Clinical utility

當前兩部都達成後,再來要看這個檢驗的結果能如何改變治療,如同companion diagnostic等等,通常此時就需要setup 臨床的試驗,來看實際在臨床場域中的效果!

screenshot.png

 

除此之外,FDA也有針對RNAseq為主的檢測項目做出一些法規的架構,可以由此深入看美國FDA對於NGS技術延生的檢測所抱持的想法

US Department of Health & Human Services. Center for Devices and Radiological Health. FDA notification and medical device reporting for laboratory developed tests (LDTs) — draft guidance. [online]  (2014).

US Department of Health & Human Services. Center for Devices and Radiological Health. Framework for regulatory oversight of laboratory developed tests (LDTs) — draft guidance. [online]  (2014)

US Department of Health & Human Services. Optimizing FDA’ s regulatory oversight of next generation sequencing diagnostic tests — preliminary discussion paper. [online]  (2014).

Evans, B. J., Burke, W. & Jarvik, G. P. The FDA and genomic tests—getting regulation right. N. Engl. J. Med. 372, 2258–2264 (2015).