Введение в расшифровку ЭКГ

Что такое ЭКГ

Электрокардиограмма (ЭКГ) – это графическая запись электрической активности сердца, отражающая, как проходят импульсы возбуждения (деполяризация) и восстановления (реполяризация) через сердечную мышцу. ЭКГ показывает работу сердца во времени: каждое сокращение регистрируется как характерная кривая. 

ЭКГ снимается несколькими электродами, которые крепятся к различным частям тела. Положения электродов называются отведениями. Как правило, при снятии ЭКГ регистрируются 12 отведений. Есть три стандартных отведения от конечностей – I, II и III, три усиленных отведения от конечностей – aVR, aVL и aVF и шесть грудных – V1, V2, V3, V4, V5, V6. Важно отметить, что сигналы, зарегистрированные в различных отведениях, отображают работу одного и того же сердца, то есть они показывают одну и ту же электрическую активность, но с разных углов обзора. Положения электродов для стандартных отведений представлены на картинке ниже.

Запись ЭКГ 

Ниже вы видите фото ЭКГ, снятого с 12 отведений. В левом нижнем углу кардиограммы можно посмотреть степень усиления сигнала, скорость движения бумаги и используемый фильтр. Степень усиления сигнала показывает масштаб записи сигнала по вертикали, обычно используется 10 мм / мВ (миниВольт), то есть 1 мм бумаги составит 0,1 мВ. От скорости движения бумаги зависит детализированность ЭКГ и соответственно расход бумаги. Чтобы правильно посчитать продолжительности элементов кардиограммы важно учитывать скорость снятия ЭКГ, так как от нее зависит, сколько секунд будет показывать один миллиметр на бумаге (маленькая клеточка). Скорость движения бумаги может составлять 25 мм/с (1 мм – 0,04 с), 50мм/с (1 мм – 0,02 с), 100мм/с (1 мм – 0,01 с) или другие значения. Фильтры используются в кардиографах для удаления шумов и устранения дрейфа базовой линии (изолинии, изоэлектрической линии). В некоторых ситуациях применяются так же фильтры для устранения возможных помех при снятии ЭКГ, например, существует специальный антитреморный фильтр. В различных моделях кардиографов используются различные фильтры, обычно их можно включить или выключить до снятия ЭКГ.

Основные элементы ЭКГ 

Для постановки диагноза врач должен найти на записи элементы кардиограммы и проанализировать их. Ключевые элементы кардиограммы: зубцы P, T, U (может отсутствовать или  находиться на волне T) комплекс QRS (состоящий из зубцов Q, R, S – зубцы Q и S могут быть малыми, большими или отсутствовать) интервалы PQ (PR), PP, RR, QT, ТР, сегмент ST и другие. На картинке ниже представлены основные элементы кардиограммы.

Все положения зубцов кардиограммы обычно ищутся на отведении II. Изоэлектрическая линия представляет собой сегмент PQ (от конца P до начала Q) или сегмент TP (от конца T до начала P), которые в норме совпадают с изолинией, однако оценка изоэлектрической линии может быть приблизительной из-за погрешностей снятия ЭКГ. Так, автор этого урока для медиков отмечает точки на нескольких сегментах TP и проводит изолинию. От изолинии считаются все амплитуды зубцов. Кроме амплитуд и продолжительностей зубцов и интервалов, врач оценивает частоту и регулярность сердечных сокращений по расстояниям между зубцами R, определяет синусовый ритм или нет по критериям. Также на кардиограмме могут быть выявлены некоторые другие дополнительные волны.

С тем, что из себя представляют все пики и интервалы вы можете ознакомиться здесь или в других статьях для медицинских работников. Нам важно только понимать, что для каждого из элементов кардиограммы установлены некоторые диапазоны значений продолжительностей и / или амплитуд, которые считаются нормальными. Отклонения этих значений от нормы могут говорить о каких-то нарушениях в работе сердца. Например, увеличение длительности комплекса QRS может указывать на блокаду ножек пучка Гиса, а снижение его амплитуды – на инфаркт, гипертрофию, электролитные нарушения. Врачам необходимо вручную на ЭКГ находить все ее элементы, а потом проводить все расчеты, поэтому автоматическая детекция всех пиков и интервалов и проведение расчетов может сократить время работы врача на анализ ЭКГ. 

Написание программы для анализа ЭКГ

