/usr/include/htslib/synced_bcf_reader.h is in libhts-dev 0.2.0~rc3-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 | /*
The synced_bcf_reader allows to keep multiple VCFs open and stream them
using the next_line iterator in a seamless matter without worrying about
chromosomes and synchronizing the sites. This is used by vcfcheck to
compare multiple VCFs simultaneously and is used also for merging,
creating intersections, etc.
The synced_bcf_reader also provides API for reading indexed BCF/VCF,
hiding differences in BCF/VCF opening, indexing and reading.
Example of usage:
bcf_srs_t *sr = bcf_sr_init();
for (i=0; i<nfiles; i++)
bcf_sr_add_reader(sr,files[i]);
while ( bcf_sr_next_line(sr) )
{
for (i=0; i<nfiles; i++)
{
bcf1_t *line = bcf_sr_get_line(sr,i);
...
}
}
bcf_sr_destroy(sr);
*/
#ifndef SYNCED_BCF_READER_H
#define SYNCED_BCF_READER_H
#include "hts.h"
#include "vcf.h"
#include "tbx.h"
// How should be treated sites with the same position but different alleles
#define COLLAPSE_NONE 0 // require the exact same set of alleles in all files
#define COLLAPSE_SNPS 1 // allow different alleles, as long as they all are SNPs
#define COLLAPSE_INDELS 2 // the same as above, but with indels
#define COLLAPSE_ANY 4 // any combination of alleles can be returned by bcf_sr_next_line()
#define COLLAPSE_SOME 8 // at least some of the ALTs must match
#define COLLAPSE_BOTH (COLLAPSE_SNPS|COLLAPSE_INDELS)
typedef struct
{
// for reading from tabix-indexed file (big data)
tbx_t *tbx; // tabix index
hts_itr_t *itr; // tabix iterator
kstring_t line; // holder of the current line
htsFile *file;
char *fname;
int is_bin; // is open in binary mode (tabix access)
char **als; // parsed alleles if targets_als set and _regions_match_alleles called
int nals, mals; // number of set alleles and the size of allocated array
int als_type; // alleles type, currently VCF_SNP or VCF_INDEL
// for in-memory regions (small data)
struct _region_t *regs; // the regions
int nregs, ireg; // number of regions; current region
// shared by both tabix-index and in-memory regions
char **snames; // sequence names
int nseqs; // number of sequences (chromosomes) in the file
char *seq; // current position: chr name
int start, end; // current position: start, end of the region (0-based)
int done; // finished the region or the whole file if streaming
}
bcf_sr_regions_t;
typedef struct
{
htsFile *file;
tbx_t *tbx_idx;
hts_idx_t *bcf_idx;
bcf_hdr_t *header;
hts_itr_t *itr;
const char *fname;
bcf1_t **buffer; // cached VCF records. First is the current record synced across the reader
int nbuffer, mbuffer; // number of cached records (including the current record); number of allocated records
int nfilter_ids, *filter_ids; // -1 for ".", otherwise filter id as returned by bcf_id2int
int type;
int *samples, n_smpl; // list of columns in the order consistent with bcf_srs_t.samples
}
bcf_sr_t;
typedef struct
{
// Parameters controlling the logic
int collapse; // How should the duplicate sites be treated. One of the COLLAPSE_* types above.
char *apply_filters; // If set, sites where none of the FILTER strings is listed
// will be skipped. Active only at the time of
// initialization, that is during the add_reader()
// calls. Therefore, each reader can be initialized with different
// filters.
int require_index; // Some tools do not need random access
int max_unpack; // When reading VCFs and knowing some fields will not be needed, boost performance of vcf_parse1
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.
// Auxiliary data
bcf_sr_t *readers;
int nreaders;
int streaming; // reading mode: index-jumping or streaming
int explicit_regs; // was the list of regions se by bcf_sr_set_regions or guessed from tabix index?
char **samples; // List of samples
bcf_sr_regions_t *regions, *targets; // see bcf_sr_set_[targets|regions] for description
int targets_als; // subset to targets not only by position but also by alleles? (todo)
char *cseq; // current sequence for targets, when regions are present
int crid; // current sequence for targets, with single VCF
kstring_t tmps;
int n_smpl;
}
bcf_srs_t;
/** Init bcf_srs_t struct */
bcf_srs_t *bcf_sr_init();
/** Destroy bcf_srs_t struct */
void bcf_sr_destroy(bcf_srs_t *readers);
/**
* bcf_sr_add_reader() - open new reader
* @readers: holder of the open readers
* @fname: the VCF file
*
* Returns 1 if the call succeeded, or 0 on error.
*
* See also the bcf_srs_t data structure for parameters controlling
* the reader's logic.
*/
int bcf_sr_add_reader(bcf_srs_t *readers, const char *fname);
/**
* bcf_sr_next_line() - the iterator
* @readers: holder of the open readers
*
* Returns the number of readers which have the current line
* (bcf_sr_t.buffer[0]) set at this position. Use the bcf_sr_has_line macro to
* determine which of the readers are set.
*/
int bcf_sr_next_line(bcf_srs_t *readers);
#define bcf_sr_has_line(readers, i) (readers)->has_line[i]
#define bcf_sr_get_line(_readers, i) ((_readers)->has_line[i] ? ((_readers)->readers[i].buffer[0]) : NULL)
/**
* bcf_sr_set_samples() - sets active samples
* @readers: holder of the open readers
* @samples: this can be one of: file name with one sample per line;
* or column-separated list of samples; or '-' for a list of
* samples shared by all files. If first character is the
* exclamation mark, all but the listed samples are included.
*
* Returns 1 if the call succeeded, or 0 on error.
*/
int bcf_sr_set_samples(bcf_srs_t *readers, const char *samples);
/**
* bcf_sr_set_targets(), bcf_sr_set_regions() - init targets/regions
* @readers: holder of the open readers
* @targets: regions can be either a comma-separated list of regions
* (chr,chr:from-to) or a name of a tabix indexed file with a list
* of regions (<chr,pos> or <chr,from,to>). Coordinates are
* one-based and inclusive.
*
* Returns 0 if the call succeeded, or -1 on error.
*
* Both functions behave the same way, unlisted positions will be skipped by
* bcf_sr_next_line(). However, there is an important difference: regions use
* index to jump to desired positions while targets streams the whole files
* and merely skip unlisted positions.
*
* Moreover, bcf_sr_set_targets() accepts an optional parameter $alleles which
* is intepreted as a 1-based column index in the tab-delimited file where
* alleles are listed. This in principle enables to perform the COLLAPSE_*
* logic also with tab-delimited files. However, the current implementation
* considers the alleles merely as a suggestion for prioritizing one of possibly
* duplicate VCF lines. It is up to the caller to examine targets->als if
* perfect match is sought after. Note that the duplicate positions in targets
* file are currently not supported.
*/
int bcf_sr_set_targets(bcf_srs_t *readers, const char *targets, int alleles);
int bcf_sr_set_regions(bcf_srs_t *readers, const char *regions);
/*
* bcf_sr_regions_init()
* @regions: regions can be either a comma-separated list of regions
* (chr,chr:from-to) or a name of a tabix indexed file with a list
* of regions (<chr,pos> or <chr,from,to>). Coordinates are
* one-based and inclusive.
*/
bcf_sr_regions_t *bcf_sr_regions_init(const char *regions);
void bcf_sr_regions_destroy(bcf_sr_regions_t *reg);
/*
* bcf_sr_regions_seek() - seek to the chromosome block
*
* Returns 0 on success or -1 on failure. Sets reg->seq appropriately and
* reg->start,reg->end to -1.
*/
int bcf_sr_regions_seek(bcf_sr_regions_t *regions, const char *chr);
/*
* bcf_sr_regions_next() - retrieves next region. Returns 0 on success and -1
* when all regions have been read. The fields reg->seq, reg->start and
* reg->end are filled with the genomic coordinates on succes or with
* NULL,-1,-1 when no region is available. The coordinates are 0-based,
* inclusive.
*/
int bcf_sr_regions_next(bcf_sr_regions_t *reg);
/*
* bcf_sr_regions_query() - checks if the interval <start,end> overlaps any of
* the regions. The sequence name is set by bcf_sr_regions_seek(). The
* coordinate queries must come in ascending order.
*
* Returns 0 if the position is in regions; -1 if the position is not in the
* regions and more regions exist; -2 if not in the regions and there are no more
* regions left.
*/
int bcf_sr_regions_query(bcf_sr_regions_t *reg, int start, int end);
#endif
|