Nuclear vs Total nominal associations

try with perm

totalandnuclear_commonassociationPerm.py

totRes=open("/project2/gilad/briana/threeprimeseq/data/perm_APAqtl_GeneLocAnno_noMP_5percUs/filtered_APApeaks_merged_allchrom_refseqGenes.GeneLocAnno.NoMP_sm_quant.Total.fixed.pheno_5perc_permRes.txt","r")
nucRes=open("/project2/gilad/briana/threeprimeseq/data/perm_APAqtl_GeneLocAnno_noMP_5percUs/filtered_APApeaks_merged_allchrom_refseqGenes.GeneLocAnno.NoMP_sm_quant.Nuclear.fixed.pheno_5perc_permRes.txt", "r")  


outfile=open("/project2/gilad/briana/threeprimeseq/data/NucSpecQTL/TotNucRes_Permoverlapassociations.txt", "w")

resDict={}
for ln in totRes:
    gene=ln.split()[0].split(":")[-1].split("_")[0]
    peak=ln.split()[0].split(":")[-1].split("_")[-1]
    id=gene + ":" + peak 
    pval=ln.split()[-2]
    resDict[id]={}
    resDict[id]["Total"]=pval

for ln in nucRes:
    gene=ln.split()[0].split(":")[-1].split("_")[0]
    peak=ln.split()[0].split(":")[-1].split("_")[-1]
    id=gene + ":" + peak
    pval=ln.split()[-2]
    if id in resDict.keys():
        resDict[id]["Nuclear"]=pval
    else:
        continue   
        
#now i have a double dictionary. i need to write it out. i want id, the total pval, then the nuc pval  
for outer, inner in resDict.items():
    id=outer
    outPval=[]
    for key in inner:
        outPval.append(inner[key])
    print(outPval)
    if(len(outPval))==2:
         outfile.write("%s\t%s\t%s\n"%(id, outPval[0], outPval[1]))


outfile.close()        
        

library(tidyverse)

── Attaching packages ───────────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──

✔ ggplot2 3.1.0       ✔ purrr   0.3.1  
✔ tibble  2.0.1       ✔ dplyr   0.8.0.1
✔ tidyr   0.8.3       ✔ stringr 1.4.0  
✔ readr   1.3.1       ✔ forcats 0.4.0  

Warning: package 'tibble' was built under R version 3.5.2

Warning: package 'tidyr' was built under R version 3.5.2

Warning: package 'purrr' was built under R version 3.5.2

Warning: package 'dplyr' was built under R version 3.5.2

Warning: package 'stringr' was built under R version 3.5.2

Warning: package 'forcats' was built under R version 3.5.2

── Conflicts ──────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

library(data.table)

Warning: package 'data.table' was built under R version 3.5.2

Attaching package: 'data.table'

The following objects are masked from 'package:dplyr':

    between, first, last

The following object is masked from 'package:purrr':

    transpose

library(cowplot)

Warning: package 'cowplot' was built under R version 3.5.2

Attaching package: 'cowplot'

The following object is masked from 'package:ggplot2':

    ggsave

file=fread("../data/NucSpecQTL/TotNucRes_Permoverlapassociations.txt", col.names =c("ID","Total", "Nuclear"))

cor_plot=ggplot(file, aes(x=Total, y=Nuclear)) + geom_point() + geom_density2d(na.rm = TRUE, size = 1, colour = 'red') + labs(title="Common Peak associations Total and Nuclear") 

cor_plot

Warning: Removed 20 rows containing missing values (geom_point).

Past versions of unnamed-chunk-7-1.png

This is similar to what was happening in the nominal case.

I think I need to plot just the total in nuclaer and nuclear in total. this way i can see the outliers.

names=c("SNP", "peak", "pval")
NucQTLinTot=read.table("../data/QTL_overlap/NucQTLs_inTotFractionRes.txt", stringsAsFactors = F, col.names = names)
TotQTLinNuc=read.table("../data/QTL_overlap/TotQTLs_inNucFractionRes.txt", stringsAsFactors = F, col.names = names)

QTL_names= c("pid", "nvar", "shape1", "shape2", "dummy", "SNP", "dist", "npval", "slope", "ppval", "bpval", "bh")
NucQTL=read.table("../data/ApaQTLs/NuclearapaQTLs.GeneLocAnno.noMP.5perc.10FDR.txt", col.names = QTL_names, stringsAsFactors = F) %>% separate(pid, into=c("chr", "start","end", "peakID"), sep=":") %>% separate(peakID, into=c("gene", "strand", "peak"), sep="_") %>% select(SNP, peak, ppval)

TotQTL=read.table("../data/ApaQTLs/TotalapaQTLs.GeneLocAnno.noMP.5perc.10FDR.txt", col.names = QTL_names, stringsAsFactors = F) %>% separate(pid, into=c("chr", "start","end", "peakID"), sep=":") %>% separate(peakID, into=c("gene", "strand", "peak"), sep="_") %>% select(SNP, peak, ppval)

Join these:

TotQTL_bothfrac= TotQTL %>% inner_join(TotQTLinNuc, by=c("SNP", "peak"))
colnames(TotQTL_bothfrac)= c("SNP", "peak", "Total", "Nuclear")

NucQTL_bothfrac= NucQTL %>% inner_join(NucQTLinTot, by=c("SNP", "peak"))
colnames(NucQTL_bothfrac)= c("SNP", "peak", "Nuclear", "Total")

All qtls together:

allQTL=TotQTL_bothfrac %>% bind_rows(NucQTL_bothfrac)

ggplot(allQTL, aes(x=Nuclear, y=Total)) + geom_point()

The pvalues are not consistent. I need to get the original pvalues from the nominal files.

nominal associations for each QTL:

getNominal4QTLs.py

nucQTLs="/project2/gilad/briana/threeprimeseq/data/ApaQTLs/Nuclear.apaQTLs.sort.bed"  
totQTLs="/project2/gilad/briana/threeprimeseq/data/ApaQTLs/Total.apaQTLs.sort.bed"  

nucNom="/project2/gilad/briana/threeprimeseq/data/nominal_APAqtl_GeneLocAnno_noMP_5percUs/filtered_APApeaks_merged_allchrom_refseqGenes.GeneLocAnno_NoMP_sm_quant.Nuclear.fixed.pheno_5perc.fc.gz.qqnorm_allNomRes.txt"
totNom="/project2/gilad/briana/threeprimeseq/data/nominal_APAqtl_GeneLocAnno_noMP_5percUs/filtered_APApeaks_merged_allchrom_refseqGenes.GeneLocAnno_NoMP_sm_quant.Total.fixed.pheno_5perc.fc.gz.qqnorm_allNomRes.txt"

outtot="/project2/gilad/briana/threeprimeseq/data/NucSpecQTL/TotQTLs_inTotFractionRes.txt"
outnuc="/project2/gilad/briana/threeprimeseq/data/NucSpecQTL/NucQTLs_inNucFractionRes.txt"


def sameFract(inRes, inQTL, out):
    fout=open(out, "w")
    qtl_dic={}
    #SNP is key, peak is value
    for ln in open(inQTL,"r"):
        snp=ln.split()[2]
        chrom=ln.split()[0]
        peak=ln.split()[3].split(":")[0]
        qtl=str(chrom) + ":" + str(snp)
        if qtl not in qtl_dic.keys():
            qtl_dic[qtl]=[peak]
        else: 
            qtl_dic[qtl].append(peak)
    #print(qtl_dic)
    for ln in open(inRes, "r"):
        pval=ln.split()[3]
        snp=ln.split()[1]
        peak=ln.split()[0].split(":")[3].split("_")[-1]
        if snp in qtl_dic.keys():
            if peak in qtl_dic[snp]:
                fout.write("%s\t%s\t%s\n"%(snp, peak, pval))   
    fout.close()
    

sameFract(nucNom, nucQTLs,outnuc)  
sameFract(totNom, totQTLs, outtot)  
                

run_getNominal4QTLs.sh

#!/bin/bash

#SBATCH --job-name=run_getNominal4QTLs
#SBATCH --account=pi-yangili1
#SBATCH --time=5:00:00
#SBATCH --output=run_getNominal4QTLs.out
#SBATCH --error=run_getNominal4QTLs.err
#SBATCH --partition=bigmem2
#SBATCH --mem=100G
#SBATCH --mail-type=END


module load Anaconda3
source activate three-prime-env

python getNominal4QTLs.py

I now want to plot the total vs nuclear nominal pvalues for the nuclear QTLs.

colnames(NucQTLinTot)=c("SNP", "peak", "totalPval")

NucQTLinNuc=read.table("../data/NucSpecQTL/NucQTLs_inNucFractionRes.txt", col.names = c("SNP", "peak","nuclearPval"), stringsAsFactors = F)


NucQTL_all= NucQTLinTot %>% inner_join(NucQTLinNuc, by=c("SNP", "peak"))

summary(lm(data=NucQTL_all,-log10(totalPval) ~ -log10(nuclearPval)))

Call:
lm(formula = -log10(totalPval) ~ -log10(nuclearPval), data = NucQTL_all)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.1638 -2.7988 -0.9846  1.5394 17.9462 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   4.1767     0.1719    24.3   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.824 on 494 degrees of freedom

bothPval_nucqtl=ggplot(NucQTL_all, aes(x=-log10(totalPval), y=-log10(nuclearPval))) + geom_point() + geom_density2d(na.rm = TRUE, size = 1, colour = 'red') + geom_smooth(aes(x=-log10(totalPval),y=-log10(nuclearPval)),method = "lm") + labs(title="Nominal Pvalues for Nuclear APA qtls")
#+geom_text(aes(label=SNP),hjust=0, vjust=0)

bothPval_nucqtl

ggplot(NucQTL_all, aes(x=totalPval, y=nuclearPval)) + geom_point() + labs(title="Nominal Pval results for Nuclear QTLs") + geom_density2d(na.rm = TRUE, size = 1, colour = 'red')  +geom_vline(xintercept=.05)

