Entering edit mode
Kripa R
▴
180
@kripa-r-4482
Last seen 10.2 years ago
Hello, I'm trying to analyze 11 samples from illumina HT12 chips. In
the lumi package I set Although I set convert NuID and input
annotation to be true, after using limma my top table does not show
the gene symbol or name.
Loading required package: annotate
ID geneSymbol geneName logFC AveExpr t
4894 2230037 NA NA 11.93955 12.13396 100.74365
1211 610112 NA NA 11.55668 11.63635 99.67799
I was wondering how I'd be able to add this information and export it.
Any help would be greatly appreciated!
> sessionInfo()
R version
2.12.2 (2011-02-25)
Platform:
i386-pc-mingw32/i386 (32-bit)
locale:
[1]
LC_COLLATE=English_Canada.1252
LC_CTYPE=English_Canada.1252
[3]
LC_MONETARY=English_Canada.1252 LC_NUMERIC=C
[5]
LC_TIME=English_Canada.1252
attached base
packages:
[1]
stats graphics grDevices utils datasets
methods base
other
attached packages:
[1]
annotate_1.28.0
lumiHumanAll.db_1.12.0 org.Hs.eg.db_2.4.6
[4]
RSQLite_0.9-2 DBI_0.2-5 AnnotationDbi_1.12.0
[7]
limma_3.6.9 lumi_2.2.1 Biobase_2.10.0
loaded via a
namespace (and not attached):
[1] affy_1.28.0 affyio_1.18.0 grid_2.12.0
[4] hdrcde_2.15 KernSmooth_2.23-4 lattice_0.19-13
[7] MASS_7.3-8 Matrix_0.999375-44 methylumi_1.6.1
[10]
mgcv_1.6-2 nlme_3.1-97 preprocessCore_1.12.0
[13]
tools_2.12.0 xtable_1.5-6
#------------------------------------------------------------#
###HT-12 illumina data.
11samples, 3 groups
#group1:AJ
#group2:BCDH
#group3:EFGIKL
#------------------------------------------------------------#
###PREAMBLE###############################
library("lumi");
date <- Sys.Date();
####################################################
setwd("\\\\rinas1p2/users/ramank/Illumina/NugenAnalysis")
#load raw data
data<-
"FinalReportA-K.txt";
#read data
x.lumi<-lumiR(fileName=data,
sep="\t",
detectionTh=0.01,
na.rm=TRUE,
convertNuID=TRUE,
lib.mapping=NULL,
dec='.',
parseColumnName=TRUE,
checkDupId=TRUE,
QC=TRUE,
columnNameGrepPattern=list(exprs='AVG_Signal',
se.exprs='BEAD_STDERR',
detection='Detection
Pval',
beadNum='Avg_NBEADS'),
inputAnnotation=TRUE,
annotationColumn=c('ACCESSION',
'SYMBOL',
'PROBE_START',
'CHROMOSOME',
'PROBE_CHR_ORIENTATION',
'PROBE_COORDINATES',
'DEFINITION'),
verbose=TRUE);
#summary of the quality
control
summary(x.lumi, 'QC');
##????how to write this to
a txt file?????##
#########2 PRE PROCESSING
METHODS (LUMI)################## ############
###no background
correction since no control,
#variance stabilization =
log2
#normalization = quantile
noBg_log_quantile <-
lumiExpresso(
x.lumi,
bg.correct
= TRUE,
bgcorrect.param
= list(method='none'),
variance.stabilize
= TRUE,
varianceStabilize.param
= list(method='log2'),
normalize
= TRUE,
normalize.param
= list(method='quantile'),
QC.evaluation
= TRUE,
QC.param = list(), verbose = TRUE)
summary(noBg_log_quantile,'QC')
#making matrix of the data
lumi.N.Q <-
noBg_log_quantile
dataMatrix <-
exprs(lumi.N.Q);
#filtering
presentCount <-
detectionCall(x.lumi);
selDataMatrix <-
dataMatrix[presentCount > 0,];
probeList <-
rownames(selDataMatrix);
####################################################################
## LIMMA
####################################################################
SampleType <-
c('1','2','2','2','3','3','3','2','3','1','3')
##????above line correct?????##
if (require(limma)) {
## compare '1'
and '2' and '3'
design <- model.matrix(~ factor(sampleType))
colnames(design) <- c('1', '2', '3')
fit <- lmFit(selDataMatrix, design)
fit <- eBayes(fit)
## Add gene symbols to gene properties
if (require(lumiHumanAll.db) &
require(annotate)) {
geneSymbol <-
getSYMBOL(probeList, 'lumiHumanAll.db')
geneName <-
sapply(lookUp(probeList, 'lumiHumanAll.db', 'GENENAME'), function(x)
x[1])
fit$genes <- data.frame(ID=
probeList, geneSymbol=geneSymbol, geneName=geneName,
stringsAsFactors=FALSE)
}
## print the top 10 genes
print(topTable(fit, coef='1', adjust='fdr',
number=10))
## get significant gene list with FDR adjusted
p.values less than 0.01
p.adj <- p.adjust(fit$p.value[,2])
sigGene.adj <- probeList[ p.adj < 0.01]
## without FDR adjustment
sigGene <- probeList[ fit$p.value[,2] < 0.01]
} ;
.kripa
[[alternative HTML version deleted]]