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 :)