sadII_lab8_ucsc_data

2282 days ago by agorska

# pakiet "rtracklayer" # Ponizej przykladowa sesja z uzyciem tego pakietu, # natomiast dane mozna takze wczytac z pliku (pare komorek nizej) 
       
library(rtracklayer) session <- browserSession("UCSC") # pakiet umożliwia przeszukiwania UCSC . 
       
Loading required package: RCurl
Loading required package: bitops
tnames<-trackNames(session) length(tnames) head(tnames,20) 
       
[1] 146
    Base Position   Chromosome Band       STS Markers       FISH
Clones       Recomb Rate 
          "ruler"        "cytoBand"          "stsMap"     
"fishClones"      "recombRate" 
     ENCODE Pilot       Map Contigs          Assembly   GRC Map
Contigs               Gap 
  "encodeRegions"          "ctgPos"            "gold"        
"ctgPos2"             "gap" 
    BAC End Pairs  Fosmid End Pairs        GC Percent GRC Patch
Release         Hg18 Diff 
    "bacEndPairs"     "fosEndPairs"         "gc5Base"
"altSeqComposite"  "hg19ContigDiff" 
     Hi Seq Depth     NCBI Incident       Short Match     Restr
Enzymes        Wiki Track 
     "hiSeqDepth"  "ncbiIncidentDb"      "oligoMatch"        
"cutters"       "wikiTrack" 
tnames[grep("miRNA",names(tnames))] # miRNA zassane z UCSC 
       
     sno/miRNA TS miRNA sites 
       "wgRna"  "targetScanS" 
%time query <- ucscTableQuery(session, "wgRna") tableNames(query) 
       
[1] "wgRna"
CPU time: 0.01 s,  Wall time: 6.91 s
tableName(query)<-"wgRna" mirna <- getTable(query) summary(mirna) 
       
      bin           chrom       chromStart           chromEnd       
name          score  
 Min.   : 110   chr14  :118   Min.   :    30143   Min.   :    30281 
14q(0)  :   1   Min.   :0  
 1st Qu.: 805   chr1   :115   1st Qu.: 28906892   1st Qu.: 28907024 
14q(I-1):   1   1st Qu.:0  
 Median :1062   chr15  :115   Median : 62622763   Median : 62622838 
14q(I-2):   1   Median :0  
 Mean   :1175   chr19  :105   Mean   : 77498437   Mean   : 77498529 
14q(I-3):   1   Mean   :0  
 3rd Qu.:1425   chrX   :101   3rd Qu.:110141514   3rd Qu.:110141589 
14q(I-4):   1   3rd Qu.:0  
 Max.   :2485   chr17  : 76   Max.   :249120575   Max.   :249120642 
14q(I-5):   1   Max.   :0  
                (Other):711                                         
(Other) :1335              
 strand    thickStart    thickEnd      type    
 -:610   Min.   :0    Min.   :0   CDBox  :269  
 +:731   1st Qu.:0    1st Qu.:0   HAcaBox:112  
         Median :0    Median :0   miRNA  :939  
         Mean   :0    Mean   :0   scaRna : 21  
         3rd Qu.:0    3rd Qu.:0                
         Max.   :0    Max.   :0                
                                               
# Wczytywanie tabeli miRNA zapisanej wczesniej w innej sesji load("data/miRNA.Rdata") ls() 
       
[1] "ind"     "logRfun" "mirna"   "query"   "session" "st"     
"tnames" 
dim(mirna) head(mirna) summary(mirna) 
       
[1] 1341   10
  bin chrom chromStart chromEnd           name score strand
thickStart thickEnd  type
1 585  chr1      30365    30503 hsa-mir-1302-2     0      +         
0        0 miRNA
2 593  chr1    1102483  1102578   hsa-mir-200b     0      +         
0        0 miRNA
3 593  chr1    1103242  1103332   hsa-mir-200a     0      +         
0        0 miRNA
4 593  chr1    1104384  1104467    hsa-mir-429     0      +         
0        0 miRNA
5 608  chr1    3044538  3044599   hsa-mir-4251     0      +         
0        0 miRNA
6 611  chr1    3477258  3477354   hsa-mir-551a     0      -         
0        0 miRNA
      bin           chrom       chromStart           chromEnd       
