مقدمه

در این مقاله به صورت خلاصه در مورد گام های اساسی برای فراخوانی واریانت های ژنتیکی و ابزارهایی که مقبول عموم هستند، صحبت می کنیم. فرض کنید هدف یافتن واریانت های کوتاه باشد. برای این منظور میزان عمق پیشنهادی برای توالی یابی در مطالعات کل ژنوم[1]   برابر 30-50 و برای مطالعات کل اگزوم[2] برابر 70-100 می باشد. مراحل کار به صورت شماتیک به سه قسمت زیر قابل تقسیم می باشد:

مراحل آنالیز فراخوانی واریانت
 مراحل فراخوانی واریانت ها variant calling))

[1] Whole Genome Sequencing (WGS)

[2] Whole Exome Sequencing (WES)

توجه: ابزار GATK  (Genome Analysis Toolkit) که توسط Broad Institute of Harvard and MIT توسعه یافته است، توسط عموم برای انجام مراحل این آنالیز استفاده می گردد. نسخه جدید این مجموعه ابزار GATK4 به صورت کاملا متن باز بوده قابلیت اجرا روی کلاسترهای سپارک و ابر گوگل را نیز علاوه بر اجرا روی سیستم لوکال را دارد.

فاز اول: پیش پردازش داده های خام

گام های زیر باید برای آماده سازی داده های خام انجام پذیرند.

گام اول- 1: حذف توالی ای آداپتور (Adapter trimming):

خروجی دستگاه های توالی یابی معمولا خوانش هایی با فرمت fastq می باشد که شامل توالی نوکلئوتیدها و کیفیت خوانش هر نوکلئوتید خوانده شده می باشد. معمولا توالی آداپتور قبل از ارائه داده خام به مشتری حذف می گردد، ولی بعضی مواقع مقداری از توالی آداپتور باقی می ماند، این توالی ممکن از در هر قسمت از خوانش و هر مقئاری از طول آداپتور باشد. لذا لازم است که این آلودگی حذف گردد. از ابزارهای رایجی که برای حذف آداپتور می توان استفاده کرد می تواند به ابزار Trimmomatic [1]، Fastx-toolkit (fastx-clipper)، Bioconductor(Short read Package) و ... استفاده نمود.

گام اول-2 : حذف خوانش ها بر اساس کیفیت (Quality trimming):

بعد از حذف آلودگی آداپتور بهتر است خوانش های با کیفیت پایین و یا نواحی از خوانش که کیفیت پایینی دارند حذف گردد. در تحلیل واریانت ها کیفیت پیشنهادی عدد بالای 30 برای Phred score می باشد. برای این منظور می توان از ابزارهایی مانند FASTQC [2] ، PrinSeq [3] برای بررسی کیفیت خوانش ها استفاده کرد. سپس با ابزارهای برش مانند Trimmomatic، FASTX-Toolkit و PrinSeq توالی های نامطلوب را حذف نمود.

گام اول- 3: حذف خوانش های بسیار کوتاه:

در این مرحله لازم است خوانش هایی که طول آن­ها بسیار کوتاه است (مثلا زیر 20   جفت باز( را حذف نمود. توجه شود که این توالی های بسیار کوتاه می توانند به قسمت های مختلفی از ژنوم به صورت تصادفی مپ گردند و باعث خطا شوند. همچنین انتخاب این طول به طول fragment ها بستگی دارد. برای حذف این خوانش ها از ابزارهایی مانند Trimmomatic و PrinSeq می توان استفاده نمود.

 

 

 

فاز دوم: فراخوانی واریانت های خام

گام دوم- 1: هم ترازی (Alignment)

برای یافتن نواحی مشابه و پلیمورفیسم ها در یک نمونه، باید خوانش های مربوط به آن نمونه به ژنوم مرجع مپ شوند. برای این منظور ژنوم Hg38 از پروژه 1000 ژنوم به عنوان مرجع پیشنهاد می گردد.  برای مپ کردن خوانش ها از الگوریتم های مختلفی می توان استفاده نمود. ولی از میان آن ها الگوریتم های BWA MEM و Bowtie2 الگوریتم های مورد اعتماد برای خوانش ها دستگاه های Illumina می باشند. همراهی این الگوریتم ها با الگوریتم های مختلف فراخوانی واریانت ها، دقت های مختلفی برای شناسایی SNP ها InDel ها و واریانت های ساختاری به همراه خواهد داشت. به عنوان یک پیشنهاد الگوریتم BWA –MEM با تنظیمات (-K 100000000  -Y  و بدون پارامتر –M) پیشنهاد می گردد.خروجی این الگوریتم یک فایل BAM می باشد که حجمی از مرتبه چند 10 گیگابایت می تواند اشغال نماید. برای تسریع این مرحله از کار پیشنهاد می شود، ژنوم رفرنس روی یک هارد SSD قرار داشته باشد. همچنین این مرحله قابلیت پردازش موازی و حتی بر روی GPU را دارد.

گام دوم -2: حذف خوانش های تکراری (De-duplication):

وجود خوانش های تکراری در توالی یابی یک مشکل است که در روش های توالی یابی بر مبنای سنتر که از PCR برای آماده سازی نمونه استفاده می کنند، اجناب ناپذیر است. توالی های تکراری نباید به عنوان شاعدی از وجود یا عدم وجود یک واریانت محسوب شوند و لذا باید قبل از مرحله فراخوانی واریانت ها حذف گردند. برای این منظور از ابزارهایی مانند samblaster [4]، ابزار تجاری Novosort از مجموعه ی Novocraft  [5]، Picard و FASTX-Toolkit(fastx-collapser) می توان استفاده نمود. در حال حاضر ابزار MarkDuplicates از مجموعه ابزار GATK4 فراخوانی می گردد، ولی در نسخه های قبلی از  مجموعه ابزار Picard [6] خوانده می شد.

گام دوم -3 : Artifact removal(Local realignment around indels)

بعضی خطا ها به در مرحله ی هم ترازی مخصوصا حوالی نواحی Indel که خوانش ها شروع یا پایان یک Indel را پوشش می دهند اتفاق می افتد. این mismatch های ایجاد شده می توانند با SNP ها اشتباه گرفته شوند. لذا در این مرحله بازهمترازی خوانش ها در این حوالی انجام می پذیرد و خوانش های این نواحی با خوانش اجماع شده از آن ها جایگزین می گردد تا فراخوانی واریانت ها به درستی انجام پذیرد. این رایند توسط ابزاری مانند GATK IndelRealigner [7] یا ابزارهای دیگری مانند Dindel [8] و SRMA [9] میتواند انجام پذیرد.

وجود مرحله ی بازهمترازی در فراخوانی واریانت ها به ابزار و الگوریتم مربوط به فراخوانی واریانت ها در مرحله ی بعدی وابسته است. در واقع این مرحله در صورتی با ارزش می باشد که الگوریتم فراخوانی واریانت ها در اصطلاح non-haplotype-aware باشد. به عنوان مثال ابزار UnifiedGenotyper از مجموعه ی GATK این حالت را دارد. ولی در صورتی که ابزار فراخوانی واریانت ها آگاه نسبت به هاپلوتایپ باشد، مثلا Platypus [10] ، FreeBayes [11]  یا HaplotypeCaller (GATK) )، در این حالت این گام لازم نبوده و توصیه نمی شود.

گام دوم -4 : Base quality score recalibration:

کیفیتی که ماشین توالی یابی به هر جفت باز هنگام خواندن توالی می دهد، معموملا می تواند غیر دقیق و یا دارای خطای بایاس باشد. لذا در این مرحله بر اساس یک مدل تجربی که از روی خصوصیات داده ساخته می شود، این خطا تصحیح می گردد. اینکار توسط ابزاری مانند GATK BQSR قابل انجام است و توصیه می شود برای مقایسه ی عملکرد الگوریتم های مختلف این گام صورت پذیرد. برای افزایش سرعت آنالیز در صورتی که از GATK4 استفاده می نمایید می توانید خروجی BaseRecalibrator را مستقیما به HaplotypeCaller بدهید  و از انجام BQSR بین این دو مرحله عبور نمایید. ابزار جایگزین برای این منظور ReQON از پکیج های Bioconductor می باشد.

گام دوم- 5: فراخوانی واریانت ها (Variant Calling):

برای فراخوانی واریانت های محتمل، یک ابزار بهینه وجود ندارد و توصیه می شود از اجماع  سه ابزار استفاده نمایید:

  1. CRISP [12]
  2. HaplotypeCaller (GATK)
  3. Mpileup(SAMtools)

همچنین اخیرا ابزار MuTect2 به عنوان ابزار شناسایی واریانت های اکتسابی (سوماتیک) برای کاربرد در مجموعه ابزارهای GATK ارائه شده است. این ابزار SNP ها و Indel های سوماتیک را با ترکیب دو ابزار Mutect [13] و HaplotypeCaller انجام می دهد. HaplotypeCaller با فرض دیپلویدی کار می کند، در حالی که Mutect2 اجازه ی نسبت های آللی متفاوت را می دهد که برای یافتن جهش های سوماتیک کاربرد دارد.  خروجی این مرحله فایل VCF نام دارد که حجم به مراتب کمتری نسبت به فایل BAM خواهد داشت ( ده ها مگابایت).

گام دوم 6: فیلتر کردن آماری واریانت ها:

نتیجه ی بدست آمده در قسمت قبل شامل تعداد بسیار زیادی واریانت غیر واقعی می باشد که در واقع نتیجه ی خطاهای توالی یابی و مراحل بعد از آن می باشد. در مطالعات کوچک معمولا استفاده از فیلتر سخت (Hard Filtering) برای جداسازی واریانت ها کافی می باشد، البته مشخص کردن آستانه های این فیلتر نیاز به تخصص مناسب دارد. در این مقاله یک سری قواعد کلی برای انتخاب این نوع فیلتر ارائه شده است [14].  برای مطالعات با تعداد نمونه های بالا (بیشتر از 30 نمونه) ، GATK ابزار VQSR را برای جداسازی واریانت های مثبت کاذب از مثبت واقعی توسط یک مدل Gaussian mixture  با ساتفاده از الگوریتم های یادگیری ماشین ارائه داده است. برای این منظور ماشین ابتدا پارامترهای این مدل را با استفاده از داده های واریانت های صحیح یادگرفته و بعد مورد استفاده قرار می دهد.

فاز سوم: تفسیر و الویت بندی واریانت ها

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

  • جستجو و الویت بندی واریانت ها در حوزه ی بیماری های انسانی با ابزار GEMINI [15].
  • واریانت های مرتبط با بیماری های مندلی : VAR-MD، KGGSeq، OMIM، FamSeq
  • پیشبینی میزان حذف کنندگی SNP  ها با ابزارهای بیوانفورماتیکی مانند: ANNOVAR, VAAST
  • شناسایی واریانت های نواحی تنظیمی: RegulomeDB