WikiSort.ru - Не сортированное

SAMtools — это набор утилит для обработки коротких фрагментов секвенированной ДНК в форматах SAM или BAM. Автор SAMtools — китайский биоинформатик Хенг Ли, который является также автором спецификаций форматов SAM и BAM. В настоящее время ведущими разработчиками SAMtools являются Petr Danecek^[4] и John Marshall^[5].

Форматы SAM, BAM и CRAM

В текстовом файле формата SAM представлена информация о выравнивании фрагментов ДНК (которые также называются прочтениями или ридами) на некоторую последовательность (которая также называется референсной). Чаще всего SAM получают в результате картирования прочтений из файла FASTQ (англ.)русск. на последовательность референсного генома (англ.)русск.. Одно выравнивание состоит из нескольких строк, каждая из которых — выравнивание одного сегмента. Строка включает 11 обязательных полей, содержащих такую информацию, как позиция и качество выравнивания, направление прочтения, указание на парное прочтение и др. Кроме того, возможно указание ряда необязательных полей в виде тэг: тип: значение.

SAM-файл также может содержать заголовок в начале файла. В нём приводится общая информация о файле, в том числе использованная для выравнивания программа и указание на то, отсортирован ли файл. Строки заголовка начинаются со знака «@».

Формат BAM представляет собой бинарный эквивалент SAM. BAM занимает меньше места и позволяет быстрее работать с информацией, чем SAM. Однако только файлы SAM доступны для чтения как текстовые файлы. Samtools позволяет эффективно работать с форматом BAM и извлекать необходимую информацию в человекочитаемом формате.

Формат CRAM являются ещё более эффективными с точки зрения занимаемого дискового пространства, чем файлы BAM. В CRAM-файле хранятся отличия прочтений от референсной последовательности (англ.)русск., поэтому для работы с ним необходимо наличие файла с референсным геномом. Спецификация^[6] формата разработана в Европейском институте биоинформатики (англ.)русск.. Samtools позволяет выполнять конвертацию между форматами SAM, BAM и CRAM.

Спецификацию^[7] формата SAM можно найти в репозитории SAMtools^[8]. Ниже рассмотрено выравнивание парных ридов в формате SAM и описание полей.

r001 99 ref 7 30 8M2I4M1D3M = 37 39 TTAGATAAAGGATACTG * r001 147 ref 37 30 9M = 7 -39 CAGCGGCAT * NM:i:1

Опциональные поля должны соответствовать формату двухбуквенный тэг: тип: значение. Например, NH:i:1 указывает число выравниваний в файле для данного прочтения как целочисленную величину, равную единице. Некоторые другие распространённые теги:

• AS — вес (score) выравнивания, рассчитанный программой-картировщиком;
• NM — редакционное расстояние от прочтения до референса;
• MD — строка с информацией о невыровненных позициях; например, 10A5^AC6’ означает 10 совпадений с референсом → A в референсе, отличное от нуклеотида в соответствующей позиции прочтения → 5 совпадений → делеция (отсутствие в прочтении) двух нуклеотидов — AC → 6 совпадений;
• CC — название референса для «следующего» выравнивания («хита») — для случая неуникального выравнивания;
• CP — координата крайней левой позиции для «следующего» выравнивания («хита»);
• HI — индекс выравнивания («хита») для данного прочтения.

Опциональные поля, теги которых начинаются с X, Y или Z, зарезервированы для использования различными программами и непосредственно пользователями.

Использование

SAMtools предназначен для работы с потоком данных. Каждая программа вызывается отдельной командой, принимает входной файл через стандартный поток ввода (stdin) и возвращает результат через стандартный поток вывода (stdout). Предупреждения и сообщения об ошибках выводятся в стандартный поток ошибок (stderr). Команды samtools могут быть скобинированы в конвейеры с другими Unix-командами.^[9]

По умолчанию поток вывода направляется на экран. Так как он может быть громоздким и сложным, используется перенаправление вывода в файл (> и >>) или следующей команде в конвейере (|).

SAMtools также может открывать BAM-файлы через FTP или HTTP.

Некоторые команды поддерживают многопоточность.

Для многих команд важно явное указание референсных последовательностей в заголовке файла (поле @SQ). Кроме того, большинство команд samtools работают только с сортированными файлами и выдают ошибку при вызове с несортированным файлом.

SAMtools написан на C и может быть использован через API. Существуют обёртки для других языков программирования:

• pysam для Python^[10],
• Bio-samtools для Ruby^[11],
• Bio-SamTools для Perl^[12],
• samtools для Haskell^[13].