name          score  
 Min.   : 110   chr14  :118   Min.   :    30143   Min.   :    30281 
14q(0)  :   1   Min.   :0  
 1st Qu.: 805   chr1   :115   1st Qu.: 28906892   1st Qu.: 28907024 
14q(I-1):   1   1st Qu.:0  
 Median :1062   chr15  :115   Median : 62622763   Median : 62622838 
14q(I-2):   1   Median :0  
 Mean   :1175   chr19  :105   Mean   : 77498437   Mean   : 77498529 
14q(I-3):   1   Mean   :0  
 3rd Qu.:1425   chrX   :101   3rd Qu.:110141514   3rd Qu.:110141589 
14q(I-4):   1   3rd Qu.:0  
 Max.   :2485   chr17  : 76   Max.   :249120575   Max.   :249120642 
14q(I-5):   1   Max.   :0  
                (Other):711                                         
(Other) :1335              
 strand    thickStart    thickEnd      type    
 -:610   Min.   :0    Min.   :0   CDBox  :269  
 +:731   1st Qu.:0    1st Qu.:0   HAcaBox:112  
         Median :0    Median :0   miRNA  :939  
         Mean   :0    Mean   :0   scaRna : 21  
         3rd Qu.:0    3rd Qu.:0                
         Max.   :0    Max.   :0                
                                               
mirna$width<-with(mirna,chromEnd-chromStart) head(mirna) 
       
  bin chrom chromStart chromEnd           name score strand
thickStart thickEnd  type width
1 585  chr1      30365    30503 hsa-mir-1302-2     0      +         
0        0 miRNA   138
2 593  chr1    1102483  1102578   hsa-mir-200b     0      +         
0        0 miRNA    95
3 593  chr1    1103242  1103332   hsa-mir-200a     0      +         
0        0 miRNA    90
4 593  chr1    1104384  1104467    hsa-mir-429     0      +         
0        0 miRNA    83
5 608  chr1    3044538  3044599   hsa-mir-4251     0      +         
0        0 miRNA    61
6 611  chr1    3477258  3477354   hsa-mir-551a     0      -         
0        0 miRNA    96
#indeksowanie wektorem true/false test <- dataframe(subset(mirna,type=="miRNA")) x<-dataframe(test$width[test$strand == "+"])) y<-dataframe(test$width[test$strand == "-"])) 
       
Error: could not find function "dataframe"
Error: unexpected ')' in "x<-dataframe(test$width[test$strand ==
"+"]))"
Error: unexpected ')' in "y<-dataframe(test$width[test$strand ==
"-"]))"
#Kolmogorov-Smirnov test ks.test(x = subset(mirna,type=="miRNA")$width[mirna$strand == "+"] , y=subset(mirna,type=="miRNA")$width[mirna$strand == "-"]) 
       
	Two-sample Kolmogorov-Smirnov test

data:  subset(mirna, type == "miRNA")$width[mirna$strand == "+"] and
subset(mirna, type == "miRNA")$width[mirna$strand == "-"] 
D = 0.0457, p-value = 0.7121
alternative hypothesis: two-sided 

