R/msaMuscle.R

In msa: Multiple Sequence Alignment

msaMuscle <- function(inputSeqs,
                      cluster="default",
                      gapOpening="default",
                      gapExtension="default",
                      maxiters="default",
                      substitutionMatrix="default",
                      type="default",
                      order=c("aligned", "input"),
                      verbose=FALSE,
                      help=FALSE,
                      ...)
{
    ##########
    ## help ##
    ##########
    if (help)
    {
        printHelp("Muscle")
        return(invisible(NULL))
    }

    if (!checkFunctionAvailable("Muscle"))
        stop("Muscle is not available via msa!")

    params <- list(...)
    ##create a copy of the parameter list which is only used to
    ##avoid params which are not checked:
    ##after every check of a parameter, the parameter is deleted in the copy(!)
    ##if all parameter were checked, the copy finally should be empty...
    ##if not, there are additional parameters
    paramsCopy <- lapply(params, function(x) TRUE)

    ##set names according to Biostrings

    ##set default values according to muscle_userguide3.8
    ##and check if the parameters are consistnt

    ########################
    ########################
    ##  Common Parameters ##
    ########################
    ########################

    #############
    # inputSeqs #
    #############
    ##check if input of inputSeqs has a valid form, stop if not;
    ##if inputSeq is a file name set Flag TRUE else FALSE
    params[["inputSeqIsFileFlag"]] <- checkInputSeq(inputSeqs)


    ########
    # type #
    ########
    ##Sequence type

    ##validation of type
    type <- checkType(type, inputSeqs, "msaMuscle")

    ##check the profile score
    ##necessary, because default-values for gapOpening and gapExtension
    ##depending on used profile score;
    ##center, smoothScoreCeiling, minBestColScore and minSmoothScore too
    temporaryHelp <- checkProfileScore(type, params)
    params[["le"]] <- temporaryHelp[["le"]]
    params[["sp"]] <- temporaryHelp[["sp"]]
    params[["sv"]] <- temporaryHelp[["sv"]]
    params[["spn"]] <- temporaryHelp[["spn"]]

    #############
    # inputSeqs #
    #############
    ##transform the input Sequences to a string vector
    inputSeqs <- transformInputSeq(inputSeqs)

    #############
    # order     #
    #############
    order <- match.arg(order)

    if (order == "input")
    {
        if (params[["inputSeqIsFileFlag"]])
            stop("msaMuscle does not support order=\"input\" for reading\n",
                 "sequences directly from a FASTA file.")
        else if (is.null(names(inputSeqs)) ||
                 length(unique(names(inputSeqs))) != length(inputSeqs))
        {
            warning("order=\"input\" requires input sequences to be named\n",
                    "uniquely! Assigning default names 'Seq1'..'Seqn'\n",
                    "to sequences.")
            names(inputSeqs) <- paste0("Seq", 1:length(inputSeqs))
        }
    }

    ###########
    # cluster #
    ###########
	##Perform fast clustering of input sequences.
	##Use the -tree1 option to save the tree.

    ##default-value
    if (identical(cluster, "default") || is.null(cluster)) {
        cluster <- "upgma"
    } else {
        possibleValues <- c("upgma", "upgmamax", "upgmamin",
                            "upgmb", "neighborjoining")

        ##check if cluster contains only one single value
        if(length(cluster) != 1) {
            stop("The parameter cluster contains more than one value!")
        }

        ##check if cluster is of type character
        if(!is.character(cluster)) {
            stop("The parameter cluster should be a string!")
        }

        ##make cluster more robust
        cluster <- tolower(cluster)

        ##check, if input of parameter is valid
        if (!(cluster %in% possibleValues)){
            ##create a string with all possible Values named text
            text <- ""
            text <- paste(possibleValues, collapse=", ")
            stop("The parameter cluster only can have the values: \n", text)
        }
    }


    ######################
    # substitutionMatrix #
    ######################
    ## Substitution matrix in NCBI or WU-BLAST format.
    ## If you specify your own matrix, you should also specify:
    ## -gapOpening <g>, -gapExtension <e> -center 0.0
    ## Note that <g> and <e> MUST be negative.
    ##default:
    ##protein & le => vtml_la
    ##protein & sp => PAM200
    ##protein & sv => vtml_sp
    ##dna|rna => nuc_sp

    ##default-value
    if (is.null(substitutionMatrix) ||
        identical(substitutionMatrix, "default")) {
        substitutionMatrix <- NULL;
    }

    ##check if substitutionMatrix is a matrix
    if ((!is.null(substitutionMatrix) && !is.matrix(substitutionMatrix)) ||
        identical(mode(substitutionMatrix), "list"))
        stop("The parameter substitutionMatrix should be a matrix!")

    if (!is.null(substitutionMatrix))
    {
        headerNames <- c("A", "C", "D", "E", "F",
                         "G", "H", "I", "K", "L",
                         "M", "N", "P", "Q", "R",
                         "S", "T", "V", "W", "Y")

        if (type == "protein")
            reqNames <- headerNames
        else if (type == "dna")
            reqNames <- c("A", "C", "G", "T")
        else
            reqNames <- c("A", "C", "G", "U")

        rowPerm <- match(reqNames, rownames(substitutionMatrix))
        if (any(is.na(rowPerm)))
            stop("substitutionMatrix does not contain all necessary rows")

        colPerm <- match(reqNames, colnames(substitutionMatrix))
        if (any(is.na(colPerm)))
            stop("substitutionMatrix does not contain all necessary columns")

        substitutionMatrix <- substitutionMatrix[rowPerm, colPerm]

        if (type == "rna")
            reqNames <- c("A", "C", "G", "T")

        auxMat <- matrix(0, length(headerNames), length(headerNames))
        rownames(auxMat) <- headerNames
        colnames(auxMat) <- headerNames
        auxMat[reqNames, reqNames] <- substitutionMatrix
        substitutionMatrix <- auxMat

        if (!isSymmetric(substitutionMatrix))
            stop("substitutionMatrix should be a symmetric matrix!")

        if (any(is.na(substitutionMatrix)) || any(is.na(substitutionMatrix)) ||
            any(is.infinite(substitutionMatrix)))
            stop("substitutionMatrix contains invalid values!")

        params[["le"]] <- FALSE
        params[["sv"]] <- FALSE

        if (type == "protein")
        {
            params[["sp"]] <- TRUE
            params[["spn"]] <- FALSE
        }
        else
        {
            params[["sp"]] <- FALSE
            params[["spn"]] <- TRUE
        }

        paramsCopy[["le"]] <- NULL
        paramsCopy[["sv"]] <- NULL
        paramsCopy[["sp"]] <- NULL
        paramsCopy[["spn"]] <- NULL
     }

    ##############
    # gapOpening #
    ##############
    ##The gap open score. Must be negative
    ##defaultValues:
    ##le: -2.9
    ##sp: -1439
    ##sv: -300
    ##spn_dna:400
    ##spn_rna:420

    if (params$le) {
        gapOpening <- checkGapOpening2(gapOpening, substitutionMatrix, 2.9)
    }
    else if (params$sp) {
        gapOpening <- checkGapOpening2(gapOpening, substitutionMatrix, 1439)
    }
    else if (params$sv) {
        gapOpening <- checkGapOpening2(gapOpening, substitutionMatrix, 300)
    }
    else if (params$spn) {
        if (identical(type,"dna")) {
        gapOpening <- checkGapOpening2(gapOpening, substitutionMatrix, 400)
        }
        if (identical(type,"rna")) {
            gapOpening <- checkGapOpening2(gapOpening, substitutionMatrix, 420)
        }
        if (identical(type,"protein")) {
           stop("If you use sequences of type \"protein\", \n",
                "you can't use the parameter \"spn\"!")
        }
    }

    ##FIXME TODO: check default-Value for type=protein
    ################
    # gapExtension #
    ################
    ##The gap extend score. Must be negative
    ##GapExtension only used if type=dna/rna
    ##default-Value:
    ##type= "dna" => gapExtension=0
    ##type= "rna" => gapExtension=0
    ##type= "protein" 0> gapExtension=???

    gapExtension <- checkGapExtension(gapExtension,
                                      type, substitutionMatrix, 0, 0)


    ############
    # maxiters #
    ############
    ##Maximum number of iterations

    maxiters <- checkMaxiters(maxiters, 16, "msaMuscle")


    ###########
    # verbose #
    ###########
    ##Write parameter settings and progress messages to log file

    verbose <- checkVerbose(FALSE, verbose)

    ######################################
    ######################################
    ######################################
    ###  ALGORITHM SPECIFIC PARAMETERS ###
    ######################################
    ######################################
    ######################################

    ####################
    ####################
    ##  VALUE-OPTIONS ##
    ####################
    ####################

    #################
    # anchorspacing #
    #################
    ##Minimum spacing between anchor columns

    params[["anchorspacing"]] <- checkIntegerParamsNew("anchorspacing", params)

    ##delete param in copy
    paramsCopy[["anchorspacing"]] <- NULL


    ##########
    # center #
    ##########
    ##Center parameter. Should be 0 or negative
    ##default-Values:
    ##le=TRUE: -0.52
    ##spn=TRUE, type="rna": -300
    ##sp=TRUE or SV=TRUE or spn=TRUE and type="dna": 0

        params[["center"]] <- checkNumericParamsNew("center", params)
        params[["center"]] <- checkNegativeParams("center", params)


    ##delete param in copy
    paramsCopy[["center"]] <- NULL

    ############
    # cluster1 #
    ############
    ##Clustering method. cluster1 is used in iteration 1 and 2,
    ##cluster2 in later iterations.

    posVal <- c("upgma", "upgmamax", "upgmamin", "upgmb", "neighborjoining")
    params[["cluster1"]] <- checkSingleValParamsNew("cluster1", params, posVal)

    ##delete param in copy
    paramsCopy[["cluster1"]] <- NULL

    ############
    # cluster2 #
    ############
    ##Clustering method. cluster1 is used in iteration 1 and 2,
    ##cluster2 in later iterations.

    posVal <- c("upgma", "upgmb", "upgmamax", "upgmamin", "neighborjoining")
    params[["cluster2"]] <- checkSingleValParamsNew("cluster2", params, posVal)

    ##delete param in copy
    paramsCopy[["cluster2"]] <- NULL

    #############
    # diagbreak #
    #############
    ##Maximum distance between two diagonals that allows
    ##them to merge into one diagonal.

    params[["diagbreak"]] <- checkIntegerParamsNew("diagbreak", params)
    params[["diagbreak"]] <- checkPositiveParams("diagbreak", params)

    ##delete param in copy
    paramsCopy[["diagbreak"]] <- NULL

    ##############
    # diaglength #
    ##############
    ##Minimum length of diagonal

    params[["diaglength"]] <- checkIntegerParamsNew("diaglength", params)
    params[["diaglength"]] <- checkPositiveParams("diaglength", params)

    ##delete param in copy
    paramsCopy[["diaglength"]] <- NULL

    ##############
    # diagmargin #
    ##############
    ## Discard this many positions at ends of diagonal

    params[["diagmargin"]] <- checkIntegerParamsNew("diagmargin", params)
    params[["diagmargin"]] <- checkPositiveParams("diagmargin", params)

    ##delete param in copy
    paramsCopy[["diagmargin"]] <- NULL

    #############
    # distance1 #
    #############
    ##Distance measure for iteration 1.

    ##possibleValues
    if (type == "protein") {
        posVal <- c("kmer6_6", "kmer20_3", "kmer20_4", "kbit20_3")
    } else {
        posVal <- c("kmer6_6", "kmer20_3", "kmer20_4", "kbit20_3", "kmer4_6")
    }
    ##check, if input of distance1 is valid
    if (!is.null(params[["distance1"]])) {
        params[["distance1"]] <- checkIsValue("distance1", params, posVal)
    }

    ##delete param in copy
    paramsCopy[["distance1"]] <- NULL

    #############
    # distance2 #
    #############
    ##Distance measure for iterations 2, 3 ...

    params[["distance2"]] <- checkSingleValParamsNew(
                      "distance2", params, c("pctidkimura", "pctidlog"))

    ##delete param in copy
    paramsCopy[["distance2"]] <- NULL


    #########
    # hydro #
    #########
    ##Window size for determining whether a region is hydrophobic

    params[["hydro"]] <- checkIntegerParamsNew("hydro", params)

    ##delete param in copy
    paramsCopy[["hydro"]] <- NULL

    ###############
    # hydrofactor #
    ###############
    ##Multiplier for gap open/close penalties in hydrophobic regions

    params[["hydrofactor"]] <- checkNumericParamsNew("hydrofactor", params)

    ##delete param in copy
    paramsCopy[["hydrofactor"]] <- NULL

    #######
    # in1 #
    #######
    if (!is.null(params[["in1"]])) {
        checkInFile("in1", params)
    }

    ##delete param in copy
    paramsCopy[["in1"]] <- NULL

    #######
    # in2 #
    #######
    if (!is.null(params[["in2"]])) {
        checkInFile("in2", params)
    }

    ##delete param in copy
    paramsCopy[["in2"]] <- NULL

    ############
    # maxhours #
    ############
    ##Maximum time to run in hours. The actual time may exceed
    ##the requested limit by a few minutes. Decimals are allowed,
    ##so 1.5 means one hour and 30 minutes.

    params[["maxhours"]] <- checkNumericParamsNew("maxhours", params)
    params[["maxhours"]] <- checkPositiveParams("maxhours", params)

    ##delete param in copy
    paramsCopy[["maxhours"]] <- NULL

    ############
    # maxtrees #
    ############
    ##Maximum number of new trees to build in iteration 2

    params[["maxtrees"]] <- checkIntegerParamsNew("maxtrees", params)
    params[["maxtrees"]] <- checkPositiveParams("maxtrees", params)

    ##delete param in copy
    paramsCopy[["maxtrees"]] <- NULL


    ###################
    # minbestcolscore #
    ###################
    ##Minimum score a column must have to be an anchor
    ##default-Values
    ##le=TRUE=2.0
    ##sp=TRUE: 300.0
    ##sv=TRUE: 130.0
    ##spn=TRUE: 90.0

    params[["minbestcolscore"]] <- checkNumericParamsNew("minbestcolscore",
                                                         params)

    ##delete param in copy
    paramsCopy[["minbestcolscore"]] <- NULL

    ##################
    # minsmoothscore #
    ##################
    ##Minimum smoothed score a column must have to be an anchor
    ##default-Values
    ##le=TRUE=1.0
    ##sp=TRUE: 125.0
    ##sv=TRUE: 40.0
    ##spn=TRUE: 90.0

    params[["minsmoothscore"]] <- checkNumericParamsNew("minsmoothscore",
                                                        params)

    ##delete param in copy
    paramsCopy[["minsmoothscore"]] <- NULL

    ############
    # objscore #
    ############
    ##Objective score used by tree dependent refinement.
    ##sp=sum-of-pairs score.
    ##spf=sum-of-pairs score (dimer approximation)
    ##spm=sp for < 100 seqs, otherwise spf
    ##dp=dynamic programming score
    ##ps=average profile-sequence score
    ##xp=cross profile score

    posVal <- c("dp", "ps", "sp", "spf", "spm", "xp")
    params[["objscore"]] <- checkSingleValParamsNew("objscore", params, posVal)


    #delete param in copy
    paramsCopy[["objscore"]] <- NULL


    ################
    # refinewindow #
    ################
    ##Length of window for -refinew

    params[["refinewindow"]] <- checkIntegerParamsNew("refinewindow", params)
    params[["refinewindow"]] <- checkPositiveParams("refinewindow", params)

    ##delete param in copy
    paramsCopy[["refinewindow"]] <- NULL


    #########
    # root1 #
    #########
    ##Method used to root tree; root1 is used in iteration 1 and 2,
    ##root2 in later iterations

    posVal <- c("pseudo", "midlongestspan", "minavgleafdist")
    params[["root1"]] <- checkSingleValParamsNew("root1", params, posVal)

    ##delete param in copy
    paramsCopy[["root1"]] <- NULL

    #########
    # root2 #
    #########
    ##Method used to root tree; root1 is used in iteration 1 and 2,
    ##root2 in later iterations

    posVal <-  c("pseudo", "midlongestspan", "minavgleafdist")
    params[["root2"]] <- checkSingleValParamsNew("root2", params, posVal)

    ##delete param in copy
    paramsCopy[["root2"]] <- NULL

    ###################
    # smoothscoreceil #
    ###################
    ##Maximum value of column score for smoothing purposes
    ##default-Values
    ##le=TRUE=3.0
    ##sp=TRUE: 200.0
    ##sv=TRUE: 90.0
    ##spn=TRUE: 999.0

    checkNumericParamsNew("smoothscoreceil", params)
    ##delete param in copy
    paramsCopy[["smoothscoreceil"]] <- NULL


    ################
    # smoothwindow #
    ################
    ##Window used for anchor column smoothing

    params[["smoothwindow"]] <- checkIntegerParamsNew("smoothwindow", params)
    params[["smoothwindow"]] <- checkPositiveParams("smoothwindow", params)

    if (!is.null(params[["smoothwindow"]]) &&
        params[["smoothwindow"]] %% 2 == 0) {
            stop("The parameter smoothwindow must be odd!")
    }

    ##delete param in copy
    paramsCopy[["smoothwindow"]] <- NULL

    #########
    # SUEFF #
    #########
    ##Constant used in UPGMB clustering. Determines the relative
    ##fraction of average linkage (SUEFF) vs. nearest-neighbor linkage
    ##(1 - SUEFF).

    params[["SUEFF"]] <- checkNumericParamsNew("SUEFF", params)

    ##check, if input is between 0 and 1
    params[["SUEFF"]] <- checkIntervalParamsNew("SUEFF", params, 0, 1)

    ##delete param in copy
    paramsCopy[["SUEFF"]] <- NULL

    #########
    # tree1 #
    #########
    ##Save tree produced in first or second iteration to given
    ##file in Newick (Phylip-compatible) format

    ##not implemented yet

    #########
    # tree2 #
    #########
    ##Save tree produced in first or second iteration to given
    ##file in Newick (Phylip-compatible) format

    ##not implemented yet

    ###########
    # usetree #
    ###########
    ##Use given tree as guide tree.
    ##Must by in Newick (Phyip-compatible) format

    ##not implemented yet

    ###########
    # weight1 #
    ###########
    ##Sequence weighting scheme. weight1 is used in iterations 1 and 2.
    ##weight2 is used for tree-dependent refinement.
    ##none=all sequences have equal weight.
    ##henikoff=Henikoff & Henikoff weighting scheme.
    ##henikoffpb=Modified Henikoff scheme as used in PSI-BLAST.
    ##clustalw=CLUSTALW method.
    ##threeway=Gotoh three-way method.

    posVal <- c("none",
                "henikoff",
                "henikoffpb",
                "gsc",
                "clustalw",
                "threeway")

    params[["weight1"]] <- checkSingleValParamsNew("weight1", params, posVal)

    ##delete param in copy
    paramsCopy[["weight1"]] <- NULL


    ###########
    # weight2 #
    ###########
    ##Sequence weighting scheme. weight1 is used in iterations 1 and 2.
    ##weight2 is used for tree-dependent refinement.
    ##none=all sequences have equal weight.
    ##henikoff=Henikoff & Henikoff weighting scheme.
    ##henikoffpb=Modified Henikoff scheme as used in PSI-BLAST.
    ##clustalw=CLUSTALW method.
    ##threeway=Gotoh three-way method.

    posVal <- c("none",
            "henikoff",
            "henikoffpb",
            "gsc",
            "clustalw",
            "threeway")
    params[["weight2"]] <- checkSingleValParamsNew("weight2", params, posVal)

    ##delete param in copy
    paramsCopy[["weight2"]] <- NULL


    ###################
    ###################
    ##  FLAG-OPTIONS ##
    ###################
    ###################

    ###########
    # anchors #
    ###########
    ##Use anchor optimization in tree dependent refinement iterations

    if (!is.null(params[["anchors"]])) {
        params[["anchors"]] <- checkLogicalParams("anchors", params, TRUE)
    }

    ##delete param in copy
    paramsCopy[["anchors"]] <- NULL

    ###########
    # brenner #
    ###########
    ##Use Steven Brenner's method for computing the root alignment

    params[["brenner"]] <- checkLogicalParams("brenner", params, FALSE)

    ##delete param in copy
    paramsCopy[["brenner"]] <- NULL

    #######
    # clw #
    #######
    ##Write output in CLUSTALW format

    ##not activated

    ##params[["clw"]] <- checkLogicalParams("clw", params, FALSE)

    ##delete param in copy
    ##paramsCopy[["clw"]] <- NULL

    #############
    # clwstrict #
    #############
    ##Write output in CLUSTALW format with the "CLUSTAL W (1.81)"
    ##header rather than the MUSCLE version.
    ##This is useful when a post-processing step is
    ##picky about the file header

    ##set by default!!!!
    ##not activated

    ##params[["clwstrict"]] <- checkLogicalParams("clwstrict", params, FALSE)

    ##delete param in copy
    ##paramsCopy[["clwstrict"]] <- NULL

    ########
    # core #
    ########
    ##Do not catch exceptions

    if (!is.null(params[["core"]])) {
        params[["core"]] <- checkLogicalParams("core", params, TRUE)
    }

    ##delete param in copy
    paramsCopy[["core"]] <- NULL

    #########
    # diags #
    #########
    ##Use diagonal optimizations. Faster, especially for closely
    ##related sequences, but may be less accurate.

    params[["diags"]] <- checkLogicalParams("diags", params, FALSE)

    ##delete param in copy
    paramsCopy[["diags"]] <- NULL

    ##########
    # diags1 #
    ##########
    ##Use diagonal optimizations in first iteration

    params[["diags1"]] <- checkLogicalParams("diags1", params, FALSE)

    ##delete param in copy
    paramsCopy[["diags1"]] <- NULL

    ##########
    # diags2 #
    ##########
    ##Use diagonal optimizations in second iteration

    params[["diags2"]] <- checkLogicalParams("diags2", params, FALSE)

    ##delete param in copy
    paramsCopy[["diags2"]] <- NULL

    #########
    # dimer #
    #########
    ##Use dimer approximation for the SP score
    ##(faster, slightly less accurate)

    params[["dimer"]] <- checkLogicalParams("dimer", params, FALSE)

    ##delete param in copy
    paramsCopy[["dimer"]] <- NULL

    #########
    # group #
    #########
    ##Group similar sequences together in the output.
    ##This is the default. See also -stable.

    ##if (!is.null(params[["group"]])){
    ##    params[["group"]] <- checkLogicalParams("group", params, TRUE)
    ##}

    ##delete param in copy
    ##paramsCopy[["group"]] <- NULL

    ########
    # html #
    ########
    ##Write output in HTML format

    ##not implemented yet
    ##not activated

    ##params[["html"]] <- checkLogicalParams("html", params, FALSE)

    ##delete param in copy
    ##paramsCopy[["html"]] <- NULL

    ######
    # le #
    ######
    ##Use log-expectation profile score (VTML240).
    ##Alternatives are to use -sp or -sv.
    ##This is the default for amino acid sequences.

    ##param already checked within function checkProfileScore

    ##delete param in copy
    paramsCopy[["le"]] <- NULL


    #######
    # msf #
    #######
    ##Write output in MSF format (default is FASTA).
    ##Designed to be compatible with the GCG package.

    ##not implemented yet
    ##not activated

    ##params[["msf"]] <- checkLogicalParams("msf", params, FALSE)

    ##delete param in copy
    ##paramsCopy[["msf"]] <- NULL

    #############
    # noanchors #
    #############
    ##Disable anchor optimization. Default is -anchors

    params[["noanchors"]] <- checkLogicalParams("noanchors", params, FALSE)

    ##check anchors <-> noanchors
    if (!is.null(params[["anchors"]])){
        ##both positive
        if (params[["anchors"]] && params[["noanchors"]]){
            stop("The parameters anchors and noanchors \n",
                 "can't be positive at the same time!")
        }
        ##both negative
        if (!params[["anchors"]] && !params[["noanchors"]]){
            stop("The parameters anchors and noanchors \n",
                 "can't be negative at the same time!")
        }
    }
    ##delete param in copy
    paramsCopy[["noanchors"]] <- NULL

    ##########
    # nocore #
    ##########
    ##Catch exceptions and give an error message if possible

    params[["nocore"]] <- checkLogicalParams("nocore", params, FALSE)
    ##check core <-> nocore
    if (!is.null(params[["core"]])){
        ##both positive
        if (params[["core"]] && params[["nocore"]]){
            stop("The parameters core and nocore \n",
                 "can't be positive at the same time!")
        }
        ##both negative
        if (!params[["core"]] && !params[["nocore"]]){
            stop("The parameters core and nocore \n",
                 "can't be negative at the same time!")
        }
    }

    ##delete param in copy
    paramsCopy[["nocore"]] <- NULL

    ########
    # phyi #
    ########
    ##Write output in Phylip interleaved format

    ##not implemented yet
    ##not activated

    ##params[["phyi"]] <- checkLogicalParams("phyi", params, FALSE)

    ##delete param in copy
    ##paramsCopy[["phyi"]] <- NULL

    ########
    # phys #
    ########
    ##Write output in Phylip sequential format

    ##not implemented yet
    ##not activated

    ##params[["phys"]] <- checkLogicalParams("phys", params, FALSE)

    ##delete param in copy
    ##paramsCopy[["phys"]] <- NULL

    ###########
    # profile #
    ###########
    ##compute profile-profile alignment. Input alignments must be given
    ##using -in1 and -in2 options.

    params[["profile"]] <- checkLogicalParams("profile", params, FALSE)

    if (params[["profile"]]) {
        if (is.null(params[["in1"]]) || is.null(params[["in2"]])) {
            stop("The parameter profile needs the following parameters: \n",
                  "in1 and in2!")
        }
    }

    ##delete param in copy
    paramsCopy[["profile"]] <- NULL

    ##########
    # refine #
    ##########
    ##Input file is already aligned, skip first two iterations
    ##and begin tree dependent refinement.

    params[["refine"]] <- checkLogicalParams("refine", params, FALSE)

    ##delete param in copy
    paramsCopy[["refine"]] <- NULL

    ###########
    # refinew #
    ###########
    ##Refine an alignment by dividing it into non-overlapping windows
    ##and re-aligning each window.
    ##Typically used for whole-genome nucleotide alignments.

    params[["refinew"]] <- checkLogicalParams("refinew", params, FALSE)

    ##delete param in copy
    paramsCopy[["refinew"]] <- NULL

    ######
    # sp #
    ######
    ##Use sum-of-pairs protein profile score (PAM200). Default is -le

    ##param already checked within function checkProfileScore

    ##delete param in copy
    paramsCopy[["sp"]] <- NULL

    #######
    # spn #
    #######
    ##Use sum-of-pairs nucleotide profile score. This is the only option
    ##for nucleotides, and is therefore the default. The substitution scores
    ##and gap penalty scores are "borrowed" from BLASTZ.

    ##param already checked within function checkProfileScore

    ##delete param in copy
    paramsCopy[["spn"]] <- NULL

    ###########
    # spscore #
    ###########
    ##Compute alignment score of profile-profile alignment.
    ##Input alignments must be given using -in1 and -in2 options.
    ##These must be pre-aligned with gapped columns as needed,
    ##i.e. must be of the same length (have same number of columns).

    params[["spscore"]] <- checkLogicalParams("spscore", params, FALSE)

    ##delete param in copy
    paramsCopy[["spscore"]] <- NULL

    ##########
    # stable #
    ##########
    ##Preserve input order of sequences in output file.
    ##Default is to group sequences by similarity (-group).
    ##WARNING THIS OPTION WAS BUGGY AND IS NOT SUPPORTED IN v3.8

    ##The termgapshalflonger-parameter is not fully supported in version 3.8
    ##not activated
    ##params[["stable"]] <- checkLogicalParams("stable", params, FALSE)

    ##delete param in copy
    ##paramsCopy[["stable"]] <- NULL

    ######
    # sv #
    ######
    ##Use sum-of-pairs profile score (VTML240). Default is -le

    ##param already checked within function checkProfileScore

    ##delete param in copy
    paramsCopy[["sv"]] <- NULL

    #############
    # termgaps4 #
    #############
    ##Use 4-way test for treatment of terminal gaps.
    ##(Cannot be disabled in this version)

    ##The termgaps4-parameter is not fully supported in version 3.8
    ##not activated

    ##params[["termgaps4"]] <- checkLogicalParams("termgaps4", params, TRUE)

    ##delete param in copy
    ##paramsCopy[["termgaps4"]] <- NULL

    ################
    # termgapsfull #
    ################
    ##Terminal gaps penalized with full penalty

    ##The termgapsfull-parameter is not fully supported in version 3.8
    ##not activated

    ##params[["termgapsfull"]] <-
    ##checkLogicalParams("termgapsfull", params, FALSE)

    ##if (identical(params[["termgapsfull"]],TRUE)) {
    ##    warning("The parameter termgapsfull is not fully
    ##    supported in version 3.8")
    ##}

    ##delete param in copy
    ##paramsCopy[["termgapsfull"]] <- NULL

    ################
    # termgapshalf #
    ################
    ##Terminal gaps penalized with half penalty

    ##The termgapshalf-parameter is not fully supported in version 3.8
    ##not activated

    ##params[["termgapshalf"]] <-
    ## checkLogicalParams("termgapshalf", params, TRUE)

    ##delete param in copy
    ##paramsCopy[["termgapshalf"]] <- NULL

    ######################
    # termgapshalflonger #
    ######################
    ##Terminal gaps penalized with half penalty if gap relative
    ##to longer sequence, otherwise with full penalty.

    ##The termgapshalflonger-parameter is not fully supported in version 3.8
    ##not activated

    ##params[["termgapshalflonger"]] <-
    ##checkLogicalParams("termgapshalflonger", params, FALSE)

    ##if (identical(params[["termgapshalflonger"]],TRUE)) {
    ##    warning("The parameter termgapshalflonger is not fully supported
    ##    in version 3.8")
    ##}

    ##delete param in copy
    ##paramsCopy[["termgapshalflonger"]] <- NULL

    ###########
    # version #
    ###########
    ##Write version string to stdout and exit

    params[["version"]] <- checkLogicalParams("version", params, FALSE)

    ##delete param in copy
    paramsCopy[["version"]] <- NULL

    #################
    #################
    ## FINAL CHECK ##
    #################
    #################
    ##check, whether there are additional parameters
    ##see line 24

    if (length(paramsCopy) != 0){
        stop("The following parameters are not known  \n",
             "(or have been specified",
             "more often than once):\n    ",
             paste(names(paramsCopy), collapse=", ", sep=""))
    }

    inputSeqNames <- names(inputSeqs)

    names(inputSeqs) <- paste0("Seq", 1:length(inputSeqs))

    result <- .Call("RMuscle", inputSeqs, cluster, -abs(gapOpening),
                    -abs(gapExtension), maxiters, substitutionMatrix, type,
                    verbose, params, PACKAGE="msa")

    out <- convertAlnRows(result$msa, type)

    if (length(inputSeqNames) > 0)
    {
        if (order == "aligned")
        {
            perm <- match(names(out@unmasked), names(inputSeqs))
            names(out@unmasked) <- inputSeqNames[perm]
        }
        else
        {
            perm <- match(names(inputSeqs), names(out@unmasked))
            out@unmasked <- out@unmasked[perm]
            names(out@unmasked) <- inputSeqNames
        }
    }
    else
        names(out@unmasked) <- NULL

    standardParams <- list(gapOpening=gapOpening,
                           gapExtension=gapExtension,
                           maxiters=maxiters,
                           verbose=verbose)

    out@params <- c(standardParams, params)
    out@call <- deparse(sys.call())
    out
}