نوع مقاله : مقاله پژوهشی

نویسندگان

1 گروه بیوتکنولوژی کشاورزی، دانشگاه تربیت مدرس

2 دانشگاه صنعتی شاهرود

3 Dept. of Plant Breeding and Biotechnolog

چکیده

با پیشرفت‌های سریع در تکنولوژی توالی‌یابی نسل جدید امروزه این تکنیک به ابزاری قدرتمند و کم‌هزینه برای مطالعات در سطح ترنسکریپتوم تبدیل شده است. سرهم‌بندی داده‌های حاصل از توالی‌یابی نسل جدید، به‌صورت de novo باعث شکل‌گیری مسیری نوین در مطالعه و شناخت گونه‌های فاقد ژنوم مرجع گردیده است. با گسترش این تکنولوژی و افزایش روز افزون نرم‌افزارهای سرهم‌بندی، انتخاب مسیر و گزینش نرم‌افزار برتر برای سرهم‌بندی داده‌های حاصل از توالی‌یابی ترنسکریپتوم به عنوان چالشی برای زیست‌شناسان در این زمینه به‌شمار می‌آید. در این پژوهش برای اولین بار داده‌های حاصل از توالی‌یابی ترنسکریپتوم زرین‌گیاه با استفاده از نرم‌افزارهای Oases-velvet، SOAPdenovo-Trans، Trans-ABySS و Trinity به دو صورت مختلف با استفاده از پارامتر K=25 و K=32 به‌منظور دستیابی به مسیر مناسب و نرم‌افزار برتر در این زمینه مورد ارزیابی و آنالیز قرار گرفت. نتایج حاصل از سرهم‌بندی براساس معیارهای متعددی مقایسه شده که گویای برتری Trinity و Trans-ABySSمی‌باشد، در پایان خروجی حاصل از بهترین سرهم‌بندی به منظور بررسی فراوانی ایزوفرم‌های مختلف و آنالیز هستی‌شناسی (Gene Ontology) مورد ارزیابی قرار گرفت. باتوجه به دارویی بودن این گیاه و بالا بودن متابولیت‌های ثانویه آن، بیشترین فراوانی در بخش فرایندهای زیستی، مربوط به فعالیت‌های متابولیتی گزارش شد.

کلیدواژه‌ها

موضوعات

عنوان مقاله [English]

The Comparison of Assembly Softwares and Gene Ontology Analysis using transcriptome sequencing data from Dracocephalum kotschyi Boiss.

نویسندگان [English]

  • Abolnaser Poursalavati 1
  • amin ebrahimi 2
  • Sajad Rashidi Monfared 3

1 Department of agricaltural Biotechnology, Faculty of Agriculture, Tarbiat modares university

2 shahrood university

3 Dept. of Plant Breeding and Biotechnolog

چکیده [English]

With fast advances in next generation sequencing technologies, they has become powerful and low-cost tools for transcriptome studies, Nowadays; de novo assembly of transcriptome data, has caused the formation of the new pathway in the study of non-model genome species. With the expansion of this technology and increasing the number of assembly softwares, choosing the best software for assembling transcriptome sequencing data has become a challenge for biologists. For the first time in this study, we used transcriptome sequencing data of Dracocephalum kotschyi in order to reach the appropriate parameters and superior software; so here we used Oases-velvet, SOAPdenovo-Trans, Trans-ABySS and Trinity softwares with two different values of K parameter; K=25 and K=32. The results of assembly by each software were compared to others in the term of several criteria. The result suggests the superiority of Trinity and Trans-ABySS softwares. Finally, the output of the best assembly was used to estimate abundance of various isoforms and Gene Ontology analysis as regards to the pharmaceutical properties of this plant and the high amount of secondary metabolites, the highest frequency of sections in the biological processes was related to the metabolic activity.

کلیدواژه‌ها [English]

  • Trinity
  • SOAPdenovo-Trans
  • Oases-velvet
  • Trans-ABySS
  • Gene Ontology

مقایسه برنامههای سرهمبندی و آنالیز هستی‌شناسی با استفاده از داده‌های حاصل از توالی یابیترنسکریپتوم زرین‌گیاه(Dracocephalum kotschyi Boiss.)

عبدالناصر پورصلواتی1، امین ابراهیمی2 و سجاد رشیدی‌منفرد1*

1 تهران، دانشگاه تربیت مدرس، گروه بیوتکنولوژی کشاورزی

2 شاهرود، دانشگاه صنعتی شاهرود، گروه زراعت و اصلاح نباتات

تاریخ دریافت: 7/12/96                تاریخ پذیرش: 25/10/97

چکیده

با پیشرفتهای سریع در تکنولوژی توالی‌یابی نسل جدید امروزه این تکنیک به ابزاری قدرتمند و کم‌هزینه برای مطالعات در سطح ترنسکریپتوم تبدیل شده است. سرهمبندی داده‌های حاصل از توالی‌یابی نسل جدید، به صورت de novo باعث شکل‌گیری مسیری نوین در مطالعه و شناخت گونه‌های فاقد ژنوم مرجع گردیده است. با گسترش این تکنولوژی و افزایش روز افزون نرم‌افزارهای سرهمبندی، انتخاب مسیر و گزینش نرم‌افزار برتر برای سرهمبندی داده‌های حاصل از توالی‌یابی ترنسکریپتوم به عنوان چالشی برای زیست‌شناسان در این زمینه به شمار می‌آید. در این پژوهش برای اولین بار داده‌های حاصل از توالی‌یابی ترنسکریپتوم زرین‌گیاه با استفاده از نرم‌افزارهای Oases-velvet، SOAPdenovo-Trans، Trans-ABySS و Trinity به دو صورت مختلف با استفاده از متغیر K=25 و K=32  به منظور دستیابی به مسیر مناسب و نرم‌افزار برتر در این زمینه مورد ارزیابی و آنالیز قرار گرفت. نتایج حاصل از سرهمبندی براساس معیارهای متعددی مقایسه شده که گویای برتری Trinity و Trans-ABySS می‌باشد، در پایان خروجی حاصل از بهترین سرهمبندی به منظور بررسی فراوانی ایزوفرمهای مختلف و آنالیز هستی‌شناسی (Gene Ontology) مورد ارزیابی قرار گرفت. باتوجه به خواص دارویی و مقادیر بالای متابولیتهای ثانویه در این گیاه؛ بیشترین فراوانی مشاهده شده در بخش فرآیندهای زیستی، مربوط به فعالیتهای متابولیتی (Metabolic Process) بود.

واژه های کلیدی: Trinity ,SOAPdenovo-Trans ,Oases-velvet ,Trans-ABySS ,Gene Ontology

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

مقدمه

 

