Qual not right after query() method

[suyanan]

irelate/parsers/bam.go

Line 166 in 1bf7c8f

func (b *BamQueryable) Query(region interfaces.IPosition) (interfaces.RelatableIterator, error) {

while this method running, first record from it.Next() then Read() method will write to interfaces.Relatable channel[0]index. Then second record also will write to interfaces.Relatable channel[1] index and update its Record.Qual, but which will unexpectedly update the Record.Qual value in nterfaces.Relatable channel[0]index.

Finally, the bam record's name and pos not coresponding to Qual.

[brentp]

Hi, can you explain the problem in more detail? I am not following.

[suyanan]

thanks.
my demo data is :

use vscode to debug it :

the 1st record is 0_1_E100005555L1C028R0333371821 and the 2nd 0_1_E100005555L1C040R0343126857.
at the second time in code's for circle, 1nd record's Qual which has stored in channel was updated to 2nd record'Qual.

[brentp]

I see. It seems that we must change here: https://github.com/brentp/irelate/blob/master/parsers/bam.go#L96

to:

			bam := Bam{Record: &(*rec), Chromosome: rec.Ref.Name(), related: nil}

it's been a while since I've done much go, but something like that should work.

[suyanan]

thanks.
i tried, its no work.
likely, channel's element and iterator's result use same pointer? before, i also use new() to get new address, no work.

[brentp]

Looks like it needs to go here as well: https://github.com/brentp/irelate/blob/master/parsers/bam.go#L215

[suyanan]

both changed in bam.go, also no work.
thanks a lot, i will get its life period to position problem.

[brentp]

If I understand correctly, this is a consequence of this change in biogo/hts
which re-uses a buffer instead of makeing a new one each time.

Perhaps it would be useful to have a sam.Record.Copy() method to do a deep copy of the record since it's now possible for the Qual to be changed in an existing record. Or perhaps there's another solution in place. @kortschak ?

[kortschak]

I'm surprised this doesn't show up in tests. @suyanan can you make a BAM available that has those two reads in it (preferably only those two reads and the necessary SQ headers)?

I'd rather not depend on a user making a clone, that is just asking for problems, but I can do use some more care when sharing a buffer.

[brentp]

Is there a test that keeps a read in memory and then reads the next from the bam? If not, I can write one in the next couple days.

[kortschak]

Yeah, I imagine not. The plan for the fix is to keep track of whether a new buffer was allocated, and only copy the Qual if it wasn't. This should fix the issue, but I'd like a test case.

[kortschak]

OK. I have a presumptive fix, but I need a test.

[kortschak]

@suyanan Would you please test your code with a replace line in your go.mod file like so:

replace github.com/biogo/hts => github.com/biogo/hts@v1.4.2-0.20210701111156-e5e9e6dcebd7

[brentp]

package main

import (
  "github.com/biogo/hts/bam"
  "fmt"
  "log"
  "os"
)

func main() {
  f, err := os.Open("x.bam")
  if err != nil {
    log.Fatalf("couldn't open file: %q", err)
  }
  defer f.Close()
  b, err := bam.NewReader(f, 1)
  if err != nil {
    log.Fatalf("couldn't read bam: %q", err)
  }
  defer b.Close()

  recA, err := b.Read()
  if err != nil {
    log.Fatalf("couldn't read first record: %q", err)
  }
  fmt.Fprintf(os.Stderr, "%q\n", recA.Qual)
  _, err = b.Read()
  if err != nil {
    log.Fatalf("couldn't read second record: %q", err)
  }
  fmt.Fprintf(os.Stderr, "%q\n", recA.Qual)
}

with bam (converted to sam of):

@HD	VN:1.5	GO:none	SO:coordinate
@SQ	SN:1	LN:249250621
A	1187	1	1000025	60	74M	=	1000064	113	AGCTGGGGGTCCACTGGGCTCAGGGAAGACCCCCTGCCAGGGAGACCCCAGGCGCCTGAATGGCCACGGGAAGG	BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
A	1107	1	1000064	60	74M	=	1000064	-113	AGCTGGGGGTCCACTGGGCTCAGGGAAGACCCCCTGCCAGGGAGACCCCAGGCGCCTGAATGGCCACGGGAAGG	DDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDD

[suyanan]

Thank you very much. I just read this, its no nessesary replying bam.
Follow you new changes, I update biogo/hts to v1.4.2 with Read() mothed updating, the irelate's Query() method and my similar method also worked.
Thanks a lot.