lamortenera/bamsignals

R package to quickly obtain count vectors from indexed bam files

bamsignals's Issues

Exclusion of fragments with read pairs having supplementary alignments

Hello,
I am interested in counting fragments (RNA-seq) that overlap a certain region.
Also I would like to exclude fragments which have supplementary alignments reported for any of its read pairs (since those are putatively chimeric reads).

As far as I understood it, bamCoverage(bampath, gr, paired.end="extend", filteredFlag=2048, tlenFilter=NULL) could help me here.
Unfortunately, the reads that have supplementary alignments are still counted (I guess, since they are marked as primary alignment).
Would it be possible to exclude those reads from counting, e.g. by filtering them out based on their "Hidden Tag: SA"?
Thanks in advance.

bamProfile and bamCoverage only has one signal in output results

Hi, I ran bamsignals in R 4.0.3. When I ran bamProfile or bamCoverage the results returned only one object, which should be a list of objects with the number of the input Grange objects. I even ran the sample data in the package, which gave the same problem.
Could you please help me out?
Many thanks,
Xiaoyong Fu

Never mind. I restart my Rstudio, and it works fine. Sorry for the post.

Extending reads in bamCoverage?

For single end data, is there a way for reads to be extended? I.e. if you know that fragments are expected to be 200 bp, can you extend the reads to be that length when computing coverage? Such extensions are common practice when analyzing ChIP-seq data. Thanks!

deal with ranges not present in the bam file

As Wolfgang pointed out, it can happen that a GRanges object contains chrY but the bamfile doesn't, and this would create an error. Maybe a better behaviour is to issue a warning and return a zero signal...

bamsignals: Write vignette

Generate a vignette that demonstrates bamsignals use cases. Bam files can be generated with the artificial read generation routine from 48ed070

bamsignals: Implement S4 object

Implement S4 bamsignal objects. This enables us to publish the package on CRAN or bioconductor

consider merging the overlapAndPileup functions

The code is very similar and hard to keep synchronized

Accept two variants for BAM index file names

Minor change to also accept index files names like this

sample.bam + sample.bai
instead of just this
sample.bam + sample.bam.bai

A user of normR asked for this feature and it might be advisable to implement it in bamsignals. It seems to be related to Rhtslib::bam_index_load.

Error in testthat routine

Sometimes tests/testthat/test_methods.R fails with the following message (seems to be an error realted to the bam file generation):

> test()
Loading bamsignals
Testing bamsignals
CountSignals class and methods : ........................................
bamsignals methods : [W::sam_parse1] mapped query cannot have zero coordinate; treated as unmapped
[W::sam_parse1] mapped mate cannot have zero coordinate; treated as unmapped
[E::hts_idx_push] NO_COOR reads not in a single block at the end 1 -1
1

1. Error: bamCount function ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
BAM indexing file is not available for file /tmp/RtmpfM1L5N/file28be18bd1e53.bam

not opening bamfiles for which path contains symlinks

I have noticed that when I try running bamCoverage on bam paths that include symbolic links the function is not able to open them; when I give the full path it opens them just fine. The paths with the symbolic links work with functions from other packages just fine. It is easy to work around this by first converting the files to full paths, but I thought I'd bring it up anyways in case you weren't aware of this issue as it may or may not be a general issue

bamsignals: Write tests

With my recent commit ( fa399a4 ), I added tests/testthat.R which should implement some naive testing procedures.

We should refactor the code from tests/compile.R, tests/count.R, tests/pileup.R, tests/depth.R because they have already some use cases implemented that uses toy data already in the repository.

what happens with negative shift values?

Installation issues - bamsignals

Hi,

I have been trying to install bamsignals package. When I install through biocLite(), it is giving the following error.

Error in system.file("lib", package = "Rhtslib", mustWork = TRUE) : 
  no file found
Calls: <Anonymous> -> system.file
Execution halted

But I already installed Rhtslib. It would be great if you could help me fixing this.

Thank you.

build problems on Windows

Nate says that there are problems building bamsignals on Windows

bamsignals: Index bamfile if no index is present

As an enhancement, if a bampath points to a bam with no index, a bam index can be created from src/bamsignals.cpp.

What do you think?

bamProfile issue with large BAM files

Hi, I really like using this package to extract coverage information from BAM files in R.

When using the function "bamProfile", I noticed that for BAM files bigger than 2GB (approximately), I always get one of the following error messages:

[W::sam_hdr_read] bgzf_check_EOF: Invalid argument
[W::sam_hdr_read] bgzf_check_EOF: No such file or directory

The funtion works normally for BAM files up to 1.8GB in size.
I was wondering if this problem is known, if it can be solved or if I am doing anything wrong.

Thank you so much for the help!

Filter for SAMFLAG (samtools view -F analogue)

Filter out reads based on certain flags, e.g. optical duplicate with MAPQ 1024 (https://broadinstitute.github.io/picard/explain-flags.html).

For paired end, we require the flag first read in a properly mapped pair to be set (SAMFLAG=66).

Many people want to count only reads which are not marked as duplicates (SAMFLAG=1024). This is true for paired end as well as single end data.

In my opinion, we should allow to pass two arguments to bamCount(), bamProfile() & bamCoverage():

1. requiredFlag to select for certain reads with a set flag
2. filterFlag to exclude certain reads with a set flag

Best,

Condition Counting Fragments on TLEN in Paired End Bam

As a further enhancement, bamsignals.cpp::overlapAndPileup() should be able to exclude fragments where TLEN exceeds/deceeds a certain user specified value, e.g. counting only mono-nucleosomal fragments.

Accounting for fragment length

Hi , I find the document in
https://bioconductor.org/packages/release/bioc/vignettes/bamsignals/inst/doc/bamsignals.html .
My question is about shift process in bamCount to extract reads count from given intervals.
In my understanding, the shift process should be as follows:

1. Shift all reads in the bam based shift size 2. Counts the reads number in the intervals.

But I also guess that your shift process may be:

1. For each interval, shift the reads only in the interval. 2. Count the reads number in the intervals.

I am not sure which one is the right for the your process.

Thanks very much.

Best regards,
Yang