SignalSiteLoc.Rmd

---
title: "Location of Signal Site"
author: "Briana Mittleman"
date: "3/6/2019"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

In the [Signal Site enrichment analysis](SignalSiteEnrich.html) I looked at the peaks to see if signal sites are enriched upstream of my peaks. I found this is true but now I want to see where the signal sites are in comparison to my peaks. I am going to use the biostrings package tool matchPWM for this analysis.
```{r}
library(workflowr)
library(tidyverse)
library(Biostrings)
library(BSgenome)
library(genomation)

```
##Make location file
I need to get the coordinates for the regions I care about. I want to look at the peak and 150bp upstream. This is probably larger than I will need  to look at but it will be good to have an inclusive look first.  


I want to use the peak file and make a file that is the peak and upstream 150:  


Upstream150Bases.py
```{python, eval=F}
#python  
def main(Fin, Fout):
  outBed=open(Fout, "w")
  chrom_lengths=open("/project2/gilad/briana/genome_anotation_data/chrom_lengths2.sort.bed","r")
  #make a dictionary with chrom lengths
  length_dic={}
  for i in chrom_lengths:
    chrom, start, end = i.split()
    length_dic[str(chrom)]=int(end)  

#write file 
  for ln in open(Fin):
    chrom, start, end, name, score, strand = ln.split()
    chrom=str(chrom)
    if strand=="+":
      start_new=int(start)-150
      if start_new <= 1:
        start_new = 0 
      end_new= int(end)
      if end_new == 0:
        end_new=1
      outBed.write("%s\t%d\t%d\t%s\t%s\t%s\n"%(chrom, start_new, end_new, name, score, strand))
    if strand == "-":
      start_new=int(start)
      end_new=int(end) + 150
      outBed.write("%s\t%d\t%d\t%s\t%s\t%s\n"%(chrom, start_new, end_new, name, score, strand))
  outBed.close()  

if __name__ == "__main__":
    import sys
    inFile = sys.argv[1]
    outFile=sys.argv[2] 
    main(inFile, outFile)
```



run_get150up.sh
```{bash,eval=F}
#!/bin/bash

#SBATCH --job-name=run_get150up
#SBATCH --account=pi-yangili1
#SBATCH --time=36:00:00
#SBATCH --output=run_get150upt.out
#SBATCH --error=run_get150up.err
#SBATCH --partition=broadwl
#SBATCH --mem=16G
#SBATCH --mail-type=END

module load Anaconda3
source activate three-prime-env

python Upstream150Bases.py /project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPeaks_5percCov_fixedStrand.bed  /project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_peakand150up.bed 
```

##Get subject (reads)

Input the regions: 


Fix chromosomes: 
```{r}
PeakRegions=read.table("../data/Signal_Loc/APAPeaks_5percCov_fixedStrand_peakand150up.bed", header=F,col.names = c("chr","start", "end", "peak", "score", "strand")) %>% mutate(Chrom=paste("chr", chr, sep="")) %>% select(Chrom, start,end,peak,score,strand)
write.table(PeakRegions, file="../data/Signal_Loc/APAPeaks_5percCov_fixedStrand_peakand150up_fixedChr.bed", quote=F, col.names = F, row.names = F, sep="\t")


```



```{r}
#convert to reads 

reads.GR= readGeneric(file="../data/Signal_Loc/APAPeaks_5percCov_fixedStrand_peakand150up_fixedChr.bed",chr =1, start = 2, end =3,  meta.cols =4, header=F, zero.based=TRUE,strand=6)


```

I need to overlap these positions with the genome

###Make motifs  

```{r}
AATAAA= PWM("AATAAA", type = c("log2probratio", "prob"), prior.params = c(A=0.25, C=0.25, G=0.25, T=0.25))
ATTAAA= PWM("ATTAAA", type = c("log2probratio", "prob"), prior.params = c(A=0.25, C=0.25, G=0.25, T=0.25))
AGTAAA= PWM("AGTAAA", type = c("log2probratio", "prob"), prior.params = c(A=0.25, C=0.25, G=0.25, T=0.25))
TATAAA= PWM("TATAAA", type = c("log2probratio", "prob"), prior.params = c(A=0.25, C=0.25, G=0.25, T=0.25))
CATAAA= PWM("CATAAA", type = c("log2probratio", "prob"), prior.params = c(A=0.25, C=0.25, G=0.25, T=0.25))
GATAAA= PWM("GATAAA", type = c("log2probratio", "prob"), prior.params = c(A=0.25, C=0.25, G=0.25, T=0.25))
AATATA= PWM("AATATA", type = c("log2probratio", "prob"), prior.params = c(A=0.25, C=0.25, G=0.25, T=0.25))
AATACA= PWM("AATACA", type = c("log2probratio", "prob"), prior.params = c(A=0.25, C=0.25, G=0.25, T=0.25))
AATAGA= PWM("AATAGA", type = c("log2probratio", "prob"), prior.params = c(A=0.25, C=0.25, G=0.25, T=0.25))
AAAAAG= PWM("AAAAAG", type = c("log2probratio", "prob"), prior.params = c(A=0.25, C=0.25, G=0.25, T=0.25))
ACTAAA= PWM("ACTAAA", type = c("log2probratio", "prob"), prior.params = c(A=0.25, C=0.25, G=0.25, T=0.25))
```



###find the mathes  

```{r}
genome.hg19 <- getBSgenome("BSgenome.Hsapiens.UCSC.hg19")

#matches <- matchPWM(pwm=AATAAA, subject = genome.hg19) %>% keepStandardChromosomes(., species= "Homo sapiens")


DNAstringSetPeaks=data.frame(seq=getSeq(genome.hg19, reads.GR))

x=DNAString(DNAstringSetPeaks[1,1])

hits <- matchPWM(AATAAA,x,with.score=T) 


start(hits)
```



Look over and make hits file for all 

```{r}
list_AATAAA_res=c()
for (i in 1:nrow(DNAstringSetPeaks)){
  x=DNAString(DNAstringSetPeaks[i,1])
  list_AATAAA_res=c(list_AATAAA_res,matchPWM(AATAAA,x,with.score=T))
}
```


Get out the start positions: 


```{r}
starts_AATAAA=c()
nsig=c()
last_oc_AATAAA=c()
for (i in list_AATAAA_res){
  nsig=c(nsig, length(start(i)))
  starts_AATAAA=c(starts_AATAAA, start(i))
  #print(length(start(i)))
  if (length(start(i)) != 0 ){
    last_oc_AATAAA=c(last_oc_AATAAA, max(start(i),na.rm =T))
  }
}
   
```




Histogram of results:

```{r}
summary(starts_AATAAA)
hist(starts_AATAAA,breaks=10000)
```

```{r}
summary(nsig)
sum(nsig==0)
sum(nsig==1)

hist(nsig,breaks=100)
```

Look at the first occurence:  

```{r}
summary(last_oc_AATAAA)
hist(last_oc_AATAAA,breaks=1000)
```

I want to get the closest occurance of the 


##Change analysis: 


