Add usage instructions to the README

README.org

@@ -1,3 +1,49 @@
 * locigenesis
 
 locigenesis is a tool that generates an immune repertoire and runs it through a sequence reader simulation tool, to generate sequencing errors.
+
+** Installation
+
+This project uses [[https://nixos.org/][Nix]] to ensure reproducible builds.
+
+1. Install Nix (compatible with MacOS, Linux and [[https://docs.microsoft.com/en-us/windows/wsl/about][WSL]]):
+
+#+begin_src shell
+curl -L https://nixos.org/nix/install | sh
+#+end_src
+
+1. Clone the repository:
+
+#+begin_src shell
+git clone https://git.coolneng.duckdns.org/coolneng/locigenesis
+#+end_src
+
+3. Change the working directory to the project:
+
+#+begin_src shell
+cd locigenesis
+#+end_src
+
+4. Enter the nix-shell:
+
+#+begin_src shell
+nix-shell
+#+end_src
+
+After running these commands, you will find yourself in a shell that contains all the needed dependencies.
+
+** Usage
+
+An execution script that accepts 2 parameters is provided, the following command invokes it:
+
+#+begin_src shell
+./generation.sh <number of sequences> <number of reads>
+#+end_src
+
+- <number of sequences>: an integer that specifies the number of different sequences to generate
+- <number of reads>: an integer that specifies the number of reads to perform on each sequence
+
+The script will generate 2 files under the data directory:
+
+| HVR.fastq         | Contains the original CDR3 sequence                             |
+| CuReSim-HVR.fastq | Contains CDR3 after the read simulation, with sequencing errors |

generation.sh

@@ -1,7 +1,7 @@
 #!/bin/sh
 
 usage() {
-	echo "usage: generation.sh <number of sequences> <number_of_reads>"
+	echo "usage: generation.sh <number of sequences> <number of reads>"
 	exit 1
 }
 
@@ -17,5 +17,6 @@ filename="sequence"
 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"
+	CuReSim -f "$data_directory$filename$fastq" -o "$data_directory$prefix$filename$fastq"
+Rscript src/alignment.r
 rm "$data_directory/log.txt"

shell.nix

@@ -2,15 +2,35 @@
 
 with pkgs;
 
-mkShell {
+let
+  CuReSim = stdenv.mkDerivation rec {
+    name = "CuReSim";
+    version = "1.3";
+    src = fetchzip {
+      url =
+        "http://www.pegase-biosciences.com/wp-content/uploads/2015/08/${name}${version}.zip";
+      sha256 = "1hvlpgy4haqgqq52mkxhcl9i1fx67kgwi6f1mijvqzk0xff77hkp";
+      stripRoot = true;
+      extraPostFetch = ''
+        chmod go-w $out
+      '';
+    };
+
+    nativeBuildInputs = [ makeWrapper ];
+
+    installPhase = ''
+      mkdir -pv $out/share/java $out/bin
+      cp -r ${src} $out/share/java/${name}
+      makeWrapper ${pkgs.jdk}/bin/java $out/bin/CuReSim --add-flags "-jar $out/share/java/${name}/${name}.jar"
+    '';
+  };
+in mkShell {
   buildInputs = [
     R
     rPackages.immuneSIM
     rPackages.Biostrings
     rPackages.stringr
     jdk
-    # Development tools
-    rPackages.languageserver
-    rPackages.lintr
+    CuReSim
   ];
 }

src/alignment.r

@@ -70,31 +70,32 @@ get_cys_coordinates <- function(alignment) {
   return(list("start" = cys_start, "end" = cys_end))
 }
 
-# TODO Refactor this mess
-get_hvr_sequences <- function(sequences, vdj_segments) {
+get_hvr_sequences <- function(sequences, vdj_segments, cores = detectCores()) {
   df <- fetch_vj_sequences(sequences, vdj_segments)
-  v_alignment <- parallel::mcmapply(sequences, df$v_seq, FUN = align_sequence)
-  j_alignment <- parallel::mcmapply(sequences, df$j_seq, FUN = align_sequence)
+  v_alignment <- parallel::mcmapply(sequences,
+    df$v_seq,
+    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))
-  )
-  hvr <- Biostrings::subseq(sequences,
-    start = unlist(cys_df$start),
-    end = unlist(j_start) + 2
+    function(x) start(Biostrings::Views(x)),
+    mc.cores = cores
   )
+  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")