Add usage instructions to the README

.gitignore

@@ -0,0 +1 @@
+*.fastq

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 |

docs/notebook.org

@@ -0,0 +1,46 @@
+#+TITLE: locigenesis
+#+AUTHOR: Amin Kasrou Aouam
+#+DATE: 2021-03-10
+* Sequence alignment
+
+Our generated sequences contain the full VJ region, but we are only interested in the CDR3 (Complementarity-determining region). We will proceed by delimiting CDR3, using the known sequences of V and J.
+
+#+begin_src R :results value silent
+v_segments <- readRDS("data/v_segments.rds")
+j_segments <- readRDS("data/j_segments_phe.rds")
+#+end_src
+
+#+begin_src R
+print(v_segments)
+print(j_segments)
+#+end_src
+
+#+RESULTS:
+#+begin_example
+  A DNAStringSet instance of length 147
+      width seq                                             names
+  [1]   326 GATACTGGAATTACCCAGACAC...ATCTCTGCACCAGCAGCCAAGA TRBV1*01_P
+  [2]   326 GATGCTGAAATCACCCAGAGCC...ATTTCTGCGCCAGCAGTGAGTC TRBV10-1*01_F
+  [3]   326 GATGCTGAAATCACCCAGAGCC...ATTTCTGCGCCAGCAGTGAGTC TRBV10-1*02_F
+  [4]   326 GATGCTGGAATCACCCAGAGCC...ATTTCTGCGCCAGCAGTGAGTC TRBV10-2*01_F
+  [5]   326 GATGCTGGAATCACCCAGAGCC...ATTTCTGCGCCAGCAGTGAGTC TRBV10-2*02_F
+  ...   ... ...
+[143]   324 GATACTGGAGTCTCCCAGAACC...GTATCTCTGTGCCAGCACGTTG TRBV7-9*06_(F)
+[144]   323 .........................TGTATCTCTGTGCCAGCAGCAG TRBV7-9*07_(F)
+[145]   325 GATTCTGGAGTCACACAAACCC...TATTTCTGTGCCAGCAGCGTAG TRBV9*01_F
+[146]   325 GATTCTGGAGTCACACAAACCC...TATTTCTGTGCCAGCAGCGTAG TRBV9*02_F
+[147]   321 GATTCTGGAGTCACACAAACCC...TTTGTATTTCTGTGCCAGCAGC TRBV9*03_(F)
+  A DNAStringSet instance of length 16
+     width seq                                              names
+ [1]    32 TGGGCGTCTGGGCGGAGGACTCCTGGTTCTGG                 TRBJ2-2P*01_ORF
+ [2]    31 TTTGGAGAGGGAAGTTGGCTCACTGTTGTAG                  TRBJ1-3*01_F
+ [3]    31 TTTGGTGATGGGACTCGACTCTCCATCCTAG                  TRBJ1-5*01_F
+ [4]    31 TTTGGCAGTGGAACCCAGCTCTCTGTCTTGG                  TRBJ1-4*01_F
+ [5]    31 TTCGGTTCGGGGACCAGGTTAACCGTTGTAG                  TRBJ1-2*01_F
+ ...   ... ...
+[12]    31 TTTGGCCCAGGCACCCGGCTGACAGTGCTCG                  TRBJ2-3*01_F
+[13]    31 TTCGGGCCAGGCACGCGGCTCCTGGTGCTCG                  TRBJ2-5*01_F
+[14]    31 TTCGGGCCAGGGACACGGCTCACCGTGCTAG                  TRBJ2-1*01_F
+[15]    31 TTCGGGCCGGGCACCAGGCTCACGGTCACAG                  TRBJ2-7*01_F
+[16]    31 GTCGGGCCGGGCACCAGGCTCACGGTCACAG                  TRBJ2-7*02_ORF
+#+end_example

generation.sh

