|
1
|
1 #! /usr/bin/Rscript
|
|
|
2
|
|
5
|
3 RetreiveAndDestroy=function(opt,root,stem,regexp,SearchNames,Out,isValue,NbFound,StockNames) {
|
|
|
4 Bool=lapply(paste(root,SearchNames,stem,sep=""),grepl,opt)
|
|
|
5 names(Bool)=StockNames
|
|
|
6 Pos=lapply(Bool,which)
|
|
|
7 names(Pos)=StockNames
|
|
|
8 disable=c()
|
|
|
9 for(i in StockNames) {
|
|
|
10 nbmatch=length(Pos[[i]])
|
|
|
11 if(nbmatch>0) {
|
|
|
12 NbFound[[i]]=NbFound[[i]]+nbmatch
|
|
|
13 disable=c(disable,-1*Pos[[i]])
|
|
|
14 if(is.null(Out[[i]])) {
|
|
|
15 if(isValue[[i]]!=0) {
|
|
|
16 if(regexp=="next") {
|
|
|
17 Out[[i]]=opt[Pos[[i]]+1]
|
|
|
18 disable=c(disable,-1*(Pos[[i]]+1))
|
|
|
19 }
|
|
|
20 else {
|
|
|
21 Out[[i]]=sub(regexp,"\\1",opt[Pos[[i]]])
|
|
|
22 }
|
|
|
23 }
|
|
|
24 else {
|
|
|
25 Out[[i]]=TRUE
|
|
1
|
26 }
|
|
5
|
27 }
|
|
|
28 else {
|
|
|
29 if(isValue[[i]]!=0) {
|
|
|
30 if(regexp=="next") {
|
|
|
31 Out[[i]]=c(Out[[i]],opt[Pos[[i]]+1])
|
|
|
32 disable=c(disable,-1*(Pos[[i]]+1))
|
|
|
33 }
|
|
|
34 else {
|
|
|
35 Out[[i]]=c(Out[[i]],sub(regexp,"\\1",opt[Pos[[i]]]))
|
|
|
36 }
|
|
1
|
37 }
|
|
5
|
38 else {
|
|
|
39 Out[[i]]=c(Out[[i]],TRUE)
|
|
|
40 }
|
|
|
41 }
|
|
|
42 }
|
|
|
43 }
|
|
|
44 if(length(disable)>0) {
|
|
|
45 opt=opt[disable]
|
|
|
46 }
|
|
|
47 Out[["ARGUMENT"]]=list()
|
|
|
48 Out[["ARGUMENT"]][["opt"]]=opt
|
|
|
49 Out[["ARGUMENT"]][["NbFound"]]=NbFound
|
|
|
50 return(Out)
|
|
1
|
51 }
|
|
|
52
|
|
5
|
53 getopt=function (spec=NULL,opt=commandArgs()) {
|
|
|
54 FindArgs=which(opt=="--args")
|
|
|
55 if(length(FindArgs)!=1) {
|
|
|
56 stop(length(FindArgs)," --args found where 1 expected.",call.=F)
|
|
|
57 }
|
|
|
58 ExecName=sub("--file=","",opt[FindArgs-1])
|
|
|
59
|
|
|
60 if(FindArgs<length(opt)) {
|
|
|
61 opt=opt[(FindArgs+1):length(opt)]
|
|
|
62 }
|
|
|
63 else {
|
|
|
64 opt=""
|
|
|
65 }
|
|
|
66
|
|
|
67 min.columns=5
|
|
|
68 colNames=c("LongName","ShortName","Flag","Mod","Default")
|
|
|
69 max.columns=6
|
|
|
70 DimSpec=dim(spec)
|
|
|
71 if(DimSpec[2]>min.columns) {
|
|
|
72 colNames=c(colNames,"Description")
|
|
|
73 }
|
|
|
74
|
|
|
75 if(is.null(spec) | !is.matrix(spec) | (DimSpec[2]<min.columns | DimSpec[2]>max.columns)) {
|
|
|
76 stop('argument "spec" is required and must be a matrix with 4|5 columns.',call.=F)
|
|
|
77 }
|
|
|
78 colnames(spec)=colNames
|
|
|
79
|
|
|
80 spec=as.data.frame(spec,stringsAsFactors=F)
|
|
|
81 #spec validation
|
|
|
82 if(length(unique(c(spec$ShortName,"ARGUMENT","args")))!=DimSpec[1]+2 | length(unique(spec$LongName))!=DimSpec[1]) {
|
|
|
83 stop('Long|Short names for flags must be unique (Long name : "ARGUMENT" and "args" forbidden).',
|
|
|
84 "\n","List of duplicated :",
|
|
|
85 "\n","Short: ",paste(spec$ShortName[duplicated(c(spec$ShortName,"ARGUMENT","args"))],collapse=" "),
|
|
|
86 "\n","Long: ",paste(spec$ShortName[duplicated(spec$LongName)],collapse=" "),call.=F)
|
|
|
87 }
|
|
|
88 if(length(which(nchar(spec$ShortName)>1))!=0) {
|
|
|
89 stop('Short names flags can\'t be longer than 1 character.')
|
|
|
90 }
|
|
|
91
|
|
|
92 #initialize
|
|
|
93 Out=list()
|
|
|
94 Short2Long=list()
|
|
|
95 NbFound=list()
|
|
|
96 isValue=list()
|
|
|
97 for(i in 1:DimSpec[1]) {
|
|
|
98 Short2Long[[spec$ShortName[i]]]=spec$LongName[i]
|
|
|
99 NbFound[[spec$LongName[i]]]=0
|
|
|
100 isValue[[spec$LongName[i]]]=spec$Flag[i]
|
|
|
101 }
|
|
|
102
|
|
|
103 #Map, retreive and suppress ARGUMENTs and arguments
|
|
|
104 #Value ARGUMENT --example=value
|
|
|
105 Out=RetreiveAndDestroy(opt,"^--","=.+$",".+=(.+)$",spec$LongName,Out,isValue,NbFound,spec$LongName)
|
|
|
106 opt=Out[["ARGUMENT"]][["opt"]]
|
|
|
107 NbFound=Out[["ARGUMENT"]][["NbFound"]]
|
|
|
108 Out[["ARGUMENT"]]=NULL
|
|
|
109 #boolean ARGUMENT --example
|
|
|
110 Out=RetreiveAndDestroy(opt,"^--","$","$",spec$LongName,Out,isValue,NbFound,spec$LongName)
|
|
|
111 opt=Out[["ARGUMENT"]][["opt"]]
|
|
|
112 NbFound=Out[["ARGUMENT"]][["NbFound"]]
|
|
|
113 Out[["ARGUMENT"]]=NULL
|
|
|
114 #short name ARGUMENT -t value OR boolean -t
|
|
|
115 Out=RetreiveAndDestroy(opt,"^-","$","next",spec$ShortName,Out,isValue,NbFound,spec$LongName)
|
|
|
116 opt=Out[["ARGUMENT"]][["opt"]]
|
|
|
117 NbFound=Out[["ARGUMENT"]][["NbFound"]]
|
|
|
118 Out[["ARGUMENT"]]=NULL
|
|
|
119 #Warn about non mapped ARGUMENTs
|
|
|
120 if(length(opt)>0) {
|
|
|
121 PosUnkArg=which(grepl("^-",opt))
|
|
|
122 if(length(PosUnkArg)) {
|
|
|
123 message("Error, argument unreconized :","\n",paste(opt[PosUnkArg],collapse="\n"),"\n\n")
|
|
|
124 }
|
|
|
125 if(length(PosUnkArg)>0) {
|
|
|
126 opt=opt[PosUnkArg*-1]
|
|
|
127 }
|
|
|
128 }
|
|
|
129 #Arguments
|
|
|
130 Out[["ARGUMENT"]]=opt
|
|
|
131
|
|
|
132 #Validation of ARGUMENTs
|
|
|
133 for(i in 1:DimSpec[1]) {
|
|
|
134 if(spec$Flag[i]=="0") {#verify boolean arguments
|
|
|
135 NbValue=length(Out[[spec$LongName[i]]])
|
|
|
136 if(NbValue>1) {
|
|
|
137 message("Warning : ",spec$LongName[i]," found ",NbValue," times")
|
|
|
138 }
|
|
|
139 }
|
|
|
140 if(length(Out[[spec$LongName[i]]])==0) {
|
|
|
141 Out[[spec$LongName[i]]]=spec$Default[i]
|
|
|
142 }
|
|
|
143 library("methods")
|
|
|
144 Out[[spec$LongName[i]]]=as(Out[[spec$LongName[i]]],spec$Mod[i])
|
|
|
145 }
|
|
|
146
|
|
|
147 return(Out)
|
|
1
|
148 }
|
|
|
149
|
|
5
|
150 ## Converts the flag numbers into binary
|
|
|
151 integer.base.b=function(x, b=2) {
|
|
|
152 xi=as.integer(x)
|
|
1
|
153 if(any(is.na(xi) | ((x-xi)!=0)))
|
|
|
154 print(list(ERROR="x not integer", x=x))
|
|
5
|
155 N=length(x)
|
|
|
156 xMax=max(c(x,1))
|
|
|
157 ndigits=11
|
|
|
158 Base.b=array(NA, dim=c(N, ndigits))
|
|
|
159 for(i in 1:ndigits) {#i=1
|
|
|
160 Base.b[, ndigits-i+1]=(x %% b)
|
|
|
161 x=(x %/% b)
|
|
1
|
162 }
|
|
|
163 Base.b
|
|
|
164 }
|
|
|
165
|
|
5
|
166 ## Checks if color used is an acceptable color
|
|
|
167 is.acceptable.color=function(character) {
|
|
1
|
168 tmp=try(col2rgb(character),TRUE)
|
|
|
169 return(class(tmp)!="try-error")
|
|
|
170 }
|
|
|
171
|
|
5
|
172 extractFromBam=function(file,which,what) {
|
|
1
|
173 return(scanBam(file, param=ScanBamParam(which=which,what=what))[[1]][[1]])
|
|
|
174 }
|
|
|
175
|
|
5
|
176 ## Returns a list from cigar expression
|
|
|
177 interpreteCIGAR=function(cigar) {
|
|
|
178 cigar_un=strsplit(unique(cigar),split="")
|
|
|
179 n_cigar_un=length(cigar_un)
|
|
|
180 taille_cigar=list()
|
|
|
181 analise_cigar=function(cigar_) {
|
|
|
182 cigar_sortie=list()
|
|
|
183 acc=""
|
|
1
|
184 for(j in 1:length(cigar_)) {
|
|
|
185 if(sum(cigar_[j]==as.character(0:9))==1) {
|
|
|
186 acc=paste(acc,cigar_[j],sep="")
|
|
|
187 }
|
|
|
188 else {
|
|
|
189 cigar_sortie[[length(cigar_sortie)+1]]=as.integer(acc)
|
|
|
190 cigar_sortie[[length(cigar_sortie)+1]]=cigar_[j]
|
|
|
191 acc=""
|
|
|
192 }
|
|
|
193 }
|
|
|
194 return(cigar_sortie)
|
|
|
195 }
|
|
5
|
196 cigar_interprete=lapply(cigar_un,analise_cigar)
|
|
|
197 names(cigar_interprete)=unique(cigar)
|
|
1
|
198
|
|
|
199 return(cigar_interprete)
|
|
|
200 }
|
|
|
201
|
|
|
202 # prend un CIGAR splités et retourne la taille occupé par le read sur la séquence génomique (introns compris)
|
|
5
|
203 calcule_longueur_cigar=function(cigar_) {
|
|
|
204 lon=0
|
|
1
|
205 N=length(cigar_)
|
|
|
206 for(j in seq(2,N,2)) {
|
|
|
207 if(cigar_[[j]]!="I") {
|
|
|
208 lon=lon+cigar_[[j-1]]
|
|
|
209 }
|
|
|
210 }
|
|
|
211 return(lon)
|
|
|
212 }
|
|
|
213
|
|
|
214 # prend un CIGAR splités et retourne les positions
|
|
5
|
215 calcule_junction_cigar=function(cigar_) {
|
|
1
|
216 retour=list()
|
|
5
|
217 lon=0
|
|
1
|
218 N=length(cigar_)
|
|
|
219 for(j in seq(2,N,2)) {
|
|
|
220 if(cigar_[[j]]!="I") {
|
|
|
221 lon=lon+cigar_[[j-1]]
|
|
|
222 }
|
|
|
223 if(cigar_[[j]]=="N") {
|
|
|
224 retour[[length(retour)+1]]=c(lon-cigar_[[j-1]]+1,lon)
|
|
|
225 }
|
|
|
226 }
|
|
|
227 return(retour)
|
|
|
228 }
|
|
|
229
|
|
5
|
230 ## Returns a list of numbers of single read with their coordinates
|
|
|
231 compresse_coordonnees=function(debut,fin) {
|
|
1
|
232 if(length(debut)==0) {
|
|
|
233 return(list(numeric(),numeric(),numeric()))
|
|
|
234 }
|
|
|
235 else {
|
|
|
236 tmp=sort(paste(debut,fin,sep="_"))
|
|
|
237 tmp_rle=rle(tmp)
|
|
|
238 poids=tmp_rle$lengths
|
|
|
239 values_split=strsplit(tmp_rle$values,split="_")
|
|
5
|
240 doit=function(j) {
|
|
1
|
241 return(as.integer(values_split[[j]][1]))
|
|
|
242 }
|
|
|
243 debut_uni=sapply(1:length(poids),doit)
|
|
5
|
244 doit=function(j) {
|
|
1
|
245 return(as.integer(values_split[[j]][2]))
|
|
|
246 }
|
|
|
247 fin_uni=sapply(1:length(poids),doit)
|
|
5
|
248 ordre_debut=order(debut_uni)
|
|
1
|
249 return(list(debut_uni[ordre_debut],fin_uni[ordre_debut],poids[ordre_debut]))
|
|
|
250 }
|
|
|
251 }
|
|
|
252
|
|
5
|
253 RDataFileName=function(file) {
|
|
1
|
254 return(paste(file,".RData",sep=""))
|
|
|
255 }
|
|
|
256
|
|
5
|
257 ## Function converts and extracts the infos from bamfile
|
|
|
258 readBam_=function(file_,insert_max_=2000,stranded_=TRUE,ncore_=1,libraryType_=c("standard","inverse"),fileNameRData_=NA,normalized_=NULL,chrName_=NULL,from_=1,to_=NULL) {
|
|
1
|
259 suppressPackageStartupMessages(require("Rsamtools"))
|
|
|
260 suppressPackageStartupMessages(require("GenomicRanges"))
|
|
5
|
261 ## Declaration of variables
|
|
|
262 flagstat=numeric(11)
|
|
1
|
263 names(flagstat)=c("total","duplicates","mapped","paired","read1","read2","properly paired","itself and mate mapped","singletons","mate mapped on a different chr","QC-failed")
|
|
|
264 genome_info=scanBamHeader(file_)[[1]]$targets
|
|
5
|
265 noms_chromosomes=names(genome_info)
|
|
|
266 longueur_chromosomes=as.integer(genome_info)
|
|
|
267 nombre_chromosomes=length(noms_chromosomes)
|
|
|
268 brin_F=list()
|
|
|
269 brin_R=list()
|
|
|
270 brin_F_junction=list()
|
|
|
271 brin_R_junction=list()
|
|
|
272 pas=c(1,2,6,7)
|
|
|
273 i_zone=0
|
|
1
|
274 if(is.null(chrName_)) {
|
|
5
|
275 chrName__=noms_chromosomes
|
|
1
|
276 }
|
|
|
277 else {
|
|
5
|
278 chrName__=chrName_
|
|
1
|
279 }
|
|
5
|
280 ## Fragments identification
|
|
1
|
281 for(i in (1:nombre_chromosomes)) {
|
|
5
|
282 i_zone=i_zone +1
|
|
|
283 nom_chromo=noms_chromosomes[i]
|
|
|
284 lon_chromo=longueur_chromosomes[i]
|
|
1
|
285
|
|
|
286 if(!(nom_chromo %in% chrName__)) {
|
|
5
|
287 brin_F[[i]]=list(numeric(),numeric(),numeric())
|
|
|
288 brin_R[[i]]=list(numeric(),numeric(),numeric())
|
|
|
289 brin_F_junction[[i]]=list(numeric(),numeric(),numeric())
|
|
|
290 brin_R_junction[[i]]=list(numeric(),numeric(),numeric())
|
|
1
|
291 }
|
|
|
292 else {
|
|
|
293 if(is.null(to_)) {
|
|
|
294 to_i=lon_chromo
|
|
|
295 }
|
|
|
296 else {
|
|
|
297 to_i=to_[min(i_zone,length(to_))]
|
|
|
298 }
|
|
|
299 from_i=from_[min(i_zone,length(from_))]
|
|
|
300
|
|
5
|
301 commande=paste("RangesList(`",nom_chromo,"`=IRanges(",from_i,",",to_i,"))",sep="")
|
|
1
|
302 expr=try(parse(text=commande),TRUE)
|
|
5
|
303 ## Function used from GenomicRanges package
|
|
|
304 which=eval(expr)
|
|
|
305 what=c("flag","mpos","cigar","mrnm","isize")
|
|
|
306 param=ScanBamParam(what=what, which=which)
|
|
|
307 ## Case of no reads on the chromosome
|
|
4
|
308 start=extractFromBam(file=file_,which=which,what="pos")
|
|
5
|
309 if(length(start)==0 ) {
|
|
|
310 brin_F[[i]]=list(numeric(),numeric(),numeric())
|
|
|
311 brin_R[[i]]=list(numeric(),numeric(),numeric())
|
|
|
312 brin_F_junction[[i]]=list(numeric(),numeric(),numeric())
|
|
|
313 brin_R_junction[[i]]=list(numeric(),numeric(),numeric())
|
|
1
|
314 }
|
|
|
315 else {
|
|
4
|
316 strand=extractFromBam(file=file_,which=which,what="strand")
|
|
|
317 flag=extractFromBam(file=file_,which=which,what="flag")
|
|
|
318 mpos=extractFromBam(file=file_,which=which,what="mpos")
|
|
|
319 cigar=extractFromBam(file=file_,which=which,what="cigar")
|
|
|
320 mrnm=extractFromBam(file=file_,which=which,what="mrnm")
|
|
|
321 isize=extractFromBam(file=file_,which=which,what="isize")
|
|
5
|
322
|
|
1
|
323 first_read=integer.base.b(flag)[,5]==1
|
|
5
|
324 strand[strand=="+" & !first_read ]="*"
|
|
|
325 strand[strand=="-" & !first_read ]="+"
|
|
|
326 strand[strand=="*" & !first_read ]="-"
|
|
1
|
327
|
|
5
|
328 ## CIGAR's interpreter
|
|
|
329 cigar_interprete=interpreteCIGAR(cigar)
|
|
|
330 longueur_cigar=lapply(cigar_interprete,calcule_longueur_cigar)
|
|
|
331 junction_cigar=lapply(cigar_interprete,calcule_junction_cigar)
|
|
4
|
332
|
|
|
333 end=start+sapply(1:length(cigar),function(j) longueur_cigar[[cigar[j]]])
|
|
|
334
|
|
5
|
335 ## Case of pairend reads
|
|
|
336 is_on_same_chr=mrnm==nom_chromo
|
|
|
337 is_on_same_chr[is.na(is_on_same_chr)]=FALSE
|
|
|
338 is_paired=is_on_same_chr & abs(isize) <=insert_max_
|
|
|
339 is_paired[first_read & strand=="+" & (isize<0 | isize>insert_max_)]=FALSE
|
|
|
340 is_paired[!first_read & strand=="+" & (isize>0 | isize < -insert_max_)]=FALSE
|
|
|
341 is_paired[first_read & strand=="-" & (isize>0 | isize < -insert_max_)]=FALSE
|
|
|
342 is_paired[!first_read & strand=="-" & (isize<0 | isize>insert_max_)]=FALSE
|
|
|
343 is_paired[is.na(is_paired)]=FALSE
|
|
|
344
|
|
|
345 debut_fragment_paired_plus=mpos[!first_read & strand=="+" & is_paired]
|
|
|
346 fin_fragment_paired_plus=end[!first_read & strand=="+" & is_paired]
|
|
|
347 debut_fragment_paired_moins=mpos[first_read & strand=="-" & is_paired]
|
|
|
348 fin_fragment_paired_moins=end[first_read & strand=="-" & is_paired]
|
|
|
349
|
|
|
350 ## Case of single reads
|
|
|
351 debut_fragment_singleton_plus=start[!is_paired & strand=="+"]
|
|
|
352 fin_fragment_singleton_plus=end[!is_paired & strand=="+"]
|
|
|
353 debut_fragment_singleton_moins=start[!is_paired & strand=="-"]
|
|
|
354 fin_fragment_singleton_moins=end[!is_paired & strand=="-"]
|
|
1
|
355
|
|
5
|
356 ## Fragments
|
|
|
357 debut_frag_plus=c(debut_fragment_paired_plus,debut_fragment_singleton_plus)
|
|
|
358 fin_frag_plus=c(fin_fragment_paired_plus,fin_fragment_singleton_plus)
|
|
1
|
359 debut_frag_moins=c(debut_fragment_paired_moins,debut_fragment_singleton_moins)
|
|
5
|
360 fin_frag_moins=c(fin_fragment_paired_moins,fin_fragment_singleton_moins)
|
|
|
361 brin_F[[i]]=compresse_coordonnees(debut_frag_plus,fin_frag_plus)
|
|
|
362 brin_R[[i]]=compresse_coordonnees(debut_frag_moins,fin_frag_moins)
|
|
1
|
363
|
|
5
|
364 ## Junction read
|
|
|
365 debut_junction=numeric()
|
|
|
366 fin_junction=numeric()
|
|
|
367 brin_junction=numeric()
|
|
|
368 i_junction=0
|
|
1
|
369 for(j in 1:length(cigar)) {
|
|
5
|
370 junctions_=junction_cigar[[cigar[j]]]
|
|
|
371 if(length(junctions_)) {
|
|
1
|
372 for(k in 1:length(junctions_)) {
|
|
5
|
373 i_junction=i_junction + 1
|
|
|
374 debut_junction[i_junction]=start[j] + junctions_[[k]][1] - 1
|
|
|
375 fin_junction[i_junction]=start[j] + junctions_[[k]][2] - 1
|
|
|
376 brin_junction[i_junction]=as.character(strand[j])
|
|
1
|
377 }
|
|
|
378 }
|
|
|
379 }
|
|
|
380 if(i_junction==0) {
|
|
5
|
381 brin_F_junction[[i]]=list(numeric(),numeric(),numeric())
|
|
|
382 brin_R_junction[[i]]=list(numeric(),numeric(),numeric())
|
|
1
|
383 }
|
|
|
384 else {
|
|
5
|
385 brin_F_junction[[i]]=compresse_coordonnees(debut_junction[brin_junction=="+"],fin_junction[brin_junction=="+"])
|
|
|
386 brin_R_junction[[i]]=compresse_coordonnees(debut_junction[brin_junction=="-"],fin_junction[brin_junction=="-"])
|
|
1
|
387 }
|
|
5
|
388
|
|
|
389 ## Flagstat interpreter
|
|
|
390 flag_bits=integer.base.b(flag)#remplie les données stat pour un flag donné
|
|
1
|
391
|
|
5
|
392 ## flagstat
|
|
|
393 ## total
|
|
|
394 flagstat[1]=flagstat[1] + sum(flag_bits[,2]==0)
|
|
|
395 ## duplicates
|
|
|
396 flagstat[2]=flagstat[2] + sum((flag_bits[,1]==1)&(flag_bits[,2]==0))
|
|
|
397 ## mapped
|
|
|
398 flagstat[3]=flagstat[3] + sum((flag_bits[,9]==0)&(flag_bits[,2]==0))
|
|
|
399 ## paired
|
|
|
400 flagstat[4]=flagstat[4] + sum((flag_bits[,11]==1)&(flag_bits[,2]==0))
|
|
|
401 ## read1
|
|
|
402 flagstat[5]=flagstat[5] + sum((flag_bits[,5]==1)&(flag_bits[,2]==0))
|
|
|
403 ## read2
|
|
|
404 flagstat[6]=flagstat[6] + sum((flag_bits[,4]==1)&(flag_bits[,2]==0))
|
|
|
405 ## iself and mate mapped
|
|
|
406 flagstat[8]=flagstat[8] + sum((flag_bits[,11]==1)&(flag_bits[,9]==0)&(flag_bits[,8]==0)&(flag_bits[,2]==0))
|
|
|
407 ## singletons
|
|
|
408 flagstat[9]=flagstat[9] + sum((flag_bits[,8]==1)&(flag_bits[,2]==0))
|
|
|
409 ## QC-failed
|
|
|
410 flagstat[11]=flagstat[11] + sum(flag_bits[,2]==1)
|
|
|
411 ## flagstat
|
|
|
412 ## mate on a different chr
|
|
|
413 flagstat[10]=flagstat[10] + sum((!is_on_same_chr)&(flag_bits[,11]==1)&(flag_bits[,9]==0)&(flag_bits[,8]==0)&(flag_bits[,2]==0))
|
|
|
414 ## flagstat
|
|
|
415 ## properly paired
|
|
|
416 flagstat[7]=flagstat[7] + sum(is_paired)
|
|
1
|
417 }
|
|
|
418 }
|
|
|
419 }
|
|
5
|
420
|
|
|
421 ## Data storing
|
|
|
422 names(brin_F)=noms_chromosomes
|
|
|
423 names(brin_R)=noms_chromosomes
|
|
|
424 names(brin_F_junction)=noms_chromosomes
|
|
|
425 names(brin_R_junction)=noms_chromosomes
|
|
1
|
426
|
|
5
|
427 bamHandler=list()
|
|
1
|
428 if(libraryType_[1]=="inverse") {
|
|
5
|
429 bamHandler[[1]]=brin_R
|
|
|
430 bamHandler[[2]]=brin_F
|
|
1
|
431 }
|
|
|
432 else {
|
|
5
|
433 bamHandler[[1]]=brin_F
|
|
|
434 bamHandler[[2]]=brin_R
|
|
1
|
435 }
|
|
5
|
436 bamHandler[[3]]=longueur_chromosomes
|
|
|
437 bamHandler[[4]]=flagstat
|
|
|
438 bamHandler[[5]]=stranded_
|
|
1
|
439
|
|
|
440 if(libraryType_[1]=="inverse") {
|
|
5
|
441 bamHandler[[6]]=brin_R_junction
|
|
|
442 bamHandler[[7]]=brin_F_junction
|
|
1
|
443 }
|
|
|
444 else {
|
|
5
|
445 bamHandler[[6]]=brin_F_junction
|
|
|
446 bamHandler[[7]]=brin_R_junction
|
|
1
|
447 }
|
|
5
|
448 bamHandler[[8]]=FALSE
|
|
|
449 if(!is.null(normalized_)) {
|
|
1
|
450 for( i in pas) {
|
|
5
|
451 for(j in 1:nombre_chromosomes) {
|
|
|
452 bamHandler[[i]][[j]][[4]]=normalized_*bamHandler[[i]][[j]][[3]]
|
|
|
453 }
|
|
1
|
454 }
|
|
5
|
455 bamHandler[[8]]=TRUE
|
|
|
456 }
|
|
1
|
457
|
|
5
|
458 names(bamHandler)=c("F","R","chrLength","flagstat","stranded","junctions_F","junctions_R","norm")
|
|
1
|
459 if((is.null(chrName_))&(from_==1)&(is.null(to_))) {
|
|
|
460 if(is.null(fileNameRData_)|is.na(fileNameRData_)) {
|
|
5
|
461 save(bamHandler,file=RDataFileName(file_))
|
|
1
|
462 }
|
|
|
463 else {
|
|
5
|
464 save(bamHandler,file=fileNameRData_)
|
|
1
|
465 }
|
|
|
466 }
|
|
|
467 return(bamHandler)
|
|
|
468 }
|
|
|
469
|
|
5
|
470 ## Returns the sum of two bamHandler objects
|
|
|
471 addBam=function(bamHandler1,bamHandler2) {
|
|
4
|
472 if(class(bamHandler1)=="try-error"|class(bamHandler2)=="try-error") {
|
|
|
473 if(class(bamHandler1)=="try-error") {
|
|
|
474 return(bamHandler2)
|
|
|
475 }
|
|
|
476 else {
|
|
|
477 return(bamHandler1)
|
|
|
478 }
|
|
|
479 }
|
|
|
480 else {
|
|
5
|
481 brin_F=list()
|
|
|
482 junctions_brin_F=list()
|
|
|
483 brin_R=list()
|
|
|
484 junctions_brin_R=list()
|
|
|
485 L=length(bamHandler1$F)
|
|
|
486 bamHandler=list()
|
|
|
487 if(bamHandler1$norm !=bamHandler2$norm ) {
|
|
4
|
488 warning(expr="Two different bam files(normalized and non normalized)!!",immediate.=TRUE)
|
|
|
489 }
|
|
|
490 for(i in 1:L) {
|
|
5
|
491 brin_F[[i]]=list(c(bamHandler1$F[[i]][[1]],bamHandler2$F[[i]][[1]]),c(bamHandler1$F[[i]][[2]],bamHandler2$F[[i]][[2]]),c(bamHandler1$F[[i]][[3]],bamHandler2$F[[i]][[3]]))
|
|
|
492 junctions_brin_F[[i]]=list(c(bamHandler1$junctions_F[[i]][[1]],bamHandler2$junctions_F[[i]][[1]]),c(bamHandler1$junctions_F[[i]][[2]],bamHandler2$junctions_F[[i]][[2]]),c(bamHandler1$junctions_F[[i]][[3]],bamHandler2$junctions_F[[i]][[3]]))
|
|
|
493 brin_R[[i]]=list(c(bamHandler1$R[[i]][[1]],bamHandler2$R[[i]][[1]]),c(bamHandler1$R[[i]][[2]],bamHandler2$R[[i]][[2]]),c(bamHandler1$R[[i]][[3]],bamHandler2$R[[i]][[3]]))
|
|
|
494 junctions_brin_R[[i]]=list(c(bamHandler1$junctions_R[[i]][[1]],bamHandler2$junctions_R[[i]][[1]]),c(bamHandler1$junctions_R[[i]][[2]],bamHandler2$junctions_R[[i]][[2]]),c(bamHandler1$junctions_R[[i]][[3]],bamHandler2$junctions_R[[i]][[3]]))
|
|
|
495 if(bamHandler1$norm & bamHandler2$norm) {
|
|
|
496 brin_F[[i]][[4]]=c(bamHandler1$F[[i]][[4]],bamHandler2$F[[i]][[4]])
|
|
|
497 junctions_brin_F[[i]][[4]]=c(bamHandler1$junctions_F[[i]][[4]],bamHandler2$junctions_F[[i]][[4]])
|
|
|
498 brin_R[[i]][[4]]=c(bamHandler1$R[[i]][[4]],bamHandler2$R[[i]][[4]])
|
|
|
499 junctions_brin_R[[i]][[4]]=c(bamHandler1$junctions_R[[i]][[4]],bamHandler2$junctions_R[[i]][[4]])
|
|
4
|
500 }
|
|
|
501 }
|
|
5
|
502 names(brin_F)=names(bamHandler1$F)
|
|
|
503 names(brin_R)=names(bamHandler2$R)
|
|
|
504 names(junctions_brin_F)=names(bamHandler1$junctions_F)
|
|
|
505 names(junctions_brin_R)=names(bamHandler2$junctions_R)
|
|
|
506 bamHandler[[1]]=brin_F
|
|
|
507 bamHandler[[2]]=brin_R
|
|
|
508 bamHandler[[3]]=bamHandler1[[3]]
|
|
|
509 bamHandler[[4]]=bamHandler1[[4]] + bamHandler2[[4]]
|
|
|
510 bamHandler[[5]]=bamHandler1[[5]] & bamHandler2[[5]]
|
|
|
511 bamHandler[[6]]=junctions_brin_F
|
|
|
512 bamHandler[[7]]=junctions_brin_R
|
|
|
513 bamHandler[[8]]=bamHandler1$norm & bamHandler2$norm
|
|
|
514 names(bamHandler)=c("F","R","chrLength","flagstat","stranded","junctions_F","junctions_R","norm")
|
|
4
|
515 return(bamHandler)
|
|
|
516 }
|
|
|
517 }
|
|
|
518
|
|
5
|
519 ## Extracts the signal from bamHandler objects
|
|
|
520 extractSignal=function(bamHandlerList,chrName,from=1, to=NULL,normalized_=FALSE) {
|
|
1
|
521 forward=list()
|
|
|
522 reverse=list()
|
|
|
523 chr=which(names(bamHandlerList[[1]]$F)==chrName)
|
|
5
|
524 if(is.null(to)) {
|
|
1
|
525 to=bamHandlerList[[1]]$chrLength[chr]
|
|
|
526 }
|
|
5
|
527 for(i in 1:length(bamHandlerList)) {
|
|
|
528 if(normalized_) {
|
|
|
529 if(bamHandlerList[[i]]$norm) {
|
|
|
530 end=4
|
|
|
531 }
|
|
|
532 else {
|
|
|
533 end=3
|
|
|
534 }
|
|
|
535 }
|
|
|
536 else {
|
|
|
537 end=3
|
|
|
538 }
|
|
4
|
539 forward_=numeric(to-from+1)
|
|
5
|
540 which_read=which((bamHandlerList[[i]]$F[[chrName]][[2]]>=from)&(bamHandlerList[[i]]$F[[chrName]][[1]]<=to))
|
|
4
|
541 n_reads=length(which_read)
|
|
|
542 if(n_reads>0) {
|
|
1
|
543 for(k in which_read) {
|
|
5
|
544 debut_read=max(1,bamHandlerList[[i]]$F[[chrName]][[1]][k]-from+1)
|
|
|
545 fin_read=min(bamHandlerList[[i]]$F[[chrName]][[2]][k]-from+1,to-from+1)
|
|
1
|
546 forward_[debut_read:fin_read]=forward_[debut_read:fin_read]+bamHandlerList[[i]]$F[[chrName]][[end]][k]
|
|
|
547 }
|
|
|
548 }
|
|
5
|
549 which_junctions=which((bamHandlerList[[i]]$junctions_F[[chrName]][[2]]>=from)&(bamHandlerList[[i]]$junctions_F[[chrName]][[1]]<=to))
|
|
4
|
550 n_junctions=length(which_junctions)
|
|
|
551 if(n_junctions>0) {
|
|
1
|
552 for(k in which_junctions) {
|
|
5
|
553 debut_junction=max(1,bamHandlerList[[i]]$junctions_F[[chrName]][[1]][k]-from+1)
|
|
|
554 fin_junction=min(bamHandlerList[[i]]$junctions_F[[chrName]][[2]][k]-from+1,to-from+1)
|
|
1
|
555 forward_[debut_junction:fin_junction]=forward_[debut_junction:fin_junction]-bamHandlerList[[i]]$junctions_F[[chrName]][[end]][k]
|
|
|
556 }
|
|
|
557 }
|
|
|
558 reverse_=numeric(to-from+1)
|
|
5
|
559 which_read=which((bamHandlerList[[i]]$R[[chrName]][[2]]>=from)&(bamHandlerList[[i]]$R[[chrName]][[1]]<=to))
|
|
|
560 n_reads=length(which_read)
|
|
4
|
561 if(n_reads>0) {
|
|
|
562 for(k in which_read) {
|
|
5
|
563 debut_read=max(1,bamHandlerList[[i]]$R[[chrName]][[1]][k]-from+1)
|
|
|
564 fin_read=min(bamHandlerList[[i]]$R[[chrName]][[2]][k]-from+1,to-from+1)
|
|
4
|
565 reverse_[debut_read:fin_read]=reverse_[debut_read:fin_read]+bamHandlerList[[i]]$R[[chrName]][[end]][k]
|
|
1
|
566 }
|
|
|
567 }
|
|
5
|
568 which_junctions=which((bamHandlerList[[i]]$junctions_R[[chrName]][[2]]>=from)&(bamHandlerList[[i]]$junctions_R[[chrName]][[1]]<=to))
|
|
4
|
569 n_junctions=length(which_junctions)
|
|
|
570 if(n_junctions>0) {
|
|
|
571 for(k in which_junctions) {
|
|
5
|
572 debut_junction=max(1,bamHandlerList[[i]]$junctions_R[[chrName]][[1]][k]-from+1)
|
|
|
573 fin_junction=min(bamHandlerList[[i]]$junctions_R[[chrName]][[2]][k]-from+1,to-from+1)
|
|
1
|
574 reverse_[debut_junction:fin_junction]=reverse_[debut_junction:fin_junction]-bamHandlerList[[i]]$junctions_R[[chrName]][[end]][k]
|
|
|
575 }
|
|
|
576 }
|
|
|
577 forward_[forward_<0]=0
|
|
|
578 reverse_[reverse_<0]=0
|
|
|
579 if(bamHandlerList[[i]]$stranded) {
|
|
|
580 forward[[i]]=forward_
|
|
|
581 reverse[[i]]=reverse_
|
|
|
582 }
|
|
|
583 else {
|
|
|
584 forward[[i]]=forward_+reverse_
|
|
|
585 reverse[[i]]=numeric(to-from+1)
|
|
|
586 }
|
|
4
|
587 }
|
|
5
|
588 chr_=list()
|
|
|
589 chr_$F=forward
|
|
|
590 chr_$R=reverse
|
|
1
|
591 return(chr_)
|
|
|
592 }
|
|
|
593
|
|
5
|
594 ## Intern function for readGff function
|
|
1
|
595 my.read.lines2=function(fname) {
|
|
5
|
596 s=file.info( fname )$size
|
|
|
597 buf=readChar( fname, s, useBytes=T)
|
|
1
|
598 strsplit( buf,"\n",fixed=T,useBytes=T)[[1]]
|
|
|
599 }
|
|
|
600
|
|
5
|
601 ## Extracts the annotation infos from Gff file
|
|
|
602 readGff=function(file_in, from=1, to=Inf, chr=NULL, infoName=c("ID","Name","Parent","gene","Alias","orf_classification","Ontology_term","Note","GO")) {
|
|
1
|
603 tmp=try(my.read.lines2(file_in))
|
|
|
604 if(!is.null(chr)) {
|
|
|
605 tmp1=grep(chr, tmp, value=TRUE,useBytes=T)
|
|
|
606 }
|
|
|
607 else {
|
|
|
608 tmp1=tmp
|
|
|
609 }
|
|
|
610 N=length(tmp1)
|
|
5
|
611 Chr=array()
|
|
|
612 Start=array()
|
|
|
613 Stop=array()
|
|
1
|
614 Strand=array()
|
|
|
615 Type=array()
|
|
|
616 info=list()
|
|
5
|
617 for(i in 1:length(infoName)) info[[i]]=array()
|
|
1
|
618 names(info)=infoName
|
|
5
|
619 j=1
|
|
1
|
620 for (i in 1:N) {
|
|
|
621 if(substr(tmp1[i],1,1)!="#") {
|
|
|
622 line_split=unlist(strsplit(tmp1[i],"\t",fixed=T,useBytes=T))
|
|
|
623 if((as.integer(line_split[4])<=to) & (as.integer(line_split[5])>=from)) {
|
|
5
|
624 Chr[j]=line_split[1]
|
|
|
625 Start[j]=as.integer(line_split[4])
|
|
|
626 Stop[j]=as.integer(line_split[5])
|
|
1
|
627 Strand[j]=line_split[7]
|
|
5
|
628 Type[j]=line_split[3]
|
|
1
|
629 ninth=unlist(strsplit(line_split[9],";",fixed=T,useBytes=T))
|
|
|
630 element_ninth_empty=rep(TRUE,length(infoName))
|
|
|
631 for(element_ninth in ninth) {
|
|
5
|
632 element_ninth_split=unlist(strsplit(element_ninth,"=",fixed=T,useBytes=T))
|
|
1
|
633 if(length(element_ninth_split)==2) {
|
|
|
634 if(element_ninth_split[1] %in% infoName) {
|
|
|
635 info[[element_ninth_split[1]]][j]=element_ninth_split[2]
|
|
|
636 element_ninth_empty[infoName==element_ninth_split[1]]=FALSE
|
|
|
637 }
|
|
|
638 }
|
|
|
639 }
|
|
|
640 for(infoName_ in infoName[element_ninth_empty]) {
|
|
|
641 info[[infoName_]][j]="."
|
|
|
642 }
|
|
5
|
643 j=j+1
|
|
1
|
644 }
|
|
|
645 }
|
|
|
646 }
|
|
5
|
647 retour=data.frame(Chr,Type,Start,Stop,Strand,info,stringsAsFactors=FALSE)
|
|
1
|
648 return(retour)
|
|
|
649 }
|
|
|
650
|
|
5
|
651 ## Returns the classic visualisation
|
|
|
652 plotRNAseq=function(forward,reverse,debut_vue=1,fin_vue=length(forward),chr=NULL,annot=NULL,style=NULL,top=NULL,bottom=NULL,x="",y="",titre="",repeated=FALSE,name_flags="",decal=0,ataxises=NULL,classic_plus_color="navyblue",classic_minus_color="mediumvioletred",stranded=TRUE) {
|
|
1
|
653 if(repeated) {
|
|
|
654 forward_=numeric(fin_vue)
|
|
|
655 forward_[debut_vue:fin_vue]=forward
|
|
|
656 reverse_=numeric(fin_vue)
|
|
|
657 reverse_[debut_vue:fin_vue]=reverse
|
|
|
658 }
|
|
|
659 else {
|
|
|
660 forward_=forward
|
|
|
661 reverse_=reverse
|
|
|
662 }
|
|
|
663 if(is.null(top)) {
|
|
|
664 top=max(forward_[debut_vue:fin_vue])
|
|
|
665 }
|
|
|
666 if(is.null(bottom)) {
|
|
|
667 bottom=max(reverse_[debut_vue:fin_vue])
|
|
|
668 }
|
|
|
669 if(is.null(ataxises)) {
|
|
|
670 plot(c(debut_vue,fin_vue)+decal,c(-bottom,top),ylim=c(-bottom,top),xlab=x,ylab=y,main=titre,col="white",xaxs="i",cex.main=2,yaxt="n",cex.lab=1.8)
|
|
|
671 }
|
|
|
672 else {
|
|
|
673 plot(c(debut_vue,fin_vue)+decal,c(-bottom,top),ylim=c(-bottom,top),xlab=x,ylab=y,main=titre,col="white",xaxs="i",xaxt="n",cex.main=2,yaxt="n",cex.lab=1.8)
|
|
|
674 ataxisesLabels=as.character(ataxises)
|
|
|
675 ataxisesLabels[((1:length(ataxises))%%2)==0]=""
|
|
5
|
676 ataxisesLabels[1]=""
|
|
|
677 ataxisesLabels[length(ataxises)]=""
|
|
1
|
678 lim=c(-bottom,top)
|
|
|
679 ataxises_y=pretty(lim,n=4)
|
|
|
680 ataxisesLabels_y=as.character(abs(ataxises_y))
|
|
|
681 axis(1,at=ataxises,labels=FALSE,cex.axis=2)
|
|
|
682 axis(1,at=ataxises,labels=ataxisesLabels,cex.axis=2,line=0.4,lwd=0)
|
|
|
683 axis(2, at=ataxises_y,labels=FALSE,cex.axis=2)
|
|
|
684 axis(2, at=ataxises_y,labels=ataxisesLabels_y,cex.axis=2,line=-0.4,lwd=0)
|
|
|
685 }
|
|
|
686 polygon(c(debut_vue,debut_vue:fin_vue,fin_vue)+decal,c(0,forward_[debut_vue:fin_vue],0),col=classic_plus_color,border=NA)
|
|
5
|
687 if(stranded) {
|
|
1
|
688 text(fin_vue+(fin_vue-debut_vue)*0.01,top/2,"+",xpd=NA,cex=3)
|
|
|
689 text(fin_vue+(fin_vue-debut_vue)*0.01,-bottom/2,"-",xpd=NA,cex=3)
|
|
|
690 text(fin_vue+(fin_vue-debut_vue)*0.025,(top-bottom)/2,"Strand",xpd=NA,cex=2,srt=-90)
|
|
|
691 polygon(c(debut_vue,debut_vue:fin_vue,fin_vue)+decal,c(0,-reverse_[debut_vue:fin_vue],0),col=classic_minus_color,border=NA)
|
|
|
692 abline(h=0,lwd=2)
|
|
|
693 abline(h=0,col="white")
|
|
|
694 }
|
|
|
695 if(name_flags!="") {
|
|
5
|
696 flags=try(get(name_flags),TRUE)
|
|
|
697 if(class(flags)!="try-error") {
|
|
|
698 f_=flags[(flags$Chr==chr)&(flags$Stop>=debut_vue)&(flags$Start<=fin_vue),]
|
|
|
699 N=dim(f_)[1]
|
|
|
700 points(f_$Start,rep(0,N),col=2,pch=19,cex=2)
|
|
|
701 }
|
|
|
702 }
|
|
|
703 if(decal<=0) {
|
|
1
|
704 lines(c(0,0),c(top,-bottom),col=1)
|
|
5
|
705 }
|
|
1
|
706 }
|
|
|
707
|
|
5
|
708 ## Intern function for heatmap visualisation
|
|
|
709 paletteFromColors=function(colMin="blue",colMax="red",n=300,method=c("hsv","rgb")) {
|
|
|
710 colMinRGB=col2rgb(colMin)[,1]/255
|
|
|
711 colMaxRGB=col2rgb(colMax)[,1]/255
|
|
|
712 seqList=list()
|
|
1
|
713 if(method[1]=="rgb") {
|
|
|
714 for(i in 1:3) {
|
|
|
715 seqList[[i]]=seq(colMinRGB[i],colMaxRGB[i],length.out=n)
|
|
|
716 }
|
|
|
717 return(rgb(seqList[[1]],seqList[[2]],seqList[[3]]))
|
|
|
718 }
|
|
|
719 else {
|
|
5
|
720 colMinHSV=rgb2hsv(colMinRGB)
|
|
|
721 colMaxHSV=rgb2hsv(colMaxRGB)
|
|
1
|
722 for(i in 1:3) {
|
|
|
723 seqList[[i]]=seq(colMinHSV[i],colMaxHSV[i],length.out=n)
|
|
|
724 }
|
|
|
725 return(hsv(seqList[[1]],s=seqList[[2]],v=seqList[[3]]))
|
|
|
726 }
|
|
|
727 }
|
|
|
728
|
|
5
|
729 ## Retuns the heatmap visualisation
|
|
|
730 myHeatMap=function(data,debut_vue,fin_vue,ataxises=NULL,lim=NULL,heatmap_max_color="#000055",heatmap_min_color="#FFFFAA",heatmap_palette_method="hsv",textOnLeftSide="") {
|
|
1
|
731 palette=paletteFromColors(heatmap_min_color,heatmap_max_color,method=heatmap_palette_method)
|
|
|
732 if(is.null(lim)) {
|
|
5
|
733 image(debut_vue:fin_vue,1:dim(data)[1],t(data), col=palette,xlab="",ylab="",xaxt="n",yaxt="n")
|
|
1
|
734 }
|
|
|
735 else {
|
|
5
|
736 image(debut_vue:fin_vue,1:dim(data)[1],t(data), col=palette,xlab="",ylab="",xaxt="n",yaxt="n",zlim=lim)
|
|
1
|
737 }
|
|
|
738 box()
|
|
|
739 if(is.null(ataxises)) {
|
|
|
740 axis(1)
|
|
|
741 }
|
|
|
742 else {
|
|
|
743 if(sum(!is.na(ataxises))!=0) {
|
|
|
744 axis(1,at=ataxises,labels=FALSE)
|
|
|
745 axis(3,at=ataxises,labels=FALSE)
|
|
|
746 }
|
|
|
747 }
|
|
|
748 if(sum(!is.na(ataxises))!=0) {
|
|
|
749 axis(2,at=1:dim(data)[1],labels=rownames(data),las=2,cex.axis=1)
|
|
|
750 }
|
|
|
751 text(fin_vue+(fin_vue-debut_vue)*0.015,(dim(data)[1]+1)/2,textOnLeftSide,xpd=NA,cex=1,srt=-90)
|
|
|
752 }
|
|
|
753
|
|
5
|
754 ## Returns the title of visualisation
|
|
1
|
755 plotTitle=function(debut_vue=1,fin_vue=length(listForward[[1]]),chr=NULL,style=NULL) {
|
|
5
|
756 plot(c(0,1),c(0,1),cex=0,ylab="",xlab="",fg="white",axes=FALSE,xaxs="i",yaxs="i")
|
|
|
757 text(0,0.5,paste(chr,":",debut_vue,"-",fin_vue,sep=""),cex=2.2,adj=0)
|
|
|
758 }
|
|
|
759
|
|
|
760 openGraphicalDevice=function(file,widthPixels,heightPixels,fileType=c("png","jpeg","tiff","bmp","pdf"),resolutionDPI=72) {
|
|
|
761 widthInches=widthPixels/72
|
|
|
762 heightInches=heightPixels/72
|
|
|
763 doit=function(x) {
|
|
|
764 switch(x,
|
|
|
765 png(file,widthInches,heightInches,units="in",res=resolutionDPI),
|
|
|
766 jpeg=jpeg(file,widthInches,heightInches,units="in",res=resolutionDPI),
|
|
|
767 tiff=tiff(file,widthInches,heightInches,units="in",res=resolutionDPI),
|
|
|
768 bmp=bmp(file,widthInches,heightInches,units="in",res=resolutionDPI),
|
|
|
769 pdf=pdf(file,widthInches,heightInches))
|
|
|
770 }
|
|
|
771 doit(fileType[1])
|
|
1
|
772 }
|
|
|
773
|
|
5
|
774 ## The main function of visualisation
|
|
|
775 plotVisu=function(file,typeVisu="classic",listForward,listReverse,which=1:length(listForward),stranded=TRUE,
|
|
1
|
776 debut_vue=1,fin_vue=length(listForward[[1]]),chr=NULL,annot=NULL,style=NULL,tops=NULL,bottoms=NULL,marks=NULL,strandMarks=NULL,
|
|
4
|
777 titres="",repeated=FALSE,name_flags="",decal=0,log=TRUE,classic_plus_color="navyblue",classic_minus_color="deeppink3",
|
|
5
|
778 heatmap_max_color="#000055",heatmap_min_color="#FFFFAA",heatmap_palette_method="hsv",heatmap_lane_height=round(10+40/(1+10^((length(which)-12)/9))),
|
|
4
|
779 lines_samples_colors=c(1,3,4,2)[((0:length(listForward))%%4)+1],lines_samples_type_line=((0:length(listForward))%/%4)+1,
|
|
5
|
780 smoothLength=trunc((fin_vue-debut_vue)/1200),annotation_color_by_strand=FALSE,annotation_placed_by_strand=FALSE,display_name=NULL,initialize_label_sizes=NULL,
|
|
|
781 fileType="png",resolutionDPI=72) {
|
|
4
|
782 if(fin_vue-debut_vue+1>1000000) {
|
|
5
|
783 openGraphicalDevice(file,1200,400,fileType=fileType,resolutionDPI=resolutionDPI)
|
|
4
|
784 plot(0:1,0:1,fg="white",col="white",axes=FALSE,frame=FALSE,xlab="",ylab="")
|
|
|
785 text(0.5,0.5,"Window too long !",cex=5)
|
|
|
786 }
|
|
|
787 else {
|
|
5
|
788 n_element_vue=length(which)
|
|
|
789 i_data=which[1]
|
|
1
|
790 forward_matrice=NULL
|
|
|
791 reverse_matrice=NULL
|
|
5
|
792 all_stranded=data[[i_data]]$stranded
|
|
1
|
793 if(!repeated) {
|
|
|
794 for(i in 1:n_element_vue) {
|
|
5
|
795 i_data=which[i]
|
|
|
796 forward_matrice=rbind(forward_matrice,listForward[[i_data]][debut_vue:fin_vue])
|
|
|
797 reverse_matrice=rbind(reverse_matrice,listReverse[[i_data]][debut_vue:fin_vue])
|
|
|
798 all_stranded=all_stranded & data[[i_data]]$stranded
|
|
1
|
799 }
|
|
|
800 }
|
|
5
|
801 else {
|
|
1
|
802 for(i in 1:n_element_vue) {
|
|
5
|
803 i_data=which[i]
|
|
|
804 forward_matrice=rbind(forward_matrice,listForward[[i_data]])
|
|
|
805 reverse_matrice=rbind(reverse_matrice,listReverse[[i_data]])
|
|
|
806 all_stranded=all_stranded & data[[i_data]]$stranded
|
|
1
|
807 }
|
|
|
808 }
|
|
|
809 rownames(forward_matrice)=titres[which]
|
|
|
810 rownames(reverse_matrice)=titres[which]
|
|
5
|
811 for( i in 1:n_element_vue) {
|
|
1
|
812 if(smoothLength>1) {
|
|
|
813 lo=smooth(forward_matrice[i,],L=smoothLength)
|
|
|
814 forward_matrice[i,]=lo
|
|
|
815 los=smooth(reverse_matrice[i,],L=smoothLength)
|
|
|
816 reverse_matrice[i,]=los
|
|
|
817 }
|
|
|
818 }
|
|
|
819 if(is.null(annot)) {
|
|
|
820 heights_ann=NULL
|
|
4
|
821 annot_selec=NULL
|
|
1
|
822 }
|
|
|
823 else {
|
|
4
|
824 annot_selec=annot[(annot$Chr==chr)&(annot$Stop>=debut_vue)&(annot$Start<=fin_vue),]
|
|
5
|
825 heights_ann=sizePlotAnnotation(annot=annot_selec,chr=chr,debut=debut_vue,fin=fin_vue,annotation_placed_by_strand=annotation_placed_by_strand,display_name=display_name,typeVisu=typeVisu,initialize_label_sizes=initialize_label_sizes)*50
|
|
1
|
826 }
|
|
|
827 ataxises=pretty((debut_vue:fin_vue)+decal,n=14)
|
|
5
|
828 if(log) {
|
|
1
|
829 label_scal="log2 tag densities"
|
|
5
|
830 }
|
|
|
831 else {
|
|
1
|
832 label_scal="tag densities"
|
|
|
833 }
|
|
|
834 if(debut_vue==1&decal==0) {
|
|
|
835 ataxises[1]=1
|
|
5
|
836 }## Classic visualisation
|
|
1
|
837 if(typeVisu=="classic") {
|
|
5
|
838 height_panels=c(40,rep(200,n_element_vue),heights_ann)
|
|
|
839 openGraphicalDevice(file,1200,sum(height_panels),fileType=fileType,resolutionDPI=resolutionDPI)
|
|
|
840 prev=par(no.readonly=TRUE)
|
|
1
|
841 n_panels=length(height_panels)
|
|
|
842 layout(matrix(1:n_panels,n_panels,1),heights=height_panels)
|
|
|
843 par(mar=c(0, 5, 0, 4)+0.1)
|
|
|
844 plotTitle(debut_vue=debut_vue,fin_vue=fin_vue,chr=chr)
|
|
|
845 par(mar=c(2.5, 5, 2.5, 4)+0.1)
|
|
|
846 for(element in 1:n_element_vue) {
|
|
5
|
847 i_data=which[element]
|
|
4
|
848 plotRNAseq(listForward[[i_data]],listReverse[[i_data]],debut_vue,fin_vue,chr,annot_selec,style,top=tops[min(element,length(tops))],bottom=bottoms[min(element,length(bottoms))],y=label_scal,titre=titres[min(i_data,length(titres))],name_flags=name_flags,repeated=repeated,decal=decal,ataxises=ataxises,classic_plus_color=classic_plus_color,classic_minus_color=classic_minus_color,stranded=stranded)
|
|
1
|
849 if(!is.null(marks)) {
|
|
|
850 if(is.null(tops)) {
|
|
4
|
851 top=max(listForward[[i_data]])
|
|
1
|
852 }
|
|
|
853 else {
|
|
|
854 top=tops[min(element,length(tops))]
|
|
|
855 }
|
|
4
|
856 if(is.null(bottoms)) {
|
|
|
857 bottom=max(listReverse[[i_data]])
|
|
|
858 }
|
|
|
859 else {
|
|
|
860 bottom=bottoms[min(element,length(bottoms))]
|
|
|
861 }
|
|
|
862 if(is.null(strandMarks)) {
|
|
|
863 segments(marks,top,y1=-bottom,col=2)
|
|
1
|
864 }
|
|
|
865 else {
|
|
4
|
866 if(strandMarks==0|!stranded) {
|
|
|
867 segments(marks,0,y1=top,col=2)
|
|
1
|
868 }
|
|
|
869 else {
|
|
4
|
870 if(strandMarks==1) {
|
|
|
871 segments(marks,0,y1=-bottom,col=2)
|
|
|
872 }
|
|
|
873 else {
|
|
|
874 segments(marks,top,y1=-bottom,col=2)
|
|
|
875 }
|
|
1
|
876 }
|
|
|
877 }
|
|
|
878 }
|
|
5
|
879 }
|
|
|
880 if(!is.null(annot_selec)) {
|
|
|
881 par(mar=c(0,5,1,4)+0.1)
|
|
|
882 plot_annotation(annot_selec,chr,debut=debut_vue,fin=fin_vue,style=style,textSize=1.5,annotation_color_by_strand=annotation_color_by_strand,annotation_placed_by_strand=annotation_placed_by_strand,display_name=display_name,typeVisu=typeVisu,initialize_label_sizes=initialize_label_sizes)
|
|
|
883 }
|
|
1
|
884 }
|
|
5
|
885 else {## Heatmap visualisation
|
|
1
|
886 if(typeVisu=="heatmap") {
|
|
|
887 if(stranded) {
|
|
5
|
888 height_panels=c(40,heatmap_lane_height*n_element_vue+28,heights_ann,heatmap_lane_height*n_element_vue+10,100)
|
|
1
|
889 }
|
|
|
890 else {
|
|
5
|
891 height_panels=c(40,heatmap_lane_height*n_element_vue+28,heights_ann,100)
|
|
1
|
892 }
|
|
5
|
893 openGraphicalDevice(file,1200,sum(height_panels),fileType=fileType,resolutionDPI=resolutionDPI)
|
|
|
894 prev=par(no.readonly=TRUE)
|
|
1
|
895 n_panels=length(height_panels)
|
|
|
896 layout(matrix(1:n_panels,n_panels,1),heights=height_panels)
|
|
|
897 par(mar=c(0, 8, 0, 2)+0.1)
|
|
|
898 plotTitle(debut_vue=debut_vue,fin_vue=fin_vue,chr=chr)
|
|
5
|
899 limIntensity=c(0,max(forward_matrice,reverse_matrice,na.rm=TRUE))
|
|
4
|
900 if(limIntensity[2]==0) {
|
|
|
901 limIntensity[2]=0.01
|
|
|
902 }
|
|
1
|
903 ataxisesLabels=as.character(ataxises)
|
|
|
904 ataxisesLabels[((1:length(ataxises))%%2)==0]=""
|
|
|
905 rownames(forward_matrice)=titres[which]
|
|
|
906 rownames(reverse_matrice)=titres[which]
|
|
|
907 par(mar=c(1, 8, 1.8, 2)+0.1,cex=1.1)
|
|
|
908 if(stranded) {
|
|
|
909 tmp="Plus strand"
|
|
|
910 }
|
|
|
911 else {
|
|
|
912 tmp="Both strands"
|
|
|
913 }
|
|
|
914 myHeatMap(forward_matrice[(dim(forward_matrice)[1]):1,],debut_vue+decal,fin_vue+decal,ataxises,lim=limIntensity,heatmap_max_color=heatmap_max_color,heatmap_min_color=heatmap_min_color,heatmap_palette_method=heatmap_palette_method,textOnLeftSide=tmp)
|
|
5
|
915 axis(3,at=ataxises,labels=FALSE,cex.axis=1.2)
|
|
|
916 axis(3,at=ataxises,labels=ataxisesLabels,cex.axis=1.2,line=-0.4,lwd=0)
|
|
1
|
917 if(!is.null(marks)) {
|
|
|
918 if(is.null(strandMarks)|!all_stranded) {
|
|
|
919 segments(marks,-2,y1=dim(forward_matrice)[1]+2,col=2,lwd=2,xpd=TRUE)
|
|
|
920 }
|
|
|
921 else {
|
|
|
922 if(strandMarks!=1) {
|
|
|
923 segments(marks,-2,y1=dim(forward_matrice)[1]+2,col=2,lwd=2,xpd=TRUE)
|
|
|
924 }
|
|
|
925 }
|
|
|
926 }
|
|
|
927 if(name_flags!="") {
|
|
|
928 flags=try(get(name_flags),TRUE)
|
|
|
929 if(class(flags)!="try-error") {
|
|
|
930 f_=flags[(flags$Chr==chr)&(flags$Stop>=debut_vue)&(flags$Start<=fin_vue),]
|
|
|
931 N=dim(f_)[1]
|
|
|
932 points(f_$Start,rep(0,N),col=2,pch=19,cex=2)
|
|
|
933 }
|
|
|
934 }
|
|
4
|
935 if(!is.null(annot_selec)) {
|
|
1
|
936 par(mar=c(0,8,0,2)+0.1)
|
|
5
|
937 plot_annotation(annot_selec,chr,debut=debut_vue,fin=fin_vue,style=style,textSize=0.9,annotation_color_by_strand=annotation_color_by_strand,annotation_placed_by_strand=annotation_placed_by_strand,display_name=display_name,typeVisu=typeVisu,initialize_label_sizes=initialize_label_sizes)
|
|
1
|
938 }
|
|
|
939 if(stranded) {
|
|
|
940 par(mar=c(1, 8, 0, 2)+0.1)
|
|
|
941 myHeatMap(reverse_matrice[(dim(reverse_matrice)[1]):1,],debut_vue+decal,fin_vue+decal,ataxises,lim=limIntensity,heatmap_max_color=heatmap_max_color,heatmap_min_color=heatmap_min_color,heatmap_palette_method=heatmap_palette_method,textOnLeftSide="Minus strand")
|
|
|
942 if(!is.null(marks)) {
|
|
|
943 if(is.null(strandMarks)|!all_stranded) {
|
|
|
944 segments(marks,-1,y1=dim(forward_matrice)[1]+1,col=2,lwd=2,xpd=TRUE)
|
|
|
945 }
|
|
|
946 else {
|
|
|
947 if(strandMarks!=0) {
|
|
|
948 segments(marks,-1,y1=dim(forward_matrice)[1]+1,col=2,lwd=2,xpd=TRUE)
|
|
|
949 }
|
|
|
950 }
|
|
|
951 }
|
|
5
|
952 axis(1,at=ataxises,labels=ataxisesLabels,cex.axis=1.2)
|
|
1
|
953 }
|
|
|
954 palette=matrix(seq(limIntensity[1],limIntensity[2],length.out=2000),1,2000)
|
|
|
955 rownames(palette)="0"
|
|
|
956 par(mar=c(2,18,3, 12)+0.1)
|
|
|
957 myHeatMap(palette,1,2000,ataxises=NA,heatmap_max_color=heatmap_max_color,heatmap_min_color=heatmap_min_color,heatmap_palette_method=heatmap_palette_method)
|
|
|
958 labelAxisHeatmapLegend=pretty(c(limIntensity[1],limIntensity[2]),n=7)
|
|
|
959 atAxisHeatmapLegend=1+((labelAxisHeatmapLegend-limIntensity[1])/(limIntensity[2]-limIntensity[1]))*1999
|
|
5
|
960 axis(1,at=atAxisHeatmapLegend,labels=FALSE,cex.axis=1.2)
|
|
|
961 axis(1,at=atAxisHeatmapLegend,labels=labelAxisHeatmapLegend,cex.axis=1.2,line=-0.4,lwd=0)
|
|
1
|
962 text(1000,2,label_scal,xpd=NA,font=2,cex=1.4)
|
|
5
|
963 }## Lines visualisation
|
|
|
964 else if(typeVisu=="lines") {
|
|
|
965 legendSize=(floor(n_element_vue/2)+n_element_vue%%2)*40
|
|
1
|
966 height_panels=c(40,legendSize,400,heights_ann)
|
|
|
967 n_panels=length(height_panels)
|
|
5
|
968 openGraphicalDevice(file,1200,sum(height_panels),fileType=fileType,resolutionDPI=resolutionDPI)
|
|
|
969 prev=par(no.readonly=TRUE)
|
|
1
|
970 par(mar=c(0, 5, 0,4)+0.1,cex=1.1)
|
|
|
971 layout(matrix(c(1:n_panels),n_panels,4),heights=height_panels)
|
|
|
972 par(mar=c(0, 5, 0,4)+0.1)
|
|
|
973 plotTitle(debut_vue=debut_vue,fin_vue=fin_vue,chr=chr)
|
|
|
974 lines_legend(n_element_vue=n_element_vue,which,titres,lines_samples_colors=lines_samples_colors,lines_samples_type_line=lines_samples_type_line)
|
|
|
975 par(mar=c(3, 5, 0,4)+0.1)
|
|
|
976 plotlines(forward_matrice,reverse_matrice,which=which,debut_vue, fin_vue , chr, annot, style, tops, bottoms,marks,strandMarks,titres, repeated,name_flags,decal,ataxises=ataxises,n_element_vue=n_element_vue,y=label_scal,lines_samples_colors=lines_samples_colors,lines_samples_type_line=lines_samples_type_line,stranded=stranded)
|
|
4
|
977 if(!is.null(annot_selec)) {
|
|
5
|
978 par(mar=c(0,5,0,4)+0.1)
|
|
|
979 plot_annotation(annot_selec,chr,debut=debut_vue,fin=fin_vue,style=style,textSize=1.5,annotation_color_by_strand=annotation_color_by_strand,annotation_placed_by_strand=annotation_placed_by_strand,display_name=display_name,typeVisu=typeVisu,initialize_label_sizes=initialize_label_sizes)
|
|
1
|
980 }
|
|
|
981
|
|
|
982 }
|
|
|
983 }
|
|
4
|
984 }
|
|
1
|
985 invisible(dev.off())
|
|
|
986 }
|
|
|
987
|
|
5
|
988 ##
|
|
|
989 smooth=function(X,L=10) {
|
|
1
|
990 x_smooth=filter(X,rep(1,L)/L)
|
|
|
991 x_smooth[is.na(x_smooth)]=0
|
|
|
992 return(x_smooth)
|
|
|
993 }
|
|
|
994
|
|
5
|
995 ## Returns lines visualisation
|
|
|
996 plotlines=function(forward_matrice,reverse_matrice,which=1:length(forward), debut_vue=1,fin_vue=length(forward), chr=NULL, annot=NULL, style=NULL, tops=NULL,bottoms=NULL,marks=NULL,strandMarks=NULL, titres="", repeated=FALSE,name_flags="",decal=0,ataxises=NULL,n_element_vue=length(which),y="",lines_samples_colors=c(1,3,4,2)[((0:length(forward))%%4)+1],lines_samples_type_line=((0:length(forward))%%4)+1,stranded=TRUE) {
|
|
|
997 limIntensity=c(-max(bottoms),max(tops))
|
|
|
998 plot(c(debut_vue,fin_vue)+decal,limIntensity,col="white",ylab=y,main="",xaxs="i",xaxt="n",yaxt="n",cex.lab=1.2,xlab="",cex.lab=1.8)
|
|
|
999 lty=lines_samples_type_line
|
|
|
1000 col=lines_samples_colors
|
|
|
1001 ataxises_y=pretty(limIntensity,n=8)
|
|
1
|
1002 for(i in ataxises) abline(v=i,lty=2,col="#808080")
|
|
|
1003 for(i in ataxises_y) abline(h=i,lty=2,col="#808080")
|
|
5
|
1004 for( i in 1:n_element_vue) {
|
|
1
|
1005 lo=forward_matrice[i,]
|
|
5
|
1006 lines((debut_vue:fin_vue)+decal,lo,type="l",lty=lty[i],col=col[i],lwd=2,xaxt="n")
|
|
1
|
1007 if(stranded) {
|
|
|
1008 los=-reverse_matrice[i,]
|
|
5
|
1009 lines((debut_vue:fin_vue)+decal,los,type="l",lty=lty[i],col=col[i],lwd=2,xaxt="n")
|
|
1
|
1010 }
|
|
|
1011 }
|
|
|
1012 ataxisesLabels=as.character(ataxises)
|
|
|
1013 ataxisesLabels[((1:length(ataxises))%%2)==0]=""
|
|
|
1014 ataxisesLabels_y=as.character(abs(ataxises_y))
|
|
|
1015 ataxisesLabels_y[((1:length(ataxises_y))%%2)==0]=""
|
|
5
|
1016 axis(1,at=ataxises,labels=FALSE,cex.axis=2)
|
|
|
1017 axis(1,at=ataxises,labels=ataxisesLabels,cex.axis=2,line=0.6,lwd=0)
|
|
1
|
1018 axis(2, at=ataxises_y,labels=FALSE,cex.axis=2)
|
|
|
1019 axis(2, at=ataxises_y,labels=ataxisesLabels_y,cex.axis=2,line=-0.4,lwd=0)
|
|
|
1020 if(stranded) {
|
|
|
1021 abline(h=0,lwd=6)
|
|
|
1022 abline(h=0,col="white",lwd=4)
|
|
|
1023 text(fin_vue+(fin_vue-debut_vue)*0.01,tops/2,"+",xpd=NA,cex=3)
|
|
|
1024 text(fin_vue+(fin_vue-debut_vue)*0.01,-bottoms/2,"-",xpd=NA,cex=3)
|
|
|
1025 text(fin_vue+(fin_vue-debut_vue)*0.025,(tops-bottoms)/2,"Strand",xpd=NA,cex=2,srt=-90)
|
|
|
1026 }
|
|
4
|
1027 if(!is.null(marks)) {
|
|
|
1028 if(is.null(tops)) {
|
|
|
1029 top=max(listForward[[i_data]])
|
|
|
1030 }
|
|
|
1031 else {
|
|
|
1032 top=max(tops)
|
|
|
1033 }
|
|
|
1034 if(is.null(bottoms)) {
|
|
|
1035 bottom=max(listReverse[[i_data]])
|
|
|
1036 }
|
|
|
1037 else {
|
|
|
1038 bottom=max(bottoms)
|
|
|
1039 }
|
|
|
1040 if(is.null(strandMarks)) {
|
|
|
1041 segments(marks,top,y1=-bottom,col=2)
|
|
|
1042 }
|
|
|
1043 else {
|
|
|
1044 if(strandMarks==0|!stranded) {
|
|
|
1045 segments(marks,0,y1=top,col=2)
|
|
|
1046 }
|
|
|
1047 else {
|
|
|
1048 if(strandMarks==1) {
|
|
|
1049 segments(marks,0,y1=-bottom,col=2)
|
|
|
1050 }
|
|
|
1051 else {
|
|
|
1052 segments(marks,top,y1=-bottom,col=2)
|
|
|
1053 }
|
|
|
1054 }
|
|
|
1055 }
|
|
|
1056 }
|
|
1
|
1057 }
|
|
|
1058
|
|
5
|
1059 ## Returns lines legend
|
|
|
1060 lines_legend=function(n_element_vue,which,titres,lines_samples_colors=c(1,3,4,2)[((0:n_element_vue)%%4)+1],lines_samples_type_line=((0:n_element_vue)%%4)+1) {
|
|
1
|
1061 lty=lines_samples_type_line
|
|
|
1062 col=lines_samples_colors
|
|
5
|
1063 n_y=floor(n_element_vue/2)+n_element_vue%%2
|
|
1
|
1064 plot(c(0,4),c(0,-(n_y+1)),col="white", ylab="",xlab="",main="",fg="white",col.axis="white",yaxs="i")
|
|
|
1065 i_style=0
|
|
|
1066 for(i in 1:n_y) {
|
|
|
1067 i_style=i_style+1
|
|
|
1068 lines(c(1.5,2),-c(i,i),col=col[i_style],lty=lty[i_style],lwd=4)
|
|
|
1069 text(1.48,-i,titres[which[i_style]],cex=2.6,adj=1)
|
|
|
1070 }
|
|
|
1071 if(n_element_vue>1) {
|
|
|
1072 for(i in (n_y+1):n_element_vue) {
|
|
|
1073 i_style=i_style+1
|
|
|
1074 lines(c(3.5,4),-c(i,i)+n_y,col=col[i_style],lty=lty[i_style],lwd=4)
|
|
|
1075 text(3.48,-i+n_y,titres[which[i_style]],cex=2.6,adj=1)
|
|
|
1076 }
|
|
|
1077 }
|
|
|
1078 }
|
|
|
1079
|
|
5
|
1080 ## Returns the shape of plain arrow for the annotation
|
|
|
1081 plain_arrow=function(left,right,y,thickness=1,pickSize=(right-left)*0.1,pickSide=c("right","left","both","none"),col="blue") {
|
|
4
|
1082 middle=(left+right)/2
|
|
1
|
1083 if(pickSide[1]=="right") {
|
|
5
|
1084 pick_point=max(right - pickSize,middle)
|
|
1
|
1085 polygon(c(left,pick_point,right,pick_point,left),c(y-thickness/2,y-thickness/2,y,y+thickness/2,y+thickness/2),col=col)
|
|
|
1086 }
|
|
|
1087 if(pickSide[1]=="left") {
|
|
5
|
1088 pick_point=min(left + pickSize,middle)
|
|
1
|
1089 polygon(c(right,pick_point,left,pick_point,right),c(y-thickness/2,y-thickness/2,y,y+thickness/2,y+thickness/2),col=col)
|
|
|
1090 }
|
|
|
1091 if(pickSide[1]=="none") {
|
|
|
1092 polygon(c(right,left,left,right),c(y-thickness/2,y-thickness/2,y+thickness/2,y+thickness/2),col=col)
|
|
|
1093 }
|
|
|
1094 if(pickSide[1]=="both") {
|
|
5
|
1095 pick_point_1=min(left + pickSize,middle)
|
|
|
1096 pick_point_2=max(right - pickSize,middle)
|
|
1
|
1097 polygon(c(left,pick_point_1,pick_point_2,right,pick_point_2,pick_point_1),c(y,y-thickness/2,y-thickness/2,y,y+thickness/2,y+thickness/2),col=col)
|
|
|
1098 }
|
|
|
1099 }
|
|
|
1100
|
|
5
|
1101 ## Returns the size of the panel of the annotation
|
|
|
1102 sizePlotAnnotation=function(annot,chr,debut,fin,annotation_placed_by_strand=FALSE,display_name=NULL,typeVisu="classic",initialize_label_sizes=NULL) {
|
|
|
1103 left=c()
|
|
|
1104 right=c()
|
|
|
1105 labels=c()
|
|
|
1106 annot_chr=annot[annot$Chr==chr,]
|
|
|
1107 N=dim(annot_chr)[1]
|
|
|
1108 all_names=annot_chr[,c(display_name,"ID")]
|
|
|
1109 unique_ID=unique(annot_chr$ID)
|
|
|
1110 for(j in 1:length(unique_ID)) {
|
|
|
1111 left[j]=min(annot_chr$Start[annot_chr$ID==unique_ID[j]])
|
|
|
1112 right[j]=max(annot_chr$Stop[annot_chr$ID==unique_ID[j]])
|
|
|
1113 all_names=unlist(annot_chr[annot_chr$ID==unique_ID[j],c(display_name,"ID")])
|
|
|
1114 all_names=all_names[all_names!="."]
|
|
|
1115 labels[j]=all_names[1]
|
|
|
1116 }
|
|
|
1117 if(annotation_placed_by_strand) {
|
|
|
1118 y_plot=parking(left,right,debut,fin,annot_chr$Strand=="+",FALSE,labels=labels,typeVisu=typeVisu,initialize_label_sizes=initialize_label_sizes)
|
|
|
1119 }
|
|
|
1120 else {
|
|
|
1121 y_plot=parking(left,right,debut,fin,biggestOnesInside=FALSE,labels=labels,typeVisu=typeVisu,initialize_label_sizes=initialize_label_sizes)
|
|
|
1122 }
|
|
|
1123 return(max(y_plot)-min(y_plot)+1)
|
|
1
|
1124 }
|
|
|
1125
|
|
5
|
1126 ## Function to organise the annotation shapes to display
|
|
|
1127 parking=function(left,right,debut,fin,plus=rep(TRUE,length(left)),labels=c(),biggestOnesInside=TRUE,typeVisu="classic",initialize_label_sizes=NULL) {
|
|
|
1128 if(length(labels)!=0) {
|
|
|
1129 if(!is.null(initialize_label_sizes)) {
|
|
|
1130 initialize_=initialize_label_sizes[[typeVisu]]
|
|
|
1131 for(i in 1:length(left)) {
|
|
|
1132 right[i]=max(right[i],left[i]+initialize_$size[initialize_$labels==labels[i]])+(fin-debut)/100
|
|
|
1133 }
|
|
|
1134 }
|
|
|
1135 }
|
|
1
|
1136 y=rep(0,length(left))
|
|
|
1137 if(sum(plus)>0) {
|
|
|
1138 left_plus=left[plus]
|
|
|
1139 right_plus=right[plus]
|
|
|
1140 y_plus=rep(0,sum(plus))
|
|
|
1141 lengths_plus=right_plus-left_plus
|
|
|
1142 for(i in order(lengths_plus,decreasing=TRUE)) {
|
|
|
1143 otherleft_plus=left_plus
|
|
|
1144 otherleft_plus[i]=NA
|
|
|
1145 otherright_plus=right_plus
|
|
|
1146 otherright_plus[i]=NA
|
|
|
1147 placed=FALSE
|
|
|
1148 y_plus[i]=0.5
|
|
|
1149 while(placed==FALSE) {
|
|
|
1150 placed=sum((right_plus[i]>otherleft_plus[y_plus==y_plus[i]])&(left_plus[i]<otherright_plus[y_plus==y_plus[i]]),na.rm=TRUE)==0
|
|
|
1151 if(placed==FALSE) {
|
|
|
1152 y_plus[i]=y_plus[i]+1
|
|
|
1153 }
|
|
|
1154 }
|
|
|
1155 }
|
|
|
1156 if(biggestOnesInside) {
|
|
|
1157 y[plus]=y_plus
|
|
|
1158 }
|
|
|
1159 else {
|
|
|
1160 y[plus]=sapply(y_plus,function(i) sort(unique(y_plus))[i==sort(unique(y_plus),decreasing=TRUE)])
|
|
|
1161 }
|
|
|
1162 }
|
|
|
1163 if(sum(!plus)>0) {
|
|
|
1164 left_minus=left[!plus]
|
|
|
1165 right_minus=right[!plus]
|
|
|
1166 y_minus=rep(0,sum(!plus))
|
|
|
1167 lengths_minus=right_minus-left_minus
|
|
|
1168 for(i in order(lengths_minus,decreasing=TRUE)) {
|
|
|
1169 otherleft_minus=left_minus
|
|
|
1170 otherleft_minus[i]=NA
|
|
|
1171 otherright_minus=right_minus
|
|
|
1172 otherright_minus[i]=NA
|
|
|
1173 placed=FALSE
|
|
|
1174 y_minus[i]=-0.5
|
|
|
1175 while(placed==FALSE) {
|
|
|
1176 placed=sum((right_minus[i]>otherleft_minus[y_minus==y_minus[i]])&(left_minus[i]<otherright_minus[y_minus==y_minus[i]]),na.rm=TRUE)==0
|
|
|
1177 if(placed==FALSE) {
|
|
|
1178 y_minus[i]=y_minus[i]-1
|
|
|
1179 }
|
|
|
1180 }
|
|
|
1181 }
|
|
|
1182 if(biggestOnesInside) {
|
|
|
1183 y[!plus]=y_minus
|
|
|
1184 }
|
|
|
1185 else {
|
|
|
1186 y[!plus]=sapply(y_minus,function(i) sort(unique(y_minus))[i==sort(unique(y_minus),decreasing=TRUE)])
|
|
|
1187 }
|
|
|
1188 }
|
|
|
1189 return(y)
|
|
|
1190 }
|
|
|
1191
|
|
5
|
1192 ## Function to
|
|
|
1193 plot_annotation=function(annot,chr,debut,fin,style=NULL,textSize=par("cex"),annotation_color_by_strand=FALSE,annotation_placed_by_strand=FALSE,display_name=NULL,typeVisu="classic",initialize_label_sizes=NULL) {
|
|
1
|
1194 left=c()
|
|
|
1195 right=c()
|
|
5
|
1196 labels=c()
|
|
1
|
1197 annot_chr=annot[annot$Chr==chr,]
|
|
5
|
1198 N=dim(annot_chr)[1]
|
|
1
|
1199 if(N>0) {
|
|
5
|
1200 all_names=annot_chr[,c(display_name,"ID")]
|
|
1
|
1201 unique_ID=unique(annot_chr$ID)
|
|
5
|
1202 for(j in 1:length(unique_ID)) {
|
|
1
|
1203 left[j]=min(annot_chr$Start[annot_chr$ID==unique_ID[j]])
|
|
|
1204 right[j]=max(annot_chr$Stop[annot_chr$ID==unique_ID[j]])
|
|
5
|
1205 all_names=unlist(annot_chr[annot_chr$ID==unique_ID[j],c(display_name,"ID")])
|
|
|
1206 all_names=all_names[all_names!="."]
|
|
|
1207 labels[j]=all_names[1]
|
|
1
|
1208 }
|
|
|
1209 if(annotation_placed_by_strand) {
|
|
5
|
1210 y_plot=parking(left,right,debut,fin,annot_chr$Strand=="+",FALSE,labels=labels,typeVisu=typeVisu,initialize_label_sizes=initialize_label_sizes)
|
|
1
|
1211 }
|
|
|
1212 else {
|
|
5
|
1213 y_plot=parking(left,right,debut,fin,biggestOnesInside=FALSE,labels=labels,typeVisu=typeVisu,initialize_label_sizes=initialize_label_sizes)
|
|
1
|
1214 }
|
|
|
1215 plot(c(debut,fin),c(min(y_plot)-0.5,max(y_plot)+0.5),col="white",ylab="",xlab="",fg="white",col.axis="white",xaxs="i",yaxs="i")
|
|
5
|
1216 for(j in 1:length(unique_ID)) {
|
|
1
|
1217 annot_ID=annot_chr[annot_chr$ID==unique_ID[j],]
|
|
5
|
1218 x_text=Inf
|
|
1
|
1219 for(i_annot_ID in 1:dim(annot_ID)[1]) {
|
|
|
1220 if(annot_ID$Strand[i_annot_ID]=="-") {
|
|
|
1221 orientation="left"
|
|
|
1222 }
|
|
|
1223 else {
|
|
|
1224 orientation="right"
|
|
|
1225 }
|
|
5
|
1226 iDraw=FALSE
|
|
1
|
1227 if(annot_ID$Strand[i_annot_ID]==".") {
|
|
|
1228 tmp="+"
|
|
|
1229 }
|
|
|
1230 else {
|
|
|
1231 tmp=annot_ID$Strand[i_annot_ID]
|
|
|
1232 }
|
|
5
|
1233 style_demande=style[style$Type==annot_ID$Type[i_annot_ID]&style$Strand==tmp,]
|
|
|
1234 x_text=min(x_text,annot_ID$Start[i_annot_ID])
|
|
1
|
1235 if(style_demande$shape=="box") {
|
|
|
1236 plain_arrow(annot_ID$Start[i_annot_ID],annot_ID$Stop[i_annot_ID],y_plot[j],thickness=0.5, pickSide=orientation,col=style_demande$col,pickSize=(fin-debut)*0.02)
|
|
|
1237 iDraw=TRUE
|
|
|
1238 }
|
|
5
|
1239 if(style_demande$shape=="line") {
|
|
1
|
1240 segments(annot_ID$Start[i_annot_ID],y_plot[j],annot_ID$Stop[i_annot_ID],y_plot[j],col=style_demande$col)
|
|
|
1241 iDraw=TRUE
|
|
|
1242 }
|
|
5
|
1243 if(style_demande$shape=="rectangle") {
|
|
1
|
1244 plain_arrow(annot_ID$Start[i_annot_ID],annot_ID$Stop[i_annot_ID],y_plot[j],thickness=0.5, pickSide="none",col=style_demande$col)
|
|
|
1245 iDraw=TRUE
|
|
|
1246 }
|
|
5
|
1247 if(style_demande$shape=="arrow") {
|
|
|
1248 arrowHeads=pretty(debut:fin,n=50)
|
|
|
1249 x=c(annot_ID$Start[i_annot_ID],arrowHeads[arrowHeads>annot_ID$Start[i_annot_ID]&arrowHeads<annot_ID$Stop[i_annot_ID]],annot_ID$Stop[i_annot_ID])
|
|
1
|
1250 if(annot_ID$Strand[i_annot_ID]=="-") {
|
|
|
1251 arrows(x[2:length(x)],y_plot[j],x[1:length(x)-1],col=style_demande$col,length=0.08)
|
|
|
1252 }
|
|
|
1253 else {
|
|
|
1254 arrows(x[1:length(x)-1],y_plot[j],x[2:length(x)],col=style_demande$col,length=0.08)
|
|
|
1255 }
|
|
|
1256 iDraw=TRUE
|
|
|
1257 }
|
|
5
|
1258 if(iDraw==FALSE) {
|
|
1
|
1259 segments(annot_ID$Start[i_annot_ID],y_plot[j],annot_ID$Stop[i_annot_ID],y_plot[j],col="black")
|
|
|
1260 }
|
|
|
1261 }
|
|
5
|
1262 text(max(x_text,debut),y_plot[j]+0.4,adj=0,labels=labels[j],cex=textSize)
|
|
1
|
1263 }
|
|
|
1264 if(annotation_placed_by_strand) {
|
|
|
1265 abline(h=0.125,lty=2,col="black")
|
|
|
1266 }
|
|
|
1267 }
|
|
|
1268 }
|
|
|
1269
|
|
5
|
1270 html2rgb=function(codeHTML) {
|
|
1
|
1271 chiffres=c(0:9,"A","B","C","D","E","F")
|
|
|
1272 codeHTMLsplit=unlist(strsplit(codeHTML,split=""))
|
|
5
|
1273 red=(((which(chiffres==codeHTMLsplit[2])-1)*16)+(which(chiffres==codeHTMLsplit[3])-1))/255
|
|
|
1274 green=(((which(chiffres==codeHTMLsplit[4])-1)*16)+(which(chiffres==codeHTMLsplit[5])-1))/255
|
|
|
1275 blue=(((which(chiffres==codeHTMLsplit[6])-1)*16)+(which(chiffres==codeHTMLsplit[7])-1))/255
|
|
1
|
1276 return(c(red,green,blue))
|
|
|
1277 }
|
|
|
1278
|
|
5
|
1279 rgb2hsv=function(vectRGB) {
|
|
1
|
1280 Cmax=max(vectRGB)
|
|
|
1281 Cmin=min(vectRGB)
|
|
|
1282 delta=Cmax-Cmin
|
|
|
1283 if(delta==0) {
|
|
|
1284 H=0
|
|
|
1285 }
|
|
|
1286 else {
|
|
|
1287 if(Cmax==vectRGB[1]) H=(1/6)*(((vectRGB[2]-vectRGB[3])/delta)%%6)
|
|
|
1288 if(Cmax==vectRGB[2]) H=(1/6)*(((vectRGB[3]-vectRGB[1])/delta)+2)
|
|
|
1289 if(Cmax==vectRGB[3]) H=(1/6)*(((vectRGB[1]-vectRGB[2])/delta)+4)
|
|
|
1290 }
|
|
|
1291 if(delta==0) {
|
|
|
1292 S=0
|
|
|
1293 }
|
|
|
1294 else {
|
|
|
1295 S=delta/Cmax
|
|
|
1296 }
|
|
|
1297 V=Cmax
|
|
|
1298
|
|
|
1299 return(c(H,S,V))
|
|
|
1300 }
|
|
|
1301
|
|
|
1302 ### =========================================================================
|
|
5
|
1303 ### Ving's Main version beta 1.1
|
|
1
|
1304 ### -------------------------------------------------------------------------
|
|
|
1305 ###
|
|
|
1306 {
|
|
5
|
1307 if(version$major !="3" & version$minor !="0.2") {
|
|
4
|
1308 write("Ving has been developped on R 3.0.2,
|
|
|
1309 therefore it might misbehave here.",stderr())
|
|
1
|
1310 }
|
|
5
|
1311 arg=commandArgs(TRUE)
|
|
1
|
1312 if(length(arg)==0) {
|
|
5
|
1313 write("Ving for visualisation of RNA seq data:
|
|
1
|
1314 Usage:
|
|
5
|
1315 ./ving.R [options] <input> [<input2> <input3> ...]
|
|
|
1316 Options:
|
|
1
|
1317 -o/--output <string> [default:./ouput.png]
|
|
5
|
1318 -F/--fileType <integer> [default: same as -o] (png,jpeg,bmp,tiff,pdf)
|
|
|
1319 -R/--resolution <integer> [default: 72 ]
|
|
1
|
1320 -v/--typeVisu <string> [default: classic ] (classic,lines,heatmap)
|
|
4
|
1321 -t/--description-data <string> [default: <filename>]
|
|
1
|
1322 -n/--normalization-coefficient <string> [default: none ]
|
|
|
1323 -i/--inverseStrand
|
|
|
1324 -u/--unstranded
|
|
|
1325 -c/--chromosome-name <string> [default: first one ]
|
|
4
|
1326 -S/--start <integer> [default: 1 ]
|
|
|
1327 -E/--end <integer> [default: chr length]
|
|
1
|
1328 -a/--annotation <filename>
|
|
4
|
1329 -r/--typeTranscript <string> [default: <all> ]
|
|
1
|
1330 -C/--annotation-colors <string> [default: fabulous ]
|
|
|
1331 -s/--annotation-shapes <string> [default: box ] (box,line,rectangle,arrow)
|
|
|
1332 -l/--scale-log
|
|
|
1333 -y/--symetric-scale
|
|
4
|
1334 --classic-plus-color <string> [default: navyblue ]
|
|
|
1335 --classic-minus-color <string> [default: deeppink3 ]
|
|
|
1336 --heatmap-max-color <string> [default: 000055 ]
|
|
|
1337 --heatmap-min-color <string> [default: FFFFAA ]
|
|
1
|
1338 --heatmap-palette-method <string> [default: hsv ] (hsv,rgb)
|
|
|
1339 --lines-samples-colors <string> [default: black ]
|
|
|
1340 --lines-samples-type-line <string> [default: plain ]
|
|
|
1341 --annotation-color-by-strand <logical> [default: FALSE ]
|
|
|
1342 --annotation-placed-by-strand <logical> [default: FALSE ]
|
|
4
|
1343 -L/--smoothLength <integer> [default: NA ]
|
|
|
1344 ",stdout())
|
|
5
|
1345 q("no")
|
|
|
1346 }
|
|
|
1347 else {
|
|
1
|
1348 tmp=suppressPackageStartupMessages(require("Rsamtools"))
|
|
|
1349 if(!tmp) {
|
|
|
1350 stop("Package Rsamtools required !!")
|
|
|
1351 }
|
|
|
1352 tmp=suppressPackageStartupMessages(require("GenomicRanges"))
|
|
|
1353 if(!tmp) {
|
|
|
1354 stop("Package GenomicRanges required !!")
|
|
|
1355 }
|
|
5
|
1356 optArgs=getopt(rbind(
|
|
|
1357 c('output','o', 1, 'character',"output.png"),
|
|
|
1358 c('fileType','F',1,'character',NA),
|
|
|
1359 c('resolution','R',1,'numeric',72),
|
|
|
1360 c('typeVisu', 'v', 1, 'character', "classic"),
|
|
|
1361 c('description-data','t',1,'character',NA),
|
|
|
1362 c('chromosome-name', 'c', 1, 'character', NA),
|
|
|
1363 c('start', 'S', 1, 'numeric',1),
|
|
|
1364 c('end', 'E', 1, 'numeric',NA),
|
|
|
1365 c('annotation','a',1,'character',NA),
|
|
|
1366 c('typeTranscript','r',1,'character',NA),
|
|
|
1367 c('annotation-colors','C',1,'character',NA),
|
|
|
1368 c('annotation-shapes','s','1','character',NA),
|
|
|
1369 c('normalization-coefficient','n','1','character',NA),
|
|
|
1370 c('classic-plus-color',1,1,'character',"navyblue"),
|
|
|
1371 c('classic-minus-color',2,1,'character',"deeppink3"),
|
|
|
1372 c('heatmap-max-color',3,1,'character',"#000055"),
|
|
|
1373 c('heatmap-min-color',4,1,'character',"#FFFFAA"),
|
|
|
1374 c('heatmap-palette-method',7,1,'character',"hsv"),
|
|
|
1375 c('lines-samples-colors',5,1,'character',NA),
|
|
|
1376 c('lines-samples-type-line',6,1,'character',NA),
|
|
|
1377 c('scale-log', 'l',0,'logical',FALSE),
|
|
|
1378 c('inverseStrand','i',0,'logical', FALSE),
|
|
|
1379 c('unstranded','u',0,'logical', FALSE),
|
|
|
1380 c('symetric-scale','y',0,'logical', FALSE),
|
|
|
1381 c('annotation-color-by-strand',8,0,'logical',FALSE),
|
|
|
1382 c('annotation-placed-by-strand',9,0,'logical',FALSE),
|
|
|
1383 c('smoothLength','L',1,'numeric',NA)
|
|
|
1384 ))
|
|
|
1385 }
|
|
1
|
1386
|
|
|
1387 ###################
|
|
|
1388 ## ARGUMENTS
|
|
|
1389 #################################################################################
|
|
|
1390 files=optArgs$ARGUMENT
|
|
5
|
1391 ## Case file doesn't exist
|
|
|
1392 for( i in 1:length(files)) {
|
|
|
1393 if(!file.exists(files[i])) {
|
|
|
1394 stop(paste(files[i],"do not exist!","\n"))
|
|
1
|
1395 }
|
|
5
|
1396 }
|
|
1
|
1397 imagefile=optArgs$output
|
|
|
1398 typeVisu=optArgs$typeVisu
|
|
|
1399 description_data=optArgs$`description-data`
|
|
5
|
1400 chrName=optArgs$`chromosome-name`[1]
|
|
1
|
1401 start=optArgs$start[1]
|
|
|
1402 end=optArgs$end[1]
|
|
|
1403 annotation=optArgs$annotation
|
|
|
1404 typeTranscript=optArgs$typeTranscript
|
|
|
1405 colors=optArgs$`annotation-colors`
|
|
|
1406 shape_data=optArgs$`annotation-shapes`
|
|
|
1407 weight=optArgs$`normalization-coefficient`
|
|
|
1408 classic_plus_color=optArgs$`classic-plus-color`
|
|
|
1409 classic_minus_color=optArgs$`classic-minus-color`
|
|
|
1410 heatmap_max_color=optArgs$`heatmap-max-color`
|
|
|
1411 heatmap_min_color=optArgs$`heatmap-min-color`
|
|
|
1412 heatmap_palette_method=optArgs$`heatmap-palette-method`
|
|
|
1413 lines_samples_colors=optArgs$`lines-samples-colors`
|
|
|
1414 lines_samples_type_line=optArgs$`lines-samples-type-line`
|
|
|
1415 log=optArgs$`scale-log`
|
|
|
1416 inverseStrand=optArgs$inverseStrand
|
|
|
1417 unstranded=optArgs$unstranded
|
|
5
|
1418 symetric=optArgs$`symetric-scale`
|
|
1
|
1419 annotation_color_by_strand=optArgs$`annotation-color-by-strand`
|
|
|
1420 annotation_placed_by_strand=optArgs$`annotation-placed-by-strand`
|
|
4
|
1421 smoothLength=optArgs$`smoothLength`
|
|
5
|
1422 resolution=optArgs$`resolution`
|
|
|
1423 fileType=optArgs$`fileType`
|
|
|
1424
|
|
1
|
1425 ###################
|
|
|
1426 ## MAIN
|
|
|
1427 ###################################################################################
|
|
|
1428
|
|
5
|
1429 genome_info=scanBamHeader(files[1])[[1]]$targets
|
|
|
1430 noms_chromosomes=names(genome_info)
|
|
|
1431 longueur_chromosomes=as.integer(genome_info)
|
|
|
1432 nombre_chromosomes=length(noms_chromosomes)
|
|
1
|
1433
|
|
5
|
1434 ## Case no chromosome specified
|
|
|
1435 if(sum(is.na(chrName))) {
|
|
|
1436 chrName=noms_chromosomes[1]
|
|
|
1437 write(paste("No chromosome specified, processing chromosome :",chrName,"\n",sep=""),stderr())
|
|
|
1438 }
|
|
|
1439
|
|
|
1440 ## Case false chromosome name
|
|
|
1441 if(!(chrName %in% noms_chromosomes)) {
|
|
|
1442 stop(paste("\"",chrName,"\" is not a proper chromosome name"))
|
|
|
1443 }
|
|
1
|
1444
|
|
5
|
1445 if(is.na(end)) {
|
|
|
1446 end=longueur_chromosomes[chrName==noms_chromosomes]
|
|
|
1447 }
|
|
|
1448
|
|
|
1449 if(start > end) {
|
|
|
1450 stop("The start is bigger than the end!")
|
|
|
1451 }
|
|
1
|
1452
|
|
5
|
1453 ## Case asked coordinates outside the chromosome
|
|
|
1454 if(start<0|end>longueur_chromosomes[chrName==noms_chromosomes]) {
|
|
|
1455 stop("You are outside the chromosome")
|
|
|
1456 }
|
|
1
|
1457
|
|
5
|
1458 if(sum(is.na(weight))>0) {
|
|
|
1459 normalized_data=NULL
|
|
|
1460 isnormalized=FALSE
|
|
|
1461 }
|
|
|
1462 else {
|
|
|
1463 isnormalized=TRUE
|
|
|
1464 normalized_data=unlist(strsplit(weight,split=","))
|
|
|
1465 if(length(files)!=length(normalized_data)) {
|
|
|
1466 stop("Different number of files and weights ")
|
|
1
|
1467 }
|
|
|
1468 else {
|
|
5
|
1469 normalized_data=as.numeric(normalized_data)
|
|
1
|
1470 }
|
|
5
|
1471 }
|
|
1
|
1472
|
|
5
|
1473 if(inverseStrand) {
|
|
|
1474 libraryType="inverse"
|
|
|
1475 }
|
|
|
1476 else {
|
|
|
1477 libraryType="standard"
|
|
|
1478 }
|
|
|
1479 ## Read the bam file and extract the infos
|
|
|
1480 doit=function(i,libraryType) {
|
|
|
1481 try(readBam_(files[i], libraryType=libraryType, chrName_=chrName, from_=start, to_=end,normalized_=normalized_data[i]))
|
|
|
1482 }
|
|
|
1483 data=lapply(1:length(files),doit,libraryType=libraryType)
|
|
|
1484 ctrl=unlist(lapply(data,class))=="try-error"
|
|
|
1485
|
|
|
1486 if(sum(ctrl)>0) {
|
|
|
1487 for(i in which(ctrl)) {
|
|
|
1488 write(paste("Problem with file :",files[i],"\n",sep=""),stderr())
|
|
|
1489 }
|
|
|
1490 stop("At least a file has problem")
|
|
|
1491 }
|
|
|
1492
|
|
|
1493 ## Read the GFF file and extract the infos
|
|
|
1494 if(sum(is.na(annotation))==0) {
|
|
|
1495 gff=try(readGff(annotation[1],from=start,to=end,chr=chrName),TRUE)
|
|
|
1496 ctrl=class(gff)=="try-error"
|
|
|
1497 if(sum(ctrl)>0) {
|
|
|
1498 stop(paste("Problem with Gff file :",annotation,"\n"))
|
|
|
1499 }
|
|
|
1500 if(length(annotation)>1) {
|
|
|
1501 for(i in 2:length(annotation)) {
|
|
|
1502 gff1=try(readGff(annotation[i],from=start,to=end,chr=chrName),TRUE)
|
|
|
1503 ctrl=class(gff1)=="try-error"
|
|
|
1504 if(sum(ctrl)>0) {
|
|
|
1505 stop(paste("Problem with Gff file :",gff1,"\n"))
|
|
|
1506 }
|
|
|
1507 if(sum(is.na(gff1))==0) {
|
|
1
|
1508 gff=rbind(gff,gff1)
|
|
|
1509 }
|
|
|
1510 }
|
|
|
1511 }
|
|
5
|
1512 }
|
|
|
1513 else {
|
|
|
1514 gff=NA
|
|
|
1515 }
|
|
|
1516
|
|
|
1517 if(sum(is.na(description_data))>0) {
|
|
|
1518 description_data=files
|
|
|
1519 }
|
|
4
|
1520
|
|
5
|
1521 ## Case of different number of files and description
|
|
|
1522 if(length(description_data)!=length(files)) {
|
|
|
1523 stop("Different number of files and description")
|
|
|
1524 }
|
|
1
|
1525
|
|
5
|
1526 ## Pooling bam files (if necessary)
|
|
|
1527 description_data_unique=unique(description_data)
|
|
|
1528 data_pooled=list()
|
|
|
1529 for(i_descri in description_data_unique) {
|
|
|
1530 i_data_voulu=which(description_data==i_descri)
|
|
|
1531 data_pooled_=data[[i_data_voulu[1]]]
|
|
|
1532 if(length(i_data_voulu)>1) {
|
|
|
1533 for(i_i_data_voulu in 2:length(i_data_voulu)) {
|
|
|
1534 data_pooled_=addBam(data_pooled_,data[[i_data_voulu[i_i_data_voulu]]])
|
|
1
|
1535 }
|
|
|
1536 }
|
|
5
|
1537 data_pooled[[i_descri]]=data_pooled_
|
|
|
1538 }
|
|
|
1539
|
|
|
1540 if(sum(is.na(typeTranscript)>0)) {
|
|
|
1541 if(sum(is.na(annotation))==0) {
|
|
|
1542 typeTranscritSplit=unique(gff$Type)
|
|
|
1543 }
|
|
1
|
1544 else {
|
|
5
|
1545 typeTranscritSplit=NA
|
|
1
|
1546 }
|
|
5
|
1547 }
|
|
|
1548 else {
|
|
|
1549 typeTranscritSplit=unlist(strsplit(typeTranscript,split=","))
|
|
|
1550 }
|
|
|
1551
|
|
|
1552 if((sum(is.na(gff))>0)|(sum(is.na(typeTranscritSplit))>0)) {
|
|
|
1553 annot=NULL
|
|
|
1554 }
|
|
|
1555 else {
|
|
|
1556 annot=gff[gff$Type %in% typeTranscritSplit,]
|
|
|
1557 }
|
|
|
1558
|
|
|
1559 ## Check the colors
|
|
|
1560 if(sum(is.na(colors)>0)) {
|
|
|
1561 nTypeTranscrit=length(typeTranscritSplit)
|
|
|
1562 if(annotation_color_by_strand) {
|
|
|
1563 colorsSplit=c(classic_plus_color,classic_minus_color)
|
|
1
|
1564 }
|
|
5
|
1565 else {
|
|
|
1566 colorsSplit=rainbow(nTypeTranscrit+1)[1:nTypeTranscrit]
|
|
|
1567 }
|
|
|
1568 }
|
|
|
1569 else {
|
|
|
1570 colorsSplit=unlist(strsplit(colors,split=","))
|
|
|
1571 }
|
|
|
1572 for(i in 1:length(colorsSplit)) {
|
|
|
1573 tmp=unlist(strsplit(colorsSplit[i],split=""))
|
|
1
|
1574 if(length(tmp)==6|length(tmp)==8) {
|
|
|
1575 if(sum(tmp %in% c(0:9,"A","B","C","D","E","F"))==length(tmp)) {
|
|
5
|
1576 colorsSplit[i]=paste("#",colorsSplit[i],sep="")
|
|
1
|
1577 }
|
|
|
1578 }
|
|
5
|
1579 }
|
|
|
1580
|
|
|
1581 tmp=unlist(strsplit(classic_plus_color,split=""))
|
|
|
1582 if(length(tmp)==6|length(tmp)==8) {
|
|
|
1583 if(sum(tmp %in% c(0:9,"A","B","C","D","E","F"))==length(tmp)) {
|
|
|
1584 classic_plus_color=paste("#",classic_plus_color,sep="")
|
|
|
1585 }
|
|
|
1586 }
|
|
|
1587 tmp=unlist(strsplit(classic_minus_color,split=""))
|
|
|
1588 if(length(tmp)==6|length(tmp)==8) {
|
|
|
1589 if(sum(tmp %in% c(0:9,"A","B","C","D","E","F"))==length(tmp)) {
|
|
|
1590 classic_minus_color=paste("#",classic_minus_color,sep="")
|
|
|
1591 }
|
|
|
1592 }
|
|
|
1593 tmp=unlist(strsplit(heatmap_max_color,split=""))
|
|
|
1594 if(length(tmp)==6|length(tmp)==8) {
|
|
|
1595 if(sum(tmp %in% c(0:9,"A","B","C","D","E","F"))==length(tmp)) {
|
|
|
1596 heatmap_max_color=paste("#",heatmap_max_color,sep="")
|
|
|
1597 }
|
|
|
1598 }
|
|
|
1599 tmp=unlist(strsplit(heatmap_min_color,split=""))
|
|
|
1600 if(length(tmp)==6|length(tmp)==8) {
|
|
|
1601 if(sum(tmp %in% c(0:9,"A","B","C","D","E","F"))==length(tmp)) {
|
|
|
1602 heatmap_min_color=paste("#",heatmap_min_color,sep="")
|
|
1
|
1603 }
|
|
5
|
1604 }
|
|
|
1605 if(sum(is.na(lines_samples_colors))>0) {
|
|
|
1606 lines_samples_colors_split=c(1,3,4,2)[((0:(length(files)-1))%%4)+1]
|
|
|
1607 }
|
|
|
1608 else {
|
|
|
1609 lines_samples_colors_split=unlist(strsplit(lines_samples_colors,split=","))
|
|
|
1610 for(i in 1:length(lines_samples_colors_split)) {
|
|
|
1611 tmp=unlist(strsplit(lines_samples_colors_split[i],split=""))
|
|
|
1612 if(length(tmp)==6|length(tmp)==8) {
|
|
|
1613 if(sum(tmp %in% c(0:9,"A","B","C","D","E","F"))==length(tmp)) {
|
|
|
1614 lines_samples_colors_split[i]=paste("#",lines_samples_colors_split[i],sep="")
|
|
|
1615 }
|
|
|
1616 }
|
|
|
1617 }
|
|
|
1618 }
|
|
|
1619 colorToCheck=c(colorsSplit,classic_plus_color,classic_minus_color,heatmap_max_color,heatmap_min_color,lines_samples_colors_split)
|
|
|
1620 ctrl=sapply(colorToCheck,is.acceptable.color)
|
|
|
1621 if(sum(!ctrl)>0) {
|
|
|
1622 for(i in which(!ctrl)) {
|
|
|
1623 write(paste("\"",colorToCheck[i],"\" is not a proper color name.\n",sep=""),stderr())
|
|
|
1624 }
|
|
|
1625 stop("At least one color has a problem")
|
|
|
1626 }
|
|
|
1627 if(annotation_color_by_strand) {
|
|
|
1628 if(length(colorsSplit)>2) {
|
|
|
1629 stop("You have to specify two and only two colors!")
|
|
|
1630 }
|
|
|
1631 }
|
|
|
1632 else {
|
|
|
1633 if(length(typeTranscritSplit)!=length(colorsSplit)) {
|
|
|
1634 stop("Please specify the same number of transcript types and colors")
|
|
|
1635 }
|
|
|
1636 }
|
|
|
1637 ## Check the line types
|
|
|
1638 if(sum(is.na(lines_samples_type_line))>0) {
|
|
|
1639 lines_samples_type_line_split=((0:(length(files)-1))%/%4)+1
|
|
|
1640 }
|
|
|
1641 else {
|
|
|
1642 lines_samples_type_line_split=unlist(strsplit(lines_samples_type_line,split=","))
|
|
|
1643 }
|
|
|
1644 if(typeVisu=="lines") {
|
|
|
1645 ctrl=sapply(lines_samples_type_line_split,function(x) {
|
|
|
1646 tmp=suppressWarnings(as.numeric(x))
|
|
|
1647 if(!is.na(tmp)) {
|
|
|
1648 return((tmp==floor(tmp))&tmp>=1&tmp<=5)
|
|
1
|
1649 }
|
|
5
|
1650 else {
|
|
|
1651 return(FALSE)
|
|
1
|
1652 }
|
|
5
|
1653 })
|
|
1
|
1654 if(sum(!ctrl)>0) {
|
|
|
1655 for(i in which(!ctrl)) {
|
|
5
|
1656 write(paste("\"",lines_samples_type_line_split[i],"\" is not a proper line style.\n",sep=""),stderr())
|
|
1
|
1657 }
|
|
5
|
1658 stop("At least one line style has problem")
|
|
1
|
1659 }
|
|
5
|
1660 lines_samples_type_line_split=as.integer(lines_samples_type_line_split)
|
|
|
1661 }
|
|
|
1662 ## Check the shapes
|
|
1
|
1663 type_shape=rep(1,length(typeTranscritSplit))
|
|
|
1664
|
|
|
1665 if(sum(is.na(shape_data)>0)) {
|
|
5
|
1666 for( i in 1:length(typeTranscritSplit)) {
|
|
|
1667 type_shape[i]="box"
|
|
|
1668 }
|
|
|
1669 }
|
|
|
1670 else {
|
|
|
1671 shape=unlist(strsplit(shape_data,split=","))
|
|
|
1672 shape=as.array(shape)
|
|
|
1673 if(length(typeTranscritSplit)!=length(shape)) {
|
|
|
1674 stop("Please specify the same number of transcript types and shapes")
|
|
|
1675 }
|
|
|
1676 else {
|
|
1
|
1677 for( i in 1:length(typeTranscritSplit)) {
|
|
5
|
1678 type_shape[i]=shape[[i]]
|
|
1
|
1679 }
|
|
|
1680 }
|
|
|
1681 }
|
|
|
1682
|
|
5
|
1683 ## Style for the annotation
|
|
1
|
1684 label=rep(1,length(typeTranscritSplit))
|
|
5
|
1685 style=data.frame(Type=c(typeTranscritSplit,typeTranscritSplit),Strand=c(rep("+",length(typeTranscritSplit)),rep("-",length(typeTranscritSplit))),col=NA,shape=NA,label,stringsAsFactors=FALSE)
|
|
1
|
1686 for(i in 1:length(typeTranscritSplit)) {
|
|
|
1687 style$shape[style$Type==typeTranscritSplit[i]]=type_shape[i]
|
|
|
1688 }
|
|
|
1689 if(annotation_color_by_strand) {
|
|
|
1690 style$col[style$Strand=="+"]=colorsSplit[1]
|
|
|
1691 style$col[style$Strand=="-"]=colorsSplit[2]
|
|
|
1692 }
|
|
|
1693 else {
|
|
|
1694 for(i in 1:length(typeTranscritSplit)) {
|
|
|
1695 style$col[style$Type==typeTranscritSplit[i]]=colorsSplit[i]
|
|
|
1696 }
|
|
|
1697 }
|
|
|
1698
|
|
5
|
1699 ## Main for visualisation
|
|
4
|
1700 databychr=extractSignal(data_pooled,chrName,from=start,to=end,normalized_=isnormalized)
|
|
1
|
1701
|
|
|
1702 {
|
|
|
1703 reverse=list()
|
|
|
1704 forward=list()
|
|
|
1705 if(log) {
|
|
5
|
1706 for(i in 1:length(databychr$F)) {
|
|
|
1707 forward_=numeric()
|
|
|
1708 tmp=log2(1+databychr$F[[i]])
|
|
|
1709 forward_[1:length(tmp)]=tmp
|
|
|
1710 forward[[i]]=forward_
|
|
|
1711 reverse_=numeric()
|
|
|
1712 tmp=log2(1+databychr$R[[i]])
|
|
|
1713 reverse_[1:length(tmp)]=tmp
|
|
|
1714 reverse[[i]]=reverse_
|
|
|
1715 }
|
|
1
|
1716 }
|
|
|
1717 else {
|
|
5
|
1718 forward=databychr$F
|
|
|
1719 reverse=databychr$R
|
|
1
|
1720 }
|
|
|
1721 }
|
|
|
1722
|
|
4
|
1723 if(sum(is.na(smoothLength))>0) {
|
|
|
1724 smoothLength=max(trunc((end-start)/1200),1)
|
|
|
1725 }
|
|
|
1726 else {
|
|
|
1727 smoothLength=smoothLength[1]
|
|
|
1728 }
|
|
1
|
1729 if(unstranded) {
|
|
5
|
1730 for(i in 1:length(databychr$F)) {
|
|
|
1731 if(smoothLength>0) {
|
|
|
1732 lo=smooth(forward[[i]]+reverse[[i]],L=smoothLength)
|
|
|
1733 }
|
|
|
1734 else {
|
|
|
1735 lo=forward[[i]]+reverse[[i]]
|
|
|
1736 }
|
|
1
|
1737 forward[[i]]=lo
|
|
|
1738 los=rep(0,length(lo))
|
|
|
1739 reverse[[i]]=los
|
|
|
1740 }
|
|
|
1741 }
|
|
|
1742 else {
|
|
5
|
1743 for(i in 1:length(databychr$F)) {
|
|
|
1744 if(smoothLength>0) {
|
|
|
1745 lo=smooth(forward[[i]],L=smoothLength)
|
|
|
1746 }
|
|
|
1747 else {
|
|
|
1748 lo=forward[[i]]
|
|
|
1749 }
|
|
1
|
1750 forward[[i]]=lo
|
|
5
|
1751 if(smoothLength>0) {
|
|
|
1752 los=smooth(reverse[[i]],L=smoothLength)
|
|
|
1753 }
|
|
|
1754 else {
|
|
|
1755 los=reverse[[i]]
|
|
|
1756 }
|
|
1
|
1757 reverse[[i]]=los
|
|
|
1758 }
|
|
|
1759 }
|
|
|
1760
|
|
5
|
1761 group_maximum=rep(1,length(databychr$F))
|
|
|
1762 max_forward=numeric(length(databychr$F))
|
|
|
1763 max_reverse=numeric(length(databychr$F))
|
|
1
|
1764 for(i_data in 1:length(databychr$F)) {
|
|
|
1765 max_forward[i_data]=max(forward[[i_data]],na.rm=TRUE)
|
|
|
1766 max_reverse[i_data]=max(reverse[[i_data]],na.rm=TRUE)
|
|
|
1767 }
|
|
|
1768 for(i_max in unique(group_maximum)) {
|
|
|
1769 max_forward[group_maximum==i_max]=max(max_forward[group_maximum==i_max],na.rm=TRUE)
|
|
|
1770 max_reverse[group_maximum==i_max]=max(max_reverse[group_maximum==i_max],na.rm=TRUE)
|
|
|
1771 }
|
|
|
1772
|
|
5
|
1773 if(symetric) {
|
|
1
|
1774 for(i_data in 1:length(databychr$F)) {
|
|
|
1775 max_forward[i_data]=max(max_forward[i_data],max_reverse[i_data],na.rm=TRUE)
|
|
|
1776 }
|
|
|
1777 max_reverse=max_forward
|
|
|
1778 }
|
|
|
1779
|
|
5
|
1780 if(is.na(fileType)) {
|
|
|
1781 imageFileSplit=unlist(strsplit(imagefile,split=".",fixed=TRUE))
|
|
|
1782 if(length(imageFileSplit)>=2) {
|
|
|
1783 fileType=imageFileSplit[length(imageFileSplit)]
|
|
|
1784 }
|
|
|
1785 else {
|
|
|
1786 fileType="png"
|
|
|
1787 }
|
|
|
1788 }
|
|
|
1789
|
|
1
|
1790 plotVisu(imagefile,typeVisu=typeVisu,listForward=forward,listReverse=reverse,
|
|
|
1791 debut_vue=start,fin_vue=end,chr=chrName,annot=annot,repeated=TRUE,
|
|
4
|
1792 titres=description_data_unique,name_flags="",style=style,log=log,stranded=!unstranded,
|
|
1
|
1793 tops=max_forward,bottoms=max_reverse,
|
|
|
1794 classic_plus_color=classic_plus_color,classic_minus_color=classic_minus_color,
|
|
|
1795 heatmap_max_color=heatmap_max_color,heatmap_min_color=heatmap_min_color,heatmap_palette_method=heatmap_palette_method,
|
|
|
1796 lines_samples_colors=lines_samples_colors_split,lines_samples_type_line=lines_samples_type_line_split,
|
|
5
|
1797 smoothLength=1,annotation_color_by_strand=annotation_color_by_strand,annotation_placed_by_strand=annotation_placed_by_strand,
|
|
|
1798 fileType=fileType,resolutionDPI=resolution)
|
|
1
|
1799
|
|
|
1800 }
|
