text:
How to calculate statistical significance of sequence motifs
I would like to show that certain types of RBP motifs are enriched in RNA editing islands (i.e. clusters of RNA editing). However, I am unsure about how to think about sequence motifs with respect to their occurrence in other genomic features.
I understand how to find the probability of a motif at a location. i.e.
P(sequence is at position i) = P(A)^[A] * P(C)^[C] * P(G)^[G] * P(T)^[T]
where (.) is the base pair and [.] is the number of bps.
Instead what I would like to find is:
P(sequence S is contained in a feature type T) = ???
where feature type T is an gene, intron, editing island, etc. I think I should incorporate length since I will mainly be comparing genes or introns vs. editing islands. Also, I am not sure what to do about the editing islands being located within genes. How can I keep from counting the same motif twice?
Any ideas would be greatly appreciated. Thank you for your time.

It's a bit unclear if you really just want enrichment in editing islands (use Fisher's exact test) or instead want to argue that there's increased enrichment in editing islands versus other features. Please clarify that.

Mainly I just want to show enrichment in the islands. I am familiar with Fisher's Exact Test, however, I wasn't sure how applicable it woudl be in this case as it is difficult to define groups i.e. in/not in editing island. Also the lengths of editing islands vary.

As far as I remember the exact probability computation is an open problem. The reason is that potential motifs can overlap, which makes probability computations for an arbitrary string non-trivial, and depends on the motif.
For example if he have a binary string of four digits, the probability of "01" will be 11/16, while the probability of "11" will be 8/16=1/2.
The simplest approximation is to assume that probability of the motif on every position of the sequence is equal and independent. In such case (this comes from this course):

(Probability of text of length N having at least t occurrences of a k-mer Pattern, A is the number of letters in the alphabet, n = N – t * k)
This approximation is very rough, there are better here.

As Iakov Davydov points out above mathematical calculation is very difficult. Particularly if you do not have uniform nucleotide, dinucleotide etc. frequencies. Thus I'd recommend using a method based on randomisation.
You'd want to randomise only within genes since the likelihood of getting an editing island is likely to be dependent on expression level of the RNA.
You might like to look at what people analysing CLIP data do for these sorts of analyses. See for example Wang et al. Or I have code for this sort of thing in my iCLIPlib. See the module kmer.py. (WARNING: this code is still in development, and is not "released" yet, so I offer no promises about bugs or how easy it would be for your to get working). However, these approaches generally only rank the significance of kmers rather than giving absolute significance values.
Genome Annotation Tester (GAT) might be of use here. You would define genes as your workspace, islands as your annotations, and motif matches as your segments of interest.


text:
Imputing missing genotypes from separate genotyping panels
What is the current standard for imputing missing genotypes between two genotyping panels?  I have two populations genotyped using two different panels (A & B), and I would like to impute all the genotypes in population B for those positions on used in panel A.
I've read the examples for impute2, and I think the closest thing to what I am looking for is this example, "Imputation with one unphased reference panel".
Simply put, I want to provide a list of SNPs, some variant file for population B, and haplotype information from 1,000 Genomes and get imputed genotypes for each SNP in the list.  Is impute2 the state of the art for this?

It depends, if are we imputing whole genome, all chromosomes then impute2 is pretty solid. Unphased panel would give a better result, but slow performance. If we are imputing a region I think beagle is better.

I don't want to impute the whole genome, just certain specified sites.

Given that you mention wanting to use 1000 Genomes as a reference panel for imputing genotypes into your two SNP chip panels, I am going to assume that you are working with human data.
In that case there are several options you can go with:

If your two panels are of European descent, then you are probably best off using the HRC reference panel together with a fast genotype imputation tool such as Beagle 4.1 to impute genotypes in each of your two SNP chip panels separately.
If your panels are not of European descent, then you will likely want to use the 1000 Genomes phase 3 reference panel with Beagle 4.1, Impute2, or Minimac3.

In either case, there are two phasing services available that will do much of the heavy lifting for you 1,2.
The second Wellcome Trust Case-Control Consortium paper performed a cross-imputation analysis as you describe.  I don't see many studies using multiple SNP chip panels.  You will need to take care in your analysis that you are not hit by batch effects from using two different SNP chip panels.
Also, none of these methods will work if the region you are imputing into has too few variants.  I'm not sure what the minimum number of variants is, but if you are using a whole genome genotyping panel of at least 500k SNPs, then you should be ok if you impute a whole chromosome at a time.

Thanks!  These seem like good options.  And yes, I'm working with human data.  I don't know exactly what ethnicity the participants are, and I suspect it is likely a diverse population, so 1000 Genomes makes the most sense to me. You mention Beagle a couple times, is there some reason you have a preference for it?

I don't have enough experience comparing the three programs Beagle 4.1, Impute2, and Minimac3 to really voice a preference.  You should get good quality imputation with any of them, but test for your self, which is easy to do by imputing held out genotypes.  The only caveat to this is that I think Impute2 will likely take a bit longer than the other programs when imputing from a reference panel the size of the HRC.

Beagle worked great and was easy to use.  I highly recommend it

source: Stack Exchange

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

sort: votes


text:
Difference between computational biology, bioinformatics and biostatistics
I find that in many contexts, the terms computational biology, bioinformatics and biostatistics are often treated as functionally equivalent, and yet for students selecting PhD programs and the like the difference could be quite significant. Is there a standard or rigorous definition of these terms and the difference between them?

Biostatistics is clear, but there's a LOT of grey area between computational biology and bioinformatics.

There are many blogs and question on this on other sites, I feel the discussions always becomes rather opinion based and semantic. Here's my personal favorite blog: http://cabbagesofdoom.blogspot.ch/2015/09/is-there-important-difference-between.html

As soon as you read "biology", you can be sure there will be no accepted standard or rigorous definition ;)

I was afraid that this would descend into opinion rather than a having a straightforward answer - but I think this is a real issue that certainly students have no clear answer on, and should be addressed by this community.

NIH has had a working definition since 2000: https://www.kennedykrieger.org/sites/default/files/research_related_files/bioinformatics-def.pdf

There's also a blog post on it here

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


text:
Are there any RepBase alternatives for genome-wide repeat element annotations?
I’m using the RepBase libraries in conjunction with RepeatMasker to get genome-wide repeat element annotations, in particular for transposable elements.
This works well enough, and seems to be the de facto standard in the field.
However, there are two issues with the use of RepBase, which is why I (and others) have been looking for alternatives (so far without success):

RepBase isn’t open data. Their academic license agreement includes a clause that explicitly forbids dissemination of data derived from RepBase. It’s unclear to what extent this is binding/enforceable, but it effectively prevents publishing at least some of the data I’m using and generating. This is unacceptable for open science.

Subordinate to this, the subscription model of RepBase also makes it impossible to integrate RepBase into fully automated pipelines, because user interaction is required to subscribe to RepBase, and to provide the login credentials.

RepBase is heavily manually curated. This is both good and bad. Good, because manual curation of sequence data is often the most reliable form of curation. On the flip side, manual curation is inherently biased; and worse, it’s hard to quantify this bias — this is acknowledged by the RepBase maintainers.

Were you only asking about defined repeat libraries? I interpreted it slightly more broadly as about tools used to build the libraries also (which becomes relevant when genomes from new taxa are sequenced)

@Chris_Rands Both (libraries and tools). Your answer is spot-on.

The goal is to build the annotated library of repetitions, or to mask repetitive parts of a genome?

@KamilSJaron I’m working with TEs, so I need the annotated library, not (merely) a repeat masked sequence.

Ouch, but for TEs just a subset of repetitive regions, there are specialised tools to annotate them (like DNApipeTE and REPET). Maybe you could specify it in the question.

@KamilSJaron Nice, and this might be worth an answer. I’ll also update teh question. That said, I am also asking for repetitive elements beyond TEs.

Also are you sure that RepBase is against publications that derive information based on their data? Could it be that they just don't want you sharing the raw data and files they provide. Since they have a good number of citations in different fields: https://scholar.google.de/scholar?um=1&ie=UTF-8&lr&cites=11574259945967474319

@story They literally say so in the academic user agreement that I link to. Here’s the relevant quote: “You agree NOT to make the Repbase (or any part thereof, including Repbase Reports, Repeat Maps and other derived materials, modified or not) available to anyone outside your research group.” Emphasis mine. In fact, another clause in the agreement technically even forbids me from signing it because my institute requires public data deposition, so I’m probably not allowed to sign such agreements.”

Ya that seems to agree with my previous statement. I guess my point is what exactly did you need to share (based on your original post) that would be considered from their database? I feel like this wouldn't include counts of features but sequences might be an issue.

@story I need to potentially share all data that was used/generated in my analysis. This particularly includes the specific repeat annotation I used, which are derived from RepBase, as well as potentially sequence data from these repeats.

This might be an old question, but someone is trying to set-up a new, open, alternative to repBase (which i snow going full commercial), or at least that is how I perceive it: https://twitter.com/TransposableMan/status/1060519887897067521

Dfam has recently launched a sister resource, Dfam_consensus, whose stated aim is to replace RepBase. From the annoucement:

Dfam_consensus provides an open framework for the community to store both seed alignments (multiple alignments of instances for a given family) and the corresponding consensus sequence model.

Both RepeatMasker and RepeatModeler have been updated to support Dfam_consensus.
I haven’t tried it yet but it looks promising.

