1 /// @file htslib/synced_bcf_reader.h 2 /// Stream through multiple VCF files. 3 /* 4 Copyright (C) 2012-2017, 2019-2021 Genome Research Ltd. 5 6 Author: Petr Danecek <pd3@sanger.ac.uk> 7 8 Permission is hereby granted, free of charge, to any person obtaining a copy 9 of this software and associated documentation files (the "Software"), to deal 10 in the Software without restriction, including without limitation the rights 11 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell 12 copies of the Software, and to permit persons to whom the Software is 13 furnished to do so, subject to the following conditions: 14 15 The above copyright notice and this permission notice shall be included in 16 all copies or substantial portions of the Software. 17 18 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 19 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 20 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 21 THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 22 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING 23 FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER 24 DEALINGS IN THE SOFTWARE. */ 25 26 /* 27 The synced_bcf_reader allows to keep multiple VCFs open and stream them 28 using the next_line iterator in a seamless matter without worrying about 29 chromosomes and synchronizing the sites. This is used by vcfcheck to 30 compare multiple VCFs simultaneously and is used also for merging, 31 creating intersections, etc. 32 33 The synced_bcf_reader also provides API for reading indexed BCF/VCF, 34 hiding differences in BCF/VCF opening, indexing and reading. 35 36 37 Example of usage: 38 39 bcf_srs_t *sr = bcf_sr_init(); 40 bcf_sr_set_opt(sr, BCF_SR_PAIR_LOGIC, BCF_SR_PAIR_BOTH_REF); 41 bcf_sr_set_opt(sr, BCF_SR_REQUIRE_IDX); 42 for (i=0; i<nfiles; i++) 43 bcf_sr_add_reader(sr,files[i]); 44 while ( bcf_sr_next_line(sr) ) 45 { 46 for (i=0; i<nfiles; i++) 47 { 48 bcf1_t *line = bcf_sr_get_line(sr,i); 49 ... 50 } 51 } 52 if ( sr->errnum ) error("Error: %s\n", bcf_sr_strerror(sr->errnum)); 53 bcf_sr_destroy(sr); 54 */ 55 56 #ifndef HTSLIB_SYNCED_BCF_READER_H 57 #define HTSLIB_SYNCED_BCF_READER_H 58 59 #include "hts.h" 60 #include "vcf.h" 61 #include "tbx.h" 62 63 #ifdef __cplusplus 64 extern "C" { 65 #endif 66 67 /* 68 When reading multiple files in parallel, duplicate records within each 69 file will be reordered and offered in intuitive order. For example, 70 when reading two files, each with unsorted SNP and indel record, the 71 reader should return the SNP records together and the indel records 72 together. The logic of compatible records can vary depending on the 73 application and can be set using the PAIR_* defined below. 74 75 The COLLAPSE_* definitions will be deprecated in future versions, please 76 use the PAIR_* definitions instead. 77 */ 78 #define COLLAPSE_NONE 0 // require the exact same set of alleles in all files 79 #define COLLAPSE_SNPS 1 // allow different alleles, as long as they all are SNPs 80 #define COLLAPSE_INDELS 2 // the same as above, but with indels 81 #define COLLAPSE_ANY 4 // any combination of alleles can be returned by bcf_sr_next_line() 82 #define COLLAPSE_SOME 8 // at least some of the ALTs must match 83 #define COLLAPSE_BOTH (COLLAPSE_SNPS|COLLAPSE_INDELS) 84 85 #define BCF_SR_PAIR_SNPS (1<<0) // allow different alleles, as long as they all are SNPs 86 #define BCF_SR_PAIR_INDELS (1<<1) // the same as above, but with indels 87 #define BCF_SR_PAIR_ANY (1<<2) // any combination of alleles can be returned by bcf_sr_next_line() 88 #define BCF_SR_PAIR_SOME (1<<3) // at least some of multiallelic ALTs must match. Implied by all the others with the exception of EXACT 89 #define BCF_SR_PAIR_SNP_REF (1<<4) // allow REF-only records with SNPs 90 #define BCF_SR_PAIR_INDEL_REF (1<<5) // allow REF-only records with indels 91 #define BCF_SR_PAIR_EXACT (1<<6) // require the exact same set of alleles in all files 92 #define BCF_SR_PAIR_BOTH (BCF_SR_PAIR_SNPS|BCF_SR_PAIR_INDELS) 93 #define BCF_SR_PAIR_BOTH_REF (BCF_SR_PAIR_SNPS|BCF_SR_PAIR_INDELS|BCF_SR_PAIR_SNP_REF|BCF_SR_PAIR_INDEL_REF) 94 95 typedef enum 96 { 97 BCF_SR_REQUIRE_IDX, 98 BCF_SR_PAIR_LOGIC, // combination of the PAIR_* values above 99 BCF_SR_ALLOW_NO_IDX, // allow to proceed even if required index is not present (at the user's risk) 100 BCF_SR_REGIONS_OVERLAP, // include overlapping records with POS outside the regions: 0=no, 1=VCF line overlap, 2=true variant overlap [1] 101 BCF_SR_TARGETS_OVERLAP // include overlapping records with POS outside the targets: 0=no, 1=VCF line overlap, 2=true variant overlap [0] 102 } 103 bcf_sr_opt_t; 104 105 struct bcf_sr_region_t; 106 107 typedef struct bcf_sr_regions_t 108 { 109 // for reading from tabix-indexed file (big data) 110 tbx_t *tbx; // tabix index 111 hts_itr_t *itr; // tabix iterator 112 kstring_t line; // holder of the current line, set only when reading from tabix-indexed files 113 htsFile *file; 114 char *fname; 115 int is_bin; // is open in binary mode (tabix access) 116 char **als; // parsed alleles if targets_als set and _regions_match_alleles called 117 kstring_t als_str; // block of parsed alleles 118 int nals, mals; // number of set alleles and the size of allocated array 119 int als_type; // alleles type, currently VCF_SNP or VCF_INDEL 120 121 // user handler to deal with skipped regions without a counterpart in VCFs 122 void (*missed_reg_handler)(struct bcf_sr_regions_t *, void *); 123 void *missed_reg_data; 124 125 // for in-memory regions (small data) 126 struct bcf_sr_region_t *regs; // the regions 127 128 // shared by both tabix-index and in-memory regions 129 void *seq_hash; // keys: sequence names, values: index to seqs 130 char **seq_names; // sequence names 131 int nseqs; // number of sequences (chromosomes) in the file 132 int iseq; // current position: chr name, index to snames 133 hts_pos_t start, end; // current position: start, end of the region (0-based) 134 int prev_seq; 135 hts_pos_t prev_start, prev_end; 136 int overlap; // see BCF_SR_REGIONS_OVERLAP/BCF_SR_TARGETS_OVERLAP 137 } 138 bcf_sr_regions_t; 139 140 typedef struct bcf_sr_t 141 { 142 htsFile *file; 143 tbx_t *tbx_idx; 144 hts_idx_t *bcf_idx; 145 bcf_hdr_t *header; 146 hts_itr_t *itr; 147 char *fname; 148 bcf1_t **buffer; // cached VCF records. First is the current record synced across the reader 149 int nbuffer, mbuffer; // number of cached records (including the current record); number of allocated records 150 int nfilter_ids, *filter_ids; // -1 for ".", otherwise filter id as returned by bcf_hdr_id2int 151 int *samples, n_smpl; // list of columns in the order consistent with bcf_srs_t.samples 152 } 153 bcf_sr_t; 154 155 typedef enum 156 { 157 open_failed, not_bgzf, idx_load_failed, file_type_error, api_usage_error, 158 header_error, no_eof, no_memory, vcf_parse_error, bcf_read_error, noidx_error 159 } 160 bcf_sr_error; 161 162 typedef struct bcf_srs_t 163 { 164 // Parameters controlling the logic 165 int collapse; // Do not access directly, use bcf_sr_set_pairing_logic() instead 166 char *apply_filters; // If set, sites where none of the FILTER strings is listed 167 // will be skipped. Active only at the time of 168 // initialization, that is during the add_reader() 169 // calls. Therefore, each reader can be initialized with different 170 // filters. 171 int require_index; // Some tools do not need random access 172 int max_unpack; // When reading VCFs and knowing some fields will not be needed, boost performance of vcf_parse1 173 int *has_line; // Corresponds to return value of bcf_sr_next_line but is not limited by sizeof(int). Use bcf_sr_has_line macro to query. 174 bcf_sr_error errnum; 175 176 // Auxiliary data 177 bcf_sr_t *readers; 178 int nreaders; 179 int streaming; // reading mode: index-jumping or streaming 180 int explicit_regs; // was the list of regions se by bcf_sr_set_regions or guessed from tabix index? 181 char **samples; // List of samples 182 bcf_sr_regions_t *regions, *targets; // see bcf_sr_set_[targets|regions] for description 183 int targets_als; // subset to targets not only by position but also by alleles? 184 int targets_exclude; 185 kstring_t tmps; 186 int n_smpl; 187 188 int n_threads; // Simple multi-threaded decoding / encoding. 189 htsThreadPool *p; // Our pool, but it can be used by others if needed. 190 void *aux; // Opaque auxiliary data 191 } 192 bcf_srs_t; 193 194 /** Allocate and initialize a bcf_srs_t struct. 195 * 196 * The bcf_srs_t struct returned by a successful call should be freed 197 * via bcf_sr_destroy() when it is no longer needed. 198 */ 199 HTSLIB_EXPORT 200 bcf_srs_t *bcf_sr_init(void); 201 202 /** Destroy a bcf_srs_t struct */ 203 HTSLIB_EXPORT 204 void bcf_sr_destroy(bcf_srs_t *readers); 205 206 HTSLIB_EXPORT 207 char *bcf_sr_strerror(int errnum); 208 209 HTSLIB_EXPORT 210 int bcf_sr_set_opt(bcf_srs_t *readers, bcf_sr_opt_t opt, ...); 211 212 213 /** 214 * bcf_sr_set_threads() - allocates a thread-pool for use by the synced reader. 215 * @n_threads: size of thread pool 216 * 217 * Returns 0 if the call succeeded, or <0 on error. 218 */ 219 HTSLIB_EXPORT 220 int bcf_sr_set_threads(bcf_srs_t *files, int n_threads); 221 222 /** Deallocates thread memory, if owned by us. */ 223 HTSLIB_EXPORT 224 void bcf_sr_destroy_threads(bcf_srs_t *files); 225 226 /** 227 * bcf_sr_add_reader() - open new reader 228 * @readers: holder of the open readers 229 * @fname: the VCF file 230 * 231 * Returns 1 if the call succeeded, or 0 on error. 232 * 233 * See also the bcf_srs_t data structure for parameters controlling 234 * the reader's logic. 235 */ 236 HTSLIB_EXPORT 237 int bcf_sr_add_reader(bcf_srs_t *readers, const char *fname); 238 239 HTSLIB_EXPORT 240 void bcf_sr_remove_reader(bcf_srs_t *files, int i); 241 242 /** 243 * bcf_sr_next_line() - the iterator 244 * @readers: holder of the open readers 245 * 246 * Returns the number of readers which have the current line 247 * (bcf_sr_t.buffer[0]) set at this position. Use the bcf_sr_has_line macro to 248 * determine which of the readers are set. 249 */ 250 HTSLIB_EXPORT 251 int bcf_sr_next_line(bcf_srs_t *readers); 252 253 #define bcf_sr_has_line(readers, i) (readers)->has_line[i] 254 #define bcf_sr_get_line(_readers, i) ((_readers)->has_line[i] ? ((_readers)->readers[i].buffer[0]) : (bcf1_t *) NULL) 255 #define bcf_sr_swap_line(_readers, i, lieu) { bcf1_t *tmp = lieu; lieu = (_readers)->readers[i].buffer[0]; (_readers)->readers[i].buffer[0] = tmp; } 256 #define bcf_sr_region_done(_readers,i) (!(_readers)->has_line[i] && !(_readers)->readers[i].nbuffer ? 1 : 0) 257 #define bcf_sr_get_header(_readers, i) (_readers)->readers[i].header 258 #define bcf_sr_get_reader(_readers, i) &((_readers)->readers[i]) 259 260 261 /** 262 * bcf_sr_seek() - set all readers to selected position 263 * @seq: sequence name; NULL to seek to start 264 * @pos: 0-based coordinate 265 */ 266 HTSLIB_EXPORT 267 int bcf_sr_seek(bcf_srs_t *readers, const char *seq, hts_pos_t pos); 268 269 /** 270 * bcf_sr_set_samples() - sets active samples 271 * @readers: holder of the open readers 272 * @samples: this can be one of: file name with one sample per line; 273 * or column-separated list of samples; or '-' for a list of 274 * samples shared by all files. If first character is the 275 * exclamation mark, all but the listed samples are included. 276 * @is_file: 0: list of samples; 1: file with sample names 277 * 278 * Returns 1 if the call succeeded, or 0 on error. 279 */ 280 HTSLIB_EXPORT 281 int bcf_sr_set_samples(bcf_srs_t *readers, const char *samples, int is_file); 282 283 /** 284 * bcf_sr_set_targets(), bcf_sr_set_regions() - init targets/regions 285 * @readers: holder of the open readers 286 * @targets: list of regions, one-based and inclusive. 287 * @is_fname: 0: targets is a comma-separated list of regions (chr,chr:from-to) 288 * 1: targets is a tabix indexed file with a list of regions 289 * (<chr,pos> or <chr,from,to>) 290 * 291 * Returns 0 if the call succeeded, or -1 on error. 292 * 293 * Both functions behave the same way, unlisted positions will be skipped by 294 * bcf_sr_next_line(). However, there is an important difference: regions use 295 * index to jump to desired positions while targets streams the whole files 296 * and merely skip unlisted positions. 297 * 298 * Moreover, bcf_sr_set_targets() accepts an optional parameter $alleles which 299 * is interpreted as a 1-based column index in the tab-delimited file where 300 * alleles are listed. This in principle enables to perform the COLLAPSE_* 301 * logic also with tab-delimited files. However, the current implementation 302 * considers the alleles merely as a suggestion for prioritizing one of possibly 303 * duplicate VCF lines. It is up to the caller to examine targets->als if 304 * perfect match is sought after. Note that the duplicate positions in targets 305 * file are currently not supported. 306 * Targets (but not regions) can be prefixed with "^" to request logical complement, 307 * for example "^X,Y,MT" indicates that sequences X, Y and MT should be skipped. 308 * 309 * API note: bcf_sr_set_regions/bcf_sr_set_targets MUST be called before the 310 * first call to bcf_sr_add_reader(). 311 */ 312 HTSLIB_EXPORT 313 int bcf_sr_set_targets(bcf_srs_t *readers, const char *targets, int is_file, int alleles); 314 315 HTSLIB_EXPORT 316 int bcf_sr_set_regions(bcf_srs_t *readers, const char *regions, int is_file); 317 318 319 320 /* 321 * bcf_sr_regions_init() 322 * @regions: regions can be either a comma-separated list of regions 323 * (chr|chr:pos|chr:from-to|chr:from-) or VCF, BED, or 324 * tab-delimited file (the default). Uncompressed files 325 * are stored in memory while bgzip-compressed and tabix-indexed 326 * region files are streamed. 327 * @is_file: 0: regions is a comma-separated list of regions 328 * (chr|chr:pos|chr:from-to|chr:from-) 329 * 1: VCF, BED or tab-delimited file 330 * @chr, from, to: 331 * Column indexes of chromosome, start position and end position 332 * in the tab-delimited file. The positions are 1-based and 333 * inclusive. 334 * These parameters are ignored when reading from VCF, BED or 335 * tabix-indexed files. When end position column is not present, 336 * supply 'from' in place of 'to'. When 'to' is negative, first 337 * abs(to) will be attempted and if that fails, 'from' will be used 338 * instead. 339 * 340 * The bcf_sr_regions_t struct returned by a successful call should be freed 341 * via bcf_sr_regions_destroy() when it is no longer needed. 342 */ 343 HTSLIB_EXPORT 344 bcf_sr_regions_t *bcf_sr_regions_init(const char *regions, int is_file, int chr, int from, int to); 345 346 HTSLIB_EXPORT 347 void bcf_sr_regions_destroy(bcf_sr_regions_t *regions); 348 349 /* 350 * bcf_sr_regions_seek() - seek to the chromosome block 351 * 352 * Returns 0 on success or -1 on failure. Sets reg->seq appropriately and 353 * reg->start,reg->end to -1. 354 */ 355 HTSLIB_EXPORT 356 int bcf_sr_regions_seek(bcf_sr_regions_t *regions, const char *chr); 357 358 /* 359 * bcf_sr_regions_next() - retrieves next region. Returns 0 on success and -1 360 * when all regions have been read. The fields reg->seq, reg->start and 361 * reg->end are filled with the genomic coordinates on success or with 362 * NULL,-1,-1 when no region is available. The coordinates are 0-based, 363 * inclusive. 364 */ 365 HTSLIB_EXPORT 366 int bcf_sr_regions_next(bcf_sr_regions_t *reg); 367 368 /* 369 * bcf_sr_regions_overlap() - checks if the interval <start,end> overlaps any of 370 * the regions, the coordinates are 0-based, inclusive. The coordinate queries 371 * must come in ascending order. 372 * 373 * Returns 0 if the position is in regions; -1 if the position is not in the 374 * regions and more regions exist; -2 if not in the regions and there are no more 375 * regions left. 376 */ 377 HTSLIB_EXPORT 378 int bcf_sr_regions_overlap(bcf_sr_regions_t *reg, const char *seq, hts_pos_t start, hts_pos_t end); 379 380 /* 381 * bcf_sr_regions_flush() - calls repeatedly regs->missed_reg_handler() until 382 * all remaining records are processed. 383 * Returns 0 on success, <0 on error. 384 */ 385 HTSLIB_EXPORT 386 int bcf_sr_regions_flush(bcf_sr_regions_t *regs); 387 388 #ifdef __cplusplus 389 } 390 #endif 391 392 #endif 393