Location of Signal Site

Make motifs

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

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)

[1] 215 219 255

Look over and make hits file for all

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:

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:

summary(starts_AATAAA)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    1.0    79.0   175.0   156.7   216.0  1612.0 

hist(starts_AATAAA,breaks=10000)

Past versions of unnamed-chunk-10-1.png

summary(nsig)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   1.000   3.000   3.857   5.000  74.000 

sum(nsig==0)

[1] 3584

sum(nsig==1)

[1] 8358

hist(nsig,breaks=100)

Past versions of unnamed-chunk-11-1.png

Look at the first occurence:

summary(last_oc_AATAAA)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    1.0   209.0   215.0   213.1   231.0  1612.0 

hist(last_oc_AATAAA,breaks=1000)

Past versions of unnamed-chunk-12-1.png

I want to get the closest occurance of the

Subset peaks less than 200bp

filterPeaks100length.py

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

#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:

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

#!/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

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

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

#!/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
python DistPAS2Sig.py AAAAAA

HEREEE

Make histograms:

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

[1] 11809

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")

Loc_AAAAAA= read.table("../data/Signal_Loc/Loc_AAAAAA_Distance2end.txt", header=T) %>% mutate(Site="AAAAAA")

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) 
  geom_histogram(bins=100,data=Loc_AAAAAA, aes(x=AAAAAA), fill="black", alpha=.3)

mapping: x = ~AAAAAA 
geom_bar: na.rm = FALSE
stat_bin: binwidth = NULL, bins = 100, na.rm = FALSE, pad = FALSE
position_stack 

dist_Loc

Past versions of unnamed-chunk-21-1.png

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

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")
colnames(Loc_AAAAAA)=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, Loc_AAAAAA)) %>% mutate(NegCount=-1*as.integer(as.character(Count)), Cononical=ifelse(Site=="AATAAA", "Yes","No"))


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

plot:

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")  

Past versions of unnamed-chunk-23-1.png

Do this as proportion:

allpeak=33432
AllsiteDF_to100_prop=AllsiteDF_to100 %>% group_by(Site,NegCount) %>% summarise(CountperPos=n()) %>% mutate(TotCount=sum(CountperPos),prop=CountperPos/allpeak)

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

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

Past versions of unnamed-chunk-25-1.png

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

Saving 7 x 5 in image

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")

Warning: Ignoring unknown parameters: binwidth, bins, pad

annoationPAS_allpeakProphist

Past versions of unnamed-chunk-26-1.png

ggsave(annoationPAS_allpeakProphist, file="../output/plots/annoationPAS_allpeakPropLocStackHist.png")

Saving 7 x 5 in image

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

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

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

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

#!/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

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

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

#!/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_nucIntron.py AAAAAA


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
python DistPAS2Sig_nuclear.py AAAAAA

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")
Loc_AAAAAA_Nuc= read.table("../data/Signal_Loc/Loc_AAAAAA_Distance2end_nuclear.txt", header=T) %>% mutate(Site="AAAAAA")


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) +
  geom_histogram(bins=100,data=Loc_AAAAAA_Nuc, aes(x=AAAAAA), fill="black", alpha=.3) 
  
  
dist_Loc_Nuc

Past versions of unnamed-chunk-34-1.png

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

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")
colnames(Loc_AAAAAA_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,Loc_AAAAAA_Nuc)) %>% mutate(NegCount=-1*as.integer(as.character(Count)), Cononical=ifelse(Site=="AATAAA", "Yes","No"))


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

plot:

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")  

Past versions of unnamed-chunk-36-1.png

Do this as proportion:

nucmore=657
AllsiteDF_to100_prop_Nuc=AllsiteDF_to100_Nuc %>% group_by(Site,NegCount) %>% summarise(CountperPos=n()) %>% mutate(TotCount=sum(CountperPos),prop=CountperPos/nucmore)

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

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")

Warning: Ignoring unknown parameters: binwidth, bins, pad

annoationPAS_Nucpeak

Past versions of unnamed-chunk-38-1.png

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

Saving 7 x 5 in image

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")

Warning: Ignoring unknown parameters: binwidth, bins, pad

annoationPAS_NucpeakProphist

Past versions of unnamed-chunk-39-1.png

ggsave(annoationPAS_NucpeakProphist, file="../output/plots/annoationPAS_NucpeakPropLocStackHist.png")

Saving 7 x 5 in image

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")
Loc_AAAAAA_Int= read.table("../data/Signal_Loc/Loc_AAAAAA_Distance2end_nuclearIntron.txt", header=T)%>% mutate(Site="AAAAAA")


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) +
  geom_histogram(bins=100,data=Loc_AAAAAA_Int, aes(x=AAAAAA), fill="black", alpha=.3) 
  
  
dist_Loc_Int

