Document the alignment script

README.org

@@ -1,6 +1,13 @@
 * locigenesis
 
-locigenesis is a tool that generates an immune repertoire and runs it through a sequence reader simulation tool, to generate sequencing errors.
+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
+
+- [[https://github.com/GreiffLab/immuneSIM/][immuneSIM]]: in silico generation of human and mouse BCR and TCR repertoires
+- [[http://www.pegase-biosciences.com/curesim-a-customized-read-simulator/][CuReSim]]: read simulator that mimics Ion Torrent sequencing
 
 ** Installation
 

docs/notebook.org

@@ -1,46 +0,0 @@
-#+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

nix/sources.json

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

shell.nix

@@ -21,16 +21,10 @@ let
     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"
+      makeWrapper ${jre}/bin/java $out/bin/CuReSim --add-flags "-jar $out/share/java/${name}/${name}.jar"
     '';
   };
 in mkShell {
-  buildInputs = [
+  buildInputs =
-    R
+    [ R rPackages.immuneSIM rPackages.Biostrings rPackages.stringr CuReSim ];
-    rPackages.immuneSIM
-    rPackages.Biostrings
-    rPackages.stringr
-    jdk
-    CuReSim
-  ];
 }

src/alignment.r

@@ -1,6 +1,10 @@
 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)
@@ -11,6 +15,10 @@ parse_data <- function(file) {
   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]
@@ -18,12 +26,24 @@ parse_metadata <- function(metadata) {
   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)
   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"])
@@ -31,6 +51,11 @@ get_vj_sequence <- function(metadata, names, vdj_segments) {
   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),
@@ -41,6 +66,11 @@ fetch_vj_sequences <- function(sequences, vdj_segments) {
   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,
@@ -50,6 +80,13 @@ align_sequence <- function(sequence, vdj_segment) {
   ))
 }
 
+#' 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]))
@@ -60,6 +97,10 @@ handle_indels <- function(insertion, deletion, cys, 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))
@@ -70,6 +111,12 @@ get_cys_coordinates <- function(alignment) {
   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,

src/repertoire.r

@@ -1,6 +1,10 @@
 library(immuneSIM)
 library(Biostrings)
 
+#' Generate the beta chain of a human T-cell receptor (TCR)
+#'
+#' @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,
@@ -10,6 +14,9 @@ generate_repertoire <- function(number_of_sequences) {
   ))
 }
 
+#' 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",
@@ -18,6 +25,11 @@ save_data <- function(data) {
   Biostrings::writeXStringSet(data$junction, "data/HVR.fastq", format = "fastq")
 }
 
+#' Applies the reverse complement and amplifies the number of sequences
+#'
+#' @param data A \code{data.frame} containing the TCR sequences and CDR3
+#' @param reads Number of times to amplify each sequence
+#' @return A \code{data.frame} with reverse complement sequences and VJ metadata
 process_data <- function(data, reads) {
   dna_sequence <- Biostrings::DNAStringSet(data$sequence)
   data$sequence <- Biostrings::reverseComplement(dna_sequence)
@@ -28,6 +40,9 @@ process_data <- function(data, reads) {
   return(amplified_data)
 }
 
+#' Checks the number of command line arguments and captures them
+#'
+#' @return A \code{vector} containing the command line arguments
 parse_cli_arguments <- function() {
   args <- commandArgs(trailingOnly = TRUE)
   if (length(args) != 2) {