@@ -1,21 +1,22 @@
 #!/bin/sh
 
 usage() {
-	echo "usage: generation.sh <number of sequences>"
+	echo "usage: generation.sh <number of sequences> <number of reads>"
 	exit 1
 }
 
-if [ $# != 1 ]; then
+if [ $# != 2 ]; then
 	usage
 fi
 
 sequences=$1
+number_of_reads=$2
 data_directory="data/"
+fastq=".fastq"
+filename="sequence"
 prefix="curesim_"
 
-Rscript src/repertoire.r "$sequences"
+Rscript src/repertoire.r "$sequences" "$number_of_reads" &&
-
+	CuReSim -f "$data_directory$filename$fastq" -o "$data_directory$prefix$filename$fastq"
-for file in "$data_directory"*.fastq; do
+Rscript src/alignment.r
-	file_name=$(echo "$file" | cut -d / -f 2)
+rm "$data_directory/log.txt"
-	java -jar tools/CuReSim.jar -f "$file" -o "$data_directory$prefix$file_name"
-done

shell.nix

@@ -2,14 +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
-    # Develoment tools
+    CuReSim
-    rPackages.languageserver
-    rPackages.lintr
   ];
 }

src/alignment.r

@@ -0,0 +1,101 @@
+library(Biostrings)
+library(parallel)
+
+parse_data <- function(file) {
+  reversed_sequences <- Biostrings::readQualityScaledDNAStringSet(file)
+  sequences <- Biostrings::reverseComplement(reversed_sequences)
+  vj_segments <- union(
+    readRDS("data/v_segments.rds"),
+    readRDS("data/j_segments_phe.rds")
+  )
+  return(list(sequences, vj_segments))
+}
+
+parse_metadata <- function(metadata) {
+  id_elements <- unlist(strsplit(metadata, split = " "))
+  v_identifier <- id_elements[2]
+  j_identifier <- id_elements[3]
+  return(list(v_id = v_identifier, j_id = j_identifier))
+}
+
+match_id_sequence <- function(names, vdj_segments, id) {
+  matches <- grep(names, pattern = id)
+  row <- matches[1]
+  return(as.character(vdj_segments[row]))
+}
+
+get_vj_sequence <- function(metadata, names, vdj_segments) {
+  identifiers <- parse_metadata(metadata)
+  v_sequence <- match_id_sequence(names, vdj_segments, id = identifiers["v_id"])
+  j_sequence <- match_id_sequence(names, vdj_segments, id = identifiers["j_id"])
+  return(list(v_seq = v_sequence, j_seq = j_sequence))
+}
+
+fetch_vj_sequences <- function(sequences, vdj_segments) {
+  vj_sequences <- sapply(names(sequences),
+    names(vdj_segments),
+    vdj_segments,
+    FUN = get_vj_sequence
+  )
+  results <- data.frame(t(vj_sequences))
+  return(results)
+}
+
+align_sequence <- function(sequence, vdj_segment) {
+  return(Biostrings::pairwiseAlignment(
+    subject = sequence,
+    pattern = vdj_segment,
+    type = "global-local",
+    gapOpening = 1
+  ))
+}
+
+handle_indels <- function(insertion, deletion, cys, alignment) {
+  ins_start <- sum(Biostrings::width(deletion[start(deletion) <= cys$start]))
+  ins_end <- sum(Biostrings::width(deletion[end(deletion) <= cys$end]))
+  shift_num <- c(0, cumsum(Biostrings::width(insertion))[-length(ins_start)])
+  shifted_ins <- IRanges::shift(insertion, shift_num)
+  gaps <- sum(width(shifted_ins[end(shifted_ins) < cys$start + ins_start])) +
+    nchar(stringr::str_extract(alignedSubject(alignment), "^-*"))
+  return(list("start" = ins_start - gaps, "end" = ins_end - gaps))
+}
+
+get_cys_coordinates <- function(alignment) {
+  cys <- list("start" = 310, "end" = 312)
+  insertion <- unlist(Biostrings::insertion(alignment))
+  deletion <- unlist(Biostrings::deletion(alignment))
+  delta_coordinates <- handle_indels(insertion, deletion, cys, alignment)
+  cys_start <- cys$start + delta_coordinates$start
+  cys_end <- cys$end + delta_coordinates$end
+  return(list("start" = cys_start, "end" = cys_end))
+}
+
+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,
+    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)),
+    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)
+}
+
+data <- parse_data(file = "data/curesim_sequence.fastq")
+hvr <- get_hvr_sequences(sequences = data[[1]], vdj_segments = data[[2]])
+Biostrings::writeXStringSet(hvr, "data/CuReSim-HVR.fastq", format = "fastq")