Past versions of unnamed-chunk-40-1.png

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

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")
colnames(Loc_AAAAAA_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,Loc_AAAAAA_Int)) %>% mutate(NegCount=-1*as.integer(as.character(Count)), Cononical=ifelse(Site=="AATAAA", "Yes","No"))


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

plot:

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")  

Past versions of unnamed-chunk-42-1.png

Do this as proportion:

nucmore_int=449
AllsiteDF_to100_prop_NucInt=AllsiteDF_to100_NucInt %>% group_by(Site,NegCount) %>% summarise(CountperPos=n()) %>% mutate(TotCount=sum(CountperPos),prop=CountperPos/nucmore_int)

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

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")

Warning: Ignoring unknown parameters: binwidth, bins, pad

annoationPAS_NucIntpeak

Past versions of unnamed-chunk-44-1.png

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

Saving 7 x 5 in image

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")

Warning: Ignoring unknown parameters: binwidth, bins, pad

annoationPAS_NucIntpeakProphist

Past versions of unnamed-chunk-45-1.png

ggsave(annoationPAS_NucIntpeakProphist, file="../output/plots/annoationPAS_NucIntronicpeakPropLocStackHist.png")

Saving 7 x 5 in image

All intronic peaks:

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

filterPeaks100length_Intron.py

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

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

#!/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

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

#!/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
python DistPAS2Sig_Intron.py AAAAAA

Loc_AATAAA_AllInt= read.table("../data/Signal_Loc/Loc_AATAAA_Distance2end_Intronic.txt", header=T) %>% mutate(Site="AATAAA")


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


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


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


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


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


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



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

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

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

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

Loc_AAAAAA_AllInt= read.table("../data/Signal_Loc/Loc_AAAAAA_Distance2end_Intronic.txt", header=T)%>% mutate(Site="AAAAAA")

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) +
  geom_histogram(bins=100,data=Loc_AAAAAA_AllInt, aes(x=AAAAAA), fill="black", alpha=.3) 
  
  
dist_Loc_AllInt

Past versions of unnamed-chunk-51-1.png

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

colnames(Loc_AATAAA_AllInt)=c("Count", "Site")
colnames(Loc_AAAAAG_AllInt)=c("Count", "Site")
colnames(Loc_AATACA_AllInt)=c("Count", "Site")
colnames(Loc_AATAGA_AllInt)=c("Count", "Site")
colnames(Loc_AATATA_AllInt)=c("Count", "Site")
colnames(Loc_ACTAAA_AllInt)=c("Count", "Site")
colnames(Loc_AGTAAA_AllInt)=c("Count", "Site")
colnames(Loc_CATAAA_AllInt)=c("Count", "Site")
colnames(Loc_ATTAAA_AllInt)=c("Count", "Site")
colnames(Loc_GATAAA_AllInt)=c("Count", "Site")
colnames(Loc_TATAAA_AllInt)=c("Count", "Site")
colnames(Loc_AAAAAA_AllInt)=c("Count", "Site")


AllsiteDF_NucAllInt=as.data.frame(rbind(Loc_AATAAA_AllInt,Loc_AAAAAG_AllInt,Loc_AATACA_AllInt,Loc_AATAGA_AllInt,Loc_AATATA_AllInt,Loc_ACTAAA_AllInt,Loc_AGTAAA_AllInt,Loc_ATTAAA_AllInt, Loc_GATAAA_AllInt,Loc_TATAAA_AllInt,Loc_CATAAA_AllInt,Loc_AAAAAA_AllInt)) %>% mutate(NegCount=-1*as.integer(as.character(Count)), Cononical=ifelse(Site=="AATAAA", "Yes","No"))


 AllsiteDF_to100_NucAllInt= AllsiteDF_NucAllInt %>% filter(Count < 100)

plot:

