Create a script that executes the pipeline

.gitignore

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

README.md

@@ -1,68 +0,0 @@
-# locigenesis
-
-locigenesis is a tool that generates a human T-cell receptor (TCR), runs
-it through a sequence reader simulation tool and extracts CDR3.
-
-The goal of this project is to generate both HVR sequences with and
-without sequencing errors, in order to create datasets for a Machine
-Learning algorithm.
-
-## Technologies
-
--   [immuneSIM](https://github.com/GreiffLab/immuneSIM/): in silico
-    generation of human and mouse BCR and TCR repertoires
--   [CuReSim](http://www.pegase-biosciences.com/curesim-a-customized-read-simulator/):
-    read simulator that mimics Ion Torrent sequencing
-
-## Installation
-
-This project uses [Nix](https://nixos.org/) to ensure reproducible
-builds.
-
-1.  Install Nix (compatible with MacOS, Linux and
-    [WSL](https://docs.microsoft.com/en-us/windows/wsl/about)):
-
-```bash
-curl -L https://nixos.org/nix/install | sh
-```
-
-2.  Clone the repository:
-
-```bash
-git clone https://git.coolneng.duckdns.org/coolneng/locigenesis
-```
-
-3.  Change the working directory to the project:
-
-```bash
-cd locigenesis
-```
-
-4.  Enter the nix-shell:
-
-```bash
-nix-shell
-```
-
-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:
-
-```bash
-./generation.sh <number of sequences> <number of reads>
-```
-
--   \<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  | curesim-HVR.fastq |
-|:----:|:-----:|
-|Contains the original CDR3 sequence|Contains CDR3 after the read simulation, with sequencing errors |

README.org

@@ -0,0 +1,3 @@
+* locigenesis
+
+locigenesis is a tool that generates an immune repertoire and runs it through a sequence reader simulation tool, to generate sequencing errors.

generation.sh

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

nix/sources.json

@@ -17,10 +17,10 @@
         "homepage": "",
         "owner": "NixOS",
         "repo": "nixpkgs",
-        "rev": "359e6542e1d41eb18df55c82bdb08bf738fae2cf",
+        "rev": "6f1ce38d0c0b1b25727d86637fd2f3baf7b0f1f6",
-        "sha256": "05v28njaas9l26ibc6vy6imvy7grbkli32bmv0n32x6x9cf68gf9",
+        "sha256": "16da722vqn96k1scls8mr8l909hl66r7y4ik6sad4ms3vmxbkbb3",
         "type": "tarball",
-        "url": "https://github.com/NixOS/nixpkgs/archive/359e6542e1d41eb18df55c82bdb08bf738fae2cf.tar.gz",
+        "url": "https://github.com/NixOS/nixpkgs/archive/6f1ce38d0c0b1b25727d86637fd2f3baf7b0f1f6.tar.gz",
         "url_template": "https://github.com/<owner>/<repo>/archive/<rev>.tar.gz"
     }
 }

shell.nix

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

src/alignment.r

@@ -1,153 +0,0 @@
-library(Biostrings)
-library(parallel)
-
-#' Import and process the TCR and VJ sequences
-#'
-#' @param file A file path with the sequences after applying a read simulator
-#' @return A \code{list} with the TCR sequences and VJ sequences
-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))
-}
-
-#' Extracts the VJ metadata from the sequences read identifier
-#'
-#' @param metadata The read identifier of a sequence
-#' @return A \code{list} with the V and J gene identifier
-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))
-}
-
-#' Fetches the sequence that matches the VJ gene identifier
-#'
-#' @param names The names of the VJ sequences
-#' @param vdj_segments A \code{DNAStringSet} containing the VJ sequences
-#' @param id The read identifier of a sequence
-#' @return A \code{character} containing the gene sequence
-match_id_sequence <- function(names, vdj_segments, id) {
-  matches <- grep(names, pattern = id)
-  if(id == "TRBJ2-2"){
-    row <- matches[2]
-  } else {
-    row <- matches[1]
-  }
-  return(as.character(vdj_segments[row]))
-}
-
-#' Gets the V and J sequences for a particular read identifier
-#'
-#' @param metadata The read identifier of a sequence
-#' @param names The names of the VJ sequences
-#' @param vdj_segments A \code{DNAStringSet} containing the VJ sequences
-#' @return A \code{list} with the V and J sequences
-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))
-}
-
-#' Obtains the VJ sequences for all the TCR sequences
-#'
-#' @param sequences A \code{QualityScaledDNAStringSet} with the TCR sequences
-#' @param vdj_segments A \code{DNAStringSet} containing the VJ sequences
-#' @return A \code{data.frame} with the V and J sequences
-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)
-}
-
-#' Perform a pairwise alignment of a sequence with the canonical V or J sequence
-#'
-#' @param sequence A \code{DNAString} containing the TCR sequences
-#' @param vdj_segment A \code{DNAString} containing the V or J sequence
-#' @return A \code{PairwiseAlignments}
-align_sequence <- function(sequence, vdj_segment) {
-  return(Biostrings::pairwiseAlignment(
-    subject = sequence,
-    pattern = vdj_segment,
-    type = "global-local",
-    gapOpening = 1
-  ))
-}
-
-#' Computes the coordinate shift of the Cysteine due to indels
-#'
-#' @param insertion An \code{IRanges} containing the insertions
-#' @param deletion An \code{IRanges} containing the deletions
-#' @param cys A \code{list} with the Cysteine coordinates
-#' @param alignment A \code{PairwiseAlignments}
-#' @return A \code{list} with the delta of the Cysteine coordinates
-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))
-}
-
-#' Find the coordinates of the first Cysteine of the HVR
-#'
-#' @param alignment A \code{PairwiseAlignments}
-#' @return A \code{list} with the Cysteine coordinates
-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)
-  read_start <- unlist(start(Biostrings::Views(alignment)))
-  cys_start <- cys$start + delta_coordinates$start + read_start - 1
-  cys_end <- cys$end + delta_coordinates$end + read_start
-  return(list("start" = cys_start, "end" = cys_end))
-}
-
-#' Delimit the hypervariable region (HVR) for each TCR sequence
-#'
-#' @param sequences A \code{QualityScaledDNAStringSet} with the TCR sequences
-#' @param vdj_segments A \code{DNAStringSet} containing the VJ sequences
-#' @param cores Number of cores to apply multiprocessing
-#' @return A \code{QualityScaledDNAStringSet} containing the HVR
-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) + 1)
-  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,57 +1,55 @@
 library(immuneSIM)
 library(Biostrings)
 