For pre-existing reliabe TE libraries it is a bit of a mess, because not everybody deposits the species-specific TE libraries to a database like RepBase. And as far as I know DFAM contains only human resources, or am I wrong?
As for de novo generation of species-specific TE libraries (which should be done for any species not already present in eg. RepBase):
There is no "gold-standard" how to tackle this best.
In principle one has to think about two main parts
-repeat detection
-annotation
For repeat detection I would recommend using a combination of two things (which is necessary, because TE copies might miss in the assemblies as repetitive regions tend to be difficult to assemble and thrown away in the final assembly).
I) Repeat detection from raw reads (as with e.g. DNApipeTE or tedna or RepeatExplorer). For me, DNAPipeTE worked quite nicely, but everything has pros and cons.
II) Repeat detection from assemblies (as with e.g. REPET or as mentioned before RepeatModeler)
Then the annotation of these repeats is tricky too, because most methods are relying on homology between the de novo TEs and the TEs from some (probably distantly) related species. But some programs take also structure into account (like REPCLASS). REPET can do both detection and annotation, but is a pain to get to run.
I would recommend using some programs to do de novo repeat detection on your species of interest on both the raw reads and assembly, clustering these libraries together (with e.g. uclust and 95% identity) and then run an annotation with homology and structural identification.
Probably the programs will not give you complete, full-length TEs but rather consensus sequences of several copies from TE families.
If you want, you could search all copies of one family, extract them from the contigs plus boundaries and align them manually and curate boundaries manually. Then extend boundaries if not hitting the surrounding (non-alignable) regions or landmarks of TEs like LTRs or TIRs or so. But this is very time consuming if you only want to compare TE abundance between species for example, I would not do this and rather compare the abundance using read coverage (as in Bast et al. 2016). Depends all on the questions you want to ask.

You could use RepeatScout, which has defined repeat libraries for a limited number of species (including human, mouse, and rat). If your taxon is not represented, you can also do de novo repeat prediction with RepeatScout to build your own library to feed to RepeatMasker. The RepeatScout publication includes some comparisons with RepBase. Another related tool is RepeatModeler, which wraps RepeatScout with RECON and some other programs, and shares authors with the RepeatMasker team.
On the plus side RepeatScout/RepeatModeler are open source and do not use manual curation, meeting your criteria. On the negative, I'm not sure exactly how RepeatModeler and the component tools are maintained. The RepeatScout web and github pages have not been updated for several years, although the RepeatModeler page shows its latest release was in 2017. Anyway, I know that some combination of RepeatScout/RepeatModeler have been used to annotate repeats for some fairly recent newly sequenced genomes, e.g. for cichlids, coelacanth, and Darwin's finch, so I think it's fair to say this kind of approach is accepted in the field, at least for vertebrate genome projects.

AFAIK Dfam and Repbase are currently the two best sources of (a variety of) TE sequences.
In my genome annotations I have used RepeatModeler+RepeatMasker and then later used Repbase+tblastx and Dfam+nhmmer to classify them.
The classification process in my pipeline PhyLTR (https://github.com/mcsimenc/PhyLTR) is based on Dfam and Repbase. The process I used for LTR identification is

Putative ID with LTRHarvest (based on structural sequence characteristics)
Classification-by-homology to Repbase and Dfam
Removal of elements without homology to sequences in Repbase or Dfam.

This results in a set of LTR-Rs which are full-length and have evidence that they are LTR-Rs.

+1 for taking issue with RepBase.
I use the annotations from the Hammell Lab GTFs that they put out with TEtoolkit. It is similar to what you described to be using so this may be a redundant and useless answer but from the digging I’ve done they do seem to be comprehensive and well curated (for Drosophila, at least).

I know this question is a bit old, but this is still an issue for a lot of researchers not being able to access RepBase. It seems now that the most recent version of RepeatMasker is dependent upon RepBase for full functionality if masking anything other than human (currently DFAM only has human models). I recently discovered a de novo repeat masking approach called REpeat Detector (Red).  This might be a solution for some looking to mask repeats on a genome assembly for annotation.  The paper is here. I also then wrote a wrapper around Red to make it a bit more easy to soft-mask a genome, which you can find here.
One of the limitations with Red is that the repeats are not classified, so they are only identified. You would have to use some of the other tools mentioned above to try to classify them.

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


text:
What's the difference between VCF spec versions 4.1 and 4.2?
What are the key differences between VCF versions 4.1 and 4.2?
 It looks like v4.3 contains a changelog (specs available here) but earlier specifications do not.
This biostar post points out one difference: the introduction of Number=R for fields with one value per allele including REF — can anyone enumerate the other changes between these two versions?

This is easy to check, you can download both specs in .tex format and do diff.
Changes to the v4.2 compared to v4.1:

Information field format: adding source and version as recommended fields.
INFO field can have one value for each possible allele (code R).
For all of the ##INFO, ##FORMAT, ##FILTER, and ##ALT metainformation, extra fields can be included after the default fields.
Alternate base (ALT) can include *: missing due to a upstream deletion.
Quality scores, a sentence removed: High QUAL scores indicate high confidence calls. Although traditionally people use integer phred scores, this field is permitted to be a floating point to enable higher resolution for low confidence calls if desired.
Examples changed a bit.

Doing a raw diff against the spec files seems pretty silly, ideally theres documented differences somewhere, but thanks for summarizing this!

Would also be worth generating a comparison with 4.3


text:
How to select high quality structures from the Protein Data Bank?
Models of structures deposited in the Protein Data Bank vary in the quality, depending both on the data quality and expertise and patience of the person who built the model. Is there a well-accepted subset of the PDB entries that has only "high quality" structures? Ideally these structures would be representative for classes of proteins in the whole PDB.
based on a real question from biology.SE

There is a very nice database, pdbcull (also known as the PISCES server in the literature). It filters the PDB for high resolution and reduced sequence identity. It also seems to be updated regularly. Depending on the cut-offs, you get between 3000 and 35000 structures.
If you are specifically interested in rotamers, you may want to look at top8000 instead, where they have checked for high resolution, and good MolProbity scores. They also provide a rotamer database.
PDB also provides their own clustering. They first cluster the sequences, and then extract a representative structure for each one, based on the quality factor (1/resolution - R_value). This has the advantage of being comprehensive, but you will have bad structures when no good ones were ever obtained.

If you choose to perform your own culling of the PDB, resolution is probably the first thing you'll want to look at, which as Davidmh mentions is the main selection criteria for PISCES. High quality structures will also have better R-factor values. You can also give preference based on experimental technique, in descending order of quality:
Neutron diffraction, X-ray diffraction, solution/solid state NMR, electron microscopy/crystallography, fiber diffraction, solution scattering.

yes, these are criteria from PISCES, although they look to be quite simplistic. R-factor is used but the gap between R-factor and Rfree is ignored. Resolution (pressumably d_min) is the only criterium of data quality (data completeness is ignored). No geometry validation.

R-factor and R-free discrepancy is a great thing to look at, but in my experience the number of structures reporting both in a standardized way is rather small. BioJava only recently made them both available, IIRC, and that's typically the tool I use.

80.7% of PDB structures reports R-free, slightly less than 86.8% that reports resolution.


text:
Tumor purity/contamination/admixture estimation
Can anyone recommend a good tool for estimating the tumor content given a matched tumor and normal file for DNA NGS whole genome sequencing data or whole exome data?
Is it possible to estimate this without a normal sample as well?

N.B., I've edited your post to include the word "admixture", which is what you should search for the find relevant tools.

Great, might as well make it purity/contamination/admixture then =). Many names for same meaning.

Yeah, but you'll get a lot more relevant hits with "cancer admixture estimate" :)

I would suggest to change "admixture", which is common in population genetics for genetic exchange between populations to "cell admixture".

I have previously estimated tumour purity with the EXPANDS an inferred tumour heterogeneity program which is designed to calculate the number of clonal subpopulations in matched tumour/normal samples.  The purity is essentially the size of the largest subpopulation identified in that sample - this is discussed in the programs FAQ.  In addition to a matched tumour / normal exome or genome sequencing you'll also need match somatic copy number as a seg file as input.  Other programs also exist for this sort of inferred heterogeneity analysis - some of which may also give you a measure of purity are: Absolute, ThetA, SciClone, CHAT, PyClone and Canopy. A more complete list looks to be here.
The only other thing I'd suggest is have an estimate via some orthogonal measure of the purity to try at least judge how the different methods are performing with you data, e.g. you might have clonality indications from cytogenetics/FISH or histology work, or perhaps FACS.  Picking samples known to be pure / impure via one or more of these might help you get a handle on how well various methods are performing.
As regards estimation without a normal sample the only program which looks like it might help is QuantumClone, but that requires more than one tumour sample form said patient either spatially or temporally distinct to perform the analysis.
Also If you have low coverage DNA sequencing CNAnorm looks like it could be useful with just a single sample.

It's usually CNV callers that make use of Tumour/Normal WGS pairs to estimate purity. It can also be done with WES (exome) Tumour/Normal pairs.
There are several tools out there, I have some experience with the one written by Illumina (public on Github):
https://github.com/Illumina/canvas
It requires realigning things with bowtie2, so I don't think it can take existing data aligned with other aligners. Not sure why.

Where are you getting the bowtie2 requirement from? The readme examples all use isaac bams

Section 5.1 of the pdf documentation (CanvasBin): Canvas has been validated using BAM files created using the Isaac and Bowtie aligners. I guess it means that it may work with other aligners, but hasn't been validated.

Thanks for your answer. Since this is primarily a CNV caller, do you think that mean it requires rearrangements to estimate admixture amount?

Celluloid emits copy number profiles and tumour purity / ploidy information.  There's a nice tutorial https://github.com/mathieu-lemire/celluloid_0.11_tutorial. One thing to keep in mind is that celluloid will find multiple solutions and taking advantage of its plotting tools is essential to determining which solution is correct. Typically the first solution is correct, but occasionally it may not be.
Another tool similar is TITAN, which may also require hand annotated the ideal solution.


text:
Using shells other than bash
As someone who's beginning to delve into bioinformatics, I'm noticing that like biology there are industry standards here, similar to Illumina in genomics and bowtie for alignment, many people use bash as shell.
Is using a shell besides bash going to cause issues for me?

I would adjust the examples you provided. Illumina is a standard for short reads, but there are many genomics labs working mostly with PacBio or Nanopore. Bowtie is hardly a standard. Even versions 1 and 2 are very different.

@burger what do you suggest then?

No suggestion. Although I agree with all the answers so far, bioinformatics is not good with standards. Even something like a SAM/BAM file that is technically a properly defined standard that almost everyone in genomics uses has many fields that are treated differently, causing issues for a lot of tools.

A statement "this is not meant to be opinionated" doesn't help much with such a broad question like this. Do you have a particular application that you would like to use a shell for, or an indication of what "industry" you're interested in?

@burger: Do you have specific problematic SAM/BAM fields in mind?  You could raise issues at https://github.com/samtools/hts-specs/issues or at least this suggests another question to ask here…

