propeQTLs_explained.Rmd

---
title: "Proportion eQTLs explained"
author: "Briana Mittleman"
date: "6/7/2019"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


```{r}
library(tidyverse)
library(workflowr)
library(reshape2)
```

I need to fix the explained_FDR10.sort.txt and unexplained_FDR10.sort.txt files because right now this file has multiple genes per snp.  

```{bash,eval=F}
python fixExandUnexeQTL.py ../data/Li_eQTLs/explained_FDR10.sort.txt ../data/Li_eQTLs/explained_FDR10.sort_FIXED.txt
python fixExandUnexeQTL.py ../data/Li_eQTLs/unexplained_FDR10.sort.txt ../data/Li_eQTLs/unexplained_FDR10.sort_FIXED.txt
```


There are 1195 explained and 814 unexplained eQTLs. I will next look at each of these in my apadata.   

Convert nominal results to have snps rather than rsids:  
```{bash,eval=F}
python convertNominal2SNPLOC.py Total
python convertNominal2SNPLOC.py Nuclear
```

```{bash,eval=F}
mkdir ../data/overlapeQTL_try2
sbatch run_getapafromeQTL.sh


```


###total  
I can group the unexplained by gene and snp then I can ask if there is at least 1 significat peak for each of these.

I will use the bonforoni correction here and multiply the pvalue by the number of peaks in the gene:snp association. 

```{r}
nomnames=c("peakID", 'snp','dist', 'pval', 'slope')
totalapaUnexplained=read.table("../data/overlapeQTL_try2/apaTotal_unexplainedQTLs.txt", stringsAsFactors = F, col.names = nomnames)
totalapaUnexplained=totalapaUnexplained %>% separate(peakID, into=c("chr","start","end","geneID"), sep=":") %>% separate(geneID, into=c("gene", "loc", "strand", "PASnum"), sep="_") %>% group_by(gene, snp)  %>% mutate(nPeaks=n(), adjPval=pval* nPeaks)%>%  dplyr::slice(which.min(adjPval))

totalapaUnexplained_sig= totalapaUnexplained %>% filter(adjPval<.05)
```

Look at distribution of these pvals:

```{r}
ggplot(totalapaUnexplained, aes(x=adjPval)) + geom_histogram(bins=50)
```

Proportion explained:

```{r}
nrow(totalapaUnexplained_sig)/nrow(totalapaUnexplained)
```




Compare to explained eQTLS:  

```{r}
totalapaexplained=read.table("../data/overlapeQTL_try2/apaTotal_explainedQTLs.txt", stringsAsFactors = F, col.names = nomnames) %>% separate(peakID, into=c("chr","start","end","geneID"), sep=":") %>% separate(geneID, into=c("gene", "loc", "strand", "PASnum"), sep="_") %>% group_by(gene, snp) %>%  mutate(nPeaks=n(), adjPval=pval* nPeaks) %>%  dplyr::slice(which.min(adjPval))

totalapaexplained_sig= totalapaexplained %>% filter(adjPval<.05)

nrow(totalapaexplained_sig)/nrow(totalapaexplained)
```

difference of proportions:

```{r}
prop.test(x=c(nrow(totalapaUnexplained_sig),nrow(totalapaexplained_sig)), n=c(nrow(totalapaUnexplained),nrow(totalapaexplained)))
```
```{r}
ggplot(totalapaUnexplained_sig,aes(x=loc)) + geom_histogram(stat="count",aes(y=..count../sum(..count..))) + labs(y="Proportion", title = "Total apaQTLs explaining eQTLs")



```


```{r}




totalapaUnexplained_sig_loc= totalapaUnexplained_sig %>% group_by(loc) %>% summarise(nLocTotalUn=n()) %>% mutate(propTotalUn=nLocTotalUn/nrow(totalapaUnexplained_sig))
totalapaexplained_sig_loc= totalapaexplained_sig %>% group_by(loc) %>% summarise(nLocTotalEx=n()) %>% mutate(propTotalEx=nLocTotalEx/nrow(totalapaexplained_sig))

BothTotalLoc=totalapaUnexplained_sig_loc %>% full_join(totalapaexplained_sig_loc,by="loc") %>%  replace_na(list(propTotalUn = 0, nLocTotalUn = 0,propTotalEx=0,nLocTotalEx=0  ))

BothTotalLoc
```


