Entering edit mode
Minyue Wang
▴
10
@minyue-wang-5197
Last seen 9.6 years ago
> Hi all,
> My name is Amy, I am a masters student in Bioinformatics at North
Carolina
> State University. I am working on a project and I am trying to use
the lumi
> R package for microarray data analysis. I have shown the sample code
here
> and have questions about modifying the sample code for my own data.
>
>
> lumi package in R, example.lumi, the sample data has 8000 features
and 4
> samples
>
> I have highlighted the code I have questions on in red, my data has
4
> different types of samples, each repeated 6 times, so a total of 24
samples
> and about 48,000 rows. how should I identify my sampleType in my
case? also
> what does colnames(design) <- c('100:0', '95:5-100:0') do, which
columns
> exactly does it take into consideration? Thanks!
>
>
> so the sample code i'm trying to follow is below:
>
> ###################################################
>
> ### code chunk number 30: filtering
>
> ###################################################
>
> presentCount <- detectionCall(example.lumi)
>
> selDataMatrix <- dataMatrix[presentCount > 0,]
>
> probeList <- rownames(selDataMatrix)
>
>
>
>
>
> ###################################################
>
> ### code chunk number 31: Identify differentially expressed genes
>
> ###################################################
>
> ## Specify the sample type
>
> sampleType <- c('100:0', '95:5', '100:0', '95:5')
>
> if (require(limma)) {
>
> ## compare '95:5' and '100:0'
>
> design <- model.matrix(~ factor(sampleType))
>
> colnames(design) <- c('100:0', '95:5-100:0')
>
> 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='95:5-100:0', 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]
>
> }
>
--
- Amy W.
[[alternative HTML version deleted]]
