Ралдугина Василиса

Студентка Факультета биоинженерии и биоинформатики

МГУ имени М.В. Ломоносова

Обо мне

Главная

Сайт ФББ МГУ

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


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


Была создана директория /nfs/srv/databases/ngs/vasidze/pr12 в которой выполнялись все задания практикума. Туда были скопированы файл с ридами (chr9_2.fastq) и с последовательностью 9 хромосомы (chr9.fasta).
Далее с помощью программы FastQC было установлено качество прочтений. Результатом работы программы был html файл и zip архив.
Команда запуска программы : fastqc chr9.1.fastq
Далее, исполюзуя программу Trimmomatiс были оставлены только чтения с длинной не менее 50 нуклеотидов.
Команда запуска программы : java -jar /usr/share/java/trimmomatic.jar SE -phred33 chr9.1.fastq out_chr9.1.fastq MINLEN:50
До чистки было 19976 чтений, после чистки стало 19267 чтений.
На рисунке 1 изображен контроль качества чтений до чистки. Синяя линия на графике - среднее качество чтений, центральные красные линии - медиана, желтые прямоугольники - интерквартальный размах (разница между верхней и нижней квартилями, диапазон значений качества, при котором качество 25% чтений на данной позиции выше нижней границы, а 75% - не выше верхней). Поле графика разделено на 3 полосы зеленого, желтого и красного цветов, попадание в которые вышеперечисленных элементов графика позволяет судить о качестве чтений.
На рисунке 2 изображен контроль качества чтений до чистки. До чистки чтения были хорошие, поэтому после нее мы не видим особых изменений. С помощью чистки мы получили самые надежные прочтения длинной меньше 50.

Рис 1. Контроль прочтений до чистки.

Рис 2. Контроль прочтений после чистки.

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


Чтения были откартированы с помощью программы hisat2 аналогично тому как это проводилось в практикуме 11. Напомним как это происходило:
Программа была импортирована с помощью команды : export PATH=${PATH}:/home/students/y06/anastaisha_w/hisat2-2.0.5
Таблицы 1. Команды, использованные при картировании последовательности.
КомандаФункцияВыдача
hisat2-build chr9.fasta chr9Индексирование референсной
последовательности		Индексированный файл chr9.fasta
hisat2 -x chr9 -U out_chr9.1.fastq --no-softclip > chr9.samВыравнивание чтений после чистки
с референсной последовательностью		Файл, который содержит
выравнивание формата SAM 1_align.sam

Мы не использовали параметр запрещающий без разрывов картировать, так как в данном случае у нас данные транскриптомного анализа, то есть мы работаем с РНК. По сравнению с ДНК с ней могли произойти модификации связанные с перегруппировкой (сплайсинг).
Далее необходимо было проанализировать, полученное выравнивание. Для этого я использовала программу Samtools. Она работает с файлами в формате SAM.
Таблица 2. Команды, использованные для анализа последовательностей в формате sam.
КомандаФункцияВыдача
samtools view chr9.sam -bo chr9.bamПрограмма переводит файл в формат bamchr9.bam
samtools sort chr9.bam -T nexact.txt -o chr9_sort.bamСортировка выравнивания
чтений и референса по
координате в референсе		chr9_sort.bam
samtools index chr9_sort.bamИндексирование отсортированного выравниванияchr9_sort.bam
samtools idxstats chr13_sort.bam > result.txtЗапись числа откартировавшихся чтенийresult.txt

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

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


Для подсчета чтений мы использовали пакет Bedtools.
Далее была использована функция bamtobed, которая позволяет из bam файла сделать bed:
Команда запуска: /P/y14/term3/block4/SNP/bedtools2/bin/bedtools bamtobed -i chr9.bam > chr9.bed
Выдача: chr9.bed
Затем нашли пересечение наших выравниваний с разметкой по генам:
Команда запуска: /P/y14/term3/block4/SNP/bedtools2/bin/bedtools intersect -a /P/y14/term3/block4/SNP/rnaseq_reads/gencode.genes.bed -b chr9.bed -c > chr.intersect.bed
Выдача: chr.intersect.bed
После этого для упрощения чтения файла была проведена сортировака:
Команда запуска: sort -k 6 -r us.bed > s.bed
Выдача: s.bed
После проведения подсчета чтений мы видим, что прочтения легли в границы двух генов. Прочтения, которые не попали скорее всего относятся к не матричной РНК.
Мои чтения относятся к GL877870 и CICP10
Ген CICP10 - псевдоген репрессора capicua 10. Длина этого гена 466 пар оснований. Координаты гена: 243,062,007-243,062,293
Ген GL877870. Длина этого гена 66021 пара оснований. Координаты гена: 243,059,660-243,189,373

Часть IV. Функции пакета BEDTOOLS


Задание 1. Получите из файла c выравниванием файл с чтениями в формате fastq
Команда для запуска:/P/y14/term3/block4/SNP/bedtools2/bin/bedtools bamtofastq -i chr9.bam -fq chr2.1.fastq
Выдача: chr2.1.fastq
Задание 4. Объедините Ваши чтения в кластеры.
Команда для запуска: /P/y14/term3/block4/SNP/bedtools2/bin/bedtools cluster -i chr9.bed -s >cluster.bed
Выдача: cluster.bed
Задание 3. Разбейте свою хромосому на фрагменты по 1 млн нуклеотидов. Какова длина хромосомы в нуклеотидах? Сколько в результате получилось интервалов?
Исходные файлы: chr9.txt. (содержание: chr9 141213431 )
Команда для запуска:/P/y14/term3/block4/SNP/bedtools2/bin/bedtools makewindows -g chr9.txt -w 1000000 > win.bed
Выдача: win.bed
Флаг -g - файл с указанием размера участков генома, w- размер промежутка. Длина хромосомы 141 213 431 п.н.. Получилось 142 промежутка.

© Raldugina Vasilisa 2016