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")) + scale_fill_brewer(palette="Paired")

Scale for 'fill' is already present. Adding another scale for 'fill',
which will replace the existing scale.

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

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")) + scale_fill_brewer(palette="Paired")

Scale for 'fill' is already present. Adding another scale for 'fill',
which will replace the existing scale.

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")) + scale_fill_brewer(palette="Paired")

Scale for 'fill' is already present. Adding another scale for 'fill',
which will replace the existing scale.

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) + scale_fill_brewer(palette="Paired")

Scale for 'fill' is already present. Adding another scale for 'fill',
which will replace the existing scale.

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)+ scale_fill_brewer(palette="Paired")

Scale for 'fill' is already present. Adding another scale for 'fill',
which will replace the existing scale.

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)+ scale_fill_brewer(palette="Paired")

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"))+ scale_fill_brewer(palette="Paired")

Scale for 'fill' is already present. Adding another scale for 'fill',
which will replace the existing scale.

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"))+ scale_fill_brewer(palette="Paired")

Scale for 'fill' is already present. Adding another scale for 'fill',
which will replace the existing scale.

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

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:]
    name=name.split("_")[0]
    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 also need to create a significant QTL snp bedfile for the total qtls. Bed files are 0 bases meaning I want the end to be the position I care about.

chrom, start (pos -1), end (pos), name, score (bh), strand

I can do this in python using the /project2/gilad/briana/threeprimeseq/data/perm_APAqtl_trans/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total_transcript_permResBH.txt. I will make the script general to use on the total or nuclear file.

QTLres2SigSNPbed.py

def main(inFile, outFile):
    fout=open(outFile, "w")
    fin=open(inFile, "r")
    for num, ln in enumerate(fin):
      if num >= 1:
          pid, nvar, shape1, shape2, dummy, sid, dist, npval, slope, ppval, bpval, bh = ln.split()
          chrom, pos= sid.split(":")
          name=sid
          start= int(pos)-1
          end=int(pos)
          strand=pid.split(":")[3].split("_")[1]
          bh=float(bh)
          if bh <= .1: 
              fout.write("%s\t%s\t%s\t%s\t%s\t%s\n"%(chrom, start, end, name, bh, strand))
    fout.close()

if __name__ == "__main__":
    import sys
    fraction=sys.argv[1]
    inFile = "/project2/gilad/briana/threeprimeseq/data/perm_APAqtl_trans/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_%s_transcript_permResBH.txt"%(fraction)
    outFile= "/project2/gilad/briana/threeprimeseq/data/perm_APAqtl_trans/sigSnps/ApaQTLsignificantSnps_10percFDR_%s.bed"%(fraction)
    main(inFile,outFile) 

I am going to try to use pybedtools instead for bedtools for this analysis. First I can add it to my conda environment.

Remove header from HMM

import pybedtools 

sigNuc=pybedtools.BedTool('/project2/gilad/briana/threeprimeseq/data/perm_APAqtl_trans/sigSnps/ApaQTLsignificantSnps_10percFDR_Nuclear.sort.bed') 

sigTot=pybedtools.BedTool('/project2/gilad/briana/threeprimeseq/data/perm_APAqtl_trans/sigSnps/ApaQTLsignificantSnps_10percFDR_Total.sort.bed')

hmm=pybedtools.BedTool("/project2/gilad/briana/genome_anotation_data/GM12878.chromHMM.sort.bed")

#map hmm to snps  
Tot_overlapHMM=sigTot.map(hmm, c=4)

Nuc_overlapHMM=sigNuc.map(hmm,c=4)

#save results  

Tot_overlapHMM.saveas("/project2/gilad/briana/threeprimeseq/data/perm_APAqtl_trans/sigSnps/Tot_overlapHMM.bed")

Nuc_overlapHMM.saveas("/project2/gilad/briana/threeprimeseq/data/perm_APAqtl_trans/sigSnps/Nuc_overlapHMM.bed")

I want to make a file that has all of the numbers for the chromatin regions for downstream analysis.

cut -f5 GM12878.chromHMM.txt | sort | uniq > chromHMM_regions.txt

I then manually seperate the numbers from the name with a tab and remove the name line.

Evaluate results for total:

chromHmm=read.table("../data/ChromHmmOverlap/chromHMM_regions.txt", col.names = c("number", "name"), stringsAsFactors = F)

TotalOverlapHMM=read.table("../data/ChromHmmOverlap/Tot_overlapHMM.bed", col.names=c("chrom", "start", "end", "sid", "significance", "strand", "number"))

TotalOverlapHMM_names=TotalOverlapHMM %>% left_join(chromHmm, by="number")

ggplot(TotalOverlapHMM_names, aes(x=name)) + geom_bar() + labs(title="ChromHMM labels for Total APAQtls" , y="Number of SNPs", x="Region")+theme(axis.text.x = element_text(angle = 90, hjust = 1))

