Главная


Практикум №13: "Предсказание генов прокариот"



Часть 1.


1. Анализ качества чтений (программа FastQC).
Команда: fastqc chr3.fastq (chr3.fastq — имя файла с чтениями).
Получен архив (.zip), который содержит отчет о программе в виде html файла

2. Очистка чтений (программа Trimmomatic).
Сначала были удалены основания с низким качеством (менее 20) с 3'-конца с помощью программы TRAILING:
java -jar /usr/share/java/trimmomatic.jar SE -phred33 chr3.fastq trailing.fastq TRAILING:20
Полученный файл: trailing.fastq
Затем были удалены риды с длиной менее 50 с помощью команды MINLEN: java -jar /usr/share/java/trimmomatic.jar SE -phred33 trailing.fastq a.fastq MINLEN:50
Полученный файл: a.fastq

Далее был выполнен анализ качества очищенных чтений с помощью FastQC:
  • 20932 чтений было до очистки, 20570 - после (362 рида удалены).

    Рис. 1. Per base quality до очистки.
    График поделён на три области: зелёную, жёлтую и красную,
    соответствующие качеству ридов.


    Рис. 2. Per base quality после очистки.

  • Из Рис. 1. и Рис. 2. видно, что Trimmomatic удалил чтения имеющие "некачественные концы". После очистки уже не наблюдается позиций, "залезающие" в красную область, а в жёлтую область попала лишь последняя позиция. Вывод: действительно, был произведён отбор наиболее качественных чтений.

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



  • Картирование чтений осуществлялось командой BWA.
    >Была проиндексирована референсная последовательность. Было проведено 100 итераций. Команда представлена на Рис. 3.

    Рис. 3. Команда index bwa.


    Далее с помощью команды bwa mem было построено выравнивание референса и ридов. Полученное выравнивание - в формате .sam ЗДЕСЬ.
  • Анализ выравнивания.
    Использовался пакет samtools.
    Для начала я перевёл формат .sam выравнивания в формат .bam. Команда samtools view представлена на Рис. 4.

    Рис. 4. Команда samtools view с опциями -b (на выходе бинарный файл) и -o ( указание имени выходного файла).

    Полученный .bam файл: ЗДЕСЬ.

    Затем выравнивание референсной последовательности с ридами было отсортировано по координате в референсе начала чтения; команда samtools sort (Рис. 5).

    Рис. 5. Команда samtools sort.

    Полученный .bam файл: out.prefix.bam.

    Далее полученный файл был проиндексирован командой samtools index (Рис. 6).

    Рис. 6. Команда samtools index.

    И вот, наконец, я узнал, сколько чтений было картировано на геном. Команда - на Рис. 7.

    Рис. 7. Команда samtools idxstats. Число картированных ридов выделено красным, некартированных - синим.


    Выдача происходит следующим образом: название последовательности, длина, число картированных ридов, число некартированных ридов.
    То есть 20569 чтений было картировано, 0 было некартировано. Результат был для меня удивительным, ведь после очистки осталось 20570 чтений, и куда делся один рид - непонятно...

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

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

    Был создан файл с полиморфизмами в формате .bcf (polymorphism.bcf); команда samtools mpileup -uf (Рис. 8). Затем командой "bcftools call -cv" пакета bcftools был создан файл со списком отличий между референсом и чтениями в формате .vcf (dif.bcf)

    Рис. 8. Команда samtools mpileup и bcftools call.

  • Описание полиморфизмов.

    Всего найдено 235 полиморфизмов, из которых 219 - замены, а 16 - индели (Рис. 9).

    Рис. 9. Три полиморфизма из файла dif.vcf.


    На Рис. 9. представлены 3 полиморфизма, каждый из которых - замена одного нуклеотида другим.
  • В первом полиморфизме (координата: 41291081) в референсе был гуанин, в чтении заменился на аденин. Качество чтения здесь 221.999 (хорошее!), глубина 25 (хорошее покрытие).
  • Во втором полиморфизме (координата: 41309334) в референсе был цитозин, в чтении заменился на аденин. Качество чтения здесь 4,12848 (не очень так), глубина 3 (так себе).
  • В третьем полиморфизме (координата: 41310211) в референсе был гуанин, в чтении заменился на аденин. Качество чтения здесь 5.46383 (тоже не ахти качество), глубина 1 (плохое).

    Аннотация SNP.


    Использовался annovar.
    Для начала файл .vcf был конвертирован в avinput формат, с которым может работать annovar. Команда:
    perl convert2annovar.pl -format vcf4 ~/dif.vcf > ~/diff.avinput

    Далее был аннотирован файл с snp с помощью предложенных баз данных с использованием скрипта annotate_variation.pl
    БД dbsnp (filter-based annotation).

    Команда:
    perl annotate_variation.pl -filter -out ~/rs -build hg19 -dbtype snp138 ~/diff.avinput humandb/

    С помощью данной команды было выявлены SNP, имеющие rs. Получено два файла:
  • dropped - здесь snp, имеющие rs, то есть те, которые аннотированы в dbSNP. Для таких SNP указываеся их идентификатор (rs) в данной БД. Таких: 180.
  • filtered - здесь snp без rs (неаннотированые в dbSNP). Таких: 40.
    БД refgene (gene-based annotation).

    Команда:
    perl annotate_variation.pl -out ~/ref.out -build hg19 ~/diff.avinput humandb/

    Было получено 3 файла, из которых наиболее важны два:
  • variant_function - в этом файле информация о том, куда попали snp - в экзоны, интроны или 5'-, 3'-некодирующие области. Так же присутствует информация о координатах snp в хромосоме, о качестве и глубине полиморфизмов и в какие гены попали snp.
  • exonic_variant_function - в данном файле информация о snp, попавших в экзоны, синонимичная или несинонимичная замена, к какой аминокислотной замене привёл данный полиморфизм, и в какие гены попали snp.

    база данных refseq в annovar делит snp на две категории: первая (файл variant_function) - все аннотированные snp, попавшие в в экзоны, интроны или 5'-, 3'-некодирующие области.
    И вторая группа (файл exonic_variant_function) - только snp, попавшие в экзоны.

    В гены попало 13 snp: 10 - в ген ULK4 (кодирует одну из серин-треониновая киназа 4, которая фосфоририрует гидроксильную группу в остатках серина или треонина), 3 полиморфизма попали - в ген FNDC3B. Продукт гена - домен фибронектина III типа, найден во множестве экстрацеллюлярных белков человека, в частности в фибронектине - гликопротеине внеклеточного матрикса, учавствует в заживлении ран, необходим для эмбриогенеза, если инактивировать ген фибронектина, это может привести к гибели эмбриона. Таким образом, данный snp может иметь большое клиническое значение.
    БД 1000 genomes (filter-based annotation).

    Команда:
    perl annotate_variation.pl -filter -dbtype 1000g2014oct_all -buildver hg19 -out ~/1000.out ~/diff.avinput humandb/

    Было получено 3 файла, из которых наиболее важны два:
  • 1000dropped - в этом файле собраны SNP, которые аннотированы в БД 1000 Genomes Project (2014 Oct). Так же для каждого такого SNP указана частота его встречаемости. Таких было 177.
  • 1000filtered - в этом файле собраны SNP, которые неаннотированы в БД 1000 Genomes Project (2014 Oct). Таких оказалось 43.
    БД Gwas (region-based annotation).

    Команда:
    perl annotate_variation.pl -regionanno -build hg19 -out ~/gwas.out -dbtype gwasCatalog ~/diff.avinput humandb/

    Было получено 2 файла, из которых наиболее важен один:
  • GWAS - в этом файле собраны SNP, которые, имеют отношение к каким-либо заболеваниям или особенностям, которые были выявлены при проведении геномных исследований.
    Таких SNP было найдено четыре, которые ассоциированы с уровнем адипонектина, давлением крови и ростом человека (два SNP).
    БД Clinvar (filter-based annotation).

    Команда:
    perl annotate_variation.pl ~/diff.avinput humandb/ -filter -dbtype clinvar_20150629 -buildver hg19 -out ~/clinvar.out

    Данная БД позволяет аннотировать SNP, которые связаны с болезнями человека.
    Выполнение данной команды иллюстрирует Рис. 10.

    Рис. 10.

    Было получено 3 файла, из которых наиболее важны два:
  • dropped - файл, в котором собраны snp, аннотированные данной БД, как snp, имеющие отношение к какому-либо заболеванию.
  • filtered - файл, в котором собраны snp, неаннотированные данной БД.

    К моему удивлению, все snp находились файл filtered, а в файл dropped был пустой. Я несколько раз перезапускал - результат тот же.
    То есть, видимо, ни один snp, по мнению Clinvar, не имеет отношения к заболеваниям человека, что довольно странно, так как по БД Gwas мы получили четыре snp, ассоциированые с уровнем адипонектина, давлением крови и ростом человека.

    Далее необходимо было составить сводную таблицу, в которую бы вошли все snp и их характеристики по использованным для аннотации базам данных.
    Но у меня нашлось слишком много SNP (220!), поэтому я составил такую таблицу только для тех snp, которые попали в экзоны, так как они, как можно предполагать, представляются наиболее интересными.
  • Таблица

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

    Команда
    Результат
    fastqc chr3.fastq анализ качества ридов
    java -jar /usr/share/java/trimmomatic.jar SE -phred33 chr3.fastq trailing.fastq TRAILING:20 Очистка чтения: удаляет с конца рида нуклеотиды с качеством ниже 20
    java -jar /usr/share/java/trimmomatic.jar SE -phred33 trailing.fastq a.fastq MINLEN:50 Очистка чтения: удалены риды с длиной менее 50 нуклеотидов.
    bwa index chr3.fasta Индексирует референс
    bwa mem chr3.fasta a.fastq > BW.sam Картирование ридов на хромосому 3: выравнивание строит чтений с референсной последовательностью
    samtools view BW.sam -b -o BW.bam Перевод файла с выравниванием в бинарный формат (.bam)
    samtools sort BW.bam out.prefix сортирует выравнивание чтений с референсом.
    samtools idxstats out.prefix.bam Определяет количество картированных чтений
    samtools mpileup -uf chr3.fasta out.prefix.bam > polymorfism.bcf Создание файла с полиморфизмами
    bcftools call -cv polymorfism.bcf > dif.vcf Перевод .bcf файла в формат .vcf.
    perl convert2annovar.pl -format vcf4 ~/dif.vcf > ~/diff.avinput Перевод файла .vcf в файл .avinput, с которым может работать Annovar.
    perl annotate_variation.pl -filter -out ~/rs -build hg19 -dbtype snp138 ~/diff.avinput humandb/ Аннотация файла со списком SNP по БД refgene.
    perl annotate_variation.pl -out ~/ref.out -build hg19 ~/diff.avinput humandb/ Аннотирует файл со списком SNP по базе данных dbsnp.
    perl annotate_variation.pl -filter -dbtype 1000g2014oct_all -buildver hg19 -out ~/1000.out ~/diff.avinput humandb/ Аннотирует файл со списком SNP по базе данных 1000genomes.
    perl annotate_variation.pl -regionanno -build hg19 -out ~/gwas.out -dbtype gwasCatalog ~/diff.avinput humandb/ Аннотация файла со списком SNP по базе данных GWAS.
    perl annotate_variation.pl ~/diff.avinput humandb/ -filter -dbtype clinvar_20150629 -buildver hg19 -out ~/clinvar.out Аннотация файла со списком SNP по базе данных Clinvar