chip_barseq.Rmd

  fdata6 = read.delim("jun10_2016_fdata6.txt",
    stringsAsFactors = F,
    check.names = F)
    
  noness = filter(fdata6, fdata6$essential == 0)
  noness = noness %>% arrange(strain)

preprocessing for low cnts – 1. count matrix 2. filter count matrix for ALL rows < 50 3. filter count matrix for all CTRLS < 50

example, SC dropouts hom samples

xbar= as.matrix(read.delim('aug6_2016_hom_barseq.txt',header = T,stringsAsFactors =F,check.names = F,strip.white = T))

p11 = read.delim("oct2_phsbar.txt",header = T,stringsAsFactors = F,check.names = F)

w11 = which(p11$type == 'ctrl')
lp11 = p11[-w11,]
#filter out essential strains
wnebar = which(noness$strain %in% rownames(xbar))

cat(paste0(' # of rows before filtering: ',nrow(xbar)))

hsbar = myall_less50(xbar[noness$strain[wnebar],p11$name])

cat(paste0(' # of rows after filtering for ALL cnts < 50: ',nrow(hsbar)))

hsbar = mymin50(hsbar,w11)

cat(paste0(' # of rows after filtering cnts less 50 in ANY control: ',nrow(hsbar)))

retrieve normalized counts from edgeR and DESeq2, and normalizing by count w/out package

group = model factor

ref = reference condition

hsbar = hsbar[,p11$name]

p11$cond = factor(p11$cond)
p11$cond = relevel(p11$cond,ref='sc')
hdds = mydds(hsbar,groups = p11$cond)

hedge = myedgeR(hsbar,group = p11$cond,ref = 'sc')

mysumtags = sum counts from up and down tags

mysumcond = sum counts from each condition

post processing normalized count matrix collapse matrix dimensions by:

1. sum up and downtags
2. summing replicate condtions

hdds2 = myproc_normcounts(hdds,p11$cond)

hedge2 = myproc_normcounts(hedge,p11$cond)

compare this to same data generated by arrays

x11= as.matrix(read.delim('x11.txt',header = T,stringsAsFactors =F,check.names = F,strip.white = T))

x111 = mymixEMctrl(x11,w11)
btst = myBesttag(x111,bgThresh = 6.3)
hbat = ComBat(btst,mod = model.matrix(~cond,p11),batch = p11$pcr)
tsw = mysweep(hbat)

hsum = mysumm(tsw[,-w11],lp11$cond)

comparison to normalizing by counts/experiment without an R package

tst = mybarnorm(hsbar)

tst2 = myproc_normcounts(tst,p11$cond)

noRpkg = tst2[,colnames(hsum)]

plot comparisions

YPD oliver and hom

get barseq ypd data

xobar= as.matrix(read.delim('count_matrix_aug4_noctrl_celnames.txt',
header = T,stringsAsFactors = F,sep = "\t",row.names = 1,check.names = F,as.is=TRUE))

perica = read.delim("jun10_2016_pfinal_erica.txt",header = T,stringsAsFactors=F, check.names = F)

w= which(perica$celfile %in% colnames(xobar))

pbar = perica[w,]

xobar = xobar[which(rownames(xobar) %in% noness$strain),pbar$celfile]

colnames(xobar)=pbar$name

#pjen = data.frame(name= colnames(xobar)[grep('Jen',colnames(xobar))])

phybar = pbar %>% filter(pool == 'hom',media == 'ypd')
poybar = pbar %>% filter(pool == 'oliver',media == 'ypd')
phybar = phybar %>% arrange(type,cond,pcr)
poybar = poybar %>% arrange(type,cond,pcr)

woybar = which(poybar$type == 'ctrl')
whybar = which(phybar$type == 'ctrl')

filter as above for low counts

dds and edge evaluation for ypd – check out the FEB17 0083 outlier for hom ypd

hydds = mydds(hybar,phybar$cond)
oydds = mydds(oybar,poybar$cond)
hyedge = myedgeR(hybar,phybar$cond,ref = 'ypd')
oyedge = myedgeR(oybar,poybar$cond,ref = 'ypd')

hyd = myproc_normcounts(hydds,phybar$cond)
oyd = myproc_normcounts(oydds,poybar$cond)