-#' Generate the beta chain of a human T-cell receptor (TCR)
+generate_repertoires <- function(number_of_sequences) {
-#'
+  a_chain <- immuneSIM(
-#' @param number_of_sequences Number of different sequences to generate
-#' @return A \code{data.frame} with the sequences, V and J genes and CDR3
-generate_repertoire <- function(number_of_sequences) {
-  return(immuneSIM(
     number_of_seqs = number_of_sequences,
     species = "hs",
     receptor = "tr",
-    chain = "b"
+    chain = "a",
-  ))
+    verbose = TRUE
-}
-
-#' Saves the sequences and CDR3 to FASTQ files
-#'
-#' @param data A \code{data.frame} with the preprocessed TCR sequences and CDR3
-save_data <- function(data) {
-  Biostrings::writeXStringSet(data$sequence,
-    "data/sequence.fastq",
-    format = "fastq"
   )
-  Biostrings::writeXStringSet(data$junction, "data/HVR.fastq", format = "fastq")
+  b_chain <- immuneSIM(
+    number_of_seqs = number_of_sequences,
+    species = "hs",
+    receptor = "tr",
+    chain = "b",
+    verbose = TRUE
+  )
+  return(list("a_chain" = a_chain, "b_chain" = b_chain))
 }
 
-#' Applies the reverse complement and amplifies the number of sequences
+process_chain <- function(repertoire) {
-#'
+  sequences <- as.character(repertoire$sequence)
-#' @param data A \code{data.frame} containing the TCR sequences and CDR3
+  counts <- as.integer(repertoire$counts)
-#' @param reads Number of times to amplify each sequence
+  reads <- Biostrings::DNAStringSet(rep(sequences, counts))
-#' @return A \code{data.frame} with reverse complement sequences and VJ metadata
+  names(reads) <- seq_len(length(reads))
-process_data <- function(data, reads) {
+  reverse_complement <- Biostrings::reverseComplement(reads)
-  dna_sequence <- Biostrings::DNAStringSet(data$sequence)
+  return(reverse_complement)
-  data$sequence <- Biostrings::reverseComplement(dna_sequence)
-  names(data$sequence) <- paste(rownames(data), data$v_call, data$j_call, " ")
-  data$junction <- Biostrings::DNAStringSet(data$junction)
-  names(data$junction) <- rownames(data)
-  amplified_data <- data[rep(seq_len(nrow(data)), reads), ]
-  return(amplified_data)
 }
 
-#' Checks the number of command line arguments and captures them
+preprocess_data <- function(repertoires) {
-#'
+  filtered_repertoires <- lapply(repertoires, process_chain)
-#' @return A \code{vector} containing the command line arguments
+  names(filtered_repertoires) <- names(repertoires)
-parse_cli_arguments <- function() {
+  return(filtered_repertoires)
-  args <- commandArgs(trailingOnly = TRUE)
+}
-  if (length(args) != 2) {
+
-    stop("usage: repertoire.r <number of sequences> <number of reads>")
+save_data <- function(repertoires) {
+  for (chain in names(repertoires)) {
+    file_name <- paste("data/", chain, ".fastq", sep = "")
+    Biostrings::writeXStringSet(repertoires[[chain]], file_name, format = "fastq")
   }
-  return(args)
 }
 
-args <- parse_cli_arguments()
+parse_cli_arguments <- function(args) {
-repertoire <- generate_repertoire(number_of_sequences = as.integer(args[1]))
+  if (length(args) != 1) {
-data <- process_data(data = repertoire, reads = args[2])
+    stop("usage: repertoire.r <number of sequences>")
-save_data(data)
+  }
+  return(as.integer(args[1]))
+}
+
+args <- commandArgs(trailingOnly = TRUE)
+number_of_sequences <- parse_cli_arguments(args)
+sim_repertoire <- generate_repertoires(number_of_sequences)
+processed_data <- preprocess_data(sim_repertoire)
+save_data(processed_data)