זיהוי שחבור חלופי ופוליאדנילציה בנתוני RNA-seq

Summary

שחבור חלופי (AS) ופוליאדנילציה חלופית (APA) מרחיבים את מגוון האיזופורמים של שעתוק ותוצרים. כאן אנו מתארים פרוטוקולים ביואינפורמטיים לניתוח מבחני RNA-seq בתפזורת ו-3' של ריצוף קצה כדי לזהות ולהמחיש AS ו-APA המשתנים בין תנאי ניסוי.

Abstract

בנוסף לניתוח הטיפוסי של RNA-Seq למדידת ביטוי גנים דיפרנציאלי (DGE) בתנאים ניסיוניים/ביולוגיים, ניתן להשתמש בנתוני RNA-seq גם כדי לחקור מנגנוני ויסות מורכבים אחרים ברמת האקסון. שחבור חלופי ופוליאדנילציה ממלאים תפקיד מכריע במגוון התפקודי של גן על ידי יצירת איזופורמים שונים כדי לווסת את ביטוי הגנים ברמה שלאחר השעתוק, והגבלת הניתוחים לרמת הגן כולה עלולה לפספס את שכבת הוויסות החשובה הזו. כאן, אנו מדגימים ניתוחים מפורטים צעד אחר צעד לזיהוי והדמיה של שימוש באקסון דיפרנציאלי ובאתר פוליאדנילציה בתנאים שונים, תוך שימוש ב- Bioconductor ובחבילות ופונקציות אחרות, כולל DEXSeq, diffSplice מחבילת לימה ו- rMATS.

Introduction

RNA-seq היה בשימוש נרחב לאורך השנים בדרך כלל להערכת ביטוי גנים דיפרנציאליים וגילוי גנים¹. בנוסף, ניתן להשתמש בו גם כדי להעריך שימוש משתנה ברמת אקסון עקב גנים המבטאים איזופורמים שונים, ובכך לתרום להבנה טובה יותר של ויסות גנים ברמה שלאחר השעתוק. רוב הגנים האאוקריוטים מייצרים איזופורמים שונים על ידי שחבור חלופי (AS) כדי להגדיל את המגוון של ביטוי mRNA. ניתן לחלק את אירועי AS לדפוסים שונים: דילוג על אקסונים שלמים (SE) שבהם אקסון ("קלטת") מוסר לחלוטין מהתמליל יחד עם האינטרונים האגפים שלו; בחירת אתר שחבור חלופי (תורם) 5' (A5SS) וחלופה 3' (מקבל) בחירת אתר שחבור (A3SS) כאשר שני אתרי שחבור או יותר נמצאים משני קצותיו של אקסון; שמירה של אינטרונים (RI) כאשר אינטרון נשמר בתוך תעתיק mRNA בוגר והדרה הדדית של שימוש באקסון (MXE) כאשר רק אחד משני האקסונים הזמינים יכול להישמר בכל פעם ^2,3. פוליאדנילציה חלופית (APA) ממלאת גם תפקיד חשוב בוויסות ביטוי גנים באמצעות אתרי פולי (A) חלופיים ליצירת איזופורמים מרובים של mRNA מתעתיק יחיד⁴. רוב אתרי הפוליאדנילציה (pAs) ממוקמים באזור 3' לא מתורגם (3' UTRs), ומייצרים איזופורמים של mRNA עם אורכי UTR מגוונים של 3'. מכיוון שה-UTR 3' הוא הרכזת המרכזית לזיהוי אלמנטים רגולטוריים, אורכי UTR שונים של 3' יכולים להשפיע על לוקליזציה, יציבות ותרגום^{mRNA 5}. ישנם סוגים של מבחני רצף קצה 3 'הממוטבים לזיהוי APA השונים בפרטי הפרוטוקול⁶. הצינור המתואר כאן מיועד ל- PolyA-seq, אך ניתן להתאים אותו לפרוטוקולים אחרים כמתואר.