hyj = myproc_normcounts(hyedge,phybar$cond)
oyj = myproc_normcounts(oyedge,poybar$cond)

ypd array datag

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

omix = xochip[,poybar$name]
omix = mymixEMctrl(omix,woybar)
omix = medShift.updn(omix)
omin = min(apply(omix[,woybar],1,median))

hmix = xochip[,phybar$name]
hmix = mymixEMctrl(hmix,whybar)
hmin = min(apply(hmix[,whybar],1,median))

bomix = myBesttag(omix,bgThresh=omin)

bhmix = myBesttag(hmix,bgThresh=hmin)

#before batch correction
mydendy(bhmix,phybar$pcr)

bhybat = ComBat(bhmix,mod=model.matrix(~cond,phybar),batch=phybar$pcr)
#after batch correction
mydendy(bhybat,phybar$cond)

#oliver before batch correction
mydendy(bomix,poybar$cond)
boybat = ComBat(bomix,mod=model.matrix(~cond,poybar),batch=poybar$pcr)

#oliver after batch correction
mydendy(boybat,poybar$cond)

hysw= mysweep(bhybat)
lphybar = phybar[-whybar,]
#summarize chip data
hysum = mysumm(hysw[,lphybar$name],lphybar$cond)

oysw= mysweep(boybat)
lpoybar = poybar[-woybar,]
oysum = mysumm(oysw[,lpoybar$name],lpoybar$cond)

normalize YPD data by count

noRhy = mybarnorm(hybar)

noRhy3 = myproc_normcounts(noRhy,phybar$cond)

noRhy3 = noRhy3[,colnames(hysum)]

noRoy = mybarnorm(oybar)

noRoy3 = myproc_normcounts(noRoy,poybar$cond)

noRoy3 = noRoy3[,colnames(oysum),drop=F]

plot YPD results

hom YPD

opar <- par(pch=19,mfrow=c(1,2),mar=c(2.2,2,1.4,1),cex=1.2)
for (i in 1:ncol(hysum)) { 
  
  p1(hysum,i)
  p1(noRhy3,i)
  p1(hyj,i)
  p1(hyd,i)
  
}

oliver YPD

opar <- par(pch=19,mfrow=c(1,2),mar=c(2.2,2,1.4,1),cex=1.2)
for (i in 1:ncol(oysum)) { 
  
  p1(oysum,i)
  p1(noRoy3,i,sig =0.7)
  p1(oyj,i,sig =0.7)
  p1(oyd,i,sig =0.7)
  
}

compare to combining biological replicates with standard edge:

1. up and downtags separately, collapsed by summation after dispersion estimates

hs_std_edge = mydfedge(hsbar,p11$cond,ref= 'sc')
mhs_std = as.matrix(hs_std_edge[,1:5])
rownames(mhs_std) = hs_std_edge$strain

hs_sumtags = mysumtags(mybarnorm(hsbar))
hs_std_sumtags = mydfedge(hs_sumtags,p11$cond,ref= 'sc')
mhs_std_sumtags = as.matrix(hs_std_sumtags[,1:5])
rownames(mhs_std_sumtags) = hs_std_sumtags$strain
mhs_std_sumtags= mhs_std_sumtags[order(rownames(mhs_std_sumtags)),]



hy_std_edge = mydfedge(hybar,phybar$cond,ref= 'ypd')
mhy_std = as.matrix(hy_std_edge[,1:4])
rownames(mhy_std) = hy_std_edge$strain


hy_sumtags = mysumtags(mybarnorm(hybar))

hy_std_sumtags = mydfedge(hy_sumtags,phybar$cond,ref= 'ypd')
mhy_std_sumtags = as.matrix(hy_std_sumtags[,1:4])
rownames(mhy_std_sumtags) = hy_std_sumtags$strain
mhy_std_sumtags= mhy_std_sumtags[order(rownames(mhy_std_sumtags)),]

2. normalize up and downtags separately before combining, predispersion

plot standard edge, pre-summed up and downtags

opar <- par(pch=19,mfrow=c(1,2),mar=c(2.2,2,1.4,1),cex=1.2)
for (i in 1:3)
psig(-mhs_std_sumtags,hedge2,i)

for (i in 1:2)
psig(-mhy_std_sumtags,hyj,i)