parklab / bamsnap Goto Github PK

Is it possible to show an empty coverage track.
I am comparing different genome capture sets, and one misses a specific region completely. Currently using the Trio example, the generated output png only shows one panel (of the bam file with the coverage). Hoping to show another panel with no coverage.

Thank you.

--disregard - figured it out.

my command was

bamsnap -bam ./NA12878.final.cram -pos chr11:12314 -out ./test.png

It worked well but suddenly it didn't work with error below.

Traceback (most recent call last):
File "/usr/lib64/python3.6/multiprocessing/process.py", line 258, in _bootstrap
self.run()
File "/usr/lib64/python3.6/multiprocessing/process.py", line 93, in run
self._target(*self._args, **self._kwargs)
File "/home/sylee/.local/lib/python3.6/site-packages/bamsnap/bamsnap.py", line 227, in run_process_drawplot_bamlist
refseq = rseq.get_refseq(pos1)
File "/home/sylee/.local/lib/python3.6/site-packages/bamsnap/bamsnap.py", line 535, in get_refseq
refseq = self.get_refseq_from_ucsc(pos1)
File "/home/sylee/.local/lib/python3.6/site-packages/bamsnap/bamsnap.py", line 575, in get_refseq_from_ucsc
refseq[gpos] = seq[i]
IndexError: string index out of range

how can I fix this problem?

Hi,
Your software seems pretty straightforward to use, but somehow I cannot make it work.
I have installed it and I can make it work on the test data from your repository. For instance, I can clone your repository, do cd bamsnap and from there run the command:

bamsnap -bam ./tests/data/test_SV1_softclipped_1.bam \
    -title "Clipped read" \
    -pos chr1:37775740 \
    -out ./out/test_SV1-7.png \
    -bamplot coverage read \
    -margin 100 \
    -no_target_line \
    -show_soft_clipped \
    -read_color_by interchrom \
    -save_image_only

And I get this:

However, when I try to do this with my own BAM file (a Drosophila RNAseq data replicate):

bamsnap -bam reads_indexed_and_sorted.bam \
    -title "my_analysis" \
    -pos 3L:747400-755492 \
    -draw bamplot \
    -bamplot read \
    -out ./my_plot.png \
    -read_group strand

I get the following errors:

Process proc 1:
Traceback (most recent call last):
  File "/Users/username/opt/anaconda3/lib/python3.7/multiprocessing/process.py", line 297, in _bootstrap self.run()
  File "/.../.../multiprocessing/process.py", line 99, in run self._target(*self._args, **self._kwargs)
  File "/.../.../site-packages/bamsnap/bamsnap.py", line 233, in run_process_drawplot_bamlist refseq = rseq.get_refseq(pos1)
  File "/.../.../site-packages/bamsnap/bamsnap.py", line 541, in get_refseq refseq = self.get_refseq_from_ucsc(pos1)
  File "/.../.../site-packages/bamsnap/bamsnap.py", line 576, in get_refseq_from_ucsc if line[0] != '<':
IndexError: string index out of range
2024-02-02 16:20:03,200 : [INFO] Total running time: 0.4 sec

I am not sure what to make of this, as the example above did not need a reference sequence, but just in case, I tried again providing the ref file, adding to the previous command -ref Drosophila_melanogaster.BDGP6.46.dna.toplevel.fa.
The FASTA file is from Ensemble (and the BAM files were produce by mapping the fastq files to the same version of the genome, and then indexing and sorting the BAM files with this reference)
But now I get this error:

Process proc 1:
Traceback (most recent call last):
  File "/.../.../multiprocessing/process.py", line 297, in _bootstrap self.run()
  File "/.../.../multiprocessing/process.py", line 99, in run self._target(*self._args, **self._kwargs)
  File "/.../.../site-packages/bamsnap/bamsnap.py", line 235, in run_process_drawplot_bamlist
    imagefname = bsplot.drawplot_bamlist(pos1, image_w, bamlist, xscale, refseq)
  File "/.../.../site-packages/bamsnap/bamsnap.py", line 473, in drawplot_bamlist
    ia_sub = self.get_bamplot_image(bam, pos1, image_w, xscale, refseq)
  File "/.../.../site-packages/bamsnap/bamsnap.py", line 357, in get_bamplot_image rset.calculate_readmap(is_strand_group=True)
  File "/.../.../site-packages/bamsnap/drawreadset.py", line 305, in calculate_readmap r = DrawRead(a, self.refseq)
  File "/.../.../site-packages/bamsnap/drawread.py", line 66, in __init__ self.set_read_variant()
  File "/.../.../site-packages/bamsnap/drawread.py", line 126, in set_read_variant refbase = self.refseq[gp].upper()
KeyError: 448823
2024-02-02 16:33:43,856 : [INFO] Total running time: 0.1 sec

I am not sure how to interpret this error message. There seems to be a problem with the reference sequence but I do not know what can that be. I have tried with the FASTA file of the individual 3L chromosome arm, but that does not help.

Is there any reason why my commands would not work?
Thanks so much in advance.

Hey @danielmsk

Is it possible to use bamsnap programmatically - I have a few use cases where I'd like to embed some plots.

Hi,
Thank you for your great work. ^_^.
Is there a parameter used to sort the aligntrack by base, just like igv?

Thanks a lot,
xiucz

command :
bamsnap -bam ./F1_dad.bam -ref new_GCA_024713975.2_ASM2471397v2_genomic.fna -ref_index_rebuild -pos CM045671.1:1-31859138

I got bam file for "olive flounder" species, and then names of the chromosomes are
"CM045671.1", "CM045672.1" ...
As I wanted to screen shot the alignment image for chromosome "CM045671.1",
I set the position like above.

And then, the error occured:
/home/jwshin0727/miniconda3/lib/python3.10/site-packages/pyfaidx-0.7.2.1-py3.10.egg/pyfaidx/init.py:523: RuntimeWarning: Index file /home/jwshin0727/CNU/Reference/Chinese/new_GCA_024713975.2_ASM2471397v2_genomic.fna.fai is older than FASTA file /home/jwshin0727/CNU/Reference/Chinese/new_GCA_024713975.2_ASM2471397v2_genomic.fna.
warnings.warn(
Process proc 1:
Traceback (most recent call last):
File "/home/jwshin0727/miniconda3/lib/python3.10/multiprocessing/process.py", line 314, in _bootstrap
self.run()
File "/home/jwshin0727/miniconda3/lib/python3.10/multiprocessing/process.py", line 108, in run
self._target(*self._args, **self._kwargs)
File "/home/jwshin0727/miniconda3/lib/python3.10/site-packages/bamsnap-0.2.19-py3.10.egg/bamsnap/bamsnap.py", line 233, in run_process_drawplot_bamlist
refseq = rseq.get_refseq(pos1)
File "/home/jwshin0727/miniconda3/lib/python3.10/site-packages/bamsnap-0.2.19-py3.10.egg/bamsnap/bamsnap.py", line 543, in get_refseq
refseq = self.get_refseq_from_localfasta(pos1)
File "/home/jwshin0727/miniconda3/lib/python3.10/site-packages/bamsnap-0.2.19-py3.10.egg/bamsnap/bamsnap.py", line 592, in get_refseq_from_localfasta
refseq[gpos+1] = seq[i]
IndexError: string index out of range
2023-05-18 23:02:10,954 : [INFO] Total running time: 0.0 sec

How can I solve the problem??

I'd like to run bamsnap with a custom reference (assembly output).
This works disabling the "gene" track with:

bamsnap -draw coordinates bamplot coverage base \ 
   -out test.png -bam aln.bam -pos NODE_1_length_31064_cov_4.851893:1000 -ref contigs.fasta

A warning could be useful to figure this out.

Hey Daniel,

I met this error when try to run bamsnap on a rna-seq bam file:

Process proc 1:
Traceback (most recent call last):
  File "/home/cc439/anaconda3/lib/python3.7/multiprocessing/process.py", line 297, in _bootstrap
    self.run()
  File "/home/cc439/anaconda3/lib/python3.7/multiprocessing/process.py", line 99, in run
    self._target(*self._args, **self._kwargs)
  File "/home/cc439/anaconda3/lib/python3.7/site-packages/bamsnap/bamsnap.py", line 436, in start_process_drawplot_bamlist
    imagefname = self.drawplot_bamlist(pos1, image_w, bamlist, xscale, refseq)
  File "/home/cc439/anaconda3/lib/python3.7/site-packages/bamsnap/bamsnap.py", line 360, in drawplot_bamlist
    ia_sub = self.get_bamplot_image(bam, pos1, image_w, xscale, refseq)
  File "/home/cc439/anaconda3/lib/python3.7/site-packages/bamsnap/bamsnap.py", line 257, in get_bamplot_image
    rset.calculate_readmap(is_strand_group=True)
  File "/home/cc439/anaconda3/lib/python3.7/site-packages/bamsnap/drawreadset.py", line 301, in calculate_readmap
    r = DrawRead(a, self.refseq)
  File "/home/cc439/anaconda3/lib/python3.7/site-packages/bamsnap/drawread.py", line 66, in __init__
    self.set_read_variant()
  File "/home/cc439/anaconda3/lib/python3.7/site-packages/bamsnap/drawread.py", line 126, in set_read_variant
    refbase = self.refseq[gp].upper()
KeyError: 32793279

Any suggestions on this? Thank you!

bamsnap runs into error when -process was used, although it still successfully writes the output file. WIthout -process, the same command runs successfully.

Traceback (most recent call last):
  File "/.conda/envs/bamsnap/bin/bamsnap", line 10, in <module>
    sys.exit(cli())
  File "/.conda/envs/bamsnap/lib/python3.9/site-packages/bamsnap/__init__.py", line 12, in cli
    bs.run()
  File "/.conda/envs/bamsnap/lib/python3.9/site-packages/bamsnap/bamsnap.py", line 211, in run
    self.start_process_drawplot(image_w, self.bamlist)
  File "/.conda/envs/bamsnap/lib/python3.9/site-packages/bamsnap/bamsnap.py", line 83, in start_process_drawplot
    image_w, bamlist, self.split_poslist[tno], self.opt, self.is_single_image_out), name='proc ' + str(tno+1))
