دوره جامع و فوق‌تخصصی بیوانفورماتیک: آنالیز داده‌های توالی‌یابی RNA (RNA-seq) با R و Bioconductor

تکنولوژی توالی‌یابی کل ترنسکریپتوم (Whole Transcriptome Sequencing) یا RNA-seq، با کنار زدن محدودیت‌های روش‌های قدیمی مبتنی بر هیبریداسیون (مانند میکروارری)، به استاندارد طلایی در ژنومیک عملکردی تبدیل شده است. این تکنولوژی با ارائه رزولوشن تک‌نوکلئوتیدی، دامنه دینامیک (Dynamic Range) فوق‌العاده وسیع و قابلیت شناسایی ایزوفرم‌های جدید (Novel Isoforms) و RNA‌های غیرکدکننده (lncRNA)، دیدگاه ما را نسبت به پیچیدگی‌های بیان ژن متحول کرده است. با این حال، تبدیل ترابایت‌ها داده خام (Raw Reads) به بینش بیولوژیکی معتبر، نیازمند تسلط بر پایپ‌لاین‌های محاسباتی چندمرحله‌ای، درک عمیق از مدل‌های آماری و مهارت در زبان برنامه‌نویسی R است. این مستند فنی، به عنوان یک مرجع آکادمیک و عملیاتی، تمامی مراحل این فرآیند پیچیده را از کنترل کیفیت اولیه تا تحلیل‌های سیستمی شبکه ژنی پوشش می‌دهد.

فصل اول: معماری داده‌های NGS و پیش‌پردازش (Upstream Analysis)

پیش از ورود به محیط R، داده‌های خام تولید شده توسط دستگاه‌های سکوئنسر (مانند Illumina NovaSeq) باید پالایش شوند. کیفیت ورودی مستقیماً تعیین‌کننده صحت خروجی است.

۱. آنالیز فایل‌های FASTQ و کنترل کیفیت (QC)

فایل‌های FASTQ حاوی میلیون‌ها توالی کوتاه (Reads) به همراه امتیاز کیفیت متناظر هر باز (Quality Score) هستند.

