Mercurial > repos > iarc > mutspec

comparison R/figs_MutSpecStat_Galaxy.r @ 7:eda59b985b1c draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author iarc
date Mon, 13 Mar 2017 08:21:19 -0400
parents
children
comparison
equal deleted inserted replaced
6:46a10309dfe2 7:eda59b985b1c
1 #!/usr/bin/Rscript
2
3 #-----------------------------------#
4 # Author: Maude #
5 # Script: figs_MutSpecStat_Galaxy.r #
6 # Last update: 18/10/16 #
7 #-----------------------------------#
8
9
10 #########################################################################################################################################
11 # Create the figures for the report and the HTML page #
12 #########################################################################################################################################
13
14
15 #-------------------------------------------------------------------------------
16 # Load library for recovering the arguments
17 #-------------------------------------------------------------------------------
18 suppressMessages(suppressWarnings(require("getopt")))
19
20
21
22 #-------------------------------------------------------------------------------
23 # Recover the arguments
24 #-------------------------------------------------------------------------------
25 spec = matrix(c(
26 "folderFigure", "folderFigure", 1, "character",
27 "folderTemp", "folderTemp", 1, "character",
28 "filename", "filename", 1, "character",
29 "help", "h", 0, "logical"
30 ),
31 byrow=TRUE, ncol=4
32 )
33
34 opt = getopt(spec)
35
36 # No argument is pass to the command line
37 if(length(opt) == 1)
38 {
39 cat(paste("Usage:\n figs_MutSpecStat_Galaxy.r --folderFigure <path_to_folder> --folderTemp <path_to_tempFolder> --filename <filename> \n",sep=""))
40 q(status=1)
41 }
42
43 # Help was asked for.
44 if ( !is.null(opt$help) )
45 {
46 # print a friendly message and exit with a non-zero error code
47 cat(paste("Usage:\n figs_MutSpecStat_Galaxy.r --folderFigure <path_to_folder> --filename <filename> \n",sep=""))
48 q(status=1)
49 }
50
51
52
53 #-------------------------------------------------------------------------------
54 # Load library
55 #-------------------------------------------------------------------------------
56 suppressMessages(suppressWarnings(library(ggplot2)))
57 suppressMessages(suppressWarnings(library(gplots)))
58 suppressMessages(suppressWarnings(library(gtable)))
59 suppressMessages(suppressWarnings(library(grid)))
60
61
62
63
64 #-------------------------------------------------------------------------------
65 # OVERALL MUTATION DISTRIBUTION
66 #-------------------------------------------------------------------------------
67 inputDistrMut <- paste0(opt$folderFigure, "/Overall_mutation_distribution/", opt$filename, "/", opt$filename, "-OverallMutationDistribution.txt")
68 outputDistrMut <- paste0(opt$folderFigure, "/Overall_mutation_distribution/", opt$filename, "/", opt$filename, "-OverallMutationDistribution.png")
69
70 # Load the input file
71 distrMut <- read.table(inputDistrMut, header=T)
72
73 # Add the count of each category in the legend
74 distrMut$Legend[[1]] <- paste0(distrMut$Variant_type[[1]], " (", distrMut$Count[[1]], ")")
75 distrMut$Legend[[2]] <- paste0(distrMut$Variant_type[[2]], " (", distrMut$Count[[2]], ")")
76 distrMut$Legend[[3]] <- paste0(distrMut$Variant_type[[3]], " (", distrMut$Count[[3]], ")")
77
78 # Base plot
79 pie <- ggplot(distrMut, aes(x=factor(""), fill=Legend, weight=Count)) + geom_bar(width=1) + coord_polar(theta="y") + scale_x_discrete("", breaks=NULL) + scale_y_continuous("", breaks=NULL) + labs(fill="")
80 # Background of the plot entire white
81 pie <- pie + theme(panel.grid.major = element_line(colour="white"), panel.grid.minor = element_line(colour="white"), panel.background = element_rect(fill="white"))
82 # Legend on right in 3 rows
83 pie <- pie + theme(legend.position="bottom") + guides(fill=guide_legend(nrow=3))
84 # Change the color and the title of the legend
85 pie <- pie + scale_fill_brewer("Variant type", palette="Set1")
86 # Remove all the margins
87 pie <- pie + theme(plot.margin=unit(c(-1, 0, -1.5, 0), "cm"))
88 # Save the pie chart for the HTML page (higher resolution)
89 options(bitmapType='cairo') # Use cairo device as isn't possible to install X11 on the server...
90 png(outputDistrMut, width=700, height=1100, res=300)
91 print(pie)
92 dev.off()
93
94
95
96 #-------------------------------------------------------------------------------
97 # SBS MUTATION DISTRIBUTION
98 #-------------------------------------------------------------------------------
99 inputDistrSBS <- paste0(opt$folderFigure, "/SBS_distribution/", opt$filename, "/", opt$filename, "-SBS_distribution.txt")
100 outputDistrSBS <- paste0(opt$folderFigure, "/SBS_distribution/", opt$filename, "/", opt$filename, "-SBS_distribution.png")
101 outputDistrSBSReport <- paste0(opt$folderTemp, "/", opt$filename, "-SBS_distribution-Report.png")
102
103
104 # Load the input file
105 distrSBS <- read.delim(inputDistrSBS)
106 distrSBS <- data.frame(distrSBS)
107 bar <- ggplot(distrSBS, aes(x=Mutation_Type, y=Percentage, fill=Mutation_Type))
108 bar <- bar + geom_bar(stat="identity", width=0.5)
109 # Theme classic
110 bar <- bar + theme_classic()
111 # Remove the axis legend
112 bar <- bar + xlab("")
113 # Set the color of the bars and Changing the labels in the legend
114 bar <- bar + scale_fill_manual(values=c("blue", "black", "red", "gray", "#00CC33", "pink"),
115 labels=c("C:G>A:T", "C:G>G:C", "C:G>T:A", "T:A>A:T", "T:A>C:G", "T:A>G:C")
116 )
117 # Remove the label in x axis
118 bar <- bar + theme(axis.text.x = element_blank())
119 # Change the name of the y label
120 bar <- bar + ylab("Percent")
121 # Save the plot for the HTML page (higher resolution)
122 png(outputDistrSBS, width=1800, height=1500, res=300)
123 print(bar)
124 dev.off()
125 # Save the plot for the report
126 bar
127 ggsave(outputDistrSBSReport)
128
129
130
131 #-------------------------------------------------------------------------------
132 # IMPACT ON PROTEIN
133 #-------------------------------------------------------------------------------
134 inputImpactProt <- paste0(opt$folderFigure, "/Impact_protein_sequence/", opt$filename, "/", opt$filename, "-DistributionExoFunc.txt")
135 outputImpactProt <- paste0(opt$folderFigure, "/Impact_protein_sequence/", opt$filename, "/", opt$filename, "-DistributionExoFunc.png")
136 outputImpactProtReport <- paste0(opt$folderTemp, "/", opt$filename, "-DistributionExoFunc-Report.png")
137
138 # Load the input file
139 impactProt <- read.table(inputImpactProt, header=T)
140 # Custom palette: black, orange, dark green, yellow, light blue, dark blue, darkslateblue, red, purple, pink, light green, turquoise, gray
141 cb_palette <- c("#000000", "#E69F00", "#006600", "#660000", "#F0E442", "#56B4E9", "#3300FF", "#483D8B", "#FF0000", "#9900CC", "#FF66CC", "#00CC00", "#66FFFF", "#C0C0C0")
142 pie <- ggplot(impactProt, aes(x=factor(""), fill=AA_Change, weight=Count)) + geom_bar(width=1) + coord_polar(theta="y") + scale_x_discrete("", breaks=NULL)+ scale_y_continuous("", breaks=NULL) + scale_fill_manual(values=cb_palette)
143 # Background of the plot entire white
144 pie <- pie + theme(panel.grid.major = element_line(colour="white"), panel.grid.minor = element_line(colour="white"), panel.background = element_rect(fill="white"))
145 # Legend in two column
146 pie <- pie + guides(fill=guide_legend(ncol=2)) + theme(legend.position="bottom")
147 # Remove the legend title
148 pie <- pie + labs(fill="")
149 # Save the plot for the HTML page (higher resolution)
150 png(outputImpactProt, width=1600, height=1800, res=300)
151 print(pie)
152 dev.off()
153 # Save the plot for the report
154 pie
155 ggsave(outputImpactProtReport)
156
157
158
159
160 #-------------------------------------------------------------------------------
161 # STRAND BIAS
162 #-------------------------------------------------------------------------------
163 inputSB <- paste0(opt$folderFigure, "/Stranded_Analysis/", opt$filename, "/", opt$filename, "-StrandBias.txt")
164 outputSB <- paste0(opt$folderFigure, "/Stranded_Analysis/", opt$filename, "/", opt$filename, "-StrandBias.png")
165 outputSBReport <- paste0(opt$folderTemp, "/", opt$filename, "-StrandBias-Report.png")
166
167 # Load the input file
168 file_sb <- read.table(inputSB, header=T)
169 # Custom palette (blue, red)
170 cb_palette_SB <- c("#0072B2", "#CC0000")
171 # Base plot
172 p_sb <- ggplot(file_sb, aes(x=Alteration, y=Count, fill=Strand)) + theme_classic() + geom_bar(stat="identity", position="dodge") + scale_fill_manual(values=cb_palette_SB) + theme(axis.text.x = element_text(angle=60, hjust=1)) + xlab("") + theme(legend.position="bottom")
173 # Save the plot for the HTML page (higher resolution)
174 png(outputSB, width=1000, height=1200, res=300)
175 print(p_sb)
176 dev.off()
177 # Save the plot for the report
178 p_sb
179 ggsave(outputSBReport)
180
181
182
183
184 #-------------------------------------------------------------------------------
185 # HEATMAP SEQUENCE CONTEXT - GENOMIC STRAND
186 #-------------------------------------------------------------------------------
187 inputHeatmapGenomic <- paste0(opt$folderFigure, "/Trinucleotide_Sequence_Context/", opt$filename, "/", opt$filename, "-HeatmapCount-Genomic.txt")
188 outputHeatmapGenomic <- paste0(opt$folderFigure, "/Trinucleotide_Sequence_Context/", opt$filename, "/", opt$filename, "-HeatmapCount-Genomic.png")
189 outputHeatmapGenomicReport <- paste0(opt$folderTemp, "/", opt$filename, "-HeatmapCount-Genomic-Report.png")
190
191 inputHeatmapGenomicPercent <- paste0(opt$folderFigure, "/Trinucleotide_Sequence_Context/", opt$filename, "/", opt$filename, "-HeatmapPercent-Genomic.txt")
192 outputHeatmapGenomicPercent <- paste0(opt$folderFigure, "/Trinucleotide_Sequence_Context/", opt$filename, "/", opt$filename, "-HeatmapPercent-Genomic.png")
193 outputHeatmapGenomicPercentReport <- paste0(opt$folderTemp, "/", opt$filename, "-HeatmapPercent-Genomic-Report.png")
194
195
196 ## COUNT
197 heatmap_C <- read.table(inputHeatmapGenomic, header=T)
198 # Save the plot for the report
199 png(filename=outputHeatmapGenomicReport, bg="transparent", width=240, height=360)
200 # Heatmap with an absolute scale
201 heatmap.2(as.matrix(heatmap_C),Rowv=F,Colv=F,col=colorpanel(384,low="yellow",high="red"),dendrogram="none",scale="none",trace="none",key=F,labRow=rownames(as.matrix(heatmap_C)),labCol=colnames(as.matrix(heatmap_C)),lmat=rbind(c(5,1,4),c(3,1,2)), lhei=c(0.75,0.75),lwid=c(0.5,1.5,0.5))
202 dev.off()
203 # Save the plot for the HTML page (higher resolution)
204 png(filename=outputHeatmapGenomic, width=1100, height=1600, res=300)
205 heatmap.2(as.matrix(heatmap_C),Rowv=F,Colv=F,col=colorpanel(384,low="yellow",high="red"),dendrogram="none",scale="none",trace="none",key=F,labRow=rownames(as.matrix(heatmap_C)),labCol=colnames(as.matrix(heatmap_C)),lmat=rbind(c(5,1,4),c(3,1,2)), lhei=c(0.75,0.75),lwid=c(0.5,1.5,0.5))
206 dev.off()
207
208 ## PERCENT
209 heatmap_P <- read.table(inputHeatmapGenomicPercent, header=T)
210 # Save the plot for the report
211 png(filename=outputHeatmapGenomicPercentReport,bg="transparent", width=240, height=360)
212 # Heatmap with an absolute scale
213 heatmap.2(as.matrix(heatmap_P),Rowv=F,Colv=F,col=colorpanel(384,low="yellow",high="red"),dendrogram="none",scale="none",trace="none",key=F,labRow=rownames(as.matrix(heatmap_P)),labCol=colnames(as.matrix(heatmap_P)),lmat=rbind(c(5,1,4),c(3,1,2)), lhei=c(0.75,0.75),lwid=c(0.5,1.5,0.5))
214 dev.off()
215 # Save the plot for the HTML page (higher resolution)
216 png(filename=outputHeatmapGenomicPercent, width=1100, height=1600, res=300)
217 heatmap.2(as.matrix(heatmap_P),Rowv=F,Colv=F,col=colorpanel(384,low="yellow",high="red"),dendrogram="none",scale="none",trace="none",key=F,labRow=rownames(as.matrix(heatmap_P)),labCol=colnames(as.matrix(heatmap_P)),lmat=rbind(c(5,1,4),c(3,1,2)), lhei=c(0.75,0.75),lwid=c(0.5,1.5,0.5))
218 dev.off()