Characterize Total ApaQTLs

Upload Libraries and 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:

I will add a column to this dataframe that will tell me if the association is significant at 10% FDR. This will help me plot based on significance later in the analysis. I am also going to seperate the PID into relevant pieces.

totalAPA=read.table("../data/perm_QTL_trans/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total_transcript_permResBH.txt", stringsAsFactors = F, header=T)  %>% mutate(sig=ifelse(-log10(bh)>=1, 1,0 )) %>%  separate(pid, sep = ":", into=c("chr", "start", "end", "id")) %>% separate(id, sep = "_", into=c("gene", "strand", "peak"))

totalAPA$sig=as.factor(totalAPA$sig)


print(names(totalAPA))

 [1] "chr"    "start"  "end"    "gene"   "strand" "peak"   "nvar"  
 [8] "shape1" "shape2" "dummy"  "sid"    "dist"   "npval"  "slope" 
[15] "ppval"  "bpval"  "bh"     "sig"   

Distance Metrics

ggplot(totalAPA, aes(x=dist, fill=sig, by=sig)) + geom_density(alpha=.5)  +  labs(title="Distance from snp to TSS", x="Base Pairs") + scale_fill_discrete(guide = guide_legend(title = "Significant QTL"))

ggplot(totalAPA, aes(x=dist, fill=sig, by=sig)) + geom_density(alpha=.5)+coord_cartesian(xlim = c(-150000, 150000))

peaks=read.table("../data/PeaksUsed/filtered_APApeaks_merged_allchrom_refseqTrans.noties_sm.fixed.bed", col.names = c("chr", "peakStart", "peakEnd", "PeakNum", "PeakScore", "Strand", "Gene")) %>% mutate(peak=paste("peak", PeakNum,sep="")) %>% mutate(PeakCenter=peakStart+ (peakEnd- peakStart))

totalAPA_peakdist= totalAPA %>%  inner_join(peaks, by="peak") %>%  separate(sid, into=c("snpCHR", "snpLOC"), by=":")
totalAPA_peakdist$snpLOC= as.numeric(totalAPA_peakdist$snpLOC)

totalAPA_peakdist= totalAPA_peakdist %>%  mutate(DisttoPeak= snpLOC-PeakCenter)

ggplot(totalAPA_peakdist, aes(x=DisttoPeak, fill=sig, by=sig)) + geom_density(alpha=.5)  +  labs(title="Distance from snp peak", x="log10 absolute value Distance to Peak") + scale_fill_discrete(guide = guide_legend(title = "Significant QTL"))

totalAPA_peakdist_sig=totalAPA_peakdist %>% filter(sig==1)
totalAPA_peakdist_notsig=totalAPA_peakdist %>% filter(sig==0)


summary(totalAPA_peakdist_sig$DisttoPeak)

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-634474.0  -26505.0   -3325.5  -23883.8     492.5  460051.0 

summary(totalAPA_peakdist_notsig$DisttoPeak)

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
-70147526   -269853     -1269     -6620    266370   5433999 

ggplot(totalAPA_peakdist, aes(y=DisttoPeak,x=sig, fill=sig, by=sig)) + geom_boxplot()  + scale_fill_discrete(guide = guide_legend(title = "Significant QTL"))

totalAPA_peakdist_filt=totalAPA_peakdist %>% filter(abs(DisttoPeak) <= 1*(10^6))

ggplot(totalAPA_peakdist_filt, aes(y=DisttoPeak,x=sig, fill=sig, by=sig)) + geom_boxplot()  + scale_fill_discrete(guide = guide_legend(title = "Significant QTL")) + facet_grid(.~strand)

ggplot(totalAPA_peakdist_filt, aes(x=DisttoPeak, fill=sig, by=sig)) + geom_density()  + scale_fill_discrete(guide = guide_legend(title = "Significant QTL")) + facet_grid(.~strand)

totalAPA_peakdist %>% filter(sig==1) %>% filter(strand=="+") %>%  filter(DisttoPeak < 0) %>% nrow()

[1] 45

totalAPA_peakdist %>% filter(sig==1) %>% filter(strand=="-") %>%  filter(DisttoPeak > 0) %>% nrow()

[1] 16

totalAPA_peakdist %>% filter(sig==1) %>% filter(strand=="+") %>%  filter(DisttoPeak > 0) %>% nrow()

[1] 17

totalAPA_peakdist %>% filter(sig==1) %>% filter(strand=="-") %>%  filter(DisttoPeak < 0) %>% nrow()

[1] 40

ggplot(totalAPA_peakdist_sig, aes(x=DisttoPeak)) + geom_density()

totalAPA_peakdist_sig %>% filter(abs(DisttoPeak) < 1000) %>% nrow()

[1] 29

totalAPA_peakdist_sig %>% filter(abs(DisttoPeak) < 10000) %>% nrow()

[1] 57

totalAPA_peakdist_sig %>% filter(abs(DisttoPeak) < 100000) %>% nrow()

[1] 98

Expression metrics

Next I want to pull in the expression values and compare the expression of genes with a total APA qtl in comparison to genes without one. I will also need to pull in the gene names file to add in the gene names from the ensg ID.

Remove the # from the file.

