perl/soap/draw/cvgstat.pl

   1 #!/bin/env perl
   2 use strict;
   3 use warnings;
   4
   5 unless (@ARGV>0){
   6         print "perl $0 depsingle\n";
   7         exit;
   8 }
   9 my ($depf)=@ARGV;
  10 my $statf=$depf;
  11 $statf =~ s/\.[^.\/\\]+$//;
  12 $statf .= '.depstat';
  13
  14 my ($ChrID,%Stat,%StatA)=('Empty');
  15 open DEPTH,'<',$depf or die "[x]Error opening $depf: $!\n";
  16 open STAT,'>',$statf or die "[x]Error opening $statf: $!\n";
  17 open STATA,'>',$statf.'all' or die "[x]Error opening ${statf}all: $!\n";
  18 while (<DEPTH>) {
  19         chomp;
  20         if (/^>/) {
  21                 if (exists $Stat{$ChrID}) {
  22                         print STAT "[Chr$ChrID]\n";     # for Reading in R, SectionName must be string
  23                         for my $v (sort {$a <=> $b} keys %{$Stat{$ChrID}}) {
  24                                 print STAT 'd',$v,"=",$Stat{$ChrID}{$v},"\n";   # for Reading in R, Key must be string
  25                                 #print $v,"=",$Stat{$ChrID}{$v},"\n";
  26                         }
  27                         print STAT "\n";
  28                         print STDERR "\b\b\bdone.\n";
  29                         delete $Stat{$ChrID};
  30                 }
  31                 print STDERR "$_ ...";
  32                 s/^>//;
  33                 $ChrID=$_;
  34                 next;
  35         } else {
  36                 for (split /\s+/) {
  37                         ++$Stat{$ChrID}{$_};
  38                         ++$StatA{$_};
  39                 }
  40         }
  41 }
  42 close DEPTH;
  43
  44                 # for the last Chr.
  45                 if (exists $Stat{$ChrID}) {
  46                         print STAT "[$ChrID]\n";
  47                         for my $v (sort {$a <=> $b} keys %{$Stat{$ChrID}}) {
  48                                 print STAT $v,"=",$Stat{$ChrID}{$v},"\n";
  49                         }
  50                         print STAT "\n";
  51                         print STDERR "\b\b\bo\n";
  52                         delete $Stat{$ChrID};
  53                 }
  54
  55 close STAT;
  56
  57 for my $v (sort {$a <=> $b} keys %StatA) {
  58         print STATA $v,"\t",$StatA{$v},"\n";
  59 }
  60
  61 close STATA;
  62
  63 __END__
  64 > x=seq(0,30,1)
  65 > plot(x,dpois(x,12),type='p')
  66 > lines(x,dpois(x,12))
  67 > lines(c(12,12),c(-1,1),col='blue')
  68
  69 HUMfcmRABDMAAPE
  70 Lib='HUMfcmRAADDAAPE'
  71 a=read.delim(paste('g',Lib,'.depstatall',sep=''),F)
  72 a[1,2]=a[1,2]-239850802
  73 EffLen=sum(as.numeric(a$V2))
  74 maxX=20
  75 maxY=0.15
  76 x=seq(0,maxX+1,1)
  77 c=sum(as.numeric(a$V1*a$V2))/EffLen
  78 pdf(paste('g',Lib,'.pdf',sep=''),width=8,height=6,paper="special")
  79 plot(a$V1,a$V2/EffLen,xlim=c(0,maxX),ylim=c(0,maxY),type='b',xlab='depth',ylab='density',main=paste('g',Lib,sep=''),
  80         sub=paste('depth=',round(c,4),', ','cvg=',round(100*(1-a[1,2]/EffLen),4),' % < ',round(100*(1-exp(-c)),4),' %',sep='')
  81         )
  82 #lines(a$V1,a$V2/EffLen)
  83 lines(x,dpois(x,c),col='blue')
  84 lines(c(c,c),c(-1,1),col='red')
  85 axis(1,at=c)
  86 dev.off()
  87
  88 https://stat.ethz.ch/pipermail/r-help/2007-June/134110.html
  89
  90 Since in my problem the structure of the INI sections is almost static and
  91 always present, I extended your example to create an in-memory list of
  92 everything in the INI file with this function:
  93
  94 # Prototype of how to read INI files to process olfactometer data
  95 # efg, 13 June 2007
  96 # Thanks to Gabor Grothendieck for helpful suggestions in the R-Help
  97 # mailing list on how to parse the INI file.
  98 Parse.INI <- function(INI.filename)
  99 {
 100   connection <- file(INI.filename)
 101   Lines  <- readLines(connection)
 102   close(connection)
 103
 104   Lines <- chartr("[]", "==", Lines)  # change section headers
 105
 106   connection <- textConnection(Lines)
 107   d <- read.table(connection, as.is = TRUE, sep = "=", fill = TRUE)
 108   close(connection)
 109
 110   L <- d$V1 == ""                    # location of section breaks
 111   d <- subset(transform(d, V3 = V2[which(L)[cumsum(L)]])[1:3],
 112                            V1 != "")
 113
 114   ToParse  <- paste("INI.list$", d$V3, "$",  d$V1, " <- '",
 115                     d$V2, "'", sep="")
 116
 117   INI.list <- list()
 118   eval(parse(text=ToParse))
 119
 120   return(INI.list)
 121 }
 122
 123
 124 Here's an example of using the above function (I'll put the sample input
 125 file below):
 126
 127 INI1 <- Parse.INI("sample.ini")
 128
 129 # Explore INI contents
 130 summary(INI1)
 131
 132 INI1$SystemSetup$OlfactometerCode
 133 INI1$DefaultLevels
 134 unlist(INI1$DefaultLevels)
 135 INI1$Map
 136
 137 INI1$Map$port1
 138 as.integer( unlist( strsplit(INI1$Map$port1, ",") ) )
 139
 140 = = = = =
 141 Sample output:
 142
 143 > INI1 <- Parse.INI("sample.ini")
 144 >
 145 > # Explore INI contents
 146 > summary(INI1)
 147               Length Class  Mode
 148 SystemSetup   1      -none- list
 149 Files         8      -none- list
 150 DefaultLevels 4      -none- list
 151 OdorNames     2      -none- list
 152 Map           3      -none- list
 153 >
 154 > INI1$SystemSetup$OlfactometerCode
 155 [1] "3"
 156 > INI1$DefaultLevels
 157 $FC00
 158 [1] "50"
 159
 160 $FC01
 161 [1] "100"
 162
 163 $FC02
 164 [1] "50"
 165
 166 $FC10
 167 [1] "50"
 168
 169 > unlist(INI1$DefaultLevels)
 170  FC00  FC01  FC02  FC10
 171  "50" "100"  "50"  "50"
 172 > INI1$Map
 173 $port0
 174 [1] "0,0,0,0,0,0,0,0,0,0,0,0"
 175
 176 $port1
 177 [1] "0,0,0,0,0,0,0,0,0,0,0,0"
 178
 179 $port2
 180 [1] "0,0,0,0,0,0,0,0,0,0,0,0"
 181
 182 >
 183 > INI1$Map$port1
 184 [1] "0,0,0,0,0,0,0,0,0,0,0,0"
 185 > as.integer( unlist( strsplit(INI1$Map$port1, ",") ) )
 186  [1] 0 0 0 0 0 0 0 0 0 0 0 0
 187
 188 = = = = =
 189 Sample input file, sample.ini:
 190
 191 [SystemSetup]
 192 OlfactometerCode=3
 193 [Files]
 194 prelog0=Part0.txt
 195 date0=2:06:27.461 PM 6/9/2007
 196 note0=group1-1
 197 name0=group1
 198 prelog1=Part1.txt
 199 date1=2:09:16.809 PM 6/9/2007
 200 note1=group1-1
 201 name1=group1-1
 202 [DefaultLevels]
 203 FC00=50
 204 FC01=100
 205 FC02=50
 206 FC10=50
 207 [OdorNames]
 208 port0=None
 209 port1=None
 210 [Map]
 211 port0=0,0,0,0,0,0,0,0,0,0,0,0
 212 port1=0,0,0,0,0,0,0,0,0,0,0,0
 213 port2=0,0,0,0,0,0,0,0,0,0,0,0
 214
 215 plot(as.matrix(INI1$Chr11),type='l')