分類: R

使用rprojroot改善麻煩的檔案路徑問題，以子資料夾結構取而代之

發表留言

screenshot.png

在Jenny Bryan登高一呼，希望有人解決這種資料夾路徑總是要重設的問題時，這個rprojroot 包便誕生了，他主要解決了路近依賴的問題，而使用sub-directory結構取而代之，本身是個簡單的package，但算是為培養好的reproducible coding踏出一步。

rprojroot要發揮作用的話，要建立在兩個使用者的好習慣：
第一個：每個專案有固定的子資料夾結構，這我個人覺得很重要，專業素養便是從檔案管理做起，可以看下面的文章，對於如何管理好專案，裡面的建議頗受用

Managing a statistical analysis project – guidelines and best practices
R Workflow: Slides from a Talk at Melbourne R Users
Dynamic documents with R and LATEX as an important part of reproducible research
How to efficiently manage a statistical analysis project?

第二個：願意系統性的寫好代碼(可能30秒能搞定的代碼，必須多花兩倍時間，為了代碼的嚴謹)

     假如使用者有上面這兩個習慣，這個包便可以幫助你的R代碼進步一點點，此包主要使用一個概念，便是為專案root資料夾設立“criteria”，藉由檢查當前的資料架內結構，是否為特定的形式如r package、r project，或是包含特定檔名的檔案。藉由這個criteria，來確定root路徑，而剩下的檔案便是用相對路徑來表示。

這邊有示範的代碼：

