На главную

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

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

    Для выполнения задания был выдан файл с одноконцевыми чтениями экзома 2 хромосомы человека chr2.fastq, и файл с последовательностью этой хромосомы из генома со сборкой версии hg19 chr2.fasta.
  Все действия выполнялись в специально созданной директории /nfs/srv/databases/ngs/p.avdiunina 

Анализ качества чтений

    Сперва был выполнен анализ качества чтений с помощью программы FastQC(использовалась версия, установленная на сервере kodomo). В командную строку была введена команда:
  
  

    Программа в той же папке выдала мне файл chr2_fastqc.html, содержащий визуализированный отчёт о работе, и архив chr2_fastqc.zip с каринками и текстовыми файлами.

Очистка чтений

 
    Далее производилась очистка чтений при помощи программы Trimmomatic, которая также была запущена с сервера kodomo следующей командой:

  

    Первая опция отрезает с конца каждого чтения нуклеотиды с качеством ниже 20, вторая удаляет риды короче 50 нуклеотидов. До чистки было 10410 чтений, после - осталось 10191 (удалено 219 чтений).
  Затем снова была запущена программа FastQC, результатом стали архив out_chr2_fastqc.zip и файл out_chr2_fastqc.html 
  С помощью команды unzip *.zip были разархивированы оба архива. Теперь выполним сравнение чтений до и после очистки.
    
    На рис. 1 представлен график качества чтений, сделанный программой FastQC до чистки. Синей линией на графике обозначено среднее качество чтений, центральными красными линиями обозначена медиана, желтыми 
  прямоугольниками - интерквартальный размах (разница между верхней и нижней квартилями, диапазон значений качества, при котором качество 25% чтений на данной позиции выше нижней границы, а 75% - не выше верхней).
  Поле графика разделено на 3 полосы зеленого, желтого и красного цветов, попадание в которые вышеперечисленных элементов графика позволяет судить о качестве чтений.
    На рис. 2 представлен график качества чтений, сделанный программой FastQC после очистки. Теперь все позиции, кроме последней, расположены в зеленой области (последняя - в желтой), из чего можно заключить, 
  что чистка прошла успешно. Еще один результат работы программы Trimmomatic состоял в удалении чтений с длиной меньше 50 нуклеотидов, теперь длины чтений колеблются от 50 до 100 нуклеотидов (раньше: от 40 до 100).

Рис. 1. "Per base sequence quality", график качества чтений, сделанный программой FastQC до чистки Рис. 2. "Per base sequence quality", график качества чтений, сделанный программой FastQC после чистки

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

Картирование чтений


    Чтения были откартированы с помощью программы  BWA - Burrows-Wheeler Alignment Tool. Использованные команды указаны в таблице. 

Команда Назначение Результат
bwa index chr2.fasta Индексирование референсной последовательности Индексированный файл chr2.fasta
bwa mem chr2.fasta out_chr2.fastq > first_align.sam Выравнивание очищенных чтений с референсной
последовательностью
Файл first_align.sam, содержащий выравнивание
в формате SAM (Sequence Alignment/Map)

Анализ выравнивания


    Для анализа полученного выравнивания применялась программа Samtools, предназначенная для обработки файлов в формате SAM. Использованные команды указаны в таблице. 


Команда Назначение Результат
samtools view first_align.sam -bo first_align.bam Перевод выравнивания в бинарный формат .bam
(опция -b меняет формат выходного файла с
установленного по умолчанию, -o обозначат имя
выходного файла).
first_align.bam
samtools sort first_align.bam -T nexact.txt -o align_sort.bam Сортировка выравнивания чтений с референсом по
координате в референсе начала чтения (опция -T
позволяет записывать временные файлы в файл
smth.txt, а не в stdout).
align_sort.bam
samtools index align_sort.bam Индексирование отсортированного выравнивания Индексированный
align_sort.bam
samtools idxstats align_sort.bam > finish.txt Выяснение числа чтений, откартированных на геном finish.txt
Оказалось, что на геном откартировалось 10191 чтение - в точности столько же, сколько осталось после очистки.

