Анализ данных RNA-seq

Задание 1: Подготовка референса

Собственно, подготовка референса заключается в том, чтобы его проиндексировать:

$ hisat2-build chr3.fna chr3 2> index.log

Собственно, на вход программе hisat2-build я подал, собственно, референс, а на выходе получил 8 бинарных файлов, каждый из которых, по всей видимости, смотря на логи, представляет собой проиндексированную часть референса.

Задание 2: Оценка качества чтений

Провёл я оценку с помощью программы fastqc:

$ fastqc RNA-seq_data/SRR2015719_1.fastq.gz

Количество чтений: 34642046;

В принципе у нас качество прочтения всех нуклеотидов находится значительно выше красной зоны (а большинство из них находится в зелёной зоне), это определённо идёт в плюс нашим чтениям. Но при этом стоит отметить, что 2-11 нуклеотиды не имеют бокса в боксплоте, остались только усы. Из этого, на мой взгляд, можно сделать несколько выводов: либо плотность качества по боксплоту распределена равномерно и поэтому там не удаётся выделить квартили, либо наоборот плотность качества прочтения сконцентрированна в одной точке, что выделять квартили бессмысленно, либо просто очень мало прочтений этих нуклеотидов, либо качества прочтения нуклеотидов слишком разрозенны, чтобы выделять среди них статистические характеристики. Также ярко видно, что посередине условного "рида" качество прочтения определённо лучше, чем с концов. Посмотрим, исправит ли это trimmomatic. В целом качества прочтения ридов можно охарактеризовать, как очень хорошие, с этими ридами можно работать.

Судя по приведённому выше графику, все наши риды имеют длину в 101 нуклеотид.

Задание 3: Картирование чтений на референс

Собственно, само картирование я буду делать с помощью программы hisat2:

$ hisat2 -p 12 -x chr3 -U SRR2015719_1.fastq.gz > RNA-seq_cartir.sam 2> cartir.log

• hisat2 - обращение к программе hisat2;
• -p 12 - уточнение, какое количество ядер сервера стоит использовать (для ускорения процесса);
• -x chr3 - подача нашего индексированного референса;
• -U SRR2015719_1.fastq.gz - подача наших одноконцевых чтений;
• > RNA-seq_cartir.sam - перенаправление всего нашего выхода от картирования в файл;
• 2> cartir.log - сохраняем логи в файл.

Исходя из логов, всего закартировалось 6.41% чтений, что логично, так как мы тотальную РНК картируем всего лишь на одну хромосому, а также РНК может сильно отличаться от генов из-за сплайсинга; всё это внесло свою лепту в столь низкую долю картированных чтений. Но при этом стоит отметить, что лишь 0.1% закартировались более 1-го раза, все остальные риды закартировались ровно один раз, что, на мой взгляд, должно говорить об однозначности картирования, и не может не радовать (меня).

Далее надо перевести sam в bam и отсортировать. Сделаю я это уже известными нам по прошлым практикумам командами:

$ samtools sort -@ 12 -o RNA-seq_cartir.bam RNA-seq_cartir.sam 2> samtoolsort.log
    
$ samtools index -@ 12 RNA-seq_cartir.bam

Далее надо отобрать чтения, которые картировались исключительно на мою хромосому. Для этого я последовательно воспользуюсь двумя командами:

$ samtools view -h RNA-seq_cartir.bam NC_000003.12 > RNA-seq_cartir_chr3.sam
    
$ samtools view -b RNA-seq_cartir_chr3.sam > RNA-seq_cartir_chr3.bam

Задание 4: Поиск экспрессирующихся генов

Внутри файл с разметкой генов выглядит таким образом:

Первые несколько строчек (шапка) содержат описание, собственно файла: формат, дата создания, контакты, провайдер и т. д. Само тело файла своим видом очень напоминает .sam, но если вглядеться, то сразу понимаешь, что это не он. Тут в каждой строчке содержится аннотация одного гена, экзона, транскрипта, CDS, UTR (untranslated regions) и др.

Далее для каждого гена из разметки я посчитал количество картированных на него чтений:

$ htseq-count -f bam -s no -t exon RNA-seq_cartir_chr3.bam gencode.chr3.gtf > htseq-count.txt

• htseq-count - обращение к программе htseq-count;
• -f bam - таким образом мы уточняем, что наши картированные чтения будут подаваться программе именно в формате .bam;
• -s no - этим параметром мы указываем, что у нас информация берётся не из цепи, с ней вообще мало что связано;
• -t exon - этим параметром мы укзываем, что из gtf-файла мы будем брать только экзонные строчки;
• RNA-seq_cartir_chr3.bam - подача в программу наших картированных на одну хромосому чтений;
• gencode.chr3.gtf - наш файл с разметкой генов;
• > htseq-count.txt - перенаправление выхода программы в отдельный файл.

Далее, полученный с помощью прошлой команды файл я обработал с помощью скрипта, который подсчитал мне, сколько всего чтений попало в границы генов. Получилось, кстати, 1541257. После этого я просто посмотрел в конец того же самого файла и увидел, что чтений, которые попали мимо границ генов (__no_feature), оказалось 527098.