library(rprojroot
#Make Criteria
root <- has_file("README.md")

#Fine the file
dir <- find_root(root)
data.path <-find_root_file("Data", criterion = root)

#True or False criteria
file.exists(find_root_file("Data", criterion = has_file("README.md")))

這邊則是我的放代碼的資料夾結構：
screenshot.png
所以使用這個包來設定檔案路徑，未來只要整個資料夾給要合作的對象，基本上，重新執行就不會發生問題。

在dplyr的pipe中，在mutate裡使用cut來將連續變數依據區間分塊，形成新的變數

發表留言

     本來在dplyr的函數中是不支持cut函數的（cut為R base函數），但將連續變數依照特定區間分塊貼標籤這使用場景太常見了，要把班上成績，依據80,85,90分別給予C,B,A等操作，且在stackover可以看到非常多類似的問題如下：
Create column with grouped values based on another colum

R dplyr – categorize numeric variable with mutate

Is cut() style binning available in dplyr?

cut function not working in dplyr, but works outside

applying cut within dplyr

screenshot.png

cars %>% mutate(type = cut(speed,c(0,15,30),c("slow","fast")))

screenshot.png

tidy data- 怎麼樣的資料才是“整齊乾淨的”

1 則迴響

在一開始使用ggplot2時，一定會被其要的data input格式所困擾，其實那是養成好習慣和觀念的開始，在Hadley Wickham 最近剛完成的R for Data Science  裡頭很重要的一件事就是使用tidy的data來“往返”整個資料分析的過程，讓各個package之間的I/O能完美的銜接，且從根源解決常因資料格式不同而要從新處理的問題。

關於怎樣的資料才算是“tidy”的呢？這篇論文中Hadley Wickham已經非常仔細地探討了，有興趣可以進一步點開。

三個重點決定了一筆資料的乾淨與否

  1. 每行就代表一個變數（量測的性質）
  2. 每一列為一個觀察點（樣本、或是病人、檢體..）
  3. 每一個table就是在講一類事情

 

 

髒的資料，無法直接filter，需要進一步整理
screenshot.png
乾淨的資料，可直接filter，甚至在excel中就可以處理
screenshot.png
相關跟處理tidy data的包，有下面這些：
  • magrittr: 主要將pipe的概念帶入R中，使用%>%符號
  • dplyr: 一整套處理data.frame的包
  • tidyr: 將資料整理成tidy form的包，其中gather, spread是最關鍵的函數
  • broom: 將R內建常見的統計輸出整理成tidy form

使用ggplot2: stat_function 簡化函數畫圖

發表留言

在R裡面有很多可以用來畫函數的方法，像是curve就是一個簡單懶人函式，其實ggplot2裡頭有一個超好用的函式，不需要像一般ggplot2必須輸入一個data.frame的data.set，只要把想要畫圖的function定義好即可。這邊有documentation的連結。

ggplot()+stat_function( aes(x=0:2), fun = 自訂函式,  args = list())

aes內的x設定要帶入函數的範圍，也是圖x軸的邊界，fun後面則放入想要畫圖的函式，args裡面可以把原本function中其他的參數指定好，這樣就可以開始畫圖了，下面用簡單的例子展現一下他的方便！

 

screenshot.png

screenshot.png

寫出漂亮的R程式碼

發表留言

就像文章會有固定的寫法，程式語言也是，練習寫得簡潔易讀，是需要多加練習的，可以參考Hadley Wicham在Advance R裡寫關於寫作style的建議，另一個則是google 的 R style rule。整體來說是很簡單的。

分成以下幾類來改善R code 的撰寫：

  • 檔案名稱
    • 將檔案名稱以底線連接，且最好檔名要讓人一眼知道裡面是什麼
      • 好的命名
        • predict_ad_revenue.R
      • 壞的命名
        • foo.R
  • 變數(identifiers)的名字
    • 可以分成兩種，函數的名稱或是其他一般的變數名稱
      • 一般變數：盡量不要用底線連接，而是用逗點來連接
        • 好的命名
          • variable.name
        • 壞的命名
          • variable_name
      • 函數名稱:加入大寫，不要用逗點或是底線
        • 好的命名
          • CalculateAvgClicks
        • 壞的命名
          • calculate_avg_clicks
  • 每行的長度
    • 不要超過80個字碼
  • 縮進（indentation）
    • 不要混用tab和space，最好的方式是使用 2 spaces
  • 空白格（spacing）
    • 在所有運算符號前面加上一個空白格如（=+<-..etc.）

RNAseq: Annotation Based Quality Control

2 則迴響

當我們將定序出來的reads使用參考基因序列對應回位置後，我們便可以進一步計算每個基因上、exon、transcripts上對應的reads數量，來定量基因的表現，除此之外我們還可以用這些資訊來做進一步的實驗或是定序品質檢驗。

 

理論上來說，計算每個基因上所對應的reads數，可以很簡單的定量出基因表現，但在真核生物中，同樣的基因會有不同的isoforms，這樣便會出現一些問題，除此之外，在製備定序使用的library其protocol也會有影響，如rRNA有無去除乾淨等等，這些可以在alignment後對bam檔做annotation-based的quality control

 

雖然我們已經對fastq檔案做過初步的fastqc，但許多重要的資訊單純的fastq檔案無法提供，因為在fastq檔案裡我們只能對每個base-pair其定序出來的品質來檢查，如有沒有被rRNA汙染、或是製備protocol的檔案有沒有造成一些bias等等。

 

以下列出一些重要的品質指標(annotation-based quality metrics)

 

定序深度是否過深Saturation of sequencing depth

定序的深度其實會決定整個RNAseq實驗後續的基因表現分析、RNA剪接或是尋找特殊transcripts，所以到底實驗設計想要回答的問題需要定序多少深度會適合，而這決定於目標轉錄體的複雜性。

定序片段的在特定功能基因區分布是否有問題Read distribution between different genomic features

 可以看定序片段（reads）在不同的genomic features上的分佈，像是外顯子區域（exonic region）、內含子區域（intronic regions）、intergenic regions（基因間隔區），看分佈是否有特殊的地方，可能會發現新穎的基因或是isoforms，當然也可能會發現有無被污染，尤其是ribosomal RNA的比例，這可以知道是否在製備library時有清除乾淨。

transcripts上是否有均一地被定序到Coverage uniformity along transcripts

  不同實驗的步驟可能會造成transcripts上某一部份location被定序到的比例不同，比如說要是在製備library的時候，有poly-A capture的步驟，則其3’ ends就會有比較高被定序到的現象．

 

可以用來做annotation-based quality control的工具

 

  • RNA-SeQC
    • 為Java撰寫的程式，採用GTF格式的annotation來分析，另外，其需要bam所使用的reference fasta且要有.fai index和sequence dictionary file(,dict)，其會輸出html格式的檔案，裡面包含mean coverage、針對不同location transcripts coverage的計算（3’或是5’ bias)，會將基因的表現和GC content分成三個程度：low、medium、high。這篇文章有進一步介紹其功能
  • Qualimap
    • 為Java撰寫的程式，裡頭會使用到R/bioconductor的package，除了使用GTF格式，也採用BED格式作為annotation，其對於saturation和biotype 分佈的分析還不錯。其對於測序深度的分析會包含不同genomic features的深度分析、且會提供每增加一百萬的定序深度會提高發現多少新的features的機會，在biotype分佈分析，會將定序片段（reads）分佈在蛋白質轉譯區域(protein coding region)、偽基因（pseudogens）、rRNA、miRNA等地方顯示出來。
  • RseQC
    • 為python撰寫的程式，其特色是在計算不同genomic features上定序片段分佈時，會把上游和下游的transcripts都一起展現出來，另外一點特色則是其對splice junction的分析會標示出是否為novel
  • Picard’s CollectRNASeqMetrics
    • Picard工具裡其中一個Java 程式，其特色為會分析出strand specific的程度

R Shiny 快速使用R搭建app應用

發表留言

Rstudio 開發的shiny基本架構

Shiny application是一款以R語言為框架，可以快速且方便搭健application，當幫忙使用R寫完清理或是視覺化資料後，最重要的是如何讓不懂R語言的人也可以享受，此時，用RStudio開發的Shiny框架便能快速有這效果，且其跟Rstudio整合得非常好，貼心宜人，即使非寫網頁應用程式出生的生物學家，只要略懂R語言，便能完成！

screenshot.png

shiny語法有非常詳細的解說，且搭配實例，只要到官方網頁就可以瀏覽，大概一兩個小時就可以了解他的架構，最棒的事其可以跟html / javascripts等等串接在一起，做出更多的效果！

做好的application可以直接使用Rstudio裡頭的快捷鍵發佈到shinyapp的網路平台，其提供後台完整的服務，付費的話可以獲得更佳的效能，但免費版已經可以玩得很開心了！

screenshot.png

 

點擊寫好的code，便能直接deploy到shinyapps的平台，一條龍服務！

screenshot.png

基本上，shiny的概念，就是直接用R code搭建網站的前端（ui.R）和資料處理的後端（server.R），server的code會處理由ui端或是使用者輸入的資料，處理後再丟回給ui!