Browse code

bugs in events reclassification

JFerrer-B authored on 02/06/2022 07:40:24
Showing 4 changed files

... ...
@@ -7641,13 +7641,29 @@ reclasify_intern <- function(SG,mievento,pp1,pp2,ppref){
7641 7641
   ppref_2 <- strsplit(unlist(strsplit(ppref,",")),"-")
7642 7642
   
7643 7643
   pp1_2 <- unlist(lapply(pp1_2, function(X){
7644
-    intersect(grep(X[1],as.vector(SG$Edges$From)),grep(X[2],as.vector(SG$Edges$To)))
7644
+    if(identical(X[1],X[2])){
7645
+      X <- paste0(X,c(".a",".b"))
7646
+    }else{
7647
+      X <- paste0(X,c(".b",".a"))
7648
+    }
7649
+    which(as.vector(SG$Edges$From)==X[1] & as.vector(SG$Edges$To)==X[2])
7650
+    
7645 7651
   }))
7646 7652
   pp2_2 <- unlist(lapply(pp2_2, function(X){
7647
-    intersect(grep(X[1],as.vector(SG$Edges$From)),grep(X[2],as.vector(SG$Edges$To)))
7653
+    if(identical(X[1],X[2])){
7654
+      X <- paste0(X,c(".a",".b"))
7655
+    }else{
7656
+      X <- paste0(X,c(".b",".a"))
7657
+    }
7658
+    which(as.vector(SG$Edges$From)==X[1] & as.vector(SG$Edges$To)==X[2])
7648 7659
   }))
7649 7660
   ppref_2 <- unlist(lapply(ppref_2, function(X){
7650
-    intersect(grep(X[1],as.vector(SG$Edges$From)),grep(X[2],as.vector(SG$Edges$To)))
7661
+    if(identical(X[1],X[2])){
7662
+      X <- paste0(X,c(".a",".b"))
7663
+    }else{
7664
+      X <- paste0(X,c(".b",".a"))
7665
+    }
7666
+    which(as.vector(SG$Edges$From)==X[1] & as.vector(SG$Edges$To)==X[2])
7651 7667
   }))
7652 7668
   
7653 7669
   Event <- list(P1=SG$Edges[pp1_2,],P2=SG$Edges[pp2_2,],Ref=SG$Edges[ppref_2,])
... ...
@@ -7712,41 +7728,32 @@ reclasify_intern <- function(SG,mievento,pp1,pp2,ppref){
7712 7728
     
7713 7729
     
7714 7730
     while(TRUE){
7715
-      torm <- which(
7716
-        rowSums(newAdj[-nrow(newAdj),])==0)
7731
+      torm <- which(rowSums(newAdj[-nrow(newAdj),])==0)
7717 7732
       if(length(torm)>0){
7718 7733
         torm_index <- c()
7719 7734
         for(ssx in 1:length(torm)){
7720 7735
           # ssx <- 2
7721
-          subexontoremove <- gsub(
7722
-            "\\..*","",names(torm)[ssx])
7723
-          torm_index <- c(torm_index,
7724
-                          grep(subexontoremove,
7725
-                               rownames(newAdj)))
7736
+          subexontoremove <- names(torm)[ssx]
7737
+          torm_index <- c(torm_index,which(rownames(newAdj)==subexontoremove))
7726 7738
         }
7727
-        newAdj <- newAdj[-torm_index,
7728
-                         -torm_index]     
7739
+        newAdj <- newAdj[-torm_index,-torm_index]     
7729 7740
       }else{
7730 7741
         break
7731 7742
       }
7732 7743
     }
7733 7744
     
7745
+    
7746
+    
7734 7747
     while(TRUE){
7735
-      torm <- which(
7736
-        colSums(newAdj[,-1])==0)
7748
+      torm <- which(colSums(newAdj[,-1])==0)
7737 7749
       if(length(torm)>0){
7738 7750
         torm_index <- c()
7739 7751
         for(ssx in 1:length(torm)){
7740 7752
           # ssx <- 2
7741
-          subexontoremove <- gsub(
7742
-            "\\..*","",names(torm)[ssx])
7743
-          torm_index <- c(
7744
-            torm_index,
7745
-            grep(subexontoremove,
7746
-                 rownames(newAdj)))
7753
+          subexontoremove <- names(torm)[ssx]
7754
+          torm_index <- c(torm_index,which(rownames(newAdj)==subexontoremove))
7747 7755
         }
7748
-        newAdj <- newAdj[
7749
-          -torm_index,-torm_index]     
7756
+        newAdj <- newAdj[-torm_index,-torm_index]     
7750 7757
       }else{
7751 7758
         break
7752 7759
       }
... ...
@@ -7767,6 +7774,12 @@ reclasify_intern <- function(SG,mievento,pp1,pp2,ppref){
7767 7774
     #possible paths:
7768 7775
     numberOfPaths <- solve(
7769 7776
       Diagonal(ncol(newAdj))-newAdj)
7777
+    
7778
+    if(numberOfPaths[
7779
+      ref_junct_ford,ref_junct_revs] > 10){
7780
+      return("Complex Event")
7781
+    }
7782
+    
7770 7783
     distances_list <- vector(
7771 7784
       mode="list",
7772 7785
       length=numberOfPaths[
... ...
@@ -22,7 +22,6 @@
22 22
 #' @param SplicingGraph A list with the splicing graph of all the genes of
23 23
 #'  a reference transcriptome. This data is returned
24 24
 #'   by the function EventDetection_transcriptome.
25
-#' @param cores Number of cores using in the parallel processing (by default = 1)
26 25
 #'
27 26
 #' @return A data.frame containing a new column with the new classification ('EventType_new'):
28 27
 #' 
... ...
@@ -39,8 +38,7 @@
39 38
 #' #this table has the information of 5 complex events.
40 39
 #' 
41 40
 #' EventTable_new <- Events_ReClassification(EventTable = EventTable,
42
-#'                                           SplicingGraph = SG_reclassify,
43
-#'                                           cores = 1)
41
+#'                                           SplicingGraph = SG_reclassify)
44 42
 #'          
45 43
 #'     
46 44
 #'
... ...
@@ -60,7 +58,7 @@
60 58
 #' @importFrom S4Vectors queryHits subjectHits split
61 59
 #' @importFrom IRanges relist disjoin %over% reduce
62 60
 
63
-Events_ReClassification <- function(EventTable,SplicingGraph,cores = 1){
61
+Events_ReClassification <- function(EventTable,SplicingGraph){
64 62
   
65 63
   if (is.null(EventTable)) {
66 64
     stop("not EventTable")
... ...
@@ -76,11 +74,13 @@ Events_ReClassification <- function(EventTable,SplicingGraph,cores = 1){
76 74
   if(length(uux)==0){
77 75
     return(EventTable)
78 76
   }
79
-  
80
-  
81
-  registerDoParallel(cores = cores)
82
-  Result <- foreach(jj = seq_len(length(uux))) %dopar%
83
-    {
77
+  cat("\nStarting reclassification\n")
78
+  pb <- txtProgressBar(min = 0, max = length(uux), 
79
+                       style = 3)
80
+  misnewtype <- c()
81
+    for(jj in seq_len(length(uux))){
82
+      setTxtProgressBar(pb, jj)
83
+      # cat(jj,"\n")
84 84
       mievento <- EventTable$ID[uux][jj]
85 85
       # mievento
86 86
       ggx <- match(EventTable$GeneName[uux][jj],names(SplicingGraph))
... ...
@@ -90,14 +90,17 @@ Events_ReClassification <- function(EventTable,SplicingGraph,cores = 1){
90 90
       pp1 <- EventTable$Path.1[uux][jj]
91 91
       pp2 <- EventTable$Path.2[uux][jj]
92 92
       ppref <- EventTable$Path.Reference[uux][jj]
93
-      return(reclasify_intern(SG = SG,mievento = mievento,
94
-                              pp1=pp1,pp2=pp2,ppref=ppref))
93
+      misnewtype <- c(misnewtype,reclasify_intern(SG = SG,mievento = mievento,
94
+                                                  pp1=pp1,pp2=pp2,ppref=ppref))
95 95
       
96 96
     }
97 97
   
98
-  Result <- unlist(Result)
98
+  close(pb)
99
+  cat("\nReclassification Finished\n")
100
+  # Result <- unlist(Result)
101
+  # EventTable$EventType_new[uux] <- Result
99 102
   
100
-  EventTable$EventType_new[uux] <- Result
103
+  EventTable$EventType_new[uux] <- misnewtype
101 104
   
102 105
   return(EventTable)
103 106
   
... ...
@@ -4,7 +4,7 @@
4 4
 \alias{Events_ReClassification}
5 5
 \title{Events_ReClassification}
6 6
 \usage{
7
-Events_ReClassification(EventTable, SplicingGraph, cores = 1)
7
+Events_ReClassification(EventTable, SplicingGraph)
8 8
 }
9 9
 \arguments{
10 10
 \item{EventTable}{Table returned by EventDetection_transcriptome. 
... ...
@@ -13,8 +13,6 @@ Can be easily loaded using the function read.delim as data.frame.}
13 13
 \item{SplicingGraph}{A list with the splicing graph of all the genes of
14 14
 a reference transcriptome. This data is returned
15 15
  by the function EventDetection_transcriptome.}
16
-
17
-\item{cores}{Number of cores using in the parallel processing (by default = 1)}
18 16
 }
19 17
 \value{
20 18
 A data.frame containing a new column with the new classification ('EventType_new'):
... ...
@@ -46,8 +44,7 @@ EventTable <- read.delim(file=inputFile)
46 44
 #this table has the information of 5 complex events.
47 45
 
48 46
 EventTable_new <- Events_ReClassification(EventTable = EventTable,
49
-                                          SplicingGraph = SG_reclassify,
50
-                                          cores = 1)
47
+                                          SplicingGraph = SG_reclassify)
51 48
          
52 49
     
53 50
 
... ...
@@ -1033,7 +1033,7 @@ Thus, EventPointer reclassifies the complex events according to how similar the
1033 1033
 
1034 1034
 The following code shows how this function works for a subset of 5 splicing events:
1035 1035
 
1036
-```{r reclassify, eval=TRUE}
1036
+```{r reclassification, eval=TRUE}
1037 1037
 #load splicing graph
1038 1038
 data("SG_reclassify")
1039 1039
 
... ...
@@ -1044,54 +1044,12 @@ EventTable <- read.delim(file=inputFile)
1044 1044
 #this table has the information of 5 complex events.
1045 1045
 
1046 1046
 EventTable_new <- Events_ReClassification(EventTable = EventTable,
1047
-                                          SplicingGraph = SG_reclassify,
1048
-                                          cores = 1)
1047
+                        SplicingGraph = SG_reclassify)
1049 1048
 ```
1050 1049
 
1051 1050
 The following figures show the structure of these events and their corresponding gene:
1052 1051
 
1053 1052
 **Example 1** 
1054
-![**Example 1** Complex Event reclassified as alternative 5' splice site and as complex casstte exon  ](example_1_alt5_ce.png)
1055
-```{r example1, eval=TRUE}
1056
-EventTable_new$EventType_new[[1]]
1057
-```
1058
-<br><br><br>
1059
-**Example 2** 
1060
-![**Example 2** Complex Event reclassified as Multiple Exon Skipping  ](example_2_multiple_exon_skipping.png)
1061
-```{r example2, eval=TRUE}
1062
-EventTable_new$EventType_new[[2]]
1063
-```
1064
-
1065
-<br><br><br>
1066
-**Example 3** 
1067
-![**Example 3** Complex Event reclassified as Complex Cassette Exon  ](example3_ce.png)
1068
-```{r example3, eval=TRUE}
1069
-EventTable_new$EventType_new[[3]]
1070
-```
1071
-
1072
-<br><br><br>
1073
-**Example 4** 
1074
-![**Example 4** Complex Event reclassified as Retained Intron and Complex Cassette Exon  ](example_4_ri_ce.png)
1075
-```{r example4, eval=TRUE}
1076
-EventTable_new$EventType_new[[4]]
1077
-```
1078
-
1079
-<br><br><br>
1080
-**Example 5** 
1081
-![**Example 5** Complex Event reclassified as Complex Cassette Exon and complex Alternative 3' Splice Site  ](example_5_ce_alt3.png)
1082
-```{r example5, eval=TRUE}
1083
-EventTable_new$EventType_new[[5]]
1084
-```
1085
-
1086
-
1087
-Further, files needed to load these pictures in IGV are available:
1088
-
1089
-```{r source_pictures,eval=FALSE}
1090
-PathFiles<-system.file("extdata",package="EventPointer")
1091
-paste0(PathFiles,"/example_recla_events.gtf")
1092
-paste0(PathFiles,"/example_recla_genes.gtf")
1093
-```
1094
-
1095 1053
 
1096 1054
 ## Statistical Tests
1097 1055