Integrated analysis of single-cell and bulk RNA-sequencing identifies a signature based on drug response genes to predict prognosis and therapeutic response in ovarian cancer

این مطالعه در حوزه بیوانفورماتیک سرطان و انکولوژی دقیق (Precision Oncology) به مسئله ناهمگونی توموری و مقاومت دارویی در سرطان تخمدان می پردازد. نویسندگان با ادغام single-cell RNA sequencing و توالی یابی RNA بالک (bulk RNA sequencing)، زیرجمعیت های سلولی مرتبط با حساسیت به پاکلیتاکسل (Paclitaxel) را شناسایی کرده و بر پایه آن ها یک امضای ژنی پنج گانه برای ارزیابی پیش آگهی بیماران طراحی کرده اند. مسئله محوری مقاله این است که پاسخ بیماران مبتلا به سرطان تخمدان به شیمی درمانی یکسان نیست و در نبود نشانگرهای زیستی (biomarkers) دقیق، انتخاب درمان مناسب و پیش بینی پیامد بالینی با دشواری همراه می شود. اهمیت این موضوع از آن رو است که سرطان تخمدان از بدخیم ترین تومورهای زنان به شمار می رود و در بسیاری از بیماران در مراحل پیشرفته تشخیص داده می شود. درمان استاندارد معمولا شامل جراحی کاهنده حجم تومور و سپس شیمی درمانی بر پایه پلاتین و تاکسان است، اما بروز مقاومت دارویی و پیچیدگی ریزمحیط تومور (tumor microenvironment; TME) موجب می شود بخش قابل توجهی از بیماران از درمان سود کامل نبرند. در همین چارچوب، مقاله نشان می دهد که ترکیب داده های چندلایه مولکولی می تواند به طراحی مدل های دقیق تر برای طبقه بندی بیماران کمک کند. رویکرد کلی مقاله بر استفاده از الگوریتم Scissor استوار است؛ روشی محاسباتی که امکان اتصال فنوتیپ های استخراج شده از داده های bulk RNAseq به سلول های منفرد را فراهم می کند. نویسندگان ابتدا سلول های اپی تلیالی توموری را از داده های تک سلولی استخراج کرده، سپس با استفاده از داده های مقاومت دارویی در سطح bulk، سلول های مرتبط با حساسیت به پاکلیتاکسل را شناسایی کرده اند. ژن های بیان شده در این زیرجمعیت به عنوان ژن های مرتبط با حساسیت به پاکلیتاکسل انتخاب شدند و در ادامه باLASSO Cox یک مدل پیش آگهی پنج ژنی ساخته شد. افزون بر این، تعاملات سلول-سلول، ویژگی های ایمنی و الگوهای جهشی نیز میان گروه های مختلف خطر بررسی شد. این مقاله برای پژوهشگران بیوانفورماتیک سرطان، زیست شناسی سامانه ها، انکولوژی مولکولی و نیز دانشجویان تحصیلات تکمیلی که به کاربردهای عملی ادغام داده های scRNA-seq و bulk RNA-seq در کشف نشانگرهای زیستی علاقه مند هستند، ارزش ویژه ای دارد.