I want to start at the end of the peak (most downstream) and look for the signal sites. For confidence a peak is one PAS,  I will look only at peaks less than 100bp long. I will extend the peak upstream 100 basepairs. I will look from the downstream end for the sites.  


###Subset peaks less than 200bp  


filterPeaks100length.py
```{bash,eval=F}
peaks=open("/project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPeaks_5percCov_fixedStrand.bed", "r")

outPeaks=open("/project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_filter4length.bed", "w")

nNotOk=0
for ln in peaks:
    start= int(ln.split()[1])
    end=int(ln.split()[2])
    length=end - start 
    if length <= 100:
        outPeaks.write(ln)
    else:
        nNotOk +=1
        
print(nNotOk)
outPeaks.close()
      
```

This filters 12105 peaks.  

###Get sequence for peaks and 100 upstream


Upstream100Bases_filteredpeaks.py
```{bash,eval=F}
#python  
def main(Fin, Fout):
  outBed=open(Fout, "w")
  chrom_lengths=open("/project2/gilad/briana/genome_anotation_data/chrom_lengths2.sort.bed","r")
  #make a dictionary with chrom lengths
  length_dic={}
  for i in chrom_lengths:
    chrom, start, end = i.split()
    length_dic[str(chrom)]=int(end)  

#write file 
  for ln in open(Fin):
    chrom, start, end, name, score, strand = ln.split()
    chrom=str(chrom)
    if strand=="+":
      start_new=int(start)-100
      if start_new <= 1:
        start_new = 0 
      end_new= int(end)
      if end_new == 0:
        end_new=1
      outBed.write("%s\t%d\t%d\t%s\t%s\t%s\n"%(chrom, start_new, end_new, name, score, strand))
    if strand == "-":
      start_new=int(start)
      end_new=int(end) + 100
      outBed.write("%s\t%d\t%d\t%s\t%s\t%s\n"%(chrom, start_new, end_new, name, score, strand))
  outBed.close()  

if __name__ == "__main__":
    import sys
    inFile = sys.argv[1]
    outFile=sys.argv[2] 
    main(inFile, outFile)
```


Run this:  

```{bash,eval=F}
python Upstream100Bases_filteredpeaks.py /project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_filter4length.bed  /project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_filter4length_upstream100.bed 
```

Run bedtools nuc for this to get the sequences:  



nucpeaksand100up.sh
```{bash,eval=F}
#!/bin/bash

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

module load Anaconda3
source activate three-prime-env

bedtools nuc -seq -fi /project2/gilad/briana/genome_anotation_data/genome/Homo_sapiens.GRCh37.75.dna_sm.all.fa -bed /project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_filter4length_upstream100.bed  > /project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_filter4length_upstream100_Seq.bed  
```


This is looking at the positive strand left to right always. I need to go from the right and look at the reverse signal sites. I can look into ways to flip a string in python 

'a string'[::-1]  




###find distance to peaks
**change region I am looking at before I do this - move to new analysis and try a new method for this analysis**  


DistPAS2Sig.py
```{bash,eval=F}
def main(Insite, out):
  sigsite=[Insite]
  
  inBed=open("/project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_filter4length_upstream100_Seq.bed", "r")
  outRes=open(out, "w")
  
  #function for reverse compliments  
  
  def ReverseComplement1(seq):
      seq_dict = {'A':'T','T':'A','G':'C','C':'G', 'a':'t', 't':'a', 'g':'c', 'c':'g'}
      bases=[seq_dict[base] for base in seq]
      bases=reversed(bases)
      return("".join(bases))
  
  
  #reverse comp each signal site 
  sigsite_revComp=[]
  for i in sigsite:
      sigsite_revComp.append(ReverseComplement1(i))
      
  #want a dictionary for each of the sites and its reverse compliment:
  sigsites_dic={}
  for i in range(len(sigsite)):
      sigsites_dic[sigsite[i]]=sigsite_revComp[i]    
      
  
  #function to get occurance: takes in sig site and sequence (give it the correct stranded stuff)
  
  #make 2 of these, this is for the pos strand
  def getOccurance(sigsite, seq):
      if sigsite in seq:
          length=len(seq)
          pos= seq.rfind(sigsite)
          posF=length-pos
          return(posF)
      else:
          return(-9)
          
  #negative strand occurance function:
  
  def getOccurance_neg(sigsite, seq):
      sigsite=sigsites_dic[sigsite]
      if sigsite in seq:
          pos= seq.find(sigsite)
          return(pos + 6)
      else:
          return(-9)
  
  #i can only addpend the value if the function does not return -9
      
  
  #function i can run on each signal site  
  
  
  #loop through peaks and check for every site, first ask stand and do the rev  
  def loop41site(site):
      resList=[]
      for ln in inBed:
          strand=ln.split()[5]
          seq= ln.split()[15]
          if strand == "+":
              loc= getOccurance(site, seq)
              if loc !=-9:
                  resList.append(loc)
          else: 
              loc=getOccurance_neg(site,seq)
              if loc !=-9:
                  resList.append(loc)
      return(resList)
          
  
  #run this for each sig site
  res_dic={}
  for i in sigsite:
      res_dic[i]=[]
  for i in sigsite:
      reslist=loop41site(i)
      res_dic[i]=reslist
  
  
  outRes.write("%s\n"%(sigsite[0]))
  for i in reslist:
      outRes.write("%d\n"%(i))
      
  
  
  outRes.close()
      
    
if __name__ == "__main__":
    import sys
    Site_in = sys.argv[1]
    outFile= "/project2/gilad/briana/threeprimeseq/data/Signal_Loc/Loc_%s_Distance2end.txt"%(Site_in)
    main(Site_in, outFile)

```


make a test with just 1 site: 


test_DistPAS2Sig.py
```{bash,eval=F}

sigsite=['ATTAAA'] 

inBed=open("/project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_filter4length_upstream100_Seq.bed", "r")
outRes=open('/project2/gilad/briana/threeprimeseq/data/Signal_Loc/test.text', "w")

#function for reverse compliments  

def ReverseComplement1(seq):
    seq_dict = {'A':'T','T':'A','G':'C','C':'G', 'a':'t', 't':'a', 'g':'c', 'c':'g'}
    bases=[seq_dict[base] for base in seq]
    bases=reversed(bases)
    return("".join(bases))


#reverse comp each signal site 
sigsite_revComp=[]
for i in sigsite:
    sigsite_revComp.append(ReverseComplement1(i))
    
#want a dictionary for each of the sites and its reverse compliment:
sigsites_dic={}
for i in range(len(sigsite)):
    sigsites_dic[sigsite[i]]=sigsite_revComp[i]    
    

#function to get occurance: takes in sig site and sequence (give it the correct stranded stuff)

#make 2 of these, this is for the pos strand
def getOccurance(sigsite, seq):
    if sigsite in seq:
        print(sigsite)
        print(seq)
        pos= seq.rfind(sigsite)
        return(pos)
    else:
        return(-9)
        
#negative strand occurance function:

def getOccurance_neg(sigsite, seq):
    sigsite=sigsites_dic[sigsite]
    if sigsite in seq:
        pos= seq.find(sigsite)
        return(pos + 6)
    else:
        return(-9)

#i can only addpend the value if the function does not return -9
    

#function i can run on each signal site  


#loop through peaks and check for every site, first ask stand and do the rev  
def loop41site(site):
    resList=[]
    for ln in inBed:
        strand=ln.split()[5]
        seq= ln.split()[15]
        if strand == "+":
            loc= getOccurance(site, seq)
            print(loc)
            if loc !=-9:
                resList.append(str(loc))
        else: 
            loc=getOccurance_neg(site,seq)
            if loc !=-9:
                resList.append(str(loc))
    return(resList)
        

#run this for each sig site
res_dic={}
for i in sigsite:
    res_dic[i]=[]
for i in sigsite:
    reslist=loop41site(i)
    res_dic[i]=reslist
    
for key, value in res_dic.items():
    valString=":".join(value)
    outRes.write("%s\t%s\n"%(key, valString))


outRes.close()
   
   
   
```



