eRNAs and Repetitive elements

overlap potential drivers of extra peaks

Does not look like there are a strong driver. I will see if any of these overlap with our peaks.I will need to look at the opposite strand overlap or use the fixed strand peaks. I will ask how many of these eRNAs or rep elements overlap a peak.

I want to run the overlap in all of the peaks as well as those that have been filtered 5%

/project2/gilad/briana/threeprimeseq/data/mergedPeaks_noMP_GeneLoc/Filtered_APApeaks_merged_allchrom_noMP.sort.named.noCHR_geneLoc.bed

/project2/gilad/briana/threeprimeseq/data/peaks4DT/Peaks_5percCov_fixedStrand.bed

fix strand for nonfiltered:

fixStrand4DTplots_allpeaks.py

peaksIn="/project2/gilad/briana/threeprimeseq/data/mergedPeaks_noMP_GeneLoc/Filtered_APApeaks_merged_allchrom_noMP.sort.named.noCHR_geneLoc.bed"
PeakOut="/project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPeaks_allPeaks_fixedStrand.bed"


def fix_strand(Fin,Fout):
    fout=open(Fout,"w")
    for ln in open(Fin, "r"):
        chrom, start, end, name, score, strand, score2, pos = ln.split()
        if strand=="+":
            nameF="peak" + name + ":" + pos
            fout.write("%s\t%s\t%s\t%s\t%s\t-\n"%(chrom,start,end,nameF,score))
        else:
            nameF="peak" + name + ":" + pos
            fout.write("%s\t%s\t%s\t%s\t%s\t+\n"%(chrom,start,end,nameF,score))
    fout.close()
    
    
fix_strand(peaksIn, PeakOut)

ernas: /project2/gilad/briana/genome_anotation_data/LCLenhancerRNA.bed repetitive elements: /project2/gilad/briana/genome_anotation_data/RepeatMask.bed

make a python script with pybedtools that will take any bed file and overlap it

overlapWFilteredPeaks.py

def main(infile, outfile):
    peak_file=open("/project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPeaks_5percCov_fixedStrand.bed","r")
    peak=pybedtools.BedTool(peak_file)
    elementFile=open(infile, "r")
    for i,ln in enumerate(elementFile):
       if i == 0:
           if len(ln.split()) > 3:
               strand= "yes"
           else:
               strand= "no"
       else:
          break
    print(strand)
    elements=pybedtools.BedTool(elementFile)
    if strand== "yes": 
        elemOverpeak=elements.intersect(peak, wa=True,wb=True, s=True)
    else:
        elemOverpeak=elements.intersect(peak, wa=True,wb=True)
    elemOverpeak.saveas(outfile)

if __name__ == "__main__":
    import sys
    import pybedtools
    infile=sys.argv[1]
    outfile=sys.argv[2]
    main(infile, outfile)

run:

python overlapWFilteredPeaks.py /project2/gilad/briana/genome_anotation_data/LCLenhancerRNA.bed  /project2/gilad/briana/threeprimeseq/data/FeatureoverlapPeaks/FilteredPeak_overeRNA.txt  

python overlapWFilteredPeaks.py /project2/gilad/briana/genome_anotation_data/RepeatMask.bed /project2/gilad/briana/threeprimeseq/data/FeatureoverlapPeaks/FilteredPeak_overRepElements.txt  

overlapWAllPeaks.py

def main(infile, outfile):
    peak_file=open("/project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPeaks_allPeaks_fixedStrand.bed","r")
    peak=pybedtools.BedTool(peak_file)
    elementFile=open(infile, "r")
    for i,ln in enumerate(elementFile):
       if i == 0:
           if len(ln.split()) > 3:
               strand= "yes"
           else:
               strand= "no"
       else:
          break
    print(strand)
    elements=pybedtools.BedTool(elementFile)
    if strand== "yes": 
        elemOverpeak=elements.intersect(peak, wa=True,wb=True, s=True)
    else:
        elemOverpeak=elements.intersect(peak, wa=True,wb=True)
    elemOverpeak.saveas(outfile)

