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

Создал папку для работы программы hisat2. Назвал её hisat2. Скопировал туда хромосому. Проиндексировал его командой:
hisat2-build chr7.fna chr7 2> log_hisat_index.txt
На вход подаётся рефересная последовательность. На выходе были получены следующие файлы:
chr7.1.ht2 chr7.2.ht2 chr7.3.ht2 chr7.4.ht2 chr7.5.ht2 chr7.6.ht2 chr7.7.ht2 chr7.8.ht2 log_hisat_index.txt
Это бинарные файлы индексирования, не считая последний - туда записывались ошибки.

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

Получим html-файлы о качестве чтений с помощью команды:
fastqc SRR2015717_1.fastq.gz
html страница полученная при работе команды
Количество чтений - 31803919. Качество чтений полность попадаёт в зеленую область, что говорит о его высоком качестве:
Картинку из раздела Sequence Length Distribution:

Длина чтений - 101 нуклетид.

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

Для картирования чтений на референс воспользуемся командой:
hisat2 -p 10 -x chr7 -k 3 -U SRR2015717_1.fastq.gz > rna_mapped.sam 2>log_er.txt
-p 10 - использовал 10 ядер
-x chr7 - референс
-k 3 - максимальное количество выравниваний (в данном случае 3), чей score больше или равен score любого другого выравнивания.
-U файл с ридами
2>log_er.txt - сюда перенапрявлялись ошибки
Относительно малая часть чтений была картированна на референсный геном - 4.64%.
Переведём sam файл в bam-формат с помощью уже известных нам команд:
samtools sort -o rna_mapped.bam rna_mapped.sam
samtools index -@ 12 rna_mapped.bamБ
Отберём чтения, которые легли на 7-ю хромосому:
samtools view -@ 10 -h rna_mapped.bam NC_000007.14 > rna_chr7.sam
samtools view -@ 10 -b rna_chr7.sam > rna_chr7.bam

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

Скопирую себе файл gencode.chr7.gtf. Он состоит из шапки и тела. В голове прописана разметка (версия), информация о базе данных, формат файла и дата публикации.
Таблица, расположенная в теле файла состоит из следующих колонок:
seqname - название последоватлеьности, где аннотирован ген
source - источник аннтации
feature - особенности гена
start & and - начало и конец гена
score - достоверность последовательности
strand - направление последовательности
frame - рамка считывания
attribute - дополнительная (но при этом важная) информация.
Для каждого гена из разметки посчитаем число картированных на этот ген чтений с помощью следующей команды:
htseq-count -f bam -s no -t exon rna_chr7.bam gencode.chr7.gtf 1> htseq
-f bam - формат входного файла
-s no - сообщаем программе, что у нас нет информации о цепи
-t exon - выбираем данные, которые мы используем из gtf файла, при использовании "exon" - информация об экзонах.
Используя команду:
wc -l htseq | head -n $linea-5 htseq | awk '{s+=$2}END{print s}'
Мы обнаружили, что в груницы генов попало 1056500 чтений. В свою очередь, 299186 чтений не попали в границы чтений.