diff --git a/deeptools/computeMatrixOperations.py b/deeptools/computeMatrixOperations.py
index 6b3272d4..0224f00a 100755
--- a/deeptools/computeMatrixOperations.py
+++ b/deeptools/computeMatrixOperations.py
@@ -546,7 +546,7 @@ def cbindMatrices(hm, args):
         # Add on additional NA initialized columns
         ncol = hm.matrix.matrix.shape[1]
         hm.matrix.matrix = np.hstack((hm.matrix.matrix, np.empty(hm2.matrix.matrix.shape)))
-        hm.matrix.matrix[:, ncol:] = np.NAN
+        hm.matrix.matrix[:, ncol:] = np.nan
 
         # Update the values
         for idx2, group in enumerate(hm2.parameters["group_labels"]):
diff --git a/deeptools/heatmapper.py b/deeptools/heatmapper.py
index f67afaf4..f86b85e2 100644
--- a/deeptools/heatmapper.py
+++ b/deeptools/heatmapper.py
@@ -376,7 +376,7 @@ def compute_sub_matrix_worker(self, chrom, start, end, score_file_list, paramete
 
         # create an empty matrix to store the values
         sub_matrix = np.zeros((len(regions), matrix_cols))
-        sub_matrix[:] = np.NAN
+        sub_matrix[:] = np.nan
 
         j = 0
         sub_regions = []
diff --git a/deeptools/plotFingerprint.py b/deeptools/plotFingerprint.py
index 3adce87a..a5a46880 100755
--- a/deeptools/plotFingerprint.py
+++ b/deeptools/plotFingerprint.py
@@ -203,10 +203,10 @@ def getCHANCE(args, idx, mat):
     """
     # Get the index of the reference sample
     if args.JSDsample not in args.bamfiles:
-        return [np.NAN, np.NAN, np.NAN]
+        return [np.nan, np.nan, np.nan]
     refIdx = args.bamfiles.index(args.JSDsample)
     if refIdx == idx:
-        return [np.NAN, np.NAN, np.NAN]
+        return [np.nan, np.nan, np.nan]
 
     subMatrix = np.copy(mat[:, [idx, refIdx]])
     subMatrix[np.isnan(subMatrix)] = 0
@@ -272,10 +272,10 @@ def getJSD(args, idx, mat):
 
     # Get the index of the reference sample
     if args.JSDsample not in args.bamfiles:
-        return np.NAN
+        return np.nan
     refIdx = args.bamfiles.index(args.JSDsample)
     if refIdx == idx:
-        return np.NAN
+        return np.nan
 
     # These will hold the coverage histograms
     chip = np.zeros(MAXLEN, dtype=int)
@@ -339,7 +339,7 @@ def signalAndBinDist(x):
 
     if abs(sum(PMFinput) - 1) > 0.01 or abs(sum(PMFchip) - 1) > 0.01:
         sys.stderr.write("Warning: At least one PMF integral is significantly different from 1! The JSD will not be returned")
-        return np.NAN
+        return np.nan
 
     # Compute the JSD from the PMFs
     M = (PMFinput + PMFchip) / 2.0
diff --git a/pyproject.toml b/pyproject.toml
index 84f14f48..1fe8eeba 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -21,7 +21,7 @@ authors = [
 ]
 requires-python = ">=3.8"
 dependencies = [
-    "numpy >= 1.9.0",
+    "numpy >= 2.0.0",
     "scipy >= 0.17.0",
     "matplotlib >= 3.5.0",
     "pysam >= 0.14.0",