expression=read.table("../data/mol_pheno/fastqtl_qqnorm_RNAseq_phase2.fixed.noChr.txt", header = T,stringsAsFactors = F)
expression_mean=apply(expression[,5:73],1,mean,na.rm=TRUE)
expression_var=apply(expression[,5:73],1,var,na.rm=TRUE)
expression$exp.mean= expression_mean 
expression$exp.var=expression_var
expression= expression %>% separate(ID, into=c("Gene.stable.ID", "ver"), sep ="[.]")

Now I can pull in the names and join the dataframes.

geneNames=read.table("../data/ensemble_to_genename.txt", sep="\t", header=T,stringsAsFactors = F) 



expression=expression %>% inner_join(geneNames,by="Gene.stable.ID") 

expression=expression %>% select(Chr, start, end, Gene.name, exp.mean,exp.var) %>%  rename("gene"=Gene.name)

Now I can join this with the qtls.

totalAPA_wExp=totalAPA %>% inner_join(expression, by="gene") 

ggplot(totalAPA_wExp, aes(x=exp.mean, by=sig, fill=sig)) + geom_density(alpha=.3)

This is not exactly what I want because there are multiple peaks in a gene so some genes are plotted multiple times. I want to group the QTLs by gene and see if there is 1 sig QTL for that gene.

gene_wQTL= totalAPA_wExp %>% group_by(gene) %>% summarise(sig_gene=sum(as.numeric(as.character(sig)))) %>% ungroup() %>% inner_join(expression, by="gene") %>% mutate(sigGeneFactor=ifelse(sig_gene>=1, 1,0))

gene_wQTL$sigGeneFactor= as.factor(as.character(gene_wQTL$sigGeneFactor))

Therea are 92 genes in this set with a QTL.

ggplot(gene_wQTL, aes(x=exp.mean, by=sigGeneFactor, fill=sigGeneFactor)) + geom_density(alpha=.3) +labs(title="Mean in RNA expression by genes with significant QTL", x="Mean in normalized expression") + scale_fill_discrete(guide = guide_legend(title = "Significant QTL"))

I can do a similar analysis but test the variance in the gene expression.

ggplot(gene_wQTL, aes(x=exp.var, by=sigGeneFactor, fill=sigGeneFactor)) + geom_density(alpha=.3) + labs(title="Varriance in RNA expression by genes with significant QTL", x="Variance in normalized expression") + scale_fill_discrete(guide = guide_legend(title = "Significant QTL"))

ggplot(totalAPA_peakdist, aes(x=PeakScore,fill=sig,by=sig)) + geom_density(alpha=.5)+ scale_x_log10() + labs(title="Peak score by significance") 

Where are the snps:

Download the GM12878 chromHMM annotation. I downleaded this from uscs and put it in:

• /Users/bmittleman1/Documents/Gilad_lab/threeprimeseq/data/GM12878.chromHMM.txt
  
• /project2/gilad/briana/genome_anotation_data/GM12878.chromHMM.txt

Column:

• bin
  
• chrom
  
• chromstart
  
• chromend
  
• name
  
• score
  
• strand
  
• thich start
  
• thick end
  
• item rgb

I can make this a bedfile to use a bedtools pipeline:

• chrom (nochr)

• start

• end

• name (txn hetero ect)

• score

• strand

fout = open("/project2/gilad/briana/genome_anotation_data/GM12878.chromHMM.bed",'w')
for ln in open("/project2/gilad/briana/genome_anotation_data/GM12878.chromHMM.txt", "r"):
    bin, chrom, start, end, name, score, strand, thSt, thE, rgb = ln.split()
    chrom=chrom[3:]
    fout.write("%s\t%s\t%s\t%s\t%s\t%s\n"%(chrom, start, end, name, score, strand))
fout.close()

fout = open("/Users/bmittleman1/Documents/Gilad_lab/threeprimeseq/data/GM12878.chromHMM.bed",'w')
for ln in open("/Users/bmittleman1/Documents/Gilad_lab/threeprimeseq/data/GM12878.chromHMM.txt", "r"):
    bin, chrom, start, end, name, score, strand, thSt, thE, rgb = ln.split()
    chrom=chrom[3:]
    fout.write("%s\t%s\t%s\t%s\t%s\t%s\n"%(chrom, start, end, name, score, strand))
fout.close()

I am gonig to try to use pybedtooksinstead for bedtools for this analysis. First I can add it to my conda environment.

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] bindrcpp_0.2.2      cowplot_0.9.3       data.table_1.11.8  
 [4] VennDiagram_1.6.20  futile.logger_1.4.3 forcats_0.3.0      
 [7] stringr_1.3.1       dplyr_0.7.6         purrr_0.2.5        
[10] readr_1.1.1         tidyr_0.8.1         tibble_1.4.2       
[13] ggplot2_3.0.0       tidyverse_1.2.1     reshape2_1.4.3     
[16] 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] 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] R.methodsS3_1.7.1    evaluate_0.11        labeling_0.3        
[25] knitr_1.20           broom_0.5.0          Rcpp_0.12.19        
[28] formatR_1.5          backports_1.1.2      scales_1.0.0        
[31] jsonlite_1.5         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          magrittr_1.5         lazyeval_0.2.1      
[40] futile.options_1.0.1 crayon_1.3.4         whisker_0.3-2       
[43] pkgconfig_2.0.2      xml2_1.2.0           lubridate_1.7.4     
[46] assertthat_0.2.0     rmarkdown_1.10       httr_1.3.1          
[49] rstudioapi_0.8       R6_2.3.0             nlme_3.1-137        
[52] git2r_0.23.0         compiler_3.5.1