| На главную |
Практикум 12
Часть 1 : Подготовка чтений
| |
|
fastqc chr21.1.fastq |
|
java -jar /nfs/srv/databases/ngs/suvorova/trimmomatic/trimmomatic-0.30.jar SE -phred33 chr21.1.fastq res_trimm_21_1.fastq MINLEN:50 |
|
fastqc res_trimm_21_1.fastq |
|
При чистке из всех чтений были оставлены те, что имеют длину не менее 50 нуклеотидов. Плохие нуклеотиды с конца не обрезались, потому что их не было в "необработанном варианте"
До чистки - 11221 чтений, после чистки - 11158 чтений. Для второй реплики до чистки 7650, после чистки - 7597.
box plot распределения качества секвенирования нуклеотидов для всех ридов из неотфильтрованного набора
box plot распределения качества секвенирования нуклеотидов для всех ридов из отфильтрованного набора
Часть 2 : Картирование чтений
| |
|
hisat2-build chr21.fasta chr21_index_base |
|
hisat2 --no-softclip -x chr21_index_base -U res_trimm_21_1.fastq -S mapped_reads_1.sam |
|
samtools view -b -o mapped_reads_1.bam mapped_reads_1.sam |
|
samtools sort -T ./tmp/sorted.tmp -o sorted_reads_1.bam -O bam mapped_reads_1.bam |
|
samtools index sorted_reads_1.bam |
|
281 ридов из 11158 не было откартировано, 10877 откартировались на геном один раз. Для второй реплики 167 не откартировались ни разу, 7430 ридов откартировались один раз (всего 7597 ридов).
| |
|
htseq-count -f bam -s no -i gene_id -m intersection-nonempty sorted_reads_2.bam /P/y14/term3/block4/SNP/rnaseq_reads/gencode.v19.chr_patch_hapl_scaff.annotation.gtf -o counted_1.txt |
|
Параметр -f bam указывает на формат входного файла с откартированными и отсортированными ридами, в моем случае он бинарный - bam
Параметр -i gene_id означает какой именно аттрибут будет использован для характеристики места картирования рида. Значение gene_id означает, что из столбца с аттрибутами в фалйле gtf для характеристики места картирования мы берем значение gene_id (так как разметка генома создана на основе ensembl, то там находятся id генов из этого геномного браузера), помимо этого можно было взять transcript_id value.
Параметр -m intersection-nonempty означает особенность присваивания риду той или иной метки. C таким значением этого параметра, программа присваивает риду gene_id, если рид хоть как-то перекрывается с геном (в том случае когда рид полностью лежит внутри гена Х и частично внутри гена У, риду присваивается gene_id гена Х). Если же рид полностью находится внутри двух проаннотированных генов одновременно, риду присваивается метка __ambiguous, а если картируется неоднозначно - метка __alignment_not_unique. Помимо такого значения параметра -m можно выбрать значение union, и он будет отличаться тем, что если рид любым образом захватывает несколько генов, ему присваивается метка __ambiguous. Также можно выбрать значение intersection-strict и тогда программа будет присваивать риду конкретный gene_id только если рид целиком лежит внутри этого гена.
Параметр -s означает цепь, в которой мы будем искать гены на которые картировались риды. Это полезно в том случае, если бы исследование было специфичным относительно цепи, с которой эти риды получены (были бы отброшены полученные по ошибке риды с другой цепи), но в нашем случае мы анализируем транскриптом и ничего не знаем о той цепи, с которой этот транскрипт получается, поэтолу нужно установит значение параметра no
Для первой реплики с данными параметрами скрипт нашел 408 ридов не откартировавшихся никуда (__no_feature), 281 неоткартированный рид (__not_aligned). Остальные риды программа приписала к трем разным генам: ENSG00000156256.10, ENSG00000156261.8, ENSG00000231125.2. (информация о том куда картированы риды записывается при помощи параметра -o counted_1.txt, подсчет ридов принадлежащих разным генам, определение и подсчет меток, присвоенных ридам производился при помощи функий bash: sort -u wc -l grep) Скорее всего те люди, которые делали это исследование хотели отсеквенировать матричную РНК гена CCT8 (chaperonin containing TCP1 subunit 8), и большая часть ридов, судя по картинке из IGV действительно ложится на экзоны этого гена. Но часть ридов легла на экзон находящегося рядом в геноме гена USP16 (ubiquitin specific peptidase 16) (но на прямой цепи, вероятно это часть мРНК CCT8 но за счет того что мы разрешали htseq-count рассматривать гены на любой из цепей, он определил риды 3'UTR CCT8 как часть гена USP16). Третий найденный ген это некодирующая РНК закодированная в интроне CCT8 - AF129075.5. Те чтения, что откартировались на геном, но не легли на границы генов просто находятся в 3'UTR гена CCT8 между его последним экзоном и геном USP16. CCT8 кодирует субъединицу шаперонина - сложного молекулярного комплекса напоминающего бочонок и имеющего внутри гидрофобную поверхность. Его функция заключается в помощи сворачиванию новосинтезированных белков.
Риды, картированные на геном
Сравнение реплик
| |
|
|
| |
|
|
| |
|
|
| |
|
|
| |
|
|
| |
|
|
| |
|
|
| |
|
|
Мы видим что в первой реплике относительно второй примерно в полтора раза больше ридов. С тем же коэффицентом пропорциональности (при увеличении числа ридов) увеличивается только число ридов, картирующихся на CCT8, что подтверждает выдвинутое ранее предположение об мРНК этого гена как целевой мРНК.
Другой способ подсчета чтений
| |
|
|
| |
|
|
| |
|
|
| |
|
|
| |
|
|
| |
|
|
| |
|
|
Так как параметр intersection-strict требует полного пересечения рида с геном, количество ридов картирующихся на гены уменьшается.
© Кристина Перевощикова, 2017