inst/doc/FindMyFriends_intro.R

In FindMyFriends: Microbial Comparative Genomics in R

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
# Unpack files
location <- tempdir()
unzip(system.file('extdata', 'Mycoplasma.zip', package='FindMyFriends'),
      exdir=location)
genomeFiles <- list.files(location, full.names=TRUE, pattern='*.fasta')

## ---- echo=TRUE, eval=TRUE, message=FALSE, warning=FALSE----------------------
library(FindMyFriends)

mycoPan <- pangenome(genomeFiles[1:9], 
                     translated = TRUE, 
                     geneLocation = 'prodigal', 
                     lowMem = FALSE)
mycoPan

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
head(geneLocation(mycoPan))

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
library(reutils)
orgMeta <- lapply(orgNames(mycoPan), function(name) {
    uid <- esearch(name, 'assembly')
    taxuid <- elink(uid, dbTo = 'taxonomy')
    reutils::content(esummary(taxuid), as = 'parsed')
})
orgMeta <- lapply(lapply(orgMeta, lapply, unlist), as.data.frame)
orgMeta <- do.call(rbind, orgMeta)

mycoPan <- addOrgInfo(mycoPan, orgMeta)
head(orgInfo(mycoPan))

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
genes(mycoPan)
genes(mycoPan, split = 'organism')

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
# Query current defaults
head(defaults(mycoPan))

# Set a new default
defaults(mycoPan)$lowerLimit <- 0.6

## ---- echo=TRUE, eval=TRUE, results='hide', message=FALSE---------------------
mycoPan <- cdhitGrouping(mycoPan, kmerSize = 5, cdhitIter = TRUE, nrep = 3)

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
mycoPan

## ---- echo=TRUE, eval=TRUE, results='hide', message=FALSE---------------------
mycoPan <- neighborhoodSplit(mycoPan, lowerLimit = 0.8)

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
mycoPan

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
mycoPan <- kmerLink(mycoPan, lowerLimit=0.8)

genes(mycoPan, split='paralogue')[[1]]

mycoPan
collapseParalogues(mycoPan, combineInfo='largest')

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
# Remove a gene by raw index
removeGene(mycoPan, ind=60)

# Remove the first organism by index
removeGene(mycoPan, organism=1)
# or by name
name <- orgNames(mycoPan)[1]
removeGene(mycoPan, organism=name)

# Remove the second member of the first gene group
removeGene(mycoPan, group=1, ind=2)

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
# Get the pangenome matrix as an ExpressionSet object
as(mycoPan, 'ExpressionSet')
# or as a regular matrix
as(mycoPan, 'matrix')[1:6, ]

# Get all genes split into gene groups
genes(mycoPan, split='group')

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
groupStat(mycoPan)[[1]]

head(orgStat(mycoPan))

## ---- echo=TRUE, eval=TRUE, fig.height=6, fig.width=9, fig.align='center'-----
plotEvolution(mycoPan)

## ---- echo=TRUE, eval=TRUE, fig.height=9, fig.width=9, fig.align='center'-----
# Pangenome matrix similarity
plotSimilarity(mycoPan)
# Kmer similarity
plotSimilarity(mycoPan, type='kmer', kmerSize=5)
# No ordering
plotSimilarity(mycoPan, ordering='none')

## ---- echo=TRUE, eval=TRUE, fig.width=9, fig.height=9, fig.align='center'-----
plotTree(mycoPan, clust='ward.D2', dist='minkowski')
plotTree(mycoPan, type='kmer', kmerSize=5, clust='ward.D2', 
         dist='cosine', circular=TRUE, info='Species') + 
    ggplot2::scale_color_brewer(type='qual', palette=6)

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
library(igraph)
getNeighborhood(mycoPan, group=15, vicinity=5)

## ---- echo=TRUE, eval=TRUE, fig.align='center', fig.height=9, fig.width=9-----
plotNeighborhood(mycoPan, group=15, vicinity=5)

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
pcGraph(mycoPan)

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
localVar <- variableRegions(mycoPan, flankSize=6)
localVar[[1]]
plot(localVar[[1]]$graph)

## ---- echo=TRUE, eval=TRUE----------------------------------------------------
sessionInfo()