inst/doc/protr.R

In protr: Generating Various Numerical Representation Schemes for Protein Sequences

## ----include=FALSE------------------------------------------------------------
knitr::opts_chunk$set(
  comment = "#>",
  collapse = TRUE
)

## -----------------------------------------------------------------------------
library("protr")

# Load FASTA files
# Note that `system.file()` is for accessing example files
# in the protr package. Replace it with your own file path.
extracell <- readFASTA(
  system.file(
    "protseq/extracell.fasta",
    package = "protr"
  )
)
mitonchon <- readFASTA(
  system.file(
    "protseq/mitochondrion.fasta",
    package = "protr"
  )
)

## ----eval = FALSE-------------------------------------------------------------
#  length(extracell)

## ----eval = FALSE-------------------------------------------------------------
#  ## [1] 325

## ----eval = FALSE-------------------------------------------------------------
#  length(mitonchon)

## ----eval = FALSE-------------------------------------------------------------
#  ## [1] 306

## ----eval = FALSE-------------------------------------------------------------
#  extracell <- extracell[(sapply(extracell, protcheck))]
#  mitonchon <- mitonchon[(sapply(mitonchon, protcheck))]

## ----eval = FALSE-------------------------------------------------------------
#  length(extracell)

## ----eval = FALSE-------------------------------------------------------------
#  ## [1] 323

## ----eval = FALSE-------------------------------------------------------------
#  length(mitonchon)

## ----eval = FALSE-------------------------------------------------------------
#  ## [1] 304

## ----eval = FALSE-------------------------------------------------------------
#  # Calculate APseAAC descriptors
#  x1 <- t(sapply(extracell, extractAPAAC))
#  x2 <- t(sapply(mitonchon, extractAPAAC))
#  x <- rbind(x1, x2)
#  
#  # Make class labels
#  labels <- as.factor(c(rep(0, length(extracell)), rep(1, length(mitonchon))))

## ----eval = FALSE-------------------------------------------------------------
#  set.seed(1001)
#  
#  # Split training and test set
#  tr.idx <- c(
#    sample(1:nrow(x1), round(nrow(x1) * 0.75)),
#    sample(nrow(x1) + 1:nrow(x2), round(nrow(x2) * 0.75))
#  )
#  te.idx <- setdiff(1:nrow(x), tr.idx)
#  
#  x.tr <- x[tr.idx, ]
#  x.te <- x[te.idx, ]
#  y.tr <- labels[tr.idx]
#  y.te <- labels[te.idx]

## ----eval = FALSE-------------------------------------------------------------
#  library("randomForest")
#  rf.fit <- randomForest(x.tr, y.tr, cv.fold = 5)
#  rf.fit

## ----eval = FALSE-------------------------------------------------------------
#  ## Call:
#  ##  randomForest(x = x.tr, y = y.tr, cv.fold = 5)
#  ##                Type of random forest: classification
#  ##                      Number of trees: 500
#  ## No. of variables tried at each split: 8
#  ##
#  ##         OOB estimate of  error rate: 25.11%
#  ## Confusion matrix:
#  ##     0   1 class.error
#  ## 0 196  46   0.1900826
#  ## 1  72 156   0.3157895

## ----eval = FALSE-------------------------------------------------------------
#  # Predict on test set
#  rf.pred <- predict(rf.fit, newdata = x.te, type = "prob")[, 1]
#  
#  # Plot ROC curve
#  library("pROC")
#  plot.roc(y.te, rf.pred, grid = TRUE, print.auc = TRUE)

## ----eval = FALSE-------------------------------------------------------------
#  ## Call:
#  ## plot.roc.default(x = y.te, predictor = rf.pred, col = "#0080ff",
#  ##                  grid = TRUE, print.auc = TRUE)
#  ##
#  ## Data: rf.pred in 81 controls (y.te 0) > 76 cases (y.te 1).
#  ## Area under the curve: 0.8697

## -----------------------------------------------------------------------------
knitr::include_graphics("figures/roc.png")

## ----extractAAC---------------------------------------------------------------
library("protr")

x <- readFASTA(system.file(
  "protseq/P00750.fasta",
  package = "protr"
))[[1]]

extractAAC(x)