سرطان تخمدان یکی از مهاجم ترین بدخیمی های دستگاه تناسلی زنان است و به علت علائم اولیه مبهم، اغلب در مراحل پیشرفته شناسایی می شود. همین ویژگی سبب شده است که این بیماری در بسیاری از منابع به عنوان یک بیماری خاموش معرفی شود. با وجود پیشرفت در جراحی و شیمی درمانی، پیش آگهی بیماران همچنان نامطلوب است و عود بیماری پس از درمان اولیه پدیده ای شایع به شمار می رود. درمان استاندارد پس از جراحی عمدتا بر پایه ترکیباتی مانند پاکلیتاکسل و کربوپلاتین (Carboplatin) است، اما تفاوت های زیستی میان بیماران و ناهمگونی درون تومور باعث می شود پاسخ به درمان به شدت متغیر باشد. مطالعات پیشین با استفاده از داده های توالی یابی bulk RNAseq توانسته اند بخشی از مبانی مولکولی مقاومت دارویی را آشکار کنند. با این حال، داده های bulk RNAseq میانگین بیان ژن را در جمعیتی مرکب از سلول های مختلف اندازه می گیرند و به همین دلیل قادر نیستند تفاوت های ظریف اما تعیین کننده میان زیرجمعیت های سلولی را به خوبی نشان دهند. در مقابل، single-cell RNAseq امکان مطالعه ناهمگونی رونویسی در سطح هر سلول را فراهم می کند و برای شناسایی زیرجمعیت های نادر اما زیست شناختی مهم بسیار سودمند است. با وجود این مزیت، استفاده از scRNA-seq به تنهایی نیز محدودیت هایی دارد؛ از جمله هزینه بالا، حجم نمونه کمتر و دشواری پیوند دادن مستقیم یافته های تک سلولی به پیامدهای بالینی در مقیاس بزرگ. از این رو، یکی از نیازهای مهم این حوزه توسعه راهبردهایی برای ادغام داده های تک سلولی و توده ای است تا بتوان هم از قدرت تفکیک سلولی scRNA-seq و هم از غنای اطلاعات بالینی داده های bulk RNA-seq بهره برد. الگوریتم Scissor دقیقا برای پاسخ به این نیاز طراحی شده است. در این مطالعه، هدف نویسندگان آن است که با ادغام داده های مقاومت به پاکلیتاکسل در سطح bulk و داده های single-cell بیماران مبتلا به سرطان تخمدان، زیرجمعیت های سلولی مرتبط با حساسیت دارویی را شناسایی کنند، از میان ژن های آن ها یک امضای ژنی پیش آگهی استخراج کنند و سپس تفاوت های زیستی، ایمنی و جهشی میان گروه های حاصل را بررسی کنند. به این ترتیب، مقاله می کوشد پلی میان زیست شناسی تک سلولی و کاربرد بالینی در پیش بینی پیامد بیماری ایجاد کند.

این پژوهش یک مطالعه محاسباتی و تحلیلی مبتنی بر داده های عمومی است که در چند مرحله اصلی انجام شده است. در مرحله نخست، داده های single-cell RNAseq از چند مجموعه داده عمومی شامل نمونه های سرطان تخمدان جمع آوری شد. همچنین داده های توالی یابی bulk RNAseq از بیماران سرطان تخمدان و نیز داده های مرتبط با مقاومت به پاکلیتاکسل در رده های سلولی تهیه شد تا مبنای ادغام تحلیلی فراهم شود (از GEO و TCGA کار دریافت دیتا صورت پذیرفت). در پردازش داده های تک سلولی، سلول ها بر اساس معیارهای کنترل کیفیت پالایش شدند. سلول های با تعداد ژن بسیار کم یا بسیار زیاد، درصد بالای ژن های میتوکندریایی، و شاخص های غیرطبیعی دیگر حذف شدند تا فقط سلول های با کیفیت مناسب برای تحلیل باقی بمانند. سپس نرمال سازی، شناسایی ژن های با تغییرپذیری بالا، و حذف اثر دسته ای (batch effect) انجام شد. برای کاهش ابعاد از تحلیل مولفه های اصلی (principal component analysis; PCA) و برای نمایش داده ها و خوشه بندی از uniform manifold approximation and projection یا همان  UMAP استفاده شد. سلول های دوتایی احتمالی نیز حذف شدند تا دقت شناسایی جمعیت های سلولی افزایش یابد. در ادامه، خوشه های سلولی با استفاده از نشانگرهای متعارف (canonical markers) تعیین هویت شدند و به جمعیت هایی مانند سلول های اپی تلیالی، فیبروبلاست ها، سلول های اندوتلیالی، سلول های T، سلول های B و سلول های میلوئیدی تفکیک شدند. پس از آن، نویسندگان بر سلول های اپی تلیالی توموری تمرکز کردند تا ناهمگونی درون تومور را با دقت بیشتری بررسی کنند. برای شناسایی سلول های مرتبط با پاسخ دارویی، الگوریتم Scissor به کار گرفته شد. این الگوریتم با اتصال الگوهای بیان ژن در سطح تک سلولی به فنوتیپ های مشاهده شده در داده های bulk RNAseq، سلول هایی را که بیشترین ارتباط را با حساسیت یا مقاومت به پاکلیتاکسل دارند مشخص می کند. سلول هایی که با الگوی حساسیت به پاکلیتاکسل هم خوان بودند به عنوان سلول های حساس به پاکلیتاکسل شناخته شدند و ژن های با بیان بالاتر در آن ها به عنوان ژن های حساسیت به پاکلیتاکسل انتخاب شدند. پس از شناسایی این ژن ها، تحلیل غنی سازی عملکردی (functional enrichment analysis) با استفاده از پایگاه های GO و KEGG انجام شد تا فرایندها و مسیرهای زیستی مرتبط مشخص شود. همچنین برای بررسی شبکه های ارتباطی در ریزمحیط تومور، تحلیل برهم کنش لیگاند-گیرنده با CellChat صورت گرفت. در کنار آن، الگوهای جهش سوماتیک، بار جهشی تومور (tumor mutational burden; TMB)، و ترکیب سلول های ایمنی نفوذ کرده به تومور نیز میان گروه های مختلف مقایسه شد. در نهایت، برای ساخت مدل پیش آگهی از رگرسیون LASSO Cox استفاده شد. این روش به انتخاب ژن های کلیدی و ساخت یک امتیاز خطر (risk score) بر پایه ضرایب مدل و میزان بیان ژن ها منجر شد. بیماران بر اساس این امتیاز به گروه های مختلف تقسیم شدند و تفاوت بقای کلی (overall survival; OS) میان آن ها با تحلیل کاپلان-مایر (Kaplan-Meier) بررسی شد. علاوه بر این، نوموگرام (nomogram) با ترکیب امتیاز خطر و متغیرهای بالینی ساخته شد تا ابزاری کاربردی تر برای برآورد بقای بیماران ارائه شود.

