diff --git a/README.rst b/README.rst
index a49db33..f93f088 100644
--- a/README.rst
+++ b/README.rst
@@ -70,7 +70,7 @@ If you have sleep EEG data in standard formats (e.g. EDF or BrainVision), you ca
   # Apply a bandpass filter from 0.1 to 40 Hz
   raw.filter(0.1, 40)
   # Select a subset of EEG channels
-  raw.pick_channels(['C4-A1', 'C3-A2'])
+  raw.pick(['C4-A1', 'C3-A2'])
 
 How do I get started with YASA?
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
diff --git a/docs/faq.rst b/docs/faq.rst
index bca44eb..adf089e 100644
--- a/docs/faq.rst
+++ b/docs/faq.rst
@@ -31,7 +31,7 @@ If you have polysomnography data in European Data Format (.edf), you can use the
   # Apply a bandpass filter from 0.1 to 40 Hz
   raw.filter(0.1, 40)
   # Select a subset of EEG channels
-  raw.pick_channels(['C4-A1', 'C3-A2'])
+  raw.pick(['C4-A1', 'C3-A2'])
 
 .. ----------------------------- VISUALIZE -----------------------------
 .. raw:: html
diff --git a/docs/index.rst b/docs/index.rst
index 6d03a8f..8499106 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -70,7 +70,7 @@ If you have sleep EEG data in standard formats (e.g. EDF or BrainVision), you ca
   # Apply a bandpass filter from 0.1 to 40 Hz
   raw.filter(0.1, 40)
   # Select a subset of EEG channels
-  raw.pick_channels(['C4-A1', 'C3-A2'])
+  raw.pick(['C4-A1', 'C3-A2'])
 
 **********
 
diff --git a/notebooks/08_bandpower.ipynb b/notebooks/08_bandpower.ipynb
index f13eeb5..224876c 100644
--- a/notebooks/08_bandpower.ipynb
+++ b/notebooks/08_bandpower.ipynb
@@ -1881,7 +1881,7 @@
    ],
    "source": [
     "# Single channel bandpower\n",
-    "yasa.bandpower(raw.copy().pick_channels(['F3']), hypno=hypno_up, include=(2, 3), bandpass=True)"
+    "yasa.bandpower(raw.copy().pick(['F3']), hypno=hypno_up, include=(2, 3), bandpass=True)"
    ]
   },
   {
diff --git a/notebooks/16_EEG-HRV_coupling.ipynb b/notebooks/16_EEG-HRV_coupling.ipynb
index 7923f65..080dc2d 100644
--- a/notebooks/16_EEG-HRV_coupling.ipynb
+++ b/notebooks/16_EEG-HRV_coupling.ipynb
@@ -139,7 +139,7 @@
     "# Use MNE to load the EDF file\n",
     "raw = mne.io.read_raw_edf(path_edf, preload=True, verbose=False)\n",
     "# Keep only one EEG channel (C4-A1) and one ECG\n",
-    "raw.pick_channels(['C4-A1', 'EKG-R-EKG-L'], ordered=True)\n",
+    "raw.pick(['C4-A1', 'EKG-R-EKG-L'], ordered=True)\n",
     "raw"
    ]
   },
diff --git a/yasa/staging.py b/yasa/staging.py
index 9a59ffa..aeca3be 100644
--- a/yasa/staging.py
+++ b/yasa/staging.py
@@ -188,7 +188,7 @@ def __init__(self, raw, eeg_name, *, eog_name=None, emg_name=None, metadata=None
         ch_names = ch_names[keep_chan].tolist()
         ch_types = ch_types[keep_chan].tolist()
         # Keep only selected channels (creating a copy of Raw)
-        raw_pick = raw.copy().pick_channels(ch_names, ordered=True)
+        raw_pick = raw.copy().pick(ch_names)
 
         # Downsample if sf != 100
         assert sf > 80, "Sampling frequency must be at least 80 Hz."
diff --git a/yasa/tests/test_detection.py b/yasa/tests/test_detection.py
index ad70ace..333e006 100644
--- a/yasa/tests/test_detection.py
+++ b/yasa/tests/test_detection.py
@@ -41,7 +41,7 @@
 # MNE Raw
 data_mne = mne.io.read_raw_fif("notebooks/sub-02_mne_raw.fif", preload=True, verbose=0)
 data_mne.pick_types(eeg=True)
-data_mne_single = data_mne.copy().pick_channels(["F3"])
+data_mne_single = data_mne.copy().pick(["F3"])
 hypno_mne = np.loadtxt("notebooks/sub-02_hypno_30s.txt", dtype=str)
 hypno_mne = hypno_str_to_int(hypno_mne)
 hypno_mne = hypno_upsample_to_data(hypno=hypno_mne, sf_hypno=(1 / 30), data=data_mne)
diff --git a/yasa/tests/test_others.py b/yasa/tests/test_others.py
index ce74c4a..478e6bc 100644
--- a/yasa/tests/test_others.py
+++ b/yasa/tests/test_others.py
@@ -31,7 +31,7 @@
 # Using MNE
 data_mne = mne.io.read_raw_fif("notebooks/sub-02_mne_raw.fif", preload=True, verbose=0)
 data_mne.pick_types(eeg=True)
-data_mne_single = data_mne.copy().pick_channels(["F3"])
+data_mne_single = data_mne.copy().pick(["F3"])
 hypno_mne = np.loadtxt("notebooks/sub-02_hypno_30s.txt", dtype=str)
 hypno_mne = hypno_str_to_int(hypno_mne)
 hypno_mne = hypno_upsample_to_data(hypno=hypno_mne, sf_hypno=(1 / 30), data=data_mne)