Return to Question

Following are the codes I used:

source("http://www.bioconductor.org/biocLite.R")
biocLite()
library(affy)
library(limma)
library(hexbin)
mydata<-ReadAffy()
eset<-rma(mydata)

## Load biomaRt package
library(biomaRt)

## Specify which "mart" (i.e., source of genetic data) that you want to use
ensembl <- useMart("ensembl")
ensembl <- useDataset("hsapiens_gene_ensembl", mart = ensembl)

## You can then ask the system what attributes are available for download
listAttributes(ensembl)

Download microRNA data
miRNA <- getBM(c("mirbase_id", "ensembl_gene_id", "start_position", "chromosome_name"), filters = c("with_mirbase"), values = list(TRUE), mart = ensembl)

##Check how much we downloaded
dim(miRNA)

##Peak at the head of our data
head(miRNA)

## Check which chromosomes are contributing to our data
table(miRNA$chromosome_name)

grep("hsa", featureNames(eset))->hsa
eset.hsa<-eset[hsa,]

group<-factor(c(rep(1,16), rep(2,20), rep(3,19)))
design<-model.matrix(~0+group)
fit<-lmFit(eset.hsa, design)
contrast.matrix <- makeContrasts(group2-group1, group3-group2, group3-group1, levels=design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
topTable(fit2, coef=1, adjust="none", sort.by="logFC", number=10)

topTable(fit2, coef=1, adjust="fdr", sort.by="B", number=100) 

write.table(topTable(fit2, coef=1, adjust="fdr", sort.by="B", number=100), file="Group1.xls", row.names=F, sep="\t")
results <- decideTests(fit2, p.value=0.05); vennDiagram(results)
x <- topTable(fit2, coef=1, adjust="fdr", sort.by="P", number=50000); y <- x[x$adj.P.Val < 0.05,]; y; print("Number of genes in this list:"); length(y$ID)

QC
library(affyQCReport)
QCReport(mydata, file="ExampleQC.pdf")

deg <- AffyRNAdeg(mydata); summaryAffyRNAdeg(deg); plotAffyRNAdeg(deg)
image(mydata[ ,1])
hist(mydata[ ,1:2])
hist(log2(pm(mydata[,1])), breaks=100, col="blue")
boxplot(mydata,col="red")
boxplot(data.frame(exprs(eset.hsa)),col="blue", main="Normalized Data")
mva.pairs(pm(mydata)[,c(1,4)])
mva.pairs(exprs(eset.hsa)[,c(1,4)])


my_fct <- function(x) hclust(x, method="complete")
heatmap(as.matrix(2^exprs(eset.hsa)[1:40,]), col = cm.colors(256), hclustfun=my_fct)
heatmap(as.matrix(2^exprs(eset.hsa)[1:40,]), col = cm.colors(256), hclustfun=my_fct)
heatmap(as.matrix(2^exprs(eset.hsa)[1:10,]), col = cm.colors(256), hclustfun=my_fct)
heatmap(as.matrix(2^exprs(eset.hsa)[1:10,]), col=topo.colors(100), hclustfun=my_fct)
heatmap(as.matrix(2^exprs(eset.hsa)[1:10,]), col=topo.colors(75), hclustfun=my_fct)
heatmap(as.matrix(2^exprs(eset.hsa)[1:10,]), col=topo.colors(75), scale="none", hclustfun=my_fct)
heatmap(as.matrix(2^exprs(eset.hsa)[1:10,]), col=topo.colors(100), scale="none", hclustfun=my_fct)
heatmap(as.matrix(2^exprs(eset.hsa)[1:10,]), col=topo.colors(100), scale="none", key=TRUE, symkey=FALSE, density.info="none", trace="none", cexRow=0.5, hclustfun=my_fct)


selected<-p.adjust(fit2$p.value[,1])<0.5
selected<-fit2$p.value[,1]<0.005
esetSel<-eset.hsa[selected,]
esetSel
heatmap(esetSel)
my_fct <- function(x) hclust(x, method="complete")
heatmap(as.matrix(2^exprs(esetSel)), col =topo.colors(256))
heatmap(as.matrix(2^exprs(esetSel)), col =redgreen(75))
library("gplots")
heatmap.2(as.matrix(2^exprs(eset.hsa)[1:10,]), col=topo.colors(100), scale="none", key=TRUE, symkey=FALSE, density.info="none", trace="none", cexRow=0.5, hclustfun=my_fct)
heatmap.2(as.matrix(2^exprs(eset.hsa)[1:10,]), col=redgreen(75), scale="row", key=TRUE, symkey=FALSE, density.info="none", trace="none", cexRow=0.5, hclustfun=my_fct)
heatmap.2(as.matrix(2^exprs(eset.hsa)[1:10,]), col=redgreen(75), scale="row", key=TRUE, symkey=FALSE, density.info="none", trace="none", cexRow=0.5, hclustfun=my_fct)

source("http://www.bioconductor.org/biocLite.R") biocLite() library(affy) library(limma) library(hexbin) mydata<-ReadAffy() eset<-rma(mydata)

Load biomaRt package

library(biomaRt)

Specify which "mart" (i.e., source of genetic data) that you want to use

ensembl <- useMart("ensembl") ensembl <- useDataset("hsapiens_gene_ensembl", mart = ensembl)

You can then ask the system what attributes are available for download

listAttributes(ensembl)

Download microRNA data miRNA <- getBM(c("mirbase_id", "ensembl_gene_id", "start_position", "chromosome_name"), filters = c("with_mirbase"), values = list(TRUE), mart = ensembl)

##Check how much we downloaded dim(miRNA)

##Peak at the head of our data head(miRNA)

Check which chromosomes are contributing to our data

table(miRNA$chromosome_name)

grep("hsa", featureNames(eset))->hsa eset.hsa<-eset[hsa,]

group<-factor(c(rep(1,16), rep(2,20), rep(3,19))) design<-model.matrix(~0+group) fit<-lmFit(eset.hsa, design) contrast.matrix <- makeContrasts(group2-group1, group3-group2, group3-group1, levels=design) fit2 <- contrasts.fit(fit, contrast.matrix) fit2 <- eBayes(fit2) topTable(fit2, coef=1, adjust="none", sort.by="logFC", number=10)

topTable(fit2, coef=1, adjust="fdr", sort.by="B", number=100)

write.table(topTable(fit2, coef=1, adjust="fdr", sort.by="B", number=100), file="Group1.xls", row.names=F, sep="\t") results <- decideTests(fit2, p.value=0.05); vennDiagram(results) x <- topTable(fit2, coef=1, adjust="fdr", sort.by="P", number=50000); y <- x[x$adj.P.Val < 0.05,]; y; print("Number of genes in this list:"); length(y$ID)

QC library(affyQCReport) QCReport(mydata, file="ExampleQC.pdf")

deg <- AffyRNAdeg(mydata); summaryAffyRNAdeg(deg); plotAffyRNAdeg(deg) image(mydata[ ,1]) hist(mydata[ ,1:2]) hist(log2(pm(mydata[,1])), breaks=100, col="blue") boxplot(mydata,col="red") boxplot(data.frame(exprs(eset.hsa)),col="blue", main="Normalized Data") mva.pairs(pm(mydata)[,c(1,4)]) mva.pairs(exprs(eset.hsa)[,c(1,4)])

my_fct <- function(x) hclust(x, method="complete") heatmap(as.matrix(2^exprs(eset.hsa)[1:40,]), col = cm.colors(256), hclustfun=my_fct) heatmap(as.matrix(2^exprs(eset.hsa)[1:40,]), col = cm.colors(256), hclustfun=my_fct) heatmap(as.matrix(2^exprs(eset.hsa)[1:10,]), col = cm.colors(256), hclustfun=my_fct) heatmap(as.matrix(2^exprs(eset.hsa)[1:10,]), col=topo.colors(100), hclustfun=my_fct) heatmap(as.matrix(2^exprs(eset.hsa)[1:10,]), col=topo.colors(75), hclustfun=my_fct) heatmap(as.matrix(2^exprs(eset.hsa)[1:10,]), col=topo.colors(75), scale="none", hclustfun=my_fct) heatmap(as.matrix(2^exprs(eset.hsa)[1:10,]), col=topo.colors(100), scale="none", hclustfun=my_fct) heatmap(as.matrix(2^exprs(eset.hsa)[1:10,]), col=topo.colors(100), scale="none", key=TRUE, symkey=FALSE, density.info="none", trace="none", cexRow=0.5, hclustfun=my_fct)

selected<-p.adjust(fit2$p.value[,1])<0.5 selected<-fit2$p.value[,1]<0.005 esetSel<-eset.hsa[selected,] esetSel heatmap(esetSel) my_fct <- function(x) hclust(x, method="complete") heatmap(as.matrix(2^exprs(esetSel)), col =topo.colors(256)) heatmap(as.matrix(2^exprs(esetSel)), col =redgreen(75)) library("gplots") heatmap.2(as.matrix(2^exprs(eset.hsa)[1:10,]), col=topo.colors(100), scale="none", key=TRUE, symkey=FALSE, density.info="none", trace="none", cexRow=0.5, hclustfun=my_fct) heatmap.2(as.matrix(2^exprs(eset.hsa)[1:10,]), col=redgreen(75), scale="row", key=TRUE, symkey=FALSE, density.info="none", trace="none", cexRow=0.5, hclustfun=my_fct) heatmap.2(as.matrix(2^exprs(eset.hsa)[1:10,]), col=redgreen(75), scale="row", key=TRUE, symkey=FALSE, density.info="none", trace="none", cexRow=0.5, hclustfun=my_fct)

source("http://www.bioconductor.org/biocLite.R")
biocLite()
library(affy)
library(limma)
library(hexbin)
mydata<-ReadAffy()
eset<-rma(mydata)

## Load biomaRt package
library(biomaRt)

## Specify which "mart" (i.e., source of genetic data) that you want to use
ensembl <- useMart("ensembl")
ensembl <- useDataset("hsapiens_gene_ensembl", mart = ensembl)

## You can then ask the system what attributes are available for download
listAttributes(ensembl)

Download microRNA data
miRNA <- getBM(c("mirbase_id", "ensembl_gene_id", "start_position", "chromosome_name"), filters = c("with_mirbase"), values = list(TRUE), mart = ensembl)

##Check how much we downloaded
dim(miRNA)

##Peak at the head of our data
head(miRNA)

## Check which chromosomes are contributing to our data
table(miRNA$chromosome_name)

grep("hsa", featureNames(eset))->hsa
eset.hsa<-eset[hsa,]

group<-factor(c(rep(1,16), rep(2,20), rep(3,19)))
design<-model.matrix(~0+group)
fit<-lmFit(eset.hsa, design)
contrast.matrix <- makeContrasts(group2-group1, group3-group2, group3-group1, levels=design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
topTable(fit2, coef=1, adjust="none", sort.by="logFC", number=10)

topTable(fit2, coef=1, adjust="fdr", sort.by="B", number=100) 

write.table(topTable(fit2, coef=1, adjust="fdr", sort.by="B", number=100), file="Group1.xls", row.names=F, sep="\t")
results <- decideTests(fit2, p.value=0.05); vennDiagram(results)
x <- topTable(fit2, coef=1, adjust="fdr", sort.by="P", number=50000); y <- x[x$adj.P.Val < 0.05,]; y; print("Number of genes in this list:"); length(y$ID)

QC
library(affyQCReport)
QCReport(mydata, file="ExampleQC.pdf")

deg <- AffyRNAdeg(mydata); summaryAffyRNAdeg(deg); plotAffyRNAdeg(deg)
image(mydata[ ,1])
hist(mydata[ ,1:2])
hist(log2(pm(mydata[,1])), breaks=100, col="blue")
boxplot(mydata,col="red")
boxplot(data.frame(exprs(eset.hsa)),col="blue", main="Normalized Data")
mva.pairs(pm(mydata)[,c(1,4)])
mva.pairs(exprs(eset.hsa)[,c(1,4)])


my_fct <- function(x) hclust(x, method="complete")
heatmap(as.matrix(2^exprs(eset.hsa)[1:40,]), col = cm.colors(256), hclustfun=my_fct)
heatmap(as.matrix(2^exprs(eset.hsa)[1:40,]), col = cm.colors(256), hclustfun=my_fct)
heatmap(as.matrix(2^exprs(eset.hsa)[1:10,]), col = cm.colors(256), hclustfun=my_fct)
heatmap(as.matrix(2^exprs(eset.hsa)[1:10,]), col=topo.colors(100), hclustfun=my_fct)
heatmap(as.matrix(2^exprs(eset.hsa)[1:10,]), col=topo.colors(75), hclustfun=my_fct)
heatmap(as.matrix(2^exprs(eset.hsa)[1:10,]), col=topo.colors(75), scale="none", hclustfun=my_fct)
heatmap(as.matrix(2^exprs(eset.hsa)[1:10,]), col=topo.colors(100), scale="none", hclustfun=my_fct)
heatmap(as.matrix(2^exprs(eset.hsa)[1:10,]), col=topo.colors(100), scale="none", key=TRUE, symkey=FALSE, density.info="none", trace="none", cexRow=0.5, hclustfun=my_fct)


selected<-p.adjust(fit2$p.value[,1])<0.5
selected<-fit2$p.value[,1]<0.005
esetSel<-eset.hsa[selected,]
esetSel
heatmap(esetSel)
my_fct <- function(x) hclust(x, method="complete")
heatmap(as.matrix(2^exprs(esetSel)), col =topo.colors(256))
heatmap(as.matrix(2^exprs(esetSel)), col =redgreen(75))
library("gplots")
heatmap.2(as.matrix(2^exprs(eset.hsa)[1:10,]), col=topo.colors(100), scale="none", key=TRUE, symkey=FALSE, density.info="none", trace="none", cexRow=0.5, hclustfun=my_fct)
heatmap.2(as.matrix(2^exprs(eset.hsa)[1:10,]), col=redgreen(75), scale="row", key=TRUE, symkey=FALSE, density.info="none", trace="none", cexRow=0.5, hclustfun=my_fct)
heatmap.2(as.matrix(2^exprs(eset.hsa)[1:10,]), col=redgreen(75), scale="row", key=TRUE, symkey=FALSE, density.info="none", trace="none", cexRow=0.5, hclustfun=my_fct)

Actually I have miRNA data from affymetrix chip. Data is classified into 3 groups. I used limma package for analysis. I did RMA normalization, eBayes etc. But I did not get any miRNA significatly expressed in any group. When I posed questions about it on net somebody suggested me to PCA analysis etc. to remove outliners and go ahead with normal analysis as I did.

Now I will appreciate it some body can guide me from the begining

Following are the codes I used:

source("http://www.bioconductor.org/biocLite.R") biocLite() library(affy) library(limma) library(hexbin) mydata<-ReadAffy() eset<-rma(mydata)

Load biomaRt package

library(biomaRt)

Specify which "mart" (i.e., source of genetic data) that you want to use

ensembl <- useMart("ensembl") ensembl <- useDataset("hsapiens_gene_ensembl", mart = ensembl)

You can then ask the system what attributes are available for download

listAttributes(ensembl)

Download microRNA data miRNA <- getBM(c("mirbase_id", "ensembl_gene_id", "start_position", "chromosome_name"), filters = c("with_mirbase"), values = list(TRUE), mart = ensembl)

##Check how much we downloaded dim(miRNA)

##Peak at the head of our data head(miRNA)

Check which chromosomes are contributing to our data

table(miRNA$chromosome_name)

grep("hsa", featureNames(eset))->hsa eset.hsa<-eset[hsa,]

group<-factor(c(rep(1,16), rep(2,20), rep(3,19))) design<-model.matrix(~0+group) fit<-lmFit(eset.hsa, design) contrast.matrix <- makeContrasts(group2-group1, group3-group2, group3-group1, levels=design) fit2 <- contrasts.fit(fit, contrast.matrix) fit2 <- eBayes(fit2) topTable(fit2, coef=1, adjust="none", sort.by="logFC", number=10)

topTable(fit2, coef=1, adjust="fdr", sort.by="B", number=100)

write.table(topTable(fit2, coef=1, adjust="fdr", sort.by="B", number=100), file="Group1.xls", row.names=F, sep="\t") results <- decideTests(fit2, p.value=0.05); vennDiagram(results) x <- topTable(fit2, coef=1, adjust="fdr", sort.by="P", number=50000); y <- x[x$adj.P.Val < 0.05,]; y; print("Number of genes in this list:"); length(y$ID)

QC library(affyQCReport) QCReport(mydata, file="ExampleQC.pdf")

deg <- AffyRNAdeg(mydata); summaryAffyRNAdeg(deg); plotAffyRNAdeg(deg) image(mydata[ ,1]) hist(mydata[ ,1:2]) hist(log2(pm(mydata[,1])), breaks=100, col="blue") boxplot(mydata,col="red") boxplot(data.frame(exprs(eset.hsa)),col="blue", main="Normalized Data") mva.pairs(pm(mydata)[,c(1,4)]) mva.pairs(exprs(eset.hsa)[,c(1,4)])

my_fct <- function(x) hclust(x, method="complete") heatmap(as.matrix(2^exprs(eset.hsa)[1:40,]), col = cm.colors(256), hclustfun=my_fct) heatmap(as.matrix(2^exprs(eset.hsa)[1:40,]), col = cm.colors(256), hclustfun=my_fct) heatmap(as.matrix(2^exprs(eset.hsa)[1:10,]), col = cm.colors(256), hclustfun=my_fct) heatmap(as.matrix(2^exprs(eset.hsa)[1:10,]), col=topo.colors(100), hclustfun=my_fct) heatmap(as.matrix(2^exprs(eset.hsa)[1:10,]), col=topo.colors(75), hclustfun=my_fct) heatmap(as.matrix(2^exprs(eset.hsa)[1:10,]), col=topo.colors(75), scale="none", hclustfun=my_fct) heatmap(as.matrix(2^exprs(eset.hsa)[1:10,]), col=topo.colors(100), scale="none", hclustfun=my_fct) heatmap(as.matrix(2^exprs(eset.hsa)[1:10,]), col=topo.colors(100), scale="none", key=TRUE, symkey=FALSE, density.info="none", trace="none", cexRow=0.5, hclustfun=my_fct)

selected<-p.adjust(fit2$p.value[,1])<0.5 selected<-fit2$p.value[,1]<0.005 esetSel<-eset.hsa[selected,] esetSel heatmap(esetSel) my_fct <- function(x) hclust(x, method="complete") heatmap(as.matrix(2^exprs(esetSel)), col =topo.colors(256)) heatmap(as.matrix(2^exprs(esetSel)), col =redgreen(75)) library("gplots") heatmap.2(as.matrix(2^exprs(eset.hsa)[1:10,]), col=topo.colors(100), scale="none", key=TRUE, symkey=FALSE, density.info="none", trace="none", cexRow=0.5, hclustfun=my_fct) heatmap.2(as.matrix(2^exprs(eset.hsa)[1:10,]), col=redgreen(75), scale="row", key=TRUE, symkey=FALSE, density.info="none", trace="none", cexRow=0.5, hclustfun=my_fct) heatmap.2(as.matrix(2^exprs(eset.hsa)[1:10,]), col=redgreen(75), scale="row", key=TRUE, symkey=FALSE, density.info="none", trace="none", cexRow=0.5, hclustfun=my_fct)