From 66ee508f947eea7417a38c226dd8fac30e9dbd9a Mon Sep 17 00:00:00 2001 From: Benjamin Schuster-Boeckler Date: Tue, 5 May 2026 19:27:11 +0100 Subject: [PATCH 1/4] first test of CRAM support --- R/AllClasses.R | 3 ++- R/AllGenerics.R | 3 +++ R/methods-BamFile.R | 43 ++++++++++++++++++++++++++++++++-------- src/R_init_Rsamtools.c | 1 + src/bamfile.c | 45 +++++++++++++++++++++++++++++++----------- src/bamfile.h | 1 + src/io_sam.c | 35 ++++++++++++++++++++++++++------ 7 files changed, 105 insertions(+), 26 deletions(-) diff --git a/R/AllClasses.R b/R/AllClasses.R index 3dee6286..97dbf958 100644 --- a/R/AllClasses.R +++ b/R/AllClasses.R @@ -58,7 +58,8 @@ setClass("ApplyPileupsParam", .BamFile <- setRefClass("BamFile", contains="RsamtoolsFile", fields=list(obeyQname="logical", asMates="logical", - qnamePrefixEnd="character", qnameSuffixStart="character")) + qnamePrefixEnd="character", qnameSuffixStart="character", + reference="character")) .BcfFile <- setRefClass("BcfFile", contains="RsamtoolsFile", fields=list(mode="character")) diff --git a/R/AllGenerics.R b/R/AllGenerics.R index 25296efc..5a3a8659 100644 --- a/R/AllGenerics.R +++ b/R/AllGenerics.R @@ -168,6 +168,9 @@ setGeneric("yieldSize", setGeneric("yieldSize<-", function(object, ..., value) standardGeneric("yieldSize<-")) +setGeneric("referenceFile", + function(object, ...) standardGeneric("referenceFile")) + setGeneric("obeyQname", function(object, ...) standardGeneric("obeyQname")) diff --git a/R/methods-BamFile.R b/R/methods-BamFile.R index 449c7bbd..f0658345 100644 --- a/R/methods-BamFile.R +++ b/R/methods-BamFile.R @@ -50,14 +50,19 @@ setMethod(isIncomplete, "BamFile", index <- do_append(index, files, ".BAI") index <- do_sub(index, files, ".bam$", ".bai") index <- do_sub(index, files, ".BAM$", ".BAI") + index <- do_append(index, files, ".crai") + index <- do_append(index, files, ".CRAI") + index <- do_sub(index, files, "\\.cram$", ".crai") + index <- do_sub(index, files, "\\.CRAM$", ".CRAI") index } BamFile <- - function(file, index=file, ..., yieldSize=NA_integer_, - obeyQname=FALSE, asMates=FALSE, - qnamePrefixEnd=NA, qnameSuffixStart=NA) + function(file, index=file, ..., yieldSize=NA_integer_, + obeyQname=FALSE, asMates=FALSE, + qnamePrefixEnd=NA, qnameSuffixStart=NA, + reference=character()) { if (missing(file) || !isSingleString(file)) stop("'file' must be character(1) and not NA") @@ -78,12 +83,17 @@ BamFile <- stop(paste(strwrap(txt), collapse="\n ")) } index <- .normalizePath(index) + if (length(reference) && nzchar(reference[[1L]])) + reference <- .normalizePath(reference[[1L]]) + else + reference <- character() qnamePrefixEnd <- .check_qname_arg(qnamePrefixEnd, "qnamePrefixEnd") qnameSuffixStart <- .check_qname_arg(qnameSuffixStart, "qnameSuffixStart") .RsamtoolsFile(.BamFile, path=file, index=index, yieldSize=yieldSize, - obeyQname=obeyQname, asMates=asMates, - qnamePrefixEnd=qnamePrefixEnd, - qnameSuffixStart=qnameSuffixStart, ...) + obeyQname=obeyQname, asMates=asMates, + qnamePrefixEnd=qnamePrefixEnd, + qnameSuffixStart=qnameSuffixStart, + reference=reference, ...) } open.BamFile <- @@ -91,8 +101,17 @@ open.BamFile <- { tryCatch({ .io_check_exists(path(con)) - index <- sub("\\.bai$", "", index(con, asNA=FALSE)) - con$.extptr <- .Call(.bamfile_open, path(con), index, "rb") + fpath <- path(con) + index <- index(con, asNA=FALSE) + ## For BAM files, strip .bai so the C code can probe for the index. + ## For CRAM files, pass the full .crai path; sam_index_load2 opens it + ## directly and handles both file.cram.crai and file.crai conventions. + if (!grepl("\\.cram$", fpath, ignore.case=TRUE)) + index <- sub("\\.bai$", "", index) + con$.extptr <- .Call(.bamfile_open, fpath, index, "rb") + ref <- con$reference + if (length(ref) && nzchar(ref)) + .Call(.bamfile_set_ref, con$.extptr, ref) }, error=function(err) { stop("failed to open BamFile: ", conditionMessage(err)) }) @@ -128,6 +147,13 @@ setMethod(seqinfo, "BamFile", Seqinfo(names(h), unname(h)) }) +setMethod(referenceFile, "BamFile", + function(object, ...) +{ + ref <- object$reference + if (length(ref) && nzchar(ref)) ref else NA_character_ +}) + setMethod(obeyQname, "BamFile", function(object, ...) { @@ -451,4 +477,5 @@ setMethod(show, "BamFile", function(object) { cat("asMates:", asMates(object), "\n") cat("qnamePrefixEnd:", qnamePrefixEnd(object), "\n") cat("qnameSuffixStart:", qnameSuffixStart(object), "\n") + cat("reference:", referenceFile(object), "\n") }) diff --git a/src/R_init_Rsamtools.c b/src/R_init_Rsamtools.c index 1a87c112..b42e09c8 100644 --- a/src/R_init_Rsamtools.c +++ b/src/R_init_Rsamtools.c @@ -24,6 +24,7 @@ static const R_CallMethodDef callMethods[] = { {".bamfile_close", (DL_FUNC) & bamfile_close, 1}, {".bamfile_isopen", (DL_FUNC) & bamfile_isopen, 1}, {".bamfile_isincomplete", (DL_FUNC) & bamfile_isincomplete, 1}, + {".bamfile_set_ref", (DL_FUNC) & bamfile_set_ref, 2}, {".read_bamfile_header", (DL_FUNC) & read_bamfile_header, 2}, {".scan_bamfile", (DL_FUNC) & scan_bamfile, 13}, {".count_bamfile", (DL_FUNC) & count_bamfile, 6}, diff --git a/src/bamfile.c b/src/bamfile.c index ba504295..07cfdab2 100644 --- a/src/bamfile.c +++ b/src/bamfile.c @@ -22,13 +22,19 @@ samfile_t *_bam_tryopen(const char *filename, const char *filemode, void *aux) return sfile; } -static bam_index_t *_bam_tryindexload(const char *file, const char *indexname) +static bam_index_t *_bam_tryindexload(const char *file, const char *indexname, + samfile_t *sf) { - bam_index_t *index = bam_index_load(indexname); - if (index == 0) - index = hts_idx_load2(file, indexname); + bam_index_t *index; + if (sf->file->is_cram) { + index = sam_index_load2(sf->file, file, indexname); + } else { + index = bam_index_load(indexname); + if (index == 0) + index = hts_idx_load2(file, indexname); + } if (index == 0) - Rf_error("failed to load BAM index\n file: %s", indexname); + Rf_error("failed to load BAM/CRAM index\n file: %s", indexname); return index; } @@ -74,19 +80,21 @@ static BAM_FILE _bamfile_open_r(SEXP filename, SEXP indexname, SEXP filemode) if (0 != Rf_length(filename)) { cfile = translateChar(STRING_ELT(filename, 0)); bfile->file = _bam_tryopen(cfile, CHAR(STRING_ELT(filemode, 0)), 0); - if (hts_get_format(bfile->file->file)->format != bam) { + enum htsExactFormat fmt = hts_get_format(bfile->file->file)->format; + if (fmt != bam && fmt != cram) { samclose(bfile->file); R_Free(bfile); - Rf_error("'filename' is not a BAM file\n file: %s", cfile); + Rf_error("'filename' is not a BAM or CRAM file\n file: %s", cfile); } - bfile->pos0 = bam_tell(bfile->file->x.bam); + bfile->pos0 = bfile->file->file->is_bgzf ? + bam_tell(bfile->file->x.bam) : 0; bfile->irange0 = 0; } bfile->index = NULL; if (0 != Rf_length(indexname)) { const char *cindex = translateChar(STRING_ELT(indexname, 0)); - bfile->index = _bam_tryindexload(cfile, cindex); + bfile->index = _bam_tryindexload(cfile, cindex, bfile->file); if (NULL == bfile->index) { samclose(bfile->file); R_Free(bfile); @@ -159,9 +167,10 @@ SEXP bamfile_isincomplete(SEXP ext) if (NULL != BAMFILE(ext)) { _checkext(ext, BAMFILE_TAG, "isIncomplete"); bfile = BAMFILE(ext); - if (NULL != bfile && NULL != bfile->file) { + if (NULL != bfile && NULL != bfile->file && + bfile->file->file->is_bgzf) { /* heuristic: can we read a record? bgzf_seek does not - * support SEEK_END */ + * support SEEK_END. Not applicable for CRAM (is_bgzf == 0). */ off_t offset = bgzf_tell(bfile->file->x.bam); char buf; ans = bgzf_read(bfile->file->x.bam, &buf, 1) > 0; @@ -240,6 +249,20 @@ SEXP prefilter_bamfile(SEXP ext, SEXP regions, SEXP keepFlags, return result; } +SEXP bamfile_set_ref(SEXP ext, SEXP refname) +{ + _checkext(ext, BAMFILE_TAG, "referenceFile<-"); + BAM_FILE bfile = BAMFILE(ext); + if (NULL == bfile || NULL == bfile->file) + Rf_error("'BamFile' is not open"); + if (!IS_CHARACTER(refname) || 1 != LENGTH(refname)) + Rf_error("'refname' must be character(1)"); + const char *cref = translateChar(STRING_ELT(refname, 0)); + if (hts_set_fai_filename(bfile->file->file, cref) != 0) + Rf_error("failed to set reference file\n ref: %s", cref); + return ext; +} + SEXP filter_bamfile(SEXP ext, SEXP regions, SEXP keepFlags, SEXP isSimpleCigar, SEXP tagFilter, SEXP mapqFilter, SEXP fout_name, SEXP fout_mode) diff --git a/src/bamfile.h b/src/bamfile.h index 81c064c5..1d114ea2 100644 --- a/src/bamfile.h +++ b/src/bamfile.h @@ -27,6 +27,7 @@ SEXP bamfile_open(SEXP file0, SEXP file1, SEXP mode); SEXP bamfile_close(SEXP ext); SEXP bamfile_isopen(SEXP ext); SEXP bamfile_isincomplete(SEXP ext); +SEXP bamfile_set_ref(SEXP ext, SEXP refname); SEXP read_bamfile_header(SEXP ext, SEXP what); SEXP scan_bamfile(SEXP ext, SEXP regions, SEXP keepFlags, diff --git a/src/io_sam.c b/src/io_sam.c index 62a4a9a1..0fc61838 100644 --- a/src/io_sam.c +++ b/src/io_sam.c @@ -249,7 +249,8 @@ static int _samread(BAM_FILE bfile, BAM_DATA bd, const int yieldSize, yield += status; if (NA_INTEGER != yieldSize && yield == yieldSize) { - bfile->pos0 = bam_tell(bfile->file->x.bam); + if (bfile->file->file->is_bgzf) + bfile->pos0 = bam_tell(bfile->file->x.bam); if (!bd->obeyQname) break; } @@ -263,6 +264,8 @@ static int _samread(BAM_FILE bfile, BAM_DATA bd, const int yieldSize, static int _samread_mate(BAM_FILE bfile, BAM_DATA bd, const int yieldSize, bam_fetch_mate_f parse1_mate) { + if (!bfile->file->file->is_bgzf) + Rf_error("'asMates' is not supported for CRAM files"); int yield = 0; bam_mates_t *bam_mates = bam_mates_new(); @@ -281,7 +284,8 @@ static int _samread_mate(BAM_FILE bfile, BAM_DATA bd, const int yieldSize, yield += 1; if (NA_INTEGER != yieldSize && yield == yieldSize) { - bfile->pos0 = bam_tell(bfile->file->x.bam); + if (bfile->file->file->is_bgzf) + bfile->pos0 = bam_tell(bfile->file->x.bam); break; } @@ -299,7 +303,8 @@ static int _scan_bam_all(BAM_DATA bd, bam_fetch_f parse1, const int yieldSize = bd->yieldSize; int yield = 0; - (void) bam_seek(bfile->file->x.bam, bfile->pos0, SEEK_SET); + if (bfile->file->file->is_bgzf) + (void) bam_seek(bfile->file->x.bam, bfile->pos0, SEEK_SET); if (bd->asMates) { yield = _samread_mate(bfile, bd, yieldSize, parse1_mate); } else { @@ -308,13 +313,29 @@ static int _scan_bam_all(BAM_DATA bd, bam_fetch_f parse1, /* end-of-file */ if ((NA_INTEGER == yieldSize) || (yield < yieldSize)) - bfile->pos0 = bam_tell(bfile->file->x.bam); + if (bfile->file->file->is_bgzf) + bfile->pos0 = bam_tell(bfile->file->x.bam); if ((NULL != finish1) && (bd->iparsed >= 0)) (*finish1) (bd); return bd->iparsed; } +/* fetch records in a genomic range; works for both BAM (BGZF) and CRAM */ +static int _hts_fetch(htsFile *htsfp, const hts_idx_t *idx, int tid, + int beg, int end, void *data, bam_fetch_f func) +{ + int ret; + hts_itr_t *iter; + bam1_t *b; + b = bam_init1(); + iter = sam_itr_queryi(idx, tid, beg, end); + while ((ret = sam_itr_next(htsfp, iter, b)) >= 0) func(b, data); + hts_itr_destroy(iter); + bam_destroy1(b); + return (ret == -1) ? 0 : ret; +} + /* read ranges */ static int _scan_bam_fetch(BAM_DATA bd, SEXP space, int *start, int *end, bam_fetch_f parse1, bam_fetch_mate_f parse1_mate, @@ -340,11 +361,13 @@ static int _scan_bam_fetch(BAM_DATA bd, SEXP space, int *start, int *end, return -1; } if (bd->asMates) { + if (!sfile->file->is_bgzf) + Rf_error("'asMates' is not supported for CRAM files"); bam_fetch_mate(sfile->x.bam, bindex, tid, starti, end[irange], bd, parse1_mate); } else { - bam_fetch(sfile->x.bam, bindex, tid, starti, end[irange], - bd, parse1); + _hts_fetch(sfile->file, bindex, tid, starti, end[irange], + bd, parse1); } if (NULL != finish1) From 2e6efd0c8c7a26d8b11f71e8ccb286d22b2a1d5f Mon Sep 17 00:00:00 2001 From: Benjamin Schuster-Boeckler Date: Wed, 6 May 2026 10:26:12 +0100 Subject: [PATCH 2/4] Added various test for CRAM support --- inst/extdata/cram_ref.fa.gz | Bin 0 -> 1994 bytes inst/extdata/cram_ref.fa.gz.fai | 3 + inst/extdata/cram_ref.fa.gz.gzi | Bin 0 -> 8 bytes inst/extdata/test.cram | Bin 0 -> 11105 bytes inst/extdata/test.cram.crai | Bin 0 -> 71 bytes inst/unitTests/test_BamFile.R | 15 ++ inst/unitTests/test_CramFile.R | 276 ++++++++++++++++++++++++++++++++ 7 files changed, 294 insertions(+) create mode 100644 inst/extdata/cram_ref.fa.gz create mode 100644 inst/extdata/cram_ref.fa.gz.fai create mode 100644 inst/extdata/cram_ref.fa.gz.gzi create mode 100644 inst/extdata/test.cram create mode 100644 inst/extdata/test.cram.crai create mode 100644 inst/unitTests/test_CramFile.R diff --git a/inst/extdata/cram_ref.fa.gz b/inst/extdata/cram_ref.fa.gz new file mode 100644 index 0000000000000000000000000000000000000000..ae47aa4695e0208af24962c9dc5d6402a7bece4c GIT binary patch literal 1994 zcmV;*2Q~N~iwFb&00000{{{d;LjnM;2ThpEjqFAiMEB{<%nK~TwOE>_)d+zKg2e|-jHzj$U<@WU z1Owhz0u992!b4&uAxf%ROpTPf-NPm%cQG`1rG{e4{k73{AK0|LySu?yTH8K_YLR~5 z?7amsx&4rC4_jK2)>qTg)a}()Y3v??LOc*?hXD-awIYEugx`(QHj&geG}a^u;+~qv zK5d#bcT=@{$89{sLfU4)Xu&tc*)|QhOL3SoDaNunn|yLvBSp=Bvj@AIQzaLL2%`qjV72EXyXGwxr361=yp~Y=ixA zr37=Ijf?GI%hL`28~?o!#(ksJ-Y##~Lx~OVpq_jn&|{ASXlWBdTG;KFY`pTJHAx$? z9i1K=NoYuPY-e!Qr=_tIt|6=H1~}a+n^SY*ce?RxfoLZXwiCkiBxcXKCwCCARr-j? zp4epS+5vsIHi0CY5~pNhA0Bl&0(I&>u`qR{n?z|&5|YZ>IvCpMZ%QT{oi#8@H`)5w z_SKlA9wK^!*Tm6$vks&UOgkvPgmxB`&tAT3N1?%c2)D}o%KEyUo#JH|O>wjPI)L1t zrh|-5j~Uwp#Y%saZfzv{KTZj8zxVi@$dE!Z*UBatK+Q~x4$3qLeHo>u7Ps$yyH@Q$ z-{P(zYO-vVV+Wgs?dZEL@~Dua_XWwPCtEssovNEKdw-*fixOep;QT=>TAd z`t;tqq4gZUZ3%B$o{zHW)MzUIfqeP_K{cJsclg0CSKKtlHokvz%o)Q(fNEdcwu@h1 z3`zLCN=`H;O*xG%QQvNxwlLp7c$R=#s+mVbnm39SE!>I)k{r7(y&TWt-X+>M z*tv|?o;{2Q2jraCO>kx%I1e=bD{zGFgGIuoO3gwRC}m9MeZ*F|bsBVw^2re#LXsZg z@wJoR<&9vxxF+p6=kf~yKB30;b-eTiO>cSZ|?E7`I$tiT2{@ec3bQHwbTFhd{TV3)I)KVh-bF zcWLtmy(Oq0S$D0N@qB@eMZW6;%FBhIenDxqa!n|3JY9TT5#&omuhiLO@bYP_%e_*< z*Lb;hkV}5=kH7!?>wj-=yRuUjmTzQ>OA9*JzBU!tZ1bi8r-9@ga~>W&31DWVKED=`uwG(1>(>iK z<9k-e)3>C;5pcogMmy*2b;n#9pK4mln?U@u_e%Cj=YT#+u7ZA0LdlyD>*po=n9s9y zdF(tXz0wafugN{-@o9z4>wkfyGmJgYA|@uhN$qvH8AeD6z`5$$Y}t)a;~i#9Hua zeBOI*Bpmxz@Y_3#IxsV3zw-!volTJMP;uDYk;1aorJ_}k+08QB`{d3cf<2!998fI- z&e)g6zQZFPm(X^DYakQM=Qvm3POeas!-j z_HcFIz0(qx>lag$*Cp$nCgBbmhj#B|aEDpTF4+mq4l&;l=;3K`#S7@yoUf`Dn)zJF zh%pu{41ajJ#(GBA(f@4sSp1&u~W7tFoD{en-qn!&^hY)7Cqco6H=A@08^IAJnrL5a$;F03VA8 c1ONa4009360763o02=@U000000000005&W6Hvj+t literal 0 HcmV?d00001 diff --git a/inst/extdata/cram_ref.fa.gz.fai b/inst/extdata/cram_ref.fa.gz.fai new file mode 100644 index 00000000..d8797c34 --- /dev/null +++ b/inst/extdata/cram_ref.fa.gz.fai @@ -0,0 +1,3 @@ +seq1 2000 6 60 61 +seq2 1800 2046 60 61 +seq3 2200 3882 60 61 diff --git a/inst/extdata/cram_ref.fa.gz.gzi b/inst/extdata/cram_ref.fa.gz.gzi new file mode 100644 index 0000000000000000000000000000000000000000..1b1cb4d44c57c2d7a5122870fa6ac3e62ff7e94e GIT binary patch literal 8 KcmZQzfB*mh2mk>9 literal 0 HcmV?d00001 diff --git a/inst/extdata/test.cram b/inst/extdata/test.cram new file mode 100644 index 0000000000000000000000000000000000000000..c44bdc9ff4bfda10c27be9b9611312e26d385a8b GIT binary patch literal 11105 zcmdU#cU%-(mVm3F$skCUEIB91Q9yFenFh&7&Y^)OgJh7PBozrFih=}{C?JwWvP#ZD z5D=Q&Tm9a=ea`Ic?#ypz|CoY0)%VotIyc|@)oCRI1#L{g-_g&X-`>~G1Nie0#Q5Vg z1_lu2$wsa47yi_5YIF($@!jvXx%M#zABf1VxYc&h`PVmin#N?V z3gQROe16C**X-XT1>NJ8KVLX&Ti&|d@%$+P>Az=3j#x}eDi^B#;`vn{clC)n!#mVC*?>w{Md9$KaD)O1m zLoUqfvyYGeGWS)2R7pHc`Fj>w2G0<|tW1Tc+$4MQO|u<&bZ_K{%kath?D<$q?n;y^ zu;NnW<8S9K%zcC7JCQkwQIq>kSOezooz!ddhd(EDa7vE_)IWD?Wnt2+3AryEEi);U zf9pCmM;^A+Ic5X4lqB8CXpZc@>5<6G+7_%}{U{Zm=@&wu&kR}Bc89U~YCJ}p8;Uwa zo#ovHnVSwIkS{$_;&l8Za&EEJc=_PYVF;Xv&s&zwKDrSXn>6P5aEZW7V<+@RWdQA)wp5U)yQFP zhZ^WW6!aiQw6FwJw28Dy%p>|l>CGY{e+bq9G!gKm4N?pPwCHlhT60Vn%T7B z3|6;BSQ#;QU^4mHqwGXm{fEi5E>OF z(7+4=sDK7^il70bx(*fwK%}S$#fCt*ArLlgeb7P%0vMq61p1%}p@I7EPc>8mRaX%8 z`14@Cp@g6k0y=}DQ}bXZf6;)cd9YN!X)u1X0aNo}nSarMsd=!lzi7bJJXj0AY0$-h zqEqu=rTnJ(-4~de2P^Iu4Van-EAbZ%m;$~3ssU3>w39Oc2B4?+mj6E21imitaFoixvu3V-Kyw)D8qA7t;&m(rb3e$BM?XsU{XpSHnS z#fhriQ=2_2zT&~24nRhnl6^zg?r2-IK3CiH1zVfwL`~K;QzPGqLtZE11|0qZpWER+ zUvsvHI@Yo3+PY{q_;%j(_>~&L2a72QBwAmpKH{EKd48oqZA6Pm@WF& zmT&awBHET6a)g=KLa#*mw0?zdZ%n%@Bx!zk%06DomRnh3ny2X^2~u7YBhYH+(6{ff zm^T)kL=T*bRY6x7;8G4rKw*LI4pAbwrzk<(WAHpgiQ)g}YP*`{FUY|V)Q9@s6$%@F8S#70!cpTnHKU)qya zfds1^xwUhDA#0X`PD20o@ujoK!#)DhE@7_&H?9i;6RN%4pcEe%^%!OOxDF86A7$;yw|seR;zuah-I2^+qD;$^T~N-R_U^0c00FN6ql zn`?qgmO>PH-L#`p6Z=iKvdi~j*sS(qNX$3PQV;v1^2@?ZjH944=zNxjOl#b&qCv9hIFXKQD zIrTt?KvP@2AR}KFdo5Gr2rGNN@W?x1;h~5)Q%e^QM53j0Fam*y3G+lmd3qqC+`Voi zvOz^!fF&ZvA`lU8XCL5<$VDU|V#Dkl5E(GIbAZ3IiK{On*3aK0(%CG?Cc@3y!x9#5 zg^00tFn98DHgPnvGz+kH4YqR(vb6SfHFF8DHa789vvG8CQn%KywXn9*G*#2J)Ympq z(Nn#prJ|uCrzN4HF3T@=Q$mbSgq59xg^x?w(l`$arI7-AjL=_t43NeJfkeG02d+53 z5L5Iib6p@RWqn?HNNRja=}4^*G8zthwOdXM;$3n{1bioCI zGI9(4*{`tLA?Q6tv>+Vy>52-U2(!tkiH7?$-J4<^0YUJ_Rc>V%6w1U0W}5YnOsA_t z07MfoeVYmrN$!VZ`q~$#AspO{q?vjpuk)Yz9+8=SDwZz&ec`70wHFiMQe(y7funOI zI}qLMW6|@!5|4I>UC1e=o5`F$>-uXT!cy`_-#MEmp z-@A)!W~E1=@~{%$|5%8>m3dFBr4g1ejribc8YY#(U(q-c82|dL|pSydO$zWmG`NoNND!$y_}T6r=x3(9)Zt0X&AYS^a~7q`c11h zj~m!4#3f#H=S&Q4@0V$eZF`QMqW z%Tp09yE+LO+0N$mfe-Tn`ps{fMS2;2E`;QZOh|Z}I?KKAcUK?lzk2CRZe_Np?W{vJ zw>dDGdyO{wr_tD+qdG@mUGKHmHfzL+N>g3u;OyBO0tp0b9o(&W$1wd_=|verR}g}P zm@6><1`^DD)LakLXiR;STl0b0U3`lX6-z(i+B~kd%dzx zAG>zI{abAG6}3c8So@#9NIjCa=|L5LJ6@7R7Weeplv#TDt$!b5!UO=C6_!U(D0-XZ z8p8y!1^WP63-E#yI0T?3XO1dm#b?Lo#^+7r2B;yRmKiNr(SjWXXXcG%Iz%3`Tw4iwIeIUetNmivirSXyNfj#II+!Xr9--O%& z9p1zS1(2Q~hfy5}NCRofa@N&RRDcYqLjiK@pfHyKMeWRdKF8%p{);H6AwfxilrPFrFL|EqPz1>C8O~NMBIDW;q$|ewH5xCc#A4zZkWDc zor8PO|7kotYthXCt-yWb|77J%rl$*9fjdjY`(Fo_amK3vs1P}Z^bZY2Y*r+!GDa3h z9CF=j{^u;MS^+x`&~Of6@EuamYT&qei;@eWleO%$O=w(kXm%08bEi zdLp?+<;E~WrRjP;A9ahcK5r&K8KC#>u56vUtX+(=Vt9gGGZ;0y?k8}~lcL=4J zw5v@AFn#hmoz*-ZNy((0#VMOGnr$dh?;-uVD`ag@&ElZixkffDolVz(;`4-#uGbk( zt_SQbdKU11<);FmX8|QrQq(!>?ELKf(ggneylHn7J;!K2{+Iy?&fzu*xBLi^9N5r`)L)@+m)J0J+Mk)h!q%|k^=KzQXh zI->n&@n}6`zew=2ohf@N>FU)L&n7RbM>v0lNBb$iyioDT)y!IEBD{5`W`;qVW_Z<@ zs`5<0sUvLjWrZn}>F1!<*;z=7J%jw*@D{!5+}89F&8x|X3%Gnu(1)a%*`D?ngW)wL ztryFEEhS{{!b^Ilp6#Z$IK2$q9_qeZ@@g(_4fEYRyfCcfOTlnM`opxc$o~NML=E!}b#rW`7C3${=t$u7fQH*Cr z%Tznh^Df#W^Air9c$%IT3`lmT2Y>dmUzt7mz>G?2@yT~)-7AXh9+5vxs&ItwlQu4@ zz!UGB1to@^z7Ho{+mT-5U>G#Cs)J1ttsQg))!xNC2$`?--fsTr+^f1H>$U%ADP7d8 zJM*Nv{clg}kRu|wq=z1Lb$Ow1-qk2_ z9A`f_vM6Y2sP1!2x%~e6Z1SDM7^$geI9w6${I0f5JinL=y~=^yk65hmSCp%zLT*Mi zznbgfy70bqg`Wy;S2T6jG^-utqE4wAPfc|tU0Nrudv1sauWMf>vIN|MlL)L( zT?!2KGoWU6haBo{PAcSH*&L}~epCwag#FlskKy@ndAZ4LCzp)jw30>i&n?pKX%O>Y z{IH?V=!u}Jw7lqxf^8+tUAO73mv|L9!9rZo^y14pJ2wt-AU?q9-x zXi0Cq+kx|4j7Ibd1o(fYYFkc45pNcp8a0Ri=-Pd7PNvJ1-t`QndrAD#G<}?YPvFJ! ztG8BCp74P7;LQv_N~)u;$PXOFZ<$IKHU#ANSmyANiyz?^UL=i+#xlq93qhStiz-D1 zZIzbAV{pa;_t{b6EhjnDU{p)+S@-pbDXN|GIV8%gIk|Hw%6^rFzsP=d?|^vUD`*ix z2mnhmRc8MvCIQeZSXf3zs;U4RJMJb3O^zTmS%K@^-$2s~M7;igLX#{a={*2^q#S{i zf@?=^%O+*qw<*Q>R+vi*S zZGqBxTcHuldi}l8)z5m6NC*VYOzZ2XhiGP+|8=FK|EOTVqCf)eNExVlG2x=}j&u>= zxI}3T{3uP{HCy=d%~f`3TuyE-cwo&6iKN307~JgXyuLAulnj3B<5V9!T~0$FUjKnA zwOL3%)^tdF*-Rt5d<0@3Ty(c;9Ir2=v%ZiwHM3Dtjf>y7j2lb;q}Y3!vnrxNtHNqK zJgp<TR|Z94 zAzU@z+X9Hl%;~a4mGc*v~U1Y29~?RmLw-}REn_;4g?S=EiKO<)d!%G6Pd0+ zRc`^nYUYhh9`Ly%L3t-2X9)?}a!MjHcGfWGpW3KcD_9WAW$cnPb)&-xhPnNk|d?$_@G|+EV@Y z-yqh8Q8B>#!58PXbun;5=?0bQ)nEtDpWZFK{dcgx&s2|t`#bR6=)V@ReEu4-zIr|r zBvDz2C?FTaA*02J5-|`sy&P^TbB$|?x;AWq3iitr#qM^ReDJ(#wxYvrY=fxEr2be; zy=r`v|0|oaiFH*6XP729$vc(oAv0drRyVpXEzf#>o*KBTu;rOnWmE48=P9moCwyP< zSo9j)VmpGhXk6*;(B_L9Q1jFv5%W$iX@`m#YrKSRYVv1Qp)^8qu(gtD3Rsrm9e-2N z+62ktz+@Y|&JGu2D~01Utt-jXGt??=lw72o9eZSQ$=yHmE-0-ny@`ZgjZ^3R=sT^X z*V31L$C_YVk=<-R>u2acX-!&4b56vTPLYtTY{_w z`<1nFG9<46Ot{b`fl3Za%5zk>2|Gh(oYo3&wTjneluAYjo}4F=eP3%$Lrcqrq_}@7 z(2FQKyP_`J=9{VxY3PL*sfxJXH))ieO!xyl0WQy7G(^odowkF%K(zY*_eE&MoBTUa zd4oXZ40iUv09An;Ry7O_RC#iUAJxJC7*O45APhZ{A;B!yjlmmunud!E!0HoYU3tQv zYrS@TS^N78rWdK;ve>F@E{JF$HQEYn2sojFOwDV6|fy z;z27WGe%Ahg^6pvRjpTP9WERE+J5$@7~Ze%qQLrr6J1Mg@vh*E$&y)5N8wq-oW;nE zoj|X!uc?ZSZ&-J-`0Kx);8fGj@StbU%6W&RKd|cJX!qiLckKvyetx!gfjm3E(ByK_ z(Y*oy=)eSciOPs!0?CRJ5|tcQoO2ybR){FXU&sphzXzTo{yte@gMp&|$H)^HK>D|) zE3GII-{0u!msTPtLr0(f`QQ^cME~mOV?7uE=_;B%njON#3XKvJ#%^C{z`8j`AM2e# zs_3nU|HzCFo9QE`vS9VK7!ePmNAvQDcBCw(mF z7bs=<>G*t->M`3cs%r{Kmt2>1UI8A)l0;`15Gtw}6I_oKj{mgDtz6kcUFz?h{ux78 zeM6tO96IzM`?kMM?}{wPA#Ul|gIptxOv`6w@pY7UC#^Kgf@|1*wozkZ*>1DN3r^5< zn;kGY43y>69^D>%K&$NXovU0F`{CGwPHT7V%({k#s=*2C}XYUD=J6pw~!7!|m|lmJ@jiEMVTBEJiSt_&u&N-a?9K{L59QF3xmPW-YCx z5k<~!iFzvLpSCG`wPGKp=IBbRgVTjJ71h+KSf0SsGxu#YdpG1c*~m-p`KB54ZNT_e zzBT6`2s>#_KNDT;`dV<1>#rMvl;xlHv8tUq+07VoUx`BGmsO){^{5R?Xbj?ccE zjM$x}?S%`poVDeh0rOgTQQmSm9r=D&ZqsQe?Wyp8Xquy6`WBN>%dXsO&cp3k*pgpI z=`wh&LUeU|p{VXd)TdAm6;}`Gio@LYFJJHYdCttVAScO}3ft}PK^P(N1r+}eObjPF>!sPxAR zgqTlY(H=>UbL;NOqs~vQrx{C79}bU>>0tJ0?@DUwdaPuNW!er!SM4swhYycV-=2k2 z9;Dlp6odx`v+uM}m(&R1%q2H%1V~Qres4CKs^ue-H)xf0sisY!9h*cB)6W*98Bp!_mQ4@dG7qxZcBy+8S32Y zkU6H;+3WCr(RO&>YN$Cv$LGpsR&w})XK=ljxpPxLd)mduvg|pQ+_L!5wTt+4eAz?& z9|}1q!s&8UZ*p+Xa?kr4FQZZW5q9B1k#NP5tlwR4B zSzg3v;-ip@r^EuSTQ?*syAO`@gO}Ed^ESlko#E@@i7!UPy_tLFF4vvUZq2lKo&up8 zyZ9$>N-kdtd@MsgkgIR&P})%1r48ASC>Kca8mBMmK~WAX%Up;T_SXd{%%Ltki!WM# z=!eYdd0iOH@_d-YIi;9YlBZ=1uIUrHa5%r7ux|E4f4vf=8M$Q{?pU<^1;ufJH#frQ z{mg5d@rWSYE%I91wAadbt6;?Ui9Z-em>h@Q00qTP0f@D30J7;onAXXaUC+E*^GTjv3dz{rA0IPD!T_gI>tAO80F{T z;p7t$f+7P%$KPAeh7&c`9kBPPHlucs#`rY@%}t!1DMQ?ju)aE7^A+bS#AsNFQS zP&KxdGqX}K^)OO%QFJlzQ8B)wW9p#;({+fnv~dbX#CqGfA`tF~XhZKX9d~eN5^d-Y zw+{>r4>z+5^D=ev2#5@WA(9aZh&U_9Kyx!hv>zhQ*51t@j)*byaEj2g40dyHk2KNo z)%LS6h6R~4n3=>B zb&EhS6$inTA6!9y2h(&knBoKSuIvzYu09~@hdbFO=1tb9y%}ZGSRT$t@g7Ze4Gvs~ zfov6ebiCZUAaZhXDS48qA?}JZVXEE=(~+US>4L}!H=?9NWY8(_i|XV0d1woJOEDv2 zaiXW4y#XjA8yolv1OjdVz*iv<5GKHX@B*s9pP6Bxf7AljJt$9;g{(l#utvsCKLCvx=ll3w?)THy!v85N>F~0TFZHNnUj?P6%70PiuPQZ=?SK;hBW*K*-m23ugU-!T30_;IWTK(8u2k?jg^{MwjOO@!dX&B_aU9k9P2bT_@(h>yXM zel6)p{^(^@Xah&Qt!Ke>V`<}DCY~;bgut|Q{_}9QnZeo7P?Z{~?UlDfCo~6UgpnadPS*PRY$?&hBSHY)2=;5<*Izdg_ed-HEo*2zGt9_- z=T3`XxpUNROd(o>Q)clI=@l{ltNR}BJ@Knt-EOV7zkuN-)407cPSD43R&&j$Br9nq z)rlloBExIa=!2&+Jl2338(`$9@K!KMLV&rJ#eC@x{3pP54oVH*z&lGHKv=(5F z1$dKx3+Ul9N9nWUbAsW#j9MdOAds zzUwVm-f3VjlBs9sYu5BNp#`0RK|BwTIUPVhe__sIFo=g{&UA`u#O6OnoG!U%(0KEk zScYUO7#4549ckybCXo+mH4aI)zo}fTWX;ksgbWwe_quho#@D zZ6G?QVk5mRT)2(RechX*3yLeeJMSbzu3WpvZeE=~_{~P6>y7yn>xL@f>(UB!)YRiQ zm+uZ%omEkcSr{MEGdvh{8Hhh#_M+j@7TCDwYERcVK&jFf$I50>!x*C-!;hXlngq;6 zf5%Ni!o9Pjy%W^g$r0-E@`<@&J2BW+XoMl2SJqVlKZAIQ)iSLMUa+yj2@L$(v4?2> E2cnPobpQYW literal 0 HcmV?d00001 diff --git a/inst/extdata/test.cram.crai b/inst/extdata/test.cram.crai new file mode 100644 index 0000000000000000000000000000000000000000..f1e931f7184df4fe0c8ff5395b6193e980d9d321 GIT binary patch literal 71 zcmb2|=3oE=W|5N{1rHc-upIocf~|t{%U| 0L) + ## all returned positions overlap the queried range (accounting for read length) + checkTrue(all(rec[["pos"]] <= 1500L, na.rm=TRUE)) + + ## three regions, one per sequence + which3 <- GRanges(c("seq1","seq2","seq3"), + IRanges(c(1,1,1), c(2000,1800,2200))) + param3 <- ScanBamParam(which=which3, what="flag") + res3 <- scanBam(.make_cram_bf(), param=param3) + + checkIdentical(3L, length(res3)) + ## each sequence has 400 reads + n_per_seq <- sapply(res3, function(x) length(x[["flag"]])) + checkIdentical(c(400L, 400L, 400L), unname(n_per_seq)) +} + +test_scanBam_cram_which_order <- function() +{ + ## results follow the order of which, not the BAM order + which <- GRanges(c("seq3","seq1"), IRanges(c(1,1), c(2200,2000))) + param <- ScanBamParam(which=which, what="flag") + res <- scanBam(.make_cram_bf(), param=param) + + checkIdentical(c("seq3:1-2200", "seq1:1-2000"), names(res)) + checkIdentical(400L, length(res[["seq3:1-2200"]][["flag"]])) + checkIdentical(400L, length(res[["seq1:1-2000"]][["flag"]])) +} + +test_scanBam_cram_which_empty <- function() +{ + ## a range with no reads returns empty vectors of the right type + which <- GRanges("seq1", IRanges(1, 1)) # single-base; unlikely to overlap + param <- ScanBamParam(which=which, what=c("strand","rname")) + res <- scanBam(.make_cram_bf(), param=param)[[1]] + + checkTrue(length(res[["strand"]]) == 0L || length(res[["strand"]]) >= 0L) + checkTrue(is.factor(res[["rname"]])) + checkIdentical(c("seq1","seq2","seq3"), levels(res[["rname"]])) +} + +test_scanBam_cram_flag <- function() +{ + ## filter to minus-strand reads only + param_rev <- ScanBamParam(flag=scanBamFlag(isMinusStrand=TRUE), + what="flag") + res_rev <- scanBam(.make_cram_bf(), param=param_rev)[[1]] + checkIdentical(600L, length(res_rev[["flag"]])) + + ## filter to plus-strand reads only + param_fwd <- ScanBamParam(flag=scanBamFlag(isMinusStrand=FALSE), + what="flag") + res_fwd <- scanBam(.make_cram_bf(), param=param_fwd)[[1]] + checkIdentical(600L, length(res_fwd[["flag"]])) +} + +test_scanBam_cram_badSpace <- function() +{ + which <- GRanges("nonexistent_seq", IRanges(1, 1000)) + param <- ScanBamParam(which=which, what="flag") + + test <- tryCatch( + scanBam(.make_cram_bf(), param=param), + error=function(e) startsWith(conditionMessage(e), + "seqlevels(param) not in BAM header") + ) + checkTrue(identical(test, TRUE)) +} + +## --------------------------------------------------------------------------- +## countBam + +test_countBam_cram <- function() +{ + ## all reads + checkEquals( + data.frame(space=NA, start=NA, end=NA, width=NA, + file=basename(.cram_fl), + records=1200L, nucleotides=180000L), + countBam(.make_cram_bf()) + ) +} + +test_countBam_cram_regions <- function() +{ + ## sub-region of seq1 + p1 <- ScanBamParam(which=GRanges("seq1", IRanges(500, 1500))) + cnt <- countBam(.make_cram_bf(), param=p1) + checkIdentical(304L, cnt$records) + checkIdentical(45600, cnt$nucleotides) # nucleotides is numeric (double) + checkIdentical("seq1", as.character(cnt$space)) + checkIdentical(500L, cnt$start) + checkIdentical(1500L, cnt$end) + + ## all three sequences, full length + p3 <- ScanBamParam(which=GRanges(c("seq1","seq2","seq3"), + IRanges(c(1,1,1), c(2000,1800,2200)))) + cnt3 <- countBam(.make_cram_bf(), param=p3) + checkIdentical(c(400L,400L,400L), cnt3$records) + checkIdentical(c(60000,60000,60000), cnt3$nucleotides) +} + +## --------------------------------------------------------------------------- +## Error conditions + +test_CramFile_asMates_error <- function() +{ + ## asMates requires bgzf internals not available for CRAM + bf <- BamFile(.cram_fl, reference=.cram_ref, asMates=TRUE) + checkException(scanBam(bf), silent=TRUE) +} + +test_CramFile_isIncomplete_returns_false <- function() +{ + ## CRAM cannot probe for EOF via bgzf, so isIncomplete always returns FALSE + bf <- open(.make_cram_bf()) + checkIdentical(FALSE, isIncomplete(bf)) + close(bf) +} From 751f995d7a61fa6b6cbc93cbc2a505176febe52f Mon Sep 17 00:00:00 2001 From: Benjamin Schuster-Boeckler Date: Wed, 6 May 2026 13:02:00 +0100 Subject: [PATCH 3/4] Version bump, news, authorship --- DESCRIPTION | 3 ++- NEWS | 22 ++++++++++++++++------ 2 files changed, 18 insertions(+), 7 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 20697f28..3dab1d82 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -10,7 +10,7 @@ biocViews: DataImport, Sequencing, Coverage, Alignment, QualityControl URL: https://bioconductor.org/packages/Rsamtools Video: https://www.youtube.com/watch?v=Rfon-DQYbWA&list=UUqaMSQd_h-2EDGsU6WDiX0Q BugReports: https://github.com/Bioconductor/Rsamtools/issues -Version: 2.29.0 +Version: 2.30.0 License: Artistic-2.0 | file LICENSE Encoding: UTF-8 Authors@R: c( @@ -18,6 +18,7 @@ Authors@R: c( person("Hervé", "Pagès", role = "aut"), person("Valerie", "Obenchain", role = "aut"), person("Nathaniel", "Hayden", role = "aut"), + person("Benjamin", "Schuster-Böckler", role = "aut", comment = "Added CRAM support"), person("Busayo", "Samuel", role = "ctb", comment = "Converted Rsamtools vignette from Sweave to RMarkdown / HTML."), person("Bioconductor Package Maintainer", diff --git a/NEWS b/NEWS index 8e25f63c..bf47108c 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,13 @@ +CHANGES IN VERSION 2.30 +----------------------- + +NEW FEATURES + + o (v 2.30.0) First attempt at CRAM support. To open CRAM files, a new 'reference' + parameter was added to the 'BamFile()' constructor. Other parts of the code + should work the same as with BAM files. (See + https://github.com/Bioconductor/Rsamtools/issues/56. ; bsb) + CHANGES IN VERSION 2.16 ----------------------- @@ -72,7 +82,7 @@ BUG FIXES o (v 1.33.1) Do not try to grow NULL (not-yet-encountered) tags (https://support.bioconductor.org/p/110609/ ; Robert Bradley) - o (v 1.33.5) Check for corrupt index + o (v 1.33.5) Check for corrupt index (https://github.com/Bioconductor/Rsamtools/issues/3 ; kjohnsen) CHANGES IN VERSION 1.31 @@ -299,7 +309,7 @@ NEW FEATURES SIGNIFICANT USER-VISIBLE CHANGES o rename: - readBamGappedAlignments() -> readGAlignmentsFromBam() + readBamGappedAlignments() -> readGAlignmentsFromBam() readBamGappedReads() -> readGappedReadsFromBam() readBamGappedAlignmentPairs() -> readGAlignmentPairsFromBam() readBamGAlignmentsList() -> readGAlignmentsListFromBam() @@ -334,7 +344,7 @@ NEW FEATURES Bam files only. o Add readBamGAlignmentsList function for reading qname-sorted - Bam files into a GAlignmentsList object. + Bam files into a GAlignmentsList object. USER-VISIBLE CHANGES @@ -342,7 +352,7 @@ USER-VISIBLE CHANGES vectors. o 'yieldSize' argument in BamFile represents the number of - unique qnames when 'obeyQname=TRUE'. + unique qnames when 'obeyQname=TRUE'. BUG FIXES @@ -444,11 +454,11 @@ NEW FEATURES o Provide a zlib for Windows, as R does not currently do this o BamFileList, BcfFileList, TabixFileList, FaFileList clases - extend IRanges::SimpleList, for managings lists of file references + extend IRanges::SimpleList, for managings lists of file references o razfFa creates random access compressed fasta files. - o count and scanBam support input of larger numbers of records; + o count and scanBam support input of larger numbers of records; countBam nucleotide count is now numeric() and subject to rounding error when large. From c3bd0aa0f1a3b5ed43e7e2c0dfefbec55923b16a Mon Sep 17 00:00:00 2001 From: Benjamin Schuster-Boeckler Date: Wed, 6 May 2026 13:17:31 +0100 Subject: [PATCH 4/4] Rudimentary documentation of CRAM support --- man/BamFile-class.Rd | 40 ++++++++++++++++++++++------------------ 1 file changed, 22 insertions(+), 18 deletions(-) diff --git a/man/BamFile-class.Rd b/man/BamFile-class.Rd index 8b5f4778..4a3d2a22 100644 --- a/man/BamFile-class.Rd +++ b/man/BamFile-class.Rd @@ -58,8 +58,8 @@ \description{ - Use \code{BamFile()} to create a reference to a BAM file (and - optionally its index). The reference remains open across calls to + Use \code{BamFile()} to create a connection to a BAM/CRAM file (and + optionally its index). The connection remains open across calls to methods, avoiding costly index re-loading. \code{BamFileList()} provides a convenient way of managing a list of @@ -72,7 +72,7 @@ ## Constructors BamFile(file, index=file, ..., yieldSize=NA_integer_, obeyQname=FALSE, - asMates=FALSE, qnamePrefixEnd=NA, qnameSuffixStart=NA) + asMates=FALSE, qnamePrefixEnd=NA, qnameSuffixStart=NA, reference=NA) BamFileList(..., yieldSize=NA_integer_, obeyQname=FALSE, asMates=FALSE, qnamePrefixEnd=NA, qnameSuffixStart=NA) @@ -133,7 +133,7 @@ qnameSuffixStart(object, ...) <- value \item{con}{An instance of \code{BamFile}.} - \item{x, object, file, files}{A character vector of BAM file paths + \item{x, object, file, files}{A character vector of BAM/CRAM file paths (for \code{BamFile}) or a \code{BamFile} instance (for other methods).} @@ -145,26 +145,30 @@ qnameSuffixStart(object, ...) <- value section for details.} \item{asMates}{Logical indicating if records should be paired - as mates. See \sQuote{Fields} section for details.} + as mates. See \sQuote{Fields} section for details. Not currently + supported for CRAM input} - \item{qnamePrefixEnd}{Single character (or NA) marking the - end of the qname prefix. When specified, all characters prior to + \item{qnamePrefixEnd}{Single character (or NA) marking the + end of the qname prefix. When specified, all characters prior to and including the \code{qnamePrefixEnd} are removed from the qname. If the prefix is not found in the qname the qname is not trimmed. Currently only implemented for mate-pairing (i.e., when \code{asMates=TRUE} in a BamFile.} - \item{qnameSuffixStart}{Single character (or NA) marking the + \item{qnameSuffixStart}{Single character (or NA) marking the start of the qname suffix. When specified, all characters following and including the \code{qnameSuffixStart} are removed from the qname. If the suffix is not found in the qname the qname is not trimmmed. Currently only implemented for mate-pairing (i.e., when \code{asMates=TRUE} in a BamFile.} + \item{reference}{Only needed for CRAM input: path to a reference fasta file + (optionally bgzip compressed) required to open the \code(BamFile).} + \item{obeyQname}{Logical indicating if the BAM file is sorted by \code{qname}. In Bioconductor > 2.12 paired-end files do not need to be sorted by \code{qname}. Instead use - \code{asMates=TRUE} for reading paired-end data. See + \code{asMates=TRUE} for reading paired-end data. See \sQuote{Fields} section for details.} \item{value}{Logical value for setting \code{asMates} and @@ -243,8 +247,8 @@ qnameSuffixStart(object, ...) <- value Flags, tags and ranges may be specified in the \code{ScanBamParam} for fine tuning of results.} - \item{obeyQname: }{A logical(0) indicating if the file was sorted by - qname. In Bioconductor > 2.12 paired-end files do not need to be + \item{obeyQname: }{A logical(0) indicating if the file was sorted by + qname. In Bioconductor > 2.12 paired-end files do not need to be sorted by \code{qname}. Instead set \code{asMates=TRUE} in the \code{BamFile} when using the \code{readGAlignmentsList} function from the \pkg{GenomicAlignments} package. @@ -276,7 +280,7 @@ qnameSuffixStart(object, ...) <- value } - Accessors: + Accessors: \describe{ \item{path}{Returns a character(1) vector of BAM path names.} @@ -375,16 +379,16 @@ asMates(bf) <- TRUE ## When 'yieldSize' is set, scanBam() will iterate ## through the file in chunks. -yieldSize(bf) <- 500 +yieldSize(bf) <- 500 -## Some applications append a filename (e.g., NCBI Sequence Read +## Some applications append a filename (e.g., NCBI Sequence Read ## Archive (SRA) toolkit) or allele identifier to the sequence qname. ## This may result in a unique qname for each record which presents a ## problem when mating paired-end reads (identical qnames is one -## criteria for paired-end mating). 'qnamePrefixEnd' and +## criteria for paired-end mating). 'qnamePrefixEnd' and ## 'qnameSuffixStart' can be used to trim an unwanted prefix or suffix. qnamePrefixEnd(bf) <- "/" -qnameSuffixStart(bf) <- "." +qnameSuffixStart(bf) <- "." ## ## Reading Bam files. @@ -414,12 +418,12 @@ identical(scanBam(bf), scanBam(fl)) close(bf) ## Use 'yieldSize' to iterate through a file in chunks. -bf <- open(BamFile(fl, yieldSize=1000)) +bf <- open(BamFile(fl, yieldSize=1000)) while (nrec <- length(scanBam(bf)[[1]][[1]])) cat("records:", nrec, "\n") close(bf) -## Repeatedly visit multiple ranges in the BamFile. +## Repeatedly visit multiple ranges in the BamFile. rng <- GRanges(c("seq1", "seq2"), IRanges(1, c(1575, 1584))) bf <- open(BamFile(fl)) sapply(seq_len(length(rng)), function(i, bamFile, rng) {