RNA-seq

В этом практикуме мы должны были проанализировать одноконцевые чтения RNA-seq.

1. Подготовка референса

Проще всего анализировать уровень экспрессии с помощью секвенирования RNA, если существует хорошо аннотированный референсный геном. В этом случае нам подойдет использованная ранее девятая хромосома человека. По традиции, проиндексируем референс:
hisat2-build -p 10 chr9.fna chr9
На вход подаем наш fasta-файл с хромосомой, указав число процессоров с помощью -p, а в итоге получили 8 .ht2 бинарных файлов.

2. Оценка качества чтений

Для анализа будем использовать одноконцевые чтения, полученные с помощью набора реагентов для выделения мРНК, без указания информации о цепи.
a. Проанализируем качество исходных чтений с помощью программы fastqc:
fastqc SRR2015719_1.fastq.gz - отчёт
b. В итоге получилось 34 642 046 чтений
c.

d. Качество чтений хорошее, однако помимо закономерного падения качества к концу, есть некоторое падение качества в начале.
e.

f. Все чтения имеют длину 101.
Меня смутило наличие адаптеров (в графе Overrepresented sequences), которых около 20%, и которые, вообще говоря, могут подпортить картину. Удалим их с помощью trimmomatic:
java -jar /usr/share/java/trimmomatic.jar SE -threads 20 SRR2015719_1.fastq.gz trimadapt.fastq.gz ILLUMINACLIP:adapter.fasta:2:30:10 MINLEN:99 2> trimadapt.err
fastqc trimadapt.fastq.gz (файл)
Теперь дела обстоят лучше.

3. Картирование чтений на референс

a. Картируем чтения на хромосому:
hisat2 -p 10 -x chr9 -k 3 -U trimadapt.fastq.gz > rna.sam 2> rna_log.txt (лог)
С помощью -p задаём число ядер, -x - задаём референс, -k - пороговое значения веса выравнивания рида на референс, -U задаёт файл с ридами (если их несколько, нужно подавать файл в формате csv).
b. Из лог файла видно, что 1 808 375 (5.26%) чтений выровнялись ровно 1 раз, а 45 084 (0.13%) - больше 1 раза;
c. Чтобы получить сортированный по самой левой координате .bam файл, переконвертируем его с помощью команды:
samtools sort -O bam -@ 10 rna.sam -o rna.bam
Проиндексируем его:
samtools index -b -@ 10 rna.bam rna.bai
И посмотрим статистику маппинга:
samtools flagstat -@ 10 rna.bam > rna_flag.txt (файл)
Изначальный sam-файл весил 8,2 Гб, полученный bam - 2,4 Гб. Закартировалось 5,53% всех ридов.
d. Отберём чтения, лёгшие только на нашу хромосому, поскольку в файле с ридами тотальная mRNA.
samtools faidx chr9.fna -o chr9.fai
samtools view -h rna.bam NC_000009.12 > rna.chr9.sam
samtools view -bS rna.chr9.sam > rna.chr9.bam
Посмотрим статистику:
samtools flagstat -@ 10 rna.chr9.bam > rna_chr9_flag.txt (файл)
Теперь все 100% ридов принадлежат нашей хромосоме. Проиндексируем файл:
samtools index -b -@ 10 rna.chr9.bam rna_chr9.bai

4. Поиск экспрессирующихся генов

a. Изучим gtf-файл с разметкой генов нашей хромосомы с помощью head gencode.chr9.gtf:

##description: evidence-based annotation of the human genome (GRCh38), version 35 (Ensembl 101)
##provider: GENCODE
##contact: gencode-help@ebi.ac.uk
##format: gtf
##date: 2020-06-03
NC_000009.12	HAVANA	gene	12134	13783	.	+	.	gene_id "ENSG00000236875.3"; gene_type "unprocessed_pseudogene"; gene_name "DDX11L5"; level 2; hgnc_id "HGNC:37106"; havana_gene "OTTHUMG00000019419.2";
NC_000009.12	HAVANA	transcript	12134	13783	.	+	.	gene_id "ENSG00000236875.3"; transcript_id "ENST00000421620.2"; gene_type "unprocessed_pseudogene"; gene_name "DDX11L5"; transcript_type "unprocessed_pseudogene"; transcript_name "DDX11L5-201"; level 2; transcript_support_level "NA"; hgnc_id "HGNC:37106"; ont "PGO:0000005"; tag "basic"; havana_gene "OTTHUMG00000019419.2"; havana_transcript "OTTHUMT00000051447.2";
NC_000009.12	HAVANA	exon	12134	12190	.	+	.	gene_id "ENSG00000236875.3"; transcript_id "ENST00000421620.2"; gene_type "unprocessed_pseudogene"; gene_name "DDX11L5"; transcript_type "unprocessed_pseudogene"; transcript_name "DDX11L5-201"; exon_number 1; exon_id "ENSE00001680583.2"; level 2; transcript_support_level "NA"; hgnc_id "HGNC:37106"; ont "PGO:0000005"; tag "basic"; havana_gene "OTTHUMG00000019419.2"; havana_transcript "OTTHUMT00000051447.2";
NC_000009.12	HAVANA	exon	12291	12340	.	+	.	gene_id "ENSG00000236875.3"; transcript_id "ENST00000421620.2"; gene_type "unprocessed_pseudogene"; gene_name "DDX11L5"; transcript_type "unprocessed_pseudogene"; transcript_name "DDX11L5-201"; exon_number 2; exon_id "ENSE00001759901.2"; level 2; transcript_support_level "NA"; hgnc_id "HGNC:37106"; ont "PGO:0000005"; tag "basic"; havana_gene "OTTHUMG00000019419.2"; havana_transcript "OTTHUMT00000051447.2";
NC_000009.12	HAVANA	exon	12726	12834	.	+	.	gene_id "ENSG00000236875.3"; transcript_id "ENST00000421620.2"; gene_type "unprocessed_pseudogene"; gene_name "DDX11L5"; transcript_type "unprocessed_pseudogene"; transcript_name "DDX11L5-201"; exon_number 3; exon_id "ENSE00002213115.2"; level 2; transcript_support_level "NA"; hgnc_id "HGNC:37106"; ont "PGO:0000005"; tag "basic"; havana_gene "OTTHUMG00000019419.2"; havana_transcript "OTTHUMT00000051447.2";

Файл формата .gtf состоит из шапки и тела. Шапка файла содержит информации о версии, базе данных и дате публикации. Тело файла состоит из таблицы со следующими колонками:

Таблица 1. Описание первых 9 столбцов тела файла в формате GTF

№	Название	Смысл
1	seqname	Название последовательности
2	source	База данных - источник информации
3	feature	Особенности гена
4	start	Координата начала гена
5	end	Координата конца гена
6	score	Вес элемента
7	strand	Направление элемента относительно цепи
8	frame	Рамка считывания (какому основанию участка соответствует первое основание кодона рамки: первому (0), второму (1) или третьему (2))
9	attribute	Дополнительная информация

b. Попробуем выяснить, сколько на нашей хромосоме аннотировано генов:
awk '{$3 == "gene"}' gencode.chr9.gtf | wc -l - всего их 2330 (из 67 924 строк таблицы)
c. Для каждого гена из разметки посчитаем число картированных на этот ген чтений с помощью программы htseq-count:
htseq-count -f bam -s no -m union -t exon rna.chr9.bam gencode.chr9.gtf -o rna_chr9.sam 1> htseq_log_cons.txt 2> htseq_log_err.txt (файл)
Параметры: -f - формат; -s no - риды не имеют информацию о цепи; -m - режим (а данном случае выдаёт количество ридов, которые легли на каждый ген); -t - из gtf-файла берутся только экзоны (используется при стандартном RNA-seq);
d. Полученим чтения, попавшие в гены:
wc -l htseq_log_cons.txt | head -n $lines-5 htseq_log_cons.txt | awk '{s+=$2}END{print s}'
Имеем 1 259 668 чтений, попавших в гены. Непопавшие в границы генов риды указаны в конце файла в строке "__no_feature". Посмотрим конец файла: tail htseq_log_cons.txt - их 398 233.