APA and Protien QTL overlap

Total

I want the overlap where I started in APA qtls and found the snp in the mol file. I am starting with the total.

totAPAinsu30=read.table("../data/mol_overlap/APA2molTotal/TotAPAqtlsPval4su30.txt", header = T, stringsAsFactors = F)
totAPAinsu60=read.table("../data/mol_overlap/APA2molTotal/TotAPAqtlsPval4su60.txt", header = T, stringsAsFactors = F)
totAPAinRNA=read.table("../data/mol_overlap/APA2molTotal/TotAPAqtlsPvalRNA.txt", header = T, stringsAsFactors = F)
totAPAinRNAg=read.table("../data/mol_overlap/APA2molTotal/TotAPAqtlsPvalRNAg.txt", header = T, stringsAsFactors = F)
totAPAinRibo=read.table("../data/mol_overlap/APA2molTotal/TotAPAqtlsPvalribo.txt", header = T, stringsAsFactors = F)
totAPAinProt=read.table("../data/mol_overlap/APA2molTotal/TotAPAqtlsPvalProtein.txt", header = T, stringsAsFactors = F)


allOverlap_T=totAPAinsu30 %>%  full_join(totAPAinsu60, by=c("Gene.name", "sid")) %>%  full_join(totAPAinRNA, by=c("Gene.name", "sid")) %>%  full_join(totAPAinRNAg, by=c("Gene.name", "sid"))  %>%  full_join(totAPAinRibo, by=c("Gene.name", "sid"))  %>%  full_join(totAPAinProt, by=c("Gene.name", "sid")) 

colnames(allOverlap_T)=c("Gene.name", "sid", "su30", "su60", "RNA", "RNAg", "ribo", "prot")

plot(sort(allOverlap_T$prot))

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

plot(allOverlap_T$RNA ~ allOverlap_T$prot)

Past versions of unnamed-chunk-16-2.png

I want to make a ggplot of these where I color them by RNA pvalue:

allOverlap_T_lowP=allOverlap_T %>% dplyr::filter(prot<.05)
ggplot(allOverlap_T_lowP, aes(x=RNA, y=prot)) + geom_point()+ geom_text(aes(label=Gene.name),hjust=0, vjust=0) + geom_vline(xintercept = .05, col="red")

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

ggplot(allOverlap_T_lowP, aes(x=RNAg, y=prot)) + geom_point()+ geom_text(aes(label=Gene.name),hjust=0, vjust=0) + geom_vline(xintercept = .05, col="red")

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

I can use these to look for examples of SNPs that are significant in prot but not in RNA.

Look at some of these:

Total RNA:
* SACM1L
* EBI3
* FBXL18
* PSMF1
* COX17

Total RNAg:
* EBI3
* FBXL18
* APBB1IP * PSMF1

Look at some examples of genes that come up in both.

EBI3 peak152751 19:4236475

Expressed in B lymphocytes in response to EB virus.

grep EBI3 /project2/gilad/briana/genome_anotation_data/ensemble_to_genename.txt
#ENSG00000105246


python createQTLsnpAPAPhenTable.py 19 19:4236475  peak152751 Total
python createQTLsnpAPAPhenTable.py 19 19:4236475  peak152751  Nuclear

sbatch run_createQTLsnpMolPhenTable.sh "19" "19:4236475" "ENSG00000105246"

