Mispriming Approach

Retreive 10 upstream bases for each read

Upstream10Bases.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)-10
      if start_new <= 1:
        start_new = 0 
      end_new= int(start)
      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(end)
      end_new=int(end) + 10
      if end_new >= length_dic[chrom]:
        end_new = length_dic[chrom]
        start_new=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))
  outBed.close()  

if __name__ == "__main__":
    import sys
    inFile = sys.argv[1]
    fileNoPath=inFile.split("/")[-1]
    fileshort=fileNoPath[:-4]
    outFile="/project2/gilad/briana/threeprimeseq/data/bed_10up/" + fileshort + "10up.bed"
    main(inFile, outFile)

I can wrap this for all of the files.
wrap_Upstream10Bases.sh

#!/bin/bash

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


module load Anaconda3  

source activate three-prime-env


for i in $(ls /project2/gilad/briana/threeprimeseq/data/bed_sort/*-combined-sort.bed); do
            python  Upstream10Bases.py  $i 
        done

I need to sort the files:

Next step is running the nuc function to get the sequences of the positions I just put in the bed files.

bedtools nuc

• -fi (fasta file) /project2/gilad/briana/genome_anotation_data/genome/Homo_sapiens.GRCh37.75.dna_sm.all.fa

• -bed results from 10up stream

• -s strand specific

• -seq print exracted sequence

• output

Nuc10BasesUp.sh

#!/bin/bash

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

for i in $(ls /project2/gilad/briana/threeprimeseq/data/bed_10up/*);do
   describer=$(echo ${i} | sed -e 's/.*YL-SP-//' | sed -e "s/-combined-sort10up.bed$//")
   bedtools nuc -s -seq -fi /project2/gilad/briana/genome_anotation_data/genome/Homo_sapiens.GRCh37.75.dna_sm.all.fa -bed $i > /project2/gilad/briana/threeprimeseq/data/nuc_10up/TenBaseUP.${describer}.txt   
done

Evaluate problem in 1 line

library(data.table)
require(ggseqlogo)

Loading required package: ggseqlogo

library(workflowr)

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

library(tidyverse)

── Attaching packages ───────────────────────────────────────────────────────────── tidyverse 1.2.1 ──

✔ ggplot2 3.0.0     ✔ purrr   0.2.5
✔ tibble  1.4.2     ✔ dplyr   0.7.6
✔ tidyr   0.8.1     ✔ stringr 1.3.1
✔ readr   1.1.1     ✔ forcats 0.3.0

── Conflicts ──────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::between()   masks data.table::between()
✖ dplyr::filter()    masks stats::filter()
✖ dplyr::first()     masks data.table::first()
✖ dplyr::lag()       masks stats::lag()
✖ dplyr::last()      masks data.table::last()
✖ purrr::transpose() masks data.table::transpose()

res_colNames=c("chrom","start", "end", "name", "score", "strand", "pctAT", "pctGC", "A", "C", "G", "T", "N", "Other", "Length", "Seq")

nuc_18486_N= fread("../data/nuc_10up/TenBaseUP.18486-N.txt", col.names = res_colNames)

#filter for full 10 bp  - removes 422 reads (too close to ends)
nuc_18486_N=nuc_18486_N %>% filter(Length==10)
seqs_18486N= nuc_18486_N$Seq

cs1 = make_col_scheme(chars=c('A', 'T', 'C', 'G', 'N'), groups=c('A', 'T', 'C', 'G', 'N'), cols=c('red', 'blue', 'green', 'yellow', 'pink'))

ggseqlogo(seqs_18486N, col_scheme=cs1,  method = 'prob')

SixT="TTTTTT"
nuc_18486_N_6Ts=nuc_18486_N %>% filter(grepl(SixT, Seq)) 

perc_Bad6T= nrow(nuc_18486_N_6Ts)/nrow(nuc_18486_N)
perc_Bad6T

[1] 0.01797875

nuc_18486_N_70perc= nuc_18486_N%>% mutate(percT=T/Length) %>% filter(percT>=.7) 

perc_Bad70= nrow(nuc_18486_N_70perc)/nrow(nuc_18486_N)

perc_Bad70

[1] 0.460071

nuc_18486_N_bad=  nuc_18486_N%>% mutate(percT=T/Length) %>% filter(percT>=.7 | grepl(SixT, Seq) )

perc_Bad=nrow(nuc_18486_N_bad)/nrow(nuc_18486_N)

perc_Bad

[1] 0.4622981

nuc_18486_N_good=nuc_18486_N%>% mutate(percT=T/Length) %>% filter(percT<.7,  !grepl(SixT, Seq) )


ggseqlogo(nuc_18486_N_good$Seq, col_scheme=cs1,  method = 'prob')

nuc_18486_T= fread("../data/nuc_10up/TenBaseUP.18486-T.txt", col.names = res_colNames)

nuc_18486_T=nuc_18486_T %>% filter(Length==10)
seqs_18486T= nuc_18486_T$Seq

ggseqlogo(seqs_18486T, col_scheme=cs1,  method = 'prob')

nuc_18486_T_6Ts=nuc_18486_T %>% filter(grepl(SixT, Seq)) 

perc_Bad6T_tot= nrow(nuc_18486_T_6Ts)/nrow(nuc_18486_T)
perc_Bad6T_tot

[1] 0.01999222

nuc_18486_T_70perc= nuc_18486_T%>% mutate(percT=T/Length) %>% filter(percT>=.7) 

perc_Bad70_tot= nrow(nuc_18486_T_70perc)/nrow(nuc_18486_T)

perc_Bad70_tot

[1] 0.2460797

nuc_18486_T_bad=  nuc_18486_T%>% mutate(percT=T/Length) %>% filter(percT>=.7, grepl(SixT, Seq) )

perc_Bad_tot=nrow(nuc_18486_T_bad)/nrow(nuc_18486_T)

perc_Bad_tot

[1] 0.01466245

nuc_18486_T_good=nuc_18486_T%>% mutate(percT=T/Length) %>% filter(percT<.7 | !grepl(SixT, Seq) )


ggseqlogo(nuc_18486_T_good$Seq, col_scheme=cs1,  method = 'prob')

For all

I may have to use python when i look at all beacuse this is not fast.

I will look at each read in a file and check if for 70% Ts or 6Ts in a row.

filterMissprimingInNuc10.py

#python  
def main(Fin, Fout):
  outBed=open(Fout, "w")
  inBed=open(Fin, "r")
  for ind, ln in enumerate(inBed):
    if ind >=1:
      chrom,start, end, name, score, strand, pctAT, pctGC, A, C, G, T, N, Other, Length, Sequence = ln.split()
      Tperc= float(T) / float(Length)
      if Tperc < .7:
          if "TTTTTT" not in Sequence:
              start_new=int(start)
              end_new=int(end)
              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]
    fileNoPath=inFile.split("/")[-1]
    sampleName=fileNoPath.split(".")[1]
    outFile="/project2/gilad/briana/threeprimeseq/data/nuc_10up_CleanReads/TenBaseUP." + sampleName + ".CleanReads.bed"
    main(inFile, outFile)

run_filterMissprimingInNuc10.sh

#!/bin/bash

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

for i in $(ls /project2/gilad/briana/threeprimeseq/data/nuc_10up/*);do
    python filterMissprimingInNuc10.py $i  
done 

I will look at these stats then move to getting rid ofthe peaks from these reads.

CleanStats=read.csv("../data/nuc_10up/CleanCount_stats.csv", header = T) %>% separate(Sample_ID, into=c("Sample", "Fraction"), by="_") %>% mutate(Perc_PostFilter=PostMPFilter/mappedReads)


cleanStatPlot=ggplot(CleanStats, aes(x=Sample, by=Fraction, fill=Fraction, y=Perc_PostFilter)) + geom_bar(stat="identity", position = "Dodge") + theme(axis.text.x = element_text(angle = 90, hjust = 1),axis.text.y = element_text(size=12),axis.title.y=element_text(size=10,face="bold"), axis.title.x=element_text(size=12,face="bold"))+ scale_fill_manual(values=c("deepskyblue3","darkviolet")) + labs(y="Percent Reads Passing Misprime Filter", title="Accounting for mispriming in 3' Seq Data")

ggsave(filename = "../output/plots/CleanStatsPlot.png", plot = cleanStatPlot)

Saving 7 x 5 in image

Plot number of clean reads per ind:

ggplot(CleanStats, aes(x=Sample, by=Fraction, fill=Fraction, y=PostMPFilter)) + geom_bar(stat="identity", position = "Dodge") + theme(axis.text.x = element_text(angle = 90, hjust = 1),axis.text.y = element_text(size=12),axis.title.y=element_text(size=10,face="bold"), axis.title.x=element_text(size=12,face="bold"))+scale_fill_manual(values=c("deepskyblue3","darkviolet")) + labs(y="Reads Passing Misprime Filter", title="Accounting for mispriming in 3' Seq Data") + scale_y_log10()

CleanStatsMelt= melt(CleanStats, id.vars=c("Sample", "Fraction")) %>% filter(variable=="PostMPFilter") %>%  group_by(Fraction) %>% summarise(mean=mean(value), sd=sd(value))

ggplot(CleanStatsMelt,aes(x=Fraction, y=mean, fill=Fraction)) + geom_col()+ geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), width=.2)+ theme(axis.text.y = element_text(size=12),axis.title.y=element_text(size=10,face="bold"), axis.title.x=element_text(size=12,face="bold"))+ scale_fill_manual(values=c("deepskyblue3","darkviolet"))+ labs(title="Clean Reads by Fraction", y="Clean Reads")

Extra not using

sort_10upbedFile.sh

#!/bin/bash

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

for i in $( ls /project2/gilad/briana/threeprimeseq/data/bed_10up/*);do
  describer=$(echo ${i} | sed -e 's/.*YL-SP-//' | sed -e "s/-combined-sort10up.bed$//")
  sort -k 1,1 -k2,2n $i > /project2/gilad/briana/threeprimeseq/data/bed_10up_sort/YL-SP-${describer}-combined-sort10up.sort.bed
  done

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS  10.14.1

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] bindrcpp_0.2.2    forcats_0.3.0     stringr_1.3.1    
 [4] dplyr_0.7.6       purrr_0.2.5       readr_1.1.1      
 [7] tidyr_0.8.1       tibble_1.4.2      ggplot2_3.0.0    
[10] tidyverse_1.2.1   workflowr_1.1.1   ggseqlogo_0.1    
[13] data.table_1.11.8

loaded via a namespace (and not attached):
 [1] tidyselect_0.2.4  reshape2_1.4.3    haven_1.1.2      
 [4] lattice_0.20-35   colorspace_1.3-2  htmltools_0.3.6  
 [7] yaml_2.2.0        rlang_0.2.2       R.oo_1.22.0      
[10] pillar_1.3.0      withr_2.1.2       glue_1.3.0       
[13] R.utils_2.7.0     modelr_0.1.2      readxl_1.1.0     
[16] bindr_0.1.1       plyr_1.8.4        munsell_0.5.0    
[19] gtable_0.2.0      cellranger_1.1.0  rvest_0.3.2      
[22] R.methodsS3_1.7.1 evaluate_0.11     labeling_0.3     
[25] knitr_1.20        broom_0.5.0       Rcpp_0.12.19     
[28] scales_1.0.0      backports_1.1.2   jsonlite_1.5     
[31] hms_0.4.2         digest_0.6.17     stringi_1.2.4    
[34] grid_3.5.1        rprojroot_1.3-2   cli_1.0.1        
[37] tools_3.5.1       magrittr_1.5      lazyeval_0.2.1   
[40] crayon_1.3.4      whisker_0.3-2     pkgconfig_2.0.2  
[43] xml2_1.2.0        lubridate_1.7.4   rstudioapi_0.8   
[46] assertthat_0.2.0  rmarkdown_1.10    httr_1.3.1       
[49] R6_2.3.0          nlme_3.1-137      git2r_0.23.0     
[52] compiler_3.5.1