coolneng
/
locigenesis
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Isolate HVR sequence and save it to a file

This commit is contained in:
coolneng 2021-05-03 21:15:40 +02:00
parent c250c139dd
commit 5afe040592
Signed by: coolneng
GPG Key ID: 9893DA236405AF57
2 changed files with 38 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
shell.nix
View File

@ -7,6 +7,7 @@ mkShell {
R R
rPackages.immuneSIM rPackages.immuneSIM
rPackages.Biostrings rPackages.Biostrings
rPackages.stringr
jdk jdk
# Development tools # Development tools
rPackages.languageserver rPackages.languageserver

49
src/alignment.r
View File

@ -50,26 +50,51 @@ align_sequence <- function(sequence, vdj_segment) {
)) ))
} }
locate_cysteine <- function(sequence, v_segment) { handle_indels <- function(insertion, deletion, cys, alignment) {
codons <- Biostrings::DNAStringSet(c("TGT", "TGC")) ins_start <- sum(Biostrings::width(deletion[start(deletion) <= cys$start]))
matches <- Biostrings::matchPDict( ins_end <- sum(Biostrings::width(deletion[end(deletion) <= cys$end]))
pdict = codons, shift_num <- c(0, cumsum(Biostrings::width(insertion))[-length(ins_start)])
subject = DNAString(toString(v_segment)) shifted_ins <- IRanges::shift(insertion, shift_num)
) gaps <- sum(width(shifted_ins[end(shifted_ins) < cys$start + ins_start])) +
position <- as.data.frame(intersect(matches[[1]], matches[[2]])) nchar(stringr::str_extract(alignedSubject(alignment), "^-*"))
return(position) return(list("start" = ins_start - gaps, "end" = ins_end - gaps))
} }
# TODO Extract CDR3 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))
}
# TODO Refactor this mess
get_hvr_sequences <- function(sequences, vdj_segments) { get_hvr_sequences <- function(sequences, vdj_segments) {
df <- fetch_vj_sequences(sequences, vdj_segments) df <- fetch_vj_sequences(sequences, vdj_segments)
v_alignment <- parallel::mcmapply(sequences, df$v_seq, FUN = align_sequence) v_alignment <- parallel::mcmapply(sequences, df$v_seq, FUN = align_sequence)
hvr_start <- parallel::mcmapply(sequences, v_alignment, FUN = locate_cysteine) j_alignment <- parallel::mcmapply(sequences, df$j_seq, FUN = align_sequence)
hvr_start_df <- as.data.frame(t(hvr_start)) cys_coordinates <- parallel::mclapply(v_alignment, FUN = get_cys_coordinates)
cys_df <- as.data.frame(do.call(rbind, cys_coordinates))
j_start <- parallel::mclapply(
j_alignment,
function(x) start(Biostrings::Views(x))
)
hvr <- Biostrings::subseq(sequences,
start = unlist(cys_df$start),
end = unlist(j_start) + 2
)
return(hvr)
}
save_data <- function(data) {
Biostrings::writeXStringSet(data, "data/CuReSim-HVR.fastq", format = "fastq")
} }
data <- parse_data(file = "data/curesim_sequence.fastq") data <- parse_data(file = "data/curesim_sequence.fastq")
hvr_sequences <- get_hvr_sequences( hvr <- get_hvr_sequences(
sequences = data[[1]], sequences = data[[1]],
vdj_segments = data[[2]] vdj_segments = data[[2]]
) )
Biostrings::writeXStringSet(hvr, "data/CuReSim-HVR.fastq", format = "fastq")