R/Henikoff_weights.R

In bgafun: BGAfun A method to identify specifity determining residues in protein families

"Henikoff_weights" <-
function(my_amino){ 
## Henikoff weights,
## For each position,
## count number of residue types then weights are 1/R
##                     AAABBC = Three residue types....
##                              1/3*(number of sequences)
##
## Input is amino acid binary array generated by create_aln_amino function
## Only want to sum over residues 



my_profile<-create_profile(my_amino)
seq_number<-dim(my_amino)[1]
seq_length<-dim(my_amino)[2]

## Total number of residues in each column
my_residues_total<-apply(my_profile,2,sum) 
my_residues_total<-rep(my_residues_total>3,each=20)

my_residue_counts<-apply(my_profile>0,2,sum)
my_residue_counts<-rep(my_residue_counts,each=20)
weight<-((1/my_profile)*(1/my_residue_counts))*as.numeric(my_residues_total)
weight[grep("Inf",weight)]<-0
weight[grep("NaN",weight)]<-0
dim(weight)<-c(1,seq_length)
weights_array<-rep(weight,each=seq_number)
dim(weights_array)<-c(seq_number,seq_length)
seq_weights<-apply(weights_array*my_amino,1,sum)
total<-sum(seq_weights)

## Normalize the weights

seq_weights<-(seq_weights/total)     
return(seq_weights)
}