Для экспериментов был использован датасет PTB-XL, содержащий 21837 клинических ЭКГ продолжительностью 10 секунд с частотой дискретизации 100 и 500 Гц / с в 12 отведениях от 18885 пациентов. Частота дискретизации (в Герцах, количество отсчетов в секунду) – частота, с которой данные отбираются из непрерывного сигнала, влияющая на плотность информации и качество сигнала. 

Детекция пиков с помощью библиотеки neurokit2

Для нахождения зубцов и их границ были использованы методы библиотеки neurokit2 на Python, предназначенной для обработки биомедицинских сигналов (не только ЭКГ). К сожалению, исходный код библиотеки не открытый, и подробности реализации алгоритмов неизвестны. Тем не менее, в документации к библиотеке можно посмотреть, какие методы там реализованы, и прочитать по ним научные статьи.

В методе ecg_clean() реализована фильтрация сигнала, заданный алгоритм ‘neurokit’ применяет высокочастотный фильтр Баттерворта частотой 0,5 Гц (порядка 5), за которым следует фильтрация по изолинии. На вход методу подается исходный сигнал и частота дискретизации в Герцах. 

После применения фильтров сначала нужно найти положение зубцов R методом ecg_peaks(). Был использован алгоритм khamis2016 (UNSW метод). Как сообщают авторы оригинальной статьи (Khamis et al., 2016), алгоритм UNSW генерирует характерный сигнал, содержащий информацию об амплитуде и производной ЭКГ, который фильтруется в соответствии с его частотным содержанием и применяется адаптивный порог. При этом эффективность алгоритма UNSW значительно превысила показатели алгоритмов Pan-Tompkins (также реализован в библиотеке neurokit2) и Gutiérrez-Rivas (в библиотеке реализована его адаптация), а чувствительность и положительная предсказуемость превысили 95%. В рамках этой работы не был проведен анализ литературы, посвященный доступным в библиотеке алгоритмам, и не были проверены все методы на размеченных данных. Вы можете самостоятельно поэкспериментировать с доступными алгоритмами.

В методе также есть возможность корректировать артефакты (ошибки снятия ЭКГ), так как выбранный датасет содержит информацию об артефактах, мы можем проверять их наличие и принимать решение об их корректировке. Метод ecg_peaks() возвращает датафрейм signals длины исходного сигнала, заполненный 0 и 1 в тех местах, где обнаружен зубец R, а также словарь info, содержащий все индексы с зубцом R и частоту дискретизации.

С помощью метода ecg_rate() находится частота сердечных сокращений. Для интерполяции скорости между пиками была использована монотонная кубическая интерполяция (значение по умолчанию).

Далее по найденным пикам R мы можем найти остальные пики и их границы в помощью метода ecg_delineate() алгоритмом «prominence». На вход метод принимает отфильтрованный сигнал, найденные зубцы R и  частоту дискретизации. Функция ниже принимает на вход сигнал (в минивольтах, отведение II) и частоту дискретизации (500) и булевое значение требуется ли корректировать артефакты. 

Нахождение изолинии и границ комплекса QRS

Приведу примеры использования полученных результатов для нахождения другой информации из ЭКГ. Изоэлектрическая линия на ЭКГ находится по сегментам PQ (или PR) и сегментам TP, которые в норме лежат на изолинии. Так как границы комплекса QRS еще не найдены, изоэлектрическую линию можно посчитать как среднее или медиану между всеми сегментами TP. Конечно, можно взять и только один найденный сегмент, но несмотря на применение фильтров, базовая линия все еще незначительно дрейфует, поэтому было решено остановиться на медиане между всеми TP. Код ниже находит положение изолинии (по отведению II).

import plotly.graph_objects as go