@JohnMarshall I don't think there is a "bug" with the SAM/BAM standard. It's just that it's open-ended and different tools require different fields. I have had to modify my BAM files many times in the past because some tool expected it in a slightly different format. Technically, it's still a valid BAM before and after, but one is compatible and one isn't. If you have a BAM, you have no idea if it will work with a tool that requires a BAM file.

@burger: If you want this situation to improve, you'll have to say which particular fields you have had to modify and what the various tools' expectations were. If you do this, specs can be clarified, tools can be modified, and everyone's bioinformatics pipelines can run a little more smoothly. Otherwise, it's just FUD.

Bioinformatics tools written in shell and other shell scripts generally specify the shell they want to use (via #!/bin/sh or e.g. #!/bin/bash if it matters), so won't be affected by your choice of user shell.
If you are writing significant shell scripts yourself, there are reasons to do it in a Bourne-style shell.  See Csh Programming Considered Harmful and other essays/polemics.
A Bourne-style shell is pretty much the industry standard, and if you choose a substantially different shell you'll have to translate some of the documentation of your bioinformatics tools.  It's not uncommon to have things like

Set some variables pointing at reference data and add the script to your PATH to run it:
export FOO_REF=/path/to/stuff
export PATH=/path/to/foo-x.y:$PATH
foo blah blah

These will typically be shown in Bourne-shell syntax.  By using a different shell you have to translate the export commands to your local syntax, and especially PATH munging is somewhat shell-dependent.
If you're experienced in Unix, this will be only a minor niggle.  If you're a beginner, IMHO this will add a non-negligible amount of friction on top of all the other things you're learning.

Do not use #!/bin/bash in the shebang. Having Bash installed in a nonstandard location is common enough that doing so will break often. Use #!/usr/bin/env bash instead, it should have no disadvantage.

TL;DR
SH adheres to an official industry standard, but it is not suitable for scientific computing. Bash is considered an informal standard (e.g., by Google). Bash 3 is preferable in most of situations in the world of bioinformatics.
Long answer
As already described in other answers, SH (/bin/sh, plain Bourne shell, the original UNIX shell) should fully adhere to POSIX which is a real industry standard. However, SH is too limited for scientific computing since many key features were incorporated later in SH successors, especially in Bash (/bin/bash, Bourne Again Shell): set -o pipefail, [[ ... ]], or process substitutions < () to name at least few.
In practice, it is much more difficult to write "safe" scripts in pure SH and only shell experts are usually capable to prevent unexpected behavior. For example, it may be hard to ensure that no command in a pipeline failed in the middle of computation. For Bash, various easy-to-follow defensive-programming recommendations have been developed and they should prevent all these problems. From this reason, many computer scientists, software engineers and companies use Bash as a kind of a standard. For instance, Google internal policy allows only Bash for writing any shell scripts.
Even though we cannot expect that Bash is present on completely every Unix machine (e.g., on mobile devices as @terdon pointed out), a vast majority of *nix machines used for scientific computation should have it. We also should be aware of the fact that Bash can be slower than SH and that it has recently suffered from  major security issues. Moreover, various Bash versions exist and scripts working on modern Linux machines with Bash 4 might not work on OS X, which is still based on Bash 3.
To sum up, Bash 3 is probably be the most reasonable choice for scientific computing.
Edit
I addressed the comments from @terdon and @John Marshall. In particular, I added an explanation why Bash is more suitable for scientific computing than SH (in my opinion).

Bash is not present in every Unix machine, sh is and that is not the same. Yes, Linux tends to have /bin/sh pointing to bash, but Linux is not Unix and, anyway, even in Linux /bin/sh is not always bash (Debian based systems use dash instead, for example). You can safely expect the Bourne shell (sh) to be present on a POSIX-compliant system, but not necessarily the Bourne again shell (bash).

@terdon Could you provide any reference, please? According to https://wiki.debian.org/Bash, bash is the default shell on Debian. Do you know about any (modern) *nix distro where bash would not be installed?

@terdon I will answer my question – e.g., FreeBSD. https://www.freebsd.org/doc/en/articles/linux-users/shells.html says that "Bash is not included in the default installation". Do you have an example of a Linux distribution without bash?

Some (all?) embedded Linux systems won't have bash and will instead have busybox sh. The main issue is that people tend to think that sh and bash are the same thing but they aren't. They are similar and bash is an extension of sh, but they are not the same.

@Karel: Asking about the “default shell” is ambiguous. As per https://wiki.debian.org/Shell, these days on Debian the default /bin/sh is dash while the default login shell (as listed in /etc/passwd) remains /bin/bash. This means that portable shell scripts that identify themselves with #!/bin/sh need to restrict themselves to POSIX shell facilities, while scripts that want to use bash extensions need to use #!/bin/bash.  This got tidied up the hard way a few years ago when various distributions switched to dash for /bin/sh…

@terdon @John Marshall Thank you for your comments. Compared to bash, I consider "pure" sh very limited and inappropriate for scientific computing, in particular because of some missing, but very important, features such as set -o pipefail or [[ ... ]]. My experience is that sh scripts can be very susceptible to unexpected behavior (unless the developer is a shell expert, which is usually not the case in bioinformatics). Several good and simple defensive programming strategies for scientific computing exist for bash.

This is why I would like to know whether /bin/bash might not return anything, or return a non-bash shell (I have seen such a problem only once with some obscure bioinformatics distribution).

I wouldn't do "scientific computing" in a shell no matter what shell it was. The shell should be used for, at most, handling the plumbing for basic utilities and applications. Computing should be handled by utilities and applications designed for those tasks.

@Kusalananda  How do you do scientific computing without shell? I believe that you use it at least for running your programs. If it is so, do you agree that the way how it handles errors matters?

@Karel I wouldn't do computing of any kind without a shell, but not in (with) a shell.

I’m a bit puzzled why this answer recommends the ancient Bash 3 instead of Bash 4, which itself is now almost 10 years old (announced in 2009). Bash 3 lacks crucial features such as associative arrays so it’s a severe restriction. True, macOS still ships with Bash 3 but so what? macOS is generally known to lag behind in its unix tools (and even for its Ruby and Python). Also, nitpick: it’s “Bash”, not “BASH”.

@KonradRudolph Thank you for the comment. I fixed the capitalization issue. Concerning Bash 4, I completely agree that it has many useful features. However, if it cannot be used on a substantial proportion of machines, it's a fatal problem. Whereas Python 3 can be easily installed (e.g., using Conda), upgrading Bash is complicated and easily results in serious problems. As for the associative arrays, the Google standard says the following: "If you find you need to use arrays for anything more than assignment of ${PIPESTATUS}, you should use Python."

@Karel I have a few choice word for the Google coding guidelines, none of which are acceptable in polite company. Anyway, upgrading Bash is actually trivial. Replacing it for the login shell may not be but in practice that’s unnecessary: on macOS you specify the shell in the terminal app, and other systems ship with Bash 4.

strongly agree with @Kusalananda, trying to write your pipelines entirely in shell is a mistake. There are plenty of workflow frameworks; I am partial to Nextflow, and a lot of my peers use Snakemake. Entirely shell-based pipelines quickly become unmanageable, overly complex, confusing to understand, and extremely difficult to debug. If you must use Bash, then you should aim for POSIX compliant implementations.

additionally, a lot of horrible Bash code for simpler scripts can be accomplished better with Makefiles. For beginners, you should aim to learn how to use these after you've become comfortable with basic shell scripting.

The Open Group Base Specifications Issue 7
IEEE Std 1003.1™-2008, 2016 Edition, or "The POSIX Standard" for short, is the standard that defines the interfaces and utilities provided by a Unix system.  Among these is the command line shell language and tools (see "Shell & Utilities" in the main index on the page linked above).
As far as I know, there is no shell that implement exactly what's specified by the standard, but both bash and ksh93 does a pretty good job of adhering to the standard along with their own, sometimes conflicting, extensions. The ksh93 shell in particular has had a big impact on the past development of the POSIX shell specification, but future POSIX specifications may borrow more from bash due to its wide use on Linux.
The bash shell is pretty much ubiquitous on Linux systems, and may be installed on all other Unices too. ksh93 is also available for most Unices but is usually not installed by default on Linux. ksh93 is available by default on at least macOS (as ksh) and Solaris.
If you are concerned about portability when writing a shell script (which is IMHO a good thing to be concerned about), you should make sure that you use only the POSIX utilities and their POSIX command line flags, as well as only use POSIX shell syntax. You should then ensure that you script is executed by /bin/sh which is supposed to be a shell that understands the POSIX specification.  /bin/sh is often implemented by bash running in "POSIX mode", but it may also be dash, ash or pdksh (or something else) depending on what Unix you are using.
For a Linux user, the most difficult bit in writing a portable script is often not the shell per se, but the multitude of non-standard command line flags provided by the GNU implementation of the many shell utilities.  The GNU coreutils (basic shell utilities) may, like bash, be installed on all Unices though.
Also note that bash, when running in POSIX mode (either when invoked as /bin/sh or with its --posix command line flag), is not strict about its POSIX conformance and may accept some syntax extensions to the POSIX standard.

I would not say bash as a "standard", but it is indeed likely to be the most widely used unix shell and available by default on most modern unix/linux distros. There are a few other more convenient shells like zsh that are broadly compatible with /bin/sh, but they are not as widely available. There is also C-shell and in particular its open-source implementation tcsh. C-shell is quite different from bash. Over ten years ago, I saw it was used from time to time, but nowadays, I rarely see its use, except by programmers from older generations.

The generic command sh is quite literally an industry standard, a POSIX standard, to be precise (IEEE 1003.2 and 1003.2a, available for purchase for hundreds of dollars at various websites). In theory, any script that starts with #!/bin/sh should conform to this standard. In practise, most Linux systems have a shell that is close to this standard, but has a few quirks and extensions.
Problems crop up when these quirks and extensions become standard practise in shell scripts. The Debian operating system changed to dash as their sh shell to encourage people to stop using "bashisms" in shell scripts that didn't specify a particular shell, i.e. those that began with #!/bin/sh. The dash shell tries to be as standards-compliant as possible:

dash is the standard command interpreter for the system.  The current version of dash is in the process of being changed to
       conform with the POSIX 1003.2 and 1003.2a specifications for the shell.  This version has many features which make it
       appear similar in some respects to the Korn shell, but it is not a Korn shell clone (see ksh(1)).  Only features designated
       by POSIX, plus a few Berkeley extensions, are being incorporated into this shell.  This man page is not intended to be a
       tutorial or a complete specification of the shell.

I'm not familiar with the differences, and generally try to stick to the sh man pages to instruct me with regards to correct standards-compliant shell scripts.

Note that sh is not bash. Even on systems whose /bin/sh points to bash, being invoked as sh changes bash's behavior and causes it to run in POSIX-compliant mode. The "real" sh shell (bourne shell) is something else again and not the same as bash (bourne again shell).

In Debian the default interactive shell, i.e. the one you'll be using on the command-line is bash https://wiki.debian.org/Shell yes /bin/sh will be symlinked to /bin/dash but the one people will be using live will be bash.


text:
What is the fastest way to calculate the number of unknown nucleotides in FASTA / FASTQ files?
I used to work with publicly available genomic references, where basic statistics are usually available and if they are not, you have to compute them only once so there is no reason to worry about performance.
Recently I started sequencing project of couple of different species with mid-sized genomes (~Gbp) and during testing of different assembly pipelines I had compute number of unknown nucleotides many times in both raw reads (in fastq) and assembly scaffolds (in fasta), therefore I thought that I would like to optimize the computation.

For me it is reasonable to expect 4-line formatted fastq files, but general solution is still prefered
It would be nice if solution would work on gzipped files as well

Q : What is the fastest way (performance-wise) to compute the number of unknown nucleotides (Ns) in fasta and fastq files?

Presumably, you might want the total number of nucleotides as well as the number of Ns (to give proportion of Ns). This is a fairly fast implementation that provides this plus a few other metrics for fasta: https://github.com/ENCODE-DCC/kentUtils/blob/master/src/utils/faSize/faSize.c

For FASTQ:
seqtk fqchk in.fq | head -2

It gives you percentage of "N" bases, not the exact count, though.
For FASTA:
seqtk comp in.fa | awk '{x+=$9}END{print x}'

This command line also works with FASTQ, but it will be slower as awk is slow.
EDIT: ok, based on @BaCH's reminder, here we go (you need kseq.h to compile):
// to compile: gcc -O2 -o count-N this-prog.c -lz
#include <zlib.h>
#include <stdio.h>
#include <stdint.h>
#include "kseq.h"
KSEQ_INIT(gzFile, gzread)

unsigned char dna5tbl[256] = {
    0, 1, 2, 3,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 5, 4, 4,
    4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 0, 4, 1,  4, 4, 4, 2,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 4, 4, 4,  3, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 0, 4, 1,  4, 4, 4, 2,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 4, 4, 4,  3, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4
};

int main(int argc, char *argv[]) {
    long i, n_n = 0, n_acgt = 0, n_gap = 0;
    gzFile fp;
    kseq_t *seq;
    if (argc == 1) {
        fprintf(stderr, "Usage: count-N <in.fa>\n");
        return 1;
    }
    if ((fp = gzopen(argv[1], "r")) == 0) {
        fprintf(stderr, "ERROR: fail to open the input file\n");
        return 1;
    }
    seq = kseq_init(fp);
    while (kseq_read(seq) >= 0) {
        for (i = 0; i < seq->seq.l; ++i) {
            int c = dna5tbl[(unsigned char)seq->seq.s[i]];
            if (c < 4) ++n_acgt;
            else if (c == 4) ++n_n;
            else ++n_gap;
        }
    }
    kseq_destroy(seq);
    gzclose(fp);
    printf("%ld\t%ld\t%ld\n", n_acgt, n_n, n_gap);
    return 0;
}

It works for both FASTA/Q and gzip'ed FASTA/Q. The following uses SeqAn:
#include <seqan/seq_io.h>

using namespace seqan;

int main(int argc, char *argv[]) {
    if (argc == 1) {
        std::cerr << "Usage: count-N <in.fastq>" << std::endl;
        return 1;
    }
    std::ios::sync_with_stdio(false);
    CharString id;
    Dna5String seq;
    SeqFileIn seqFileIn(argv[1]);
    long i, n_n = 0, n_acgt = 0;
    while (!atEnd(seqFileIn)) {
        readRecord(id, seq, seqFileIn);
        for (i = beginPosition(seq); i < endPosition(seq); ++i)
            if (seq[i] < 4) ++n_acgt;
            else ++n_n;
    }
    std::cout << n_acgt << '\t' << n_n << std::endl;
    return 0;
}

On a FASTQ with 4-million 150bp reads:

The C version: ~0.74 sec
The C++ version: ~2.15 sec
An older C version without a lookup table (see the previous edit): ~2.65 sec

based on Devon Ryan's benchmarking you are the fastest so far...

Can you add #include <inttypes.h> to you first C example?

I can’t not upvote SeqAn, but I’m seriously concerned that your C++ code is so much slower than your C code. They should be virtually identical. That said, you didn’t specify std::ios::sync_with_stdio(false), which might account for some of the difference; and the loop could be simplified by using iterators, though I’d be surprised if this impacted the speed on any setting above -O1.

@KonradRudolph sync_with_stdio helps little when reading from a normal file. I don't know the internal of seqan, so I can't explain why it is slower. That said, I wouldn't worry too much about seqan performance. Gzip decompression is much slower than parsing fastq. @DevonRyan added.

The difference is in readfq. That beast is fast. I know, I tried to replace Heng's code by something less ugly in C++ and couldn't get anywhere near it in simple read performance. Best I could come up with in clean C++ was a speed factor of 2x slower.

you need to add -lz to your compile flags (at least for the current Ubuntu LTS) gcc -O2 -o count-N this-prog.c -lz

Are your numbers based on flushing the cache between samples? Otherwise, those numbers will almost certainly be incorrect.

Out of curiosity I just checked the performance with a straightforward implementation in seqan3 and to my great disappointment it’s still a factor 3 slower than kseq.h. It’s of course entirely possible that I did something wrong but I don’t think so.

5 hours and no benchmarks posted? I am sorely disappointed.
I'll restrict the comparison to just be fasta files, since fastq will end up being the same. So far, the contenders are:

R with the ShortRead package (even if not the fastest, certainly a super convenient method).
A pipeline of grep -v "^>" | tr -cd A | wc -c
A pipeline of grep -v "^>" | grep -io A | wc -l
A pipeline of seqtk comp | awk
Something custom in C/C++

The R code can't actually count N records for some reason when I try it, so I counted A for everything.
I'll exclude all python solutions, because there's no conceivable universe in which python is fast. For the C/C++ solution I just used what @user172818 posted (everything I tried writing took the same amount of time).
Repeating 100 times and taking the average (yes, one should take the median):

1.7 seconds
0.65 seconds
15 seconds
1.2 seconds
0.48 seconds

Unsurprisingly, anything in straight C or C++ is going to win. grep with tr is surprisingly good, which surprises me since even though grep is very very optimized I still expected it to have more overhead. Piping to grep -io is a performance killer. The R solution is surprisingly good given the typical R overhead.
Update1: As suggested in the comments

sum(x.count('A') for x in open('seqs.fa','r') if x[0] != '>') in python

This yields a time of

1.6 seconds (I'm at work now, so I've adjusted the time as best I can to account for the different computers)

Update2: More benchmarks from the comments

awk -F A '!/^>"/ {cnt+=NF-1}END{print cnt}' foo.fa
perl -ne 'if(!/^>/){$count += tr/A//} END{print "$count\n"}' foo.fa

These yield times of:

5.2 seconds
1.18 seconds

The awk version is the best hack I could come up with, there are likely better ways.
Update 3: Right, so regarding fastq files, I'm running this on a 3.3GB gzipped file with 10 repetitions (this takes a bit of time to run), so I'm not going to initially limit tests to commands that I can trivially modify to handle compressed files (after all, uncompressed fastq files are an abomination).

seqtk fqchk
sed -n '2~4p' <(zcat foo.fastq.gz) | grep -io A | wc -l
sed -n '1d;N;N;N;P;d' <(zcat Undetermined_S0_R1_001.fastq.gz) | tr -cd A | wc -c
The C example from user1728181 (the comparison to C++ is already provided there, so no need to include it).
bioawk -c fastx '{n+=gsub(/A/, "", $seq)} END {print n}' foo.fastq.gz

The average times are:

1 minute 44 seconds
2 minutes 6 seconds
1 minute 52 seconds
1 minute 15 seconds
3 minutes 8 seconds

N.B., these will often not handle fastq files with entries spanning more than 4 lines. However, such files are the dodo birds of bioinformatics, so that's not practically important. Also, I couldn't get ShortRead to work with the compressed fastq file. I suspect there's a way to fix that.

Nice benchmark, what is the size of your test fasta file? Also, why not add Python if you have time, yes it will not be fastest, but interesting to see for comparison, I would try: sum(x.count('A') for x in open('seqs.fa','r') if x[0] != '>'), probably BioPython will be slower for this

@Chris_Rands I used dm6 (uncompressed), so it's not terribly large. I figured that python would be so self-evidently slower that I didn't want to bother, but I suppose it wouldn't hurt to add.

@Chris_Rands I've updated with a python benchmark. It does as well as R, which is better than I would have guessed.

Thanks! I'm not so surprised since python's str.count is implemented in C. I also feel there should be a fast awk only solution, but I'm not able to figure it out

Perl: perl -ne 'if(!/^>/){$count +=  tr/A//} END{print "$count\n"}' dm6.fa

Okay, now that the question has been updated, you can modify these to work for the general FASTQ case (multi-line sequence + quality string, with '@' possibly appearing in the quality line(s)).

My awk-fu is sufficient to come up with a working example, but not good enough to come up with one that's not painfully slow.

@gringer There goes me being productive today...

Can you try my bioawk -c fastx '{n+=gsub(/A/, "", $seq)} END {print n}' ? I'm curious to see how it compares with the other solutions on your data.

@bli: Yup, I have 5 others running at the moment...

The R code can count Ns once you write it correctly. ;-) (My bad)