در سالهای اخیر آنالیز ترنسکریپتوم به عنوان یکی از مراحل بنیادی در مطالعات زیستی مطرح شده که به‌این منظور روشهای متعددی از جمله؛ نورترن بلات(Northern blot) ، واکنش زنجیره ای پلیمراز رونوشت برداری معکوس
(RT-PCR)، ریزآرایه‌ها(Microarray)  و توالی‌یابی به روشهای سنتی را می توان نام برد (30) و با‌توجه به پیشرفتهای سریع در توالی‌یابی نسل جدید
(Next-generation sequencing: NGS)، این مورد به ابزاری قدرتمند در آنالیز ترنسکریپتوم تبدیل شده‌است (30 و 46) امروزه با استفاده از تکنولوژی NGS، امکان توالی‌یابی و شناسایی طیف بسیار گسترده‌ای از ژنها در مدتی کوتاه و در میان ژنوم پیچیده و عظیم گیاهی فراهم‌آمده است (40). ازسوی دیگر نرم‌افزارهای متعددی برای سرهمبندی خوانشهای خام (Raw Read) حاصل از توالی‌یابی به وسیله همردیفی خوانشها بر روی ژنوم‌مرجع(Reference assembly)  معرفی شده‌ که عمل شناسایی ژن را در محیط نرم‌افزاری انجام می‌دهند (11)، اما ژنوم‌مرجع برای بسیاری از گونه‌ها از‌جمله گونه مورد مطالعه در این پژوهش )زرین‌گیاه) در دسترس نیست. با این وجود و با بهره‌گیری از نرم‌افزارهای به روز و معرفی تکنیکهای NGS امکان مطالعه در سطوح «امیک»(Omics) برای گونه‌های فاقد نقشه ژنومی نیز فراهم آمده (7)، که در این حالت سرهمبندی de novo، این امکان را ایجاد کرده‌ تا توالی کاملی از ترنسکریپتوم موجود بازسازی و ژنهای بیان شده در یک بافت خاص شناسایی، مقدار سنجی و دسته‌بندی گردد (46).

اولین تلاشها در مطالعات نوین ترنسکریپتوم با استفاده از داده‌های RNA-Seq در سال 2009 به وسیله ونگ و همکاران تحت عنوان ابزاری انقلابی در مطالعه ترنسکریپتوم مطرح شد (40). نرم‌افزارهای مربوط در این‌روش، برای سرهمبندی از دو الگوریتم مختلف De Bruijn graph (45) و Overlap layout-consensus (13) پیروی می‌کنند. در نرم‌افزارهای نسل جدید از الگوریتم De Bruijn graph استفاده شده که در این روش سرهمبندی از طریق شکستن خوانش خام اولیه، به توالیهای کوچک‌تر که K-mer نامیده می‌شوند انجام شده و یافتن همپوشانی میان K-mer ها صورت می‌گیرد (14). نرم‌افزارهای مورد نظر برای این پژوهش شامل (v2.4.0) Trinity (14)، SOAPdenovo-Trans (v1.03) (42)، Oases-Velvet (v0.2.08) (37) و Trans-ABySS (v1.5.5) (34) بوده که همگی از الگوریتم De Bruijn graph پیروی کرده و با‌توجه به گستردگی این دست نرم‌افزارها و پارامترهای مربوط؛ انتخاب نرم‌افزار برتر و روش بهینه برای سرهمبندی داده‌ها امری ضروری است. در سال 2012 مقایسه‌ای میان سه ابزار مختلف سرهمبندی یعنی SOAPdenovo ، Oases و Trinity بر روی خوانشهای حاصل از سیب زمینی شیرین به منظور یافتن پوشش ژنومی قوی‌تر انجام گرفت (38). همچنین بر روی خوانشهای مگس سرکه موجود در دیتابیس، مقایسه میان نرم‌افزارهای سرهمبندی صورت گرفته و نتایج حاصل از سرهمبندی با ژنهای شناخته شده این موجود مقایسه گردید (46). در تحقیقی دیگر برای سرهم بندی داده‌های حاصل از توالی‌یابی با سیستم 454 برای اولین‌بار از ابزارهای مبتنی بر الگوریتم De Bruijn graph استفاده شد (33). با این همه در تحقیقات انجام شده برای مقایسه نرم‌افزارها به طور معمول از داده‌های موجودات مدل استفاده شده و این در حالی‌است که این پژوهش خوانشهای زرین‌گیاه را به صورت de novo و بر اساس مقادیر مختلف K برای سرهمبندی درنظر گرفته است. بایستی به این نکته نیز توجه نمود که الگوریتمهای این برنامه­ها درحال به روز‌رسانی است و مقایسه آخرین ورژن برنامه­ها برای به دست آوردن نتیجه‌ مطلوب ضروری می‌نماید.

زرین‌گیاه با نام علمی Dracocephalum kotschyi یکی از گونه‌های ایندمیک ایران است که در شمال، غرب و مرکز ایران یافت شده (32) و با نام بادرنجبویه‌دنایی نیز شناخته می‌شود (31). با توجه به تعداد کم و خطر انقراض بالقوه،  ضرورت حفاظت، اصلاح و اهلی کردن این گیاه بیش از هر زمانی احساس می‌شود (12 و 20). زرین‌گیاه در طب‌سنتی و داروسازی نیز مورد توجه بوده و تحقیقات متعددی روی آن صورت گرفته‌است (4، 21، 22، 36 و 44).  با این‌ حال، میزان پژوهشهای بیوانفورماتیکی صورت گرفته در زمینه گیاهان ارزشمند دارویی بسیار اندک بوده و به گواه بانک‌ اطلاعاتی NCBI به جز یک توالی rRNA و سه توالی کلروپلاستی، هیچ‌گونه فعالیت و پژوهشی در راستای آنالیزهای مولکولی بر روی زرین‌گیاه صورت نگرفته است (48).

اکنون با استفاده از روشهای نوین با کارآیی بالا و بهره‌گیری از نرم‌افزارهای بیوانفورماتیکی برای تبدیل داده‌های‌خام به اطلاعات مفید و آنالیزهای In silico (37) می توان اقدام به شناسایی هرچه بهتر این گیاه نمود و برای دست‌یابی مؤثر و کارآمد به این هدف، اولین گام انتخاب ابزار و نرم‌افزار مناسب با کارآیی و دقت بالاست. هدف از این پژوهش بررسی مهمترین نرم‌افزارهای موجود در این زمینه و دستیابی به مسیری قابل اعتماد و کارآمد برای سرهمبندی خوانشها به صورت de novo ‌جهت آنالیزهای پایین‌دست می‌باشد. امید است نتایج حاصل از این پژوهش بتواند در موارد مشابه راه‌گشای محققین در تحقیقات آتی قرار گیرد.

مواد و روشها

جمع‌آوری نمونه‌های گیاهی و استخراج mRNA: در این پژوهش، از زرین‌گیاه موجود در ارتفاعات رشته‌کوههای زاگرس در استان لرستان، شهرستان الشتر (با مشخصات، ارتفاع از سطح دریا: 3585 متر، عرض جغرافیایی: 33.955996 و طول جغرافیایی: 48.320011) نمونه‌های برگ جمع‌آوری و در ازت مایع به آزمایشگاه منتقل گردید. مراحل استخراج RNA باکیفیت بالا از 100 میلی‌گرم بافت منجمد شده گیاهی با استفاده از کیت استخراج RNA کیاژن RNeasy Plant Mini Kit (QIAGEN., Cat No.: 74904) انجام گرفت. کیفیت و کمیت RNA استخراج شده توسط ژل آگارز 1 درصد و نانودراپ بررسی شده و مجموعه RNA حاصل از استخراج جهت توالی‌یابی دوطرفه (paired end) با 20میلیون خوانش به طول 150جفت‌باز به وسیله دستگاه Illumina HiSeqTM 2000 برای شرکت توپاز ارسال شد.

