annotate pre_process_protein_name_set.R @ 32:63008bdf576e draft

Uploaded
author bornea
date Tue, 26 Apr 2016 17:16:03 -0400
parents e6e456d3ac14
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
28
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
1 #######################################################################################
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
2 # R-code: Protein Name and Tukey's Normalization
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
3 # Author: Adam L Borne
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
4 # Contributers: Paul A Stewart, Brent Kuenzi
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
5 #######################################################################################
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
6 # Assigns uniprot id from MaxQuant peptides file. Filters and normalizes the
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
7 # intensities of each proteins. Resulting in a one to one list of intensities to
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
8 # uniprot id.
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
9 #######################################################################################
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
10 # Copyright (C) Adam Borne.
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
11 # Permission is granted to copy, distribute and/or modify this document
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
12 # under the terms of the GNU Free Documentation License, Version 1.3
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
13 # or any later version published by the Free Software Foundation;
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
14 # with no Invariant Sections, no Front-Cover Texts, and no Back-Cover Texts.
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
15 # A copy of the license is included in the section entitled "GNU
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
16 # Free Documentation License".
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
17 #######################################################################################
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
18 ## REQUIRED INPUT ##
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
19
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
20 # 1) peptides_file: MaxQuant peptides file.
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
21 #######################################################################################
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
22
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
23
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
24 ins_check_run <- function() {
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
25 if ("affy" %in% rownames(installed.packages())){}
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
26 else {
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
27 source("https://bioconductor.org/biocLite.R")
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
28 biocLite(c('mygene','affy'))
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
29 }
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
30 if ('data.table' %in% rownames(installed.packages())){}
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
31 else {
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
32 install.packages('data.table', repos='http://cran.us.r-project.org')
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
33 }
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
34 if ('stringr' %in% rownames(installed.packages())){}
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
35 else {
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
36 install.packages('stringr', repos='http://cran.us.r-project.org')
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
37 }
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
38 if ('VennDiagram' %in% rownames(installed.packages())){}
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
39 else {
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
40 install.packages('VennDiagram', repos='http://cran.us.r-project.org')
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
41 }
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
42 }
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
43
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
44 ins_check_run()
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
45 library(data.table)
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
46 library(affy)
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
47 library(stringr)
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
48 library(mygene)
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
49 library(VennDiagram)
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
50
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
51
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
52 main <- function(peptides_file, db_path) {
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
53 peptides_file = read.delim(peptides_file,header=TRUE,stringsAsFactors=FALSE,fill=TRUE)
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
54 peptides_txt = peptides_file
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
55 intensity_columns = names(peptides_txt[,str_detect(names(peptides_txt),"Intensity\\.*")])
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
56 # Pulls out all lines with Intensity in them.
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
57 intensity_columns = intensity_columns[2:length(intensity_columns)]
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
58 # Removes the first column that does not have a bait.
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
59 peptides_txt_mapped = as.data.frame(map_peptides_proteins(peptides_txt))
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
60 # This function as below sets every line to a 1 to 1 intensity to each possible protein.
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
61 peptides_txt_mapped$Uniprot = str_extract(peptides_txt_mapped$mapped_protein, "[OPQ][0-9][A-Z0-9]{3}[0-9]|[A-NR-Z][0-9]([A-Z][A-Z0-9]{2}[0-9]){1,2}") #Pulls out just Uniprot id from the script.
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
62 peptides_txt_mapped = subset(peptides_txt_mapped,!is.na(Uniprot))
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
63 # Removes reverse sequences and any that didn't match a uniprot accession.
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
64 columns_comb = c("Uniprot", intensity_columns)
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
65 peptides_mapped_intensity = subset(peptides_txt_mapped, select = columns_comb)
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
66 # Subsets out only the needed cloumns for Tukeys (Uniprot IDS and baited intensities)/
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
67 swissprot_fasta = scan(db_path, what="character")
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
68 peptides_txt_mapped_log2 = peptides_mapped_intensity
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
69 # Takes the log2 of the intensities.
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
70 for (i in intensity_columns) {
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
71 peptides_txt_mapped_log2[,i] = log2(subset(peptides_txt_mapped_log2, select = i))
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
72 }
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
73 # Get the minimum from each column while ignoring the -Inf; get the min of these mins for the
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
74 # global min; breaks when there's only one intensity column.
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
75 global_min = min(apply(peptides_txt_mapped_log2[,2:ncol(peptides_txt_mapped_log2)],2,function(x) {
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
76 min(x[x != -Inf])
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
77 }))
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
78 peptides_txt_mapped_log2[peptides_txt_mapped_log2 == -Inf] <- NA
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
79 #uniprot accessions WITHOUT isoforms; it looks like only contaminants contain isoforms anyways.
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
80 mapped_protein_uniprotonly = str_extract(peptides_txt_mapped_log2$Uniprot,"[OPQ][0-9][A-Z0-9]{3}[0-9]|[A-NR-Z][0-9]([A-Z][A-Z0-9]{2}[0-9]){1,2}")
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
81 mapped_protein_uniprot_accession = str_extract(peptides_txt_mapped_log2$Uniprot,"[OPQ][0-9][A-Z0-9]{3}[0-9](-[0-9]+)?|[A-NR-Z][0-9]([A-Z][A-Z0-9]{2}[0-9]){1,2}(-[0-9]+)?|[OPQ][0-9][A-Z0-9]{3}[0-9]|[A-NR-Z][0-9]([A-Z][A-Z0-9]{2}[0-9]){1,2}")
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
82 peptides_txt_mapped_log2$mapped_protein = mapped_protein_uniprotonly
32
63008bdf576e Uploaded
bornea
parents: 30
diff changeset
83 names_db = str_extract(swissprot_fasta,"[OPQ][0-9][A-Z0-9]{3}[0-9]|[A-NR-Z][0-9]([A-Z][A-Z0-9]{2}[0-9]){1,2}")
63008bdf576e Uploaded
bornea
parents: 30
diff changeset
84 names_db = names_db[!is.na(names_db)]
28
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
85 # Runs the Tukey function returning completed table.
32
63008bdf576e Uploaded
bornea
parents: 30
diff changeset
86 peptides_txt_mapped_log2 = subset(peptides_txt_mapped_log2,mapped_protein %in% names_db)
28
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
87 if (nrow(peptides_txt_mapped_log2) == 0) {
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
88 print("Uniprot Database does not have any of the proteins in the peptides file")
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
89 quit()
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
90 }
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
91 protein_intensities_tukeys = get_protein_values(peptides_txt_mapped_log2,intensity_columns)
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
92 protein_intensities_tukeys[protein_intensities_tukeys == 1] <- 0
30
e6e456d3ac14 Uploaded
bornea
parents: 28
diff changeset
93 # Eliminates all pure NA rows.
e6e456d3ac14 Uploaded
bornea
parents: 28
diff changeset
94 protein_intensities_tukeys = protein_intensities_tukeys[rowSums(is.na(protein_intensities_tukeys)) != (ncol(protein_intensities_tukeys)-1),]
e6e456d3ac14 Uploaded
bornea
parents: 28
diff changeset
95 write.table(protein_intensities_tukeys, "./tukeys_output.txt", row.names = FALSE, col.names = TRUE, quote = FALSE, sep = "\t")
28
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
96
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
97 }
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
98
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
99 map_peptides_proteins = function(peptides_in) {
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
100 peptides_in = subset(peptides_in,peptides_in$Proteins != "")
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
101 results_list = list()
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
102 k = 1
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
103 for (i in 1:nrow(peptides_in)) {
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
104 protein_names = peptides_in[i,"Proteins"]
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
105 protein_names_split = unlist(str_split(protein_names,";"))
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
106 for (j in 1:length(protein_names_split)) {
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
107 peptides_mapped_proteins = data.frame(peptides_in[i,],mapped_protein=protein_names_split[j],stringsAsFactors=FALSE)
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
108 results_list[[k]] = peptides_mapped_proteins
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
109 k = k+1
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
110
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
111 }
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
112 }
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
113 return(rbindlist(results_list))
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
114 }
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
115
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
116 get_protein_values = function(mapped_peptides_in,intensity_columns_list) {
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
117 unique_mapped_proteins_list = unique(mapped_peptides_in$mapped_protein)
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
118 # Gets list of all peptides listed.
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
119 # Generates a blank data frame with clomns of Intensities and rows of Uniprots.
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
120 Tukeys_df = data.frame(mapped_protein = unique_mapped_proteins_list, stringsAsFactors = FALSE )
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
121 for (q in intensity_columns_list) {Tukeys_df[,q] = NA}
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
122 for (i in 1:length(unique_mapped_proteins_list)) {
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
123 mapped_peptides_unique_subset = subset(mapped_peptides_in, mapped_protein == unique_mapped_proteins_list[i])
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
124 # Calculate Tukey's Biweight from library(affy); returns a single numeric.
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
125 # Results_list[[i]] = data.frame(Protein=unique_mapped_proteins_list[i],Peptides_per_protein=nrow(mapped_peptides_unique_subset)).
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
126 for (j in intensity_columns_list) {
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
127 # Populates with new Tukeys values.
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
128 Tukeys_df[i,j] = 2^(tukey.biweight(na.omit(mapped_peptides_unique_subset[,j])))
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
129 }
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
130 }
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
131 return(Tukeys_df)
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
132 }
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
133
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
134
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
135 args <- commandArgs(trailingOnly = TRUE)
dbd1af88f060 Uploaded
bornea
parents: 3
diff changeset
136 main(args[1], args[2])