#ct2015 - part 2. Data analysis from JMP genomics
library(reshape2)
library(ggplot2)
rm(list=ls())
setwd("~/Library/Mobile Documents/com~apple~CloudDocs/cistranstest")
options(warn=-1)
ct<-read.csv("ct2015_jmpg_pipe_cistrans.csv")
des<-read.csv("ct2015_design4.csv")
kdml<-read.csv("ct2015_counts4_kml.csv")
kd<-data.frame(t(kdml))
des$generation<-as.character(des$Species)
des$generation[grep("f1",des$generation)]<-"F1"
des$generation[grep("f0",des$generation)]<-"F0"
colnames(kd)<-kdml[,1]
kd<-kd[-1,]
kd<-cbind(des,kd)
# annot<-read.csv("ph2015_v1.1annot.edited.csv")
# annot.cull<-annot[annot$id %in% genes.parents.f1,c("chr","start","id")]
# colnames(ct)
# head(des)
#make a dataframe with the effects and direction of each term.
#figure out what the effects mean
pvals<-ct[,c("gene",colnames(ct)[grep("PrF", colnames(ct))])]
#pdf("ct2015_pvaluesummary.pdf")
par(mfrow=c(2,2))
hist(ct$Chi_Sq_Div_DF, main="chi2 model fit")
top.ps<-data.frame()
for(i in 2:8) {
hist(pvals[,i], main=colnames(pvals)[i])
best<-pvals[pvals[,i]==min(pvals[,i]),]
top.ps<-rbind(top.ps,best)
}
par(mfrow=c(1,1))
for(i in 2:8) {
best.genes<-as.character(pvals$gene[order(pvals[,i])][1:4])
stats<-ct[ct$gene %in% best.genes,grep("PrF", colnames(ct))]
test<-kd[,c("Treatment","Allele","generation",best.genes)]
test$unique<-paste(test$Treatment,test$Allele,sep="_")
test<-melt(test,id.vars=c("Treatment","Allele","generation","unique"))
colnames(test)[5:6]<-c("gene","log.knmm.counts")
ps<-paste(stats[,i-1], collapse = ' .. ')
print(
ggplot(test, aes(x=generation, y=as.numeric(log.knmm.counts), col=unique, group=unique, size=unique))+
theme_bw()+
stat_summary(fun.y = mean,
fun.ymin = function(x) mean(x) - sd(x),
fun.ymax = function(x) mean(x) + sd(x),
geom = "pointrange",
position=position_dodge(0.1)) +
stat_summary(fun.y = mean,
geom = "line",
position=position_dodge(0.1))+
scale_size_manual(values=c(1,1.5,1,1.5))+
scale_color_manual(values=c("pink","darkred","lightblue","darkblue"))+
ggtitle(paste("top", colnames(pvals)[i], "\n",ps))+
facet_wrap(~gene, scale="free")
)
}
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
## ymax not defined: adjusting position using y instead
#dev.off()
#qvalue transformations for pvalues
library(qvalue)
##
## Attaching package: 'qvalue'
##
## The following object is masked from 'package:ggplot2':
##
## qplot
pqvals<-pvals
sig.genes<-list()
sig.out<-matrix(nrow=dim(pvals)[1], ncol=dim(pvals)[2]-1)
count<-0
for ( i in colnames(pvals)[grep("PrF",colnames(pvals))]){
count<-count+1
dat<-pvals[,i]
q.out<-data.frame(qvalue(dat,robust=T)$qvalue)
colnames(q.out)<-gsub("PrF_T3_","q_",i)
sig<-data.frame(ifelse(q.out[,1]<0.05,"sig","not.sig"))
colnames(sig)<-gsub("PrF_T3_","sig_",i)
pqvals<-cbind(pqvals,q.out,sig)
print(gsub("PrF_T3_","q_",i))
print(table(sig))
sig.genes[[gsub("PrF_T3_","sig_",i)]]<-as.character(pvals$gene[sig[,1]=="sig"])
sig.out[,count]<-q.out[,1]<0.05
}
## [1] "q_Treatment"
## sig
## not.sig sig
## 3295 13183
## [1] "q_Cis"
## sig
## not.sig sig
## 5835 10643
## [1] "q_Trans_fil"
## sig
## not.sig sig
## 11995 4483
## [1] "q_Trans_hal"
## sig
## not.sig sig
## 12222 4256
## [1] "q_Treatment_Cis"
## sig
## not.sig sig
## 15852 626
## [1] "q_Treatment_Trans_fil"
## sig
## not.sig sig
## 16163 315
## [1] "q_Treatment_Trans_hal"
## sig
## not.sig sig
## 16407 71
length(intersect(sig.genes[["sig_Treatment"]],sig.genes[["sig_Cis"]]))
## [1] 8470
head(pqvals)
## gene PrF_T3_Treatment PrF_T3_Cis PrF_T3_Trans_fil
## 1 Pahalv11b000009m.g 1.3e-02 2.4e-01 0.630
## 2 Pahalv11b000011m.g 4.7e-08 7.1e-08 0.810
## 3 Pahalv11b000012m.g 9.0e-21 4.3e-07 0.066
## 4 Pahalv11b000014m.g 5.5e-02 3.0e-11 0.490
## 5 Pahalv11b000017m.g 2.9e-04 1.5e-01 0.920
## 6 Pahalv11b000020m.g 4.0e-12 8.5e-07 0.190
## PrF_T3_Trans_hal PrF_T3_Treatment_Cis PrF_T3_Treatment_Trans_fil
## 1 0.5600 0.0083 0.002
## 2 0.3600 0.1200 0.270
## 3 0.0021 0.3500 0.320
## 4 0.2200 0.1500 0.420
## 5 0.6100 0.4200 0.570
## 6 0.0043 0.2100 0.990
## PrF_T3_Treatment_Trans_hal q_Treatment sig_Treatment q_Cis
## 1 0.77 0.0041790245 sig 1.040199e-01
## 2 0.20 0.0000365399 sig 7.318932e-05
## 3 0.45 0.0000365399 sig 7.318932e-05
## 4 0.12 0.0156588971 sig 7.318932e-05
## 5 0.12 0.0001274116 sig 6.982220e-02
## 6 0.89 0.0000365399 sig 7.318932e-05
## sig_Cis q_Trans_fil sig_Trans_fil q_Trans_hal sig_Trans_hal
## 1 not.sig 0.38872695 not.sig 0.369384339 not.sig
## 2 sig 0.44848033 not.sig 0.282781974 not.sig
## 3 sig 0.09100879 not.sig 0.009116491 sig
## 4 sig 0.33545158 not.sig 0.206323886 not.sig
## 5 not.sig 0.47850837 not.sig 0.388455716 not.sig
## 6 sig 0.18190294 not.sig 0.014834021 sig
## q_Treatment_Cis sig_Treatment_Cis q_Treatment_Trans_fil
## 1 0.1291912 not.sig 0.06174627
## 2 0.5453309 not.sig 0.47393878
## 3 0.8249442 not.sig 0.50121276
## 4 0.5938799 not.sig 0.55476380
## 5 0.8669049 not.sig 0.61965920
## 6 0.6904278 not.sig 0.72790122
## sig_Treatment_Trans_fil q_Treatment_Trans_hal sig_Treatment_Trans_hal
## 1 not.sig 0.7543327 not.sig
## 2 not.sig 0.4974981 not.sig
## 3 not.sig 0.6499453 not.sig
## 4 not.sig 0.4163191 not.sig
## 5 not.sig 0.4163191 not.sig
## 6 not.sig 0.7810281 not.sig
library(limma)
par(mfrow=c(1,1))
colnames(sig.out)<- gsub("sig_","",names(sig.genes))
sig.out2<-data.frame(sig.out)
sig.out2$Trans<-apply(sig.out[,c("Trans_fil","Trans_hal")],1, function(x) ifelse(any(x),TRUE,FALSE))
sig.out2$Treatment_Trans<-apply(sig.out[,c("Treatment_Trans_fil","Treatment_Trans_hal")],1, function(x) ifelse(any(x),TRUE,FALSE))
sig.out2<-sig.out2[,-c(3,4,6,7)]
colnames(sig.out2)<-gsub("_","*",colnames(sig.out2))
a <- vennCounts(sig.out2)
vennDiagram(a,circle.col=rainbow(5))
