text:
What methods exist to calculate RNA expression profile similarity
Some of the work in our lab requires a comparison of a strain across several experimental conditions. We are looking to identify most similar experimental conditions based on the gene transcription response similarity from the cell.
While we could easily invent and create home-grown methods to do it, their implementation and testing are a laborious project in themselves and are outside the scope of our current work.
Are there any methods for RNA expression profile similarity computation that has been already published? If yes, what is your experience using them?

There's a new paper that's just appeared on the subject of "proportionality", including a method by which RNA expression might be compared.
This is a new concept for me, and the article is not easy enough for me to read that I can write up a quick summary; the authors don't seem to devote a section in the paper to defining "proportionality". However, here's an interesting chunk from the article:

We graphed the network of relationships between these mRNAs (S5 Fig.), an approach similar to gene co-expression network [12] or weighted gene co-expression analysis [13] but founded on proportionality and therefore valid for relative data. The network revealed one cluster of 96, and many other smaller clusters of mRNAs behaving proportionally across conditions.
  ...
  We are also keen to raise awareness that correlation (and other statistical methods that assume measurements come from real coordinate space) should not be applied to relative abundances. This is highly relevant to gene coexpression networks [12]. Correlation is at the heart of methods like Weighted Gene Co-expression Network Analysis [13] and heatmap visualization [14]. These methods are potentially misleading if applied to relative data.

proportionality is defined as (1+ beta^2 - 2betaabs(r))
where r is the Pearson correlation of log x and log y, and beta^2 is var(log(y))/var(log(x)).
The authors also note: "Proportionality is appropriate, but ϕ does not satisfy the properties of a distance—most obviously, it is not symmetric unless β = 1". They go on to say "Hence, our approach to forming a dissimilarity matrix is simply to work with ϕ(log xi, log xj) where i < j, in effect, the lower triangle of the matrix of ϕ values between all pairs of components."

But why is it called "proportionality" (in non-formulaic terms)? What is the effect of that formula on expression values? Can you fill in the remainder of this sentence: "Proportionality is a measure that represents how the expression of two ..."

No, not really - I didn't quite understand the use of proportionality. Although I do note that phi is not a measure of proportionality, but a measure of goodness-of-fit to proportionality.

Proportionalty implies the ratio of expression values of two genes will be a constant. Instead of fitting $y=mx+c$ in regression, proportionalty fits $y=mx$

The suggestion to use proportionality is probably the correct one if you are interested in similar patterns between samples. But not if you are interested in absolute differences.
For example: the following two samples are similar in pattern, but not similar on absolute levels:
        Sample 1    Sample 2    Sample 3
Gene A    10          100           80
Gene B     8           80          100
Gene C    12          120          120

Samples 1 and 2 have a perfect proportionality (phi is 0) and also a perfect correlation (as a side note any pair with a perfect proportionality will always have a perfect correlation). However, in terms of logfold changes, samples 2 and 3 are more like each other.
Of course in real life you would never see a comparison like the sample 1 - sample 2 one because normalisation would have removed the scale difference. This was exactly the point brought up by the proportionality paper. But normalisation methods don't normally guarantee that the sum of expression for each sample is identical, and such differences can still occur.
An alternative that might suit more in the second case is either the euclidean distance between the samples, or the euclidean distance on the first two components of a principle component or multi-dimensional scaling. The latter is effectively using mean logFC between samples.
Assuming that x is a matrix containing normalised, log transformed expression values, you could use R and limma to calculate distance in multi-dimensional scaled space as follows:
library(limma)
mds <- plotMDS(x, plot=FALSE)
mds <- data.frame(mds$x, mds$y)
distances <- dist(mds)

Normalisation would only remove the scale difference if it were consistent across a large proportion of genes. In the example you have given, this would indeed be the case if only the indicated genes were assayed, but with this low number of genes it would be more common to normalise based on a housekeeping gene rather than the total gene expression.

Yes, this was my point - hence why there might be reason to use some distance measure (such as euclidean distance) that measured scale difference rather than pattern differences.

If you have a lot of samples with very different environmental conditions, then a Weighted Gene Correlation Network Analysis (WGCNA) might be appropriate.
This type of analysis looks for genes that trace similar (or opposing) expression patterns throughout different conditions (e.g. high-medium-medium-low-absent-high would be highly negatively correlated with low-medium-medium-high-very high-low).
That particular paper introduces the concept of "modules", which are groups of genes that share similar expression patterns. Functions are available for plotting how the expression of canonical module members changes throughout the different conditions, and for identifying which module (or modules) a particular gene is likely to be a member of.
WGCNA works best when there's a lot of different expression changes in the different conditions, which sounds like it would fit well with your project. However, it concentrates more on the genes, rather than the conditions (which seems like it would be less useful for you).

The euclidean distance is probably the simplest, both conceptually and in terms of implementation. It's by no means an elegant solution, and it may not perform well in certain circumstances.
Euclidean distance is easiest to conceptualize as the distance between two points in a two dimensional space.
Y
^
|
|
|
|            * p = (3, 3)
|
|
|
|        * q = (2, 1)
|
----------------------------> X

In this example, the distance between the two points is
d(p, q) = sqrt( (p_x-q_x)^2 + (p_y-q_y)^2 )
        = sqrt( (2-3)^2 + (1-3)^2 )
        = sqrt(5)
        ≈ 2.24

For a gene expression profile with two genes, this is exactly how the euclidean distance would be calculated, using expression values from one gene as the X axis and expression values from the other gene as the Y axis. Realistically, though, gene expression profiles typically contain thousands or tens of thousands of genes, so instead we use the generalization of the distance calculation for N dimensions.
d(p, q) = sqrt( (p_1-q_1)^2 + (p_2-q_2)^2 + ... + (p_N-q_N)^2 )

Packages for R and Python make these types of calculations trivial once you have the data loaded into the correct data structure. See Ian's answer for some example R code.


text:
What causes the difference in total length of assembled contigs and scaffolds in SOAPdenovo2?
I use SOAPdenovo2 to assemble a large genome (4.8G) using ~20X paired-end reads. The  total length of contig sizes is 6.3G while total length of scaffolds is 2.7G. Note that this is a haploid genome, so there is no issue of heterozygosity for scaffolding.
I am wondering what happened during the scaffolding procedure? How comes that sum of scaffolds is so much smaller than sum of contigs?

Scaffolding probably throws out small contigs. Without being an expert in soapdenovo2, I would first bwasw the contigs to he scaffolds to see which sequence looks "lost".

Something else that could be happening is that contigs that are being collapsed into "heterozygous" groups. This would be a particular problem when a genome has a substantial amount of repeated sequence. Digging deep into the supplementary information of the SOAPdenovo2 paper, I've found the following information:

In SOAPdenovo2, heterozygous contig pairs are recognized by utilizing the information of contig depth and the locality of contig. The recognized heterozygous contig pairs should obey the following rules: 1) the similarity between contigs should be high enough, for example, ≥ 95%; 2) the depth of both contigs should be near half of the average depth or all contigs, complying Poisson distribution; 3) the two contigs should be located adjacently in a scaffold and have no relationship to each other inferred by paired-end reads information. The normal contigs neighboring the heterozygous regions, if they exist, could be connected to both of the heterozygous contig pairs (H1 and H2). Only the contig with relatively higher depth in a heterozygous contig pair were kept for scaffolding. The method reduces the influence of genome heterozygosity on final scaffold length. All heterozygous contig pairs were outputted to a file to facilitate further analysis. However, the trade-off of this method is that it might incorrectly remove paralogous contigs. This problem could be relieved by a gap-filling procedure while the removed copy of paralogous contigs would be represented by gaps during scaffolding.

Nice. However, this can't be a complete answer, because of "2) the depth of both contigs should be near half of the average depth or all contigs, complying Poisson distribution;". Because in the OP's case the reduction of size is from 6.3G to 2.7G, which means that more than a half of contigs would be coverage outliers.

That would be true for a diploid genome that has a low amount of repeats, but if there's substantial gene duplication (as happens with very large genomes), then the collapsed sequence will be a lot smaller. However, I'm not excluding other factors, just pointing out that this is one of the possible explanations.

I just wanted to point out, that there got to be something else as well given those numbers...

Maybe also the filtering of small contigs as @manuel says.

Nice. I see there are at least a few factors we agreed: low coverage of contigs, small contigs, also possible high gene duplication levels. I also have 3 mate-pair libraries with 10X each, will use them for scaffolding and gap closer. Hopefully this will tell us a bit more on this issue. i will keep you update.

source: Stack Exchange

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

sort: votes


text:
How can I speed up INDEL calling/correction on BAM files?
The samtools mpileup command has quite a neat feature that it is able to correct mapping errors associated with misalignment of INDELs. By default, the mpileup command will not work for reads that have more than 250X coverage of the reference genome. While this limit can be increased, very high coverage causes the mpileup program to grind to a halt, so it'd be nice to know if there's some easy way to make that faster.
To add a bit more context, I've been doing this with mitochondrial genome reads that were extracted from both Illumina whole-genome sequencing (coverage ~1000X), and from targeted amplicon sequencing done on the IonTorrent (coverage up to ~4000X).
I see that @rightskewed has mentioned the downsampling ability of samtools with samtools view -s <float> (see here), which seems like it might work as a solution for this if used prior to the mpileup operation.

