Top APA QTL in other phenotypes

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()

library(VennDiagram)

Loading required package: grid

Loading required package: futile.logger

library(data.table)

Attaching package: 'data.table'

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

    between, first, last

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

    transpose

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

    dcast, melt

library(cowplot)

Attaching package: 'cowplot'

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

    ggsave

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)  

I want to use a buzz swarm plot to plot peak usage for some of the top QTLs. I can use the examples I gave Tony.

Nuclear:
* peak305794, sid: 7:128635754

• peak: 164036, sid: 2:3502035

Total:

• Peak: peak228606, SID 3:150302010

• Peak: peak152751, SID 19:4236475

I need to pull out the genotypes for each snp and the corresponding phenotype. I want to make a python script that I can give a snp and a peak and it will make a table with the genotypes and phenotypes for the necessary gene snp pair.

geno3_m=fread("../data/apaExamp/geno3_150302010.txt", header = T) %>% dplyr::select(starts_with("NA")) %>% t
geno3df= data.frame(geno3_m) %>% separate(geno3_m, into=c("geno", "dose", "extra"), sep=":") %>% dplyr::select(dose) %>% rownames_to_column(var="ind")
apaphen228606_m= fread("../data/apaExamp/Total.peak228606.txt", header = T) %>% dplyr::select(starts_with("NA")) %>% t
apaphen228606_df=data.frame(apaphen228606_m) %>% rownames_to_column(var="ind")

toplotAPA=geno3df %>% inner_join(apaphen228606_df, by="ind")
toplotAPA$dose= as.factor(toplotAPA$dose)
colnames(toplotAPA)= c("ind", "Genotype", "APA")
EIF2A_APAex=ggplot(toplotAPA, aes(y=APA, x=Genotype, by=Genotype, fill=Genotype)) + geom_boxplot() + geom_jitter() + labs(y="APA phenotype", title="Total APA: Peak 228606, EIF2A")
ggsave("../output/plots/EIF2a_APA.png", EIF2A_APAex)

Saving 7 x 5 in image

This is in the gene EIF2A, I need to find this in the eQTL data. The ensg id for this gene is ENSG00000144895.

RNAseqEIF2A_m=read.table("../data/apaExamp/RNAseq.phenoEIF2A.txt", header=T) %>% dplyr::select(starts_with("NA")) %>% t 
RNAseqEIF2A_df= data.frame(RNAseqEIF2A_m) %>% rownames_to_column("ind")

plotRNA=geno3df %>% inner_join(RNAseqEIF2A_df, by="ind")
plotRNA$dose= as.factor(plotRNA$dose)
colnames(plotRNA)= c("ind", "Genotype", "Expression")

EIF2A_RNAex=ggplot(plotRNA, aes(y=Expression, x=Genotype, by=Genotype, fill=Genotype)) + geom_boxplot() + geom_jitter() + labs(y="Normalized Expression", title="Gene Expression: EIF2A")

ggsave("../output/plots/EIF2a_RNA.png", EIF2A_RNAex)

Saving 7 x 5 in image

Try this in protein:

ProtEIF2A_m=read.table("../data/apaExamp/ProtEIF2A.txt", header=T) %>% dplyr::select(starts_with("NA")) %>% t 
ProtEIF2A_df= data.frame(ProtEIF2A_m) %>% rownames_to_column("ind")

plotProt=geno3df %>% inner_join(ProtEIF2A_df, by="ind")
plotProt$dose= as.factor(plotProt$dose)
colnames(plotProt)= c("ind", "Genotype", "Prot_level")

IF2A_Protex= ggplot(plotProt, aes(y=Prot_level, x=Genotype, by=Genotype, fill=Genotype)) + geom_boxplot() + geom_jitter() + labs(y="Normalized Protein Level", title="Protein Level: EIF2A")

ggsave("../output/plots/EIF2a_Prot.png", IF2A_Protex)

Saving 7 x 5 in image

multphenoEIF2a=plot_grid(EIF2A_APAex,IF2A_Protex,EIF2A_RNAex,nrow=1)
ggsave("../output/plots/EIF2a_multpheno.png", multphenoEIF2a, width=15, height=5)

Do this with 4su 60:

have to remove the #

su60_EIF2A_m=read.table("../data/apaExamp/Foursu60EIF2A.txt", header=T) %>% dplyr::select(starts_with("NA")) %>% t 
su60_EIF2A_df= data.frame(su60_EIF2A_m) %>% rownames_to_column("ind")

plot4su60=geno3df %>% inner_join(su60_EIF2A_df, by="ind")
plot4su60$dose= as.factor(plot4su60$dose)
colnames(plot4su60)= c("ind", "Genotype", "su60")

EIF2A_4su60ex=ggplot(plot4su60, aes(y=su60, x=Genotype, by=Genotype, fill=Genotype)) + geom_boxplot() + geom_jitter() + labs(y="4su60", title="4su 60min Value: EIF2A")

ggsave("../output/plots/EIF2a_4su60.png", EIF2A_4su60ex)

Saving 7 x 5 in image

Geuvadis RNA

rnaG_EIF2A_m=read.table("../data/apaExamp/RNA_GEU_EIF2A.txt", header=T) %>% dplyr::select(starts_with("NA")) %>% t 
rnaG_EIF2A_df= data.frame(rnaG_EIF2A_m) %>% rownames_to_column("ind")

plotRNAg=geno3df %>% inner_join(rnaG_EIF2A_df, by="ind")
plotRNAg$dose= as.factor(plotRNAg$dose)
colnames(plotRNAg)= c("ind", "Genotype", "RNAg")

EIF2A_RNAgex=ggplot(plotRNAg, aes(y=RNAg, x=Genotype, by=Genotype, fill=Genotype)) + geom_boxplot() + geom_jitter() + labs(y="RNA expression", title="RNA Expression Geuvadis: EIF2A")

ggsave("../output/plots/EIF2a_RNAg.png", EIF2A_RNAgex)

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

other attached packages:
 [1] cowplot_0.9.3       data.table_1.11.8   VennDiagram_1.6.20 
 [4] futile.logger_1.4.3 forcats_0.3.0       stringr_1.3.1      
 [7] dplyr_0.7.6         purrr_0.2.5         readr_1.1.1        
[10] tidyr_0.8.1         tibble_1.4.2        ggplot2_3.0.0      
[13] 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] lambda.r_1.2.3       modelr_0.1.2         readxl_1.1.0        
[16] bindrcpp_0.2.2       bindr_0.1.1          plyr_1.8.4          
[19] munsell_0.5.0        gtable_0.2.0         cellranger_1.1.0    
[22] rvest_0.3.2          R.methodsS3_1.7.1    evaluate_0.11       
[25] labeling_0.3         knitr_1.20           broom_0.5.0         
[28] Rcpp_0.12.19         formatR_1.5          backports_1.1.2     
[31] scales_1.0.0         jsonlite_1.5         hms_0.4.2           
[34] digest_0.6.17        stringi_1.2.4        rprojroot_1.3-2     
[37] cli_1.0.1            tools_3.5.1          magrittr_1.5        
[40] lazyeval_0.2.1       futile.options_1.0.1 crayon_1.3.4        
[43] whisker_0.3-2        pkgconfig_2.0.2      xml2_1.2.0          
[46] lubridate_1.7.4      assertthat_0.2.0     rmarkdown_1.10      
[49] httr_1.3.1           rstudioapi_0.8       R6_2.3.0            
[52] nlme_3.1-137         git2r_0.23.0         compiler_3.5.1