---
title: "plasmid_pool_barcode_dialOut"
author: "yaara"
date: "6/22/2021"
output: html_document
---

This script confirms barcode diversity based on PCR dial out of a watermelon plasmid pool
Input: Starcode table that was generated by clustering the plasmid dialout fastq file (34 bases comprising 30bp of vatiable barcode and 4 of the backbone) 

Output: plots showing library diversity (barcode diversity for the watermelon plasmid pool should be ~1 million).  

Additional files:
(1) plasmid pool dialout protocol
(2) plasmid pool dialout primers


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

#1. collapse barcodes using stracode
Starcode: https://github.com/gui11aume/starcode
the result of this step should be a 2 column starcode table 

#2. read stracode data
```{r read_data}
starecode_file="/Volumes/BreastSeq2/Me/lineageTracing/1.7.20_wm-lib-23-3_for_addGene/starcode_table" #change this to your local starcode location
Barcode.data=read.table(starecode_file) 
names(Barcode.data)=c("seq","count")
```

#3. keep only sequences that end with GAAG (assuming 34 bp were sequnced this would be 4bp from the plasmid)
```{r match backbone}
Barcode.data=Barcode.data[grepl(".+GAAG$",Barcode.data$seq),]
```


#4. check barocdes follow the semi-random sequence pattern ( should be ~0.9531849 )
```{r barcode_validity}
Barcode.data$semiRandom=ifelse(grepl('^([CG][AT]){15}',Barcode.data$seq,perl=TRUE),1,0)
```


#5. plot library complexity
```{r library_complexity}
library(scales)

#compute cumulative frequencies 
Bardode_feqTable=as.data.frame(table(Barcode.data$count))
Bardode_feqTable$Var1=as.numeric(as.character(Bardode_feqTable$Var1))
Bardode_feqTable$cumulative=cumsum(Bardode_feqTable$Freq/sum(Bardode_feqTable$Freq))
Bardode_feqTable$Ecdf_scaled=Bardode_feqTable$cumulative*max(Bardode_feqTable$Freq)
ymax <- max(Bardode_feqTable$Freq)

#plot
ggplot(data=Bardode_feqTable) +geom_col(aes(Var1, Freq),fill = "red",colour = "red",width = 0.5)+coord_trans(x="log10")+scale_x_continuous(breaks=c(1,2,3,4,5,10,20,40,100,400))+theme_bw() +theme(panel.border = element_blank(), panel.grid.major = element_blank(),panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),axis.text=element_text(size=14,colour ="black"))+geom_line(aes(x=Var1,y=Ecdf_scaled), col="blue")+geom_point(aes(x=Var1,y=Ecdf_scaled), col="blue", lwd=1)+xlab("Counts per barcode")+scale_y_continuous(label=comma,sec.axis = sec_axis(trans=~.*100/max(ymax), name = "Percentage of total unique barcodes"))+ylab("Number of unique barcodes")
```



#6. create consensus logo
```{r logo}
library(Biostrings)
library("seqLogo") #this is a Bioconductor package. See https://www.bioconductor.org/install/ 

#remove sequences that include N's
Barcode.data=Barcode.data[!grepl('.*N.*',Barcode.data$seq),]

#calc PSSM
string_set <- AAStringSet(Barcode.data$seq)
PSSM <- consensusMatrix(string_set, as.prob = TRUE)
p <- makePWM(PSSM)
l=seqLogo(p,xfontsize=10,ic.scale=F)

#only barcode area
PSSM=PSSM[,1:30]
p <- makePWM(PSSM)
l=seqLogo(p,xfontsize=10,ic.scale=F)
```

