annotate bin/Rcall.pm @ 0:198009598544 draft

Uploaded
author romaingred
date Wed, 11 Oct 2017 09:57:58 -0400
parents
children ee6b4b2072a9
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
198009598544 Uploaded
romaingred
parents:
diff changeset
1 package Rcall;
198009598544 Uploaded
romaingred
parents:
diff changeset
2
198009598544 Uploaded
romaingred
parents:
diff changeset
3 use strict;
198009598544 Uploaded
romaingred
parents:
diff changeset
4 use warnings;
198009598544 Uploaded
romaingred
parents:
diff changeset
5 use Statistics::R;
198009598544 Uploaded
romaingred
parents:
diff changeset
6
198009598544 Uploaded
romaingred
parents:
diff changeset
7 use Exporter;
198009598544 Uploaded
romaingred
parents:
diff changeset
8 our @ISA = qw(Exporter);
198009598544 Uploaded
romaingred
parents:
diff changeset
9 our @EXPORT_OK = qw( &histogram &pie_chart &bg_to_png );
198009598544 Uploaded
romaingred
parents:
diff changeset
10
198009598544 Uploaded
romaingred
parents:
diff changeset
11 sub histogram
198009598544 Uploaded
romaingred
parents:
diff changeset
12 {
198009598544 Uploaded
romaingred
parents:
diff changeset
13 my ($size_hashR, $out_png, $size) = @_;
198009598544 Uploaded
romaingred
parents:
diff changeset
14 my (@abs, @ord);
198009598544 Uploaded
romaingred
parents:
diff changeset
15 my $i = 0;
198009598544 Uploaded
romaingred
parents:
diff changeset
16 foreach my $k (sort {$a <=> $b} keys %{$size_hashR})
198009598544 Uploaded
romaingred
parents:
diff changeset
17 {
198009598544 Uploaded
romaingred
parents:
diff changeset
18 my $percentage = $size_hashR->{$k} * 100 / $size;
198009598544 Uploaded
romaingred
parents:
diff changeset
19 $abs[$i] = $k ; $ord[$i] = $percentage; $i++;
198009598544 Uploaded
romaingred
parents:
diff changeset
20 }
198009598544 Uploaded
romaingred
parents:
diff changeset
21 my $abs = join (",", @abs );
198009598544 Uploaded
romaingred
parents:
diff changeset
22 my $ord = join (",", @ord );
198009598544 Uploaded
romaingred
parents:
diff changeset
23 if (scalar(@abs) != 0)
198009598544 Uploaded
romaingred
parents:
diff changeset
24 {
198009598544 Uploaded
romaingred
parents:
diff changeset
25
198009598544 Uploaded
romaingred
parents:
diff changeset
26 my $R = Statistics::R->new();
198009598544 Uploaded
romaingred
parents:
diff changeset
27 $R->startR;
198009598544 Uploaded
romaingred
parents:
diff changeset
28 $R->send(
198009598544 Uploaded
romaingred
parents:
diff changeset
29 qq`library(ggplot2)
198009598544 Uploaded
romaingred
parents:
diff changeset
30 percentage = c($ord)
198009598544 Uploaded
romaingred
parents:
diff changeset
31 size =c($abs)
198009598544 Uploaded
romaingred
parents:
diff changeset
32 min = min(size)
198009598544 Uploaded
romaingred
parents:
diff changeset
33 max = max(size)
198009598544 Uploaded
romaingred
parents:
diff changeset
34 dat = data.frame(size,percentage)
198009598544 Uploaded
romaingred
parents:
diff changeset
35 png(filename=\"$out_png\", width = 800, height = 480)
198009598544 Uploaded
romaingred
parents:
diff changeset
36 c = ggplot(dat,aes(size,percentage))
198009598544 Uploaded
romaingred
parents:
diff changeset
37 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") )
198009598544 Uploaded
romaingred
parents:
diff changeset
38 dev.off()`);
198009598544 Uploaded
romaingred
parents:
diff changeset
39 $R->stopR();
198009598544 Uploaded
romaingred
parents:
diff changeset
40
198009598544 Uploaded
romaingred
parents:
diff changeset
41 }
198009598544 Uploaded
romaingred
parents:
diff changeset
42 }
198009598544 Uploaded
romaingred
parents:
diff changeset
43
198009598544 Uploaded
romaingred
parents:
diff changeset
44 sub bg_to_png
198009598544 Uploaded
romaingred
parents:
diff changeset
45 {
198009598544 Uploaded
romaingred
parents:
diff changeset
46 my ( $fai, $bgP, $bgM, $dir, $sb ) = @_;
198009598544 Uploaded
romaingred
parents:
diff changeset
47 my $R = Statistics::R->new();
198009598544 Uploaded
romaingred
parents:
diff changeset
48 $R->startR;
198009598544 Uploaded
romaingred
parents:
diff changeset
49 $R->send(
198009598544 Uploaded
romaingred
parents:
diff changeset
50 qq`library('Sushi')
198009598544 Uploaded
romaingred
parents:
diff changeset
51 fai =read.table("$fai")
198009598544 Uploaded
romaingred
parents:
diff changeset
52 bgP = read.table("$bgP")
198009598544 Uploaded
romaingred
parents:
diff changeset
53 bgM = read.table("$bgM")
198009598544 Uploaded
romaingred
parents:
diff changeset
54 fun = function(chr,end) {
198009598544 Uploaded
romaingred
parents:
diff changeset
55 jpeg( paste0("$dir",as.character(chr),".png"), quality=100)
198009598544 Uploaded
romaingred
parents:
diff changeset
56 par(mfrow=c(2,1),mar=c(1,10,1,3))
198009598544 Uploaded
romaingred
parents:
diff changeset
57 plotBedgraph(bgP, chrom=chr,chromstart=0,chromend=end,transparency=.50, color=SushiColors(2)(2)[1])
198009598544 Uploaded
romaingred
parents:
diff changeset
58 axis(side=2,las=2,tcl=.2)
198009598544 Uploaded
romaingred
parents:
diff changeset
59 mtext("Scaled Read Depth",side=2,line=4,cex=1,font=2)
198009598544 Uploaded
romaingred
parents:
diff changeset
60 plotBedgraph(bgM, chrom=chr,chromstart=0,chromend=end,transparency=.50, flip=TRUE, color=SushiColors(2)(2)[2])
198009598544 Uploaded
romaingred
parents:
diff changeset
61 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)
198009598544 Uploaded
romaingred
parents:
diff changeset
62 axis(side=2,las=2,tcl=.2,at=pretty(par("yaxp")[c(1,2)]),labels=-1*pretty(par("yaxp")[c(1,2)]))
198009598544 Uploaded
romaingred
parents:
diff changeset
63 mtext("Scaled Read Depth",side=2,line=4.5,cex=1,font=2)
198009598544 Uploaded
romaingred
parents:
diff changeset
64 dev.off()
198009598544 Uploaded
romaingred
parents:
diff changeset
65 }
198009598544 Uploaded
romaingred
parents:
diff changeset
66 mapply( fun, fai\$V1, fai\$V2)`);
198009598544 Uploaded
romaingred
parents:
diff changeset
67 $R->stopR();
198009598544 Uploaded
romaingred
parents:
diff changeset
68 }
198009598544 Uploaded
romaingred
parents:
diff changeset
69
198009598544 Uploaded
romaingred
parents:
diff changeset
70 sub pie_chart
198009598544 Uploaded
romaingred
parents:
diff changeset
71 {
198009598544 Uploaded
romaingred
parents:
diff changeset
72 my $dir = shift;
198009598544 Uploaded
romaingred
parents:
diff changeset
73 my $in = $dir.'repartition.txt';
198009598544 Uploaded
romaingred
parents:
diff changeset
74 my $out = $dir.'pie_chart.png';
198009598544 Uploaded
romaingred
parents:
diff changeset
75
198009598544 Uploaded
romaingred
parents:
diff changeset
76 my $R = Statistics::R->new();
198009598544 Uploaded
romaingred
parents:
diff changeset
77 $R->startR;
198009598544 Uploaded
romaingred
parents:
diff changeset
78 $R->send(
198009598544 Uploaded
romaingred
parents:
diff changeset
79 qq`
198009598544 Uploaded
romaingred
parents:
diff changeset
80 library(plotrix)
198009598544 Uploaded
romaingred
parents:
diff changeset
81 library(RColorBrewer)
198009598544 Uploaded
romaingred
parents:
diff changeset
82 R =read.table("$in",header=T)
198009598544 Uploaded
romaingred
parents:
diff changeset
83 values = round(R\$percentage)
198009598544 Uploaded
romaingred
parents:
diff changeset
84 keys = R\$type
198009598544 Uploaded
romaingred
parents:
diff changeset
85 lab = paste(values, "%", sep="")
198009598544 Uploaded
romaingred
parents:
diff changeset
86 png("$out")
198009598544 Uploaded
romaingred
parents:
diff changeset
87 colors <- brewer.pal(7,"Paired")
198009598544 Uploaded
romaingred
parents:
diff changeset
88 pie(values, col=colors, labels=lab, clockwise=TRUE)
198009598544 Uploaded
romaingred
parents:
diff changeset
89 legend("bottom", legend = keys, fill=colors, bty="n", ncol = 3)
198009598544 Uploaded
romaingred
parents:
diff changeset
90 par(mai = c(0,0,0,0))
198009598544 Uploaded
romaingred
parents:
diff changeset
91 layout(c(1,2),heights=c(0.3,1))
198009598544 Uploaded
romaingred
parents:
diff changeset
92 plot.new()
198009598544 Uploaded
romaingred
parents:
diff changeset
93 legend("bottom", legend = keys, fill=colors, bty="n",ncol = 3)
198009598544 Uploaded
romaingred
parents:
diff changeset
94 pie(values, col=colors, labels=lab, clockwise=TRUE)
198009598544 Uploaded
romaingred
parents:
diff changeset
95 dev.off()`
198009598544 Uploaded
romaingred
parents:
diff changeset
96 );
198009598544 Uploaded
romaingred
parents:
diff changeset
97 $R->stopR();
198009598544 Uploaded
romaingred
parents:
diff changeset
98 }
198009598544 Uploaded
romaingred
parents:
diff changeset
99
198009598544 Uploaded
romaingred
parents:
diff changeset
100 1;