В данном практикуме необходимо проанализировать одноконцевые чтения RNA-seq.
1. Подготовка референса
Для начала референс, 8-я хромосома человека, был проиндексирован командой:
hisat2-build -p 10 chr8.fna chr8
Подавая на вход fasta-файл, получили на выход 8 бинарных файлов с расширением .ht2.
2. Оценка качества чтений
Я анализировал одноконцевые чтнеия, которые были мне выданы (SRR2015717). При помощи команды:
fastqc SRR2015717_1.fasta.gz
получили страничку (ссылка). Всего получилось 31803919 чтений.
Качество чтений достаточно неплохое, но почему-то есть провал по качеству не только в конце, но и в начале.
Как видно, все чтения имеют длину 101 (так нет таких чтений, которые бы имели длину 100 или 102)
3. Картирование чтений на референс
Откартируем наши чтения на хромосому:
hisat2 -p 10 -x chr8 -k 3 -U SRR2015719_1.fastq.gz > rna.sam 2> log.txt (файл)
Здесь -p задаёт число ядер, -х - референс, -k - пороговое значение веса выравнивания рида на референс, -U задаёт файл с ридами. Из файла: 1180890 (3.71%) aligned exactly 1 time и 21718 (0.07%) aligned >1 times (простите, надеюсь, что можно просто вставлять английский текст, так как я так больше не могу).
Теперь отсортируем 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
Как видно из файла, закартировалось 3.87% ридов
Теперь надо отбрать только те чтения, что легли на нашу хромосому (в файле с ридами тотальная мРНК). Команды:
samtools faidx chr8.fna -o chr8.fai
samtools view -h rna.bam NC_000008.11 > rna.chr8.sam
samtools view -bS rna.chr8.sam > rna.chr8.bam
samtools flagstat -@ 10 rna.chr8.bam > rna.chr8_flag.txt
samtools index -b -@ 10 rna.chr8.bam rna.chr8.bai
Из файла видим, что теперь все чтения принадлежат нашей хромосоме, значит, всё хорошо.
4. Поиск экспрессирующихся генов
При помощи команды head изучим файл с разметкой.
##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_000008.11 HAVANA gene 64091 64320 . - . gene_id "ENSG00000253620.2"; gene_type "processed_pseudogene"; gene_name "AC144568.1"; level 2; havana_gene "OTTHUMG00000163938.2"; NC_000008.11 HAVANA transcript 64091 64320 . - . gene_id "ENSG00000253620.2"; transcript_id "ENST00000520139.2"; gene_type "processed_pseudogene"; gene_name "AC144568.1"; transcript_type "processed_pseudogene"; transcript_name "AC144568.1-201"; level 2; transcript_support_level "NA"; ont "PGO:0000004"; tag "basic"; havana_gene "OTTHUMG00000163938.2"; havana_transcript "OTTHUMT00000376500.2"; NC_000008.11 HAVANA exon 64269 64320 . - . gene_id "ENSG00000253620.2"; transcript_id "ENST00000520139.2"; gene_type "processed_pseudogene"; gene_name "AC144568.1"; transcript_type "processed_pseudogene"; transcript_name "AC144568.1-201"; exon_number 1; exon_id "ENSE00002132061.2"; level 2; transcript_support_level "NA"; ont "PGO:0000004"; tag "basic"; havana_gene "OTTHUMG00000163938.2"; havana_transcript "OTTHUMT00000376500.2"; NC_000008.11 HAVANA exon 64091 64175 . - . gene_id "ENSG00000253620.2"; transcript_id "ENST00000520139.2"; gene_type "processed_pseudogene"; gene_name "AC144568.1"; transcript_type "processed_pseudogene"; transcript_name "AC144568.1-201"; exon_number 2; exon_id "ENSE00002126516.2"; level 2; transcript_support_level "NA"; ont "PGO:0000004"; tag "basic"; havana_gene "OTTHUMG00000163938.2"; havana_transcript "OTTHUMT00000376500.2"; NC_000008.11 HAVANA gene 72601 79775 . + . gene_id "ENSG00000253896.3"; gene_type "lncRNA"; gene_name "AC144568.2"; level 2; havana_gene "OTTHUMG00000163939.3";
Gene Transfer Format (GTF) - одна из трёх версий GFF. Каждая строка в файле формата GFF содержит 9 колонок, разделенных знаком табуляции. Каждая колонка называется полем и имеет своё назначение. Список названий полей и их содержание приведены ниже.
| № | Название поля | Описание (спасибо, вики) |
|---|---|---|
| 1 | seqid | Название (идентификатор) последовательности, где находится данный элемент. |
| 2 | source | Источник определения элемента |
| 3 | type | Тип элемента |
| 4 и 5 | start | Начальные и конечные координаты элемента. |
| 6 | score | Вес элемента |
| 7 | strand | Направление элемента относительно цепи, на которой располагается. + 5'-->3' и - 3'-->5' |
| 8 | frame | Рамка считывания или фаза для белок-кодирующих последовательностей. Указывает какому основанию участка соответствует первое основание кодона раки |
| 9 | attribute | Поле для дополнительной информации |
При помощи htseq-count посчитаем для каждого гена число картированных на этот ген чтений:
htseq-count -f bam -s no -m union -t exon rna.chr8.bam gencode.chr8.gtf -o rna.chr8.sam 1> log_cons.txt 2> log_err.txt (файл)
А теперь поподробнее о команде. -а - задаёт формат (bam), -s no - чтение считается перекрывающимся с функцией, независимо от того, сопоставлено ли оно с той же или противоположной цепью, что и объект, -m - режим для обработки операций чтения, перекрывающих более одной функции, -t - здесь (exon) задаёт, чтобы из файла брались только экзоны.
Осталось только получить чтения, попавшие в гены. Команда:
wc -l log_cons.txt | head -n $lines-5 log_cons.txt | awk '{s+=$2}END{print s}'
Получили 840337 чтений, попавших в гены. Чтобы узнать, сколько чтений не попали в границы генов, надо прочитать конец файла: tail log_cons.txt.
Получаем 288871 ридов.