def visualize_ecg(result, fs=500, save_path='ecg.json'):

    fig = go.Figure()

    time = result['timeline_x']

    x_tick = 0.02 # 1мм по x

    y_tick = 0.1 # 1мм по y

    colors = {'p_peaks': 'purple',

              'q_peaks': 'green',

              'r_peaks': 'red',

              's_peaks': 'blue',

              't_peaks': 'orange',

              'u_peaks': 'magenta',

              'j_points': 'black'}

    peaks = {

        'p_peaks': result['p_peaks'],

        'q_peaks': result['q_peaks'],

        'r_peaks': result['r_peaks'],

        's_peaks': result['s_peaks'],

        't_peaks': result['t_peaks'],

        'u_peaks': result['u_peaks'],

        'j_points': result['j_points']

    }

    baseline = result['baseline']

    for i, (lead_name, lead_result) in enumerate(result['leads_results'].items()):

        offset = -i * 2 # расстояние между отведениями

        y = np.array(lead_result['clean_signal']) + offset

        fig.add_trace(go.Scatter(

            x=time,

            y=y,

            mode='lines',

            line=dict(color='black', width=1),

            showlegend=False

        ))

        for key, color in colors.items():

            lead_peaks = [p/fs for p in peaks[key] if p is not None]

            fig.add_trace(go.Scatter(

                x=lead_peaks,

                y=[y[int(p*fs)] for p in lead_peaks if int(p*fs) < len(y)],

                mode='markers',

                marker=dict(color=color, size=6, symbol="circle"),

                name=f"{key}",

                showlegend=(i == 0)

            ))

        fig.add_hline(y=baseline + offset, line_color='gray', line_dash='dash',

                annotation_text="baseline", annotation_position="bottom left")

        fig.add_annotation(

            x=0,

            y=y[0] + 0.1,

            text=f"{lead_name}",

            showarrow=False,

            xanchor='left',

            yanchor='bottom',

            font=dict(size=10)

        )

    for key, color in colors.items():

        lead_peaks = [p/fs for p in peaks[key] if p is not None]

        for peak in lead_peaks:

            fig.add_vline(x=peak, line_color=color)

    aspect_ratio = x_tick / y_tick

    fig.update_xaxes(

        showgrid=True,

        gridcolor="red",

        gridwidth=0.4,

        dtick=x_tick*5,

        minor=dict(dtick=x_tick, showgrid=True, gridcolor="lightcoral", gridwidth=0.2),

        zeroline=False,

        showticklabels=False

    )

    fig.update_yaxes(

        showgrid=True,

        gridcolor="red",

        gridwidth=0.4,

        dtick=y_tick*5,

        minor=dict(dtick=y_tick, showgrid=True, gridcolor="lightcoral", gridwidth=0.2),

        zeroline=False,

        scaleanchor="x",

        scaleratio=aspect_ratio,

        showticklabels=False

    )

    fig.update_layout(

        height=800,

        width=1200,

        plot_bgcolor="white",

        paper_bgcolor="white",

    )

    graph_json = fig.to_json()

    with open(save_path, "w", encoding="utf-8") as f:

        f.write(graph_json)

    fig.show()

    return graph_json

Если вы захотите разработать любое приложение с использованием plotly, вы можете сохранить график в json или html и передать на фронтенд. Вызов метода show() открывает страницу в браузере с полученным графиком (был визуализирован сигнал файла 0003_hr датасета PTB-XL). 

Мы можем приближать и отдалять график, двигать его в разные стороны и наводить на конкретные точки. Так, на рисунке ниже точно видны погрешности в детекции первого R-пика (первая красная линия должна быть левее). Вероятно, и второй зубец T должен быть чуть левее. Тем не менее, методы библиотеки neurokit2 более-менее хорошо справляются с обработкой сигнала. Найденные изолиния (серый пунктир) и точка j (черная) тоже в целом выглядят корректными. 
.

Обучение модели классификации болезней сердца на обнаруженных зубцах в отведении II 

Датасет PTB-XL содержит разметку для каждого ЭКГ с поставленными диагнозами или определением ЭКГ как нормальной. Проверим, что мы сможем предсказать, используя пики, частоту сердечных сокращений и характеристики сигнала на отведении II (среднее и стандартное отклонение). Обучение модели проводилось на kaggle, поэтому пути к файлам начинаются с kaggle/input. 

Импортируем библиотеки:

import os
import wfdb
import pandas as pd
import numpy as np
from tqdm import tqdm
import neurokit2 as nk
from sklearn.preprocessing import StandardScaler, MultiLabelBinarizer
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import classification_report

Загружаем данные по сигналам ЭКГ с более высокой частотой дискретизации 500 (в датасете ЭКГ хранятся в файлах .dat и .hea).

base_path = "/kaggle/input/ptb-xl-dataset/ptb-xl-a-large-publicly-available-electrocardiography-dataset-1.0.1/records500"

