-
Notifications
You must be signed in to change notification settings - Fork 0
/
main3.R
68 lines (50 loc) · 1.25 KB
/
main3.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
library(Rsamtools)
library(GenomicAlignments)
gr <- GRanges("chr14",IRanges(104029449,104029449),"*")
pathToFiles <- "data"
files <- list.files(pathToFiles, full.names=TRUE)
###########
#first version
###########
#part 1 - import to a galignments object
which <- gr
what <- scanBamWhat()
flag <- scanBamFlag(isUnmappedQuery = FALSE)
param <- ScanBamParam(flag = flag, which = which, what = what)
bf <- BamFile(files[1])
open(bf)
ga <- readGAlignments(bf, param = param)
close(bf)
#part 2 - summarize allele counts
my_IGPOI <- gr
seqlevels(ga) <- seqlevels(my_IGPOI)
qseq <- mcols(ga)$seq
nuclpiles <- pileLettersAt(qseq, seqnames(ga), start(ga), cigar(ga), my_IGPOI)
nstr <- table(strsplit(as.character(nuclpiles), ""))
nstr
###########
#second version
###########
countF <-
function(x){
x[["seq"]][-5,,1]
}
which <- gr
flag <- scanBamFlag(isUnmappedQuery = FALSE)
p1 <- ApplyPileupsParam(flag=flag,
which=which,
minBaseQuality = 0L,
what="seq",
yieldBy = "position",
yieldAll=TRUE,
maxDepth=.Machine$integer.max,
)
pf <- PileupFiles(bf)
res <- applyPileups(pf, countF, param=p1)
###########
# investigate
###########
#why do these indicate different read depth?
length(ga)
nstr
res