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

Дано

• Чтения, картирующиеся на участок хромосомы человека (получены путем секвенирования РНК), две реплики - chr2.1.fastq, chr2.2.fastq

• Хромосомы человеческого генома (сборка версии hg19) - chr2.fasta

Часть I: подготовка чтений

Анализ качества чтений

С помощью программы FastQC был проведён контроль качества чтений для обеих реплик. Команды:

fastqc chr2.1.fastq
fastqc chr2.2.fastq 

Получив на вход файл в формате .fastq, программа выдала архив (.zip) с отчетом по работе в виде текстового файла и рисунков, в котором помимо прочей информации находились файлы chr2.1_fastqc.html и chr2.2_fastqc.html, содержащие визуализацию качества чтений. В Таблице 1 приведены различия по чтениям первой и второй реплик.

На Рис. 1-2 представлена информация, полученная с помощью программы FastQC, а именно - характеристика качества каждого нуклеотида в ридах. По абсциссе расположены позиции в ридах, а по ординате - Quality score, причем данные представлены в виде диаграмм boxplot, где границами желтых столбцов служат первый и третий квартили, красная линия - медиана, синяя линия - среднее значение, а отходящие от столбцов "усы" демонстрируют размах между 10-м и 90-м процентилями. Помимо этого само поле графика покрашено в три цвета: нахождение диаграмм в зеленой плотности говорит о том, что качество данной позиции хорошее и она достоверна (выше 28), в оранжевой - среднее (от 20 до 28), в красной - плохое, таким последовательностям верить нельзя (от 0 до 20).

Можно заметить, что качество последовательностей ридов примерно одинаковое в обеих репликах, все диаграммы входят входят в зелёную область, следовательно, для всех позиций в чтениях quality score > 28, что я является очень хорошим показателем. Помимо этого, в html-файлах содержится ещё много полезной информации, ниже приведен анализ некоторых графиков, выдаваемых программой. На Рис. 3-4 представлена тепловая карта качества последовательностей, поступающих из проточных ячеек в Illumina - такую информацию можно получить только для этого метода секвенирования. Здесь по вертикальной оси обозначены номера ячеек, а по горизонтальной - позиции в риде. Синий цвет характеризует хорошее качество, а градация к красному демонстрирует падение качества позиции.

Как можно заметить, у второй реплики большее количество ячеек с позициями с низким качеством и встречаются даже красные, означающие очень низкое качество позиции в чтении, ячейки. В следствие чего программа указала, что нужно обратить внимание на эти данные для первой реплики, а для второй, что эти данные могут вносить какие-то помехи.
Также программа FastQC выдает частоту встречаемости азотистых оснований для каждой позиции в ридах. В хорошем секвенировании частота встречаемости оснований должны быть постоянными из позиции в позицию. На Рис. 5-6 представлена такая характеристика, можно заметить, что в обеих репликах на протяжении всех позиций наблюдаются скачки частот, возможно, это может быть объяснено наличием оверэкспрессированных участков, в частности 200 таких участков найдены в первой реплике и 179 - во второй.

Часть II: картирование чтений

Картирование чтений

С помощью программы Hisat2 было осуществлено картирование чтений. Сначала была проиндексирована референсная последовательность с помощью команды:

hisat2-build chr2.fasta chr2

В результате чего были получены 8 файлов в формате .ht2. Затем построено выравнивание обеих реплик и референса в формате .sam с помощью команды, где параметр --no-softclip запрещает подрезать чтения:

hisat2 -x chr2 -U chr2.1.fastq  --no-softclip > align.sam    
hisat2 -x chr2 -U chr2.2.fastq  --no-softclip > align2.sam    

Примечание: В отличие от предыдущего практикума при выполненнии этой команды был убран параметр --no-spliced-alignmen, так как анализируется последовательности РНК, подвергающиеся сплайсингу, поэтому в данном случае разрывы можно разрешить.
В результате выдача сохранена в файл align.sam - для первой реплике, align2.sam - для второй реплики.

