تحديد الربط البديل والبولي أدينيل في بيانات RNA-seq

ملخص

يعمل الربط البديل (AS) والبولي أدينيل البديل (APA) على توسيع تنوع الأشكال المتماثلة ومنتجاتها. هنا ، نصف بروتوكولات المعلوماتية الحيوية لتحليل مقايسات تسلسل RNA-seq السائبة و 3 'للكشف عن وتصور AS و APA المتغيرة عبر الظروف التجريبية.

الملخص

بالإضافة إلى التحليل النموذجي ل RNA-Seq لقياس التعبير الجيني التفاضلي (DGE) عبر الظروف التجريبية / البيولوجية ، يمكن أيضا استخدام بيانات RNA-seq لاستكشاف آليات تنظيمية معقدة أخرى على مستوى exon. يلعب الربط البديل و polyadenylation دورا حاسما في التنوع الوظيفي للجين من خلال توليد أشكال متساوية مختلفة لتنظيم التعبير الجيني على مستوى ما بعد النسخ ، ويمكن أن يؤدي قصر التحليلات على مستوى الجينات بالكامل إلى تفويت هذه الطبقة التنظيمية المهمة. هنا ، نوضح تحليلات مفصلة خطوة بخطوة لتحديد وتصور استخدام موقع exon و polyadenylation التفاضلي عبر الظروف ، باستخدام Bioconductor والحزم والوظائف الأخرى ، بما في ذلك DEXSeq و diffSplice من حزمة Limma و rMATS.

المقدمة

تم استخدام RNA-seq على نطاق واسع على مر السنين عادة لتقدير التعبير الجيني التفاضلي واكتشاف الجينات¹. بالإضافة إلى ذلك ، يمكن استخدامه أيضا لتقدير الاستخدام المتغير لمستوى إكسون بسبب تعبير الجينات عن أشكال متساوية مختلفة ، وبالتالي المساهمة في فهم أفضل لتنظيم الجينات على مستوى ما بعد النسخ. تولد غالبية الجينات حقيقية النواة أشكالا متساوية مختلفة عن طريق الربط البديل (AS) لزيادة تنوع تعبير mRNA. يمكن تقسيم أحداث AS إلى أنماط مختلفة: تخطي exons كاملة (SE) حيث تتم إزالة exon ("كاسيت") تماما من النص جنبا إلى جنب مع الإنترونات المرافقة له ؛ اختيار موقع لصق بديل (مانح) 5 '(A5SS) واختيار موقع لصق بديل 3 '(متقبل) (A3SS) عند وجود موقعين أو أكثر من مواقع لصق على طرفي إكسون ؛ الاحتفاظ بالإنترونات (RI) عند الاحتفاظ ب intron ضمن نسخة mRNA الناضجة والاستبعاد المتبادل لاستخدام exon (MXE) حيث يمكن الاحتفاظ بواحد فقط من الإكسونات المتاحين في وقت ^2,3. يلعب polyadenylation البديل (APA) أيضا دورا مهما في تنظيم التعبير الجيني باستخدام مواقع poly البديلة (A) لتوليد أشكال متماثلة متعددة من mRNA من نسخة واحدة⁴. تقع معظم مواقع polyadenylation (pAs) في المنطقة غير المترجمة 3 '(3' UTRs) ، مما يولد أشكالا متماثلة من mRNA بأطوال UTR متنوعة 3 بوصات. نظرا لأن 3 'UTR هو المحور المركزي للتعرف على العناصر التنظيمية ، يمكن أن تؤثر أطوال UTR المختلفة 3 'على توطين mRNA واستقراره وترجمته⁵. هناك فئة من 3 'مقايسات تسلسل نهاية محسنة للكشف عن APA التي تختلف في تفاصيل البروتوكول⁶. تم تصميم خط الأنابيب الموصوف هنا ل PolyA-seq ، ولكن يمكن تكييفه مع البروتوكولات الأخرى كما هو موضح.

في هذه الدراسة ، نقدم مجموعة من طرق تحليل exon التفاضلية ^7,8 (الشكل 1) ، والتي يمكن تقسيمها إلى فئتين عريضتين: القائمة على exon (DEXSeq⁹ ، diffSplice 10) والقائمة على الحدث (تكرار التحليل متعدد المتغيرات لربط النص (rMATS)¹¹). تقارن الطرق القائمة على الإكسون تغير الطي عبر ظروف الإكسونات الفردية ، مقابل مقياس للتغير الكلي في طية الجينات لاستدعاء استخدام إكسون المعبر عنه بشكل تفاضلي ، ومن ذلك تحسب مقياسا على مستوى الجينات لنشاط AS. تستخدم الطرق المستندة إلى الأحداث قراءات تقاطع exon-intron-spanning لاكتشاف وتصنيف أحداث الربط المحددة مثل تخطي exon أو الاحتفاظ بالإنترونات ، وتمييز أنواع AS هذه في الإخراج³. وبالتالي ، توفر هذه الطرق وجهات نظر تكميلية لتحليل كامل لمعيار المحاسبة^12,13. اخترنا DEXSeq (استنادا إلى حزمة DESeq2¹⁴ DGE) و diffSplice (استنادا إلى حزمة Limma¹⁰ DGE) للدراسة لأنهما من بين الحزم الأكثر استخداما لتحليل الربط التفاضلي. تم اختيار rMATS كطريقة شائعة للتحليل القائم على الأحداث. طريقة أخرى شائعة قائمة على الأحداث هي MISO (خليط من الأشكال المتساوية)¹. بالنسبة إلى APA ، نقوم بتكييف النهج القائم على exon.

الشكل 1. خط أنابيب التحليل. مخطط انسيابي للخطوات المستخدمة في التحليل. تشمل الخطوات: الحصول على البيانات ، وإجراء فحوصات الجودة ومحاذاة القراءة متبوعة بحساب القراءات باستخدام التعليقات التوضيحية لمواقع exons و introns و pA المعروفة ، والتصفية لإزالة الأعداد المنخفضة والتطبيع. تم تحليل بيانات PolyA-seq لمواقع pA البديلة باستخدام طرق diffSplice / DEXSeq ، وتم تحليل RNA-Seq السائب للربط البديل على مستوى exon باستخدام طرق diffSplice / DEXseq ، وتم تحليل أحداث AS باستخدام rMATS. الرجاء الضغط هنا لعرض نسخة أكبر من هذا الشكل.

