في الجزء الأول من هذا المقال عَرَفْنا بحزمة Biostrings وأعطينا نظرة عامة عن بعض امكانياتها. في هذا الجزء سوف نواصل من حيث توقفنا في الجزء الأول وسوف نتحدث عن الحزم والدوال المستعملة للتعامل مع الإحداثيات الجينومية وقرائة الملفات مثل BAM أو FASTQ.
التعامل مع الإحداثيات الجينومية
لدى التعامل مع الجينات أو المناطق الجينومية كالمحفزات والمحسنات فاننا بطريقة أو باخرى نتعامل مع احداثيات جينومية. فمثلا يُعَرف كل جين برقم الكروموزوم الذي يتواجد فيه و احداثيات اكسوناته. تكون هذه الاحداثيات إما متصلة (تمثل منطقة واحدة من البداية إلى النهاية) او منفصلة (عبارة عن مناقص متقطعة, كاحداثيات اكزونات جين معين).
للتعامل مع هذه البيانات يمكننا استعمال حزمة GenomicRanges المبنية فوق حزمة IRanges. الفرق بين هذين الحزمتين هو أن حزمة IRanges تختص بتمثيل المجالات الرقمية بصفة عامة وليس بالضرورة المجالات الجينية, حيث يُمثل كل مجال بنقطتي بداية ونهاية, كما انها تُوفر دوال للقيام بعمليات تقاطع و دمج المجالات.عند التعامل مع المجالات الجينومية نأخذ أيضا بعين الاعتبار رقم الكروموزوم والسلسلة التي يتواجد فيها (الموجبة أو السالبة), حيث ان بداية المجال في السلسلة الموجبة تكون من الشمال إلى اليمين و العكس في حالة السلسلة السالبة.
حزمة IRanges
قبل المرور إلى حزمة GenomicRanges أرى أنه من الظروري التعريف بحزمة IRanges المصممة للتعامل مع المجالات الرقمية. تَم تصميم حزمة IRanges بطريقة تسمح بالقيام بعدة عمليات خاصة بالمجالات الرقمية كالتقاطع بين المجالات, ايجاد اقرب مجال,... إلخ باستعمال خوازميات فعالة مبرمجة في العادة بلغة الـ C++ وبطريقة تسمح باستعمال اقل للذاكرة.
تحتوي الحزمة على فئة IRanges لتمثيل مجموعة من المجالات, حيث يُمَثَلُ كل مجال بنقطتي بداية ونهاية أو بنقطة بداية و عُرض المجال. يمكننا قرائة هذه المعلومات باستعمال الدوال end, start و width كما هو موضح في المثال التالي:
library(IRanges) # يمكن انشاء مجموعة مجالات بتحديد نقطة البداية والنهاية range1 <- IRanges(start=c(10,50,300),end =c(30,90,456)) range1 #IRanges of length 3 # start end width #[1] 10 30 21 #[2] 50 90 41 #[3] 300 456 157 # يمكن انشاء مجموعة مجالات بتحديد نقطة البداية وعُرض المجال range2 <- IRanges(start=c(10,50,300),width=c(21,41,157)) range2 #IRanges of length 3 # start end width #[1] 10 30 21 #[2] 50 90 41 #[3] 300 456 157 # Start, end , widht يمكن قرائة هذه الاحداثيا باستعمال دوال start(range2) #[1] 10 50 300 end(range2) #[1] 30 90 456 width(range2) #[1] 21 41 157
من بين الدوال المُعرفة نجد دالة reduce التي تسمح بدمج المجالات المتقاطعة ودالةdisjoin التي تقوم بفصل المجالات وتعريف مناطق التقاطع على أنها مجالات جديدة. من بين الدوال المهمة أيضا لدينا دالة findOverlaps لايجاد المجالات المتقاطعة و دالة nearest لإيجاد أقرب مجال من مجال آخر, مثلا لايجاد اقرب جين من موقع ارتباط بروتين.
library(IRanges) # يمكن انشاء مجموعة مجالات بتحديد نقطة البداية والنهاية range1 <- IRanges(start=c(10,50,300),end =c(30,90,456)) range1 #IRanges of length 3 # start end width #[1] 10 30 21 #[2] 50 90 41 #[3] 300 456 157 # يمكن انشاء مجموعة مجالات بتحديد نقطة البداية وعُرض المجال range2 <- IRanges(start=c(10,50,300),width=c(21,41,157)) range2 #IRanges of length 3 # start end width #[1] 10 30 21 #[2] 50 90 41 #[3] 300 456 157 # Start, end , widht يمكن قرائة هذه الاحداثيا باستعمال دوال start(range2) #[1] 10 50 300 end(range2) #[1] 30 90 456 width(range2) #[1] 21 41 157
بالاظافة إلى الدوال والكائنات المختصة بتمثيل المجالات توفر حزمة IRanges دوال للتعامل مع الأشعة (الجداول) بطريقة احسن لايسع المجال لذكرها هنا, لكن من بين اهم الفئات المعرفة هي فئة Rleالتي تسمح بضعط البيانات باستعمال خوارزمية Run Length Encoding التي تتفادى حفظ البيانات المتكررة المتتالية مثلا إذا كان لدينا الجدول (1,1,1,2,3,3) فستقوم بحفظ الجدول (1,2,3) وعدد مرات تكرار كل عدد (3,1,2).
# ننشئ شعاع يحتوي على أرقام من 1 إلى 10 بطول 40 x<- sort(sample(1:10,40,replace=T)) head(x) #[1] 1 1 1 1 1 1 # Rle لحفضه x <- Rle(x) x #integer-Rle of length 40 with 10 runs # Lengths: 6 2 1 3 4 5 6 9 2 2 # Values : 1 2 3 4 5 6 7 8 9 10
حزمة GenomicRanges
تعتبر حزمة GenomicRanges هي الأساس الذي تقوم عليه الكثير من الحزم التي تتعامل مع الاحداثيات الجينية. تعتبر هذه الحزمة الطبقة الوسط ما بين حزمة IRanges والحزم الأخرى مثل BSgenome و ShortReads. من الجانب البرمجي توفر الحزمة الثلاث فئات الممثلة في الجدول التالي:
| GRanges | تستعمل لتمثيل احداثيات مجموعة من المناطق الجينومية (مثلا مجموعة اكسونات جين). |
| GRangesList | سِلْسلة من كائنات GRanges, مفيدة لتجميع كائنات GRanges في مجموعات (مثلا مجموعة جينات كل مع اكسوناته). |
| GAlignments | تستعمل لتمثيل نتائج المطابقة الجينومية مع جينوم مرجعي (تقبل الفراغات Gap). |
كانئات GRanges
تعتبر فئة GRanges امتدادا لفئة IRanges بحيث تسمح بتمثيل احداثيات مجموعة من المناطق الجينومية بالاظافة إلى انها تمكن المبرمج من إظافة معلومات إظافية خاصة بكل مجال. مثلا إذا كنا نتعامل مع بيانات Chip-Seq يمكننا اظافة معلومة عن عدد السلاسل المتواجدة في كل قمة,... إلخ.
يمكن استعمال دالة البناء GRanges من أجل تحديد الاحداثيات الجينومية والبيانات الوصفية (Metadata) المتعلقة بهذه البيانات. في حالة عدم تحديد البيانات الوصفية عند الانشاء يمكن استعمال دالة mcols كما هو موضح في المثال التالي:
# بعدة طرق GRanges يمكن انشاء
# مثلا يمكننا تحديد فقط المجالات الجينومية
gr <- GRanges(seqnames = c("chr1","chr1","chr2","chrX"),
ranges = IRanges(start = c(130,30050,4509,69098),
width= c(250,1300,400,590)),
strand = c("+","+","-","*"))
gr
#GRanges object with 4 ranges and 0 metadata columns:
# seqnames ranges strand
# <Rle> <IRanges> <Rle>
# [1] chr1 [ 130, 379] +
# [2] chr1 [30050, 31349] +
# [3] chr2 [ 4509, 4908] -
# [4] chrX [69098, 69687] *
# -------
# seqinfo: 3 sequences from an unspecified genome; no seqlengths
# mcols يمكن اظافة البيانات الوصفية باستعمال دالة
mcols(gr) <- DataFrame(score = runif(4))
gr
#GRanges object with 4 ranges and 1 metadata column:
# seqnames ranges strand | score
# <Rle> <IRanges> <Rle> | <numeric>
# [1] chr1 [ 130, 379] + | 0.779228801373392
# [2] chr1 [30050, 31349] + | 0.325288600986823
# [3] chr2 [ 4509, 4908] - | 0.609552363865077
# [4] chrX [69098, 69687] * | 0.229161157971248
# -------
# seqinfo: 3 sequences from an unspecified genome; no seqlengths
# او يمكننا اظافتها عند الانشاء, يمكن اعطاء اي اسم للحقول المضافة
gr <- GRanges(seqnames = c("chr1","chr1","chr2","chrX"),
ranges = IRanges(start = c(130,30050,4509,69098),
width= c(250,1300,400,590)),
strand = c("+","+","-","*"),
score = runif(4))
gr
#GRanges object with 4 ranges and 1 metadata column:
# seqnames ranges strand | score
# <Rle> <IRanges> <Rle> | <numeric>
# [1] chr1 [ 130, 379] + | 0.0297835641540587
# [2] chr1 [30050, 31349] + | 0.42395940516144
# [3] chr2 [ 4509, 4908] - | 0.112338428385556
# [4] chrX [69098, 69687] * | 0.800393168348819
# -------
# seqinfo: 3 sequences from an unspecified genome; no seqlengths
يمكن التعامل مع هذا الكائن كأنك تتعامل مع جدول, فيمكنك مثلا قراءة اي مجال باستعمال [] أما اذا اردت معرفة نقاط البداية أو النهاية أو اسم الكروموزومات التي تتواجد فيها هذه المجالات يمكنك استعمال دوالend , start و seqnames ولقراءة البيانات الوصفية يمكنك استعمال دالة elementMetadata أو values :
gr[1:2] #GRanges object with 2 ranges and 1 metadata column: # seqnames ranges strand | score # <Rle> <IRanges> <Rle> | <numeric> # [1] chr1 [ 130, 379] + | 0.0297835641540587 # [2] chr1 [30050, 31349] + | 0.42395940516144 # ------- # seqinfo: 3 sequences from an unspecified genome; no seqlengths start(gr) #[1] 130 30050 4509 69098 end(gr) #[1] 379 31349 4908 69687 width(gr) #[1] 250 1300 400 590 seqnames(gr) #factor-Rle of length 4 with 3 runs # Lengths: 2 1 1 # Values : chr1 chr2 chrX #Levels(3): chr1 chr2 chrX values(gr) #DataFrame with 4 rows and 1 column # score # <numeric> #1 0.02978356 #2 0.42395941 #3 0.11233843 #4 0.80039317
في العادة تكون هذه البيانات في ملفات BED او GTF يمكنك مثلا قراءة الـ BED على شكل data.frame ثم انشاء كائنGRangesيدويا أو يمكن استعمال دالة import في حزمة rtracklayer. في المثال التالي سوف نقوم بقراءة ملف BED x باستعمال دالة import و باستعمال الطريقة الكلاسيكية بدون الحاجة إلى حزم اخرى (للفضولين البيانات المستعملة هي الـ 1 بالمئة من الجينوم البشري التي تم دراستها في المشروع التجريبي لـ ENCODE, لمزيد من المعلومات هنا)
# rtracklayer يمكن استعمال الدوال الجاهزة في حزمة
encodePilotRegions <- import.bed("encodePilotRegion.bed")
encodePilotRegions
#GRanges object with 44 ranges and 1 metadata column:
# seqnames ranges strand | name
# <Rle> <IRanges> <Rle> | <character>
# [1] chr1 [151158061, 151658060] * | ENr231
# [2] chr2 [ 51658705, 52158704] * | ENr112
# [3] chr2 [118294574, 118794573] * | ENr121
# [4] chr2 [220277346, 220777345] * | ENr331
# [5] chr2 [234491825, 234991888] * | ENr131
# ... ... ... ... ... ...
# [40] chr19 [54331773, 55332648] * | ENm007
# [41] chr20 [33841513, 34341514] * | ENr333
# [42] chr21 [33746366, 35442351] * | ENm005
# [43] chr21 [40322597, 40822596] * | ENr133
# [44] chr22 [31803954, 33503953] * | ENm004
# -------
# seqinfo: 21 sequences from an unspecified genome; no seqlengths
# أو يمكن تحضير كل شيء يدويا
df <- read.table("encodePilotRegion.bed",sep="\t")
head(df)
# V1 V2 V3 V4
#1 chr1 151158060 151658060 ENr231
#2 chr2 51658704 52158704 ENr112
#3 chr2 118294573 118794573 ENr121
#4 chr2 220277345 220777345 ENr331
#5 chr2 234491824 234991888 ENr131
#6 chr4 118246655 118746655 ENr113
encodePilotRegions2
#GRanges object with 44 ranges and 1 metadata column:
# seqnames ranges strand | name
# <Rle> <IRanges> <Rle> | <factor>
# [1] chr1 [151158060, 151658060] * | ENr231
# [2] chr2 [ 51658704, 52158704] * | ENr112
# [3] chr2 [118294573, 118794573] * | ENr121
# [4] chr2 [220277345, 220777345] * | ENr331
# [5] chr2 [234491824, 234991888] * | ENr131
# ... ... ... ... ... ...
# [40] chr19 [54331772, 55332648] * | ENm007
# [41] chr20 [33841512, 34341514] * | ENr333
# [42] chr21 [33746365, 35442351] * | ENm005
# [43] chr21 [40322596, 40822596] * | ENr133
# [44] chr22 [31803953, 33503953] * | ENm004
# -------
# seqinfo: 21 sequences from an unspecified genome; no seqlengths
يمكن القيام بعدة عمليات تقاطع أو اظافة مجالات أو غيرها من العمليات على كائناتGRangesالتي يمكن اختصارها في البيان اسفله. هنا لدينا المتغير Gene model الذي يمثل جين ذو منسوخين, حيث مثلنا كل منسوخ بكائنGRanges ووضعناهم في كائنGRangesList. كل العمليات الموضحة في البيان يمكن تطبيقها على كائناتGRanges
رسم تخطيطي يوضح الدوال التي يمكن استعمالها للتعامل مع كائنات GRanges (المصدر)
كانئات GAlignments
الهدف من هذه الفئة هو توفير بيئة للتعامل مع نتائج المطابقة (mapping أو alignment) التي تكون في العادة على شكل ملفات BAM او SAM بطريقة سَلِسَة واستعمال ذكي لذاكرة الكمبيوتر. هناك عدة حزم لقرائة ملفات نتائج المطابقة في هذا المثال سوف نستعمل دالة readGAlignmentsFromBam من حزمة GenomicAlignments. هنا قمنا بتحميل ملف BAM صغير (يحتوي على 204355 سِلْسِلة) بالكامل إلى الذاكرة لكن في العادة نقوم بقرائته على أجزاء أو نقراء منطقة معينة, سنتطرق إلى هذه التفاصيل في بقية المقال. استعملنا الحزمة pasillaBamSubset هنا لاحتوائها على امثلة ملفات BAM ماخوذة من نتائج البحث من طرف Brooks حول ذبابة الفاكهة.
require(pasillaBamSubset) require(GenomicAlignments) # BAM نقراء الملف bamFile <- untreated1_chr4() gReads <- readGAlignmentsFromBam(bamFile) #GAlignments object with 204355 alignments and 0 metadata columns: # seqnames strand cigar qwidth start end width njunc # <Rle> <Rle> <character> <integer> <integer> <integer> <integer> <integer> # [1] chr4 - 75M 75 892 966 75 0 # [2] chr4 - 75M 75 919 993 75 0 # [3] chr4 + 75M 75 924 998 75 0 # ... ... ... ... ... ... ... ... ... # [204353] chr4 + 75M 75 1348268 1348342 75 0 # [204354] chr4 - 75M 75 1348449 1348523 75 0 # [204355] chr4 - 75M 75 1350124 1350198 75 0 # ------- # seqinfo: 8 sequences from an unspecified genome
يحتوي كل سطر من كائناتGAlignments على معلومات حول بعض الحقول المهمة من ملفات SAM أو BAM كحقل CIGAR الذي يعطي معلومات حول عدد الجزئيات القاعدية التي تَم أو لم يتم مطابقتها وعدد الفراغات,... إلخ من المعلومات. يمكن الاطلاع على معنى كل حقل في الرابط التالي.
التعامل مع ملفات المطابقة
تعتبر ملفات BAM و SAM من بين الملفات التي يُتَعَامل معها بصفة روتينية في مجال المعلوماتية الحيوية, حيث أنها تحتوي على معلومات حول احداثيات المناطق المُطابقة وجودتها وكيفيتها. في الكثيرمن الاحيان يمكن ان نكتفي بالأدوات الموجودة لكن في بعض الاحيان وخاصة عند التعامل مع بيانات جديدة نحاول كتابة سكريبتات للقيام بذلك. يتوفر مستودع Bioconductor على مجموعة من الحزم التي تساعدنا في التعامل مع هذا النوع من الملفات. في هذا المقال سوف نكتفي بتقديم بعض الحزم المهمة مثل حزمة GenomicAlignments و حزمتي Rsamtools و ShortRead.
حزمة Rsamtools
تعتبر حزمة Rsamtools واجهة للدوال البرامج samtools , bcftool و tabix للتعامل مع ملفات BAM, SAM, BCF و TABIX. من اهم الفئات التي توفرها حزمة Rsamtools هي فئة ScanBamParam التي تسمح بتحديد الحقول والمجالات الجينومية التي نريد قرائتها, ففي حالة توفر ملفات bam indexذات الامتداد bai يمكننا قرائة ملفات كبيرة بسرعة و بدون استهلاك للذاكرة. مثلا في هذا المثال سوف نقوم بقرائة منطقة جينومية من ملف BAM.
require(pasillaBamSubset)
require(Rsamtools)
bamFile <- untreated1_chr4()
region <- GRanges("chr4",IRanges(100000L, 110000L))
param <- ScanBamParam(which=region, what=scanBamWhat())
mappedReads <- scanBam(bamFile, param=param)
يعتبر حقل FLAG من بين اهم الحقول حيث انه يخبرنا عن نوعيةكل مطابقة (هل تطابقت السلسلة مع منطقة واحدة ام لا, هل السِلْسِلة احادية أم نناية الحد,...إلخ ), يمكنا طباعة ملخص عن بعض الاحصائيات حقل FLAG باستعمال الدالة quickBamFlagSummary :
require(RNAseqData.HNRNPC.bam.chr14) require(pasillaBamSubset) require(Rsamtools) bamFile <- untreated1_chr4() # نلاحظ أن المجموعة اس هي المتواجدة هنا مما يدل على ان السلاسل احادية النهاية quickBamFlagSummary(bamFile, main.groups.only = TRUE) # group | nb of | nb of | mean / max # of | records | unique | records per # records | in group | QNAMEs | unique QNAME #All records........................ A | 204355 | 190770 | 1.07 / 10 # o template has single segment.... S | 204355 | 190770 | 1.07 / 10 <--- # o template has multiple segments. M | 0 | 0 | NA / NA # - first segment.............. F | 0 | 0 | NA / NA # - last segment............... L | 0 | 0 | NA / NA # - other segment.............. O | 0 | 0 | NA / NA # #Note that (S, M) is a partitioning of A, and (F, L, O) is a partitioning of M. #Indentation reflects this. # نلاحظ هنا أن السلاسل ثنائية النهاية bamfile2 <- RNAseqData.HNRNPC.bam.chr14_BAMFILES[1] quickBamFlagSummary(bamfile2, main.groups.only = TRUE) # group | nb of | nb of | mean / max # of | records | unique | records per # records | in group | QNAMEs | unique QNAME #All records........................ A | 800484 | 393300 | 2.04 / 10 # o template has single segment.... S | 0 | 0 | NA / NA # o template has multiple segments. M | 800484 | 393300 | 2.04 / 10 <-- # - first segment.............. F | 400242 | 393300 | 1.02 / 5 # - last segment............... L | 400242 | 393300 | 1.02 / 5 # - other segment.............. O | 0 | 0 | NA / NA # #Note that (S, M) is a partitioning of A, and (F, L, O) is a partitioning of M. #Indentation reflects this.
من بين الدوال المفيدة أيضا لدينا دلة BamSampler التي تقوم بالاختيار العشوائي للسلاسل من ملف BAM, يمكن استعمالها مثلا للتأكد من مدى كفاية عمق عملية السَلْسَلَة (sequencing depth) ورسم بيان مشابه للرسم اسفله الذي يين عدد الجينات التي تم مطابقتها في كل مرة ناخذ جزء من نتائج المطابقة. في حالة ما إذا كانت عملية السَلْسَلة بكمية كافية فانه عند حد معين (مثلا 80%) لن نلاحظ زيادة ملوحة في عدد الجينات المكتشفة مما يدل على اننا غطينا الجينوم بشكل كامل, اما حالة الحصول على منحنى متصاعد إلى غاية 100% فتدل ان السَلْسَلَة لم تكن بالعمق الكاففي.
رسم تخطيطي يوضح مخططات التشبع (Saturation plot) (المصدر)
إذا كان حجم الملف كبيرا يمكنك قرائته على اجزاء بايستعمال دالة BamFile وتحديد عدد السلاسل التي تريد قرائها في الحقل: yieldSize
حزمة GenomicAlignments
كما نوهنا إليه سابقا, توفر هذه الحزمة مجموعة من الدوال لقرائة ملفات BAM, تحتوي هذه الجزمة على ثلات فئات رئيسية:
| GAlignments | تستعمل لتخرين نتائج السلسلة على شكل جدول. في حالة السلاسلة المتعددة مناطق المطابقة, تُمثل كل عن نتيجة مُطَابقة على انها مستقلة عن الاخرى لكن يشتركون في الاسم |
| GAlignmentPairs | تستعمل لتخزين نتائج السلسلة ثنائية النهاية (pair-end). يمكن استعمال دالتيleft و right لقرائة النهايتين. |
| GAlignmentsList | جدول من كائنات GAlignments . |
يمكن ايضا استعمال ScanBamParam لتحديد المنطقة التي تريد دراستها. مثلا يمكنك قرائة السلاسل المطابقة لكل كروموزموم في كائن وتحليلهابطريقة متوازية.
البيان التالي يلخص أهم الحزم المستعملة في مختلف مراحل التعامل مع نتائج السَلْسَلة.
رسم تخطيطي يوضح اهم حزم Bioconductor المستعملة في تحليل نتائج السَلْسَلة (المصدر)
وفي الأخير نتمنى ان يكون هذا المقال مفيدا في التعريف ببعض خصائص الحزم واعطاء نظرة عن امكانيات آر في التعامل مع السلاسل الجينومية. بطبيعة الحال المقال على طوله لم يتطرق إلى الكثير من الدوال ونترك المجال للقارئ في اكتشافا أو سؤالنا عبر التعليقات.
مصادر:
- Lawrence M, Huber W, Pagès H, Aboyoun P, Carlson M, Gentleman R, Morgan M and Carey V (2013). “Software for Computing and Annotating Genomic Ranges.” PLoS Computational Biology, 9.
- Morgan M, Pagès H, Obenchain V and Hayden N. Rsamtools: Binary alignment (BAM), FASTA, variant call (BCF), and tabix file import. R package version 1.18.2,
- Bioconductor for Sequence Analysis
- SAM/BAM and related specifications
- Angela R W et al, "Quantitative assessment of single-cell RNA-sequencing methods", 2013, Nature Methods.