shaunpwilkinson / aphid

Analysis with Profile Hidden Markov Models
21 stars 4 forks source link

aphid crashes on sequences of length 1 #9

Open leonjessen opened 3 years ago

leonjessen commented 3 years ago

Hi @shaunpwilkinson,

When including sequences of length 1 in the pHMM training, Aphid crashes with this error:

Error in checkForRemoteErrors(val) : 
  3 nodes produced errors; first error: Invalid model for AA, residue alphabet does not correspond to
            20-letter amino acid alphabet

Perhaps it has to do with lacking dimensions for the sequence length 1:

> str(x_AAbin)
List of 25
 $ E                        : 'AAbin' raw E
 $ IA                       : 'AAbin' raw [1:2] I A
 $ KQL                      : 'AAbin' raw [1:3] K Q L
 $ ICHW                     : 'AAbin' raw [1:4] I C H W
 $ QPCRH                    : 'AAbin' raw [1:5] Q P C R ...
 $ GVHWFV                   : 'AAbin' raw [1:6] G V H W ...
 $ LNRSNHF                  : 'AAbin' raw [1:7] L N R S ...
 $ TTARHMHY                 : 'AAbin' raw [1:8] T T A R ...
 $ GGEMLKSMS                : 'AAbin' raw [1:9] G G E M ...
 $ PYYYKCVNGH               : 'AAbin' raw [1:10] P Y Y Y ...
 $ KQMNMMSLLKG              : 'AAbin' raw [1:11] K Q M N ...
 $ IHVQMAQWQRLP             : 'AAbin' raw [1:12] I H V Q ...
 $ SWTGEYIGIAELI            : 'AAbin' raw [1:13] S W T G ...
 $ AWWNTWAWWTRKVY           : 'AAbin' raw [1:14] A W W N ...
 $ HEKPAMPYFWRNYMA          : 'AAbin' raw [1:15] H E K P ...
 $ HKMEEPFTNTNVPSMV         : 'AAbin' raw [1:16] H K M E ...
 $ GNKRWGCQLWCQEWHGK        : 'AAbin' raw [1:17] G N K R ...
 $ RWKPWQDELKMCYQKTDH       : 'AAbin' raw [1:18] R W K P ...
 $ GFVDRYQIANAAVINIDTQ      : 'AAbin' raw [1:19] G F V D ...
 $ VHIRTPRNYHSCFSDMYHHF     : 'AAbin' raw [1:20] V H I R ...
 $ ECTVFWPKEWTTPSHYDPCCL    : 'AAbin' raw [1:21] E C T V ...
 $ DCCRVNMFYPCGWEMGVDKKCI   : 'AAbin' raw [1:22] D C C R ...
 $ IIMWWEMGPNFLAISLMWKLWWN  : 'AAbin' raw [1:23] I I M W ...
 $ HPQYPLKSCPNGIHLRILNDGGNG : 'AAbin' raw [1:24] H P Q Y ...
 $ FTSSSSWTSNPWCCILTKWTWCNVA: 'AAbin' raw [1:25] F T S S ...
 - attr(*, "class")= chr "AAbin"

Reproducible example:

# Load libraries ----------------------------------------------------------
library("aphid")
library("ape")

# Define example data -----------------------------------------------------
set.seed(24643)
amino_acids = unlist(strsplit("ARNDCQEGHILKMFPSTWYV" , ""))
x = sapply(seq(1, 25),
           function(i){
             paste0(
               sample(amino_acids, size = i, replace = TRUE),
               collapse = "") })
names(x) = x

# Train pHMM --------------------------------------------------------------

# Set parameters
n_cores = 3
quiet = TRUE

# Binary encode sequences
x_AAbin = ape::as.AAbin(x = strsplit(x = x, split = ""))

# Derive initial pHMM
x_phmm = aphid::derivePHMM.AAbin(x = x_AAbin,
                                 seqweights = NULL,
                                 cores = n_cores,
                                 residues = "AMINO",
                                 maxiter = 1e3,
                                 deltaLL = 0.01,
                                 quiet = quiet)

# Train pHMM
x_phmm_trained = aphid::train.PHMM(x = x_phmm,
                                   y = x_AAbin,
                                   method = "BaumWelch",
                                   seqweights = NULL,
                                   cores = n_cores,
                                   maxiter = 1e3,
                                   deltaLL = 0.01,
                                   quiet = quiet)

If I change

sapply(seq(1, 25)...

to

sapply(seq(2, 25)...

The error disappears