record_paths = []
for root, dirs, files in os.walk(base_path):
    for file in files:
        if file.endswith(".dat"):
            record_paths.append(os.path.join(root, file).replace(".dat", "")) 

Загружаем датасет с метками ЭКГ в датафрейм meta.

db_csv = "/kaggle/input/ptb-xl-dataset/ptb-xl-a-large-publicly-available-electrocardiography-dataset-1.0.1/ptbxl_database.csv"
meta = pd.read_csv(db_csv)

meta["record_name"] = meta["filename_hr"].apply(lambda x: os.path.basename(x))
meta["record_name"] = meta["record_name"].str.replace(".dat", "")

Проходимся по всем файлам и находим пики в отведении II, считаем ЧСС и добавляем признаки об исходном сигнале. На этот раз был использован метод ecg_process(), который совмещает в себе и методы фильтрации и методы нахождения зубцов и их границ. 

fs = 500  # частота дискретизации
records = []
for path in tqdm(record_paths, desc="Обработка сигналов"):
    record = wfdb.rdrecord(path)
    signal = record.p_signal[:, 1]  # возьмём отведение II
    record_name = record.record_name

    try:
        signals, info = nk.ecg_process(signal, sampling_rate=fs)
        cleaned_signal = signals['ECG_Clean']

        q_peaks_idx = np.where(signals['ECG_Q_Peaks'] == 1)[0]
        r_peaks_idx = np.where(signals['ECG_R_Peaks'] == 1)[0]
        t_peaks_idx = np.where(signals['ECG_T_Peaks'] == 1)[0]
        s_peaks_idx = np.where(signals['ECG_S_Peaks'] == 1)[0]
        p_peaks_idx = np.where(signals['ECG_P_Peaks'] == 1)[0]

        q_peaks = np.mean(cleaned_signal[q_peaks_idx]) if len(q_peaks_idx) > 0 else 0
        r_peaks = np.mean(cleaned_signal[r_peaks_idx]) if len(r_peaks_idx) > 0 else 0
        s_peaks = np.mean(cleaned_signal[s_peaks_idx]) if len(s_peaks_idx) > 0 else 0
        p_peaks = np.mean(cleaned_signal[p_peaks_idx]) if len(p_peaks_idx) > 0 else 0
        t_peaks = np.mean(cleaned_signal[t_peaks_idx]) if len(t_peaks_idx) > 0 else 0

        mean_signal = np.mean(cleaned_signal)
        std_signal = np.std(cleaned_signal)

        heart_rate = np.mean(signals['ECG_Rate'])

        records.append({
            "record_name": record_name,
            "mean_signal": mean_signal, 
            "std_signal": std_signal, 
            "Q_peaks": q_peaks,
            "R_peaks": r_peaks,
            "T_peaks": t_peaks,
            "P_peaks": p_peaks,
            "S_peaks": s_peaks,
            "heart_rate": heart_rate
        })
    except Exception as e:
        print(f'{e}')
        continue

features_df = pd.DataFrame(records)
print("Извлечено признаков:", features_df.shape)

Без ошибок получилось обработать 21795 файлов из исходных 21837, в датафрейме содержится 9 признаков. Объединяем датафреймы features_df и meta по полю record_name (имя записи).

df = pd.merge(features_df, meta, on="record_name", how="inner")

df["scp_codes"] = df["scp_codes"].apply(eval)
df["labels"] = df["scp_codes"].apply(lambda d: list(d.keys()))

Выбираем самые частые метки (встречающиеся более 500 раз).

from collections import Counter

all_labels = [lbl for lst in df["labels"] for lbl in lst]
label_counts = Counter(all_labels)

valid_labels = {lbl for lbl, count in label_counts.items() if count >= 500}
df["labels"] = df["labels"].apply(lambda lst: [x for x in lst if x in valid_labels])
df = df[df["labels"].map(len) > 0]
print("Осталось меток:", len(valid_labels)) 

Осталось 24 метки (из исходных 71).

Применяем стандартизацию к признакам.

scaler = StandardScaler()
X = scaler.fit_transform(df[["mean_signal", "std_signal", "Q_peaks", "R_peaks", "T_peaks", "P_peaks", "S_peaks", "heart_rate"]]) 

Обучаем модель RandomForest для решения задачи мультилейбловой классификации (предсказывание нескольких меток для одного объекта). MultiLabelBinarizer сводит задачу мультилейбловой классификации к множеству задач бинарной классификации.