The awk and perl versions do not give the same results on my test fasta file: The awk version find one more "A" in C. elegans genome than the perl version. I guess it depends whether one of the sequence starts or ends with "A".

I've added a fastq comparison with some methods (including bioawk).

@bli There should just be a -1 in the awk code. Among the many reasons to use bioawk instead!

I think the zcat is not necessary with bioawk. Does it improve the speed?

@bli It's not needed and that's a mistake on my part (I didn't actually use zcat in the benchmark). Thanks for catching that!

Why are you counting the As and not the Ns? Or am I missing something?

Because the original ShortRead variant couldn't handle N and the performance comparisons are the same regardless. "The R code can't actually count N records for some reason when I try it, so I counted A for everything."

Yes, sorry, I missed that when I first read the answer.

You say: "Repeating 100 times and taking the average (yes, one should take the median):" But the average is the mean not median - so which one did you use? I appreciate the median is less susceptible to skew but in this case surely run time is normally distributed (in which case mean is fine) unless you're running on a system with very high instantaneous CPU or IO load from other users which is creating contention.

The mean, not the median, thus the caveat. Yes, run time should be normally distributed, but occasional weird IO spikes and such do happen.

Did you flush system cache between all these tests?

Nah, in practice the files were in the cache for all of the runs.

You should flush cache between tests. Your numbers will be misleading, otherwise, (well, wrong, really) regardless of whether you draw a mean or a median from test results.

I think that this should be pretty fast:
FASTA:
grep -v "^>" seqs.fa | tr -cd N | wc -c

FASTQ:
sed -n '1d;N;N;N;P;d' seqs.fq | tr -cd N | wc -c

See this answer on SO about how to count characters in BASH using different approaches.

Your FASTQ answer unfortunately has the same problem, fundamentally, as Iakov’s. That is, FASTQ records have a variable number of lines.

I know that the original standard admits FASTQ line splitting. However, have you ever seen such a file in practice?

I haven’t, but I’ve never myself worked with Sanger sequencing/PacBio data. I think you might encounter them there. For short read data I guess the “record = 4 lines” assumption is safe.

Long read data do not use fastq. Currently PacBio uses bam (yes, really bam) and ONT fast5 (h5 archive). I think it is very reasonable to assume 4 line fastq.

FASTQ
As it was pointed out, fastq can be complicated. But in a simple case when you have four lines per record, one possible solution in bash is:
sed -n '2~4p' seqs.fastq | grep -io N | wc -l

sed -n '2~4p' will print every fourth line
grep -o N will output a line with N for every matching symbol
wc -l will count the lines

I suspect this python approach will work faster:
cat seqs.fastq | python3 -c "import sys; print(sum(line.upper().count('N') for line in sys.stdin))"

FASTA
Coreutils:
grep -v "^>" seqs.fasta | grep -io N | wc -l

Python:
cat seqs.fasta | python3 -c "import sys; print(sum(line.upper().count('N') for line in sys.stdin if not line.startswith('>')))"

Now you’re counting lines containing Ns, not Ns. Also, FASTQ records can span more than four lines. :-(

@KonradRudolph with -o grep outputs individual symbols, I will think how to correct fastq code

FASTQ is tricky because @ can also start the quality string line; see http://chat.stackexchange.com/transcript/message/37809200#37809200 and following lines. For reference, here’s a Perl script that parses FASTQ and trims the mate information off of the identifier. It can probably serve as a starting point; I think this is very close to the minimal code: https://gist.github.com/klmr/81a2b86cd93c706d28611f40bdfe2702

@KonradRudolph, I updated the answer to indicate that this works only in the simple case

Using bioawk
With bioawk (counting "A" in the C. elegans genome, because there seem to be no "N" in this file), on my computer:
$ time bioawk -c fastx '{n+=gsub(/A/, "", $seq)} END {print n}' genome.fa
32371810

real    0m1.645s
user    0m1.548s
sys     0m0.088s

bioawk is an extension of awk with convenient parsing options. For instance -c fastx makes the sequence available as $seq in fasta and fastq format (including gzipped versions), so the above command should also handle fastq format robustly.
The gsub command is a normal awk function. It returns the number of substituted characters (got it from https://unix.stackexchange.com/a/169575/55127). I welcome comments about how to count inside awk in a less hackish way.
On this fasta example it is almost 3 times slower than the grep, tr, wc pipeline:
$ time grep -v "^>" genome.fa | tr -cd "A" | wc -c
32371810

real    0m0.609s
user    0m0.744s
sys     0m0.184s

(I repeated the timings, and the above values seem representative)
Using python
readfq
I tried the python readfq retrieved from the repository mentioned in this answer: https://bioinformatics.stackexchange.com/a/365/292
$ time python -c "from sys import stdin; from readfq import readfq; print sum((seq.count('A') for _, seq, _ in readfq(stdin)))" < genome.fa
32371810

real    0m0.976s
user    0m0.876s
sys     0m0.100s

"pure python"
Comparing with something adapted from this answer: https://bioinformatics.stackexchange.com/a/363/292
$ time python -c "from sys import stdin; print sum((line.count('A') for line in stdin if not line.startswith('>')))" < genome.fa
32371810

real    0m1.143s
user    0m1.100s
sys     0m0.040s

(It is slower with python 3.6 than with python 2.7 on my computer)
pyGATB
I just learned (08/06/2017) that GATB includes a fasta/fastq parser and has recently released a python API. I tried to use it yesterday to test another answer to the present question and found a bug. This bug is now fixed, so here is a pyGATB-based answer:
$ time python3 -c "from gatb import Bank; print(sum((seq.sequence.count(b\"A\") for seq in Bank(\"genome.fa\"))))"
32371810

real    0m0.663s
user    0m0.568s
sys     0m0.092s

(You can also do sequence.decode("utf-8").count("A") but this seems a little slower.)
Although I used python3.6 here (pyGATB seems python3-only), this is faster than the other two python approaches (for which the reported timings are obtained with python 2.7). This is even almost as fast as the grep, tr, wc pipeline.
Biopython
And, to have even more comparisons, here is a solution using SeqIO.parse from Biopython (with python2.7):
$ time python -c "from Bio import SeqIO; print(sum((rec.seq.count(\"A\") for rec in SeqIO.parse(\"genome.fa\", \"fasta\"))))"
32371810

real    0m1.632s
user    0m1.532s
sys     0m0.096s

This is a bit slower than the "pure python" solution, but perhaps more robust.
There seems to be a slight improvement with @peterjc's suggestion to use the lower level SimpleFastaParser:
$ time python -c "from Bio.SeqIO.FastaIO import SimpleFastaParser; print(sum(seq.count('A') for title, seq in SimpleFastaParser(open('genome.fa'))))"
32371810

real    0m1.618s
user    0m1.500s
sys     0m0.116s

(I did a series of timings and tried to take one that seemed representative, but there's a lot of overlap with the higher-level parser's timings.)
scikit-bio
I read about scikit-bio today (23/06/2017), which has a sequence reader.
$ time python3 -c  "import skbio; print(sum(seq.count('A') for seq in skbio.io.read('genome.fa', format='fasta', verify=False)))"
32371810

real    0m3.643s
user    0m3.440s
sys     0m1.228s

pyfastx
I read about pyfastx today (27/08/2020), which has a sequence reader.
$ time python3 -c "from pyfastx import Fasta; print(sum(seq.count('A') for (_, seq) in Fasta('genome.fa', build_index=False)))"
32371810

real    0m0.781s
user    0m0.740s
sys     0m0.040s

(The timing is done using python 3.6, on the same machine as previous timings from 3 years ago, so it should be comparable.)
fastx_read from mappy
I realize (27/08/2020) that I hadn't added a benchmark for fastx_read from mappy which I regularly use since quite some time already:
$ time python3 -c "from mappy import fastx_read; print(sum(seq.count('A') for (_, seq, _) in fastx_read('genome.fa')))"
32371810

real    0m0.731s
user    0m0.648s
sys     0m0.080s

Benchmarks on a fastq.gz file
I tested some of the above solutions (or adaptations thereof), counting "N" in the same file that was used here: https://bioinformatics.stackexchange.com/a/400/292
bioawk
$ time bioawk -c fastx '{n+=gsub(/N/, "", $seq)} END {print n}' SRR077487_2.filt.fastq.gz
306072

real    1m9.686s
user    1m9.376s
sys     0m0.304s

pigz + readfq python module
readfq doesn't complain and is very fast when I pass directly the compressed fastq, but returns something wrong, so don't forget to manually take care of the decompression.
Here I tried with pigz:
$ time pigz -dc SRR077487_2.filt.fastq.gz | python -c "from sys import stdin; from readfq import readfq; print sum((seq.count('N') for _, seq, _ in readfq(stdin)))"
306072

real    0m52.347s
user    1m40.716s
sys     0m8.604s

The computer has 16 cores, but I suspect the limiting factor for pigz is reading from the disk: the processors are very far from running full speed.
And with gzip:
$ time gzip -dc SRR077487_2.filt.fastq.gz | python -c "from sys import stdin; from readfq import readfq; print sum((seq.count('N') for _, seq, _ in readfq(stdin)))"
306072

real    0m49.448s
user    1m31.984s
sys     0m2.312s

Here gzip and python both used a full processor. Resource usage balance was slightly better. I work on a desktop computer. I suppose a computing server would take advantage of pigz better.
pyGATB
$ time python3 -c "from gatb import Bank; print(sum((seq.sequence.count(b\"N\") for seq in Bank(\"SRR077487_2.filt.fastq.gz\"))))"
306072

real    0m46.784s
user    0m46.404s
sys     0m0.296s

biopython
$ time python -c "from gzip import open as gzopen; from Bio.SeqIO.QualityIO import FastqGeneralIterator; print(sum(seq.count('N') for (_, seq, _) in FastqGeneralIterator(gzopen('SRR077487_2.filt.fastq.gz'))))"
306072

real    3m18.103s
user    3m17.676s
sys     0m0.428s

pyfastx (added 27/08/2020)
$ time python3 -c "from pyfastx import Fastq; print(sum(seq.count('N') for (_, seq, _) in Fastq('SRR077487_2.filt.fastq.gz', build_index=False)))"
306072

real    0m51.140s
user    0m50.784s
sys     0m0.352s

fastx_read from mappy (added 27/08/2020)
$ time python3 -c "from mappy import fastx_read; print(sum(seq.count('N') for (_, seq, _) in fastx_read('SRR077487_2.filt.fastq.gz')))"
306072

real    0m52.336s
user    0m52.024s
sys     0m0.300s

Conclusion
pyGATB and bioawk both handle transparently compression (gzipped or not) and format differences (fasta or fastq). pyGATB is quite a new tool, but seems more efficient compared to the other python modules I tested.
(Update 27/08/2020) pyfastx and mappy seem quite convenient too, and are only slightly slower than pyGATB.

Could you add the Biopython low-level string FASTA parser numbers? time python -c "from Bio.SeqIO.FastaIO import SimpleFastaParser; print(sum(seq.count('A') for title, seq in SimpleFastaParser(open('genome.fa'))))"

Why do all your examples count the letter A, rather than N (and perhaps n) as per the question?

I count letter "A" because the first large file I thought about for my testing happened to have no "N".

I would think that all the Python examples will score badly on this small example file due to the overhead of loading Python and importing libraries. But if the user only wants to count one genome, that is a fair test all the same.

Decompression of gzipped FASTQ is the main issue
If we take real world gzipped FASTQ (which as the OP suggested would be beneficial) rather than trivial FASTA as the starting point then the real issue is actually decompressing the file not counting the Ns and in this case the C program count-N is no longer the fastest solution.
Additionally it would be good to use a specific file for benchmarking which actually has Ns, because you'll get some quite interesting execution time differences with some methods counting the more frequently occurring As rather than Ns.
One such file is:
http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG00096/sequence_read/SRR077487_2.filt.fastq.gz
It's also worth checking the various solutions return the correct answer there should be 306072 Ns in the above file.
Next note that decompression of this file redirected to /dev/null is slower with zcat and gzip (which are both gzip 1.6 on my system) than say a parallel implementation of gzip like Mark Adler's pigz, which appears to use 4 threads for decompression.  All timings represent an average of 10 runs reporting real (wall clock time).
time pigz -dc SRR077487_2.filt.fastq.gz > /dev/null

real    0m29.0132s

time gzip -dc SRR077487_2.filt.fastq.gz > /dev/null

real    0m40.6996s

There is an ~11.7 second difference between the two.  Next if I then try to benchmark a one-liner which performs on FASTQ and gives the correct answer (Note I've yet to encounter FASTQ which is not 4 line, and seriously who generates these files!)
time pigz -dc SRR077487_2.filt.fastq.gz | awk 'NR%4==2{print $1}' | tr -cd N | wc -c
306072

real    0m34.793s

As you can see the counting adds ~5.8 second to the total run time versus pigz based decompression.  Additionally this time delta is higher when using gzip ~6.7 seconds above gzip decompression alone.
time gzip -dc SRR077487_2.filt.fastq.gz | awk 'NR%4==2{print $1}' | tr -cd N | wc -c
306072

real    0m44.399s

The pigz awk tr wc based solution is however ~4.5 seconds faster than the count-N based C code solution:
time count-N SRR077487_2.filt.fastq.gz
2385855128      306072  0

real    0m39.266s

This difference appears to be robust to re-running as many times as I like.  I expect if you could use pthread in the C based solution or alter it to take the standard out from pigz it would also show an increase in performance.
Benchmarking another alternative pigz grep variant appears to take more or less the same time as the tr based variant:
time pigz -dc SRR077487_2.filt.fastq.gz | awk 'NR%4==2{print $1}' | grep -o N | wc -l
306072

real    0m34.869s

Note that the seqtk based solution discussed above is noticeably slower,
time seqtk comp  SRR077487_2.filt.fastq.gz | awk '{x+=$9}END{print x}'
306072

real    1m42.062s

However it's worth noting seqtk comp is doing a bit more than the other solutions.

Honestly, the easiest way (especially for FASTQ) is probably to use a dedicated parsing library, such as R/Bioconductor:
suppressMessages(library(ShortRead))

seq = readFasta(commandArgs(TRUE)[1]) # or readFastq
cat(colSums(alphabetFrequency(sread(seq))[, 'N', drop = FALSE]), '\n')

This may not be the fastest, but it’s pretty fast, since the relevant methods are highly optimised. The biggest overhead is actually probably from loading ShortRead, which is an unjustifiable slog.

If it is raw speed you're after, then writing an own little C/C++ program is probably what you need to do.
Fortunately, the worst part (a fast and reliable parser) has already been tackled: the readfq from Heng Li is probably the fastest FASTA/FASTQ parser around.
And it's easy to use, the example on GitHub can easily be expanded to do what you need. Just add in a parameter parser (for the filename you want to analyse), program a simple 'N'-counter loop and have the results printed to stdout. Done.
Else, if it is simplicity you are after, see Konrad's answer for an R based solution. Or a very simple script using Biopython.

God that’s absolutely atrocious code. ☹ I think that for raw speed, SeqAn (C++) is probably on par, and a lot higher quality.

It's not the prettiest around. Then again, it does what it needs to do and doesn't come with a whole library dependency. Do you have a benchmark with SeqAn & readfq?

I'm curious to see how readfq ans SeqAn compare with the fasta/fastq parser from GATB: http://gatb-core.gforge.inria.fr/doc/api/classgatb_1_1core_1_1bank_1_1impl_1_1BankFasta.html#details
I'm not really able to code in C/C++ though.

For FASTA files, I've implemented a relatively efficient method in pyfaidx v0.4.9.1. This post made me realize that my previous code was quite slow and easy to replace:
$ pip install pyfaidx
$ time faidx -i nucleotide ~/Downloads/hg38.fa
name    start   end A   T   C   G   N   others
chr1    1   248956422   67070277    67244164    48055043    48111528    18475410
chr10   1   133797422   38875926    39027555    27639505    27719976    534460
chr11   1   135086622   39286730    39361954    27903257    27981801    552880
chr11_KI270721v1_random 1   100316  18375   19887   31042   31012   0
...
chrX    1   156040895   46754807    46916701    30523780    30697741    1147866
chrY    1   57227415    7886192 7956168 5285789 5286894 30812372
chrY_KI270740v1_random  1   37240   8274    12978   7232    8756    0

real    2m28.251s
user    2m15.548s
sys 0m9.951s

Here are perl 1-liners for fasta and fastq files that are fast enough if you just want totals:
$ perl -ne 'BEGIN { my $n=0 }; /^>/ or $n += s/N//gi; END { print "N=$n\n" }' contigs.fasta
For each line of a FASTA file, if not header (/^>/) then it's a sequence line.
 Substitute Ns, which returns the number of substitutions made, which is added to the total.
$ zcat reads.fastq | perl -ne 'BEGIN { my $n=0 }; /^[ATCGN]+$/i and $n += s/N//gi; END { print "N=$n\n" }'
For each sequence line of a FASTQ file, matched by /^[ATCGN]+$/i (i.e. all NTs), then count substitutions made in the same way.
You may use pigz instead of zcat (or gunzip -c) if you have pigz installed.
An uncompressed fastq file of 14M 250bp Illumina reads completed in 40sec.

If you could explain what each perl line do, it would improve the answer, specially commenting the differences between both approaches


text:
How exactly is "effective length" used in FPKM calculated?
According to this famous blog post, the effective transcript length is:
$\tilde{l}_i = l_i - \mu$
where $l_i$ is the length of transcript and $\mu$ is the average fragment length. However, typically fragment length is about 300bp. What if when the transcript $l_i$ is smaller than 300? How do you compute the effective length in this case?
A related question: when computing the FPKM of a gene, how to choose a transcript? Do we choose a "canonical" transcript (how?) or combine the signals from all transcripts to a gene-level FPKM?

I have a blog post that describes the effective length (as well as these different relative abundance units).  The short explanation is that what people refer to as the "effective length" is actually the expected effective length (i.e., the expectation, in a statistical sense, of the effective length).  The notion of effective length is actually a property of a transcript, fragment pair, and is equal to the number of potential starting locations for a fragment of this length on the given transcript.  If you take the average, over all fragments mapping to a transcript (potentially weighted by the conditional probability of this mapping), this quantity is the expected effective length of the transcript.  This is often approximated as simply $l_i - \mu$, or $l_i - \mu_{l_i}$ --- where $\mu_{l_i}$ is the mean of the conditional fragment length distribution (conditioned on the fragment length being < $l_i$ to account for exactly the issue that you raise).

The effective length is $\tilde{l}_i = l_i - \mu + 1$ (note the R code at the bottom of Harold's blog post), which in the case of $\mu < l_i$ should be 1. Ideally, you'd use the mean fragment length mapped to the particular feature, rather than a global $\mu$, but that's a lot more work for likely 0 benefit.
Regarding choosing a particular transcript, ideally one would use a method like salmon or kallisto (or RSEM if you have time to kill). Otherwise, your options are (A) choose the major isoform (if it's known in your tissue and condition) or (B) use a "union gene model" (sum the non-redundant exon lengths) or (C) take the median transcript length. None of those three options make much of a difference if you're comparing between samples, though they're all inferior to a salmon/kallisto/etc. metric.
Why are salmon et al. better methods? They don't use arbitrary metrics that will be the same across samples to determine the feature length. Instead, they use expectation maximization (or similarish, since at least salmon doesn't actually use EM) to quantify individual isoform usage. The effective gene length in a sample is then the average of the transcript lengths after weighting for their relative expression (yes, one should remove $\mu$ in there). This can then vary between samples, which is quite useful if you have isoform switching between samples/groups in such a way that methods A-C above would miss (think of cases where the switch is to a smaller transcript with higher coverage over it...resulting in the coverage/length in methods A-C to be tamped down).

But \tilde{l} is a denominator. Setting it to 1 would dramatically increase the value for short transcripts. This sounds dangerous to me... Also, could you clarify what is the advantage of salmon/kallisto over A/B/C? Thanks.

Short transcripts do have absurdly high FPKMs, it's one of the useless things about FPKMs. I'll update the question regarding salmon/kallisto/etc.

For the effective length part please see to Devons answer. I just have a small addition: Kallisto/Salmon/RSEM incorporate all bias estimates into the effective length meaning the effective length not only represent the length bias if you take the values from those tools (given that they were run with the bias algorithms enabled naturally).
With regards to getting gene level estimates you should not choose a specific transcript. Instead you should extract/calculate the RPKM/FPKM/TxPM (transcript per million that Kallisto/Salmon/RSEM outputs) for each transcript and sum them up to get the gene level estimate.

authors: กชพร รัชตวุฒิมงคล


text:
Missing genes and normalisation of RSEM output using EBSeq
Without going into too much background, I just joined up with a lab as a bioinformatics intern while I'm completing my masters degree in the field. The lab has data from an RNA-seq they outsourced, but the only problem is that the only data they have is preprocessed from the company that did the sequencing: filtering the reads, aligning them, and putting the aligned reads through RSEM. I currently have output from RSEM for each of the four samples consisting of: gene id, transcript id(s), length, expected count, and FPKM. I am attempting to get the FASTQ files from the sequencing, but for now, this is what I have, and I'm trying to get something out of it if possible.
I found this article that talks about how expected read counts can be better than raw read counts when analyzing differential expression using EBSeq; it's just one guy's opinion, and it's from 2014, so it may be wrong or outdated, but I thought I'd give it a try since I have the expected counts.
However, I have just a couple of questions about running EBSeq that I can't find the answers to:
1: In the output RSEM files I have, not all genes are represented in each, about 80% of them are, but for the ones that aren't, should I remove them before analysis with EBSeq? It runs when I do, but I'm not sure if it is correct.
2: How do I know which normalization factor to use when running EBSeq? This is more of a conceptual question rather than a technical question.
Thanks!

Yes, that blog post does represent just one guy's opinion (hi!) and it does date all the way back to 2014, which is, like, decades in genomics years. :-) By the way, there is quite a bit of literature discussing the improvements that expected read counts derived from an Expectation Maximization algorithm provide over raw read counts. I'd suggest reading the RSEM papers for a start[1][2].
But your main question is about the mechanics of running RSEM and EBSeq. First, RSEM was written explicitly to be compatible with EBSeq, so I'd be very surprised if it does not work correctly out-of-the-box. Second, EBSeq's MedianNorm function worked very well in my experience for normalizing the library counts. Along those lines, the blog you mentioned above has another post that you may find useful.
But all joking aside, these tools are indeed dated. Alignment-free RNA-Seq tools provide orders-of-magnitude improvements in runtime over the older alignment-based alternatives, with comparable accuracy. Sailfish was the first in a growing list of tools that now includes Salmon and Kallisto. When starting a new analysis from scratch (i.e. if you ever get the original FASTQ files), there's really no good reason not to estimate expression using these much faster tools, followed by a differential expression analysis with DESeq2, edgeR, or sleuth.

1Li B, Ruotti V, Stewart RM, Thomson JA, Dewey CN (2010) RNA-Seq gene expression estimation with read mapping uncertainty. Bioinformatics, 26(4):493–500, doi:10.1093/bioinformatics/btp692.
2Li B, Dewey C (2011) RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics, 12:323, doi:10.1186/1471-2105-12-323.

"There's really no good reason not to estimate expression using these much faster tools" -- unless you don't have the raw reads, as is the case here

Oh wow. Big oversight on my part!

Wow I never expected to get a response from the actual author!  I did use R to process the data frames into a single matrix of all the expected counts for each gene of each sample.

I was just as surprised to see a link to my old blog on StackExchange! :)

As a small follow-up question, I'm ultimately trying to obtain the fold change and associated p-value for each gene per condition.  I've found GetMultiFC() to get the fold changes, but I'm unclear about the exact difference between the fold change and the posterior fold change.  I'm thinking that the posterior fold change is just the fold change for the normalized values, so it's the one I should be using, but I'm not sure about that.  Also, is it possible to get associated p-values?

Include all genes/transcripts in your analysis.

A transcript that is not detected could be undetected through sampling error (i.e. the sequencer / library prep just happened to miss that transcript), or it could be because the transcript isn't generated in a particular sample. It's not uncommon for genes to be switched off in response to different biological factors, so zero-count genes shouldn't be ignored. I can't speak from experience with EBSeq, but as long as the analysis package treats a zero count as "unobserved" rather than "absent" (and makes relevant corrections), it's a good idea to keep them in.


text:
Finding the location and unit length of repetitive sequences within a long read
After discovering a few difficulties with genome assembly, I've taken an interest in finding and categorising repetitive DNA sequences, such as this one from Nippostrongylus brasiliensis [each base is colour-coded as A: green; C: blue; G: yellow; T: red]:

[FASTA file associated with this sequence can be found here]
These sequences with large repeat unit sizes are only detectable (and assembleable) using long reads (e.g. PacBio, nanopore) because any subsequence smaller than the unit length will not be able to distinguish between sequencing error and hitting a different location within the repeat structure. I have been tracking these sequences down in a bulk fashion by two methods:

Running an all-vs-all mapping, and looking for sequences that map to themselves lots of times
Carrying out a compression of the sequence (e.g. bzip2), and finding sequences that have a compression rate that is substantially higher than normal

After I've found a suspicious sequence, I then want to be able to categorise the repeat (e.g. major repeat length, number of tandem repeats, repetitive sequence). This is where I'm getting stuck.
For doing a "look, shiny" demonstration, I currently have a very manual process of getting these sequences into a format that I can visualise. My process is as follows:

Use LAST to produce a dot plot of self-mapping for the mapping
Visually identify the repetitive region, and extract out the region from the sequence
Use a combination of fold -w <width> and less -S to visually inspect the sequence with various potential repeat unit widths to find the most likely repeat unit size
Display the sequence in a rectangular and circular fashion using my own script, wrapping at the repeat unit length

But that process is by no means feasible when I've got thousands of potential repetitive sequences to fish through.
Is there any better way to do this? Given an arbitrary DNA sequence of length >10kb, how can I (in an automated fashion) find both the location of the repeat region, and also the unit length (bearing in mind that there might be multiple repeat structures, with unit lengths from 30bp to 10kb)?
An example sequence can be found here, which has a ~21kb repeat region with ~171bp repeat units about 1/3 of the way into the sequence.
A Kmer-based Analysis
I've now seen human sequences with repetitive regions in excess of 10kb (i.e. out of the range of most linked-read applications). My current idea is centred around creating hash tables of short sequences (currently 13-mers) and tracking their location:

Process the sequence, storing the location of each kmer
For each kmer, find out how many times it appears in the sequence
For repeated kmers, find out how much of a gap there is between the next time that kmer occurs
Report the median and modal gap length of repeated kmers, and statistics associated with their frequency in the sequence

Some local repetitive regions may be lost in the statistics with this approach, it's hard to tell if there are multiple repetitive regions within a single sequence, and if the repeat units are themselves slightly repetitive (enough that a kmer is duplicated within a repeat unit), then the algorithm will under-report the repetitiveness (see step 3).

Would RepeatMasker help you here? Or does it only find known repetitive elements and regions of low complexity? I've never tried to use it to identify simple repeats.

These repeats have high-complexity repeat units; I haven't found the regions being marked by tantan or dust, and suspect that RepeatMasker will be similar. The web service at least isn't able to find the repeat region in the test sequence that I've linked.

I know about this blog post from Gene Myers where he is using strategy of masking repeats in long reads to speed up genome assembly. I am not sure if it is exactly what you want, but I thought it might be useful.

Thank you. TANmask (at the very end of the post) does look very much like the sort of algorithm I'm looking for. It will at least solve the boundary problem, if not the repeat unit length statistic. Can you write that up as an answer?

I've sorted the visualisation out. Here are three alternative representations of repetitive structures for the same sequence:

These were generated using the same R script, callable from the command line:
 $ repaver.r -style dotplot -prefix dotplot MAOA_chrX_43645715-43756016.fa
 $ repaver.r -style profile -prefix profile MAOA_chrX_43645715-43756016.fa
 $ repaver.r -style semicircular -prefix semi MAOA_chrX_43645715-43756016.fa

More details about this are in the presentation I gave at Queenstown Research Week, 2018. I also wrote a chapter in a peer-reviewed eBook here.
This is fast enough that I can run it on the nanopore C. elegans genome in about half an hour, producing these plots for each contig. I don't quite have a method to iterate through this plot and pick out the dominant repeat length at each location, but that's a [relatively] simple extension of what I've already done.
With a lot of optimisations for speed and memory consumption, I've now been able to run it on the full human genome. It takes a couple of days on my desktop computer (64GB RAM + SSD swap) to categorise the 100bp repeat structure of the assembled T2T/CHM13-v1.0 chromosomes.
Code here.

It could be an idea to fragment the long reads into small sequences, like simulating Illumina reads of 150 bp, and then map these small sequences against the original long reads and extract regions with a high coverage?

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


text:
Generic HMM solvers in bioinformatics?
Hidden Markov models (HMMs) are used extensively in bioinformatics, and have been adapted for gene prediction, protein family classification, and a variety of other problems. Indeed, the treatise by Durbin, Eddy and colleagues is one of the defining volumes in this field.
Although the details of each of these different applications of HMMs differ, the core mathematical model remains unchanged, and there are efficient algorithms for computing the probability of the observed sequence given the model, or (perhaps more useful) the most likely hidden sequence given the sequence of observed states.
Accordingly, it seems plausible that there could be a generic software library for solving HMMs. As far as I can tell that's not the case, and most bioinformaticians end up writing HMMs from scratch. Perhaps there's a good reason for this? (Aside from the obvious fact that it's already difficult, nigh impossible, to get funding to build and provide long-term support for open source science software. Academic pressures incentivize building a new tool that you can publish a paper on much more than building on and extending existing tools.)
Do any generic HMM solver libraries exist? If so, would this be tempting enough for bioinformaticians to use rather than writing their own from scratch?

You make it sound like “writing HMMs from scratch” is a huge endeavour. ;-) That said, of course having existing implementations of the relevant algorithms makes sense.

I would also recommend to take a look at pomegranate, a nice Python package for probabilistic graphical models. It includes solvers for HMMs and much more. Under the hood it uses Cythonised code, so it's also quite fast.

There are certainly software libraries for working with HMMs. For a general-purpose implementation in C++, take a look at the SeqAn HMM algorithms.
For your purposes, i.e. “computing … the most likely hidden sequence given the sequence of observed states”, you’d invoke viterbiAlgorithm with your observed sequence and the HMM graph.
More fundamentally I think that most existing, mature implementations are probably found in the domain of signal processing, which has been using them longer than biology, and where most of the underlying theory was developed.

If I remember correctly Ewan Birney's Dynamite (a compiler-compiler) as presented at ISMB 1997 had this functionality, there is also some code here on GitHub https://github.com/birney/wise3 which at least mentions Dynamite.  Suspect Ewan is too busy these days to work on this, although he has tweeted about blowing dust of his old Dynamite, sorry Dynamite code: https://twitter.com/ewanbirney/status/788121636973142016


text:
Are there any databases of templates for common bioinformatic file formats?
I want some templates of different file formats that I can use to test my scripts and identify possible bugs in my code.
For example, consider nucleotide FASTA, a simple but often abused format, I would want templates to capture regular and irregular formats, like I have seen all of these:
1) Single line sequence
>1
ATG

2) Multi-line sequence
>1
AT
G

3) Upper and lower case letters in sequence
>1
Atg

4) Ns and Xs (and possibly other letters) in sequence
>1
ANnxX