###nuclear  

```{r}
nuclearapaUnexplained=read.table("../data/overlapeQTL_try2/apaNuclear_unexplainedQTLs.txt", stringsAsFactors = F, col.names = nomnames) %>% separate(peakID, into=c("chr","start","end","geneID"), sep=":") %>% separate(geneID, into=c("gene", "loc", "strand", "PASnum"), sep="_") %>% group_by(gene, snp)  %>%  mutate(nPeaks=n(), adjPval=pval* nPeaks) %>% dplyr::slice(which.min(adjPval))

nuclearapaUnexplained_sig= nuclearapaUnexplained %>% filter(adjPval<.05)

nrow(nuclearapaUnexplained_sig)/nrow(nuclearapaUnexplained)
```



```{r}
nuclearapaexplained=read.table("../data/overlapeQTL_try2/apaNuclear_explainedQTLs.txt", stringsAsFactors = F, col.names = nomnames) %>% separate(peakID, into=c("chr","start","end","geneID"), sep=":") %>% separate(geneID, into=c("gene", "loc", "strand", "PASnum"), sep="_") %>% group_by(gene, snp) %>%  mutate(nPeaks=n(), adjPval=pval* nPeaks) %>%  dplyr::slice(which.min(adjPval))

nuclearapaexplained_sig= nuclearapaexplained %>% filter(adjPval<.05)

nrow(nuclearapaexplained_sig)/nrow(nuclearapaexplained)
```


```{r}
prop.test(x=c(nrow(nuclearapaUnexplained_sig),nrow(nuclearapaexplained_sig)), n=c(nrow(nuclearapaUnexplained),nrow(nuclearapaexplained)))
```


```{r}
ggplot(nuclearapaUnexplained_sig,aes(x=loc))  + geom_histogram(stat="count",aes(y=..count../sum(..count..))) + labs(title = "Nuclear apaQTLs explaining eQTLs", y="Proportion")

```
```{r}
nuclearapaUnexplained_sig_loc= nuclearapaUnexplained_sig %>% group_by(loc) %>% summarise(nLocnuclearUn=n()) %>% mutate(propnuclearUn=nLocnuclearUn/nrow(nuclearapaUnexplained_sig))
nuclearapaexplained_sig_loc= nuclearapaexplained_sig %>% group_by(loc) %>% summarise(nLocnuclearEx=n()) %>% mutate(propnuclearEx=nLocnuclearEx/nrow(nuclearapaexplained_sig))

BothnuclearLoc=nuclearapaUnexplained_sig_loc %>% full_join(nuclearapaexplained_sig_loc,by="loc") %>%  replace_na(list(propnuclearUn = 0, nLocnuclearUn = 0,propnuclearEx=0,nLocnuclearEx=0  ))

BothnuclearLoc
```



###total v nuclear  

```{r}
prop.test(x=c(nrow(nuclearapaUnexplained_sig),nrow(totalapaUnexplained_sig)), n=c(nrow(nuclearapaUnexplained),nrow(totalapaUnexplained)))
```

Differences in proportion by location
```{r}
allLocProp=BothnuclearLoc %>% full_join(BothTotalLoc, by="loc") %>% select(loc,propnuclearUn,propnuclearEx,propTotalUn,propTotalEx )

allLocPropmelt= melt(allLocProp, id.vars = "loc") %>% mutate(Fraction=ifelse(grepl("Total", variable), "Total", "Nuclear"),eQTL=ifelse(grepl("Un", variable), "Unexplained", "Explained"))


ggplot(allLocPropmelt,aes(x=loc, fill=eQTL, y=value)) + geom_histogram(stat="identity", position = "dodge") + facet_grid(~Fraction)+ labs(y="Proportion of PAS", title="apaQTLs overlaping eQTLs by PAS location")  + scale_fill_manual(values=c("orange", "blue"))
```

This is a very stringent test. A less stringent way to get an upper bound would be to make an informed decision about which peak to use. This will make it so I am only testing one PAS per gene.  

