Hi guys, I know this question mqybe look silly but it is taking me a lot of time. I have downoaded some miRNA expression data from the AML samples available on the webstie TCGA (Arround 18 samples). I calsified my samples to groups according to some criteria (aproximatly 4 samples per group), and did a simple diffrent expression analysis sing DESeq Here is the used code:
get.DiffrentlyExpressed<-function(counts,groups){
library("DESeq")
my.design <- data.frame(row.names = colnames( counts ),condition = groups)
conds <- factor(my.design$condition)
cds <- newCountDataSet( counts, conds )
cds <- estimateSizeFactors( cds )
cds <- estimateDispersions( cds )
sizeFactors( cds )
groups<-factor(groups);
result <- nbinomTest( cds, levels(groups)[1],levels(groups)[2] );
result = result[order(result$padj), ]
return(result)
}
I my data I have 705 miRNA, the problem is after doing the DE analysis I didn't ant DE miRNA (lmfit
model didn't converge)
I tried to write my own t-test analysis but I got similar results (when using the Counts, none was DE, but when I used RPKM about 2 or 3 miRNA where DE between the diffrent groups). I used the following code for the t-test and multiple test correction,
get.DiffrentlyExpressedMiRNA<-function(counts, groups){
require(qvalue);
Pvals<- rep(0,nrow(counts));
FC<-rep(0,nrow(counts));
logFC<- rep(0,nrow(counts));
groups<-factor(groups)
grp1<-which(groups == levels(groups)[1]);
grp2<-which(groups == levels(groups)[2]);
for(i in 1:nrow(counts)){
Pvals[i]<-t.test(counts[i, grp1], counts[i, grp2])$p.value;
FC[i]<- mean(counts[i,grp1]) / mean(counts[i,grp2]);
logFC[i]<-log2(FC[i]);
}
nas<-which(!is.nan(Pvals));
inf<-which(is.infinite(logFC))
logFC[inf]<-0;
results<-data.frame(miRNA_ID = rownames(counts)[nas], P_val= Pvals[nas], FC = FC[nas], log2FC= logFC[nas], q = qvalue(Pvals[nas])$qvalues);
ord<-order(results$log2FC,decreasing=T);
results<-results[ord,];
return(results);
}
How do you generally guys, analyse the miRNA data?
NB: I just have access to the level 3 data on TCGA (which means I just have the counts and RPKM for each miRNA)
Any help is apreciated. Thanks,
Hi Jeremy, I removed the weakly expressed miRNA, I got some improvement but still no diffrently expressed miRNA. I think the problem is in the set that I selected. I will check further :)
another strategy is to normalize the mirna counts against housekeeping genes (are there housekeeping mirnas?) , but i've never tried that. http://pubs.rsc.org/en/content/articlelanding/2013/tx/c3tx50034a/unauth
Thanks, Jeremy I will take a look at that :)