text:
What's the most efficient file format for the storage of DNA sequences?
I'd like to learn which format is most commonly used for storing the full human genome sequence (4 letters without a quality score) and why.
I assume that storing it in plain-text format would be very inefficient. I expect a binary format would be more appropriate (e.g. 2 bits per nucleotide).
Which format is most common in terms of space efficiency?

see: https://www.biostars.org/p/75178/ "Why Don't We Use Binary Format?"

Also important question: is the goal to create the smallest footprint on disk for an isolated single genome, or a multiple genomes?

@GWW Even if you had 5 letters (i.e. with N) you could get away with 3 bits per nucleotide and still have room for 3 more nucleotide encodings, maybe for like U, mC, hmC.

Genomes are commonly stored as either fasta files (.fa) or twoBit (.2bit) files. Fasta files store the entire sequence as text and are thus not particularly compressed.
twoBit files store each nucleotide in two bits and contain additional metadata that indicates where there's regions containing N (unknown) bases.
For more information, see the documentation on the twoBit format at the UCSC genome browser.
You can convert between twoBit and fasta format using the faToTwoBit and twoBitToFa utilities.
For the human genome, you can download it in either fasta or twoBit format here: http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/

The standard formats for storing sequence data are fasta and fastq.  Fasta is used if you only need the raw sequence data, fastq is used if you want to store the sequence data along with the quality information from base calling.  Each of these can be compressed using gzip or another standard compression algorithm.
Typically we want to keep the quality information along with the raw sequence data, but the quality information accounts for half the storage space required.  Some people have developed algorithms for lossy compression of the quality data that allow us to reduce the storage requirements.
If you are interested in storing variant calling data, the standard format for that is VCF.  VCF is useful if you want to store quality information of the variant calls, genomic positions, and any annotations you might have about the position.  VCFs can be compressed and indexed using bgzip and tabix.  Many tools require variant data to be compressed and indexed using these tools.

fastq is generally used for storing read information garnered from sequencing. The information in a fastq is certainly not in chromosome or positional order and can be duplicated with multiple reads representing a position.

I can confirm with FASTA. Even Matlab's Bioinformatics Toolbox has a function (fastaread) for the import of such data. Plus, all of the genome data (at least that I have used) on NCBI is available in fasta.

I read the question as being about storing the reference sequence, not sequencing data.

Yes, this answer is about raw sequence data from sequencing projects which is a very inefficient format for storing stable, large sequences. FASTA is indeed very popular, but not FASTQ, not for things like genomes.

The standard and the most common sequence format is FASTA for sure. You can compress it with a compressor. For the ~3GB human genome, gzip reduces the size to ~900MB, depending on the option in use.
Another often used format is UCSC's 2-bit format. This format keeps each A/C/G/T with 2 bits. As I remember, it keeps non-A/C/G/T bases and lowercases in two separate lists. These lists basically tell you that bases between offset x and y are all "N"/lowercase. The 2-bit format loses IUB codes that GRCh37 has. UCSC's hg19 differs from GRCh37 at a few bases.
BWA also produces its own 2-bit format with indexing. You can generate it separately with:
bwa fa2pac -f hg19.fa

Unlike UCSC, BWA keeps all IUB codes but loses letter cases. BWA does not provide utilities to convert its 2-bit representation to FASTA, either.
The 2-bit format typically reduces the file size down to 1/4 of its original size, unless there are too many scattered ambiguous bases. For human genome, you get a file ~784MB in size. You can compress it further with gzip, but that actually doesn't work well. A gzip'd 2-bit file is only ~5-10% smaller.
If you want to achieve an even smaller file size, you can compress the BWT of 2-bit file. This gives you a ~633MB file:
bwa pac2bwtgen hg19.fa.pac tmp.bwt && gzip tmp.bwt

A bit-aware compression algorithm may achieve an even higher compression ratio. However, such BWT-based compression prevents you from extracting subsequences. It is probably of little use in practice.

BWA replaces ambiguous bases by random nucleotides. See the original BWA paper:
Section 2.7.1 ... "Non-A/C/G/T
bases on the reference genome are converted to random nucleotides. Doing
so may lead to false hits to regions full of ambiguous bases. Fortunately,
the chance that this may happen is very small given relatively long reads.
We tried 2 million 32 bp reads and did not see any reads mapped to poly-N
regions by chance."

@Karel All the ambiguous bases are kept in the .amb file. BWA can reconstruct every base, well, at least in principle.

Thank you very much for this information!! I have never used .amb files and they seem to be very useful for me. I wish I were aware of them earlier. Btw. I think we have some code to reconstruct the original sequences from BWA .bwt files. During our work on ProPhyle, we played a little bit with this type of compression. Maybe we will create a separate program for bwt2fa.

There are several things to consider when asking for "the most efficient" way to store data, it all depends on your use case. Do you just need ACGT, or are there also IUPAC codings for combinations? Do you need additional data (like quality values)? What kind of application are you using the data for (does it need to load all at once or in in chunks? Once or multiple times? Sequential or random access? etc.pp)?
E.g., most efficient for:

Lowest footprint on disk, without a lot of hassle: use either FASTA or 2bit, but run through standard compressor (gzip, bzip2, others). The literature you want to consult here is that of standard text compression I think. Also of interest Large Text Compression Benchmark
Keeping the file on disk, but ultra-fast loading small subsets into memory, being able to work in memory with character sized entities: a simple dump of the DNA as characters to disk, maybe combined with an index file to know which chromosome starts where. Then use mmap
Storing quality values: See papers like Compression of FASTQ and SAM Format Sequencing Data or Sequence squeeze: an open contest for sequence compression
Any combination of the above use cases + a lot more

I would classify 2bit format as 1b. Lowest footprint on disk, permitting some hassle. In order to use it for anything other than storage, it will need to be converted back into plaintext (fasta), or compressed plaintext (e.g. fasta.gz).

Your case 2 doesn’t require storing in uncompressed character data format. In fact, A compressed index (for instance) can make access faster by avoiding cache thrashing.

I think the question is a bit ambiguous so please excuse this answer that's a bit redundant from the rest of the ones provided.
As others have mentioned, if you want to store a full genome, FASTA and 2bit formats are appropriate.  For some context, hg19 is about 900Mb compressed for the FASTA file and about 780Mb compressed for the 2bit file .  hg19 is a reference and is haploid so doesn't represent a "full" human genome that would normally have two alleles for the autosome (non-sex chromosomes).
A common format for representing variant information is Variant Call Format (VCF). The VCF format represent differences from a reference (hg19, say) that can be used to recover the original full sequence by using the reference and the differences encoded in the VCF file.  I've seen VCF files in the range of 100Mb, but a reference file is still needed to recover the full genome sequence which is the range of 800Mb+, as mentioned above.
If you're considering just one "whole genome" in isolation, then the answer is pretty clear: 2bit format is probably approaching the entropy limit of the human genome and you probably won't be able to do much better.
The reason why your question is a bit ambiguous is that as soon as you start encoding more than one genome, a population of genomes, say, then you can start exploiting the redundancy of the genome as shared by the population.
For example, say you want to store two "whole genomes".  You could download the hg19 reference and download two VCF files which would give around 1Gb worth of data (around 800Mb for the 2bit file and around 200Mb for both of the VCF files).  Now you've been able to represent a "whole genome" in 500Mb instead of the 800Mb.  You can see a similar argument for downloading 3 VCF files and more.
The minimum amount of information needed to represent a population of genomes is, as far as I know, unknown, but I would guess in the 2.5Mb-5Mb range.  For example, see "Human genomes as email attachments" by Christley, Lu, Li and Xie which claims a 4Mb encoding of a genome.
Things get tricky because you have to ask what you're claiming as a "whole genome".  VCF files are notoriously bad because older versions of the specification only store high quality differences from reference, throwing away high quality called sections.  If you want to store low quality information, the encoding is now going to depend on the sequencing technology in weird ways.
Insertions, deletions, mobile insertion elements, copy number variants, other structural variants, etc. all complicate this matter further.  Genome Graphs are trying to tackle at least some of these problems but the focus is on variant calling rather than efficient individual whole genome representation, though perhaps can be adapted in the future.

How can there be 2 orders of magnitude differences between the 2bit compression (hundreds of Mb) and the "e-mail attachment" (a few Mb)? Is the "e-mail attachment" case a really impractical storage format so that no one actually uses it? The abstract doesn't say it, but it seems that what they store is actually variations than the full standalone info.

@bli, the order of magnitudes difference comes from exploiting the redundancy in a population. Storing one persons (encoded) DNA requires ~800Mb. Storing thousands of peoples takes (probably) a couple megabytes (each). When dealing with a population of DNA data that you want to encode, there are many ways of doing it.  One way is to store variants from a reference. Another is to store a library of short "reads" and then reference that library. The paper is meant as a proof of concept to answer the underlying question of "what is the theoretical minimum space required to store a whole genome".

It is not yet standardized, but graph format has the potential for being the most space-efficient method for storing genomes.  The idea is this: rather than store a genome as a linear string of sequenced nucleotides, genomes are stored as overlapping graphs, where sequence variants branch off from the reference genome, and then rejoin when the alignment continues.  Basically, you start with a reference genome, and for every subsequent genome added to the graph, only the differences are stored.  This could allow for an enormous gain in space efficiency.

In terms of raw storage capacity 2 bits per nucleotide, and then further compressed with standard compression techniques would be the most efficient. However, you'd still have other storage considerations. Like what to do about non-standard bases: like if you want to indicate a gap or ambiguity.
I'd also query if it is really necessary to store them as binary since it reduces the readability of the data. It is quite convenient having a whole bunch of unix and programming tools that can operate on the string level in text files.

In terms of storage space '2 bits per nucleotide' and compressed would be a way to save resources.

In all seriousness, the most efficient way to store DNA sequence data is...you guessed it...in DNA. (Church, Gao, and Kasuri, 2012) and others have used DNA synthesis and sequencing as an information writing/reading mechanism.
Practical? Not yet.
Storage efficiency? Unparalleled!

Genozip is lossless compression software that has been optimized for genomic data. If the benchmarks are to be believed, Illumina FASTQ files compressed with Genozip are nearly 5× smaller than FASTQ files compressed with gzip. Equally impressive gains over gzip are seen for FASTA, VCF, and GFF3 files. See the publication in Bioinformatics for more details:
Genozip: a universal extensible genomic data compressor

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/

