7
$\begingroup$

I have this read in my BAM file. It maps on chromosome 1.

I open this BAM file in IGV, and I can see the alignment on chromosome 1.

But when I open this file in R with Rsamtools:

bamContigsCel <- Rsamtools::scanBam('output/alignment/pacbio/bwa/ref     /bristolAssemblySorted.bam', param = Rsamtools::ScanBamParam(what = Rsamtools::scanBamWhat(), flag = Rsamtools::scanBamFlag(isMinusStrand = FALSE), tag = bamTags))[[1]]

I then check if the read maps to chromosome 1 in my R object but I cannot find it.

bamContigsCel$rname[bamContigsCel$qname == '000000F|arrow']
[1] II II IV
Levels: I II III IV MtDNA V X

But if I look at the BAM file, it's there. Why isn't Rsamtools importing my read into R?

$ samtools view bristolAssemblySorted.bam | grep -n '000000F' | head  -c80
000000F|arrow   2064    I   336331  7   2926310H260M1774267H    *   0 0   GAAGCTGTCTAAACTTTGGC

Cross-posted on biostars

$\endgroup$
1
  • $\begingroup$ Careful, your code is using an almost certainly wrong path to the BAM file. $\endgroup$ Jul 7, 2017 at 12:51

1 Answer 1

4
$\begingroup$

Note the flag; that read is mapped in a reverse-complemented manner, so isMinusStrand = FALSE is filtering it out. I tested this by making a BAM file with only that read:

@SQ SN:I    LN:1000000
000000F|arrow   2064    I   336331  7   2926310H260M1774267H    *   00  *   *

Then in R:

> library(Rsamtools)
> blah = scanBam("foo.bam", param=ScanBamParam(what=scanBamWhat(), flag=scanBamFlag(isMinusStrand = TRUE)))[[1]]
> blah$qname
[1] "000000F|arrow"
> blah = scanBam("foo.bam", param=ScanBamParam(what=scanBamWhat(), flag=scanBamFlag(isMinusStrand = FALSE)))[[1]]
> blah$qname
character(0)
$\endgroup$
2
  • $\begingroup$ Thanks! I don't know why I was filtering out reverse-complemented alignments! Can I just use bam <- Rsamtools::scanBam('output/alignment/pacbio/bwa/ref/bristolAssemblySorted.bam')[[1]] to get all the information in the BAM file? $\endgroup$ Jul 7, 2017 at 15:18
  • $\begingroup$ Yup, that should work! $\endgroup$
    – Devon Ryan
    Jul 7, 2017 at 15:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.