Hi,

I am using the ASpli R package to do a project on Chlamydomonas reinhardtii when I ran into a problem using one of the commands.

Here is my R session and everything I have done so far:

library(GenomicFeatures)
library(ASpli)
aTxDb <- loadDb(file="Cours/Stage/Cre.v5.5.genome.txdb")

features <- binGenome( aTxDb )
geneCoord <- featuresg( features )
binCoord <- featuresb( features )
junctionCoord <- featuresj( features )
binMetadata <- mcols( binCoord )
symbols <- data.frame( row.names = unique(genes( aTxDb )),
                       symbol = paste( 'This is symbol of gene:', unique(genes( aTxDb ) ) ))
features <- binGenome( aTxDb, geneSymbols = symbols )
features

This is what I get with "Features": 

Object of class ASpliFeatures 
Genes: GRangesList of length  17741 Access using featuresg(object) 
Bins: GRanges of length 320712 Access using featuresb(object) 
Junctions: GRanges of length 169749 Access using featuresj(object) 

bamFiles <- c( "bam1.bam","bam2.bam", "bam3.bam", "bam4.bam", "bam5.bam","bam6.bam","bam7.bam" )
targets <- data.frame( row.names = c("SRR685","SRR686","SRR698","SRR699", "SRR700", "SRR701","SRR702"),
                         bam = bamFiles,
                         genotype = c("ControlLight","ControlLight","mt-vegetative","mt+gamete","mt-gamete","zygote","post-germination") ,
                         stringsAsFactors = FALSE )
targets

library(GenomicAlignments)
bam <- loadBAM(targets)
counts <- readCounts (
  features,
  bam,
  targets,
  cores = 1,
  readLength = 51L,
  maxISize = 50000,
  minAnchor = 10 )
counts

This is what I get with counts:

Object of class ASpliCounts 
Gene counts: 17741 genes analysed. Access using countsg(object) 
Gene RD: 17741 genes analysed. Access using rdsg(object) 
Bin counts: 320712 bins analysed. Access using countsb(object) 
Bin RD: 320712 bins analysed. Access using rdsb(object) 
Junction counts: 0 junctions analysed. Access using countsj(object) 

As you can see, I somehow end up with 0 Junctions analysed and that prevents me from going any further. I have tried re-creating the bam files and loading them again, but I end up with the same result. I have made some changes to the maxISize, minAnchor but it doens't seem to change anything. Of course, when using the data from the examples of the package, the function works and junctions are indeed analyzed. So it makes me wonder what could be the problem with my files...

Is it coming from my Bam files ? The sam files were made with the Barracuda technique and then converted into bam files with samtools view. 

What bothers me is that I get Junctions: GRanges of length 169749 Access using featuresj(object) from "features", but then I get nothing analysed by the "readCounts" function.

If anyone has any idea of where the problem could be coming from and how to solve it, it would be very much appreciated !

Nicolas