تم الحصول على بيانات RNA-seq المستخدمة في هذا المسح من التعبير الجيني الجامع (GEO) (GSE138691)¹⁵. استخدمنا بيانات RNA-seq للفأر من هذه الدراسة مع مجموعتين من الحالات: النوع البري (WT) والضربة القاضية من النوع 1 الشبيه بالعضلات (Mbnl1 KO) مع ثلاث نسخ متماثلة لكل منهما. لإثبات تحليل استخدام موقع polyadenylation التفاضلي ، حصلنا على بيانات الخلايا الليفية الجنينية للفأر (MEFs) PolyA-seq (GEO Accession GSE60487)¹⁶. تحتوي البيانات على أربع مجموعات شروط: النوع البري (WT) ، النوع 1 / النوع 2 الشبيه بالعضلات بالضربة القاضية المزدوجة (Mbnl1/2 DKO) ، Mbnl 1/2 DKO مع ضربة قاضية Mbnl3 (KD) و Mbnl1/2 DKO مع التحكم Mbnl3 (Ctrl). تتكون كل مجموعة شرط من نسختين متماثلتين.

	انضمام GEO	رقم تشغيل SRA	اسم العينة	شرط	تكرار	نسيج	التسلسل	طول القراءة
RNA-Seq	جي إس إم 4116218	SRR10261601	Mbnl1KO_Thymus_1	Mbnl1 بالضربة القاضية	ممثل 1	الغدة الصعترية	نهاية مزدوجة	100 نقطة أساس
	جي إس إم 4116219	SRR10261602	Mbnl1KO_Thymus_2	Mbnl1 بالضربة القاضية	ممثل 2	الغدة الصعترية	نهاية مزدوجة	100 نقطة أساس
	جي إس إم 4116220	SRR10261603	Mbnl1KO_Thymus_3	Mbnl1 بالضربة القاضية	مندوب 3	الغدة الصعترية	نهاية مزدوجة	100 نقطة أساس
	جي إس إم 4116221	SRR10261604	WT_Thymus_1	نوع البرية	ممثل 1	الغدة الصعترية	نهاية مزدوجة	100 نقطة أساس
	جي إس إم 4116222	SRR10261605	WT_Thymus_2	نوع البرية	ممثل 2	الغدة الصعترية	نهاية مزدوجة	100 نقطة أساس
	جي إس إم 4116223	SRR10261606	WT_Thymus_3	نوع البرية	مندوب 3	الغدة الصعترية	نهاية مزدوجة	100 نقطة أساس
3P-Seq	جي إس إم 1480973	ريال1553129	WT_1	النوع البري (WT)	ممثل 1	الخلايا الليفية الجنينية للفأر (MEFs)	نهاية واحدة	40 نقطة أساس
	جي إس إم 1480974	ريال1553130	WT_2	النوع البري (WT)	ممثل 2	الخلايا الليفية الجنينية للفأر (MEFs)	نهاية واحدة	40 نقطة أساس
	جي إس إم 1480975	SRR1553131	DKO_1	Mbnl 1/2 ضربة قاضية مزدوجة (DKO)	ممثل 1	الخلايا الليفية الجنينية للفأر (MEFs)	نهاية واحدة	40 نقطة أساس
	جي إس إم 1480976	SRR1553132	DKO_2	Mbnl 1/2 ضربة قاضية مزدوجة (DKO)	ممثل 2	الخلايا الليفية الجنينية للفأر (MEFs)	نهاية واحدة	40 نقطة أساس
	جي إس إم 1480977	ريال1553133	DKOsiRNA_1	Mbnl 1/2 ضربة قاضية مزدوجة مع Mbnl 3 siRNA (دينار كويتي)	ممثل 1	الخلايا الليفية الجنينية للفأر (MEFs)	نهاية واحدة	40 نقطة أساس
	جي إس إم 1480978	ريال1553134	DKOsiRNA_2	Mbnl 1/2 ضربة قاضية مزدوجة مع Mbnl 3 siRNA (دينار كويتي)	ممثل 2	الخلايا الليفية الجنينية للفأر (MEFs)	نهاية واحدة	36 نقطة أساس
	جي إس إم 1480979	ريال1553135	DKONTsiRNA_1	Mbnl 1/2 ضربة قاضية مزدوجة مع siRNA غير مستهدف (Ctrl)	ممثل 1	الخلايا الليفية الجنينية للفأر (MEFs)	نهاية واحدة	40 نقطة أساس
	جي إس إم 1480980	SRR1553136	DKONTsiRNA_2	Mbnl 1/2 ضربة قاضية مزدوجة مع siRNA غير مستهدف (Ctrl)	ممثل 2	الخلايا الليفية الجنينية للفأر (MEFs)	نهاية واحدة	40 نقطة أساس

الجدول 1. ملخص مجموعات بيانات RNA-Seq و PolyA-seq المستخدمة في التحليل.

البروتوكول

1. تركيب الأدوات وحزم R المستخدمة في التحليل

1. Conda هو مدير حزم شائع ومرن يسمح بالتثبيت المريح للحزم مع تبعياتها عبر جميع الأنظمة الأساسية. استخدم "Anaconda" (مدير حزم conda) لتثبيت "conda" والذي يمكن استخدامه لتثبيت الأدوات / الحزم المطلوبة للتحليل.
2. قم بتنزيل "Anaconda" وفقا لمتطلبات النظام من https://www.anaconda.com/products/individual#Downloads وقم بتثبيته باتباع المطالبات في برنامج التثبيت الرسومي. قم بتثبيت جميع الحزم المطلوبة باستخدام "conda" عن طريق كتابة ما يلي في سطر أوامر Linux.
   conda install -c daler sratoolkit
conda install -c conda-forge parallel
conda install -c bioconda star bowtie fastqc rmats rmats2sashimiplot samtools fasterq-dump cutadapt bedtools deeptools
3. لتنزيل جميع حزم R المستخدمة في البروتوكول ، اكتب التعليمة البرمجية التالية في وحدة تحكم R (التي بدأت في سطر أوامر Linux بكتابة "R") أو وحدة تحكم Rstudio.
   bioc_packages<- c("DEXSeq", "Rsubread", "EnhancedVolcano", "edgeR", "limma", "maser","GenomicRanges")
packages<- c("magrittr", "rtracklayer", "tidyverse", "openxlsx", "BiocManager")
#Install if not already installed
installed_packages<-packages%in% rownames(installed.packages())
installed_bioc_packages<-bioc_packages%in% rownames(installed.packages())
if(any(installed_packages==FALSE)) {
install.packages(packages[!installed_packages],dependencies=TRUE)
BiocManager::install(packages[!installed_bioc_packages], dependencies=TRUE)
}
   ملاحظة: في هذا البروتوكول الحسابي ، سيتم إعطاء الأوامر إما كملفات R Notebook (ملفات ذات امتداد ". Rmd") ، ملفات رمز R (ملفات ذات امتداد ". R") ، أو البرامج النصية Linux Bash shell (الملفات ذات الامتداد ".sh"). يجب فتح ملفات R Notebook (Rmd) في RStudio باستخدام File| افتح ملف...، ثم يتم تشغيل أجزاء التعليمات البرمجية الفردية (التي قد تكون أوامر R أو أوامر Bash shell) بشكل تفاعلي بالنقر فوق السهم الأخضر في الجزء العلوي الأيمن. يمكن تشغيل ملفات التعليمات البرمجية R عن طريق الفتح في RStudio ، أو على سطر أوامر Linux عن طريق التقديم ب "Rscript" ، على سبيل المثال Rscript example.R. يتم تشغيل البرامج النصية Shell على سطر أوامر Linux عن طريق تمهيد البرنامج النصي بالأمر "sh" على سبيل المثال .sh example.sh.