5) Unusual headers (sometimes non-ASCI characters, need to consider the encoding)
>ATG >汉字
ATG

6) Whitespace between records
>1
ATG

>2
ATG

7) Duplicated headers
>1
ATG
>1
ATC

8) Empty headers or sequences (valid FASTA?)
>
>

9) No new line '\n' character on last line (can mess-up file concatenation)
>1
A# < no new-line here

10) Different newline characters depending on the OS
>1
A# \r\n vs \n

etc.
There should be separate templates for nucleotide and protein FASTA, and separate ones for aligned FASTA.
It would ideally include other aspects too, like different compression formats (such as .gz, .bzip2) and different file extensions (such as .fa, .fasta).
I have never seen resources that provide templates covering these, but I think it would be useful. Of course I could build my own templates but it would take time to capture all the likely variations of the formats, particularly for more complex file formats.
Note, I am not just interested in FASTA format, it was an example.
Also note, I know about tools (such as BioPython) that should handle many formats well, but they may also have bugs. Anyway, in practice sometimes I end up parsing files myself directly because I don't want the overhead or dependency of an external package.
EDIT: Please don't answer this question to say you don't know of any such resources, me neither, hence the question. bli's helpful answer shows that there is at least one test suite that could be used as a starting point. I know that it is normally easy to look up the specification of any particular file format.

