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)

# 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 • 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

ADD REPLY
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

ADD REPLY
0
Entering edit mode

yeah I am getting no miRNAs.

ADD REPLY
0
Entering edit mode

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

ADD REPLY
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,
  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
ADD REPLY
0
Entering edit mode

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

ADD REPLY
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?

ADD REPLY
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)
ADD REPLY

Login before adding your answer.

Traffic: 3172 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6