nNucQTLNotSigTot= NucQTL_all %>% filter(totalPval>.05) %>% nrow()
nNucQTLNotSigTot

[1] 118

For Total QTLs:

colnames(TotQTLinNuc)=c("SNP", "peak", "nuclearPval")

TotQTLinTot=read.table("../data/NucSpecQTL/TotQTLs_inTotFractionRes.txt", col.names = c("SNP", "peak","totalPval"), stringsAsFactors = F)


TotQTL_all= TotQTLinNuc %>% inner_join(TotQTLinTot, by=c("SNP", "peak"))

summary(lm(data=NucQTL_all,-log10(nuclearPval) ~ -log10(totalPval)))

Call:
lm(formula = -log10(nuclearPval) ~ -log10(totalPval), data = NucQTL_all)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.4270 -1.9985 -1.1978  0.9371 15.8122 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   7.1208     0.1385   51.41   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.082 on 494 degrees of freedom

bothPval_totqtl=ggplot(TotQTL_all, aes(x=-log10(nuclearPval), y=-log10(totalPval))) + geom_point() + geom_density2d(na.rm = TRUE, size = 1, colour = 'red') + geom_smooth(aes(x=-log10(nuclearPval),y=-log10(totalPval)),method = "lm") + labs(title="Nominal Pvalues for Total APA qtls")
#+geom_text(aes(label=SNP),hjust=0, vjust=0)

bothPval_totqtl

bothFracPlots=plot_grid(bothPval_nucqtl,bothPval_totqtl)
bothFracPlots

ggsave(bothFracPlots, file="../output/plots/NominalPvalcondQTL.png", height=7, width=12 )

ggplot(TotQTL_all, aes(x=nuclearPval, y=totalPval)) + geom_point() + labs(title="Nominal Pval results for Total QTLs") + geom_density2d(na.rm = TRUE, size = 1, colour = 'red') 

nTotQTLNotSigNuc= TotQTL_all %>% filter(nuclearPval>.05) %>% nrow()
nTotQTLNotSigNuc

[1] 39

Is the difference significant. More likely to have nuclear specific qlts?

phyper(118, 495, 272, 157, lower.tail = F)

[1] 0.0005207115

Plot for this:

compareSpecQTL=as.data.frame(cbind(fraction=c("Total QTL", "Total QTL","Nuclear QTL","Nuclear QTL"), SigInOther=c("No","Yes","No","Yes"),Value=c(.14,.86 ,.24,.76)))

compareSpecQTL$Value=as.numeric(as.character(compareSpecQTL$Value))

qtls_oppfrac=ggplot(compareSpecQTL,aes(x=fraction, fill=SigInOther, y=Value))  + geom_bar(stat="identity") +scale_fill_manual(values=c("red","grey")) + labs(y="Proportion", title="QTLs by nominal significance in other fraction")  + annotate("text", label="p < .0005", y=1.2, x=1.5) + geom_segment(aes(x=1.1, xend=1.9, y=1.1, yend=1.1))

qtls_oppfrac

ggsave(qtls_oppfrac, file="../output/plots/QTLsSiginOppFraction.png")

Saving 7 x 5 in image

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS  10.14.1

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] cowplot_0.9.4     data.table_1.12.0 forcats_0.4.0    
 [4] stringr_1.4.0     dplyr_0.8.0.1     purrr_0.3.1      
 [7] readr_1.3.1       tidyr_0.8.3       tibble_2.0.1     
[10] ggplot2_3.1.0     tidyverse_1.2.1  

loaded via a namespace (and not attached):
 [1] tidyselect_0.2.5 xfun_0.5         haven_2.1.0      lattice_0.20-38 
 [5] colorspace_1.4-0 generics_0.0.2   htmltools_0.3.6  yaml_2.2.0      
 [9] rlang_0.3.1      pillar_1.3.1     glue_1.3.0       withr_2.1.2     
[13] modelr_0.1.4     readxl_1.3.0     plyr_1.8.4       munsell_0.5.0   
[17] gtable_0.2.0     workflowr_1.2.0  cellranger_1.1.0 rvest_0.3.2     
[21] evaluate_0.13    labeling_0.3     knitr_1.21       broom_0.5.1     
[25] Rcpp_1.0.0       scales_1.0.0     backports_1.1.3  jsonlite_1.6    
[29] fs_1.2.6         hms_0.4.2        digest_0.6.18    stringi_1.3.1   
[33] grid_3.5.1       rprojroot_1.3-2  cli_1.0.1        tools_3.5.1     
[37] magrittr_1.5     lazyeval_0.2.1   crayon_1.3.4     whisker_0.3-2   
[41] pkgconfig_2.0.2  MASS_7.3-51.1    xml2_1.2.0       lubridate_1.7.4 
[45] assertthat_0.2.0 rmarkdown_1.11   httr_1.4.0       rstudioapi_0.9.0
[49] R6_2.4.0         nlme_3.1-137     git2r_0.24.0     compiler_3.5.1