Mercurial > repos > iuc > limma_voom
comparison limma_voom.R @ 4:a61a6e62e91f draft
planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/limma_voom commit 6a458881c0819b75e55e64b3f494679d43bb9ee8
| author | iuc |
|---|---|
| date | Sun, 29 Apr 2018 17:36:42 -0400 |
| parents | 38aab66ae5cb |
| children | d8a55b5f0de0 |
comparison
equal
deleted
inserted
replaced
| 3:38aab66ae5cb | 4:a61a6e62e91f |
|---|---|
| 347 # Extract counts and annotation data | 347 # Extract counts and annotation data |
| 348 print("Extracting counts") | 348 print("Extracting counts") |
| 349 data <- list() | 349 data <- list() |
| 350 data$counts <- counts | 350 data$counts <- counts |
| 351 if (haveAnno) { | 351 if (haveAnno) { |
| 352 data$genes <- geneanno | 352 # order annotation by genes in counts (assumes gene ids are in 1st column of geneanno) |
| 353 annoord <- geneanno[match(row.names(counts), geneanno[,1]), ] | |
| 354 data$genes <- annoord | |
| 353 } else { | 355 } else { |
| 354 data$genes <- data.frame(GeneID=row.names(counts)) | 356 data$genes <- data.frame(GeneID=row.names(counts)) |
| 355 } | 357 } |
| 356 | 358 |
| 357 # If filter crieteria set, filter out genes that do not have a required cpm/counts in a required number of | 359 # If filter crieteria set, filter out genes that do not have a required cpm/counts in a required number of |
| 429 | 431 |
| 430 if (wantTrend) { | 432 if (wantTrend) { |
| 431 # limma-trend approach | 433 # limma-trend approach |
| 432 logCPM <- cpm(data, log=TRUE, prior.count=opt$trend) | 434 logCPM <- cpm(data, log=TRUE, prior.count=opt$trend) |
| 433 fit <- lmFit(logCPM, design) | 435 fit <- lmFit(logCPM, design) |
| 436 fit$genes <- data$genes | |
| 434 fit <- contrasts.fit(fit, contrasts) | 437 fit <- contrasts.fit(fit, contrasts) |
| 435 if (wantRobust) { | 438 if (wantRobust) { |
| 436 fit <- eBayes(fit, trend=TRUE, robust=TRUE) | 439 fit <- eBayes(fit, trend=TRUE, robust=TRUE) |
| 437 } else { | 440 } else { |
| 438 fit <- eBayes(fit, trend=TRUE, robust=FALSE) | 441 fit <- eBayes(fit, trend=TRUE, robust=FALSE) |
| 457 | 460 |
| 458 plotData <- logCPM | 461 plotData <- logCPM |
| 459 | 462 |
| 460 # Save normalised counts (log2cpm) | 463 # Save normalised counts (log2cpm) |
| 461 if (wantNorm) { | 464 if (wantNorm) { |
| 462 write.table(logCPM, file=normOut, row.names=TRUE, sep="\t") | 465 write.table(logCPM, file=normOut, row.names=TRUE, sep="\t", quote=FALSE) |
| 463 linkData <- rbind(linkData, c((paste0(deMethod, "_", "normcounts.tsv")), (paste0(deMethod, "_", "normcounts.tsv")))) | 466 linkData <- rbind(linkData, c((paste0(deMethod, "_", "normcounts.tsv")), (paste0(deMethod, "_", "normcounts.tsv")))) |
| 464 } | 467 } |
| 465 } else { | 468 } else { |
| 466 # limma-voom approach | 469 # limma-voom approach |
| 467 voomOutPdf <- makeOut("voomplot.pdf") | 470 voomOutPdf <- makeOut("voomplot.pdf") |
| 508 } | 511 } |
| 509 | 512 |
| 510 # Save normalised counts (log2cpm) | 513 # Save normalised counts (log2cpm) |
| 511 if (wantNorm) { | 514 if (wantNorm) { |
| 512 norm_counts <- data.frame(vData$genes, vData$E) | 515 norm_counts <- data.frame(vData$genes, vData$E) |
| 513 write.table(norm_counts, file=normOut, row.names=FALSE, sep="\t") | 516 write.table(norm_counts, file=normOut, row.names=FALSE, sep="\t", quote=FALSE) |
| 514 linkData <- rbind(linkData, c((paste0(deMethod, "_", "normcounts.tsv")), (paste0(deMethod, "_", "normcounts.tsv")))) | 517 linkData <- rbind(linkData, c((paste0(deMethod, "_", "normcounts.tsv")), (paste0(deMethod, "_", "normcounts.tsv")))) |
| 515 } | 518 } |
| 516 | 519 |
| 517 # Fit linear model and estimate dispersion with eBayes | 520 # Fit linear model and estimate dispersion with eBayes |
| 518 voomFit <- contrasts.fit(voomFit, contrasts) | 521 voomFit <- contrasts.fit(voomFit, contrasts) |
| 552 downCount[i] <- sumStatus["Down", i] | 555 downCount[i] <- sumStatus["Down", i] |
| 553 flatCount[i] <- sumStatus["NotSig", i] | 556 flatCount[i] <- sumStatus["NotSig", i] |
| 554 | 557 |
| 555 # Write top expressions table | 558 # Write top expressions table |
| 556 top <- topTable(fit, coef=i, number=Inf, sort.by="P") | 559 top <- topTable(fit, coef=i, number=Inf, sort.by="P") |
| 557 if (wantTrend) { | 560 write.table(top, file=topOut[i], row.names=FALSE, sep="\t", quote=FALSE) |
| 558 write.table(top, file=topOut[i], row.names=TRUE, sep="\t") | |
| 559 } else { | |
| 560 write.table(top, file=topOut[i], row.names=FALSE, sep="\t") | |
| 561 } | |
| 562 | 561 |
| 563 linkName <- paste0(deMethod, "_", contrastData[i], ".tsv") | 562 linkName <- paste0(deMethod, "_", contrastData[i], ".tsv") |
| 564 linkAddr <- paste0(deMethod, "_", contrastData[i], ".tsv") | 563 linkAddr <- paste0(deMethod, "_", contrastData[i], ".tsv") |
| 565 linkData <- rbind(linkData, c(linkName, linkAddr)) | 564 linkData <- rbind(linkData, c(linkName, linkAddr)) |
| 566 | 565 |
| 625 | 624 |
| 626 cata("<html>\n") | 625 cata("<html>\n") |
| 627 | 626 |
| 628 cata("<body>\n") | 627 cata("<body>\n") |
| 629 cata("<h3>Limma Analysis Output:</h3>\n") | 628 cata("<h3>Limma Analysis Output:</h3>\n") |
| 630 cata("PDF copies of JPEGS available in 'Plots' section.<br />\n") | 629 cata("Links to PDF copies of plots are in 'Plots' section below />\n") |
| 631 if (wantWeight) { | 630 if (wantWeight) { |
| 632 HtmlImage(imageData$Link[1], imageData$Label[1], width=1000) | 631 HtmlImage(imageData$Link[1], imageData$Label[1], width=1000) |
| 633 } else { | 632 } else { |
| 634 HtmlImage(imageData$Link[1], imageData$Label[1]) | 633 HtmlImage(imageData$Link[1], imageData$Label[1]) |
| 635 } | 634 } |