Стоит отметить, что существуют независимые программы для работы с форматами SAM и BAM, написанные на других языках:

• BamTools для C++^[14],
• Picard для Java^[15],
• cl-sam для Common Lisp^[16].

Команды SAMtools

Вызов команд осуществляется в виде «samtools название_команды». Далее указываются опции вызова и необходимые файлы (если файл не был передан через конвейер).

view

Команда view позволяет проводить различные операции с форматами SAM, BAM и CRAM. Например, можно конвертировать файл из одного формата в другой, фильтровать по региону, качеству или тэгу, просматривать содержимое через stdout.

Конвертация из BAM в SAM формат представляет собой простую команду с направлением в файл (Внимание! здесь информация заголовка bam файла не будет включена в sam файл):

samtools view sample.bam > sample.sam

Конвертация из SAM в BAM требует явного указания входного (-S) и выходного (-b) форматов:

samtools view -bS sample.sam > sample.bam

Часто в заголовке SAM-файла отсутствует информация о референсных последовательностях (поле @SQ). В этом случае нужно явно указать файл, который использовался для картирования (-T):

samtools view -bT ref.fa sample.sam > sample.bam

Конвертация из BAM в CRAM с производится с обязательным указанием локальной референсной последовательности:

samtools view -C -T ref.fa sample.bam > sample.cram

Команда также позволяет извлекать картирования, перекрывающихся с указанным регионом. Для этого входной файл должен быть сортирован и содержать поле @SQ в заголовке. Пример извлечения картирований, перекрывающихся с основаниями 10-13 на хромосоме 1:

samtools view sample_sorted.bam chr1:10-13

Чтобы записать извлеченные картирования в файл, требуется указать формат выходного файла (опция -b для BAM), флаг -h для включения заголовка в результирующий файл, а также формат файла:

samtools view -h -b sample_sorted.bam chr1:10-13 > tiny_sorted.bam

faidx

faidx индексирует входной FASTA-файл и создает индекс-файл fasta.fai. Команду можно применять как для полной последовательности, так и для заданного региона.

Например, если в заголовке SAM-файла не указано, по какой последовательности проводилось выравнивание (поле @SQ), можно дополнить файл этой информацией:

samtools faidx ref.fa

samtools view -bt ref.fa.fai sample.sam > sample.bam

Интересно, что при вызове команды view с опцией -T и указанием неиндексированного fastq-файла, соответствующий индекс будет сгенерирован автоматически.

tview

tview вызывает интерактивный консольный текстовый просмотрщик. Он позволяет просматривать выравнивание картирований на небольшой фрагмент референсного генома. С помощью клавиши ‘g’ можно перемещаться на разные позиции выравнивания. ‘?’ вызывает помощник с перечислением возможных горячих клавиш.

Использование:

samtools tview sample_sorted.bam ref.fasta

reheader

Команда reheader предназначена для переименования заголовка входного BAM на заголовок другого входного SAM-файла. Такая команда работает быстрее, чем конвертация BAM→ SAM → BAM.

cat

cat позволяет быстро объединять BAM-файлы, порядок записей во входных файлах сохраняется.

sort

sort сортирует прочтения во входном файле. По умолчанию сортировка производится по левой координате прочтений (то есть по первой координате прочтений на прямой цепи и по второй координате прочтений на обратной цепи). С использованием опции -n возможна сортировка картирований по имени.

Сортировка требует больших затрат памяти, поэтому для экономии ресурсов программа записывает промежуточные BAM-файлы. По умолчанию используется порог в 768 MiB на один поток выполнения. С помощью опции -m можно менять этот порог, а с помощью -@ можно увеличивать количество потоков.

Пример использования, в котором явно указан префикс для временных файлов (-T) и выходной файл (-o):

samtools sort -T /tmp/sample_sorted -o sample_sorted.bam aln.bam

Однако до сих пор часто встречается устаревшая форма использования. В ней после названия входного файла обязательно указывается префикс для выходного файла и временных файлов. Например, команда принимает на вход файл unsorted_in.bam, записывает промежуточные файлы вида sorted_out.%d.out (%d обозначает номер временного файла) и записывает результирующий файл sorted_out.bam:

samtools sort unsorted_in.bam sorted_out

index

Команда index формирует индекс-файл на основе отсортированного BAM или CRAM файла. Индекс-файл имеет формат .cram.crai для CRAM и .bam.bai или .bam.csi для BAM. Например, он позволяет быстрее работать командам view и tview в случае, если требуется извлечь какой-то регион из ассоциированного файла.

