Hi,

I have been trying to use Rsamtools’ applyPileups function to compare two files position-by-position. In order to test it out, I simply ran:
minBaseQuality <- 20
minMapQuality <- 30
minDepth <- 10
maxDepth <- 1000
testparam <- PileupParam(what="seq",
                         which=GRanges(“chr20", IRanges(1, 1000000)),
                         minBaseQuality=minBaseQuality,
                         minMapQuality=minMapQuality,
                         minDepth=minDepth,
                         maxDepth=maxDepth,
)
fl1 <- "10696X1_140408_SN141_0782_BC4J3NACXX_2_STP1N90A.bam"
fl2 <- "10696X4_140408_SN141_0782_BC4J3NACXX_2_STP1N90A.bam"
r3 = applyPileups(PileupFiles(c(fl1, fl2)), function(x) x, param=testparam)

My understanding is that this should result in a three-dimensional array with “ACTGN Counts” in the first dimension, files in the second and position in the third. Positions with overlapping reads in both files should thus be collapsed into a single line in the third dimension. However, selecting one of these positions shows that they are duplicated:

r3[[1]][["seq"]][ , , r3[[1]][["pos"]] == 135003]
yields:

, , 1

  10696X1_140408_SN141_0782_BC4J3NACXX_2_STP1N90A.bam 10696X4_140408_SN141_0782_BC4J3NACXX_2_STP1N90A.bam
A                                                   0                                                   0
C                                                   0                                                   0
G                                                  10                                                   0
T                                                   0                                                   0
N                                                   0                                                   0

, , 2

  10696X1_140408_SN141_0782_BC4J3NACXX_2_STP1N90A.bam 10696X4_140408_SN141_0782_BC4J3NACXX_2_STP1N90A.bam
A                                                   0                                                   0
C                                                   0                                                   0
G                                                   0                                                  13
T                                                   0                                                   0
N                                                   0                                                   0

Even though the position is the same, it is showing up twice. Each time, one of the files shows zeroes. This is not consistent with what happens if the files are identical (as in the example from the help docs).

For example,

r3 = applyPileups(PileupFiles(c(fl1, fl1)), function(x) x, param=testparam) #file 1 entered twice
 r3[[1]][["seq"]][ , , r3[[1]][["pos"]] == 135003]

yields:

  10696X1_140408_SN141_0782_BC4J3NACXX_2_STP1N90A.bam 10696X1_140408_SN141_0782_BC4J3NACXX_2_STP1N90A.bam
A                                                   0                                                   0
C                                                   0                                                   0
G                                                  10                                                  10
T                                                   0                                                   0
N                                                   0                                                   0

Is this the expected behavior? It seems like each position should only show up once in the output. Is there something I am missing?

Thanks,

Jonathon Hill
Postdoc
Yost Lab, University of Utah



	[[alternative HTML version deleted]]

