Merge = read.table('./Figure3/Cross.tissue.comparison.haQTL.txt',head=T) 

Merge$logP.tis2<- -log10(Merge$nominalP_tis2) * Merge$slope_tis2/abs(Merge$slope_tis2)
Merge$logP.tis1<- -log10(Merge$nominalP_tis1) * Merge$slope_tis1/abs(Merge$slope_tis1)
Merge<-Merge[order(Merge$ifQTL_tiss2),]

Merge$Tiss2 <- factor(Merge$Tiss2,levels=c("Brain","Muscle","Heart","Lung"))

ggplot(Merge, aes(x=slope_tis1, y=slope_tis2, colour=abs(logP.tis2),shape=ifQTL_tiss2,size=ifQTL_tiss2)) + geom_point() + theme_bw() + geom_hline(yintercept = 0) + geom_vline(xintercept = 0)+xlim(-2,2) +ylim(-2,2) +
  theme(aspect.ratio=1,axis.text = element_text(face="bold", size=rel(1.4)),axis.title = element_text(face="bold", size=rel(1.4)),legend.text = element_text(size=14, face="bold"),legend.title =element_text(size=14, face="bold") ) + scale_colour_gradient(low="grey", high="darkred")+
  facet_grid(Tiss1 ~ Tiss2)+scale_size_discrete(range = c(0.5,1.2)) + ggtitle("Tissue1 haQTL effect in Tissue 2, top QTLs")

# group the haQTLs into p-value bins based on the p-values in the replication tissues
Merge$P_tiss2_bin = 0
for(i in c(0.5:4.5)){
  Merge[abs(Merge$logP.tis2)>= i-0.5 & abs(Merge$logP.tis2)< i+0.5,]$P_tiss2_bin = i
}
Merge[abs(Merge$logP.tis2)>=5,]$P_tiss2_bin = 5.5

mycolor=c("cadetblue3","pink","pink2","pink3","pink4","darkred")

ggplot(Merge, aes(x=slope_tis1, y=slope_tis2, colour=as.character(P_tiss2_bin),shape=ifQTL_tiss2,size=ifQTL_tiss2)) + geom_point() + theme_bw() + geom_hline(yintercept = 0) + geom_vline(xintercept = 0)+xlim(-2,2) +ylim(-2,2) +
  theme(aspect.ratio=1,axis.text = element_text(face="bold", size=rel(1.4)),axis.title = element_text(face="bold", size=rel(1.4)),legend.text = element_text(size=14, face="bold"),legend.title =element_text(size=14, face="bold") ) +
  facet_grid(Tiss1 ~ Tiss2)+scale_size_discrete(range = c(0.5,0.8)) + ggtitle("Tissue1 haQTL effect in Tissue 2, top QTLs")+scale_color_manual(values = mycolor)

####### further check effect size correlation, partitioned by p-value in tissue2
Merge$Tiss2 <- factor(Merge$Tiss2,levels=c("Brain","Muscle","Heart","Lung"))
Merge$Tiss1 <- factor(Merge$Tiss1,levels=c("Brain","Muscle","Heart","Lung"))
Merge.diffTis <- Merge[Merge$Tiss1 != Merge$Tiss2,]
dim(Merge.diffTis)
# exclude lung
Merge.diffTis = Merge.diffTis[Merge.diffTis$Tiss1 !="Lung" & Merge.diffTis$Tiss2 !="Lung",]
dim(Merge.diffTis)

Merge.diffTis$P_tiss2_bin = 0
for(i in c(0.5:4.5)){
  Merge.diffTis[abs(Merge.diffTis$logP.tis2)>= i-0.5 & abs(Merge.diffTis$logP.tis2)< i+0.5,]$P_tiss2_bin = i
}
Merge.diffTis[abs(Merge.diffTis$logP.tis2)>=5,]$P_tiss2_bin = 5.5

df=data.frame()
for(i in c(0.5:5.5)){
  Merge.diffTis.bin = Merge.diffTis[Merge.diffTis$P_tiss2_bin == i,]
  out=summary(lm(Merge.diffTis.bin$slope_tis2~Merge.diffTis.bin$slope_tis1))
  pcc=cor.test(Merge.diffTis.bin$slope_tis1,Merge.diffTis.bin$slope_tis2)
  df0=data.frame(pval.bin=i,beta=out$coefficients[2,1],beta.sd=out$coefficients[2,2],R2=out$adj.r.squared,pearson=pcc$estimate,pearson.ci.low=pcc$conf.int[1],pearson.ci.high=pcc$conf.int[2])
  df=rbind(df,df0)
}
print(df)

df$cate = "non_sharing"
df[df$pval.bin==2.5|df$pval.bin==3.5|df$pval.bin==4.5,]$cate = "sub_sharing"
df[df$pval.bin>4.5,]$cate = "Shared"
df$cate = factor(df$cate,levels = c("non_sharing","sub_sharing","Shared"))
mycolor=c("cadetblue3","pink3","darkred")