##Vary the pvalue cuttoff  

To test if .05 is a good cuttoff for this analysis I will create a function that computes the overlap at different cutoffs. I will go from .01 to .5 by .05  


totalapaUnexplained
totalapaexplained


nuclearapaUnexplained
nuclearapaexplained

```{r}
prop_overlap=function(status, fraction, cutoff){
  if (fraction=="Total"){
    if (status=="Explained"){
      file=totalapaexplained
      sig=file %>% filter(adjPval<=cutoff)
      proportion=round(nrow(sig)/nrow(file),digits=2)
    }else {
      file=totalapaUnexplained
      sig=file %>% filter(adjPval<=cutoff)
      proportion=round(nrow(sig)/nrow(file),digits=2)
    }
  } else{
    if (status=="Explained"){
      file=nuclearapaexplained
      sig=file %>% filter(adjPval<=cutoff)
      proportion=round(nrow(sig)/nrow(file),digits=2)
     }else {
      file=nuclearapaUnexplained
      sig=file %>% filter(adjPval<=cutoff)
      proportion=round(nrow(sig)/nrow(file),digits=2)
     }
  }
  return(proportion)
}

```

```{r}
cutoffs=c(0.001,0.01,0.02,0.03,0.04,0.05,0.1,0.2,0.3,0.4,0.5)

TotalExplained_Proportions=c()
for(i in cutoffs){
  TotalExplained_Proportions=c( TotalExplained_Proportions, prop_overlap("Explained", "Total", i))
}
TotalExplained_ProportionsDF=as.data.frame(cbind(cutoffs,Prop=TotalExplained_Proportions, Status=rep("Explained", 11), Fraction=rep("Total", 11)))

TotalUnexplained_Proportions=c()
for(i in cutoffs){
  TotalUnexplained_Proportions=c(TotalUnexplained_Proportions, prop_overlap("Unexplained", "Total", i))
}
TotalUnexplained_ProportionsDF=as.data.frame(cbind(cutoffs,Prop=TotalUnexplained_Proportions, Status=rep("Unexplained", 11), Fraction=rep("Total", 11)))

NuclearExplained_Proportions=c()
for(i in cutoffs){
  NuclearExplained_Proportions=c( NuclearExplained_Proportions, prop_overlap("Explained", "Nuclear", i))
}
NuclearExplained_ProportionsDF=as.data.frame(cbind(cutoffs,Prop=NuclearExplained_Proportions, Status=rep("Explained", 11), Fraction=rep("Nuclear", 11)))


NuclearUnexplained_Proportions=c()
for(i in cutoffs){
  NuclearUnexplained_Proportions=c( NuclearUnexplained_Proportions, prop_overlap("Unexplained", "Nuclear", i))
}
NuclearUnexplained_ProportionsDF=as.data.frame(cbind(cutoffs,Prop=NuclearUnexplained_Proportions, Status=rep("Unexplained", 11), Fraction=rep("Nuclear", 11)))



AllPropDF=bind_rows(TotalExplained_ProportionsDF,TotalUnexplained_ProportionsDF,NuclearExplained_ProportionsDF,NuclearUnexplained_ProportionsDF)

AllPropDF$Prop=as.numeric(AllPropDF$Prop)
```
Plot this:

```{r}
ggplot(AllPropDF, aes(x=cutoffs, y=Prop, fill=Status)) + geom_bar(position = "dodge", stat="identity") + facet_grid(~Fraction) + labs(title="Proportion of eQTLs explained by apaQTLs", y="Proportion", "P-Value cut off") + scale_fill_manual(values=c("orange", "blue"))
```

##Concordance of directions for intronic pas usage and eQTL  

I want to look at the intronic pas and the eQTLs. To do this I want to look at a correaltion of effect sizes for the eQTLs and and intronic PAS. 

The eQTL information is in ../data/molQTLs/fastqtl_qqnorm_RNAseq_phase2.fixed.nominal.AllNomRes.GeneName.txt. I need to converte the RSID into snp loc.  

```{bash,eval=F}
python eQTL_switch2snploc.py
```


prepare eQTL:

