2015-12-06 82 views
0

我如何知道函數(x),以及函數的座標,從edgeR軟件包生成散佈圖(BCV)?從邊緣R的BCV圖中找出座標R

我需要知道「趨勢」色散的「頂點座標」(每個頂點的最大值和最小值),以及「趨勢線(藍色)」與公共色散線(紅色線)相交的位置的座標(

enter image description here

爲了方便起見,我做一個例子:

data <- "chr start end depth depth1 depth2 depth3 depth4 depth5 depth6 
chr1 3273 3273 7 200 35 1 200 850 0 
chr1 3274 3274 3 50 25 5 300 1500 2 
chr1 3275 3275 8 600 15 8 100 300 5 
chr1 3276 3276 4 30 2 10 59 20 0 
chr1 3277 3277 25 20 7 4 600 45 0" 
data <- read.table(text=data, header=T) 

datamatrix <- data [, c(4:10)] 
library (edgeR) 

#grouping factor 
group <- c(1, 2, 2, 2, 2, 2, 2) #groups of each sample) 

#create a DGEList 
y <- DGEList(counts=datamatrix,group=group) 

根據軋邊機包,我可以利用估計我的數據集的分散體:

y <- estimateDisp(y) 

y$common.dispersion 

常見色散的平方根給出生物變異係數。而且,分散估計可以在BCV情節來看:

plotBCV(y) 

回答

0

我認爲這解決了招

#finding dispersion coordinates 

    afun <- approxfun(y$AveLogCPM, y$trended.dispersion) 

#ranging between -2 and 12 the "par" number to check the intersection 
    optim(par=-1, function(x) { (afun(x) - y$common.dispersion)^2 }) #-0.9056 
    optim(par=0, function(x) { (afun(x) - y$common.dispersion)^2 }) #-0.4955 
    optim(par=2, function(x) { (afun(x) - y$common.dispersion)^2 }) #1.8237 

#with the plotBCV function code, I have the access to tagwise, trended, common dispersions, and the abundance values, for each CpG as: y$tagwise.dispersion, y$trended.dispersion, y$common.dispersion, and y$AveLogCPM, respectively. 

min (y$trended.dispersion) #[1] 1.09009 
max (y$trended.dispersion) #[1] 1.50691 

#to the vertex, you can subset a table according to the AveLogCPM range that you want to know the max and min like this: 

afuntable<- data.table(y$AveLogCPM, y$trended.dispersion) 
afuntablemod <-subset (afuntable, V1 >= 0) 
afuntablemod <-subset (afuntablemod, V1 <= 2) 
min (afuntablemod$V2) #[1] 1.318002 

#and again 
afuntable<- data.table(y$AveLogCPM, y$trended.dispersion) 
afuntablemod <-subset (afuntable, V1 >= 2) 
afuntablemod <-subset (afuntablemod, V1 <= 12) 
max (afuntablemod$V2) #[1] 1.50691