Remove redundant directories

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,56 @@
+* 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
+
+- [[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
+
+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 |

nix/sources.json

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

src/alignment.r

@@ -1,10 +1,6 @@
 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)
@@ -15,10 +11,6 @@ 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]
@@ -26,28 +18,12 @@ 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)
-  if(id == "TRBJ2-2"){
-    row <- matches[2]
-  } else {
-    row <- matches[1]
-  }
+  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"])
@@ -55,11 +31,6 @@ 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),
@@ -70,11 +41,6 @@ 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,
@@ -84,13 +50,6 @@ 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]))
@@ -101,27 +60,16 @@ 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))
   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
+  cys_start <- cys$start + delta_coordinates$start
+  cys_end <- cys$end + delta_coordinates$end
   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,
@@ -131,7 +79,7 @@ get_hvr_sequences <- function(sequences, vdj_segments, cores = detectCores()) {
   )
   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)
+  remaining <- Biostrings::subseq(sequences, start = unlist(cys_df$end))
   j_alignment <- parallel::mcmapply(remaining,
     df$j_seq,
     FUN = align_sequence,
@@ -150,4 +98,4 @@ get_hvr_sequences <- function(sequences, vdj_segments, cores = detectCores()) {
 
 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")
+Biostrings::writeXStringSet(hvr, "data/CuReSim-HVR.fastq", format = "fastq")

src/repertoire.r

@@ -1,10 +1,6 @@
 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,
@@ -14,9 +10,6 @@ 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",
@@ -25,11 +18,6 @@ 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)
@@ -40,9 +28,6 @@ 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) {