Практикум №12 и №13

Ресеквенирование. Поиск полиморфизмов у человека

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

В данном практикуме необходимо было провести картирование почищенных геномов. Картируем только парные, так как из прошлого практикума стало понятно, что в непарных просто отвратительное качество. Используемая команда:

bwa men -t 10 chr8.fna /mnt/scratch/NGS/vitbuev/bwa/chr8.fna /mnt/scratch/NGS/vitbuev/dna_reads/dna_f_p.fastq.gz /mnt/scratch/NGS/vitbuev/dna_reads/dna_r_p.fastq.gz > dna.sam 2> dna.txt (log)

Для исполнения команды я задействовал 10 ядер, в команде указан референсный файл и 2 последовательности. Выход был перенаправлен в файл dna.sam, также был создан лог-файл с ошибками.

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

Формат SAM (англ. Sequence Alignment Map) — это текстовый формат для хранения биологических последовательностей, выровненных по эталонной последовательности, также называемой референсной. SAM-файл может содержать заголовок, строки которого всегда начинаются с символа «@», за которым следует один из двухбуквенных кодов типа заголовка. В заголовке каждая строка разделена символом табуляции, и, кроме строк @CO, каждое поле данных соответствует формату тэг: значение, где тэг представляет собой двухсимвольную строку, которая определяет формат и содержимое значения. Разделы выравнивания содержат 11 обязательных полей и некоторое количество дополнительных. При помощи команды head были просмотрены первые 10 строк файла dna.sam (ссылка) 

 
@SQ	SN:NC_000008.11	LN:145138636
@PG	ID:bwa	PN:bwa	VN:0.7.17-r1188	CL:bwa mem -t 10 ../bwa/chr8.fna ../12pr/dna_f_p.fastq.gz ../12pr/dna_r_p.fastq.gz
SRR10720414.1	77	*	0	0	*	*	0	0	GAAGAAATCTGATGGTGNNNNNNANNNTCCANNNNNTGACNNNNNNNNNTCCAGTCAGATTGCTCAGGACCATCA	FHIIIHDIIIIIIGI8@######1###5>>8#####3<7;#########7>?776IIIDFHIIHIIIIIHIIDIH	AS:i:0	XS:i:0
SRR10720414.1	141	*	0	0	*	*	0	0	CAAAATCATGTCTACATAAAAAATNNNNCATTTAGTTGCTAGAGAACCCTAAAAGATGTGTTCAT	HDGFHHHGHHHHHHHHHHHH1100####<4>><5;HHHHHGGGFFGGGGGDGDGGDDGBGFHHHH	AS:i:0	XS:i:0
SRR10720414.2	77	*	0	0	*	*	0	0	TTTGCAGTGTTTCCCNANNNNNNNNNNTGCCNNNNNATCANNNNNNNNNATCTTCTGCCTTTTCTTTCNTTCATT	IIIDIIIIIIDDDEE#?##########6>?8#####<75;#########<<:<;;IIIGIFFHHI>>=#BAA?AB	AS:i:0	XS:i:0
SRR10720414.2	141	*	0	0	*	*	0	0	AATGTTTTCCCTTAGCCAATGCCCNNNNATAAGTGAACTTACCTTGGTCTTTGCTGACTCCTCCC	HIHIIIHIIIIHIIGIIIHI=@B@####<99:469HIIIIGHIHIHGEIHIDDIIFFIFHHFGGI	AS:i:0	XS:i:0
SRR10720414.3	77	*	0	0	*	*	0	0	CATGTTTTTACTCATAANNNNNNGNNNACAGNNNNNGAAANNNNNNNNNCGAATCCTCCACCGTCCCTTGCAATT	IIIIIIIIIHIIIII>>######:###:<:<#####<7;<#########699799IIGGHHGIDIGHEGIIGEEG	AS:i:0	XS:i:0
SRR10720414.3	141	*	0	0	*	*	0	0	TTCCTCTAAATAAGGGAGATCTGTNNNNTTTTCTTCTTTTTGCTTAACTTCTTCAGACTTGTTGA	IIIIIIIIIHIFIIIIHHII@=>9####:>>>8::IIIIIIHIIIFIIIIIIIHHIHIIIGHIHE	AS:i:0	XS:i:0
SRR10720414.5	83	NC_000008.11	118928912	48	75M	=	118928795	-192	TATTTCGCTCTGGGGTTCCTACGGAANNNNNNNNNCAATNNNNNACCTNNNCNNNNNNTCGTTTCCCAGCAGATG	GGCGGGEGGGDIHIBIIBII666966#########6693#####,8>6###:######:BIIHHIIDIHIIIIHI	NM:i:24	MD:Z:22A3A0A0T0A0C0C0A0A0G4T0T0A0G0T4T0C0T1C0T0C0A0A0A17	MC:Z:65M	AS:i:24	XS:i:0
SRR10720414.5	163	NC_000008.11	118928795	60	65M	=	118928912	192	CATTTCCTTTCTGAGTTAGCAGGANNNNAAAGACACTGCAATTTGTGTGTTTTCTACAGGGTGCT	IHIIIIIIIIIIIIIIIIGI@@@;####<>8CEEEHICHEIIDIAIIGHID=GAD	NM:i:4	MD:Z:24G0A0C0C37	MC:Z:75M	AS:i:57	XS:i:20
Номер Поле Описание
1 QNAME Название рида (парным будет присвоено одинаковое)
2 FLAG Представляет собой комбинацию нескольких бинарных флагов, обозначающих парность, картированность и т. д.
3 RNAME Название референса
4 POS Один минус самая левая позиция выравнивания
5 MAPQ Качество маппинга. MAPQ = −10 log10P, где P - вероятность того, что выравнивание неверно
6 CIGAR Описание выравнивания, в записи которого используется набор операций (совпадение, инсерция и др.)
7 RNEXT Название референса пары
8 PNEXT Позиция в референсе, которой соответствует начало пары. Аналогично четвёртому полю.
9 TLEN Расстояние между крайними точками выравнивания
10 SEQ Последовательность рида
11 QUAL Качество в кодировке ASCII по шкале Phred+33
12 TAGs Специальные тэги, предназначенные для хранения информации о прочтении или выравнивании

