perl/etc/justonce/sperm/getlm.r

   1 #!/usr/bin/env littler
   2
   3 # R --save < lm.R --args bamrsplot.tsv tsv/nrss1m.tsv read.lm uncover.lm base.lm
   4
   5 # gzcat bamrsplot8.tsv.gz|perl -lane 'print $_ unless $F[0] =~ /[XYM]/i' > bamrsplot8.auto.tsv
   6 # gzcat rss1m.tsv.gz|perl -lane 'print $_ unless $F[0] =~ /[XYM]/i' > rss1m.auto.tsv
   7 # mv base.lm read.lm uncover.lm t/
   8 names <- c("MDA-Blood","MALBAC-Blood","MDA-Sperm23","MDA-Sperm24","MDA-Sperm28","MALBAC-SpermS01","MALBAC-SpermS02","MALBAC-SpermS03")
   9 theLen <- length(names)
  10
  11 DATa <- read.table("bamrsplot8.auto.tsv",skip=1)
  12 DATb <- read.table("rss1m.auto.tsv",skip=3)
  13
  14 for (i in seq(3,1+theLen)) {
  15         for (j in seq(i+1,2+theLen)) {
  16                 ta <- cor.test(DATa[,i],DATa[,j],method='pearson')
  17                 tb <- cor.test(DATa[,i],DATa[,j],method='kendall')
  18
  19                 cat(as.character(names[i-2]), as.character(names[j-2]), as.character(ta$p.value), as.character(ta$estimate), as.character(tb$p.value), as.character(tb$estimate),"\n",file='read.lm',append=T,sep="\t")
  20         }
  21 }
  22
  23 for (i in seq(3,2+theLen)){
  24         var=sd(DATa[,i])/mean(DATa[,i])
  25         cat(as.character(names[i-2]), as.character(var),"\n",file='read.lm',append=T,sep="\t")
  26 }
  27
  28 for (i in seq(3+theLen,1+2*theLen)) {
  29         for (j in seq(i+1,2+2*theLen)) {
  30                 ta <- cor.test(DATb[,i],DATb[,j],method='pearson')
  31                 tb <- cor.test(DATb[,i],DATb[,j],method='kendall')
  32                 cat(as.character(names[i-2-theLen]), as.character(names[j-2-theLen]), as.character(ta$p.value), as.character(ta$estimate), as.character(tb$p.value), as.character(tb$estimate),"\n",file='uncover.lm',append=T,sep="\t")
  33         }
  34
  35 }
  36 for(i in seq(3+theLen,2+2*theLen)){
  37         var=sd(DATb[,i])/mean(DATb[,i])
  38         cat(as.character(names[i-2-theLen]), as.character(var),"\n",file='uncover.lm',append=T,sep="\t")
  39 }
  40
  41 for (i in seq(3+3*theLen,1+4*theLen)) {
  42         for (j in seq(i+1,2+4*theLen)) {
  43                 ta <- cor.test(DATb[,i],DATb[,j],method='pearson')
  44                 tb <- cor.test(DATb[,i],DATb[,j],method='kendall')
  45                 cat(as.character(names[i-2-3*theLen]), as.character(names[j-2-3*theLen]), as.character(ta$p.value), as.character(ta$estimate), as.character(tb$p.value), as.character(tb$estimate),"\n",file='base.lm',append=T,sep="\t")
  46         }
  47 }
  48
  49 for(i in seq(3+3*theLen,2+4*theLen)){
  50         var=sd(DATb[,i])/mean(DATb[,i])
  51         cat(as.character(names[i-2-3*theLen]),as.character(var),"\n",file='base.lm',append=T,sep="\t")
  52 }