Анализ транскриптомов.
I часть. Подготовка чтений.
Перед началом работы, как и в практикуме 11, мною в директорию /nfs/srv/databases/ngs/yudna были скопированы фыйлы:
референсный файл chr4.fasta (сборка версии hg19) и chr4.1.fastq, chr4.2.fastq.
1. Проведение анализа качества чтения. Необходимо сделать контроль качества чтения с помощью
программы fastqc. Я использовала программу, установленную на kodomo.
Команда: fastqc chr4.1.fastq и fastqc chr4.2.fastq
Выдача: отчет о программе в виде html файла chr4.1_fastq.html, chr4.2_fastq.html.
Все прочтения лежат в зеленой зоне на графикакх качества прочтений, длина прочтений в соответствующих сборках
32-51 и 39-51 нуклеотид, поэтому было решено нецелесообразно проводить очистку чтений программой trimmomatic.
 
Рис.1 График качества прочтений для chr4.1_fastq.html и chr4.2_fastq.html соответственно.
II часть. Картирование чтений.
3. Картирование чтений. При выполнении практикума 11 уже было выполнено индексирование
референсной последовательности. В таблице эти команды отмечены зеленым цветом.
Таблица 1
| Команда |
Назначение/Выдача |
| export PATH=${PATH}:/home/students/y06/anastaisha_w/hisat2-2.0.5 |
Вызывает программу, лежащую в указанной директории. |
| hisat2-build chr4.fasta chr4 |
Индексирует референснуб последовательность, выдает множество файло с расширением .ht2 |
hisat2 -x chr4 -U chr4.1.fastq --no-softclip > align1.sam;
hisat2 -x chr4 -U chr4.2.fastq --no-softclip > align2.sam |
Строит выранивание прочтения и референса, сохраняет результаты в отдельный файл align1.sam, align2.sam |
В отличие от практикума 11 был убран параметр --no-spliced-alignment при построении выравниваний прочтений.
Так следует сделать, так как мы анализирыум последовательности РНК, уже прошедшие вырезания интронов, поэтому относительно последовательности генома в них могут быть разрывы.
4. Анализ выравнивания. Для анализа выравнивания использовалась программа sаmtools. Результаты ее работы
сведены в таблицу.
Таблица 2
| Команда |
Назначение/Выдача |
samtools view align1.sam -bo align1.bam
samtools view align2.sam -bo align2.bam |
Переводит выравние с референсом в бинарный формат |
samtools sort align1.bam -T myfile.txt -o sorted1.bam
samtools sort align2.bam -T myfile.txt -o sorted2.bam |
Сортирует выравнивание чтений с референсом по координате в референсе начала чтения. |
samtools index sorted1.bam
samtools index sorted2.bam |
Индексирует отсортированный файл .bam |
Количесво откартированных на референсную последовательность (программa Hisat2 выдает на stdout по окончании работы).
Таблица 3
| chr4.1 |
chr4.2 |
72(2,63%) - откартированы 0 раз
2658(97,18) - ровно 1 раз
5(0,18) - более 1 раза |
33(1,45%) - откартированы 0 раз
2236(98,50) - ровно 1 раз
1(0,04) - более 1 раза |
III часть. Подсчет прочтений.
Для подсчета глубины прочтений и выяснения их координат была использована программа bedtools .
Программе задавлись следующие команды:
Таблица 4
| Команда |
Назначение/Выдача |
| export PATH=${PATH}:/P/y14/term3/block4/SNP/bedtools2/bin |
Вызывает программу bedtools из указанной директории. |
| bedtools bamtobed -i sorted1.bam > b.bed |
Переводит файл формата .bam в формат .bed и записывает в файл b.bed. |
| bedtools intersect -a gencode.genes.bed -b b.bed -u > b2.bed |
Пересекает общий файл с разметкой генов с чтениями. Параметр -u оставляет только те чтения, которые хотя бы раз пересеклись с генами из общего файла. |
ПРИМЕЧАНИЕ. Праматетр -с у команды intersect сразу считает количество чтений, отложившихся на какой-либо ген (выдает отдельный столбик значений, большая часть из которых 0 и только для 2-4 генов будут совпадения).
Однако в файле gencode.genes.bed собраны гены со всех хромосом, поэтому такой поиск перекрываний чтений только для одной хромосомы не удобен.
Для подсчета количества чтений я использовала команду подсчета количества строк в файле, содержащих заданное слово: grep (название гена) b2.bed -c, где b2.bed - файл, в который собраны все, хотя бы раз пересекающиеся, прочтения.
По итогам анализа с помощью некоторых подпрограмм программы bedtools были получены следующие результаты:
| Хромосома |
Кординаты |
Ген |
Функция |
Направление цепи |
Кол-во чтений |
| Chr4 |
10069713:10074643 |
RP11-448G15.3 |
Кодирует фермент - глицерол-3-фосфатацетилтрансферазу, участвующий в регуляции липидного обмена [3]. |
+ |
3 |
| Chr4 |
10075963:10118573 |
WDR1 |
Кодирует белок, состоящий из 9 доменов, осуществляющих белок-белковые взаимодействия. Участвует в разборке актиновых филаментов в клетке.[1] |
+ |
226 |
| Chr4 |
10080235:10080316 |
MIR3138 |
МикроРНК 3138 [2] |
+ |
3 |
| Chr4 |
10117380:10117508 |
RNA5SP155 |
Рибосомальная 5S РНК [2] |
+ |
3 |
IV часть. Работа с bedtools.
Для освоения программы bedtools были опробованы следующие команды:
#1
Команда:bedtools bamtofastq -i sorted1.bam -fq sorted1.fq
Назначение:Переводит файл формата .bam с выравниванием в файл формата .fastq.
#2
1)Команда: bedtools intersect -a gencode.genes.bed -b sorted1.bed -u | head -1 > test.bed
Назначение:Извлекли 1 строчку, содержащую координаты 1-ого гена, покрытого чтениями.
2)Команда:bedtools getfasta -fi chr4.fasta -bed test.bed
Назначение:Из файла, содержащего целиком нуклеотидную последовательность хромосомы, извлекли последовательность нужного нам гена.
#4
Команда:bedtools cluster -i b.bed -d 50 > clusted.bed
Назначение:объединяет полученные прочтения в кластеры. По умолчанию в клатеры объединяются перекрывающиеся прочтения, параметр -d позволяет объединять прочтения, находящиеся в моем случае на расстоянии 50 нуклеотидов.
#7
Команда:bedtools slop -i test.bed -g chr4.genome -b 1000
Назначение:удлиняет координаты указанного в файле test.bed гена на 1000 нуклеотидов с каждой стороны. В файле chr4.genome указана длина хромосомы в нуклеотидах, чтобы запрашиваемое расширение координат не вышло за допустимые пределы.
Источники
[1] wikipedia.org
[2] genenames.org
[3] Биохимия
|