Actually, the FASTA format definition is very simple. It only has two constraints and the second is often ignored: i) header lines must start with > and can contain anything they like except a \n. ii) sequence lines should be 60 characters per line. Nothing else is relevant, not the number of sequence lines, not the characters in the header, not the extension (which is usually irrelevant outside the Windows world anyway), there is no constraint on what characters a sequence can have (hence no difference between protein and nucleotide fasta specifications).

Sure, the specification is simple, but as you say in practice it's not always followed and when you parse a file all these variations becomes relevant. I have had many tools trip up over strange FASTA headers or sequences and clearly it's important for the tools to distinguish nucleotides from amino acids. You may glob files with just 1 particular suffix for example.

If the tools choke on that, it's the tools that are ignoring the standard. The FASTA standard is one of the very few that is actually clear and easy to understand. It is simply a very free format. And yes, of course you can glob on extensions. I, for example, tend to name my protein files .pep and nt .fa but there is no standard for this. All your tools need to do is i) take the entire line after a > as a sequence name and ii) everything else as the sequence. That's all, and that's all it take to be fasta.

Sure, the tools may ignore the standard, and that's why I think it's good to have varied templates to test the tools. And it's not always simple, e.g. if you use FASTA headers to name subsequent files, special characters can be a nuisance. Anyway, as I said FASTA was just an example, the logic extends for other more complex formats

