Mercurial > repos > brasset_jensen > srnapipe
diff bin/Rcall.pm @ 40:4bc00caa60b4 draft
Uploaded
author | pierre.pouchin |
---|---|
date | Thu, 12 Jul 2018 08:27:08 -0400 |
parents | 6d9f127da28f |
children | 9185ca0a7b43 |
line wrap: on
line diff
--- a/bin/Rcall.pm Thu Jul 12 08:24:21 2018 -0400 +++ b/bin/Rcall.pm Thu Jul 12 08:27:08 2018 -0400 @@ -1,136 +1,136 @@ -package Rcall; - -use strict; -use warnings; -use Statistics::R; - -use Exporter; -our @ISA = qw(Exporter); -our @EXPORT_OK = qw( &histogram &pie_chart &bg_to_png ); - -sub histogram -{ - my ($size_hashR, $out_png, $size) = @_; - my (@abs, @ord); - my $i = 0; - foreach my $k (sort {$a <=> $b} keys %{$size_hashR}) - { - my $percentage = 0; - $percentage = $size_hashR->{$k} * 100 / $size if $size != 0; - $abs[$i] = $k ; $ord[$i] = $percentage; $i++; - } - my $abs = join (",", @abs ); - my $ord = join (",", @ord ); - if (scalar(@abs) != 0) - { - - my $R = Statistics::R->new(); - $R->startR; - $R->send( - qq`library(ggplot2) - percentage = c($ord) - size =c($abs) - min = min(size) - max = max(size) - dat = data.frame(size,percentage) - png(filename=\"$out_png\", width = 640, height = 640) - c = ggplot(dat,aes(size,percentage)) - 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") ) - dev.off()`); - $R->stopR(); - - } -} - -sub bg_to_png -{ - my ( $fai, $bgP, $bgM, $dir, $sb ) = @_; - my $R = Statistics::R->new(); - $R->startR; - $R->send( - qq`library('Sushi') - fai =read.table("$fai") - if ( file.info("$bgP")\$size !=0 ) - { - bgP = read.table("$bgP") - } else { bgP = data.frame(factor(),integer()) } - - if ( file.info("$bgM")\$size !=0 ) - { - bgM = read.table("$bgM") - } else { bgM = data.frame(factor(),integer()) } - - f_both = function(chr,end) { - jpeg( paste0("$dir",as.character(chr),".png"), quality=100) - par(mfrow=c(2,1),mar=c(1,10,1,3)) - plotBedgraph(bgP, chrom=chr,chromstart=0,chromend=end,transparency=.50, color=SushiColors(2)(2)[1]) - axis(side=2,las=2,tcl=.2) - mtext("Scaled Read Depth",side=2,line=4,cex=1,font=2) - plotBedgraph(bgM, chrom=chr,chromstart=0,chromend=end,transparency=.50, flip=TRUE, color=SushiColors(2)(2)[2]) - 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) - axis(side=2,las=2,tcl=.2,at=pretty(par("yaxp")[c(1,2)]),labels=-1*pretty(par("yaxp")[c(1,2)])) - mtext("Scaled Read Depth",side=2,line=4.5,cex=1,font=2) - dev.off() - } - - f_plus = function(chr,end) { - jpeg( paste0("$dir",as.character(chr),".png"), quality=100) - plotBedgraph(bgP, chrom=chr,chromstart=0,chromend=end,transparency=.50, color=SushiColors(2)(2)[1]) - 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) - axis(side=2,las=2,tcl=.2) - mtext("Scaled Read Depth",side=2,line=4,cex=1,font=2) - dev.off() - } - - f_minus = function(chr,end) { - jpeg( paste0("$dir",as.character(chr),".png"), quality=100) - plotBedgraph(bgM, chrom=chr,chromstart=0,chromend=end,transparency=.50, flip=TRUE, color=SushiColors(2)(2)[2]) - 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) - axis(side=2,las=2,tcl=.2,at=pretty(par("yaxp")[c(1,2)]),labels=-1*pretty(par("yaxp")[c(1,2)])) - mtext("Scaled Read Depth",side=2,line=4.5,cex=1,font=2) - dev.off() - } - - fai_b = fai[fai\$V1 %in% intersect(bgM\$V1,bgP\$V1), ] - mapply( f_both, fai_b\$V1, fai_b\$V2) - - fai_p = fai[fai\$V1 %in% setdiff(bgP\$V1,bgM\$V1), ] - mapply( f_plus, fai_p\$V1, fai_p\$V2) - - fai_m = fai[fai\$V1 %in% setdiff(bgM\$V1,bgP\$V1), ] - mapply( f_minus, fai_m\$V1, fai_m\$V2) `); - - $R->stopR(); -} - -sub pie_chart -{ - my $dir = shift; - my $in = $dir.'repartition.txt'; - my $out = $dir.'pie_chart.png'; - - my $R = Statistics::R->new(); - $R->startR; - $R->send( - qq` - library(plotrix) - library(RColorBrewer) - R =read.table("$in",header=T) - values = round(R\$percentage) - keys = R\$type - lab = paste(values, "%", sep="") - png("$out") - colors <- brewer.pal(7,"Paired") - pie(values, col=colors, labels=lab, clockwise=TRUE) - legend("bottom", legend = keys, fill=colors, bty="n", ncol = 3) - par(mai = c(0,0,0,0)) - layout(c(1,2),heights=c(0.3,1)) - plot.new() - legend("bottom", legend = keys, fill=colors, bty="n",ncol = 3) - pie(values, col=colors, labels=lab, clockwise=TRUE) - dev.off()` - ); - $R->stopR(); -} - -1; +package Rcall; + +use strict; +use warnings; +use Statistics::R; + +use Exporter; +our @ISA = qw(Exporter); +our @EXPORT_OK = qw( &histogram &pie_chart &bg_to_png ); + +sub histogram +{ + my ($size_hashR, $out_png, $size) = @_; + my (@abs, @ord); + my $i = 0; + foreach my $k (sort {$a <=> $b} keys %{$size_hashR}) + { + my $percentage = 0; + $percentage = $size_hashR->{$k} * 100 / $size if $size != 0; + $abs[$i] = $k ; $ord[$i] = $percentage; $i++; + } + my $abs = join (",", @abs ); + my $ord = join (",", @ord ); + if (scalar(@abs) != 0) + { + + my $R = Statistics::R->new(); + $R->startR; + $R->send( + qq`library(ggplot2) + percentage = c($ord) + size =c($abs) + min = min(size) + max = max(size) + dat = data.frame(size,percentage) + png(filename=\"$out_png\", width = 640, height = 640) + c = ggplot(dat,aes(size,percentage)) + 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") ) + dev.off()`); + $R->stopR(); + + } +} + +sub bg_to_png +{ + my ( $fai, $bgP, $bgM, $dir, $sb ) = @_; + my $R = Statistics::R->new(); + $R->startR; + $R->send( + qq`library('Sushi') + fai =read.table("$fai") + if ( file.info("$bgP")\$size !=0 ) + { + bgP = read.table("$bgP") + } else { bgP = data.frame(factor(),integer()) } + + if ( file.info("$bgM")\$size !=0 ) + { + bgM = read.table("$bgM") + } else { bgM = data.frame(factor(),integer()) } + + f_both = function(chr,end) { + jpeg( paste0("$dir",as.character(chr),".png"), quality=100) + par(mfrow=c(2,1),mar=c(1,10,1,3)) + plotBedgraph(bgP, chrom=chr,chromstart=0,chromend=end,transparency=.50, color=SushiColors(2)(2)[1]) + axis(side=2,las=2,tcl=.2) + mtext("Scaled Read Depth",side=2,line=4,cex=1,font=2) + plotBedgraph(bgM, chrom=chr,chromstart=0,chromend=end,transparency=.50, flip=TRUE, color=SushiColors(2)(2)[2]) + 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) + axis(side=2,las=2,tcl=.2,at=pretty(par("yaxp")[c(1,2)]),labels=-1*pretty(par("yaxp")[c(1,2)])) + mtext("Scaled Read Depth",side=2,line=4.5,cex=1,font=2) + dev.off() + } + + f_plus = function(chr,end) { + jpeg( paste0("$dir",as.character(chr),".png"), quality=100) + plotBedgraph(bgP, chrom=chr,chromstart=0,chromend=end,transparency=.50, color=SushiColors(2)(2)[1]) + 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) + axis(side=2,las=2,tcl=.2) + mtext("Scaled Read Depth",side=2,line=4,cex=1,font=2) + dev.off() + } + + f_minus = function(chr,end) { + jpeg( paste0("$dir",as.character(chr),".png"), quality=100) + plotBedgraph(bgM, chrom=chr,chromstart=0,chromend=end,transparency=.50, flip=TRUE, color=SushiColors(2)(2)[2]) + 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) + axis(side=2,las=2,tcl=.2,at=pretty(par("yaxp")[c(1,2)]),labels=-1*pretty(par("yaxp")[c(1,2)])) + mtext("Scaled Read Depth",side=2,line=4.5,cex=1,font=2) + dev.off() + } + + fai_b = fai[fai\$V1 %in% intersect(bgM\$V1,bgP\$V1), ] + mapply( f_both, fai_b\$V1, fai_b\$V2) + + fai_p = fai[fai\$V1 %in% setdiff(bgP\$V1,bgM\$V1), ] + mapply( f_plus, fai_p\$V1, fai_p\$V2) + + fai_m = fai[fai\$V1 %in% setdiff(bgM\$V1,bgP\$V1), ] + mapply( f_minus, fai_m\$V1, fai_m\$V2) `); + + $R->stopR(); +} + +sub pie_chart +{ + my $dir = shift; + my $in = $dir.'repartition.txt'; + my $out = $dir.'pie_chart.png'; + + my $R = Statistics::R->new(); + $R->startR; + $R->send( + qq` + library(plotrix) + library(RColorBrewer) + R =read.table("$in",header=T) + values = round(R\$percentage) + keys = R\$type + lab = paste(values, "%", sep="") + png("$out") + colors <- brewer.pal(7,"Paired") + pie(values, col=colors, labels=lab, clockwise=TRUE) + legend("bottom", legend = keys, fill=colors, bty="n", ncol = 3) + par(mai = c(0,0,0,0)) + layout(c(1,2),heights=c(0.3,1)) + plot.new() + legend("bottom", legend = keys, fill=colors, bty="n",ncol = 3) + pie(values, col=colors, labels=lab, clockwise=TRUE) + dev.off()` + ); + $R->stopR(); +} + +1;