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 дублированных чтений