microarray miRNA expression data analysis
1
0
Entering edit mode
10 weeks ago

I wrote a script on how to analyze the microarray-based miRNA expression data. Here is my code:

# general config
baseDir <- '.'
annotfile <- 'mirbase_genelist.tsv'
setwd(baseDir)
options(scipen = 99)
require(limma)

# retain information about background via gIsWellAboveBG
project <- read.maimages(targets,source = 'agilent.median',green.only = TRUE, other.columns = c("gIsWellAboveBG","gIsSaturated","gIsFeatPopnOL","gIsFeatNonUnifOL")
)
# annotate the probes (gene annotation)
annotfile,
sep = '\t',
stringsAsFactors = FALSE)
colnames(annotLookup)[1] <- 'ProbeID'
annotLookup <- annotLookup[which(annotLookup$ProbeID %in% project$genes$ProbeName),] annotLookup <- annotLookup[match(project$genes$ProbeName, annotLookup$ProbeID),]
table(project$genes$ProbeName == annotLookupProbeID) # check that annots are aligned projectgenes$ProbeID <- annotLookup$ProbeID
project$genes$wikigene_description <- annotLookup$wikigene_description project$genes$ensembl_gene_id <- annotLookup$ensembl_gene_id
project$genes$targetID <- annotLookup$TargetID project$genes$mirbase_accession <- annotLookup$mirbase_accession
project$genes$external_gene_name <- annotLookup$external_gene_name # perform background correction on the fluorescent intensities (subtract the background) project.bgcorrect <- backgroundCorrect(project, method = 'normexp', offset=16) # normalize the data with the 'quantile' method project.bgcorrect.norm <- normalizeBetweenArrays(project.bgcorrect, method = 'quantile') # filter out control probes, those with no symbol, and those that fail (probe filter) Control <- project.bgcorrect.norm$genes$ControlType==1L NoSymbol <- is.na(project.bgcorrect.norm$genes$ProbeName) project.bgcorrect.norm.filt <- project.bgcorrect.norm[!Control & !NoSymbol,] # for replicate probes, replace values with the mean # ID is used to identify the replicates project.bgcorrect.norm.filt.mean <- avereps(project.bgcorrect.norm.filt,ID = project.bgcorrect.norm.filt$genes$ProbeName) #The normalised expression matrix is then contained in project.bgcorrect.norm.filt.mean$E

#Create the study design and comparison model(DEG filter)
Group <- factor(targets$Treatment, levels=c("Normal","Diseased")) design <- model.matrix(~0+Group) colnames(design) <- c("Normal","DiseasedvsNormal") #Checking the experimental design design #Applying the empirical Bayes method with t-test fit <- lmFit(project.bgcorrect.norm.filt.mean, design) contrast.matrix <- makeContrasts(Normal-DiseasedvsNormal, levels=design) fit2 <- contrasts.fit(fit, contrast.matrix) fit2 <- eBayes(fit2) summary(decideTests(fit2)) probeset.list <- topTable(fit2,coef=1, number=Inf, adjust.method ="BH", lfc=1, p.value = 0.05) results <- decideTests(fit2, method="global") vennDiagram(results, include=c("up","down"))  After running this, I am not getting any differentially expressed miRNAs. Can anyone please tell me whether this code is correct or having some error? I will fix it, but it would be very helpful if someone can help me regarding this. Thank you. microarray R • 458 views ADD COMMENT 0 Entering edit mode 10 weeks ago Hi smrutimayipanda, The code looks okay. The only key parts to verify would be the design matrix and the coefficient being selected when running topTable(). Could you please share the output of model.matrix(~0+Group)? Need to see this because the contrast names that you later select seem unusual, i.e., 'Normal' and 'DiseasedvsNormal' When running topTable(), I also output everything by selecting number = Inf and lfc = 0 and p.value = 1 Kevin ADD COMMENT 0 Entering edit mode  GroupNormal GroupDiseased 1 1 0 2 1 0 3 1 0 4 1 0 5 0 1 6 0 1 7 0 1 8 0 1 attr(,"assign") [1] 1 1 attr(,"contrasts") attr(,"contrasts")$Group
[1] "contr.treatment"


I got this after running model.matrix

1
Entering edit mode

Cool, in this case, I would change this line:

colnames(design) <- c("Normal","DiseasedvsNormal")


to:

colnames(design) <- c("Normal","Diseased")


It looks like this will not make any difference to the results, though. It can be that there are genuinely no statistically significantly differentially expressed miRNAs

0
Entering edit mode

yeah I am getting no miRNAs.

0
Entering edit mode

then what should be the contrast matrix- normal-diseased?

0
Entering edit mode

Yes, probably then need:

contrast.matrix <- makeContrasts(
Disease_vs_Normal = Disease - Normal,
levels = design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
summary(decideTests(fit2))
res <- topTable(fit2, coef = 'Disease_vs_Normal',
number = Inf,
lfc = 0,
p.value = 1)



I never use the decideTests() function.

Directionality is important:

• Disease - Normal == Numerator, Disease; Denominator, Normal
• Normal - Disease == Numerator, Normal; Denominator, Disease
0
Entering edit mode

Thank you so much kevin. It helped me a lot.

0
Entering edit mode

This code gave me all upregulated miRNAs. Why is it? It should give me some up and downregulated miRNAs. When I am using Disease - Normal, it is giving up and Normal - Disease, it is giving me down. How it can be rectified?

0
Entering edit mode

This code gave me all upregulated miRNAs. Why is it? It should give me some up and downregulated miRNAs.

Perhaps it is the genuine result. Please show some of the output table, if you can?

When I am using Disease - Normal, it is giving up and Normal - Disease, it is giving me down. How it can be rectified?

There is no problem here. This is expected behaviour.

• Disease - Normal equals Disease divided by Normal (for fold change calculations)
• Normal - Disease equals Normal divided by Disease (for fold change calculations)