Часть III: анализ SNP

Поиск SNP и инделей


Команда Назначение Результат
samtools mpileup -uf chr2.fasta
align_sort.bam > snp.bcf
Создание файла с полиморфизмами
в формате .bcf
snp.bcf
bcftools call -cv snp.bcf -o snp.vcf Определение отличий между референсом
и чтениями
snp.vcf
Всего было найдено 78 полиморфизмов, из них 9 инделей и 69 замен. Три из них были выбраны и описаны в таблице ниже. У первого из описанных полиморфизмов хорошие покрытие и качество, у других двух - плохие.
Координата Тип
полиморфизма
Было в
референсе
Найдено
в чтениях
Глубина покрытия
данного места
Качество прочтения
на участке
55516588 Замена G C 23 184.999
55523307 Индель ACC AC 9.04379
234160243 Замена G C 7 81.0075

Аннотация SNP


    Далее необходимо было аннотировать полученные SNP с помощью программы ANNOVAR. Для этого сначала нужно было подготовить входной файл, чтобы программа могла работать с ним. В первую очередь из файла snp.vcf были удалены 
  все индели, полученный файл - snp_noindel.vcf. Затем был использован скрипт convert2annovar.pl с помощью команды perl /nfs/srv/databases/annovar/convert2annovar.pl -format vcf4 /nfs/srv/databases/ngs/p.avdiunina/snp.vcf >
  /nfs/srv/databases/ngs/p.avdiunina/snp.avinput.
    Далее были проведены аннотирования по БД, указанным выше. Команды:
  • dbnsp: perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -out rs.snp -build hg19 -dbtype snp138 snp.avinput /nfs/srv/databases/annovar/humandb/
  • 1000genomes: perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -out rs.1000g -buildver hg19 -dbtype 1000g2014oct_all snp.avinput /nfs/srv/databases/annovar/humandb/
  • GWAS: perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -out rs.gwas -build hg19 -dbtype gwasCatalog snp.avinput /nfs/srv/databases/annovar/humandb/
  • Clinvar: perl /nfs/srv/databases/annovar/annotate_variation.pl -filter -out rs.clinvar -dbtype clinvar_20150629 -buildver hg19 snp.avinput /nfs/srv/databases/annovar/humandb/
  • Refgene: perl /nfs/srv/databases/annovar/annotate_variation.pl -out rs.refgene -build hg19 snp.avinput /nfs/srv/databases/annovar/humandb/
Из всех SNP 70 имеют rs, 9 не имеют. Наименьшая частота SNP = 0.0227636, наибольшая частота = 0.994209. (rs.snp.hg19_snp138_dropped). В клинической аннотации SNP представлено 4 замены. Первая влияет на рост, вторая и третья увеличивает вероятность возникновения болезни Крона, а четвертая - рака простаты. (rs.snp.hg19_snp138_filtered). В БД RefSeq SNP подразделяется на exonic (variant overlaps a coding), splicing (variant is within 2-bp of a splicing junction), ncRNA (variant overlaps a transcript without coding annotation in the gene definition), UTR5 (variant overlaps a 5' untranslated region), UTR3 (variant overlaps a 3' untranslated region), intronic (variant overlaps n intron), upstream (variant overlaps 1-kb region upstream of transcription start site), downstream (variant overlaps 1-kb region downstream of transcription start site), intergenic (variant is in integric region). БД RefSeq выделила 35 het замен и 43 hom. Из них в UTR3 попало 7 замен, в intronic - 62, в UTR5 - 2, в exonic - 7 и в ncRNA - одна замена. Замены затронули гены CCDC88A (UTR3, intronic SNP), ATG16L1 (UTR5, intronic, UTR3, exonic SNP) и MLPH (intronic, exonic, ncRNA SNP).

Источники:


© Avdiunina Polina, 2015