run_DistPAS2Sig.sh

```{bash,eval=F}
#!/bin/bash

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

module load Anaconda3
source activate three-prime-env


python DistPAS2Sig.py AATAAA
python DistPAS2Sig.py ATTAAA
python DistPAS2Sig.py AGTAAA
python DistPAS2Sig.py TATAAA
python DistPAS2Sig.py CATAAA
python DistPAS2Sig.py GATAAA
python DistPAS2Sig.py AATATA
python DistPAS2Sig.py AATACA
python DistPAS2Sig.py AATAGA
python DistPAS2Sig.py AAAAAG
python DistPAS2Sig.py ACTAAA
```

###Make histograms:

```{r}
Loc_AATAAA= read.table("../data/Signal_Loc/Loc_AATAAA_Distance2end.txt", header=T) %>% mutate(Site="AATAAA")
nrow(Loc_AATAAA)

Loc_AAAAAG= read.table("../data/Signal_Loc/Loc_AAAAAG_Distance2end.txt", header=T) %>% mutate(Site="AAAAAG")


Loc_AATACA= read.table("../data/Signal_Loc/Loc_AATACA_Distance2end.txt", header=T) %>% mutate(Site="AATACA")


Loc_AATAGA= read.table("../data/Signal_Loc/Loc_AATAGA_Distance2end.txt", header=T) %>% mutate(Site="AATAGA")


Loc_AATATA= read.table("../data/Signal_Loc/Loc_AATATA_Distance2end.txt", header=T) %>% mutate(Site="AATATA")


Loc_ACTAAA= read.table("../data/Signal_Loc/Loc_ACTAAA_Distance2end.txt", header=T) %>% mutate(Site="ACTAAA")


Loc_AGTAAA= read.table("../data/Signal_Loc/Loc_AGTAAA_Distance2end.txt", header=T) %>% mutate(Site="AGTAAA")



Loc_ATTAAA= read.table("../data/Signal_Loc/Loc_ATTAAA_Distance2end.txt", header=T) %>% mutate(Site="ATTAAA")

Loc_CATAAA= read.table("../data/Signal_Loc/Loc_CATAAA_Distance2end.txt", header=T) %>% mutate(Site="CATAAA")

Loc_GATAAA= read.table("../data/Signal_Loc/Loc_GATAAA_Distance2end.txt", header=T) %>% mutate(Site="GATAAA")

Loc_TATAAA= read.table("../data/Signal_Loc/Loc_TATAAA_Distance2end.txt", header=T) %>% mutate(Site="TATAAA")


```


```{r}
dist_Loc= ggplot(Loc_AATAAA, aes(x=AATAAA)) + 
  geom_histogram(bins=100, fill="red", alpha=.3) + 
  labs(title="Distance from PAS", x="Number of PAS", y="N basepair from PAS") +
  geom_histogram(bins=100,data=Loc_AAAAAG, aes(x=AAAAAG), fill="orange", alpha=.3) +
  geom_histogram(bins=100,data=Loc_AATACA, aes(x=AATACA), fill="yellow", alpha=.3) +
  geom_histogram(bins=100,data=Loc_AATAGA, aes(x=AATAGA), fill="green", alpha=.3) +
  geom_histogram(bins=100,data=Loc_AATATA, aes(x=AATATA), fill="blue", alpha=.3) +
  geom_histogram(bins=100,data=Loc_ACTAAA, aes(x=ACTAAA), fill="purple", alpha=.3) +
  geom_histogram(bins=100,data=Loc_AGTAAA, aes(x=AGTAAA), fill="firebrick3", alpha=.3) +
  geom_histogram(bins=100,data=Loc_AGTAAA, aes(x=AGTAAA), fill="darksalmon", alpha=.3) +
  geom_histogram(bins=100,data=Loc_CATAAA, aes(x=CATAAA), fill="darkslategray", alpha=.3) +
  geom_histogram(bins=100,data=Loc_GATAAA, aes(x=GATAAA), fill="deeppink1", alpha=.3) +
  geom_histogram(bins=100,data=Loc_TATAAA, aes(x=TATAAA), fill="lightcyan1", alpha=.3) 
  
  
dist_Loc

```
 
Make a long dataframe for all of this to make it easier to manipulate.  

```{r}
colnames(Loc_AATAAA)=c("Count", "Site")
colnames(Loc_AAAAAG)=c("Count", "Site")
colnames(Loc_AATACA)=c("Count", "Site")
colnames(Loc_AATAGA)=c("Count", "Site")
colnames(Loc_AATATA)=c("Count", "Site")
colnames(Loc_ACTAAA)=c("Count", "Site")
colnames(Loc_AGTAAA)=c("Count", "Site")
colnames(Loc_CATAAA)=c("Count", "Site")
colnames(Loc_ATTAAA)=c("Count", "Site")
colnames(Loc_GATAAA)=c("Count", "Site")
colnames(Loc_TATAAA)=c("Count", "Site")



AllsiteDF=as.data.frame(rbind(Loc_AATAAA,Loc_AAAAAG,Loc_AATACA,Loc_AATAGA,Loc_AATATA,Loc_ACTAAA,Loc_AGTAAA,Loc_ATTAAA, Loc_GATAAA,Loc_TATAAA,Loc_CATAAA)) %>% mutate(NegCount=-1*as.integer(as.character(Count)), Cononical=ifelse(Site=="AATAAA", "Yes","No"))


 AllsiteDF_to100= AllsiteDF %>% filter(Count < 100)

```


plot: 

```{r}
ggplot(AllsiteDF_to100, aes(group=Site, x=NegCount, fill=Site)) + geom_histogram(position="stack",bins=50 ) + labs(x="Distance from PAS", y="N annotated Sites", title="Location of annotated signal sites")  
```


Do this as proportion: 

```{r}
AllsiteDF_to100_prop=AllsiteDF_to100 %>% group_by(Site,NegCount) %>% summarise(CountperPos=n()) %>% mutate(TotCount=sum(CountperPos),prop=CountperPos/TotCount)

#%>% ungroup() %>% group_by(Site) %>% mutate(nType=sum(Count), prop=CountperPos/nType)
```