במחקר זה אנו מציגים צינור של שיטות ניתוח אקסון דיפרנציאליות^7,8 (איור 1), שניתן לחלק לשתי קטגוריות רחבות: מבוססות אקסון (DEXSeq^9, diffSplice 10) ומבוססות אירועים (ניתוח רב-משתני משוכפל של שחבור תעתיק (rMATS)¹¹). השיטות המבוססות על אקסון משוות את שינוי הקיפול על פני תנאים של אקסונים בודדים, כנגד מידה של שינוי קיפול גנים כולל כדי לקרוא לשימוש באקסון המתבטא באופן דיפרנציאלי, ומתוך כך מחשבים מדידה ברמת הגן של פעילות AS. שיטות מבוססות אירועים משתמשות בקריאות צומת exon-intron-spanning כדי לזהות ולסווג אירועי שחבור ספציפיים כגון דילוג אקסון או שמירה של אינטרונים, ולהבחין בין סוגי AS אלה בפלט³. לפיכך, שיטות אלה מספקות תצוגות משלימות לניתוח מלא של AS^12,13. בחרנו ב- DEXSeq (בהתבסס על חבילת DESeq2¹⁴ DGE) ו- diffSplice (בהתבסס על חבילת Limma¹⁰ DGE) למחקר מכיוון שהם בין החבילות הנפוצות ביותר לניתוח שחבור דיפרנציאלי. rMATS נבחרה כשיטה פופולרית לניתוח מבוסס אירועים. שיטה פופולרית נוספת המבוססת על אירועים היא MISO (תערובת של איזופורמים)¹. עבור APA אנו מתאימים את הגישה מבוססת exon.

איור 1. צינור ניתוח. תרשים זרימה של השלבים המשמשים בניתוח. השלבים כוללים: קבלת הנתונים, ביצוע בדיקות איכות ויישור קריאה ולאחר מכן ספירת קריאות באמצעות ביאורים עבור exons, introns ואתרי pA ידועים, סינון כדי להסיר ספירות נמוכות ונורמליזציה. נתוני PolyA-seq נותחו עבור אתרי pA חלופיים באמצעות שיטות diffSplice/DEXSeq, RNA-Seq בתפזורת נותח עבור שחבור חלופי ברמת האקסון בשיטות diffSplice/DEXseq, ואירועי AS נותחו עם rMATS. אנא לחץ כאן כדי להציג גרסה גדולה יותר של נתון זה.

נתוני ה-RNA-seq ששימשו בסקר זה נרכשו מביטוי גנים אומניבוס (GEO) (GSE138691)¹⁵. השתמשנו בנתוני RNA-seq של עכברים ממחקר זה עם שתי קבוצות מצב: נוקאאוט מסוג פראי (WT) ונוקאאוט דמוי שריר מסוג 1 (Mbnl1 KO) עם שלושה שכפולים כל אחד. כדי להדגים ניתוח שימוש באתר פוליאדנילציה דיפרנציאלית, השגנו נתוני פוליA-seq של עוברי עכברים (MEFs) (GEO Accession GSE60487)¹⁶. לנתונים יש ארבע קבוצות מצבים: סוג פראי (WT), דמוי שרירים מסוג 1/סוג 2 נוקאאוט כפול (Mbnl1/2 DKO), Mbnl 1/2 DKO עם נוקאאוט Mbnl3 (KD) ו- Mbnl1/2 DKO עם בקרת Mbnl3 (Ctrl). כל קבוצת תנאים מורכבת משני עותקים משוכפלים.

	הצטרפות גיאוגרפית	מספר הפעלה של SRA	שם לדוגמה	תנאי	לשכפל	טישו	רצף	אורך קריאה