I wasn't aware of the samtools subsampling when I had this problem a couple of years ago, so ended up writing my own digital normalisation method to deal with mapped reads. This method reduces the genome coverage, but preserves reads where coverage is low.
Because I was working with IonTorrent reads (which have variable length), I came up with the idea of selecting the longest read that mapped to each location in the genome (assuming such a read existed). This meant that the highly variable coverage for different samples (sometimes as low as 200X, sometimes as high as 4000X) was flattened out to a much more consistent coverage of about 100-200X. Here's the core of the Perl code that I wrote:
  if(($F[2] ne $seqName) || ($F[3] != $pos) || (length($bestSeq) <= length($F[9]))){
    if(length($bestSeq) == length($F[9])){
      ## reservoir sampling with a reservoir size of 1
      ## See https://en.wikipedia.org/wiki/Reservoir_sampling
      ## * with probability 1/i, keep the new item instead of the current item
      $seenCount++;
      if(!rand($seenCount)){
        ## i.e. if rand($seenCount) == 0, then continue with replacement
        next;
      }
    } else {
      $seenCount = 1;
    }
    if(($F[2] ne $seqName) || ($F[3] != $pos)){
      if($output eq "fastq"){
        printSeq($bestID, $bestSeq, $bestQual);
      } elsif($output eq "sam"){
        print($bestLine);
      }
    }
    $seqName = $F[2];
    $pos = $F[3];
    $bestLine = $line;
    $bestID = $F[0];
    $bestFlags = $F[1];
    $bestSeq = $F[9];
    $bestQual = $F[10];

The full code (which works as a filter on uncompressed SAM files) can be found here.

To clarify: this is Perl code?

yes, Perl code. I can tell by the dollar symbols everywhere, but I guess it's not obvious to people who are not familiar with Perl.

Well I do speak Perl but it’s hardly the only language with dollars before identifiers. ;-)

Indeed. My mind was blank, but my wife reminded me of bash/sh.

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


text:
How do I efficiently subset a very large line-based file?
This has come up repeatedly recently: I have a very large text file (in the order of several GiB) and I need to perform line-based subsetting for around 10,000 lines. There exist solutions for specific scenarios (e.g. samtools view -s for randomly sampling BAM files) but sometimes my use-case doesn’t fit into these categories.
Unfortunately a naïve sed-based solution is extremely slow:
time sed -n -f <(awk -vOFS='' '{print $0, "p"}' line_numbers.txt) input_file > selected_lines.txt

Where line_numbers.txt is a file containing one line number per line.
Forget running this for 10,000 lines; it’s already grinding to a halt for a mere 1000.
How can I speed this up, ideally so that it scales only with the size of the input file, and has more or less constant runtime n the number of lines that I subset?

I think this question belongs more to stackoverflow than here, because it wants the general case to be addressed, not some specific bioinformatics file format. I tried to vote to close, but somehow, I don't have the choice to select stackoverflow as a more appropriate site.

@bli I disagree, I think we should be open to any question that can arise in the course of a bioinformatician's work. Many, many bioinformatics questions can be boiled down to simple text parsing operations (consider translating between different sequence formats, for example) but I still feel those should be on topic. Also, you can't migrate to another site unless a specific migration path has been set up and we don't have any of those yet.

@bli I disagree as well, because Stack Overflow questions are pretty much always tied to a language (or the opposite, tied to none, and no code solution is required). By contrast, I'm interested in a solution but I don't care about the technology.

It would be helpful to include a bit more bioinformatics context / story into this question. You've given an example of what wouldn't be a use case, but not what would. I can't think off the top of my head about a situation where I've needed to subset based on the line number. It's more common for me to be subsetting based on a value in one of the columns.

With awk only: awk 'BEGIN{while((getline<"line_num.txt")>0)l[$1]=1}NR in l' input_file.

Related post on SO, using awk, Subset a file by row and column numbers. Agreed this comes up frequently, and I think this post is on-topic for bioinformatics.

It would be nice to see benchmarking results for top answers.

If you don't mind reading the whole files into memory, this should be quick in Python operator.itemgetter(*[int(line) for line in open(line_file)])(open(input_file).readlines())

@Chris_Rands Why not use itertools.islice on a generator? This will avoid loading the whole file at once, and have minimal overhead.

@KonradRudolph you mean like islice each continuous range of line numbers separately, using the relative line number as the generator is consumed?

@Chris_Rands Damn, I’ve been using too much R. I thought you could pass non-contiguous index ranges to islice.

Turns out, simply keeping track of the next candidate line (after sorting the sample line numbers) fixes the performance issue, and most of the remaining slowness seems to be due to the overhead of actually reading the file so there’s not very much to improve.
Since I don’t know how how to do this in sed, and it’s not trivial in awk either, here’s a Perl script:
#!/usr/bin/env perl

use strict;
use warnings;

my $file = $ARGV[0];
my $lines_file = $ARGV[1];

open my $lines_fh, '<', $lines_file or die "Cannot read file $lines_file";
chomp (my @lines = <$lines_fh>);
close $lines_fh;

@lines = sort {$a <=> $b} @lines;

open my $fh, '<', $file or die "Cannot read file $file";
my $line = 1;
my $next_line = 0;
while (<$fh>) {
    last if $next_line == scalar @lines;
    if ($line++ == $lines[$next_line]) {
        $next_line++;
        print;
    }
}
close $fh;

I’ve implemented a similar function in C++ for an R package, that's only slightly longer than the Perl script. It is ~3 times faster than the Perl script on my test file.

Wouldn't this be far faster if you used a hash instead? $line{$_}++ for @lines and then, in the while loop: print if $line{$.}.

@terdon No. Hash lookup time is constant on average, but slower than a direct index lookup in a contiguous array (plus the occasional index increment). So much for theory. I also tested it, and it holds in practice. (And, just as expected, the actual difference is almost negligible.)