Например, получение индекс-файла sorted.bam.bai:

samtools index sorted.bam

Некоторые команды SAMtools требуют обязательного индексирования файлов (например, idxstats).

merge

merge объединяет несколько отсортированных файлов BAM в один выходной файл, также отсортированный. Если заголовки исходных файлов отличаются, merge формирует новый, обобщённый. Если в файлах присутствуют записи в одинаковыми ID, по умолчанию к ним добавляется суффикс, позволяющий их различить.

Пример, который позволяет слить отсортированные файлы in1.bam, in2.bam и in3.bam в один out.bam:

samtools merge out.bam in1.bam in2.bam in3.bam

idxstats

idxstats выводит короткую аннотацию по входному BAM-файлу в формате значений, разделённых табуляций (TSV). В аннотации содержится информация о названиях референсных последовательностей, их длинах, количестве картированных и некартированных ридов.

Важно, что для использования этой команды входной файл должен быть отсортирован и индексирован:

samtools idxstats sorted.bam

mpileup

mpileup принимает на вход один или несколько BAM-файлов и генерирует форматы pileup, en:VCF или его бинарный вариант BCF. Они предназначены для хранения или поиска однокулеотидных полиморфизмов и инделов (вставок или делеций).

В формате pileup каждая строка представляет одну геномную позицию и содержит информацию о её покрытии, нуклеотидном составе, качестве выравнивания и других характеристиках. В формате VCF (или BCF) хранится информация только о позициях, измененных по сравнению с референсным геномом.

Пример использования для региона 1,000-2,000 третьей хромосомы референсного генома ref.fasta:

samtools mpileup -gf ref.fasta -r chr3:1,000-2,000 in1.bam in2.bam

rmdup

rmdup удаляет потенциальные дупликации ПЦР. По умолчанию команда работает для парных прочтений. Если несколько пар прочтений имеет идентичные координаты, программа оставляет только одну пару с наилучшим качеством. С помощью опции -s и -S команда работает и с одиночными прочтениями.

Пример использования:

samtools rmdup -s sorted.bam sorted_nodup.bam

fixmate

fixmate позволяет дополнять файл с парными прочтениями: заполняет координаты соответствующих пар, добавляет тэг ISIZE (отвечает за расстояние между двумя парами) и другие тэги, отвечающие за информацию о парах. С помощью опции -r можно удалять пары с некартированными прочтениями.

Пример использования:

samtools fixmate in.namesorted.sam out.bam

calmd

Команда calmd дополняет записи входного файла тэгом MD, в котором содержится информация о полиморфизмах и инделах.

phase

Команда phase позволяет находить и фазировать гетерозиготные полиморфизмы.

flags

flags конвертирует флаги картирования из численного в текстовое представление и наоборот.

Например, команда

samtools flags PAIRED,UNMAP,MUNMAP

выведет на экран:

0xd 13 PAIRED,UNMAP,MUNMAP

bam2fq

Конвертирует BAM в FASTQ.

Пример:

samtools bam2fq input.bam > output.fastq

Примеры использования

Преобразование SAM в отсортированный BAM:

samtools view -bS sample.sam | samtools sort - file_sorted

Извлечение только картированных прочтений:

samtools view -h -F 4 -b sample.bam > sample_only_mapped.bam

Извлечь только первые картирования из пар можно с помощью фильтрации по соответствующему флагу:

samtools view -h -f 0x0040 sample.bam > sample_first_pair.sam

Случайное семплирование файла, или выбор случайных прочтений, можно произвести с помощью параметра -s команды view. Параметр имеет формат числа с плавающей точкой, причём до точки указывается зерно генератора случайных чисел, а после точки — доля прочтений исходного файла, определяющая объём результата:

samtools view -s 0.5 -b sample.bam > random_half_of_file.bam

Подсчет всех прочтений:

samtools idxstats in.bam | awk '{s+=$3+$4} END {print s}'

Подсчет только картированных прочтений:

samtools idxstats in.bam | awk '{s+=$3} END {print s}'

Примечания

Литература

• Li, Heng, et al. «The sequence alignment/map format and SAMtools». Bioinformatics 25.16 (2009): 2078—2079.
• Li, Heng. «Improving SNP discovery by base alignment quality». Bioinformatics 27.8 (2011): 1157—1158.
• Ramirez-Gonzalez, Ricardo, et al. «Bio-samtools: Ruby bindings for SAMtools, a library for accessing BAM files containing high-throughput sequence alignments». Source code for biology and medicine 7.6 (2012).