7

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

Devon Ryan
  • 19,602
  • 2
  • 29
  • 60
Biomagician
  • 2,459
  • 16
  • 30

1 Answers1

4

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)
llrs
  • 4,693
  • 1
  • 18
  • 42
Devon Ryan
  • 19,602
  • 2
  • 29
  • 60
  • 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? – Biomagician Jul 07 '17 at 15:18
  • Yup, that should work! – Devon Ryan Jul 07 '17 at 15:25