2. تحليل الربط البديل (AS) باستخدام RNA-seq

1. تنزيل البيانات ومعالجتها مسبقا
   ملاحظة: تتوفر مقتطفات التعليمات البرمجية المشروحة أدناه في ملف التعليمات البرمجية التكميلية "AS_analysis_RNASeq.Rmd"، لاتباع الخطوات الفردية بشكل تفاعلي ، ويتم توفيرها أيضا كبرنامج نصي bash ليتم تشغيله دفعة واحدة على سطر أوامر Linux (ش downloading_data_preprocessing.sh).
   1. تنزيل البيانات الأولية.
      1. قم بتنزيل البيانات الأولية من أرشيف قراءة التسلسل (SRA) باستخدام الأمر "الجلب المسبق" من مجموعة أدوات SRA (الإصدار 2.10.8)¹⁷. أعط معرفات انضمام SRA بالتسلسل في الأمر التالي لتنزيلها بالتوازي باستخدام الأداة المساعدة المتوازية^{GNU 18}. لتنزيل ملفات SRA لمعرفات الانضمام من SRR10261601 إلى SRR10261606 بالتوازي ، استخدم ما يلي في سطر أوامر Linux.
         ​seq 10261601 10261606 | parallel prefetch SRR{}
      2. استخرج ملفات fastq من الأرشيف باستخدام وظيفة "fastq-dump" من مجموعة أدوات SRA. استخدم غنو بالتوازي وأعط أسماء جميع ملفات SRA معا.
         ​parallel -j 3 fastq-dump --gzip --skip-technical --read-filter pass --dumpbase --split-e --clip --origfmt {} :::
      3. قم بتنزيل الجينوم المرجعي والتعليقات التوضيحية للماوس (تجميع الجينوم GRCm39) من www.ensembl.org باستخدام ما يلي في سطر أوامر Linux.
         wget -nv -O annotation.gtf.gz http://ftp.ensembl.org/pub/release-103/gtf/mus_musculus/Mus_musculus.GRCm39.103.gtf.gz \ && gunzip -f annotation.gtf.gz
