annotate baseline/Baseline_Functions.r @ 77:58d2377b507d draft

Uploaded
author davidvanzessen
date Wed, 19 Jun 2019 04:31:44 -0400
parents ba33b94637ca
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
67
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1 #########################################################################################
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2 # License Agreement
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
3 #
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
4 # THIS WORK IS PROVIDED UNDER THE TERMS OF THIS CREATIVE COMMONS PUBLIC LICENSE
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
5 # ("CCPL" OR "LICENSE"). THE WORK IS PROTECTED BY COPYRIGHT AND/OR OTHER
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
6 # APPLICABLE LAW. ANY USE OF THE WORK OTHER THAN AS AUTHORIZED UNDER THIS LICENSE
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
7 # OR COPYRIGHT LAW IS PROHIBITED.
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
8 #
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
9 # BY EXERCISING ANY RIGHTS TO THE WORK PROVIDED HERE, YOU ACCEPT AND AGREE TO BE
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
10 # BOUND BY THE TERMS OF THIS LICENSE. TO THE EXTENT THIS LICENSE MAY BE CONSIDERED
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
11 # TO BE A CONTRACT, THE LICENSOR GRANTS YOU THE RIGHTS CONTAINED HERE IN
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
12 # CONSIDERATION OF YOUR ACCEPTANCE OF SUCH TERMS AND CONDITIONS.
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
13 #
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
14 # BASELIne: Bayesian Estimation of Antigen-Driven Selection in Immunoglobulin Sequences
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
15 # Coded by: Mohamed Uduman & Gur Yaari
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
16 # Copyright 2012 Kleinstein Lab
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
17 # Version: 1.3 (01/23/2014)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
18 #########################################################################################
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
19
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
20 # Global variables
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
21
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
22 FILTER_BY_MUTATIONS = 1000
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
23
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
24 # Nucleotides
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
25 NUCLEOTIDES = c("A","C","G","T")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
26
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
27 # Amino Acids
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
28 AMINO_ACIDS <- c("F", "F", "L", "L", "S", "S", "S", "S", "Y", "Y", "*", "*", "C", "C", "*", "W", "L", "L", "L", "L", "P", "P", "P", "P", "H", "H", "Q", "Q", "R", "R", "R", "R", "I", "I", "I", "M", "T", "T", "T", "T", "N", "N", "K", "K", "S", "S", "R", "R", "V", "V", "V", "V", "A", "A", "A", "A", "D", "D", "E", "E", "G", "G", "G", "G")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
29 names(AMINO_ACIDS) <- c("TTT", "TTC", "TTA", "TTG", "TCT", "TCC", "TCA", "TCG", "TAT", "TAC", "TAA", "TAG", "TGT", "TGC", "TGA", "TGG", "CTT", "CTC", "CTA", "CTG", "CCT", "CCC", "CCA", "CCG", "CAT", "CAC", "CAA", "CAG", "CGT", "CGC", "CGA", "CGG", "ATT", "ATC", "ATA", "ATG", "ACT", "ACC", "ACA", "ACG", "AAT", "AAC", "AAA", "AAG", "AGT", "AGC", "AGA", "AGG", "GTT", "GTC", "GTA", "GTG", "GCT", "GCC", "GCA", "GCG", "GAT", "GAC", "GAA", "GAG", "GGT", "GGC", "GGA", "GGG")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
30 names(AMINO_ACIDS) <- names(AMINO_ACIDS)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
31
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
32 #Amino Acid Traits
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
33 #"*" "A" "C" "D" "E" "F" "G" "H" "I" "K" "L" "M" "N" "P" "Q" "R" "S" "T" "V" "W" "Y"
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
34 #B = "Hydrophobic/Burried" N = "Intermediate/Neutral" S="Hydrophilic/Surface")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
35 TRAITS_AMINO_ACIDS_CHOTHIA98 <- c("*","N","B","S","S","B","N","N","B","S","B","B","S","N","S","S","N","N","B","B","N")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
36 names(TRAITS_AMINO_ACIDS_CHOTHIA98) <- sort(unique(AMINO_ACIDS))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
37 TRAITS_AMINO_ACIDS <- array(NA,21)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
38
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
39 # Codon Table
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
40 CODON_TABLE <- as.data.frame(matrix(NA,ncol=64,nrow=12))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
41
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
42 # Substitution Model: Smith DS et al. 1996
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
43 substitution_Literature_Mouse <- matrix(c(0, 0.156222928, 0.601501588, 0.242275484, 0.172506739, 0, 0.241239892, 0.586253369, 0.54636291, 0.255795364, 0, 0.197841727, 0.290240811, 0.467680608, 0.24207858, 0),nrow=4,byrow=T,dimnames=list(NUCLEOTIDES,NUCLEOTIDES))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
44 substitution_Flu_Human <- matrix(c(0,0.2795596,0.5026927,0.2177477,0.1693210,0,0.3264723,0.5042067,0.4983549,0.3328321,0,0.1688130,0.2021079,0.4696077,0.3282844,0),4,4,byrow=T,dimnames=list(NUCLEOTIDES,NUCLEOTIDES))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
45 substitution_Flu25_Human <- matrix(c(0,0.2580641,0.5163685,0.2255674,0.1541125,0,0.3210224,0.5248651,0.5239281,0.3101292,0,0.1659427,0.1997207,0.4579444,0.3423350,0),4,4,byrow=T,dimnames=list(NUCLEOTIDES,NUCLEOTIDES))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
46 load("FiveS_Substitution.RData")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
47
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
48 # Mutability Models: Shapiro GS et al. 2002
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
49 triMutability_Literature_Human <- matrix(c(0.24, 1.2, 0.96, 0.43, 2.14, 2, 1.11, 1.9, 0.85, 1.83, 2.36, 1.31, 0.82, 0.52, 0.89, 1.33, 1.4, 0.82, 1.83, 0.73, 1.83, 1.62, 1.53, 0.57, 0.92, 0.42, 0.42, 1.47, 3.44, 2.58, 1.18, 0.47, 0.39, 1.12, 1.8, 0.68, 0.47, 2.19, 2.35, 2.19, 1.05, 1.84, 1.26, 0.28, 0.98, 2.37, 0.66, 1.58, 0.67, 0.92, 1.76, 0.83, 0.97, 0.56, 0.75, 0.62, 2.26, 0.62, 0.74, 1.11, 1.16, 0.61, 0.88, 0.67, 0.37, 0.07, 1.08, 0.46, 0.31, 0.94, 0.62, 0.57, 0.29, NA, 1.44, 0.46, 0.69, 0.57, 0.24, 0.37, 1.1, 0.99, 1.39, 0.6, 2.26, 1.24, 1.36, 0.52, 0.33, 0.26, 1.25, 0.37, 0.58, 1.03, 1.2, 0.34, 0.49, 0.33, 2.62, 0.16, 0.4, 0.16, 0.35, 0.75, 1.85, 0.94, 1.61, 0.85, 2.09, 1.39, 0.3, 0.52, 1.33, 0.29, 0.51, 0.26, 0.51, 3.83, 2.01, 0.71, 0.58, 0.62, 1.07, 0.28, 1.2, 0.74, 0.25, 0.59, 1.09, 0.91, 1.36, 0.45, 2.89, 1.27, 3.7, 0.69, 0.28, 0.41, 1.17, 0.56, 0.93, 3.41, 1, 1, NA, 5.9, 0.74, 2.51, 2.24, 2.24, 1.95, 3.32, 2.34, 1.3, 2.3, 1, 0.66, 0.73, 0.93, 0.41, 0.65, 0.89, 0.65, 0.32, NA, 0.43, 0.85, 0.43, 0.31, 0.31, 0.23, 0.29, 0.57, 0.71, 0.48, 0.44, 0.76, 0.51, 1.7, 0.85, 0.74, 2.23, 2.08, 1.16, 0.51, 0.51, 1, 0.5, NA, NA, 0.71, 2.14), nrow=64,byrow=T)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
50 triMutability_Literature_Mouse <- matrix(c(1.31, 1.35, 1.42, 1.18, 2.02, 2.02, 1.02, 1.61, 1.99, 1.42, 2.01, 1.03, 2.02, 0.97, 0.53, 0.71, 1.19, 0.83, 0.96, 0.96, 0, 1.7, 2.22, 0.59, 1.24, 1.07, 0.51, 1.68, 3.36, 3.36, 1.14, 0.29, 0.33, 0.9, 1.11, 0.63, 1.08, 2.07, 2.27, 1.74, 0.22, 1.19, 2.37, 1.15, 1.15, 1.56, 0.81, 0.34, 0.87, 0.79, 2.13, 0.49, 0.85, 0.97, 0.36, 0.82, 0.66, 0.63, 1.15, 0.94, 0.85, 0.25, 0.93, 1.19, 0.4, 0.2, 0.44, 0.44, 0.88, 1.06, 0.77, 0.39, 0, 0, 0, 0, 0, 0, 0.43, 0.43, 0.86, 0.59, 0.59, 0, 1.18, 0.86, 2.9, 1.66, 0.4, 0.2, 1.54, 0.43, 0.69, 1.71, 0.68, 0.55, 0.91, 0.7, 1.71, 0.09, 0.27, 0.63, 0.2, 0.45, 1.01, 1.63, 0.96, 1.48, 2.18, 1.2, 1.31, 0.66, 2.13, 0.49, 0, 0, 0, 2.97, 2.8, 0.79, 0.4, 0.5, 0.4, 0.11, 1.68, 0.42, 0.13, 0.44, 0.93, 0.71, 1.11, 1.19, 2.71, 1.08, 3.43, 0.4, 0.67, 0.47, 1.02, 0.14, 1.56, 1.98, 0.53, 0.33, 0.63, 2.06, 1.77, 1.46, 3.74, 2.93, 2.1, 2.18, 0.78, 0.73, 2.93, 0.63, 0.57, 0.17, 0.85, 0.52, 0.31, 0.31, 0, 0, 0.51, 0.29, 0.83, 0.54, 0.28, 0.47, 0.9, 0.99, 1.24, 2.47, 0.73, 0.23, 1.13, 0.24, 2.12, 0.24, 0.33, 0.83, 1.41, 0.62, 0.28, 0.35, 0.77, 0.17, 0.72, 0.58, 0.45, 0.41), nrow=64,byrow=T)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
51 triMutability_Names <- c("AAA", "AAC", "AAG", "AAT", "ACA", "ACC", "ACG", "ACT", "AGA", "AGC", "AGG", "AGT", "ATA", "ATC", "ATG", "ATT", "CAA", "CAC", "CAG", "CAT", "CCA", "CCC", "CCG", "CCT", "CGA", "CGC", "CGG", "CGT", "CTA", "CTC", "CTG", "CTT", "GAA", "GAC", "GAG", "GAT", "GCA", "GCC", "GCG", "GCT", "GGA", "GGC", "GGG", "GGT", "GTA", "GTC", "GTG", "GTT", "TAA", "TAC", "TAG", "TAT", "TCA", "TCC", "TCG", "TCT", "TGA", "TGC", "TGG", "TGT", "TTA", "TTC", "TTG", "TTT")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
52 load("FiveS_Mutability.RData")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
53
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
54 # Functions
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
55
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
56 # Translate codon to amino acid
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
57 translateCodonToAminoAcid<-function(Codon){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
58 return(AMINO_ACIDS[Codon])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
59 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
60
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
61 # Translate amino acid to trait change
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
62 translateAminoAcidToTraitChange<-function(AminoAcid){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
63 return(TRAITS_AMINO_ACIDS[AminoAcid])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
64 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
65
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
66 # Initialize Amino Acid Trait Changes
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
67 initializeTraitChange <- function(traitChangeModel=1,species=1,traitChangeFileName=NULL){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
68 if(!is.null(traitChangeFileName)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
69 tryCatch(
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
70 traitChange <- read.delim(traitChangeFileName,sep="\t",header=T)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
71 , error = function(ex){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
72 cat("Error|Error reading trait changes. Please check file name/path and format.\n")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
73 q()
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
74 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
75 )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
76 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
77 traitChange <- TRAITS_AMINO_ACIDS_CHOTHIA98
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
78 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
79 TRAITS_AMINO_ACIDS <<- traitChange
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
80 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
81
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
82 # Read in formatted nucleotide substitution matrix
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
83 initializeSubstitutionMatrix <- function(substitutionModel,species,subsMatFileName=NULL){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
84 if(!is.null(subsMatFileName)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
85 tryCatch(
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
86 subsMat <- read.delim(subsMatFileName,sep="\t",header=T)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
87 , error = function(ex){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
88 cat("Error|Error reading substitution matrix. Please check file name/path and format.\n")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
89 q()
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
90 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
91 )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
92 if(sum(apply(subsMat,1,sum)==1)!=4) subsMat = t(apply(subsMat,1,function(x)x/sum(x)))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
93 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
94 if(substitutionModel==1)subsMat <- substitution_Literature_Mouse
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
95 if(substitutionModel==2)subsMat <- substitution_Flu_Human
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
96 if(substitutionModel==3)subsMat <- substitution_Flu25_Human
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
97
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
98 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
99
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
100 if(substitutionModel==0){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
101 subsMat <- matrix(1,4,4)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
102 subsMat[,] = 1/3
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
103 subsMat[1,1] = 0
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
104 subsMat[2,2] = 0
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
105 subsMat[3,3] = 0
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
106 subsMat[4,4] = 0
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
107 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
108
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
109
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
110 NUCLEOTIDESN = c(NUCLEOTIDES,"N", "-")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
111 if(substitutionModel==5){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
112 subsMat <- FiveS_Substitution
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
113 return(subsMat)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
114 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
115 subsMat <- rbind(subsMat,rep(NA,4),rep(NA,4))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
116 return( matrix(data.matrix(subsMat),6,4,dimnames=list(NUCLEOTIDESN,NUCLEOTIDES) ) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
117 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
118 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
119
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
120
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
121 # Read in formatted Mutability file
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
122 initializeMutabilityMatrix <- function(mutabilityModel=1, species=1,mutabilityMatFileName=NULL){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
123 if(!is.null(mutabilityMatFileName)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
124 tryCatch(
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
125 mutabilityMat <- read.delim(mutabilityMatFileName,sep="\t",header=T)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
126 , error = function(ex){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
127 cat("Error|Error reading mutability matrix. Please check file name/path and format.\n")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
128 q()
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
129 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
130 )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
131 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
132 mutabilityMat <- triMutability_Literature_Human
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
133 if(species==2) mutabilityMat <- triMutability_Literature_Mouse
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
134 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
135
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
136 if(mutabilityModel==0){ mutabilityMat <- matrix(1,64,3)}
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
137
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
138 if(mutabilityModel==5){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
139 mutabilityMat <- FiveS_Mutability
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
140 return(mutabilityMat)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
141 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
142 return( matrix( data.matrix(mutabilityMat), 64, 3, dimnames=list(triMutability_Names,1:3)) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
143 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
144 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
145
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
146 # Read FASTA file formats
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
147 # Modified from read.fasta from the seqinR package
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
148 baseline.read.fasta <-
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
149 function (file = system.file("sequences/sample.fasta", package = "seqinr"),
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
150 seqtype = c("DNA", "AA"), as.string = FALSE, forceDNAtolower = TRUE,
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
151 set.attributes = TRUE, legacy.mode = TRUE, seqonly = FALSE,
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
152 strip.desc = FALSE, sizeof.longlong = .Machine$sizeof.longlong,
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
153 endian = .Platform$endian, apply.mask = TRUE)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
154 {
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
155 seqtype <- match.arg(seqtype)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
156
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
157 lines <- readLines(file)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
158
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
159 if (legacy.mode) {
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
160 comments <- grep("^;", lines)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
161 if (length(comments) > 0)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
162 lines <- lines[-comments]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
163 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
164
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
165
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
166 ind_groups<-which(substr(lines, 1L, 3L) == ">>>")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
167 lines_mod<-lines
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
168
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
169 if(!length(ind_groups)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
170 lines_mod<-c(">>>All sequences combined",lines)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
171 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
172
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
173 ind_groups<-which(substr(lines_mod, 1L, 3L) == ">>>")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
174
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
175 lines <- array("BLA",dim=(length(ind_groups)+length(lines_mod)))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
176 id<-sapply(1:length(ind_groups),function(i)ind_groups[i]+i-1)+1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
177 lines[id] <- "THIS IS A FAKE SEQUENCE"
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
178 lines[-id] <- lines_mod
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
179 rm(lines_mod)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
180
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
181 ind <- which(substr(lines, 1L, 1L) == ">")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
182 nseq <- length(ind)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
183 if (nseq == 0) {
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
184 stop("no line starting with a > character found")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
185 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
186 start <- ind + 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
187 end <- ind - 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
188
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
189 while( any(which(ind%in%end)) ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
190 ind=ind[-which(ind%in%end)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
191 nseq <- length(ind)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
192 if (nseq == 0) {
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
193 stop("no line starting with a > character found")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
194 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
195 start <- ind + 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
196 end <- ind - 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
197 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
198
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
199 end <- c(end[-1], length(lines))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
200 sequences <- lapply(seq_len(nseq), function(i) paste(lines[start[i]:end[i]], collapse = ""))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
201 if (seqonly)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
202 return(sequences)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
203 nomseq <- lapply(seq_len(nseq), function(i) {
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
204
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
205 #firstword <- strsplit(lines[ind[i]], " ")[[1]][1]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
206 substr(lines[ind[i]], 2, nchar(lines[ind[i]]))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
207
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
208 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
209 if (seqtype == "DNA") {
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
210 if (forceDNAtolower) {
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
211 sequences <- as.list(tolower(chartr(".","-",sequences)))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
212 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
213 sequences <- as.list(toupper(chartr(".","-",sequences)))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
214 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
215 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
216 if (as.string == FALSE)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
217 sequences <- lapply(sequences, s2c)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
218 if (set.attributes) {
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
219 for (i in seq_len(nseq)) {
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
220 Annot <- lines[ind[i]]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
221 if (strip.desc)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
222 Annot <- substr(Annot, 2L, nchar(Annot))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
223 attributes(sequences[[i]]) <- list(name = nomseq[[i]],
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
224 Annot = Annot, class = switch(seqtype, AA = "SeqFastaAA",
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
225 DNA = "SeqFastadna"))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
226 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
227 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
228 names(sequences) <- nomseq
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
229 return(sequences)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
230 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
231
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
232
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
233 # Replaces non FASTA characters in input files with N
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
234 replaceNonFASTAChars <-function(inSeq="ACGTN-AApA"){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
235 gsub('[^ACGTNacgt[:punct:]-[:punct:].]','N',inSeq,perl=TRUE)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
236 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
237
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
238 # Find the germlines in the FASTA list
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
239 germlinesInFile <- function(seqIDs){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
240 firstChar = sapply(seqIDs,function(x){substr(x,1,1)})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
241 secondChar = sapply(seqIDs,function(x){substr(x,2,2)})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
242 return(firstChar==">" & secondChar!=">")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
243 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
244
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
245 # Find the groups in the FASTA list
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
246 groupsInFile <- function(seqIDs){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
247 sapply(seqIDs,function(x){substr(x,1,2)})==">>"
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
248 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
249
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
250 # In the process of finding germlines/groups, expand from the start to end of the group
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
251 expandTillNext <- function(vecPosToID){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
252 IDs = names(vecPosToID)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
253 posOfInterests = which(vecPosToID)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
254
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
255 expandedID = rep(NA,length(IDs))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
256 expandedIDNames = gsub(">","",IDs[posOfInterests])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
257 startIndexes = c(1,posOfInterests[-1])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
258 stopIndexes = c(posOfInterests[-1]-1,length(IDs))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
259 expandedID = unlist(sapply(1:length(startIndexes),function(i){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
260 rep(i,stopIndexes[i]-startIndexes[i]+1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
261 }))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
262 names(expandedID) = unlist(sapply(1:length(startIndexes),function(i){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
263 rep(expandedIDNames[i],stopIndexes[i]-startIndexes[i]+1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
264 }))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
265 return(expandedID)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
266 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
267
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
268 # Process FASTA (list) to return a matrix[input, germline)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
269 processInputAdvanced <- function(inputFASTA){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
270
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
271 seqIDs = names(inputFASTA)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
272 numbSeqs = length(seqIDs)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
273 posGermlines1 = germlinesInFile(seqIDs)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
274 numbGermlines = sum(posGermlines1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
275 posGroups1 = groupsInFile(seqIDs)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
276 numbGroups = sum(posGroups1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
277 consDef = NA
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
278
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
279 if(numbGermlines==0){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
280 posGermlines = 2
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
281 numbGermlines = 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
282 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
283
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
284 glPositionsSum = cumsum(posGermlines1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
285 glPositions = table(glPositionsSum)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
286 #Find the position of the conservation row
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
287 consDefPos = as.numeric(names(glPositions[names(glPositions)!=0 & glPositions==1]))+1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
288 if( length(consDefPos)> 0 ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
289 consDefID = match(consDefPos, glPositionsSum)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
290 #The coservation rows need to be pulled out and stores seperately
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
291 consDef = inputFASTA[consDefID]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
292 inputFASTA = inputFASTA[-consDefID]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
293
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
294 seqIDs = names(inputFASTA)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
295 numbSeqs = length(seqIDs)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
296 posGermlines1 = germlinesInFile(seqIDs)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
297 numbGermlines = sum(posGermlines1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
298 posGroups1 = groupsInFile(seqIDs)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
299 numbGroups = sum(posGroups1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
300 if(numbGermlines==0){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
301 posGermlines = 2
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
302 numbGermlines = 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
303 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
304 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
305
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
306 posGroups <- expandTillNext(posGroups1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
307 posGermlines <- expandTillNext(posGermlines1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
308 posGermlines[posGroups1] = 0
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
309 names(posGermlines)[posGroups1] = names(posGroups)[posGroups1]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
310 posInput = rep(TRUE,numbSeqs)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
311 posInput[posGroups1 | posGermlines1] = FALSE
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
312
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
313 matInput = matrix(NA, nrow=sum(posInput), ncol=2)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
314 rownames(matInput) = seqIDs[posInput]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
315 colnames(matInput) = c("Input","Germline")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
316
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
317 vecInputFASTA = unlist(inputFASTA)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
318 matInput[,1] = vecInputFASTA[posInput]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
319 matInput[,2] = vecInputFASTA[ which( names(inputFASTA)%in%paste(">",names(posGermlines)[posInput],sep="") )[ posGermlines[posInput]] ]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
320
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
321 germlines = posGermlines[posInput]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
322 groups = posGroups[posInput]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
323
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
324 return( list("matInput"=matInput, "germlines"=germlines, "groups"=groups, "conservationDefinition"=consDef ))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
325 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
326
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
327
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
328 # Replace leading and trailing dashes in the sequence
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
329 replaceLeadingTrailingDashes <- function(x,readEnd){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
330 iiGap = unlist(gregexpr("-",x[1]))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
331 ggGap = unlist(gregexpr("-",x[2]))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
332 #posToChange = intersect(iiGap,ggGap)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
333
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
334
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
335 seqIn = replaceLeadingTrailingDashesHelper(x[1])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
336 seqGL = replaceLeadingTrailingDashesHelper(x[2])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
337 seqTemplate = rep('N',readEnd)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
338 seqIn <- c(seqIn,seqTemplate[(length(seqIn)+1):readEnd])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
339 seqGL <- c(seqGL,seqTemplate[(length(seqGL)+1):readEnd])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
340 # if(posToChange!=-1){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
341 # seqIn[posToChange] = "-"
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
342 # seqGL[posToChange] = "-"
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
343 # }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
344
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
345 seqIn = c2s(seqIn[1:readEnd])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
346 seqGL = c2s(seqGL[1:readEnd])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
347
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
348 lenGL = nchar(seqGL)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
349 if(lenGL<readEnd){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
350 seqGL = paste(seqGL,c2s(rep("N",readEnd-lenGL)),sep="")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
351 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
352
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
353 lenInput = nchar(seqIn)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
354 if(lenInput<readEnd){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
355 seqIn = paste(seqIn,c2s(rep("N",readEnd-lenInput)),sep="")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
356 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
357 return( c(seqIn,seqGL) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
358 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
359
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
360 replaceLeadingTrailingDashesHelper <- function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
361 grepResults = gregexpr("-*",x)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
362 grepResultsPos = unlist(grepResults)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
363 grepResultsLen = attr(grepResults[[1]],"match.length")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
364 #print(paste("x = '", x, "'", sep=""))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
365 x = s2c(x)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
366 if(x[1]=="-"){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
367 x[1:grepResultsLen[1]] = "N"
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
368 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
369 if(x[length(x)]=="-"){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
370 x[(length(x)-grepResultsLen[length(grepResultsLen)]+1):length(x)] = "N"
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
371 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
372 return(x)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
373 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
374
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
375
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
376
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
377
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
378 # Check sequences for indels
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
379 checkForInDels <- function(matInputP){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
380 insPos <- checkInsertion(matInputP)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
381 delPos <- checkDeletions(matInputP)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
382 return(list("Insertions"=insPos, "Deletions"=delPos))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
383 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
384
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
385 # Check sequences for insertions
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
386 checkInsertion <- function(matInputP){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
387 insertionCheck = apply( matInputP,1, function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
388 inputGaps <- as.vector( gregexpr("-",x[1])[[1]] )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
389 glGaps <- as.vector( gregexpr("-",x[2])[[1]] )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
390 return( is.finite( match(FALSE, glGaps%in%inputGaps ) ) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
391 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
392 return(as.vector(insertionCheck))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
393 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
394 # Fix inserstions
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
395 fixInsertions <- function(matInputP){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
396 insPos <- checkInsertion(matInputP)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
397 sapply((1:nrow(matInputP))[insPos],function(rowIndex){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
398 x <- matInputP[rowIndex,]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
399 inputGaps <- gregexpr("-",x[1])[[1]]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
400 glGaps <- gregexpr("-",x[2])[[1]]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
401 posInsertions <- glGaps[!(glGaps%in%inputGaps)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
402 inputInsertionToN <- s2c(x[2])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
403 inputInsertionToN[posInsertions]!="-"
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
404 inputInsertionToN[posInsertions] <- "N"
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
405 inputInsertionToN <- c2s(inputInsertionToN)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
406 matInput[rowIndex,2] <<- inputInsertionToN
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
407 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
408 return(insPos)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
409 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
410
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
411 # Check sequences for deletions
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
412 checkDeletions <-function(matInputP){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
413 deletionCheck = apply( matInputP,1, function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
414 inputGaps <- as.vector( gregexpr("-",x[1])[[1]] )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
415 glGaps <- as.vector( gregexpr("-",x[2])[[1]] )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
416 return( is.finite( match(FALSE, inputGaps%in%glGaps ) ) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
417 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
418 return(as.vector(deletionCheck))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
419 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
420 # Fix sequences with deletions
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
421 fixDeletions <- function(matInputP){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
422 delPos <- checkDeletions(matInputP)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
423 sapply((1:nrow(matInputP))[delPos],function(rowIndex){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
424 x <- matInputP[rowIndex,]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
425 inputGaps <- gregexpr("-",x[1])[[1]]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
426 glGaps <- gregexpr("-",x[2])[[1]]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
427 posDeletions <- inputGaps[!(inputGaps%in%glGaps)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
428 inputDeletionToN <- s2c(x[1])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
429 inputDeletionToN[posDeletions] <- "N"
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
430 inputDeletionToN <- c2s(inputDeletionToN)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
431 matInput[rowIndex,1] <<- inputDeletionToN
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
432 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
433 return(delPos)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
434 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
435
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
436
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
437 # Trim DNA sequence to the last codon
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
438 trimToLastCodon <- function(seqToTrim){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
439 seqLen = nchar(seqToTrim)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
440 trimmedSeq = s2c(seqToTrim)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
441 poi = seqLen
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
442 tailLen = 0
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
443
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
444 while(trimmedSeq[poi]=="-" || trimmedSeq[poi]=="."){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
445 tailLen = tailLen + 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
446 poi = poi - 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
447 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
448
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
449 trimmedSeq = c2s(trimmedSeq[1:(seqLen-tailLen)])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
450 seqLen = nchar(trimmedSeq)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
451 # Trim sequence to last codon
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
452 if( getCodonPos(seqLen)[3] > seqLen )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
453 trimmedSeq = substr(seqToTrim,1, ( (getCodonPos(seqLen)[1])-1 ) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
454
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
455 return(trimmedSeq)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
456 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
457
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
458 # Given a nuclotide position, returns the pos of the 3 nucs that made the codon
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
459 # e.g. nuc 86 is part of nucs 85,86,87
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
460 getCodonPos <- function(nucPos){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
461 codonNum = (ceiling(nucPos/3))*3
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
462 return( (codonNum-2):codonNum)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
463 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
464
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
465 # Given a nuclotide position, returns the codon number
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
466 # e.g. nuc 86 = codon 29
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
467 getCodonNumb <- function(nucPos){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
468 return( ceiling(nucPos/3) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
469 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
470
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
471 # Given a codon, returns all the nuc positions that make the codon
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
472 getCodonNucs <- function(codonNumb){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
473 getCodonPos(codonNumb*3)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
474 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
475
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
476 computeCodonTable <- function(testID=1){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
477
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
478 if(testID<=4){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
479 # Pre-compute every codons
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
480 intCounter = 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
481 for(pOne in NUCLEOTIDES){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
482 for(pTwo in NUCLEOTIDES){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
483 for(pThree in NUCLEOTIDES){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
484 codon = paste(pOne,pTwo,pThree,sep="")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
485 colnames(CODON_TABLE)[intCounter] = codon
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
486 intCounter = intCounter + 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
487 CODON_TABLE[,codon] = mutationTypeOptimized(cbind(permutateAllCodon(codon),rep(codon,12)))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
488 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
489 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
490 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
491 chars = c("N","A","C","G","T", "-")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
492 for(a in chars){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
493 for(b in chars){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
494 for(c in chars){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
495 if(a=="N" | b=="N" | c=="N"){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
496 #cat(paste(a,b,c),sep="","\n")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
497 CODON_TABLE[,paste(a,b,c,sep="")] = rep(NA,12)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
498 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
499 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
500 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
501 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
502
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
503 chars = c("-","A","C","G","T")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
504 for(a in chars){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
505 for(b in chars){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
506 for(c in chars){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
507 if(a=="-" | b=="-" | c=="-"){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
508 #cat(paste(a,b,c),sep="","\n")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
509 CODON_TABLE[,paste(a,b,c,sep="")] = rep(NA,12)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
510 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
511 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
512 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
513 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
514 CODON_TABLE <<- as.matrix(CODON_TABLE)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
515 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
516 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
517
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
518 collapseClone <- function(vecInputSeqs,glSeq,readEnd,nonTerminalOnly=0){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
519 #print(length(vecInputSeqs))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
520 vecInputSeqs = unique(vecInputSeqs)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
521 if(length(vecInputSeqs)==1){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
522 return( list( c(vecInputSeqs,glSeq), F) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
523 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
524 charInputSeqs <- sapply(vecInputSeqs, function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
525 s2c(x)[1:readEnd]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
526 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
527 charGLSeq <- s2c(glSeq)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
528 matClone <- sapply(1:readEnd, function(i){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
529 posNucs = unique(charInputSeqs[i,])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
530 posGL = charGLSeq[i]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
531 error = FALSE
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
532 if(posGL=="-" & sum(!(posNucs%in%c("-","N")))==0 ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
533 return(c("-",error))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
534 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
535 if(length(posNucs)==1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
536 return(c(posNucs[1],error))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
537 else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
538 if("N"%in%posNucs){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
539 error=TRUE
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
540 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
541 if(sum(!posNucs[posNucs!="N"]%in%posGL)==0){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
542 return( c(posGL,error) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
543 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
544 #return( c(sample(posNucs[posNucs!="N"],1),error) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
545 if(nonTerminalOnly==0){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
546 return( c(sample(charInputSeqs[i,charInputSeqs[i,]!="N" & charInputSeqs[i,]!=posGL],1),error) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
547 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
548 posNucs = charInputSeqs[i,charInputSeqs[i,]!="N" & charInputSeqs[i,]!=posGL]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
549 posNucsTable = table(posNucs)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
550 if(sum(posNucsTable>1)==0){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
551 return( c(posGL,error) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
552 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
553 return( c(sample( posNucs[posNucs%in%names(posNucsTable)[posNucsTable>1]],1),error) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
554 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
555 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
556
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
557 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
558 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
559 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
560
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
561
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
562 #print(length(vecInputSeqs))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
563 return(list(c(c2s(matClone[1,]),glSeq),"TRUE"%in%matClone[2,]))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
564 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
565 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
566
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
567 # Compute the expected for each sequence-germline pair
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
568 getExpectedIndividual <- function(matInput){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
569 if( any(grep("multicore",search())) ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
570 facGL <- factor(matInput[,2])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
571 facLevels = levels(facGL)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
572 LisGLs_MutabilityU = mclapply(1:length(facLevels), function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
573 computeMutabilities(facLevels[x])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
574 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
575 facIndex = match(facGL,facLevels)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
576
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
577 LisGLs_Mutability = mclapply(1:nrow(matInput), function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
578 cInput = rep(NA,nchar(matInput[x,1]))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
579 cInput[s2c(matInput[x,1])!="N"] = 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
580 LisGLs_MutabilityU[[facIndex[x]]] * cInput
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
581 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
582
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
583 LisGLs_Targeting = mclapply(1:dim(matInput)[1], function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
584 computeTargeting(matInput[x,2],LisGLs_Mutability[[x]])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
585 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
586
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
587 LisGLs_MutationTypes = mclapply(1:length(matInput[,2]),function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
588 #print(x)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
589 computeMutationTypes(matInput[x,2])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
590 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
591
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
592 LisGLs_Exp = mclapply(1:dim(matInput)[1], function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
593 computeExpected(LisGLs_Targeting[[x]],LisGLs_MutationTypes[[x]])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
594 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
595
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
596 ul_LisGLs_Exp = unlist(LisGLs_Exp)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
597 return(matrix(ul_LisGLs_Exp,ncol=4,nrow=(length(ul_LisGLs_Exp)/4),byrow=T))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
598 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
599 facGL <- factor(matInput[,2])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
600 facLevels = levels(facGL)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
601 LisGLs_MutabilityU = lapply(1:length(facLevels), function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
602 computeMutabilities(facLevels[x])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
603 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
604 facIndex = match(facGL,facLevels)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
605
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
606 LisGLs_Mutability = lapply(1:nrow(matInput), function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
607 cInput = rep(NA,nchar(matInput[x,1]))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
608 cInput[s2c(matInput[x,1])!="N"] = 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
609 LisGLs_MutabilityU[[facIndex[x]]] * cInput
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
610 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
611
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
612 LisGLs_Targeting = lapply(1:dim(matInput)[1], function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
613 computeTargeting(matInput[x,2],LisGLs_Mutability[[x]])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
614 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
615
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
616 LisGLs_MutationTypes = lapply(1:length(matInput[,2]),function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
617 #print(x)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
618 computeMutationTypes(matInput[x,2])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
619 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
620
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
621 LisGLs_Exp = lapply(1:dim(matInput)[1], function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
622 computeExpected(LisGLs_Targeting[[x]],LisGLs_MutationTypes[[x]])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
623 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
624
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
625 ul_LisGLs_Exp = unlist(LisGLs_Exp)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
626 return(matrix(ul_LisGLs_Exp,ncol=4,nrow=(length(ul_LisGLs_Exp)/4),byrow=T))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
627
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
628 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
629 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
630
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
631 # Compute mutabilities of sequence based on the tri-nucleotide model
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
632 computeMutabilities <- function(paramSeq){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
633 seqLen = nchar(paramSeq)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
634 seqMutabilites = rep(NA,seqLen)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
635
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
636 gaplessSeq = gsub("-", "", paramSeq)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
637 gaplessSeqLen = nchar(gaplessSeq)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
638 gaplessSeqMutabilites = rep(NA,gaplessSeqLen)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
639
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
640 if(mutabilityModel!=5){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
641 pos<- 3:(gaplessSeqLen)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
642 subSeq = substr(rep(gaplessSeq,gaplessSeqLen-2),(pos-2),(pos+2))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
643 gaplessSeqMutabilites[pos] =
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
644 tapply( c(
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
645 getMutability( substr(subSeq,1,3), 3) ,
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
646 getMutability( substr(subSeq,2,4), 2),
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
647 getMutability( substr(subSeq,3,5), 1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
648 ),rep(1:(gaplessSeqLen-2),3),mean,na.rm=TRUE
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
649 )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
650 #Pos 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
651 subSeq = substr(gaplessSeq,1,3)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
652 gaplessSeqMutabilites[1] = getMutability(subSeq , 1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
653 #Pos 2
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
654 subSeq = substr(gaplessSeq,1,4)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
655 gaplessSeqMutabilites[2] = mean( c(
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
656 getMutability( substr(subSeq,1,3), 2) ,
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
657 getMutability( substr(subSeq,2,4), 1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
658 ),na.rm=T
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
659 )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
660 seqMutabilites[which(s2c(paramSeq)!="-")]<- gaplessSeqMutabilites
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
661 return(seqMutabilites)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
662 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
663
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
664 pos<- 3:(gaplessSeqLen)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
665 subSeq = substr(rep(gaplessSeq,gaplessSeqLen-2),(pos-2),(pos+2))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
666 gaplessSeqMutabilites[pos] = sapply(subSeq,function(x){ getMutability5(x) }, simplify=T)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
667 seqMutabilites[which(s2c(paramSeq)!="-")]<- gaplessSeqMutabilites
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
668 return(seqMutabilites)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
669 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
670
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
671 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
672
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
673 # Returns the mutability of a triplet at a given position
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
674 getMutability <- function(codon, pos=1:3){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
675 triplets <- rownames(mutability)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
676 mutability[ match(codon,triplets) ,pos]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
677 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
678
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
679 getMutability5 <- function(fivemer){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
680 return(mutability[fivemer])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
681 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
682
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
683 # Returns the substitution probabilty
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
684 getTransistionProb <- function(nuc){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
685 substitution[nuc,]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
686 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
687
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
688 getTransistionProb5 <- function(fivemer){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
689 if(any(which(fivemer==colnames(substitution)))){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
690 return(substitution[,fivemer])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
691 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
692 return(array(NA,4))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
693 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
694 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
695
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
696 # Given a nuc, returns the other 3 nucs it can mutate to
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
697 canMutateTo <- function(nuc){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
698 NUCLEOTIDES[- which(NUCLEOTIDES==nuc)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
699 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
700
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
701 # Given a nucleotide, returns the probabilty of other nucleotide it can mutate to
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
702 canMutateToProb <- function(nuc){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
703 substitution[nuc,canMutateTo(nuc)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
704 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
705
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
706 # Compute targeting, based on precomputed mutatbility & substitution
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
707 computeTargeting <- function(param_strSeq,param_vecMutabilities){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
708
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
709 if(substitutionModel!=5){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
710 vecSeq = s2c(param_strSeq)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
711 matTargeting = sapply( 1:length(vecSeq), function(x) { param_vecMutabilities[x] * getTransistionProb(vecSeq[x]) } )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
712 #matTargeting = apply( rbind(vecSeq,param_vecMutabilities),2, function(x) { as.vector(as.numeric(x[2]) * getTransistionProb(x[1])) } )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
713 dimnames( matTargeting ) = list(NUCLEOTIDES,1:(length(vecSeq)))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
714 return (matTargeting)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
715 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
716
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
717 seqLen = nchar(param_strSeq)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
718 seqsubstitution = matrix(NA,ncol=seqLen,nrow=4)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
719 paramSeq <- param_strSeq
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
720 gaplessSeq = gsub("-", "", paramSeq)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
721 gaplessSeqLen = nchar(gaplessSeq)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
722 gaplessSeqSubstitution = matrix(NA,ncol=gaplessSeqLen,nrow=4)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
723
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
724 pos<- 3:(gaplessSeqLen)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
725 subSeq = substr(rep(gaplessSeq,gaplessSeqLen-2),(pos-2),(pos+2))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
726 gaplessSeqSubstitution[,pos] = sapply(subSeq,function(x){ getTransistionProb5(x) }, simplify=T)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
727 seqsubstitution[,which(s2c(paramSeq)!="-")]<- gaplessSeqSubstitution
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
728 #matTargeting <- param_vecMutabilities %*% seqsubstitution
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
729 matTargeting <- sweep(seqsubstitution,2,param_vecMutabilities,`*`)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
730 dimnames( matTargeting ) = list(NUCLEOTIDES,1:(seqLen))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
731 return (matTargeting)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
732 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
733 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
734
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
735 # Compute the mutations types
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
736 computeMutationTypes <- function(param_strSeq){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
737 #cat(param_strSeq,"\n")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
738 #vecSeq = trimToLastCodon(param_strSeq)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
739 lenSeq = nchar(param_strSeq)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
740 vecCodons = sapply({1:(lenSeq/3)}*3-2,function(x){substr(param_strSeq,x,x+2)})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
741 matMutationTypes = matrix( unlist(CODON_TABLE[,vecCodons]) ,ncol=lenSeq,nrow=4, byrow=F)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
742 dimnames( matMutationTypes ) = list(NUCLEOTIDES,1:(ncol(matMutationTypes)))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
743 return(matMutationTypes)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
744 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
745 computeMutationTypesFast <- function(param_strSeq){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
746 matMutationTypes = matrix( CODON_TABLE[,param_strSeq] ,ncol=3,nrow=4, byrow=F)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
747 #dimnames( matMutationTypes ) = list(NUCLEOTIDES,1:(length(vecSeq)))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
748 return(matMutationTypes)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
749 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
750 mutationTypeOptimized <- function( matOfCodons ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
751 apply( matOfCodons,1,function(x){ mutationType(x[2],x[1]) } )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
752 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
753
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
754 # Returns a vector of codons 1 mutation away from the given codon
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
755 permutateAllCodon <- function(codon){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
756 cCodon = s2c(codon)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
757 matCodons = t(array(cCodon,dim=c(3,12)))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
758 matCodons[1:4,1] = NUCLEOTIDES
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
759 matCodons[5:8,2] = NUCLEOTIDES
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
760 matCodons[9:12,3] = NUCLEOTIDES
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
761 apply(matCodons,1,c2s)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
762 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
763
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
764 # Given two codons, tells you if the mutation is R or S (based on your definition)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
765 mutationType <- function(codonFrom,codonTo){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
766 if(testID==4){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
767 if( is.na(codonFrom) | is.na(codonTo) | is.na(translateCodonToAminoAcid(codonFrom)) | is.na(translateCodonToAminoAcid(codonTo)) ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
768 return(NA)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
769 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
770 mutationType = "S"
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
771 if( translateAminoAcidToTraitChange(translateCodonToAminoAcid(codonFrom)) != translateAminoAcidToTraitChange(translateCodonToAminoAcid(codonTo)) ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
772 mutationType = "R"
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
773 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
774 if(translateCodonToAminoAcid(codonTo)=="*" | translateCodonToAminoAcid(codonFrom)=="*"){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
775 mutationType = "Stop"
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
776 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
777 return(mutationType)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
778 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
779 }else if(testID==5){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
780 if( is.na(codonFrom) | is.na(codonTo) | is.na(translateCodonToAminoAcid(codonFrom)) | is.na(translateCodonToAminoAcid(codonTo)) ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
781 return(NA)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
782 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
783 if(codonFrom==codonTo){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
784 mutationType = "S"
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
785 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
786 codonFrom = s2c(codonFrom)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
787 codonTo = s2c(codonTo)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
788 mutationType = "Stop"
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
789 nucOfI = codonFrom[which(codonTo!=codonFrom)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
790 if(nucOfI=="C"){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
791 mutationType = "R"
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
792 }else if(nucOfI=="G"){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
793 mutationType = "S"
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
794 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
795 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
796 return(mutationType)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
797 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
798 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
799 if( is.na(codonFrom) | is.na(codonTo) | is.na(translateCodonToAminoAcid(codonFrom)) | is.na(translateCodonToAminoAcid(codonTo)) ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
800 return(NA)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
801 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
802 mutationType = "S"
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
803 if( translateCodonToAminoAcid(codonFrom) != translateCodonToAminoAcid(codonTo) ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
804 mutationType = "R"
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
805 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
806 if(translateCodonToAminoAcid(codonTo)=="*" | translateCodonToAminoAcid(codonFrom)=="*"){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
807 mutationType = "Stop"
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
808 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
809 return(mutationType)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
810 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
811 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
812 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
813
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
814
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
815 #given a mat of targeting & it's corresponding mutationtypes returns
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
816 #a vector of Exp_RCDR,Exp_SCDR,Exp_RFWR,Exp_RFWR
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
817 computeExpected <- function(paramTargeting,paramMutationTypes){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
818 # Replacements
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
819 RPos = which(paramMutationTypes=="R")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
820 #FWR
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
821 Exp_R_FWR = sum(paramTargeting[ RPos[which(FWR_Nuc_Mat[RPos]==T)] ],na.rm=T)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
822 #CDR
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
823 Exp_R_CDR = sum(paramTargeting[ RPos[which(CDR_Nuc_Mat[RPos]==T)] ],na.rm=T)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
824 # Silents
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
825 SPos = which(paramMutationTypes=="S")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
826 #FWR
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
827 Exp_S_FWR = sum(paramTargeting[ SPos[which(FWR_Nuc_Mat[SPos]==T)] ],na.rm=T)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
828 #CDR
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
829 Exp_S_CDR = sum(paramTargeting[ SPos[which(CDR_Nuc_Mat[SPos]==T)] ],na.rm=T)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
830
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
831 return(c(Exp_R_CDR,Exp_S_CDR,Exp_R_FWR,Exp_S_FWR))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
832 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
833
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
834 # Count the mutations in a sequence
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
835 # each mutation is treated independently
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
836 analyzeMutations2NucUri_website <- function( rev_in_matrix ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
837 paramGL = rev_in_matrix[2,]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
838 paramSeq = rev_in_matrix[1,]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
839
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
840 #Fill seq with GL seq if gapped
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
841 #if( any(paramSeq=="-") ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
842 # gapPos_Seq = which(paramSeq=="-")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
843 # gapPos_Seq_ToReplace = gapPos_Seq[paramGL[gapPos_Seq] != "-"]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
844 # paramSeq[gapPos_Seq_ToReplace] = paramGL[gapPos_Seq_ToReplace]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
845 #}
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
846
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
847
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
848 #if( any(paramSeq=="N") ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
849 # gapPos_Seq = which(paramSeq=="N")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
850 # gapPos_Seq_ToReplace = gapPos_Seq[paramGL[gapPos_Seq] != "N"]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
851 # paramSeq[gapPos_Seq_ToReplace] = paramGL[gapPos_Seq_ToReplace]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
852 #}
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
853
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
854 analyzeMutations2NucUri( matrix(c( paramGL, paramSeq ),2,length(paramGL),byrow=T) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
855
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
856 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
857
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
858 #1 = GL
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
859 #2 = Seq
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
860 analyzeMutations2NucUri <- function( in_matrix=matrix(c(c("A","A","A","C","C","C"),c("A","G","G","C","C","A")),2,6,byrow=T) ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
861 paramGL = in_matrix[2,]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
862 paramSeq = in_matrix[1,]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
863 paramSeqUri = paramGL
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
864 #mutations = apply(rbind(paramGL,paramSeq), 2, function(x){!x[1]==x[2]})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
865 mutations_val = paramGL != paramSeq
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
866 if(any(mutations_val)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
867 mutationPos = {1:length(mutations_val)}[mutations_val]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
868 mutationPos = mutationPos[sapply(mutationPos, function(x){!any(paramSeq[getCodonPos(x)]=="N")})]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
869 length_mutations =length(mutationPos)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
870 mutationInfo = rep(NA,length_mutations)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
871 if(any(mutationPos)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
872
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
873 pos<- mutationPos
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
874 pos_array<-array(sapply(pos,getCodonPos))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
875 codonGL = paramGL[pos_array]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
876
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
877 codonSeq = sapply(pos,function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
878 seqP = paramGL[getCodonPos(x)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
879 muCodonPos = {x-1}%%3+1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
880 seqP[muCodonPos] = paramSeq[x]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
881 return(seqP)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
882 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
883 GLcodons = apply(matrix(codonGL,length_mutations,3,byrow=TRUE),1,c2s)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
884 Seqcodons = apply(codonSeq,2,c2s)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
885 mutationInfo = apply(rbind(GLcodons , Seqcodons),2,function(x){mutationType(c2s(x[1]),c2s(x[2]))})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
886 names(mutationInfo) = mutationPos
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
887 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
888 if(any(!is.na(mutationInfo))){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
889 return(mutationInfo[!is.na(mutationInfo)])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
890 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
891 return(NA)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
892 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
893
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
894
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
895 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
896 return (NA)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
897 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
898 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
899
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
900 processNucMutations2 <- function(mu){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
901 if(!is.na(mu)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
902 #R
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
903 if(any(mu=="R")){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
904 Rs = mu[mu=="R"]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
905 nucNumbs = as.numeric(names(Rs))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
906 R_CDR = sum(as.integer(CDR_Nuc[nucNumbs]),na.rm=T)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
907 R_FWR = sum(as.integer(FWR_Nuc[nucNumbs]),na.rm=T)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
908 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
909 R_CDR = 0
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
910 R_FWR = 0
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
911 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
912
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
913 #S
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
914 if(any(mu=="S")){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
915 Ss = mu[mu=="S"]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
916 nucNumbs = as.numeric(names(Ss))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
917 S_CDR = sum(as.integer(CDR_Nuc[nucNumbs]),na.rm=T)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
918 S_FWR = sum(as.integer(FWR_Nuc[nucNumbs]),na.rm=T)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
919 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
920 S_CDR = 0
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
921 S_FWR = 0
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
922 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
923
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
924
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
925 retVec = c(R_CDR,S_CDR,R_FWR,S_FWR)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
926 retVec[is.na(retVec)]=0
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
927 return(retVec)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
928 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
929 return(rep(0,4))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
930 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
931 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
932
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
933
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
934 ## Z-score Test
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
935 computeZScore <- function(mat, test="Focused"){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
936 matRes <- matrix(NA,ncol=2,nrow=(nrow(mat)))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
937 if(test=="Focused"){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
938 #Z_Focused_CDR
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
939 #P_Denom = sum( mat[1,c(5,6,8)], na.rm=T )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
940 P = apply(mat[,c(5,6,8)],1,function(x){(x[1]/sum(x))})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
941 R_mean = apply(cbind(mat[,c(1,2,4)],P),1,function(x){x[4]*(sum(x[1:3]))})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
942 R_sd=sqrt(R_mean*(1-P))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
943 matRes[,1] = (mat[,1]-R_mean)/R_sd
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
944
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
945 #Z_Focused_FWR
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
946 #P_Denom = sum( mat[1,c(7,6,8)], na.rm=T )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
947 P = apply(mat[,c(7,6,8)],1,function(x){(x[1]/sum(x))})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
948 R_mean = apply(cbind(mat[,c(3,2,4)],P),1,function(x){x[4]*(sum(x[1:3]))})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
949 R_sd=sqrt(R_mean*(1-P))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
950 matRes[,2] = (mat[,3]-R_mean)/R_sd
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
951 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
952
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
953 if(test=="Local"){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
954 #Z_Focused_CDR
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
955 #P_Denom = sum( mat[1,c(5,6,8)], na.rm=T )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
956 P = apply(mat[,c(5,6)],1,function(x){(x[1]/sum(x))})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
957 R_mean = apply(cbind(mat[,c(1,2)],P),1,function(x){x[3]*(sum(x[1:2]))})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
958 R_sd=sqrt(R_mean*(1-P))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
959 matRes[,1] = (mat[,1]-R_mean)/R_sd
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
960
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
961 #Z_Focused_FWR
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
962 #P_Denom = sum( mat[1,c(7,6,8)], na.rm=T )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
963 P = apply(mat[,c(7,8)],1,function(x){(x[1]/sum(x))})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
964 R_mean = apply(cbind(mat[,c(3,4)],P),1,function(x){x[3]*(sum(x[1:2]))})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
965 R_sd=sqrt(R_mean*(1-P))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
966 matRes[,2] = (mat[,3]-R_mean)/R_sd
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
967 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
968
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
969 if(test=="Imbalanced"){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
970 #Z_Focused_CDR
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
971 #P_Denom = sum( mat[1,c(5,6,8)], na.rm=T )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
972 P = apply(mat[,5:8],1,function(x){((x[1]+x[2])/sum(x))})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
973 R_mean = apply(cbind(mat[,1:4],P),1,function(x){x[5]*(sum(x[1:4]))})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
974 R_sd=sqrt(R_mean*(1-P))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
975 matRes[,1] = (mat[,1]-R_mean)/R_sd
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
976
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
977 #Z_Focused_FWR
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
978 #P_Denom = sum( mat[1,c(7,6,8)], na.rm=T )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
979 P = apply(mat[,5:8],1,function(x){((x[3]+x[4])/sum(x))})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
980 R_mean = apply(cbind(mat[,1:4],P),1,function(x){x[5]*(sum(x[1:4]))})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
981 R_sd=sqrt(R_mean*(1-P))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
982 matRes[,2] = (mat[,3]-R_mean)/R_sd
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
983 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
984
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
985 matRes[is.nan(matRes)] = NA
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
986 return(matRes)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
987 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
988
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
989 # Return a p-value for a z-score
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
990 z2p <- function(z){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
991 p=NA
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
992 if( !is.nan(z) && !is.na(z)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
993 if(z>0){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
994 p = (1 - pnorm(z,0,1))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
995 } else if(z<0){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
996 p = (-1 * pnorm(z,0,1))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
997 } else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
998 p = 0.5
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
999 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1000 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1001 p = NA
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1002 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1003 return(p)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1004 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1005
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1006
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1007 ## Bayesian Test
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1008
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1009 # Fitted parameter for the bayesian framework
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1010 BAYESIAN_FITTED<-c(0.407277142798302, 0.554007336744485, 0.63777155771234, 0.693989162719009, 0.735450014674917, 0.767972534429806, 0.794557287143399, 0.816906816601605, 0.83606796225341, 0.852729446430296, 0.867370424541641, 0.880339760590323, 0.891900995024999, 0.902259181289864, 0.911577919359,0.919990301665853, 0.927606458124537, 0.934518806350661, 0.940805863754375, 0.946534836475715, 0.951763691199255, 0.95654428191308, 0.960920179487397, 0.964930893680829, 0.968611312149038, 0.971992459313836, 0.975102110004818, 0.977964943023096, 0.980603428208439, 0.983037660179428, 0.985285800977406, 0.987364285326685, 0.989288037855441, 0.991070478823525, 0.992723699729969, 0.994259575477392, 0.995687688867975, 0.997017365051493, 0.998257085153047, 0.999414558305388, 1.00049681357804, 1.00151036237481, 1.00246080204981, 1.00335370751909, 1.0041939329768, 1.0049859393417, 1.00573382091263, 1.00644127217376, 1.00711179729107, 1.00774845526417, 1.00835412715854, 1.00893143010366, 1.00948275846309, 1.01001030293661, 1.01051606798079, 1.01100188771288, 1.01146944044216, 1.01192026195449, 1.01235575766094, 1.01277721370986)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1011 CONST_i <- sort(c(((2^(seq(-39,0,length.out=201)))/2)[1:200],(c(0:11,13:99)+0.5)/100,1-(2^(seq(-39,0,length.out=201)))/2))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1012
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1013 # Given x, M & p, returns a pdf
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1014 calculate_bayes <- function ( x=3, N=10, p=0.33,
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1015 i=CONST_i,
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1016 max_sigma=20,length_sigma=4001
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1017 ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1018 if(!0%in%N){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1019 G <- max(length(x),length(N),length(p))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1020 x=array(x,dim=G)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1021 N=array(N,dim=G)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1022 p=array(p,dim=G)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1023 sigma_s<-seq(-max_sigma,max_sigma,length.out=length_sigma)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1024 sigma_1<-log({i/{1-i}}/{p/{1-p}})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1025 index<-min(N,60)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1026 y<-dbeta(i,x+BAYESIAN_FITTED[index],N+BAYESIAN_FITTED[index]-x)*(1-p)*p*exp(sigma_1)/({1-p}^2+2*p*{1-p}*exp(sigma_1)+{p^2}*exp(2*sigma_1))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1027 if(!sum(is.na(y))){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1028 tmp<-approx(sigma_1,y,sigma_s)$y
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1029 tmp/sum(tmp)/{2*max_sigma/{length_sigma-1}}
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1030 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1031 return(NA)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1032 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1033 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1034 return(NA)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1035 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1036 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1037 # Given a mat of observed & expected, return a list of CDR & FWR pdf for selection
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1038 computeBayesianScore <- function(mat, test="Focused", max_sigma=20,length_sigma=4001){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1039 flagOneSeq = F
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1040 if(nrow(mat)==1){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1041 mat=rbind(mat,mat)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1042 flagOneSeq = T
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1043 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1044 if(test=="Focused"){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1045 #CDR
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1046 P = c(apply(mat[,c(5,6,8)],1,function(x){(x[1]/sum(x))}),0.5)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1047 N = c(apply(mat[,c(1,2,4)],1,function(x){(sum(x))}),0)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1048 X = c(mat[,1],0)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1049 bayesCDR = apply(cbind(X,N,P),1,function(x){calculate_bayes(x=x[1],N=x[2],p=x[3],max_sigma=max_sigma,length_sigma=length_sigma)})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1050 bayesCDR = bayesCDR[-length(bayesCDR)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1051
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1052 #FWR
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1053 P = c(apply(mat[,c(7,6,8)],1,function(x){(x[1]/sum(x))}),0.5)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1054 N = c(apply(mat[,c(3,2,4)],1,function(x){(sum(x))}),0)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1055 X = c(mat[,3],0)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1056 bayesFWR = apply(cbind(X,N,P),1,function(x){calculate_bayes(x=x[1],N=x[2],p=x[3],max_sigma=max_sigma,length_sigma=length_sigma)})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1057 bayesFWR = bayesFWR[-length(bayesFWR)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1058 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1059
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1060 if(test=="Local"){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1061 #CDR
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1062 P = c(apply(mat[,c(5,6)],1,function(x){(x[1]/sum(x))}),0.5)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1063 N = c(apply(mat[,c(1,2)],1,function(x){(sum(x))}),0)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1064 X = c(mat[,1],0)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1065 bayesCDR = apply(cbind(X,N,P),1,function(x){calculate_bayes(x=x[1],N=x[2],p=x[3],max_sigma=max_sigma,length_sigma=length_sigma)})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1066 bayesCDR = bayesCDR[-length(bayesCDR)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1067
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1068 #FWR
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1069 P = c(apply(mat[,c(7,8)],1,function(x){(x[1]/sum(x))}),0.5)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1070 N = c(apply(mat[,c(3,4)],1,function(x){(sum(x))}),0)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1071 X = c(mat[,3],0)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1072 bayesFWR = apply(cbind(X,N,P),1,function(x){calculate_bayes(x=x[1],N=x[2],p=x[3],max_sigma=max_sigma,length_sigma=length_sigma)})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1073 bayesFWR = bayesFWR[-length(bayesFWR)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1074 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1075
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1076 if(test=="Imbalanced"){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1077 #CDR
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1078 P = c(apply(mat[,c(5:8)],1,function(x){((x[1]+x[2])/sum(x))}),0.5)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1079 N = c(apply(mat[,c(1:4)],1,function(x){(sum(x))}),0)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1080 X = c(apply(mat[,c(1:2)],1,function(x){(sum(x))}),0)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1081 bayesCDR = apply(cbind(X,N,P),1,function(x){calculate_bayes(x=x[1],N=x[2],p=x[3],max_sigma=max_sigma,length_sigma=length_sigma)})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1082 bayesCDR = bayesCDR[-length(bayesCDR)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1083
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1084 #FWR
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1085 P = c(apply(mat[,c(5:8)],1,function(x){((x[3]+x[4])/sum(x))}),0.5)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1086 N = c(apply(mat[,c(1:4)],1,function(x){(sum(x))}),0)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1087 X = c(apply(mat[,c(3:4)],1,function(x){(sum(x))}),0)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1088 bayesFWR = apply(cbind(X,N,P),1,function(x){calculate_bayes(x=x[1],N=x[2],p=x[3],max_sigma=max_sigma,length_sigma=length_sigma)})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1089 bayesFWR = bayesFWR[-length(bayesFWR)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1090 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1091
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1092 if(test=="ImbalancedSilent"){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1093 #CDR
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1094 P = c(apply(mat[,c(6,8)],1,function(x){((x[1])/sum(x))}),0.5)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1095 N = c(apply(mat[,c(2,4)],1,function(x){(sum(x))}),0)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1096 X = c(apply(mat[,c(2,4)],1,function(x){(x[1])}),0)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1097 bayesCDR = apply(cbind(X,N,P),1,function(x){calculate_bayes(x=x[1],N=x[2],p=x[3],max_sigma=max_sigma,length_sigma=length_sigma)})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1098 bayesCDR = bayesCDR[-length(bayesCDR)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1099
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1100 #FWR
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1101 P = c(apply(mat[,c(6,8)],1,function(x){((x[2])/sum(x))}),0.5)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1102 N = c(apply(mat[,c(2,4)],1,function(x){(sum(x))}),0)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1103 X = c(apply(mat[,c(2,4)],1,function(x){(x[2])}),0)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1104 bayesFWR = apply(cbind(X,N,P),1,function(x){calculate_bayes(x=x[1],N=x[2],p=x[3],max_sigma=max_sigma,length_sigma=length_sigma)})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1105 bayesFWR = bayesFWR[-length(bayesFWR)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1106 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1107
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1108 if(flagOneSeq==T){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1109 bayesCDR = bayesCDR[1]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1110 bayesFWR = bayesFWR[1]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1111 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1112 return( list("CDR"=bayesCDR, "FWR"=bayesFWR) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1113 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1114
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1115 ##Covolution
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1116 break2chunks<-function(G=1000){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1117 base<-2^round(log(sqrt(G),2),0)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1118 return(c(rep(base,floor(G/base)-1),base+G-(floor(G/base)*base)))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1119 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1120
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1121 PowersOfTwo <- function(G=100){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1122 exponents <- array()
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1123 i = 0
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1124 while(G > 0){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1125 i=i+1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1126 exponents[i] <- floor( log2(G) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1127 G <- G-2^exponents[i]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1128 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1129 return(exponents)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1130 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1131
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1132 convolutionPowersOfTwo <- function( cons, length_sigma=4001 ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1133 G = ncol(cons)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1134 if(G>1){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1135 for(gen in log(G,2):1){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1136 ll<-seq(from=2,to=2^gen,by=2)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1137 sapply(ll,function(l){cons[,l/2]<<-weighted_conv(cons[,l],cons[,l-1],length_sigma=length_sigma)})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1138 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1139 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1140 return( cons[,1] )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1141 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1142
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1143 convolutionPowersOfTwoByTwos <- function( cons, length_sigma=4001,G=1 ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1144 if(length(ncol(cons))) G<-ncol(cons)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1145 groups <- PowersOfTwo(G)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1146 matG <- matrix(NA, ncol=length(groups), nrow=length(cons)/G )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1147 startIndex = 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1148 for( i in 1:length(groups) ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1149 stopIndex <- 2^groups[i] + startIndex - 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1150 if(stopIndex!=startIndex){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1151 matG[,i] <- convolutionPowersOfTwo( cons[,startIndex:stopIndex], length_sigma=length_sigma )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1152 startIndex = stopIndex + 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1153 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1154 else {
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1155 if(G>1) matG[,i] <- cons[,startIndex:stopIndex]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1156 else matG[,i] <- cons
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1157 #startIndex = stopIndex + 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1158 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1159 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1160 return( list( matG, groups ) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1161 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1162
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1163 weighted_conv<-function(x,y,w=1,m=100,length_sigma=4001){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1164 lx<-length(x)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1165 ly<-length(y)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1166 if({lx<m}| {{lx*w}<m}| {{ly}<m}| {{ly*w}<m}){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1167 if(w<1){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1168 y1<-approx(1:ly,y,seq(1,ly,length.out=m))$y
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1169 x1<-approx(1:lx,x,seq(1,lx,length.out=m/w))$y
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1170 lx<-length(x1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1171 ly<-length(y1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1172 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1173 else {
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1174 y1<-approx(1:ly,y,seq(1,ly,length.out=m*w))$y
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1175 x1<-approx(1:lx,x,seq(1,lx,length.out=m))$y
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1176 lx<-length(x1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1177 ly<-length(y1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1178 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1179 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1180 else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1181 x1<-x
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1182 y1<-approx(1:ly,y,seq(1,ly,length.out=floor(lx*w)))$y
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1183 ly<-length(y1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1184 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1185 tmp<-approx(x=1:(lx+ly-1),y=convolve(x1,rev(y1),type="open"),xout=seq(1,lx+ly-1,length.out=length_sigma))$y
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1186 tmp[tmp<=0] = 0
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1187 return(tmp/sum(tmp))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1188 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1189
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1190 calculate_bayesGHelper <- function( listMatG,length_sigma=4001 ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1191 matG <- listMatG[[1]]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1192 groups <- listMatG[[2]]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1193 i = 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1194 resConv <- matG[,i]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1195 denom <- 2^groups[i]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1196 if(length(groups)>1){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1197 while( i<length(groups) ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1198 i = i + 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1199 resConv <- weighted_conv(resConv, matG[,i], w= {{2^groups[i]}/denom} ,length_sigma=length_sigma)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1200 #cat({{2^groups[i]}/denom},"\n")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1201 denom <- denom + 2^groups[i]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1202 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1203 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1204 return(resConv)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1205 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1206
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1207 # Given a list of PDFs, returns a convoluted PDF
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1208 groupPosteriors <- function( listPosteriors, max_sigma=20, length_sigma=4001 ,Threshold=2 ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1209 listPosteriors = listPosteriors[ !is.na(listPosteriors) ]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1210 Length_Postrior<-length(listPosteriors)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1211 if(Length_Postrior>1 & Length_Postrior<=Threshold){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1212 cons = matrix(unlist(listPosteriors),length(listPosteriors[[1]]),length(listPosteriors))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1213 listMatG <- convolutionPowersOfTwoByTwos(cons,length_sigma=length_sigma)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1214 y<-calculate_bayesGHelper(listMatG,length_sigma=length_sigma)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1215 return( y/sum(y)/(2*max_sigma/(length_sigma-1)) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1216 }else if(Length_Postrior==1) return(listPosteriors[[1]])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1217 else if(Length_Postrior==0) return(NA)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1218 else {
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1219 cons = matrix(unlist(listPosteriors),length(listPosteriors[[1]]),length(listPosteriors))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1220 y = fastConv(cons,max_sigma=max_sigma, length_sigma=length_sigma )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1221 return( y/sum(y)/(2*max_sigma/(length_sigma-1)) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1222 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1223 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1224
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1225 fastConv<-function(cons, max_sigma=20, length_sigma=4001){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1226 chunks<-break2chunks(G=ncol(cons))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1227 if(ncol(cons)==3) chunks<-2:1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1228 index_chunks_end <- cumsum(chunks)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1229 index_chunks_start <- c(1,index_chunks_end[-length(index_chunks_end)]+1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1230 index_chunks <- cbind(index_chunks_start,index_chunks_end)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1231
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1232 case <- sum(chunks!=chunks[1])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1233 if(case==1) End <- max(1,((length(index_chunks)/2)-1))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1234 else End <- max(1,((length(index_chunks)/2)))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1235
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1236 firsts <- sapply(1:End,function(i){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1237 indexes<-index_chunks[i,1]:index_chunks[i,2]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1238 convolutionPowersOfTwoByTwos(cons[ ,indexes])[[1]]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1239 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1240 if(case==0){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1241 result<-calculate_bayesGHelper( convolutionPowersOfTwoByTwos(firsts) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1242 }else if(case==1){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1243 last<-list(calculate_bayesGHelper(
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1244 convolutionPowersOfTwoByTwos( cons[ ,index_chunks[length(index_chunks)/2,1]:index_chunks[length(index_chunks)/2,2]] )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1245 ),0)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1246 result_first<-calculate_bayesGHelper(convolutionPowersOfTwoByTwos(firsts))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1247 result<-calculate_bayesGHelper(
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1248 list(
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1249 cbind(
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1250 result_first,last[[1]]),
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1251 c(log(index_chunks_end[length(index_chunks)/2-1],2),log(index_chunks[length(index_chunks)/2,2]-index_chunks[length(index_chunks)/2,1]+1,2))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1252 )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1253 )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1254 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1255 return(as.vector(result))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1256 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1257
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1258 # Computes the 95% CI for a pdf
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1259 calcBayesCI <- function(Pdf,low=0.025,up=0.975,max_sigma=20, length_sigma=4001){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1260 if(length(Pdf)!=length_sigma) return(NA)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1261 sigma_s=seq(-max_sigma,max_sigma,length.out=length_sigma)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1262 cdf = cumsum(Pdf)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1263 cdf = cdf/cdf[length(cdf)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1264 return( c(sigma_s[findInterval(low,cdf)-1] , sigma_s[findInterval(up,cdf)]) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1265 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1266
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1267 # Computes a mean for a pdf
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1268 calcBayesMean <- function(Pdf,max_sigma=20,length_sigma=4001){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1269 if(length(Pdf)!=length_sigma) return(NA)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1270 sigma_s=seq(-max_sigma,max_sigma,length.out=length_sigma)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1271 norm = {length_sigma-1}/2/max_sigma
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1272 return( (Pdf%*%sigma_s/norm) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1273 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1274
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1275 # Returns the mean, and the 95% CI for a pdf
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1276 calcBayesOutputInfo <- function(Pdf,low=0.025,up=0.975,max_sigma=20, length_sigma=4001){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1277 if(is.na(Pdf))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1278 return(rep(NA,3))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1279 bCI = calcBayesCI(Pdf=Pdf,low=low,up=up,max_sigma=max_sigma,length_sigma=length_sigma)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1280 bMean = calcBayesMean(Pdf=Pdf,max_sigma=max_sigma,length_sigma=length_sigma)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1281 return(c(bMean, bCI))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1282 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1283
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1284 # Computes the p-value of a pdf
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1285 computeSigmaP <- function(Pdf, length_sigma=4001, max_sigma=20){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1286 if(length(Pdf)>1){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1287 norm = {length_sigma-1}/2/max_sigma
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1288 pVal = {sum(Pdf[1:{{length_sigma-1}/2}]) + Pdf[{{length_sigma+1}/2}]/2}/norm
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1289 if(pVal>0.5){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1290 pVal = pVal-1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1291 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1292 return(pVal)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1293 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1294 return(NA)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1295 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1296 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1297
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1298 # Compute p-value of two distributions
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1299 compareTwoDistsFaster <-function(sigma_S=seq(-20,20,length.out=4001), N=10000, dens1=runif(4001,0,1), dens2=runif(4001,0,1)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1300 #print(c(length(dens1),length(dens2)))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1301 if(length(dens1)>1 & length(dens2)>1 ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1302 dens1<-dens1/sum(dens1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1303 dens2<-dens2/sum(dens2)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1304 cum2 <- cumsum(dens2)-dens2/2
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1305 tmp<- sum(sapply(1:length(dens1),function(i)return(dens1[i]*cum2[i])))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1306 #print(tmp)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1307 if(tmp>0.5)tmp<-tmp-1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1308 return( tmp )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1309 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1310 else {
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1311 return(NA)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1312 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1313 #return (sum(sapply(1:N,function(i)(sample(sigma_S,1,prob=dens1)>sample(sigma_S,1,prob=dens2))))/N)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1314 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1315
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1316 # get number of seqeunces contributing to the sigma (i.e. seqeunces with mutations)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1317 numberOfSeqsWithMutations <- function(matMutations,test=1){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1318 if(test==4)test=2
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1319 cdrSeqs <- 0
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1320 fwrSeqs <- 0
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1321 if(test==1){#focused
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1322 cdrMutations <- apply(matMutations, 1, function(x){ sum(x[c(1,2,4)]) })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1323 fwrMutations <- apply(matMutations, 1, function(x){ sum(x[c(3,4,2)]) })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1324 if( any(which(cdrMutations>0)) ) cdrSeqs <- sum(cdrMutations>0)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1325 if( any(which(fwrMutations>0)) ) fwrSeqs <- sum(fwrMutations>0)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1326 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1327 if(test==2){#local
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1328 cdrMutations <- apply(matMutations, 1, function(x){ sum(x[c(1,2)]) })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1329 fwrMutations <- apply(matMutations, 1, function(x){ sum(x[c(3,4)]) })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1330 if( any(which(cdrMutations>0)) ) cdrSeqs <- sum(cdrMutations>0)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1331 if( any(which(fwrMutations>0)) ) fwrSeqs <- sum(fwrMutations>0)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1332 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1333 return(c("CDR"=cdrSeqs, "FWR"=fwrSeqs))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1334 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1335
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1336
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1337
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1338 shadeColor <- function(sigmaVal=NA,pVal=NA){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1339 if(is.na(sigmaVal) & is.na(pVal)) return(NA)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1340 if(is.na(sigmaVal) & !is.na(pVal)) sigmaVal=sign(pVal)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1341 if(is.na(pVal) || pVal==1 || pVal==0){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1342 returnColor = "#FFFFFF";
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1343 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1344 colVal=abs(pVal);
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1345
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1346 if(sigmaVal<0){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1347 if(colVal>0.1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1348 returnColor = "#CCFFCC";
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1349 if(colVal<=0.1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1350 returnColor = "#99FF99";
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1351 if(colVal<=0.050)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1352 returnColor = "#66FF66";
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1353 if(colVal<=0.010)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1354 returnColor = "#33FF33";
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1355 if(colVal<=0.005)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1356 returnColor = "#00FF00";
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1357
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1358 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1359 if(colVal>0.1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1360 returnColor = "#FFCCCC";
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1361 if(colVal<=0.1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1362 returnColor = "#FF9999";
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1363 if(colVal<=0.05)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1364 returnColor = "#FF6666";
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1365 if(colVal<=0.01)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1366 returnColor = "#FF3333";
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1367 if(colVal<0.005)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1368 returnColor = "#FF0000";
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1369 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1370 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1371
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1372 return(returnColor)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1373 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1374
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1375
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1376
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1377 plotHelp <- function(xfrac=0.05,yfrac=0.05,log=FALSE){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1378 if(!log){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1379 x = par()$usr[1]-(par()$usr[2]-par()$usr[1])*xfrac
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1380 y = par()$usr[4]+(par()$usr[4]-par()$usr[3])*yfrac
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1381 }else {
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1382 if(log==2){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1383 x = par()$usr[1]-(par()$usr[2]-par()$usr[1])*xfrac
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1384 y = 10^((par()$usr[4])+((par()$usr[4])-(par()$usr[3]))*yfrac)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1385 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1386 if(log==1){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1387 x = 10^((par()$usr[1])-((par()$usr[2])-(par()$usr[1]))*xfrac)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1388 y = par()$usr[4]+(par()$usr[4]-par()$usr[3])*yfrac
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1389 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1390 if(log==3){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1391 x = 10^((par()$usr[1])-((par()$usr[2])-(par()$usr[1]))*xfrac)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1392 y = 10^((par()$usr[4])+((par()$usr[4])-(par()$usr[3]))*yfrac)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1393 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1394 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1395 return(c("x"=x,"y"=y))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1396 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1397
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1398 # SHMulation
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1399
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1400 # Based on targeting, introduce a single mutation & then update the targeting
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1401 oneMutation <- function(){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1402 # Pick a postion + mutation
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1403 posMutation = sample(1:(seqGermlineLen*4),1,replace=F,prob=as.vector(seqTargeting))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1404 posNucNumb = ceiling(posMutation/4) # Nucleotide number
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1405 posNucKind = 4 - ( (posNucNumb*4) - posMutation ) # Nuc the position mutates to
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1406
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1407 #mutate the simulation sequence
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1408 seqSimVec <- s2c(seqSim)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1409 seqSimVec[posNucNumb] <- NUCLEOTIDES[posNucKind]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1410 seqSim <<- c2s(seqSimVec)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1411
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1412 #update Mutability, Targeting & MutationsTypes
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1413 updateMutabilityNTargeting(posNucNumb)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1414
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1415 #return(c(posNucNumb,NUCLEOTIDES[posNucKind]))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1416 return(posNucNumb)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1417 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1418
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1419 updateMutabilityNTargeting <- function(position){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1420 min_i<-max((position-2),1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1421 max_i<-min((position+2),nchar(seqSim))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1422 min_ii<-min(min_i,3)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1423
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1424 #mutability - update locally
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1425 seqMutability[(min_i):(max_i)] <<- computeMutabilities(substr(seqSim,position-4,position+4))[(min_ii):(max_i-min_i+min_ii)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1426
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1427
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1428 #targeting - compute locally
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1429 seqTargeting[,min_i:max_i] <<- computeTargeting(substr(seqSim,min_i,max_i),seqMutability[min_i:max_i])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1430 seqTargeting[is.na(seqTargeting)] <<- 0
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1431 #mutCodonPos = getCodonPos(position)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1432 mutCodonPos = seq(getCodonPos(min_i)[1],getCodonPos(max_i)[3])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1433 #cat(mutCodonPos,"\n")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1434 mutTypeCodon = getCodonPos(position)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1435 seqMutationTypes[,mutTypeCodon] <<- computeMutationTypesFast( substr(seqSim,mutTypeCodon[1],mutTypeCodon[3]) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1436 # Stop = 0
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1437 if(any(seqMutationTypes[,mutCodonPos]=="Stop",na.rm=T )){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1438 seqTargeting[,mutCodonPos][seqMutationTypes[,mutCodonPos]=="Stop"] <<- 0
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1439 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1440
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1441
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1442 #Selection
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1443 selectedPos = (min_i*4-4)+(which(seqMutationTypes[,min_i:max_i]=="R"))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1444 # CDR
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1445 selectedCDR = selectedPos[which(matCDR[selectedPos]==T)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1446 seqTargeting[selectedCDR] <<- seqTargeting[selectedCDR] * exp(selCDR)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1447 seqTargeting[selectedCDR] <<- seqTargeting[selectedCDR]/baseLineCDR_K
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1448
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1449 # FWR
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1450 selectedFWR = selectedPos[which(matFWR[selectedPos]==T)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1451 seqTargeting[selectedFWR] <<- seqTargeting[selectedFWR] * exp(selFWR)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1452 seqTargeting[selectedFWR] <<- seqTargeting[selectedFWR]/baseLineFWR_K
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1453
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1454 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1455
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1456
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1457
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1458 # Validate the mutation: if the mutation has not been sampled before validate it, else discard it.
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1459 validateMutation <- function(){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1460 if( !(mutatedPos%in%mutatedPositions) ){ # if it's a new mutation
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1461 uniqueMutationsIntroduced <<- uniqueMutationsIntroduced + 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1462 mutatedPositions[uniqueMutationsIntroduced] <<- mutatedPos
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1463 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1464 if(substr(seqSim,mutatedPos,mutatedPos)==substr(seqGermline,mutatedPos,mutatedPos)){ # back to germline mutation
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1465 mutatedPositions <<- mutatedPositions[-which(mutatedPositions==mutatedPos)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1466 uniqueMutationsIntroduced <<- uniqueMutationsIntroduced - 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1467 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1468 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1469 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1470
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1471
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1472
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1473 # Places text (labels) at normalized coordinates
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1474 myaxis <- function(xfrac=0.05,yfrac=0.05,log=FALSE,w="text",cex=1,adj=1,thecol="black"){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1475 par(xpd=TRUE)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1476 if(!log)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1477 text(par()$usr[1]-(par()$usr[2]-par()$usr[1])*xfrac,par()$usr[4]+(par()$usr[4]-par()$usr[3])*yfrac,w,cex=cex,adj=adj,col=thecol)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1478 else {
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1479 if(log==2)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1480 text(
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1481 par()$usr[1]-(par()$usr[2]-par()$usr[1])*xfrac,
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1482 10^((par()$usr[4])+((par()$usr[4])-(par()$usr[3]))*yfrac),
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1483 w,cex=cex,adj=adj,col=thecol)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1484 if(log==1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1485 text(
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1486 10^((par()$usr[1])-((par()$usr[2])-(par()$usr[1]))*xfrac),
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1487 par()$usr[4]+(par()$usr[4]-par()$usr[3])*yfrac,
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1488 w,cex=cex,adj=adj,col=thecol)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1489 if(log==3)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1490 text(
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1491 10^((par()$usr[1])-((par()$usr[2])-(par()$usr[1]))*xfrac),
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1492 10^((par()$usr[4])+((par()$usr[4])-(par()$usr[3]))*yfrac),
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1493 w,cex=cex,adj=adj,col=thecol)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1494 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1495 par(xpd=FALSE)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1496 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1497
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1498
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1499
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1500 # Count the mutations in a sequence
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1501 analyzeMutations <- function( inputMatrixIndex, model = 0 , multipleMutation=0, seqWithStops=0){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1502
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1503 paramGL = s2c(matInput[inputMatrixIndex,2])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1504 paramSeq = s2c(matInput[inputMatrixIndex,1])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1505
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1506 #if( any(paramSeq=="N") ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1507 # gapPos_Seq = which(paramSeq=="N")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1508 # gapPos_Seq_ToReplace = gapPos_Seq[paramGL[gapPos_Seq] != "N"]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1509 # paramSeq[gapPos_Seq_ToReplace] = paramGL[gapPos_Seq_ToReplace]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1510 #}
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1511 mutations_val = paramGL != paramSeq
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1512
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1513 if(any(mutations_val)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1514 mutationPos = which(mutations_val)#{1:length(mutations_val)}[mutations_val]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1515 length_mutations =length(mutationPos)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1516 mutationInfo = rep(NA,length_mutations)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1517
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1518 pos<- mutationPos
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1519 pos_array<-array(sapply(pos,getCodonPos))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1520 codonGL = paramGL[pos_array]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1521 codonSeqWhole = paramSeq[pos_array]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1522 codonSeq = sapply(pos,function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1523 seqP = paramGL[getCodonPos(x)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1524 muCodonPos = {x-1}%%3+1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1525 seqP[muCodonPos] = paramSeq[x]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1526 return(seqP)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1527 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1528 GLcodons = apply(matrix(codonGL,length_mutations,3,byrow=TRUE),1,c2s)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1529 SeqcodonsWhole = apply(matrix(codonSeqWhole,length_mutations,3,byrow=TRUE),1,c2s)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1530 Seqcodons = apply(codonSeq,2,c2s)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1531
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1532 mutationInfo = apply(rbind(GLcodons , Seqcodons),2,function(x){mutationType(c2s(x[1]),c2s(x[2]))})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1533 names(mutationInfo) = mutationPos
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1534
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1535 mutationInfoWhole = apply(rbind(GLcodons , SeqcodonsWhole),2,function(x){mutationType(c2s(x[1]),c2s(x[2]))})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1536 names(mutationInfoWhole) = mutationPos
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1537
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1538 mutationInfo <- mutationInfo[!is.na(mutationInfo)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1539 mutationInfoWhole <- mutationInfoWhole[!is.na(mutationInfoWhole)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1540
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1541 if(any(!is.na(mutationInfo))){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1542
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1543 #Filter based on Stop (at the codon level)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1544 if(seqWithStops==1){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1545 nucleotidesAtStopCodons = names(mutationInfoWhole[mutationInfoWhole!="Stop"])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1546 mutationInfo = mutationInfo[nucleotidesAtStopCodons]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1547 mutationInfoWhole = mutationInfo[nucleotidesAtStopCodons]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1548 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1549 countStops = sum(mutationInfoWhole=="Stop")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1550 if(seqWithStops==2 & countStops==0) mutationInfo = NA
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1551 if(seqWithStops==3 & countStops>0) mutationInfo = NA
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1552 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1553
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1554 if(any(!is.na(mutationInfo))){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1555 #Filter mutations based on multipleMutation
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1556 if(multipleMutation==1 & !is.na(mutationInfo)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1557 mutationCodons = getCodonNumb(as.numeric(names(mutationInfoWhole)))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1558 tableMutationCodons <- table(mutationCodons)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1559 codonsWithMultipleMutations <- as.numeric(names(tableMutationCodons[tableMutationCodons>1]))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1560 if(any(codonsWithMultipleMutations)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1561 #remove the nucleotide mutations in the codons with multiple mutations
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1562 mutationInfo <- mutationInfo[!(mutationCodons %in% codonsWithMultipleMutations)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1563 #replace those codons with Ns in the input sequence
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1564 paramSeq[unlist(lapply(codonsWithMultipleMutations, getCodonNucs))] = "N"
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1565 matInput[inputMatrixIndex,1] <<- c2s(paramSeq)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1566 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1567 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1568
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1569 #Filter mutations based on the model
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1570 if(any(mutationInfo)==T | is.na(any(mutationInfo))){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1571
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1572 if(model==1 & !is.na(mutationInfo)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1573 mutationInfo <- mutationInfo[mutationInfo=="S"]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1574 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1575 if(any(mutationInfo)==T | is.na(any(mutationInfo))) return(mutationInfo)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1576 else return(NA)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1577 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1578 return(NA)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1579 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1580 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1581 return(NA)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1582 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1583
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1584
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1585 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1586 return(NA)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1587 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1588
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1589
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1590 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1591 return (NA)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1592 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1593 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1594
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1595 analyzeMutationsFixed <- function( inputArray, model = 0 , multipleMutation=0, seqWithStops=0){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1596
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1597 paramGL = s2c(inputArray[2])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1598 paramSeq = s2c(inputArray[1])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1599 inputSeq <- inputArray[1]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1600 #if( any(paramSeq=="N") ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1601 # gapPos_Seq = which(paramSeq=="N")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1602 # gapPos_Seq_ToReplace = gapPos_Seq[paramGL[gapPos_Seq] != "N"]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1603 # paramSeq[gapPos_Seq_ToReplace] = paramGL[gapPos_Seq_ToReplace]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1604 #}
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1605 mutations_val = paramGL != paramSeq
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1606
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1607 if(any(mutations_val)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1608 mutationPos = which(mutations_val)#{1:length(mutations_val)}[mutations_val]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1609 length_mutations =length(mutationPos)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1610 mutationInfo = rep(NA,length_mutations)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1611
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1612 pos<- mutationPos
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1613 pos_array<-array(sapply(pos,getCodonPos))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1614 codonGL = paramGL[pos_array]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1615 codonSeqWhole = paramSeq[pos_array]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1616 codonSeq = sapply(pos,function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1617 seqP = paramGL[getCodonPos(x)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1618 muCodonPos = {x-1}%%3+1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1619 seqP[muCodonPos] = paramSeq[x]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1620 return(seqP)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1621 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1622 GLcodons = apply(matrix(codonGL,length_mutations,3,byrow=TRUE),1,c2s)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1623 SeqcodonsWhole = apply(matrix(codonSeqWhole,length_mutations,3,byrow=TRUE),1,c2s)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1624 Seqcodons = apply(codonSeq,2,c2s)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1625
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1626 mutationInfo = apply(rbind(GLcodons , Seqcodons),2,function(x){mutationType(c2s(x[1]),c2s(x[2]))})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1627 names(mutationInfo) = mutationPos
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1628
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1629 mutationInfoWhole = apply(rbind(GLcodons , SeqcodonsWhole),2,function(x){mutationType(c2s(x[1]),c2s(x[2]))})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1630 names(mutationInfoWhole) = mutationPos
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1631
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1632 mutationInfo <- mutationInfo[!is.na(mutationInfo)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1633 mutationInfoWhole <- mutationInfoWhole[!is.na(mutationInfoWhole)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1634
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1635 if(any(!is.na(mutationInfo))){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1636
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1637 #Filter based on Stop (at the codon level)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1638 if(seqWithStops==1){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1639 nucleotidesAtStopCodons = names(mutationInfoWhole[mutationInfoWhole!="Stop"])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1640 mutationInfo = mutationInfo[nucleotidesAtStopCodons]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1641 mutationInfoWhole = mutationInfo[nucleotidesAtStopCodons]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1642 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1643 countStops = sum(mutationInfoWhole=="Stop")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1644 if(seqWithStops==2 & countStops==0) mutationInfo = NA
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1645 if(seqWithStops==3 & countStops>0) mutationInfo = NA
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1646 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1647
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1648 if(any(!is.na(mutationInfo))){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1649 #Filter mutations based on multipleMutation
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1650 if(multipleMutation==1 & !is.na(mutationInfo)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1651 mutationCodons = getCodonNumb(as.numeric(names(mutationInfoWhole)))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1652 tableMutationCodons <- table(mutationCodons)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1653 codonsWithMultipleMutations <- as.numeric(names(tableMutationCodons[tableMutationCodons>1]))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1654 if(any(codonsWithMultipleMutations)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1655 #remove the nucleotide mutations in the codons with multiple mutations
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1656 mutationInfo <- mutationInfo[!(mutationCodons %in% codonsWithMultipleMutations)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1657 #replace those codons with Ns in the input sequence
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1658 paramSeq[unlist(lapply(codonsWithMultipleMutations, getCodonNucs))] = "N"
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1659 #matInput[inputMatrixIndex,1] <<- c2s(paramSeq)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1660 inputSeq <- c2s(paramSeq)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1661 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1662 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1663
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1664 #Filter mutations based on the model
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1665 if(any(mutationInfo)==T | is.na(any(mutationInfo))){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1666
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1667 if(model==1 & !is.na(mutationInfo)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1668 mutationInfo <- mutationInfo[mutationInfo=="S"]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1669 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1670 if(any(mutationInfo)==T | is.na(any(mutationInfo))) return(list(mutationInfo,inputSeq))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1671 else return(list(NA,inputSeq))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1672 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1673 return(list(NA,inputSeq))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1674 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1675 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1676 return(list(NA,inputSeq))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1677 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1678
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1679
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1680 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1681 return(list(NA,inputSeq))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1682 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1683
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1684
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1685 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1686 return (list(NA,inputSeq))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1687 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1688 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1689
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1690 # triMutability Background Count
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1691 buildMutabilityModel <- function( inputMatrixIndex, model=0 , multipleMutation=0, seqWithStops=0, stopMutations=0){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1692
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1693 #rowOrigMatInput = matInput[inputMatrixIndex,]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1694 seqGL = gsub("-", "", matInput[inputMatrixIndex,2])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1695 seqInput = gsub("-", "", matInput[inputMatrixIndex,1])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1696 #matInput[inputMatrixIndex,] <<- cbind(seqInput,seqGL)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1697 tempInput <- cbind(seqInput,seqGL)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1698 seqLength = nchar(seqGL)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1699 list_analyzeMutationsFixed<- analyzeMutationsFixed(tempInput, model, multipleMutation, seqWithStops)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1700 mutationCount <- list_analyzeMutationsFixed[[1]]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1701 seqInput <- list_analyzeMutationsFixed[[2]]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1702 BackgroundMatrix = mutabilityMatrix
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1703 MutationMatrix = mutabilityMatrix
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1704 MutationCountMatrix = mutabilityMatrix
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1705 if(!is.na(mutationCount)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1706 if((stopMutations==0 & model==0) | (stopMutations==1 & (sum(mutationCount=="Stop")<length(mutationCount))) | (model==1 & (sum(mutationCount=="S")>0)) ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1707
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1708 fivermerStartPos = 1:(seqLength-4)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1709 fivemerLength <- length(fivermerStartPos)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1710 fivemerGL <- substr(rep(seqGL,length(fivermerStartPos)),(fivermerStartPos),(fivermerStartPos+4))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1711 fivemerSeq <- substr(rep(seqInput,length(fivermerStartPos)),(fivermerStartPos),(fivermerStartPos+4))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1712
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1713 #Background
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1714 for(fivemerIndex in 1:fivemerLength){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1715 fivemer = fivemerGL[fivemerIndex]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1716 if(!any(grep("N",fivemer))){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1717 fivemerCodonPos = fivemerCodon(fivemerIndex)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1718 fivemerReadingFrameCodon = substr(fivemer,fivemerCodonPos[1],fivemerCodonPos[3])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1719 fivemerReadingFrameCodonInputSeq = substr(fivemerSeq[fivemerIndex],fivemerCodonPos[1],fivemerCodonPos[3])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1720
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1721 # All mutations model
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1722 #if(!any(grep("N",fivemerReadingFrameCodon))){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1723 if(model==0){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1724 if(stopMutations==0){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1725 if(!any(grep("N",fivemerReadingFrameCodonInputSeq)))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1726 BackgroundMatrix[fivemer] <- (BackgroundMatrix[fivemer] + 1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1727 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1728 if( !any(grep("N",fivemerReadingFrameCodonInputSeq)) & translateCodonToAminoAcid(fivemerReadingFrameCodon)!="*" ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1729 positionWithinCodon = which(fivemerCodonPos==3)#positionsWithinCodon[(fivemerCodonPos[1]%%3)+1]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1730 BackgroundMatrix[fivemer] <- (BackgroundMatrix[fivemer] + probNonStopMutations[fivemerReadingFrameCodon,positionWithinCodon])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1731 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1732 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1733 }else{ # Only silent mutations
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1734 if( !any(grep("N",fivemerReadingFrameCodonInputSeq)) & translateCodonToAminoAcid(fivemerReadingFrameCodon)!="*" & translateCodonToAminoAcid(fivemerReadingFrameCodonInputSeq)==translateCodonToAminoAcid(fivemerReadingFrameCodon)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1735 positionWithinCodon = which(fivemerCodonPos==3)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1736 BackgroundMatrix[fivemer] <- (BackgroundMatrix[fivemer] + probSMutations[fivemerReadingFrameCodon,positionWithinCodon])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1737 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1738 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1739 #}
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1740 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1741 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1742
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1743 #Mutations
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1744 if(stopMutations==1) mutationCount = mutationCount[mutationCount!="Stop"]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1745 if(model==1) mutationCount = mutationCount[mutationCount=="S"]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1746 mutationPositions = as.numeric(names(mutationCount))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1747 mutationCount = mutationCount[mutationPositions>2 & mutationPositions<(seqLength-1)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1748 mutationPositions = mutationPositions[mutationPositions>2 & mutationPositions<(seqLength-1)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1749 countMutations = 0
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1750 for(mutationPosition in mutationPositions){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1751 fivemerIndex = mutationPosition-2
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1752 fivemer = fivemerSeq[fivemerIndex]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1753 GLfivemer = fivemerGL[fivemerIndex]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1754 fivemerCodonPos = fivemerCodon(fivemerIndex)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1755 fivemerReadingFrameCodon = substr(fivemer,fivemerCodonPos[1],fivemerCodonPos[3])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1756 fivemerReadingFrameCodonGL = substr(GLfivemer,fivemerCodonPos[1],fivemerCodonPos[3])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1757 if(!any(grep("N",fivemer)) & !any(grep("N",GLfivemer))){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1758 if(model==0){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1759 countMutations = countMutations + 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1760 MutationMatrix[GLfivemer] <- (MutationMatrix[GLfivemer] + 1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1761 MutationCountMatrix[GLfivemer] <- (MutationCountMatrix[GLfivemer] + 1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1762 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1763 if( translateCodonToAminoAcid(fivemerReadingFrameCodonGL)!="*" ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1764 countMutations = countMutations + 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1765 positionWithinCodon = which(fivemerCodonPos==3)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1766 glNuc = substr(fivemerReadingFrameCodonGL,positionWithinCodon,positionWithinCodon)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1767 inputNuc = substr(fivemerReadingFrameCodon,positionWithinCodon,positionWithinCodon)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1768 MutationMatrix[GLfivemer] <- (MutationMatrix[GLfivemer] + substitution[glNuc,inputNuc])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1769 MutationCountMatrix[GLfivemer] <- (MutationCountMatrix[GLfivemer] + 1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1770 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1771 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1772 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1773 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1774
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1775 seqMutability = MutationMatrix/BackgroundMatrix
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1776 seqMutability = seqMutability/sum(seqMutability,na.rm=TRUE)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1777 #cat(inputMatrixIndex,"\t",countMutations,"\n")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1778 return(list("seqMutability" = seqMutability,"numbMutations" = countMutations,"seqMutabilityCount" = MutationCountMatrix, "BackgroundMatrix"=BackgroundMatrix))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1779
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1780 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1781 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1782
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1783 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1784
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1785 #Returns the codon position containing the middle nucleotide
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1786 fivemerCodon <- function(fivemerIndex){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1787 codonPos = list(2:4,1:3,3:5)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1788 fivemerType = fivemerIndex%%3
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1789 return(codonPos[[fivemerType+1]])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1790 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1791
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1792 #returns probability values for one mutation in codons resulting in R, S or Stop
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1793 probMutations <- function(typeOfMutation){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1794 matMutationProb <- matrix(0,ncol=3,nrow=125,dimnames=list(words(alphabet = c(NUCLEOTIDES,"N"), length=3),c(1:3)))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1795 for(codon in rownames(matMutationProb)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1796 if( !any(grep("N",codon)) ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1797 for(muPos in 1:3){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1798 matCodon = matrix(rep(s2c(codon),3),nrow=3,ncol=3,byrow=T)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1799 glNuc = matCodon[1,muPos]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1800 matCodon[,muPos] = canMutateTo(glNuc)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1801 substitutionRate = substitution[glNuc,matCodon[,muPos]]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1802 typeOfMutations = apply(rbind(rep(codon,3),apply(matCodon,1,c2s)),2,function(x){mutationType(c2s(x[1]),c2s(x[2]))})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1803 matMutationProb[codon,muPos] <- sum(substitutionRate[typeOfMutations==typeOfMutation])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1804 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1805 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1806 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1807
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1808 return(matMutationProb)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1809 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1810
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1811
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1812
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1813
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1814 #Mapping Trinucleotides to fivemers
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1815 mapTriToFivemer <- function(triMutability=triMutability_Literature_Human){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1816 rownames(triMutability) <- triMutability_Names
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1817 Fivemer<-rep(NA,1024)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1818 names(Fivemer)<-words(alphabet=NUCLEOTIDES,length=5)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1819 Fivemer<-sapply(names(Fivemer),function(Word)return(sum( c(triMutability[substring(Word,3,5),1],triMutability[substring(Word,2,4),2],triMutability[substring(Word,1,3),3]),na.rm=TRUE)))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1820 Fivemer<-Fivemer/sum(Fivemer)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1821 return(Fivemer)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1822 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1823
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1824 collapseFivemerToTri<-function(Fivemer,Weights=MutabilityWeights,position=1,NUC="A"){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1825 Indices<-substring(names(Fivemer),3,3)==NUC
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1826 Factors<-substring(names(Fivemer[Indices]),(4-position),(6-position))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1827 tapply(which(Indices),Factors,function(i)weighted.mean(Fivemer[i],Weights[i],na.rm=TRUE))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1828 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1829
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1830
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1831
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1832 CountFivemerToTri<-function(Fivemer,Weights=MutabilityWeights,position=1,NUC="A"){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1833 Indices<-substring(names(Fivemer),3,3)==NUC
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1834 Factors<-substring(names(Fivemer[Indices]),(4-position),(6-position))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1835 tapply(which(Indices),Factors,function(i)sum(Weights[i],na.rm=TRUE))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1836 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1837
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1838 #Uses the real counts of the mutated fivemers
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1839 CountFivemerToTri2<-function(Fivemer,Counts=MutabilityCounts,position=1,NUC="A"){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1840 Indices<-substring(names(Fivemer),3,3)==NUC
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1841 Factors<-substring(names(Fivemer[Indices]),(4-position),(6-position))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1842 tapply(which(Indices),Factors,function(i)sum(Counts[i],na.rm=TRUE))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1843 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1844
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1845 bootstrap<-function(x=c(33,12,21),M=10000,alpha=0.05){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1846 N<-sum(x)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1847 if(N){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1848 p<-x/N
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1849 k<-length(x)-1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1850 tmp<-rmultinom(M, size = N, prob=p)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1851 tmp_p<-apply(tmp,2,function(y)y/N)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1852 (apply(tmp_p,1,function(y)quantile(y,c(alpha/2/k,1-alpha/2/k))))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1853 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1854 else return(matrix(0,2,length(x)))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1855 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1856
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1857
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1858
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1859
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1860 bootstrap2<-function(x=c(33,12,21),n=10,M=10000,alpha=0.05){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1861
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1862 N<-sum(x)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1863 k<-length(x)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1864 y<-rep(1:k,x)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1865 tmp<-sapply(1:M,function(i)sample(y,n))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1866 if(n>1)tmp_p<-sapply(1:M,function(j)sapply(1:k,function(i)sum(tmp[,j]==i)))/n
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1867 if(n==1)tmp_p<-sapply(1:M,function(j)sapply(1:k,function(i)sum(tmp[j]==i)))/n
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1868 (apply(tmp_p,1,function(z)quantile(z,c(alpha/2/(k-1),1-alpha/2/(k-1)))))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1869 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1870
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1871
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1872
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1873 p_value<-function(x=c(33,12,21),M=100000,x_obs=c(2,5,3)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1874 n=sum(x_obs)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1875 N<-sum(x)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1876 k<-length(x)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1877 y<-rep(1:k,x)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1878 tmp<-sapply(1:M,function(i)sample(y,n))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1879 if(n>1)tmp_p<-sapply(1:M,function(j)sapply(1:k,function(i)sum(tmp[,j]==i)))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1880 if(n==1)tmp_p<-sapply(1:M,function(j)sapply(1:k,function(i)sum(tmp[j]==i)))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1881 tmp<-rbind(sapply(1:3,function(i)sum(tmp_p[i,]>=x_obs[i])/M),
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1882 sapply(1:3,function(i)sum(tmp_p[i,]<=x_obs[i])/M))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1883 sapply(1:3,function(i){if(tmp[1,i]>=tmp[2,i])return(-tmp[2,i])else return(tmp[1,i])})
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1884 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1885
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1886 #"D:\\Sequences\\IMGT Germlines\\Human_SNPless_IGHJ.FASTA"
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1887 # Remove SNPs from IMGT germline segment alleles
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1888 generateUnambiguousRepertoire <- function(repertoireInFile,repertoireOutFile){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1889 repertoireIn <- read.fasta(repertoireInFile, seqtype="DNA",as.string=T,set.attributes=F,forceDNAtolower=F)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1890 alleleNames <- sapply(names(repertoireIn),function(x)strsplit(x,"|",fixed=TRUE)[[1]][2])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1891 SNPs <- tapply(repertoireIn,sapply(alleleNames,function(x)strsplit(x,"*",fixed=TRUE)[[1]][1]),function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1892 Indices<-NULL
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1893 for(i in 1:length(x)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1894 firstSeq = s2c(x[[1]])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1895 iSeq = s2c(x[[i]])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1896 Indices<-c(Indices,which(firstSeq[1:320]!=iSeq[1:320] & firstSeq[1:320]!="." & iSeq[1:320]!="." ))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1897 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1898 return(sort(unique(Indices)))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1899 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1900 repertoireOut <- repertoireIn
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1901 repertoireOut <- lapply(names(repertoireOut), function(repertoireName){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1902 alleleName <- strsplit(repertoireName,"|",fixed=TRUE)[[1]][2]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1903 geneSegmentName <- strsplit(alleleName,"*",fixed=TRUE)[[1]][1]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1904 alleleSeq <- s2c(repertoireOut[[repertoireName]])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1905 alleleSeq[as.numeric(unlist(SNPs[geneSegmentName]))] <- "N"
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1906 alleleSeq <- c2s(alleleSeq)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1907 repertoireOut[[repertoireName]] <- alleleSeq
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1908 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1909 names(repertoireOut) <- names(repertoireIn)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1910 write.fasta(repertoireOut,names(repertoireOut),file.out=repertoireOutFile)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1911
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1912 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1913
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1914
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1915
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1916
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1917
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1918
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1919 ############
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1920 groupBayes2 = function(indexes, param_resultMat){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1921
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1922 BayesGDist_Focused_CDR = calculate_bayesG( x=param_resultMat[indexes,1], N=apply(param_resultMat[indexes,c(1,2,4)],1,sum,na.rm=T), p=apply(param_resultMat[indexes,5:8],1,function(x){x[1]/(x[1]+x[2]+x[4])}))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1923 BayesGDist_Focused_FWR = calculate_bayesG( x=param_resultMat[indexes,3], N=apply(param_resultMat[indexes,c(3,2,4)],1,sum,na.rm=T), p=apply(param_resultMat[indexes,5:8],1,function(x){x[3]/(x[3]+x[2]+x[4])}))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1924 #BayesGDist_Local_CDR = calculate_bayesG( x=param_resultMat[indexes,1], N=apply(param_resultMat[indexes,c(1,2)],1,sum,na.rm=T), p=apply(param_resultMat[indexes,5:8],1,function(x){x[1]/(x[1]+x[2])}))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1925 #BayesGDist_Local_FWR = calculate_bayesG( x=param_resultMat[indexes,3], N=apply(param_resultMat[indexes,c(3,4)],1,sum,na.rm=T), p=apply(param_resultMat[indexes,5:8],1,function(x){x[3]/(x[3]+x[4])}))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1926 #BayesGDist_Global_CDR = calculate_bayesG( x=param_resultMat[indexes,1], N=apply(param_resultMat[indexes,c(1,2,3,4)],1,sum,na.rm=T), p=apply(param_resultMat[indexes,5:8],1,function(x){x[1]/(x[1]+x[2]+x[3]+x[4])}))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1927 #BayesGDist_Global_FWR = calculate_bayesG( x=param_resultMat[indexes,3], N=apply(param_resultMat[indexes,c(1,2,3,4)],1,sum,na.rm=T), p=apply(param_resultMat[indexes,5:8],1,function(x){x[3]/(x[1]+x[2]+x[3]+x[4])}))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1928 return ( list("BayesGDist_Focused_CDR"=BayesGDist_Focused_CDR,
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1929 "BayesGDist_Focused_FWR"=BayesGDist_Focused_FWR) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1930 #"BayesGDist_Local_CDR"=BayesGDist_Local_CDR,
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1931 #"BayesGDist_Local_FWR" = BayesGDist_Local_FWR))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1932 # "BayesGDist_Global_CDR" = BayesGDist_Global_CDR,
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1933 # "BayesGDist_Global_FWR" = BayesGDist_Global_FWR) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1934
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1935
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1936 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1937
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1938
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1939 calculate_bayesG <- function( x=array(), N=array(), p=array(), max_sigma=20, length_sigma=4001){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1940 G <- max(length(x),length(N),length(p))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1941 x=array(x,dim=G)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1942 N=array(N,dim=G)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1943 p=array(p,dim=G)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1944
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1945 indexOfZero = N>0 & p>0
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1946 N = N[indexOfZero]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1947 x = x[indexOfZero]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1948 p = p[indexOfZero]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1949 G <- length(x)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1950
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1951 if(G){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1952
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1953 cons<-array( dim=c(length_sigma,G) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1954 if(G==1) {
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1955 return(calculate_bayes(x=x[G],N=N[G],p=p[G],max_sigma=max_sigma,length_sigma=length_sigma))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1956 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1957 else {
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1958 for(g in 1:G) cons[,g] <- calculate_bayes(x=x[g],N=N[g],p=p[g],max_sigma=max_sigma,length_sigma=length_sigma)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1959 listMatG <- convolutionPowersOfTwoByTwos(cons,length_sigma=length_sigma)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1960 y<-calculate_bayesGHelper(listMatG,length_sigma=length_sigma)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1961 return( y/sum(y)/(2*max_sigma/(length_sigma-1)) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1962 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1963 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1964 return(NA)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1965 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1966 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1967
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1968
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1969 calculate_bayesGHelper <- function( listMatG,length_sigma=4001 ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1970 matG <- listMatG[[1]]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1971 groups <- listMatG[[2]]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1972 i = 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1973 resConv <- matG[,i]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1974 denom <- 2^groups[i]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1975 if(length(groups)>1){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1976 while( i<length(groups) ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1977 i = i + 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1978 resConv <- weighted_conv(resConv, matG[,i], w= {{2^groups[i]}/denom} ,length_sigma=length_sigma)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1979 #cat({{2^groups[i]}/denom},"\n")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1980 denom <- denom + 2^groups[i]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1981 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1982 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1983 return(resConv)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1984 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1985
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1986 weighted_conv<-function(x,y,w=1,m=100,length_sigma=4001){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1987 lx<-length(x)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1988 ly<-length(y)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1989 if({lx<m}| {{lx*w}<m}| {{ly}<m}| {{ly*w}<m}){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1990 if(w<1){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1991 y1<-approx(1:ly,y,seq(1,ly,length.out=m))$y
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1992 x1<-approx(1:lx,x,seq(1,lx,length.out=m/w))$y
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1993 lx<-length(x1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1994 ly<-length(y1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1995 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1996 else {
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1997 y1<-approx(1:ly,y,seq(1,ly,length.out=m*w))$y
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1998 x1<-approx(1:lx,x,seq(1,lx,length.out=m))$y
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
1999 lx<-length(x1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2000 ly<-length(y1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2001 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2002 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2003 else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2004 x1<-x
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2005 y1<-approx(1:ly,y,seq(1,ly,length.out=floor(lx*w)))$y
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2006 ly<-length(y1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2007 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2008 tmp<-approx(x=1:(lx+ly-1),y=convolve(x1,rev(y1),type="open"),xout=seq(1,lx+ly-1,length.out=length_sigma))$y
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2009 tmp[tmp<=0] = 0
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2010 return(tmp/sum(tmp))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2011 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2012
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2013 ########################
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2014
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2015
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2016
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2017
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2018 mutabilityMatrixONE<-rep(0,4)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2019 names(mutabilityMatrixONE)<-NUCLEOTIDES
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2020
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2021 # triMutability Background Count
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2022 buildMutabilityModelONE <- function( inputMatrixIndex, model=0 , multipleMutation=0, seqWithStops=0, stopMutations=0){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2023
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2024 #rowOrigMatInput = matInput[inputMatrixIndex,]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2025 seqGL = gsub("-", "", matInput[inputMatrixIndex,2])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2026 seqInput = gsub("-", "", matInput[inputMatrixIndex,1])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2027 matInput[inputMatrixIndex,] <<- c(seqInput,seqGL)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2028 seqLength = nchar(seqGL)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2029 mutationCount <- analyzeMutations(inputMatrixIndex, model, multipleMutation, seqWithStops)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2030 BackgroundMatrix = mutabilityMatrixONE
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2031 MutationMatrix = mutabilityMatrixONE
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2032 MutationCountMatrix = mutabilityMatrixONE
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2033 if(!is.na(mutationCount)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2034 if((stopMutations==0 & model==0) | (stopMutations==1 & (sum(mutationCount=="Stop")<length(mutationCount))) | (model==1 & (sum(mutationCount=="S")>0)) ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2035
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2036 # ONEmerStartPos = 1:(seqLength)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2037 # ONEmerLength <- length(ONEmerStartPos)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2038 ONEmerGL <- s2c(seqGL)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2039 ONEmerSeq <- s2c(seqInput)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2040
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2041 #Background
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2042 for(ONEmerIndex in 1:seqLength){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2043 ONEmer = ONEmerGL[ONEmerIndex]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2044 if(ONEmer!="N"){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2045 ONEmerCodonPos = getCodonPos(ONEmerIndex)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2046 ONEmerReadingFrameCodon = c2s(ONEmerGL[ONEmerCodonPos])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2047 ONEmerReadingFrameCodonInputSeq = c2s(ONEmerSeq[ONEmerCodonPos] )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2048
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2049 # All mutations model
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2050 #if(!any(grep("N",ONEmerReadingFrameCodon))){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2051 if(model==0){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2052 if(stopMutations==0){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2053 if(!any(grep("N",ONEmerReadingFrameCodonInputSeq)))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2054 BackgroundMatrix[ONEmer] <- (BackgroundMatrix[ONEmer] + 1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2055 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2056 if( !any(grep("N",ONEmerReadingFrameCodonInputSeq)) & translateCodonToAminoAcid(ONEmerReadingFrameCodonInputSeq)!="*"){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2057 positionWithinCodon = which(ONEmerCodonPos==ONEmerIndex)#positionsWithinCodon[(ONEmerCodonPos[1]%%3)+1]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2058 BackgroundMatrix[ONEmer] <- (BackgroundMatrix[ONEmer] + probNonStopMutations[ONEmerReadingFrameCodon,positionWithinCodon])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2059 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2060 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2061 }else{ # Only silent mutations
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2062 if( !any(grep("N",ONEmerReadingFrameCodonInputSeq)) & translateCodonToAminoAcid(ONEmerReadingFrameCodonInputSeq)!="*" & translateCodonToAminoAcid(ONEmerReadingFrameCodonInputSeq)==translateCodonToAminoAcid(ONEmerReadingFrameCodon) ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2063 positionWithinCodon = which(ONEmerCodonPos==ONEmerIndex)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2064 BackgroundMatrix[ONEmer] <- (BackgroundMatrix[ONEmer] + probSMutations[ONEmerReadingFrameCodon,positionWithinCodon])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2065 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2066 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2067 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2068 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2069 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2070
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2071 #Mutations
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2072 if(stopMutations==1) mutationCount = mutationCount[mutationCount!="Stop"]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2073 if(model==1) mutationCount = mutationCount[mutationCount=="S"]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2074 mutationPositions = as.numeric(names(mutationCount))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2075 mutationCount = mutationCount[mutationPositions>2 & mutationPositions<(seqLength-1)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2076 mutationPositions = mutationPositions[mutationPositions>2 & mutationPositions<(seqLength-1)]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2077 countMutations = 0
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2078 for(mutationPosition in mutationPositions){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2079 ONEmerIndex = mutationPosition
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2080 ONEmer = ONEmerSeq[ONEmerIndex]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2081 GLONEmer = ONEmerGL[ONEmerIndex]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2082 ONEmerCodonPos = getCodonPos(ONEmerIndex)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2083 ONEmerReadingFrameCodon = c2s(ONEmerSeq[ONEmerCodonPos])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2084 ONEmerReadingFrameCodonGL =c2s(ONEmerGL[ONEmerCodonPos])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2085 if(!any(grep("N",ONEmer)) & !any(grep("N",GLONEmer))){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2086 if(model==0){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2087 countMutations = countMutations + 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2088 MutationMatrix[GLONEmer] <- (MutationMatrix[GLONEmer] + 1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2089 MutationCountMatrix[GLONEmer] <- (MutationCountMatrix[GLONEmer] + 1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2090 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2091 if( translateCodonToAminoAcid(ONEmerReadingFrameCodonGL)!="*" ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2092 countMutations = countMutations + 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2093 positionWithinCodon = which(ONEmerCodonPos==ONEmerIndex)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2094 glNuc = substr(ONEmerReadingFrameCodonGL,positionWithinCodon,positionWithinCodon)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2095 inputNuc = substr(ONEmerReadingFrameCodon,positionWithinCodon,positionWithinCodon)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2096 MutationMatrix[GLONEmer] <- (MutationMatrix[GLONEmer] + substitution[glNuc,inputNuc])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2097 MutationCountMatrix[GLONEmer] <- (MutationCountMatrix[GLONEmer] + 1)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2098 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2099 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2100 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2101 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2102
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2103 seqMutability = MutationMatrix/BackgroundMatrix
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2104 seqMutability = seqMutability/sum(seqMutability,na.rm=TRUE)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2105 #cat(inputMatrixIndex,"\t",countMutations,"\n")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2106 return(list("seqMutability" = seqMutability,"numbMutations" = countMutations,"seqMutabilityCount" = MutationCountMatrix, "BackgroundMatrix"=BackgroundMatrix))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2107 # tmp<-list("seqMutability" = seqMutability,"numbMutations" = countMutations,"seqMutabilityCount" = MutationCountMatrix)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2108 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2109 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2110
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2111 ################
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2112 # $Id: trim.R 989 2006-10-29 15:28:26Z ggorjan $
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2113
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2114 trim <- function(s, recode.factor=TRUE, ...)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2115 UseMethod("trim", s)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2116
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2117 trim.default <- function(s, recode.factor=TRUE, ...)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2118 s
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2119
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2120 trim.character <- function(s, recode.factor=TRUE, ...)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2121 {
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2122 s <- sub(pattern="^ +", replacement="", x=s)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2123 s <- sub(pattern=" +$", replacement="", x=s)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2124 s
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2125 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2126
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2127 trim.factor <- function(s, recode.factor=TRUE, ...)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2128 {
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2129 levels(s) <- trim(levels(s))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2130 if(recode.factor) {
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2131 dots <- list(x=s, ...)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2132 if(is.null(dots$sort)) dots$sort <- sort
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2133 s <- do.call(what=reorder.factor, args=dots)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2134 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2135 s
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2136 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2137
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2138 trim.list <- function(s, recode.factor=TRUE, ...)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2139 lapply(s, trim, recode.factor=recode.factor, ...)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2140
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2141 trim.data.frame <- function(s, recode.factor=TRUE, ...)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2142 {
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2143 s[] <- trim.list(s, recode.factor=recode.factor, ...)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2144 s
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2145 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2146 #######################################
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2147 # Compute the expected for each sequence-germline pair by codon
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2148 getExpectedIndividualByCodon <- function(matInput){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2149 if( any(grep("multicore",search())) ){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2150 facGL <- factor(matInput[,2])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2151 facLevels = levels(facGL)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2152 LisGLs_MutabilityU = mclapply(1:length(facLevels), function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2153 computeMutabilities(facLevels[x])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2154 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2155 facIndex = match(facGL,facLevels)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2156
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2157 LisGLs_Mutability = mclapply(1:nrow(matInput), function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2158 cInput = rep(NA,nchar(matInput[x,1]))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2159 cInput[s2c(matInput[x,1])!="N"] = 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2160 LisGLs_MutabilityU[[facIndex[x]]] * cInput
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2161 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2162
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2163 LisGLs_Targeting = mclapply(1:dim(matInput)[1], function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2164 computeTargeting(matInput[x,2],LisGLs_Mutability[[x]])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2165 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2166
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2167 LisGLs_MutationTypes = mclapply(1:length(matInput[,2]),function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2168 #print(x)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2169 computeMutationTypes(matInput[x,2])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2170 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2171
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2172 LisGLs_R_Exp = mclapply(1:nrow(matInput), function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2173 Exp_R <- rollapply(as.zoo(1:readEnd),width=3,by=3,
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2174 function(codonNucs){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2175 RPos = which(LisGLs_MutationTypes[[x]][,codonNucs]=="R")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2176 sum( LisGLs_Targeting[[x]][,codonNucs][RPos], na.rm=T )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2177 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2178 )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2179 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2180
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2181 LisGLs_S_Exp = mclapply(1:nrow(matInput), function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2182 Exp_S <- rollapply(as.zoo(1:readEnd),width=3,by=3,
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2183 function(codonNucs){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2184 SPos = which(LisGLs_MutationTypes[[x]][,codonNucs]=="S")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2185 sum( LisGLs_Targeting[[x]][,codonNucs][SPos], na.rm=T )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2186 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2187 )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2188 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2189
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2190 Exp_R = matrix(unlist(LisGLs_R_Exp),nrow=nrow(matInput),ncol=readEnd/3,T)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2191 Exp_S = matrix(unlist(LisGLs_S_Exp),nrow=nrow(matInput),ncol=readEnd/3,T)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2192 return( list( "Expected_R"=Exp_R, "Expected_S"=Exp_S) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2193 }else{
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2194 facGL <- factor(matInput[,2])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2195 facLevels = levels(facGL)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2196 LisGLs_MutabilityU = lapply(1:length(facLevels), function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2197 computeMutabilities(facLevels[x])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2198 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2199 facIndex = match(facGL,facLevels)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2200
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2201 LisGLs_Mutability = lapply(1:nrow(matInput), function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2202 cInput = rep(NA,nchar(matInput[x,1]))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2203 cInput[s2c(matInput[x,1])!="N"] = 1
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2204 LisGLs_MutabilityU[[facIndex[x]]] * cInput
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2205 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2206
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2207 LisGLs_Targeting = lapply(1:dim(matInput)[1], function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2208 computeTargeting(matInput[x,2],LisGLs_Mutability[[x]])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2209 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2210
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2211 LisGLs_MutationTypes = lapply(1:length(matInput[,2]),function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2212 #print(x)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2213 computeMutationTypes(matInput[x,2])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2214 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2215
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2216 LisGLs_R_Exp = lapply(1:nrow(matInput), function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2217 Exp_R <- rollapply(as.zoo(1:readEnd),width=3,by=3,
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2218 function(codonNucs){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2219 RPos = which(LisGLs_MutationTypes[[x]][,codonNucs]=="R")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2220 sum( LisGLs_Targeting[[x]][,codonNucs][RPos], na.rm=T )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2221 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2222 )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2223 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2224
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2225 LisGLs_S_Exp = lapply(1:nrow(matInput), function(x){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2226 Exp_S <- rollapply(as.zoo(1:readEnd),width=3,by=3,
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2227 function(codonNucs){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2228 SPos = which(LisGLs_MutationTypes[[x]][,codonNucs]=="S")
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2229 sum( LisGLs_Targeting[[x]][,codonNucs][SPos], na.rm=T )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2230 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2231 )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2232 })
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2233
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2234 Exp_R = matrix(unlist(LisGLs_R_Exp),nrow=nrow(matInput),ncol=readEnd/3,T)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2235 Exp_S = matrix(unlist(LisGLs_S_Exp),nrow=nrow(matInput),ncol=readEnd/3,T)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2236 return( list( "Expected_R"=Exp_R, "Expected_S"=Exp_S) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2237 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2238 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2239
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2240 # getObservedMutationsByCodon <- function(listMutations){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2241 # numbSeqs <- length(listMutations)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2242 # obsMu_R <- matrix(0,nrow=numbSeqs,ncol=readEnd/3,dimnames=list(c(1:numbSeqs),c(1:(readEnd/3))))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2243 # obsMu_S <- obsMu_R
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2244 # temp <- mclapply(1:length(listMutations), function(i){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2245 # arrMutations = listMutations[[i]]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2246 # RPos = as.numeric(names(arrMutations)[arrMutations=="R"])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2247 # RPos <- sapply(RPos,getCodonNumb)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2248 # if(any(RPos)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2249 # tabR <- table(RPos)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2250 # obsMu_R[i,as.numeric(names(tabR))] <<- tabR
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2251 # }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2252 #
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2253 # SPos = as.numeric(names(arrMutations)[arrMutations=="S"])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2254 # SPos <- sapply(SPos,getCodonNumb)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2255 # if(any(SPos)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2256 # tabS <- table(SPos)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2257 # obsMu_S[i,names(tabS)] <<- tabS
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2258 # }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2259 # }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2260 # )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2261 # return( list( "Observed_R"=obsMu_R, "Observed_S"=obsMu_S) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2262 # }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2263
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2264 getObservedMutationsByCodon <- function(listMutations){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2265 numbSeqs <- length(listMutations)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2266 obsMu_R <- matrix(0,nrow=numbSeqs,ncol=readEnd/3,dimnames=list(c(1:numbSeqs),c(1:(readEnd/3))))
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2267 obsMu_S <- obsMu_R
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2268 temp <- lapply(1:length(listMutations), function(i){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2269 arrMutations = listMutations[[i]]
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2270 RPos = as.numeric(names(arrMutations)[arrMutations=="R"])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2271 RPos <- sapply(RPos,getCodonNumb)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2272 if(any(RPos)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2273 tabR <- table(RPos)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2274 obsMu_R[i,as.numeric(names(tabR))] <<- tabR
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2275 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2276
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2277 SPos = as.numeric(names(arrMutations)[arrMutations=="S"])
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2278 SPos <- sapply(SPos,getCodonNumb)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2279 if(any(SPos)){
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2280 tabS <- table(SPos)
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2281 obsMu_S[i,names(tabS)] <<- tabS
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2282 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2283 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2284 )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2285 return( list( "Observed_R"=obsMu_R, "Observed_S"=obsMu_S) )
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2286 }
ba33b94637ca Uploaded
davidvanzessen
parents: 0
diff changeset
2287