sort: votes


text:
Accuracy of the original human DNA datasets sequenced by Human Genome Project?
The Human Genome Project was the project of 'determining the sequence of nucleotide base pairs that make up human DNA, and of identifying and mapping all of the genes of the human genome'. It was declared complete in 2003, i.e. 99% of the euchromatic human genome completed with 99.99% accuracy.
Are the datasets provided by HGP still accurate, or as accurate as was claimed in 2003?
Given the technology in the past (such as using old techniques), or any other reason (newer research studies), is it possible that the datasets are not as accurate as originally expected?

It should be noted that the reason why a technology is replaced by a newer one is not always because the old one was not accurate. Sometimes the newer technology can be a bit less accurate on some aspects, but so much cheaper that it is deemed worth switching. The sequencing and assembly techniques used to generate the data at the basis of the first release of the human genome may have been quite accurate. That said, using newer techniques to improve on this basis will of course make newer versions better.

The HGP developed the first "reference" human genome - a genome that other genomes could be compared to, and was actually a composite of multiple human genome sequences.
The standard human reference genome is actually continually updated with major and minor revisions, a bit like software. The latest major version is called GRCh38, was released in 2013, and has since had a number of minor updates.

Are the datasets provided by HGP still accurate?

Yes, in a sense, but we certainly have better information now. One way to measure the quality of the assembly is that the initial release from the HGP had hundreds of thousands of gaps - sequences that could not be resolved (this often occurs because of repetitive sequences). The newest reference genome has less than 500 gaps.

While the quality of the reference human assembly keeps improving, there are still misassemblies in it. A common problem is recent segmental duplications are occasionally collapsed into one sequence in the reference. Another issue is that the centromeric sequences in the reference are computationally generated, which are probably different from real sequences. Issues like these often complicate data analyses.
You should also beware that each human has a different genome. A large region having one copy in the reference genome may have two copies in a specific sample. While this is not really the problem with the reference, such copy-number changes will have the same effect as reference errors and mess up your pipelines.
There is still room for improvement to the human reference genome. In some regions, the CHM1 and CHM13 PacBio assemblies are better than the current reference genome at the larger scale. Illumina population data can produce better consensus at the single base level. GRC is continuously releasing new patches to the latest assembly.

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
Mapping drug names to ATC codes
I'm interested working with the medication information provided by the UK Biobank.  In order to get these into a usable form I would like to map them to ATC codes.  Since many of the drugs listed in the data showcase include dosage information, doing an exact string match between drug names is not very effective.  I've considered using something like fuzzywuzzy to do string matching between the medications in the data showcase and the ATC drug names but validating the matches could still be a laborious process.  Does anyone know of a tool that can match drug names to ATC codes or some other drug ontology?  If not, maybe there's a better way to do it that I haven't thought of.

Can you add a couple of examples for drug names and corresponding ATC codes (in case the link stops working)? In general, matching a non-controlled to a controlled vocabulary will always require manual checking.

The CART tool let's you upload a set of names and map them (optionally in a fuzzy way) to STITCH 4 identifiers, and then use those to map to ATC codes (using the chemicals sources download file). It's a bit indirect, and I'm not sure what CART will do with the dosage info you mention.

DrugBank seems to have a tool to map ATC codes to drug names and DrugBank IDs.
A quick look in the XSD schema on the release page suggests the complete database includes ATC codes for drugs, you could then do a fuzzy match of the BioBank names against all of DrugBank's synonyms, or match on some other data (e.g. canonicalised SMILES) if available.
The downside is that there may not be complete overlap between UK BioBank and DrugBank. Additionally, DrugBank is under licence for commercial use.


text:
What are the optimal parameters for docking a large ligand using Hex?
I'm looking to dock a large ligand (~90kDa) to a receptor slightly larger receptor  (~125kDa) using Hex. If anyone is familiar with docking large structures, are there any recommended parameters for finding the best docking solution?
Parameters in particular:

Number of Solutions
N order of correlation for initial and final
searches
Receptor Range
Ligand Range

While I have no experience with this specific question you have, according to Feinstein & Brylinski (2015), a "fully automated procedure", i.e. a Perl script, can be used to optimize the box size itself, and it can be found here.


text:
What is the difference between FASTA, FASTQ, and SAM file formats?
I'd like to learn the differences between 3 common formats such as FASTA, FASTQ and SAM. How they are different? Are there any benefits of using one over another?
Based on Wikipedia pages, I can't tell the differences between them.

Let’s start with what they have in common: All three formats store

sequence data, and
sequence metadata.

Furthermore, all three formats are text-based.
However, beyond that all three formats are different and serve different purposes.
Let’s start with the simplest format:
FASTA
FASTA stores a variable number of sequence records, and for each record it stores the sequence itself, and a sequence ID. Each record starts with a header line whose first character is >, followed by the sequence ID. The next lines of a record contain the actual sequence.
The Wikipedia artice gives several examples for peptide sequences, but since FASTQ and SAM are used exclusively (?) for nucleotide sequences, here’s a nucleotide example:
>Mus_musculus_tRNA-Ala-AGC-1-1 (chr13.trna34-AlaAGC)
GGGGGTGTAGCTCAGTGGTAGAGCGCGTGCTTAGCATGCACGAGGcCCTGGGTTCGATCC
CCAGCACCTCCA
>Mus_musculus_tRNA-Ala-AGC-10-1 (chr13.trna457-AlaAGC)
GGGGGATTAGCTCAAATGGTAGAGCGCTCGCTTAGCATGCAAGAGGtAGTGGGATCGATG
CCCACATCCTCCA

The ID can be in any arbitrary format, although several conventions exist.
In the context of nucleotide sequences, FASTA is mostly used to store reference data; that is, data extracted from a curated database; the above is adapted from GtRNAdb (a database of tRNA sequences).
FASTQ
FASTQ was conceived to solve a specific problem arising during sequencing: Due to how different sequencing technologies work, the confidence in each base call (that is, the estimated probability of having correctly identified a given nucleotide) varies. This is expressed in the Phred quality score. FASTA had no standardised way of encoding this. By contrast, a FASTQ record contains a sequence of quality scores for each nucleotide.
A FASTQ record has the following format:

A line starting with @, containing the sequence ID.
One or more lines that contain the sequence.
A new line starting with the character +, and being either empty or repeating the sequence ID.
One or more lines that contain the quality scores.

Here’s an example of a FASTQ file with two records:
@071112_SLXA-EAS1_s_7:5:1:817:345
GGGTGATGGCCGCTGCCGATGGCGTC
AAATCCCACC
+
IIIIIIIIIIIIIIIIIIIIIIIIII
IIII9IG9IC
@071112_SLXA-EAS1_s_7:5:1:801:338
GTTCAGGGATACGACGTTTGTATTTTAAGAATCTGA
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII6IBI

FASTQ files are mostly used to store short-read data from high-throughput sequencing experiments. The sequence and quality scores are usually put into a single line each, and indeed many tools assume that each record in a FASTQ file is exactly four lines long, even though this isn’t guaranteed.
As with FASTA, the format of the sequence ID isn’t standardised, but different producers of FASTQ use fixed notations that follow strict conventions.
SAM
SAM files are so complex that a complete description [PDF] takes 15 pages. So here’s the short version.
The original purpose of SAM files is to store mapping information for sequences from high-throughput sequencing. As a consequence, a SAM record needs to store more than just the sequence and its quality, it also needs to store information about where and how a sequence maps into the reference.
Unlike the previous formats, SAM is tab-based, and each record, consisting of either 11 or 12 fields, fills exactly one line. Here’s an example (tabs replaced by fixed-width spacing):
r001  99  chr1  7 30  17M         =  37  39  TTAGATAAAGGATACTG   IIIIIIIIIIIIIIIII
r002  0   chrX  9 30  3S6M1P1I4M  *  0   0   AAAAGATAAGGATA      IIIIIIIIII6IBI    NM:i:1

For a description of the individual fields, refer to the documentation. The relevant bit is this: SAM can express exactly the same information as FASTQ, plus, as mentioned, the mapping information. However, SAM is also used to store read data without mapping information.
In addition to sequence records, SAM files can also contain a header, which stores information about the reference that the sequences were mapped to, and the tool used to create the SAM file. Header information precede the sequence records, and consist of lines starting with @.
SAM itself is almost never used as a storage format; instead, files are stored in BAM format, which is a compact, gzipped, binary representation of SAM. It stores the same information, just more efficiently. And, in conjunction with a search index, allows fast retrieval of individual records from the middle of the file (= fast random access). BAM files are also much more compact than compressed FASTQ or FASTA files.

The above implies a hierarchy in what the formats can store: FASTA ⊂ FASTQ ⊂ SAM.
In a typical high-throughput analysis workflow, you will encounter all three file types:

FASTA to store the reference genome/transcriptome that the sequence fragments will be mapped to.
FASTQ to store the sequence fragments before mapping.
SAM/BAM to store the sequence fragments after mapping.

Why is there a '+' sign in FASTQ format?

@charlesdarwin I have no idea. The line with the plus sign is completely redundant. The original developers of the FASTQ format probably intended it as a redundancy to simplify error checking (= to see if the record was complete) but it fails at that. In hindsight it shouldn’t have been included. Unfortunately we’re stuck with it for now.

@KonradRudolph as far as I know fastq is a combination of fasta and qual files, see also https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2847217/
This explains the header of the quality part. It, however, doesn't make sense we're stuck with it...

In a nutshell,
FASTA file format is a DNA sequence format for specifying or representing DNA sequences and was first described by Pearson (Pearson,W.R. and Lipman,D.J. (1988) Improved tools for biological sequence comparison. Proc. Natl Acad. Sci. USA, 85, 2444–2448)
FASTQ is another DNA sequence file format that extends the FASTA format with the ability to store the sequence quality. The quality scores are often represented in  ASCII characters which correspond to a phred score)
Both FASTA and FASTQ are common sequence representation formats and have emerged as key data interchange formats for molecular biology and bioinformatics.
SAM is format for representing sequence alignment information from a read aligner. It represents sequence information in respect to a given reference sequence. The information is stored in a series of tab delimited ascii columns. The full SAM format specification is available at http://samtools.sourceforge.net/SAM1.pdf

On a historical note, the Sanger Institute first used the FASTQ format.

