evfratov (evfratov) wrote,
evfratov
evfratov

Запиливаем rnaSeqMap

Сначала ставим R-пакет rnaSeqMap, пока оно ставится - оно скажет что ему именно нужно доставлять. Обычно нужно дополнительно нексколько девелоперских вариантов базовых пакетов. Потом ставим R-пакет xmapcore (основа системы работы с данными чипОв), а для него организуем базу mysql xmapcore_{organizm}_{NN} и в него-же будут потом импортироваться данные из SAM файлов. Непосредственно для работы только с SAM база не должна быть нужна, но в xmap-базу нужно импортировать SAM файлы после их обработки специальным awk скриптом.
Заводим mysql: запускаем от mysql рута: mysql -uroot -p , создаём пользователя > create '{USER}'@'localhost'; позволяем ему всё grant all on *.* to '{USER}'@'localhost'; и выходим из рута mysql. Входим mysql просто: там create database xmapcore_{organism}_{NN}; выходим.
Формируем базу xmapcore: качаем с сайта xmapcore нужную базу, распаковываем, cd <куда мы её распаковали> открываем ./importdb.sh. В нём меняем имя пользователя (строка 17) на своё. И запускаем - должны лишь раз 30 спросить пароль к базе (себе оставил пустым). По окончании база должна быть готова.
Составляем скрипт-преобразователь SAM файлов для подгтовки к импорту:
В оригинале (руководстве) скрипт выглядит так -

# For mouse 19 chromosomes + 20=X,21=Y,22=MT
{
    OFS = "\t"
    strand = rshift($2,4) % 2
    if (strand==0) strand = -1
    cc = substr($3,4)
    if (cc!=cc+0)
    {
        if (cc=="X") cc=20
        else if (cc=="Y") cc=21
        else if (cc=="M") cc=22
        else cc=30;
    }
print(6, cc, $4, $4+ length($10),strand)
}
Ориентация цепи выполняется с помощью интерпретации 2 колонки SAM ("байтовое" представление цепи в  SAM) при помощи хитрой операции байтового сдвига. Имя хромосомы заменяется на числовую величину, предварительно забирая в "сс" подстроку с 4й буквы (вероятно, отсечение "chr" из "chr{N}") или заменой с буквы (X, Y или MT) на число, для этого выполняется проверка корректности числовой операции с не числом. В самом скрипте на выводе первый числовой параметр - номер сэмпла, задаётся руками.

