Mercurial > repos > yhoogstrate > edger_with_design_matrix
changeset 6:149a52c74f39 draft
Uploaded
| author | yhoogstrate |
|---|---|
| date | Thu, 09 Jan 2014 04:30:24 -0500 |
| parents | c0cd78c1c10d |
| children | 61e42740b13a |
| files | edgeR_DGE.xml |
| diffstat | 1 files changed, 11 insertions(+), 11 deletions(-) [+] |
line wrap: on
line diff
--- a/edgeR_DGE.xml Thu Jan 09 02:55:12 2014 -0500 +++ b/edgeR_DGE.xml Thu Jan 09 04:30:24 2014 -0500 @@ -46,7 +46,7 @@ library(edgeR) ## Fetch commandline arguments -args <- commandArgs(trailingOnly = TRUE) +args <- commandArgs(trailingOnly = TRUE) designmatrix = args[1] contrast = args[2] @@ -67,7 +67,7 @@ library(edgeR) -raw_data <- read.delim(designmatrix,header=T,stringsAsFactors=T) +raw_data <- read.delim(designmatrix,header=T,stringsAsFactors=T) ## Obtain read-counts read_counts = read.delim(as.character(raw_data[1,1]),header=F,stringsAsFactors=F,row.names=1) @@ -101,10 +101,10 @@ colnames(read_counts) = raw_data[,2] dge = DGEList(counts=read_counts,genes=rownames(read_counts)) -design_tmp <- raw_data[3:length(raw_data)] -rownames(design_tmp) <- colnames(dge) +design_tmp <- raw_data[3:length(raw_data)] +rownames(design_tmp) <- colnames(dge) formula = paste(c("~0",colnames(design_tmp)),collapse = " + ") -design <- model.matrix(as.formula(formula),design_tmp) +design <- model.matrix(as.formula(formula),design_tmp) prefixes = colnames(design_tmp)[attr(design,"assign")] avoid = nchar(prefixes) == nchar(colnames(design)) @@ -144,10 +144,10 @@ fit = glmFit(dge,design) print(paste("Performing likelihood ratio test: ",contrast,sep="")) -cont <- c(contrast) -cont <- makeContrasts(contrasts=cont, levels=design) +cont <- c(contrast) +cont <- makeContrasts(contrasts=cont, levels=design) -lrt <- glmLRT(fit, contrast=cont[,1]) +lrt <- glmLRT(fit, contrast=cont[,1]) print(paste("Exporting to file: ",output_1,sep="")) write.table(file=output_1,topTags(lrt,n=nrow(read_counts))$table,sep="\t",row.names=T) write.table(file=output_2,cpm(dge,normalized.lib.sizes=TRUE),sep="\t") @@ -160,11 +160,11 @@ print("Creating MA plots...") - etable <- topTags(lrt, n=nrow(dge))$table - etable <- etable[order(etable$FDR), ] + etable <- topTags(lrt, n=nrow(dge))$table + etable <- etable[order(etable$FDR), ] pdf(output_7) with(etable, plot(logCPM, logFC, pch=20, main="edgeR: Fold change vs abundance")) - with(subset(etable, FDR<0.05), points(logCPM, logFC, pch=20, col="red")) + with(subset(etable, FDR<0.05), points(logCPM, logFC, pch=20, col="red")) abline(h=c(-1,1), col="blue") dev.off() }