پیشگویی دقیق ساختار دوم RNA مبتنی بر الگوریتم ژنتیک 

سهیلا منتصری¹، نصرالله مقدم-چرکری^2* و فاطمه زارع میرک آباد³

¹ تهران، دانشگاه تربیت مدرس، دانشکده علوم ریاضی، گروه علوم کامپیوتر

² تهران، دانشگاه تربیت مدرس، دانشکده مهندسی برق و کامپیوتر، گروه مهندسی کامپیوتر

³ تهران، دانشگاه صنعتی امیرکبیر، دانشکده ریاضی و علوم کامپیوتر، گروه علوم کامپیوتر

تاریخ دریافت: 1/9/90                 تاریخ پذیرش: 4/5/91

چکیده

مولکول RNA نقش مهم و اساسی در فرآیندهای زیستی ایفاء می­کند. در بیشتر مواقع، عملکرد RNAها توسط ساختار آنها مشخص می­شود. با توجه به پیچیدگی و هزینه بر بودن روشهای آزمایشگاهی برای پیشگویی ساختار RNAها، از روشهای محاسباتی استفاده می­گردد. الگوریتمهای متنوعی جهت پیشگویی ساختار دوم مولکول RNA وجود دارد. در این مقاله، یک الگوریتم  ژنتیک به نام RNAG جهت پیشگویی ساختار دوم مولکول RNA براساس حداقل انرژی آزاد ارائه می­شود. در این الگوریتم، هر فرد از جمعیت شامل تعدادی ساقه می­باشد. افراد براساس مقدار برازندگی حداقل انرژی آزاد شده از ساقه­ها و حلقه­ها به­ترتیب صعودی رتبه­بندی شده و در ادامه به­ترتیب عملگرهای تقاطع و جهش روی آنها برای ایجاد نسل بعد اجرا می گردد. فرآیند تولید نسل تا زمان تولید یک فرد با حداقل انرژی آزاد مناسب ادامه می­یابد. در پایان این فرد به­عنوان ساختار دوم بهینه در نظر گرفته می­شود. الگوریتم پیشنهادی روی تعدادی از RNAها در باکتریها اجرا می­گردد. نتایج حاصل از این تحقیق نشان می­دهد که الگوریتم RNAG در مقایسه با سایر روشهای مشابه دارای دقت بسیار بالا است.

واژه های کلیدی: حداقل انرژی آزاد، ساقه، مقدار برازندگی.

* نویسنده مسئول، تلفن: 82883301، پست الکترونیکی: charkari@modares.ac.ir

مقدمه