Well, TIL. Thanks, I had always assumed hashes were faster by definition (that's what you get when you let biologists write code :).

Perl should be fairly fast with this when using a hash set to store the list of lines. A structure like this also works for subsetting based on a field value, where the comparison would be with the field rather than "$.":
#!/usr/bin/perl

use strict;
use warnings;

my $lines_file = $ARGV[0];
my %include_lines = ();

open my $lines_fh, '<', $lines_file or die "Cannot read file $lines_file";
while(<$lines_fh>){
  chomp;
  $include_lines{$_} = 1;
}
close $lines_fh;

while(<>){
  if($include_lines{$.}){ # "$." -- line number of current file
    print;
  }
}

Note that according to this SO answer, the "$." operator is not strictly the current line number, and can be influenced by different file operations or other settings.
Edit: just saw your comment about speed in your answer, comparing hash sets to a sorted list. The $lines[$next_line] bit feels a bit odd to me. Have you tried out using shift or pop on a sorted list to fetch the next line:
#!/usr/bin/perl

use strict;
use warnings;

my $lines_file = $ARGV[0];

open my $lines_fh, '<', $lines_file or die "Cannot read file $lines_file";
chomp (my @lines = <$lines_fh>);
close $lines_fh;

@lines = sort {$a <=> $b} @lines;
my $next_line = shift(@lines);

while (<>) {
    if ($. == $next_line) {
        $next_line = shift(@lines);
        print;
        last if (!@lines);
    }
}

I changed the shift to a pop (reversing the sort order to match), and got times of 4.8s, 5.8s and 4.1s for Konrad's original code, my hash code, and my pop code respectively, fetching 10,000 lines from /usr/share/dict/british-english-insane 25 times over (after copying input files to /tmp). They're all the same from my perspective: quick enough that it would take me longer to type out the command than it would to run it. Using shift instead of pop doesn't seem to change the time noticeably.

shift/pop modify the array which may be quite slow if memory is moved as a consequence (but I admit I don’t know whether Perl does that). It should definitely never be faster than indexed access; though to be honest it’s less code.

Sorry about that; I wrote these up quickly without testing. I've fixed those bugs (changed a % to a $ in line 12), and put a "my" at the front of line 18.

Some related questions appear in other sites, with potentially interesting solutions, which I report here:
To sample approximately 1% of the non-empty lines:
awk 'BEGIN {srand()} !/^$/ { if (rand() <= .01) print $0}' input_file

(from https://stackoverflow.com/a/692321/1878788)
To select 1000 random lines:
shuf -n 1000 input_file

(from https://stackoverflow.com/a/15065490/1878788, and https://unix.stackexchange.com/a/108604/55127)
Edit: Python solutions using a list of lines
Using a set of line indices and selecting lines by testing set membership:
#!/usr/bin/env python3

import sys

with open(sys.argv[2], "r") as line_numbers_file:
    line_indices = set(int(line) - 1 for line in line_numbers_file)

with open(sys.argv[1], "r") as input_file:
    print(*(line.strip() for (idx, line) in enumerate(input_file)
            if idx in line_indices), sep="\n")

Using a numpy boolean array together with itertools.compress:
#!/usr/bin/env python3

import sys
from itertools import compress
from numpy import zeros

with open(sys.argv[2], "r") as line_numbers_file:
    line_indices = [int(line) - 1 for line in line_numbers_file]

selector = zeros(max(line_indices) + 1, dtype=bool)
selector[line_indices] = 1

with open(sys.argv[1], "r") as input_file:
    print(*(line.strip() for line in compress(input_file, selector)), sep="\n")

I did some tests on a file containing 15774756 sam records and a list of 10000 pre-generated random line numbers.
The perl script proposed by Konrad Rudolph (https://bioinformatics.stackexchange.com/a/454/292) runs in about 5.3 seconds.
The set membership testing python solution runs in about 4.45 seconds.
The compress based solution runs in about 3.4 seconds.
I suspect this may vary a lot depending on the highest line number we want, since the number of iterations will depend on the length of the boolean array. Here the highest line number was 15773768, so pretty high compared with the total number of lines.
I tried with python 3.6. I suspect that python 2.7 could be slightly faster, but haven't tested.

I think the question is about selecting a specific set of lines from the file, not just a random sample.

@terdon is right, the user cases I've come across recently required specific, nonrandom lines. Nevertheless, this is a good complement.

@KonradRudolph Since you mentioned samtools view -s, I thought that you were thinking of random selection for the rest of the question and that your list of line numbers was randomly generated.

@bli My intention was to show this as a counter-example how I can't solve the problem.

int(line.strip()) can be simplified to int(line); int ignores leading/trailing whitespace anyway

I wrote a command-line (C++14) tool called subset which is up on Github: https://github.com/alexpreynolds/subset
This should be reasonably memory efficient and fast. The subset tool does not store input lines in a table, but instead streams through the file once, storing a 4 or 8k buffer chunk of the input file (depending on OS).
It stores line numbers in an array, but eight bytes per integer * 100k is 800kB, for that use case — not very much memory there.
There's an O(nlogn) sort penalty on the line number array, but again this list will be much smaller than the query file, and integer sorting is fairly optimized, so the hit should be small.
If your line number list is already sorted, I could add an option to skip sorting; let me know if that would be useful.
The filtering step walks through the line number array and input file linearly, printing lines where there are index matches, and skipping over the rest.
Indeed, subset will quit early in parsing the input file if there are no more line numbers to query. So this feature is especially useful for speeding up filtering of very large query files. (If your query file has 1M rows, say, and your last line number of interest is 12345, there's no reason to read through the rest of the file.)
You can grab, build and install it like so:
$ git clone https://github.com/alexpreynolds/subset.git
$ cd subset
$ make
$ cp subset /usr/local/bin

Once the binary is in your path, there are a couple ways to use it.
For example, you can specify a start index and length value. The following grabs seven lines starting with the 33rd line (32 as a 0-indexed value):
$ subset --prefix-with-indices -s 32 -n 7 -i query.txt > answer.txt

Or you can specify a text file containing line numbers, each on a separate line. The following reads in a file called line-numbers.txt and uses that to filter query.txt:
$ subset --prefix-with-indices -l line-numbers.txt -i query.txt > answer.txt

The indices in line-numbers.txt should be positive, 0-indexed integers. The list of numbers does not need to be sorted, as subset will sort the list of numbers for you. This is so that an efficient single pass through the input/query file can be done.
You can leave out --prefix-with-indices to leave out the debug prefix. This is there so that you can do a sanity check on the result.
The test/makefile tests demonstrate options and usages for the two types of filtering.


text:
R packages for data analyses of pooled CRISPR screens
We are designing a CRISPR/Cas9 experiment and thinking of the down-stream data analyses.
Are there any R packages for analysing raw NGS read count data from pooled genetic screens using CRIPSR/Cas9 to disrupt gene expression in a population of cells? I guess we will need to start with basics, i.e. sequences processing, data exploration, visualisation etc.

How is this different from other expression analysis? I.e. why can’t you use one of the standard workflows for differential expression analysis?

I would just like to highlight that using CRISPR/Cas9 for editing and just estimating the expression might not be a great idea per se. If you have really a lot number of samples you might be interested first to check the mutational burden introduced by this editing. Then preferably go with expression analysis.  I would go with the below tools as mentioned as well but will not do anything without even considering if my editing is throwing a lot of SNVs or not. If its not throwing then you should be fine doing the downstream analysis with linear mixed effect.

Assuming you mean CRISPR screens targeting many loci (for example using the GeCKO library) there is an R package here.
It uses a linear mixed effect model to compare guide counts before and after a selection step, allowing for multiple guides/gene and multiple replicates.
It can also do the initial read alignment AFAIK and the author should be quite responsive in case you have questions.
Another option could be the python program MAGeCK

From my experience, there is few about CRISPR pooled screenings in R. One of them is ScreenBEAM. It uses a linear mixture model to determine significance at gene level by considering all the guides together at the same step in the analysis. However, it does not work very well for me.
There are other ways of analyze these screenings outside of R that are easy to use such the suggested MAGeCK or others such HiTSelect, BAGEL or CRISPR-Analyzer.
CRISPR-Analyzer is an online tool that could be useful and that recapitulates some of the mentioned methods in a single analysis. However, if you want to work from command line I would recommend you BAGEL for gene essentiality experiments and MAGeCK for other experimental setup (i.e, control vs treatment)
Sorry for not putting all the links but I can only put two at the moment!

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


text:
confidence ellipses for MDS plot in edgeR?
Is it possible to draw e.g. 95% confidence ellipses around samples from the same group on the results from the plotMDS function under edgeR? If so, how?

What does the 95% confidence correspond to? MDS doesn’t directly give confidence intervals.

That depends on whether or not you have a frequentist or bayesian view about confidence intervals. I would expect (from a bayesian interpretation) that the confidence ellipses on a 2D plot would contain 95% of the points that would be expected in a given sample group.

@gringer Hmm. For a typical experimental setup in NGS, a sample group contains < 10 points. Often < 5. So there’s a lot of uncertainty about the placement of the ellipses based on this. ;-)

I agree to what @KonradRudolph said. If your group has less < 10 or  < 5 , in that case the ellipse criterion might not even hold true or I would not believe in it. First of all MDSplot might not be needed as well for doing the 95% ellipse as well. This link explains pretty well about its use. I am just thinking some times we do a lot of over-estimation or over-calculation of certain statistical measures which for the first choice might not be required. Do you even want to do that? https://stats.stackexchange.com/questions/217374/real-meaning-of-confidence-ellipse .

The question assumes that the position on an MDS plot is a population parameter for a population from which each of the samples in the group are drawn from. It assumes that each sample is a point estimate of the population for that group. Under this interpretation, a confidence ellipse is a valid thing. However, i'm not sure position on an MDS plot can be called a population parameter due to the nature of how an MDS plot is constructed (specifcally the position of each sample is not independent of the position of the others). But perhaps some suitable parameter does exist?

I can see that there would be a three step process to doing this:

Merge counts from all samples in the group and then resample pseudo-replicates from this. If x is a matrix with samples in the group being columns and genes being rows
s <- sample(row.names(x), n = mean(colSums(x)), probs=rowSums(x)/sum(rowSums(x))
stab <- table(s)
s <- as.vector(stab)
names(s) <- names(stab)

Do this thousands of times.
You would then want to project those onto your MDS plot space - i'm less clear about how you could do this without perturbing the space itself.
Calculate the parameter for an ellipise that would contain 95% of these points. Again, I'm not so sure about how to do this, but I think the car pacakge might be a good place to start looking.


text:
When to account for the blacklisted genomic regions in ChIP-seq data analyses?
We have heard in the group that it is important to keep track of and to filter artifact regions when analysing data from functional genomics experiments, especially ChIP-seq.
Here, we have seen pipelines that remove the ENCODE tracks i) before cross-correlations QCs, ii) after cross-correlations QC but before peak calling and iii) after peak calling.
We have noticed that removal of the tracks does not affect significantly cross-correlation and peak-independent QCs. However, we are not sure whether peak calling should be done on the filtered tracks or not?

Aside: Cross-correlation is largely meaningless, regardless of what some of the ENCODE folks might argue. When we process our DEEP samples we don't even look at that value.
Regardless, if you're using SPP/phantomPeakQual for cross-correlation then note that it already removes the highest peaks from your dataset before computing the cross-correlation (in fact, it can remove most of the actual peaks too, which makes one further wonder what it's actually telling you). I don't know that this is actually documented anywhere, it's something I noticed when going through the code while pondering whether to implement it in deepTools. But at least it's ignoring these regions already :)
In general, it's most convenient to just remove peaks overlapping blacklisted regions. In an ideal world you'd filter out the blacklisted reads before peak calling, but (1) this is really inconvenient (more time and disk required) and (2) I've never seen an appreciable gain in peak calling performance. In theory at least you should be losing sensitivity right around blacklisted regions if you don't remove reads overlapping blacklisted regions, but you have to ask yourself whether you want to trust such peaks anyway. For other QC steps, at least with deepTools we provide a parameter with every tool to specify a BED file of blacklisted regions to skip.
As an aside, there are many fewer blacklisted regions in more recent genome builds (GRCh38 and GRCm38 at least), so this is less of an issue in general with them.

I filtered my bams out of curiosity. I had around 500 peaks before filtering. When i called on the filtered bams, i lost 100 (from the blacklist) and gained 38. The peak call was using macs2. Looking at the gained peaks, they are legit and a binding partner of my TF had peaks there.
Note however i ran out of ram for bam > 8GB using bedtools intersect to filter. And if you do it, sort the bam and make sure the chromosome names are ordered using the same tool between blacklist bed and your bams.

It is much more efficient to use samtools view. It accepts a BED file and only keeps reads overlapping this. So you take a BED file of the entire genome, make the complement of this with your blacklist and what is left is the regions to keep. Provide that to samtools view, this is fast with almost no memory footprint.


text:
How to compare groups using WGS data?
We have whole genome sequencing data for patients (not-cancer) (n=60) and for healthy controls (n=20). The sequencing centre has provided us with the best practice bioinformatics analyses including reads mapping (.BAM) and variant calling using GATK (.vcf) as well as annotation (annotated .vcf and .gVCF).
What should be our next steps? We are interested to see if there are any differences (global and/or specific) between the groups.

Is that really what you're interested in, or would you rather find likely causative differences?

Hi Jessica, thanks for writing a question up for Bioinformatics Stack Exchange. Your question fits well with the format we like on this site, in that it is specific to a particular task, and includes a bit of a story behind it.
If you would like to improve the quality of answers for this question, my suggestion would be to add in information about what types of analyses you are already familiar or comfortable with. Have you used GATK or R before? Is this a diagnostic setting, where a structured, documented, and replicable analysis is important?

Sounds like you would like to do GWAS, but for that you would probably need way more samples. Can you subset your variants to test only relevant regions of genome? That would help you to get a statistical power (less tests -> better).

I'm not familiar with the program, but apparently Hail is setting itself up as a swiss-army chainsaw project for doing downstream analysis on variant-called datasets.
An overview of Hail can be found here:
http://blog.cloudera.com/blog/2017/05/hail-scalable-genomics-analysis-with-spark/
A tutorial on association testing can be found here:
https://hail.is/hail/tutorial.html#Association-testing

If you have gVCFs, the first thing you should try is joint variant calling. According to GATK, joint variant calling "empowers variant discovery by providing the ability to leverage population-wide information from a cohort of multiple sample[sic], allowing us to detect variants with great sensitivity and genotype samples as accurately as possible." source.
Once you have two sets of high-quality variants, what you do next depends on your research question. Are you looking for druggable mutations? Underlying causes? Biomarkers? Patient prognosis? You can look at the most frequently mutated positions in the patient cohort and compare them to your reference cohort, check for co-occurring mutations, cluster them, do principal component analysis, do some machine learning to stratify them, etc.


text:
How do I carry out an ancestry/admixture test on a single VCF file?
This is a question from /u/beneficii9 on reddit. The original post can be found here.

Through the Personal Genome Project, I have had my whole genome sequenced by Veritas, and have it in the form of a single VCF file for the whole genome and one BAS file for each chromosome. The reference genome associated with the VCF file is hg19. It has been helpful in health data; for example, I discovered I'm homozygous for the non-functional variant CYP-2D6 gene (rs3892097), which can render several common medications useless, and helps explain why some medicines didn't really work for me. My doctor has found this information very helpful.
Unfortunately, I can't find any way of looking at admixture or ancestry. I've tried setting everything up using a combination of VCFTools, Plink1.9, and ADMIXTURE, but I can't get it to work. I think for ADMIXTURE you have to have a bunch of genomes sorted by geographical origin to compare your genome against, but I'm not sure how to do that, and what's online isn't very clear to me. So scratch that one off.
I've tried converting the file to 23andme format (and at this /u/psychosomaticism has been very helpful). I did that (though it seems there were problems because of the way the VCF file was set up). But the websites that take the data want you to point them to your 23andme account, and that doesn't really work if you only have the file. 23andme doesn't provide for people who had their whole genomes sequenced. They want you to give them a saliva sample like everyone else.
So, what can I do?

A modified implementation of Vivek's answer.
peddy is a Python package that samples an input .vcf at ~25000 sites and projects onto a principal component space built on 2504 thousand genome samples. The author has extensive documentation of the tool's features and a link to the preprint.
I downloaded the .vcf and .vcf.tbi for the NA12878 sample from Genome in a Bottle's ftp here. Then, created a custom .ped file, NA12878.ped, with the contents:
NA12878 HG001   0       0       2       0
At the command line:
$ peddy --plot --prefix myvcf HG001_GRCh37_GIAB_highconf_CG-IllFB-IllGATKHC-Ion-10X-SOLID_CHROM1-X_v.3.3.2_highconf_PGandRTGphasetransfer.vcf.gz NA12878.ped

The output files all have the prefix myvcf., here's myvcf.pca_check.png

I think it might be worth adding: peddy requires coloredlogs and cyvcf2 the later of which will required you to have gcc&g++ and the headers (e.g. -dev versions in ubuntu) of libbz2, liblzma and curl. This knowledge would have saved a precious 15 precious moments of my life :)

