Filter reads comming from double strand sequences from a bam File

Filter putative double strand DNA from a strand specific RNA-seq using a window sliding across the genome.

Usage

filterDNA(
  file,
  destination,
  statFile,
  sequences,
  mapqFilter = 0,
  paired,
  yieldSize = 1e+06,
  winWidth = 1000L,
  winStep = 100L,
  readProp = 0.5,
  threshold = 0.7,
  pvalueThreshold = 0.05,
  useCoverage = FALSE,
  mustKeepRanges,
  getWin = FALSE,
  minCov = 0,
  maxCov = 0,
  errorRate = 0.01
)

Arguments

file: the input bam file to be filterd. Your bamfile should be sorted and have an index file located at the same path.
destination: the file path where the filtered output will be written
statFile: the file to write the summary of the results
sequences: the list of sequences to be filtered
mapqFilter: every read that has mapping quality below mapqFilter will be removed before any analysis. If missing, the entire bam file will be read.
paired: if TRUE then the input bamfile will be considered as paired-end reads. If missing, 100 thousands first reads will be inspected to test if the input bam file in paired-end or single-end.
yieldSize: by default is 1e6, i.e. the bam file is read by block of reads whose size is defined by this parameter. It is used to pass to same parameter of the scanBam function.
winWidth: the length of the sliding window, 1000 by default.
winStep: the step length to sliding the window, 100 by default.
readProp: a read is considered to be included in a window if at least readProp of it is in the window. Specified as a proportion. 0.5 by default.
threshold: the strand proportion threshold to test whether to keep a window or not. 0.7 by default
pvalueThreshold: the threshold for the p-value in the test of keeping windows. 0.05 by default
useCoverage: if TRUE, then the strand information in each window corresponds to the sum of coverage coming from positive/negative reads; and not the number of positive/negative reads as default.
mustKeepRanges: a GRanges object; all reads that map to those ranges will be kept regardless the strand proportion of the windows containing them.
getWin: if TRUE, the function will not only filter the bam file but also return a data frame containing the information of all windows of the original and filtered bam file.
minCov: if useCoverage=FALSE, every window that has less than minCov reads will be rejected regardless the strand proportion. If useCoverage=TRUE, every window has max coverage least than minCov will be rejected. 0 by default
maxCov: if useCoverage=FALSE, every window that has more than maxCov reads will be kept regardless the strand proportion. If useCoverage=TRUE, every window with max coverage more than maxCov will be kept. If 0 then it doesn't have effect on selecting window. 0 by default.
errorRate: the probability that an RNA read takes the false strand. 0.01 by default.

Value

if getWin is TRUE: a DataFrame object which could also be obtained by the function getStrandFromBamFile

Details

filterDNA reads a bam file containing strand specific RNA reads, and filter reads coming from putative double strand DNA. Using a window sliding across the genome, we calculate the positive/negative proportion of reads in each window. We then use logistic regression to estimate the strand proportion of reads in each window, and calculate the p-value when comparing that to a given threshold. Let \(\pi\) be the strand proportion of reads in a window.

Null hypothesis for positive window: \(\pi \le threshold\).

Null hypothesis for negative window: \(\pi \ge 1-threshold\).

Only windows with p-value <= pvalueThreshold are kept. For a kept positive window, each positive read in this window is kept with the probability (P-M)/P where P be the number of positive reads, and M be the number of negative reads. That is because those M negative reads are supposed to come from double-strand DNA, then there should be also M postive reads among the P positive reads come from double-strand DNA. In other words, there are only (P-M) positive reads come from RNA. Each negative read is kept with the probability equalling the rate that an RNA read of your sample has wrong strand, which is errorRate. Similar for kept negative windows.

Since each alignment can be belonged to several windows, then the probability of keeping an alignment is the maximum probability defined by all windows that contain it.

Examples

file <- system.file('extdata','s2.sorted.bam',package = 'strandCheckR')
filterDNA(file,sequences='10',destination='out.bam')
#> Testing paired end by checking the first 1e+05 reads of file /__w/_temp/Library/strandCheckR/extdata/s2.sorted.bam
#> Your bam file is single end
#> Summary will be written to out.stat
#> Read sequences 10