You mention "I am not just interested in FASTA format, it was an example", but bear in mind that It's best to keep questions as specific as possible.

BioJulia maintains BioFmtSpecimens. From the README:

Bioinformatics is rife with formats and parsers for those formats. These parsers don't always agree on the definitions of these formats, since many lack any sort of formal standard.
This repository aims to consolidate a collection of format specimens to create a unified set of data with which to test software against. Testing against the same cases is a first step towards agreeing on the details and edge cases of a format.

You mention Biopython, which contains tests: https://github.com/biopython/biopython/tree/master/Tests.
Some of the tests consist in reading files present in the folders listed in the above link. These files could be a starting point for a database of test files. Whenever one comes across a test case not covered with these files, one could construct a new test file and contribute it to Biopython, along with a test, or at least file an issue: https://github.com/biopython/biopython/issues.
That would be a way to contribute to Biopython while constituting a database of test files.

Thanks, I like this idea; I wonder if BioPerl and other projects also contain similar test suites that could be used

A "me too" here: I've sometimes used the example or test data for BioPython as test data for other programs. You at least have the assurance of it being valid. You might find some similar useful data with the Galaxy tools.

Not that I am aware. It is best to go with format specifications when coding.
Also it may be good to look at the example files that come together with various tools performing file conversions and handling. E.g.