```{r}
annoationPAS_allpeak=ggplot(AllsiteDF_to100_prop, aes(fill=Site, y=prop, x=NegCount)) + geom_density(stat="identity") + facet_wrap(~Site) + labs(x="Distance from PAS", y="Proportion of Signal Site", title="Location of annotated signal sites")

annoationPAS_allpeak

ggsave(annoationPAS_allpeak, file="../output/plots/annoationPAS_allpeakPropLocFacet.png")
```

```{r}
annoationPAS_allpeakProphist=ggplot(AllsiteDF_to100_prop, aes(fill=Site, by=Site, y=prop, x=NegCount)) + geom_histogram(stat="identity", position="stack") + labs(x="Distance from PAS", y="Proportion of Signal Site", title="Location of annotated signal sites")
annoationPAS_allpeakProphist
ggsave(annoationPAS_allpeakProphist, file="../output/plots/annoationPAS_allpeakPropLocStackHist.png")
```


###Run for nuclear specific PAS:  


*Used more in NUclear* 
*  /project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPeaks_5percCov_fixedStrand_SigUsageNuc.bed

*Used more in nucelar, in intron* 
* /project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPeaks_5percCov_fixedStrand_SigUsageNuc_Intron.bed


Filter out ones that are too long:  

filterPeaks100length_nuc.py
```{bash,eval=F}
peaks=open(" /project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPeaks_5percCov_fixedStrand_SigUsageNuc.bed", "r")

outPeaks=open("/project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_SigUsageNuc_filter4length.bed", "w")

nNotOk=0
for ln in peaks:
    start= int(ln.split()[1])
    end=int(ln.split()[2])
    length=end - start 
    if length <= 100:
        outPeaks.write(ln)
    else:
        nNotOk +=1
        
print(nNotOk)
outPeaks.close()
      
```

filterPeaks100length_nucintron.py
```{bash,eval=F}
peaks=open("/project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPeaks_5percCov_fixedStrand_SigUsageNuc_Intron.bed", "r")

outPeaks=open("/project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_SigUsageNuc_Intron_filter4length.bed", "w")

nNotOk=0
for ln in peaks:
    start= int(ln.split()[1])
    end=int(ln.split()[2])
    length=end - start 
    if length <= 100:
        outPeaks.write(ln)
    else:
        nNotOk +=1
        
print(nNotOk)
outPeaks.close()
      
```

```{bash,eval=F}
python Upstream100Bases_filteredpeaks.py /project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_SigUsageNuc_Intron_filter4length.bed /project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_SigUsageNuc_Intron_filter4length_upstream100.bed

python Upstream100Bases_filteredpeaks.py /project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_SigUsageNuc_filter4length.bed /project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_SigUsageNuc_filter4length_upstream100.bed

```


nucpeaksand100up_nuc.sh
```{bash,eval=F}
#!/bin/bash

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

module load Anaconda3
source activate three-prime-env

bedtools nuc -seq -fi /project2/gilad/briana/genome_anotation_data/genome/Homo_sapiens.GRCh37.75.dna_sm.all.fa -bed /project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_SigUsageNuc_Intron_filter4length_upstream100.bed  > /project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_SigUsageNuc_Intron_filter4length_upstream100_seq.bed

bedtools nuc -seq -fi /project2/gilad/briana/genome_anotation_data/genome/Homo_sapiens.GRCh37.75.dna_sm.all.fa -bed /project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_SigUsageNuc_filter4length_upstream100.bed  > /project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_SigUsageNuc_filter4length_upstream100_seq.bed



```



DistPAS2Sig_nuclear.py
```{bash,eval=F}
def main(Insite, out):
  sigsite=[Insite]
  
  inBed=open("/project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_SigUsageNuc_filter4length_upstream100_seq.bed", "r")
  outRes=open(out, "w")
  
  #function for reverse compliments  
  
  def ReverseComplement1(seq):
      seq_dict = {'A':'T','T':'A','G':'C','C':'G', 'a':'t', 't':'a', 'g':'c', 'c':'g'}
      bases=[seq_dict[base] for base in seq]
      bases=reversed(bases)
      return("".join(bases))
  
  
  #reverse comp each signal site 
  sigsite_revComp=[]
  for i in sigsite:
      sigsite_revComp.append(ReverseComplement1(i))
      
  #want a dictionary for each of the sites and its reverse compliment:
  sigsites_dic={}
  for i in range(len(sigsite)):
      sigsites_dic[sigsite[i]]=sigsite_revComp[i]    
      
  
  #function to get occurance: takes in sig site and sequence (give it the correct stranded stuff)
  
  #make 2 of these, this is for the pos strand
  def getOccurance(sigsite, seq):
      if sigsite in seq:
          length=len(seq)
          pos= seq.rfind(sigsite)
          posF=length-pos
          return(posF)
      else:
          return(-9)
          
  #negative strand occurance function:
  
  def getOccurance_neg(sigsite, seq):
      sigsite=sigsites_dic[sigsite]
      if sigsite in seq:
          pos= seq.find(sigsite)
          return(pos + 6)
      else:
          return(-9)
  
  #i can only addpend the value if the function does not return -9
      
  
  #function i can run on each signal site  
  
  
  #loop through peaks and check for every site, first ask stand and do the rev  
  def loop41site(site):
      resList=[]
      for ln in inBed:
          strand=ln.split()[5]
          seq= ln.split()[15]
          if strand == "+":
              loc= getOccurance(site, seq)
              if loc !=-9:
                  resList.append(loc)
          else: 
              loc=getOccurance_neg(site,seq)
              if loc !=-9:
                  resList.append(loc)
      return(resList)
          
  
  #run this for each sig site
  res_dic={}
  for i in sigsite:
      res_dic[i]=[]
  for i in sigsite:
      reslist=loop41site(i)
      res_dic[i]=reslist
  
  
  outRes.write("%s\n"%(sigsite[0]))
  for i in reslist:
      outRes.write("%d\n"%(i))
      
  
  
  outRes.close()
      
    
if __name__ == "__main__":
    import sys
    Site_in = sys.argv[1]
    outFile= "/project2/gilad/briana/threeprimeseq/data/Signal_Loc/Loc_%s_Distance2end_nuclear.txt"%(Site_in)
    main(Site_in, outFile)

```


