Overlap Molecular QTLs

Upload Data:

Library

library(workflowr)

This is workflowr version 1.1.1
Run ?workflowr for help getting started

library(reshape2)
library(tidyverse)

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

✔ ggplot2 3.0.0     ✔ purrr   0.2.5
✔ tibble  1.4.2     ✔ dplyr   0.7.6
✔ tidyr   0.8.1     ✔ stringr 1.3.1
✔ readr   1.1.1     ✔ forcats 0.3.0

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

Permuted Results from APA:

nuclearAPA=read.table("../data/perm_QTL_trans/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Nuclear_transcript_permResBH.txt", stringsAsFactors = F, header = T)
totalAPA=read.table("../data/perm_QTL_trans/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total_transcript_permResBH.txt", stringsAsFactors = F, header=T)  

Permuted results for other QTLs

perm_names=c("pid" ,"nvar","shape1" ,"shape2", "dummy","sid" ,"dist","npval", "slope" , "ppval" ,"bpval")
su30=read.table("../data/other_qtls/fastqtl_qqnorm_4su30.fixed.perm.out", stringsAsFactors = F,col.names = perm_names)
su60=read.table("../data/other_qtls/fastqtl_qqnorm_4su60.fixed.perm.out", stringsAsFactors = F,col.names = perm_names)

rna=read.table("../data/other_qtls/fastqtl_qqnorm_RNAseq_phase2.fixed.perm.out", stringsAsFactors = F,col.names = perm_names)
rib=read.table("../data/other_qtls/fastqtl_qqnorm_ribo_phase2.fixed.perm.out", stringsAsFactors = F,col.names = perm_names)
prot=read.table("../data/other_qtls/fastqtl_qqnorm_prot.fixed.perm.out", stringsAsFactors = F,col.names = perm_names)

Create overlap plot

I will write this in multiple functions and put them together. The first function will take in the permuted results and return the significant snps at a given FDR.

#returns significant snps given a file and a cutoff 
sigsnp=function(file, cutoff){
  file$bh=p.adjust(file$bpval, method="fdr")
  file_sig=file %>% filter(-log10(bh)> cutoff) %>% select(sid)
  return(file_sig)
}

testsigsnp=sigsnp(rna,1 )

Next step is to choose a random subset with the same number of snps as were found significant.

#takes the file and the list of sig snps, returns a df with the same number of random snps  
randomsnps=function(file, SigSnpList){
  nsnp=nrow(SigSnpList)
  randomSnpDF= file %>% sample_n(nsnp) %>% arrange(sid) %>% select(sid)
  return(randomSnpDF)
}
testrandomsnps=randomsnps(rna, testsigsnp)

The next step is to filter permuted file by the snp id. To do this I will join on the snpIDs then group by the snp ids. I should then be able to take the lowest Pvalue from each group and count how many are in each group to multiply by the number of tests. I will practice this with a small set then make the general function.

#filter and fix pvals
filt_tot= totalAPA %>% semi_join(testrandomsnps, by="sid") %>% group_by(sid) %>% add_tally() %>% ungroup() %>% mutate(corrPval=bpval * n)
#take top snp
filt_tot_top= filt_tot %>% group_by(sid) %>% top_n(-1, corrPval)

Make this into a function for the total and nuclear:

#takes a list of snps and filters the top corrected snp for each one, returns df
top_Total=function(snp_list){
  filt_tot=totalAPA %>% semi_join(snp_list, by="sid") %>% group_by(sid) %>% add_tally() %>% ungroup() %>% mutate(corrPval=bpval)
  filt_tot_top= filt_tot %>% group_by(sid) %>% top_n(-1, corrPval)
  return(filt_tot_top)
}

#same for nuclear:  
top_Nuclear=function(snp_list){
  filt_nuc=nuclearAPA %>% semi_join(snp_list, by="sid") %>% group_by(sid) %>% add_tally() %>% ungroup() %>% mutate(corrPval=bpval)
  filt_nuc_top= filt_nuc %>% group_by(sid) %>% top_n(-1, corrPval)
  return(filt_nuc_top)
}

In the full script I will run this on the real QTLs and the random snps.

The next function will make the plots. I will make one that takes the results of the top_total or top_Nuclear snps.

#function returns a QQplot when given the results of the top_X functions. One will be the test set (real QTLs) and 1 will be the baseline snps.

makeQQ=function(test, baseline, Mol, Fraction){
  plot=qqplot(-log10(runif(nrow(baseline))), -log10(baseline$corrPval), ylab="Observed", xlab="Expected", main=paste("Overlap QTLs:", Mol, "with APA", Fraction, sep=" "))
  points(sort(-log10(runif(nrow(test)))), sort(-log10(test$corrPval)), col= alpha("Red"))
  abline(0,1)
  return(plot)
}

Put these together in a function: I want to give the function the molQTL file and it will make the total and nuclear plots. This means I need to give it the file to write the png files to.

