generation.sh

@@ -18,5 +18,4 @@ prefix="curesim_"
 
 Rscript src/repertoire.r "$sequences" "$number_of_reads" &&
 	java -jar tools/CuReSim.jar -f "$data_directory$filename$fastq" -o "$data_directory$prefix$filename$fastq"
-Rscript src/alignment.r
 rm "$data_directory/log.txt"

src/alignment.r

@@ -70,32 +70,31 @@ get_cys_coordinates <- function(alignment) {
   return(list("start" = cys_start, "end" = cys_end))
 }
 
-get_hvr_sequences <- function(sequences, vdj_segments, cores = detectCores()) {
+# TODO Refactor this mess
+get_hvr_sequences <- function(sequences, vdj_segments) {
   df <- fetch_vj_sequences(sequences, vdj_segments)
-  v_alignment <- parallel::mcmapply(sequences,
+  v_alignment <- parallel::mcmapply(sequences, df$v_seq, FUN = align_sequence)
-    df$v_seq,
+  j_alignment <- parallel::mcmapply(sequences, df$j_seq, FUN = align_sequence)
-    FUN = align_sequence,
-    mc.cores = cores
-  )
   cys_coordinates <- parallel::mclapply(v_alignment, FUN = get_cys_coordinates)
   cys_df <- as.data.frame(do.call(rbind, cys_coordinates))
-  remaining <- Biostrings::subseq(sequences, start = unlist(cys_df$end))
-  j_alignment <- parallel::mcmapply(remaining,
-    df$j_seq,
-    FUN = align_sequence,
-    mc.cores = cores
-  )
   j_start <- parallel::mclapply(
     j_alignment,
-    function(x) start(Biostrings::Views(x)),
+    function(x) start(Biostrings::Views(x))
-    mc.cores = cores
+  )
+  hvr <- Biostrings::subseq(sequences,
+    start = unlist(cys_df$start),
+    end = unlist(j_start) + 2
   )
-  hvr_start <- unlist(cys_df$start)
-  hvr_end <- unlist(cys_df$start) + unlist(j_start) + 2
-  hvr <- Biostrings::subseq(sequences, start = hvr_start, end = hvr_end)
   return(hvr)
 }
 
+save_data <- function(data) {
+  Biostrings::writeXStringSet(data, "data/CuReSim-HVR.fastq", format = "fastq")
+}
+
 data <- parse_data(file = "data/curesim_sequence.fastq")
-hvr <- get_hvr_sequences(sequences = data[[1]], vdj_segments = data[[2]])
+hvr <- get_hvr_sequences(
+  sequences = data[[1]],
+  vdj_segments = data[[2]]
+)
 Biostrings::writeXStringSet(hvr, "data/CuReSim-HVR.fastq", format = "fastq")