R中有很多的模組與套裝軟體都是PHP所沒有的, 若能夠將這兩個平台結合起來將會將很多工作簡化起來, 以下是結合PHP與R的一個範例, 但通常我們使用到網頁服務(web service)都要考量到TIMEOUT的問題. 若R需要大量的computational time, 則表示PHP勢必會造成TIMEOUT; 所以解決方法可能就要再另外寫一個LOADING BAR, 或者用FLASH, SILVERLIGHT等RIA的方式來操控了.

以下範例應該是以LINUX裝R的機器為主, 若要在windows Server的R上去執行exec()可能需要設定其他的環境變數, 由於尚未實際操作, 所以先略過這個部分了.

http://www.r-bloggers.com/integrating-php-and-r/

台灣R的官方網站載點: http://cran.csie.ntu.edu.tw/

Preprocessing includes
– Image analysis and data import
– Background adjustment
– Normalization
– Summarization (for Affymetrix)
– Quality assessment

Practise Data:

使用limmaGUI匯入two-channel microarray data

參考: http://bioinf.wehi.edu.au/limmaGUI/

A. 下載 Swirl Zebrafish資料

http://bioinf.wehi.edu.au/limmaGUI/Swirl/swirl.zip

B. 安裝 Bioconductor

> source("http://bioconductor.org/biocLite.R")

> biocLite()