بررسی کیفیت نتایج توالی‌یابی، Trimming و Normalization: نتایج حاصل از توالی‌یابی شامل 47 میلیون خوانش به طول 150 جفت‌باز، در قالب فایل فست‌کیو (FASTQ) از طرف شرکت ارسال گردید. این خوانشها Trim شده بودند و توالی آداپتور از ابتدای آنها برش خورده بود. کیفیت داده‌های خام با استفاده از نرم‌افزار FastQC (Version 0.11.5; Simon Andrews) مورد سنجش قرار گرفت و برای رسیدن به صحت و دقت بالاتر در سرهمبندی، توالیهایی با کیفیت کمتر با استفاده از نرم‌افزارهای(v0.36)  Trimmomatic (6) و AfterQC (8) حذف گردید. در ادامه عمل  Normalization بر روی داده‌های خام اعمال گردید که این مرحله با استفاده از ابزار In Silico Normalization که در بسته نرم‌افزاری Trinity قرار دارد صورت گرفت. فرآیند Normalization  براساس فرمول P=min(1 T/C) صورت گرفته (17) که در این پژوهش پارامتر T (Max Target Coverage) در تنظیمات نرم‌افزار برابر با 30 (-max_cov 30) درنظر گرفته شد. در این حالت برای خوانشهای پرتکرار حداقل 30 تکرار از هر خوانش حفظ شده و مابقی حذف خواهد شد (18). همچنین این مقدار توسط Haas و همکاران به عنوان مقدار بهینه در هنگام سرهمبندی خوانشها توصیه شده است (17). 

سرهمبندی خوانشهای حاصل از توالی‌یابی: نتایج حاصل از مراحل قبل به منظور سرهمبندی و تشکیل توالیهای ترنسکریپتوم با استفاده از نرم‌افزارهای (v2.4.0) Trinity، SOAPdenovo-Trans (v1.03)، Oases-Velvet (v0.2.08) و Trans-ABySS (v1.5.5) مورد آنالیز و ارزیابی قرار گرفت. به منظور بررسی تأثیر پارامترهای دخیل در سرهمبندی، این فرآیند با دو مقدار مختلف 25 و 32 برای پارامتر K در تمامی نرم‌افزارها (به جز نرم‌افزار Oases-Velvet که فقط متغیرهای فرد را پذیرفته و مقدار 33 برای آن در نظر گرفته شد) استفاده گردید. براساس مطالعات پیشین؛ بالاترین اندازه برای N50 هنگامی ایجاد می‌شود که بازه عددی 25 تا 35 برای K-mer در نظر گرفته شده و مقادیر 25 و32 در اکثر مطالعات دیگر نیز مورد استفاده قرار گرفته است (5، 9، 41 و 47). همچنین مقدار 25 به عنوان پیش‌فرض، در برخی از نرم‌افزارهای مورد بررسی، استفاده شده و از سوی دیگر در نرم‌افزار Trinity نیز بیشترین عدد مورد پذیرش برای متغیر K-mer، مقدار 32 بود (15). نرم‌افزارهای مورد استفاده در این پژوهش تا زمان نوشتن این مقاله، آخرین نسخه منتشر شده از سوی توسعه‌دهندگان بوده و برای بی‌اثر کردن شرایط محیطی، تمامی آنالیزها در سیستم‌عامل لینوکس و در شرایط مشابه از نظر مشخصات سیستمی (Ram=40GB, CPU=2.4GHz×16, OS=Ubuntu Linux 16.04 LTS) صورت گرفت. در ادامه همه فایلهای خروجی سرهمبندی شده از تمامی نرم‌افزارهای مورد مطالعه، بر‌اساس حداقل طول توالی، با استفاده از ابزار فیلتر طول (fasta_filter_by_min_length) از مجموعه Trinity به ‌مقدار 200 جفت‌باز فیلتر گردید. مقدار 200 جفت باز، با توجه به طول خوانشهای اولیه 150 و مقادیر درنظر گرفته شده برای متغیر K(25 و 32)، می‌تواند مبنای مناسبی برای فیلتر کردن نتایج سرهمبندی باشد. همچنین با توجه به این نکته که در این حالت کوچکترین پروتئین حاصل از سرهمبندی، طولی در حدود 66 آمینواسید خواهد داشت. به‌این ترتیب همه‌ فایلهای خروجی از نظر حداقل طول بازسازی شده، یکسان گردید و از توالیهای کوتاه‌تر و بی‌معنی که به طور معمول حاصل از ضعف الگوریتمها و ابزارهای سرهمبندی بوده، چشم‌پوشی شد. در مرحله بعد فایلهای حاصل از سرهمبندی بر اساس پارامتر‌های تعداد توالی، N50، طول بیشینه و کمینه، میانگین طول توالیهای ایجاد شده، زمان‌ مورد نیاز برای اجرا(RunTime)  بررسی گردید.

همردیفی خوانشهای اولیه روی خروجیهای سرهمبندی: در ادامه خوانشهای خام اولیه (قبل از فرآیند Normalization) با استفاده از نرم‌افزار(v2.3.2)  Bowtie2 (23) بر روی تمامی فایلهای خروجی شامل؛ سرهمبندی با مقادیرK=25  و K=32 همچنین قبل از اعمال فیلتر 200 و بعد از آن، همردیف گردید. به‌این ترتیب نرخ همردیفی (Alignment rate) و میزان پوشش (Coverage values) ایجاد شده میان نرم‌افزارهای مختلف مورد بررسی و مقایسه قرار گرفت.

برآورد میزان فراوانی رونوشتهای حاصل از سرهمبندی: برای بررسی میزان فراوانی رونوشتهای حاصل از سرهمبندی، از ابزار RSEM (24) بهره برده که این فرآیند با استفاده از همردیفی و شمارش تعداد خوانشها برای هر رونوشت با در نظر گرفتن طول رونوشت محاسبه گردید. سپس با استفاده از اسکریپت filter_low_expr_transcripts از مجموعه Trinity، به ازای هر ژن تنها یک ایزوفرم که دارای بیشترین فراوانی بود برای ادامه مسیر انتخاب شد.

