R version 2.5.0 (2007-04-23) Copyright (C) 2007 The R Foundation for Statistical Computing ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > if(!exists("baseenv", mode="function")) baseenv <- function() NULL > options(STERM='iESS', editor='emacsclient') > setwd("/usr/local/data/supplements/sfp/data/barley") > library(limma) > source("../SFP_function.R") > load("barleyRNA.pm.RData") load("pheno.RData") > > allSNP <- read.csv("SNPtable.csv",as.is=T) > temp <- which(duplicated(allSNP$Probe)) > allSNP <- allSNP[-temp,] > str(allSNP <- allSNP[allSNP$SNPgenotype != "GP,MX",]) 'data.frame': 2601 obs. of 5 variables: $ Probe : chr "Contig10521_at3" "Contig10838_at11" "Contig11037_at3" "Contig12503_s_at10" ... $ SNPgenotype : chr "GP" "GP" "GP" "GP" ... $ dataset : chr "EST" "EST" "EST" "EST" ... $ SNPscore : chr "12" NA "9+1" "6" ... $ Probe.sequence: chr "GAAGCCGGTTTCGATCTGAGCTAAT" "GATTGTGTTGCCGTTCATTGAATAA" "GGTACTACCAGCAGGATAATCCCAA" "ACAGGAAGTGGTGCTTTGTTGTGCC" ... > pm.medpolishres <- get_medpolish_res(barleyRNA.pm,ps.n=11) Loading required package: preprocessCore Attaching package: 'preprocessCore' The following object(s) are masked from package:affy : normalize.quantiles, normalize.quantiles.robust colnames(pm.medpolishres) <- pheno[,2] 4870(cl <- factor(substring(colnames(pm.medpolishres),1,1))) 5541design <- model.matrix(~0+cl);colnames(design) <- sub('cl','',colnames(design));contrast.matrix <- makeContrasts(M-G,levels=design) 6212fit <- do_limma(pm.medpolishres,design,contrast.matrix) 7911(fixnum <- summary(fit.data <- decideTests(fit,p=0.001)));table(fit.data!=0) 19730save(fit.data,file="fit.data.RData") 22801 probeSNP <- match(as.character(allSNP$Probe),sfpName) > fixnum <- 10504 str(my.pqvalue <- p.adjust(fit$p.value[,1],'fdr')) (thresh <-quantile(my.pqvalue, fixnum/length(my.pqvalue))) > [1] G G G M M M G G G M M M G G G M M M G G G M M M G G G M M M G G M M M M Levels: G M > > fit.data <- matrix(rep(0,length(my.pqvalue))) fit.data <- new("TestResults", array(0, dim(fit$t))) dimnames(fit.data) <- dimnames(fit$t) rownames(fit.data) <- sub('--.*$','',rownames(fit.data)) fit.data[my.pqvalue0] <- 1 fit.data[my.pqvalue M - G -1 4950 0 241317 1 4544 FALSE TRUE 241317 9494 > > Error in inherits(x, "factor") : object "sfpName" not found > > Named num [1:250811] 0.977 0.990 0.995 0.828 0.789 ... - attr(*, "names")= chr [1:250811] "1200459_Reg_88-1740_at1--80--545" "1200459_Reg_88-1740_at2--108--171" "1200459_Reg_88-1740_at3--51--35" "1200459_Reg_88-1740_at4--154--287" ... > 4.188014% 0.002805159 > > > > > > > M - G -1 5440 0 240307 1 5064 > save(fit.data,file="fit.data.medianpolish_10504.RData",compress=T) > table(snp.geno <- sfp.check <- fit.data[rownames(fit.data)%in%allSNP$Probe,]) -1 0 1 142 2224 235 > snp.geno[snp.geno!=0] <- 0 > snp.geno[names(snp.geno)%in%allSNP$Probe[allSNP$SNPgenotype=='GP']] <- 1 > snp.geno[names(snp.geno)%in%allSNP$Probe[allSNP$SNPgenotype=='MX']] <- -1 > sfp.check <- fit.data[match(names(snp.geno),rownames(fit.data)),] > print(y <- table(snp.geno,sfp.check)) sfp.check snp.geno -1 0 1 -1 123 44 11 0 11 2129 60 1 8 51 164 > print(summary(y)) Number of cases in table: 2601 Number of factors: 2 Test for independence of all factors: Chisq = 2758.6, df = 4, p-value = 0 > cat("Sensitivity", sum(c(y[c(1,3),-2]))/sum(c(y[c(1,3),])), "\n" ) Sensitivity 0.7630923 > cat("False Negative Rate", 1-sum(c(y[c(1,3),-2]))/sum(c(y[c(1,3),])), "\n") False Negative Rate 0.2369077 > cat("False positive Rate", sum(c(y[2,-2]))/sum(c(y[2,])), "\n") False positive Rate 0.03227273 > cat("False sequence polymorphism Rate", sum(c(y[2,-2]))/sum(c(y[,c(1,3)])), "\n") False sequence polymorphism Rate 0.1883289 > q() Save workspace image? [y/n/c]: n Process R finished at Sat May 3 17:03:27 2008