רנ"א-סק	GSM4116218	SRR10261601	Mbnl1KO_Thymus_1	נוקאאוט Mbnl1	נציג 1	התימוס	סוף מזווג	100 כ"ס
	GSM4116219	SRR10261602	Mbnl1KO_Thymus_2	נוקאאוט Mbnl1	חזרה 2	התימוס	סוף מזווג	100 כ"ס
	GSM4116220	SRR10261603	Mbnl1KO_Thymus_3	נוקאאוט Mbnl1	חזרה 3	התימוס	סוף מזווג	100 כ"ס
	GSM4116221	SRR10261604	WT_Thymus_1	סוג פראי	נציג 1	התימוס	סוף מזווג	100 כ"ס
	GSM4116222	SRR10261605	WT_Thymus_2	סוג פראי	חזרה 2	התימוס	סוף מזווג	100 כ"ס
	GSM4116223	SRR10261606	WT_Thymus_3	סוג פראי	חזרה 3	התימוס	סוף מזווג	100 כ"ס
3P-Seq	GSM1480973	SRR1553129	WT_1	סוג פראי (WT)	נציג 1	פיברובלסטים עובריים לעכבר (MEFs)	קצה יחיד	40 bp
	GSM1480974	SRR1553130	WT_2	סוג פראי (WT)	חזרה 2	פיברובלסטים עובריים לעכבר (MEFs)	קצה יחיד	40 bp
	GSM1480975	SRR1553131	DKO_1	Mbnl 1/2 נוקאאוט כפול (DKO)	נציג 1	פיברובלסטים עובריים לעכבר (MEFs)	קצה יחיד	40 bp
	GSM1480976	SRR1553132	DKO_2	Mbnl 1/2 נוקאאוט כפול (DKO)	חזרה 2	פיברובלסטים עובריים לעכבר (MEFs)	קצה יחיד	40 bp
	GSM1480977	SRR1553133	DKOsiRNA_1	Mbnl 1/2 נוקאאוט כפול עם Mbnl 3 siRNA (KD)	נציג 1	פיברובלסטים עובריים לעכבר (MEFs)	קצה יחיד	40 bp
	GSM1480978	SRR1553134	DKOsiRNA_2	Mbnl 1/2 נוקאאוט כפול עם Mbnl 3 siRNA (KD)	חזרה 2	פיברובלסטים עובריים לעכבר (MEFs)	קצה יחיד	36 bp
	GSM1480979	SRR1553135	DKONTsiRNA_1	Mbnl 1/2 נוקאאוט כפול עם siRNA ללא מיקוד (Ctrl)	נציג 1	פיברובלסטים עובריים לעכבר (MEFs)	קצה יחיד	40 bp
	GSM1480980	SRR1553136	DKONTsiRNA_2	Mbnl 1/2 נוקאאוט כפול עם siRNA ללא מיקוד (Ctrl)	חזרה 2	פיברובלסטים עובריים לעכבר (MEFs)	קצה יחיד	40 bp

טבלה 1. סיכום של מערכי נתונים של RNA-Seq ו- PolyA-seq המשמשים לניתוח.

Protocol

1. התקנת כלים וחבילות R המשמשים בניתוח

1. קונדה הוא מנהל חבילות פופולרי וגמיש המאפשר התקנה נוחה של חבילות עם התלות שלהן בכל הפלטפורמות. השתמש ב- '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 (קבצים עם סיומת ". Rmd"), קבצי קוד R (קבצים עם סיומת ". R"), או לינוקס Bash מעטפת סקריפטים (קבצים עם סיומת ".sh"). יש לפתוח קבצי מחברת R (Rmd) ב- RStudio באמצעות קובץ | פתח קובץ..., וגושי קוד בודדים (שעשויים להיות פקודות R או פקודות מעטפת Bash) ולאחר מכן הפעל באופן אינטראקטיבי על-ידי לחיצה על החץ הירוק בפינה השמאלית העליונה. ניתן להפעיל קבצי קוד R על ידי פתיחה ב- RStudio, או בשורת הפקודה של לינוקס על ידי הקדמה עם "Rscript", למשל Rscript example.R. סקריפטים של Shell מופעלים בשורת הפקודה של לינוקס על ידי הקדמת הסקריפט עם הפקודה "sh" למשל.sh example.sh.