Warning message:
In ks.test(x = subset(mirna, type == "miRNA")$width[mirna$strand == 
:
  p-values will be approximate in the presence of ties
#Kolmogorov-Smirnov test ks.test(x = mirna$width[mirna$chrom == "chr3"] , y=mirna$width[mirna$chrom == "chr5"]) 
       
	Two-sample Kolmogorov-Smirnov test

data:  mirna$width[mirna$chrom == "chr3"] and
mirna$width[mirna$chrom == "chr5"] 
D = 0.1432, p-value = 0.6425
alternative hypothesis: two-sided 

Warning message:
In ks.test(x = mirna$width[mirna$chrom == "chr3"], y =
mirna$width[mirna$chrom ==  :
  cannot compute exact p-values with ties
 
       
library(lattice) histogram(~(chromEnd-chromStart)|type,data=mirna,breaks=20) dev.off() 
       
Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf
null device 
          1 
histogram(~width|type,data=subset(mirna,type=="miRNA"),breaks=50) dev.off() #KS test -> prównuje deystrbuanty -> wiadomo czy dane pochodzą z tego samego rozkładu 
       
Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf
null device 
          1 
png(width=800,height=1000) histogram(~width|chrom,data=subset(mirna,type=="miRNA"),breaks=20) dev.off() 
       
Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf
null device 
          1 
histogram(~width|strand,data=subset(mirna,type=="miRNA"),breaks=50) dev.off() 
       
Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf
null device 
          1 
# Czy rozklady wielkosci na obydwu wstegach sa takie same? Zastosuj odpowiedni test. 
       
## Jaki to moze byc rozklad? Wykorzystujac wiedze z wykladu postaraj sie dopasowac najlepszy rozklad. 
       
# Ilosc miRNA na chromosomach mirna.chr.counts<-table(subset(mirna,type=="miRNA")$chrom) mirna.chr.counts barchart(~sort(mirna.chr.counts)) dev.off() 
       
 chr1 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 
chr2 chr20 chr21 chr22  chr3 
   75    36    41    32    21    69    33    28    45    14    91   
48    30    10    21    44 
 chr4  chr5  chr6  chr7  chr8  chr9  chrX 
   34    39    23    40    39    37    89 

null device 
          1 
# Pytanie: czy ilosc miRNA zalezy tylko od wielkosci chromosomu? query <- ucscTableQuery(session, "wgRna") tableNames(query) # Czy jakis chromosom jest wyjatkowy pod tym wzgledem? 
       
[1] "wgRna"
# Wielkosci chromosomu: sciagamy wspolrzedne prazkow od kariotypu # Sciagnij query "cytoBand" i tabele "cytoBand" 
       
%time library(rtracklayer) session <- browserSession("UCSC") query <- ucscTableQuery(session, "cytoBand") tableNames(query) 
       
[1] "cytoBand"
CPU time: 0.00 s,  Wall time: 7.15 s
tableName(query)<-"cytoBand" cband <- getTable(query) 
       
head(cband) summary(cband) with[cband, tapply(chromEnd, chrom)] 
       
  chrom chromStart chromEnd   name gieStain
1  chr1          0  2300000 p36.33     gneg
2  chr1    2300000  5400000 p36.32   gpos25
3  chr1    5400000  7200000 p36.31     gneg
4  chr1    7200000  9200000 p36.23   gpos25
5  chr1    9200000 12700000 p36.22     gneg
6  chr1   12700000 16200000 p36.21   gpos50
     chrom       chromStart           chromEnd              name    
gieStain  
 chr1   : 63   Min.   :        0   Min.   :  2200000   p11.1  : 20  
gneg   :414  
 chr2   : 62   1st Qu.: 29000000   1st Qu.: 32475000   q11.1  : 17  
gpos50 :121  
 chr3   : 62   Median : 62950000   Median : 66650000   q22.2  : 17  
gpos75 : 89  
 chr6   : 48   Mean   : 72987007   Mean   : 76578280   q21.2  : 16  
gpos25 : 87  
 chr4   : 47   3rd Qu.:107475000   3rd Qu.:111650000   q22.1  : 16  
gpos100: 81  
 chr5   : 45   Max.   :243700000   Max.   :249250621   p11.2  : 15  
acen   : 48  
 (Other):535                                           (Other):761  
(Other): 22  
Error in tapply(chromEnd, chrom) : object 'chromEnd' not found
 
       
chr.sizes<-with(cband, tapply(chromEnd, chrom,max)) library(lattice) barchart(~sort(chr.sizes)) dev.off() 
       
null device 
          1 
# Zadanie: czy ilosc miRNA zalezy tylko od wielkosci chromosomu? # Ktore chromosomy sa wyjatkowe pod tym wzgledem? 
       
mirna.size.data<-data.frame(mirna.cnt=c(mirna.chr.counts,0), chr.size=as.vector(chr.sizes), chr=names(chr.sizes), row.names=names(chr.sizes)) mirna.size.data 
       
      mirna.cnt  chr.size   chr
chr1         75 249250621  chr1
chr10        36 135534747 chr10
chr11        41 135006516 chr11
chr12        32 133851895 chr12
chr13        21 115169878 chr13
chr14        69 107349540 chr14
chr15        33 102531392 chr15
chr16        28  90354753 chr16
chr17        45  81195210 chr17
chr18        14  78077248 chr18
chr19        91  59128983 chr19
chr2         48 243199373  chr2
chr20        30  63025520 chr20
chr21        10  48129895 chr21
chr22        21  51304566 chr22
chr3         44 198022430  chr3
chr4         34 191154276  chr4
chr5         39 180915260  chr5
chr6         23 171115067  chr6
chr7         40 159138663  chr7
chr8         39 146364022  chr8
chr9         37 141213431  chr9
chrX         89 155270560  chrX
chrY          0  59373566  chrY
fit<-lm( chr.size~mirna.cnt, data=mirna.size.data ) summary(fit) 
       
Call:
lm(formula = chr.size ~ mirna.cnt, data = mirna.size.data)

Residuals:
       Min         1Q     Median         3Q        Max 
-116535888  -37935374    3412699   35033864  106226886 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 93780989   22979129   4.081 0.000495 ***
mirna.cnt     899823     511168   1.760 0.092250 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 55440000 on 22 degrees of freedom
Multiple R-squared: 0.1235,	Adjusted R-squared: 0.08362 
F-statistic: 3.099 on 1 and 22 DF,  p-value: 0.09225 
cor(mirna.size.data$chr.size, mirna.size.data$mirna.cnt) 
       
[1] 0.351372
xyplot(mirna.cnt~chr.size, data=mirna.size.data, type=c("p","smooth","r"), auto.key=T) dev.off() 
       
null device 
          1 
png(h=800,w=800) layout(matrix(1:4,2,2)) plot(fit) dev.off() 
       
null device 
          1 
fit19<-lm(mirna.cnt~chr.size,data=subset(mirna.size.data,!chr%in%paste("chr",c(19),sep="") ) ) summary(fit19) fit1<-lm(mirna.cnt~chr.size,data=subset(mirna.size.data,!chr%in%paste("chr",c(1),sep="") ) ) summary(fit1) fit14<-lm(mirna.cnt~chr.size,data=subset(mirna.size.data,!chr%in%paste("chr",c(14),sep="") ) ) summary(fit14) fitX<-lm(mirna.cnt~chr.size,data=subset(mirna.size.data,!chr%in%paste("chr",c(X),sep="") ) ) summary(fitX) fitOST<-lm(mirna.cnt~chr.size,data=subset(mirna.size.data,!chr%in%paste("chr",c("X",14,1,19,17),sep="") ) ) summary(fitOST) 
       
Call:
lm(formula = mirna.cnt ~ chr.size, data = subset(mirna.size.data, 
    !chr %in% paste("chr", c(19), sep = "")))

Residuals:
    Min      1Q  Median      3Q     Max 
-22.384 -10.589  -1.702   2.773  47.495 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept) 1.055e+01  9.098e+00   1.159  0.25942   
chr.size    1.994e-07  6.344e-08   3.143  0.00491 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 17.03 on 21 degrees of freedom
Multiple R-squared: 0.3199,	Adjusted R-squared: 0.2875 
F-statistic: 9.877 on 1 and 21 DF,  p-value: 0.004914 



Call:
lm(formula = mirna.cnt ~ chr.size, data = subset(mirna.size.data, 
    !chr %in% paste("chr", c(1), sep = "")))

Residuals:
    Min      1Q  Median      3Q     Max 
-31.254  -9.817  -2.485  -0.813  59.770 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) 2.543e+01  1.165e+01   2.183   0.0406 *
chr.size    9.802e-08  8.683e-08   1.129   0.2717  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 21.63 on 21 degrees of freedom
Multiple R-squared: 0.05722,	Adjusted R-squared: 0.01232 
F-statistic: 1.275 on 1 and 21 DF,  p-value: 0.2717 



Call:
lm(formula = mirna.cnt ~ chr.size, data = subset(mirna.size.data, 
    !chr %in% paste("chr", c(14), sep = "")))

Residuals:
    Min      1Q  Median      3Q     Max 
-27.463  -9.641  -3.828   0.599  63.573 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)  
(Intercept) 1.874e+01  1.074e+01   1.746   0.0955 .
chr.size    1.469e-07  7.551e-08   1.945   0.0653 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 20.91 on 21 degrees of freedom
Multiple R-squared: 0.1527,	Adjusted R-squared: 0.1123 
F-statistic: 3.784 on 1 and 21 DF,  p-value: 0.06526 

Error in paste("chr", c(X), sep = "") : object 'X' not found
Error in summary(fitX) : 
  error in evaluating the argument 'object' in selecting a method
for function 'summary': Error: object 'fitX' not found


Call:
lm(formula = mirna.cnt ~ chr.size, data = subset(mirna.size.data, 
    !chr %in% paste("chr", c("X", 14, 1, 19, 17), sep = "")))

Residuals:
    Min      1Q  Median      3Q     Max 
-17.971  -6.342   1.938   4.803  11.394 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 7.655e+00  4.908e+00   1.560 0.137241    
chr.size    1.738e-07  3.524e-08   4.931 0.000127 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 8.215 on 17 degrees of freedom
Multiple R-squared: 0.5885,	Adjusted R-squared: 0.5643 
F-statistic: 24.32 on 1 and 17 DF,  p-value: 0.0001266 
fit2<-lm(mirna.cnt~chr.size,data=subset(mirna.size.data,!chr%in%paste("chr",c(19,1,14,"X"),sep="") ) ) summary(fit2) 
       
Call:
lm(formula = mirna.cnt ~ chr.size, data = subset(mirna.size.data, 
    !chr %in% paste("chr", c(19, 1, 14, "X"), sep = "")))

Residuals:
    Min      1Q  Median      3Q     Max 
-20.372  -7.128   1.994   4.388  21.241 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)   
(Intercept) 1.116e+01  5.502e+00   2.028  0.05768 . 
chr.size    1.552e-07  4.017e-08   3.864  0.00114 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 9.547 on 18 degrees of freedom
Multiple R-squared: 0.4534,	Adjusted R-squared: 0.423 
F-statistic: 14.93 on 1 and 18 DF,  p-value: 0.001137 
png(h=800,w=800) layout(matrix(1:4,2,2)) plot(fit2) dev.off() 
       