```{r}
eQTLeffect=read.table("../data/molQTLs/fastqtl_qqnorm_RNAseq_phase2.fixed.nominal.AllNomRes.GeneName_snploc.txt", stringsAsFactors = F, col.names = c("gene","snp","dist", "pval", "eQTL_es")) %>% select(gene, snp, eQTL_es)
```

total: 
```{r}

#totalunex_all=read.table("../data/overlapeQTL_try2/apaTotal_unexplainedQTLs.txt", stringsAsFactors = F, col.names = nomnames) %>% separate(peakID, into=c("chr","start","end","geneID"), sep=":") %>% separate(geneID, into=c("gene", "loc", "strand", "PASnum"), sep="_")

#totalex_all=read.table("../data/overlapeQTL_try2/apaTotal_explainedQTLs.txt", stringsAsFactors = F, col.names = nomnames) %>%  separate(peakID, into=c("chr","start","end","geneID"), sep=":") %>% separate(geneID, into=c("gene", "loc", "strand", "PASnum"), sep="_")


alleQTLS_total=bind_rows(totalapaUnexplained, totalapaexplained) %>% filter(loc=="intron") %>% inner_join(eQTLeffect, by=c("gene","snp"))


ggplot(alleQTLS_total,aes(x=eQTL_es, y=slope)) + geom_point() + geom_smooth(method = "lm") +geom_text(y=-1, x=-1.5, label="slope: -0.22 p-value: 0.00002, r2=0.08") + labs(title="Total apa effect sizes vs eQTL eqtl effect sizes", y="Total apaQTL effect size",x="eQTL effect size")

summary(lm(alleQTLS_total$slope ~alleQTLS_total$eQTL_es))
```

Nuclear:

```{r}
alleQTLS_nuclear=bind_rows(nuclearapaUnexplained,nuclearapaexplained) %>% filter(loc=="intron") %>% inner_join(eQTLeffect, by=c("gene","snp"))


ggplot(alleQTLS_nuclear,aes(x=eQTL_es, y=slope)) + geom_point() + geom_smooth(method = "lm") +geom_text(y=1.5, x=-1, label="slope: -0.20 p-value: 9.0 * 10 ^ -9, r2=0.08") + labs(title="", y="apaQTL effect size",x="eQTL effect size")


summary(lm(alleQTLS_nuclear$slope ~alleQTLS_nuclear$eQTL_es))
```

```{r figure3B, include=FALSE, dev="pdf", fig.height=10, fig.width=10, crop=FALSE}
ggplot(alleQTLS_nuclear,aes(x=eQTL_es, y=slope)) + geom_point() + geom_smooth(method = "lm") +geom_text(y=1.5, x=-1, label="slope: -0.20 p-value: 9.0 * 10 ^ -9, r2=0.08") + labs(title="", y="apaQTL effect size",x="eQTL effect size")
```

##Examples for overlap:  

```{r}
unexplained_snps=read.table("../data/Li_eQTLs/unexplained_FDR10.sort_FIXED.txt", col.names = c("chr", "loc", "gene"),stringsAsFactors = F)
totQTL=read.table("../data/apaQTLs/Total_apaQTLs4pc_5fdr.bed", header = T, stringsAsFactors = F, col.names = c("chr", "bedstart","loc","ID", "score", "strand"))
nucQTL=read.table("../data/apaQTLs/Nuclear_apaQTLs4pc_5fdr.bed", stringsAsFactors = F, header = T, col.names = c("chr", "bedstart","loc","ID", "score", "strand"))
```

Overlap:

```{r}
totQTL_unex=totQTL %>% inner_join(unexplained_snps, by=c("chr", "loc"))
nucQTL_unex=nucQTL %>% inner_join(unexplained_snps, by=c("chr", "loc"))
```


```{r}
totQTL_unex
```
```{r}
nucQTL_unex
```


Make a plot for KDELR2 7:6497501  


