brentp / bigly Goto Github PK

Hi,
I am using mpileup for supporting reads from forward or reverse strand, but it is very slow. If bigly can track it as well, I would like to give it a try.

Thanks!
Tommy

Hi Brent,

I was trying to fetch records from a bam file with golang and bamat seems to be the most straight-forward solution I have ever seen. But when I did some tests, biogo/hts or bamat seemed to give unexpected results. Did I use biogo/hts or bamat correctly?

Here is my code:

First I tested with samtools:

# this command gives only one record
# samtools/htslib version: 1.11
./samtools view ./test.bam  chr1:100000-100200

Then I tried biogo/hts & bamat

import (
	"bufio"
	"github.com/biogo/hts/bam"
	"github.com/biogo/hts/sam"
	"os"
)

/*
The entire bamat source codes are copied here
*/

// this program gives 2000+ records
// the position of the first record starts from the beginning of bam
func main() {

	bamat, _ := New("test.bam")
	defer bamat.Close()

	iter, _ := bamat.Query("chr1", 100000, 100200)
	for iter.Next() {
		rec := iter.Record()
		println(rec.Name)
	}
}

Additional information:
system: centos7

1. test.bam is the alignment result of a paired-end sequencing experient.
2. test.bam is sorted by chromosome and position
3. test.bam.bai is built with samtools index command

I also printed the information of offsets and chunks:

// Query the BamAt with 0-base half-open interval.
func (b *BamAt) Query(chrom string, start int, end int) (*bam.Iterator, error) {
	if chrom == "" { // stdin
		return bam.NewIterator(b.Reader, nil)
	}
	ref := b.Refs[chrom]
	if end <= 0 {
		end = ref.Len() - 1
	}
	chunks, err := b.idx.Chunks(ref, start, end)
	for _, c := range chunks {
		fmt.Printf("Chunk, c.Begin.File: %v, c.Begin.Block: %v \n", c.Begin.File, c.Begin.Block)
		fmt.Printf("Chunk, c.End.File: %v, c.End.Block: %v \n", c.End.File, c.End.Block)
	}
	if err != nil {
		return nil, err
	}
	return bam.NewIterator(b.Reader, chunks)
}

output:

Chunk, c.Begin.File: 1157, c.Begin.Block: 0
Chunk, c.End.File: 152443, c.End.Block: 61430
Chunk, c.Begin.File: 305256, c.Begin.Block: 45799
Chunk, c.End.File: 369546, c.End.Block: 49973
Chunk, c.Begin.File: 4495641, c.Begin.Block: 6375
Chunk, c.End.File: 4495641, c.End.Block: 7201
Chunk, c.Begin.File: 4678143, c.Begin.Block: 22248
Chunk, c.End.File: 4678143, c.End.Block: 22521
Chunk, c.Begin.File: 5535256, c.Begin.Block: 32237
Chunk, c.End.File: 5535256, c.End.Block: 33167
Chunk, c.Begin.File: 8076117, c.Begin.Block: 22784
Chunk, c.End.File: 8076117, c.End.Block: 24440
Chunk, c.Begin.File: 8319574, c.Begin.Block: 22863
Chunk, c.End.File: 8319574, c.End.Block: 24508
Chunk, c.Begin.File: 54542702, c.Begin.Block: 55586
Chunk, c.End.File: 54542702, c.End.Block: 62429
Chunk, c.Begin.File: 59486419, c.Begin.Block: 10427
Chunk, c.End.File: 59486419, c.End.Block: 10975

This is my python code which gives the correct record:

import pysam

def query_bam(region, bamfile):

    samfile = pysam.AlignmentFile(bamfile, "rb")
    iterable = samfile.fetch(region=region)

    record_count = 0
    for record in iterable:
        print(record.query_name)
        record_count += 1

    samfile.close()
    return record_count

if __name__ == "__main__":
    bam_file = "test.bam"
    region = "chr1:100000-100200"
    total_records = query_bam(region, bam_file)
    print(f"Total records found: {total_records}")

Hi,

nice, I like the plotting function very much.

Any chance the plotting function could take an argument to save an image directly to disk instead of requiring x11 ?

Then it would be a nice way to script lots of images without manual work for example, and supply them to users.

Thanks,
Colin

• coordinate y-axes across samples

• add per-sample spinner.

• cancel previous request when new one is made

• make legend click work on all plots (currently just top plot)

• add "weird orientation line"

• bigly: split weird count into +/+, -/-, -/+

• don't show soft clips < 4.

• add simple repeats + micro satellites track

• make IGV bounds match those from the plot exactly.

• report insert size as number of std-devs above/below the mean calculated by sampling the bam.

• change bigly to accept a bam.Reader