Анализ выравнивания

С помощью команды, приведенной ниже, выравнивания из формата .sam было переведено в бинарный формат .bam:

samtools view align.sam -bo align1.bam
samtools view align2.sam -bo align2.bam

Для этого использовалась команда view из пакета samtools. Затем выравнивания чтений с референсом было отсортировано по координате в референсе начала чтения с помощью следующей команды:

samtools sort align1.bam -T 1.txt -o chr2.1.bam
samtools sort align2.bam -T 1.txt -o chr2.2.bam

И далее отсортированные файлы chr2.1.bam и chr2.2.bam были проиндексированы командой:

samtools index chr2.1.bam
samtools index chr2.2.bam

После чего получены файлы chr2.1.bam.bai и chr2.2.bam.bai. Затем выяснялось, сколько чтений откартировано на геном, эта информация выдавалась сразу после работы программы Hisat2, в Таблицу 2 сведены эти данные.

12507 reads; of these:
  12507 (100.00%) were unpaired; of these:
    82 (0.66%) aligned 0 times
    12377 (98.96%) aligned exactly 1 time
    48 (0.38%) aligned >1 times
99.34% overall alignment rate

7126 reads; of these:
  7126 (100.00%) were unpaired; of these:
    46 (0.65%) aligned 0 times
    7056 (99.02%) aligned exactly 1 time
    24 (0.34%) aligned >1 times
99.35% overall alignment rate

samtools idxstats chr2.1.bam > map1.txt
samtools idxstats chr2.2.bam > map2.txt

В результате в файлах map1.txt и map2.txt была записана следующая информация: длина референсной последовательности - 243199373 и число откартированных чтений, однако никаких пояснений не даётся.

С помощью программы IGV - Integrative Genomics Viewer, были визуализированы выравнивания обеих реплик. На вход программа получает отсортированный файл с выравниванием в формате .bam. Рис. 7 демонстрируют результат работы этой программы, причем она сразу способна отобразить оба выравнивания относительно одного референса. Верхней линией на рисунках показана картированная хромосома, что позволяет довольно легко найти нужный участок.

Можно заметить, что отличие такого выравнивания от геномных чтений, рассматриваемых в предыдущем практикуме в том, что длина ридов намного меньше длины ридов для геномов. Также видно, что гены намного более полно покрываются ридами, причем чтения в основном покрывают кодирующие участки генов. Полосами показывается то, что рид приходится на кодирующие участки, получающиеся на мРНК в результате вырезания участка, изображенного полосой.

Часть III: Работа с Bedtools

Подсчет чтений

С помощью программы Bedtools был выполнен перевод полученных выравниваний из формата .bam в формат .bed, с которым способна легко работать данная программа, это выполнено для обеих реплик. Команды:

bedtools bamtobed -i chr2.1.bam >2.1.bed
bedtools bamtobed -i chr2.2.bam >2.2.bed

Далее с помощью программы intersect производилось пересечение реплик с известной разметкой генов с помощью команд, где параметр -a отражает разметку, с которой выполняется пересечение, а -b - реплику, -с проводит подсчет количества ридов, перекрывающих ген:

bedtools intersect -a gencode.genes.bed -b 2.1.bed -c > coun2.1.bed
bedtools intersect -a gencode.genes.bed -b 2.2.bed -c > coun2.2.bed

В результате были получены два файла coun2.1.bed и coun2.2.bed, которые затем были отсортированы по 6 колонке, где и содержалась информация о присутствии чтений, покрывающих гены. Команды:

sort -k 6 -r coun2.1.bed > sort2.1.bed
sort -k 6 -r coun2.2.bed > sort2.2.bed

Информация о генах, в которые вошли чтения из обеих реплик, сведена в Таблицу 3.

Изучение опций Bedtools

Был взят файл, полученный выше в формате .bam и переведен в fastq-формат. Команда:

bedtools bamtofastq -i chr2.1.bam -fq chr2.1.fq

2. Получите файл с нуклеотидной последовательностью (.fasta) для одного из покрытых Вашими чтениями генов

Для выполнения данного задания был создан файл 1.bed с координатами одного из генов, который покрывается ридами, информация об этом приведена выше. А затем из полной последовательности хромосомы 2 человека была получена его последовательность и выдача записана в файл gene.fasta в fasta-формате. Команда:

bedtools getfasta -fi chr2.fasta -bed 1.bed > gene.fasta

3. Разбейте свою хромосому на фрагменты по 1 млн нуклеотидов

В данном задании на вход программа получала файл chr2.txt, в котором была записана длина хромосомы, равняющаяся 242193529 п.н., а затем программа разбивала хромосому на интервалы длиной 1000000 п.н. и записывала свою выдачу в файл wind.txt. В результате работы программы получилось 243 интервала. Команда:

bedtools makewindows -g chr2.txt -w 1000000 > wind.txt

4. Объедините Ваши чтения в кластеры

При выполнении данной команды на вход поступает файл с сортированными выравненными чтениями в формате .bed, параметр -d позволяет объединять чтения, находящиеся на расстоянии 5 п.н. и ближе, в результате получен файл с объединением в кластеры cluster.txt. Команда:

bedtools cluster -i 2.1.bed -d 5 > cluster.txt

5. Наберите из Вашей хромосомы 1000 случайных фрагментов по 200 нуклеотидов

В данном задании на вход программа получала файл chr2.txt, в котором была записана длина хромосомы, равняющаяся 242193529 п.н., а затем программа сгенерировала 1000 рандомных последовательностей длиной 200 п.н., а выдача программы записывалась random.txt. Команда:

bedtools random -g chr2.txt -n 1000 -l 200 > random.txt

7. Получите координаты одного из покрытых Вашими чтениями генов, расширенные на 1000 нуклеотидов в обе стороны

Для выполнения данного задания был создан файл 1.bed с координатами одного из генов, который покрывается ридами, информация об этом приведена выше. С помощью программы slop был получен файл slop.txt, в котором записаны координаты гена на 1000 п.н. больше слева и справа от него. Команда:

bedtools slop -i 1.bed -g chr2.txt -b 1000 > slop.txt

8. Получите координаты одного из покрытых Вашими чтениями генов, сдвинутые на 500 нуклеотидов ближе к началу хромосомы

На вход программа получала всё те же файлы с координатами гена и с указанным числом пар нуклеотидов в хромосоме, а её выдача записывалась в файл shift.txt, в котором указывались координаты гена, сдвинутые на 500 п.н. ближе к началу хромосомы. Команда:

bedtools shift -i 1.bed -g chr2.txt -s -500 > shift.txt

9. Получите непересекающиеся фрагменты, соответствующие области, покрытой Вашими чтениями

Если подразумевается получение фрагментов, которые не пересекаются с другими чтениями, то после получения файла cov.txt (см. Пункт 10), была выполнена команда, достающая из этого файла строчки, где пересечение равно 1, и получаем файл с координатами таких участков cov1.txt. Команда:

grep -w 1 cov.txt > cov1.txt

Если же подразумевается получение участков, непересеченных ридами, то с помощью программы complement можно получить координаты участков хромосомы, непересечнной ридами, выдача программы записывается в файл compl.txt:

bedtools complement -i 2.1.bed -g chr2.txt > compl.txt

10. Получите файл с координатами интервалов, покрытых Вашими чтениями, с информацией о покрытии в любом формате

На вход программа получала файл с сортированным выравниванием чтений и длиной хромосомы, а на выходе получется файл с координатами интервалов, покрытых чтениями - cov.txt, где в 4 столбце представлено количество чтений, покрывающий данный промежуток. Команда:

bedtools genomecov -i 2.1.bed -g chr2.txt -bg > cov.txt

Источники: