Эта статья расскажет о том, как обработать данные, полученные после pipeline, выходом которого будет sam/bam файл[1], создать несложный bed graph файл (http://genome.ucsc.edu/FAQ/FAQformat.html) и просмотреть его с помощью UCSC genome browser[2]. Очень сложно решиться, на чем писать программы, ибо уже есть огромное количество чужих наработок и совсем не хочется сочинять колесо там, где этот этап уже пройден. Долго мучаясь, я решил остановиться на C++, хотя Python и R рассматривались на равных. Также сохранилась идея, что может понадобиться графика, да ещё и под Linux, поэтому к С++ прибавилось Qt. Надеюсь, в этой статье я расскажу достаточно подробно о всем выше перечисленном, чтобы ответить на вопрос, заданный мне в начале пути и озвученный в первой части повествования.
Для работы с sam/bam файлами нам понадобится собранный пакет samtools[1]. Скачиваем пакет с сайта samtools.sourceforge.net, разворачиваем в директорию, например, /usr/src заходим в создавшуюся директорию и набираем make. У меня в системе не было установлено XCurses и я заменил строчку “LIBCURSES= -lXCurses” на “LIBCURSES= -lncurses” и все собралось. Результатом работы программы make стала собранная программа samtools и библиотека libbam.a.
Нам нужны C++ классы, которые будут хранить информацию о ридах, генах, интронах, эксонах и т.д. Для организации таких классов я воспользовался boost.intervals, хотя и не совсем в том виде, как хотелось. Ни boost.intervals, ни boost.icl не позволяют сохранять полную информацию об отрезках. Мне нужна следующая информация о множестве отрезков:
В частности, эта информация ответит на вопрос, какова высота покрытия отрезками текущего участка координатной оси.
С каждой новой статьей я буду пополнять и изменять эти классы. Сейчас их достаточно, чтобы собрать программу и прикинуть чего не хватает для будующего. Вот пример получившихся классов.
Reads.hpp
Теперь приступим к обертке над функциями для работы с sam файлами. Отличие между sam и bam заключается в том, что первый — это обычный текстовый файл, хоть и форматированный, а второй — сжатый и бинарный, т.е. структурированный. Я предпочитаю работать с bam при большом количестве ридов, sam файл может достигать гигабайтов.
BEDHandler.hpp
BEDHandler.cpp
Оставшуюся часть кода я приводить в статье не буду, вот ссылка на архив, структура довольно простая: в корне две директории thirdparty и src, в первой лежит samtools, во второй — sam2bedgraph и global. Для того, чтобы собрать проект, надо в директории sam2bedgraph запустить qmake и затем make. Проверял под openSUSE 12.1 x64 с родной Qt (4.7.4) и boost (1.46.1).
Ссылка на описание структуры bedgraph файла приведена в первом абзаце, но вкратце упомяну. Первая строчка задает характеристики файла, если присутствует. Тело файла содержит как минимум 4 столбца. В первом столбце указывается название хромосомы, во втором и третьем — координаты начала и конца отрезка («окна»), соответственно. В последнем столбце — «высота», в нашем случае это количество начал ридов на это окно. Размеры «окна» для каждого файла лучше подбирать экспериментально. Я для ChIP-seq использую окно в 200 bp, для RNA-seq окно в 20 bp.
У программы может быть два параметра: входной .bam файл и выходной файл, куда запишутся данные о bedgraph. Параметры можно и не задавать, тогда программа попытается открыть файл по умолчанию input.bam и создать выходной output.data. В домашней директории пользователя будет создана директория с конфигурационным файлом .config/CCHMC/sam2bedgraph.ini, в котором можно поменять значения имен файлов по умолчанию, изменить длину «окна», для которого считаем высоты и задать вид первой информационной строчки.
Получившийся файл можно загрузить в genome browser. Идем на сайт genome.ucsc.edu[2] нажимаем на ссылку Genomes.
Попадаем на следующую страницу, где необходимо задать условия, при которых был получен наш bam файл: Clade(Эволюционная линия), genome (Геном), assembly (Аннотация, Сборка). И затем нажать кнопку “add custom tracks”.
Теперь можно нажать кнопку «Choose file» и загрузить наш файл, нажав кнопку «Submit». В результате попадаем на страницу с аннотацией, на которую будет выведена информация из базы данных и из нашего bedgraph. Если установить genome browser локально, то можно добавлять bedgraph на сайт не временно в рамках сессии, а постоянно. Можно организовывать специальные линки на сайт, в которых будет указано, где брать bedgraph файл. На сайте приведена инструкция, как скопировать genome browser.
Можно загрузить несколько файлов, что я и сделал. Результат вы видите на следующем скриншоте.
Как видно на этом скриншоте, есть участки с нулевым уровнем, а есть участки с обогащением, там где пики. Иногда виден фон или шум. Также наглядно видно, чем отличается ChIP-seq от RNA-seq. Определение этих участков из программы — отдельный вопрос, однозначного ответа на который нет. Дело в том, что пиков может быть столько, сколько генов в геноме или, и того хуже, сколько экзонов. И отделить обогащённые участки в рамках эксперимента очень трудоемко.
1. Li, H., et al., The Sequence Alignment/Map format and SAMtools. Bioinformatics, 2009. 25(16): p. 2078-9.
2. Kent, W.J., et al., The human genome browser at UCSC. Genome Res, 2002. 12(6): p. 996-1006.
3. Dreszer, T.R., et al., The UCSC Genome Browser database: extensions and updates 2011. Nucleic Acids Res, 2012. 40(Database issue): p. D918-23.
Для работы с sam/bam файлами нам понадобится собранный пакет samtools[1]. Скачиваем пакет с сайта samtools.sourceforge.net, разворачиваем в директорию, например, /usr/src заходим в создавшуюся директорию и набираем make. У меня в системе не было установлено XCurses и я заменил строчку “LIBCURSES= -lXCurses” на “LIBCURSES= -lncurses” и все собралось. Результатом работы программы make стала собранная программа samtools и библиотека libbam.a.
Нам нужны C++ классы, которые будут хранить информацию о ридах, генах, интронах, эксонах и т.д. Для организации таких классов я воспользовался boost.intervals, хотя и не совсем в том виде, как хотелось. Ни boost.intervals, ни boost.icl не позволяют сохранять полную информацию об отрезках. Мне нужна следующая информация о множестве отрезков:
- любая пара начало/конец, начало/длина;
- все отрезки закрытые, т.е. включают концы;
- сколько отрезков начинается в заданной точке;
- сколько отрезков пересекается в точке;
- сколько отрезков начинается в заданном отрезке или пересекается с заданным отрезком.
В частности, эта информация ответит на вопрос, какова высота покрытия отрезками текущего участка координатной оси.
С каждой новой статьей я буду пополнять и изменять эти классы. Сейчас их достаточно, чтобы собрать программу и прикинуть чего не хватает для будующего. Вот пример получившихся классов.
Reads.hpp
#ifndef _READS_29122011_HPP_ #define _READS_29122011_HPP_ #include <boost/numeric/interval.hpp> #include <QString> #include <QMap> #include <QtDebug> namespace genome { namespace bni = boost::numeric; typedef bni::interval<int> read_position; /********************************************************************************** **********************************************************************************/ class Read { public: Read(): multiplying(1), length(0){}; Read(Read const &r): multiplying(r.multiplying), length(r.length), position(r.position){ sentenceRepresentation=r.sentenceRepresentation; qualityRepresentation=r.qualityRepresentation;}; Read(int start,int len,QString sr="",QString qr=""): multiplying(1), length(len), position(start,start+len-1), sentenceRepresentation(sr), qualityRepresentation(qr){}; int getLevel() {return multiplying;}; void plusLevel() {++multiplying;}; int getStart() {return position.lower();}; int getLength() {return length;}; void operator+= (const int& c) {this->multiplying+=c;}; bool operator== (const Read& r) const {return this->position==r.position;}; bool operator!= (const Read& r) const {return this->position!=r.position;}; void operator++ (int) {this->multiplying++;}; private: int multiplying; int length; read_position position; QString sentenceRepresentation; QString qualityRepresentation; }; /********************************************************************************** **********************************************************************************/ typedef QMap<int,Read> cover_map; class Cover { public: Cover():max_len(0){}; void add(Read&); int getHeight(int); int getHeight(int,int); int getStarts(int); int getStarts(int,int); QList<int> getStarts(); cover_map::iterator getBeginIterator(){return covering.begin();}; cover_map::iterator getEndIterator(){return covering.end();}; bool operator== (const Cover& c) const {return this==&c;}; bool isEmpty(){return covering.size()==0;}; // static Cover empty(){ return Cover();}; private: cover_map covering; int max_len; }; /********************************************************************************** **********************************************************************************/ class Lines { public: Lines(){}; Lines(Lines&){}; void addLine(QString, Read&); Cover& getLineCover(QString); QList<QString> getLines(void) { return lines.keys(); }; /* */ private: QMap<QString,Cover> lines; }; /********************************************************************************** **********************************************************************************/ class GenomeDescription: public Lines { public: quint64 notAligned; // number of reads (ussualy form sam/bam file) that are not aligned quint64 total; /* */ void setGene(const QChar &sense,const QString &chrName,const qint32 &pos,const qint32 &,const qint32 &len) { Read r(pos,len); addLine(chrName+sense,r); }; /* */ GenomeDescription():Lines(), notAligned(0), total(0) {}; }; } #endif
Теперь приступим к обертке над функциями для работы с sam файлами. Отличие между sam и bam заключается в том, что первый — это обычный текстовый файл, хоть и форматированный, а второй — сжатый и бинарный, т.е. структурированный. Я предпочитаю работать с bam при большом количестве ридов, sam файл может достигать гигабайтов.
BEDHandler.hpp
#ifndef BEDHandler_H #define BEDHandler_H #include <config.hpp> #include <Reads.hpp> template <class Storage, class Result> class BEDHandler : public QState { private: Storage *sam_input; Result *output; QSettings setup; QFile _outFile; public: BEDHandler(Storage &sam,Result &output,QState *parent=0); ~BEDHandler(); protected: virtual void onEntry(QEvent* event); }; #include <BEDHandler.cpp> #endif
BEDHandler.cpp
//------------------------------------------------------------- //------------------------------------------------------------- template <class Storage,class Result> BEDHandler<Storage,Result>::~BEDHandler() { } //------------------------------------------------------------- //------------------------------------------------------------- template <class Storage,class Result> BEDHandler<Storage,Result>::BEDHandler(Storage& sam,Result &_output,QState * parent): QState(parent), sam_input(&sam), output(&_output) { _outFile.setFileName(setup.value("outFileName").toString()); _outFile.open(QIODevice::WriteOnly|QIODevice::Truncate); } //------------------------------------------------------------- //------------------------------------------------------------- template <class Storage,class Result> void BEDHandler<Storage,Result>::onEntry(QEvent*) { if(!setup.contains("graphWindow")) setup.setValue("graphWindow",200); if(!setup.contains("siteShift")) setup.setValue("siteShift",75); if(!setup.contains("separateStrand")) setup.setValue("separateStrand",false); if(!setup.contains("HeaderString")) setup.setValue("HeaderString","track type=bedGraph name=%1"); if(setup.value("HeaderString").toString()!="") { _outFile.write((setup.value("HeaderString").toString().arg(_outFile.fileName())+"\n").toLocal8Bit()); _outFile.flush(); } quint32 window=setup.value("graphWindow").toUInt(); quint32 shift= setup.value("siteShift").toUInt(); QString line; QString chrome; foreach(line,sam_input->getLines()) { if(line.endsWith("-")) continue; chrome=line; chrome.truncate(line.length()-1); QMap <int,int> bed; { genome::cover_map::iterator i=sam_input->getLineCover(chrome+QChar('+')).getBeginIterator(); genome::cover_map::iterator e=sam_input->getLineCover(chrome+QChar('+')).getEndIterator(); while(i!=e) { int val=i.key()+shift; bed[val-val%window]+=i.value().getLevel(); ++i; } } { genome::cover_map::iterator i=sam_input->getLineCover(chrome+QChar('-')).getBeginIterator(); genome::cover_map::iterator e=sam_input->getLineCover(chrome+QChar('-')).getEndIterator(); while(i!=e) { int val=i.key()-shift; if(val<0) val=0; bed[val-val%window]+=i.value().getLevel(); ++i; } } QMap<int,int>::iterator i = bed.begin(); for(;i!=bed.end();i++) { _outFile.write(QString(chrome+"\t%1\t%2\t%3\n"). arg(i.key()).arg(i.key()+window).arg(i.value()).toLocal8Bit()); _outFile.flush(); } } }//end of function
Оставшуюся часть кода я приводить в статье не буду, вот ссылка на архив, структура довольно простая: в корне две директории thirdparty и src, в первой лежит samtools, во второй — sam2bedgraph и global. Для того, чтобы собрать проект, надо в директории sam2bedgraph запустить qmake и затем make. Проверял под openSUSE 12.1 x64 с родной Qt (4.7.4) и boost (1.46.1).
Ссылка на описание структуры bedgraph файла приведена в первом абзаце, но вкратце упомяну. Первая строчка задает характеристики файла, если присутствует. Тело файла содержит как минимум 4 столбца. В первом столбце указывается название хромосомы, во втором и третьем — координаты начала и конца отрезка («окна»), соответственно. В последнем столбце — «высота», в нашем случае это количество начал ридов на это окно. Размеры «окна» для каждого файла лучше подбирать экспериментально. Я для ChIP-seq использую окно в 200 bp, для RNA-seq окно в 20 bp.
У программы может быть два параметра: входной .bam файл и выходной файл, куда запишутся данные о bedgraph. Параметры можно и не задавать, тогда программа попытается открыть файл по умолчанию input.bam и создать выходной output.data. В домашней директории пользователя будет создана директория с конфигурационным файлом .config/CCHMC/sam2bedgraph.ini, в котором можно поменять значения имен файлов по умолчанию, изменить длину «окна», для которого считаем высоты и задать вид первой информационной строчки.
Получившийся файл можно загрузить в genome browser. Идем на сайт genome.ucsc.edu[2] нажимаем на ссылку Genomes.
Попадаем на следующую страницу, где необходимо задать условия, при которых был получен наш bam файл: Clade(Эволюционная линия), genome (Геном), assembly (Аннотация, Сборка). И затем нажать кнопку “add custom tracks”.
Теперь можно нажать кнопку «Choose file» и загрузить наш файл, нажав кнопку «Submit». В результате попадаем на страницу с аннотацией, на которую будет выведена информация из базы данных и из нашего bedgraph. Если установить genome browser локально, то можно добавлять bedgraph на сайт не временно в рамках сессии, а постоянно. Можно организовывать специальные линки на сайт, в которых будет указано, где брать bedgraph файл. На сайте приведена инструкция, как скопировать genome browser.
Можно загрузить несколько файлов, что я и сделал. Результат вы видите на следующем скриншоте.
Как видно на этом скриншоте, есть участки с нулевым уровнем, а есть участки с обогащением, там где пики. Иногда виден фон или шум. Также наглядно видно, чем отличается ChIP-seq от RNA-seq. Определение этих участков из программы — отдельный вопрос, однозначного ответа на который нет. Дело в том, что пиков может быть столько, сколько генов в геноме или, и того хуже, сколько экзонов. И отделить обогащённые участки в рамках эксперимента очень трудоемко.
1. Li, H., et al., The Sequence Alignment/Map format and SAMtools. Bioinformatics, 2009. 25(16): p. 2078-9.
2. Kent, W.J., et al., The human genome browser at UCSC. Genome Res, 2002. 12(6): p. 996-1006.
3. Dreszer, T.R., et al., The UCSC Genome Browser database: extensions and updates 2011. Nucleic Acids Res, 2012. 40(Database issue): p. D918-23.
