---
title: "Peptide Annotation and Protein Inference"
author: "Witold Wolski"
date: "March 16, 2017"
output: html_document
vignette: >
  %\VignetteIndexEntry{Peptide Annotation and Protein Inference}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---
<!--
%\VignetteEngine{knitr::rmarkdown}
%\VignetteIndexEntry{Spectra Annotation and Protein Inference}
-->



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


```{r}
library(prozor)
#library(reshape2)

rm(list=ls())

file = system.file("extdata/IDResults.txt.gz" , package = "prozor")
specMeta <- readr::read_tsv(file)

```


```{r}

nrow(specMeta)
hist(specMeta$score, breaks=100)

```

Annotate peptide sequences with protein sequences

```{r}
length(unique(specMeta$peptideSeq))
upeptide <-unique(specMeta$peptideSeq)

resAll <- prozor::readPeptideFasta(system.file("extdata/Annotation_allSeq.fasta.gz" , package = "prozor"))
resRev <- prozor::readPeptideFasta(system.file("extdata/Annotation_canSeq.fasta.gz" , package = "prozor"))

annotAll <- prozor::annotatePeptides(upeptide,resAll)
pcAll <- length(unique(annotAll$peptideSeq))/ length(upeptide)

annotRev <- prozor::annotatePeptides(upeptide, resRev)
pcCan <-length(unique(annotRev$peptideSeq))/ length(upeptide)

barplot(c(Canonical = pcCan, All = pcAll))

```


```{r, fig="Nr of protein ID's per peptide"}

par(mfrow=c(1,3))
plot(sort(table(annotAll$peptideSeq)),axes=F, ylab="Nr protein IDs")
axis(2)
PCProteotypic_all <- sum(table(annotAll$peptideSeq) == 1)/ length(table(annotAll$peptideSeq)) * 100

#plot(sort(table(annotIso$peptideSeq)),axes=F, ylab="Nr protein IDs")
#axis(2)
#PCProteotypic_iso <- sum(table(annotIso$peptideSeq) == 1)/ length(table(annotIso$peptideSeq)) * 100

plot(sort(table(annotRev$peptideSeq)),axes=F, ylab="Nr protein IDs")
axis(2)

PCProteotypic_canonical <- sum(table(annotRev$peptideSeq) == 1)/ length(table(annotRev$peptideSeq)) * 100
```


```{r}
barplot(c(All = PCProteotypic_all, canonical =  PCProteotypic_canonical),las=2, ylab="% proteotypic" )

```



## Do protein inference



```{r}
library(Matrix)
precursors <- unique(subset(specMeta,select = c(peptideModSeq,precursorCharge,peptideSeq )))

```

### For trembl and swissprot 

Annotate the precursors with protein ID's

```{r}

annotatedPrecursors <-merge(precursors ,
                            subset(annotAll, select= c(peptideSeq,proteinID)),
                            by.x="peptideSeq", 
                            by.y="peptideSeq")

annotatedPrecursors$precursorCharge <- annotatedPrecursors$precursorCharge
annotatedPrecursors$peptideModSeq <- annotatedPrecursors$peptideModSeq
head(annotatedPrecursors)
xx<-prepareMatrix(annotatedPrecursors,proteinID = "proteinID", peptideID = "peptideSeq")

library(Matrix)
image(xx)
xxAll <- greedy(xx)

```

### For swissprot no isoforms

```{r}

annotatedPrecursors <-merge(precursors , subset(annotRev, select= c(peptideSeq,proteinID)), by.x="peptideSeq", 
                            by.y="peptideSeq")

annotatedPrecursors$precursorCharge <- annotatedPrecursors$precursorCharge
annotatedPrecursors$peptideModSeq <- annotatedPrecursors$peptideModSeq

xx<-prepareMatrix(annotatedPrecursors ,proteinID = "proteinID", peptideID = "peptideSeq")
image(xx)
xxCAN <- greedy(xx)
```


```{r fig.cap="Number of Proteins after protein inference."}
barplot(c(All = length(unique(unlist(xxAll))) , canonical = length(unique(unlist(xxCAN)))    ))
```


# TODO

[Protein Grouping and Clustering in Scaffold](https://www.dropbox.com/s/a4kqyc4gxln8sfj/scaffold_protein_grouping_clustering.pdf?dl=1)