DistPAS2Sig_nucIntron.py
```{bash,eval=F}
def main(Insite, out):
  sigsite=[Insite]
  
  inBed=open("/project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_SigUsageNuc_Intron_filter4length_upstream100_seq.bed", "r")
  outRes=open(out, "w")
  
  #function for reverse compliments  
  
  def ReverseComplement1(seq):
      seq_dict = {'A':'T','T':'A','G':'C','C':'G', 'a':'t', 't':'a', 'g':'c', 'c':'g'}
      bases=[seq_dict[base] for base in seq]
      bases=reversed(bases)
      return("".join(bases))
  
  
  #reverse comp each signal site 
  sigsite_revComp=[]
  for i in sigsite:
      sigsite_revComp.append(ReverseComplement1(i))
      
  #want a dictionary for each of the sites and its reverse compliment:
  sigsites_dic={}
  for i in range(len(sigsite)):
      sigsites_dic[sigsite[i]]=sigsite_revComp[i]    
      
  
  #function to get occurance: takes in sig site and sequence (give it the correct stranded stuff)
  
  #make 2 of these, this is for the pos strand
  def getOccurance(sigsite, seq):
      if sigsite in seq:
          length=len(seq)
          pos= seq.rfind(sigsite)
          posF=length-pos
          return(posF)
      else:
          return(-9)
          
  #negative strand occurance function:
  
  def getOccurance_neg(sigsite, seq):
      sigsite=sigsites_dic[sigsite]
      if sigsite in seq:
          pos= seq.find(sigsite)
          return(pos + 6)
      else:
          return(-9)
  
  #i can only addpend the value if the function does not return -9
      
  
  #function i can run on each signal site  
  
  
  #loop through peaks and check for every site, first ask stand and do the rev  
  def loop41site(site):
      resList=[]
      for ln in inBed:
          strand=ln.split()[5]
          seq= ln.split()[15]
          if strand == "+":
              loc= getOccurance(site, seq)
              if loc !=-9:
                  resList.append(loc)
          else: 
              loc=getOccurance_neg(site,seq)
              if loc !=-9:
                  resList.append(loc)
      return(resList)
          
  
  #run this for each sig site
  res_dic={}
  for i in sigsite:
      res_dic[i]=[]
  for i in sigsite:
      reslist=loop41site(i)
      res_dic[i]=reslist
  
  
  outRes.write("%s\n"%(sigsite[0]))
  for i in reslist:
      outRes.write("%d\n"%(i))
      
  
  
  outRes.close()
      
    
if __name__ == "__main__":
    import sys
    Site_in = sys.argv[1]
    outFile= "/project2/gilad/briana/threeprimeseq/data/Signal_Loc/Loc_%s_Distance2end_nuclearIntron.txt"%(Site_in)
    main(Site_in, outFile)

```

Run both of these:    
run_DistPAS2Sig_nuc.sh
```{bash,eval=F}
#!/bin/bash

#SBATCH --job-name=run_DistPAS2Sig_nuc
#SBATCH --account=pi-yangili1
#SBATCH --time=36:00:00
#SBATCH --output=un_DistPAS2Sig_nuc.out
#SBATCH --error=un_DistPAS2Sig_nuc.err
#SBATCH --partition=broadwl
#SBATCH --mem=16G
#SBATCH --mail-type=END

module load Anaconda3
source activate three-prime-env


python DistPAS2Sig_nucIntron.py AATAAA
python DistPAS2Sig_nucIntron.py ATTAAA
python DistPAS2Sig_nucIntron.py AGTAAA
python DistPAS2Sig_nucIntron.py TATAAA
python DistPAS2Sig_nucIntron.py CATAAA
python DistPAS2Sig_nucIntron.py GATAAA
python DistPAS2Sig_nucIntron.py AATATA
python DistPAS2Sig_nucIntron.py AATACA
python DistPAS2Sig_nucIntron.py AATAGA
python DistPAS2Sig_nucIntron.py AAAAAG
python DistPAS2Sig_nucIntron.py ACTAAA


python DistPAS2Sig_nuclear.py AATAAA
python DistPAS2Sig_nuclear.py ATTAAA
python DistPAS2Sig_nuclear.py AGTAAA
python DistPAS2Sig_nuclear.py TATAAA
python DistPAS2Sig_nuclear.py CATAAA
python DistPAS2Sig_nuclear.py GATAAA
python DistPAS2Sig_nuclear.py AATATA
python DistPAS2Sig_nuclear.py AATACA
python DistPAS2Sig_nuclear.py AATAGA
python DistPAS2Sig_nuclear.py AAAAAG
python DistPAS2Sig_nuclear.py ACTAAA

```


```{r}
Loc_AATAAA_Nuc= read.table("../data/Signal_Loc/Loc_AATAAA_Distance2end_nuclear.txt", header=T) %>% mutate(Site="AATAAA")


Loc_AAAAAG_Nuc= read.table("../data/Signal_Loc/Loc_AAAAAG_Distance2end_nuclear.txt", header=T) %>% mutate(Site="AAAAAG")


Loc_AATACA_Nuc= read.table("../data/Signal_Loc/Loc_AATACA_Distance2end_nuclear.txt", header=T) %>% mutate(Site="AATACA")


Loc_AATAGA_Nuc= read.table("../data/Signal_Loc/Loc_AATAGA_Distance2end_nuclear.txt", header=T) %>% mutate(Site="AATAGA")


Loc_AATATA_Nuc =read.table("../data/Signal_Loc/Loc_AATATA_Distance2end_nuclear.txt", header=T) %>% mutate(Site="AATATA")


Loc_ACTAAA_Nuc= read.table("../data/Signal_Loc/Loc_ACTAAA_Distance2end_nuclear.txt", header=T) %>% mutate(Site="ACTAAA")


Loc_AGTAAA_Nuc= read.table("../data/Signal_Loc/Loc_AGTAAA_Distance2end_nuclear.txt", header=T) %>% mutate(Site="AGTAAA")



Loc_ATTAAA_Nuc=read.table("../data/Signal_Loc/Loc_ATTAAA_Distance2end_nuclear.txt", header=T) %>% mutate(Site="ATTAAA")

Loc_CATAAA_Nuc= read.table("../data/Signal_Loc/Loc_CATAAA_Distance2end_nuclear.txt", header=T) %>% mutate(Site="CATAAA")

Loc_GATAAA_Nuc= read.table("../data/Signal_Loc/Loc_GATAAA_Distance2end_nuclear.txt", header=T) %>% mutate(Site="GATAAA")

Loc_TATAAA_Nuc= read.table("../data/Signal_Loc/Loc_TATAAA_Distance2end_nuclear.txt", header=T) %>% mutate(Site="TATAAA")



dist_Loc_Nuc= ggplot(Loc_AATAAA_Nuc, aes(x=AATAAA)) + 
  geom_histogram(bins=100, fill="red", alpha=.3) + 
  labs(title="Distance from PAS in Nuclear Used peaks", x="Number of PAS", y="N basepair from PAS") +
  geom_histogram(bins=100,data=Loc_AAAAAG_Nuc, aes(x=AAAAAG), fill="orange", alpha=.3) +
  geom_histogram(bins=100,data=Loc_AATACA_Nuc, aes(x=AATACA), fill="yellow", alpha=.3) +
  geom_histogram(bins=100,data=Loc_AATAGA_Nuc, aes(x=AATAGA), fill="green", alpha=.3) +
  geom_histogram(bins=100,data=Loc_AATATA_Nuc, aes(x=AATATA), fill="blue", alpha=.3) +
  geom_histogram(bins=100,data=Loc_ACTAAA_Nuc, aes(x=ACTAAA), fill="purple", alpha=.3) +
  geom_histogram(bins=100,data=Loc_AGTAAA_Nuc, aes(x=AGTAAA), fill="firebrick3", alpha=.3) +
  geom_histogram(bins=100,data=Loc_AGTAAA_Nuc, aes(x=AGTAAA), fill="darksalmon", alpha=.3) +
  geom_histogram(bins=100,data=Loc_CATAAA_Nuc, aes(x=CATAAA), fill="darkslategray", alpha=.3) +
  geom_histogram(bins=100,data=Loc_GATAAA_Nuc, aes(x=GATAAA), fill="deeppink1", alpha=.3) +
  geom_histogram(bins=100,data=Loc_TATAAA_Nuc, aes(x=TATAAA), fill="lightcyan1", alpha=.3) 
  
  
dist_Loc_Nuc



```