Also, now that it's installed I seem to get the error: AssertionError: error loading tabix index for b'HG001_GRCh37_GIAB_highconf_CG-IllFB-IllGATKHC-Ion-10X-SOLID_CHROM1-X_v.3.3.2_highconf_PGandRTGphasetransfer.vcf.gz' (Also based on my understanding shouldn't the 2 be 0 or left to the user to set based on his gender ?)

A suggestion from the UK Biobank QC methods white paper:

Create a set of SNPs common to your VCF and the 1000 genomes phase 3 call set.
Perform PCA of the 1000 genomes samples using eigenstrat smartpca. You might have to convert to binary plink format.
Project your genotypes onto this pre-computed PCA space and visualize using ggplot to see which cluster you fall in.

Ancestry testing is a tricky subject. I spent a good chunk of my PhD project on questions similar to this, and haven't really found a good answer for how to detect unmodelled ancestry.
The general idea of how I tried to approach ancestry determination was to create a model set containing well-defined groups of individuals with a known, specific ancestral background. A query individual (or individuals) were then added to this group, and an ancestry-estimating program was run (specifically Structure) to work out what proportion of the unknown individual was attributable to each of the known groups. I would expect that most genetic ancestry tests follow a similar approach to this, although perhaps with a little bit less care about the definition of the model population groups.
Problems arise when individuals are incorrectly assigned to a particular group, when groups are present in the model set that represent a larger proportion of individuals than most other groups, when groups are closely-related to other groups, and when an ancestral history is present in a test individual that doesn't match any of the model groups. And all that assumes that the marker set used for ancestry determination is perfect: no bias towards any particular group, and no systematic genotyping error.
This doesn't mean that ancestry testing won't work, but it's a good idea to take the results with a large grain of salt. There's a good example of a media personality in New Zealand who was told that she had a very high probability of being 100% Māori, despite having a good knowledge of her own family history that indicated one European ancestor a few generations back on both sides of her family.

Adding another answer that, while much more involved, offers flexibility, walks through minimal steps, and provides visualizations.
I created a Dash app and repository that explores this type of analysis using published ancestry informative SNPs and 1000 Genomes Project Data.

Identify published ancestry informative SNPs (AISNPs).

Kidd et al.
Kosoy et al.

Download 1000 Genomes Project genotype bcf data.
Limit genotype data from (2) to AISNPs loci from (1).
One-hot encode the genotypes.
Perform dimensionality reduction (PCA, t-SNE, or UMAP)
Plot components.

The original post seemed to want to project their genome against a reference sample. The plot_walkthrough.ipynb in the tgviz repo shows how this can be achieved.
Output from the notebook on my own data:

Related: ethnicity check either from bam or vcf files
The missing piece of your puzzle is an ancestry-labeled reference dataset.  1000 Genomes phase 3 usually works well enough here, and it now has variant callsets for both GRCh37 and GRCh38.


text:
Accidental mapping of eukaryotic reads in a metagenomic dataset
This is a question from /u/wipeyourmit on reddit. The original post can be found here.

If I have a metagenomic dataset that contains reads from both
  eukaryotes and prokaryotes and then I annotate by running DIAMOND or
  HMMER against a bacterial database how much of a risk do I run of
  eukaryotic reads being annotated in the process?
I was hoping to use the eggNOG mapper to search against the bacterial
  and archaeal databases and to exclude the eukaryotic portion of my
  dataset. Is the eukaryotic filtering something that I would need to do
  in a step prior to this?

Can you be more specific; by metagenomic data, do you mean shotgun metagenomic data or 16S? I assume without depletion of host DNA or an enrichment protocol? What is the sample site and species? For example, the ratio of prokaryotic:eukaryotic reads is much higher in the human gut than most other human body sites

If it includes eukaryotes, it's probably not 16s. I guess I would assume shotgun metagenomic data.

I guess this is a disadvantage of cross-posting, one can't really clarify the original question. Note cross-posting is being discussed now on meta https://bioinformatics.meta.stackexchange.com/questions/78/please-stop-cross-posting-behavior

Please don't do this. If you want to repeat a question, it must be one for which you can provide clarification. Otherwise, the exercise is futile.

If you know the eukaryotic contaminant present you could use bbsplit.sh from BBMap suite to split/bin reads first using a reference for that contaminant (into one or as many bins as reference sequences you provide).

The contaminants are probably numerous and/or unknown in a metagenome?

My comment was with reference to the eukaryote being the contaminant. The original question refers to a vague "both eukaryotes and prokaryotes" specification. Hopefully there is a small number of known "contaminants". In absence of concrete information none of the tools/answers are going to be a perfect solution.

Sure, I was just saying that typically in metagenomics, you will capture many eukaryotes (that you may deem contaminants) in your sample. Often these environments are lacking in good references, too. In scenarios such as this, we often have to rely on less effective solutions than BBMap.

@SamNicholls if the sample is from a body site (like the human gut) then there is only likely 1 major eukaryotic contaminant (i.e. human DNA; ignoring rare microbial eukaryotes or food), but if you're trawling the ocean or something will be a lot of eukaryotes; think the question needs clarifying

You can use centrifuge with the NT library to profile the taxonomy, and then remove the reads from eukaryotes.

Answer from /u/Romanticon on reddit. The original answer can be found here

Depending on the length of your reads, the level of error in the
  reads, the quality of the database, and the specificity settings you
  give DIAMOND, you can change the level of "bleed" that you get - but
  you'll always have a couple reads, especially at metagenome dataset
  sizes, that will be annotated incorrectly.
But you'll always have this happen, whenever you annotate any large
  dataset. The best way to handle it is to set a threshold after
  annotation that removes wrongly annotated reads (looking for 'clearly
  wrong' organisms in your annotated dataset can help you figure out
  where to set the threshold level - there's probably no Bos taurus
  reads in a mouse metagenome).
You can do things like enable the sensitive flag on the DIAMOND
  annotation search to help improve annotation accuracy.


text:
Converting gene names from one public database format to another
This is a question from /u/apivan19 on reddit. The original post can be found here.
I have some proteomics data that was given to me with the UniProt gene identifiers in column 1. I've been trying to convert these to normal gene symbols using various programs, but it is proving to be difficult.
The Uniprot website does it fairly decently, but it is not able to convert all of them and then adds some unknown genes into my list.
For example, I will give it 5439 genes in UniProt notation, and will say "5420 of 5439 UniProt identifiers have been converted to 5450 gene symbols"... which is ridiculous.
I tried using David to change the symbols, but it returns them to me in some ridiculous, random order and there's no way I can sort... actually there might be but it'll take a second.
What are some of the easiest ways to do this? It's already very time consuming and am looking for simpler solutions

I tend to use Ensembl Biomart for such queries since there are APIs for various programming languages, e.g. biomaRt, and, maybe more interestingly, via a REST API (although it’s a pretty terrible one).
To translate identifiers from different databases, proceed as follows:

Choose database “Ensembl genes”
Choose dataset your desired oganism
Go on “Filters” › “Gene:” › “Input external reference ID list”

Select the chosen source database
Provide a list of IDs, delimited by newline

Go to “Attributes” › “Gene:” › untick “Transcript stable ID”

If Ensembl IDs are desired, leave “Gene stable ID” ticked …
Otherwise untick it; go to “External:”, tick your desired identifier format

Click “Results” at the top left. This gives a preview that can be exported into various formats; alternatively the top centre buttons “XML” and “Perl” provide the query in XML (for SOAP/REST requests) and as a (horrendously formatted) executable Perl script.

The genenames biomart can be very useful for this - it's got a nice user interface for bulk translation, covers a wide variety of identifiers, can access synonyms or defunct names, etc. However, the service is sometimes up and down like a yoyo and some wrapper client libraries are difficult about versions of biomart.

