made changes to take into account the new cases in the test_data_2

deletion/EJR_junction_auxfns.py

@@ -436,6 +436,8 @@ def get_deln_indices(ori_seq, deleted_seq):
         j += 1
     if j < len(ori_seq)-1:
         print "original seq and deletion not properly aligned(C)"
+        print(ori_seq)
+        print(deleted_seq)
         return None
     return (del_index1, del_index2)
 

deletion/SDMMEJDeletionProgram_cli.py

@@ -19,6 +19,7 @@ def parse_args():
     return args.h_in, args.d_in, args.dir_path, int(args.nick), int(args.search_radius), int(args.min_consist_rep_size)
 
 def parse_ref(h_in, nick):
+    print(h_in)
     hifi_input = pd.read_csv(h_in)
 
     reference = hifi_input.loc[hifi_input['CLASS'] == 'exact']
@@ -45,6 +46,7 @@ def parse_ref(h_in, nick):
     return original_seq, left_dna, mid_dna, right_dna, rightnick
 
 h_in, d_in, dir_path, nick, search_radius, min_consist_rep_size = parse_args()
+print(h_in)
 original_seq, left_dna, mid_dna, right_dna, rightnick = parse_ref(h_in, nick)
 
 # get and process data file
@@ -55,6 +57,10 @@ df = open(data_file)
 base_name = os.path.basename(data_file).strip('.txt')
 
 
+# make output dir
+if not os.path.exists(dir_path):
+    os.makedirs(dir_path)
+
 # get output filename
 outfile = dir_path + "/" + base_name + "_consistency_log.txt"
 print(outfile)

insertion/INSERTION_PROGRAM.R

@@ -16,26 +16,33 @@ suppressPackageStartupMessages(library(stringr))
 suppressPackageStartupMessages(library(Biostrings))
 
 ##test if there is at least one argument: if not, return an error
-# if (length(args)<4) {
-#   print(length(args))
-#   stop("Usage: Rscript INSERTION_PROGRAM.R hifibr_reclassified.csv insertion_file.txt nick_location outdir", call.=FALSE)
-# } else if (length(args)==4) {
-#   hifi_in = args[1]
-#   insertion_in = args[2]
-#   nick = as.integer(args[3])
-#   out_dir = args[4]
-# }
+if (length(args)<5) {
+  print(length(args))
+  stop("Usage: Rscript INSERTION_PROGRAM.R hifibr_reclassified.csv insertion_file.txt nick_location outdir", call.=FALSE)
+} else if (length(args)==5) {
+  hifi_in = args[1]
+  insertion_in = args[2]
+  out_dir = args[3]
+  nick = as.integer(args[4])
+  search_radius=as.integer(args[5])
+  
+}
 
 #out_dir = "~/Documents/git/sdmmej/test_data/TestData_HiFiBR_Output_mod_output/"
 #hifi_in = paste0(out_dir,"TestData_HiFiBR_Output_mod_reclassified.csv")
 #insertion_in = paste0(out_dir,"TestData_HiFiBR_Output_mod_insertion.txt")
 #nick=71
 
-out_dir = "~/Documents/git/sdmmej/test_data/expected_results/"
-hifi_in = "~/Documents/git/sdmmej/test_data/TestData_HiFiBR_Output_mod_output/TestData_HiFiBR_Output_mod_reclassified.csv"
-insertion_in = paste0(out_dir,"TestData_Insertions.csv")
-nick=71
+# out_dir = "~/Documents/git/sdmmej/test_data/expected_results/"
+# hifi_in = "~/Documents/git/sdmmej/test_data/TestData_HiFiBR_Output_mod_output/TestData_HiFiBR_Output_mod_reclassified.csv"
+# insertion_in = paste0(out_dir,"TestData_Insertions.csv")
+# nick=71
 
+# out_dir = "~/Box/rebecca_documents/sdmmej_data/"
+# hifi_in = "~/Box/rebecca_documents/sdmmej_data/PolyA1Seq_reclassified.csv"
+# insertion_in = "~/Box/rebecca_documents/sdmmej_data/PolyA1Seq_insertion.txt"
+# nick=161
+# 
 a<-read.csv(insertion_in)
 
 ## get reference
