Use csi index in GenomicAlignmets (Rsamtools )

Question

[SSK]

Hi. I'm trying to do random sampling with R. When I run

library(GenomicAlignments)
library(rtracklayer)

## your bam file
bam.file <- 'your_bam_file.bam'

## is it paired end or not
paired <- TRUE

## read in the bam file
if(paired){
  bam.input <- readGAlignmentPairs(bam.file)
}else{
  bam.input <- readGAlignments(bam.file)
}

## select your number of subsets
num.samples <- 100

## subsample WITH replacement
bam.subsample <- sample(bam.input,num.samples,replace=TRUE)

## export new bam file
export(bam.subsample,BamFile('bam_subsampled.bam'))

I get the following error message

[E::hts_idx_push] Region 554468587..554468738 cannot be stored in a bai index. Try using a csi index with min_shift = 14, n_lvls >= 6
 value[[3L]](cond) でエラー: 'asBam' failed to build index
  file: bam_subsampled.bam
  SAM file: 'bam_subsampled.sam'
 追加情報:  警告メッセージ: 
 .make_GAlignmentPairs_from_GAlignments(gal, strandMode = strandMode,  で: 
    79328 alignments with ambiguous pairing were dumped.
    Use 'getDumpedAlignments()' to retrieve them from the dump environment.

I searched for csi option in Rsamtools or GenomicAlignmets, but I wasn't able to find that.

How can I solve this error ?? Thanks.

[SSK]

I got this code from C: Sampling random reads from BAM files with allowing same reads

Answer 1

[ATpoint]

Simply add index = FALSE to avoid building the index.

export(bam.subsample,BamFile('bam_subsampled.bam'), index = FALSE)

If you need the index you can create it with samtools at any time (samtools index -c bam_subsampled.bam)

Also asked for some details over at the BioC forum, see what they say. https://support.bioconductor.org/p/126603/

[SSK]

ATpoint It worked! Thanks!!