bwa mem -t 10 ../../bwa/chr6.fna ../trimmomatic_filtration/SRR10720416_1_paired.fq.gz
../trimmomatic_filtration/SRR10720416_2_paired.fq.gz > SRR10720416.sam 2> bwa_mem.log
Версия 0.7.17-r1188, log-файл.
Подача на вход двух файлов (_1 и _2) автоматически распознаётся как парные чтения. Параметр -t зачаёт число используемых потоков выполнения (threads), файл chr6.fna содержит референс.
Файл .sam весит 16G.
2.1 Устройство файла
Первая строка содержит название референсной последовательности и её длину. Вторая строка содержит информацию о программе: название, версию, все параметры, с которыми она была запущена.
Далее тело файла:
Первый столбец содержит название временное последовательности.
Второй - флаг, шифрующий свойства чтения и картирования (парное ли чтение, первое или второе из пары, картированы ли оно и его пара и др.).
Третий столбец - название последовательности референса, на которую картирован фрагмент.
Четвёртый - позиция картирования (координата первого нуклеотида, самая левая, если картировался на несколько мест).
Пятый - качество картирования.
Шестой - информацию о выравнивании в формате CIGAR (для блока подряд идущих совпадающих/несовпадающих/делетированных нуклеотидов указывается тип совпадения и длина блока, но не сами нуклеотиды).
Седьмой - название пары/следующего чтения.
Восьмой - позицию пары/следующего чтения.
Девятый - длину временно рассматриваемого сегмента: если несколько картирований на одном референсе, то это длина от самого левого картированного нуклеотида до самого правого, иначе 0.
Десятый - сама последовательность.
Одиннадцатый - строка, описывающая качество прочтения понуклеотидно. Дальше идёт всякая дополнительная информация.
Если какого-либо значения нет, то строки заменяются на "*", а числа на "0".
2.2 Анализ части файла
Чтения в теле фала расположены парами (с совпадающими названиями). В той части файла, которую я посмотрела (по 25 строк с начала и конца) довольно много некартированных пар чтений, у которых заполнены только название, флаг(77 для первого, 141 для второго), последовательность и качество прочтения. Есть чтения, у которых из пары картировано только одно, например, флаги 117 и 185
Картированные встречаются среди них двойками, описывающими парные чтения: они имеют соответствующие флаги (например, 117 и 185), тогда нули и * будут только для одного из чтений, но название референса и позиция у него будут заполнены согласно значению для его пары.
Есть и те, которые картировались оба, например, с флагами 99 и 147, в этом случае нормально заполнены все поля, а TLEN у них одинаковые со знаком + у левого и знаком - у правого.
Во всех случаях, которые я видела, название пары совпадает с названием самого чтения, поэтому в седьмом столбце стоит "=".
3.1 Получение
samtools sort -o SRR10720416.bam SRR10720416.sam &> samtools_sort.log
Версия samtools (здесь и ниже) 1.9 (using htslib 1.9), log-файл.
-o задаёт имя итогового файла (иначе выдача производилась бы в STDOUT).
Файл .bam весит 4.6G
3.2 Индексирование
samtools index SRR10720416.bam &> samtools_index.log
log-файл (пуст)
Дополнительные параметры не указаны, поэтому индекс будет в формате bai.
samtools flagstat SRR10720416.bam &> samtools_flagstat.log
| Картировано | 11 516 782 (14.8%) |
| Картировано корректно | 8 475 774 (10.9%) |
Процент считается в первой строке от чсех чтений в файле, во второй - от чтений, для которых есть пара (paired in sequencing).
Числа отличаются, потому что в некоторых парах чтений могло не картироваться одно из пары или картироваться далеко - например, если участок, на который картировалось первое, на самом деле не является тем нужным, а просто содержит случайно похожий фрагмент. Или если рядом случайно оказался участок, похожий на пару, но в другом направлении, поэтому они картировались близко, но неправильно.
5.1 Выделение нужной хромосомы
samtools view -b SRR10720416.bam NC_000006.12 > SRR10720416_chr6.bam 2> samtools_view.log
log-файл(пуст)
NC_000006.12 в данном случае - название моей хромосомы. Поскольку в качестве референса я использовала фал только с её последовательностью, то по сути эта команда уберёт чтения, не картировавшиеся никуда (хотя, как видно из анализа sam файла, вероятно останется второе чтение из пары, если первое картировалось на хромосому), а все картировавшиеся корректно гарантированно находятся на хромосоме и не исчезнут при фильтрации.
5.2 Анализ чтений на хромосоме
samtools flagstat SRR10720416_chr6.bam &> samtools_flagstat_chr6.log
| Картировано | 11 516 782 (14.8%) |
| Картировано корректно | 8 475 774 (10.9%) |
Видно, что поменялось только общее число, но не число картированных или картированных корректно.
5.3 Фильтрация по корректности
samtools view -f 0x2 -bS SRR10720416_chr6.bam > SRR10720416_chr6_f.bam 2> samtools_view_f.log
log-файл(пуст)
5.4 Анализ после фильтрации
samtools flagstat SRR10720416_chr6_f.bam &> samtools_flagstat_chr6_f.log
| Картировано | 8 479 089 (100%) |
| Картировано корректно | 8 475 774 (100%) |
Почему-то общее число картированных корректно увеличилось: на предыдущих шагах было 8 475 774, а теперь стало 8 479 089.
5.5 Индексация итогового файла
samtools index SRR10720416_chr6_f.bam &> samtools_index_chr6_f.log
log-файл(пуст)
picard MarkDuplicates -M metrix.file.txt -I SRR10720416_chr6_f.bam
-O marked_SRR10720416_chr6_f.bam &> picard_MarkDuplicates.log
Версия 2.23.3, log-файл.
В log-файле есть предупреждение о том, что название чтения не начинается с числа, но выполнить программу это ему не помешало.
samtools flagstat marked_SRR10720416_chr6_f.bam &> samtools_flagstat_marked.log
Дубликатов найдено 511 848.