مولکولهای RNA در تمام موجودات زنده دارای نقش حیاتی هستند. شناخت ساختار RNA در درک فعالیت آن اهمیت فراوانی دارد )8(. ساختار مولکولهای RNA در بیان ژن، پیرایش RNAهای پیک (Messenger RNA)، ساخت پروتئین و عملکردهای زیستی دیگر مؤثر است (1،2و11). به­عنوان مثال، خاتمه رونویسی (Transcription) بعضی از ژنها در باکتری، براساس ساختار سنجاق­سری انتهایی RNA پیک انجام می­شود )14(. از نکات مهم در پیشگویی ساختار دوم RNA می­توان به دنباله­هایی اشاره کرد که هنوز ساختار آنها از راه آزمایش مشخص نشده و در نتیجه هیچ نظیری در پایگاه داده­ برای آنها نمی­توان یافت. از اهداف مهم پیشگویی ساختارها حل این مشکل است (16). موضوع دیگری که کاربرد مهمی جهت طراحی ساختارهای RNA دارد، برهمکنش دو مولکول RNA است )13(. پیشگویی ساختارهای RNA مقدمه­ای در تعیین ساختار برهمکنش دو RNA می­باشد.

تلاشهایی جهت پیشگویی ساختار دوم مولکولهای RNA با بیشینه کردن تعداد جفت­بازها (Base pairs)، با استفاده از برنامه­نویسی پویا انجام شد که در آن بهترین ساختار برای هر زیردنباله­ محاسبه می­گردد )10(. پس از آن الگوریتم مشابهی ارائه شد که در آن از مقادیر انرژی آزاد­ جفت­بازها برای محاسبه ساختاری با کمترین انرژی آزاد (Minimum free energy) استفاده می­شود (9 و 18). در یک روش دیگر برای پیشگویی ساختار دوم، توابع تسهیم (Partition function) مولکولهای RNA براساس برنامه نویسی پویا محاسبه­ می­گردند )7(. ابزار MFold )19( توسط پارامترهای موجود برای محاسبه ساختار دوم RNA )15(، به پیشگویی ساختار دوم می­پردازد. در تعدادی از رویکردها، انرژی آزاد با استفاده از مدل ترمودینامیکی نزدیکترین همسایه تعیین می­شود. در این مدلها، انرژی آزاد ساختار به­عنوان مجموع انرژیهای آزاد شده از هر ساقه (Stem) و حلقه با استفاده از داده­های ترمودینامیکی محاسبه می­گردد )6 و 17(. روشی براساس گرامرهای مستقل ازمتن ارائه شد که در آن از الگوریتمهای آماری برای ایجاد ساختار دوم استفاده می­شود )12(. ابزار RNAFold )3( با استفاده از پارامترهای انرژی ایجاد شده )5( ساختار دوم RNA را پیشگویی می­کند.

در این مقاله، یک الگوریتم ژنتیک به نام RNAG جهت پیشگویی دقیق ساختار دوم مولکول RNA ارائه می­شود. در این الگوریتم، یک ماتریس­ نقطه­ای ایجاد می­گردد که نشان­دهنده تمام جفت­بازهای ممکن در RNA است. هر زیرقطر در ماتریس نقطه ای را می­توان به­عنوان یک ساقه در نظر گرفت. سپس جمعیتی از ساقه­هایی که به طور تصادفی انتخاب می­شوند، ایجاد شده و مقدار برازندگی (Fitness value) حداقل انرژی آزاد شده از ساقه­ها و حلقه­ها برای هر فرد موجود در جمعیت محاسبه می­گردد. برای ایجاد نسل جدید، عملگرهای تقاطع (Crossover) و جهش (Mutation) به­ترتیب روی تعدادی از افراد نسل جاری انجام می­شود. فرآیند تولید نسل ادامه می­یابد تا زمانی که انرژی آزاد فردی به حد مطلوب برسد. در نهایت، این فرد با حداقل انرژی آزاد جهت تشکیل ساختار دوم RNA انتخاب می­شود. الگوریتم پیشنهادی روی تعدادی از داده­ها شامل CopA، CopT، R1inv، R2inv، Tar، Tar*، DIS، IncRNA₅₄ و RepZ در باکتریها به کار رفته است. نتایج حاصل از این تحقیق نشان می­دهد که الگوریتم RNAG در مقایسه با سایر روشهای مشابه دقت بالایی دارد.

مواد و روشها 

پایگاه داده: RNAهای مورد بررسی در این مقاله شامل CopA، CopT، R1inv، R2inv، Tar، Tar*، DIS، IncRNA₅₄ و RepZ هستند (4).

تعاریف پایه: دنباله RNA از چهار نوع نوکلئوتید تشکیل می­شود که شامل آدنین (A)، گوانین (G)، سیتوزین (C) و یوراسیل(U)  است. هر RNA دارای دو انتهای مجزا است که به­عنوان ¢3 و ¢5 شناخته می­شوند. یک دنباله RNA به نام R ، |R| = n در جهت ¢5 به ¢3 به صورت زیر تعریف می­گردد:

R = r₁r₂ ... r_n : "i (1£i£n)r_i Î{A, C, G, U}. 

معکوسR  باr_n r_n-1... r₁. در جهت ¢3 به ¢5 مشخص می شود. بنابراین r_{i j} =r_i r_i+1... r_j زیردنباله­ای از R است که از موقعیت i شروع شده و به موقعیت j ختم می­گردد. دنباله RNA با تشکیل پیوند هیدروژنی بین بازهای آن تشکیل ساختار می­دهد. بیشتر پیوندها بین بازهای مکمل واتسون-کریک روی می­دهند که در آنها G با C و A با U جفت می­شوند و برعکس. این پیوندها می­توانند ساختار دوم RNA را تشکیل دهند.

ساختار دوم RNA از ساقه­ها و نواحی­ منفرد (Single regions) تشکیل می­شود. هر ساقه مجموعه­ای از جفت­بازهای مجاور مانند (r_i , r_j) و (r_i¢ , r_j¢) است به طوری که
j¢ i, i¢, j, توسط یکی از شرایط زیر ارضاء می گردند:

i < i¢ < j¢ < j

i¢ < i < j < j¢

به فرض اینکه r_{i j} و r_{k l} دو زیردنباله از RNA باشند که تشکیل ساقه می­دهند. بنابراین زیردنباله r_{i j} به معکوس r_{k l} متصل می­شود. به­عبارت دیگر، بین هریک از بازهای r_{i j} و معکوس r_{k l} به­ترتیب پیوند هیدروژنی برقرار می­گردد. برای نشان دادن ساقه در ساختار دوم RNA، هر باز در r_{i j} با ¢(¢ و هر باز در r_{k l} با ¢)¢ مشخص می­شود. نواحی منفرد به­عنوان حلقه­ یا تک­رشته­ای شده (Single-stranded) در نظر گرفته می­شوند. لازم به ذکر است که دو انتهای هر ناحیه حلقه به ساقه­ها متصل می گردند در حالی که تنها یک انتهای هر تک­رشته­ای شده به یک ساقه پیوند می خورد. هر باز در نواحی منفرد با ¢.¢ نشان داده می­شود. بنابراین S = s₁ ... s_n ساختار دوم RNA را نشان می­دهد که در آن برای هر باز i ، این فرمول 1£ i £ n ،
s_i Î{¢(¢,¢)¢,¢.¢} در نظر گرفته می شود.

الگوریتم ژنتیک روشی ​​است که در حل مسائل بهینه­سازی مورد استفاده قرار می­گیرد و براساس فرآیندهای ژنتیکی موجودات زنده است. الگوریتم ژنتیک با جمعیتی از افراد بیان می­شود که هر فرد نشان­دهنده راه­حلی برای مسئله است. مقدار برازندگی به هر فرد با توجه به میزان مناسب بودن آن به­عنوان یک راه­حل، اختصاص داده می­شود. افراد براساس مقدار برازندگی جفت­گیری کرده و نسل جدید تشکیل می­گردد. فرآیند تولید نسل تا زمانی ادامه می­یابد که راه­حل بهینه برای مسئله یافت شود.

روش پیشنهادی: در مسئله پیشگویی ساختار دوم RNA، یک RNA به­عنوان ورودی در نظر گرفته می­شود و هدف یافتن ساختار دوم RNA است که دارای حداقل انرژی آزاد باشد. تعریف دقیق­تر مسئله به شرح زیر است:

ورودی: یک RNA با دنباله R = r₁ r₂... r_n. در جهت ¢5 به ¢3.

خروجی: ساختار دوم  که با دنباله­ای از کاراکترهای  ¢(¢،  ¢)¢ و ‘.’ نشان داده می­شود.

در این مقاله، روش پیشنهادی برای حل مسئله ساختار دوم RNA براساس یک الگوریتم ژنتیک به نام RNAG است که شامل مراحل ایجاد جمعیت اولیه، عملگرهای تقاطع و جهش و شرط خاتمه الگوریتم می­باشد که در ادامه به توضیح آنها پرداخته می شود.

ایجاد جمعیت اولیه: چگونگی تولید جمعیت اولیه به­ترتیب زیر انجام می­شود:

3)   ماتریس نقطه­ای   برای دنباله R براساس بازهای مکمل واتسون-کریک ایجاد می­گردد که مقدار آن در موقعیت  به صورت زیر تعریف می­شود:

M^R[i,j] =  

به طوری که   و  به­ترتیب نشان­دهنده باز  ام و  ام در دنباله  برای هر  و ، ، هستند. (توجه کنید که در ایجاد این ماتریس جفت­­باز  در نظر گرفته نمی­شود چون دقت پیشگویی را کاهش می­دهد.)

4)   در ماتریس نقطه­ای ، تمام مقادیر مورب متوالی 1 که روی قطر اصلی یا موازی آن قرار داشته باشند به­عنوان یک زیرقطر در نظر گرفته می­شوند. مجموعه­ای از زیرقطرهای   به صورت زیر تعریف می­گردد:

D^R = { < i, j, k, l > |1£ i £ k £ n & 1 £ j £ l £ n }

به طوری که  و  به­ترتیب موقعیت­ شروع و پایان یک زیرقطر را مشخص می­کنند. فرض کنید  و  است. این زیرقطر نشان دهنده این است که زیردنباله  در  به زیردنباله  در معکوس  متصل می­شود. اگر  و  با شرایط  ،   و   موجود باشند، آنگاه  بایستی از مجموعه  حذف و دو زیرقطر   و   به مجموعه مورد نظر اضافه شوند. زیرا در این حالت تعدادی از بازهای تشکیل دهنده زیرقطر  دو مرتبه جهت تشکیل پیوند محاسبه می­گردند.