## ----extractDC----------------------------------------------------------------
dc <- extractDC(x)
head(dc, n = 30L)

## ----extractTC----------------------------------------------------------------
tc <- extractTC(x)
head(tc, n = 36L)

## -----------------------------------------------------------------------------
knitr::include_graphics("figures/AAindex.png")

## ----extractMoreau1-----------------------------------------------------------
moreau <- extractMoreauBroto(x)
head(moreau, n = 36L)

## ----extractMoreau2-----------------------------------------------------------
# Define 3 custom properties
myprops <- data.frame(
  AccNo = c("MyProp1", "MyProp2", "MyProp3"),
  A = c(0.62, -0.5, 15), R = c(-2.53, 3, 101),
  N = c(-0.78, 0.2, 58), D = c(-0.9, 3, 59),
  C = c(0.29, -1, 47), E = c(-0.74, 3, 73),
  Q = c(-0.85, 0.2, 72), G = c(0.48, 0, 1),
  H = c(-0.4, -0.5, 82), I = c(1.38, -1.8, 57),
  L = c(1.06, -1.8, 57), K = c(-1.5, 3, 73),
  M = c(0.64, -1.3, 75), F = c(1.19, -2.5, 91),
  P = c(0.12, 0, 42), S = c(-0.18, 0.3, 31),
  T = c(-0.05, -0.4, 45), W = c(0.81, -3.4, 130),
  Y = c(0.26, -2.3, 107), V = c(1.08, -1.5, 43)
)

# Use 4 properties in the AAindex database and 3 customized properties
moreau2 <- extractMoreauBroto(
  x,
  customprops = myprops,
  props = c(
    "CIDH920105", "BHAR880101",
    "CHAM820101", "CHAM820102",
    "MyProp1", "MyProp2", "MyProp3"
  )
)

head(moreau2, n = 36L)

## ----extractMoran-------------------------------------------------------------
# Use the 3 custom properties defined before
# and 4 properties in the AAindex database
moran <- extractMoran(
  x,
  customprops = myprops,
  props = c(
    "CIDH920105", "BHAR880101",
    "CHAM820101", "CHAM820102",
    "MyProp1", "MyProp2", "MyProp3"
  )
)

head(moran, n = 36L)

## ----extractGeary-------------------------------------------------------------
# Use the 3 custom properties defined before
# and 4 properties in the AAindex database
geary <- extractGeary(
  x,
  customprops = myprops,
  props = c(
    "CIDH920105", "BHAR880101",
    "CHAM820101", "CHAM820102",
    "MyProp1", "MyProp2", "MyProp3"
  )
)

head(geary, n = 36L)

## -----------------------------------------------------------------------------
knitr::include_graphics("figures/CTD.png")

## ----extractCTDC--------------------------------------------------------------
extractCTDC(x)

## ----extractCTDT--------------------------------------------------------------
extractCTDT(x)

## ----extractCTDD--------------------------------------------------------------
extractCTDD(x)

## -----------------------------------------------------------------------------
knitr::include_graphics("figures/ctriad.png")

## ----extractCTriad------------------------------------------------------------
ctriad <- extractCTriad(x)
head(ctriad, n = 65L)

## ----extractSOCN--------------------------------------------------------------
extractSOCN(x)

## -----------------------------------------------------------------------------
knitr::include_graphics("figures/QSO.png")

## ----extractQSO---------------------------------------------------------------
extractQSO(x)

## -----------------------------------------------------------------------------
knitr::include_graphics("figures/PAAC.png")

## ----extractPAAC--------------------------------------------------------------
extractPAAC(x)

## -----------------------------------------------------------------------------
knitr::include_graphics("figures/APAAC.png")

## ----extractAPAAC-------------------------------------------------------------
extractAPAAC(x)

## ----extractDescScales--------------------------------------------------------
x <- readFASTA(system.file(
  "protseq/P00750.fasta",
  package = "protr"
))[[1]]

descscales <- extractDescScales(
  x,
  propmat = "AATopo",
  index = c(37:41, 43:47),
  pc = 5, lag = 7, silent = FALSE
)

## ----extractDescScales2-------------------------------------------------------
length(descscales)
head(descscales, 15)