null device 
          1 
mirna.size.data$model.fit <- NA mirna.size.data[names(fitted(fit2)),"model.fit"] <- fitted(fit2) 
       
mirna.size.data 
       
      mirna.cnt  chr.size   chr model.fit
chr1         75 249250621  chr1        NA
chr10        36 135534747 chr10  32.19372
chr11        41 135006516 chr11  32.11172
chr12        32 133851895 chr12  31.93250
chr13        21 115169878 chr13  29.03260
chr14        69 107349540 chr14        NA
chr15        33 102531392 chr15  27.07081
chr16        28  90354753 chr16  25.18070
chr17        45  81195210 chr17  23.75892
chr18        14  78077248 chr18  23.27494
chr19        91  59128983 chr19        NA
chr2         48 243199373  chr2  48.90584
chr20        30  63025520 chr20  20.93855
chr21        10  48129895 chr21  18.62639
chr22        21  51304566 chr22  19.11918
chr3         44 198022430  chr3  41.89330
chr4         34 191154276  chr4  40.82720
chr5         39 180915260  chr5  39.23786
chr6         23 171115067  chr6  37.71663
chr7         40 159138663  chr7  35.85761
chr8         39 146364022  chr8  33.87468
chr9         37 141213431  chr9  33.07518
chrX         89 155270560  chrX        NA
chrY          0  59373566  chrY  20.37168
library(lattice) xyplot(mirna.cnt + model.fit ~ chr.size,data=mirna.size.data, type=c("p","smooth"),auto.key=T, panel = function(x, y, ...) { panel.xyplot(x, y, ...) ltext(x, y-3, labels = mirna.size.data$chr, cex=0.7) }) dev.off() 
       