Используя команду ls -sh, я узнал, что файл весит 16 Гб

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

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

Итого, при помощи команды samtools sort -O bam -@ 10 dna.sam -o dna.bam был получен bam файл. -O bam отвечает за формат выходного файла, а -o - за файл выхода. Чтобы не ждать очень долго, количество ядер было задано "@", что существенно ускорило процесс. Файл весит 4.6 Гб. Чтобы далее работать с bam файлом, надо его проиндексировать. Команда: samtools index -b -@ 10 dna.bam dna.bai, где -b отвечает за формат файла.

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

Как я уже упоминал, файл бинарен, а значит просто так его не проанализировать. Команда: samtools flagstat -@ 10 dna.bam > dna_binary.txt (файл)

В файле находится 11.7% картированных чтений (правда не все из них картированы нормально, хорошо только 8.78%). Корректно откартировались чтения из одной пары, так как они должны картироваться недалеко друг от друга и должны быть направлены друг к другу. Число маленькое, но мы же картируем полногеномное секвенирование, картированное на одну из хромосом.

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

Чтобы получить чтения, картированные только на 8 хромосому, необходимо сначала проиндексировать файл: samtools faidx chr9.fna -o chr8.fai . Из файла .fai я узнал, что имя хромосомы - NC_000008.11.

Далее используем команды:

samtools view -h dna.bam NC_000008.11 > chr8.sam

samtools view -bS chr8.sam > chr8.bam

samtools flagstat -@ 10 dna.chr9.bam > chr8.txt (файл)