Make a long dataframe for all of this to make it easier to manipulate.  

```{r}
colnames(Loc_AATAAA_Nuc)=c("Count", "Site")
colnames(Loc_AAAAAG_Nuc)=c("Count", "Site")
colnames(Loc_AATACA_Nuc)=c("Count", "Site")
colnames(Loc_AATAGA_Nuc)=c("Count", "Site")
colnames(Loc_AATATA_Nuc)=c("Count", "Site")
colnames(Loc_ACTAAA_Nuc)=c("Count", "Site")
colnames(Loc_AGTAAA_Nuc)=c("Count", "Site")
colnames(Loc_CATAAA_Nuc)=c("Count", "Site")
colnames(Loc_ATTAAA_Nuc)=c("Count", "Site")
colnames(Loc_GATAAA_Nuc)=c("Count", "Site")
colnames(Loc_TATAAA_Nuc)=c("Count", "Site")



AllsiteDF_Nuc=as.data.frame(rbind(Loc_AATAAA_Nuc,Loc_AAAAAG_Nuc,Loc_AATACA_Nuc,Loc_AATAGA_Nuc,Loc_AATATA_Nuc,Loc_ACTAAA_Nuc,Loc_AGTAAA_Nuc,Loc_ATTAAA_Nuc, Loc_GATAAA_Nuc,Loc_TATAAA_Nuc,Loc_CATAAA_Nuc)) %>% mutate(NegCount=-1*as.integer(as.character(Count)), Cononical=ifelse(Site=="AATAAA", "Yes","No"))


 AllsiteDF_to100_Nuc= AllsiteDF_Nuc %>% filter(Count < 100)

```


plot: 

```{r}
ggplot(AllsiteDF_to100_Nuc, aes(group=Site, x=NegCount, fill=Site)) + geom_histogram(position="stack",bins=50 ) + labs(x="Distance from PAS", y="N annotated Sites", title="Location of annotated signal sites- Used more in Nuclear")  
```


Do this as proportion: 

```{r}
AllsiteDF_to100_prop_Nuc=AllsiteDF_to100_Nuc %>% group_by(Site,NegCount) %>% summarise(CountperPos=n()) %>% mutate(TotCount=sum(CountperPos),prop=CountperPos/TotCount)

#%>% ungroup() %>% group_by(Site) %>% mutate(nType=sum(Count), prop=CountperPos/nType)
```


```{r}
annoationPAS_Nucpeak=ggplot(AllsiteDF_to100_prop_Nuc, aes(fill=Site, y=prop, x=NegCount)) + geom_histogram(stat="identity") + facet_wrap(~Site) + labs(x="Distance from PAS", y="Proportion of Signal Site", title="Location of annotated signal sites for PAS used more in Nuclear")

annoationPAS_Nucpeak

ggsave(annoationPAS_Nucpeak, file="../output/plots/annoationPAS_NuclearpeaksPropLocFacet.png")
```

```{r}
annoationPAS_NucpeakProphist=ggplot(AllsiteDF_to100_prop_Nuc, aes(fill=Site, by=Site, y=prop, x=NegCount)) + geom_histogram(stat="identity", position="stack") + labs(x="Distance from PAS", y="Proportion of Signal Site", title="Location of annotated signal sites used more in Nuc")
annoationPAS_NucpeakProphist
ggsave(annoationPAS_NucpeakProphist, file="../output/plots/annoationPAS_NucpeakPropLocStackHist.png")
```


```{r}
Loc_AATAAA_Int= read.table("../data/Signal_Loc/Loc_AATAAA_Distance2end_nuclearIntron.txt", header=T)%>% mutate(Site="AATAAA")


Loc_AAAAAG_Int= read.table("../data/Signal_Loc/Loc_AAAAAG_Distance2end_nuclearIntron.txt", header=T)%>% mutate(Site="AAAAAG")


Loc_AATACA_Int= read.table("../data/Signal_Loc/Loc_AATACA_Distance2end_nuclearIntron.txt", header=T)%>% mutate(Site="AATACA")


Loc_AATAGA_Int= read.table("../data/Signal_Loc/Loc_AATAGA_Distance2end_nuclearIntron.txt", header=T)%>% mutate(Site="ATAGA")


Loc_AATATA_Int =read.table("../data/Signal_Loc/Loc_AATATA_Distance2end_nuclearIntron.txt", header=T)%>% mutate(Site="AATATA")

Loc_ACTAAA_Int= read.table("../data/Signal_Loc/Loc_ACTAAA_Distance2end_nuclearIntron.txt", header=T)%>% mutate(Site="ACTAAA")


Loc_AGTAAA_Int= read.table("../data/Signal_Loc/Loc_AGTAAA_Distance2end_nuclearIntron.txt", header=T)%>% mutate(Site="AGTAAA")



Loc_ATTAAA_Int=read.table("../data/Signal_Loc/Loc_ATTAAA_Distance2end_nuclearIntron.txt", header=T)%>% mutate(Site="ATTAAA")

Loc_CATAAA_Int= read.table("../data/Signal_Loc/Loc_CATAAA_Distance2end_nuclearIntron.txt", header=T)%>% mutate(Site="CATAAA")

Loc_GATAAA_Int= read.table("../data/Signal_Loc/Loc_GATAAA_Distance2end_nuclearIntron.txt", header=T)%>% mutate(Site="GATAAA")

Loc_TATAAA_Int= read.table("../data/Signal_Loc/Loc_TATAAA_Distance2end_nuclearIntron.txt", header=T)%>% mutate(Site="TATAAA")



dist_Loc_Int= ggplot(Loc_AATAAA_Int, aes(x=AATAAA)) + 
  geom_histogram(bins=100, fill="red", alpha=.3) + 
  labs(title="Distance from PAS in Intronic Nuclear Used peaks", x="Number of PAS", y="N basepair from PAS") +
  geom_histogram(bins=100,data=Loc_AAAAAG_Int, aes(x=AAAAAG), fill="orange", alpha=.3) +
  geom_histogram(bins=100,data=Loc_AATACA_Int, aes(x=AATACA), fill="yellow", alpha=.3) +
  geom_histogram(bins=100,data=Loc_AATAGA_Int, aes(x=AATAGA), fill="green", alpha=.3) +
  geom_histogram(bins=100,data=Loc_AATATA_Int, aes(x=AATATA), fill="blue", alpha=.3) +
  geom_histogram(bins=100,data=Loc_ACTAAA_Int, aes(x=ACTAAA), fill="purple", alpha=.3) +
  geom_histogram(bins=100,data=Loc_AGTAAA_Int, aes(x=AGTAAA), fill="firebrick3", alpha=.3) +
  geom_histogram(bins=100,data=Loc_AGTAAA_Int, aes(x=AGTAAA), fill="darksalmon", alpha=.3) +
  geom_histogram(bins=100,data=Loc_CATAAA_Int, aes(x=CATAAA), fill="darkslategray", alpha=.3) +
  geom_histogram(bins=100,data=Loc_GATAAA_Int, aes(x=GATAAA), fill="deeppink1", alpha=.3) +
  geom_histogram(bins=100,data=Loc_TATAAA_Int, aes(x=TATAAA), fill="lightcyan1", alpha=.3) 
  
  
dist_Loc_Int



```
Make a long dataframe for all of this to make it easier to manipulate.  