null device 
          1 
# Uzyj biblioteki mcmc by dopasowac model liniowy z jednoczesna estymacja outlierow library(mcmc) #?metrop 
       
fit19<-lm(mirna.cnt~chr.size,data=mirna.size.data[ind,] ) summary(fit19) 
       
Call:
lm(formula = mirna.cnt ~ chr.size, data = mirna.size.data[ind, 
    ])

Residuals:
       Min         1Q     Median         3Q        Max 
-4.102e-15 -4.102e-15 -4.102e-15 -4.102e-15  4.513e-14 

Coefficients: (1 not defined because of singularities)
             Estimate Std. Error   t value Pr(>|t|)    
(Intercept) 7.500e+01  4.102e-15 1.828e+16   <2e-16 ***
chr.size           NA         NA        NA       NA    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 1.421e-14 on 11 degrees of freedom
# ind -> Wektor indykacji wl/wyl czy dana zmienna jest outlierem # co to jest a i b ? Rfun<-function(ind){ level <- 0.05 ind_tmp<- ind>=1 #print(ind_tmp) f<-lm(mirna.cnt~chr.size,data=mirna.size.data[ind_tmp,] ) f.sum<-summary(f) f.stat<-f.sum$fstatistic p.value <- 1-pf(f.stat["value"],f.stat["numdf"],f.stat["dendf"]) if(p.value<level){ return(1) } else{ return(-Inf) } } 
       