## ----extractBLOSUM------------------------------------------------------------
x <- readFASTA(system.file(
  "protseq/P00750.fasta",
  package = "protr"
))[[1]]

blosum <- extractBLOSUM(
  x,
  submat = "AABLOSUM62",
  k = 5, lag = 7, scale = TRUE, silent = FALSE
)

## ----extractBLOSUM2-----------------------------------------------------------
length(blosum)
head(blosum, 15)

## ----eval = FALSE-------------------------------------------------------------
#  s1 <- readFASTA(system.file("protseq/P00750.fasta", package = "protr"))[[1]]
#  s2 <- readFASTA(system.file("protseq/P08218.fasta", package = "protr"))[[1]]
#  s3 <- readFASTA(system.file("protseq/P10323.fasta", package = "protr"))[[1]]
#  s4 <- readFASTA(system.file("protseq/P20160.fasta", package = "protr"))[[1]]
#  s5 <- readFASTA(system.file("protseq/Q9NZP8.fasta", package = "protr"))[[1]]
#  plist <- list(s1, s2, s3, s4, s5)
#  psimmat <- parSeqSim(plist, cores = 4, type = "local", submat = "BLOSUM62")
#  psimmat

## ----eval = FALSE-------------------------------------------------------------
#  ##            [,1]       [,2]       [,3]       [,4]       [,5]
#  ## [1,] 1.00000000 0.11825938 0.10236985 0.04921696 0.03943488
#  ## [2,] 0.11825938 1.00000000 0.18858241 0.12124217 0.06391103
#  ## [3,] 0.10236985 0.18858241 1.00000000 0.05819984 0.06175942
#  ## [4,] 0.04921696 0.12124217 0.05819984 1.00000000 0.05714638
#  ## [5,] 0.03943488 0.06391103 0.06175942 0.05714638 1.00000000

## ----eval = FALSE-------------------------------------------------------------
#  # By GO Terms
#  go1 <- c(
#    "GO:0005215", "GO:0005488", "GO:0005515",
#    "GO:0005625", "GO:0005802", "GO:0005905"
#  ) # AP4B1
#  go2 <- c(
#    "GO:0005515", "GO:0005634", "GO:0005681",
#    "GO:0008380", "GO:0031202"
#  ) # BCAS2
#  go3 <- c(
#    "GO:0003735", "GO:0005622", "GO:0005840",
#    "GO:0006412"
#  ) # PDE4DIP
#  golist <- list(go1, go2, go3)
#  parGOSim(golist, type = "go", ont = "CC", measure = "Wang")

## ----eval = FALSE-------------------------------------------------------------
#  ##       [,1]  [,2] [,3]
#  ## [1,] 1.000 0.344 0.17
#  ## [2,] 0.344 1.000 0.24
#  ## [3,] 0.170 0.240 1.00

## ----eval = FALSE-------------------------------------------------------------
#  # By Entrez gene id
#  genelist <- list(c("150", "151", "152", "1814", "1815", "1816"))
#  parGOSim(genelist, type = "gene", ont = "BP", measure = "Wang")

## ----eval = FALSE-------------------------------------------------------------
#  ##        150   151   152  1814  1815  1816
#  ## 150  1.000 0.702 0.725 0.496 0.570 0.455
#  ## 151  0.702 1.000 0.980 0.456 0.507 0.551
#  ## 152  0.725 0.980 1.000 0.460 0.511 0.538
#  ## 1814 0.496 0.456 0.460 1.000 0.687 0.473
#  ## 1815 0.570 0.507 0.511 0.687 1.000 0.505
#  ## 1816 0.455 0.551 0.538 0.473 0.505 1.000

## ----eval = FALSE-------------------------------------------------------------
#  ids <- c("P00750", "P00751", "P00752")
#  prots <- getUniProt(ids)
#  prots

