This is a test version of Biostars. For the public version, visit https://www.biostars.org.
microarray miRNA expression data analysis

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)

# read in the data
targets <- read.csv("/media/mdrcubuntu/46B85615B8560439/microarray_text_files/targets.txt", sep="")

# 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)
annotLookup <- read.csv(
  annotfile,
  header = TRUE,
  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 == annotLookup$ProbeID)  # check that annots are aligned
project$genes$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
• 2,572 views • link
updated by Kevin Blighe Kevin Blighe • written by smrutimayipanda • smrutimayipanda

1 answer

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

• 0 views • link
written by Kevin Blighe Kevin Blighe 
 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

• 0 views • link
written by smrutimayipanda • smrutimayipanda

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 views • link
written by Kevin Blighe Kevin Blighe

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

• 0 views • link
written by smrutimayipanda • smrutimayipanda

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,
  adjust.method = 'BH',
  lfc = 0,
  p.value = 1)

head(res)

I never use the decideTests() function.

Directionality is important:

  • Disease - Normal == Numerator, Disease; Denominator, Normal
  • Normal - Disease == Numerator, Normal; Denominator, Disease
• 0 views • link
written by Kevin Blighe Kevin Blighe

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 views • link
written by smrutimayipanda • smrutimayipanda

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)
• 0 views • link
written by Kevin Blighe Kevin Blighe

Log in to answer this question.

More posts like this