Mercurial > repos > romaingred > pirna_pipeline

comparison bin/Rcall.pm @ 11:ee6b4b2072a9 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author romaingred
date Fri, 13 Oct 2017 08:46:31 -0400
parents 198009598544
children fcaea740c23a
comparison
equal deleted inserted replaced
10:42bc59c7db3a 11:ee6b4b2072a9
47 my $R = Statistics::R->new(); 47 my $R = Statistics::R->new();
48 $R->startR; 48 $R->startR;
49 $R->send( 49 $R->send(
50 qq`library('Sushi') 50 qq`library('Sushi')
51 fai =read.table("$fai") 51 fai =read.table("$fai")
52 bgP = read.table("$bgP") 52 if ( file.info("$bgP")\$size !=0 )
53 bgM = read.table("$bgM") 53 {
54 fun = function(chr,end) { 54 bgP = read.table("$bgP")
55 } else { bgP = data.frame(factor(),integer()) }
56
57 if ( file.info("$bgM")\$size !=0 )
58 {
59 bgM = read.table("$bgM")
60 } else { bgM = data.frame(factor(),integer()) }
61
62 f_both = function(chr,end) {
55 jpeg( paste0("$dir",as.character(chr),".png"), quality=100) 63 jpeg( paste0("$dir",as.character(chr),".png"), quality=100)
56 par(mfrow=c(2,1),mar=c(1,10,1,3)) 64 par(mfrow=c(2,1),mar=c(1,10,1,3))
57 plotBedgraph(bgP, chrom=chr,chromstart=0,chromend=end,transparency=.50, color=SushiColors(2)(2)[1]) 65 plotBedgraph(bgP, chrom=chr,chromstart=0,chromend=end,transparency=.50, color=SushiColors(2)(2)[1])
58 axis(side=2,las=2,tcl=.2) 66 axis(side=2,las=2,tcl=.2)
59 mtext("Scaled Read Depth",side=2,line=4,cex=1,font=2) 67 mtext("Scaled Read Depth",side=2,line=4,cex=1,font=2)
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) 69 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)
62 axis(side=2,las=2,tcl=.2,at=pretty(par("yaxp")[c(1,2)]),labels=-1*pretty(par("yaxp")[c(1,2)])) 70 axis(side=2,las=2,tcl=.2,at=pretty(par("yaxp")[c(1,2)]),labels=-1*pretty(par("yaxp")[c(1,2)]))
63 mtext("Scaled Read Depth",side=2,line=4.5,cex=1,font=2) 71 mtext("Scaled Read Depth",side=2,line=4.5,cex=1,font=2)
64 dev.off() 72 dev.off()
65 } 73 }
66 mapply( fun, fai\$V1, fai\$V2)`); 74
75 f_plus = function(chr,end) {
76 jpeg( paste0("$dir",as.character(chr),".png"), quality=100)
77 plotBedgraph(bgP, chrom=chr,chromstart=0,chromend=end,transparency=.50, color=SushiColors(2)(2)[1])
78 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)
79 axis(side=2,las=2,tcl=.2)
80 mtext("Scaled Read Depth",side=2,line=4,cex=1,font=2)
81 dev.off()
82 }
83
84 f_minus = function(chr,end) {
85 jpeg( paste0("$dir",as.character(chr),".png"), quality=100)
86 plotBedgraph(bgM, chrom=chr,chromstart=0,chromend=end,transparency=.50, flip=TRUE, color=SushiColors(2)(2)[2])
87 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)
88 axis(side=2,las=2,tcl=.2,at=pretty(par("yaxp")[c(1,2)]),labels=-1*pretty(par("yaxp")[c(1,2)]))
89 mtext("Scaled Read Depth",side=2,line=4.5,cex=1,font=2)
90 dev.off()
91 }
92
93 fai_b = fai[fai\$V1 %in% intersect(bgM\$V1,bgP\$V1), ]
94 mapply( f_both, fai_b\$V1, fai_b\$V2)
95
96 fai_p = fai[fai\$V1 %in% setdiff(bgP\$V1,bgM\$V1), ]
97 mapply( f_plus, fai_p\$V1, fai_p\$V2)
98
99 fai_m = fai[fai\$V1 %in% setdiff(bgM\$V1,bgP\$V1), ]
100 mapply( f_minus, fai_m\$V1, fai_m\$V2) `);
101
67 $R->stopR(); 102 $R->stopR();
68 } 103 }
69 104
70 sub pie_chart 105 sub pie_chart
71 { 106 {