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

Часть 1. Подготовка чтений

В данном практикуме мы работали с чтениями, картирующимися на участок хромосомы человека (получены путем секвенирования РНК).
Файл chr10.1.fastq с одноконцевыми чтениями в формате fastq был взят мной из директории /P/y14/term3/block4/SNP/rnaseq_reads.
Подготовка чтений осуществлялась в несколько этапов:
1. Анализ качества прочтений
2. Очистка чтений с последующим анализом и сравнением параметров прочтений до и после очистки

Пункт 1. Анализ качества прочтений

Анализ прочтений был осуществлен с помощью прогрммы FastQC, установленной на kodomo. Результатом работы программы является zip архив и html файл, содержащий отчет о качестве прочтений.
Программа была запущена с помощью команды:
fastqc chr10.1.fastq

Пункт 2. Очистка чтений с последующим анализом и сравнением параметров прочтений до и после очистки

Чтения не нуждаются в чистке, так как имеют высокое качество и все попадпют в зеленую область.
Синяя линия - среднее качество чтений, центральные красные линии - медианы, желтые прямоугольники - интерквартальные размахи (разница между верхн. и нижн. квартилями; диапазон значений качества, при котором качество 25% чтений на данной позиции выше нижней границы, а 75% - не выше верхней).
Поле графика разделено на 3 полосы - зеленую, желтую и красную; попадание в данные полосы вышеперечисленных элементов графика позволяет сделать вывод о качестве чтений. После чистки все риды располагаются в зеленой области, следовательно, с помощью чистки мы получили надежные прочтения. Также после чистки стало меньше ридов, так как были удалены риды длинной меньше 50.

Контроль качества чтений

Часть 2. Картирование чтений

Чтения были откартированы с помощью программы hisat2. Программа была импортирована с помощью команды:
export PATH=${PATH}:/home/students/y06/anastaisha_w/hisat2-2.0.5

Таблица 1. Команды, использованные при картировании последовательности.

Команда	Функция	Выдача
hisat2-build chr10.fasta chr10.1	Индексирование референсной последовательности	Индексированный файл chr10.1.fasta
hisat2 -x chr10.1 -U chr10.1_out.fastq --no-softclip >chr10.1.sam	Выравнивание чтений после чистки с референсной последовательностью	Файл, который содержит выравнивание формата SAM chr10.1.sam

Программа картирования была запущена с параметрами -x (путь к индексу) -U (путь к чтениям) --no-softclip (кртирование без подрезания чтений).
Далее необходимо было проанализировать, полученное выравнивание. Для этого я использовала программу Samtools. Она работает с файлами в формате SAM.

Мы не использовали параметр запрещающий без разрывов картировать, так как в данном случае у нас данные транскриптомного анализа, то есть мы работаем с РНК. По сравнению с ДНК с ней могли произойти модификации связанные с перегруппировкой (сплайсинг).
Далее необходимо было проанализировать, полученное выравнивание. Для этого я использовала программу Samtools. Она работает с файлами в формате SAM.

Таблица 2. Команды, использованные для анализа последовательностей в формате sam.

Команда	Функция	Выдача
samtools view chr10.1.sam -bo chr10.1.bam	Программа переводит файл в формат bam	chr10.1.bam
samtools sort chr10.1.bam -T nexact.txt -o chr10.1_sort.bam	Сортировка выравнивания чтений и референса по координате в референсе	chr10.1_sort.bam
samtools index chr10.1_sort.bam	Индексирование отсортированного выравнивания	chr10.1_sort.bam
samtools idxstats chr10.1_sort.bam > result.1.txt	Запись числа откартировавшихся чтений	result.1.txt

Выяснилось, что на хромосому откартировалось 15108 ридов. 205 ридов не были откартированы.

Часть 3. Подсчет чтений.

Команда:
htseq-count -f 'bam' -s 'no' -i 'gene_id' -m 'union' chr10.1_sort.bam /P/y14/term3/block4/SNP/rnaseq_reads/gencode.v19.chr_patch_hapl_scaff.annotation.gtf >> counter_out

И далее: grep -vw 0 counter_out

Выдача:

ENSG00000165732.8       14501
ENSG00000266122.1       3
__no_feature    604
__not_aligned   205

Опции:
-f - Формат файла выравнивания: .bam или .sam
-s - Направление цепи: прямое, обратное, нет направления
-i - Атрибут GFF как feature ID
-m - Режим обработки считывания, перекрывающий более чем одну особенность

14501 чтений легли в границы гена ENSG00000165732.8 и 3 чтения легли в границы гена ENSG00000266122.1, 604 не легли.
Мне достался ген ENSG00000165732.8 (DDX21) и ген ENSG00000266122.1 (RN7SL373P). Описание: DEAD (Asp-Glu-Ala-Asp) box helicase 21 и Misc RNA.