#into apaExamp
scp brimittleman@midway2.rcc.uchicago.edu:/project2/gilad/briana/threeprimeseq/data/ApaQTL_examples/*19:4236475* .

plotQTL_func(SNP="19:4236475", peak="peak152751", gene="ENSG00000105246")

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

FBXL18 peak291746 7:5524129

“The protein encoded by this gene is a member of a family of proteins that contain an approximately 40-amino acid F-box motif. This motif is important for interaction with SKP1 and for targeting some proteins for degradation.” genecards.org

grep FBXL18 /project2/gilad/briana/genome_anotation_data/ensemble_to_genename.txt
#ENSG00000155034


python createQTLsnpAPAPhenTable.py 7 7:5524129  peak291746 Total
python createQTLsnpAPAPhenTable.py 7 7:5524129  peak291746  Nuclear

sbatch run_createQTLsnpMolPhenTable.sh "7" "7:5524129" "ENSG00000155034"

#into apaExamp
scp brimittleman@midway2.rcc.uchicago.edu:/project2/gilad/briana/threeprimeseq/data/ApaQTL_examples/*7:5524129* .

plotQTL_func(SNP="7:5524129", peak="peak291746", gene="ENSG00000155034")

Warning: Removed 4 rows containing non-finite values (stat_boxplot).

Warning: Removed 4 rows containing non-finite values (stat_compare_means).

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

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

• PSMF1 peak193648 20:1131308

This gene codes the 26S proteasome.

grep PSMF1 /project2/gilad/briana/genome_anotation_data/ensemble_to_genename.txt
#ENSG00000125818


python createQTLsnpAPAPhenTable.py 20 20:1131308   peak193648 Total
python createQTLsnpAPAPhenTable.py 20 20:1131308   peak193648  Nuclear

sbatch run_createQTLsnpMolPhenTable.sh "20" "20:1131308" "ENSG00000125818"

#into apaExamp
scp brimittleman@midway2.rcc.uchicago.edu:/project2/gilad/briana/threeprimeseq/data/ApaQTL_examples/*20:1131308* .

plotQTL_func(SNP="20:1131308", peak="peak193648", gene="ENSG00000125818")

Warning: Removed 2 rows containing non-finite values (stat_boxplot).

Warning: Removed 2 rows containing non-finite values (stat_compare_means).

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

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

I want to know the number of these that are <.05 in protien and above .1 in RNA

allOverlap_T_lowP_highRNA=allOverlap_T %>% dplyr::filter(prot<.05) %>% dplyr::filter(RNA>.05)
allOverlap_T_lowP_highRNAg=allOverlap_T %>% dplyr::filter(prot<.05) %>% dplyr::filter(RNAg>.05)

8 snps with < .05 for protein. Of those 6 have RNA pvalues greater than .05

6/8

[1] 0.75

Nuclear

nucAPAinsu30=read.table("../data/mol_overlap/APA2molNuclear/NucAPAqtlsPval4su30.txt", header = T, stringsAsFactors = F)
nucAPAinsu60=read.table("../data/mol_overlap/APA2molNuclear/NucAPAqtlsPval4su60.txt", header = T, stringsAsFactors = F)
nucAPAinRNA=read.table("../data/mol_overlap/APA2molNuclear/NucAPAqtlsPvalRNA.txt", header = T, stringsAsFactors = F)
nucAPAinRNAg=read.table("../data/mol_overlap/APA2molNuclear/NucAPAqtlsPvalRNAg.txt", header = T, stringsAsFactors = F)
nucAPAinRibo=read.table("../data/mol_overlap/APA2molNuclear/NucAPAqtlsPvalribo.txt", header = T, stringsAsFactors = F)
nucAPAinProt=read.table("../data/mol_overlap/APA2molNuclear/NucAPAqtlsPvalProtein.txt", header = T, stringsAsFactors = F)


allOverlap_N=nucAPAinsu30 %>%  full_join(nucAPAinsu60, by=c("Gene.name", "sid")) %>%  full_join(nucAPAinRNA, by=c("Gene.name", "sid")) %>%  full_join(nucAPAinRNAg, by=c("Gene.name", "sid"))  %>%  full_join(nucAPAinRibo, by=c("Gene.name", "sid"))  %>%  full_join(nucAPAinProt, by=c("Gene.name", "sid")) 

colnames(allOverlap_N)=c("Gene.name", "sid", "su30", "su60", "RNA", "RNAg", "ribo", "prot")

#subset by prot < .05  

allOverlap_N_lowP=allOverlap_N %>% dplyr::filter(prot<.05)
ggplot(allOverlap_N_lowP, aes(x=RNA, y=prot)) + geom_point()+ geom_text(aes(label=Gene.name),hjust=0, vjust=0)+ geom_vline(xintercept = .05, col="red")

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

ggplot(allOverlap_N_lowP, aes(x=RNAg, y=prot)) + geom_point()+ geom_text(aes(label=Gene.name),hjust=0, vjust=0)+ geom_vline(xintercept = .05, col="red")

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

allOverlap_N_lowP_highRNA=allOverlap_N %>% dplyr::filter(prot<.05) %>% dplyr::filter(RNA>.05)
allOverlap_N_lowP_highRNAg=allOverlap_N %>% dplyr::filter(prot<.05) %>% dplyr::filter(RNAg>.05)

39 snps with < .05 for protein. Of those 28 have RNA pvalues greater than .05

28/39

[1] 0.7179487

inBothN= allOverlap_N_lowP_highRNAg %>%  inner_join(allOverlap_N_lowP_highRNA, by=c("Gene.name", "sid", "su30", "su60", "RNA", "RNAg", "ribo", "prot")) %>% arrange(desc(RNA))

inBothN$Gene.name[1:5]

[1] "MSMO1"  "LYAR"   "CD2BP2" "KDM2A"  "RINT1" 

• MSMO1 4:166260601 peak249109

contains metal binding motifs, known alternative splice isoforms

grep MSMO1 /project2/gilad/briana/genome_anotation_data/ensemble_to_genename.txt
#ENSG00000052802


python createQTLsnpAPAPhenTable.py 4 4:166260601  peak249109 Total
python createQTLsnpAPAPhenTable.py 4 4:166260601  peak249109  Nuclear

sbatch run_createQTLsnpMolPhenTable.sh "4" "4:166260601" "ENSG00000052802"

#into apaExamp
scp brimittleman@midway2.rcc.uchicago.edu:/project2/gilad/briana/threeprimeseq/data/ApaQTL_examples/*4:166260601* .

plotQTL_func(SNP="4:166260601", peak="peak249109", gene="ENSG00000052802")

Warning: Removed 3 rows containing non-finite values (stat_boxplot).

Warning: Removed 3 rows containing non-finite values (stat_compare_means).

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

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

• LYAR peak235215 4:4196045

involved in processing pre-rRNA

grep LYAR /project2/gilad/briana/genome_anotation_data/ensemble_to_genename.txt
#ENSG00000145220


python createQTLsnpAPAPhenTable.py 4 4:4196045  peak235215 Total
python createQTLsnpAPAPhenTable.py 4 4:4196045  peak235215  Nuclear

sbatch run_createQTLsnpMolPhenTable.sh "4" "4:4196045" "ENSG00000145220"

#into apaExamp
scp brimittleman@midway2.rcc.uchicago.edu:/project2/gilad/briana/threeprimeseq/data/ApaQTL_examples/*4:4196045* .

plotQTL_func(SNP="4:4196045", peak="peak235215", gene="ENSG00000145220")

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

CD2BP2 peak122237 16:29898001

From genecards: “in the cytoplasm, the encoded protein binds the cytoplasmic tail of human surface antigen CD2 via its C-terminal GYF domain, and regulate CD2-triggered T lymphocyte activation. In the nucleus, this protein is a component of the U5 small nuclear ribonucleoprotein complex and is involved in RNA splicing.”

grep CD2BP2 /project2/gilad/briana/genome_anotation_data/ensemble_to_genename.txt
#ENSG00000169217


python createQTLsnpAPAPhenTable.py 16 16:29898001 peak122237 Total
python createQTLsnpAPAPhenTable.py 16 16:29898001 peak122237  Nuclear

sbatch run_createQTLsnpMolPhenTable.sh "16" "16:29898001" "ENSG00000169217"

#into apaExamp
scp brimittleman@midway2.rcc.uchicago.edu:/project2/gilad/briana/threeprimeseq/data/ApaQTL_examples/*16:29898001* .

plotQTL_func(SNP="16:29898001", peak="peak122237", gene="ENSG00000169217")

Warning: Removed 5 rows containing non-finite values (stat_boxplot).

Warning: Removed 5 rows containing non-finite values (stat_compare_means).

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

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

KDM2A peak55622 11:66851583

grep KDM2A /project2/gilad/briana/genome_anotation_data/ensemble_to_genename.txt
#ENSG00000173120


python createQTLsnpAPAPhenTable.py 11 11:66851583 peak55622 Total
python createQTLsnpAPAPhenTable.py 11 11:66851583 peak55622  Nuclear

sbatch run_createQTLsnpMolPhenTable.sh "11" "11:66851583" "ENSG00000173120"

#into apaExamp
scp brimittleman@midway2.rcc.uchicago.edu:/project2/gilad/briana/threeprimeseq/data/ApaQTL_examples/*11:66851583* .

plotQTL_func(SNP="11:66851583", peak="peak55622", gene="ENSG00000173120")

Warning: Removed 8 rows containing non-finite values (stat_boxplot).

Warning: Removed 8 rows containing non-finite values (stat_compare_means).

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

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

RINT1 peak303436 7:105155320

Interacts with double strand break repair protiens, regulates cell cycle and telomere length

grep RINT1 /project2/gilad/briana/genome_anotation_data/ensemble_to_genename.txt
#ENSG00000135249


python createQTLsnpAPAPhenTable.py 7 7:105155320 peak303436 Total
python createQTLsnpAPAPhenTable.py 7 7:105155320 peak303436  Nuclear

sbatch run_createQTLsnpMolPhenTable.sh "7" "7:105155320" "ENSG00000135249"

#into apaExamp
scp brimittleman@midway2.rcc.uchicago.edu:/project2/gilad/briana/threeprimeseq/data/ApaQTL_examples/*7:105155320* .

plotQTL_func(SNP="7:105155320", peak="peak303436", gene="ENSG00000135249")

Warning: Removed 2 rows containing non-finite values (stat_boxplot).

Warning: Removed 2 rows containing non-finite values (stat_compare_means).

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

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

Significance

In the next step I need to add significance to the boxplots and think more about the significance cutoffs.

Maybe I can compare 2 other phenotypes for <.05 and >.05 to see if the percentage is less than what I see for RNA and protein.

sigProt=allOverlap_N_lowP %>% nrow()
paste("Sig Prot", sigProt, sep=" ")

[1] "Sig Prot 39"

allOverlap_N_lowP_higrna=allOverlap_N %>% dplyr::filter(prot<.05) %>% dplyr::filter(RNA>.05) %>%nrow()
paste("Sig Prot not RNA", allOverlap_N_lowP_higrna, sep=" ")

[1] "Sig Prot not RNA 28"

allOverlap_N_lowP_higribo=allOverlap_N %>% dplyr::filter(prot<.05) %>% dplyr::filter(ribo>.05) %>%nrow()
paste("Sig Prot not ribo", allOverlap_N_lowP_higribo, sep=" ")

[1] "Sig Prot not ribo 23"

allOverlap_N_lowP_higsu30=allOverlap_N %>% dplyr::filter(prot<.05) %>% dplyr::filter(su30>.05) %>% nrow()
paste("Sig Prot not 4su30", allOverlap_N_lowP_higsu30, sep=" ")

[1] "Sig Prot not 4su30 31"

allOverlap_N_lowP_higsu60=allOverlap_N %>% dplyr::filter(prot<.05) %>% dplyr::filter(su60>.05) %>% nrow()
paste("Sig Prot not 4su60", allOverlap_N_lowP_higsu60, sep=" ")

[1] "Sig Prot not 4su60 27"

sigProtT=allOverlap_T_lowP %>% nrow()
paste("Sig Prot", sigProtT, sep=" ")

[1] "Sig Prot 8"

allOverlap_T_lowP_higrna=allOverlap_T %>% dplyr::filter(prot<.05) %>% dplyr::filter(RNA>.05) %>%nrow()
paste("Sig Prot not RNA", allOverlap_T_lowP_higrna, sep=" ")

[1] "Sig Prot not RNA 6"

allOverlap_T_lowP_higribo=allOverlap_T %>% dplyr::filter(prot<.05) %>% dplyr::filter(ribo>.05) %>%nrow()
paste("Sig Prot not ribo", allOverlap_T_lowP_higribo, sep=" ")

[1] "Sig Prot not ribo 5"

allOverlap_T_lowP_higsu30=allOverlap_T %>% dplyr::filter(prot<.05) %>% dplyr::filter(su30>.05) %>% nrow()
paste("Sig Prot not 4su30", allOverlap_T_lowP_higsu30, sep=" ")

[1] "Sig Prot not 4su30 6"

allOverlap_T_lowP_higsu60=allOverlap_T %>% dplyr::filter(prot<.05) %>% dplyr::filter(su60>.05) %>% nrow()
paste("Sig Prot not 4su60", allOverlap_T_lowP_higsu60, sep=" ")

[1] "Sig Prot not 4su60 5"

CLoser to a global scale

In order to think about this more broadly I am going to plot overlaps:

Given all of these are APA qtls, how many are significant in:

• protein
  
• RNA and protein
• RNA but not protein
• Not in RNA or Protein

Do this first for total:

p_N=allOverlap_N %>% filter(prot<.05) %>% nrow()
RP_N=allOverlap_N %>% filter(prot<.05 & RNA <.05) %>% nrow()
RnotP_N=allOverlap_N %>% filter(prot>.05 & RNA <.05) %>% nrow()
notRnotP_N=allOverlap_N %>% filter(prot>.05 & RNA >.05) %>% nrow()
notRP_N= allOverlap_N %>% filter(prot<.05 & RNA >.05) %>% nrow()


overlapNames=c("Protein", "RNAandProtein", "RNAnotProt", "Neither", "ProteinnotRNA")
NucRPoverlap=c(p_N/nrow(allOverlap_N),RP_N/nrow(allOverlap_N), RnotP_N/nrow(allOverlap_N), notRnotP_N/nrow(allOverlap_N),notRP_N/nrow(allOverlap_N) )


NucOverlapDF=data.frame(overlapNames,NucRPoverlap)
ggplot(NucOverlapDF, aes(x=overlapNames, y=NucRPoverlap)) + geom_bar(stat="identity")

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

DO the same for total:

p_T=allOverlap_T %>% filter(prot<.05) %>% nrow()
RP_T=allOverlap_T %>% filter(prot<.05 & RNA <.05) %>% nrow()
RnotP_T=allOverlap_T %>% filter(prot>.05 & RNA <.05) %>% nrow()
notRnotP_T=allOverlap_T %>% filter(prot>.05 & RNA >.05) %>% nrow()
notRP_T= allOverlap_T %>% filter(prot<.05 & RNA >.05) %>% nrow()

TotRPoverlap=c(p_T/nrow(allOverlap_T),RP_T/nrow(allOverlap_T), RnotP_T/nrow(allOverlap_T), notRnotP_T/nrow(allOverlap_T),notRP_T/nrow(allOverlap_T))
TotOverlapDF=data.frame(overlapNames,TotRPoverlap)
ggplot(TotOverlapDF, aes(x=overlapNames, y=TotRPoverlap)) + geom_bar(stat="identity")

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

allOverlapDF=TotOverlapDF %>% full_join(NucOverlapDF, by="overlapNames")

allOverlapDF_melt=melt(allOverlapDF, id.vars="overlapNames")

ggplot(allOverlapDF_melt, aes(x=overlapNames, y=value, by=variable, fill=variable)) + geom_bar(stat="identity", position="dodge") + scale_fill_manual(values=c("darkviolet", "deepskyblue3")) + labs(y="Proportion of APA QTLs", x="Category")

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

oT=overlapplotT=draw.pairwise.venn(area1=.08, area2=.07, cross.area = .02, category = c("Protein", "RNA"),lty = rep("solid", 2), fill = c("Blue", "Orange"), alpha = rep(0.5, 2), cat.pos = c(0, 0), cat.dist = rep(0.025, 2), euler.d = T, scaled=T)

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

png("../output/plots/overlapProtRNAvenT.png") 
grid.arrange(gTree(children=oT), top="Total: Protein and RNA QTL overlap", bottom="Neither=.37")

dev.off()

quartz_off_screen 
                2 

oN=overlapplotN=draw.pairwise.venn(area1=.056, area2=.096, cross.area = .016, category = c("Protein", "RNA"),lty = rep("solid", 2), fill = c("Blue", "Orange"), alpha = rep(0.5, 2), cat.pos = c(0, 0), cat.dist = rep(0.025, 2), euler.d = T, scaled=T)

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

png("../output/plots/overlapProtRNAvenN.png") 

grid.arrange(gTree(children=oN), top="Nuclear:Protein and RNA QTL overlap", bottom="Neither=.38")
dev.off()

quartz_off_screen 
                2 

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] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] bindrcpp_0.2.2      gridExtra_2.3       VennDiagram_1.6.20 
 [4] futile.logger_1.4.3 ggpubr_0.1.8        magrittr_1.5       
 [7] cowplot_0.9.3       data.table_1.11.8   forcats_0.3.0      
[10] stringr_1.4.0       dplyr_0.7.6         purrr_0.2.5        
[13] readr_1.1.1         tidyr_0.8.1         tibble_1.4.2       
[16] ggplot2_3.0.0       tidyverse_1.2.1     workflowr_1.2.0    

loaded via a namespace (and not attached):
 [1] tidyselect_0.2.4     reshape2_1.4.3       haven_1.1.2         
 [4] lattice_0.20-35      colorspace_1.3-2     htmltools_0.3.6     
 [7] yaml_2.2.0           rlang_0.2.2          pillar_1.3.0        
[10] glue_1.3.0           withr_2.1.2          RColorBrewer_1.1-2  
[13] lambda.r_1.2.3       modelr_0.1.2         readxl_1.1.0        
[16] bindr_0.1.1          plyr_1.8.4           munsell_0.5.0       
[19] gtable_0.2.0         cellranger_1.1.0     rvest_0.3.2         
[22] evaluate_0.13        labeling_0.3         knitr_1.20          
[25] broom_0.5.0          Rcpp_0.12.19         formatR_1.5         
[28] scales_1.0.0         backports_1.1.2      jsonlite_1.6        
[31] fs_1.2.6             hms_0.4.2            digest_0.6.17       
[34] stringi_1.2.4        rprojroot_1.3-2      cli_1.0.1           
[37] tools_3.5.1          lazyeval_0.2.1       futile.options_1.0.1
[40] crayon_1.3.4         whisker_0.3-2        pkgconfig_2.0.2     
[43] xml2_1.2.0           lubridate_1.7.4      assertthat_0.2.0    
[46] rmarkdown_1.11       httr_1.3.1           rstudioapi_0.9.0    
[49] R6_2.3.0             nlme_3.1-137         git2r_0.24.0        
[52] compiler_3.5.1