Bioinformatics Mapping Quality Assessment Pipeline

Описание проекта

Данный проект демонстрирует полный цикл биоинформатического анализа: от скачивания данных секвенирования до создания воспроизводимого пайплайна оценки качества картирования ридов с использованием фреймворка Common Workflow Language (CWL).

Цель: Реализовать алгоритм оценки качества картирования прочтений на референсный геном с автоматизацией через CWL.

Организм: Escherichia coli (геном ~5 МБ)

---

Установка

Установка биоинформатических инструментов

# Обновление пакетов
sudo apt update

# Установка основных инструментов
sudo apt install -y \
    bwa \
    samtools \
    sra-toolkit \
    python3-pip \
    python3-full \
    graphviz \
    pipx

# Добавление pipx в PATH
pipx ensurepath
source ~/.bashrc

# Установка CWL фреймворка
pipx install cwltool

# Проверка установки
bwa 2>&1 | head -1
samtools --version | head -1
cwltool --version

Инструкция по установке CWL: Фреймворк CWL устанавливается через pipx install cwltool. Для визуализации графов необходим системный пакет graphviz. Альтернативно можно использовать pip install cwltool --break-system-packages для быстрой установки без виртуального окружения.

---

Структура проекта

bioinfo-cwl-project/
├── data/                          # Данные секвенирования (FASTQ)
│   ├── SRR2584863_1.fastq        # Forward reads
│   └── SRR2584863_2.fastq        # Reverse reads
├── reference/                     # Референсный геном
│   ├── ecoli_ref.fasta           # Геном E. coli
│   ├── ecoli_ref.fasta.amb       # BWA индексы
│   ├── ecoli_ref.fasta.ann
│   ├── ecoli_ref.fasta.bwt
│   ├── ecoli_ref.fasta.pac
│   └── ecoli_ref.fasta.sa
├── tools/                         # CWL инструменты (CommandLineTool)
│   ├── bwa_mem.cwl               # Картирование ридов
│   ├── samtools_view_sort.cwl    # Конвертация SAM → BAM + сортировка
│   ├── samtools_flagstat.cwl     # Статистика картирования
│   └── evaluate_quality.cwl      # Оценка качества (bash-скрипт)
├── logs_and_results/              # Результаты выполнения
│   ├── flagstat.txt              # Вывод samtools flagstat
│   ├── quality.txt               # Результат оценки качества (OK/not OK)
│   └── cwl_execution.log         # Лог работы CWL
├── evaluate_mapping.sh            # Bash-скрипт оценки качества
├── hello.cwl                      # Тестовый пайплайн "Hello World"
├── mapping_qc_workflow.cwl        # Основной CWL пайплайн
├── inputs.yml                     # Входные параметры для пайплайна
├── workflow_visualization.png     # Визуализация DAG пайплайна
└── README.md                      # Этот файл

---

Данные

Источники данных

1. Reads (SRA):

• Accession: SRR2584863
• Организм: Escherichia coli
• Тип: Paired-end Illumina sequencing
• Размер: ~150 МБ (сжатый)

Скачивание данных:

mkdir -p data && cd data
prefetch SRR2584863
fasterq-dump SRR2584863
rm SRR2584863.sra

2. Reference Genome:

• Source: NCBI Assembly GCF_000005845.2
• Strain: E. coli str. K-12 substr. MG1655
• Размер: ~4.6 МБ

Скачивание референса:

mkdir -p reference && cd reference
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.fna.gz
gunzip GCF_000005845.2_ASM584v2_genomic.fna.gz
mv GCF_000005845.2_ASM584v2_genomic.fna ecoli_ref.fasta

---

Локальный анализ (Bash)

Индексация референса

bwa index reference/ecoli_ref.fasta

Ручной запуск пайплайна

# 1. Картирование ридов
bwa mem reference/ecoli_ref.fasta \
    data/SRR2584863_1.fastq \
    data/SRR2584863_2.fastq > results/aligned.sam

# 2. Конвертация SAM → BAM и сортировка
samtools view -bS results/aligned.sam | samtools sort -o results/aligned.sorted.bam

# 3. Статистика картирования
samtools flagstat results/aligned.sorted.bam > results/flagstat.txt

