# 01 Feb 2011 # G. Robertson # Testing edgeR on miRNA-seq consensus clusters: pairwise # k=6 consensus clusters for 173 LAML samples, q80 miRs # # Robinson MD, McCarthy DJ, Smyth GK. # edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. # Bioinformatics. 2010 Jan 1;26(1):139-40. # source("/projects/remc_bigdata/TCGA/DE/20110421/DE_LAML-STAD.edgeR.20110421.R") ################################################### ### chunk number 1: poisson ################################################### library(edgeR) set.seed(101) setwd("/projects/remc_bigdata/TCGA/DE/20110421") cat("\nread expression matrix\n") mydata <- read.table(file="LAMLSTAD_miRNA_expn.aboveQMedianRecs.20110420.txt", header = TRUE, row.names = 1, sep="\t") dim(mydata) nbr.cols <- ncol(mydata) cat("calculate library.sizes\n") # Need a library.sizes list. mydata.lib.sizes <- colSums(mydata) # lib.sizes[[1]] is the sum of the first data column of 'mydata' cat("read list of cluster assignments across the sample dendrogram\n") # Need to assign consensus cluster numbers to each column of mydata # From ConsensusClusterPlus, have cluster assignments. Use 'header=F' to avoid skipping the first row in the file. # "TCGA.AB.2803",1 # "TCGA.AB.2807",2 # "TCGA.AB.2963",3 consensusClass <- read.csv( file="LAMLSTAD_miRNA_expn.edgeR-sample-types.20110420.txt", header=F) dim(consensusClass) # 173, 2 consensusClassClusterNbrs <- consensusClass[[2]] cat("set the cluster of interest\n") ######################### clusterOfInterest <- 1 ######################### # Use the next line to do the analysis on the cluster of interest vs. all other clusters mydata.groups <- as.integer( consensusClassClusterNbrs==clusterOfInterest ) # Use the next line to do the analysis on the other clusters vs. the cluster-of-interest #mydata.groups <- as.integer( consensusClassClusterNbrs!=clusterOfInterest ) ######################### length(mydata.groups) #pairwiseLabelsFn <- function(k, lis) { as.integer(lis==k) } cat("load the DGEList data structure\n") # This finishes almost instantly, because it's just setting up a data structure d <- DGEList( counts=mydata, group=mydata.groups, lib.size=mydata.lib.sizes ) #dP$samples # group lib.size norm.factors #TCGA.AB.2803 1 993556.1 1 #TCGA.AB.2807 0 991049.3 1 #... cat("calculate normalization factor(s)\n") # This takes only a second to run d <- calcNormFactors(d) #dP$samples # group lib.size norm.factors #TCGA.AB.2803 1 993556.1 0.5241402 #TCGA.AB.2807 0 991049.3 1.5994537 #TCGA.AB.2963 0 995427.2 0.8469909 #TCGA.AB.2826 0 993389.5 0.8838622 #... cat("write MA plot\n") par( mfrow = c(1, 2) ) maPlot( d$counts[, 1], d$counts[, 2], normalize = TRUE, pch = 19, cex = 0.4, ylim = c(-8, 8) ) grid( col = "blue" ) abline( h = log2(d$samples$norm.factors[2]/d$samples$norm.factors[1]), col = "red", lwd = 4 ) eff.libsize <- d$samples$lib.size * d$samples$norm.factors maPlot( d$counts[, 1]/eff.libsize[1], d$counts[, 2]/eff.libsize[2], normalize = FALSE, pch = 19, cex = 0.4, ylim = c(-8, 8) ) grid( col = "blue" ) cat("estimate common dispersion\n") # common dispersion d <- estimateCommonDisp(d) names(d) # [1] "samples" "common.dispersion" "counts" "pseudo.alt" "genes" "all.zeros" "conc" # [8] "common.lib.size" d$samples$lib.size d$common.lib.size # 993484.2 d$common.dispersion # [1] 1.010957 cat("run exactTest(d)\n") # Takes ~5 min to run # Using common dispersion? ######################## de.com <- exactTest(d) ######################## names(de.com) # "table" "comparison" "genes" names(de.com$table) options(digits = 4) topTags(de.com) #Comparison of groups: 1-0 # logConc logFC PValue FDR #hsa-mir-1269 -13.776 6.870 1.132e-140 1.596e-138 #hsa-mir-136 -13.415 6.633 1.095e-132 7.718e-131 #hsa-mir-431 -14.516 6.517 7.859e-125 3.694e-123 #hsa-mir-127 -9.553 6.346 4.816e-123 1.697e-121 #hsa-mir-381 -12.398 6.236 8.383e-118 2.364e-116 #hsa-mir-493 -13.745 6.246 1.782e-116 4.187e-115 #hsa-mir-654 -12.480 6.019 4.821e-109 9.710e-108 #hsa-mir-337 -13.968 6.016 3.126e-107 5.510e-106 #hsa-mir-379 -9.855 5.919 6.223e-106 9.749e-105 #hsa-mir-409 -13.555 5.865 4.194e-102 5.914e-101 cat("add BH-corrected p-val to de.com$table\n") #xx <- cbind(anno[k,-5], counts[k,], id=rownames(de$table), de$table, adjp=p.adjust(de$table$p.value, method="BH") ) de.com$table <- cbind( de.com$table, adjp=p.adjust( de.com$table$p.value, method="BH") ) #topTags(de.com, sort.by = "logFC") tt <- topTags(de.com, sort.by = "p.value") #> d$counts[detags.com, ] #Error in d$counts : $ operator is invalid for atomic vectors sum(de.com$table$p.value < 0.01) # 43 topN <- topTags(de.com, n = 500) #names(topN) # "table" "adjust.method" "comparison" cat("write out the pval-sorted results table\n") write.table( topN$table, "top500.LAML-vs-STAD.20110421.csv",sep="\t", row.names=T, quote=FALSE ) nbrPosLogFC <- sum(topN$table$logFC > 0) # 44 nbrMegLogFC <- sum(topN$table$logFC < 0) # 97 sum( p.adjust(de.com$table$p.value, method = "BH") < 0.05 ) # 45 (38 have p.adj < 0.01) mean( p.adjust(de.com$table$p.value, method = "BH") < 0.05 ) * 100 # 13.91 detagsN <- rownames( topTags(de.com, n = 141)$table ) cat("generate smear plot") plotSmear(d, de.tags = detagsN, main = "FC plot using common dispersion") abline(h = c(-2, 2), col = "dodgerblue") sessionInfo() cat("done") #> top141 #Comparison of groups: 1-0 # logConc logFC PValue FDR #hsa-mir-1269 -13.776 6.86965 1.132e-140 1.596e-138 #hsa-mir-136 -13.415 6.63290 1.095e-132 7.718e-131 #hsa-mir-431 -14.516 6.51740 7.859e-125 3.694e-123 #hsa-mir-127 -9.553 6.34624 4.816e-123 1.697e-121 #hsa-mir-381 -12.398 6.23639 8.383e-118 2.364e-116 #hsa-mir-493 -13.745 6.24586 1.782e-116 4.187e-115 #hsa-mir-654 -12.480 6.01904 4.821e-109 9.710e-108 #hsa-mir-337 -13.968 6.01605 3.126e-107 5.510e-106 #hsa-mir-379 -9.855 5.91860 6.223e-106 9.749e-105 #hsa-mir-409 -13.555 5.86460 4.194e-102 5.914e-101 #hsa-mir-370 -14.163 5.54903 5.005e-90 6.416e-89 #hsa-mir-134 -9.943 5.39349 3.079e-86 3.618e-85 #hsa-mir-224 -12.140 5.24387 1.983e-80 2.151e-79 #hsa-mir-452 -13.142 4.95408 4.358e-70 4.389e-69 #hsa-mir-100 -7.285 4.11249 4.984e-46 4.685e-45 #hsa-mir-125b-1 -9.671 3.37997 4.853e-29 4.276e-28 #hsa-mir-10a -5.660 -5.26918 6.185e-15 5.130e-14 #hsa-mir-542 -11.236 2.49335 1.082e-14 8.478e-14 #hsa-mir-196b -9.147 -4.69266 6.452e-13 4.788e-12 #hsa-mir-203 -12.339 2.22524 1.435e-11 1.012e-10 #hsa-mir-99a -13.514 -3.71031 1.681e-09 1.109e-08 #hsa-mir-424 -10.100 2.01046 1.731e-09 1.109e-08 #hsa-mir-196a-1 -16.274 -2.82832 2.024e-06 1.241e-05 #hsa-let-7c -12.223 -2.69335 2.715e-06 1.595e-05 #hsa-mir-153-2 -15.088 -2.69661 3.539e-06 1.996e-05 #hsa-mir-20b -12.551 -2.23536 6.090e-05 3.303e-04 #hsa-mir-9-2 -10.485 -2.22037 6.630e-05 3.462e-04 #hsa-mir-363 -13.193 -2.19741 8.019e-05 3.946e-04 #hsa-mir-9-1 -10.451 -2.18921 8.115e-05 3.946e-04 #hsa-mir-23a -13.177 -2.00138 2.783e-04 1.308e-03 #hsa-mir-551b -15.364 -2.01673 3.183e-04 1.448e-03 #hsa-mir-151 -11.068 -1.96356 3.385e-04 1.492e-03 #hsa-mir-181a-1 -4.263 1.21951 4.845e-04 2.070e-03 #hsa-let-7b -5.721 -1.73959 1.299e-03 5.386e-03 #hsa-mir-155 -8.011 -1.72099 1.449e-03 5.839e-03 #hsa-mir-181b-1 -6.975 1.10863 1.568e-03 6.142e-03 #hsa-mir-574 -12.021 1.09577 1.823e-03 6.946e-03 #hsa-mir-143 -7.948 1.08760 1.933e-03 7.174e-03 #hsa-mir-584 -11.011 -1.59429 3.017e-03 1.091e-02 #hsa-mir-500 -15.551 -1.59991 3.549e-03 1.251e-02 #hsa-mir-145 -10.654 0.98005 5.285e-03 1.817e-02 #hsa-mir-30c-2 -11.845 0.96976 5.800e-03 1.947e-02 #hsa-mir-27a -10.642 -1.45794 6.375e-03 2.090e-02 #hsa-mir-378 -10.683 -1.30417 1.423e-02 4.559e-02 #hsa-mir-199b -10.177 0.85253 1.506e-02 4.719e-02 #hsa-mir-532 -10.178 -1.27109 1.679e-02 5.147e-02 #hsa-mir-30b -11.702 0.78062 2.568e-02 7.611e-02 #hsa-mir-582 -9.730 -1.18270 2.591e-02 7.611e-02 #hsa-mir-20a -8.842 -1.08870 4.025e-02 1.158e-01 #hsa-mir-183 -10.941 -1.07210 4.362e-02 1.230e-01 #hsa-mir-223 -5.102 -1.05499 4.688e-02 1.296e-01 #hsa-mir-1201 -12.681 -1.03145 5.261e-02 1.426e-01 #hsa-mir-1977 -12.822 -0.99223 6.264e-02 1.659e-01 #hsa-mir-425 -8.077 -0.98650 6.355e-02 1.659e-01 #hsa-mir-181a-2 -8.791 -0.97643 6.642e-02 1.703e-01 #hsa-mir-326 -12.526 -0.96991 6.889e-02 1.735e-01 #hsa-mir-1978 -12.639 -0.94040 7.801e-02 1.909e-01 #hsa-mir-7-1 -11.920 -0.93851 7.859e-02 1.909e-01 #hsa-mir-146a -7.786 0.60173 7.987e-02 1.909e-01 #hsa-mir-19a -11.627 -0.91253 8.755e-02 2.057e-01 #hsa-mir-210 -11.604 -0.90377 9.043e-02 2.090e-01 #hsa-mir-150 -8.121 -0.88813 9.628e-02 2.190e-01 #hsa-mir-15b -8.132 -0.86684 1.050e-01 2.351e-01 #hsa-mir-28 -7.900 -0.83920 1.174e-01 2.585e-01 #hsa-mir-126 -6.115 -0.83455 1.195e-01 2.592e-01 #hsa-mir-34a -12.345 -0.82358 1.254e-01 2.679e-01 #hsa-mir-374a -9.817 0.50947 1.312e-01 2.762e-01 #hsa-mir-182 -9.589 -0.78573 1.447e-01 2.991e-01 #hsa-mir-455 -12.649 0.48825 1.464e-01 2.991e-01 #hsa-mir-1259 -12.253 -0.77826 1.493e-01 2.991e-01 #hsa-mir-30a -10.419 -0.77604 1.506e-01 2.991e-01 #hsa-mir-629 -11.527 -0.76618 1.562e-01 3.059e-01 #hsa-mir-128-2 -11.656 0.44479 1.799e-01 3.475e-01 #hsa-mir-152 -12.987 -0.72418 1.842e-01 3.511e-01 #hsa-mir-24-2 -7.525 -0.70421 1.958e-01 3.680e-01 #hsa-mir-93 -6.224 -0.69552 2.019e-01 3.703e-01 #hsa-mir-1974 -6.389 -0.69515 2.022e-01 3.703e-01 #hsa-mir-505 -11.776 -0.69035 2.066e-01 3.735e-01 #hsa-mir-17 -10.166 -0.66722 2.238e-01 3.994e-01 #hsa-mir-92a-2 -3.660 -0.65061 2.365e-01 4.120e-01 #hsa-mir-19b-2 -8.669 -0.65053 2.367e-01 4.120e-01 #hsa-mir-335 -10.975 0.37752 2.418e-01 4.159e-01 #hsa-mir-10b -8.971 -0.63957 2.460e-01 4.178e-01 #hsa-mir-16-1 -6.806 -0.63365 2.507e-01 4.209e-01 #hsa-mir-148a -4.785 0.36320 2.564e-01 4.254e-01 #hsa-mir-25 -5.170 -0.61262 2.692e-01 4.414e-01 #hsa-mir-21 -2.814 0.34671 2.743e-01 4.445e-01 #hsa-mir-1975 -7.309 -0.59709 2.835e-01 4.522e-01 #hsa-mir-185 -11.905 -0.59463 2.870e-01 4.522e-01 #hsa-let-7f-2 -7.068 0.33390 2.887e-01 4.522e-01 #hsa-mir-128-1 -11.501 0.32816 2.953e-01 4.576e-01 #hsa-let-7g -8.372 -0.57213 3.077e-01 4.716e-01 #hsa-mir-342 -10.846 -0.56459 3.155e-01 4.774e-01 #hsa-mir-101-1 -6.115 0.30869 3.183e-01 4.774e-01 #hsa-mir-628 -12.309 -0.55771 3.229e-01 4.792e-01 #hsa-mir-486 -10.675 -0.53842 3.428e-01 5.035e-01 #hsa-mir-766 -12.630 -0.52215 3.609e-01 5.246e-01 #hsa-mir-140 -8.874 -0.51044 3.736e-01 5.293e-01 #hsa-mir-146b -7.204 0.26396 3.752e-01 5.293e-01 #hsa-mir-92a-1 -10.282 -0.50905 3.754e-01 5.293e-01 #hsa-mir-338 -11.425 -0.50315 3.833e-01 5.351e-01 #hsa-mir-99b -6.140 0.24128 4.061e-01 5.614e-01 #hsa-mir-106b -8.441 -0.46085 4.330e-01 5.927e-01 #hsa-let-7e -10.058 0.20663 4.560e-01 6.183e-01 #hsa-mir-15a -9.776 -0.43351 4.685e-01 6.291e-01 #hsa-mir-29b-1 -11.794 -0.43054 4.734e-01 6.297e-01 #hsa-mir-484 -10.724 -0.42619 4.782e-01 6.301e-01 #hsa-mir-197 -8.396 -0.41698 4.904e-01 6.402e-01 #hsa-mir-1826 -6.609 -0.41019 4.996e-01 6.463e-01 #hsa-mir-451 -8.503 -0.40506 5.068e-01 6.496e-01 #hsa-mir-130a -11.536 -0.38545 5.353e-01 6.800e-01 #hsa-mir-192 -10.667 -0.37906 5.440e-01 6.836e-01 #hsa-mir-30e -4.962 -0.37166 5.543e-01 6.836e-01 #hsa-let-7a-2 -6.336 -0.36656 5.618e-01 6.836e-01 #hsa-mir-125a -11.176 0.13857 5.618e-01 6.836e-01 #hsa-mir-29a -7.004 0.13819 5.624e-01 6.836e-01 #hsa-mir-191 -6.143 -0.33820 6.044e-01 7.283e-01 #hsa-mir-30d -5.613 -0.32578 6.236e-01 7.451e-01 #hsa-mir-625 -9.615 -0.28601 6.873e-01 8.093e-01 #hsa-let-7d -8.095 0.06287 6.901e-01 8.093e-01 #hsa-mir-423 -9.896 -0.28105 6.953e-01 8.093e-01 #hsa-mir-29b-2 -11.073 -0.27807 7.002e-01 8.093e-01 #hsa-mir-26b -7.977 -0.27342 7.078e-01 8.113e-01 #hsa-mir-345 -10.879 -0.25240 7.432e-01 8.451e-01 #hsa-mir-130b -12.063 -0.24205 7.612e-01 8.527e-01 #hsa-mir-144 -10.569 -0.24142 7.620e-01 8.527e-01 #hsa-mir-221 -8.378 -0.23702 7.691e-01 8.539e-01 #hsa-mir-328 -11.233 0.01009 7.843e-01 8.613e-01 #hsa-mir-181c -11.439 -0.22668 7.880e-01 8.613e-01 #hsa-mir-22 -5.282 -0.19875 8.356e-01 9.012e-01 #hsa-mir-148b -10.345 -0.01894 8.373e-01 9.012e-01 #hsa-mir-186 -9.114 -0.02670 8.512e-01 9.093e-01 #hsa-let-7i -8.241 -0.05088 8.956e-01 9.495e-01 #hsa-mir-1307 -8.657 -0.14921 9.242e-01 9.625e-01 #hsa-mir-29c -9.967 -0.06899 9.286e-01 9.625e-01 #hsa-mir-26a-2 -7.769 -0.14329 9.350e-01 9.625e-01 #hsa-mir-361 -9.627 -0.07648 9.426e-01 9.625e-01 #hsa-let-7a-1 -6.743 -0.08409 9.567e-01 9.625e-01 #hsa-let-7a-3 -6.738 -0.08421 9.569e-01 9.625e-01 #hsa-mir-142 -3.377 -0.12886 9.612e-01 9.625e-01 #hsa-mir-222 -9.098 -0.08716 9.625e-01 9.625e-01