Peak To Gene Assignment

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

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

library(cowplot)

Attaching package: 'cowplot'

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

    ggsave

Prepare referece

#!/bin/bash

#SBATCH --job-name=intGenes_combfilterPeaksOppStrand
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=intGenes_combfilterPeaksOppStrand.out
#SBATCH --error=intGenes_combfilterPeaksOppStrand.err
#SBATCH --partition=broadwl
#SBATCH --mem=12G
#SBATCH --mail-type=END

module load Anaconda3
source activate three-prime-env

bedtools intersect -wa -wb -sorted -S -a /project2/gilad/briana/threeprimeseq/data/mergedPeaks_comb/filtered_APApeaks_merged_allchrom.named.fixed.bed -b /project2/gilad/briana/genome_anotation_data/ncbiRefSeq_sm_noChr.sort.mRNA.bed > /project2/gilad/briana/threeprimeseq/data/mergedPeaks_comb/filtered_APApeaks_merged_allchrom_refseqGenes.OppStrand.bed

This results in peaks being mapped to multiple genes. I want to use a method where I look for the closest end of transcript to each peak then use that gene for the assignment. This would mean each peak is assigned to one gene.

Create a python script to process the NCBI file. I want protien coding transcript ends with the associated gene names. Original file: ncbiRefSeq.txt

• Column 2 transcript name
• Column 13 gene name
• NM is protein coding

EndOfProCodTrans.py

def main(inF, outF):
  infile= open(inF, "r")
  fout = open(outF,'w')
  for line in infile:
      linelist=line.split()
      transcript=linelist[1]
      transcript_id=transcript.split("_")[0]
      if transcript_id=="NM":
          chr=linelist[2][3:]
          strand=linelist[3] 
          gene= linelist[12]
          if strand == "+" :
              end = int(linelist[7])
              end2= end - 1
              fout.write("%s\t%d\t%d\t%s:%s\t.\t%s\n"%(chr, end2, end, transcript,gene, strand))
          if strand == "-":
              end= int(linelist[4])
              end2= end + 1
              fout.write("%s\t%d\t%d\t%s:%s\t.\t%s\n"%(chr, end, end2, transcript,gene, strand))


if __name__ == "__main__":
    inF = "/project2/gilad/briana/genome_anotation_data/ncbiRefSeq.txt"
    outF= "/project2/gilad/briana/genome_anotation_data/ncbiRefSeq_endProtCodGenes.txt"
    main(inF, outF)

Find closest gene to each peak

bedtools closest

-A peaks /project2/gilad/briana/threeprimeseq/data/mergedPeaks_comb/filtered_APApeaks_merged_allchrom.named.fixed.bed -B transcript file /project2/gilad/briana/genome_anotation_data/ncbiRefSeq_endProtCodGenes_sort.txt -S (opposite strand) -D b (give distance wrt to gene strand)

#!/bin/bash

#SBATCH --job-name=TransClosest2End
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=TransClosest2End.out
#SBATCH --error=TransClosest2End.err
#SBATCH --partition=broadwl
#SBATCH --mem=12G
#SBATCH --mail-type=END

module load Anaconda3 

source activate three-prime-env


bedtools closest -S -D b -a /project2/gilad/briana/threeprimeseq/data/mergedPeaks_comb/filtered_APApeaks_merged_allchrom.named.fixed.bed -b  /project2/gilad/briana/genome_anotation_data/ncbiRefSeq_endProtCodGenes_sort.txt > /project2/gilad/briana/threeprimeseq/data/mergedPeaks_comb/filtered_APApeaks_merged_allchrom_refseqTrans.closest2End.bed

I will take a look at this file in R then I will process the file in python.

names=c("PeakChr", "PeakStart", "PeakEnd", "PeakName","PeakScore", "PeakStrand", "GeneChr", "GeneStart", "GeneEnd", "Transcript", "GeneScore", "GeneStrand", "Distance" )
peak2transDist=read.table("../data/filtered_APApeaks_merged_allchrom_refseqTrans.closest2End.bed", col.names = names, stringsAsFactors = F, header=F)

ggplot(peak2transDist, aes(x=abs(Distance)))+ geom_density() + scale_x_log10()

Warning: Transformation introduced infinite values in continuous x-axis

Warning: Removed 4362 rows containing non-finite values (stat_density).

Expand here to see past versions of unnamed-chunk-6-1.png:

Version	Author	Date
aaed5fd	Briana Mittleman	2018-09-26

peak2transDist0=peak2transDist %>% filter(Distance==0)
nrow(peak2transDist0)

[1] 4362

peak2transDist200=peak2transDist %>% filter(abs(Distance)<200)
nrow(peak2transDist200)

[1] 23778

summary(peak2transDist$Distance)

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-5523243   -57698   -12830   -23711     3373  5592124 

try adding the no ties flag -t first.

peak2transDist_noties=read.table("../data/filtered_APApeaks_merged_allchrom_refseqTrans.closest2End.noties.bed", col.names = names, stringsAsFactors = F, header=F)
ggplot(peak2transDist_noties, aes(x=abs(Distance)))+ geom_density() + scale_x_log10()

Warning: Transformation introduced infinite values in continuous x-axis

Warning: Removed 2044 rows containing non-finite values (stat_density).

Expand here to see past versions of unnamed-chunk-8-1.png:

Version	Author	Date
1501499	Briana Mittleman	2018-09-26

peak2transDist0_noT=peak2transDist_noties%>% filter(Distance==0)
nrow(peak2transDist0_noT)

[1] 2044

peak2transDist200_noT=peak2transDist_noties %>% filter(abs(Distance)<200)
nrow(peak2transDist200_noT)

[1] 10488

summary(peak2transDist$Distance)

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-5523243   -57698   -12830   -23711     3373  5592124 

ggplot(peak2transDist_noties, aes(x=abs(Distance)))+ geom_histogram(binwidth = .5) + scale_x_log10()

Warning: Transformation introduced infinite values in continuous x-axis

Warning: Removed 2044 rows containing non-finite values (stat_bin).

Expand here to see past versions of unnamed-chunk-9-1.png:

Version	Author	Date
1501499	Briana Mittleman	2018-09-26

Looking at this visually suggests that we have way too many peaks. I want to compare the peak score which is related to the coverage to the abs(distace)

ggplot(peak2transDist_noties, aes(y=PeakScore, x=abs(Distance + 1))) + geom_point() + scale_x_log10() + scale_y_log10() + geom_density2d(na.rm = TRUE, size = 1, colour = 'red') 

Warning: Transformation introduced infinite values in continuous x-axis

Warning: Transformation introduced infinite values in continuous x-axis

Alternatively let me try to remove low peak score values.

allPeakplot=ggplot(peak2transDist_noties, aes(x=abs(Distance + 1)))+ geom_density() + scale_x_log10() + labs(title="Distance all peaks to gene end") + annotate("text", label=nrow(peak2transDist_noties), x=10, y=.4)

peak2transDist_score500=peak2transDist_noties%>% filter(PeakScore>500)
score500plot=ggplot(peak2transDist_score500, aes(x=abs(Distance + 1)))+ geom_density() + scale_x_log10() + labs(title="Peak Score > 500") + annotate("text", label=nrow(peak2transDist_score500), x=10, y=.4)


peak2transDist_score200=peak2transDist_noties%>% filter(PeakScore>200)
score200plot=ggplot(peak2transDist_score200, aes(x=abs(Distance + 1)))+ geom_density() + scale_x_log10() + labs(title="Peak Score > 200") + annotate("text", label=nrow(peak2transDist_score200), x=10, y=.4)


peak2transDist_score100=peak2transDist_noties%>% filter(PeakScore>100)
score100plot=ggplot(peak2transDist_score100, aes(x=abs(Distance + 1)))+ geom_density() + scale_x_log10() + labs(title="Peak Score > 100") + annotate("text", label=nrow(peak2transDist_score100), x=10, y=.4)

peak2transDist_score50=peak2transDist_noties%>% filter(PeakScore>50)
score50plot=ggplot(peak2transDist_score50, aes(x=abs(Distance + 1)))+ geom_density() + scale_x_log10() + labs(title="Peak Score > 50")+ annotate("text", label=nrow(peak2transDist_score50), x=10, y=.4)

peak2transDist_score20=peak2transDist_noties%>% filter(PeakScore>20)
score20plot=ggplot(peak2transDist_score20, aes(x=abs(Distance + 1)))+ geom_density() + scale_x_log10() + labs(title="Peak Score > 10")+ annotate("text", label=nrow(peak2transDist_score20), x=10, y=.4)

plot_grid(allPeakplot,score20plot,score50plot,score100plot,score200plot, score500plot)

Warning: Transformation introduced infinite values in continuous x-axis

Warning: Removed 662 rows containing non-finite values (stat_density).

Warning: Transformation introduced infinite values in continuous x-axis

Warning: Removed 431 rows containing non-finite values (stat_density).

Warning: Transformation introduced infinite values in continuous x-axis

Warning: Removed 327 rows containing non-finite values (stat_density).

Warning: Transformation introduced infinite values in continuous x-axis

Warning: Removed 234 rows containing non-finite values (stat_density).

Warning: Transformation introduced infinite values in continuous x-axis

Warning: Removed 150 rows containing non-finite values (stat_density).

Warning: Transformation introduced infinite values in continuous x-axis

Warning: Removed 78 rows containing non-finite values (stat_density).

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  cowplot_0.9.3   workflowr_1.1.1 forcats_0.3.0  
 [5] stringr_1.3.1   dplyr_0.7.6     purrr_0.2.5     readr_1.1.1    
 [9] tidyr_0.8.1     tibble_1.4.2    ggplot2_3.0.0   tidyverse_1.2.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     labeling_0.3      knitr_1.20       
[25] broom_0.5.0       Rcpp_0.12.18      scales_1.0.0     
[28] backports_1.1.2   jsonlite_1.5      hms_0.4.2        
[31] digest_0.6.16     stringi_1.2.4     grid_3.5.1       
[34] rprojroot_1.3-2   cli_1.0.0         tools_3.5.1      
[37] magrittr_1.5      lazyeval_0.2.1    crayon_1.3.4     
[40] whisker_0.3-2     pkgconfig_2.0.2   MASS_7.3-50      
[43] xml2_1.2.0        lubridate_1.7.4   assertthat_0.2.0 
[46] rmarkdown_1.10    httr_1.3.1        rstudioapi_0.7   
[49] R6_2.2.2          nlme_3.1-137      git2r_0.23.0     
[52] compiler_3.5.1