• امتیاز کیفیت Phred (Q-Score): مقیاس لگاریتمی که احتمال خطا را نشان می‌دهد. Q30 به معنای دقت ۹۹.۹٪ است (یک خطا در هر ۱۰۰۰ باز). در آنالیز استاندارد، میانگین Q-score باید در تمام طول خوانش بالای ۳۰ باشد. افت کیفیت در انتهای خوانش‌ها (3' end) امری رایج در تکنولوژی Illumina است که نیازمند تریم کردن (Trimming) است.
• محتوای GC (GC Content): توزیع غیرنرمال درصد GC می‌تواند نشان‌دهنده آلودگی باکتریایی یا سوگیری (Bias) در مرحله PCR باشد. نمودار توزیع GC باید منطبق بر توزیع نظری ژنوم گونه مورد مطالعه باشد.
• تکرار توالی‌ها (Sequence Duplication): سطوح بالای دوپلیکیشن می‌تواند ناشی از غنی‌سازی بیش از حد PCR (PCR Over-amplification) در کتابخانه‌های با تنوع کم باشد که منجر به کژتابی در کمی‌سازی بیان ژن می‌شود.

۲. استراتژی‌های هم‌ترازسازی (Alignment vs. Pseudo-alignment)

برای شمارش تعداد رونوشت‌های هر ژن، خوانش‌ها باید به ژنوم یا ترنسکریپتوم مرجع نگاشت شوند. دو پارادایم اصلی وجود دارد:

• هم‌ترازسازی مبتنی بر ژنوم (Genome-based Alignment): ابزارهایی مانند STAR و HISAT2 که Splice-aware هستند؛ یعنی می‌توانند اینترون‌ها را شناسایی کرده و خوانش‌هایی که روی مرز اگزون-اگزون قرار می‌گیرند (Gapped reads) را مپ کنند. این روش برای کشف ایزوفرم‌های جدید و واریانت‌ها ضروری است اما به منابع سخت‌افزاری سنگین (بیش از ۳۰ گیگابایت RAM) نیاز دارد. فایل خروجی این مرحله فرمت BAM است.
• شبه‌هم‌ترازسازی (Pseudo-alignment): ابزارهای نوین مانند Salmon و Kallisto. این الگوریتم‌ها به جای مپ کردن دقیق باز-به-باز، از گراف‌های دی‌بروین (De Bruijn Graphs) و k-mers برای تخمین سریع و دقیق تعلق هر خوانش به ترنسکریپت‌ها استفاده می‌کنند. این روش‌ها بسیار سریع‌تر هستند و خطای ناشی از نگاشت‌های چندگانه (Multi-mapping reads) را با مدل‌های احتمالات بیزی (Bayesian inference) اصلاح می‌کنند.

فصل دوم: ورود به R و ساختار داده‌های Bioconductor

پروژه Bioconductor اکوسیستم اصلی آنالیز ژنومیک در R است که ساختارهای داده‌ای اختصاصی را برای مدیریت یکپارچه داده‌های حجیم ارائه می‌دهد.

۱. کلاس SummarizedExperiment

این کلاس شیء (Object Class)، کانتینر مرکزی برای ذخیره‌سازی داده‌هاست و از سه مولفه هماهنگ تشکیل شده است:

• Assay: ماتریس اصلی شمارش‌ها (Counts Matrix) که در آن سطرها نمایانگر ژن‌ها و ستون‌ها نمایانگر نمونه‌ها هستند. دسترسی به این داده‌ها با تابع assay() امکان‌پذیر است.
• ColData: دیتافریمی که اطلاعات فنوتیپی نمونه‌ها (Metadata) مانند گروه تیمار، زمان، جنسیت و بچ آزمایشگاهی را نگه می‌دارد. این بخش معادل pData در ExpressionSet های قدیمی است.
• RowData (GenomicRanges): اطلاعات ساختاری ژن‌ها شامل کروموزوم، موقعیت شروع و پایان و رشته (Strand). این ساختار امکان انجام عملیات جبری روی بازه‌های ژنومی (مانند یافتن همپوشانی با پروموتورها) را فراهم می‌کند.

۲. پکیج tximport: پل ارتباطی بین Salmon و DESeq2

هنگام استفاده از ابزارهایی مانند Salmon، خروجی در سطح ترنسکریپت است، اما آنالیز بیان افتراقی معمولاً در سطح ژن انجام می‌شود تا توان آماری افزایش یابد. پکیج tximport در R وظیفه حیاتی تجمیع (Summarization) شمارش‌های ترنسکریپت به سطح ژن را با استفاده از فایل‌های Annotation (مانند GTF/GFF3) بر عهده دارد. همچنین این پکیج "طول موثر" (Effective Length) ژن‌ها را محاسبه می‌کند که برای نرمال‌سازی دقیق و تصحیح سوگیری طول ژن (Gene Length Bias) ضروری است.

فصل سوم: نرمال‌سازی داده‌ها و مدل‌سازی آماری

داده‌های خام RNA-seq به دلیل تفاوت در عمق توالی‌یابی (Sequencing Depth) و ترکیب RNA (RNA Composition) قابل مقایسه مستقیم نیستند.

۱. چرا RPKM و FPKM منسوخ شده‌اند؟

معیارهای قدیمی FPKM و RPKM تنها برای مقایسه ژن‌های مختلف در "یک نمونه" طراحی شده بودند. در مقایسه بین نمونه‌ها، اگر یک ژن بسیار پربیان در یک نمونه وجود داشته باشد، سهم نسبی سایر ژن‌ها کاهش می‌یابد، حتی اگر بیان واقعی آن‌ها تغییر نکرده باشد. بنابراین استفاده از آن‌ها در آنالیز بیان افتراقی (DGE) ممنوع است. واحد TPM (Transcripts Per Million) جایگزین بهتری برای مقایسه فراوانی نسبی است، اما برای DGE همچنان نیاز به روش‌های نرمال‌سازی مبتنی بر فاکتور مقیاس (Scaling Factor) داریم.

۲. روش نرمال‌سازی DESeq2: میانه نسبت‌ها (Median of Ratios)

این الگوریتم استاندارد فعلی است که فرض می‌کند اکثر ژن‌ها تغییر بیان ندارند (Housekeeping genes predominant):

• ابتدا یک "نمونه مرجع مجازی" (Pseudo-reference sample) با محاسبه میانگین هندسی هر ژن در تمام نمونه‌ها ساخته می‌شود.
• برای هر نمونه، نسبت بیان هر ژن به این مرجع محاسبه می‌شود.
• میانه (Median) این نسبت‌ها به عنوان "فاکتور سایز" (Size Factor) آن نمونه در نظر گرفته می‌شود.
• این روش در برابر ژن‌های پرت (Outliers) و تفاوت‌های شدید در ترکیب کتابخانه بسیار مقاوم است.

۳. توزیع دوجمله‌ای منفی (Negative Binomial Distribution)

شمارش‌های RNA-seq ماهیت گسسته دارند و واریانس آن‌ها معمولاً بزرگتر از میانگین است (پدیده‌ای به نام Overdispersion). توزیع پواسون (Poisson) که در آن واریانس برابر میانگین است، برای این داده‌ها مناسب نیست.

• تخمین پراکندگی (Dispersion Estimation): از آنجا که تعداد تکرارهای بیولوژیکی معمولاً کم است (n=3)، تخمین واریانس برای هر ژن به تنهایی دقیق نیست. DESeq2 با استفاده از روش "اشتراک اطلاعات" (Sharing information across genes)، برآورد پراکندگی ژن‌ها را به سمت یک منحنی کلی سوق می‌دهد (Shrinkage). ژن‌هایی که واریانس بسیار بالاتری نسبت به منحنی مورد انتظار دارند، به عنوان Outlier شناسایی شده و ممکن است از آنالیز حذف شوند (قانون کوک - Cook's distance).

فصل چهارم: آنالیز بیان افتراقی (Differential Expression Analysis)

هدف نهایی، شناسایی ژن‌هایی است که بیان آن‌ها تحت تأثیر متغیر آزمایشی تغییر معنادار داشته است.

۱. آزمون والد (Wald Test)

برای مقایسه دو گروه (مثلاً تیمار در برابر کنترل)، آزمون والد استفاده می‌شود. این آزمون نسبت ضریب تغییرات (Log2 Fold Change) به خطای استاندارد آن را محاسبه می‌کند. هرچه این نسبت بزرگتر باشد، احتمال معناداری بیشتر است.

۲. فرضیه صفر و تصحیح تست‌های چندگانه

در RNA-seq، ما برای ۲۰,۰۰۰ ژن به طور همزمان آزمون آماری انجام می‌دهیم. اگر سطح معناداری (alpha) را ۰.۰۵ در نظر بگیریم، به طور تصادفی ۵٪ از ژن‌ها (یعنی ۱۰۰۰ ژن) ممکن است به اشتباه معنادار شوند (False Positives).

• نرخ کشف کاذب (FDR - False Discovery Rate): برای حل این مشکل، از روش Benjamini-Hochberg استفاده می‌شود تا P-value های خام به Adjusted P-values تبدیل شوند. معمولاً padj < 0.05 به عنوان آستانه پذیرش ژن‌های افتراقی (DEGs) در نظر گرفته می‌شود.
• Independent Filtering: الگوریتمی که به طور خودکار ژن‌هایی با میانگین بیان بسیار پایین را از تست‌های آماری حذف می‌کند تا قدرت آزمون (Statistical Power) برای سایر ژن‌ها افزایش یابد و جریمه اصلاح FDR کاهش یابد.

۳. جمع‌بندی تغییرات (LFC Shrinkage)

ژن‌هایی که بیان بسیار کمی دارند، اغلب نویز بالایی نشان می‌دهند و ممکن است به طور تصادفی Fold Change های بزرگی داشته باشند. DESeq2 از روش‌هایی مانند apeglm یا ashr برای تعدیل (Shrinkage) تغییرات لاگاریتمی این ژن‌ها استفاده می‌کند. این کار P-value را تغییر نمی‌دهد اما رتبه‌بندی ژن‌ها بر اساس Fold Change را برای آنالیزهای بعدی (مانند GSEA) دقیق‌تر می‌کند.

فصل پنجم: تحلیل اکتشافی داده‌ها (EDA) و مصورسازی

پیش از نتیجه‌گیری نهایی، کیفیت تفکیک نمونه‌ها باید بصری‌سازی شود.

۱. تبدیل‌های تثبیت واریانس (VST rlog)

برای مصورسازی (مانند PCA و Heatmap)، داده‌های شمارشی باید تبدیل شوند تا وابستگی واریانس به میانگین حذف شود (Homoscedasticity). روش VST (Variance Stabilizing Transformation) برای دیتاست‌های بزرگ (n > 30) سریع‌تر است، در حالی که rlog (Regularized Logarithm) برای دیتاست‌های کوچک حساس‌تر عمل می‌کند.

۲. تحلیل مؤلفه‌های اصلی (PCA)

PCA یک تکنیک کاهش ابعاد است که واریانس کلی داده‌ها را در دو محور اصلی خلاصه می‌کند. در یک آزمایش ایده‌آل:

• PC1 (محور X): باید بیشترین واریانس را داشته باشد و نمونه‌ها را بر اساس متغیر اصلی آزمایش (مثلاً تیمار و کنترل) جدا کند.
• PC2 (محور Y): باید واریانس‌های ثانویه را نشان دهد. اگر جداسازی نمونه‌ها در PCA بر اساس متغیرهای مزاحم (مانند روز آزمایش یا اپراتور) باشد، نشان‌دهنده وجود "اثر دسته‌ای" (Batch Effect) است که باید اصلاح شود.

۳. نقشه‌های حرارتی (Heatmaps)

استفاده از پکیج‌های pheatmap یا ComplexHeatmap برای نمایش الگوی بیان ژن‌های منتخب (مثلاً ۵۰ ژن با بیشترین واریانس). استفاده از خوشه‌بندی سلسله‌مراتب (Hierarchical Clustering) با متریک‌های فاصله اقلیدسی یا همبستگی پیرسون، به شناسایی گروه‌های هم‌بیان (Co-expressed genes) کمک می‌کند. داده‌های ورودی به هیت‌مپ باید حتماً نرمال‌سازی و تبدیل شده باشند (Z-score scaling).

فصل ششم: آنالیز غنی‌سازی عملکردی (Functional Enrichment Analysis)

لیستی از نام ژن‌ها به تنهایی فاقد ارزش بیولوژیکی است. برای درک مکانیسم‌ها، باید ژن‌ها را به عملکردها ترجمه کنیم.

۱. Gene Ontology (GO) Analysis

پایگاه داده GO ژن‌ها را در سه دسته طبقه‌بندی می‌کند: فرآیندهای بیولوژیکی (BP)، اجزای سلولی (CC) و عملکرد مولکولی (MF). آنالیز ORA (Over-Representation Analysis) با استفاده از آزمون فیشر یا فوق‌هندسی بررسی می‌کند که آیا ژن‌های یک مسیر خاص در لیست DEGs ما بیش از حد انتظار تصادفی حضور دارند یا خیر. ساختار GO به صورت گراف جهت‌دار بدون دور (DAG) است، بنابراین نتایج ممکن است همپوشانی داشته باشند که با متد simplify در پکیج clusterProfiler قابل کاهش است.

۲. تحلیل مسیرهای KEGG

پایگاه داده KEGG ژن‌ها را در شبکه‌های تعاملی (Pathway Maps) قرار می‌دهد. این آنالیز به ما می‌گوید که تغییر بیان ژن‌ها کدام مسیرهای سیگنالینگ (مثلاً مسیر NF-kappaB یا Apoptosis) را فعال یا مهار کرده است. استفاده از پکیج pathview اجازه می‌دهد تا تغییرات بیان (Fold Change) را مستقیماً روی نقشه‌های رنگی KEGG مشاهده کنیم.

۳. آنالیز GSEA (Gene Set Enrichment Analysis)

برخلاف ORA که نیاز به یک آستانه برش (Cut-off) دلخواه دارد، GSEA از تمام ژن‌های آنالیز شده استفاده می‌کند. ژن‌ها بر اساس میزان تغییر بیان رتبه‌بندی می‌شوند و الگوریتم بررسی می‌کند که آیا اعضای یک مجموعه ژنی خاص در بالا یا پایین این لیست رتبه‌بندی شده تجمع یافته‌اند. این روش بسیار قدرتمند است و می‌تواند تغییرات جزئی اما هماهنگ در یک مسیر بیولوژیکی را شناسایی کند که ممکن است در روش ORA نادیده گرفته شوند.

فصل هفتم: مدیریت اثرات دسته‌ای (Batch Effects)

یکی از بزرگترین چالش‌ها در پروژه‌های بزرگ RNA-seq، اثرات تکنیکال ناخواسته است که ناشی از انجام آزمایش در زمان‌های مختلف، توسط افراد مختلف یا با سری‌های متفاوت معرف‌هاست.

• تشخیص: اگر در نمودار PCA نمونه‌ها بر اساس تاریخ آزمایش یا سری ساخت کیت جدا شوند، اثر دسته‌ای وجود دارد.
• اصلاح در مدل‌سازی (توصیه شده): افزودن متغیر Batch به فرمول طراحی DESeq2 (مثلاً: design = ~ batch + condition). این کار به مدل اجازه می‌دهد تا اثر Batch را تخمین زده و آن را از اثر تیمار تفکیک کند.
• حذف مستقیم (ComBat): استفاده از تابع ComBat-seq از پکیج sva. این روش داده‌های شمارشی را مستقیماً تغییر می‌دهد و تنها باید برای مصورسازی (PCA/Heatmap) استفاده شود. استفاده از داده‌های تغییر یافته توسط ComBat به عنوان ورودی برای آنالیز DGE توصیه نمی‌شود زیرا فرضیات توزیع آماری را مخدوش می‌کند.