getSnpMats10X: Helper function to get coverage and allele count matrices...

In JEFworks/HoneyBADGER: HMM-integrated Bayesian approach for detecting CNV and LOH events using single-cell RNA-seq data

Description

Helper function to get coverage and allele count matrices given a set of putative heterozygous SNP positions specific for 10X data

Usage

getSnpMats10X(snps, bamFiles, indexFiles, barcodes, n.cores = 1,
  verbose = FALSE)

Arguments

snps	GenomicRanges object for positions of interest
bamFiles	list of bam file
indexFiles	list of bai index file
barcodes	Cell barcodes
n.cores	number of cores
verbose	Boolean of whether or not to print progress and info

Value

refCount reference allele count matrix for each cell and each position of interest altCount alternative allele count matrix for each cell and each position of interest cov total coverage count matrix for each cell and each position of interest

Examples

## Not run: 
# Get putative hets from ExAC
vcfFile <- "../data-raw/ExAC.r0.3.sites.vep.vcf.gz"
testRanges <- GRanges(chr, IRanges(start = 1, width=1000))
param = ScanVcfParam(which=testRanges)
vcf <- readVcf(vcfFile, "hg19", param=param)
## common snps by MAF
info <- info(vcf)
if(nrow(info)==0) {
    if(verbose) {
        print("ERROR no row in vcf")
    }
    return(NA)
}
maf <- info[, 'AF'] # AF is Integer allele frequency for each Alt allele
if(verbose) {
    print(paste0("Filtering to snps with maf > ", maft))
}
vi <- sapply(maf, function(x) any(x > maft))
if(verbose) {
    print(table(vi))
}
snps <- rowData(vcf)
snps <- snps[vi,]
## get rid of non single nucleotide changes
vi <- width(snps@elementMetadata$REF) == 1
snps <- snps[vi,]
## also gets rid of sites with multiple alt alleles though...hard to know which is in our patient
vi <- width(snps@elementMetadata$ALT@partitioning) == 1
snps <- snps[vi,]
## Get bams
files <- list.files(path = "../data-raw")
bamFiles <- files[grepl('.bam$', files)]
bamFiles <- paste0(path, bamFiles)
indexFiles <- files[grepl('.bai$', files)]
indexFiles <- paste0(path, indexFiles)
barcodes <- read.table("filtered_gene_bc_matrices/hg19/barcodes.tsv", header=FALSE)
results <- getSnpMats10X(snps, bamFiles, indexFiles, barcodes)

## End(Not run)

JEFworks/HoneyBADGER documentation built on July 24, 2021, 3:01 p.m.