Monthly Archives: November 2012

Подготовка данных полногеномного секвенирования (Illumina) к сборке геномаPre-processing steps for paired-end Illumina before genome assembly

Первыми шагами сборки генома или картирования заново отсеквенированного генома к уже существующему является оценка качества сырых данных и их предварительная обработка.

В данной статье мы начнем обсуждать процеду применимую в полном объеме к данным полученным с прибора Illumina по протоколу >Casava 1.8, paired-end для одного организма (то есть разбивка файлов на отдельные виды по баркодам – так называемый demultiplexing – не требуется). В случае с другими приборами и более ранними протоколами Illumina предобработка может отличатся в деталях, но принципиальные шаги останутся теми же. Короткие отсеквенированные последовательности будем называть “риды” (от англ. reads) для краткости.

Какое потребуется аппаратное обеспечение?
На аппаратном уровне рекомендованными являются сервера с мощностями от 48 Гб оперативной (!) памяти и выше, хотя именно на этапе предварительной обработки данных количество оперативной памяти не так существенно. Для большинства программ и пайплайнов (от англ pipelines) участвующих в сборке генома или картировании ридов нужен еще более мощный сервер с 512 или 1Тб оперативной памяти. Объем жесткого дискового пространства >5 ТБ. Рекомендуемая система на сервере – основанная на Unix.

Какое потребуется програмное обеспечение?
FastQC
FastX-Toolkit
Trimmomatic
На сервере должны быть установлены Perl/Python
Кроме того, очень желательно умение писать скрипты на bash и/или знание ключвых команд командной строки UNIX

Какова общая схема действий в процессе предварительной обработки данных Illumina?

  1. Разархивировать полученные файлы *.fastq.gz. С одной полной линейки IlluminaHiSeq 2000 вам будет передано порядка 50 архивных файлов, каждый файл в разахивированном виде будет занимать примерно 1 Гб (в одном файле – 4 миллиона ридов). Так как мы обсуждаем протокол paired-end это значит, что уникальные риды будут разбиты на два файла, с названиями, например, XL_001_1.fastq и XL_001_2.fastq. Каждый рид в файле XL_001_1.fastq должен иметь пару в файле XL_001_2.fastq. Пару можно определить по заголовку рида: в нём будет отличатся только одна цифра (1 или 2) обозначающая первую или вторую половину рида:

    @I-HWUSI-EAS1826:5:70N3AAAXX_FL:8:4:16707:8219 1:N:0:CGATGT
    @I-HWUSI-EAS1826:5:70N3AAAXX_FL:8:4:16707:8219 2:N:0:CGATGT

  2. Отфильтровать (удалить или записать в отдельный файл) риды на предмет наличия в них в строк с заголовком не прошедшим контроль Illumina по показателям chastity/purity (буква Y или N в подчеркнутой части заголовка рида ниже). Исключать нужно риды с индексом Y(!)

    @I-HWUSI-EAS1826:5:70N3AAAXX_FL:8:4:16707:8219 1:N:0:CGATGT
    @I-HWUSI-EAS1826:5:70N3AAAXX_FL:8:4:16707:8219 2:Y:0:CGATGT

  3. Каждый файл визуально проанализировать в FastQC и проверить такие параметры как per base sequence quality, per base GC content, наличие адаптеров и другие значения
  4. На основании визуальной оценки в FastQC определить стратегию дальнейшей фильтрации по следующим критериям: удаление адаптера с 3’-конца (clipping), обрезка 3’ конца (trimming), оценка по среднему значения phred-score и фильтрация (filtering)
  5. Удалить адаптеры, если они есть
  6. Обрезать концы (или один из концов) рида по значеням phred-score, наличию неопределенных нуклеотидов (N) или содержанию GC
  7. Отфильтровать риды с низким phred-score (либо фильтром средним по риду, либо с помощью плавающего окна)
  8. Удалить риды меньшей длины, чем некоторый порог (например, 30 bp)
  9. Если покрытие секвенирование было достаточно большим (более 10-20 на рид), исключить все риды, в которых встречается хотя бы одна N (в некоторых случаях, используют иные пороговые значения 2N, 3N)
  10. Если на одном из предыдущих этапов использовалась процедура фильтрации без учета формата paired-end – отфильтровать все файлы на наличие полных пар и отделить так называемые риды-сироты (orphan, рид без пары) в отдельный файл
  11. Объединить все риды попарно в один файл (XL_001.fastq) через разделитель между ридами (табуляция, точка с запятой, и т.д) для каждой пары файлов (XL_001_1.fastq и XL_001_2.fastq). Разделитель нужен, чтобы потом можно было сделать обратную разбивку на отдельные файлы, если потребуется
  12. Объединить все файлы в один большой файл
  13. Отфильтровать синглтоны – риды встречающиеся всего один. Этот этап возможен, разумеется, только при приличном покрытии (более 10 повторов на рид), и важен, потому что позволяет избавиться от «шума» в данных.
  14. Если следующим шагом является сборка генома, сделать файл содержащий только уникальные риды (речь идет именно о уже попарно объединенных ридах)
  15. Разбить обратно на два файла (forward и reverse complement) по разделителю

В связи с обширностью вопроса, каждый пункт будем рассматривать в деталях в отдельных заметках

В целом, важно понимать, что процедура предварительной обработки данных полногеномного секвенирования довольно творческий процесс, в котором нужно ориентироваться на исходное качество данных, метод секвенирования, протокол использованный при секвенировании и покрытие. Кроме того, можно попробовать несколько стратегий обработки данных и выбрать ту, которая в результате даст оптимальный результат при сборке генома и картировании ридов.

Картинка для featured image взята отсюда.