Trimmomatic comes with few .fa files (http://www.usadellab.org/cms/?page=trimmomatic)
Samtools with a toy.sam (https://github.com/lh3/samtools-legacy/blob/master/examples/toy.sam)

As far as I know, there is no single repository that collects all of the common data formats used in bioinformatics.  Typically, you have to go to the source to find the specifications for each format.  There are a few places that collect descriptions of file formats, though:

IGV File Formats, coveres all of the formats usable in the Broad Institute's Integrative Genomics Viewer software (which is a lot).
NCI File Formats, mostly formats used by TGCA (including MAF and VCF).
UCSC Genomics, covers BED, MAF, and a few others.
GenePattern, covers many of the file formats related to microarray data.
GSEA, Broad Gene Set Enrichment Analysis documentation.

No. At least none that I've heard of and I doubt there ever will be. There is no central repository for formats and each tool, community, field etc have their own.
The best you can do is look up the official standard for each format and hope they include examples. Having a truly comprehensive collection of all possible variations is basically impossible for horrible complex formats like VCF. Just consider the simple fact that it allows user-defined INFO fields with arbitrary contents.
So the best you can do is make sure your scripts conform with the standard and hope your input does as well.

There could be a central repository though with templates for the most common formats. I don't see why not? It can't be completely comprehensive of course. To me this seems like a good idea...

@Chris_Rands as a general rule, there are sadly very few centralized anything in biology and that's a general issue with the field. That said, I admit I don't really see the use for this. What you need is the standard. Multiple examples of the standard being applied in different ways doesn't seem all that helpful to me; at least not if the standard is well defined (which is not always the case). Come to think of it, I can't think of any such repository of format examples in general, not only in biology.

Thanks, well I agree about the lack of standardisation/centerlisation. The main use for me would be test cases/unit testing. Knowing the specification is one thing, but spotting all the edge cases can be tricky without testing

@Chris_Rands s/tricky/impossible :) I've now been working in industry after a decade and a half in academia and you wouldn't believe the crap we've been sent as "vcf files". There's no way to cover all the edge cases, some people will give you an excel file with a .vcf extension and expect it to work :(

I think a github repo with examples of edge cases for each format is a great idea - everytime any of us came across a new way some one had done something werid with a format we could issue a PR on repo.

These may not be exactly what your are looking for, but they do contain a wide range of formats with examples. If you want more, you might try searching other major databases hosting other types of into e.g. Uniprot, PDB, NCBI.
https://genome.ucsc.edu/FAQ/FAQformat.html
http://www.ensembl.org/info/website/upload/bed.html

Thanks, I'm familiar with UCSC and Ensembl and their format descriptions but that's not really what I meant by templates


text:
validating identified sub-populations of cells in scRNA-seq
In the analyses of single-cell RNA-seq data there are different unsupervised approaches to identify putative subpopulations (e.g. as available with Suerat or SCDE packages).
Is there a good way of computationally validating the cluster solutions? Different methods may results in slightly different clustering results. How to know which one is the best i.e. representative of biological sub-populations?

There's a bit of folly in trying to computationally validate computational methods without a known truth dataset...

I changed scRNA to single-cell RNA (since scRNA also means small conditional RNA).

@Kamil “scRNA-seq” is an established name for single-cell RNA-seq.

A SC3, single-cell consensus clustering, approach could be helpful here. It aims at achieving "high accuracy and robustness by combining multiple clustering solutions through a consensus approach" https://www.nature.com/nmeth/journal/v14/n5/full/nmeth.4236.html

While better methods of evaluating your clusters would be to use an external dataset or a dataset with known truth, there are a variety of internal validation metrics that can be used to compare clustering solutions without another dataset.
Here are a few metrics:

Davies-Bouldin Index
Calinski-Harabasz Index
Root-Mean-Square Standard Deviation

Many more can be found in this clustering review: http://stke.sciencemag.org/content/9/432/re6
These internal validation metrics grade your clustering solution based on three measures: compactness, connectedness, and separation. When using these metrics to compare clustering solutions, be sure to consider which metric is appropriate for your results as some algorithms work by optimizing certain measures.


text:
What are all the reference files produced by bwa index, and are these dependent upon whether the reference is zipped?
I have indexed a gzipped reference with bwa: bwa index reference.fa.gz, which produces a series of other files reference.fa.gz.{amb,ann,bwt,pac,sa}. These are working fine with bwa alignment.
I have discovered that samtools does not take a gzipped reference, so I am planning to use an unzipped version of the reference for my workflow instead of dealing with two separate representations of the reference.
Are the files reference.fa.gz.{amb,ann,bwt,pac,sa} dependent upon the reference being gzipped? Do I need to reindex the unzipped reference.fa with bwa, or can I just rename the current files to remove the .gz portion of the filename?

You'll get the exact same index (the amb, ann, bwt, pac and sa files) whether the reference is gzipped or not. BWA also makes its own packed reference sequence (the .pac file) so you don't even need the genome around after you index.