مستند‌سازی نتایج و آنالیز GO (Gene Ontology): تمامی رونوشتهای مورد مطالعه، در ابتدا با استفاده از BLASTX  (3) در برابر پایگاه‌ داده گیاه آرابیدوپسیس تالیانا (Arabidopsis thaliana) به شماره اختصاصی (taxid:3702) مورد ارزیابی قرار گرفت. هدف از این کار؛ مستندسازی رونوشتها با حساسیت بالا (1e-5) در برابر اطلاعات فراوان پروتئینی در گیاه آرابیدوپسیس تالیانا به عنوان گیاه مدل می‌باشد که در غیر این صورت امکان بررسی هستی‌شناسی در ادامه فراهم نخواهد بود، زیرا برای انجام تجزیه‌ وتحلیل هستی‌شناسی با استفاده از ابزارهایی مانند PANTHER یا AgriGO تنها بررسی ژنهای شناخته شده و دارای شناسه GO در گیاهان محدودی از جمله آرابیدوپسیس تالیانا، برنج، کاهو، گندم و غیره... امکان‌پذیر بوده که در این میان گیاه آرابیدوپسیس تالیانا بیشترین شناسه GO منحصر به فرد را در میان سایر گیاهان به خود اختصاص داده است. این فرآیند پیش از این، در پژوهشهای مشابه برای مستند‌سازی و هستی‌شناسی رونوشتها صورت پذیرفته است (16 و 29). برای افزایش هرچه بیشتر دقت و اطمینان از نتایج حاصل از این ارزیابی مقدار 5-10 E-value ≤ در نظر گرفته شد. سپس با طراحی اسکریپت مجزایی نتایج حاصل از BLASTX به منظور دست‌یابی و جداسازی معتبرترین پاسخ برای هر رونوشت مورد بررسی قرار گرفت؛ که در طراحی این اسکریپت ابتدا تمامی نتایج حاصل از BLASTX برای هر رونوشت به طور مجزا بر اساس مقدار E-value پایین‌تر و امتیاز بالاتر مرتب گردیده و سپس بهترین نتیجه به ازای هر رونوشت در فایل خروجی نهایی ذخیره و ثبت گردید. در آخر با استفاده از شماره اختصاصی مرتبط به هریک از نتایج، فرآیند هستی‌شناسی ‌ژنها (Gene Ontology)در ابزار آنلاین PANTHER  (27) انجام شد.

نتایج

پس از دریافت خوانشهای حاصل از توالی‌یابی دوطرفه با استفاده از دستگاه Illumina HiSeqTM 2000 ، کیفیت نتایج

توالی‌یابی توسط FastQC بررسی گردید (شکل 1) و با توجه به کیفیت قابل قبول خوانشها، سه نوکلئوتید از انتهای ′3 برخی خوانشها با استفاده از ابزار  Trimmomaticو AfterQC حذف شد.

 

 

