hr {BioSeqClass} | R Documentation

## Homolog Reduction

### Description

Filter homolog sequences by sequence similarity.

### Usage

      hr(seq, method, identity, cdhit.path)

      cdhitHR(seq, identity=0.3, cdhit.path)
      aligndisHR(seq, identity=0.6)
      distance(seq1,seq2)

      getTrain(seqfile, posfile, aa, w, identity, balance=T)
      getNegSite(posSite, seq, aa)

### Arguments

`seq` | a list with one element for each protein/gene sequence. The elements are in two parts, one the description ("desc") and the second is a character string of the biological sequence ("seq").
`identity` | a numeric value ranged from 0 to 1. It is used as a maximum identity cutoff among input sequences.
`method` | a string for the method of homolog redunction. This must be one of the strings "cdhit" or "aligndis".
`cdhit.path` | a string for the path of cdhit program directory. eg: "/people/hongli/cd-hit". It is necessary when method="cdhit".
`seq1` | a string for the protein or gene sequence.
`seq2` | a string for the protein or gene sequence. seq1 and seq2 must have same length.
`seqfile` | a string for the name of FASTA file.
`posfile` | a string for the name of file which contains the positive site dataset. It has two columns: 1st column is the protein name; 2st column is the positive site. Protein name should be consistent with the name used in seqfile.
`aa` | a character for the interested amino acid. eg: "C".
`w` | an integer for the window size of flanking peptide sequence. Window size is 2*w+1, and the central residues are the positive sites in posfile.
`balance` | a logical value indicating whether negative sites will be random selected to have the same number with positive sites.
`posSite` | a string vector for the positive sites. It is consisted of protein description and positive site, eg: "P278168:952".

### Details

`hr` employs `cdhitHR` and `aligndisHR` to filter homolog sequences. It supported following methods:

"cdhit": Use cd-hit program to quickly filter sequences by given identity. It is designed to filter full-length protein or gene sequences. "formatdb" and "blastall" are required for running cd-hit program. (http://www.bioinformatics.org/download.php/cd-hit/cd-hit-2007-0131.tar.gz or http://www.bioinformatics.org/download.php/cd-hit/cd-hit-2007-0131-win32.tar.gz)

"aligndis": Use the number of different residues to meature the identity between two sequences. It is designed to filter aligned seuqnces with equal length.

`getTrain` extract 2*w+1 flanking peptides of positive sites and filter homolog sequences. Negative sites are non-positive sites in the same proteins.

`distance` calculate the number of positions with different residues between two sequences.

### Value

`hr` return a list of reduced sequences.

### Author(s)

Hong Li

### Examples

      distance("AABD","ACBD")
      distance("AABD","ECBD")
      if(interactive()){
        file = file.path(path.package("BioSeqClass"), "example", "acetylation_K.fasta")
        library(Biostrings)
        seq = as.character(readAAStringSet(file))
        ## Homolog reduction of whole-length sequence by cd-hit
        # need cd-hit program;
        reducSeq50 = hr(seq, method="cdhit", identity=0.5, cdhit.path="/people/hongli/cd-hit")

        file = file.path(path.package("BioSeqClass"), "example", "acetylation_K.site")
        tmp = as.matrix(read.csv(file, sep="\t",header=F))
        logical = apply(tmp,1,function(x){ l=nchar(seq[x[1]]); (l>=as.numeric(x[2])+7 & as.numeric(x[2])-7>0) })
        fragment = sub.seq(seq[tmp[logical,1]], as.numeric(tmp[logical,2])-7, as.numeric(tmp[logical,2])+7)
        ## Homolog reduction of short sequence fragment
        # It may be slow.
        reducSeq = hr(fragment, method="aligndis", identity=0.4)

        ## produce train set based on given positive sites and fasta sequences.
        file = file.path(path.package("BioSeqClass"), "example", "acetylation_K.fasta")
        posfile = file.path(path.package("BioSeqClass"), "example", "acetylation_K.site")
        ## "getTrain" integrate negative set construction and homolog reduction. It is designed for site level training data.
        # It may be very slow.
        data = getTrain(file, posfile, aa="K", w=7, identity=0.4)
      }

[Package BioSeqClass version 1.38.0 Index]