if __name__ == "__main__":
    import sys
    import pybedtools
    infile=sys.argv[1]
    outfile=sys.argv[2]
    main(infile, outfile)
    

run:

python overlapWAllPeaks.py /project2/gilad/briana/genome_anotation_data/LCLenhancerRNA.bed  /project2/gilad/briana/threeprimeseq/data/FeatureoverlapPeaks/AllPeak_overeRNA.txt  

python overlapWAllPeaks.py /project2/gilad/briana/genome_anotation_data/RepeatMask.bed /project2/gilad/briana/threeprimeseq/data/FeatureoverlapPeaks/AllPeak_overRepElements.txt  

How long are each of these

Full eRNA file: 1167 All peak eRNA: 128 Filt peak eRNA: 14

Full rep file: 5298130 All peak rep:52965 Filt peak rep: 9542

For Repetitive seq. I can pull this in here and look at the proportion in each type.

repEl=read.table("../data/FeatureoverlapPeaks/FilteredPeak_overRepElements.txt", col.names = c("repCHR", "repStart", "repEnd", "Type", "length", "repStrand", "peakChr", "peakStart", "peakEnd", "peak", "score", "peakStrand"),stringsAsFactors = F)

repEl_sum= repEl %>% group_by(Type) %>% summarise(NperType=n()) %>% mutate(Proppertype=NperType/nrow(repEl))

ggplot(repEl_sum, aes(x=Type,y=Proppertype) ) + geom_bar(stat="identity") +  theme(axis.text.x = element_text(angle = 90, hjust = 1))

summary(repEl_sum$Proppertype)

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.0001048 0.0001048 0.0003144 0.0014025 0.0008384 0.0552295 

To get an expectation for this I will shuffle my peaks around the genome and overlap again.

Do this interactively

import pybedtools 
peaks= pybedtools.BedTool("/project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPeaks_5percCov_fixedStrand.bed")
peaks_shuf=peaks.shuffle(genome='hg19')
peaks_shuf.saveas("/project2/gilad/briana/threeprimeseq/data/FeatureoverlapPeaks/shuffled_FilterPeaks.bed")

sort -k1,1 -k2,2n /project2/gilad/briana/threeprimeseq/data/FeatureoverlapPeaks/shuffled_FilterPeaks.bed | sed 's/^chr//'  > /project2/gilad/briana/threeprimeseq/data/FeatureoverlapPeaks/shuffled_FilterPeaks.sort.bed

overlapWShuffledPeaks.py

def main(infile, outfile):
    peak_file=open("/project2/gilad/briana/threeprimeseq/data/FeatureoverlapPeaks/shuffled_FilterPeaks.sort.bed","r")
    peak=pybedtools.BedTool(peak_file)
    elementFile=open(infile, "r")
    for i,ln in enumerate(elementFile):
       if i == 0:
           if len(ln.split()) > 3:
               strand= "yes"
           else:
               strand= "no"
       else:
          break
    print(strand)
    elements=pybedtools.BedTool(elementFile)
    if strand== "yes": 
        elemOverpeak=elements.intersect(peak, wa=True,wb=True, s=True)
    else:
        elemOverpeak=elements.intersect(peak, wa=True,wb=True)
    elemOverpeak.saveas(outfile)

if __name__ == "__main__":
    import sys
    import pybedtools
    infile=sys.argv[1]
    outfile=sys.argv[2]
    main(infile, outfile)

python overlapWShuffledPeaks.py /project2/gilad/briana/genome_anotation_data/RepeatMask.bed /project2/gilad/briana/threeprimeseq/data/FeatureoverlapPeaks/ShuffledPeak_overRepElements.txt  

repEl_shuf=read.table("../data/FeatureoverlapPeaks/ShuffledPeak_overRepElements.txt", col.names = c("repCHR", "repStart", "repEnd", "Type", "length", "repStrand", "peakChr", "peakStart", "peakEnd", "peak", "score", "peakStrand"),stringsAsFactors = F)

