1. Картирование чтений на референсный геном.

Для картирования чтения была использованна следуюшая команда:
bwa mem -t 15 ./bwa/chr7.fna ./DNA_reading/paired_1.fastq.gz ./DNA_reading/paired_2.fastq.gz > reads.sam 2> mem_error.txt
Описание параметров:
-t 15 - использованно 15 ядер ( это было ночью, поэтому ядра никем не использовались и я посчитал, что могу взять побольше для ускорения вычислений)
./bwa/chr7.fna - референс
./DNA_reading/paired_1.fastq.gz & ./DNA_reading/paired_2.fastq.gz- файлы с триммированными чтениями
reads.sam - название файла, полученного на выходе
2> mem_error.txt - вывод потока ошибок в файл

2. Анализ sam файла с выравниванием.

Файл содержит заголовки с тегами и таблицу. Заголовки заключают в себе информацию о чтениях, самом референсе и программу, использованную для картирования. Тогда как непосредственно в теле файла находится информация о качестве секвенированния и картирования для каждого рида.
То, что содержится в первых 11 столбцах файла представленно в таблице:
Название столбца Тип данных Описание
QNAME String ID рида
FLAG Int Характеристика выравнивания, выраженная суммой "флагов", характеризующие свойства выравнивания
RNAME String Имя референса
POS Int 1 - самая левая позиция выравнивания
MAPQ Int Качество картирования
CIGAR String Строка операций CIGAR
RNEXT String Референсное название следующего чтения
PNEXT Int Позиция следующего чтения
TLEN Int Длина участка, покрытая двумя соседними чтениями
SEQ String Последовательность рида
QUAL String Phred-score каждого нуклеотида в риде
Опциональные столбцы (Optional fields):
PG - программа
AS - score выравнивания
Вес файла reads.sam - 16GB.

3. Получение и индексирование bam файла.

Необходимо конвертировать sam файл в bam, из соображений о том, что sam файл очень много весит. После чего можно удалить sam файл, чтобы он не занимал место.
Команда для конвертирования:
samtools sort -@ 16 -o reads.bam reads.sam 2> bam_error.txt
-@ 16 - увеличивает количество тредов(я делал это в 10 утра, предварительно посмотрев командов w, что ядра свободны.
-о определяет название файла bam, а ошибки перенаправляются в файл bam_error.txt.
Вес файла bam - 4,7GB.
Команда для индексирования bam файла:
samtools index -@ 16 reads.bam

4. Анализ bam файла.

Для просмотра информации о bam файле использовал следующую команду:
samtools flagstat -@ 16 reads.bam > flagstat.txt
Информация их файла:
Чтений картированно: 10982930 (13,76%)
Картированы в корректных парах чтений: 8618528 (10.80%)
Эти числа могут отличаться, из-за того, что чтение может картироваться на разные позиции в референсе, а его парное чтение только на одно.

5. Получение картированных чтений.

Для получения чтений, картированных на хромосому используем следующую команду:
samtools view -h reads.bam NC_000007.14 > read_chr7.sam
И команда для коневертирования:
samtools view -b read_chr7.sam > read_chr7.bam
Получим информацию о файле с помощью уже известной команды:
samtools flagstat -@ 16 read_chr7.bam > flagstat_chr7.txt
Картировано 10982930 чтений (88.87%). Картировано в корректно картированных парах чтений 8618528 (70.03%) чтений.
Сравнение: численные значения не изменились, но изменился процент. ЭТо связано с тем, что тепеьр мы картируем не на весь геном, а только на определнную хромосому.
Теперь выберем только правильно картированные чтения. Команда:
samtools view -@ 16 -f 0x2 -bS read_chr7.bam > read_chr7_cor.bam
Посмотрим информацию о нём:
samtools flagstat -@ 10 read_chr7_cor.bam > flagstat_chr7_cor.txt
Картировано 8621082 чтений (100%). Картировано в корректно картированных парах чтений 8618528 (100%) чтений.
Теперь все чтения картированны. Мы добились желаемого.
Проиндексируем файл:
samtools index -@ 12 read_chr7_cor.bam

6. Подготовка выравнивания к поиску вариантов.

Команда для маркировки дублированных чтений:
picard MarkDuplicates -M metrix.file.txt -I read_chr7_cor.bam -O marked_chr7.bam
Узнаем информацию о файле с помощью flagstat:
samtools flagstat -@ 10 marked_chr7.bam > marked_flagstat.txt
Было найдено 590768 дублированных чтений