Mercurial > repos > brasset_jensen > srnapipe

comparison bin/Rcall.pm @ 61:9185ca0a7b43 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Updated package according to recommendations.
author pierre.pouchin
date Wed, 16 Jan 2019 08:18:13 -0500
parents 4bc00caa60b4
children
comparison
equal deleted inserted replaced
60:9645d995fb3c 61:9185ca0a7b43
8 our @ISA = qw(Exporter); 8 our @ISA = qw(Exporter);
9 our @EXPORT_OK = qw( &histogram &pie_chart &bg_to_png ); 9 our @EXPORT_OK = qw( &histogram &pie_chart &bg_to_png );
10 10
11 sub histogram 11 sub histogram
12 { 12 {
13 my ($size_hashR, $out_png, $size) = @_; 13 my ($size_hashR, $out_png, $size) = @_;
14 my (@abs, @ord); 14 my (@abs, @ord);
15 my $i = 0; 15 my $i = 0;
16 foreach my $k (sort {$a <=> $b} keys %{$size_hashR}) 16 foreach my $k (sort {$a <=> $b} keys %{$size_hashR})
17 { 17 {
18 my $percentage = 0; 18 my $percentage = 0;
19 $percentage = $size_hashR->{$k} * 100 / $size if $size != 0; 19 $percentage = $size_hashR->{$k} * 100 / $size if $size != 0;
20 $abs[$i] = $k ; $ord[$i] = $percentage; $i++; 20 $abs[$i] = $k ; $ord[$i] = $percentage; $i++;
21 } 21 }
22 my $abs = join (",", @abs ); 22 my $abs = join (",", @abs );
23 my $ord = join (",", @ord ); 23 my $ord = join (",", @ord );
24 if (scalar(@abs) != 0) 24 if (scalar(@abs) != 0)
25 { 25 {
26
27 my $R = Statistics::R->new();
28 $R->startR;
29 $R->send(
30 qq`library(ggplot2)
31 percentage = c($ord)
32 size =c($abs)
33 min = min(size)
34 max = max(size)
35 dat = data.frame(size,percentage)
36 png(filename=\"$out_png\", width = 640, height = 640)
37 c = ggplot(dat,aes(size,percentage))
38 c + geom_bar(stat="identity") + scale_x_continuous(breaks=min:max)+theme( axis.text.x = element_text(angle=90, hjust=0.5, size=20), axis.text.y = element_text( size=20 ), axis.title.x = element_text( size=25, face="bold"), axis.title.y = element_text( size=25, face="bold") )
39 dev.off()`);
40 $R->stopR();
41
42 }
43 }
44
45 sub bg_to_png
46 {
47 my ( $fai, $bgP, $bgM, $dir, $sb ) = @_;
48 my $R = Statistics::R->new();
49 $R->startR;
50 $R->send(
51 qq`library('Sushi')
52 fai =read.table("$fai")
53 if ( file.info("$bgP")\$size !=0 )
54 {
55 bgP = read.table("$bgP")
56 } else { bgP = data.frame(factor(),integer()) }
57
58 if ( file.info("$bgM")\$size !=0 )
59 {
60 bgM = read.table("$bgM")
61 } else { bgM = data.frame(factor(),integer()) }
62
63 f_both = function(chr,end) {
64 jpeg( paste0("$dir",as.character(chr),".png"), quality=100)
65 par(mfrow=c(2,1),mar=c(1,10,1,3))
66 plotBedgraph(bgP, chrom=chr,chromstart=0,chromend=end,transparency=.50, color=SushiColors(2)(2)[1])
67 axis(side=2,las=2,tcl=.2)
68 mtext("Scaled Read Depth",side=2,line=4,cex=1,font=2)
69 plotBedgraph(bgM, chrom=chr,chromstart=0,chromend=end,transparency=.50, flip=TRUE, color=SushiColors(2)(2)[2])
70 labelgenome(chrom=chr,chromstart=0,chromend=end,side=3,n=3,scale="$sb", line=0, chromline = 0.5, scaleline = 0.5, scaleadjust =1.05, chromadjust = -0.4)
71 axis(side=2,las=2,tcl=.2,at=pretty(par("yaxp")[c(1,2)]),labels=-1*pretty(par("yaxp")[c(1,2)]))
72 mtext("Scaled Read Depth",side=2,line=4.5,cex=1,font=2)
73 dev.off()
74 }
75
76 f_plus = function(chr,end) {
77 jpeg( paste0("$dir",as.character(chr),".png"), quality=100)
78 plotBedgraph(bgP, chrom=chr,chromstart=0,chromend=end,transparency=.50, color=SushiColors(2)(2)[1])
79 labelgenome(chrom=chr,chromstart=0,chromend=end,n=3,scale="$sb", line=0, chromline = 0.5, scaleline = 0.5, scaleadjust =1.05, chromadjust = -0.4)
80 axis(side=2,las=2,tcl=.2)
81 mtext("Scaled Read Depth",side=2,line=4,cex=1,font=2)
82 dev.off()
83 }
84
85 f_minus = function(chr,end) {
86 jpeg( paste0("$dir",as.character(chr),".png"), quality=100)
87 plotBedgraph(bgM, chrom=chr,chromstart=0,chromend=end,transparency=.50, flip=TRUE, color=SushiColors(2)(2)[2])
88 labelgenome(chrom=chr,chromstart=0,chromend=end,n=3,scale="$sb", line=0, chromline = 0.5, scaleline = 0.5, scaleadjust =1.05, chromadjust = -0.4)
89 axis(side=2,las=2,tcl=.2,at=pretty(par("yaxp")[c(1,2)]),labels=-1*pretty(par("yaxp")[c(1,2)]))
90 mtext("Scaled Read Depth",side=2,line=4.5,cex=1,font=2)
91 dev.off()
92 }
93
94 fai_b = fai[fai\$V1 %in% intersect(bgM\$V1,bgP\$V1), ]
95 mapply( f_both, fai_b\$V1, fai_b\$V2)
96
97 fai_p = fai[fai\$V1 %in% setdiff(bgP\$V1,bgM\$V1), ]
98 mapply( f_plus, fai_p\$V1, fai_p\$V2)
99
100 fai_m = fai[fai\$V1 %in% setdiff(bgM\$V1,bgP\$V1), ]
101 mapply( f_minus, fai_m\$V1, fai_m\$V2) `);
102
103 $R->stopR();
104 }
105
106 sub pie_chart
107 {
108 my $dir = shift;
109 my $in = $dir.'repartition.txt';
110 my $out = $dir.'pie_chart.png';
26 111
27 my $R = Statistics::R->new(); 112 my $R = Statistics::R->new();
28 $R->startR; 113 $R->startR;
29 $R->send( 114 $R->send(
30 qq`library(ggplot2) 115 qq`
31 percentage = c($ord) 116 library(plotrix)
32 size =c($abs) 117 library(RColorBrewer)
33 min = min(size) 118 R =read.table("$in",header=T)
34 max = max(size) 119 values = round(R\$percentage)
35 dat = data.frame(size,percentage) 120 keys = R\$type
36 png(filename=\"$out_png\", width = 640, height = 640) 121 lab = paste(values, "%", sep="")
37 c = ggplot(dat,aes(size,percentage)) 122 png("$out")
38 c + geom_bar(stat="identity") + scale_x_continuous(breaks=min:max)+theme( axis.text.x = element_text(angle=90, hjust=0.5, size=20), axis.text.y = element_text( size=20 ), axis.title.x = element_text( size=25, face="bold"), axis.title.y = element_text( size=25, face="bold") ) 123 colors <- brewer.pal(7,"Paired")
39 dev.off()`); 124 pie(values, col=colors, labels=lab, clockwise=TRUE)
125 legend("bottom", legend = keys, fill=colors, bty="n", ncol = 3)
126 par(mai = c(0,0,0,0))
127 layout(c(1,2),heights=c(0.3,1))
128 plot.new()
129 legend("bottom", legend = keys, fill=colors, bty="n",ncol = 3)
130 pie(values, col=colors, labels=lab, clockwise=TRUE)
131 dev.off()`
132 );
40 $R->stopR(); 133 $R->stopR();
41
42 }
43 }
44
45 sub bg_to_png
46 {
47 my ( $fai, $bgP, $bgM, $dir, $sb ) = @_;
48 my $R = Statistics::R->new();
49 $R->startR;
50 $R->send(
51 qq`library('Sushi')
52 fai =read.table("$fai")
53 if ( file.info("$bgP")\$size !=0 )
54 {
55 bgP = read.table("$bgP")
56 } else { bgP = data.frame(factor(),integer()) }
57
58 if ( file.info("$bgM")\$size !=0 )
59 {
60 bgM = read.table("$bgM")
61 } else { bgM = data.frame(factor(),integer()) }
62
63 f_both = function(chr,end) {
64 jpeg( paste0("$dir",as.character(chr),".png"), quality=100)
65 par(mfrow=c(2,1),mar=c(1,10,1,3))
66 plotBedgraph(bgP, chrom=chr,chromstart=0,chromend=end,transparency=.50, color=SushiColors(2)(2)[1])
67 axis(side=2,las=2,tcl=.2)
68 mtext("Scaled Read Depth",side=2,line=4,cex=1,font=2)
69 plotBedgraph(bgM, chrom=chr,chromstart=0,chromend=end,transparency=.50, flip=TRUE, color=SushiColors(2)(2)[2])
70 labelgenome(chrom=chr,chromstart=0,chromend=end,side=3,n=3,scale="$sb", line=0, chromline = 0.5, scaleline = 0.5, scaleadjust =1.05, chromadjust = -0.4)
71 axis(side=2,las=2,tcl=.2,at=pretty(par("yaxp")[c(1,2)]),labels=-1*pretty(par("yaxp")[c(1,2)]))
72 mtext("Scaled Read Depth",side=2,line=4.5,cex=1,font=2)
73 dev.off()
74 }
75
76 f_plus = function(chr,end) {
77 jpeg( paste0("$dir",as.character(chr),".png"), quality=100)
78 plotBedgraph(bgP, chrom=chr,chromstart=0,chromend=end,transparency=.50, color=SushiColors(2)(2)[1])
79 labelgenome(chrom=chr,chromstart=0,chromend=end,n=3,scale="$sb", line=0, chromline = 0.5, scaleline = 0.5, scaleadjust =1.05, chromadjust = -0.4)
80 axis(side=2,las=2,tcl=.2)
81 mtext("Scaled Read Depth",side=2,line=4,cex=1,font=2)
82 dev.off()
83 }
84
85 f_minus = function(chr,end) {
86 jpeg( paste0("$dir",as.character(chr),".png"), quality=100)
87 plotBedgraph(bgM, chrom=chr,chromstart=0,chromend=end,transparency=.50, flip=TRUE, color=SushiColors(2)(2)[2])
88 labelgenome(chrom=chr,chromstart=0,chromend=end,n=3,scale="$sb", line=0, chromline = 0.5, scaleline = 0.5, scaleadjust =1.05, chromadjust = -0.4)
89 axis(side=2,las=2,tcl=.2,at=pretty(par("yaxp")[c(1,2)]),labels=-1*pretty(par("yaxp")[c(1,2)]))
90 mtext("Scaled Read Depth",side=2,line=4.5,cex=1,font=2)
91 dev.off()
92 }
93
94 fai_b = fai[fai\$V1 %in% intersect(bgM\$V1,bgP\$V1), ]
95 mapply( f_both, fai_b\$V1, fai_b\$V2)
96
97 fai_p = fai[fai\$V1 %in% setdiff(bgP\$V1,bgM\$V1), ]
98 mapply( f_plus, fai_p\$V1, fai_p\$V2)
99
100 fai_m = fai[fai\$V1 %in% setdiff(bgM\$V1,bgP\$V1), ]
101 mapply( f_minus, fai_m\$V1, fai_m\$V2) `);
102
103 $R->stopR();
104 }
105
106 sub pie_chart
107 {
108 my $dir = shift;
109 my $in = $dir.'repartition.txt';
110 my $out = $dir.'pie_chart.png';
111
112 my $R = Statistics::R->new();
113 $R->startR;
114 $R->send(
115 qq`
116 library(plotrix)
117 library(RColorBrewer)
118 R =read.table("$in",header=T)
119 values = round(R\$percentage)
120 keys = R\$type
121 lab = paste(values, "%", sep="")
122 png("$out")
123 colors <- brewer.pal(7,"Paired")
124 pie(values, col=colors, labels=lab, clockwise=TRUE)
125 legend("bottom", legend = keys, fill=colors, bty="n", ncol = 3)
126 par(mai = c(0,0,0,0))
127 layout(c(1,2),heights=c(0.3,1))
128 plot.new()
129 legend("bottom", legend = keys, fill=colors, bty="n",ncol = 3)
130 pie(values, col=colors, labels=lab, clockwise=TRUE)
131 dev.off()`
132 );
133 $R->stopR();
134 } 134 }
135 135
136 1; 136 1;