repEl_shuf_sum= repEl_shuf %>% group_by(Type) %>% summarise(NperType_shuf=n()) %>% mutate(Proppertype_shuf=NperType_shuf/nrow(repEl_shuf)) 

Join both:

repEl_both=repEl_sum %>% full_join(repEl_shuf_sum, by="Type")

#fill NAs with 0  

repEl_both$Proppertype= repEl_both$Proppertype %>% replace_na(0)
repEl_both$NperType= repEl_both$NperType %>% replace_na(0)
repEl_both$Proppertype_shuf= repEl_both$Proppertype_shuf %>% replace_na(0)
repEl_both$NperType_shuf= repEl_both$NperType_shuf %>% replace_na(0)

Hyper geometric:

repEl_both_sig= repEl_both %>% mutate(Peak=sum(repEl_both$NperType), Peak_Shuf=sum(repEl_both$NperType_shuf), Chosen=NperType+NperType_shuf) %>% mutate(hyper=phyper(NperType, Peak, Peak_Shuf, Chosen, lower.tail = F)) %>% mutate(sig=ifelse(hyper<.05, "yes", "no")) %>% filter(sig=="yes") %>% arrange(desc(Proppertype))


repEl_both_sig

# A tibble: 227 x 10
   Type  NperType Proppertype NperType_shuf Proppertype_shuf  Peak
   <chr>    <dbl>       <dbl>         <dbl>            <dbl> <dbl>
 1 AT_r…      527      0.0552           278          0.0191   9542
 2 AluY       396      0.0415           345          0.0237   9542
 3 AluJb      348      0.0365           323          0.0222   9542
 4 AluJr      245      0.0257           186          0.0128   9542
 5 AluJo      188      0.0197           185          0.0127   9542
 6 MIRc       148      0.0155           165          0.0113   9542
 7 L1PA4      132      0.0138           154          0.0106   9542
 8 GA-r…      117      0.0123            24          0.00165  9542
 9 L1PA5      117      0.0123           130          0.00892  9542
10 L1PA3      116      0.0122           138          0.00946  9542
# ... with 217 more rows, and 4 more variables: Peak_Shuf <dbl>,
#   Chosen <dbl>, hyper <dbl>, sig <chr>

write.table(repEl_both_sig, "../data/FeatureoverlapPeaks/RepElementsWShuf.txt", col.names = T, row.names = F, quote=F)

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

other attached packages:
 [1] bindrcpp_0.2.2  reshape2_1.4.3  cowplot_0.9.3   forcats_0.3.0  
 [5] stringr_1.4.0   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
[13] workflowr_1.2.0

loaded via a namespace (and not attached):
 [1] tidyselect_0.2.4 haven_1.1.2      lattice_0.20-35  colorspace_1.3-2
 [5] htmltools_0.3.6  yaml_2.2.0       utf8_1.1.4       rlang_0.2.2     
 [9] pillar_1.3.0     glue_1.3.0       withr_2.1.2      modelr_0.1.2    
[13] readxl_1.1.0     bindr_0.1.1      plyr_1.8.4       munsell_0.5.0   
[17] gtable_0.2.0     cellranger_1.1.0 rvest_0.3.2      evaluate_0.13   
[21] labeling_0.3     knitr_1.20       fansi_0.4.0      broom_0.5.0     
[25] Rcpp_0.12.19     scales_1.0.0     backports_1.1.2  jsonlite_1.6    
[29] fs_1.2.6         hms_0.4.2        digest_0.6.17    stringi_1.2.4   
[33] grid_3.5.1       rprojroot_1.3-2  cli_1.0.1        tools_3.5.1     
[37] magrittr_1.5     lazyeval_0.2.1   crayon_1.3.4     whisker_0.3-2   
[41] pkgconfig_2.0.2  xml2_1.2.0       lubridate_1.7.4  assertthat_0.2.0
[45] rmarkdown_1.11   httr_1.3.1       rstudioapi_0.9.0 R6_2.3.0        
[49] nlme_3.1-137     git2r_0.24.0     compiler_3.5.1