```{r}
genohead=as.data.frame(read.table("../data/ExampleQTLPlots/genotypeHeader.txt", stringsAsFactors = F, header = F)[,10:128] %>% t())
colnames(genohead)=c("header")
genotype=as.data.frame(read.table("../data/ExampleQTLPlots/KDELR2_TotalPeaksGenotype.txt", stringsAsFactors = F, header = F) [,10:128] %>% t())

full_geno=bind_cols(Ind=genohead$header, dose=genotype$V1) %>% mutate(numdose=round(dose), genotype=ifelse(numdose==0, "TT", ifelse(numdose==1, "TG", "GG")))

RNAhead=as.data.frame(read.table("../data/molPhenos/RNAhead.txt", stringsAsFactors = F, header = F)[,5:73] %>% t())

RNApheno=as.data.frame(read.table("../data/molPhenos/RNA_KDELr2.txt", stringsAsFactors = F, header = F) [,5:73] %>% t())

full_pheno=bind_cols(Ind=RNAhead$V1, Expression=RNApheno$V1)

allRNA=full_geno %>% inner_join(full_pheno, by="Ind")

allRNA$genotype=as.factor(allRNA$genotype)
```

Ref,T	Alt= G
```{r}
ggplot(allRNA, aes(x=genotype, y=Expression,group=genotype, fill=genotype)) + geom_boxplot() + geom_jitter()+scale_fill_brewer(palette = "YlOrRd") + labs(title="Unexplained eQTL: KDELR2 - rs6962012")
```

Make locus zoom 

```{bash,eval=F}
mkdir ../data/locusZoom
```

peak119699  KDELR2 ENSG00000136240.5
```{bash,eval=F}

grep peak119699  ../data/apaQTLNominal_4pc/APApeak_Phenotype_GeneLocAnno.Total.5perc.fc.gz.qqnorm_AllChrom.txt > ../data/locusZoom/TotalAPA.peak119699.KDELR2.nomNuc.txt

grep ENSG00000136240.5 ../data/molQTLs/fastqtl_qqnorm_RNAseq_phase2.fixed.nominal.AllNomRes.txt > ../data/locusZoom/RNA.KDELR2.txt

```


```{r}
APATotal_KDELR2_LZ=read.table("../data/locusZoom/TotalAPA.peak119699.KDELR2.nomNuc.txt", stringsAsFactors = F, col.names = c("PeakID", "SNP", "Dist", "P","slope"))  %>% select( SNP, P)

write.table(APATotal_KDELR2_LZ,"../data/locusZoom/apaTotalKDELR2_LZ.txt", col.names = T, row.names = F, quote = F)

RNA_KDELR2_LZ=read.table("../data/locusZoom/RNA.KDELR2.txt", stringsAsFactors = F, col.names = c("PeakID", "SNP", "Dist", "P","slope"))  %>% select( SNP, P)

write.table(RNA_KDELR2_LZ,"../data/locusZoom/RNAKDELR2_LZ.txt", col.names = T, row.names = F, quote = F)
```
Use locuszoom.org  

locus zoom plot for C10ofr88 variant in nuclear: 


peak19881
```{bash,eval=F}

grep peak19881  ../data/apaQTLNominal_4pc/APApeak_Phenotype_GeneLocAnno.Nuclear.5perc.fc.gz.qqnorm_AllChrom.txt > ../data/locusZoom/NuclearAPA.peak119699.C10ofr88.nomNuc.txt

grep ENSG00000119965 ../data/molQTLs/fastqtl_qqnorm_RNAseq_phase2.fixed.nominal.AllNomRes.txt > ../data/locusZoom/RNA.C10ofr88.txt

```


```{r}
APATNuclear_orf_LZ=read.table("../data/locusZoom/NuclearAPA.peak119699.C10ofr88.nomNuc.txt", stringsAsFactors = F, col.names = c("PeakID", "SNP", "Dist", "P","slope"))  %>% select( SNP, P)

write.table(APATNuclear_orf_LZ,"../data/locusZoom/apaNuclearC10orf88_LZ.txt", col.names = T, row.names = F, quote = F)

RNA_orf_LZ=read.table("../data/locusZoom/RNA.C10ofr88.txt", stringsAsFactors = F, col.names = c("PeakID", "SNP", "Dist", "P","slope"))  %>% select( SNP, P)

write.table(RNA_orf_LZ,"../data/locusZoom/RNAC10orf88_LZ.txt", col.names = T, row.names = F, quote = F)
```