overlapQTLplot=function(mol_file, cut, mol_type, totQQfile, nucQQfile){

    #helper functions  
    sigsnp=function(file, cutoff){
      file$bh=p.adjust(file$bpval, method="fdr")
      file_sig=file %>% filter(-log10(bh)> cutoff) %>% select(sid)
      print(paste("Sig snps=", nrow(file_sig), sep=" "))
      return(file_sig)
    }
    randomsnps=function(file, SigSnpList){
      nsnp=nrow(SigSnpList)
      randomSnpDF= file %>% sample_n(nsnp) %>% arrange(sid) %>% select(sid)
      return(randomSnpDF)
    }
    top_Total=function(snp_list){
      filt_tot=totalAPA %>% semi_join(snp_list, by="sid") %>% group_by(sid) %>% add_tally() %>% ungroup() %>% mutate(corrPval=bpval)
      filt_tot_top= filt_tot %>% group_by(sid) %>% top_n(-1, corrPval)
      print(paste("Total overlap=", nrow(filt_tot_top), sep=" "))
      return(filt_tot_top)
    }
    top_Nuclear=function(snp_list){
      filt_nuc=nuclearAPA %>% semi_join(snp_list, by="sid") %>% group_by(sid) %>% add_tally() %>% ungroup() %>% mutate(corrPval=bpval)
      filt_nuc_top= filt_nuc %>% group_by(sid) %>% top_n(-1, corrPval)
      print(paste("Nuclear overlap=", nrow(filt_nuc_top), sep=" "))
      return(filt_nuc_top)
    }
    makeQQ=function(test, baseline, Mol, Fraction){
      plot=qqplot(-log10(runif(nrow(baseline))), -log10(baseline$corrPval), ylab="Observed", xlab="Expected", main=paste("Overlap QTLs:", Mol, "with APA", Fraction, sep=" "))
      points(sort(-log10(runif(nrow(test)))), sort(-log10(test$corrPval)), col= alpha("Red"))
      abline(0,1)
      return(plot)
    }
    TL=sigsnp(mol_file, cut)
    BL=randomsnps(mol_file, TL)
    #top snps test and base total
    topT_T=top_Total(TL)
    topT_B=top_Total(BL)
    #top snps test and base total
    topN_T=top_Nuclear(TL)
    topN_B=top_Nuclear(BL)
    
    #plot Total
    png(totQQfile)
    totalPlot=makeQQ(topT_T,topT_B, mol_type, "Total")
    dev.off()
    #plot Nuc
    png(nucQQfile)
    totalPlot=makeQQ(topN_T,topN_B, mol_type, "Nuclear")

}

overlapQTLplot(su30, 1, "4su 30", "../output/plots/APAoverlap4su30_Total.png","../output/plots/APAoverlap4su30_Nuclear.png")

[1] "Sig snps= 711"
[1] "Total overlap= 85"
[1] "Total overlap= 112"
[1] "Nuclear overlap= 87"
[1] "Nuclear overlap= 100"

overlapQTLplot(su60, 1, "4su 60", "../output/plots/APAoverlap4su60_Total.png","../output/plots/APAoverlap4su60_Nuclear.png")

[1] "Sig snps= 666"
[1] "Total overlap= 91"
[1] "Total overlap= 90"
[1] "Nuclear overlap= 76"
[1] "Nuclear overlap= 104"

overlapQTLplot(rna, 1, "RNAseq", "../output/plots/APAoverlapRNA_Total.png","../output/plots/APAoverlapRNA_Nuclear.png")

[1] "Sig snps= 847"
[1] "Total overlap= 107"
[1] "Total overlap= 116"
[1] "Nuclear overlap= 109"
[1] "Nuclear overlap= 142"

overlapQTLplot(rib, 1, "RiboSeq", "../output/plots/APAoverlapRibo_Total.png","../output/plots/APAoverlapRibo_Nuclear.png")

[1] "Sig snps= 483"
[1] "Total overlap= 57"
[1] "Total overlap= 72"
[1] "Nuclear overlap= 61"
[1] "Nuclear overlap= 79"

overlapQTLplot(rib, 1, "Protein", "../output/plots/APAoverlapProtein_Total.png","../output/plots/APAoverlapProtein_Nuclear.png")

[1] "Sig snps= 483"
[1] "Total overlap= 57"
[1] "Total overlap= 67"
[1] "Nuclear overlap= 61"
[1] "Nuclear overlap= 77"

Insert plots

Session information

sessionInfo()

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

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] bindrcpp_0.2.2  forcats_0.3.0   stringr_1.3.1   dplyr_0.7.6    
 [5] purrr_0.2.5     readr_1.1.1     tidyr_0.8.1     tibble_1.4.2   
 [9] ggplot2_3.0.0   tidyverse_1.2.1 reshape2_1.4.3  workflowr_1.1.1

loaded via a namespace (and not attached):
 [1] tidyselect_0.2.4  haven_1.1.2       lattice_0.20-35  
 [4] colorspace_1.3-2  htmltools_0.3.6   yaml_2.2.0       
 [7] rlang_0.2.2       R.oo_1.22.0       pillar_1.3.0     
[10] glue_1.3.0        withr_2.1.2       R.utils_2.7.0    
[13] modelr_0.1.2      readxl_1.1.0      bindr_0.1.1      
[16] plyr_1.8.4        munsell_0.5.0     gtable_0.2.0     
[19] cellranger_1.1.0  rvest_0.3.2       R.methodsS3_1.7.1
[22] evaluate_0.11     knitr_1.20        broom_0.5.0      
[25] Rcpp_0.12.18      backports_1.1.2   scales_1.0.0     
[28] jsonlite_1.5      hms_0.4.2         digest_0.6.16    
[31] stringi_1.2.4     grid_3.5.1        rprojroot_1.3-2  
[34] cli_1.0.0         tools_3.5.1       magrittr_1.5     
[37] lazyeval_0.2.1    crayon_1.3.4      whisker_0.3-2    
[40] pkgconfig_2.0.2   xml2_1.2.0        lubridate_1.7.4  
[43] assertthat_0.2.0  rmarkdown_1.10    httr_1.3.1       
[46] rstudioapi_0.7    R6_2.2.2          nlme_3.1-137     
[49] git2r_0.23.0      compiler_3.5.1