Mercurial > repos > devteam > cummerbund

comparison cummeRbund.R @ 0:587c425b4e76 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial commit with version 1.0.0 of the cummeRbund wrapper.
author devteam
date Tue, 23 Dec 2014 15:58:27 -0500
parents
children 78fcfc04fcfe
comparison
equal deleted inserted replaced
-1:000000000000 0:587c425b4e76
1 ## Feature Selection ##
2 get_features <- function(myGenes, f="gene") {
3 if (f == "isoforms")
4 return(isoforms(myGenes))
5 else if (f == "tss")
6 return(TSS(myGenes))
7 else if (f == "cds")
8 return(CDS(myGenes))
9 else
10 return(myGenes)
11 }
12
13 ## Main Function ##
14
15 library(argparse)
16
17 parser <- ArgumentParser(description='Create a plot with cummeRbund')
18
19 parser$add_argument('--type', dest='plotType', default='Density', required=TRUE)
20 parser$add_argument('--height', dest='height', type='integer', default=960, required=TRUE)
21 parser$add_argument('--width', dest='width', type='integer', default=1280, required=TRUE)
22 parser$add_argument('--outfile', dest='filename', default="plot-unknown-0.png", required=TRUE)
23 parser$add_argument('--input', dest='input_database', default="cuffData.db", required=TRUE)
24 parser$add_argument('--smooth', dest='smooth', action="store_true", default=FALSE)
25 parser$add_argument('--gene_selector', dest='gene_selector', action="store_true", default=FALSE)
26 parser$add_argument('--replicates', dest='replicates', action="store_true", default=FALSE)
27 parser$add_argument('--labcol', dest='labcol', action="store_true", default=FALSE)
28 parser$add_argument('--labrow', dest='labrow', action="store_true", default=FALSE)
29 parser$add_argument('--border', dest='border', action="store_true", default=FALSE)
30 parser$add_argument('--summary', dest='summary', action="store_true", default=FALSE)
31 parser$add_argument('--count', dest='count', action="store_true", default=FALSE)
32 parser$add_argument('--error_bars', dest='error_bars', action="store_true", default=FALSE)
33 parser$add_argument('--log10', dest='log10', action="store_true", default=FALSE)
34 parser$add_argument('--features', dest='features', action="store", default="genes")
35 parser$add_argument('--clustering', dest='clustering', action="store", default="both")
36 parser$add_argument('--iter_max', dest='iter_max', action="store")
37 parser$add_argument('--genes', dest='genes', action="append")
38 parser$add_argument('--k', dest='k', action="store")
39 parser$add_argument('--x', dest='x', action="store")
40 parser$add_argument('--y', dest='y', action="store")
41
42 args <- parser$parse_args()
43
44 ## Load cummeRbund library
45 library("cummeRbund")
46
47 ## Initialize cuff object
48 cuff <- readCufflinks(dir = "", dbFile = args$input_database, rebuild = FALSE)
49
50 ## Print out info
51 print(cuff)
52 sink("cuffdb_info.txt")
53 print(cuff)
54 print("SAMPLES:")
55 samples(cuff)
56 print("REPLICATES:")
57 replicates(cuff)
58 print("FEATURES:")
59 print(annotation(genes(cuff)))
60 cat(annotation(genes(cuff))[[1]],sep=",")
61 sink()
62
63 png(filename = args$filename, width = args$width, height = args$height, type=c('cairo-png'))
64 tryCatch({
65 if (args$plotType == 'density') {
66 csDensity(genes(cuff), replicates=args$replicates, logMode=args$log10)
67 }
68 else if (args$plotType == 'boxplot') {
69 csBoxplot(genes(cuff), replicates=args$replicates, logMode=args$log10)
70 }
71 else if (args$plotType == 'mds') {
72 MDSplot(genes(cuff), replicates=args$replicates)
73 }
74 else if (args$plotType == 'pca') {
75 PCAplot(genes(cuff), "PC1", "PC2", replicates=args$replicates)
76 }
77 else if (args$plotType == 'dendrogram') {
78 csDendro(genes(cuff), replicates=args$replicates)
79 }
80 else if (args$plotType == 'scatter') {
81 if (args$gene_selector) {
82 myGenes <- getGenes(cuff, args$genes)
83 csScatter(get_features(myGenes, args$features), args$x, args$y, smooth=args$smooth, logMode=args$log10)
84 }
85 else {
86 csScatter(genes(cuff), args$x, args$y, smooth=args$smooth, logMode=args$log10)
87 }
88 }
89 else if (args$plotType == 'volcano') {
90 if (args$gene_selector) {
91 myGenes <- get_features(getGenes(cuff, args$genes), args$features)
92 }
93 else {
94 myGenes <- genes(cuff)
95 }
96 csVolcano(myGenes, args$x, args$y)
97 }
98 else if (args$plotType == 'heatmap') {
99 if (args$gene_selector) {
100 myGenes <- getGenes(cuff, args$genes)
101 }
102 else {
103 myGenes <- getGenes(cuff,annotation(genes(cuff))[[1]])
104 }
105 csHeatmap(get_features(myGenes, args$features), clustering=args$clustering, labCol=args$labcol, labRow=args$labrow, border=args$border, logMode=args$log10)
106 }
107 else if (args$plotType == 'cluster') {
108 myGenes <- getGenes(cuff, args$genes)
109 csCluster(get_features(myGenes, args$features), k=args$k)
110 }
111 else if (args$plotType == 'dispersion') {
112 dispersionPlot(genes(cuff))
113 }
114 else if (args$plotType == 'fpkmSCV') {
115 fpkmSCVPlot(genes(cuff))
116 }
117 else if (args$plotType == 'scatterMatrix') {
118 csScatterMatrix(genes(cuff))
119 }
120 else if (args$plotType == 'expressionplot') {
121 myGenes <- getGenes(cuff, args$genes)
122 expressionPlot(get_features(myGenes, args$features), drawSummary=args$summary, showErrorbars=args$error_bars, replicates=args$replicates)
123 }
124 else if (args$plotType == 'expressionbarplot') {
125 myGeneId <- args$genes
126 myGenes <- getGenes(cuff, myGeneId)
127 expressionBarplot(get_features(myGenes, args$features), showErrorbars=args$error_bars, replicates=args$replicates)
128 }
129 else if (args$plotType == 'mds') {
130 MDSplot(genes(cuff),replicates=args$replicates)
131 }
132 else if (args$plotType == 'pca') {
133 PCAplot(genes(cuff),"PC1","PC2", replicates=args$replicates)
134 }
135 else if (args$plotType == 'maplot') {
136 MAplot(genes(cuff), args$x, args$y, useCount=args$count)
137 }
138 else if (args$plotType == 'genetrack') {
139 myGene <- getGene(cuff, args$genes)
140 plotTracks(makeGeneRegionTrack(myGene))
141 }
142 },error = function(e) {
143 write(paste("Failed:", e, sep=" "), stderr())
144 q("no", 1, TRUE)
145 })
146 devname = dev.off()