wget -nv -O genome.fa.gz http://ftp.ensembl.org/pub/release-103/fasta/mus_musculus/dna/Mus_musculus.GRCm39.dna.primary_assembly.fa.gz \ && gunzip -f genome.fa.gz
GTF=$(readlink -f annotation.gtf)
GENOME=$(readlink -f genome.fa)
   2. تتم قراءة المعالجة المسبقة ورسم الخرائط لتجميع الجينوم
      1. مراقبة الجودة. قم بتقييم جودة القراءات الأولية باستخدام FASTQC (فحص جودة FASTQ v0.11.9)¹⁹. قم بإنشاء مجلد إخراج وتشغيل fastqc بالتوازي على ملفات fasta متعددة الإدخال. ستنشئ هذه الخطوة تقرير جودة لكل عينة. افحص التقارير للتأكد من أن جودة القراءات مقبولة قبل إجراء مزيد من التحليل. (راجع دليل المستخدم لفهم التقارير في https://www.bioinformatics.babraham.ac.uk/projects/fastqc/)
         mkdir fastqc_out
parallel "fastqc {} -o fastqc_out" ::: $RAW_DATA/*.fastq.gz
         ملاحظة: إذا لزم الأمر ، قم بإجراء تشذيب المحول باستخدام "cutadapt'20 أو "trimmomatic"²¹ لإزالة التسلسل إلى محولات جانبية ، والتي تختلف بناء على حجم جزء RNA وطول القراءة. في هذا التحليل ، تخطينا هذه الخطوة لأن جزء القراءات المتأثرة كان ضئيلا.
      2. قراءة المحاذاة. تتضمن الخطوة التالية في المعالجة المسبقة تعيين القراءات إلى الجينوم المرجعي. أولا ، قم ببناء فهرس الجينوم المرجعي باستخدام وظيفة "genomeGenerate" ل STAR²² ثم قم بمحاذاة القراءات الأولية مع المرجع (بدلا من ذلك ، تتوفر الفهارس المعدة مسبقا من موقع STAR ويمكن استخدامها مباشرة للمحاذاة). قم بتشغيل الأوامر التالية في سطر أوامر Linux.
         #Build STAR index
GDIR=STAR_indices
mkdir $GDIR
STAR --runMode genomeGenerate --genomeFastaFiles $GENOME --sjdbGTFfile $GTF --runThreadN 8 --genomeDir $GDIR
ODIR=results/mapping
mkdir -p $ODIR
#Align reads to the genome
for fq1 in $RAW_DATA/*R1.fastq.gz;
do
fq2=$(echo $fq1| sed 's/1.fastq.gz/2.fastq.gz/g');
OUTPUT=$(basename ${fq1}| sed 's/R1.fastq.gz//g');
STAR --genomeDir $GDIR \
--runThreadN 12 \
--readFilesCommand zcat \
--readFilesIn ${fq1}${fq2}\
--outFileNamePrefix $ODIR\/${OUTPUT} \
--outSAMtype BAM SortedByCoordinate \
--outSAMunmapped Within \
--outSAMattributes Standard
Done
         ملاحظة: سيقوم محاذاة STAR بإنشاء وفرز ملفات BAM (خريطة المحاذاة الثنائية) لكل عينة بعد محاذاة القراءة. يجب فرز ملفات Bam قبل المتابعة إلى خطوات أخرى.
2. إعداد شروح إكسون.
   1. قم بتشغيل ملف التعليمات البرمجية التكميلية "prepare_mm_exon_annotation. R" مع التعليق التوضيحي الذي تم تنزيله بتنسيق GTF (تنسيق نقل الجينات) لإعداد التعليقات التوضيحية. للتشغيل، اكتب ما يلي في سطر أوامر Linux.
      Rscript prepare_mm_exon_annotation.R annotation.gtf
      ملاحظة: يحتوي ملف GTF إدخالات إكسون متعددة لأشكال isoforms مختلفة. يستخدم هذا الملف "لطي" معرفات النصوص المتعددة لكل exon. إنها خطوة مهمة لتحديد صناديق عد إكسون.
3. العد يقرأ. الخطوة التالية هي حساب عدد القراءات المعينة إلى نصوص / إكسونات مختلفة. انظر الملف التكميلي: "AS_analysis_RNASeq.Rmd".
   1. تحميل المكتبات المطلوبة:
      packages<- c("Rsubread","tidyverse", "magrittr", "EnhancedVolcano", "edgeR","openxlsx")
invisible(lapply(packages, library, character.only=TRUE))
   2. قم بتحميل ملف التعليقات التوضيحية الذي تمت معالجته والذي تم الحصول عليه من الخطوة السابقة (2.2).
      load("mm_exon_anno.RData")
   3. اقرأ جميع ملفات bam التي تم الحصول عليها في الخطوة 2.2.2 كإدخال ل "featureCounts" لحساب القرعات. اقرأ المجلد الذي يحتوي على ملفات bam عن طريق سرد كل ملف أولا من الدليل الذي ينتهي ب .bam. استخدم "featureCounts" من حزمة Rsubread التي تأخذ ملفات bam والتعليق التوضيحي GTF المعالج (المرجع) كمدخل لإنشاء مصفوفة من الأعداد المرتبطة بكل ميزة مع صفوف تمثل exons (الميزات) والأعمدة التي تمثل العينات.
      countData <- dir("bams", pattern=".bam$", full.names=T) %>%
featureCounts(annot.ext=anno,
isGTFAnnotationFile=FALSE,
minMQS=0,useMetaFeatures=FALSE,
allowMultiOverlap=TRUE,
largestOverlap=TRUE,
countMultiMappingReads=FALSE,
primaryOnly=TRUE,
isPairedEnd=TRUE,
nthreads=12)
   4. بعد ذلك ، قم بإجراء تصفية غير محددة لإزالة exons المعبر عنها بشكل متواضع ("غير محدد" يشير إلى أن معلومات الحالة التجريبية لا تستخدم في التصفية ، لتجنب تحيزات الاختيار). قم بتحويل البيانات من مقياس أولي إلى عدد لكل مليون (cpm) باستخدام وظيفة cpm من حزمة 'edgeR'²³ واحتفظ ب exons بأعداد أكبر من عتبة قابلة للتعيين (يتم استخدام نسخة واحدة لكل ألف ظهور لمجموعة البيانات هذه) في ثلاث عينات على الأقل. أيضا إزالة الجينات مع إكسون واحد فقط.
      # Non-specific filtering: Remove the exons with low counts
isexpr<- rownames(countData$counts)[rowSums(cpm(countData$counts)>1) >=3]
countData$counts<-countData$counts[rownames(countData$counts) %in%isexpr, ]
anno<-anno%>% filter(GeneID%in% rownames(countData$counts))
# Remove genes with only 1 site and NA in geneIDs
dn<-anno%>%group_by(GeneID)%>%summarise(nsites=n())%>% filter(nsites>1&!is.na(GeneID))
anno<-anno%>% filter(GeneID%in%dn$GeneID)
countData$counts<-countData$counts[rownames(countData$counts) %in%anno$GeneID, ]
      ملاحظة: تحقق من المعلمات المطلوبة ل featureCounts عند استخدام بيانات مختلفة ، على سبيل المثال ، للقراءات أحادية الطرف ، قم بتعيين 'isPairedEnd = FALSE'. ارجع إلى دليل مستخدم RSubread لاختيار خيارات بياناتك، وراجع قسم المناقشة أدناه.
4. الربط التفاضلي وتحليل استخدام إكسون. وصفنا بديلين لهذه الخطوة: DEXSeq و DiffSplice. إما يمكن استخدامها وإعطاء نتائج مماثلة. للتناسق ، حدد DEXSeq إذا كنت تفضل حزمة DESeq2 ل DGE واستخدم DiffSplice لتحليل DGE المستند إلى Limma. انظر الملف التكميلي: "AS_analysis_RNASeq.Rmd".
   1. استخدام حزمة DEXSeq لتحليل إكسون التفاضلي.
      1. قم بتحميل المكتبة وإنشاء جدول عينة لتحديد التصميم التجريبي.
         library(DEXSeq)
sampleTable<-data.frame(row.names= c("Mbnl1KO_Thymus_1", "Mbnl1KO_Thymus_2", "Mbnl1KO_Thymus_3", "WT_Thymus_1", "WT_Thymus_2", "WT_Thymus_3"), condition= rep(c("Mbnl1_KO", "WT"),c(3,3)), libType= rep(c("paired-end")))
         ملاحظة: يجب أن تكون أسماء الصفوف متناسقة مع أسماء ملفات bam المستخدمة من قبل featureCounts لحساب القروءات. يتكون جدول العينة من تفاصيل كل عينة والتي تشمل: نوع المكتبة والحالة. هذا مطلوب لتحديد التباين أو مجموعة الاختبار للكشف عن الاستخدام التفاضلي.
      2. قم بإعداد ملف معلومات exon. معلومات Exon في شكل كائنات GRanges (النطاقات الجينومية) (https://bioconductor.org/packages/release/bioc/html/GenomicRanges.html) مطلوبة كمدخل لإنشاء كائن DEXSeq في الخطوة التالية. قم بمطابقة معرفات الجينات مع أعداد القراءة لإنشاء كائن exoninfo.
         exoninfo<-anno[anno$GeneID%in% rownames(countData$counts),]
exoninfo<-GRanges(seqnames=anno$Chr,
ranges=IRanges(start=anno$Start, end=anno$End, width=anno$Width),strand=Rle(anno$Strand))
mcols(exoninfo)$TranscriptIDs<-anno$TranscriptIDs
mcols(exoninfo)$Ticker<-anno$Ticker
mcols(exoninfo)$ExonID<-anno$ExonID
mcols(exoninfo)$n<-anno$n
mcols(exoninfo)$GeneID<-anno$GeneID
transcripts_l= strsplit(exoninfo$TranscriptIDs, "\\,")
save(countData, sampleTable, exoninfo, transcripts_l, file="AS_countdata.RData")
      3. قم بإنشاء كائن DEXSeq باستخدام وظيفة DEXSeqDataSet. يجمع كائن DEXSeq عدد مرات القراءة ومعلومات معالم exon ومعلومات العينة. استخدم أعداد القراءة التي تم إنشاؤها في الخطوة 3 ومعلومات exon التي تم الحصول عليها من الخطوة السابقة لإنشاء كائن DEXSeq من مصفوفة العد. تأخذ الوسيطة sampleData إدخال إطار بيانات يحدد العينات (وسماتها: نوع المكتبة وشرطها) ، ويستخدم "التصميم" sampleData لإنشاء مصفوفة تصميم للاختبار التفاضلي باستخدام تدوين صيغة النموذج. لاحظ أن مصطلح التفاعل المهم ، condition: exon ، يشير إلى أن جزء القراءات على جين يسقط على إكسون معين يعتمد على الحالة التجريبية ، أي أن هناك AS. راجع وثائق DEXSeq للحصول على وصف كامل لإعداد صيغة النموذج للتصميمات التجريبية الأكثر تعقيدا. للحصول على معلومات الميزة ، يلزم وجود معرفات exon والجينات المقابلة والنسخ.
         ​dxd<-DEXSeqDataSet(countData$counts,sampleData=sampleTable, design=~sample+exon+condition:exon,featureID=exoninfo$ExonID,groupID=exoninfo$GeneID,featureRanges=exoninfo, transcripts=transcripts_l)
      4. التطبيع وتقدير التشتت. بعد ذلك ، قم بإجراء التطبيع بين العينات وتقدير تباين البيانات ، بسبب كل من ضوضاء حساب بواسون من الطبيعة المنفصلة ل RNA-seq والتباين البيولوجي ، باستخدام الأوامر التالية.
         ​dxd %<>% estimateSizeFactors %>% estimateDispersions %T>% plotDispEsts
      5. اختبار الاستخدام التفاضلي. بعد تقدير التباين ، اختبر استخدام الإكسون التفاضلي لكل جين وقم بإنشاء النتائج.
         dxd%<>%testForDEU%>%estimateExonFoldChanges(fitExpToVar=
"condition")#Estimate fold changes
dxr=DEXSeqResults(dxd)
      6. تصور أحداث الربط للجينات المحددة باستخدام الأمر التالي.
         plotDEXSeq(dxr,"Wnk1", displayTranscripts=TRUE, splicing=TRUE,legend
=TRUE,cex.axis=1.2,cex=1.3,lwd=2)
         افحص ملف R Notebook "AS_analysis_RNASeq.Rmd" لإنشاء قطع إضافية للجينات ذات الأهمية ولإنشاء قطع أراضي بركانية عند عتبات مختلفة.
   2. استخدام diffSplice من Limma لتحديد الربط التفاضلي. اتبع ملف R Notebook "AS_analysis_RNASeq.Rmd". تأكد من اتباع الخطوات 2.1-2.3 لإعداد ملفات الإدخال قبل المضي قدما.
      1. تحميل المكتبات
         library(limma)
library(edgeR)
      2. تصفية غير محددة. استخراج مصفوفة عدد القراءة التي تم الحصول عليها في 2.3. قم بإنشاء قائمة بالميزات باستخدام وظيفة "DGEList" من حزمة edgeR ، حيث تمثل الصفوف الجينات وتمثل الأعمدة العينات.
         mycounts=countData$counts
#Change the rownames of the countdata to exon Ids instead of genes for unique rownames.
rownames(mycounts) = exoninfo$ExonID
dge<-DGEList(counts=mycounts)
#Filtering
isexpr<- rowSums(cpm(dge)>1) >=3
dge<-dge[isexpr,,keep.lib.sizes=FALSE]
#Extract the exon annotations for only transcripts meeting non-specific filter
exoninfo=anno%>% filter(ExonID%in% rownames(dge$counts))
#Convert the exoninfo into GRanges object
exoninfo1<-GRanges(seqnames=exoninfo$Chr,
ranges=IRanges(start=exoninfo$Start, end=exoninfo$End, width=exoninfo$Width),strand=Rle(exoninfo$Strand))
mcols(exoninfo1)$TranscriptIDs<-exoninfo$TranscriptIDs
mcols(exoninfo1)$Ticker<-exoninfo$Ticker
mcols(exoninfo1)$ExonID<-exoninfo$ExonID
mcols(exoninfo1)$n<-exoninfo$n
mcols(exoninfo1)$GeneID<-exoninfo$GeneID
transcripts_l= strsplit(exoninfo1$TranscriptIDs, "\\,")
         ملاحظة: كخطوة تصفية غير محددة، تتم تصفية الأعداد بواسطة cpm < 1 في x من أصل n عينات، حيث x هو الحد الأدنى لعدد النسخ المتماثلة في أي شرط. n = 6 و x = 3 لبيانات المثال هذه.
      3. تطبيع الأعداد عبر العينات ، باستخدام وظيفة "calcNormFactors" من حزمة "edgeR" باستخدام المتوسط المقتطع لقيم M (طريقة تسوية TMM)²⁴ سيحسب عوامل القياس لضبط أحجام المكتبة.
         ​dge<-calcNormFactors(dge)
      4. استخدم sampleTable كما تم إنشاؤه في الخطوة 2.4.1.1 وقم بإنشاء مصفوفة التصميم. مصفوفة التصميم تميز التصميم. راجع دليل مستخدم Limma (https://www.bioconductor.org/packages/devel/bioc/vignettes/limma/inst/doc/usersguide.pdf) الفصلين 8 و 9 للحصول على تفاصيل حول مصفوفات التصميم للحصول على تصميمات تجريبية أكثر تقدما.
         Treat<- factor(sampleTable$condition)
design<- model.matrix(~0+Treat)
colnames(design) <- levels(Treat)
      5. تناسب نموذج خطي لكل إكسون. قم بتشغيل وظيفة "voom" لحزمة "limma" لمعالجة بيانات RNA-seq لتقدير التباين وإنشاء أوزان دقيقة لتصحيح ضوضاء عدد Poisson ، وتحويل الأعداد على مستوى exon إلى log2-counts لكل مليون (logCPM). ثم قم بتشغيل النمذجة الخطية باستخدام وظيفة "lmfit" لملاءمة النماذج الخطية مع بيانات التعبير لكل إكسون. احسب إحصائيات بايز التجريبية للنموذج المجهز باستخدام وظيفة "eBayes" للكشف عن تعبير إكسون التفاضلي. بعد ذلك ، حدد مصفوفة تباين للمقارنات التجريبية ذات الأهمية. استخدم "contrasts.fit" للحصول على المعاملات والأخطاء القياسية لكل زوج من المقارنة.
         v<-voom(dge,design,plot=FALSE)
fit<-lmFit(v,design)
fit<-eBayes(fit)
colnames(fit)
cont.matrix<-makeContrasts(
Mbnl1_KO_WT=Mbnl1_KO-WT,
levels=design)
fit2<-contrasts.fit(fit,cont.matrix)
      6. تحليل الربط التفاضلي. قم بتشغيل "diffSplice" على النموذج المجهز لاختبار الاختلافات في استخدام exon للجينات بين النوع البري والضربة القاضية واستكشاف النتائج الأعلى مرتبة باستخدام وظيفة "topSplice": test = "t" يعطي ترتيبا ل AS exons ، test = "simes" يعطي ترتيبا للجينات.
         ex<-diffSplice(fit2,geneid=exoninfo$GeneID,exonid=exoninfo$ExonID)
ts<-topSplice(ex,n=Inf,FDR=0.1, test="t", sort.by="logFC")
​tg<-topSplice(ex,n=Inf,FDR=0.1, test="simes")
      7. التصور. ارسم النتائج باستخدام وظيفة "plotSplice" ، مع إعطاء الجين محل الاهتمام في حجة geneid. احفظ أفضل النتائج التي تم فرزها حسب السجل قم بطي التغيير إلى كائن وقم بإنشاء مخطط بركان لعرض exons.
         plotSplice(ex,geneid="Wnk1", FDR=0.1)
#Volcano plot
EnhancedVolcano(ts,lab=ts$ExonID,selectLab= head((ts$ExonID),2000), xlab= bquote(~Log[2]~'fold change'), x='logFC', y='P.Value', title='Volcano Plot', subtitle='Mbnl1_KO vs WT (Limma_diffSplice)', FCcutoff=2, labSize=4,legendPosition="right", caption= bquote(~Log[2]~"Fold change cutoff, 2; FDR 10%"))
   3. باستخدام rMATS
      1. تأكد من تثبيت أحدث إصدار من rMATS v4.1.1 (المعروف أيضا باسم rMATS turbo نظرا لانخفاض وقت المعالجة ومتطلبات الذاكرة الأقل) إما باستخدام conda أو github (https://github.com/Xinglab/rmats-turbo/releases/download/v4.1.1/rmats_turbo_v4_1_1.tar.gz) في دليل العمل. اتبع القسم 4.3 في "AS_analysis_RNASeq.Rmd".
      2. انتقل إلى المجلد الذي يحتوي على ملفات bam التي تم الحصول عليها بعد التعيين وقم بإعداد الملفات النصية ، كما هو مطلوب بواسطة rMATS ، للشرطين عن طريق نسخ اسم ملفات bam (جنبا إلى جنب مع المسار) مفصولة ب ','. يجب تشغيل الأوامر التالية في سطر أوامر Linux:
         mkdir rMATS_analysis
cd bams/
ls -pd "$PWD"/*| grep "WT"| tr '\n'','> Wt.txt
ls -pd "$PWD"/*| grep "Mb"| tr '\n'','> KO.txt
mv *.txt ../rMATS_analysis
      3. قم بتشغيل rmats.py باستخدام ملفي الإدخال اللذين تم إنشاؤهما في الخطوة السابقة ، جنبا إلى جنب مع ملف GTF الذي تم الحصول عليه في 2.1.1.3. سيؤدي هذا إلى إنشاء مجلد إخراج "rmats_out" يحتوي على ملفات نصية تصف الإحصائيات (قيم p ومستويات التضمين) لكل حدث ربط على حدة.
         python rmats-turbo/rmats.py --b1 KO.txt --b2 Wt.txt --gtf annotation.gtf -t paired --readLength 50 --nthread 8 --od rmats_out/ --tmp rmats_tmp --task pos
         ملاحظة: التعليق التوضيحي المرجعي في شكل ملف GTF مطلوب أيضا. تحقق من المعلمات إذا كانت البيانات أحادية الطرف ، وقم بتغيير الخيار -t وفقا لذلك.
      4. استكشاف نتائج rMATS. استخدم حزمة الموصلات الحيوية "^maser'25 لاستكشاف نتائج rMATS. قم بتحميل الملفات النصية ل Junction و exon counts (JCEC) في كائن "maser" وقم بتصفية النتيجة بناء على التغطية من خلال تضمين خمسة قراءات متوسطة على الأقل لكل حدث ربط.
         library(maser)
mbnl1<-maser("/rmats_out/", c("WT","Mbnl1_KO"), ftype="JCEC")
#Filtering out events by coverage
mbnl1_filt<-filterByCoverage(mbnl1,avg_reads=5)
      5. تصور نتائج rMATS. حدد أحداث الربط المهمة بمعدل الاكتشاف الخاطئ (FDR) 10٪ والحد الأدنى للتغيير بنسبة 10٪ في النسبة المئوية للتقسيم (deltaPSI) باستخدام وظيفة "topEvents" من حزمة "maser". بعد ذلك ، تحقق من أحداث الجينات للجينات الفردية ذات الأهمية (عينة الجين Wnk1) وارسم قيم PSI لكل حدث ربط لهذا الجين. قم بإنشاء مخطط بركان عن طريق تحديد نوع الحدث.
         #Top splicing events at 10% FDR
mbnl1_top<-topEvents(mbnl1_filt,fdr=0.1, deltaPSI=0.1)
mbnl1_top
#Check the gene events for a particular gene
mbnl1_wnk1<-geneEvents(mbnl1_filt,geneS="Wnk1", fdr=0.1, deltaPSI=0.1)
maser::display(mbnl1_wnk1,"SE")
plotGenePSI(mbnl1_wnk1,type="SE", show_replicates
=TRUE)
​volcano(mbnl1_filt,fdr=0.1, deltaPSI=0.1,type="SE")
+xlab("deltaPSI")+ylab("Log10 Adj. Pvalue")+ggtitle("Volcano Plot of exon skipping events")
      6. قم بإنشاء مخططات Sashimi لنتائج أحداث الربط التي تم الحصول عليها باستخدام rMATS في شكل ملفات نصية باستخدام حزمة "rmats2shahimiplot". قم بتشغيل البرنامج النصي python في سطر أوامر Linux.
         python ./src/rmats2sashimiplot/rmats2sashimiplot.py --b1 ../bams/WT_Thymus_1.bam,../bams/WT_Thymus_2.bam,../bams/WT_Thymus_3.bam --b2 ../bams/Mbnl1KO_Thymus_1.bam,../bams/Mbnl1KO_Thymus_2.bam,../bams/Mbnl1KO_Thymus_3.bam -t SE -e ../rMATS_analysis/rmats_out/SE.MATS.JC.txt --l1 WT --l2 Mbnl1_KO --exon_s 1 --intron_s 5 -o ../rMATS_analysis/rmats2shasmi_output
         ملاحظة: يمكن أن تستغرق هذه العملية وقتا طويلا لأنها ستنشئ مؤامرة الساشيمي لجميع النتائج في ملف الأحداث. اختر أفضل النتائج (أسماء الجينات والإكسونات) كما هو معروض بواسطة وظيفة topEvents من "maser" وتصور مؤامرة الساشيمي المقابلة.

3. تحليل polyadenylation البديل (APA) باستخدام تسلسل نهاية 3 '

1. تنزيل البيانات ومعالجتها مسبقا
   ملاحظة: راجع ملف R Notebook التكميلي "APA_analysis_3PSeq_notebook. Rmd" للأوامر الكاملة لتنزيل البيانات وخطوات المعالجة المسبقة ، أو قم بتشغيل ملف bash التكميلي "APA_data_downloading_preprocessing.sh" في سطر أوامر Linux.
   1. قم بتنزيل البيانات من SRA باستخدام معرفات الانضمام (1553129 إلى 1553136).
   2. تقليم المحولات وتكملة عكسية للحصول على تسلسل حبلا الإحساس.
      ملاحظة: هذه الخطوة خاصة بمقايسة PolyA-seq المستخدمة.
   3. تقرأ الخريطة لتجميع جينوم الماوس باستخدام مصفف ربطة العنق²⁶.
2. إعداد التعليقات التوضيحية لمواقع pA.
   ملاحظة: تتم معالجة ملف التعليقات التوضيحية لموقع pA أولا باستخدام ملف R Notebook التكميلي "APA_analysis_3PSeq_notebook. Rmd" (2.1 - 2.6) ، ثم استخدام ملف bash "APA_annotation_preparation.sh".
   1. قم بتنزيل التعليق التوضيحي لمواقع pA من قاعدة بيانات PolyASite 2.0⁶.
   2. حدد التعليقات التوضيحية لموقع pA للاحتفاظ بمواقع pA للمنطقة غير المترجمة (UTR) ، والتي يتم شرحها ك Terminal Exon (TE) أو 1000 nt في اتجاه مجرى النهر من exon الطرفي المشروح (DS) لتحليل المصب.
   3. الحصول على قمم موقع pA. مرساة في كل موقع انشقاق pA ، وتصور متوسط تغطية القراءة باستخدام أدوات السرير والأدوات العميقة^27,28. أظهرت النتائج أن قمم القراءات المعينة كانت مشتتة بشكل أساسي داخل ~ 60 نقطة أساس في المنبع لمواقع الانقسام (الشكل 5 والشكل التكميلي 5). لذلك ، تم توسيع إحداثيات مواقع pA من ملف التعليقات التوضيحية إلى 60 نقطة أساس في المنبع لمواقع الانقسام الخاصة بها. اعتمادا على بروتوكول التسلسل النهائي المحدد 3 'المستخدم ، ستحتاج هذه الخطوة إلى التحسين للمقايسات بخلاف PolyA-seq.
3. يقرأ العد
   1. قم بإعداد ملف التعليقات التوضيحية لمواقع pA.
      anno<- read.table(file= "flanking60added.pA_annotation.bed",
stringsAsFactors=FALSE, check.names=FALSE, header=FALSE, sep="")
colnames(anno) <- c("chrom", "chromStart", "chromEnd", "name", "score", "strand", "rep", "annotation", "gene_name", "gene_id")
anno<- dplyr::select(anno,name,chrom, chromStart,chromEnd, strand,gene_id,gene_name,rep)
colnames(anno) <- c("GeneID", "Chr", "Start", "End", "Strand", "Ensembl", "Symbol", "repID")
   2. قم بتطبيق "featureCounts" للحصول على الأعداد الأولية. احفظ جدول العد كملف RData "APA_countData.Rdata" لتحليل APA باستخدام أدوات مختلفة.
      countData<- dir("bamfiles", pattern="sorted.bam$", full.names=TRUE) %>%
# Read all bam files as input for featureCounts
featureCounts(annot.ext=anno, isGTFAnnotationFile= FALSE,minMQS=0,useMetaFeatures= TRUE,allowMultiOverlap=TRUE, largestOverlap= TRUE,strandSpecific=1, countMultiMappingReads =TRUE,primaryOnly= TRUE,isPairedEnd= FALSE,nthreads=12)%T>%
save(file="APA_countData.Rdata")
      ملاحظة: كن واعيا لتغيير أي من المعلمات المدرجة في الدالة 'featureCounts'. قم بتعديل المعلمة "strandSpecific" للتأكد من توافقها مع اتجاه التسلسل لمقايسة التسلسل النهائي 3 'المستخدمة (تجريبيا ، فإن تصور البيانات في متصفح الجينوم على الجينات على خيوط زائد وناقص سيوضح ذلك).
   3. تطبيق تصفية غير محددة ل countData. يمكن أن تؤدي التصفية إلى تحسين المتانة الإحصائية بشكل كبير في اختبارات استخدام موقع pA التفاضلية. أولا ، أزلنا تلك الجينات مع موقع pA واحد فقط ، والذي لا يمكن تحديد استخدام موقع pA بشكل تفاضلي. ثانيا ، نطبق ترشيحا غير محدد بناء على التغطية: يتم تصفية الأعداد بواسطة cpm أقل من 1 في x من عينات n ، حيث x هو الحد الأدنى لعدد النسخ المتماثلة في أي حالة. N = 8 و x = 2 لبيانات المثال هذه.
      load(file= "APA_countData.Rdata")# Skip this step if already loaded
# Non-specific filtering: Remove the pA sites not differentially expressed in the samples
countData<-countData$counts%>%as.data.frame%>% .[rowSums(edgeR::cpm(.)>1) >=2, ]
anno%<>% .[.$GeneID%in% rownames(countData), ]
# Remove genes with only 1 site and NA in geneIDs
dnsites<-anno%>%group_by(Symbol)%>%summarise(nsites=n())%>% filter(nsites>1&!is.na(Symbol))
anno<-anno%>% filter(Symbol%in%dnsites$Symbol)
countData<-countData[rownames(countData) %in%anno$GeneID, ]
4. تحليل استخدام موقع البولي أدينيل التفاضلي باستخدام خطوط أنابيب DEXSeq و diffSplice.
   1. باستخدام حزمة DEXSeq
      ملاحظة: نظرا لأنه لا يمكن تحديد مصفوفة تباين لخط أنابيب DEXSeq ، يجب إجراء تحليل APA التفاضلي لكل شرطين تجريبيين بشكل منفصل. يتم إجراء تحليل APA التفاضلي للحالة WT والحالة DKO كمثال لشرح الإجراء. الرجوع إلى الملف التكميلي "APA_analysis_3PSeq_notebook. رمد" لسير العمل خطوة بخطوة لهذا القسم وتحليل APA التفاضلي للتناقضات الأخرى.
      1. قم بتحميل المكتبة وإنشاء جدول عينة لتحديد التصميم التجريبي.
         c("DEXSeq", "GenomicRanges") %>% lapply(library, character.only=TRUE) %>%invisible
sampleTable1<- data.frame(row.names= c("WT_1","WT_2","DKO_1","DKO_2"),
condition= c(rep("WT", 2), rep("DKO", 2)),
​libType= rep("single-end", 4))
      2. قم بإعداد ملف معلومات مواقع pA باستخدام حزمة الموصلات الحيوية GRanges.
         # Prepare the GRanges object for DEXSeqDataSet object construction
PASinfo <- GRanges(seqnames = anno$Chr,
ranges = IRanges(start = anno$Start, end = anno$End),strand = Rle(anno$Strand))
mcols(PASinfo)$PASID<-anno$repID
mcols(PASinfo)$GeneEns<-anno$Ensembl
mcols(PASinfo)$GeneID<-anno$Symbol
# Prepare the new feature IDs, replace the strand information with letters to match the current pA site clusterID
new.featureID <- anno$Strand %>% as.character %>% replace(. %in% "+", "F") %>% replace(. %in% "-", "R") %>% paste0(as.character(anno$repID), .)
      3. استخدم عدد مرات القراءة التي تم إنشاؤها في الخطوة 3.3 ومعلومات موقع pA التي تم الحصول عليها من الخطوة السابقة لإنشاء كائن DEXSeq.
         # Select the read counts of the condition WT and DKO
countData1<- dplyr::select(countData, SRR1553129.sorted.bam, SRR1553130.sorted.bam, SRR1553131.sorted.bam, SRR1553132.sorted.bam)
# Rename the columns of countData using sample names in sampleTable
colnames(countData1) <- rownames(sampleTable1)
dxd1<-DEXSeqDataSet(countData=countData1,
sampleData=sampleTable1,
design=~sample+exon+condition:exon,
featureID=new.featureID,
groupID=anno$Symbol,
featureRanges=PASinfo)
      4. حدد زوج التباين من خلال تحديد مستويات الظروف في كائن DEXSeq.
         dxd1$condition<- factor(dxd1$condition, levels= c("WT", "DKO"))
# The contrast pair will be "DKO - WT"
      5. التطبيع وتقدير التشتت. على غرار بيانات RNA-seq ، بالنسبة لبيانات التسلسل النهائي 3 بوصات ، قم بإجراء التطبيع بين العينات (متوسط النسب حسب العمود لكل عينة) باستخدام وظيفة "estimateSizeFactors" ، وتقدير تباين البيانات باستخدام وظيفة "estimateDispersions" ، ثم تصور نتيجة تقدير التشتت باستخدام وظيفة "plotDispEsts".
         ​dxd1 %<>% estimateSizeFactors %>% estimateDispersions %T>% plotDispEsts
      6. اختبار استخدام موقع pA التفاضلي لكل جين باستخدام الوظيفة "testForDEU" ، ثم تقدير تغيير أضعاف استخدام موقع pA باستخدام الوظيفة "estimateExonFoldChanges". تحقق من النتائج باستخدام الوظيفة "DEXSeqResults" وقم بتعيين "FDR < 10٪" كمعيار لمواقع pA التفاضلية بشكل كبير.
         dxd1 %<>% testForDEU %>% estimateExonFoldChanges(fitExpToVar = "condition")
dxr1 <- DEXSeqResults(dxd1)
dxr1
mcols(dxr1)$description
table(dxr1$padj<0.1) # Check the number of differential pA sites (FDR < 0.1)
table(tapply(dxr1$padj<0.1, dxr1$groupID, any)) # Check the number of gene overlapped with differential pA site
      7. تصور نتائج استخدام موقع pA التفاضلي باستخدام مخططات APA التفاضلية الناتجة عن الوظيفة "plotDEXSeq" ومخطط البركان بواسطة الوظيفة "EnhancedVolcano".
         # Select the top 100 significant differential pA sites ranked by FDR
topdiff.PAS<- dxr1%>%as.data.frame%>%rownames_to_column%>%arrange(padj)%$%groupID[1:100]

# Apply plotDEXSeq for the visualization of differential polyA usage
plotDEXSeq(dxr1,"S100a7a", legend=TRUE, expression=FALSE,splicing=TRUE, cex.axis=1.2, cex=1.3,lwd=2)

# Apply perGeneQValue to check the top genes with differential polyA site usage
dxr1%<>% .[!is.na(.$padj), ]
dgene<- data.frame(perGeneQValue= perGeneQValue (dxr1)) %>%rownames_to_column("groupID")

dePAS_sig1<-dxr1%>% data.frame() %>%
dplyr::select(-matches("dispersion|stat|countData|genomicData"))%>%
inner_join(dgene)%>%arrange(perGeneQValue)%>%distinct()%>%
filter(padj<0.1)

# Apply EnhancedVolcano package to visualise differential polyA site usage
"EnhancedVolcano"%>% lapply(library, character.only=TRUE) %>%invisible
EnhancedVolcano(dePAS_sig1, lab=dePAS_sig1$groupID, x='log2fold_DKO_WT',
y='pvalue',title='Volcano Plot',subtitle='DKO vs WT',
FCcutoff=1,labSize=4, legendPosition="right",
caption= bquote(~Log[2]~"Fold change cutoff, 1; FDR 10%"))
   2. باستخدام حزمة diffSplice. راجع ملف R Notebook التكميلي "APA_analysis_3PSeq_notebook. Rmd" لسير العمل خطوة بخطوة لهذا القسم.
      1. تحديد تباينات الاهتمام لتحليل استخدام pA التفاضلي.
         ملاحظة: يجب تنفيذ هذه الخطوة بعد إنشاء ومعالجة كائن DGEList ، والذي تم تضمينه في ملف R Notebook "APA_analysis_3PSeq_notebook. رمد".
         contrast.matrix<-makeContrasts(DKO_vs_WT=DKO-
WT,Ctrl_vs_DKO=Ctrl-DKO,
KD_vs_Ctrl=KD-Ctrl,KD_vs_DKO=KD-DKO,levels=design)
fit2<-fit%>%contrasts.fit(contrast.matrix)%>%eBayes
summary(decideTests(fit2))
ex<-diffSplice(fit2,geneid=anno$Symbol,exonid=new.featureID)
topSplice(ex) #Check the top significant results with topSplice
      2. تصور نتيجة تناقضات الاهتمام (هنا "DKO - WT") باستخدام مخططات APA التفاضلية بواسطة الوظيفة "plotSplice" ومؤامرات البركان مع الوظيفة "EnhancedVolcano". راجع ملف R Notebook "APA_analysis_3PSeq_notebook. Rmd" 4.2.7 - 4.2.9 لتصور أزواج التباين الأخرى.
         sig1<-topSplice(ex,n=Inf,FDR=0.1,coef=1, test="t", sort.by="logFC")
sig1.genes<-topSplice(ex,n=Inf,FDR=0.1,coef=1, test="simes")
plotSplice(ex, coef=1,geneid="S100a7a", FDR = 0.1)
plotSplice(ex,coef=1,geneid="Tpm1", FDR = 0.1)
plotSplice(ex,coef=1,geneid="Smc6", FDR = 0.1)
EnhancedVolcano(sig1, lab=sig1$GeneID,xlab= bquote(~Log[2]~'fold change'),
x='logFC', y='P.Value', title='Volcano Plot', subtitle='DKO vs WT',
FCcutoff=1, labSize=6, legendPosition="right")

النتائج

بعد تشغيل سير العمل خطوة بخطوة أعلاه ، تكون مخرجات تحليل AS و APA والنتائج التمثيلية في شكل جداول ومخططات بيانات ، يتم إنشاؤها على النحو التالي.

مثل:
الناتج الرئيسي لتحليل AS (الجدول التكميلي 1 ل diffSplice; الجدول 2 ل DEXSeq) عبارة عن قائمة بالإكسونات التي تظهر الاستخدام التفاضلي عبر الظروف ، وقائمة بالجينات التي تظهر نشاط الربط الكلي الكبير لواحد أو أكثر من الإكسونات المكونة لها ، مرتبة حسب الدلالة الإحصائية. الجدول التكميلي 1

المناقشة

في هذه الدراسة ، قمنا بتقييم النهج القائمة على exon والقائمة على الأحداث للكشف عن AS و APA في بيانات تسلسل RNA-Seq السائبة و 3 'التسلسل النهائي. تنتج مناهج AS القائمة على exon كلا من قائمة الإكسونات المعبر عنها بشكل تفاضلي وترتيب على مستوى الجينات مرتبة حسب الأهمية الإحصائية لنشاط الربط التفاضلي على مستوى الجينات الإجمالي (الجداول 1-2 ، 4-5). بالنسبة لحزمة diffSplice ، يتم تحديد الاستخدام التفاضلي عن طريق تركيب نماذج خطية مرجحة على مستوى exon لتقدير تغير طية اللوغاريتم التفاضلية ل exon مقابل متوسط تغير ط...

المواد

قائمة المواد المستخدمة في هذه المقالة

الاسم	الشركة	رقم فهرسي	التعليقات
غير ذات صلة بالدراسة