```{r}
colnames(Loc_AATAAA_Int)=c("Count", "Site")
colnames(Loc_AAAAAG_Int)=c("Count", "Site")
colnames(Loc_AATACA_Int)=c("Count", "Site")
colnames(Loc_AATAGA_Int)=c("Count", "Site")
colnames(Loc_AATATA_Int)=c("Count", "Site")
colnames(Loc_ACTAAA_Int)=c("Count", "Site")
colnames(Loc_AGTAAA_Int)=c("Count", "Site")
colnames(Loc_CATAAA_Int)=c("Count", "Site")
colnames(Loc_ATTAAA_Int)=c("Count", "Site")
colnames(Loc_GATAAA_Int)=c("Count", "Site")
colnames(Loc_TATAAA_Int)=c("Count", "Site")



AllsiteDF_NucInt=as.data.frame(rbind(Loc_AATAAA_Int,Loc_AAAAAG_Int,Loc_AATACA_Int,Loc_AATAGA_Int,Loc_AATATA_Int,Loc_ACTAAA_Int,Loc_AGTAAA_Int,Loc_ATTAAA_Int, Loc_GATAAA_Int,Loc_TATAAA_Int,Loc_CATAAA_Int)) %>% mutate(NegCount=-1*as.integer(as.character(Count)), Cononical=ifelse(Site=="AATAAA", "Yes","No"))


 AllsiteDF_to100_NucInt= AllsiteDF_NucInt %>% filter(Count < 100)

```


plot: 

```{r}
ggplot(AllsiteDF_to100_NucInt, aes(group=Site, x=NegCount, fill=Site)) + geom_histogram(position="stack",bins=50 ) + labs(x="Distance from PAS", y="N annotated Sites", title="Location of annotated signal sites- Intronic Used more in Nuclear")  
```


Do this as proportion: 

```{r}
AllsiteDF_to100_prop_NucInt=AllsiteDF_to100_NucInt %>% group_by(Site,NegCount) %>% summarise(CountperPos=n()) %>% mutate(TotCount=sum(CountperPos),prop=CountperPos/TotCount)

#%>% ungroup() %>% group_by(Site) %>% mutate(nType=sum(Count), prop=CountperPos/nType)
```


```{r}
annoationPAS_NucIntpeak=ggplot(AllsiteDF_to100_prop_NucInt, aes(fill=Site, y=prop, x=NegCount)) + geom_histogram(stat="identity") + facet_wrap(~Site) + labs(x="Distance from PAS", y="Proportion of Signal Site", title="Location of annotated signal sites for Intronic PAS used more in Nuclear")

annoationPAS_NucIntpeak

ggsave(annoationPAS_NucIntpeak, file="../output/plots/annoationPAS_NuclearIntronicpeaksPropLocFacet.png")
```

```{r}
annoationPAS_NucIntpeakProphist=ggplot(AllsiteDF_to100_prop_NucInt, aes(fill=Site, by=Site, y=prop, x=NegCount)) + geom_histogram(stat="identity", position="stack") + labs(x="Distance from PAS", y="Proportion of Signal Site", title="Location of annotated signal sites Intronic PAS used more in Nuclear")
annoationPAS_NucIntpeakProphist
ggsave(annoationPAS_NucIntpeakProphist, file="../output/plots/annoationPAS_NucIntronicpeakPropLocStackHist.png")
```

###All intronic peaks:  

*/project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPeaks_5percCov_fixedStrand_INTRON.bed  (17854)


filterPeaks100length_Intron.py
```{bash,eval=F}
peaks=open("/project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPeaks_5percCov_fixedStrand_INTRON.bed", "r")

outPeaks=open("/project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_INTRON_filter4length.bed", "w")

nNotOk=0
for ln in peaks:
    start= int(ln.split()[1])
    end=int(ln.split()[2])
    length=end - start 
    if length <= 100:
        outPeaks.write(ln)
    else:
        nNotOk +=1
        
print(nNotOk)
outPeaks.close()
      
```

12250 peaks  


```{bash,eval=F}
python Upstream100Bases_filteredpeaks.py /project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_INTRON_filter4length.bed /project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_INTRON_filter4length_upstream100.bed
```


nucpeaksand100up_intron.sh
```{bash,eval=F}
#!/bin/bash

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

module load Anaconda3
source activate three-prime-env

bedtools nuc -seq -fi /project2/gilad/briana/genome_anotation_data/genome/Homo_sapiens.GRCh37.75.dna_sm.all.fa -bed /project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_INTRON_filter4length_upstream100.bed > /project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_INTRON_filter4length_upstream100_seq.bed


```

DistPAS2Sig_Intron.py
```{bash,eval=F}
def main(Insite, out):
  sigsite=[Insite]
  
  inBed=open("/project2/gilad/briana/threeprimeseq/data/Signal_Loc/APAPeaks_5percCov_fixedStrand_INTRON_filter4length_upstream100_seq.bed", "r")
  outRes=open(out, "w")
  
  #function for reverse compliments  
  
  def ReverseComplement1(seq):
      seq_dict = {'A':'T','T':'A','G':'C','C':'G', 'a':'t', 't':'a', 'g':'c', 'c':'g'}
      bases=[seq_dict[base] for base in seq]
      bases=reversed(bases)
      return("".join(bases))
  
  
  #reverse comp each signal site 
  sigsite_revComp=[]
  for i in sigsite:
      sigsite_revComp.append(ReverseComplement1(i))
      
  #want a dictionary for each of the sites and its reverse compliment:
  sigsites_dic={}
  for i in range(len(sigsite)):
      sigsites_dic[sigsite[i]]=sigsite_revComp[i]    
      
  
  #function to get occurance: takes in sig site and sequence (give it the correct stranded stuff)
  
  #make 2 of these, this is for the pos strand
  def getOccurance(sigsite, seq):
      if sigsite in seq:
          length=len(seq)
          pos= seq.rfind(sigsite)
          posF=length-pos
          return(posF)
      else:
          return(-9)
          
  #negative strand occurance function:
  
  def getOccurance_neg(sigsite, seq):
      sigsite=sigsites_dic[sigsite]
      if sigsite in seq:
          pos= seq.find(sigsite)
          return(pos + 6)
      else:
          return(-9)
  
  #i can only addpend the value if the function does not return -9
      
  
  #function i can run on each signal site  
  
  
  #loop through peaks and check for every site, first ask stand and do the rev  
  def loop41site(site):
      resList=[]
      for ln in inBed:
          strand=ln.split()[5]
          seq= ln.split()[15]
          if strand == "+":
              loc= getOccurance(site, seq)
              if loc !=-9:
                  resList.append(loc)
          else: 
              loc=getOccurance_neg(site,seq)
              if loc !=-9:
                  resList.append(loc)
      return(resList)
          
  
  #run this for each sig site
  res_dic={}
  for i in sigsite:
      res_dic[i]=[]
  for i in sigsite:
      reslist=loop41site(i)
      res_dic[i]=reslist
  
  
  outRes.write("%s\n"%(sigsite[0]))
  for i in reslist:
      outRes.write("%d\n"%(i))
      
  
  
  outRes.close()
      
    
if __name__ == "__main__":
    import sys
    Site_in = sys.argv[1]
    outFile= "/project2/gilad/briana/threeprimeseq/data/Signal_Loc/Loc_%s_Distance2end_Intronic.txt"%(Site_in)
    main(Site_in, outFile)

```