SAM can also (and is increasingly used for it, see PacBio) store unaligned sequence information, and in this regard equivalent to FASTQ.

Note that fasta is also often use for protein data, not just DNA.

FASTA (officially) just stores the name of a sequence and the sequence, unofficially people also add comment fields after the name of the sequence. FASTQ was invented to store both sequence and associated quality values (e.g. from sequencing instruments). SAM was invented to store alignments of (small) sequences (e.g. generated from sequencing) with associated quality values and some further data onto a larger sequences, called reference sequences, the latter being anything from a tiny virus sequence to ultra-large plant sequences.

FASTA and FASTQ formats are both file formats that contain sequencing reads while SAM files are these reads aligned to a reference sequence. In other words, FASTA and FASTQ are the "raw data" of sequencing while SAM is the product of aligning the sequencing reads to a refseq.
A FASTA file contains a read name followed by the sequence. An example of one of these reads for RNASeq might be:
>Flow cell number: lane number: chip coordinates etc.
ATTGGCTAATTGGCTAATTGGCTAATTGGCTAATTGGCTAATTGGCTAATTGGCTAATTGGCTA

The FASTQ version of this read will have two more lines, one + as a place holder and then a line of quality scores for the base calls. The qualities are given as characters with '!' being the lowest and '~' being the highest, in increasing ASCII value. It would look something like this
@Flow cell number: lane number: chip coordinates etc.
ATTGGCTAATTGGCTAATTGGCTAATTGGCTAATTGGCTAATTGGCTAATTGGCTAATTGGCTA
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65

A SAM file has many fields for each alignment, the header begins with the @ character. The alignment contains 11 mandatory fields and various optional ones. You can find the spec file here: https://samtools.github.io/hts-specs/SAMv1.pdf .
Often you'll see BAM files which are just compressed binary versions of SAM files. You can view these alignment files using various tools, such as SAMtools, IGV or USCS Genome browser.
As to the benefits, FASTA/FASTQ vs. SAM/BAM is comparing apples and oranges. I do a lot of RNASeq work so generally we take the FASTQ files and align them the a refseq using an aligner such as STAR which outputs SAM/BAM files. There's a lot you can do with just these alignment files, looking at expression, but usually I'll use a tool such as RSEM to "count" the reads from various genes to create an expression matrix, samples as columns and genes as rows. Whether you get FASTQ or FASTA files just depends on your sequencing platform. I've never heard of anybody really using the quality scores.

Careful, the FASTQ format description is wrong: a FASTQ record can span more than four lines; also, + isn’t a placeholder, it’s a separator between the sequence and the quality score, with an optional repetition of the record ID following it. Finally, the quality score string has to be the same length as the sequence.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Difference between BWA-backtrack and BWA-MEM
Many of my colleagues recommend I use BWA-MEM instead of regular old BWA. The problem is I don't understand why and reading the BWA man page doesn't seem to help the matter.
What is the difference between BWA and BWA-MEM? And, in which instances would you employ one over the other?

You should read the BWA and BWA mem publications.

This paper: https://arxiv.org/abs/1303.3997

Possible duplicate: https://bioinformatics.stackexchange.com/questions/418/why-is-bwa-mem-the-standard-algorithm-when-using-bwa

TL;DR:
BWA-backtrack is based on backtracking.
This approach is appropriate only when the dissimilarity between the reads and the reference is low,
or when you want to find all best hits or enumerate all possible alignments up to a specified
number of errors.
In all other situations, BWA-MEM is preferable as it can,
thanks to its sophisticated strategy based on maximum exact matches,
deal with errors better and also automatically switch between local and
global alignment modes.
Long description:
I would like to provide some algorithmic insight since I believe that it
can be very useful in this case.
Both BWA-backtrack and BWA-MEM use the same indexing strategy
(heavily relying on BWT-index), but the actual algorithms are quite different.
BWT-index (and also other full-text indexes such as
suffix arrays or suffix trees)
can easily find exact matches (imagine Ctrl+F-like search in a text editor),
but any differences between the read and the reference,
such as sequencing errors or genomic variants,
make the situation complicated. One then needs to
somehow transform inexact matching to exact matching,
and all three BWA mappers (note that there exists also BWA-SW, but it is deprecated)
use quite different strategies.
BWA-backtrack looks for substrings of the reference, which would be similar
to the entire read (end-to-end) using an algorithm
called backtracking.
First, it searches occurences of the read without any "corrections".
If nothing found, it consideres all possible single edits; then two edits, etc.
To make mapping efficient, one usually wants to stop with the first
found alignment as this would be the best one.
When required, it is also possible to find the other equally good alignments
or to enumerate all alignments up to some edit distance or withing some divergance rate
(see the -N option of BWA-backtrack).
It turns out that
the time required for finding an alignment
can be exponential in the number of errors, which is probably the main problem of
backtracking-based approaches.
To prevent huge overheads due to dissimilar reads,
one needs to limit the number of allowed errors to some reasonable number
(see maxDiff in the BWA man page)
and consider the other reads unaligned.
In the case of BWA-backtrack,
the minimum required identity level is ~97% with the default options (see the -n option).
In fact, the algorithm is more complicated and uses various heuristics such as
seed-and-extend or Z-dropoff in order to make the computation fast enough (at the price of lower accuracy). If you are interested
in more details, all these tricks are well described in the paper.
BWA-MEM uses quite a different strategy. It detects long exact matches
between the read and the reference, and then chains them into local or global alignments,
based on what is more appropriate in that specific case.
Such an automatic local-global switching can be very powerful and
BWA-MEM works well with various types of data
(short reads, long reads, low error rates, high error rates, etc.).

To quote the Introduction to BWA on sourceforge:

BWA is a software package for mapping low-divergent sequences against a large reference genome, such as the human genome. It consists of three algorithms: BWA-backtrack, BWA-SW and BWA-MEM. The first algorithm is designed for Illumina sequence reads up to 100bp, while the rest two for longer sequences ranged from 70bp to 1Mbp. BWA-MEM and BWA-SW share similar features such as long-read support and split alignment, but BWA-MEM, which is the latest, is generally recommended for high-quality queries as it is faster and more accurate. BWA-MEM also has better performance than BWA-backtrack for 70-100bp Illumina reads.

In short, for anything where you have read lengths over 70bp BWA-MEM is faster, and more accurate.

One important addition is that BWA-MEM will split reads more aggressively. Whereas BWA-ALN requires (if I remember correctly) one end-to-end alignment for anchoring, BWA-MEM allows for partial alignments.
Practicable, use the -M flag for bwa mem if you are processing the alignments with older or less mainstream programs if they don't understand the BAM supplementary read flag yet.

@Manuel: What specific programs do you have in mind that need -M?  That flag was added for use with Picard < 1.96, which is now ancient — Picard itself hasn't benefited from or needed this flag for four years now.  See also https://github.com/lh3/bwa/pull/26


text:
Are there any rolling hash functions that can hash a DNA sequence and its reverse complement to the same value?
A common bioinformatics task is to decompose a DNA sequence into its constituent k-mers and compute a hash value for each k-mer. Rolling hash functions are an appealing solution for this task, since they can be computed very quickly. A rolling hash does not compute each the hash value from scratch with each k-mer: rather it updates a running hash value using an update strategy and a sliding window over the data.
It's also very useful for many applications to have a k-mer hash to the same value as its reverse complement. Unless the data were generated using a strand-specific sample prep, it's impossible to distinguish a k-mer from its reverse complement, and they should be treated as the same sequence.
Are there any rolling hashes that will map reverse complements to the same value? If not, how would we develop such an algorithm?
UPDATE: Ideally the hash function would be able to support k > 32, which would be lossy unless using something larger than a 64-bit integer.
ANOTHER UPDATE: I don't think it's necessary to store both the running k-mer and its reverse complement in a single value. If storing two k-mer strings and/or two hash values makes this easier, I'm totally cool with that.

I don't understand. Surely if you have the kmer you also know the reverse complement of that kmer. It would be a matter of mapping the forward hash to its complement.

See related question on biostars

@morgantaschuk Indeed, it's trivial to compute a k-mer's reverse complement. But you don't want to do this from scratch for each k-mer, since that would negate any performance benefit gained from the rolling hash.

Can you please specify in your question whether the rolling-hash function you envisage is a lossless one (trivial) or a lossy function? E.g, for lossy: mapping a 64 kmer (128 bit info) to a 32 bit hash value.

see also: https://www.biostars.org/p/184993/

A rolling hash function for DNA sequences called ntHash has recently been published in Bioinformatics and the authors dealt with reverse complements:

Using this table, we can easily compute the hash value for the reverse-complement (as well as the canonical form) of a sequence efficiently, without actually reverse- complementing the input sequence, as follows:
  ...

EDIT (by @user172818): I will add more details about how ntHash works. The notations used in its paper are somewhat uncommon. The source code is more informative.
Let's first define rotation functions for 64-bit integers:
rol(x,k) := x << k | x >> (64-k)
ror(x,k) := x >> k | x << (64-k)

We then define a hash function h() for each base. In the implementation, the authors are using:
h(A) = 0x3c8bfbb395c60474
h(C) = 0x3193c18562a02b4c
h(G) = 0x20323ed082572324
h(T) = 0x295549f54be24456
h(N) = 0

The rolling hash function of a forward k-mer s[i,i+k-1] is:
f(s[i,i+k-1]) := rol(h(s[i]),k-1) ^ rol(h(s[i+1]),k-2) ^ ... ^ h(s[i+k-1])

where ^ is the XOR operator. The hash function of its reverse complement is:
r(s[i,i+k-1]) := f(~s[i,i+k-1])
               = rol(h(~s[i+k-1]),k-1) ^ rol(h(~s[i+k-2]),k-2) ^ ... ^ h(~s[i])

where ~ gives the reverse complement of a DNA sequence. Knowing f(s[i,i+k-1]) and r(s[i,i+k-1]), we can compute their values for the next k-mer:
f(s[i+1,i+k]) = rol(f(s[i,i+k-1]),1) ^ rol(h(s[i]),k)  ^ h(s[i+k])
r(s[i+1,i+k]) = ror(r(s[i,i+k-1]),1) ^ ror(h(~s[i]),1) ^ rol(h(~s[i+k]),k-1)

This works because rol, ror and ^ can all be switched in order. Finally, for a k-mer s, the hash function considering both strands is the smaller between f(s) and r(s):
h(s) = min(f(s),r(s))