2. ניתוח שחבור חלופי (AS) באמצעות RNA-seq

1. הורדה ועיבוד מראש של נתונים
   הערה: מקטעי הקוד המבוארים להלן זמינים בקובץ הקוד המשלים "AS_analysis_RNASeq.RMD", כדי לבצע את השלבים הבודדים באופן אינטראקטיבי, והם מסופקים גם כסקריפט bash להיות מופעל באצווה על שורת הפקודה לינוקס (SH downloading_data_preprocessing.sh).
   1. הורדת הנתונים הגולמיים.
      1. הורד את הנתונים הגולמיים מארכיון קריאת רצף (SRA) באמצעות הפקודה 'prefetch' מתוך ערכת הכלים של SRA (v2.10.8)¹⁷. תן את מזהי הגישה של SRA ברצף בפקודה הבאה כדי להוריד אותם במקביל באמצעות כלי השירות המקבילי של גנו¹⁸. כדי להוריד קבצי 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 באמצעות הדברים הבאים בשורת הפקודה של לינוקס.
         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'21 כדי להסיר את הריצוף למתאמי אגף, המשתנה בהתאם לגודל מקטע ה- 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 מכיל ערכי אקסון מרובים עבור איזופורמים שונים. קובץ זה משמש כדי "לכווץ" את מזהי התמלילים המרובים עבור כל אקסון. זהו צעד חשוב להגדרת פחי ספירת אקסון.
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 (הפניה) כקלט כדי ליצור מטריצה של ספירות המשויכות לכל תכונה עם שורות המייצגות אקסונים (תכונות) ועמודות המייצגות דוגמאות.
      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. לאחר מכן, בצע סינון לא ספציפי כדי להסיר אקסונים בעלי ביטוי נמוך ("לא ספציפי" מציין שמידע מצב הניסוי אינו משמש בסינון, כדי למנוע הטיות בחירה). המר את הנתונים מקנה מידה גולמי לספירה למיליון (cpm) באמצעות הפונקציה cpm מחבילת 'edgeR'²³ ושמור אקסונים עם ספירות הגדולות מסף ניתן להגדרה (עבור ערכת נתונים זו נעשה שימוש ב- cpm אחד) בשלוש דגימות לפחות. גם להסיר גנים עם אקסון אחד בלבד.
      # 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 מבוסס לימה. ראה קובץ משלים: "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. מידע אקסון בצורה של 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 מקבל קלט מסגרת נתונים המגדיר את הדוגמאות (ואת התכונות שלהן: סוג ספרייה ותנאי), 'design' משתמש ב- 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 "AS_analysis_RNASeq.Rmd" כדי ליצור חלקות נוספות עבור גנים מעניינים וליצור חלקות הר געש בסף שונה.
   2. שימוש ב-diffSplice מלימה לזיהוי שחבור דיפרנציאלי. עקוב אחר קובץ המחברת R "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. השתמש בטבלה לדוגמה כפי שנוצרה בשלב 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' של חבילת 'לימה' כדי לעבד נתוני RNA-seq כדי להעריך שונות וליצור משקלים מדויקים לתיקון רעש ספירת פואסון, ולהפוך את הספירות ברמת אקסון ל- 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' על המודל המותאם כדי לבחון את ההבדלים בשימוש באקסון של גנים בין סוג בר לנוקאאוט ולחקור את התוצאות המדורגות ביותר באמצעות פונקציית 'topSplice': test="t" נותן דירוג של אקסונים AS, 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', המעניקה גן בעל עניין בארגומנט הגנטי. שמור את התוצאות העליונות ממוינות לפי יומן שינוי קיפול לאובייקט וצור חלקת הר געש כדי להציג את האקסונים.
         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. צור עלילות סשימי עבור תוצאת אירועי השחבור המתקבלת עם rMATS בצורה של קבצי טקסט באמצעות חבילת 'rmats2shahimiplot'. הפעל את סקריפט הפייתון בשורת הפקודה של לינוקס.
         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. ניתוח פוליאדנילציה חלופית (APA) באמצעות ריצוף קצה 3'