run_DistPAS2Sig_Intron.sh
```{bash,eval=F}
#!/bin/bash

#SBATCH --job-name=run_DistPAS2Sig_Intronc
#SBATCH --account=pi-yangili1
#SBATCH --time=36:00:00
#SBATCH --output=run_DistPAS2Sig_Intron.out
#SBATCH --error=run_DistPAS2Sig_Intron.err
#SBATCH --partition=broadwl
#SBATCH --mem=16G
#SBATCH --mail-type=END

module load Anaconda3
source activate three-prime-env


python DistPAS2Sig_Intron.py AATAAA
python DistPAS2Sig_Intron.py ATTAAA
python DistPAS2Sig_Intron.py AGTAAA
python DistPAS2Sig_Intron.py TATAAA
python DistPAS2Sig_Intron.py CATAAA
python DistPAS2Sig_Intron.py GATAAA
python DistPAS2Sig_Intron.py AATATA
python DistPAS2Sig_Intron.py AATACA
python DistPAS2Sig_Intron.py AATAGA
python DistPAS2Sig_Intron.py AAAAAG
python DistPAS2Sig_Intron.py ACTAAA
```






```{r}
Loc_AATAAA_AllInt= read.table("../data/Signal_Loc/Loc_AATAAA_Distance2end_Intronic.txt", header=T)


Loc_AAAAAG_AllInt= read.table("../data/Signal_Loc/Loc_AAAAAG_Distance2end_Intronic.txt", header=T)


Loc_AATACA_AllInt= read.table("../data/Signal_Loc/Loc_AATACA_Distance2end_Intronic.txt", header=T)


Loc_AATAGA_AllInt= read.table("../data/Signal_Loc/Loc_AATAGA_Distance2end_Intronic.txt", header=T)


Loc_AATATA_AllInt =read.table("../data/Signal_Loc/Loc_AATATA_Distance2end_Intronic.txt", header=T)


Loc_ACTAAA_AllInt= read.table("../data/Signal_Loc/Loc_ACTAAA_Distance2end_Intronic.txt", header=T)


Loc_AGTAAA_AllInt= read.table("../data/Signal_Loc/Loc_AGTAAA_Distance2end_Intronic.txt", header=T)



Loc_ATTAAA_AllInt=read.table("../data/Signal_Loc/Loc_ATTAAA_Distance2end_Intronic.txt", header=T)

Loc_CATAAA_AllInt= read.table("../data/Signal_Loc/Loc_CATAAA_Distance2end_Intronic.txt", header=T)

Loc_GATAAA_AllInt= read.table("../data/Signal_Loc/Loc_GATAAA_Distance2end_Intronic.txt", header=T)

Loc_TATAAA_AllInt= read.table("../data/Signal_Loc/Loc_TATAAA_Distance2end_Intronic.txt", header=T)



dist_Loc_AllInt= ggplot(Loc_AATAAA_AllInt, aes(x=AATAAA)) + 
  geom_histogram(bins=100, fill="red", alpha=.3) + 
  labs(title="Distance from PAS in Intronic peaks", x="Number of PAS", y="N basepair from PAS") +
  geom_histogram(bins=100,data=Loc_AAAAAG_AllInt, aes(x=AAAAAG), fill="orange", alpha=.3) +
  geom_histogram(bins=100,data=Loc_AATACA_AllInt, aes(x=AATACA), fill="yellow", alpha=.3) +
  geom_histogram(bins=100,data=Loc_AATAGA_AllInt, aes(x=AATAGA), fill="green", alpha=.3) +
  geom_histogram(bins=100,data=Loc_AATATA_AllInt, aes(x=AATATA), fill="blue", alpha=.3) +
  geom_histogram(bins=100,data=Loc_ACTAAA_AllInt, aes(x=ACTAAA), fill="purple", alpha=.3) +
  geom_histogram(bins=100,data=Loc_AGTAAA_AllInt, aes(x=AGTAAA), fill="firebrick3", alpha=.3) +
  geom_histogram(bins=100,data=Loc_AGTAAA_AllInt, aes(x=AGTAAA), fill="darksalmon", alpha=.3) +
  geom_histogram(bins=100,data=Loc_CATAAA_AllInt, aes(x=CATAAA), fill="darkslategray", alpha=.3) +
  geom_histogram(bins=100,data=Loc_GATAAA_AllInt, aes(x=GATAAA), fill="deeppink1", alpha=.3) +
  geom_histogram(bins=100,data=Loc_TATAAA_AllInt, aes(x=TATAAA), fill="lightcyan1", alpha=.3)
  
  
dist_Loc_AllInt



```


###Compare cononical vs other:  

```{r}
AllsiteDF_to100_con=AllsiteDF_to100 %>% group_by(NegCount,Cononical) %>% summarise(PerSite=n()) %>% ungroup() %>% group_by(Cononical) %>% mutate(NCon=sum(PerSite), PropPerSite=PerSite/NCon)


ggplot(AllsiteDF_to100_con,aes(x=NegCount, by=Cononical, fill=Cononical, y=PropPerSite)) +geom_histogram(stat="identity", alpha=.5,bins=50) 

#+ facet_grid(~Cononical)
```
Look in nuclear 


```{r}
AllsiteDF_to100_Nuc_con=AllsiteDF_to100_Nuc %>% group_by(NegCount,Cononical) %>% summarise(PerSite=n()) %>% ungroup() %>% group_by(Cononical) %>% mutate(NCon=sum(PerSite), PropPerSite=PerSite/NCon)


ggplot(AllsiteDF_to100_Nuc_con,aes(x=NegCount, by=Cononical, fill=Cononical, y=PropPerSite)) +geom_histogram(stat="identity",alpha=.5,bins=50)
```