اطلس سلولی سرطان تخمدان پس از اعمال مراحل کنترل کیفیت و یکپارچه سازی داده ها، مجموعه ای گسترده از سلول های با کیفیت بالا از بیماران مبتلا به سرطان تخمدان به دست آمد. این سلول ها در چندین خوشه مجزا سازمان دهی شدند که در نهایت به شش رده اصلی شامل سلول های اپی تلیالی، فیبروبلاست ها، سلول های اندوتلیالی، سلول های T، سلول های B و سلول های میلوئیدی نسبت داده شدند. این اطلس سلولی نشان داد که بافت سرطان تخمدان از نظر ترکیب سلولی و برنامه های رونویسی به شدت ناهمگون است. تمرکز بر سلول های اپی تلیالی توموری نشان داد که این جمعیت نیز یکنواخت نیست و خود از چند زیرخوشه تشکیل شده است. این زیرخوشه ها الگوهای متفاوتی از بیان ژن، وضعیت چرخه سلولی و امضای گذار اپی تلیال به مزانشیمی (epithelial-mesenchymal transition; EMT) را نشان دادند. چنین یافته ای بر این نکته تاکید دارد که حتی در درون سلول های توموری نیز جمعیت هایی با ویژگی های زیستی متفاوت حضور دارند که می توانند در پاسخ به درمان نقش تعیین کننده داشته باشند. شناسایی سلول های حساس به پاکلیتاکسل و ویژگی های مولکولی آن ها با به کارگیری الگوریتم Scissor، زیرمجموعه ای از سلول های توموری شناسایی شد که بیشترین ارتباط را با فنوتیپ حساسیت به پاکلیتاکسل داشتند. این سلول ها به عنوان سلول های حساس به پاکلیتاکسل معرفی شدند و ژن هایی که در آن ها بیان بالاتری داشت، به عنوان ژن های مرتبط با حساسیت به پاکلیتاکسل در نظر گرفته شد. بررسی پیامد بالینی نشان داد که افزایش بیان این ژن ها با بقای کلی مطلوب تر همراه است و این موضوع ارزش بالقوه آن ها را به عنوان نشانگر پیش آگهی تقویت می کند. تحلیل های عملکردی نشان داد که این ژن ها عمدتا در فرایندهای مرتبط با پاسخ به استرس سلولی، تاشدگی پروتئین، تنظیم رونویسی و برخی مسیرهای متابولیک نقش دارند. در تحلیل KEGG نیز غنی سازی در مسیرهایی مانند پردازش پروتئین در شبکه آندوپلاسمی (endoplasmic reticulum)، متابولیسم لیپید و برخی مسیرهای پیام رسانی مشاهده شد. این نتایج نشان می دهد که حساسیت به پاکلیتاکسل صرفا به یک ژن منفرد وابسته نیست، بلکه بازتاب یک حالت زیستی چندوجهی در سلول های توموری است. ریزمحیط تومور و ارتباطات سلول-سلول نمونه ها بر اساس امتیاز حاصل از ژن های حساسیت به پاکلیتاکسل به گروه های با امتیاز بالا و پایین تقسیم شدند و سپس تفاوت های ریزمحیط تومور میان آن ها بررسی شد. یکی از یافته های مهم این بود که در گروه دارای امتیاز بالاتر، حضور سلول های T پررنگ تر بود و شبکه ارتباطی میان سلول های توموری و سلول های ایمنی الگوی متفاوتی نشان می داد. تحلیل CellChat نشان داد که اگرچه شمار برخی برهم کنش ها در گروه های مختلف متفاوت بود، شدت ارتباطات خاص در گروه دارای امتیاز بالاتر افزایش می یافت. این تغییر به ویژه در مسیرهای مرتبط با MIF و گیرنده های آن، از جمله CD74+CXCR4 و CD74+CD44، میان سلول های اپی تلیالی و سلول های T برجسته بود. این یافته ها از نقش مهم ارتباطات میان سلول های توموری و اجزای ایمنی در تعیین پاسخ به شیمی درمانی پشتیبانی می کند و نشان می دهد که حساسیت دارویی صرفا ویژگی درون سلولی نیست، بلکه در تعامل با ریزمحیط تومور شکل می گیرد. چشم انداز جهشی و ویژگی های ایمنی بررسی الگوهای جهش سوماتیک نشان داد که برخی ژن ها در هر دو گروه با فراوانی بالا دچار جهش هستند، اما تفاوت هایی نیز در مسیرهای جهشی غالب دیده می شود. در کنار آن، مقایسه بار جهشی تومور تفاوت چشمگیری میان گروه ها نشان نداد، که حاکی از آن است که تفاوت پیامد بالینی بیشتر ممکن است به برنامه های رونویسی و تعاملات ریزمحیطی وابسته باشد تا صرفا تعداد کلی جهش ها. تحلیل نفوذ ایمنی نیز بیانگر آن بود که گروه های مختلف خطر از نظر ترکیب سلول های ایمنی یکسان نیستند. برخی زیرجمعیت های لنفوسیتی و سلول های ارائه کننده آنتی ژن در گروه ها به نسبت های متفاوت حضور داشتند. این تفاوت ها نشان می دهد که وضعیت ایمنی تومور می تواند در کنار امضای ژنی، به توضیح تفاوت در بقای بیماران و پاسخ به درمان کمک کند. امضای ژنی پنج گانه و مدل پیش آگهی از میان ژن های مرتبط با حساسیت به پاکلیتاکسل، مجموعه ای از ژن های دارای ارتباط معنی دار با بقا شناسایی شد و در نهایت با استفاده از LASSO Cox پنج ژن HSP90AA1، JUN، CALM1، RBP1 و GLRX5 برای ساخت مدل نهایی انتخاب شدند. بر اساس بیان این ژن ها، برای هر بیمار امتیاز خطر محاسبه شد و بیماران به دو گروه تقسیم شدند. این مدل توانست تفاوت بقای کلی را میان گروه ها نشان دهد و در مجموعه داده های مستقل نیز تا حد قابل قبولی اعتبارسنجی شد. همچنین تحلیل های چندمتغیره نشان داد که امتیاز خطر در کنار برخی ویژگی های بالینی، یک عامل مستقل برای پیش بینی بقا است. در ادامه، نویسندگان با تلفیق این امتیاز و اطلاعات بالینی، نوموگرامی برای پیش بینی بقای یک، سه و پنج ساله طراحی کردند که از نظر کالیبراسیون نیز عملکرد مطلوبی داشت.

