bfilereader
(short for big file reader) is a MATLAB tool for reading and parsing big delimited files (can also be GNU zip gz
compressed) in a fast and efficient manner. bfilereader
takes advantages of a dependent Java class (bFileReaderDep.class
) with multiple methods implemented for reading, mapping and filtering big delimited files.
There are few requirements to be met prior to use bfilereader
:
- MATLAB version must be R2019b or newer.
bFileReaderDep
was writtent and compiled in Java 8. For more information see MATLAB documents on Java Environment.- To avoid Java out of memory issues when working with large files, configure and increase your MATLAB Java Heap Memory Preferences.
Add the both bfilereader.m
and bFileReaderDep.class
to your MATLAB path.
addpath(pwd);
savepath;
bfilereader
can be used for different purposes:
- Reading the whole content or multiple columns of a delimited file (similar behavior but limited functionality to MATLAB
readtable
). - Text pattern matching (with regular expression).
- Filtering data of numeric type based on different criteria.
Different uses of bfilereader
have been covered in the following examples. The working file for these examples is publicly available GWAS summary statistics on iron metabolism disorder which can be freely download in GZIP format (phenocode-275.1.tsv.gz
).
info = dir('phenocode-275.1.tsv.gz');
fprintf('file size: %.2f mb\n', info.bytes/1e6)
file size: 641.33 mb
To see what sort of data this file contains, we can simply call the function with summary
option set to only
to fetch only first few lines. It's also useful to see the number of rows by setting verbose
flag to on
.
out = bfilereader('phenocode-275.1.tsv.gz', 'summary', 'only', 'verbose', 'on');
Elapsed time is 11.828086 seconds.
file has 28336915 lines and 13 columns
file first 6 rows:
Var1 Var2 Var3 Var4 Var5 Var6 Var7 Var8 Var9 Var10 Var11 Var12 Var13
_______ _______ _____ _____ _____________ _______________ _______________________ ______ _______ ________ _________ __________ ________
"chrom" "pos" "ref" "alt" "rsids" "nearest_genes" "consequence" "pval" "beta" "sebeta" "af" "ac" "tstat"
"1" "16071" "G" "A" "rs541172944" "OR4F5" "intron_variant" "0.71" "-2.8" "7.6" "5e-05" "40.9" "-0.048"
"1" "16280" "T" "C" "rs866639523" "OR4F5" "intron_variant" "0.56" "-2.3" "3.9" "0.00015" "125.6" "-0.15"
"1" "49298" "T" "C" "rs10399793" "OR4F5" "upstream_gene_variant" "0.66" "0.042" "0.097" "0.62" "508337.5" "-20.0"
"1" "54353" "C" "A" "rs140052487" "OR4F5" "intron_variant" "0.8" "0.81" "3.3" "0.00036" "289.1" "0.076"
"1" "54564" "G" "T" "rs558796213" "OR4F5" "intron_variant" "0.98" "0.091" "3.1" "0.00015" "121.2" "0.0095"
This file has ~28,300,000 rows and 13 columns, which may not fit into the memory. From the first row, we easily notice that the file has one header row (variable names). We can use this summary table to extract additional data.
Function can apply pattern matching to extract desired information. For instance, to extract variants on either PNPLA3 or GPAM genes (column nearest_genes
):
out = bfilereader('phenocode-275.1.tsv.gz', 'header', true, 'pattern', ["PNPLA3", "GPAM"], 'patternCol', "nearest_genes");
Elapse time: 28.334 sec
size(out)
ans =
7776 13
unique(out.nearest_genes)
ans =
3ร1 string array
"GPAM"
"PNPLA3"
"PNPLA3,SAMM50"
We can also use regular expression (regex) and further add another filter. This time, we want to find 1)all variants on PNPLA family genes and 2)missense variants on other genes. We also would like to only extract frist 10 columns (by using extractCol
option).
patterns = ["PNPLA*", "missense"];
cols = ["nearest_genes", "consequence"];
out = bfilereader('phenocode-275.1.tsv.gz', 'header', true, 'pattern', patterns, 'patternCol', cols, 'extractCol', 1:10);
Elapse time: 31.639 sec
size(out)
171838 10
head(out, 4)
chrom pos ref alt rsids nearest_genes consequence pval beta sebeta
_____ __________ ___ ___ _____________ _____________ __________________ _____ _____ ______
1 1.3903e+05 "G" "A" "rs751110858" "AL627309.1" "missense_variant" 0.79 -1.4 5.4
1 1.3906e+05 "G" "A" "rs568513188" "AL627309.1" "missense_variant" 0.039 4.6 2.2
1 7.3854e+05 "T" "C" "rs147999235" "OR4F16" "missense_variant" 0.95 0.045 0.78
1 8.6135e+05 "C" "T" "rs200686669" "SAMD11" "missense_variant" 0.7 -0.81 2.1
% check variants on PNPLA2 gene
pnpla2 = out(out.nearest_genes == "PNPLA2", :);
head(pnpla2, 4)
chrom pos ref alt rsids nearest_genes consequence pval beta sebeta
_____ __________ ___ ___ _____________ _____________ _______________________ ____ ______ ______
11 8.1591e+05 "C" "T" "rs191403268" "PNPLA2" "upstream_gene_variant" 0.94 -0.054 0.78
11 8.1602e+05 "A" "G" "rs11246321" "PNPLA2" "upstream_gene_variant" 0.33 0.2 0.2
11 8.1615e+05 "G" "A" "rs546654103" "PNPLA2" "upstream_gene_variant" 0.64 0.4 0.85
11 8.1619e+05 "C" "T" "rs755386980" "PNPLA2" "upstream_gene_variant" 0.78 -0.68 2.4
In above example, we fetched all missense variants except for PNPLA family, for which we extracted all variants regardless of their consequences. But, what if we were only interested in finding missense variants on PNPLA family genes? To do so, we need to set multiCol
flag to true
to tell bfilereader
to apply each pattern to its corresponding column (i.e. pattern 1 to column 1, pattern 2 to column 2, ...):
out = bfilereader('phenocode-275.1.tsv.gz', 'header', true, 'pattern', patterns, 'patternCol', cols,...
'extractCol', 1:10, 'multiCol', true);
Elapse time: 40.223 sec
size(out)
107 10
head(out, 4)
chrom pos ref alt rsids nearest_genes consequence pval beta sebeta
_____ __________ ___ ___ _____________ _____________ __________________ ____ ______ ______
6 3.6259e+07 "C" "T" "rs140585347" "PNPLA1" "missense_variant" 0.7 -0.95 2.5
6 3.6262e+07 "G" "A" "rs74946910" "PNPLA1" "missense_variant" 0.82 -1.5 6.6
6 3.6263e+07 "G" "A" "rs45524833" "PNPLA1" "missense_variant" 0.95 -0.029 0.49
6 3.627e+07 "A" "G" "rs371888522" "PNPLA1" "missense_variant" 0.61 -1.2 2.3
unique(out.nearest_genes)
"AC008878.2,PNPLA6"
"PNPLA1"
"PNPLA2"
"PNPLA3"
"PNPLA3,SAMM50"
"PNPLA5"
"PNPLA6"
"PNPLA7"
"PNPLA8"
"PNPLA8,THAP5"
In addition to pattern matching for strings, data of numeric type can be filtered as well. For instance, if we would like to only see how many variants pass genome-wide significance threshold (5e-8
), we can set filter
and filterCol
options:
out = bfilereader('phenocode-275.1.tsv.gz', 'header', true, 'extractCol', 1:10,...
'filter', 5e-8, 'filterCol', "pval", 'operator', '<='); % any pval <= 5e-8
Elapse time: 26.312 sec
size(out)
10598 10
fprintf('pval range: %.3g - %.3g\n', min(out.pval), max(out.pval))
pval range: 0 - 5e-08
Similar to pattern matching, we can filter for other columns. However, we don't need to set multiCol
option in this case since every filtering value can only be applied to one column. This time we want to find variants passing genome-wide significance threshold and having a negative effect size (beta):
out = bfilereader('phenocode-275.1.tsv.gz', 'header', true, 'extractCol', 1:10, ...
'filter', [5e-8, 0], 'filterCol', ["pval", "beta"], 'operator', '<='); % any pval <= 5e-8
Elapse time: 27.426 sec
size(out)
7806 10
fprintf('pval range: %.3g - %.3g\n', min(out.pval), max(out.pval))
fprintf('beta range: %.3g to %.3g\n', min(out.beta), max(out.beta))
pval range: 1.5e-197 - 5e-08
beta range: -2.9 to -0.3
We can also use a mixture of examples 2 and 3 by including both filtering and pattern options. Suppose we want to get missense variants on MHC class I genes passing 5e-8 threshold and having a negative beta:
patterns = ["HLA-\w{1}$", "missense"];
cols = ["nearest_genes", "consequence"];
out = bfilereader('phenocode-275.1.tsv.gz', 'header', true, 'extractCol', 1:10, 'filter', [5e-8, 0], 'filterCol',...
["pval", "beta"], 'operator', "<=", 'pattern', patterns, 'patternCol', cols, 'multiCol', true);
Elapse time: 41.277 sec
disp(out)
chrom pos ref alt rsids nearest_genes consequence pval beta sebeta
_____ __________ ___ ___ ___________ _____________ __________________ _______ _____ ______
6 2.9692e+07 "C" "G" "rs2072895" "HLA-F" "missense_variant" 5.5e-09 -0.32 0.056
6 2.9693e+07 "C" "T" "rs1736924" "HLA-F" "missense_variant" 4.3e-52 -1.2 0.08
6 2.991e+07 "C" "G" "rs1143146" "HLA-A" "missense_variant" 6.5e-18 -0.48 0.056
6 2.9911e+07 "T" "G" "rs1059542" "HLA-A" "missense_variant" 1.2e-66 -1.5 0.088
6 3.0458e+07 "G" "A" "rs1264457" "HLA-E" "missense_variant" 4.4e-08 -0.31 0.056
bfilereader
can also parse the input file in parallel; however, whether using this option positively or negatively influences the performance depends on several factors (size of file in hand, the available memory, memory overhead and more). For a good discussion, see here. To show how it may affect the file processing, we consider 3 scenarios using different delimited files and we compare sequential and parallel bfilereader
with MATLAB tall datastore. Under all scenarios, we will use only uncompressed files.
We begin with a similar but smaller GWAS summary statistics file.
file = "phenocode-286.81.tsv";
fprintf('file size: %.2f mb\n', dir(file).bytes/1e6)
ile size: 490.55 mb
patt = ["LDLR$", "missense"];
col = ["nearest_genes", "consequence"];
out = bfilereader(file, 'header', true, 'pattern', patt, 'patternCol', col, 'multiCol', true); % sequential
out = bfilereader(file, 'header', true, 'pattern', patt, 'patternCol', col, 'multiCol', true, 'parallel', true); % parallel
% check with MATLAB tall datastore
parpool('local', 8); % max available cores
ds = tabularTextDatastore(file, 'FileExtensions', '.tsv', 'TextType', 'string');
ds.SelectedFormats{1} = '%q'; % for chromosome "X"
ds = tall(ds);
idx = endsWith(ds.(col(1)), "LDLR") & contains(ds.(col(2)), patt(2)); % regexp and contains cannot be applied to tall arrays
out2 = gather(ds(idx, :));
% benchmark: mean elapsed time over 3 repetitions
tbench.sequential = [7.3670 7.2930 7.3010];
tbench.parallel = [3.7540 3.7330 3.8020];
tbench.tall = [5.6000 5.4000 5.2000];
disp(table(structfun(@mean, tbench), 'VariableNames', {'mean elapsed time'}, 'RowNames', fieldnames(tbench)))
mean elapsed time
_________________
sequential 7.3203
parallel 3.763
tall 5.4
Next, we use the same delimited but uncompressed file (phenocode-275.1.tsv
) we used in examples above. This file is relatively bigger (~4.6 times) than the file we used in scenario 1.
file = "phenocode-275.1.tsv";
fprintf('file size: %.2f GB\n', dir(file).bytes/1e9)
file size: 2.43 GB
patterns = ["HLA-", "missense"];
cols = ["nearest_genes", "consequence"];
out = bfilereader(file, 'header', true, 'extractCol', 1:10, 'filter', [5e-8, 0], 'filterCol',["pval", "beta"], 'operator', "<=",...
'pattern', patterns, 'patternCol', cols, 'multiCol', true);
out = bfilereader(file, 'header', true, 'extractCol', 1:10, 'filter', [5e-8, 0], 'filterCol',["pval", "beta"], 'operator', "<=",...
'pattern', patterns, 'patternCol', cols, 'multiCol', true, 'parallel', true);
ds = tabularTextDatastore(file, 'FileExtensions', '.tsv', 'TextType', 'string');
ds.SelectedFormats{1} = '%q';
tt = tall(ds);
idx = startsWith(tt.(cols(1)), patterns(1)) & startsWith(tt.(cols(2)), patterns(2)) & tt.pval <= 5e-8 & tt.beta <= 0;
out2 = gather(ds(idx, :));
% benchmark: mean elapsed time over 3 repetitions
tbench.sequential = [30.593, 30.447, 30.371];
tbench.parallel = [19.740 19.907 20.116];
tbench.tall = [23.361, 25.508, 24.617];
disp(table(structfun(@mean, tbench), 'VariableNames', {'mean elapsed time'}, 'RowNames', fieldnames(tbench)))
mean elapsed time
_________________
sequential 30.47
parallel 19.921
tall 24.495
Lastly, we use a much bigger file (I used chromosome 1 from dbNSFP project version 4.1). In this case, parallel
would throw out of memory error, so we only benchmark with sequential bfilereader
and tall datastore.
file = "dbNSFP4.1a_variant.chr1";
fprintf('file size: %.2f GB\n', dir(file).bytes/1e9)
file size: 18.44 GB
filterCol = "CADD_phred";
patternCol = "genename";
patt = "^MARC1$";
filter = 15;
out = bfilereader(file, 'header', true, 'pattern', patt, 'patternCol', patternCol,...
'filter', filter, 'filterCol', filterCol, 'operator', '>=');
parpool('local', 8); % max available cores
ds = tabularTextDatastore(file, 'FileExtensions', '.chr1', 'TextType', 'string', 'TreatAsMissing', {'.', '-'});
ds.SelectedFormats = repmat({'%q'}, 1, numel(ds.SelectedFormats)); % to avoid conversion to double error
ds.SelectedFormats(ismember(ds.SelectedVariableNames, filterCol)) = {'%f'};
tt = tall(ds);
idx = ismember(tt.(patternCol), "MARC1") & tt.(filterCol) >= filter;
out2 = gather(ds(idx, :));
% benchmark: mean elapsed time over 3 repetitions
tbench.sequential = [113.885, 113.662, 113.943];
tbench.tall = [389.2, 392.66 405.336];
disp(table(structfun(@mean, tbench), 'VariableNames', {'mean elapsed time'}, 'RowNames', fieldnames(tbench)))
mean elapsed time
_________________
sequential 113.83
tall 395.73
In scenarios 1 and 2, parallel bfilereader
performed better than both MATLAB tall datastore and sequential bfilereader
. However, when file does not fit into the memory like the case in scenario 3, memory overhead can be a serious issue. Under such circumstances, parallel computing can negatively affect the performance. Therefore, don't apply parallel
flag just because it seems cool! proper benchmarking can show if your task really benefits from parallel computing or not.