Bash-скрипт оценки качества

Файл evaluate_mapping.sh реализует алгоритм оценки качества картирования. Он принимает файл flagstat.txt в качестве аргумента, вычисляет процент картированных ридов и сравнивает его с пороговым значением (80%):

#!/bin/bash
FLAGSTAT_FILE=$1
TOTAL_READS=$(head -n 1 "$FLAGSTAT_FILE" | awk '{print $1}')
MAPPED_READS=$(grep "mapped (" "$FLAGSTAT_FILE" | head -n 1 | awk '{print $1}')

if [ -z "$TOTAL_READS" ] || [ -z "$MAPPED_READS" ] || [ "$TOTAL_READS" -eq 0 ]; then
    echo "Ошибка чтения"
    exit 1
fi

PERCENT=$(awk "BEGIN {printf \"%.2f\", ($MAPPED_READS/$TOTAL_READS)*100}")
IS_OK=$(awk "BEGIN {print ($PERCENT >= 80.0)}")

if [ "$IS_OK" -eq 1 ]; then echo "OK"; else echo "not OK"; fi```

**Запуск:**
```bash
chmod +x evaluate_mapping.sh
./evaluate_mapping.sh results/flagstat.txt

---

CWL Пайплайн

Тест "Hello World"

Файл hello.cwl демонстрирует базовую структуру CWL инструмента:

#!/usr/bin/env cwl-runner
cwlVersion: v1.0
class: CommandLineTool
baseCommand: echo
inputs:
  message:
    type: string
    inputBinding:
      position: 1
outputs: []

Запуск:

cwltool hello.cwl --message "Hello CWL World!"

Основной пайплайн оценки качества

Структура пайплайна:

1. map_step (bwa_mem.cwl) — картирование ридов на референс
2. convert_step (samtools_view_sort.cwl) — конвертация SAM → BAM + сортировка
3. flagstat_step (samtools_flagstat.cwl) — сбор статистики картирования
4. evaluate_step (evaluate_quality.cwl) — оценка качества (OK/not OK)

Входные параметры (inputs.yml):

ref_genome:
  class: File
  path: reference/ecoli_ref.fasta
reads_fwd:
  class: File
  path: data/SRR2584863_1.fastq
reads_rev:
  class: File
  path: data/SRR2584863_2.fastq
eval_script:
  class: File
  path: evaluate_mapping.sh

Запуск пайплайна:

cwltool mapping_qc_workflow.cwl inputs.yml > logs_and_results/cwl_execution.log 2>&1

---

Результаты

1. Статистика картирования (samtools flagstat)

Файл: logs_and_results/flagstat.txt

Пример вывода:

133336 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
130000 + 0 mapped (97.50% : N/A)
133336 + 0 paired in sequencing
66668 + 0 read1
66668 + 0 read2
...

2. Оценка качества картирования

Файл: logs_and_results/quality.txt

Результат:

OK

Интерпретация:

• Всего ридов: 133,336
• Картировано: 130,000 (97.50%)
• Порог качества: 80%
• Статус: OK

---

Визуализация пайплайна

Генерация DAG (Directed Acyclic Graph)

# Генерация DOT-описания графа
cwltool --print-dot mapping_qc_workflow.cwl > workflow.dot

# Конвертация в PNG (требует graphviz)
dot -Tpng workflow.dot > workflow_visualization.png

Результат: Файл workflow_visualization.png

Описание способа визуализации

Использованные инструменты:

• cwltool --print-dot — встроенная функция CWL для генерации графа в формате DOT
• Graphviz (dot) — утилита для рендеринга графов из DOT-описания в графические форматы

---

Решение проблем

Ошибка: "externally-managed-environment"

Используйте pipx install cwltool или pip install --break-system-packages

Ошибка: "fail to locate the index files"

Убедитесь, что BWA индексы присутствуют в reference/ и указаны в secondaryFiles

Ошибка: "JavaScript evaluation error"

Проверьте экранирование $() в bash-скриптах внутри YAML

---

Ссылки на документацию

• CWL Specification
• cwltool Documentation
• BWA Manual
• Samtools Documentation
• NCBI SRA

---