# NOT RUN {
path <- tempdir()
dir.create (path)
# load blast results with subject accession numbers
data(blast.fill)
#load metadata of all Chlamydia pneumoniae sequences - they are subjects that
# do not count as nonspecific and may be aligned
data(meta.all)
# load metadata with target Chlamydia pneumoniae sequences - they are specific subjects
# that must be aligned
# make new accession numbers to count all WGS sequences as one (see unite_NCBI_ac.nums ())
meta.target.new.ids <- unite_NCBI_ac.nums (data = meta.target,
ac.num.var = meta.target$GB_AcNum,
title.var = meta.target$title,
db.var = meta.target$source_db,
type = "shotgun", order = TRUE,
new.titles = TRUE)
# summarize blast results, count aligned specific subjects with "switch ids" option
# (WGS sequences are counted as one). Add query cover information.
blast.sum.sp <- summarize_blast_result (sum.aligned = "sp",
blast.probe.id.var = blast.fill$Qid,
blast.res.id.var = blast.fill$Racc,
blast.res.title.var = blast.fill$Rtitle,
reference.id.var = meta.target.new.ids$new.id,
reference.title.var = meta.target.new.ids$new.title,
titles = TRUE,
add.blast.info = TRUE,
data.blast.info = data.frame(Qcover = blast.fill$Qcover),
switch.ids = TRUE, switch.table = meta.target.new.ids,
temp.db = paste0 (path, "/temp.db"), delete.temp.db = TRUE,
return = "summary", write.alignment = "DB",
alignment.db = paste0 (path, "/alig.db"))
# summarize nonspecific alignments (that are not in meta.all dataframe)
blast.sum.nonsp <- summarize_blast_result (sum.aligned = "nonsp",
blast.probe.id.var = blast.fill$Qid,
blast.res.id.var = blast.fill$Racc,
blast.res.title.var = blast.fill$Rtitle,
reference.id.var = meta.all$GB_AcNum,
reference.title.var = meta.all$title,
titles = TRUE, switch.ids = FALSE,
add.blast.info = TRUE,
data.blast.info = data.frame(Qcover = blast.fill$Qcover),
temp.db = paste0 (path, "/temp.db"),
delete.temp.db = TRUE,
return = "summary", write.alignment = "DB",
alignment.db = paste0 (path, "/alig.db"))
# all specific targets are aligned
sp.aligned <- read_from_DB(database = paste0 (path, "/alig.db"), table = "sp_aligned")
# no targets that are not aligned
sp.not.aligned <- read_from_DB(database = paste0 (path, "/alig.db"), table = "sp_not_aligned")
# No nonspecific alignments
nonsp <- read_from_DB(database = paste0 (path, "/alig.db"), table = "nonsp")
file.remove (paste0 (path, "/alig.db"))
# }
Run the code above in your browser using DataLab