Теперь имеем 87.12% картированных чтений, 65.65% корректно картированы. Видим, что процент картированых чтений сильно вырос, хоть и не получили 100%. Скорее всего, здесь дело в работе samtools. Теперь получим только правильно картированые чтения при помощи команды samtools view -@ 10 -f 0x2 -bS chr8.bam > chr8_fixed.bam, а затем посмотрим статистику samtools flagstat -@ 10 chr8_fixed.bam > chr8_fixed.txt. Как и ожидалось, получили 100% корректно картированных чтений, значит всё сделали правильно.

Как и в прошлый раз, осталось только проиндексировать файл samtools index -b -@ 10 chr8_fixed.bam chr8_fixed.bai.

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

picard MarkDuplicates -M metrix.txt -I chr8_fixed.bam -O markdup.bam - маркировка дублированых чтений

samtools flagstat -@ 10 markdup.bam > markdup.txt (файл). Всего 603736 дублированых чтений.


Далее идёт практикум 13. В нём необходимо было получить, отфильтровать и минимально исследовать варианты, а также поработать с покрытием.

1. Получение вариантов

При помощи bcftools: bcftools mpileup --threads 10 -f /mnt/scratch/NGS/vitbuev/bwa/chr8.fna /mnt/scratch/NGS/vitbuev/12pr/markdup.bam | bcftools call -mv -o chr8.vcf был получен файл с вариантами. Mpileup создаёт файл с расширение .vcf, в котором находятся частоты встречаемости вариантов во множественном выравнивании. -threads как обычно задаёт количество ядер, -f задаёт индексированный референс хромомосомы, а затем передаём всё на выход команде call, которая и производит собственно вызов вариантов (эта команда заменяет старую view caller и вроде даже уступает ей по функционалу), где m - multiallelic-caller, v - даёт нам на вывод только варианты, ну и o - собственно выход.

При помощи знакомой команды head -40 были рассмотрены первые 40 строчек файла. 