## ----eval = FALSE-------------------------------------------------------------
#  ## [[1]]
#  ## [1] "MDAMKRGLCCVLLLCGAVFVSPSQEIHARFRRGARSYQVICRDEKTQMIYQQHQSWLRPVLRSNRVEYCWCN
#  ## SGRAQCHSVPVKSCSEPRCFNGGTCQQALYFSDFVCQCPEGFAGKCCEIDTRATCYEDQGISYRGTWSTAESGAECT
#  ## NWNSSALAQKPYSGRRPDAIRLGLGNHNYCRNPDRDSKPWCYVFKAGKYSSEFCSTPACSEGNSDCYFGNGSAYRGT
#  ## HSLTESGASCLPWNSMILIGKVYTAQNPSAQALGLGKHNYCRNPDGDAKPWCHVLKNRRLTWEYCDVPSCSTCGLRQ
#  ## YSQPQFRIKGGLFADIASHPWQAAIFAKHRRSPGERFLCGGILISSCWILSAAHCFQERFPPHHLTVILGRTYRVVP
#  ## GEEEQKFEVEKYIVHKEFDDDTYDNDIALLQLKSDSSRCAQESSVVRTVCLPPADLQLPDWTECELSGYGKHEALSP
#  ## FYSERLKEAHVRLYPSSRCTSQHLLNRTVTDNMLCAGDTRSGGPQANLHDACQGDSGGPLVCLNDGRMTLVGIISWG
#  ## LGCGQKDVPGVYTKVTNYLDWIRDNMRP"
#  ##
#  ## [[2]]
#  ## [1] "MGSNLSPQLCLMPFILGLLSGGVTTTPWSLARPQGSCSLEGVEIKGGSFRLLQEGQALEYVCPSGFYPYPVQ
#  ## TRTCRSTGSWSTLKTQDQKTVRKAECRAIHCPRPHDFENGEYWPRSPYYNVSDEISFHCYDGYTLRGSANRTCQVNG
#  ## RWSGQTAICDNGAGYCSNPGIPIGTRKVGSQYRLEDSVTYHCSRGLTLRGSQRRTCQEGGSWSGTEPSCQDSFMYDT
#  ## PQEVAEAFLSSLTETIEGVDAEDGHGPGEQQKRKIVLDPSGSMNIYLVLDGSDSIGASNFTGAKKCLVNLIEKVASY
#  ## GVKPRYGLVTYATYPKIWVKVSEADSSNADWVTKQLNEINYEDHKLKSGTNTKKALQAVYSMMSWPDDVPPEGWNRT
#  ## RHVIILMTDGLHNMGGDPITVIDEIRDLLYIGKDRKNPREDYLDVYVFGVGPLVNQVNINALASKKDNEQHVFKVKD
#  ## MENLEDVFYQMIDESQSLSLCGMVWEHRKGTDYHKQPWQAKISVIRPSKGHESCMGAVVSEYFVLTAAHCFTVDDKE
#  ## HSIKVSVGGEKRDLEIEVVLFHPNYNINGKKEAGIPEFYDYDVALIKLKNKLKYGQTIRPICLPCTEGTTRALRLPP
#  ## TTTCQQQKEELLPAQDIKALFVSEEEKKLTRKEVYIKNGDKKGSCERDAQYAPGYDKVKDISEVVTPRFLCTGGVSP
#  ## YADPNTCRGDSGGPLIVHKRSRFIQVGVISWGVVDVCKNQKRQKQVPAHARDFHINLFQVLPWLKEKLQDEDLGFL"
#  ##
#  ## [[3]]
#  ## [1] "APPIQSRIIGGRECEKNSHPWQVAIYHYSSFQCGGVLVNPKWVLTAAHCKNDNYEVWLGRHNLFENENTAQF
#  ## FGVTADFPHPGFNLSLLKXHTKADGKDYSHDLMLLRLQSPAKITDAVKVLELPTQEPELGSTCEASGWGSIEPGPDB
#  ## FEFPDEIQCVQLTLLQNTFCABAHPBKVTESMLCAGYLPGGKDTCMGDSGGPLICNGMWQGITSWGHTPCGSANKPS
#  ## IYTKLIFYLDWINDTITENP"

## ----protcheck----------------------------------------------------------------
x <- readFASTA(system.file("protseq/P00750.fasta", package = "protr"))[[1]]
# A real sequence
protcheck(x)
# An artificial sequence
protcheck(paste(x, "Z", sep = ""))

## ----protseg------------------------------------------------------------------
protseg(x, aa = "M", k = 5)

## -----------------------------------------------------------------------------
knitr::include_graphics("figures/protrweb.png")