library(knitr)
knitr::opts_chunk$set(echo = T)
library(mixtools)

## mixtools package, version 1.0.4, Released 2016-01-11
## This package is based upon work supported by the National Science Foundation under Grant No. SES-0518772.

library(reshape2)
library(ggplot2)
library(VennDiagram)

## Loading required package: grid

## 
## Attaching package: 'grid'

## The following object is masked from 'package:mixtools':
## 
##     depth

## Loading required package: futile.logger

## 
## Attaching package: 'VennDiagram'

## The following object is masked from 'package:mixtools':
## 
##     ellipse

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(sva)

## Loading required package: mgcv

## Loading required package: nlme

## 
## Attaching package: 'nlme'

## The following object is masked from 'package:dplyr':
## 
##     collapse

## This is mgcv 1.8-11. For overview type 'help("mgcv-package")'.

## 
## Attaching package: 'mgcv'

## The following object is masked from 'package:futile.logger':
## 
##     scat

## Loading required package: genefilter

## 
## Attaching package: 'genefilter'

## The following object is masked from 'package:base':
## 
##     anyNA

library(Biobase)

## Loading required package: BiocGenerics

## Loading required package: parallel

## 
## Attaching package: 'BiocGenerics'

## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB

## The following objects are masked from 'package:dplyr':
## 
##     combine, intersect, setdiff, union

## The following objects are masked from 'package:stats':
## 
##     IQR, mad, xtabs

## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, as.vector, cbind,
##     colnames, do.call, duplicated, eval, evalq, Filter, Find, get,
##     grep, grepl, intersect, is.unsorted, lapply, lengths, Map,
##     mapply, match, mget, order, paste, pmax, pmax.int, pmin,
##     pmin.int, Position, rank, rbind, Reduce, rownames, sapply,
##     setdiff, sort, table, tapply, union, unique, unlist, unsplit

## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.

library(grid)

mycolors = c("darkorange" ,"dodgerblue","limegreen","navy","hotpink",    "plum4","blue" , "magenta2","mediumpurple"  ,  "firebrick" ,"darkolivegreen4","royalblue3")
palette(mycolors)
source("081916_funcs_erica.R")
opar <- par(pch=19,mfrow=c(1,1),mar=c(2.2,2,1.4,1),cex=1)

#source("aug15_working_chip_processing.R")

SC media

post background tag removal, log ratio summarys

functions

myLRsig: keeps only significant log ratios if they are greater than threshold (> 1) in all three replicates
mylogLM_sum: call significance on log ratio matrix using linear model

bhymix= as.matrix(read.delim(
"aug18_2016_bhymix.txt",header = T,stringsAsFactors = F,check.names = F))


boymix= as.matrix(read.delim(
"aug18_2016_boymix.txt",header = T,stringsAsFactors = F,check.names = F))

pypd = read.delim(
"aug19_2016_allpools_pYPD.txt",
header = T,
stringsAsFactors = F,
check.names = F
)

purpose: correct for level of inhibition in fitness assays – YPD

put this in text

For the drug conditions an additional normalization step was required to separate condition-specific from growth-dependent fitness defects. Tag intensities that correlated with the trend line of growth inhibition were shifted to remove these growth-rate dependent effects.

steps

sort the pdata file in order of increasing growth inhibition fit a loess model to strains with high degree of correlation

phy = filter(pypd,pool == "hom")
poy = filter(pypd,pool == "oliver")

phy = phy %>% arrange(type,gen,cond)
poy = poy %>% arrange(type,gen,cond)

bhymix = bhymix[,phy$name]
boymix = boymix[,poy$name]

these tags are annoying – they are in barseq data either –

clc = grep('CHC1',rownames(bhymix))
if(length(clc > 0)) bhymix = bhymix[-clc,]

clc = grep('CHC1',rownames(boymix))
if(length(clc > 0)) boymix = boymix[-clc,]

clc = grep('CLC1',rownames(bhymix))
if(length(clc > 0)) bhymix = bhymix[-clc,]

clc = grep('CLC1',rownames(boymix))
if(length(clc > 0)) boymix = boymix[-clc,]

checkout a couple different correlation thresholds – strains that correlate with inhibition level –

opar <- par(pch=19,mfrow=c(1,2),mar=c(2.2,2,1.4,1),cex=1.2)

why = which(phy$type == 'ctrl')
woy = which(poy$type == 'ctrl')
plot(phy$gen,col=color(phy$cond))
plot(poy$gen,col=color(poy$cond))

# relative inhibtion = 1 - (doubling time of ctrl / by doubling time of treatment)

phy$inb = 1 - (mean(phy$gen[why])/phy$gen)

xcor<- cor(phy$gen, t(bhymix), use="pairwise.complete.obs", method="s")

### checking to see what threshold works best -- 

hadj5 = myinhibCts(bhymix,ctsVal = phy$gen,corThresh=0.5,nCores = 2)

## Loading required package: multicore

## WARNING: multicore has been superseded and will be removed shortly

## 
## Attaching package: 'multicore'

## The following object is masked from 'package:dplyr':
## 
##     collect

## Loading required package: limma

## 
## Attaching package: 'limma'

## The following object is masked from 'package:BiocGenerics':
## 
##     plotMA

oadj5 = myinhibCts(boymix,ctsVal = poy$gen,corThresh=0.5,nCores = 2)

hadj6 = myinhibCts(bhymix,ctsVal = phy$gen,corThresh=0.6,nCores = 2)
oadj6 = myinhibCts(boymix,ctsVal = poy$gen,corThresh=0.6,nCores = 2)


lhadj5= getLR(hadj5,why)
lhadj6 = getLR(hadj6,why)
loadj5  = getLR(oadj5 ,woy)
loadj6  = getLR(oadj6 ,woy)

lbhymix = getLR(bhymix,why)
lboymix = getLR(boymix,woy)

check how well correction worked; focus on log ratio < -1

remember its only going to correct for strains that correlate with #levels of inhibition – not for otherwise flaky tags

and will only be as good as the our growth rate estimates

left panel NO correction, middle cor = 0.5, right panel 0.6

YPG works best –

opar <- par(pch=19,mfrow=c(1,3),mar=c(2.2,2,1.4,1),cex=1.2)


for (i in 1:17){
  psig(lboymix,loadj5,i)
  p1(loadj6,i)
  
  
}

for (i in 1:15){
psig(lbhymix,lhadj5,i)
p1(lhadj6,i)
}

conclusion: 0.5 threshold seems to work best for oliver and hom

anna used +/- correlation threshold, but for our small dataset this is too extreme;

removes too many true positive log ratios:

otst = medShift.strainCts(boymix,ctsVal = poy$gen,corThresh=0.5,nCores = 2)

htst = medShift.strainCts(bhymix,ctsVal = phy$gen,corThresh=0.5,nCores = 2)

lotst = getLR(otst,woy)

lhtst = getLR(htst,why)

opar <- par(pch=19,mfrow=c(1,2),mar=c(2.2,2,1.4,1),cex=1)


for (i in c(1,4,8,12,15)) psig(lotst,loadj5,i)

for (i in c(1,4,8,12,15)) psig(lhtst,lhadj5,i)

inhibition_correction mito drugs

ggiaever

08/19/2016