C. 安裝以下packages: convert, limma, arrayQuality, marray, mclust, hexbin, limmaGUI, sma, tkrplot, R2HTML, statmod

  1. 程式套件 > 選擇存放處 > CRAN, CRAN(extra), BioC software
  2. 程式套件 > 安裝程式套件 (或用install.packages("convert")

利用affylmGUI匯入two-channel microarray data:

準備動作:

A 安裝estrogen package以取得estrogene dataset:

  1. 程式套件 > 選擇存放處 > BioC experiment
  2. 程式套件 > 安裝程式套件

B. 安裝affylmGUI packages:

  1. 程式套件 > 選擇存放處 > BioC Software
  2. 程式套件 > 安裝程式套件

Normalization

To identify and remove systematic sources of variation in the measured fluorescence intensities, such as

- different labelling efficiencies of the dyes;

- different scanning parameters;

- sector/print-tip, spatial, or plate effects, etc.

不要做過多的normalization, 因為可能把原本是signal的data弄不見.

Self-self hybridization: 去探討dye的系統性誤差的典型例子; 使用相同的sample, 不同的染劑, 去做染色就可以知道不同染劑的差異性, 理論上呈現的pattern應該是如下左圖, 但實際上會是長得像右圖的:

所以從以上的例子來看, 矯正(Normalization)的動作變得很重要.

可以用: Within-Array Normalization.

實作:

使用limmaGUI來作Normalization:

Normalization > Select Within-Array Normalization

Posted in R.

講師: 陳倩瑜

Feature Selection 課程組織
屬性的選擇

Clustering:分群
(Unsupervised Learning)

Classification: 分類
(Supervised Learning)

這堂課的大鋼會比較涵蓋在演算法上, 一些函數的背景理論與應用.

首先假設我們拿到的資料已經經過整理,
來自不同病人的資料整理成一個二維的陣列.

縱軸可能是 時間的差異, Condition,
----------------------------------------
| sample |
----------------------------------------
| gene1  |
----------------------------------------
| gene2  |
----------------------------------------

若作clustering, 屬性相同的會被歸類在一起.Data Clustering concerns how to group similar objects together while spearating dissimilar objects.如何判斷行不行? 當你拿到一群資料的時候, 就可以從一大群資料去判斷有甚麼關連性, 以及如何分群.

這樣的方法在很多地方都有用到, 如: machine learning, Data mining, Pattern recognition, Image Analysis, Bioinformatics.

Hierarchical

http://en.wikipedia.org/wiki/Cluster_analysis

http://nlp.stanford.edu/IR-book/html/htmledition/hierarchical-agglomerative-clustering-1.html

巢狀式,
階層式 (應用在基因的概念上)

Agglomerative
Divisive
HAC (hierarchical agglomerative clustering) 先把像的東西放在一起, 決定第一刀切在哪裡.

方法如下:

1. 先決定兩個人的相似程度, 比如: A有25000個FEATURE, B有25000個FEATURE

Proximity matrix containing the distance between each pair of objects. Treat each object as  cluster.

2. 把兩個最像的東西放在一起, UPDATE 剛剛的MATRIX

> class(iris)
[1] "data.frame"
> dim(iris)
[1] 150   5
> class(iris[1])
[1] "data.frame"
> class(iris[2])
[1] "data.frame"
> class(iris[4])
[1] "data.frame"
> class(iris[,1])
[1] "numeric"
> class(iris[,2])
[1] "numeric"
> class(iris[,6])
錯誤在`[.data.frame`(iris, , 6) : undefined columns selected
> class(iris[,5])
[1] "factor"
> table(iris[,5])

setosa versicolor  virginica
50         50         50

用 Iris 來繪圖 > plot(iris[1:4], col=iris[,5]) ,

顏色表示以第五個欄位來分類(因為有三個category,所以有三種顏色)

分群常用指令:

> h <- agnes(iris[3:4])
> h

Call:    agnes(x = iris[3:4])
Agglomerative coefficient:  0.9791964
Order of objects:
[1]   1   2   5   9  29  34  48  50   3  37  39  43   7  18  46  41  42   4   8  11
[21]  28  35  40  49  20  12  26  30  31  47  21  10  33  13  38  17   6  27  24  19
[41]  16  22  32  44  14  15  36  23  25  45  51  64  77  59  92  57  87  52  67  69
[61]  79  85  55  86 107  74  56  88  66  76  91  75  98  95  97  96  62  54  72  90
[81]  89 100  93  60  63  68  70  83  81  58  94  61  80  82  65  99  53  73  84 120
[101] 134  71 127 139 124 128  78 102 143 147 150 111 148 114 122 112 101 110 136 103
[121] 105 121 144 137 141 145 113 140 125 129 133 115 142 116 146 149 104 117 138 109
[141] 130 135 106 123 118 119 108 132 126 131
Height (summary):
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
0.0000  0.0000  0.1000  0.1864  0.2038  3.7370

Available components:
[1] "order"  "height" "ac"     "merge"  "diss"   "call"   "method" "data"
> h$order
[1]   1   2   5   9  29  34  48  50   3  37  39  43   7  18  46  41  42   4   8  11
[21]  28  35  40  49  20  12  26  30  31  47  21  10  33  13  38  17   6  27  24  19
[41]  16  22  32  44  14  15  36  23  25  45  51  64  77  59  92  57  87  52  67  69
[61]  79  85  55  86 107  74  56  88  66  76  91  75  98  95  97  96  62  54  72  90
[81]  89 100  93  60  63  68  70  83  81  58  94  61  80  82  65  99  53  73  84 120
[101] 134  71 127 139 124 128  78 102 143 147 150 111 148 114 122 112 101 110 136 103
[121] 105 121 144 137 141 145 113 140 125 129 133 115 142 116 146 149 104 117 138 109
[141] 130 135 106 123 118 119 108 132 126 131
> h$height
[1] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.1000000
[9] 0.0000000 0.0000000 0.0000000 0.1193300 0.0000000 0.0000000 0.1000000 0.0000000
[17] 0.1884517 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.1000000
[25] 0.1214732 0.0000000 0.0000000 0.0000000 0.0000000 0.1000000 0.1599476 0.0000000
[33] 0.1000000 0.0000000 0.2577194 0.2986122 0.1000000 0.1207107 0.1471405 0.1868034
[41] 0.0000000 0.0000000 0.2217252 0.3709459 0.1414214 0.0000000 0.1804738 0.4836443
[49] 0.2000000 3.7365080 0.0000000 0.1000000 0.1510749 0.1000000 0.1869891 0.1000000
[57] 0.2335410 0.0000000 0.0000000 0.0000000 0.0000000 0.1000000 0.1554190 0.1000000
[65] 0.3195579 0.3582523 0.1000000 0.1207107 0.0000000 0.1603553 0.2038276 0.0000000
[73] 0.1000000 0.0000000 0.1207107 0.2598380 0.5662487 0.0000000 0.0000000 0.1000000
[81] 0.0000000 0.1165685 0.1825141 0.2663756 0.1000000 0.2352187 0.1000000 0.1207107
[89] 0.9780526 0.0000000 0.2666667 0.0000000 0.2000000 0.3594423 0.4986370 1.5284638
[97] 0.0000000 0.1745356 0.1207107 0.1000000 0.3927946 0.0000000 0.0000000 0.1000000
[105] 0.0000000 0.1907326 0.2947898 0.0000000 0.1000000 0.1138071 0.1589347 0.1000000
[113] 0.1898273 0.1000000 0.2521467 0.8001408 0.1000000 0.2118034 0.4796864 0.2080880
[121] 0.1414214 0.1707107 0.2979830 0.0000000 0.1414214 0.3207243 0.1000000 0.1941714
[129] 0.1207107 0.1000000 0.5593867 0.1000000 0.2135427 0.1000000 0.1500000 0.6629414
[137] 0.1000000 0.0000000 0.3006588 0.2000000 0.3909355 1.1817205 0.1414214 0.1707107
[145] 0.3149057 0.6020066 0.2236068 0.3217620 0.1414214
> iris[h$order[1:50],5] //找出1到50已經分群的資料
[1] setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa
[12] setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa
[23] setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa
[34] setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa
[45] setosa setosa setosa setosa setosa setosa
Levels: setosa versicolor virginica

> k <- agnes(iris[1:4])
> k
Call:    agnes(x = iris[1:4])
Agglomerative coefficient:  0.9300174
Order of objects:
[1]   1  18  41  28  29   8  40  50   5  38  36  24  27  44  21  32  37   6  19  11  49  20  22  47  17  45   2  46  13  10  35  26   3   4
[35]  48  30  31   7  12  25   9  39  43  14  23  15  16  33  34  42  51  53  87  77  78  55  59  66  76  52  57  86  64  92  79  74  72  75
[69]  98  69  88 120  71 128 139 150  73  84 134 124 127 147 102 143 114 122 115  54  90  70  81  82  60  65  80  56  91  67  85  62  89  96
[103]  97  95 100  68  83  93  63 107  58  94  99  61 101 121 144 141 145 125 116 137 149 104 117 138 112 105 129 133 111 148 113 140 142 146
[137] 109 135 103 126 130 108 131 136 106 123 119 110 118 132
Height (summary):
Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
0.0000  0.2189  0.3317  0.4377  0.5081  4.0630

Available components:
[1] "order"  "height" "ac"     "merge"  "diss"   "call"   "method" "data"
> table(iris[k$order[101:150],5])

setosa versicolor  virginica
0         13         37
> table(iris[k$order[1:50],5])

setosa versicolor  virginica
50          0          0

Measures of Dissimilarity/Similarity Between Objects

  • Euclidean Distance
  • Pearson Correlation coefficient

Normalization

  • Decimal Scaling
  • Min-Max normalization
    for normalized interval [0,1]
  • Standard deviation normalization
  • 在R作normalization可使用scale, 範例如下:

> j <- data.frame(scale(iris[1:4])) //先轉換為data.frame (因為SCALE之後就是matrix)
> class(j) //查看j 的資料結構
[1] "data.frame"
> plot(j, col=iris[,5]) //畫圖

經過Normalization的Iris圖

TIPS:

1. 太高維度的OBJECT就不再相像

2. 分群其實沒有標準答案, 怎麼分是由自己決定; 應該是由自己分群後再向別人解釋自己的方法以及結果.

3. 東西如果像, 怎麼分都會在一起(如使用Single-Link與 Complete Link 演算法).

4. 點跟點的距離, 群跟群的距離決定了不同的結果. 這就是為甚麼不同的Algorithm會造成不同的樹狀圖.

> plot(agnes(iris[1:4], method="single"))
//使用single-link演算法畫出的樹狀圖, 會發現分得不錯, 有兩大族群

//UPGMA (Average Link)

Practise:

1. x <- read.delim(filechoose(), row.names = 1, header= TRUE)

2. mean(as.numeric(x[1:27]))

3. a <- function(y) {

(mean(as.numeric(x[y,1:27])) – mean(as.numeric(x[y,28:38])) ) / (sd(as.numeric(x[y,1:27])) – sd(as.numeric(x[y,28:38])) )

}

4. for(i in 1:7129) {

x[i] = a(i);

}

Partitional

K-Means

Posted in R.

Class Comparison: 比較不同Class的類別, 屬性

Class Discovery: 從一群Data中去歸類, 相似性

Class Prediction: 有新的data, 能不能用來歸類於現有的Class


Class Comparison

找出顯著差異的.

不考慮基因與基因的Interactions.

用FoldChange(FC) = Expression of Experriment al Sample / Expresssion of Reference Sample

可用Scatter plots 及 MA-plots Visualized, 但Fold-Change沒辦法找出顯著差異, 不具有統計意義,但還是很多人喜歡使用, why? 因為統計上需要大量的SAMPLE, 而一般實驗沒有那麼多的sample, 只能用Fold Change 來表示.

P-value: 機率的概念. 在正常情況下, 我看到這樣的CASE的機率到底有多少? 例如: 在大家都買樂透的情況下. p-value越小越好, 表示越顯著. alpha-value是臨界值, 若p-value小於alpha-value, 則有顯著. 所以一般上習慣把alpha-value設得比較小.

分析一般的 Microarray data方法:

Parametric  Hypothesis Testing (有常態分布)

  • Paired Data :z-test, t-test
  • Unpaired Data :two-sample, t-test
  • Complex Data (more than 2 groups) :One-way Analysis of Variance(ANOVA)

Non-Parametric Hypothesis Testing (不合常態分布)

  • Paired data: Sign test, Wilcoxon signed-rank test
  • Unpaired Data :Wilcoxon rank-sum test, (Mann-Whitney U test).
  • Complex Data :Kruskal-Wallis test

如何知道有無常態分布? 用 QQplot, Shapiro….

Parametric Statistics

Paired Data: 如化療前, 化療後的比較, 吃藥前, 吃藥後.

t.test(x,y,paired =TRUE, alternative = c("two.sided", "less", "greater"))

一定要先對好兩組數據, 然後參數paired = TRUE要寫.

Unpaired Data: sample來自不同的地方

先用var.test 測試p-value是否有顯著, 若越大就越不顯著.

if(sigmaa = sigmab),

var.test(x,y, var.equal = TRUE)

if(sigmaa != sigmab),

var.test(x,y, var.equal=FALSE)

Complex Data: (Using ANOVA)

Example: All dataset

y = drop(exprs(eset[1,]))

out = lm(y~factor(c1)) //lm 變方分析

anova(out)

Non-Parametric Statistics

基本上microarray 是很雜亂的, 不適合用以上的parametric 方式, 故使用這種統計方式.

Paired

wilcox.test(x, y, paired= TRUE, alternative = c("two.sided", "less", "greater"))

Unpaired

wilcox.test(x, y, alternative = c("two.sided", "less", "greater"))

Anova

參考投影片


Permutation Tests

What’s this? 排列. 隨機排列.

比如說, 有10個基因, 每個作2次 就有20次. 把10個結果作隨機排列, 分前五個, 後五個, 計算每個出現的機率.

會形成一個常態分布. Steps可參考教學投影片.

資料數太少的時候不適合用Permutation test, 應該使用"Bootstrap".

用R作permutation test:

"multtest" package

mt.sample.teststat: to compute permuted statistcs

mt.sample.rawp: to compute the p-values

Note: "test" includes:

t, t.equalvar, pairt, wilcoxon, f

Bootstrap

boot.out = boot(englife, function(x, id) {median(x[id])}, 1000}

how many bootstraps ? rule of thumb: 100 times. 可以再試試看1000次, 10000次.

Posted in R.