logRfun<-function(ind){ if (sum(ind[-c(1,2)]>1)>1){ return(c(-Inf,-Inf,-Inf)) } ind_tmp<- ind[-c(1,2)]>0 n<- sum(ind_tmp) a<-ind[1] b<-ind[2] num.out<- NROW(mirna.size.data) -NROW(mirna.size.data[ind_tmp,]) #na liczbę out prior <- 4 pois.prior <- dpois(num.out, 4) loglik<-with(mirna.size.data[ind_tmp,],{ err2_tmp <- ((a+chr.size*b) - mirna.cnt)^2 #wektor bledow kwadratowych err2 <- sum(err2_tmp) # suma błędów kwadratowych sigma2 <- mean(err2) #Estymator wariancji bledu z danych # http://en.wikipedia.org/wiki/Normal_distribution loglik <- -n/2 *(log(2*pi) - log(sigma2)) - 1/(2*sigma2)*err2 loglik }) c(loglik=loglik, pois.prior=pois.prior, score=loglik/n + pois.prior) } 
       
#wektor stanu 0 #na pczątku a i b ind<-c(0.5,1,rep(1,NROW(mirna.size.data))) logRfun(ind) 
       
      loglik   pois.prior        score 
465.90692603   0.01831564  19.43110422 
#autokorelacja (wektor z przesuniętym samym sobą, im dalej tym powinna być większa), wariancja z ogona, 
       
# Unormujmy dane, by wspolczynniki mialy podobna skale mirna.size.data.mcmc<-transform(mirna.size.data, mirna.cnt.sd=sd(mirna.cnt), chr.size.sd=sd(chr.size)) mirna.size.data.mcmc$mirna.cnt<-with(mirna.size.data.mcmc,mirna.cnt/sd(mirna.cnt)) mirna.size.data.mcmc$chr.size<-with(mirna.size.data.mcmc,chr.size/sd(chr.size)) head(mirna.size.data.mcmc) nrep <-10000 
       
      mirna.cnt chr.size   chr model.fit mirna.cnt.sd chr.size.sd
chr1  3.3164470 4.303860  chr1        NA     22.61456    57913269
chr10 1.5918945 2.340306 chr10  32.19372     22.61456    57913269
chr11 1.8129910 2.331185 chr11  32.11172     22.61456    57913269
chr12 1.4150174 2.311247 chr12  31.93250     22.61456    57913269
chr13 0.9286051 1.988661 chr13  29.03260     22.61456    57913269
chr14 3.0511312 1.853626 chr14        NA     22.61456    57913269
f<-logRfun(ind) obj<-function(v)logRfun(v)[3] initial<- ind # na pierwszych dwóch chcemy skakać rzadziej scale<-rep(1,length(ind)) scale[1:2]<-0.3 m<-metrop(obj, initial, nbatch=1e3, scale = scale) #jako wynik to srednia z ostatnich wektorow m 
       
WARNING: Output truncated!  
full_output.txt










$accept
[1] 0

$batch
        [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
[,12] [,13] [,14] [,15] [,16]
   [1,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
   [2,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
   [3,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
   [4,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
   [5,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
   [6,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
   [7,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
   [8,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
   [9,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [10,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [11,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [12,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [13,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [14,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [15,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [16,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [17,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [18,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [19,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [20,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [21,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [22,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [23,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [24,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [25,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [26,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [27,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [28,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [29,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [30,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [31,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [32,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [33,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [34,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [35,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [36,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [37,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [38,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [39,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [40,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [41,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [42,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [43,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [44,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [45,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [46,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1
  [47,]  0.5    1    1    1    1    1    1    1    1     1     1    
1     1     1     1     1

...

[400]  1514423388   754153947  -337335091    35509159 -1609139483  
902399760 -1289942241
[407]  -634623016   565599116   426890807  2094180376 -2041014560 
-550475444   990784105
[414]   602470381  1103627607  -941488942  1128144886  -875528650
-1529081563  1068397160
[421]   295690892    -4488576  1186079057   -89460382  -942361441 
1617078916  1338520547
[428]   865915675  1844917748  -116579212 -1877071614  1448043089 
1189366797   784080155
[435] -2086151449  1039970898    26837919  1254891780 -1090788486  
420959006  -287435594
[442]   418862200  1731099228  -667864411  1493872748 -1800967883   
6648552  -821135074
[449]   612252499  1995174881  1750823484  1116022162 -1809229874 
-492084435 -1630600875
[456]  1954533399  1840165957  -544206688 -1585978056  1617717337 
1831366529  1911574333
[463]     5401146  1998383101  1840087288   488198813 -2088134016 
-476074050    81701000
[470]  -228229566   610310266   640751378  1890058231  1929959392 
1830412356  -314626738
[477]  1466366430   429817077 -1189036763   434244247  1820929198
-1981249149  2055510413
[484]   829652602  -275437991  1213602425  -814845769  1978418249  
834314832   973229233
[491]  -907895720   719974299   428749954 -1022570285  1342674437
-1920950523   985403140
[498]  2060666101 -1338458898   170505345  2073972413  -781802103 
1701148931   826316992
[505]   335748702    85114042 -1838217251  1847070030  1829158664 
1992683755 -1988829453
[512] -1388412686  1484466753  2109810917 -1671110170  1437457306
-1995165221  1592339817
[519]  1667755188  -222148279   329506236   821030402 -1096691124  
697579288  -350331886
[526] -1135119086  1927521087 -1488196456  -262389853   625299014 
-488005693  -110177065
[533] -1163842693  2117066280  -916237568    63678482   749170857
-1327998973  -132897161
[540]  -319079354 -1549161754  -145082720 -1553017440  -771956447  
262176470    28679897
[547] -1523719328  -508918566  1971567435 -1806833956 -1514420041
-1341742435    38907926
[554]  1137912476 -1292563226  -854949098  1808782709  -722487613
-2009411850 -2005320695
[561]  -545230157 -1674208259  -826233660 -1026341203 -1831197874
-1882611626  1284978126
[568]   -63305283 -1602726284  -448300324 -1999046464 -1533232526  
271640415  1207929447
[575] -2114100352   420804678  1728315896  2127857257 -1299862800   
33579269   208761925
[582]  1811626833 -1417904901  -306059299 -1262202409 -1765350233 
1726310475  1676832772
[589]  1617082933    23146513  1885462982  1944700965  1815863652  
684490749   977156848
[596] -1903945500  -298019227  -715558505  -886128516 -2099676749  
272742186 -1811237146
[603]  2082983249  1797119374  1549993910 -1013220620  1078523434 
1241271924  -466687060
[610] -1946398232   632752012   -54852437   738514957   999524903 
-228518140  1917132767
[617]   742531821  1751852589  -594077636  1663972964 -1112395495 
-864188704   -15940703
[624] -1753514552 -1863330220  -175383952

$time
   user  system elapsed 
  0.030   0.000   0.034 

$lud
function (v) 
logRfun(v)[3]

$nbatch
[1] 1000

$blen
[1] 1

$nspac
[1] 1

$scale
 [1] 0.3 0.3 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
[25] 1.0 1.0

$debug
[1] FALSE

attr(,"class")
[1] "mcmc"       "metropolis"
#zastosowac metrop # czy da mu się kazać chodzić tylko po jakichś wartościach ? nrep <-10000 out <- metrop(Rfun, ind, nrep) t <- out$final t>=1 
       
Error in system.time(out <- .Call("metrop", func1, initial,
nbatch, blen,  : 
  log unnormalized density -Inf at initial state
Timing stopped at: 0.01 0 0.007 
Error: object 'out' not found
Error in t >= 1 : 
  comparison (5) is possible only for atomic and list types