Анализ транскриптомов

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

Команда	Функция
Анализ чтений
fastqc chr2.1.fastq	Выдает информацию о ридах, находящихся в файле chr2.1.fastq, в графическом виде; в частности график с "ящиками с усами", позволяющий дать оценку качеству нуклеотидов в ридах
java -jar /nfs/srv/databases/ngs/suvorova/trimmomatic/trimmomatic-0.30.jar SE -phred33 chr2.1.fastq chr2.1edited.fastq TRAILING:20 MINLEN:50	Отрезает с конца каждого рида в файле chr2.1.fastq нуклеотиды с качеством ниже 20 и оставляет риды длиной не менее 50 нуклеотидов; выходные данные записывает в chr2.1edited.fastq

FastQC: "Per base quality" до чистки

Число ридов: 12507

FastQC: "Per base quality" после чистки

Число ридов: 12399

Из рисунков видно, что после чистки среднее качество нуклеотидов незначительно улучшилось, причем преймуществено во второй половине для ридов.

Часть II: картирование чтений

Команда	Функция
Картирование чтений
hisat2 --no-softclip ../hisat2-build/hisat2-build ../chr2.1edited.fastq -S hisat2-build.1.sam	Строит выравнивание (выходной файл - hisat2-build.1.sam) ридов, прошедших очистку (chr2.1edited.fastq), и индексированной последовательности второй хромосомы (=референсная последовательность); --no-softclip запрещает подрезание чтений; --no-spliced-alignment, обеспечивающий картирование без разрывов, был использован для аналогичного задания в предыдущем практикуме, но при картировании транскриптома разрывы нужны, так как они обозначают места интронов, поэтому этот параметр здесь не использовался.
Анализ выравнивания
samtools view hisat2-build.1.sam -b -o hisat2-build.1.bam	Переводит выравнивания из формата .sam в бинарный формат .bam; -b указывает на формат .bam выходного файла (hisat2-build.1.bam), -o - на сам выходной файл.
samtools sort hisat2-build.1.bam hisat2-build.1sorted	Сортирует выравнивания чтений с референсной последовательностью; входной файл - hisat2-build.1.bam, выходной - hisat2-build.1sorted.bam.
samtools index hisat2-build.1sorted.bam	Индексирует выравнивания в формате .bam; входной файл - hisat2-build.1sorted.bam; выходной - hisat2-build.1sorted.bam.bai.

Вывод hisat2:

12399 reads; of these: 12399 (100.00%) were unpaired; of these: 54 (0.44%) aligned 0 times 12345 (99.56%) aligned exactly 1 time 0 (0.00%) aligned >1 times 99.56% overall alignment rate

Видим, что из 12399 чтений было откартировано 12345, причем все по одному разу. 54 ридов не было откартировано.

Часть III: анализ чтений

Команда	Функция
Подсчет чтений
htseq-count -f bam -s no -i gene_id -m union hisat2-build.1sorted.bam /P/y14/term3/block4/SNP/rnaseq_reads/gencode.v19.chr_patch_hapl_scaff.annotation.gtf	Программа считает сколько ридов из выравнивания в файле hisat2-build.1sorted.bam откартировалось на каждый ген из файла с разметкой человеческого генома по версии Gencode19 gencode.v19.chr_patch_hapl_scaff.annotation.gtf. -f обозначает формат входного файла (здесь: .bam). -s показывает, взят ли транскрипт с какой-то конкретной цепи ДНК (в данном случае мы не знаем, поэтому no). -i показывает, какой атрибут файла .gtf будет использоваться, чтобы присвоить ему риды, если они откартировались на месте, описываемым этим атрибутом (здесь: gene_id, из базы данных Ensembl). -m может принимать три значения: 1)union - либо присваивает рид определенный gene_id, если он откартировался на него, либо считает его как ambiguous, если рид перекрывается с двумя генами одновременно, либо alignment_not_unique, если рид выравнивается на больше чем одно место в референсе, либо not_aligned, если в .bam файле с выравниванием они не были выравнены, либо too_low_aQual, если рид был выкинут параметром -a (выкидывает риды с качеством выравнивания ниже указанного, по умолчанию - 10), либо no_feature если рид не смог присвоиться ни к одному gene_id; 2)intersection-nonempty - все то же самое что для union, но если рид целиком попал в один ген, а с другим перекрылся частично, то он считаетсч как gene_id того гена, с которым полностью перекрылся (для union он ambiguous); 3)intersection-strict - то же самое что для intersection-nonempty, но если рид только частично перекрылся с геном, а частично не перекрылся ни с каким другим или же часть этого рида выравнялась с одним фрагментом гена, а другая часть - с другим фрагментом этого же гена, то этот рид считается как no_feature (intersection-nonempty присваивал такие риду gene_id того гена, с которым частично этот рид перекрывается).

Вывод htseq-count (приведены только те gene_id, для которых значение найденных ридов не ноль):

ENSG00000115053.11 12058 ENSG00000206885.1 12 ENSG00000202400.1 5 ENSG00000233538.1 2 __no_feature 268 __ambiguous 0 __too_low_aQual 0 __not_aligned 54 __alignment_not_unique 0

Видим, что 54 чтений не были выровнены еще командой hisat2, а 268 не легли в границы генов. В число невыровненных могли войти переведенные в ходе обратной транскрипции в ДНК-последовательность полиА-хвосты мРНК.

Большинство ридов перекрываются с геном NCL (ENSG00000115053.11, координаты гена - chr2:231,453,531-231,483,641), кодирующим белок нуклеолин. Координаты двух других генов лежат в пределах координат гена NCL, они кодируют некодирующие РНК: ENSG00000206885 (chr2:455,800-231,455,936) и ENSG00000202400 (chr2:231,460,371-231,460,440) кодируют малые ядрышковые РНК SNORA75 и SNORD82 соответственно. Ген ENSG00000233538 лежит до гена NCL (chr2:231,452,195-231,453,153) кодирует AC017104.3, lincRNA. Скорее всего те риды, что были отмечены как no_feature, попали в область chr2:231,453,153-231,453,531 - это между координатами конца гена AC017104.3 и начала гена NCL.