mlb = MultiLabelBinarizer()
y = mlb.fit_transform(df["labels"])

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

model = RandomForestClassifier(n_estimators=200, random_state=42)
model.fit(X_train, y_train)

Протестируем модель и посчитаем метрики. Посмотрим на Accuracy и Hamming Loss.

from sklearn.metrics import classification_report, roc_auc_score, hamming_loss, accuracy_score, roc_curve, auc

hamming = hamming_loss(y_test, y_pred)
subset_acc = accuracy_score(y_test, y_pred)

print(f"Hamming loss: {hamming:.3f}")
print(f"Subset accuracy: {subset_acc:.3f}")

Hamming loss: 0.061
Subset accuracy: 0.324 

Hamming loss измеряет долю неверно предсказанных меток по отношению к общему числу меток, то есть вычисляет среднее количество ошибок в предсказании меток. Accuracy показывает долю всех верных предсказаний (когда все метки были предсказаны верно). Получились очень хорошие показатели, но они мало о чем говорят с учетом дисбаланса классов. Посмотрим на ROC-кривые по всем меткам.

Как видно на изображении, значения AUC (площадь под ROC-кривой) для различных классов ожидаемо сильно отличаются друг от друга, но при этом они все более-менее адекватны (большая часть > 0.7). Посмотрим на f1. Макро-среднее f1 показывает средний f1 по всем классам без учёта частоты, из-за дисбаланса классов результаты ожидаемо не очень хорошие. Модель хорошо распознает частые диагнозы, но не очень хорошо справляется с редкими. Метрика микро-среднего f1 довольно хорошая для решения задачи мультилейбловой классификации (0,64), в целом, модель попадает в правильные метки. 

По полученным метрикам можно сделать вывод о том, что обучение модели классификации на выявленных признаках подходит для решения задачи предсказывания болезней сердца. Необходимы дальнейшие эксперименты с моделями и с признаками, их расширение и добавление сигналов других отведений. Для сравнения, для получения отличных метрик на тесте, модель HuBERT-ECG была обучена на 9 миллионах сигналов ЭКГ по II отведению (эмбеддингах самих исходных сигналов, отфильтрованных). При этом использование сигналов вместо признаков лишает модель интерпретируемости. 

Расшифруем некоторые наши хорошо предсказанные метки: 

SR – синусовый ритм – F1 = 0.84, recall = 0.98;

STACH – синусовая тахикардия – F1 = 0.84;

NORM – нормальная электрокардиограмма – F1 = 0.77;

LAFB – блокада передней ветви левой ножки пучка Гиса – F1 = 0.73;

AFIB – фибрилляция предсердий – F1 = 0.46.

Пример интерпретации результатов обучения модели. Синусовый ритм определяется по постоянству P-пиков, а также их положительности или отрицательности в некоторых отведениях. Хотя мы подавали на вход модели только признаки по II отведению, она все равно смогла хорошо распознать синусовый ритм (меток было много). По синусовой тахикардии, хотя меток было в разы меньше, модель тоже хорошо справилась, так как синусовая тахикардия по сути должна быть всегда сгруппирована с синусовым ритмом и мы подавали на вход модели ЧСС (это синусовый ритм с большой ЧСС). Фибрилляцию предсердий модель тоже способна распознать, вероятно, по нерегулярности ритма, комплекса QRS и пиков P, хотя для постановки такого диагноза тоже нужно смотреть на все отведения. 

Выводы 

• Для решения такой важной задачи как создание программы для автоматического анализа ЭКГ, необходимо глубоко погрузиться в предметную область (на данный момент вообще нет хорошо работающего решения, по крайней мере в России);
• Для проверки корректности детекции всех пиков и интервалов, все еще требуется размеченный датасет;
• Постановка диагноза (естественно, предварительного) полностью аналитическим путем по рассчитанным зубцам, интервалам и т.д. потребует серьезного анализа всех возможных сочетаний признаков, однако предварительные расчеты и описание ЭКГ текстом и графиком с отмеченными зубцами, сегментами, интервалами в любом случае ускорят работу врача;
• Постановка предварительного диагноза с помощью моделей машинного обучения (или нейросетей) на заранее определенных из сигнала ЭКГ его признаков может иметь успех и может добавить модели интерпретируемости.