شکل 1- کیفیت خوانشها برای هر توالی (a)، کیفیت خوانشها برای هر نوکلئوتید )امتیاز (Phred (b)، میزان توالیهای تکراری در خوانشها (c)، مقدار GC در توالی خوانشها (d)

 

در ادامه با‌توجه به نتایج FastQC و وجود خوانشهای تکراری فراوان در داده‌های خام؛ به منظور سهولت و سرعت بالاتر در فرآیند سرهمبندی و از سوی دیگر به علت محدودیت در منابع سیستمی برای آنالیز داده‌ها، لازم است فرآیند Normalization بر روی خوانشها اعمال شود. به این ترتیب برای فایل خام اولیه که دارای 23,599,495 خوانش به طول 3,539,924,250 نوکلئوتید با درصد GC 48 و حجم 5/8 گیگابایت می باشد؛ بدون از دست دادن خوانشها با تکرار کمتر از 30، خوانشهای تکراری از فایل اولیه حذف شده و نتیجه حاصل از آن، ایجاد فایل نرمالایز شده نهایی با 6,529,853 خوانش به طول 979,477,950 نوکلئوتید و درصد GC 45 به حجم 4/2 گیگابایت بود )تمامی فایلهای ذکر شده در قالب FASTQ قرار دارند(.

سپس فرآیند سرهمبندی با استفاده از خوانشهای نرمالایز شده از مرحله قبل و نرم‌افزارهای ذکر شده صورت گرفت. تمامی نرم‌افزارهای مورد استفاده به جز Trinity همگی براساس نسخه ژنومی آنها ایجاد شده‌اند، به‌این صورت که SOAPdenovo-Trans تحت چهارچوب (Framework) SOAPdenovo، Trans-ABySS در چهارچوب ABySS و Oases با استفاده از نسخه ژنومی آن یعنی Velvet نوشته شده، در حالی که Trinity به طور اختصاصی برای سرهمبندی داده‌های حاصل از توالی‌یابی ترنسکریپتوم توسعه داده شده است. از سوی دیگر درحالی‌که تمامی این نرم‌افزارها از فرآیند ایجاد K-mer و الگوریتم De Bruijn graph برای سرهمبندی خوانشها استفاده می‌کنند، ولی با‌توجه به گزینه‌ها و متغیرهای خود، فرآیند متفاوتی برای ایجاد فایل رونوشت نهایی و سرهمبندی K-mer ها درنظر می‌گیرند.  در Trinity تنها یک فایل سرهمبندی به عنوان خروجی ایجاد شده و شامل ایزوفرمها و یونی ‌ژنها می‌باشد؛ این درحالی است که سایر نرم‌افزارها علاوه بر ایجاد فایل کانتیگ (Contig)‌ اولیه؛ قادر به ایجاد سطح بالاتری از سرهمبندی به نام اسکفولد (Scafold) براساس الگوریتم Overlap layout-consensus بوده که اگرچه ممکن است براساس نواحی همپوشان با اختصاصیت پایین تشکیل شده باشد ولی در مراحل بعد و آنالیزهای پایین‌دست می تواند مورد استفاده قرار گیرد.

در SOAPdenovo-Trans فایلهای خروجی شامل Contig و ScafSeq و نرم‌افزار Oases-Velvet نیز دارای دو فایل خروجی با نامهای Contigs و Transcripts بود. اما Trans-AbySS سه فایل با نامهای  Transabyss.jn، Transabyss.ref و  Transabyss.final تولید نمود. با‌توجه به حداقل طول‌ توالیهای مختلف در فایلهای خروجی هریک از نرم‌افزارها و به منظور بی‌اثر کردن تأثیر توالیهای کوتاه در این پژوهش، فایلهای حاصل از سرهمبندی برای ادامه کار تحت تأثیر فیلتر طول به میزان 200 جفت‌باز قرار گرفت. در نرم‌افزار Trinity، حداقل طول توالی بازسازی شده به طور پیش‌فرض 200 جفت‌باز بوده و به این ترتیب خروجیهای سایر نرم‌افزارها با استفاده از ابزار فیلتر طول، به میزان 200 نوکلئوتید فیلتر گردید.

در ادامه برای رسیدن به میزان صحت و پوشش خوانشها بر روی فایلهای خروجی، خوانشهای اولیه (قبل از نرمالایز شدن) بر روی تمامی فایلهای سرهمبندی شده حاصل از نرم‌افزارهای مختلف توسط Bowtie2 همردیف شده و نرخ همردیفی برای هر نرم‌افزار در هر فایل خروجی بررسی گردید (39). خروجیهای مختلف تمام نرم‌افزارها؛ قبل و بعد از فیلتر 200 که در مجموع شامل 30 فایل سرهمبندی مختلف می‌باشد در جدول 1 باهم مقایسه شده است. به این ترتیب با در نظر گرفتن K=25، متغیر N50 که بیانگر حداقل طول 50 درصد از توالیهای سرهمبندی شده می‌باشد (28) برای خروجی transcripts از نرم‌افزار Oases-Velvet ( bp1976 (، خروجی‌ transabyss.ref از Trans-ABySS (bp1193( و همچنین خروجی Trinity ( bp1585 (در مقدار قابل قبولی قرار دارند و با‌توجه به تعداد توالیهای ایجاد شده، خروجی Trinity، فایل transcripts از نرم‌افزار Oases-Velvet در وضعیت بهتری قرار دارند. همچنین در حالتی که متغیر K=32 در نظر گرفته شد، فایل transcripts از نرم‌افزار Oases-Velvet ( bp2147 (، خروجی transabyss.ref از نرم‌افزار Trans-ABySS ( bp1312 ( و خروجی Trinity ( bp1647) بیشترین مقدار N50 را به خود اختصاص داده (جدول 1) و در مطالعات قبلی که روی Trinity ، SOAPdenovo و Trans-ABySS صورت گرفت نیز نتایج مشابهی به دست آمد؛ همچنین در اثر افزایش مقدار K-mer، میزان N50 به طور محسوس افزایش یافته که در این پژوهش نیز نتایج به دست آمده بیانگر این مورد است (9 و 25).

 بر اساس تعداد توالی پیش‌‌بینی شده با‌توجه به میزان N50 هریک از خروجیها، در هر دو حالت K=32 و  K=25، فایل  transabyss-final تعداد بیشتری توالی را بازسازی نموده (K=25 شامل 390072 توالی و  K=32شامل 286004 توالی) که در مطالعه Wang و Gribskov نیز این مورد بالاترین تعداد (1070887) را ایجاد نمود (41). پارامتر طول توالیهای بازسازی شده در این پژوهش در حالت K=25 نشان‌دهنده برتری Oases-Velvet در فایل transcripts با میانگین طول bp 1051 و در حالت K=32 نیز با میانگین bp 1363 دارای برتری بوده ولی با‌توجه به فواصل خالی (gaps) (43) در توالیهای تشکیل شده توسط Oases-Velvet، این امر می‌تواند امتیازی منفی به شمار آید. در این‌صورت فایل خروجی Trinity با میانگین طول bp 1036 در K=25 و  طول bp 1066 در K=32، از این نظر موفقیت بیشتری به همراه داشته (جدول 1) و این مورد در پژوهشهای پیشین نیز مورد توجه قرار گرفته است (Zhao و همکاران میانگین طول  bp604 و Wang و Gribskov نیز برتری Oases-velvet را برای میانگین طول بهتر به میزان  bp1593 برای K=25 و bp 1710 برای K=31  گزارش کردند) (41 و 47).  بر‌اساس میزان منابع سیستمی، از جمله مقدار رم و زمان مورد نیاز برای اجرا، SOAPdenovo-Trans کمترین مقدار رم و سریع‌ترین زمان (35 دقیقه) برای سرهمبندی را به خود اختصاص داده و در مقابل، Trinity علاوه بر مقدار بالاتر رم (یک گیگابایت به ازای هر یک میلیون خوانش (19))، زمان بالاتری (15 ساعت) نیز برای سرهمبندی خوانشها صرف می‌کند. البته باید به این نکته توجه نمود که کارآیی آنالیز در هر نرم‌افزار برای سرهمبندی خوانشها، ارتباط مستقیمی با دقت و صحت سرهمبندی ندارد (47).

همردیف کردن خوانشهای خام اولیه بر روی خروجی هر نرم‌افزار موجب شکل‌گیری مبنای مناسبی برای مقایسه میان نتایج سرهمبندی گردیده که براساس نتایج جدول 1 با افزایش مقدار K نرخ همردیفی به طور محسوس افزایش یافته و همه‌ی نرم‌افزارها به غیر از SOAPdenovo-Trans افزایش قابل توجهی را در مقدار همردیفی و تعداد خوانشهای همردیف شده نشان می‌دهند، به‌این صورت که با در نظر گرفتن K=25 بیشترین همردیفی مربوط به transabyss.final (87/70 درصد) و خروجی Trinity (71/67 درصد) همچنین در حالت K=32 نیز خروجی transabyss.final (49/75 درصد) و Trinity (98/72 درصد) بیشترین مقادیر را به خود اختصاص داده (جدول 1) که در مطالعه سال 2012 روی دروزوفیلا نتایج مشابهی برای پوشش ژنومی خوانشها در برابر ژنوم مرجع موجود در پایگاه‌داده حاصل آمد (Trinity با 6/78 درصد و Trans-ABySS با 3/64 درصد) (47). همچنین در مطالعه اخیر بر روی نرم‌افزارهای مطرح در این زمینه و سرهمبندی خوانشهای آرابیدوپسیس و همردیفی آنها روی ژنوم، نتیجه مشابهی به دست آمد (Trans-ABySS 44/93 درصد و Trinity 21/90 درصد) (41).

همچنین با افزایش مقدار K-mer، تعداد توالیهای تشکیل شده در نرم‌افزار Trinity از 165125 به 165697و در نرم‌افزار SOAPdenovo-Trans (فایل scafseq) از 68979 به 76223 عدد افزایش یافته و در دو نرم‌افزار دیگر یعنی Oases-Velvet (فایل transcripts) از 124813 به 100048و در نرم‌افزار Trans-ABySS‌ (فایل transabyss.final) از 390072 به 286004 عدد کاهش یافته که این امر می‌تواند با الگوریتمهای متفاوت این نرم‌افزار‌ها در سرهمبندی در ارتباط باشد. از سوی دیگر طول توالیهای بازسازی شده با افزایش مقدار K-mer به طور محسوس افزایش داشته و تنها در مورد SOAPdenovo-Trans با افزایش مقدار K-mer، میانگین طول توالیها از bp 441 به bp 405 کاهش یافته و با‌توجه به کوتاه شدن توالیهای حاصل از سرهمبندی، این امر می‌تواند علت کاهش نرخ همردیفی در وضعیت K=32 در خروجی scafseq را توجیه کند. ایجاد ترنسکریپت بلندتر و بدون فاصله، نقش مهمی در تشکیل نواحی همپوشان برای آنالیزهای پایین‌دست به منظور شناسایی و پیش‌بینیORF ها (Open Reading Frame) و ژنهای این گیاه دارد (جدول 1).

در مطالعه‌ای که با استفاده از خوانشهای شبیه‌سازی شده در کنار خوانشهای واقعی روی کروموزوم شماره‌22 انسان جهت بررسی کارآیی چند نرم‌افزار سرهمبندی از جمله Trinity، Oases و ABySS صورت گرفت، نتایج از برتری Trinity  در زمینه نرخ پوشش ژنومی و پارامترهای آماری مرتبط از جمله N50 حکایت داشته و از طرفی در مقایسه میان دو ابزار دیگر یعنی Oases و ABySS، به برتری نرم‌افزار Oases اشاره شده که علت این امر می‌تواند با به روزرسانیهای متعدد از زمان انتشار مقاله مذکور در ارتباط باشد (10).

 

جدول 1- مقایسه نرم‌افزارهای Trinity،  SOAPdenovo-Trans، Oases-Velvet و Trans-ABySS در سرهمبندی خوانشهای حاصل از توالی‌یابی ترنسکریپتوم زرین‌گیاه

 

Oases-Velvet

SOAPdenovo-Trans

Trans-ABySS

Trinity

File names

contigs

Transcripts

contig

scafSeq

Transabyss.jn

Transabyss.ref

Transabyss.final

Trinity

K=25

N50 bp

264

(98)

2028

(1976)

606

(355)

1273

(999)

654

(642)

1302

(1193)

955

(676)

1585

Sequence number

99102

(1354622)

94910

(124813)

89318

(304772)

68979

(144075)

55382

(60063)

57038

(90319)

112242

(390072)

165125

Average Length bp

9/282

(9/91)

8/1336

(1051)

8/505

(2/219)

778

(2/441)

5671

(536)

7/833

(1/536)

8/684

(248)

88/1036

Max length bp

Min length bp

2400

200 (49)

15429

200 (100)

8552

200 (26)

15035

200 (100)

5276

200 (48)

14730

200 (26)

14730

200 (25)

13250

201

Alignment Rate %

72/0

(79/0)

41/58

(43/58)

55/39

(59/39)

45/44

(48/44)

36/15

(38/15)

56/67

(58/67)

87/70

(91/70)

71/67

Mapped Reads

170815

(185459)

13784614

(13790068)

9332805

(9344184)

10490889

(10496146)

3624599

(3628502)

15944148

(15949403)

16725570

(16733450)

15979185

RunTime

00:05'

(00:13')

00:58'

(2:20')

00:22'

(00:29')

00:59'

(00:29')

00:39'

(00:41')

00:32'

(00:33')

00:39'

(00:45')

1:50'

K=32

N50 bp

264

(123)

2170

(2147)

581

(332)

1288

(955)

791

(784)

1365

(1312)

1050

(830)

1647

Sequence number

114503

(1010046)

89623

(100048)

95826

(305136)

76223

(182956)

55597

(57949)

59322

(72971)

115274

(286004)

165597

Average Length bp

6/268

(117)

8/1501

(1363)

6/490

(2/232)

1/791

(5/405)

6/662

(7/642)

7/881

(8/739)

2/741

(8/353)

4/1066

Max length bp

Min length bp

2562

200 (65)

15857

120 (200)

6858

200 (34)

15075

200 (100)

4877

200 (73)

16150

200 (33)

16150

200 (32)

12331

201

Alignment Rate %

05/1

(14/1)

92/65

(94/65)

17/40

(22/40)

53/43

(56/43)

59/21

(60/21)

82/71

(83/71)

59/75

(52/75)

98/72

Mapped Reads

248768

(269193)

15557224

(15561648)

9479314

(9491360)

10273837

(10280241)

5094996

(5097866)

16949044

(16951093)

17814817

(17822511)

17223876

RunTime

00:25'

(00:36')

1:28'

(1:20')

00:29'

(00:48')

00:32'

(00:36')

00:11'

(00:16')

00:50'

(00:56')

00:48'

(00:50')

2:08'

* مقادیر داخل پرانتز بدون فیلتر 200 برای حداقل طول توالی محاسبه شده‌اند.

 

براساس مشاهدات این تحقیق، Trinity و Oases-Velvet در K-merهای بزرگ‌تر عملکرد مناسب‌تری داشته و همان گونه که در جدول 1 مشاهده می‌شود با افزایش مقدار K، میزان پوشش همردیفی بالاتری برای خروجیهای این دو نرم‌افزار ایجاد گردیده‌است. به این‌صورت که Trinity از 71/67 درصد به 98/72 درصد و Oases-Velvet از 43/58 درصد به 94/65 درصد در K=32 رسیده است و این درحالی است که SOAPdenovo-Trans درحالت K=25 نتایج مناسب‌تری تولید نمود (48/44 درصد و در K=32 نرخ همردیفی به 56/43 درصد کاهش یافت).

ابزار Trinity، نرم‌افزاری رایگان و متن‌باز بوده و از سه ماژول مختلف Butterfly) و  (Inchworm, Chrisalis برای انجام فرآیند سرهمبندی استفاده کرده که این عمل از طریق: سرهمبندی خوانشهای توالی‌یابی به صورت رونوشتهای اولیه در Inchworm، دسته‌بندی این رونوشتها و تشکیل گراف دی‌بروین در Chrisalisو درنهایت پردازش گراف به منظور گزارش تمام رونوشتها با طول کامل و ایزوفرمهای حاصل از پیرایش ثانویه در Butterfly صورت می‌گیرد. فایل نهایی حاصل از سرهمبندی در قالب FASTA (.fa) و با اندازه تقریبی 240 مگابایت؛ شامل 67859 ژن و 165597 رونوشت با درصد GC 64/42 بود. میزان N50 نیز برای رونوشتها و ژنها به ترتیب برابر با 1647 و 1353 بوده، همچنین اطلاعات آماری بیشتر در ارتباط با ترنسکریپتوم سرهمبندی شده در جدول 2 قابل مشاهده است که با توجه به نرخ همردیفی (Mapping) 98/72 (جدول 1) و اطلاعات آماری مناسب (جدول 2)، خروجی Trinity (K-mer=32) در وضعیت خوبی قرار داشته و نقطه شروع مناسبی برای آنالیزهای پایین دست می‌باشد. در ادامه  با استفاده از ابزار RSEM فراوانی هریک از رونوشتها بر مبنای دو پارامتر TPM و FPKM محاسبه شد که نتایج این بررسی در فایل S1 ضمیمه شده است.

جدول 2- اطلاعات آماری فایل خروجی سرهمبندی با ابزار Trinity

تمام ژنهای Trinity

67859

تمام رونوشتهای Trinity

165597

درصد GC

64/42

اطلاعات آماری بر اساس تمام رونوشتها

Contig N10

 bp3325

Contig N20

bp 2621

Contig N30

bp 2188

Contig N40

bp 1867

Contig N50

bp 1647

طول میانه

bp 742

طول متوسط

bp 88/1036

تعداد بازهای سرهمبندی شده

171215531

اطلاعات آماری بر اساس بلندترین ایزوفرم در هر ژن

Contig N10

bp 3275

Contig N20

bp 2521

Contig N30

bp 2065

Contig N40

bp 1693

Contig N50

bp 1353

طول میانه

bp 446

طول متوسط

bp 62/795

تعداد بازهای سرهمبندی شده

53989972

 

سپس به منظور جداسازی رونوشت دارای بیشترین بیان به ازای هر ژن، از ابزار filter_low_expr_transcripts استفاده کرده که در نتیجه فایل fasta با حجم 65 مگابایت و دارای 67859 توالی منحصر به فرد با میزان N50 به طول bp 1021، میانگین bp 683 و حداکثر طول  bp11820 برای انجام آنالیز BLASTX ایجاد شد. پس از انجام BLASTX، تعداد 454399 خروجی برای 24887 رونوشت از مجموع 67859 عدد رونوشت اولیه توسط فرآیند BLASTX از میان پروتئینهای گیاه آرابیدوپسیس تالیانا استخراج شد. سپس با استفاده از اسکریپت طراحی شده برای این پژوهش، معتبرترین نتایج براساس میزان امتیاز بالاتر و مقدار 5-10 E-value ≤ به ازای هر رونوشت جدا سازی گردید که شماره اختصاصی این پروتئینها در پایگاه داده آرابیدوپسیس تالیانا، در فایل S2 ضمیمه شده است. بدین‌ ترتیب فرآیند هستی‌شناسی با استفاده از 24887 نتیجه نهایی در ابزار PANTHER صورت گرفته و 7684 عدد پاسخ GO منحصر به فرد ثبت گردید که نتایج این بررسی در شکل 2 مشاهده می‌شود. خروجی نهایی بررسی GO در فایل S3 قرار دارد.

 

 

 

 

اجزای تشکیل‌دهنده سلول(Cellular Component)

 

فعالیتهای مولکولی (Mollecular Function)

 

فرآیندهای زیستی(Biological Process)

 

دسته‌بندی پروتئینها

شکل2- نتایج هستی‌شناسی (GO) 24888 عدد از رونوشتهای سرهمبندی شده در سیستم PANTHER: اجزای تشکیل‌دهنده سلول (Cellular Component)، فعالیتهای مولکولی (Mollecular Function)، فرآیندهای زیستی (Biological Process) و نتایج حاصل از دسته‌بندی پروتئینهای شناسایی شده


بحث

براساس نتایج این پژوهش خروجی scafSeq از نرم‌افزار SOAPdenovo-Trans و خروجی transcripts از نرم‌افزار Oases-Velvet همچنین در نرم‌افزار Trans-ABySS دو فایل خروجی Transabyss.ref و Transabyss.final با توجه به پوشش بهتر و میزان N50 بالاتر نسبت به سایر خروجیها در وضعیت مناسب‌تری قرار دارند. از سوی دیگر با افزایش میزان K تغییرات محسوسی در نرخ همردیفی و پوشش خوانشها ایجاد شده که K=32 می تواند مقدار مناسبی برای انجام این‌دست پژوهشها به شمار آید. در نهایت می توان خروجی نرم‌افزار Trinity و همچنین Transabyss.final از نرم‌افزار  Trans-ABySS را بهترین گزینه برای انجام آنالیزهای پایین‌دست معرفی نمود. که باتوجه به سرعت بالاتر و زمان کوتاه‌تر در نرم‌افزار Trans-ABySS، این مورد می تواند امتیاز دیگری برای آن محسوب شود. با‌ این حال خروجی Trinity تنها موردی است که علاوه بر تشکیل ترنسکریپتها، می‌تواند توالیهای ایجاد شده و مشابه به هر یونی ‌ژن را تحت عنوان یک ایزوفرم دسته‌بندی و مشخص نماید که این عمل تنها مخصوص به این نرم‌افزار بوده و برای بررسی میزان بیان ژنها و رونوشتها و همچنین تغییرات ژنتیکی در محیط نرم‌افزاری مفید باشد. در ادامه فرآیند هستی‌شناسی در میان رونوشتهای دارای بیشترین فراوانی نشان دهنده فعالیت بالای کاتالیتیک و حضور بسیار بالای پروتئینهای دخیل در اتصال و همچنین پروتئینهای دارای فعالیتهای آنتی اکسیدانی در میان فعالیتهای مولکولی این گیاه بوده و باتوجه به دارویی بودن این گیاه و اهمیت متابولیتهای ثانویه در آن؛ از جمله متوکسی فلاون‌ها و رزمارینیک اسید (1)، بیشترین فراوانی در بخش فرآیندهای زیستی با بیش از 2500 ژن (شکل 2)، مربوط به فرآیندهای متابولیتی در این گیاه دارویی می‌باشد که با یافته‌های مشابه در گونه‌های نزدیک مطابقت دارد (2 و 26). کلاس‌بندی پروتئینهای شناسایی شده نیز حاکی از بیان بالای پروتئینهای ترنسفراز و اکسیدو ردوکتازها بوده که با توجه به استخراج نمونه اولیه از برگ گیاه توجیه می‌گردد. با این همه؛ توسعه و گسترش نرم‌افزارها و الگوریتمهای سرهمبندی برای رسیدن به مسیر و ابزاری با قابلیت بسیار بالاتر و دقیق‌تر همچنان به عنوان یک چالش مهم در دنیای بیوانفورماتیک مطرح  بوده و از سوی دیگر برای بررسی دقیق ژنها و پروفایل بیان این گیاه نیز، مستندسازی تمامی رونوشتها و بررسی میزان بیان برخی از فراوان‌ترین آنها توسط روشهای معمول آزمایشگاهی و مقایسه با روشهای نرم‌افزاری می‌تواند در یافتن مسیرهای متابولیتی با اهداف افزایش بیان متابولیت سودمند، راهگشای محققین و علاقه‌مندان قرار گیرد.

1-      ایوبی، ن.، حسینی، ب. و فتاحی، م. (1396). اثر القایی کیتوزان و کلشی‌سین بر تولید رزمارینیک اسید در ریشه مویین زرین گیاه  .(Dracocephalum kotschyi Boiss)مجله پژوهشهای سلولی و مولکولی، 30(1)، 1-13.

2-      میرزایی، س. (1396). پروفایل ترانسکریپتوم نوک شاخه و ریشه گیاه سویا. مجله پژوهشهای سلولی و مولکولی، 30(4)، 499-511.

 

3-      Altschul SF, Madden TL, Schäffer AA, Zhang J, Zhang Z, Miller W, Lipman DJ (1997) Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res 25(17):3389–3402

4-      Amirghofran Z, Azadbakht M, Karimi MH (2000) Evaluation of the immunomodulatory effects of five herbal plants. J Ethnopharmacol 72(1):167–172

5-      Blande D, Halimaa P, Tervahauta AI, Aarts MGM, Kärenlampi SO (2017) De novo transcriptome assemblies of four accessions of the metal hyperaccumulator plant Noccaea caerulescens. Sci Data 4:160131

6-      Bolger AM, Lohse M, Usadel B (2014) Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30(15):2114–2120

7-      Bräutigam A, Mullick T, Schliesky S, Weber APM (2011) Critical assessment of assembly strategies for non-model species mRNA-Seq data and application of next-generation sequencing to the comparison of C3 and C4 species. J Exp Bot 62(9):3093–3102

8-      Chen S, Huang T, Zhou Y, Han Y, Xu M, Gu J (2017) AfterQC: automatic filtering, trimming, error removing and quality control for fastq data. BMC Bioinformatics 18(Suppl 3):80

9-      Chopra R, Burow G, Farmer A, Mudge J, Simpson CE, Burow MD (2014) Comparisons of de novo transcriptome assemblers in diploid and polyploid species using peanut (Arachis spp.) RNA-seq data. PLoS One 9(12):e115055

10-   Clarke K, Yang Y, Marsh R, Xie L, Zhang KK (2013) Comparative analysis of de novo transcriptome assembly. Sci China Life Sci 56(2):156

11-   Duan J, Xia C, Zhao G, Jia J, Kong X (2012) Optimizing de novo common wheat transcriptome assembly using short-read RNA-Seq data. BMC Genomics 13(1):392

12-   Fattahi M, Nazeri V, Torras-Claveria L, Sefidkon F, Cusido RM, Zamani Z, Palazon J (2013) Identification and quantification of leaf surface flavonoids in wild-growing populations of Dracocephalum kotschyi by LC–DAD–ESI-MS. Food Chem 141(1):139–146

13-   Flicek P, Birney E (2009) Sense from sequence reads: methods for alignment and assembly. Nat Methods 6:S6–S12

14-   Grabherr MG, Haas BJ, Yassour M, et al (2011) Trinity: recontructing a full-length transcriptome assembly without a genome from RNA-Seq data. Nat Biotechnol 29(7):644–652

15-   Grabherr M, Haas BJ, Yassour M, et al (2011) Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol 29(7):644–52

16-   Guo L, Allen KS, Deiulio G, Zhang Y, Madeiras AM, Wick RL, Ma L (2016) A De Novo -Assembly Based Data Analysis Pipeline for Plant Obligate Parasite Metatranscriptomic Studies. 7(July):1–9

17-   Haas BJ, Papanicolaou A, Yassour M, et al (2013) De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat Protoc 8(8):1494–1512

18-   Haas BJ, Papanicolaou A, Yassour M, et al (2013) De novo transcript sequence reconstruction from RNA-Seq: reference generation and analysis with Trinity. Nat Protoc 8(8):10.1038/nprot.2013.084

19-   Haas BJ, Papanicolaou A, Yassour M, et al (2013) De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat Protoc 8(8):1494–1512

20-   Jalili A, Jamzad Z (1999) Red data book of Iran: A preliminary survey of endemic, rare and endangered plant species in Iran.

21-   Javidnia K, Miri R, Kamalinejad M, Khoshneviszadeh M (2006) Constituents of the volatile oils of Dracocephalum kotschyi Boiss. from Iran. J Essent Oil Res 18(3):342–344

22-   Kamali M, Khosroyar S, Jalilvand MR (2014) Evaluation of phenolic, flavonoids, anthocyanin contents and antioxidant capacities of different extracts of aerial parts of Dracocephalum kotschyi. JN Khorasan Univ Med Sci 6:627–634

23-   Langmead B, Salzberg SL (2012) Fast gapped-read alignment with Bowtie 2. Nat Methods 9(4):357–359

24-   Li B, Dewey CN (2011) RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 12(1):323

25-   Li G, Du J, Li Y, Wu J (2015) Identification of putative olfactory genes from the oriental fruit moth grapholita molesta via an antennal transcriptome analysis. PLoS One 10(11):1–30

26-   Li H, Fu Y, Sun H, Zhang Y, Lan X (2017) Transcriptomic analyses reveal biosynthetic genes related to rosmarinic acid in Dracocephalum tanguticum. Sci Rep 7(1):74

27-   Mi H, Muruganujan A, Casagrande JT, Thomas PD (2013) Large-scale gene function analysis with the panther classification system. Nat Protoc 8(8):1551–1566

28-   Miller JR, Koren S, Sutton G (2010) Assembly Algorithms for Next-Generation Sequencing Data. Genomics 95(6):315–327

29-   Mizrachi E, Hefer CA, Ranik M, Joubert F, Myburg AA (2010) De novo assembled expressed gene catalog of a fast-growing Eucalyptus tree produced by Illumina mRNA-Seq.

30-   Morozova O, Marra MA (2008) Applications of next-generation sequencing technologies in functional genomics. Genomics 92(5):255–264

31-   Mozaffarian V (1996) A dictionary of Iranian plant names: Latin, English, Persian. Farhang Mo’aser

32-   Rechinger KH (1982) Salvia in Flora Iranica, Labiatae, No. 150, edits. Rechinger KH Hedge IC, Akad. Druck Vertagsanstalt, Graz, Austria. , p 477

33-   Ren X, Liu T, Dong J, Sun L, Yang J, Zhu Y, Jin Q (2012) Evaluating de Bruijn Graph Assemblers on 454 Transcriptomic Data. PLoS One. doi: 10.1371/journal.pone.0051188

34-   Robertson G, Schein J, Chiu R, et al (2010) De novo assembly and analysis of RNA-seq data. Nat Meth 7(11):909–912

35-   Schulz MH, Zerbino DR, Vingron M, Birney E (2012) Oases: robust de novo RNA-seq assembly across the dynamic range of expression levels. Bioinformatics 28(8):1086–1092

36-   Sharafi A, Sohi HH, Azadi P, Sharafi AA (2014) Hairy root induction and plant regeneration of medicinal plant Dracocephalum kotschyi. Physiol Mol Biol Plants 20(2):257–262

37-   Tang Q, Ma X, Mo C, Wilson IW, Song C, Zhao H, Yang Y, Fu W, Qiu D (2011) An efficient approach to finding Siraitia grosvenorii triterpene biosynthetic genes by RNA-seq and digital gene expression analysis. BMC Genomics 12(1):343

38-   Tao X, Gu Y-H, Wang H-Y, Zheng W, Li X, Zhao C-W, Zhang Y-Z (2012) Digital gene expression analysis based on integrated de novo transcriptome assembly of sweet potato [Ipomoea batatas (L.) Lam.]. PLoS One 7(4):e36234

39-   Trapnell C, Salzberg SL (2009) How to map billions of short reads onto genomes. Nat Biotechnol 27(5):455–457

40-   Wang Z, Gerstein M, Snyder M (2009) RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet 10(1):57–63

41-   Wang S, Gribskov M (2016) Comprehensive evaluation of de novo transcriptome assembly programs and their effects on differential gene expression analysis. Bioinformatics 215:btw625

42-   Xie Y, Wu G, Tang J, et al (2014) SOAPdenovo-Trans: de novo transcriptome assembly with short RNA-Seq reads. Bioinformatics 30(12):1660–1666

43-   Xu D-L, Long H, Liang J-J, Zhang J, Chen X, Li J-L, Pan Z-F, Deng G-B, Yu M-Q (2012) De novo assembly and characterization of the root transcriptome of Aegilops variabilis during an interaction with the cereal cyst nematode. BMC Genomics 13(1):133

44-   Yaghmai MS, Taffazoli R (1988) The essential oil of Dracocephalum kotschyi Boiss. Flavour Fragr J 3(1):33–36

45-   Zerbino DR, Birney E (2008) Velvet: algorithms for de novo short read assembly using de Bruijn graphs. Genome Res 18(5):821–829

46-   Zhao Q-Y, Wang Y, Kong Y-M, Luo D, Li X, Hao P (2011) Optimizing de novo transcriptome assembly from short-read RNA-Seq data: a comparative study. BMC Bioinformatics 12(14):S2

47-   ZHAO L, Zachary L-R, CHEN S-Y, GUO Z-H (2012) comparing De Novo Transcriptome Assemblers Using Illumina RNA-Seq Reads. Plant Divers Resour 34(5):487

48-   NCBI database search. https://www.ncbi.nlm.nih.gov/gquery/?term=Dracocephalum+kotschyi.