Mercurial > repos > rlegendre > ving

comparison ving.R @ 1:5feff08a06ea draft

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