This is a linear algorithm regardless of the k-mer length. It only uses simple arithmetic operations, so should be fairly fast. I have briefly tested its randomness. It seems comparable to murmur. ntHash is probably the best algorithm so far if you want to hash an arbitrarily long k-mer into 64 bits.

I had a look at the paper and did a small experiment. The idea behind ntHash is interesting, simple, efficient and effective. If @DanielS wants to get a 64-bit hash value for a k-mer of arbitrary length, this should be the accepted answer. This hash function however is not the best for k-mer counting as in theory, it may hash two <=32-mers into the same 64-bit integer by a tiny chance.

@user172818 Thank you for the update. Do you think that it would be possible to easily modify the code to make it work with __int128? It would be just a quick and dirty hack, not a systematic solution, but it could be sufficient for many use cases.

Collisions in hash functions are expected, and should not be discouraged (as long as the function is sufficiently random). The counting of kmers of arbitrary length can be done by storing the actual value (and count) in a bucket that is addressed by the hash.

@Karel I believe many 128-bit operations can't be done in one CPU cycle. Using __int128 will be slower. In addition, nthash incurs collisions anyway. Using longer words just reduce the chance, but won't solve the root problem.

@gringer a popular approach to k-mer counting is sorting. We don't need to hash at all. When using a hash table, a better approach is to use invertible hash functions (e.g. jellyfish & bfc). This way, you don't need to compute hashes as often.

Let f and r be two integers. They always keep the k-mer on the forward and reverse strand, respectively. At a new base c in a proper 2-bit encoding, we update the two integers as follows:
f = (f<<2|c) & ((1ULL<<2*k) - 1)
r = r>>2 | (3ULL-c)<<2*(k-1)

With this updating rule, f keeps the forward strand k-mer ending at c and r is f's reverse complement. The hash of the k-mer can be min(f,r). If integer operations take constant time, computing the hashes of all k-mers of a sequence of length L takes O(L) time.
The above only works if k-mer fits a word, which means k<=32 on x86. For long k-mers, you can overload bit operators in C++. In case of 32<k<=64, a simpler solution is to split into lower and higher bits:
c1 = c&1, c2 = c>>1
f1 = (f1<<1|c1) & ((1ULL<<k) - 1)
f2 = (f2<<1|c2) & ((1ULL<<k) - 1)
r1 = r1>>1 | (1ULL-c1)<<(k-1)
r2 = r2>>1 | (1ULL-c2)<<(k-1)

This is twice as slow as the k<=32 version. It is possible to implement these lines with SSE2 intrinsics, but that is overkilling.
When you want to hash a long k-mer into fewer bits, there are a few options:
xor:      f^r
min:      min(f,r)
min+hash: hash(min(f,r))

where hash() is a generic randomization hash function, which could be Murmur3, Wang's integer hash function, etc.
The following is a microbenchmark. Here, we squeeze all k-mers into a 32-bit hash table. collision_rate = 1 - #distinct_hash/#distinct_kmer. This is not the best metric, but it somehow measures the randomness, which is better than nothing.
====================================================
 Algorithm    data  max-k  k   %collision  CPU time
----------------------------------------------------
 xor          chr11  32    31     8.5%       0.9s
 xor+Wang     chr11  32    31     9.6%       1.3s
 min+Wang     chr11  32    31     1.4%       1.3s
 min+Wang     chr11  64    31     1.4%       1.9s
 min+murmur3  chr11  inf   31     1.4%       5.7s
 min+Wang     chr11  64    51     1.5%       2.1s
 min+murmur3  chr11  inf   51     1.5%       6.8s
 min+Wang     chr11  64    63     1.5%       2.1s
 min+murmur3  chr11  inf   63     1.5%       7.5s
 radix sort   chr11  -     -      -        5.6-7.3s
----------------------------------------------------
 min+Wang     hg38   64    63    26.2%        37s
 min+murmur3  hg38   inf   63    26.2%       192s
 radix sort   hg38   -     -      -         ~138s
====================================================

Observations:

Don't use XOR. In fact, XOR hashes all palindromes to 0xffff.... This is already worrying enough.
Unless you want to work with very long k-mers, don't use generic string hash functions like FNV and Murmur. Murmur is even slower than radix sorting all hashes.
For k<=64, min+Wang is the best here. It is fast and simple to compute and has randomness nearly as good as murmur.

No, your solution does not work if the poster wanted, e.g., compute all hashes for a genome sequence.

Lol, 2-bit encoding is certainly popular but I wouldn't go so far as to call it "proper". But even still, your answer is informative. Add next base to the end of the forward k-mer, add the complement of the next base to the beginning of the reverse k-mer, recompute both hashes.

Damn, wish I could upvote this again. Yet it still has fewer votes than the answer that … recommends Murmur?!

Yeah, wish I could +1 again for the follow through!

