hr {BioSeqClass} | R Documentation |
Filter homolog sequences by sequence similarity.
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)
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". |
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.
hr
return a list of reduced sequences.
Hong Li
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) }