perl/etc/justonce/sperm/kmercmp.r

   1 #!/usr/bin/env r
   2 if (is.null(argv) | length(argv)<2) {
   3   cat("Usage: kmerfreq.r blood_hist hist_file\n")
   4   q()
   5 }
   6 bloodfile <- argv[1]
   7 infile <- argv[2]
   8
   9 if(!exists("D__myfit"))
  10         source("myfit.R")
  11
  12 maxXrange <- 30
  13 testXrange <- 15
  14 #infile <- 'mlbacDonor.k25.lz4.hist'
  15 #avgcvg <- 6.85
  16
  17 #library('vcd')
  18
  19 blood=read.table(bloodfile,skip=8)
  20 #mda=read.table('mdaS23.k25.lz4.hist',skip=8)
  21 #malbac=read.table('S01.k25.lz4.hist',skip=8)
  22 indata=read.table(infile,skip=8)
  23
  24 getdata <- function(x,range) {
  25         thefreq <- x[,1]
  26         cnt <- as.numeric(x[,2])
  27         theratio <- x[,3]
  28         thedata <- cbind(thefreq,cnt)
  29         flag <- thedata[,1] <= range
  30         #RVAL <- cbind(theratio[flag],thefreq[flag])
  31         RVAL <- cbind(cnt[flag],thefreq[flag])
  32         RVAL
  33 }
  34
  35 usedata1 <- getdata(blood,testXrange)
  36 usedata2 <- getdata(indata,testXrange)
  37 #usedata3 <- getdata(malbac,testXrange)
  38
  39 #xx=cbind(usedata2[,1],usedata1[,1])
  40 xx=usedata2[,1];
  41 yy=usedata1[,1]
  42 #yy=cbind(usedata1[,1],usedata2[,1])
  43 #aa=chisq.test(usedata2[,1], p = usedata1[,1], rescale.p = TRUE)
  44 aa=chisq.test(cbind(xx,yy))
  45 #bb=chisq.test(yy)
  46
  47 bb<-wilcox.test(xx,yy,alternative='two',paired=TRUE,conf.int=TRUE)
  48
  49 print(usedata2[,1])
  50 print(usedata1[,1])
  51 #print(xx)
  52 print(c(sum(xx),sum(yy)))
  53 print(aa)
  54 print(bb)
  55 "
  56 scaleres <- myfit(usedata2, par = list(lambda = 3))
  57 scaleres$fitted <- as.vector(usedata1[, 1])
  58 print('===scaledmyfit===')
  59 print(scaleres)
  60 dbg1 <- function(x) {
  61         arg <- deparse(substitute(x))
  62         RVAL <- cbind(x$observed,x$fitted,(x$observed-x$fitted)/x$fitted)
  63         #RVAL <- round(RVAL,digits = 6)
  64         colnames(RVAL) <- c(paste0(arg,'.ob'),paste0(arg,'.fi'),'(O-E)/E')
  65         rownames(RVAL) <- 1+ 1:nrow(RVAL)
  66         RVAL
  67 }
  68 print(cbind(dbg1(scaleres)))
  69 scalesum <- summary(scaleres)
  70 "
  71 mm=pdf2cdf(xx,T,T)
  72 nn=pdf2cdf(yy,T,T)
  73
  74 ksres <- ks.test(xx,yy,alternative='t')
  75 print(ksres)
  76 ksres <- ks.test(mm,nn,alternative='t')
  77 print(ksres)
  78
  79 "
  80 ./kmercmp.r blood.fq.k25.gz.hist.fixed S1.fq.k25.gz.hist.fixed
  81 ./kmercmp.r blood.fq.k25.gz.hist.fixed S2.fq.k25.gz.hist.fixed
  82 ./kmercmp.r blood.fq.k25.gz.hist.fixed S3.fq.k25.gz.hist.fixed
  83 ./kmercmp.r blood.fq.k25.gz.hist.fixed S23.fq.k25.gz.hist.fixed
  84 ./kmercmp.r blood.fq.k25.gz.hist.fixed S24.fq.k25.gz.hist.fixed
  85 ./kmercmp.r blood.fq.k25.gz.hist.fixed S28.fq.k25.gz.hist.fixed
  86
  87 ./kmercmp.r S2.fq.k25.gz.hist.fixed S1.fq.k25.gz.hist.fixed
  88 ./kmercmp.r blood.fq.k25.gz.hist.fixed mlbacDonor.k25.lz4.hist.fixed
  89 ./kmercmp.r mlbacDonor.k25.lz4.hist.fixed S23.fq.k25.gz.hist.fixed
  90
  91 "