@agapow Absolutely, agreed. :-(

If you're comfortable doing a little programming, check out mygene.info (web services for gene annotations of all sorts).  ID translation is specifically one of the use cases addressed in the bioconductor client (see the vignette), and there is a python client as well available through pypi. The documentation for mygene can be found here.

You can do the same using AnnotationDbi Package from Bioconductor. Download the organism specific annotation file like org.Mm.eg.db for mouse and map current gene ids to the gene names/gene symbols.

My favourite gene database conversion site is db2db. You provide a list of IDs in one of a large number of different public formats, and can select one or more IDs as translation targets. It will then walk through various known paths to do the translation, picking what it determines to be the most reliable route to get the information you have asked for. Results appear in the browser as a table, but can also be exported as an Excel file, or as a tab-separated text file.
Note that the mapping of genes from one database to another is not a one-to-one mapping. It is likely the case that there will be some genes in the source database that map to multiple genes in the target database (and vice versa), and some genes that aren't present in the target database. These phenomena probably account for the "ridiculous" results that have been seen here.

I'm not a huge fan of the Ensembl BioMart system because I find it difficult to use. The Synergizer has a very straightforward interface and works pretty well for most lists. Note: it hasn't been updated in a while.

By far my preferred option for doing this manually is PICR:
http://www.ebi.ac.uk/Tools/picr/
BTW it is not "ridiculous" to get different numbers of genes reported for a given set of proteins. For several reasons:

Uniprot IDs can disappear, merge or be split
not all uniprot and gene IDs have a 1-to-1 relationship
depending on species some gene symbols can be ambiguous or synonymous.


text:
How are snakemake's --cluster and --drmaa options implemented?
I'm fairly new to snakemake and I'm trying to understand the difference between the --cluster and --drmaa flags, both of which are used to submit jobs to compute clusters/nodes.
The docs give a few hints about the advantages of using --drmaa here:

If your cluster system supports DRMAA, Snakemake can make use of that
  to increase the control over jobs. E.g. jobs can be cancelled upon
  pressing Ctrl+C, which is not possible with the generic --cluster
  support.

And here:

If available, DRMAA is preferable over the generic cluster modes
  because it provides better control and error handling.

So I have a conceptual understanding of the advantages of using --drmaa. However, I don't consider the above a very complete explanation, and I don't know how these flags are implemented in snakemake under the hood, can someone elaborate?
Note, that although this could be considered a more general programming question, snakemake is predominantly used in bioinformatics; I was persuaded that this question would be considered on-topic because of this meta-post and this answer.

Hi Chris! What's the actual question here? Seems that DRMAA is an open API describing submission and control of cluster/grid jobs. So if you have a cluster or submission system that has support for it (such as SLURM), snakemake will afford you the additional controls as supported by that API -- including as per your question, the ability to pass on a Ctrl+C at your console to the grid -- a scenario that would require use of the qdel command on Sun Grid Engine, for example. So if you've infrastructure with a scheduler that has DRMAA support, I can see why you might as well use --drmaa.

Maybe wait until we've actually poached the snakemake developers :)

The code is open source and fairly easy to follow, I found it quite easy to work out how DRMAA jobs were submitted by nosing around the source for 10 minutes. I'm not entirely sure what the question is looking for, though.

I'd always kind of wondered how this worked too, so I took this as an excuse to look into the snakemake code. At the end of the day this becomes a question of (1) how are jobs actually submitted and (2) how is it determined if jobs are done (and then whether they failed)?
For DRMAA, python has a module (appropriately named "drmaa") that wraps around the libdrmaa library that comes with most schedulers. This is a very popular route, for example the Galaxy project uses this for dealing with most clusters (e.g., I use it to connect our internal Galaxy instance to our slurm cluster). The huge benefit here is that DRMAA does magic to allow you to simply submit commands to your cluster without needing to know whether you should execute qsub or srun or something else. Further, it then provides methods to simply query whether a job is running or not and what its exit status was.
Using the --cluster command requires a LOT more magic on the side of snakemake. At the end of the day it creates a shell script that it then submits using the command you provided. Importantly, it includes some secret files in that script that it can then watch for (ever noticed that .snakemake directory where you execute a script? That appears to be among the things it's used for.). These are named {jobid}.jobfinished and {jobid}.jobfailed and one of them will get touched, depending on the exit status of your command/script. Once one of those is there, then snakemake can move on in its DAG (or not, if there's a failure). This is obviously a LOT more to keep track of and it's then not possible for snakemake to cancel running jobs, which is something one can do easily with DRMAA.

Having written a package for submitting jobs to Sun Grid Engine without DRMAA, can definitely confirm the magic must be strong here. Though instead of monitoring whether files have been touched, I chose to fire a POST request at the monitoring server. I can see that DRMAA must expose some nice API to get job status without having to do some odd workarounds.

Yup, it's more convenient than what I would have come up with too :)

DRMAA (Distributed Resource Management Application API) appears to be an open API that describes a specification for submission and management of work submitted to some grid/cluster. If your scheduler is DRMAA compliant, my guess would be that using snakemake's --drmaa flag will afford you the additional controls exposed by that API.
As mentioned in your question with DRMAA support, a Ctrl+C on one's console would be passed on to the grid to kill or stop submitted jobs. A task that would otherwise require use of the qdel command on Sun Grid Engine, for example.
I could see why the suggestion would be to enable --drmaa where possible as it's likely to expose functionality (such as control of jobs, and availability of logs and errors) that makes the submission and management of jobs a little easier for you and snakemake. Perhaps someone with snakemake experience can explain a little more, but I hope this helps?
For what it's worth, I use Sun Grid Engine without DRMAA and it's quite a pain to get good information out of qstat, and its qacct log format is quite possibly the worst file format I have ever encountered. I suspect DRMAA provides a nice (or at least more reasonable) API for snakemake (and others) to get at this information more readily.
Appended:
Under the hood, the flags are parsed by Python's argparse module in the __init__.py. A mutually exclusive group of options forces a selection of --cluster or --drmaa (or --cluster-sync).
--cluster or -c requires you to specify an argument; the command for which to submit a job on your cluster (the example names qsub). The --drmaa flag just seems to indicate to snakemake that DRMAA is to be used, which doesn't change a lot in terms of how the command you run from your console is processed in the __init__.py or workflow.py.
However, when it is time to interface with your scheduler (as specified in scheduler.py), an elif statement checks whether you're in some form of cluster mode and the else catches the case where you've raised the --drmaa flag instead.
Here's where the magic happens, as now your job is submitted with the DRMAAExecutor as specified in executors.py, instead of the default executor.
Without investigating much more, I can see that the DRMAAExecutor features some more class attributes and exposes additional functions when compared to the others.
Hope this is of more use? Again, I've only just looked at the code myself but in lieu of a snakemake user chiming in, I thought I'd have a nose.

Thanks, but this tells me what I already know; i have a conceptual understanding but i'd like the details on how the 2 different approaches are implemented in snakemake, i'll try to clarify the question a bit

Sorry I pre-empted what you might have been asking for as I didn't have much more room to write in the comments. I'll try and expand to address your edited question if I can help!

@Chris_Rands I've expanded, is this relevant?

Thanks, upvoted, sorry if I wasn't very clear with my question, I just thought this might be useful question for others,i was planning on answering myself if i didn't get any answers, but i've not delved into the source yet

@Chris_Rands No problem (and thanks)! I enjoy digging around code bases to see how they work. Between Devon and I, have we found what you were looking for? Happy to chat on this with you if not.


text:
Library for computing BWT-based alignments
I am writing a software tool to which I would like to add the ability to compute alignments using the efficient Burrows-Wheeler Transform (BWT) approach made popular by tools such as BWA and Bowtie. As far as I can tell, though, both of these tools and their derivatives are invoked strictly via a command-line interface. Are there any libraries that implement BWT-based read alignment with a C/C++ API?
Python bindings would also be great, but probably too much to expect.

Not C/C++, but just to mention that rust-bio seems to have BWT-related capabilities: https://rust-bio.github.io/

First, let us remark that there exist several hundred read mappers, most of which have been even published (see, e.g., pages 25-29 of this thesis). Developing a new mapper probably makes sense only as a programming exercise. Whereas developing a quick proof-of-concept read mapper is usually easy, turning it into a real competitor of existing and well-tuned mappers can last for years.
It is not clear from the provided description how long is your reference, how many alignments you need to compute, etc. In certain situations, it may be useful to write a wrapper over existing mappers (e.g., using subprocess.Popen for running the mapper and PySam for parsing the output); while in some other situations, standard dynamic programming may be sufficient (e.g., using SSW with its Python binding).
Let assume that you want to develop a toy read mapper. Most of read mappers are based on a so called seed-and-extend paradigm. Simply speaking, first you detect candidates for alignments, usually as exact matches between a read and the reference (using either a hash table or some full-text index – e.g., BWT-index). Then you would need to compute alignments for these candidates, typically using some algorithm based on dynamic programming, and report the best obtained alignments (e.g., the ones with the highest alignment score).
There exist two big, powerful and well debugged libraries implementing BWT-indexes which can be easily used for building a read mapper:

SeqAn. See the FMIndex tutorial and the Pairwise Sequence Alignment tutorial for quick examples of how to detect exact matches and how to do pairwise alignments. Also, they provide a tutorial about a quick development of a read mapper, but the resulting mapper seems to be hash-table based, not BWT-based. SeqAn is used, e.g., in YARA.
SDSL-Lite. This is a general library for succinct data structures. See the tutorial slides for an idea of how it works. For instance, GramTools are based on SDSL.

Never heard of YARA … better-known mappers also use SeqAn, e.g. both TopHat and Bowtie (though not for the mapping itself, as far as I know).

Sigh, I was hoping to avoid writing BWA's SAM output to disk, but PySAM treats input strings as file names and explicitly doesn't support StringIO objects. Guess I'll have to look into named pipes or something like that.

PySAM supports also file objects – see http://pysam.readthedocs.io/en/latest/api.html#pysam.AlignmentFile. AlignmentFile(filepath_or_object,....: "If filepath_or_object is a python File object, the already opened file will be used."

Yeah, but running proc.communicate gives me BWA's SAM output as a string. It looks like none of the usual string trickery will work, and somehow I'll have to go through the file system to get PySAM to read it.

BWA-MEM can be used as a library. File bwa/example.c shows the basic functionality for single-end mapping. It should give identical mapping to the bwa-mem command line. Header bwa/bwamem.h contains basic documentation. Paired-end mapping is doable, but not well exposed. Several teams, including GATK and MS genomics, have been using bwa-mem as a library.

That is extremely cool.

SeqAn supports BWT tables for use with their parametrisable alignment algorithms.
To use it, follow the general outline for building a SeqAn short-read aligner, and use the FMIndex specialisation instead of — as in the example — the IndexQGram.

If you're looking for a python API, nanopore have created one (although I am not sure how "maintained" it is) https://github.com/nanoporetech/bwapy


text:
What are some good practices to follow during EPIC DNA methylation data analysis?
I recently got some EPIC DNA methylation data and I was wondering what are some good practices to follow?
I am interested in knowing about normalization and differential analysis. Thank you.

Hi deepseas, welcome and thank you for posting a question to Bioinformatics StackExchange. I've added a link to the Infinium EPIC web page for people who might be curious about what EPIC DNA methylation data is. Your question is quite broad, and would be substantially improved by adding more of a story behind why you want to do this analysis (e.g. is this a single individual study? do you want to compare with RNASeq data?). Bioinformatics is a very large subject area, and there are many different ways to do [almost] the same thing, so any help with clarifying questions is greatly appreciated.

Thanks for the response! This data has 16 samples from different individuals with 5 groups (5+6+1+2+2). It will be difficult to get power for groups with 1 and 2. but the main aim is to get an idea about certain known genes of interest.

It will be difficult to make contrasts like edgeR,limma or deseq2 makes with low number of replicates per condition. 2 might still work but 1 in a single condition will not hold good. Ideally when you work with EPIC arrays or 450k both will work with any standard 450k analysis tool. Try tools like DimMER/Rnbeads. However things hold better with higher replicates per groups. All the best!

EPIC data can be processed in the same manner as the previous iteration of methylation array data from Illumina (450k). This means that starting with .idat files, normalization should be performed (for example, via the minfi package). A recent paper from the creators of minfi is particularly helpful because it makes clear that normalized EPIC data from their package can be immediately compared against, for example, level 3 TCGA data.
After that, I suggest using the manifest to attach genomic coordinates to your probes and segregating them into functional regions. By testing differential methylation in only regulatory regions, for example, you can increase the statistical power by reducing the overall number of tests to the ones you expect to yield major differences.
There are existing packages out there for differential methylation analysis, but without knowing your replicate structure or aims, it is difficult to point you in the right direction.


text:
Obtaining uniquely mapped reads from BWA mem alignment
This is based on a question from betsy.s.collins on BioStars. The original post can be found here.

Does anyone have any suggestions for other tags or filtering steps on BWA-generated BAM files that can be used so reads only map to one location? One example application would be to find seeds for the TULIP assembler/scaffolder, which works best for reads that map to unique genomic locations.

The questioner is referring to genomically-unique alignments, which are something different from multiple possible alignments at a single location. If the only multiple alignment is at the same genomic location, those alignments are still unique in the genomic sense. The impossibility of finding the exact, correct, alignment is a well-known problem, and there are a few downstream methods (e.g. left normalisation) to make sure that multiple local alignments and/or sequencing errors have reduced impact.
The concept of "uniquely mapped reads" is a loaded term, and most sources suggest filtering by MAPQ should do the trick. However, this approach doesn't seem to work when using BWA as a read mapper.
Uniquely mapped reads (i.e. a read that maps to a single location in the genome) are sometimes preferred when running analyses that depend on the quantification of reads, rather than just the coverage (e.g. RNASeq). Duplicated reads require additional filtering or normalisation on the analysis side, and most downstream programs won't consider concepts like "a fifth of the read maps somewhere here on chromosome 5, and four-fifths of the read maps somewhere here on chromosome 14."

Update - as of January 2021, samtools can now do filtering based on an expression that includes tag variables. In this case, this expression can be used to exclude any reads that have either an XA or SA tag:
samtools view -b mapped.bam -e '!([XA] | [SA])' > unique_mapped.bam

For more details on the samtools expression parameter, see the samtools documentation:
https://www.htslib.org/doc/samtools.html#FILTER_EXPRESSIONS

Original answer follows....
To exclude all possible multi-mapped reads from a BWA-mapped BAM file, it looks like you need to use grep on the uncompressed SAM fields:
samtools view -h mapped.bam | grep -v -e 'XA:Z:' -e 'SA:Z:' | samtools view -b > unique_mapped.bam

Explanation follows...
I'm going to assume a situation in which a bioinformatician is presented with a mapped BAM file produced by BWA, and has no way of getting the original reads. One high-effort solution would be to extract the mapped reads from the BAM file and re-map with a different mapper that uses the MAPQ score to indicate multiple mappings.
... but what if that were not possible?
My understanding of BWA's output is that if a read maps perfectly to multiple genomic locations, it will be given a high mapping quality (MAPQ) score for both locations. Many people expect that a read that maps to at least two locations can have (at best) a 50% probability of mapping to one of those locations (i.e. MAPQ = 3). Because BWA doesn't do this, it makes it difficult to filter out multiply-mapped reads from BWA results using the MAPQ filter that works for other aligners; this is likely to be why the current answer on Biostars [samtools view -bq 1] probably won't work.
Here is an example line from a BWA mem alignment that I've just made. These are Illumina reads mapped to a parasite genome that has a lot of repetitive sequence:
ERR063640.7     16      tig00019544     79974   21      21M2I56M1I20M   *       0       0       TATCACATATCATCCGACTCAGCTCGACGAGTACAATGCTAATTTAACACTTAGAATGCCCGGCAATGAAATTCGTTTTCCGTCAATTCTTGAAAATTTC    <AABBEGABFJKKKIM7GHKKJK>JLKLDGMHLIMIHHCGIJKKLJKLNJGLLLKLILKLMFNDLKGHJEKMKKMIJHGLOJLLLKIJLKKJEJLIGG>D    NM:i:4  MD:Z:83A13      AS:i:77 XS:i:67 XA:Z:tig00019544,-78808,21M2I56M1I20M,6;tig00019544,-84624,79M1I20M,6;tig00019544,-79312,33M4I42M1I20M,8;

BWA mem has found that this particular read, ERR063460.7, maps to at least three different locations: tig00019544, tig00019544, and tig00019544. Note that the MAPQ for this read is 21, so even though the read maps to multiple locations, MAPQ can't be used to determine that.
However, the alternative locations are shown by the presence of the XA tag in the custom fields section of the SAM output. Perhaps just filtering on lines that contain the XA tag will be able to exclude multiply-mapped reads. The samtools view man page suggests that -x will filter out a particular tag:
$ samtools view -x XA output.bam | grep '^ERR063640\.7[[:space:]]'
ERR063640.7 16  tig00019544 79974   21  21M2I56M1I20M   *   0   0   TATCACATATCATCCGACTCAGCTCGACGAGTACAATGCTAATTTAACACTTAGAATGCCCGGCAATGAAATTCGTTTTCCGTCAATTCTTGAAAATTTC    <AABBEGABFJKKKIM7GHKKJK>JLKLDGMHLIMIHHCGIJKKLJKLNJGLLLKLILKLMFNDLKGHJEKMKKMIJHGLOJLLLKIJLKKJEJLIGG>D    NM:i:4  MD:Z:83A13  AS:i:77 XS:i:67

... so it filtered out the tag (i.e. the tag no longer exists in the SAM output), but not the read. There are no useful bits in the FLAG field to indicate multiple genomic mappings (which I know can be filtered to exclude the read as well), so I have to resort to other measures.
In this particular case, I can use grep -v on the uncompressed SAM output to exclude alignment lines that have the XA tag (and re-compress to BAM afterwards, just to be tidy):
$ samtools view -h output.bam | grep -v 'XA:Z:' | samtools view -b > output_filtered.bam
$ samtools view output_filtered.bam | grep '^ERR063640\.7[[:space:]]'
<no output>

Hurray! reads filtered. As a little aside, this grep search has a fairly substantial computational load: it's looking for some string with the text XA:Z: somewhere in the line, and doesn't actually capture every situation. Some masochistic person might come along at a later date and decide that they're going to call all their reads HAXXA:Z:AWESOME!:<readNumber>, in which case a tweak to this grep search would be needed to make sure that there's a space (or more specifically a tab character) prior to the XA:Z:.
Now I do a check for any duplicated read names, just to be sure:
$ samtools view output_filtered.bam | awk '{print $1}' | sort | uniq -d
ERR063640.1194
ERR063640.1429
ERR063640.1761
ERR063640.2336
ERR063640.2825
ERR063640.3458
ERR063640.4421
ERR063640.4474
ERR063640.4888
ERR063640.49
ERR063640.4974
ERR063640.5070
ERR063640.5130
ERR063640.5300
ERR063640.5868
ERR063640.6116
ERR063640.6198
ERR063640.6468
ERR063640.6717
ERR063640.6797
ERR063640.7322
ERR063640.750
ERR063640.7570
ERR063640.7900
ERR063640.8115
ERR063640.8405
ERR063640.911
ERR063640.9206
ERR063640.9765
ERR063640.9986

Oh... damn. I wonder what they are:
$ samtools view output_filtered.bam | grep '^ERR063640.3458[[:space:]]'
ERR063640.3458  16  tig00002961 5402    60  58S38M  *   0   0   AGGTACCATTCGATAGAGGGAGAAAGGCACTACTAAAGATTTTGCCACATTTGCTATATCCGTATCGCGAAGATCAGGACTTACTCCGCAGAAGAA    DD6HFFJBKFH=KDILKLGLJEKLKGFJIH8IKHLLMJEK:L:HBGJIHJKFLLKIHJDHLNKCK;KMKGMFKJILIIIMKI9JLKKHEJFII?CC    NM:i:0  MD:Z:38 AS:i:38 XS:i:0  SA:Z:tig00002377,202353,-,14M3I5M1I35M38S,19,5;
ERR063640.3458  2064    tig00002377 202353  19  14M3I5M1I35M38H *   0   0   AGGTACCATTCGATAGAGGGAGAAAGGCACTACTAAAGATTTTGCCACATTTGCTATA  DD6HFFJBKFH=KDILKLGLJEKLKGFJIH8IKHLLMJEK:L:HBGJIHJKFLLKIHJ  NM:i:5  MD:Z:5G48   AS:i:35 XS:i:27 SA:Z:tig00002961,5402,-,58S38M,60,0;

Aha! Supplemental alignments, which use the official SA tag [other canonical alignments in a chimeric alignment]. These appear to be situations where a single read has been split up, and maps to multiple locations. Note that in this case, both of the alignments still have MAPQ scores of over 3. It sounds like the questioner would also want to get rid of these situations as well. This time, there are standard flag fields as well to deal with these situations (0x800: secondary alignment). Except it's not enough to just filter the supplemental alignment, because both read mappings should be removed, rather than just the one (or ones) that happened to be tagged as secondary.
Luckily, BWA appears to put the SA tag into all reads containing supplementary alignments (if this is not the case, I'm sure someone will correct me on that). So, I add in the SA search as an additional grep filter:
$ samtools view -h output.bam | grep -v -e 'XA:Z:' -e 'SA:Z:' | samtools view -b > output_filtered.bam
$ samtools view output_filtered.bam | awk '{print $1}' | sort | uniq -d
<no output>

Done. Easy peasy! </sarcasm>
... that original "high-effort" solution of using a different aligner doesn't look so bad now.

that's a neat solution, and well argued. I don't think it's that high effort if you condense the whole thing to a few lines of shell script either. :)

The thing with using a MAPQ to filter out multimappers is that it uses the presumption of a multimapper being "equally well mapping", which for even the examples presented here won't be the case. I've set your reply as the accepted answer on biostars but asked op there to clarify a bit on what exactly is wanted.

Yes, it's not hard if you know what to expect, which is why I've tried to run through all the problems that might happen.

About the read names possibly containing "XA:Z:", a way to achieve robustness would be to use a proper SAM format parser, such as pysam, if one programs in python, rather than grep. I don't know how much slower this would be.

Or like samtools view, which I tried, but didn't do what I had expected. I can imagine that samtools might implement something like a tag filter in the future. Maybe it already has, and I'm just using the wrong sub-tool.

Okay, so it's not yet available in samtools, but the idea of doing tag filtering has already been raised in at least one issue

Worst case, you can write a couple lines of python or likely use Pierre's jvarkit to filter things out.

For much faster and stable$ implementation, sambamba can handle this using following one-liner:
sambamba view -t 12 -h -f bam -F "mapping_quality >= 1 and not (unmapped or secondary_alignment) and not ([XA] != null or [SA] != null)" test.bwa-mem.bam -o test-uniq.bam

See manpage and sambamba synatx for filters for details on further usage.
$: Besides grep based search being slow, it may return false positive negative hits (@gringer's reply above). I also do not find awk based approach reliable because fields like XA, SA, etc. are optional fields in SAM format, and not positionally fixed fields.

In what way were the hits false positive?

Hmm.. I guess it should be false negative because you were excluding those search patterns with grep -v, right?

The only likely way grep could produce the wrong results is if someone modified the read names or a read group to contain one of those strings, which seems quite unlikely.

agree, it is quite unlikely but nonetheless, grep based approach is inherently slow.

The slow part is in the conversion from BAM to SAM by samtools view. The grep program is very quick at what it does.

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


text:
State of the art in predicting Translation Initiation Sites
I'm working on a university project of predicting Translation Initiation Sites in human DNA. I searched the net for papers and documentation to get guidelines and inspiration, but I feel uncertain that I was able to find the state of the art in solving this problem, so I thought of asking the community here: What's the state of the art method of predicting Translation Initiation Sites?
(I'm also posting the question to help this nascent Q&A site by participating, please close it if it completely misses the mark of being a good question.)

Ribosome Profiling experiments profile the positions of ribosome-protected mRNA. a cycloheximide treatment fixes the translating ribosomes. A naive approach to define Translation Initiation Sites (TIS) would be to look at positions where you see distinct 'peaks' in the distribution of these fragments since the ribosome density will be higher at the start of the CDS (at the TIS).
TISDB is a database comprised of GTI-Seq global mapping of TIS codons at nearly single-nucleotide resolution. GTI-Seq uses two 'fixers' that fix ribosomes where one of them specifically fixes the initiating ribosomes.


text:
What's the scaling for HOMER metagenes?
I'm trying to use HOMER to make a metagene profile over gene bodies using a bedgraph file I've generated. The problem is that every time I do, I get really weird scaling on the y-axis. I should be getting average values across the gene body on the order of 5-10, but instead I'm getting values on the order of 0.03-0.05 or less. The weird thing is, when I don't use metagene script or when I don't use histogram with -size given, I get perfectly normal values -- but that's not what I want, unfortunately, for the metagene profile.
What am I missing?

Figured it out! For anyone who finds this thread in the future and is wondering what is going on: HOMER, when using the -size given option, normalizes the y-axis to the number of basepairs in the intervals -- in this case, the length of the gene body. To get the "raw" values, as I wanted, you need to multiply by the length of the gene body to "undo" the normalization.


text:
Building STAR Genome Index for nanopore RNA sequencing
I am aligning a dataset of 1,000,000 reads oh human mRNA sequenced on Oxford Nanopore Technologies' MinION, and would like to use the STAR aligner, using the parameters recommended by Pacific Biosciences for long reads.
According to this Google Groups thread, in setting up the genome index for short reads, the parameter sjdbOverhang should be set to 1 less than the read length. Obviously, with long (mean 1.7Kb, max >50Kb) reads, this doesn't make sense.
Running STARlong --runMode genomeGenerate without setting sjdbOverhang sets the parameter to 0. Does anyone know what this means, how it might affect mapping, and what I should set it to for my long reads?

Not an answer, but for full-length transcripts, it is worth considering traditional cDNA mappers such as gmap and blat. How well they perform highly depends on the error rate. They are probably ok at ~5% error rate. If you are using the older chemistry, they may not work, though.

A very good point. We're using GMAP as our baseline aligner, but I wanted to try out STAR to confirm the rather surprising results in this bioRxiv preprint: http://biorxiv.org/content/early/2017/04/11/126656

The parameter is used to determine how much sequence STAR indexes on each side of a splice junction to improve its alignment accuracy.  For very long reads, this may not be ideal. I am not sure if STAR is capable of including multiple splice junctions since a long read is more than likely to span more than one.
It may be worthwhile to consider aligning your reads directly to transcripts and taking the reads that align poorly and using them to look for novel transcripts.

I've found the STAR manual to be incredibly helpful with trying to navigate my way round all the STAR command line parameters. Here's the section on the sjdbOverhang parameter:

--sjdbOverhang specifies the length of the genomic sequence around the annotated junction to be used in constructing the splice junctions database. Ideally, this length should be equal to the ReadLength-1, where ReadLength is the length of the reads. For instance, for Illumina 2x100b paired-end reads, the ideal value is 100-1=99. In case of reads of varying length, the ideal value is max(ReadLength)-1. In most cases, the default value of 100 will work as well as the ideal value.

[emphasis preserved from the original text]
I would recommend trying out running STAR with this value unset, or set to 100, and seeing how it goes. If you want to be particularly adventurous, you could set it to 1000 and see if anything changes.


text:
Entrez or Ensembl gene IDs?
I'm using several datasets that are encoded using either Entrez or Ensembl IDs to specify genes, and need to decide on which to standardise on.
Are there any major reasons to use one over the other?
Which tools and conversion tables are best to use for this?

It may be useful to specify some more about your particular use case. I don't have an answer to your question, but it's feasible that the pros and cons of the two databases differ depending on what you're trying to do.

There is a discussion on Biostar in the same topic. https://www.biostars.org/p/16505/

Hi einonm, thanks for your question, and welcome to Bioinformatics Stack Exchange. Bioinformatics is a big subject area with a vast array of tools and databases, many which have almost the same purpose but are designed to attack a particular problem. As @roblanf has pointed out, it would be helpful if you could add a bit more story / context to your question so that the community can better understand the problem. A more specific question means the chance of a good, relevant answer is increased.

It's not just a matter of what IDs to use, but of what gene models to use. Even if an Ensembl gene has a cross reference to an Entrez gene, that's only because the Ensembl xref pipeline has decided that they are similar enough.

I think there is no general answer to your question.
Which type of ID you will use depends on how practical it will be for your downstream processing. This depends on what other source of information you have, how much conversion work you will need to do in order to match your data to these information if you use one or the other.
I think the choice may emerge by trial and error, as you realize that a given choice is actually not convenient in your particular use case. Try to set up easy to use conversion procedures that you can easily insert in some steps of your workflow, so that you can easily switch from one choice to the other.
Concerning the tools to do the conversion, see the following question and its answers: Converting gene names from one public database format to another


text:
How can HISAT2/StringTie report decimal coverage values
I have performed RNA-seq analysis using HISAT2 & StringTie workflow suggested in: Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown.
Some of the transcripts/exons have decimal coverage values (eg., 1.1 or 2.59).
My question is: How can these tools report decimal coverage? If only half of the read overlaps exon will this exon have 0.5 coverage?

I suspect this is average coverage across the transcript or exon.
Formula for average coverage:

Where:

N: number of reads mapping/aligning to your transcript/exon
L: average read length
T: Transcript/exon length


text:
BAM to BigWig without intermediary BedGraph
I have a pipeline for generating a BigWig file from a BAM file:
BAM -> BedGraph -> BigWig

Which uses bedtools genomecov for the BAM -> BedGraph part and bedGraphToBigWig for the BedGraph -> BigWig part.
The use of bedGraphToBigWig to create the BigWig file requires a BedGraph file to reside on disk in uncompressed form as it performs seeks. This is problematic for large genomes and variable coverage BAM files when there are more step changes/lines in the BedGraph file. My BedGraph files are in the order of 50 Gbytes in size and all that IO for 10-20 BAM files seems unnecessary.
Are there any tools capable generate BigWig without having to use an uncompressed BedGraph file on disk? I'd like this conversion to hapen as quickly as possible.
I have tried the following tools, but they still create/use a BedGraph intermediary file:

deepTools

Some Benchmarks, Ignoring IO
Here are some timings I get for creating a BigWig file from a BAM file using 3 different pipelines. All files reside on a tmpfs i.e. in memory.
BEDTools and Kent Utils
This is the approach taken by most.
time $(bedtools genomecov -bg -ibam test.bam -split -scale 1.0 > test.bedgraph \
  && bedGraphToBigWig test.bedgraph test.fasta.chrom.sizes kent.bw \
  && rm test.bedgraph)

real    1m20.015s
user    0m56.608s
sys     0m27.271s

SAMtools and Kent Utils
Replacing bedtools genomecov with samtools depth and a custom awk script (depth2bedgraph.awk) to output bedgraph format has a significant performance improvement:
time $(samtools depth -Q 1 --reference test.fasta test.bam \
  | mawk -f depth2bedgraph.awk \
  > test.bedgraph \
  && bedGraphToBigWig test.bedgraph test.fasta.chrom.sizes kent.bw \
  && rm test.bedgraph)

real    0m28.765s
user    0m44.999s
sys     0m1.166s

Although it is has less features, we used mawk here as it's faster than gawk (we don't need those extra features here).
Parallelising with xargs
If you want to have a BigWig file per chromosome, you can easily parallelise this across chromosome/reference sequences. We can use xargs to run 5 parallel BAM->BedGraph->BigWig pipelines, each using the tmpfs mounted /dev/shm for the intermediary BedGraph files.
cut -f1 test.fasta.chrom.sizes \
  | xargs -I{} -P 5  bash -c 'mkdir /dev/shm/${1} \
  && samtools depth -Q 1 --reference test.fasta -r "${1}" test.bam \
  | mawk -f scripts/depth2bedgraph.awk \
  > "/dev/shm/${1}/test.bam.bedgraph" \
  && mkdir "./${1}" \
  && bedGraphToBigWig \
    "/dev/shm/${1}/test.bam.bedgraph" \
    test.fasta.chrom.sizes \
    "./${1}/test.bam.bw" \
  && rm "/dev/shm/${1}/test.bam.bedgraph"' -- {}

deepTools
Let's see how deepTools performs.
time bamCoverage --numberOfProcessors max \
  --minMappingQuality 1 \
  --bam test.bam --binSize 1 --skipNonCoveredRegions \
  --outFileName deeptools.bw

real    0m40.077s
user    3m56.032s
sys     0m9.276s

It's been on my to do list to remove the "make a single bedGraph file" intermediate step from deepTools. Given how that's all implemented there's no good way to avoid at least some intermediate file creation (without blowing up the memory at least).

I needed to generate the intermediate bedgraph files when I was creating BigWig files for GBrowse in 2011; it's disappointing to see that it hasn't changed much in the many years since then. The file seeking is a killer -- it excludes most of the shell tricks that can be used to work around programs that don't explicitly support the common UNIX pipeline (e.g. <(), FIFO buffers, /dev/stdin).

@DevonRyan when you say blow up memory, will that be equivilent to the size of the intermediary bedgraph file? I have a 720Gbyte fat memory compute node...so for me, speed is a bigger issue than memory. Perhaps an option to state the max memory usage, at which point you start going to disk?

@NathanS.Watson-Haigh Yeah, it'd be equivalent more or less to the file size. deepTools itself was never engineered to make this easier to implement and the generally crappy nature of multithreading in python doesn't help that :(

@DevonRyan how about using the system's temp dir as the default location for the intermediary bedgraph file and expose this as a command line argument. That way I could use my raided SSD's for more convenient and faster IO. I was doing this with my bedtools and bedGraphToBigWig pipeline.

The system temp dir (or shared memory) is already what's used for all of the tmp files. Honestly, the current implementation is a hold-over from before I wrote pyBigWig (or even joined the group) and deepTools was using bedGraphToBigwig in a system call. It's mostly a matter of me finding time to change how bamCoverage does this final step.

@DevonRyan is there any reason a named pipe as suggested by my answer wouldn't work for this? I have never used these tools so I can't be sure.

@terdon: Sorry, I was referring specifically to deepTools, which was mentioned in the post. One can certainly do what you've described, it'll just be slower than what deeptools can do (but better with IO, so maybe it's all a wash at that point!).

@terdon I don't see your named pipe based answer

@bli no I deleted it since it won't work here because the question requires the ability to seek the file. FIFOs won't work with seeking AFAIK, so the answer wasn't useful. I hadn't realized the OP required a seekable file.

This can be done in R very easily from an indexed .bam file.
Given single-end file for sample1.
library(GenomicAlignments)
library(rtracklayer)

## read in BAM file (use readGAlignmentPairs for paired-end files)
gr <- readGAlignments('sample1.bam')

## convert to coverages
gr.cov <- coverage(gr)

## export as bigWig
export.bw(gr.cov,'sample1.bigwig')

Be aware that this method doesn't include normalization steps (such as normalizing to total coverage). Most of these additional steps can be added if necessary.

I had hoped for a command line tool but I'm not totally against going into R if needs be. I need to test more extensively with regards to speed so I need to have more of a play. If tests show that it's slow for the number and sizes of BAM files I have, or a better command line answer materialises, I will accept this answer.

Just for reference, there is a not very in depth, but first-hand reasoning (Jim Kent) given here about why bedGraphToBigWig does not support streaming http://genome.soe.ucsc.narkive.com/2S1Z3VpG/bedgraphtobigwig-reading-in-from-stdin
Edit: as noted wigToBigWig allows bedgraph streaming input via stdin but takes (possibly much) more memory than bedGraphToBigWig itself does
cat file.bedgraph | wigToBigWig stdin chrom.sizes output.bigwig

Note: You can also use wiggletools bin 10 for example to bin rows and downsample

Ah! "However if you do have enough memory, you can still use wigToBigWig on bedGraphs, and it will fit in a pipe."

That's sort of a weird decision. The zoom levels are written after the raw entries, so one can write all of the raw interval entries to disk and then do another pass over them when making the zoom levels. Then you don't need to store everything in memory. That's how libBigWig works, in fact.

The QuasR package in R has a useful function called qExportWig.
To use it you first need to create a sample sheet and a qProject object - this can be done starting with bam files or fastq files (in which case it also  performs alignments). You can then create wig or bigWig files without writing a bedGraph file to disk.
An example would look something like this:
library(QuasR)

# Create qProject object using qAlign
# Sample file should contain tab-separated file names (BAM or fastq) and sample names
myproj <- qAlign("samples.txt", genome=mygenome, paired="no")

# Export all samples in project as bigWig, scaling to mean number of alignments by default
qExportWig(myproj, binsize=50, createBigWig=T)

It has a lot of flexible options, including scaling or log-transforming counts, shifting reads and changing the bin size.

It would be helpful to provide some example code to demonstrate how one can go from BAM to BigWig. Like the answer provided by @story


text:
Where is an up to date miRNA database and what happened to miRBase?
miRBase 21 was published June 26, 2014 and was still in its growth phase. Why is it not being updated anymore or the project declared officially dead? ENSEMBL also uses miRBase as a starting point (http://www.ensembl.org/info/genome/genebuild/ncrna.html).
What do people use now as a good reference for miRNAs especially for less well annotated/sequenced species?

Hi Ido Tamir, welcome to Bioinformatics Stack Exchange, and thanks for your question. The Stack Exchange community tries to encourage the asking of specific questions that can have definite, objective answers. It would be appreciated if you could modify your question to fit this format, either by adding more context/story around why you're asking the question, or by reducing the open-endedness of the question. Bioinformatics is a very large subject area, and there are frequently many different ways to do the same thing.

what part is unspecific in asking for an up to date miRNA database? You mean I should specify that the last update should have been on the 5. june 2016?

A quick internet search for me returns at least 4 different databases in the first page: 1 -- miRBase; 2; 3; 4. This is a very common occurrence in bioinformatics, where multiple solutions are created to solve very similar problems. The particular database that is used will depend on your particular application, as each will have their own strengths and weaknesses (e.g. a tradeoff between sensitivity and specificity).

miRbase is a database of miRNA sequences, their genomic locations and the evidence for their existence. It is not a database of their targets, so I am assuming that this is the information you are asking for.
Despite not having been updated since 2014, as far as I am aware miRBase is still the go to place for this data. For example, courses at the European Bioinformatics Institute were still teaching the use of miRBase in April this year (2017).
RNACentral is repository for all non-coding RNAs, but as it is a synthesis of other databases, including miRBase, I suspect data their will still be largely the same as miRBase.

If you are still interested, last year miRBase generated new updates. Currently, according to ftp site the last release is 22.1. So, it is not a dead project and for more specific information you should reference the miRBase blog.
In order to get the best set of miRNA annotations, you have to define your target species or even if exists some related annotated species and then, retrieve the current annotation along different miRNA databases, but most importantly you have to be really clear about the definition of miRNA families, because it is clearly different in miRBase and also in RFAM.

I mostly use TargetScan predictions. But here are few curated databases:

miRGate

miRGate is a curated database of human, mouse and rat miRNAs/mRNAs targets. It is designed to analyze miRNA and gene isoforms lists under a common and consistent space of annotations. Including all existing 3 UTR and the entirely known miRNAs. All Havana biotypes and ENCODE principal isoforms for the three organisms are also included. Publication

miR2Disease :

a manually curated database, aims at providing a comprehensive resource of miRNA deregulation in various human diseases. Publication

PASmiR

PASmiR, a literature-curated and web-accessible database, was thus developed to provide a solid platform for the collection, standardization, and searching of these miRNA-stress regulation data in plants. Publication

For miRNA precursor predictions in species in Ensembl but not in miRBase, you can have a go to MapMi (http://www.ebi.ac.uk/enright-srv/MapMi/) .
You can easily run it locally for other species or if you have smallRNA sequencing done on something new.
Let me know if you run into difficulties or have ideas to improve it.
Ensembl itself has its own miRNA predictions.

I would guess DIANA-Tarbase ( evidences, https://www.ncbi.nlm.nih.gov/pubmed/25416803 ) and TargetScan (predictions, https://www.ncbi.nlm.nih.gov/pubmed/26267216 ) are good starting points.
Especially for the less annotated species, TargetScan might be helpful.

thank you very much, but I am not interested in miRNA-target predictions or interactions but a curated database of miRNAs (-5p/-3p , precursor, genomic location, families). targetScan - i think - also relies on miRBase. Still one vote for the effort.

It seems that TargetScan uses more miRNAs than those in miRBase (see FAQ number 11: http://www.targetscan.org/faqs.Release_7.html).