##fileformat=VCFv4.2
##FILTER=
##bcftoolsVersion=1.9+htslib-1.9
##bcftoolsCommand=mpileup --threads 10 -f /mnt/scratch/NGS/vitbuev/bwa/chr8.fna /mnt/scratch/NGS/vitbuev/12pr/markdup.bam
##reference=file:///mnt/scratch/NGS/vitbuev/bwa/chr8.fna
##contig=
##ALT=
##INFO=
---
##FORMAT=
##FORMAT=
---
##INFO=
##bcftools_callVersion=1.9+htslib-1.9
##bcftools_callCommand=call -mv -o chr8.vcf; Date=Mon Dec 21 20:53:29 2020
#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	/mnt/scratch/NGS/vitbuev/12pr/markdup.bam
NC_000008.11	60089	.	A	AC	15.2883	.	INDEL;IDV=2;IMF=0.666667;DP=3;VDB=0.131948;SGB=-0.453602;MQSB=1;MQ0F=0;ICB=1;HOB=0.5;AC=1;AN=2;DP4=0,1,1,1;MQ=50	GT:PL	0/1:48,0,25
NC_000008.11	60111	.	A	G	47.6405	.	DP=4;VDB=0.5;SGB=-0.453602;RPB=1;MQB=1;MQSB=1;BQB=1;MQ0F=0;ICB=1;HOB=0.5;AC=1;AN=2;DP4=0,1,1,1;MQ=50	GT:PL	0/1:81,0,29
NC_000008.11	60119	.	T	G	47.721	.	DP=5;VDB=0.52;SGB=-0.453602;RPB=1;MQB=1;MQSB=1;BQB=1;MQ0F=0;ICB=1;HOB=0.5;AC=1;AN=2;DP4=0,1,1,1;MQ=50	GT:PL	0/1:81,0,31
NC_000008.11	60129	.	G	A	25.3564	.	DP=5;VDB=0.64;SGB=-0.453602;RPB=0.666667;MQB=1;MQSB=1;BQB=0.666667;MQ0F=0;ICB=1;HOB=0.5;AC=1;AN=2;DP4=1,2,0,2;MQ=41	GT:PL	0/1:58,0,80
NC_000008.11	60181	.	G	A	13.6078	.	DP=3;VDB=0.28;SGB=-0.453602;MQSB=1;MQ0F=0;AC=2;AN=2;DP4=0,0,1,1;MQ=33	GT:PL	1/1:43,6,0
NC_000008.11	60323	.	C	G	55	.	DP=7;VDB=0.703527;SGB=-0.556411;RPB=0.810265;MQB=0.607698;MQSB=1;BQB=0.810265;MQ0F=0;ICB=1;HOB=0.5;AC=1;AN=2;DP4=2,1,4,0;MQ=50	GT:PL	0/1:88,0,91
NC_000008.11	60345	.	A	G	136	.	DP=11;VDB=0.311141;SGB=-0.616816;RPB=0.783809;MQB=0.727822;MQSB=0.523791;BQB=0.503877;MQ0F=0;ICB=1;HOB=0.5;AC=1;AN=2;DP4=1,3,4,2;MQ=53	GT:PL	0/1:169,0,106
NC_000008.11	60371	.	C	T	9.80985	.	DP=13;VDB=0.154358;SGB=-0.511536;RPB=0.382304;MQB=0.857404;MQSB=0.506703;BQB=0.485672;MQ0F=0;ICB=1;HOB=0.5;AC=1;AN=2;DP4=4,5,1,2;MQ=54	GT:PL	0/1:44,0,211
NC_000008.11	60408	.	G	A	105	.	DP=9;VDB=0.551573;SGB=-0.616816;RPB=0.346719;MQB=0.346719;MQSB=0.89338;BQB=0.924584;MQ0F=0.111111;ICB=1;HOB=0.5;AC=1;AN=2;DP4=1,2,3,3;MQ=45	GT:PL	0/1:138,0,85
NC_000008.11	60422	.	C	T	158	.	DP=8;VDB=0.410943;SGB=-0.636426;MQSB=1.01283;MQ0F=0.125;AC=2;AN=2;DP4=0,0,4,3;MQ=41	GT:PL	1/1:188,21,0
NC_000008.11	60460	.	A	G	24.316	.	DP=9;VDB=0.220755;SGB=-0.511536;RPB=0.28717;MQB=0.28717;MQSB=0.5;BQB=0.28717;MQ0F=0;ICB=1;HOB=0.5;AC=1;AN=2;DP4=1,2,1,2;MQ=39	GT:PL	0/1:57,0,83
На месте прочерков стоят ещё строки INFO, но, раз на страничке это не отображается, то я их убрал. VCF формат предназначен для записи вариантов последовательностей. Сначала идут данные о файле, а затем таблица с вариантами, в которой представлено имя нашей хромосомы, ID варианта, референсный нуклеотид (если имеется), список аллельных вариантов и другое.

При помощи команды bcftools stats chr8.vcf > chr8_stats.txt была получена статистика (файл).

Получили 342714 вариантов, из которых 336648 - однонуклеотидные замены, а 6066 - индели.

2. Фильтрация вариантов

Проведём фильтрацию по качеству больше 30 и по покрытию больше 50. Команда bcftools filter -i'%QUAL>30 && DP>50' chr8.vcf >chr8_filtered.vcf, а затем посмотрим статистику bcftools stats chr8_filtered.vcf > chr8_filtered_stats.txt (файл). "Выжили" только 19758 вариантов, из которых 19512 - SNP, а 246 - индели.

3. Изучение покрытия

Теперь посмотрим на участки генома покрытые нашими чтениями хотя бы 1 раз. Команда: bedtools genomecov -bg -ibam ../12pr/markdup.bam > genomecov.bed, где параметр -bg отвечает за выход файла в формате BedGraph. В конечном файле есть таблица из 4 колонок, в которой присутсвуют название хромосомы (референса), начало и конец варианта и покрытие.

Теперь отфильтруем файл при помощи awk: awk '$4 > 50' genomecov.bed >genomecov_50.bed. В данном файле 1121075 срок.

Найдём варианты, попавшие в хорошо покрытые участки: bedtools intersect -a chr8_filtered.vcf -b genomecov.bed > final.txt. (файл) В данном файле всего 20406 строк.