src/repertoire.r

@@ -1,55 +1,42 @@
 library(immuneSIM)
 library(Biostrings)
 
-generate_repertoires <- function(number_of_sequences) {
+generate_repertoire <- function(number_of_sequences) {
-  a_chain <- immuneSIM(
+  return(immuneSIM(
     number_of_seqs = number_of_sequences,
     species = "hs",
     receptor = "tr",
-    chain = "a",
+    chain = "b"
-    verbose = TRUE
+  ))
+}
+
+save_data <- function(data) {
+  Biostrings::writeXStringSet(data$sequence,
+    "data/sequence.fastq",
+    format = "fastq"
   )
-  b_chain <- immuneSIM(
+  Biostrings::writeXStringSet(data$junction, "data/HVR.fastq", format = "fastq")
-    number_of_seqs = number_of_sequences,
-    species = "hs",
-    receptor = "tr",
-    chain = "b",
-    verbose = TRUE
-  )
-  return(list("a_chain" = a_chain, "b_chain" = b_chain))
 }
 
-process_chain <- function(repertoire) {
+process_data <- function(data, reads) {
-  sequences <- as.character(repertoire$sequence)
+  dna_sequence <- Biostrings::DNAStringSet(data$sequence)
-  counts <- as.integer(repertoire$counts)
+  data$sequence <- Biostrings::reverseComplement(dna_sequence)
-  reads <- Biostrings::DNAStringSet(rep(sequences, counts))
+  names(data$sequence) <- paste(rownames(data), data$v_call, data$j_call, " ")
-  names(reads) <- seq_len(length(reads))
+  data$junction <- Biostrings::DNAStringSet(data$junction)
-  reverse_complement <- Biostrings::reverseComplement(reads)
+  names(data$junction) <- rownames(data)
-  return(reverse_complement)
+  amplified_data <- data[rep(seq_len(nrow(data)), reads), ]
-}
+  return(amplified_data)
-
-preprocess_data <- function(repertoires) {
-  filtered_repertoires <- lapply(repertoires, process_chain)
-  names(filtered_repertoires) <- names(repertoires)
-  return(filtered_repertoires)
-}
-
-save_data <- function(repertoires) {
-  for (chain in names(repertoires)) {
-    file_name <- paste("data/", chain, ".fastq", sep = "")
-    Biostrings::writeXStringSet(repertoires[[chain]], file_name, format = "fastq")
-  }
-}
-
-parse_cli_arguments <- function(args) {
-  if (length(args) != 1) {
-    stop("usage: repertoire.r <number of sequences>")
-  }
-  return(as.integer(args[1]))
 }
 
+parse_cli_arguments <- function() {
   args <- commandArgs(trailingOnly = TRUE)
-number_of_sequences <- parse_cli_arguments(args)
+  if (length(args) != 2) {
-sim_repertoire <- generate_repertoires(number_of_sequences)
+    stop("usage: repertoire.r <number of sequences> <number of reads>")
-processed_data <- preprocess_data(sim_repertoire)
+  }
-save_data(processed_data)
+  return(args)
+}
+
+args <- parse_cli_arguments()
+repertoire <- generate_repertoire(number_of_sequences = as.integer(args[1]))
+data <- process_data(data = repertoire, reads = args[2])
+save_data(data)