व्यावहारिक जैव सूचना विज्ञान, भाग 2

यह आलेख आपको बताएगा कि पाइपलाइन के बाद प्राप्त डेटा को कैसे संसाधित किया जाए, जिसके आउटपुट में एक सैम / बैम फ़ाइल होगी [1], एक साधारण बेड ग्राफ फ़ाइल बनाएं (http://genome.ucsc.edu/FAQ/FAQformat.html) और इसे देखें UCSC जीनोम ब्राउज़र का उपयोग कर [2]। यह तय करना बहुत मुश्किल है कि प्रोग्राम को क्या लिखना है, क्योंकि पहले से ही अन्य लोगों की उपलब्धियों की एक बड़ी संख्या है और मैं वास्तव में एक पहिया की रचना नहीं करना चाहता हूं जहां यह चरण पहले ही पूरा हो चुका है। लंबे समय से परेशान, मैंने सी ++ पर रहने का फैसला किया, हालांकि पायथन और आर को एक समान पायदान पर माना जाता था। इस विचार को भी संरक्षित किया गया था कि ग्राफिक्स की आवश्यकता हो सकती है, और लिनक्स के तहत भी, इसलिए क्यूटी को C ++ में जोड़ा गया था। मुझे उम्मीद है कि इस लेख में मैं आपको यात्रा की शुरुआत में मेरे समक्ष रखे गए प्रश्न का उत्तर देने के लिए उपरोक्त सभी के बारे में पर्याप्त विस्तार से बताऊंगा और कहानी के पहले भाग में आवाज दी जाएगी।

Sam / bam फ़ाइलों के साथ काम करने के लिए, हमें इकट्ठे samtools पैकेज [1] की आवश्यकता है। Samtools.sourceforge.net से पैकेज डाउनलोड करें, इसे एक निर्देशिका में विस्तारित करें, उदाहरण के लिए, / usr / src, बनाई गई निर्देशिका पर जाएं और टाइप करें। मेरे पास XCurses मेरे सिस्टम पर स्थापित नहीं था और मैंने "LIBCURSES = -Lncurses" के साथ "LIBCURSES = -LXCurses" पंक्ति को प्रतिस्थापित किया और यह सब एक साथ आया। मेक प्रोग्राम का परिणाम संकलित समतुल्य कार्यक्रम और लिबाम.का पुस्तकालय था।
हमें C ++ क्लासेस की आवश्यकता है जो रीड, जीन, इंट्रोन्स, एक्सॉन, आदि के बारे में जानकारी स्टोर करेगी। इस तरह की कक्षाओं को व्यवस्थित करने के लिए, मैंने boost.intervals का उपयोग किया, हालांकि मैं जिस रूप में चाहता था, उस रूप में नहीं। न तो बूस्ट.इनवॉटर और न ही बूस्टर ।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 files के साथ काम करने के लिए रैपिंग फंक्शन शुरू करते हैं। सैम और बैम के बीच का अंतर यह है कि पहला एक सादा पाठ फ़ाइल है, हालांकि स्वरूपित है, और दूसरा संकुचित और द्विआधारी है, अर्थात्। संरचित। मैं बड़ी संख्या में रीड्स के साथ bam के साथ काम करना पसंद करता हूं, सैम फाइल गीगाबाइट तक पहुंच सकती है।
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 

मैं लेख में बाकी कोड नहीं दूंगा, यहां आर्काइव का लिंक दिया गया है , संरचना काफी सरल है: रूट में दो थर्डपार्टी और src निर्देशिकाएं हैं, पहले में samtools और दूसरे में sam2bedgraph और ग्लोबल। प्रोजेक्ट बनाने के लिए, आपको sam2bedgraph डायरेक्टरी में qmake को रन करना होगा और फिर बनाना होगा। देशी Qt (4.7.4) और बूस्ट (1.46.1) के साथ ओपनएसयूएसई 12.1 x64 के तहत जांच की गई।

बेडग्राफ फ़ाइल संरचना के विवरण के लिए एक लिंक पहले पैराग्राफ में दिया गया है, लेकिन मैं संक्षेप में इसका उल्लेख करूंगा। पहली पंक्ति फ़ाइल विशेषताओं को सेट करती है, यदि मौजूद है। फ़ाइल के शरीर में कम से कम 4 कॉलम होते हैं। पहले कॉलम में, क्रोमोसोम के नाम को इंगित किया गया है, दूसरे और तीसरे में - क्रमशः खंड ("विंडो") की शुरुआत और अंत के निर्देशांक। अंतिम कॉलम में - "ऊंचाई", हमारे मामले में, यह इस विंडो पर रीड की संख्या है। प्रत्येक फ़ाइल के लिए "विंडो" के आयाम प्रयोगात्मक रूप से चुनना बेहतर है। मैं Chip-seq के लिए 200 bp विंडो और RNA-seq के लिए 20 bp विंडो का उपयोग करता हूं।
एक प्रोग्राम के दो पैरामीटर हो सकते हैं: इनपुट .bam फ़ाइल और आउटपुट फ़ाइल, जहाँ बेडग्राफ के बारे में डेटा को लिखा जाता है। पैरामीटर सेट नहीं किया जा सकता है, तो प्रोग्राम डिफ़ॉल्ट फ़ाइल input.bam को खोलने और आउटपुट output.data बनाने का प्रयास करेगा। कॉन्फ़िगरेशन फ़ाइल के साथ उपयोगकर्ता के घर निर्देशिका में एक निर्देशिका बनाई जाएगी ।config / CCHMC / sam2bedgraph.ini, जिसमें आप डिफ़ॉल्ट फ़ाइल नाम मानों को बदल सकते हैं, "विंडो" की लंबाई बदल सकते हैं, जिसके लिए हम ऊंचाई पर विचार करते हैं और पहली सूचना रेखा की उपस्थिति सेट करते हैं।
परिणामी फ़ाइल को जीनोम ब्राउज़र में डाउनलोड किया जा सकता है। हम साइट genome.ucsc.edu पर जाते हैं [2] लिंक जीनोम पर क्लिक करें।


हम अगले पृष्ठ पर पहुंचते हैं जहां उन परिस्थितियों को सेट करना आवश्यक है जिनके तहत हमारी बैम फ़ाइल प्राप्त हुई थी: क्लेड (इवोल्यूशन लाइन), जीनोम (जीनोम), असेंबली (सार, असेंबली)। और फिर "कस्टम ट्रैक जोड़ें" बटन पर क्लिक करें।


अब आप “फ़ाइल चुनें” बटन पर क्लिक करें और “सबमिट” बटन पर क्लिक करके हमारी फ़ाइल डाउनलोड कर सकते हैं। नतीजतन, हम एनोटेशन के साथ पृष्ठ पर पहुंचते हैं, जो डेटाबेस और हमारे बेडग्राफ से जानकारी प्रदर्शित करेगा। यदि आप जीनोम ब्राउज़र को स्थानीय रूप से स्थापित करते हैं, तो आप अस्थायी रूप से सत्र के भीतर साइट पर बेडग्राफ जोड़ सकते हैं, लेकिन लगातार। आप साइट के लिए विशेष लिंक व्यवस्थित कर सकते हैं, जो इंगित करेगा कि बेडग्राफ फ़ाइल कहाँ प्राप्त करें। साइट जीनोम ब्राउज़र को कॉपी करने के लिए निर्देश प्रदान करती है।
आप कई फाइलें अपलोड कर सकते हैं, जो मैंने किया। आप निम्न स्क्रीनशॉट में परिणाम देखें।


जैसा कि आप इस स्क्रीनशॉट में देख सकते हैं, शून्य स्तर वाले क्षेत्र हैं, और संवर्धन वाले क्षेत्र हैं, जहाँ चोटियाँ हैं। कभी-कभी पृष्ठभूमि या शोर दिखाई देता है। यह भी स्पष्ट रूप से देखा गया है कि ChIP-seq RNA-seq से कैसे भिन्न होता है। कार्यक्रम से इन क्षेत्रों की परिभाषा एक अलग सवाल है, जिसके लिए एक भी जवाब नहीं है। तथ्य यह है कि जीन के रूप में जीन के रूप में कई चोटियां हो सकती हैं या, और इससे भी बदतर, कितने एक्सन हैं। और प्रयोग के ढांचे में समृद्ध क्षेत्रों को अलग करना बहुत श्रमसाध्य है।

1. ली, एच।, एट अल।, द सीक्वेंस एलाइनमेंट / मैप फॉर्मेट और सैटम्स। जैव सूचना विज्ञान, 2009.25 (16): पी। 2078-9।
2. केंट, डब्ल्यूजे, एट अल।, यूसीएससी में मानव जीनोम ब्राउज़र। जीनोम रेस, 2002.12 (6): पी। 996-1006।
3. ड्रेज़र, टीआर, एट अल।, यूसीएससी जीनोम ब्राउज़र डेटाबेस: एक्सटेंशन और अपडेट 2011। न्यूक्लिक एसिड रेज, 2012। 40 (डेटाबेस मुद्दा): पी। D918-23।