@@ -43,13 +50,14 @@ hifibr_input = read.csv(hifi_in)
 ref = hifibr_input %>% dplyr::filter(., CLASS_final == 'exact')
 ref = ref$ALIGNED_SEQ
 
-n <- 30  # number of bases to the left and right of the break you want to search of repeated motifs
+n <- search_radius  # number of bases to the left and right of the break you want to search of repeated motifs
 p <- 10 # number of bases to the left and right of the initial repeat motif to search for homology
 #plasmid <- "test_data_m6"
 plasmid <- gsub("_reclassified.csv","",basename(hifi_in))
 
 #L <- "GGAAAAAATTCGTACTTTGGAGTACGAAATGCGTCGTTTAGAGCAGCAGCCGAATTCGGTACATTACCCTGTTAT"
 L = substr(ref,0,nick)
+
 #R <- "TTATCCCTAGCTATGGTCTGCGCTACTAGTGGATCTGGGGCCGCATAGGCCATCCTCTAGAGTCGACCTCGAACGTAAACGTTAACGTAACGTTAACTCG"
 l=length(ref)
 R <- substr(ref,nick+1,nchar(ref)) 
@@ -57,6 +65,7 @@ R <- substr(ref,nick+1,nchar(ref))
 if (substr(R,0,4) == 'TTAT'){
  L =  paste0(L,'TAT')
 }
+L
 l <- nchar(L) # number of nucleotides of the left hand sequence
 r <- nchar(R) # number of nucleotides of the right hand sequence
 k1 <- 30 # how far you want to cut back to search, this needs to to be adjusted based on how large the deletions are, but i think this covers up to 30 bp of deletion on either side.
@@ -71,6 +80,15 @@ for (i in a[, 1]){
   lbb <- lb[nrow(lb),2]
   a2[i] = lbb
 }
