diff --git a/docs/src/graph_creation.md b/docs/src/graph_creation.md
index e6cdaa5..76ea06c 100644
--- a/docs/src/graph_creation.md
+++ b/docs/src/graph_creation.md
@@ -6,7 +6,7 @@ are intrinsically compatible with _nauty_, performing isomorphism checks or grap
 ### Creating `NautyGraphs`
 `NautyGraphs` and `NautyDiGraphs` can be created in the same way as graphs from `Graphs.jl`.  As an example, here are three different ways to define the same graph:
 
-```jldoctest default
+```jldoctest default; output=false
 using NautyGraphs, Graphs
 
 A = [0 1 0 0;
@@ -25,7 +25,7 @@ end
 
 g3 = NautyGraph(edges)
 
-g1 == g2 == g3
+g1 == g2 == g3 # true
 
 # output
 true
@@ -40,13 +40,29 @@ If labels are not explicitly provided, they are set to zero. Here is an example
 ```jldoctest default
 julia> g4 = NautyGraph(edges; vertex_labels=[4, 3, 2, 1])
 {4, 4} undirected NautyGraph
+```
 
+After graph creation, vertex labels can be accessed and modified using
+```jldoctest default
 julia> labels(g4)
 4-element Vector{Int64}:
  4
  3
  2
  1
+
+julia> label(g4, 2) # returns the second vertex label
+4
+
+julia> setlabel!(g4, 1, 20) # sets the first label == 20
+20
+
+julia> setlabels!(g4, [1, 2, 3, 4])
+4-element Vector{Int64}:
+ 1
+ 2
+ 3
+ 4
 ```
 
 ### Adding or removing vertices and edges
diff --git a/src/NautyGraphs.jl b/src/NautyGraphs.jl
index f27a804..f3de474 100644
--- a/src/NautyGraphs.jl
+++ b/src/NautyGraphs.jl
@@ -30,7 +30,10 @@ export
     NautyDiGraph,
     DenseNautyGraph,
     AutomorphismGroup,
-    labels,
+    labels, 
+    label, 
+    setlabels!,
+    setlabel!,
     iscanon,
     nauty,
     canonize!,
diff --git a/src/densenautygraph.jl b/src/densenautygraph.jl
index 32a10e4..f6ff4d9 100644
--- a/src/densenautygraph.jl
+++ b/src/densenautygraph.jl
@@ -7,7 +7,7 @@ unsigned integer type that holds the individual bits of the graph's adjacency ma
 """
 mutable struct DenseNautyGraph{D,W<:Unsigned} <: AbstractNautyGraph{Int}
     graphset::Graphset{W}
-    labels::Vector{Int}
+    _labels::Vector{Int}
     ne::Int
     iscanon::Bool
 end
@@ -38,6 +38,17 @@ function DenseNautyGraph{D,W}(n::Integer; vertex_labels=nothing) where {D,W<:Uns
 end
 DenseNautyGraph{D}(n::Integer; vertex_labels=nothing) where {D} = DenseNautyGraph{D,UInt}(n; vertex_labels)
 
+"""
+    DenseNautyGraph{D}(; vertex_labels) where {D}
+
+Construct a vertex-labeled `DenseNautyGraph` on `length(vertex_labels)` vertices and 0 edges.
+Can be directed (`D = true`) or undirected (`D = false`).
+"""
+function DenseNautyGraph{D,W}(; vertex_labels) where {D,W<:Unsigned}
+    return DenseNautyGraph{D}(length(vertex_labels); vertex_labels)
+end
+DenseNautyGraph{D}(; vertex_labels) where {D} = DenseNautyGraph{D,UInt}(; vertex_labels)
+
 """
     DenseNautyGraph{D}(A::AbstractMatrix; [vertex_labels]) where {D}
 
@@ -63,7 +74,7 @@ function (::Type{G})(g::AbstractGraph) where {G<:AbstractNautyGraph}
 end
 function (::Type{G})(g::AbstractNautyGraph) where {G<:AbstractNautyGraph}
     h = invoke(G, Tuple{AbstractGraph}, g)
-    @views h.labels .= g.labels
+    @views h._labels .= g._labels
     h.iscanon = g.iscanon
     return h
 end
@@ -91,10 +102,10 @@ end
 DenseNautyGraph{D}(edge_list::Vector{<:AbstractEdge}; vertex_labels=nothing) where {D} = DenseNautyGraph{D,UInt}(edge_list; vertex_labels)
 
 
-Base.copy(g::G) where {G<:DenseNautyGraph} = G(copy(g.graphset), copy(g.labels), g.ne, g.iscanon)
+Base.copy(g::G) where {G<:DenseNautyGraph} = G(copy(g.graphset), copy(g._labels), g.ne, g.iscanon)
 function Base.copy!(dest::G, src::G) where {G<:DenseNautyGraph}
     copy!(dest.graphset, src.graphset)
-    copy!(dest.labels, src.labels)
+    copy!(dest._labels, src._labels)
     dest.ne = src.ne
     dest.iscanon = src.iscanon
     return dest
@@ -103,12 +114,51 @@ end
 Base.show(io::Core.IO, g::DenseNautyGraph{false}) = print(io, "{$(nv(g)), $(ne(g))} undirected NautyGraph")
 Base.show(io::Core.IO, g::DenseNautyGraph{true}) = print(io, "{$(nv(g)), $(ne(g))} directed NautyGraph")
 
-Base.hash(g::DenseNautyGraph, h::UInt) = hash(g.labels, hash(g.graphset, h))
-Base.:(==)(g::DenseNautyGraph, h::DenseNautyGraph) = (g.graphset == h.graphset) && (g.labels == h.labels)
+Base.hash(g::DenseNautyGraph, h::UInt) = hash(g._labels, hash(g.graphset, h))
+Base.:(==)(g::DenseNautyGraph, h::DenseNautyGraph) = (g.graphset == h.graphset) && (g._labels == h._labels)
 
 # BASIC GRAPH API
-labels(g::AbstractNautyGraph) = g.labels
-iscanon(g::AbstractNautyGraph) = g.iscanon
+
+"""
+    labels(g::AbstractNautyGraph)
+
+Return the vertex labels of `g`. 
+
+Do not modify the vector of labels returned. Use [`setlabels!`](@ref) or [`setlabel!`](@ref) instead.
+"""
+@inline labels(g::AbstractNautyGraph) = g._labels
+
+"""
+    label(g::AbstractNautyGraph, i::Integer)
+
+Return the label of vertex `i` of `g`.
+"""
+@inline label(g::AbstractNautyGraph, index::Integer) = g._labels[index]
+
+"""
+    setlabels!(g::AbstractNautyGraph, vertex_labels)
+
+Set the vertex labels of `g` equal to `vertex_labels`.
+"""
+@inline setlabels!(g::AbstractNautyGraph, vertex_labels) = (g.iscanon = false; g._labels .= vertex_labels)
+
+"""
+    setlabel!(g::AbstractNautyGraph, i::Integer, vertex_label)
+
+Set the label of vertex `i` of `g` equal to `vertex_label`.
+"""
+@inline setlabel!(g::AbstractNautyGraph, index::Integer, vertex_label) = (g.iscanon = false; g._labels[index] = vertex_label)
+
+"""
+    iscanon(g::AbstractNautyGraph)
+
+Return true if `g` has previously been canonized.
+
+`iscanon(g) == false` does not necessarily imply that `g` is not in canonical form, it just means `g` has never
+been explicitly canonized. This function should be considered internal and may be removed in future versions.
+"""
+@inline iscanon(g::AbstractNautyGraph) = g.iscanon
+
 Graphs.nv(g::DenseNautyGraph) = g.graphset.n
 Graphs.ne(g::DenseNautyGraph) = g.ne
 Graphs.vertices(g::DenseNautyGraph) = Base.OneTo(nv(g))
@@ -208,7 +258,8 @@ Base.zero(::Type{G}) where {G<:AbstractNautyGraph} = G(0)
 
 function _induced_subgraph(g::DenseNautyGraph, iter)
     h, vmap = invoke(Graphs.induced_subgraph, Tuple{AbstractGraph,typeof(iter)}, g, iter)
-    @views h.labels .= g.labels[vmap]
+    @views h._labels .= g._labels[vmap]
+    h.iscanon = false
     return h, vmap
 end
 Graphs.induced_subgraph(g::DenseNautyGraph, iter::AbstractVector{<:Integer}) = _induced_subgraph(g::DenseNautyGraph, iter)
@@ -246,18 +297,19 @@ function Graphs.add_vertices!(g::DenseNautyGraph, n::Integer; vertex_labels=0)
     vertex_labels isa Number || n != length(vertex_labels) && throw(ArgumentError("Incompatible length: trying to add `n=$n` vertices, but`vertex_labels` has length $(length(vertex_labels))."))
     ng = nv(g)
     _add_vertices!(g.graphset, n)
-    resize!(g.labels, ng + n)
-    g.labels[ng+1:end] .= vertex_labels
+    resize!(g._labels, ng + n)
+    g._labels[ng+1:end] .= vertex_labels
     g.iscanon = false
     return n
 end
+Graphs.add_vertices!(g::DenseNautyGraph; vertex_labels) = Graphs.add_vertices!(g, length(vertex_labels); vertex_labels)
 Graphs.add_vertex!(g::DenseNautyGraph; vertex_label::Integer=0) = Graphs.add_vertices!(g, 1; vertex_labels=vertex_label) > 0
 
 function Graphs.rem_vertices!(g::DenseNautyGraph, inds)
     all(i->has_vertex(g, i), inds) || return false
 
     _rem_vertices!(g.graphset, inds)
-    deleteat!(g.labels, inds)
+    deleteat!(g._labels, inds)
 
     g.ne = is_directed(g) ? sum(g.graphset) : (sum(g.graphset) + tr(g.graphset)) ÷ 2
     g.iscanon = false
@@ -272,5 +324,5 @@ function Graphs.blockdiag(g::DenseNautyGraph{D1,W}, h::DenseNautyGraph{D2}) wher
     gset[1:ng, 1:ng] .= g.graphset
     gset[ng+1:end, ng+1:end] .= h.graphset
     D = D1 || D2
-    return DenseNautyGraph{D,W}(gset; vertex_labels=vcat(g.labels, h.labels))
+    return DenseNautyGraph{D,W}(gset; vertex_labels=vcat(labels(g), labels(h)))
 end
\ No newline at end of file
diff --git a/src/nauty.jl b/src/nauty.jl
index 502d3b0..ff0c5d8 100644
--- a/src/nauty.jl
+++ b/src/nauty.jl
@@ -80,7 +80,7 @@ function _densenauty(g::DenseNautyGraph{D,W}, options::NautyOptions=default_opti
     # TODO: allow the user to pass pre-allocated arrays for lab, ptn, orbits, canong in a safe way.
     n, m = g.graphset.n, g.graphset.m
 
-    lab, ptn = vertexlabels2labptn(g.labels)
+    lab, ptn = vertexlabels2labptn(labels(g))
     orbits = zeros(Cint, n)
     canong = Graphset{W}(n, m)
 
@@ -123,7 +123,10 @@ function nauty(g::DenseNautyGraph, options::NautyOptions=default_options(g); can
     # generators = Vector{Cint}[] # TODO: extract generators from nauty call
     autg = AutomorphismGroup(statistics.grpsize1 * 10^statistics.grpsize2, orbits)
 
-    canonize && _copycanon!(g, canong, canonperm)
+    if canonize
+        _copycanon!(g, canong, canonperm)
+        g.iscanon = true
+    end
     return canonperm, autg
 end
 
@@ -135,14 +138,14 @@ Reorder `g`'s vertices into canonical order and return the permutation used.
 function canonize!(::AbstractNautyGraph) end
 
 function canonize!(g::DenseNautyGraph)
-    iscanon(g) && return collect(Cint(1):Cint(nv(g))) # to be type stable, this needs to be Cints
+    iscanon(g) && return canonical_permutation(g)
     canong, canonperm, _ = _densenauty(g)
     _copycanon!(g, canong, canonperm)
     return canonperm
 end
 function _copycanon!(g, canong, canonperm)
     copy!(g.graphset, canong)
-    permute!(g.labels, canonperm)
+    permute!(g._labels, canonperm)
     g.iscanon = true
     return
 end
@@ -171,7 +174,7 @@ function is_isomorphic(g::DenseNautyGraph, h::DenseNautyGraph)
     iscanon(g) && iscanon(h) && return g == h
     canong, permg, _ = _densenauty(g)
     canonh, permh, _ = _densenauty(h)
-    return canong == canonh && view(g.labels, permg) == view(h.labels, permh)
+    return canong == canonh && view(g._labels, permg) == view(h._labels, permh)
 end
 ≃(g::AbstractNautyGraph, h::AbstractNautyGraph) = is_isomorphic(g, h)
 
@@ -187,10 +190,10 @@ function canonical_id end
 
 function canonical_id(g::DenseNautyGraph)
     if iscanon(g)
-        return _SHAhash(g.graphset, g.labels)
+        return _SHAhash(g.graphset, g._labels)
     else
         canong, canonperm, _ = _densenauty(g)
-        return _SHAhash(canong, @view g.labels[canonperm])
+        return _SHAhash(canong, @view g._labels[canonperm])
     end
 end
 
diff --git a/test/densenautygraph.jl b/test/densenautygraph.jl
index 52ce37c..532cfe2 100644
--- a/test/densenautygraph.jl
+++ b/test/densenautygraph.jl
@@ -197,10 +197,10 @@
 
     g4 = copy(g)
     add_vertices!(g4, 5; vertex_labels=1:5)
-    @test g4.labels == vcat(g.labels, 1:5)
+    @test labels(g4) == vcat(labels(g), 1:5)
 
     g5 = copy(g)
-    g5.labels = [1, 4, 5, 10]
+    setlabels!(g5, [1, 4, 5, 10])
     g6 = copy(g5)
     @test labels(g6) == [1, 4, 5, 10]
 
@@ -213,8 +213,8 @@
     @test h.graphset.n == g.graphset.n
     @test h.ne == g.ne
     @test h.graphset.m == g.graphset.m
-    @test h.labels == g.labels
-    @test h.iscanon == g.iscanon
+    @test labels(h) == labels(g)
+    @test iscanon(h) == iscanon(g)
 
     glab = NautyGraph(5; vertex_labels=1:5)
     add_edge!(glab, 1, 2)
@@ -224,10 +224,10 @@
     add_edge!(glab, 2, 5)
 
     gind1 = glab[[1, 5, 2]]
-    @test gind1.labels == [1, 5, 2]
+    @test labels(gind1) == [1, 5, 2]
 
     gind2 = glab[[Edge(1, 2), Edge(1, 4)]]
-    @test gind2.labels == [1, 2, 4]
+    @test labels(gind2) == [1, 2, 4]
 
     gb = NautyGraph(g0)
     vg = DiGraph(g)
@@ -235,4 +235,53 @@
     bb_ng = blockdiag(gb, g)
     bb_g = NautyDiGraph(blockdiag(DiGraph(g0), vg))
     @test bb_ng == bb_g
+
+    gl1 = NautyGraph(3; vertex_labels=[1,2,3])
+    gl2 = NautyGraph(3)
+    setlabels!(gl2, [1,2,3])
+    gl3 = NautyGraph(3)
+    foreach(1:3) do i
+        setlabel!(gl3, i, i)
+    end
+
+    @test labels(gl1) == labels(gl2) == labels(gl3)
+    @test label(gl1, 1) == label(gl2, 1)
+    @test label(gl1, 2) == label(gl2, 2)
+
+    gl4 = copy(gl1)
+    add_edge!(gl4, 1, 2)
+    gl4_id1 = canonical_id(gl4)
+    setlabels!(gl4, [3, 4, 5])
+    gl4_id2 = canonical_id(gl4)
+
+    @test gl4_id1 != gl4_id2
+
+    canonize!(gl4)
+    @test NautyGraphs.iscanon(gl4)
+
+    setlabels!(gl4, [3, 4, 5])
+    @test !NautyGraphs.iscanon(gl4)
+
+    canonize!(gl4)
+    @test NautyGraphs.iscanon(gl4)
+
+    setlabel!(gl4, 3, 3)
+    @test !NautyGraphs.iscanon(gl4)
+
+    gl5 = copy(gl1)
+    add_edge!(gl5, 1, 2)
+    gl5_id1 = canonical_id(gl5)
+    setlabel!(gl5, 3, 6)
+    gl5_id2 = canonical_id(gl5)
+    
+    @test gl5_id1 != gl5_id2
+
+    gls1 = NautyGraph(; vertex_labels=[1, 2, 3, 4])
+    gls2 = NautyGraph(4; vertex_labels=[1, 2, 3, 4])
+    @test gls1 == gls2
+
+    add_vertices!(gls1, 2; vertex_labels=[5, 6])
+    add_vertices!(gls2; vertex_labels=[5, 6])
+
+    @test gls1 == gls2
 end
\ No newline at end of file
diff --git a/test/nauty.jl b/test/nauty.jl
index 698f671..2c5eebc 100644
--- a/test/nauty.jl
+++ b/test/nauty.jl
@@ -118,7 +118,7 @@
     canonize!(canon5)
 
     canonperm5 = canonical_permutation(g5)
-    @test canon5.labels == g5.labels[canonperm5]
+    @test labels(canon5) == labels(g5)[canonperm5]
 
     # Just test that multithreading doesnt lead to errors
     thread_gs = [copy(g4) for i in 1:10]