| Команда | Что делает |
|---|---|
| fastqc chr2.1.fastq, fastqc chr2.1_2.fastq | Анализ программой FastQC.Получаем отчет в виде html файла |
| /nfs/srv/databases/ngs/denwer02$ java -jar /nfs/srv/databases/ngs/suvorova/trimmomatic/trimmomatic-0.30.jar SE -phred33 chr2.1.fastq chr2.1_2.fastq TRAILING:20 MINLEN:50 | Очистка при помощи Trimmomatic. Чтения одноконцевые. Необходимо обрезать с конца нуклеотды с качеством ниже 20 и убрать все последовательности длиной менее 50. |
| hisat2 -x chr2_denwer02 -U chr2.1_2.fastq -S chr2_align.sam --no-softclip | Построение выравнивания референсной последовательности и прочтения в формате sam. Параметр --no-spliced-alignment убран (работаем с последовательностью РНКтраскриптов, вырезаются интроны, поэтому они могут картироваться на геном с разрывами).-X - индексированная референсная последовательность, после -U указываем на файл с прочтениями, после S - название выходного файла. |
| samtools view -b chr2_align.sam -o chr2_align.bam | Перевод в bam формат |
| samtools sort chr2_align.bam chr2.1_sorted | Сортировка выравниваний по координате в референсе |
| samtools index chr2.1_sorted.bam | Индексирование |
| htseq-count -f bam -s no -i gene_id -m union chr2.1_sorted.bam gencode.v19.chr_patch_hapl_scaff.annotation.gtf > htseqc_results | Команда для подсчета чтений. Опция -f - формат файла с выравниванием (bam, sam). По умолчанию sam. Опция -s указывает на цепь, по которой были выравнены риды. Из трех возможных вариантов (yes, reverse, no) был выбран вариант no, так как нам точно неизвестно, по какой из цепей строилось выравнивание. -i - GFF атрибут, используемый как feature-ID (нам подходит вариант по умолчанию). Последняя опция -m определяет режим работы команды для неоднозначных выравниваний. |
| grep -w 0 -v htseqc_results > allres.txt | Опция -w 0, чтобы выбрать строки с нулем, ограниченным пробелами, и -v, чтобы вывести все остальные. |
Рис. 1. Оценка качества чтений до триммирования
Рис. 2. Оценка качества чтений после триммирования
Файл chr2.1.fastq содержал 12507 ридов, длина которых составляла от 28 до 51 нуклеотидов, а среднее качество прочтения - около 39%. После очистки чтений их количество практически не уменьшилось: вырезалось 108 ридов (0,86% от всех), а длина стала колебаться от 50 до 51 нуклеотида. На рисунках 1 и 2 можно видеть вывод программы FastQC. По ним видно, что данное триммирование не было необходимым, так как графики практически не отличаются.
Индексные файлы были просто скопированы в директорию, так как в этом практикуме работа идет с той же хромосомой. Потом было произведено картирование чтений. По выводу программы Hisat2:
12399 (100.00%) were unpaired; of these:
54 (0.44%) aligned 0 times
12345 (99.56%) aligned exactly 1 time
0 (0.00%) aligned >1 times
99.56% overall alignment rate
видно, что 99,56% чтений выравнились ровно один раз, всего 54 рида не выравнились ни разу. Можно сказать, что качество картирования высокое.
Программа htseq-count выдает файл со строками, в которых указаны названия генов и количество прочтений, которые легли на этот ген. Была использована команда grep, которая записала в файл allres.txt только строчки, в которых не встречается 0. Как видно:
ENSG00000115053.11 12058
ENSG00000202400.1 5
ENSG00000206885.1 12
ENSG00000233538.1 2
__no_feature 268
__not_aligned 54
большинство ридов легли на один ген. Есть 54 прочтения, которые на предыдущих шагах не картировались на хромосому ни разу, и 268 прочтений, для которых границы генов не определились. Ген ENSG00000115053.11 имеет 14 экзонов и кодирует фосфопротеин, который участвует в синтезе и созревании рибосом и расположен в ядрышке. Ген ENSG00000202400.1, на который легло пять прочтений, кодирует белок SNORD82 - малую ядрышковую РНК¬ как и ген ENSG00000206885.1 - продукт (SNORA75). Транскрипт последнего гена относится к некодирующим РНК.2