ggplot(AllsiteDF_to100_NucAllInt, 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 PAS")  

Past versions of unnamed-chunk-53-1.png

Do this as proportion:

allIntron=12251 
AllsiteDF_to100_prop_NucAllInt=AllsiteDF_to100_NucAllInt %>% group_by(Site,NegCount) %>% summarise(CountperPos=n()) %>% mutate(TotCount=sum(CountperPos),prop=CountperPos/allIntron)

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

annoationPAS_NucAllIntpeak=ggplot(AllsiteDF_to100_prop_NucAllInt, 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 all Intronic PAS")

Warning: Ignoring unknown parameters: binwidth, bins, pad

annoationPAS_NucAllIntpeak

Past versions of unnamed-chunk-55-1.png

ggsave(annoationPAS_NucAllIntpeak, file="../output/plots/annoationPAS_AllIntronicpeaksPropLocFacet.png")

Saving 7 x 5 in image

annoationPAS_NucAllIntpeakProphist=ggplot(AllsiteDF_to100_prop_NucAllInt, 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 for all Intronic PAS")

Warning: Ignoring unknown parameters: binwidth, bins, pad

annoationPAS_NucAllIntpeakProphist

Past versions of unnamed-chunk-56-1.png

ggsave(annoationPAS_NucAllIntpeakProphist, file="../output/plots/annoationPAS_AllIntronicpeakPropLocStackHist.png")

Saving 7 x 5 in image

Subsample from AllsiteDF_NucAllInt and look at it in comparison to the nuclear specific peaks:

subsample 220

AllsiteDF_to100_NucAllInt_sampl= AllsiteDF_to100_NucAllInt %>% sample_n(138)


AllsiteDF_to100_prop_NucAllInt_samp=AllsiteDF_to100_NucAllInt_sampl %>% group_by(Site,NegCount) %>% summarise(CountperPos=n()) %>% mutate(TotCount=sum(CountperPos),prop=CountperPos/TotCount)


annoationPAS_NucAllIntpeak_samp=ggplot(AllsiteDF_to100_prop_NucAllInt_samp, 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 sample of Intronic PAS")

Warning: Ignoring unknown parameters: binwidth, bins, pad

annoationPAS_NucAllIntpeak_samp

Past versions of unnamed-chunk-57-1.png

ggsave(annoationPAS_NucAllIntpeak_samp, file="../output/plots/annoationPAS_SampleofIntronicpeaksPropLocFacet.png")

Saving 7 x 5 in image

Compare cononical vs other:

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) 

Warning: Ignoring unknown parameters: binwidth, bins, pad

Past versions of unnamed-chunk-58-1.png

#+ facet_grid(~Cononical)

Look in nuclear

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)

Warning: Ignoring unknown parameters: binwidth, bins, pad

Past versions of unnamed-chunk-59-1.png

Make plots with all sites as demoniator

I want to look at percentages but for all sites:

AllsiteDF_to100_fuklprop=AllsiteDF_to100  %>%  group_by(Site,NegCount) %>% summarise(CountperPos=n()) %>% mutate(TotCount=nrow(AllsiteDF_to100),prop=CountperPos/33432)

allSite=ggplot(AllsiteDF_to100_fuklprop, aes(group=Site, x=NegCount, y=prop, fill=Site)) + geom_histogram(stat="identity",position="stack" ) + labs(x="Distance from PAS", y="Proportion of all annotated sites", title="Location of annotated signal sites") 

Warning: Ignoring unknown parameters: binwidth, bins, pad

AllsiteDF_to100_nuc_fullprop=AllsiteDF_to100_Nuc  %>%  group_by(Site,NegCount) %>% summarise(CountperPos=n()) %>% mutate(TotCount=nrow(AllsiteDF_to100_Nuc),prop=CountperPos/657)

nucSite=ggplot(AllsiteDF_to100_nuc_fullprop, aes(group=Site, x=NegCount, y=prop, fill=Site)) + geom_histogram(stat="identity",position="stack" ) + labs(x="Distance from PAS", y="Proportion of all annotated sites", title="Location of annotated signal sites- Nuclear") 

Warning: Ignoring unknown parameters: binwidth, bins, pad

AllsiteDF_to100_nucInt_fullprop=AllsiteDF_to100_NucInt  %>%  group_by(Site,NegCount) %>% summarise(CountperPos=n()) %>% mutate(TotCount=nrow(AllsiteDF_to100_NucInt),prop=CountperPos/449)

nucIntsite=ggplot(AllsiteDF_to100_nucInt_fullprop, aes(group=Site, x=NegCount, y=prop, fill=Site)) + geom_histogram(stat="identity",position="stack" ) + labs(x="Distance from PAS", y="Proportion of all annotated sites", title="Location of annotated signal sites- Nuclear Intronic") 

Warning: Ignoring unknown parameters: binwidth, bins, pad

AllsiteDF_to100_prop_NucAllInt

AllsiteDF_to100_nucAllInt_fullprop=AllsiteDF_to100_NucAllInt  %>%  group_by(Site,NegCount) %>% summarise(CountperPos=n()) %>% mutate(TotCount=nrow(AllsiteDF_to100_NucAllInt),prop=CountperPos/12251)

AllIntsite=ggplot(AllsiteDF_to100_nucAllInt_fullprop, aes(group=Site, x=NegCount, y=prop, fill=Site)) + geom_histogram(stat="identity",position="stack" ) + labs(x="Distance from PAS", y="Proportion of all annotated sites", title="Location of annotated signal sites- All Intronic") 

Warning: Ignoring unknown parameters: binwidth, bins, pad

PropSignalSite_ofall=plot_grid(allSite, AllIntsite, nucSite,nucIntsite)

ggsave(PropSignalSite_ofall, file="../output/plots/PropofAllSigSite_bySet.png", height = 10, width=10)