If your goal is to minimise storage by just having one hash per kmer and its reverse complement, there is a simple solution for non-rolling hashes. For any sequence S, you compute and store the hash of the smaller of S and its reverse complement. A simple lexicographical comparison for "smaller" is enough. In programming terms:
hash=computeHash(min(S,rev(S));

If your goal is to minimise computing time via rolling hashes AND store only one kmer value, this will be hard. Hashes are normally designed to minimise collisions and you are asking it not only to collide, but also to collide in very, very specific circumstances.
Question is: is the computation of a hash a bottleneck for an application? Maybe, maybe not.
This post on StackExchange seems to imply that, in 2012, MurmurHash was able to compute de-novo a hash for a UUID (similar to a standard kmer in size and complexity) in ~250ns, i.e., ~4 million hashes per second. That's already quite fast.
My guess is that most applications will use a hash to look up things in memory or disk. And here, even with look-ups in memory, real life data will let you run into CPU cache miss speed penalties very quickly and I expect the impact of this to be higher than the speed of the hash computation itself.

Murmur is good but conventional DNA k-mer hashes are much, much faster, since they essentially constitute an identity function (the k-mer is its own hash value). And yes, that performance difference is absolutely noticeable in relevant applications.

At its easiest, you just store the forward ($F$) and reverse ($R$) hash value.
You update the forward hash value by conventional means, e.g. bit-shifting the base value into its lower bits:
$$
F_{n + 1} = ((F_n \ll B \mid x) \mathop\& M,
$$
$B$ is the bit size of the base encoding, $M$ is the word mask for a word of length $W$ bits, and can be omitted if the width of the hash is the width of the machine word; and the reverse hash value by doing the mathematical inverse, e.g. bit-shifting the value of the  reverse complement of the base into the upper bits:
$$
R_{n + 1} = ((R_n \gg B) \mid (\operatorname{compl}x \ll (W - B))) \mathop\& M.
$$
And then you compute a combined hash value by xoring the two hashes, $H  = F \oplus R$.

I was hoping for something that would work more generically (2-bit encoding is limited to k <= 32), but this is very insightful. And for k <= 32, it will likely be much more performant than an adapted generic hash function.

@DanielStandage, I have answered what to do for k>32: you just use more integers. That will still be fast, times faster than any alternatives. Also, this answer is not perfectly correct. In 2-bit encoding, you need to shift by 2.

Furthermore, XOR actually doesn't work well. Suppose we only have 0 base and 1 base. 0 and 1 are complement to each other. You can enumerate all 3-mers and their complements, you will find that there are only two possible hash values: 111 or 010. XOR is more likely to cause hash collisions.

@DanielStandage You say “2-bit encoding is limited to k <= 32” but that’s not actually true, you just need bigger machine words (std::bitset comes to mind, and SeqAn has its own type, AFAIR). But at any rate my scheme isn’t limited to 2-bit encoding; all it requires is that the rolling hash can be computed by shifting in a value at either end.

In our testing, using std::bitset (or rolling our own large int using int arrays) led to a >10x hit to performance. I'd prefer to stick with native int types.

@DanielStandage Well of course there’s a performance hit. That said, you can of course truncate the k-mer bit representation. This is then of course equivalent to a shorter k-mer (i.e. collisions on the full k-mer) but the pigeonhole lemma trivially proves that it’s not possible to preserve all information without collisions, whatever you do …

Don't use bitset. It has a loop inside! The right way is to set a max k-mer length and manually write bit operations. Supporting up to 64bp, which is often sufficient, should only be twice as slow. SSE2 can shift two 64-bit integers in parallel and is likely to give you a faster implementation.

@user172818 And have you checked that the loop isn’t unrolled? It has a compile-time constant width. I haven’t checked this particular instance but such loops in general can be unrolled. And modern compilers usually also emit the correct SSE intrinsics if you target a supported architecture.

A good catch. No, I haven't. I was deducing from the 10x performance hit. 64bp should only be twice as slow than 32bp. The 10x hit wouldn't make sense if the compiler were generating the optimal assembly.


text:
Feature annotation: RefSeq vs Ensembl vs Gencode, what's the difference?
What are the actual differences between different annotation databases?
My lab, for reasons still unknown to me, prefers Ensembl annotations (we're working with transcript/exon expression estimation), while some software ship with RefSeq annotations. Are there significant differences between them today, or are they, for all intents and purposes, interchangeable (e.g., are exon coordinates between RefSeq and Ensembl annotations interchangeable)?

To add to rightskewed answer:
While it is true that:
Gencode is an additive set of annotation (the manual one done by Havana and an automated one done by Ensembl),
the annotation (GTF) files are quite similar for a few exceptions involving the X chromosome and Y par and additional remarks in the Gencode file (see more at FAQ - Gencode).

What are the actual differences between different annotation databases?

They are a few differences, but the main one  for me (and it could be stupid) is
that Refseq is developed by the American NCBI and
the ENSEMBL is mainly developed by the European EMBL-EBI.
Often, labs or people will just start using what is the best known to them (because of a course or workshop) or because they start working with one of the databases with one specific tool and keep with it later.

My lab, for reasons still unknown to me, prefers Ensembl annotations (we're working with transcript/exon expression estimation), while some software ship with RefSeq annotations.

Your lab might be mostly European based people or they might also have read papers like the one from Frankish et al.  Comparison of GENCODE and RefSeq gene annotation and the impact of reference geneset on variant effect prediction. BMC Genomics 2015; 16(Suppl 8):S2 - DOI: 10.1186/1471-2164-16-S8-S2
From the Frankish et al. paper paper:

The GENCODE Comprehensive transcripts contain more exons, have greater genomic coverage and capture many more variants than RefSeq in both genome and exome datasets, while the GENCODE Basic set shows a higher degree of concordance with RefSeq and has fewer unique features.

As for:

Are there significant differences between them today, or are they, for all intents and purposes, interchangeable (e.g., are exon coordinates between RefSeq and Ensembl annotations interchangeable)?

No. I don't think they are great differences between them as that the global picture should stay the same (although you will see different results if you are interested in a small set of genes). However, they are not directly interchangeable. Particularly as there are many versions of Ensembl and Refseq based on different genome annotations (and those won't be interchangeable between themselves either in most cases).
However, you can easily translate most[1] of your Refseq IDs to ENSEMBL IDs and vice-versa with tools as http://www.ensembl.org/biomart/martview for example (there are devoted libraries/API as well like Biocondutor: biomaRt
[1] Most as sometimes, they might be annotated in one of the database but haven't (yet) an equivalent in the other.
EDIT
In fine, even if people tends to keep to what they are used to (and that the annotations are constantly expanded and corrected) depending on the research subject one might be interested in using one database over another:
From Zhao S, Zhang B. A comprehensive evaluation of ensembl, RefSeq, and UCSC annotations in the context of RNA-seq read mapping and gene quantification. BMC Genomics. 2015;16: 97. paper:

When choosing an annotation database, researchers should keep in mind that no database is perfect and some gene annotations might be inaccurate or entirely wrong. [..] Wu et al. [27] suggested that when conducting research that emphasizes reproducible and robust gene expression estimates, a less complex genome annotation, such as RefGene, might be preferred. When conducting more exploratory research, a more complex genome annotation, such as Ensembl, should be chosen.
[..]
[27] Wu P-Y, Phan JH, Wang MD. Assessing the impact of human genome annotation choice on RNA-seq expression estimates. BMC Bioinformatics. 2013;14(Suppl 11):S8. doi: 10.1186/1471-2105-14-S11-S8.

Ensembl vs Gencode
https://www.gencodegenes.org/faq.html

The GENCODE annotation is made by merging the Havana manual gene annotation and the Ensembl automated gene annotation. [...] In practical terms, the GENCODE annotation is identical to the Ensembl annotation.

Further, for the GTF file differences:

The only exception is that the genes which are common to the human chromosome X and Y PAR regions can be found twice in the GENCODE GTF, while they are shown only for chromosome X in the Ensembl file.

Gencode(Ensembl) vs RefSeq
Gencode is in almost all cases more comprehensive. For example, this is NCBI RefSeq vs Ensembl (v24, release 83) for BRCA gene:

RefSeq and Gencode are not interchangeable in most cases, though RefSeq annotations will often be a subset of the Gencode ones.

How is the BRCA screenshot making your point? It's not very obvious at first sight (at least for me)

To add practical advice to what others have said:
In a practical sense, I think the biggest difference between RefSeq and Ensembl/GENCODE is in the sensitivity/specificity trade off.
Ensembl aims more towards the inclusive end, including a far larger number of transcript variants, many of which are only weakly supported.
RefSeq trades some of this sensitivity for specificity - you can be more confident that a RefSeq transcript exists, but less confident that the ReqSeq annotation includes all of the real transcripts for a gene.

That's why I prefer the Ensembl annotation as you can query for a most confident set by selecting only the Havana (Havana or Ensembl/Havana) transcripts. See: http://www.ensembl.org/Help/Faq?id=152

While annotations between RefSeq and Gencode are not so different is the coding regions (genes), Gencode is much more rich in the intergenic regions. This could be very advantagous for epigenetic studies, where regulation is of interest.

Hi! Is there any publication or other material you could link regarding your claim? That would be very interesting...

The UCSC Genome Browser Genes FAQ discusses this question in detail:
https://genome.ucsc.edu/FAQ/FAQgenes.html#ens

Officially, the Ensembl and GENCODE gene models are the same. On the
  latest human and mouse genome assemblies (hg38 and mm10), the
  identifiers, transcript sequences, and exon coordinates are almost
  identical between equivalent Ensembl and GENCODE versions (excluding
  alternative sequences or fix sequences).
GENCODE uses the UCSC convention of prefixing chromosome names with
  "chr", e.g. "chr1" and "chrM", but Ensembl calls these "1" or "MT". At
  the time of writing (Ensembl 89), a few transcripts differ due to
  conversion issues. In addition, around 160 PAR genes are duplicated in
  GENCODE but only once in Ensembl. The differences affect fewer than 1%
  of the transcripts. Apart from gene annotation itself, the links to
  external databases differ.
The GENCODE Release History shows the release dates and can be linked
  to corresponding Ensembl releases. You can download the gene
  transcript models from the website https://gencodegenes.org or from
  http://ensembl.org. For most applications, the files distributed on
  the GENCODE website should be easier to use, as the third party
  database links are easier to parse and the sequence identifiers match
  the UCSC genome files, at least for the primary chromosomes.
Additional information on this question can be found on the GENCODE
  FAQ page.

Note that the pDot of a certain mutation might differ based on which transcript you are using, and so Ensembl and RefSeq cannot be used interchangeably.

This is a result of transcript coordinates not aligning 100% (already addressed elsewhere) and does not merit being its own answer.

It is a real concern where I work so I thought it would helpful to make it clear to users, but I suppose it is a very niche concern (for example, with the MET gene)

Yes, the protein coordinates changing can be a problem but the solution is to use MANE transcripts so the annotation is database agnostic.


text:
Convert a BAM file from one reference to another?
I have a set of BAM files that are aligned using the NCBI GRCh37 human genome reference (with the chromosome names as NC_000001.10) but I want to analyze it using a BED file that has the UCSC hg19 chromosome names (e.g. chr1). I want to use bedtools to pull out all the on-target and off-target reads.

Are NCBI and UCSC directly comparable? Or do I need to re-align the BAM/lift-over the BED to the UCSC reference?
Should I convert the BED file or the BAM file? Everyone here uses the UCSC chromosome names/positions so I'll need to convert the eventual files to UCSC anyway.

You're the second person I have ever seen using NCBI "chromosome names" (they're more like supercontig IDs). Normally I would point you to a resource providing mappings between chromosome names, but since no one has added NCBI names (yet, maybe I'll add them now) you're currently out of luck there.
Anyway, the quickest way to do what you want is to samtools view -H foo.bam > header to get the BAM header and then change each NCBI "chromosome name" to its corresponding UCSC chromosome name. DO NOT REORDER THE LINES! You can then use samtools reheader and be done.
Why, you might ask, would this work? The answer is that chromosome/contig names in BAM files aren't stored in each alignment. Rather, the names are stored in a list in the header and each alignment just contains the integer index into that list (read group IDs are similar, for what it's worth). This also leads to the warning above against reordering entries, since that's a VERY convenient way to start swapping alignments between chromosomes.
As an aside, you'd be well served switching to Gencode or Ensembl chromosome names, they're rather more coherent than the something_random mess that's present in hg19 from UCSC.
Update: Because I'm nice, here is the conversion between NCBI and UCSC. Note that if you have any alignments to patches that there is simply no UCSC equivalent. One of the many reasons not to use UCSC (avoid their annotations too).

How well does this work? Have you run any benchmarks? I ask because I tried to convert various bed files between hg and GRC genomes and the three tools I used all gave very different results. This sort of mapping really should be a simple thing but it appears not to be that straightforward at all.

For cases where it's just a change of name (most cases) there's nothing to benchmark. For cases where you additionally have position changes then you'd need a different resource (namely, liftOver or crossmap).

Yeah, it's tolls like liftOver and crossmap I've used and found issues with. I was expecting this to be a solved problem but each of the three tools I used gave different results, unfortunately. Which makes me leery of using the results.

The results are deterministic, they should be the same regardless of tool provided you use the same settings.

You'd think so, yes. But they're not. I'd be happy to chat ([chat]) about this (this is one of things I need to sort out at work) if you like, and I'd be very pleased to learn I was just doing it wrong, but my preliminary tests and various threads I've read online (see here for example) suggest it isn't as simple as you'd think.

I suspect this will come down to deficiencies in the chain files, but I'd like to hear about that.

The "right" solution would be realignment, but that's expensive and most of us would not go that route. My preferred solution would be to convert the bed file, as opposed to the bam. Here's why:
1) Reheadering the bam means that you may have reads aligned to contigs without a corresponding entry in UCSC (see Devon's list for the mappings). This is a problem because:

Some of those reads would likely have been mapped elsewhere if a reference without those contigs was used.
I'm not even sure what happens to those reads after reheadering - I guess they would need to be marked as unmapped? Lots of potential for screwiness there.

2) It seems cleaner to convert the bed file from UCSC->NCBI, where you are guaranteed that every entry has a "home". Then, after you pull your info from the bam, you can always convert chromosome names back if you need to.

Alignments that end up losing their chromosome end up unmapped, though they don't get the SAM unmapped flag (rather, they end up represented as being on the * chromosome). But yeah, samtools reheader should be used with extreme caution.

Several tools will throw a wobbly if you pass reads on * that are not marked as unmapped. pysam is is one. I know this from bitter (bug fixing ) experience.

Just to point out that if you want to follow @Devon Ryan's answer for a different organism/assembly, that is not in his very useful linked resource, you can download NCBI to UCSC contig to chromosome number mappings from https://www.ncbi.nlm.nih.gov/assembly.
To the site and search for your assembly name. At the bottom of the page is a box called "Global assembly definition" containing a link titled "Download full sequence report".
The downloaded file contains a table with:

Chromosome numbers in "Sequence-Name"/"Assigned-molecule"
NCBI names in Refseq-Accn
UCSC contig name in "UCSC-style-name"


text:
Pseudo-temporal ordering in heterogeneous populations
What's the most widely accepted tool for doing pseudo-temporal ordering from scRNAseq data? Also is there away to separate differential expression that occurs based on "cell identity" or maybe more accurately cell type fate from that which arises from cells being in different stages differentiation.
To be more concrete, lets say there's a population of cells, some of which were born at time 1, time 2, and time 3. The progression along the temporal trajectory can be described via a set of genes that are fluctuating as the cell matures. So you might have the same cell type which was born at time 1 be transcriptionally distinct from a younger one born at time 3.  On the other hand, within this population there are subpopulations which will have different cell fates and are transcriptionally distinct. Is there away to reliably separate the temporal axis from the "cell fate axes". If not is this something people are working on or is flawed logic to think this kind of thing is possible?

I think this question might be too speculative for this platform. In particular, this involves defining cell types and states, which in many ways is philosophical.

What does it mean the tag? it stands for small conditional RNA? I think then that this question has two different issues, one about scRNA and another about separating subpopulations from a pool. Could you clarify it?

I think the tag stands for single-cell RNA-seq (correct me if wrong?). Although the question is broad, I do know people working on this problem and developing tools - unfortunately I'm not an expert.

@SarahCarl Maybe you can become one researching an answer for this question :D. If it is too broad, I think you can vote it too, I think it doesn't show research effort, but I am unsure if it is too broad or not

Without being an expert on scRNA-seq it sounds to me like you are looking for something like the monocle package.

@Llopis I think it is too broad - but if edited could work well. If this question were instead a question about pros/cons of a few different tools, that might work better.

For reference, here is a list with most (but not all) methods available: https://github.com/agitter/single-cell-pseudotime . I don't think this is a question which can have a clear answer without major research. See also this review of strategies: http://onlinelibrary.wiley.com/wol1/doi/10.1002/eji.201646347/abstract . Monocle is probably the most well known, but has been criticized often, leading to people making new tools when it doesn't work.

@Valentine thanks for your response, I've seen the github page before but it's good to link it here. And thanks for the review I'll read through it and see what I can glean from it. I wanted some opinions on what method/tool seems like it does the best job at this at this point in time. If anybody has opinions or experiences with this I'd love to hear them.

Also we'd like to be able to incorporate known expression patterns into the models.

@Alon Your question is interesting, but it's sounding like a research project rather than something with a readily available tool :) . If you know your expression patterns are switch-like, have a look at Ouija: https://github.com/kieranrcampbell/ouija. If the patterns are different, try finding a statistician nearby who is familiar with latent variable models.

Regarding the first part of the question: scRNA-seq is a rapidly developing field so it may be hard to talk about "widely accepted tool for doing pseudo-temporal ordering from scRNAseq data". Few of the tools aiming to do just that include Monocole, Waterfall or Sincell (see this paper for references https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4122333/)
The second part is a bit more complex. Many aspects have come together here, e.g. cell-cycles, sub-populaitons identification and pseudo-temporal ordering, to get a true reflection of the biological processes. There are efforts on all (and more) these fronts (see again e.g. the above paper, though not the latest) and there are most likely people working on their integration. I'm not aware of any studies published yet at this depth though


text:
Compare alignment quality of multiple sequencing runs aligned against the same reference genome
I have run Oxford Nanopore Technologies' MinION sequencing on the same DNA sample using three flowcells, each aligned against the same reference genome (E.coli K12 MG1655) using both BWA MEM and GraphMap and stored as BAM files.
How can I quantitatively and efficiently analyse the quality of alignment (percentage identity, insertion rate, deletion rate) of each of these files?

Qualimap will do this for you.

Go to qualimap.bioinfo.cipf.es
Run qualimap (default params are fine) on each BAM file
Open up the HTML output, and you can read off the %identity (they measure the opposite, i.e. mismatch rate, but 100% - mismatch rate is %identity of course), indel rate, etc.

One thing to watch out for (you don't mention it in your question, but just in case) is that you cannot directly compare Q scores - these are a bit of a mess and calculated very differently in each piece of software.
Unsolicited suggestion: you might also try NGM-LR for mapping MinION data. We've found it beats the others for our data (though we map to a distant reference).

Thanks for the suggestion. To clarify, is this the aligner you were referring to? https://github.com/philres/ngmlr

Yep. Poster is here: http://schatzlab.cshl.edu/publications/posters/2015/2015.GI.NextGenMap-LR.pdf
You could check out Sniffles too - detects structural variants from 3rd gen data. https://github.com/fritzsedlazeck/Sniffles


text:
Read length distribution from FASTA file
I have a single ~10GB FASTA file generated from an Oxford Nanopore Technologies' MinION run, with >1M reads of mean length ~8Kb. How can I quickly and efficiently calculate the distribution of read lengths?
A naive approach would be to read the FASTA file in Biopython, check the length of each sequence, store the lengths in a numpy array and plot the results using matplotlib, but this seems like reinventing the wheel.
Many solutions that work for short reads are inadequate for long reads. If I'm hey output a single (text) line per 1/10 bases, which would lead to a text output of upwards of 10,000 lines (and potentially more than 10x that) for a long read fasta.

If you want something quick and dirty you could rapidly index the FASTA with samtools faidx and then put the lengths column through R (other languages are available) on the command line.
samtools faidx $fasta
cut -f2 $fasta.fai | Rscript -e 'data <- as.numeric (readLines ("stdin")); summary(data); hist(data)'

This outputs a statistical summary, and creates a PDF in the current directory called Rplots.pdf, containing a histogram.

You can go even simpler with something like: faidx --transform chromsizes | cut -f2 | Rscript -e 'data <- as.numeric (readLines ("stdin")); summary(data); hist(data)'. This requires pyfaidx: pip install pyfaidx.

How can the Rscript be modified so the pdf file has a name based on the sample name?

Statistics for nanopore reads are tricky because of the huge range of read lengths that can be present in a single run. I have found that the best way to display lengths is by using a log scale on both the x axis (length) and the y axis (sequenced bases, or counts, depending on preference).
I have written my own scripts for doing this: one for generating the read lengths, and another for plotting the length distribution in various ways. The script that generates read lengths also spits out basic length summary statistics to standard error:
$ ~/scripts/fastx-length.pl > lengths_mtDNA_called.txt
Total sequences: 2110
Total length: 5.106649 Mb
Longest sequence: 107.414 kb
Shortest sequence: 219 b
Mean Length: 2.42 kb
Median Length: 1.504 kb
N50: 336 sequences; L50: 3.644 kb
N90: 1359 sequences; L90: 1.103 kb

$ ~/scripts/length_plot.r lengths_mtDNA_called.txt
lengths_mtDNA_called.txt ... done
Number of sequences: 2110
Length quantiles:
      0%      10%      20%      30%      40%      50%      60%      70%
   219.0    506.9    724.4    953.0   1196.2   1503.0   1859.2   2347.3
     80%      90%     100%
  3128.2   4804.7 107414.0

Here are a couple of of the produced graphs:

The scripts to generate these can be found here:

https://github.com/gringer/bioinfscripts/blob/master/fastx-length.pl
https://github.com/gringer/bioinfscripts/blob/master/length_plot.r

Using Biopython and matplotlib would seem like the way to go, indeed.
It really just boils down to three lines of code to get that graph:
import Bio, pandas
lengths = map(len, Bio.SeqIO.parse('/path/to/the/seqs.fasta', 'fasta'))
pandas.Series(lengths).hist(color='gray', bins=1000)

Of course you might want to make a longer script that's callable from the command line, with a couple options. You are welcome to use mine:
#!/usr/bin/env python2

"""
A custom made script to plot the distribution of lengths
in a fasta file.

Written by Lucas Sinclair.
Kopimi.

You can use this script from the shell like this:
$ ./fastq_length_hist --input seqs.fasta --out seqs.pdf
"""

###############################################################################
# Modules #
import argparse, sys, time, getpass, locale
from argparse import RawTextHelpFormatter
from Bio import SeqIO
import pandas

# Matplotlib #
import matplotlib
matplotlib.use('Agg', warn=False)
from matplotlib import pyplot

################################################################################
desc = "fasta_length_hist v1.0"
parser = argparse.ArgumentParser(description=desc, formatter_class=RawTextHelpFormatter)

# All the required arguments #
parser.add_argument("--input", help="The fasta file to process", type=str)
parser.add_argument("--out", type=str)

# All the optional arguments #
parser.add_argument("--x_log", default=True, type=bool)
parser.add_argument("--y_log", default=True, type=bool)

# Parse it #
args        = parser.parse_args()
input_path  = args.input
output_path = args.out
x_log       = bool(args.x_log)
y_log       = bool(args.y_log)

################################################################################
# Read #
lengths = map(len, SeqIO.parse(input_path, 'fasta'))

# Report #
sys.stderr.write("Read all lengths (%i sequences)\n" % len(lengths))
sys.stderr.write("Longest sequence: %i bp\n" % max(lengths))
sys.stderr.write("Shortest sequence: %i bp\n" % min(lengths))
sys.stderr.write("Making graph...\n")

# Data #
values = pandas.Series(lengths)

# Plot #
fig   = pyplot.figure()
axes  = values.hist(color='gray', bins=1000)
fig   = pyplot.gcf()
title = 'Distribution of sequence lengths'
axes.set_title(title)
axes.set_xlabel('Number of nucleotides in sequence')
axes.set_ylabel('Number of sequences with this length')
axes.xaxis.grid(False)

# Log #
if x_log: axes.set_yscale('symlog')
if y_log: axes.set_xscale('symlog')

# Adjust #
width=18.0; height=10.0; bottom=0.1; top=0.93; left=0.07; right=0.98
fig.set_figwidth(width)
fig.set_figheight(height)
fig.subplots_adjust(hspace=0.0, bottom=bottom, top=top, left=left, right=right)

# Data and source #
fig.text(0.99, 0.98, time.asctime(), horizontalalignment='right')
fig.text(0.01, 0.98, 'user: ' + getpass.getuser(), horizontalalignment='left')

# Nice digit grouping #
sep = ('x','y')
if 'x' in sep:
    locale.setlocale(locale.LC_ALL, '')
    seperate = lambda x,pos: locale.format("%d", x, grouping=True)
    axes.xaxis.set_major_formatter(matplotlib.ticker.FuncFormatter(seperate))
if 'y' in sep:
    locale.setlocale(locale.LC_ALL, '')
    seperate = lambda x,pos: locale.format("%d", x, grouping=True)
    axes.yaxis.set_major_formatter(matplotlib.ticker.FuncFormatter(seperate))

# Save it #
fig.savefig(output_path, format='pdf')

EDIT - an example output:

There are several potential approaches. For example:

histogram of sequence lengths in the Biopython tutorial
plot_distribution from the Ruby-based biopieces framework
various solutions to get sequence length including bioawk and EMBOSS infoseq

As to which of these are "quick and efficient" using a 10 GB file...it's hard to say in advance. You may have to try and benchmark a few of them.

You specifically asked about FASTA files, but it's important to always consider read length and quality jointly when assessing high-error long-read data. FASTA files do not provide the quality. This information will help you determine how successful the run was, how many reads were 'high quality', etc.
I originally posted a full answer here, suggesting pauvre, but I decided it was a little off topic since you seem to only have the FASTA files.
I recommend generating FASTQ files, but I'm not sure whether you have the original base-called fast5 files. If so, generate FASTQs using poretools as follows (poretools doc for generating FASTQ files):
poretools fastq fast5/

Then I recommend generating a heat map and histogram margin plot using pauvre with both read length and quality like shown below.

[Note: I am not the original author for pauvre, but I am now contributing to this project]

I would appreciate any feedback on why my answer may be disagreeable.

Dear Mark, I am having problems with pauvre. Can you help me?

bioawk could be reasonably efficient for this kind of task.
$ bioawk -c fastx '{histo[length($seq)]++} END {for (l in histo) print l,histo[l]}' \
    | sort -n
0   33270
1   1542
2   1132
3   3397
4   8776
5   11884
6   12474
7   14341
8   13165
9   15467
10  21089
11  30469
12  45204
13  62311
14  88744
15  115767
16  140770
17  191810
18  313088
19  518111
20  1097867
21  4729715
22  6575557
23  2734062
24  1015476
25  493323
26  323827
27  164419
28  107120
29  72487
30  40120
31  24538
32  22295
33  13121
34  9382
35  4847
36  3858
37  3161
38  2852
39  2388
40  1639
41  961
42  686
43  377
44  199
45  114
46  78
47  59
48  50
49  52
50  48
51  42
52  39
53  28
54  49
55  59
56  55
57  51
58  55
59  43
60  52
61  56
62  48
63  67
64  95
65  488

The -c fastx tells the program to parse the data as fastq or fasta. This gives access to the different parts of the records as $name, $seq (and $qual in the case of fastq format) in the awk code (bioawk is based on awk, so you can use whatever language features you want from awk).
Between the single quotes come a series of <condition> {<action>} blocks.
The first one has no <condition> part, which mean it is executed for each record. Here, it updates the lengths counts in a table which I named "histo". length is a predefined function in awk.
In the second block, the END condition means we want it to be executed after all the input has been processed. The action part consists in looping over the recorded length values and print them together with the associated count.
The output is piped to sort -n in order to sort the results numerically.
On my workstation, the above code took 20 seconds to execute for a 1.2G fasta file.

I realize that the output will not be convenient when dealing with sparse length values, because there is no binning (or, equivalently, the bins have width 1).

It is not exactly what you asked, but you can generate a histogram of read length distribution of your nanopore data directly from the HDF5 files using poretools

Due to disk space constraints, default behaviour for nanopore basecallers no longer creates the FAST5 (HDF5) files as output. While I actually have them, most users will not, and additionally this method would not generalize to other sequencers like PacBio.


text:
What is the difference between a transcriptome and a genome?
I have a computer engineering background, not biology.
I started working on a bioinformatics project recently, which involves de-novo assembly. I came to know the terms Transcriptome and Genome, but I cannot identify the difference between these two.
I know a transcriptome is the set of all messenger RNA molecules in a cell, but am not sure how this is different from a genome.

What information have you researched ? Did Wikipedia pages of Transcriptome and Genome helped you? If not where they didn't help you

This is a rather basic question and I don't think it's a good one for the site, in the message from stack exchange it says "if you ask beginner questions, survey questions, or social-conversation questions, experts and pros will not be interested"

This is not a bioinformatics question, this is a terminology question and needs to be be better framed. The approaches to transcriptome and genome de novo assembly would be accceptable, but this is a question for Google.

You can vote to close the question

In brief, the  “genome”  is the collection of all  DNA  present  in  the  nucleus  and  the  mitochondria of a  somatic  cell. The initial product of genome expression is the “transcriptome”, a collection of RNA molecules derived from those genes.

If you want to learn more about relationships between these terms, see this Wikipedia page.

This is a good concise answer.

They are two very different things. Your genome is a large section of about 3 billion DNA nucleotide bases. It has no concept of exon and introns.
Transcriptome is a study of transcriptions. You have introns and exons. We can now talk about alternative splicing and gene expression.
You can think your genome is like a cooking recipe. While it's good to have a good recipe, you can't do much if you don't use it for cooking.

I like to compare the genome with the cookbook, the transcriptome with your shopping lists, the proteome with your groceries.

wouldn't the proteome be the cooked food?


text:
Alignment based vs reference-free (transcriptome analysis)?
I want to focus on transcriptome analysis. We know it's possible to analyze RNA-Seq experiment based on alignment or k-mers.
Possible alignment workflow:

Align sequence reads with TopHat2
Quantify the gene expression with Cufflinks

Possible reference-free workflow:

Quantify sequence reads with Kallisto reference-free index

Both strategy generate gene expression table.
Q: What are pros and cons for each of the approach? Can you give guideline?

Reference-free quantification doesn’t really exist (and Kallisto isn’t it).

Here is a project that tries to build tools for high throughput sequencing data without a reference: https://colibread.inria.fr/project/

First of all, I would emphasize that "alignment-free" quantification tools like Salmon and Kallisto are not reference-free. The basic difference between them and more traditional aligners is that they do not report a specific position (either in a genome or transcriptome) to which a read maps. However, their overall purpose is still to quantify the expression levels (or differences) of a known set of transcripts; hence, they require a reference (which could be arbitrarily defined).
The most important criterion for deciding which approach to use (and this is true of almost everything in genomics) is exactly what question you would like to answer. If you are primarily interested in quantifying and comparing expression of mature mRNA from known transcripts, then a transcriptome-based alignment may be fastest and best. However, you may miss potentially interesting features outside of those known transcripts, such as new isoforms, non-coding RNAs, or information about pre-mRNA levels, which can often be gleaned from intronic reads (see the EISA method).
This paper also has some good considerations about which tools may work best depending on the question you want to answer.
Finally, another fast and flexible aligner (which can be used with or without a reference transcriptome) is STAR.

I wouldn't say Kallisto (or Salmon) are reference-free. They use a transcriptome as reference anda concept called pseudo-alignment which greatly speed up the process of assigning your reads to a transcript.
That said, both approaches of (i) mapping against a reference genome (what you called alignment workflow ) and (ii) mapping against a reference transcriptome will serve different purposes
Transcriptome mapping using pseudoalignent is becoming the method of choice for gene/transcript quantification and differential expression analysis. The drawback is that you only focus on known transcripts
Two typical workflows are:

Kallisto followed by sleuth
Salmon, followed by tximport and DESeq2/EdgeR

Genome mapping is useful for, per example discovery of new isoforms. You shouldn't use TopHat anymore as it has been discontinued by the author.
A typical workflow would be:

Hisat2 (alignment)
StringTie (transcript assembly and abundance estimation)
Ballgown (differential expression)

It's just worth noting (in addition to this very useful answer) that Salmon can also be used with sleuth by means of the wasabi package.


text:
Normalization methods with RNA-Seq ERCC spike in?
ERCC spike-in is a set of synthetic controls developed for RNA-Seq. I'm interested in using it to normalize my RNA-Seq samples. In particular, I'd like to use the spike-ins to remove technical bias and any variation that should not be part of my analysis.
The site doesn't give any details on how I can do that.
Q: What are the possible normalization strategies? Can you briefly describe them?

Are you interested in bulk or single-cell RNA-seq? The value of spike-ins is vastly different depending on that

You may consider using RUVSeq. Here is an excerpt from the 2013 Nature Biotechnology publication:

We evaluate the performance of the External RNA Control Consortium (ERCC) spike-in controls and investigate the possibility of using them directly for normalization. We show that the spike-ins are not reliable enough to be used in standard global-scaling or regression-based normalization procedures. We propose a normalization strategy, called remove unwanted variation (RUV), that adjusts for nuisance technical effects by performing factor analysis on suitable sets of control genes (e.g., ERCC spike-ins) or samples (e.g., replicate libraries).

RUVSeq essentially fits a generalized linear model (GLM) to the expression data, where your expression matrix $Y$ is a $m$ by $n$ matrix, where $m$ is the number of samples and $n$ the number of genes. The model boils down to
$Y = X*\beta + Z*\gamma + W*\alpha + \epsilon$
where $X$ describes the conditions of interest (e.g., treatment vs. control), $Z$ describes observed covariates (e.g., gender) and $W$ describes unobserved covariates (e.g., batch, temperature, lab). $\beta$, $\gamma$ and $\alpha$ are parameter matrices which record the contribution of $X$, $Z$ and $W$, and $\epsilon$ is random noise. For subset of carefully selected genes (e.g., ERCC spike-ins, housekeeping genes, or technical replicates) we can assume that $X$ and $Z$ are zero, and find $W$ - the "unwanted variation" in your sample.

We have added ERCC spike-ins to all our RNASeq data, just in case other people might find it useful in the future. However, I have never used it in my own analyses because I can't think of a reasonable way that it could be used.
The typical recommendation for ERCC is to add it in proportion to the input RNA amount, but that makes an assumption that total cell RNA counts are similar across different cells (which is demonstrably false by looking at single cell RNASeq results).
I have yet to think of a situation where ERCC would provide better results than a "housekeeping" gene set sampled from the original reads.

Why would you spike-in ERCC if you don't have any use?

We do the same thing, the sequencing depth needed is really small, so it's cheap and "better safe than sorry".

We tried to find things to add in that would mean we wouldn't need to re-do runs in the future.

One potential benefit is that if you seem extreme biases in certain samples (e.g. a HUGE % of reads going to ERCC) then you know that something went wrong with the nucleic acids (e.g. bad extraction, too low input, etc...).

Yes, I guess ERCCs are a reasonable positive control for the sample prep. Those sample prep issues tend to come out through other means as well (e.g. high ribosomal mapping proportion, low mapping rate, GC difference, number of transcripts expressed above X level, PCA).


text:
Linear models of complex diseases
A popular framework to analyze differences between groups, either experiments or diseases, in transcriptomics is using linear models (limma is a popular choice).
For instance we have a disease D with three stages as defined by clinicians, A, B and C. 10 samples each stage and the healthy H to compare with is RNA-sequenced. A typical linear model would be to observe the three stages~A+B+C independently. The data of each stage is not from the same person. (but for the question assume it isn't)
My understanding is that such a model would not take into account that stage C appears only on 30% of patients in stage B. And that a healthy patient upon external factors can jump to stage B.
If we want to find the role of a gene in the disease we should include somehow this information in the model. Which makes me think about mixing linear models and hidden Markov chains.
How can such a disease be described in terms of linear models with such data and information?

Of course, your model is not taking into account that observations (i.e 30% of patients from B evolve to C). Here, limma would look for genes with most significant differences among states based on F-statistic. If you are more interested in disease progression, you should read more about time-course patterns (states follow a progression, so positions are not transposable).

But time-course patterns wouldn't handle the jump from healthy to B. A quick look for disease time course patterns show the top hit with an article from 2001. Could you point me to a review or an introductory title

Could you clarify (probably in the question) if stages A, B & C come from the same individuals?

@IakovDavydov Done, the stages are from the disease, the data might come from the same individuals or not. But we can assume the data from each stage doesn't come from the same individuals

Why do you think it's useful to include the progression information in the model ("we should include somehow this information in the model"). You're only able to discern when/if the gene is changing and no statistical model will change that.

What about saying that if you progressed to B from C, then you are positive for both those variables, rather than that you just have B?

@DevonRyan My rationale behind is the following: imagine that a gene between B and C is found to change the expression. Is the change in the gene expression a consequence or a reason of changing the stage?

@IanSudbery I thought about that possibility but clinicians don't think that way. But I would actually like to test this model with the dataset I have at hands.

@Llopis You're not going to get causality from this sort of data.

I don't understand your description of the linear model as ~A+B+C. If these are indeed stages of a disease - then the model would rightly just be ~Disease. And the levels of the data column or vector would be the factors (H,A,B,C) - where H is reset to the reference level.

@StephenHenderson Doing that way the model I couldn't compare between stages.

@Llopis I believe you are confused and have that the wrong way round. I would refer you to section 9.3 page 42 of the limma manual on dealing with several groups. Note the actual congrats of interest are specified separately.. https://www.bioconductor.org/packages/devel/bioc/vignettes/limma/inst/doc/usersguide.pdf

@StephenHenderson AFAIK, section 9.3 use ~0+f (where f is the factor with the groups) which is the same as ~A+B+C (that's why later the columns are renamed to RNA1, RNA2, RNA3 in the tutorial, if It had one Disease/Group column how could it rename three columns?)

@Llopis that's the design matrix - 3 columns for 3 groups in the single covariate (f). ~A+B+C is like a model with 3 covariates Disease+Sex+Age... see for instance section 9.42 on the next page 43.

@Stephen The same design matrix can be build of different ways. I think we are talking about the same model but we build the design matrix differently. But that would be another question IMHO

Like others, I am confused by your model. You say A, B, C are not all observed in a given sample - some individuals never make it to C. Whether C is binary or a quantitative trait, if it never happens for some individuals, then no data is observed that would enter the model.  These sample’s estimates would be functions of their observations of A, B.  If you wish to more strongly impact the model estimation process, you can look at like a general likelihood, where you add priors to model coefficients, as well as encode other prior information.

There are two potential sources of bias in this design.

We cannot distinguish correlation from causation.

Imagine two cases. In the first, the disease progression is inducing immune response. Later stages will be associated with the higher gene expression levels. In the second scenario, the disease is caused by overexpression of a gene. Later stages will be also associated with the higher expression.
This is typical for observational studies. But I just want to mention that a special care should be taken during interpretation of the results.

If we are not following our individuals, we cannot distinguish correlation and avoidance.

Let's say the disease is lethal in certain cases; the survival is negatively correlated with the disease stage. Now imagine there is a gene, which causes severe symptoms when highly expressed. On the later stages you will only observe those patients in which the gene was not highly expressed. From that you would conclude that gene expression is decreasing with the disease progression. In reality this gene is very important and causal, you just do not have patients which are alive to see this.
This is similar to Wald's studies of aircrafts (Survival bias).

Researchers from the Center for Naval Analyses had conducted a study
  of the damage done to aircraft that had returned from missions, and
  had recommended that armor be added to the areas that showed the most
  damage.
  Wald proposed that the Navy instead reinforce the areas where the
  returning aircraft were unscathed, since those were the areas that, if
  hit, would cause the plane to be lost.

I think the second point is crucial, and can and will lead to false conclusions. I do not see immediately how to avoid this bias if you study different individuals.
Therefore I suggest that the same individuals are followed for a long time.
There are different approaches you can use later. For example you can have two-stage procedure:

Identify genes which are differentially expressed between healthy (H) and sick (A, B or C).
Build a linear model of disease stage stage ~ gene1 + gene2 + ..., using genes identified at step 1.
Similarly build a linear model of survival as a function gene expression.
It is possible you can use logistic regression to infer probability of transition from stage B to stage C.

This is not the only possible approach of course. As you suggested, Markov model is also applicable.
While following the same individuals it is possible to use continuous-time Markov model for such a study. In this case discrete or continuous measurements are recorded at certain time points, and the model parameters are inferred using maximum likelihood.
I'm not an expert in this field, but I think this paper describing the msm package for R will be useful. It also supports hidden Markov models.

You explained better than myself why this model is not good to find causation :). However in your proposed solutions I have two problems. 1) In certain diseases is not easy to follow the same individuals for a long time (for instance you might have to operate to extract the data, so he/she can't be operated every year). 2) Your two models +survival analysis doesn't take into account that stage C appears on only 30% of patients in stage B.

I agree that it is not always possible, although I do not immediately see a way how to overcome this bias. 2) I expanded question on this (see point #4 and msm).

Ok, I guess there isn't a way to mix the data in linear models :(


text:
How to run HYPHY on multiple files non-interactively?
I would like to test for positive selection in a large set of genes.
I want to have a yes/no answer to the question if gene was evolving under positive selection. Therefore I chose model BUSTED for my analysis, which is implemented in HYPHY.
When I launch HYPHY in the command-line, I get a series of question I  have to answer in the interactive form.
Is there way to perform a batch analysis involving multiple datasets, i.e. many sequence alignments?

What is BUSTED and HYPHY? What have you tried to do that. Does the manual help you?

@Llopis HYPHY is a software for positive seleciton analysis, which is popular but complicated to use.
BUSTED is a model for gene-wide identification of positive selection.
https://doi.org/10.1093/molbev/msv035
I added some background to the question, so it is more clear.
This question is brought up quite regularly by people trying to perform dN/dS analysis. I shared my knowledge here in Q & A form.

This is not a yet an official Q&A if the site doesn't work it will be closed and deleted. Adding ([edit]ing) the context to the question would make the question better for future readers if the site finally lunchs

@Llopis I corrected the question. Do you think it is better now?

It is better but what king of series of questions and answers do you get? Which version of the program are you using ? Under which OS (if that matters)? (I am still not convinced it is on-topic... but let's see)

TL;DR
I made the following bash script to generate an input file for HYPHY.
#!/bin/bash
cat << EOF
inputRedirect = {};
inputRedirect["01"]="Universal"; // genetic code
inputRedirect["02"]="$(readlink -f $1)"; // codon data
inputRedirect["03"]="$(readlink -f $2)"; // tree
inputRedirect["04"]="${3:-All}"; // Test for selection on a branch
inputRedirect["05"]=""; // complete selection

ExecuteAFile (HYPHY_LIB_DIRECTORY+"TemplateBatchFiles/BUSTED.bf", inputRedirect);
EOF

If you want to test all the branches jointly run:
./gen_busted.sh alignment.phy tree.nwk > script.bf

Or for a particular branch:
./gen_busted.sh alignment.phy tree.nwk branchLabel > script.bf

And then run:
HYPHYMP script.bf

Now you can use the script to generate a .bf file for every sequence alignment you have, and run HYPHY for every file.
Longer version
I will describe how to make a similar script for any model, such as MEME, RELAX or PARRIS.
First, find the .bf script which performs the analysis. On GNU/Linux it should be at HYPHY_INSTALLATION_PATH/lib/hyphy/TemplateBatchFiles (it might be different between different OS). In case of BUSTED your script has a name BUSTED.bf.
Now run this model in the interactive mode in order to record all the inputs required by the model. Make sure to specify full paths to all the files.
HYPHYMP HYPHY_INSTALLATION_PATH/lib/hyphy/TemplateBatchFiles/BUSTED.bf

In case of BUSTED the inputs are:

Genetic code.
Codon data (sequence alignment in phylip format).
Tree (in newick format).
Which branches to test.

Now for every input alignment you need to generate a batch file which looks like this:
inputRedirect = {};
inputRedirect["01"]="Universal"; // genetic code
inputRedirect["02"]="/path/to/alignment.phy"; // codon data
inputRedirect["03"]="/path/to/tree.nwk"; // tree
inputRedirect["04"]="All"; // Test for selection on all branches
inputRedirect["05"]="BRANCH1"; // Test for selection on branch1
inputRedirect["06"]="BRANCH2"; // Test for selection on branch2
inputRedirect["07"]=""; // complete selection

ExecuteAFile (HYPHY_LIB_DIRECTORY+"TemplateBatchFiles/BUSTED.bf", inputRedirect);

Now you can create a script which generates the file. The most important thing, don't forget to use full paths in the .bf file.

I am not a moderator, but we should be asking expert questions that would define the on-topicness of the site at this stage of the site. If we start asking about how to do X in each program. We'll end up without expert question.

@Llopis I was sharing similar (but less detailed) information with the three colleagues of mine from my department, some of them where postdocs.
This particular question is brought up over and over again, because despite the fact the the software is used a lot (almost 2k citations), the documentation is very preliminary and scattered.

To document that you could use your private site, or a blog, not necessarily this site. Or you could even contact with the software providers to address that problem if they want people to use their software.

It probably is still up for debate, but what is wrong with answering commonly asked question like: 'How to make tool X do task Y'? The bioinformatics field is probably littered with poorly documented tools abandoned by the authors due to lack of funding. For sure Stack Overflow is full of these kind of questions: https://stackoverflow.com/search?q=how+to

@holmrenser nothing wrong, but the private beta compels us to write great questions for experts: "The private beta gives you the opportunity to get the site off to a great start with expert questions and answers.". I am  saying the answers doesn't seem to me an expert question. But I might be totally wrong :D