یافته های این مطالعه نشان می دهد که ادغام داده های scRNA-seq و bulk RNA-seq می تواند به شناسایی زیرجمعیت های سلولی مرتبط با پاسخ دارویی در سرطان تخمدان منجر شود و این زیرجمعیت ها مبنای مناسبی برای ساخت مدل های پیش آگهی فراهم می کنند. ارزش اصلی مقاله در این است که حساسیت به پاکلیتاکسل را نه فقط به صورت یک پدیده بالینی، بلکه به عنوان نتیجه تعامل میان برنامه های درون سلولی، وضعیت رونویسی و ریزمحیط ایمنی بررسی می کند. نویسندگان به ویژه بر اهمیت مسیرهای ارتباطی میان سلول های اپی تلیالی و سلول های T تاکید می کنند. افزایش برخی برهم کنش های لیگاند-گیرنده در گروه دارای امتیاز بالاتر ژن های حساسیت به پاکلیتاکسل، نشان می دهد که پاسخ مطلوب تر به درمان ممکن است با وضعیت ایمنی فعال تر یا دست کم با معماری متفاوت تری از ارتباطات ایمنی-توموری همراه باشد. در نتیجه، این مطالعه از این دیدگاه حمایت می کند که شیمی درمانی و ایمنی ضدتوموری دو حوزه جدا از هم نیستند، بلکه به صورت پویا بر یکدیگر اثر می گذارند. در عین حال، مقاله چند محدودیت مهم نیز دارد. نخست آن که بیشتر تحلیل ها بر داده های عمومی و گذشته نگر تکیه دارد و اطلاعات بالینی کامل، به ویژه درباره رژیم های درمانی و پاسخ واقعی بیماران، در همه مجموعه داده ها در دسترس نیست. دوم آن که پیوند میان سلول های شناسایی شده با الگوریتم Scissor و پاسخ واقعی بیماران به پاکلیتاکسل بیشتر غیرمستقیم است، زیرا بخشی از داده های مرجع از رده های سلولی به دست آمده است. سوم آن که یافته های مربوط به ارتباطات سلول-سلول و تفاوت های ایمنی در سطح محاسباتی باقی مانده و هنوز نیازمند اعتبارسنجی تجربی هستند. نکته قابل توجه دیگر این است که جهت امتیاز خطر در این مطالعه تا حدی غیرمعمول به نظر می رسد، زیرا در برخی بخش ها گروه با امتیاز بالاتر پیامد بقایی بهتری نشان می دهد. برای چنین الگویی، تبیین زیستی دقیق تر می تواند به فهم بهتر مدل کمک کند. با این حال، این مسئله از ارزش کلی مطالعه در معرفی یک چارچوب تحلیلی نوآورانه نمی کاهد.

این مقاله نشان می دهد که ادغام داده های single-cell RNAseq و توده ای، همراه با استفاده از الگوریتم Scissor، راهبردی موثر برای کشف زیرجمعیت های سلولی مرتبط با حساسیت دارویی در سرطان تخمدان است. بر پایه این راهبرد، نویسندگان یک امضای ژنی پنج گانه شامل HSP90AA1، JUN، CALM1، RBP1 و GLRX5 معرفی کردند که توانست در چند مجموعه داده، بیماران را از نظر پیش آگهی تفکیک کند. افزون بر آن، تحلیل های ریزمحیط تومور، شبکه های ارتباطی سلولی و الگوهای جهشی نشان داد که پاسخ به پاکلیتاکسل پدیده ای چندعاملی است و باید در بستری گسترده تر از ویژگی های مولکولی و ایمنی تومور تفسیر شود. در مجموع، این مطالعه چارچوبی علمی و کاربردی برای توسعه نشانگرهای زیستی در سرطان تخمدان ارائه می دهد و می تواند زمینه ساز مطالعات آینده نگر و چندمرکزی برای ارزیابی قابلیت استفاده بالینی این مدل باشد.

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

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