Скрипт на выходе даёт 5 столбчатый поток "Номер сэмпла Хромосома Координата где лёг рид Координата ГЛР + длина рида Цепь". Проверяем - можно-ли загнать в базу буквенные имена хромосом (поскольку у меня гены из библиотеки): формируем по инструкции файлы (в моём случае 12 со сквозной нумерацией), пытаеся затолкать в базу xmap - оно отказывается, т.к. не существует соответствующей таблицы (поскольку mysql загружает данные в таблицу по имени файла с отсечением расширения). Созданием в ручную таблицы "n chr start end strand" с параметрами int(11) varchar(8) int(11) int(11) int(11) и последующи слиянием файлов экспериентов в тестовой базе удалось как-то загрузить и повыполнять запросы (напоинаем себе что оно было загружено в базу test в таблицу seq_reads).
Пытаемся подключиться к базе через R: не забываем загрузить rnaSeqMap, пытаемся включить xmap.connect - оказывается что Can't find database config file '/home/evfr/.xmapcore/databases.txt'. See package installation instructions for help setting this up. В установочном мане расписано как форировать файл баз.
Формируем файл баз xmapcore: пишем из мана табулированный файл с кучей очевидных дефолтов, убеждаемся в том что в mxapcore просто так новую базу не добавить и формировать базы ридов надо как дополнительные таблицы в xmapcore_{organism}_{NN}.
Формируем таблицу в xmapcore для ридов: TINYINT [UNSIGNED] для num, VARCHAR(80) для chr, INT [UNSIGNED] для start и end, BOOL для chain. В общем create table mtor(num TINYINT UNSIGNED, chr TINYINT UNSIGNED, start INT UNSIGNED, end INT UNSIGNED, chain BOOL); Таблица есть, выходим, называем слитый файл ридов mtor. Дальше  mysqlimport -v -u evfr -p xmapcore_homo_sapiens_66 --local mtor.txt xmapcore_homo_sapiens_66.mtor: Records: 431384  Deleted: 0  Skipped: 0  Warnings: 96 - хм.. ну что-то работает.
Проверяем функционирование через R: видим что "> data(mtor) Предупреждение In data(mtor) : данные ‘mtor’ не найдены", ищем наконец скрипты snaReqMap rnaSeqMap.sql и indices.sql, находим их в папке home/R/<пакеты>/rnaSeqMap/Scripts/*, читаем, перечёркиваем половину вышестоящего текста, удаляем mtor, кидаем mtor.txt в папку с этими скриптами, понимаем что ничего не работает и начинаем разбираться.
Формируем таблицу ридов в xmapcore: скрипт importdb.sh работает через последовательные исполнения
rnaSeqMap.sql - creates tables, partitioning
bio_sample.sql - data for bio_sample table
seq_read.sql - example data from sequencing for seq_read table
procedures.sql - creates stored procedures
indices.sql - creates indices
при этом в исполняемом скрипте и indices.sql стоит ALTER TABLE `xmapcore_homo_sapiens_58`.`seq_read` , а база уже новая (66) и ещё это строго привязывает имя таблицы и файла. А в самом исполняемом скрипте вообще умудрились опечататься:
$MYS    -h $SVR -u $USR $PWD xmapcore_homo_sapiens_66 < rnaSeqMap.sql
$MYS    -h $SVR -u $USR $PWD xmapcore_homo_sapiens_66 < bio_sample.sql
$MYS    -h $SVR -u $USR $PWD xmapcore_homo_sapiens_66 < seq_read.sql
$MYS    -h $SVR -u $USR $PWD xmapcore_homo_sapiens_66 < procedures.sql
$MYS    -h $SVR -u $USR $PWD xmapcore_homo_sapiens_66 < indeces.sql
    indeces вместо indices, в результате чего база не формировалась.
При этом в ходе работы какие-то проблемы с bio_sample - ERROR 1054 (42S22) at line 1: Unknown column 'B' in 'field list' где insert into bio_sample values ((1,B),(2,T),(3,B),(4,B),(5,T),(6,T));
Проверяем в R: коннектимся, запрашиваем пресловутый data(sample_data_rnaSeqMap) и наконец-то не получам ошибок. Пытаемся что-то сделать, видим какую-то лажу, понимаем что импортировался только набор примеров, перечитываем руководтсво и видим что myslqimport -u mysqlexampleuser -p xmapcore_homo_sapiens_58 -local seq_reads , а The seq_reads table has then to be indexed as in the script inst scripts indices.sql , т.е. для импорта данных (из файла seq_read) надо сначала создать таблицу для seq_reads `xmapcore_homo_sapiens_66`.`seq_read`  и после этого уже загонять в базу. Её должен потом индексировать скрипт indices.sql. Табличка  bio_sample должна хранить где и кто. Обе таблицы создаёт rnaSeqMap.sql, в принципе через его всё настраиваемо должно быть.
Формируем таблицу из экспериментальных данных: удаляем seq_read и bio_sample, сбрасывая изменения, понимаем что таким образом система в норме не совместима с библиотеками генов (только хромосомы, неприменимы имена генов), придётся адаптировать таблицу - например, заменить имена на номера. Заменяем, обнаруживаем что awk пропустил отсечение SAM заголовков (вероятно, они и провоцировали предупреждения, т.к. в них совсем другого типа данные), удаляем руками (простые идентификаторы). Кидаем в рабочую папку indices.sql и на всякий случай rnaSeqMap.sql. Начинаем: mysql -v -u evfr -p xmapcore_homo_sapiens_66 < ./rnaSeqMap.sql  , проверяем: mysql - таблицы те есть, описание выглядит нормально. Загружаем: mysqlimport -v -u <USER> -p xmapcore_homo_sapiens_66 --local ./seq_read.txt  - загрузились без багов, индексируемся: mysql -v -u evfr -p xmapcore_homo_sapiens_66 < ./indices.sql  -- вроде без багов. Проверяем: семплы пустые, риды полны данных - наконец-то загрузились в MySQL.
Проверяем в R: наблюдаем как обычно - в rs только примеры, в мане указано что эксперименты идентифицируются по таблице bio_sample - там всего лишь номер эксперимента (число 2 байта) и описание (строка 40 байт). Запишем вручную: mysqlimport -v -u evfr -p xmapcore_homo_sapiens_66 --local ../PROC/RNAmod/mTOR/bio_sample.tab  без лагов. Проверяем: > str(rs)
Formal class 'SeqReads' [package "rnaSeqMap"] with 12 slots  ну хоть что-то.
Пытаемся работать: экспорт с лагами, но есть, картинки строятся, считать что работает.
Tags: bioinformatics, Учёба
Subscribe
  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic

    Your IP address will be recorded 

  • 1 comment