Using AWK and R to parse 25TB
livefreeordichotomize.com
livefreeordichotomize.com
https://adamdrake.com/command-line-tools-can-be-235x-faster-...
Right tool for the right job, as always. For a 2-3GB dataset size you don't need to bother with Hadoop just as for a 2-3PB dataset size you probably don't need to bother with awk.
In OP's article there is a link to a book "Data Science at the Command Line" which sounds quite relevant: https://www.datascienceatthecommandline.com
[1] https://github.com/datascienceworkshops/dockerfiles/blob/mas...
One of the most eye-opening aspects of awk (goes for other pipeable commandline tools too), was how they support iterative development of pipelines regardless of data size.
We tried SQLite at some point for some of the stages because of pretty complicated selections sometimes, but it often won't give back a single result in minutes. Switching to AWK, I could immediately get some output, so I could quickly validate and iterate on the awk code until I got what I expected. The actual execution was likeways always very fast.
I'm not sure what your sample was but seems like it would have been better to use one of the special binary file formats for genetic data. You wrote SNP chips, But in order to get to 25 Tb I assume there must be imputed calls, so it seems like a BGEN might have been a lot easier.
This is speculation of course, I'm not sure exactly what your situation was.
My main practical recommendation for Nick is to become familiar with bgzip and zstd. bgzip sacrifices a little bit of compression efficiency relative to plain gzip, but in exchange it solves the more important problems of (i) letting you take advantage of all your cores when decompressing and (ii) supporting random-access reads with an appropriate index, while remaining compatible with all .gz-reading programs. When backward compatibility is unimportant, zstd tends to have much better compression/decompression speed for the same compression ratio than gzip.
In the end CPUs are fast enough and compression algorithms good enough that I would argue it doesn't really matter what format you use for storage, as long as it's reasonably easy to read back.
I work as a data analyst and I never got to worry about big data as the DWH takes care of the aggregation for us, plus I only work in Windows.
I see now that it would be very useful to learn *nix tools in general, as it seems that the skills to process (not to predict/analyze) terabytes+ data are very valuable and expensive to acquire and could be one's butter and bread.
https://github.com/bigdatagenomics/adam
We've fallen off the first Google hit the past few years but are still quite relevant (e.g. Databricks' commercial offering uses ADAM under the hood). Drop in our Gitter some time!
You would have to add some additional reliability to your pipeline so it could continue processing when instances are unexpectedly terminated. But this might be well worth it as it sounds like your research group is cost-constrained.
AWS made some changes this year so the spot prices are more stable and instances don't get shut down as frequently.
But there were bugs in mawk and it seemed basically unmaintained. So you'd run into something and have to use gawk or perl instead.
That's no longer the case, the xterm guy adopted it, ten years ago, and now I know!
Mawk is my go-to version because of speed. GNU Awk when its extensions are needed, or the task is over "small data" and the system default version is sufficient.
Do you know why?
Initially I was trying to use SQLite for it but I kept running out of memory and crashing the system. Turned out using grep, join, sort, and paste got the job done in seconds.
I think using 'make' with the -j parameter (# of parallel jobs) is more useful than using gnu 'parallel'. The reason is that if one of the job fails for some reason, you just re-run 'make' and only the required jobs are started instead of restarting the entire computation.
A couple of quick questions:
Was the 25TB raw data gathered from a single human genome?
What would be the size in bytes of a unique genomic fingerprint once raw data is all fully processed into high confidence base values? (including non-coding regions)
If we just look at coding regions and further compress by only looking at SNPs, how many bytes is that?
Considering that each base has ~2B of information... it would be super interesting to know how much space it takes to describe our uniqueness!
Essay: "Each row contained a data for a single SNP for a single person." ... "There were ~2.5 million SNPS and ~60 thousand people"
Statement: "it would be super interesting to know how much space it takes to describe our uniqueness"
Genome size is "3,234.83 Mb (Mega-basepairs) per haploid genome" says https://en.wikipedia.org/wiki/Human_genome .
However, question of uniqueness depends on your model. We are all unique.
If you have human genomes, then use a reference template and only list those which differ. Eg, https://en.wikipedia.org/wiki/Compression_of_Genomic_Sequenc... .
Describing your total unique genetic profile would obviously require a lot more space, and wouldn't be constant across individuals/ancestral backgrounds (e.g. there's more genetic diversity in people of African descent).
BTW not sure, but is it OK to post other's article here? Maybe I should add a short commentary in the title.
Just counterparts to all the favourites that utilise the GPU ... imagine GPU awk.
My intuition tells me that awk and other text processing tools won’t scale well to a GPGPU. I might be wrong though. Is there any example of something like grep etc working well on a GPU?
https://www.cs.cmu.edu/afs/cs/academic/class/15418-s12/www/c...
I would have expected a GPU regex too perform much worse, given that regex matching is probably very branchy code. Especially since computation is generally way faster than IO.
What's at play here is that the needle in a haystack search of regex is going to spend almost all its time 0 or 1 deep in the state machine, so the threads skip the branches and the penalty is not large.
But as to your question, a lot of traditional tools like grep, sed, aren't really suited for the GPU unless you are running them on a lot of files at once.
AMA: Explaining my 750 line compiler+runtime designed to GPU self-host APL (youtube.com) https://news.ycombinator.com/item?id=13797797
If you only need a single SNP or a group of SNPs within a region, you can use tabix[1] to index gzipped TSVs and query by genomic position. The position of SNPs can be obtained from a lookup table (2.5M is not very big even for R) or from an API if you were to say - query by rsid.
tabix also works over http (and s3) and can utilize RANGE queries to select a subset of a file...so you only wind up downloading a or reading a small portion once it is indexed and can do something like this:
tabix https://www.file.url.tsv chr1:1-1000
The command above would return variants on chromosome 1 between 1 and 1000.
The following variant browser works in this way: https://elegansvariation.org/data/browser/ - Theres no formal database (e.g. MySQL) running here, just tabix (actually bcftools which uses tabix) to select variants in a particular region, wrap them in JSON, and return to the client.
Setting this up on S3 requires configuring CORS... the igv browser also uses tabix indexes and provides guidance on how to set this up [2]
[1] https://www.htslib.org/doc/tabix.html [2] https://github.com/igvteam/igv.js/wiki/Data-Server-Requireme...
I created a similar solution to what you have done using this alternative approach, writing a wrapper in R that invoked bcftools under the hood. The dataset I was working with was a lot smaller (1.6M Snps x 252 individuals), but should work with larger genotype sets as well.
Since I never dealt with big data before, I'm wondering what would you do in this situation?
If you have such an specific task, why not just write an "actual program" (as opposed to pipeline of scripts)? From the looks of it, it sounds like this problem could have been solved with, say, 50 lines of Java, C, Go, etc, etc. Maybe a bit more verbose but it would give you full control, you wouldn't need to lookup how to use command line parameters on S/O, and would probably give you a bit more performance.
Dealing with exceptions is extremely error prone if these exceptions are not mapped beforehand. Thus it can be very costly.
Similarly doing stuff in parallel is extremely error prone due to race conditions: What does not happen when running on your 1 GB test set, may very well happen when running on your 25 TB production data.
In a program I'd try/catch defensively "just in case", if missing one line out of 25TB is not a bit deal.
For parallel processing I'd reach for the nearest standard library at hand on the language of choice.
That is a good example of what I mean: The nearest standard library is likely to either buffer output in memory or not buffer at all (in which case you can have the start of one line ending with another line). This means you cannot deal with output bigger than physical RAM. And your test set will often be so small that this problem will not show up.
GNU Parallel buffers on disk. It checks whether the disk runs full during a run and exits with a failure if that happens. It also removes the temporary files immediately, so if GNU Parallel is killed, you do not have to clean up any mess left behind.
You _could_ do all that yourself, but then we are not talking 50 lines of code. Parallelizing is hard to get right for all the corner cases - even with a standard library.
And while you would not have to look up how to use command line parameters on S/O you _would_ be doing exactly the same for the standard libraries.
Assuming you can get better performance is also not given: GNU Sort has built-in parallel sorting. So you clearly would not want to use a standard non-parallelized sort.
Basically I see you have 2 choices: Built it yourself from libraries, or build it as a shell script from commands.
You would have to spend time understanding how to use the libraries and the commands in both cases, and you are limited by whatever the library or the command can do in both cases.
I agree that if you need tighter control than a shell script will give you, then you need to switch to another language.
With the pipeline you get free parallelism and it's much easier to iterate over the individual steps, run step manually, check the output, add it to the script. You can also trivially break bits into a make file for improved parellelism and incrementalism. Performance wise these tools have had a lot of work put into them, even when they're not the most efficient tool on paper they'll often beat out the most naive versions in "real" languages.
Mastering the usual command line interface (terminal emulator, interactive shell, maybe a terminal multiplexer) is non-optional if you want to use CLI tools at or close to peak effectiveness.
One very widely used option is Galaxy https://galaxyproject.org
You might also want to check out:
- Arvados https://arvados.org
- Chipster https://chipster.csc.fi
- Knime https://www.knime.com
They each have their strenghts and weaknesses.
I'd probably be going with Galaxy as the default option, unless you find specific reasons to go with some of the others.
Task-splitting itself is inherently recoverable, as incomplete work units don't produce the final output you're looking for, and can. be retried, re-run, or re-entered into the task pool.
A GUI quite frankly simply gets in the way.
paste - merge lines of files1. I wrote a plain SQL mapper program, which spawns a worker-thread pool, where each worker opens its own "part" file for each SQL table, such that a document consumed by worker N gets records written to "tableA/partN.pgcopy".
2. And then, after the mapper is done, to do the reduce step, I just spawned a `sort -m -u -k n1` invocation to collate the row-files of each table together into a single `.sql` file. This not only efficiently merge-sorts the (presorted) "part" files into one file (without needing to sort the files themselves), but also blows away any rows with duplicate primary-keys [i.e. duplicate first columns in the row's TSV representation]—meaning I can restart the mapper program in the middle of a job (causing it to create a new set of "parts") and sort(1) will take care of cleaning up the result!
I honestly don't think anything could be done to make this system more optimal for its use-case. sort(1) goes crazy fast when it can mmap(2) files on disk.
(Also, I'm pretty sure that even the framework part of the mapper program—the "N persistent workers that each greedily consume a document-at-a-time from a shared pipe, as if they were accept(2)ing connections on a shared socket"—could be created with Unix tooling as well, though I'm not sure how. GNU parallel(1), maybe?)
Bonus: once you have SQL rowsets in TSV form like this, you can calculate a "differential" rowset (against a rowset you've already inserted) using `comm -23 $new $old`. No need for a staging table in your data warehouse; you can dedup your data-migrations at the source.
On a mid-1990s iMac running Debian, and a rusty residential DSL connection.
R played a role in other community-related analysis and reporting.
1) Assuming there is no joins/merges requirement, read in chunks and output GB dumps.
2) If joins/merges are required, use external merge sort.
Is this correct? Actually I'm wondering whether I could earn some bread and butter by focusing on the big data processing problems (e.g. sort/filter Terabytes+ dumps, do transformation for each line for Terabytes+ dumps, those kind of things) without actually knowing how to implement math algorithms (required for data science).
If so what kind of tools I need to master? I'm thinking about basic *nix tools like mentioned above, and also Python and maybe some compiled language for optimization (someone managed to speed up a Python external merge algorithm on 500GB file by 50% by implementing in Go), then maybe some easy algorithms (merge join, heap, etc.)
http://shop.oreilly.com/product/9780596003302.do
Python is slower compared to some of it’s compiled cousins, but it’s quick to write and a great skill to have when bash scripting can’t handle some of the complexity or you need dB access. We use it sometimes to call c programs to do DNA sequence alignments and process the returns.
I can support “misguided” Big Data projects from a political perspective if they help fund fixing the fundamental problems (similar to Agile consultants) that plague an organization, but most consultants are not going to do very well by suggesting going back and fixing something unrelated to their core value proposition itself. For example, if you hire a bunch of machine learning engineers and they all say “we need to spend months or even years cleaning up and tagging your completely unstructured data slop because nothing we have can work without clean data” that’ll probably frustrate the people paying them $1MM+ / year each to get some results ASAP. The basics are missing by default and it’s why the non-tech companies are falling further and further behind despite massive investments in technology - technology is not a silver bullet to crippling organizational and business problems (this is pretty much the TL;DR of 15+ years of “devops” for me at least).
In the end it doesn't matter if you wind up with a multi-TB copy of some large database or a handful of small XML files - it's all in one place, it gets updated, there are usable ACL in place, and it can be accessed and worked with. That's the point where you think about running a Spark job or the above AWK magic.
I would go further and even call long or at least not instant report generation a perceived feature. Similar to flight and hotel booking sites that show some kind of loading screen even if they could give instant search results, the duration of the generation itself seems to add trust to the reports.
Absolutely. I really want to see advanced AI/ML tools developed to address THIS problem. Don’t make me solve the data before I use ML, give me ML to fix my data!
That’s hard though, because data chaos is unbounded and computers are still dumb. I think there’s still tons of room for improvement though.
Well - Alex Ratner [stanford], for one: https://ajratner.github.io/
And several of Chris Re's other students have as well: https://cs.stanford.edu/~chrismre/
Trifacta is Joseph Hellerstein's [berkeley] startup for data wrangling: https://www.trifacta.com/
Sanjay Krishnan [berkeley]: http://sanjayk.io/
There are a lot of people tackling bits and pieces of the problem. Tom Mitchell's NELL project was an early one, using the web in all its messy glory...http://rtw.ml.cmu.edu/rtw/
Lots of other folks here (CMU). Particularly if you add an active learning. Hard messy problem that crosses databases and ML.
https://developer.ibm.com/tutorials/l-awk1/
From 2001.
https://developer.ibm.com/tutorials/l-awk2/ https://developer.ibm.com/tutorials/l-awk3/
These were mentioned but not linked to in the previous form of the article/blog, had a quick look at the newer version you linked to and that may still be the case.
A lot of the Awk info I had found prior to stumbling on these articles was focused on command line one-liners. So the sections on defining Awk scripts as files and multiline records were a great help to me.