5)   جمعیت اولیه  براساس  به این صورت ساخته می­شود که برای هر ، ، مراحل زیر انجام می­پذیرد:

الف) ،  به طور تصادفی از مجموعه  ایجاد می­گردد. به بیان دقیق­تر، برای هر فرد (که ابتدا تهی است) در جمعیت، ابتدا یک زیرقطر تصادفی از  انتخاب شده و در آن قرار می­گیرد. زیرقطرهای بعدی نیز به طور تصادفی انتخاب شده و در صورتی در فرد قرار می­گیرند که با هیچ یک از زیرقطرهای موجود در آن همپوشانی نداشته باشند. اگر همپوشانی وجود داشته باشد، قسمتهای همپوشان از زیرقطر جدا شده و زیرقطر حاصل به مجموعه زیرقطرهای قبلی اضافه می­شود. فرض کنید  که  و . همپوشانی دو زیرقطر  و   به صورت زیر تعریف می­گردد:

Overlap(d₁, d₂) =  

ب) مقدار برازندگی فرد   به صورت زیر محاسبه می شود:

Fitness (C[i]) =

به طوری که   نشان­دهنده یک زیرقطر در مجموعه  است و  حلقه­ای در مجموعه­ حلقه­های هیرپین، بالج، داخلی و چندحلقه­ای در فرد را نشان می­دهد. حداقل انرژی آزاد شده از  به صورت زیر محاسبه می­گردد:

MFE (d¢ ) =

که  انرژی آزاد شده از دو جفت­باز مجاور  و  می­باشد و  نشان­دهنده مجموعه­ا­ی از جفت­بازها در زیرقطر است. حداقل انرژی آزاد شده از تمام دو جفت­بازهای مجاور که تشکیل ساقه می­دهند، در جدول1 نشان داده شده است. انرژی آزاد حلقه­های هیرپین، بالج و داخلی مجموع دو مقدار زیر می باشد (20):

3)   حداقل انرژی آزاد شده از این نوع حلقه­ها براساس اندازه حلقه در جدول2 مشخص شده است. برای حلقه­های با طول بیشتر از 30، حداقل انرژی آزاد به صورت زیر محاسبه می­شود:

MFE(l) = MFE(30) + 1.75*RT* ln (size/30)

به طوری که  ثابت جهانی گاز،  دمای خالص و  اندازه حلقه است.

جدول 1 - حداقل انرژی آزاد شده از تمام دو جفت­بازهای مجاور در ساقه.

جدول 2- حداقل انرژی آزاد شده از حلقه­های داخلی، بالج و هیرپین براساس اندازه حلقه.

4)   تعدادی از حلقه­ها انرژی مازادی براساس نوکلئوتیدهای حلقه دارند. این نوع حلقه­ها و انرژی مازاد آنها در جدول 3 نشان داده شده است.

توجه کنید که انرژی دیگری بین جفت­باز انتهایی حلقه­ها و دو باز جفت­نشده مجاور آن وجود دارد که در اینجا به آن پرداخته نشده است.

در انتها، افراد براساس مقدار برازندگی حداقل انرژی آزاد به­ترتیب صعودی مرتب می­شوند.

جدول 3- حداقل انرژی آزاد شده از حلقه­ها با توجه به نوکلئوتیدهای حلقه.

