R/make_snpsRNA.R

In joshstolz/brainstorm: Validating Sample Identiy Through RNA-seq and DNA Genotyping Correlation

#' make_snpsRNA
#'
#' Prep rna snps and create DNA & RNA snps (snpsGeno2) and RNA snps (snpsCalled)
#'
#' @param genotyped CollapsedVCF containing snp info for genotype data
#' @param called CollapsedVCF containing filtered snp info for rna data (generated by filtered_called)
#'
#' @return Named list containing snpsGeno2 and snpsCalled, both representing snp
#' arrays as integers.
#' @export
#'
#' @examples
#'
#' snpsCalled_filter <- filter_called(snpsCalled_VCF)
#' snpsRNA <- make_snpsRNA(snpsGeno_VCF, snpsCalled_filter)
#' snpsRNA$snpsGeno2[1:5, 1:5]
#' snpsRNA$snpsCalled[1:5, 1:5]
#' @importFrom SummarizedExperiment rowRanges
#' @importFrom VariantAnnotation alt geno
#' @importFrom Biostrings nchar
make_snpsRNA <- function(genotyped, called) {

    ## subset to intersection
    n <- intersect(rownames(called), names(genotyped))
    genotyped <- genotyped[n, ]
    called <- called[n, ]
    message("# matching snps: ", table(unlist(Biostrings::nchar(VariantAnnotation::alt(called)))))

    ## make snpGeno2
    snpsGeno2 <- VariantAnnotation::geno(genotyped)$GT
    snpsGeno2[snpsGeno2 == "./."] <- NA
    snpsGeno2[snpsGeno2 == "0|0"] <- 0
    snpsGeno2[snpsGeno2 == "1|0"] <- 1
    snpsGeno2[snpsGeno2 == "0|1"] <- 1
    snpsGeno2[snpsGeno2 == "1|1"] <- 2
    class(snpsGeno2) <- "numeric"

    ## flip
    table(rowRanges(called)$REF == rowRanges(genotyped)$REF |
        rowRanges(called)$REF == as.character(unlist(VariantAnnotation::alt(genotyped))))
    toFlip <- which(rowRanges(called)$REF != rowRanges(genotyped)$REF)
    snpsGeno2[toFlip, ] <- 2 - snpsGeno2[toFlip, ]

    ## VCF string to int - RNA style
    snpsCalled <- VariantAnnotation::geno(called)$GT
    snpsCalled[snpsCalled == "./."] <- 0
    snpsCalled[snpsCalled == "0/1"] <- 1
    snpsCalled[snpsCalled == "1/1"] <- 2
    class(snpsCalled) <- "numeric"

    snpsRNA <- list(snpsGeno2 = snpsGeno2, snpsCalled = snpsCalled)
    return(snpsRNA)
}