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

Часть 1. Подготовка чтений

В данном практикуме мы работали с чтениями экзома, картирующимися на участок хромосомы человека.
Файл chr10.fastq с одноконцевыми чтениями в формате fastq был взят мной из директории /P/y14/term3/block4/SNP/reads.
Подготовка чтений осуществлялась в несколько этапов:
1. Анализ качества прочтений
2. Очистка чтений с последующим анализом и сравнением параметров прочтений до и после очистки

Пункт 1. Анализ качества прочтений

Анализ прочтений был осуществлен с помощью прогрммы FastQC, установленной на kodomo. Результатом работы программы является zip архив и html файл, содержащий отчет о качестве прочтений.
Программа была запущена с помощью команды:
fastqc chr10.fastq

Пункт 2. Очистка чтений с последующим анализом и сравнением параметров прочтений до и после очистки

Очистка чтений проводилась с использованием программы Trimmomatic. Были удалены чтения с длиной менее 50 нуклеотидов и с концов каждого чтения были отрезаны нуклеотиды с качеством ниже 20.
Ниже приведена команда:
java -jar /nfs/srv/databases/ngs/suvorova/trimmomatic/trimmomatic-0.30.jar SE -phred33 chr10.fastq chr10_out.fastq TRAILING:20 MINLEN:50

Сравним параметры чтений до и после чистки. До чистки было 10666 чтений, после чистки стало 10526 чтений.
Ниже на рисунке слева изображен контроль качества чтений до чистки, а на рисунке справа - после чистки.
Синяя линия - среднее качество чтений, центральные красные линии - медианы, желтые прямоугольники - интерквартальные размахи (разница между верхн. и нижн. квартилями; диапазон значений качества, при котором качество 25% чтений на данной позиции выше нижней границы, а 75% - не выше верхней).
Поле графика разделено на 3 полосы - зеленую, желтую и красную; попадание в данные полосы вышеперечисленных элементов графика позволяет сделать вывод о качестве чтений. После чистки все риды располагаются в зеленой области, следовательно, с помощью чистки мы получили надежные прочтения. Также после чистки стало меньше ридов, так как были удалены риды длинной меньше 50.

Контроль качества чтений до чистки

Контроль качества чтений после чистки

Часть 2. Картирование чтений

Чтения были откартированы с помощью программы hisat2. Программа была импортирована с помощью команды:
export PATH=${PATH}:/home/students/y06/anastaisha_w/hisat2-2.0.5

Таблица 1. Команды, использованные при картировании последовательности.

Команда	Функция	Выдача
hisat2-build chr10.fasta chr10	Индексирование референсной последовательности	Индексированный файл chr10.fasta
hisat2 -x chr10 -U chr10_out.fastq --no-spliced-alignment --no-softclip >align1.sam	Выравнивание чтений после чистки с референсной последовательностью	Файл, который содержит выравнивание формата SAM align1.sam

Программа картирования была запущена с параметрами -x (путь к индексу) -U (путь к чтениям) --no-spliced-alignment (запрет на картирование с разрывами) --no-softclip (кртирование без подрезания чтений).
Далее необходимо было проанализировать, полученное выравнивание. Для этого я использовала программу Samtools. Она работает с файлами в формате SAM.

Таблица 2. Команды, использованные для анализа последовательностей в формате sam.

Команда	Функция	Выдача
samtools view align1.sam -bo align1.bam	Программа переводит файл в формат bam	align1.bam
samtools sort align1.bam -T sorted.txt -o sorted.bam	Сортировка выравнивания чтений и референса по координате в референсе	sorted.bam
samtools index sorted.bam	Индексирование отсортированного выравнивания	sorted.bam
samtools idxstats sorted.bam > totalread.txt	Запись числа откартировавшихся чтений	totalread.txt

Выяснилось, что на хромосому откартировался 10398 рид. 128 ридов не были откартированы вообще.

Часть 3. Анализ SNP

Таблица 3. Команды, использованные при поиске SNP и инделей.

Команда	Функция	Выдача
samtools mpileup -uf chr10.fasta sorted.bam > snp.bcf	Создание файла с полиморфизмами	snp.bcf
bcftools call -cv snp.bcf -o snp.vcf	Определения различий	snp.vcf

В результате были найдены 66 snp. Из них были выбраны три и описаны в таблице. У всех высокое качество прочтения и количсетво покрытий.

Таблица 4. Найденные snp и индели.

Координата	Тип полиморфизма	В референсе	В чтениях	Покрытие в этом месте	Качество покрытия
chr10:5781628	замена	T	G	21	117.008
chr10:5766152	индель	AGTAT	AGTATGTAT		92.4668
chr10:5766337	замена	G	A	93	225.009

Следующий этап нашей работы - аннотация SNP:
Нам необходимо проаннотировать полученные snp по 5 базам данных - refgene, dbsnp, 1000 genomes, GWAS, Clinvar

convert2annovar.pl -format vcf4 snp.vcf > chr10.avinput	Переводим файл .vcf формат, удобный для работы annovar
annotate_variation.pl -filter -out SR_SNP -build hg19 -dbtype snp138 chr10.avinput /nfs/srv/databases/annovar/humandb.old/	Аннотация по Dbsnp
annotate_variation.pl -out refgen -build hg19 chr10.avinput /nfs/srv/databases/annovar/humandb.old/	Аннотация по Refgene
annotate_variation.pl -filter -dbtype 1000g2014oct_all -buildver hg19 -out 1000Genomes chr10.avinput /nfs/srv/databases/annovar/humandb.old/	Аннотация по 1000 Genomes
annotate_variation.pl -regionanno -build hg19 -out GWAS -dbtype gwasCatalog chr10.avinput /nfs/srv/databases/annovar/humandb.old/	Аннотация по Gwas
annotate_variation.pl chr10.avinput -filter -dbtype clinvar_20150629 -buildver hg19 -out CLINVAR /nfs/srv/databases/annovar/humandb.old/	Аннотация по Clinvar

При аннотации SNP по Refgen получаем три файла - refgen.exonic_variant_function, refgen.log
(системные сообщения, информация об ошибках и комментарии), refgen.variant_function
Из файла refgen.variant_function мы понимаем, что база данных Refseq делит snp по
их локализации (exonic - 10, intronic - 50, UTR3 - 4)

Гены, в которые попали наши snp:

                                                                                                                                       
FAM208B	27                                                                                                                                  
RTKN2	28 (+3)                                                                                                                                  
CASP7	6

Помимимо всего прочего в файлах есть информация о синониминости (несинонимичности) замен
аминокислотных остатков
У 9 есть rs.

Аннотация 1000Genomes показывает частоты аллелей в output - файле
Среднее значение частоты аллелей - 20,08333333 (подсчитано с помощью эксорта данных в Excel)
Также, в output - файле есть информация о характере замены (что на что было заменено - в форме гэпов и/или букв-нуклеотидов)
Также указаны координаты каждой такой замены

По аннотации Gwas видим, что 3 snp связаны с каким-либо заболеванием, либо физиологическим параметром
связь геномных признаков с фенотипическими:

                                                                                                                                       
Остеосаркома      		        5804531 	5804531 	G	A	70                                                   					                                                                                                                                                                                   
Ревматоидный артрит              	63958112	63958112	T	C	75	                                            
Vitiligo                 		115481018	115481018	C	T	91

ClinVar объединяет информацию о геномных вариациях (полиморфизмах), их отношении к здоровью человека.
В базе данных ClinVar не было какой-либо клинической информации по интересующим нас полиморфизмам.