+
+### testing
+#sigma <- nucleotideSubstitutionMatrix(match = 2, mismatch = -1, baseOnly = FALSE)
+#seq=i
+#align <- pairwiseAlignment(seq, ref, substitutionMatrix = sigma, gapOpening = -2,
+#                           gapExtension = -8, scoreOnly = FALSE)
+### 
+
+
 a3=NULL    # create empty vector to insert RIGHT del boundary
 for (i in a[, 1]){
   rb <- str_locate(as.character(i), sR1)

process_hifibr.R

@@ -4,21 +4,41 @@ args = commandArgs(trailingOnly=TRUE)
 suppressPackageStartupMessages(library(tidyverse))
 suppressPackageStartupMessages(library(Biostrings))
 
-# test if there is at least one argument: if not, return an error
-if (length(args)<2) {
-  stop("Usage: Rscript process_hifibr.R input_file.csv outdir", call.=FALSE)
-} else if (length(args)==2) {
+# test if there is at least 5 argument: if not, return an error
+if (length(args)<5) {
+  stop("Usage: Rscript process_hifibr.R input_file.csv outdir search_radius debug", call.=FALSE)
+} else if (length(args)==5) {
   input_file = args[1]
   out_dir = args[2]
+  search_radius=as.integer(args[3])
+  break_location=as.integer(args[4])
+  debug=as.integer(args[5])
 }
 
+print(paste0("Input file: " ,input_file))
+print(paste0("Output dir: ", out_dir))
+print(paste0("Search radius: ", search_radius))
+print(paste0("Break location: ", break_location))
+print(paste0("Debug: ", debug))
+
+## functions
+
+## 
+# input_file="~/Documents/git/sdmmej/test_data_2/PolyA1Seq_questions.csv"
+# out_dir="~/Documents/git/sdmmej/test_data_2/"
+# search_radius=30
+# break_location=161
+# debug=1
+
 #input_file="~/Documents/git/sdmmej/test_data/TestData_HiFiBR_Output_mod.csv"
+
+
 hifibr_input = read.csv(input_file)
 
 ## add an id col and col for reclassified CLASS
 hifibr_input <- tibble::rowid_to_column(hifibr_input, "ID")
 hifibr_input$CLASS_final = ifelse(hifibr_input$CLASS %in% c('deletion','insertion','exact') , hifibr_input$CLASS, NA)
-hifibr_input$CLASS_final = ifelse(hifibr_input$`X._OF_READS` < 10 & hifibr_input$CLASS != 'exact', 'filtered',hifibr_input$CLASS)
+hifibr_input$CLASS_final = ifelse(hifibr_input$READS < 10 & hifibr_input$CLASS != 'exact', 'filtered',hifibr_input$CLASS)
 
 ## names for output
 in_string  = basename(input_file)
@@ -34,11 +54,10 @@ if ( n_ref != 1){
   quit(save = "no", status = 1, runLast = FALSE)
 }
 
-
 # separate sequences and filter for length
-hifibr_input_filter_del = hifibr_input %>% dplyr::filter(., `X._OF_READS` > 10 & CLASS == 'deletion')
-hifibr_input_filter_ins = hifibr_input %>% dplyr::filter(., `X._OF_READS` > 10 & CLASS == 'insertion')
-hifibr_input_filter_complex = hifibr_input %>% dplyr::filter(., `X._OF_READS` > 10 & CLASS == 'complex')
+hifibr_input_filter_del = hifibr_input %>% dplyr::filter(.,READS > 10 & CLASS == 'deletion')
+hifibr_input_filter_ins = hifibr_input %>% dplyr::filter(.,READS > 10 & CLASS == 'insertion')
+hifibr_input_filter_complex = hifibr_input %>% dplyr::filter(.,READS > 10 & CLASS == 'complex')
 
 # process deletions
 deletions = hifibr_input_filter_del$ALIGNED_SEQ
@@ -46,7 +65,7 @@ deletions = hifibr_input_filter_del$ALIGNED_SEQ
 aligned_deletions = c()
 
 for (seq in deletions){
-  sigma <- nucleotideSubstitutionMatrix(match = 2, mismatch = -1, baseOnly = TRUE)
+  sigma <- nucleotideSubstitutionMatrix(match = 2, mismatch = -1, baseOnly = FALSE)
   align <- pairwiseAlignment(seq, ref_seq, substitutionMatrix = sigma, gapOpening = -2,
                                        gapExtension = -8, scoreOnly = FALSE)
   align_string <- toString(align@pattern)
@@ -65,51 +84,188 @@ insertions = hifibr_input_filter_ins$ALIGNED_SEQ
 # process complex
 comp_unknown = c()
 
+# ### debug
+# seq = hifibr_input_filter_complex[2, 'ALIGNED_SEQ']
+# id = hifibr_input_filter_complex[2,'ID']
+# print(seq)
+# print(id)
+# sigma <- nucleotideSubstitutionMatrix(match = 2, mismatch = -1, baseOnly = FALSE)
+# align <- pairwiseAlignment(seq, ref_seq, substitutionMatrix = sigma, gapOpening = -2,
+#                            gapExtension = -8, scoreOnly = FALSE)
+# 
+# ins_info = insertion(align)[[1]]
+# del_info = deletion(align)[[1]]
+# mis_info = mismatchTable(align)
+# 
+# ins_n = length(ins_info)
+# ins_pos = start(ins_info)
+# ins_r = abs(break_location - ins_pos)
+# 
+# del_n = length(del_info)
+# del_pos = start(del_info)
+# del_r = abs(break_location - del_pos)
+# 
+# mis_n = nmismatch(align)
+# mis_pos = mis_info$SubjectStart
+# mis_r = abs(break_location - mis_pos)
+# 
+# 
+# print(c(ins_n,ins_pos, ins_r))
+# print(c(del_n,del_pos, del_r))
+# print(c(mis_n,mis_pos, mis_r))
+# 
+# ## it seems they are always IRanges objects of length 1
+# label = ""
+# 
+# ## check number of deletions
+# deletion=deletion(align)[[1]]
+# n_del=length(deletion)
+# print(n_del)
+# # check
+# if (n_del >1){
+#   label = "complex"  
+# }
+# 
+# # check position of deletions
+# pos_del=start(deletion)[1]
+# r_del=abs(break_location-pos_del)
+# print(r_del)
+# if (r_del > search_radius){
+#   label = "complex"  
+# }
+# 
+# # check number of insertions
+# insertion=insertion(align)[[1]]
+# n_ins=length(insertion)
+# print(n_ins)
+# if (n_ins >1){
+#   label="complex"  
+# }
+# 
+# # check position of insertions
+# pos_ins=start(insertion)[1]
+# print(pos_ins)
+# 
+# r_ins=abs(break_location-pos_ins)
+# print(r_ins)
+# # check
+# if (r_ins > search_radius){
+#   label="complex" 
+# }
+# 
+# # check position of mismatch, assuming it is only one
+# pos_mis=mismatchTable(align)$SubjectStart
+# r_mis=abs(break_location-pos_mis)
+# 
+# if (r_mis > search_radius){
+#   label="complex" 
+# }
+# 
+# 
+# print(label)
+
+### 
 for (row in 1:nrow(hifibr_input_filter_complex)){
   seq = hifibr_input_filter_complex[row, 'ALIGNED_SEQ']
   id = hifibr_input_filter_complex[row,'ID']
-  #print(id)
+  if(debug == 1){
+    print(id)
+  }
+  sigma <- nucleotideSubstitutionMatrix(match = 2, mismatch = -1, baseOnly = FALSE)
   align <- pairwiseAlignment(seq, ref_seq, substitutionMatrix = sigma, gapOpening = -2,
                              gapExtension = -8, scoreOnly = FALSE)
   
-  ins_width = length(width(insertion(align))[[1]])
-  del_width = length(width(deletion(align))[[1]])
-  nmismatch = nmismatch(align)
+  if(debug == 1){
+    print(insertion(align))
+    print(deletion(align))
+    print(mismatchTable(align))
+  }
   
-  if (ins_width == 0 & del_width > 0 & nmismatch == 0){
-    align_string <- toString(align@pattern)
-    aligned_deletions = c(aligned_deletions, align_string)
-    hifibr_input[hifibr_input$ID == id,]$CLASS_final = 'deletion'
-    
-    #print("found del")
-    #print("ref")
-    #print(toString(align@subject))
-    #print("seq")
-    #print(toString(align@pattern))
-  }else if (ins_width > 0 & del_width == 0 & nmismatch == 0){
-    insertions = c(insertions, seq)
-    hifibr_input[hifibr_input$ID == id,]$CLASS_final = 'insertion'
-    
-    # print("found ins")
-    # print("ref")
-    # print(toString(align@subject))
-    # print("seq")
-    # print(toString(align@pattern))
-  }else if (ins_width == 0 & del_width == 0 & nmismatch == 1){
-    insertions = c(insertions, seq)
-    hifibr_input[hifibr_input$ID == id,]$CLASS_final = 'insertion'
+  ins_info = insertion(align)[[1]]
+  del_info = deletion(align)[[1]]
+  mis_info = mismatchTable(align)
+  
+  ins_n = length(ins_info)
+  ins_pos = start(ins_info)
+  ins_r = abs(break_location - ins_pos)
     
-    # print("found single mismatch ins")
-    # print("ref")
-    # print(toString(align@subject))
-    # print("seq")
-    # print(toString(align@pattern))
+  del_n = length(del_info)
+  del_pos = start(del_info)
+  del_r = abs(break_location - del_pos)
+  
+  mis_n = nmismatch(align)
+  mis_pos = mis_info$SubjectStart
+  mis_r = abs(break_location - mis_pos)
+  
+  # check if it's a single deletion within the search radius
+  if (ins_n == 0 & mis_n == 0 & del_n ==1 ){
+    if (del_r < search_radius){
+      align_string <- toString(align@pattern)
+      aligned_deletions = c(aligned_deletions, align_string)
+      hifibr_input[hifibr_input$ID == id,]$CLASS_final = 'deletion'
+      
+      if(debug == 1){
+        print("found del")
+        print("ref")
+        print(toString(align@subject))
+        print("seq")
+        print(toString(align@pattern))
+      }
+    }
+  }else if (mis_n == 0 & del_n ==0 & ins_n == 1 ){
+    if (ins_r < search_radius){
+      insertions = c(insertions, seq)
+      hifibr_input[hifibr_input$ID == id,]$CLASS_final = 'insertion'
+      
+      if(debug == 1){
+        print("found ins")
+        print("ref")
+        print(toString(align@subject))
+        print("seq")
+        print(toString(align@pattern))
+      }
+    }else{
+      if(debug == 1){
+        print("found complex (del outside search radius)")
+        print("ref")
+        print(toString(align@subject))
+        print("seq")
+        print(toString(align@pattern))
+      }
+      comp_unknown = c(comp_unknown, seq)
+      hifibr_input[hifibr_input$ID == id,]$CLASS_final = 'complex'
+    }
+  }else if (del_n ==0 & ins_n == 0 & mis_n == 1){
+    if (mis_r < search_radius){
+      insertions = c(insertions, seq)
+      hifibr_input[hifibr_input$ID == id,]$CLASS_final = 'insertion'
+      
+      if(debug == 1){
+        print("found single mismatch ins")
+        print("ref")
+        print(toString(align@subject))
+        print("seq")
+        print(toString(align@pattern))
+      }
+    }else{
+      if(debug == 1){
+        print("found complex (ins outside search radius)")
+        print("ref")
+        print(toString(align@subject))
+        print("seq")
+        print(toString(align@pattern))
+      }
+      comp_unknown = c(comp_unknown, seq)
+      hifibr_input[hifibr_input$ID == id,]$CLASS_final = 'complex'
+    }
   }else{
-    # print("found complex")
-    # print("ref")
-    # print(toString(align@subject))
-    # print("seq")
-    # print(toString(align@pattern))
+    if(debug == 1){
+      print("found complex")
+      print("ref")
+      print(toString(align@subject))
+      print("seq")
+      print(toString(align@pattern))
+    }
     comp_unknown = c(comp_unknown, seq)
     hifibr_input[hifibr_input$ID == id,]$CLASS_final = 'complex'
   }
@@ -130,3 +286,4 @@ write(comp_unknown, file = paste0(out_dir,"/", out_string,"_complex.txt"),append
 
 ## output file
 write.csv(hifibr_input, file = paste0(out_dir,"/", out_string,"_reclassified.csv"), row.names = FALSE)
+

run_script.sh

+#!/bin/bash
+
+infile=/Users/rbator01/Box/rebecca_documents/sdmmej_data/PolyGSeq.csv
+indir=/Users/rbator01/Box/rebecca_documents/sdmmej_data/
+search_radius=30
+breakpoint=161
+debug=0
+
+reclass=${infile%.csv}_reclassified.csv
+deletion=${infile%.csv}_deletion.txt
+insertion=${infile%.csv}_insertion.txt
+
+Rscript process_hifibr.R $infile $indir $search_radius $breakpoint $debug
+
+cd deletion/
+python SDMMEJDeletionProgram_cli.py -hi $reclass -del $deletion -n $breakpoint -out $indir
+
+cd ../insertion/
+Rscript INSERTION_PROGRAM.R $reclass $insertion $indir $breakpoint $search_radius
+
+

test_data_2/PolyA1Seq_questions.csv

+UMI,CIGAR_STRING,READ_LENGTH,SPLIT_CIGAR_STRING,MATCH_LEFT,MATCH_RIGHT,DISTANCE_FROM_BREAK_LEFT,DISTANCE_FROM_BREAK_RIGHT,DELETION_FROM_LEFT,DELETION_FROM_RIGHT,TOTAL_DELETION,INSERTION_START,INSERTION_END,INSERTION_LENGTH,INSERTED_SEQ,CLASS,ALIGNED_SEQ,READS,MICROHOMOLOGY,MH_Length,NUMBER_OF_ALIGNMENTS,MISMATCH_PERCENTAGE_TO_RECONSTRUCTED
+AATTGGAAAA_CTCATGTCCT,325M,325,325M,325,325,164,161,0,0,0,161,164,0,,exact,GATCCTAGGATCCTAGGATCCTAGGATCCTAGGATCCTAGGATCCTAGGATCCTAGGATCCTAGGATCCTAGGATCCTAGGAGGGAAAAAATTCGTACTTTGGAGTACGAAAATTGGAAAATAGAGCAGCACTCACCTTATTGTCATTACCCTGTTATCCCTAGAGGCCAAACAGGCCGGCGCCTCCTAACGATCCTCTAGCTCATGTCCTGAACGTTAACGTTAACGTAACGTTAACTCGAGGCCCAGGACCCCAGGACCCCAGGACCCCAGGACCCCAGGACCCCAGGACCCCAGGACCCCAGGACCCCAGGACCCCAGGACC,214208,,NA,1,1.847269943
+AATTGGAAAA_CTCATGTCCT,325M,325,325M,150,164,-11,0,-11,0,-11,150,164,11,NCTGTTATCCC,complex,GATCCTAGGATCCTAGGATCCTAGGATCCTAGGATCCTAGGATCCTAGGATCCTAGGATCCTAGGATCCTAGGATCCTAGGAGGGAAAAAATTCGTACTTTGGAGTACGAAAATTGGAAAATAGAGCAGCACTCACCTTATTGTCATTACNCTGTTATCCCTAGAGGCCAAACAGGCCGGCGCCTCCTAACGATCCTCTAGCTCATGTCCTGAACGTTAACGTTAACGTAACGTTAACTCGAGGCCCAGGACCCCAGGACCCCAGGACCCCAGGACCCCAGGACCCCAGGACCCCAGGACCCCAGGACCCCAGGACCCCAGGACC,24,,,1,8.333333333
+AATTGGAAAA_CTCATGTCCT,160M3D162M,322,160M3D162M,160,161,-1,-3,-1,-3,-4,160,161,1,A,complex,GATCCTAGGATCCTAGGATCCTAGGATCCTAGGATCCTAGGATCCTAGGATCCTAGGATCCTAGGATCCTAGGATCCTAGGAGGGAAAAAATTCGTACTTTGGAGTACGAAAATTGGAAAATAGAGCAGCACTCACCTTATTGTCATTACCCTGTTATCCAAGGCCAAACAGGCCGGCGCCTCCTAACGATCCTCTAGCTCATGTCCTGAACGTTAACGTTAACGTAACGTTAACTCGAGGCCCAGGACCCCAGGACCCCAGGACCCCAGGACCCCAGGACCCCAGGACCCCAGGACCCCAGGACCCCAGGACCCCAGGACC,40,,NA,1,2.5
+AATTGGAAAA_CTCATGTCCT,117M1I208M,326,117M1I208M,121,208,-40,44,-40,0,-40,121,164,41,ATAGAGCAGCACTCACCTTATTGTCATTACCCTGTTATCCC,complex,GATCCTAGGATCCTAGGATCCTAGGATCCTAGGATCCTAGGATCCTAGGATCCTAGGATCCTAGGATCCTAGGATCCTAGGAGGGAAAAAATTCGTACTTTGGAGTACGAAAATTGGAAAAATAGAGCAGCACTCACCTTATTGTCATTACCCTGTTATCCCTAGAGGCCAAACAGGCCGGCGCCTCCTAACGATCCTCTAGCTCATGTCCTGAACGTTAACGTTAACGTAACGTTAACTCGAGGCCCAGGACCCCAGGACCCCAGGACCCCAGGACCCCAGGACCCCAGGACCCCAGGACCCCAGGACCCCAGGACCCCAGGACC,15,,,1,0
\ No newline at end of file