pd <- position_dodge(width = 0.2)
mycolor=c("cadetblue3","pink","pink2","pink3","pink4","darkred")
ggplot(df,aes(x=pval.bin,y=beta,fill=as.character(pval.bin))) + geom_bar(stat="identity") +mytheme6+scale_fill_manual(values = mycolor)+
  geom_errorbar(aes(ymin=beta-beta.sd, ymax=beta+beta.sd), colour="black" , width=.15,position=pd,size=1)

####### p val in Tis2 vs. slope in Tis2, separated by the direction in Tissue 1
Merge$Tiss2 <- factor(Merge$Tiss2,levels=c("Brain","Muscle","Heart","Lung"))
Merge$Tiss1 <- factor(Merge$Tiss1,levels=c("Brain","Muscle","Heart","Lung"))
Merge.diffTis <- Merge[Merge$Tiss1 != Merge$Tiss2,]
dim(Merge.diffTis)
Merge.diffTis$logP.tis2.dirInTis1  <- -log10(Merge.diffTis$nominalP_tis2) * Merge.diffTis$slope_tis1/abs(Merge.diffTis$slope_tis1)
ggplot(Merge.diffTis,aes(x=logP.tis2.dirInTis1,y=slope_tis2,colour=abs(logP.tis2),shape=ifQTL_tiss2)) + geom_point() + theme_bw() + mytheme5 + 
  theme(aspect.ratio = .6)+ scale_colour_gradient(high="darkred", low="grey")+
  xlim(-22,18)+ylim(-2,2) + facet_grid(Tiss1 ~ Tiss2) #+geom_hline(yintercept = 0) + geom_vline(xintercept = 0) +   geom_vline(xintercept = -5) + geom_vline(xintercept = 5)+  geom_vline(xintercept = log10(0.005)) + geom_vline(xintercept = -log10(0.005)) + ggtitle("Tissue1 haQTL effect in Tissue 2, top QTLs")

mycolor=c("cadetblue3","pink","pink2","pink3","pink4","darkred")
ggplot(Merge.diffTis,aes(x=logP.tis2.dirInTis1,y=slope_tis2,colour=as.character(P_tiss2_bin),shape=ifQTL_tiss2)) + geom_point(size=1) + theme_bw() + mytheme5 + 
  theme(aspect.ratio = .6)+ geom_vline(xintercept = 0,color="grey",linetype="dashed")+geom_hline(yintercept = 0,color="grey",linetype="dashed")+
  xlim(-22,18)+ylim(-2,2) + facet_grid(Tiss1 ~ Tiss2) +scale_color_manual(values = mycolor)


ggplot(Merge.diffTis,aes(x=logP.tis2.dirInTis1,y=slope_tis2,colour=abs(logP.tis2),shape=ifQTL_tiss2)) + geom_point() + theme_bw() + mytheme5 + 
  theme(aspect.ratio = .6)+ scale_colour_gradient(high="darkred", low="grey")+
  xlim(-22,18)+ylim(-2,2)  #+geom_hline(yintercept = 0) + geom_vline(xintercept = 0) +   geom_vline(xintercept = -5) + geom_vline(xintercept = 5)+  geom_vline(xintercept = log10(0.005)) + geom_vline(xintercept = -log10(0.005)) + ggtitle("Tissue1 haQTL effect in Tissue 2, top QTLs")

# cumulative directionality changes
Merge.diffTis = Merge.diffTis[order(Merge.diffTis$nominalP_tis2),]
Merge.diffTis$If.consistent = "No"
Merge.diffTis$If.consistent[Merge.diffTis$slope_tis1*Merge.diffTis$slope_tis2>0] = "Yes"
Merge.diffTis$Consist_ratio = 1
for(i in 1:nrow(Merge.diffTis)){
  Merge.diffTis.sub = Merge.diffTis[1:i,]
  Merge.diffTis[i,"Consistency"] = nrow(Merge.diffTis.sub[Merge.diffTis.sub$If.consistent=="Yes",])/nrow(Merge.diffTis.sub)
}

# define tissue sharing haQTLs by letting consistency > 0.95
tis_share_cut = tail(Merge.diffTis[Merge.diffTis$Consistency>=0.95,],1)[1,"nominalP_tis2"]
ggplot(Merge.diffTis,aes(x= log10(nominalP_tis2),y=Consistency)) + geom_point() + 
  mytheme10 + geom_hline(yintercept = 0.95,linetype="dashed",color="red") +
  geom_vline(xintercept = log10(tis_share_cut),linetype="dashed",color="red")

ggplot(Merge.diffTis,aes(x= log10(nominalP_tis2))) +  geom_histogram(fill="cornsilk", colour="grey60", size=.5,bins=20) + 
  mytheme10 + geom_hline(yintercept = 0.95,linetype="dashed",color="red") +
  geom_vline(xintercept = log10(tis_share_cut),linetype="dashed",color="red")