KeyError: 1
2021-08-09 13:57:15,362 : [INFO] (proc 1) Saved test.png
2021-08-09 13:57:15,362 : [INFO] (proc 1) Saved test.png : 3.18431 sec

Versions used:

• python - 3.9.4
• bamsnap - bamsnap ver0.2.19 (2021-03-02)
• OS - Red Hat Enterprise Linux Server release 7.8 (Maipo)
• Conda - conda 4.10.1

Thanks for this useful tool. Straightforward to use.

Some fasta files contain lowercase characters (e.g. 'a' instead of 'A') and these characters cause the rendering to stop. The errors are due to lookups of color hex codes for the lowercase characters.

    dr.line([(x, y11), (x, y21)], fill=COLOR[base], width=self.xscale.base_width)
KeyError: 'a'

I looked into the code and can submit a pull request with a fix.

Hi,

first of all, thanks for this very handy software!

I'm using bamsnap 2.7 and I was looking for a way to disable writing the .txt and .log output files. I tried the -save_image_only option, but if that's supposed to prevent those additional files from being written it does not seem to be working for me.

This is how I'm invoking the software:

  bamsnap \
    -bam ${input.bam} \
    -ref ${params.genome} \
    -refversion ${params.refversion} \
    -title ${name} \
    -out ${outfile}.png \
    -pos ${chr}:${start}-${end} \
    -save_image_only \
    -silence

Is it possible to accomplish this with version 2.7?

I was eager to try your tool and the examples in the repository work well.
I tried directly plotting from your examples as well as re-alignment (minimap2) on hg38 with subsequent plotting.
In that case both without and with provided fasta sequence worked.

Whenever I try though to use a non-human reference genome, I get an error :

2021-01-18 10:51:57,140 : [INFO] Total running time: 0.0 sec
Process proc 1:
Traceback (most recent call last):
  File "/usr/local/lib/python3.7/multiprocessing/process.py", line 297, in _bootstrap
    self.run()
  File "/usr/local/lib/python3.7/multiprocessing/process.py", line 99, in run
    self._target(*self._args, **self._kwargs)
  File "/usr/local/lib/python3.7/site-packages/bamsnap/bamsnap.py", line 227, in run_process_drawplot_bamlist
    refseq = rseq.get_refseq(pos1)
  File "/usr/local/lib/python3.7/site-packages/bamsnap/bamsnap.py", line 537, in get_refseq
    refseq = self.get_refseq_from_localfasta(pos1)
  File "/usr/local/lib/python3.7/site-packages/bamsnap/bamsnap.py", line 586, in get_refseq_from_localfasta
    refseq[gpos+1] = seq[i]
IndexError: string index out of range

In order to assure that is should work I took really headers such as ">sequence1",">sequence2" and so on and verified that reads mapped actually on provided positions but it always gave the same error.
Any hints how to fix that, or plans to add to your test repository some user reference sequences to make it easier to understand what is the proper syntax in that case ?

I had this issue with basic script

bamsnap -bam 15L136Fb6001F9.s.md.ir.bam -pos 1:22312455 -out test.png

Process proc 1: Traceback (most recent call last): File "/home/kds/.conda/envs/py311/lib/python3.10/multiprocessing/process.py", line 314, in _bootstrap self.run() File "/home/kds/.conda/envs/py311/lib/python3.10/multiprocessing/process.py", line 108, in run self._target(*self._args, **self._kwargs) File "/home/kds/.conda/envs/py311/lib/python3.10/site-packages/bamsnap/bamsnap.py", line 235, in run_process_drawplot_bamlist imagefname = bsplot.drawplot_bamlist(pos1, image_w, bamlist, xscale, refseq) File "/home/kds/.conda/envs/py311/lib/python3.10/site-packages/bamsnap/bamsnap.py", line 473, in drawplot_bamlist ia_sub = self.get_bamplot_image(bam, pos1, image_w, xscale, refseq) File "/home/kds/.conda/envs/py311/lib/python3.10/site-packages/bamsnap/bamsnap.py", line 377, in get_bamplot_image ia_sub = covplot.get_image(image_w, self.opt['coverage_height'], self.opt['coverage_bgcolor']) File "/home/kds/.conda/envs/py311/lib/python3.10/site-packages/bamsnap/drawreadset.py", line 26, in get_image self.draw_coverage(dr, w, h) File "/home/kds/.conda/envs/py311/lib/python3.10/site-packages/bamsnap/drawreadset.py", line 61, in draw_coverage dr.rectangle([(x1, y1), (x2, y2)], fill=getrgb(self.coverage_color), outline=getrgb(self.coverage_color, 15), width=1) File "/home/kds/.conda/envs/py311/lib/python3.10/site-packages/PIL/ImageDraw.py", line 294, in rectangle self.draw.draw_rectangle(xy, fill, 1) ValueError: y1 must be greater than or equal to y0

Anyone had the same problem?

How can I solve this?

Thank you so much for this great tool. I recently replace my visualization tools with bamsnap on our pipeline. However, today I just figured out that bamsnap didn't support mm10. Is it possible to support it?

I was wondering if bamsnap supports PacBio CCS reads. When I use bamsnap with my CCS bam files, it produces the following error:

Traceback (most recent call last):
  File "/Users/sl17/Languages/python3/lib/python3.8/multiprocessing/process.py", line 315, in _bootstrap
    self.run()
  File "/Users/sl17/Languages/python3/lib/python3.8/multiprocessing/process.py", line 108, in run
    self._target(*self._args, **self._kwargs)
  File "/Users/sl17/Languages/python3/lib/python3.8/site-packages/bamsnap/bamsnap.py", line 235, in run_process_drawplot_bamlist
    imagefname = bsplot.drawplot_bamlist(pos1, image_w, bamlist, xscale, refseq)
  File "/Users/sl17/Languages/python3/lib/python3.8/site-packages/bamsnap/bamsnap.py", line 473, in drawplot_bamlist
    ia_sub = self.get_bamplot_image(bam, pos1, image_w, xscale, refseq)
  File "/Users/sl17/Languages/python3/lib/python3.8/site-packages/bamsnap/bamsnap.py", line 383, in get_bamplot_image
    ia_sub = rset.get_image(
  File "/Users/sl17/Languages/python3/lib/python3.8/site-packages/bamsnap/drawreadset.py", line 350, in get_image
    self.draw_read(dr, group, readcolor, readcolorby)
  File "/Users/sl17/Languages/python3/lib/python3.8/site-packages/bamsnap/drawreadset.py", line 357, in draw_read
    r.yidx = self.get_yidx(r, group)
  File "/Users/sl17/Languages/python3/lib/python3.8/site-packages/bamsnap/drawreadset.py", line 199, in get_yidx
    self.readmap[gpos][group] = []
KeyError: 2725160
2021-08-05 14:22:24,744 : [INFO] Total running time: 3.1 sec

The command that I have used is:

bamsnap -bam input.bam -pos 20:2740569  -refversion hg19  -out chr20_2740569

Is there a limit on the size of genomic region that can be plotted using bamsnap? My attempts to plot 1.5MB region failed with this error:

Process proc 1:
Traceback (most recent call last):
  File ".conda/envs/bamsnap/lib/python3.9/multiprocessing/process.py", line 315, in _bootstrap
    self.run()
  File ".conda/envs/bamsnap/lib/python3.9/multiprocessing/process.py", line 108, in run
    self._target(*self._args, **self._kwargs)
  File ".conda/envs/bamsnap/lib/python3.9/site-packages/bamsnap/bamsnap.py", line 235, in run_process_drawplot_bamlist
    imagefname = bsplot.drawplot_bamlist(pos1, image_w, bamlist, xscale, refseq)
  File ".conda/envs/bamsnap/lib/python3.9/site-packages/bamsnap/bamsnap.py", line 467, in drawplot_bamlist
    ia = self.append_coordinates_image(ia, pos1, image_w, xscale)
  File ".conda/envs/bamsnap/lib/python3.9/site-packages/bamsnap/bamsnap.py", line 430, in append_coordinates_image
    ia_sub = coord.get_image()
  File ".conda/envs/bamsnap/lib/python3.9/site-packages/bamsnap/coordinates.py", line 106, in get_image
    self.draw(dr)
  File ".conda/envs/bamsnap/lib/python3.9/site-packages/bamsnap/coordinates.py", line 62, in draw
    self.cal_axis()
  File ".conda/envs/bamsnap/lib/python3.9/site-packages/bamsnap/coordinates.py", line 52, in cal_axis
    unitlen = 10 ** (len(str(axis_base_unit)) - 1)
UnboundLocalError: local variable 'axis_base_unit' referenced before assignment

Error is occuring due to conditional used with hardcoded list here:

Awesome tool. Could you add cram support? Having to go through a cram > bam conversion step isn't ideal.

Thanks!

Hi! This is a very nice tool that I just learnt about. Looking in the docs, this is offered as a command line tool, which is nice. But, given that this tool is written in python, is there a library API one can use to make and configure these plots from within python or jupyter notebooks?

Hello, is it possible to specify a custom annotation file instead of relying on the shipped ones? Because, the gene track plotting runs into error when using custom references, or references other than the ones for which the GTFs are shipped with this package. For example, there is no way to reliably use BAMs generated using T2T reference genome CHMv13.

Hi,
I am trying to use bamsnap with long read data from Nanopore aligned reads using minimap. I do not manage to get any output and the exit error is a number that I cannot debug (when I Google "KeyError: 28747748" it does not return any meaningful results). Here is the code that I am using and the output:

bamsnap -bam ../aln_bam/dRNA_aln.sorted.bam -pos chr16:28932034 -out ./test.png

2021-01-26 23:07:53,894 : [INFO] Total running time: 0.0 sec Process proc 1: Traceback (most recent call last): File "/home/mcortesl/.conda/envs/myenv/lib/python3.7/multiprocessing/process.py", line 297, in _bootstrap self.run() File "/home/mcortesl/.conda/envs/myenv/lib/python3.7/multiprocessing/process.py", line 99, in run self._target(*self._args, **self._kwargs) File "/home/mcortesl/.conda/envs/myenv/lib/python3.7/site-packages/bamsnap/bamsnap.py", line 229, in run_process_drawplot_bamlist imagefname = bsplot.drawplot_bamlist(pos1, image_w, bamlist, xscale, refseq) File "/home/mcortesl/.conda/envs/myenv/lib/python3.7/site-packages/bamsnap/bamsnap.py", line 467, in drawplot_bamlist ia_sub = self.get_bamplot_image(bam, pos1, image_w, xscale, refseq) File "/home/mcortesl/.conda/envs/myenv/lib/python3.7/site-packages/bamsnap/bamsnap.py", line 351, in get_bamplot_image rset.calculate_readmap(is_strand_group=True) File "/home/mcortesl/.conda/envs/myenv/lib/python3.7/site-packages/bamsnap/drawreadset.py", line 305, in calculate_readmap r = DrawRead(a, self.refseq) File "/home/mcortesl/.conda/envs/myenv/lib/python3.7/site-packages/bamsnap/drawread.py", line 66, in __init__ self.set_read_variant() File "/home/mcortesl/.conda/envs/myenv/lib/python3.7/site-packages/bamsnap/drawread.py", line 126, in set_read_variant refbase = self.refseq[gp].upper() KeyError: 28747748

Could this be a problem because of the formatting of the minimap2 aligned file?

When using the -process option with multiple threads, no html output or zip is created.

For example:

bamsnap -bamlist list_bam.txt -bed batch_sites_50.bed \
-ref ../angsd_calling/resources/ref_genome/mgal_01.fa \
-out ./batch_sites -draw coordinates bamplot coverage base \
-separated_bam -process 4 -zipout

It would be great if it was possible to sort reads by tag and color then by tag (similar to IGV).

I have been struggling with KeyErrors in bamsnap

Process proc 1: Traceback (most recent call last): File "/home/jaesoon/miniconda3/envs/gunhee/lib/python3.6/multiprocessing/process.py", line 258, in _bootstrap self.run() File "/home/jaesoon/miniconda3/envs/gunhee/lib/python3.6/multiprocessing/process.py", line 93, in run self._target(*self._args, **self._kwargs) File "/home/jaesoon/miniconda3/envs/gunhee/lib/python3.6/site-packages/bamsnap/bamsnap.py", line 235, in run_process_drawplot_bamlist imagefname = bsplot.drawplot_bamlist(pos1, image_w, bamlist, xscale, refseq) File "/home/jaesoon/miniconda3/envs/gunhee/lib/python3.6/site-packages/bamsnap/bamsnap.py", line 473, in drawplot_bamlist ia_sub = self.get_bamplot_image(bam, pos1, image_w, xscale, refseq) File "/home/jaesoon/miniconda3/envs/gunhee/lib/python3.6/site-packages/bamsnap/bamsnap.py", line 377, in get_bamplot_image ia_sub = covplot.get_image(image_w, self.opt['coverage_height'], self.opt['coverage_bgcolor']) File "/home/jaesoon/miniconda3/envs/gunhee/lib/python3.6/site-packages/bamsnap/drawreadset.py", line 26, in get_image self.draw_coverage(dr, w, h) File "/home/jaesoon/miniconda3/envs/gunhee/lib/python3.6/site-packages/bamsnap/drawreadset.py", line 81, in draw_coverage dr.line([(x, y11), (x, y21)], fill=COLOR[base], width=self.xscale.base_width) KeyError: **'a'**

and I assume that the KeyError is because the reference genome contains lower case letters.

When I tried with other species with reference genomes without any lower case letters, it worked

I would be grateful if you could consider this opinion.

Thank you

HTML created for multi-sample snaps can't load the snaps. The template is not updated and expects one snap for each sample while the program creates a single snap with all the samples. As a result, the HTML can't load pictures properly.

It would be fantastic to have bamsnap available via bioconda!

all the best,
Peter