عملگر تقاطع: عملگر تقاطع با نرخ 9/0  روی  افراد  انجام می­گیرد، به این صورت که ابتدا 5 درصد از بهترین افراد و 5 درصد از افراد با برازندگی متوسط به نسل بعد منتقل می­شوند. باقی­مانده افراد، به­ترتیب میزان برازندگی دو به دو برای جفت­گیری انتخاب می­گردند. به عبارت دیگر دو فرد  و  از جمعیت به عنوان والدین در نظر گرفته می­شوند، سپس از یک موقعیت تصادفی جفت­گیری می­کنند و در پایان دو فرزند ایجاد می­گردد. با توجه به اینکه طول هر فرد  است، احتمال انتخاب هر موقعیت تصادفی  می­باشد. در این موقعیت، والدین به دو بخش تقسیم شده و فرزند اول و دوم به­ترتیب زیرقطرهای موجود در طرف چپ والدین اول و دوم را به خود اختصاص می دهند. هریک از زیرقطرهای طرف راست والدین اول و دوم در صورتی به­ترتیب در فرزندان دوم و اول قرار می­گیرند که با هیچ یک از زیرقطرهای موجود در فرزند همپوشانی نداشته باشند. اگر همپوشانی وجود داشته باشد، در فرزند قرار داده نمی­شود یا در صورت امکان بخشی که همپوشانی ندارد انتخاب شده و در فرزند قرار می­گیرد. این فرآیند نسل بعد را با فرزندان جدید تشکیل می­دهد.

عملگر جهش: با توجه به اینکه نرخ جهش 1/0 است، 10 درصد از ضعیف­ترین افراد جمعیت گزینش می­گردند تا عملگر جهش روی آنها اجرا شود. برای این افراد، یک زیرقطر تصادفی از ماتریس ­نقطه­ای انتخاب شده و تنها در صورتی با یک زیرقطر تصادفی از فرد جایگزین می­گردد که با هیچ یک از زیرقطرهای موجود در فرد (جز زیرقطر انتخابی از فرد) همپوشانی نداشته باشد. اگر همپوشانی وجود داشته باشد، انتخاب زیرقطر تصادفی از ماتریس ­نقطه­ای ادامه می­یابد تا زمانی که همپوشانی موجود نباشد یا زمان خاتمه یابد.

خاتمه الگوریتم: فرآیند تولید نسل هنگامی متوقف می­شود که شرط  برقرار گردد. در انتها فردی جهت تشکیل ساختار دوم گزینش می­شود که انرژی آزاد کمتری داشته باشد.

نتایج

الگوریتم پیشنهادی، RNAG، روی تعدادی از RNAها جهت پیشگویی ساختار دوم آنها اجرا شده است. مجموعه داده­ها شامل CopA، CopT، R1inv، R2inv، Tar، Tar*، DIS، IncRNA₅₄ و RepZ است. به­عنوان مثال CopA را در نظر بگیرید. شکل1 ساختار دوم پیشگویی شده CopA را نشان می­دهد که به ساختار واقعی بسیار نزدیک است. برای ارزیابی دقت پیشگویی RNAG از دو معیار حساسیت و برجستگی ویژه استفاده می­شود که به صورت زیر محاسبه می­گردند:

=        (1)

=   (2)

معیارF با در نظر گرفتن هر دو مقدار حساسیت و برجستگی ویژه به صورت زیر تعیین می­شود:

معیارF =      (3)

5’-GUGGGCCCCGGUAAUCUUUUCGUACUCGCCAAAGUUGAAGAAGAUUAUCGGGGUUU-3’

                . . . . .((((((((((((( …………………. ))))))))))))). . .

ساختار دوم واقعی CopA

5’-GUGGGCCCCGGUAAUCUUUUCGUACUCGCCAAAGUUGAAGAAGAUUAUCGGGGUUU-3’

                . . . . .(((((.((((((((………………….))))))).))))). . .    

ساختار دوم پیشگویی شده CopA

شکل1

جدول 4 - دقت پیشگویی RNAG روی مجموعه­ای از RNAها.

جدول 5 - مقایسه حساسیت RNAG با تعدادی از رویکردها.

جدول6- مقایسه برجستگی ویژه RNAG با تعدادی از رویکردها.

جدول7- مقایسه معیارF روش RNAG با تعدادی از رویکردها.