1 Commits

Author SHA1 Message Date
848a239f2c Convert org mode README to markdown 2021-05-05 12:38:13 +02:00
2 changed files with 7 additions and 12 deletions

View File

@@ -63,6 +63,6 @@ command invokes it:
The script will generate 2 files under the data directory:
|HVR.fastq | curesim-HVR.fastq |
|HVR.fastq |Contains the original CDR3 sequence |
|:----:|:-----:|
|Contains the original CDR3 sequence|Contains CDR3 after the read simulation, with sequencing errors |
|CuReSim-HVR.fastq | Contains CDR3 after the read simulation, with sequencing errors |

View File

@@ -34,11 +34,7 @@ parse_metadata <- function(metadata) {
#' @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]))
}
@@ -110,9 +106,8 @@ get_cys_coordinates <- function(alignment) {
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))
}
@@ -131,7 +126,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 +145,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")