This is the count but I want enrichemnt. I need to randomly choose 188 snps from the ones I tested (nominal res) and get the same inforamtion on where they are.

shuf -n 118 /project2/gilad/briana/threeprimeseq/data/nominal_APAqtl_trans/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total_NomRes.txt > /project2/gilad/briana/threeprimeseq/data/nominal_APAqtl_trans/randomSnps/ApaQTL_total_Random118.txt

Now I need to make these into the snp bed format.

QTLNOMres2SigSNPbed.py give this * total or nuclear
* number of snps

def main(inFile, outFile):
    fout=open(outFile, "w")
    fin=open(inFile, "r")
    for ln in fin:
          pid, sid, dist, pval, slope = ln.split()
          chrom, pos= sid.split(":")
          name=sid
          start= int(pos)-1
          end=int(pos)
          strand=pid.split(":")[3].split("_")[1]
          pval=float(pval)
          fout.write("%s\t%s\t%s\t%s\t%s\t%s\n"%(chrom, start, end, name, pval, strand))
    fout.close()

if __name__ == "__main__":
    import sys
    fraction=sys.argv[1]
    number=sys.argv[2]
    inFile = "/project2/gilad/briana/threeprimeseq/data/nominal_APAqtl_trans/randomSnps/ApaQTL_%s_Random%s.txt"%(fraction,number)
    outFile= "/project2/gilad/briana/threeprimeseq/data/nominal_APAqtl_trans/randomSnps/ApaQTL_%s_Random%s.bed"%(fraction,number)
    main(inFile,outFile) 

Sort output

import pybedtools 

RANDtot=pybedtools.BedTool('/project2/gilad/briana/threeprimeseq/data/nominal_APAqtl_trans/randomSnps/ApaQTL_total_Random118.sort.bed') 



hmm=pybedtools.BedTool("/project2/gilad/briana/genome_anotation_data/GM12878.chromHMM.sort.bed")

#map hmm to snps  
TotRnad_overlapHMM=RANDtot.map(hmm, c=4)


#save results  

TotRnad_overlapHMM.saveas("/project2/gilad/briana/threeprimeseq/data/nominal_APAqtl_trans/randomSnps/ApaQTL_total_Random_overlapHMM.bed")

TotalRandOverlapHMM=read.table("../data/ChromHmmOverlap/ApaQTL_total_Random_overlapHMM.bed", col.names=c("chrom", "start", "end", "sid", "significance", "strand", "number"))

TotalRandOverlapHMM_names=TotalRandOverlapHMM %>% left_join(chromHmm, by="number")

ggplot(TotalRandOverlapHMM_names, aes(x=name)) + geom_bar() + labs(title="ChromHMM labels for Total APAQtls (Random)" , y="Number of SNPs", x="Region")+theme(axis.text.x = element_text(angle = 90, hjust = 1))

To put this on the same plot I can count the number in each then plot them next to eachother.

random_perChromHMM=TotalRandOverlapHMM_names %>%  group_by(name) %>% summarise(Random=n())
sig_perChromHMM= TotalOverlapHMM_names %>%  group_by(name) %>%  summarise(Total_QTLs=n())

perChrommHMM=random_perChromHMM %>%  full_join(sig_perChromHMM, by="name", ) %>% replace_na(list(Random=0,Total_QTLs=0))  

perChrommHMM_melt=melt(perChrommHMM, id.vars="name")
names(perChrommHMM_melt)=c("Region","Set", "N_Snps" )

chromenrichTotalplot=ggplot(perChrommHMM_melt, aes(x=Region, y=N_Snps, by=Set, fill=Set)) + geom_bar(position="dodge", stat="identity") +theme(axis.text.x = element_text(angle = 90, hjust = 1)) + labs(title="Enrichment of Total QTLs by chromatin region", y="Number of Snps", x="Chromatin Region") + scale_fill_brewer(palette="Paired")
chromenrichTotalplot

ggsave("../output/plots/ChromHmmEnrich_Total.png", chromenrichTotalplot)

Saving 7 x 5 in image

I want to make a plot with the enrichment by fraction. I am first going to get an enrichemnt score for each bin naively by looking at the QTL-random in each category.

perChrommHMM$Random= as.integer(perChrommHMM$Random)
perChrommHMM$Total_QTLs= as.integer(perChrommHMM$Total_QTLs)
perChrommHMM_enr=perChrommHMM %>%  mutate(Total=Total_QTLs- Random)

Write this file so I can put it in the nuclear analysis and compare between groups.

write.table(perChrommHMM_enr, file="../data/ChromHmmOverlap/perChrommHMM_Total_enr.txt", quote=F, sep="\t", col.names = T, row.names = F)

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] RColorBrewer_1.1-2   lambda.r_1.2.3       modelr_0.1.2        
[16] readxl_1.1.0         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