consensus               package:bio3d               R Documentation

_S_e_q_u_e_n_c_e _C_o_n_s_e_n_s_u_s _f_o_r _a_n _A_l_i_g_n_m_e_n_t

_D_e_s_c_r_i_p_t_i_o_n:

     Determines the consensus sequence for a given alignment at a given
     identity cutoff value.

_U_s_a_g_e:

       consensus(alignment, cutoff = 0.6)

_A_r_g_u_m_e_n_t_s:

alignment: an 'alignment' object created by the 'read.fasta' function
          or an alignment character matrix. 

  cutoff: a numeric value beteen 0 and 1, indicating the minimum
          sequence identity threshold for determining a consensus amino
          acid. Default is 0.6, or 60 percent residue identity. 

_V_a_l_u_e:

     A vector containing the consensus sequence, where '-' represents
     positions with no consensus (i.e. under the 'cutoff')

_A_u_t_h_o_r(_s):

     Barry Grant

_R_e_f_e_r_e_n_c_e_s:

     Grant, B.J. et al. (2006) _Bioinformatics_ *22*, 2695-2696.

_S_e_e _A_l_s_o:

     'read.fasta'

_E_x_a_m_p_l_e_s:

     #-- Read kinesin alignment
     aln <- read.fasta(system.file("examples/kinesin_xray.fa",package="bio3d"))

     # Generate consensus
     con <- consensus(aln)
     print(con$seq)

     # Generate consensus for sub-sequence
     con <- consensus(aln$ali[,330:474])
     print(con$seq)

     #-- Read HIV protease alignment
     aln <- read.fasta(system.file("examples/hivp_xray.fa",package="bio3d"))
     con <- consensus(aln$ali)

     # Plot residue frequency matrix
     ##png(filename = "freq.png", width = 1500, height = 780)
     col <- mono.colors(32)
     aa  <- rev(rownames(con$freq))

     image(x=1:ncol(con$freq),
           y=1:nrow(con$freq),
           z=as.matrix(rev(as.data.frame(t(con$freq)))),
           col=col, yaxt="n", xaxt="n",
           xlab="Alignment Position", ylab="Residue Type")

     # Add consensus along the axis
     axis(side=1, at=seq(0,length(con$seq),by=5))
     axis(side=2, at=c(1:22), labels=aa)
     axis(side=3, at=c(1:length(con$seq)), labels =con$seq)
     axis(side=4, at=c(1:22), labels=aa)
     grid(length(con$seq), length(aa))
     box()

     # Add consensus sequence
     for(i in 1:length(con$seq)) {
       text(i, which(aa==con$seq[i]),con$seq[i],col="white")
     }

     # Add lines for residue type separation
     abline(h=c(2.5,3.5, 4.5, 5.5, 3.5, 7.5, 9.5,
              12.5, 14.5, 16.5, 19.5), col="gray")