1. הורדה ועיבוד מראש של נתונים
   הערה: עיין בקובץ מחברת R המשלים "APA_analysis_3PSeq_notebook. Rmd" עבור הפקודות המלאות להורדת נתונים ושלבי עיבוד מראש, או הפעל את קובץ bash המשלים "APA_data_downloading_preprocessing.sh" בשורת הפקודה של לינוקס.
   1. הורד נתונים מ-SRA עם מזהי ההצטרפות (1553129 עד 1553136).
   2. חתוך מתאמים והשלמה הפוכה כדי לקבל את רצף גדילי החישה.
      הערה: שלב זה הוא ספציפי לבדיקת PolyA-seq שבה נעשה שימוש.
   3. המפה קוראת להרכבת גנום העכבר באמצעות קשת²⁶.
2. הכנת ביאורים לאתרי pA.
   הערה: העיבוד של קובץ הביאור של אתר ה- pA מתבצע תחילה באמצעות קובץ מחברת R משלים "APA_analysis_3PSeq_notebook. Rmd" (2.1 - 2.6), ולאחר מכן באמצעות קובץ bash "APA_annotation_preparation.sh".
   1. הורד ביאור אתרי pA ממסד הנתונים PolyASite 2.0⁶.
   2. בחר ביאורים של אתרי pA כדי לשמור על אתרי pA של 3'-אזור לא מתורגם (UTR), המבוארים כ- Terminal Exon (TE) או 1000 nt במורד הזרם של אקסון מסוף מבואר (DS) לניתוח במורד הזרם.
   3. השג פסגות אתר pA. עוגן בכל אתר מחשוף pA, והצג באופן חזותי את כיסוי הקריאה הממוצע באמצעות מצעים^{ו-deeptools 27,28}. התוצאות הראו כי הפסגות של הקריאות הממופות פוזרו בעיקר בטווח של ~60 bp במעלה הזרם של אתרי הבקיעה (איור 5 ואיור משלים 5). לכן, הקואורדינטות של אתרי pA הורחבו מקובץ הביאור ל-60 bp במעלה הזרם של אתרי הבקיעה שלהם. בהתאם לפרוטוקול הריצוף הסופי הספציפי בגודל 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. RMD" עבור זרימת העבודה שלב אחר שלב של סעיף זה וניתוח 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 באמצעות חבילת Bioconductor 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 המשלים "APA_analysis_3PSeq_notebook. Rmd" עבור זרימת העבודה שלב אחר שלב של סעיף זה.
      1. הגדר ניגודים מעניינים לניתוח שימוש דיפרנציאלי ב- pA.
         הערה: שלב זה צריך להתבצע לאחר בנייה ועיבוד של אובייקט DGEList, אשר כלול בקובץ מחברת R "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 "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) היא רשימה של אקסונים המציגים שימוש דיפרנציאלי ...

Discussion

במחקר זה, הערכנו גישות מבוססות אקסון ומבוססות אירועים לזיהוי AS ו-APA בנתוני RNA-Seq בתפזורת ובנתוני ריצוף קצה של 3'. גישות ה-AS המבוססות על אקסון מייצרות הן רשימה של אקסונים המבוטאים באופן דיפרנציאלי והן דירוג ברמת הגן המסודר לפי המובהקות הסטטיסטית של פעילות השחבור הדיפרנציאלית הכוללת ברמת הגן (...

Materials

Name	Company	Catalog Number	Comments
Not relevent for computational study