text:
Is it possible to merge scRNAseq data from experiments with different design?
I have 4 different single-cell RNAseq experiments, each one representing a different sample of cell type population. I'd like to merge them to a single dataset. However, different cell types are enriched in each sample, and different samples can represent also different conditions (control and disease).
If I run a PCA on the merged dataset (without any processing of data) I find that the different sample's cells cluster separately. This could be due to biological and technical (i.e. batch effects) variation.
Is it possible to merge these different samples while retaining a good amount of information or should they be analyzed separately?

Is there some samples (controls perhaps) in both experiments? What is your goal by merging them: compare a cell line of experiment A with a cell line in experiment B?

No controls unfortunately. My goal by merging them is basically to find different sub-populations in the data and markers by exploiting a larger dataset. This data was presented to me after sequencing, and I think that it is only possible to use each experiment separately, but I'd like to be proven wrong - at least in some specific analysis.

I'm not familiar enough with the methods for scRNA, sorry. But why do you think that increasing the number of cells (of different experiments) will improve your statistical power to find new/different sub-populations? (Perhaps is specific of the methods, but I think they should already be present in each individual dataset)

I am thinking of sampling cells from different part of the same organ as a proxy for spatial resolution, for instance. The ability of merging different datasets may give the opportunity to compare the different cells populations with added spatial heterogeneity distribution. But however, this is just an example.. in my case I was wondering if merging could add more information.

How different is the experimental design? cellranger supports “reanalyse” across samples but I’m not sure this applies with different platforms generating the data.

@TomKelly they are generated with the same platform on 10X, but they are enriched differently using markers from FACS, so I think it is not possible.

@TomKelly which is the tool within cellranger to merge different samples? I was looking at https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/reanalyze but it seems to me it cannot be used to merge samples

I remember attending a seminar where scientists merged data from different types of sc experiments (i think from the same platform though). Its long ago, so my memory might not be correct on this, but I do recall that it would be a functionality planned to be implemented in Seurat - although I haven't yet had any luck finding documentation on https://satijalab.org/seurat/. In addition I think that someone at the seminar mentioned using WGCNA to map batch effects. Since my memory is not clear on this subject I don't dare post this as an answer though!

This can possibly be of interest:
https://www.ncbi.nlm.nih.gov/pubmed/29608177 - Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors

It’s right there on the cellranger manual:
#aggregate results of counts for separate samples
cellranger aggr
#analyse the combined results
cellranger reanalyze

Note: I’m not sure this applies to the above case where batch effects could be an issue. However, this is how the analysis is performed for 10X data from multiple runs in principle.
cellranger aggr
cellranger reanalyze

Adding to this answer months later: it does not work with different designs, it just runs depth normalization but it can be severely affected by batch effects.

Yes batch effects must be considered. Since this question was specifically about combining datasets, that is what was described here. This is the first step to adjust for batch effects, there are several methods to do this on an aggregate dataset and there are separate questions on this already (these assume some common cells between samples or it will overfit). I recommend not to use downsampling but aggregate raw. Depth normalisation does not remove batch effects or background noise in each sample.

I agree. As you said, methods like ComBat are not suitable for experiments with different design because they require a balanced design (some common cells between samples). Thanks.

The 3.0 release of the Seurat R package (currently in alpha) will have extensive support for combining data from different experiments.

See this paper from the Marioni group, where they propose a method for correcting batch effects between single cell sequencing experiments when each experiment contains different sub-populations:

Large-scale single-cell RNA sequencing (scRNA-seq) data sets that are
  produced in different laboratories and at different times contain
  batch effects that may compromise the integration and interpretation
  of the data. Existing scRNA-seq analysis methods incorrectly assume
  that the composition of cell populations is either known or identical
  across batches. We present a strategy for batch correction based on
  the detection of mutual nearest neighbors (MNNs) in the
  high-dimensional expression space. Our approach does not rely on
  predefined or equal population compositions across batches; instead,
  it requires only that a subset of the population be shared between
  batches. We demonstrate the superiority of our approach compared with
  existing methods by using both simulated and real scRNA-seq data sets.
  Using multiple droplet-based scRNA-seq data sets, we demonstrate that
  our MNN batch-effect-correction method can be scaled to large numbers
  of cells

Have you tried seurat3? Here is the reference:
https://satijalab.org/seurat/v3.0/pancreas_integration_label_transfer.html


text:
When working with Bioconductor and the curatedOvarianData dataset, what are the unit values of the expression data?
I am wondering how to interpret values inside the gene expression levels contained within Bioconductor. For example, if we have the following commands via R, we get:
source("https://bioconductor.org/biocLite.R")
biocLite("curatedOvarianData")
library(curatedOvarianData)
data(TCGA_eset)
head(exprs(TCGA_eset)["CXCL12",])

TCGA.20.0987 TCGA.23.1031 TCGA.24.0979 TCGA.23.1117 TCGA.23.1021 TCGA.04.1337
    5.437923     6.307250     4.193748     9.071293     5.101532     6.697770

The above calls out gene expression levels for the gene CXCL12. However, I am wondering how to interpret the values 5.437923, 6.307250, etc.
I know they are expression levels, but do they come with units? Do larger values mean greater cancer possibility?

These values represent the expression of a particular gene in each sample. While expression of an immune system gene like CXLC12 might be associated with the aggressiveness of a cancer, that's not going to be the case for the majority of genes.
Looking at the numbers, these look like they are probably log$_2$-transformed values. In order to know for sure, plotting a density distribution of the dataset is a good idea:
plot(density(exprs(TCGA_eset)));

If the values are in the range of around 0..15, then it'll be log$_2$-transformed values, whereas if they're in the range of 0..100000, it'll be "raw" expression values. The log-transformed values are preferred for data analysis because gene expression tends to be close to a normal distribution in log-space.
These numbers are best treated as relative counts (i.e. the most appropriate unit is "expression"), as variation in the preparation of the sample (or the sequencer, or the microarray) can influence the results that come out the other end. Converting them to something more meaningful (e.g. "transcripts per cell") requires normalising against a known reference sample.

The values are rma normalized expression values from affy microarrays. They are indeed log2 transformed (as @gringer already expected). You can find more info about the data here, under heading TCGA_eset.
To get more info about rma, check this out.
Larger values are not necessarily greater cancer possibility, I wish biology was that simple.

source: Stack Exchange

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

sort: votes


text:
What is the fastest way to get the reverse complement of a DNA sequence in python?
I am writing a python script that requires a reverse complement function to be called on DNA strings of length 1 through around length 30. Line profiling programs indicate that my functions spend a lot of time getting the reverse complements, so I am looking to optimize.
What is the fastest way to get the reverse complement of a sequence in python? I am posting my skeleton program to test different implementations below with DNA string size 17 as an example.
#!/usr/bin/env python
import random
import timeit

global complement
complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}

DNAlength = 17

#randomly generate 100k bases
int_to_basemap = {1: 'A', 2: 'C', 3: 'G', 4: 'T'}
num_strings = 500000
random.seed(90210)
DNAstrings = ["".join([int_to_basemap[random.randint(1,4)] for i in range(DNAlength)])
              for j in range(num_strings)]
#get an idea of what the DNAstrings look like
print(DNAstrings[0:5])

def reverse_complement_naive(seq):
    this_complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}
    return "".join(this_complement.get(base, base) for base in reversed(seq))

def reverse_complement(seq):
    return "".join(complement.get(base, base) for base in reversed(seq))

tic=timeit.default_timer()
rcs = [reverse_complement_naive(seq) for seq in DNAstrings]
toc=timeit.default_timer()
baseline = toc - tic
namefunc = {"naive implementation": reverse_complement_naive,
            "global dict implementation": reverse_complement}

for function_name in namefunc:
    func = namefunc[function_name]
    tic=timeit.default_timer()
    rcs = [func(seq) for seq in DNAstrings]
    toc=timeit.default_timer()
    walltime = toc-tic
    print("""{}
       {:.5f}s total,
       {:.1f} strings per second
       {:.1f}% increase over baseline""".format(
           function_name,
           walltime,
           num_strings/walltime,
           100- ((walltime/baseline)*100)   ))

By the way, I get output like this. It varies by the call, of course!
naive implementation
       1.83880s total,
       271916.7 strings per second
       -0.7% increase over baseline
global dict implementation
       1.74645s total,
       286294.3 strings per second
       4.3% increase over baseline

Edit: Great answers, everyone! When I get a chance in a day or two I will add all of these to a test file for the final run. When I asked the question, I had not considered whether I would allow for cython or c extensions when selecting the final answer. What do you all think?
Edit 2: Here are the results of the final simulation with everyone's implementations. I am going to accept the highest scoring pure python code with no Cython/C. For my own sake I ended up using user172818's c implementation. If you feel like contributing to this in the future, check out the github page I made for this question.
the runtime of reverse complement implementations.
10000 strings and 250 repetitions
╔══════════════════════════════════════════════════════╗
║ name                       %inc   s total  str per s ║
╠══════════════════════════════════════════════════════╣
║ user172818 seqpy.c        93.7%  0.002344  4266961.4 ║
║ alexreynolds Cython (v2)  93.4%  0.002468  4051583.1 ║
║ alexreynolds Cython (v1)  90.4%  0.003596  2780512.1 ║
║ devonryan string          86.1%  0.005204  1921515.6 ║
║ jackaidley bytes          84.7%  0.005716  1749622.2 ║
║ jackaidley bytesstring    83.0%  0.006352  1574240.6 ║
║ global dict                5.4%  0.035330   283046.7 ║
║ revcomp_translateSO       45.9%  0.020202   494999.4 ║
║ string_replace            37.5%  0.023345   428364.9 ║
║ revcom from SO            28.0%  0.026904   371694.5 ║
║ naive (baseline)           1.5%  0.036804   271711.5 ║
║ lambda from SO           -39.9%  0.052246   191401.3 ║
║ biopython seq then rc    -32.0%  0.049293   202869.7 ║
╚══════════════════════════════════════════════════════╝

You might also get some good advice from CodeReview.SE.

Your implementation of my approach is not doing what I suggested. I suggested working with bytes instead of strings throughout. This means you need your DNAStrings to be in bytes instead of a string and so it would need a separate generation function.

Ah, you meant use them for the entire program. I didn't think to do that. If I were to test that then I would need to convert the entire list of strings to bytestrings before testing, correct?

I don't know if it's the fastest, but the following provides an approximately 10x speed up over your functions:
import string
tab = string.maketrans("ACTG", "TGAC")

def reverse_complement_table(seq):
    return seq.translate(tab)[::-1]

The thing with hashing is that it adds a good bit of overhead for a replacement set this small.
For what it's worth, I added that to your code as "with a translation table" and here is what I got on my workstation:
global dict implementation
       1.37599s total,
       363374.8 strings per second
       3.3% increase over baseline
naive implementation
       1.44126s total,
       346919.4 strings per second
       -1.3% increase over baseline
with a translation table
       0.16780s total,
       2979755.6 strings per second
       88.2% increase over baseline

If you need python 3 rather than python 2, then substitute tab = str.maketrans("ACTG", "TGAC") for tab = string.maketrans("ACTG", "TGAC"), since maketrans is now a static method on the str type.
For those wondering, using biopython is slower for this (~50% slower than the naive implementation), presumably due to the overhead of converting the strings to Seq objects. If one were already reading sequences in using biopython, though, I wouldn't be surprised if the performance was much different.

@Chris_Rands True, updated with the change needed for python3 (thankfully, it's only a single line difference).

Yeh, it''s str.makestrans though, Biopython's implementation is pretty much the same but with some more alphabet and sanity checks https://github.com/biopython/biopython/blob/master/Bio/Seq.py#L860

Here's a Cython approach that might suggest a generic approach to speeding up Python work.
If you're manipulating (ASCII) character strings and performance is a design consideration, then C or Perl are probably preferred options to Python.
In any case, this Cython test uses Python 3.6.3:
$ python --version
Python 3.6.3 :: Anaconda custom (64-bit)

To install Cython, e.g.:
$ conda install -c anaconda cython

The Cython code below seems to offer about the same speed bump as the translation table — perhaps similar code is run under the hood of that.
Two files are needed, starting with setup.py:
from distutils.core import setup
from Cython.Build import cythonize

setup(
  name = 'Reverse complement C test',
  ext_modules = cythonize("revcomp_c.pyx"),
)

And then a second file called revcomp_c.pyx:
from libc.stdlib cimport malloc

cdef int seq_len = 17
cdef char *seq_dest = <char *>malloc(seq_len + 1)
seq_dest[seq_len] = '\0'

def reverse_complement_c_v1(str seq):
    cdef bytes py_bytes = seq.encode('ascii')
    cdef char *seq_src = py_bytes
    cdef int i = 0
    cdef int d = 0
    for i in range(seq_len):
        d = seq_len - i - 1
        if seq_src[i] == 'A':
            seq_dest[d] = 'T'
        elif seq_src[i] == 'G':
            seq_dest[d] = 'C'
        elif seq_src[i] == 'C':
            seq_dest[d] = 'G'
        elif seq_src[i] == 'T':
            seq_dest[d] = 'A'
    return seq_dest[:seq_len].decode('ascii')

This can be compiled into a Python module like so:
$ python setup.py build_ext --inplace

Then we can modify the test bench to include this Cython module and the relevant test method:
#!/usr/bin/env python
import random
import timeit
from revcomp_c import reverse_complement_c_v1

global complement
complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}

DNAlength = 17

#randomly generate 100k bases
int_to_basemap = {1: 'A', 2: 'C', 3: 'G', 4: 'T'}
num_strings = 500000
random.seed(90210)
DNAstrings = ["".join([int_to_basemap[random.randint(1,4)] for i in range(DNAlength)])
              for j in range(num_strings)]
#get an idea of what the DNAstrings look like
print(DNAstrings[0:5])

def reverse_complement_naive(seq):
    this_complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}
    return "".join(this_complement.get(base, base) for base in reversed(seq))

def reverse_complement(seq):
    return "".join(complement.get(base, base) for base in reversed(seq))

tic=timeit.default_timer()
rcs = [reverse_complement_naive(seq) for seq in DNAstrings]
toc=timeit.default_timer()
baseline = toc - tic
namefunc = {"naive implementation": reverse_complement_naive,
            "global dict implementation": reverse_complement,
            "Cython implementation (v1)": reverse_complement_c_v1}

for function_name in namefunc:
    func = namefunc[function_name]
    tic=timeit.default_timer()
    rcs = [func(seq) for seq in DNAstrings]
    toc=timeit.default_timer()
    walltime = toc-tic
    print("""{}
       {:.5f}s total,
       {:.1f} strings per second
       {:.1f}% increase over baseline""".format(
           function_name,
           walltime,
           num_strings/walltime,
           100- ((walltime/baseline)*100)   ))

A sample run:
$ python ./revcomp_test.py
['ACTGCAATAGACCACGC', 'CAGCTTGAGCCATTAAC', 'GGCCCAAGAGTTCGAAC', 'CGACTGTCTCGAATTGT', 'ATCCGCTATGCGCCGTC']
naive implementation
       1.73295s total,
       288525.0 strings per second
       -2.5% increase over baseline
global dict implementation
       1.70896s total,
       292575.5 strings per second
       -1.0% increase over baseline
Cython implementation (v1)
       0.21458s total,
       2330158.7 strings per second
       87.3% increase over baseline

One easy way to speed this up is to use a static const unsigned char array as an ASCII lookup table, which maps a residue directly to its complement. This would replace the nest of if statements and probably give a nice little boost (and it appears it does, making it among the best performers so far!).
Here is my fast implementation of a reverse complement function in C:
https://gist.github.com/alexpreynolds/4f75cab4350e9d937f4a
You might be able to use this directly in Python via the subprocess library. Outsourcing the reverse complement step to a utility written in C will almost always beat the best that Python can do, and you can do nice and important things like bounds checking etc. without losing much speed.
It's unclear how "pure" the answer needs to be, but making a system call from Python seems fair if you're processing strings and your goal is performance.
Another direction to take may be to look at multithreading, if you don't need ordered output. If you have many thousands of sequences stored in memory, you could split an array of sequences up into smaller arrays by use of offsets or array indices. Each thread would work on "rc"-ing sequences in its own piece of the array.

Nice, +1 for Cython

Nice! Can you compare this to the .translate() version in my code? It'd be interesting to see how much Cython is able to speed things up in comparison given the general python overhead.

@Devon_Ryan: With this test bench, the "Cython implementation (v2)" on my Python 3 setup gave a 91.1% increase over baseline and "table" (translate) gave a 84.6% increase. I'm not sure how a Python 2 Cython setup compares.

The most reliable and simplest way is probably using Biopython:
from Bio.Seq import Seq

def my_reverse_complement(seq):
    return Seq(seq).reverse_complement()

print(my_reverse_complement('ATGCGTA'))
# TACGCAT

As Devon has already said here using Biopython isn't as fast as the naive Python solution, and I also tested that shown here with ipython. However, this is because Biopython's implementation, although similar to the naive approach, includes other features; it can reverse complement RNA as well as DNA and it will tell you if you're mixing DNA and RNA
Note that if you really want a fast way you could look at Cython or another python extension. As I edit this now, there are several nice answers taking this approach from user172818 and Alex Reynolds.

@JackAidley I mentioned in my own reply that biopython is ~50% slower than the naive code in the original post. It's good that this one actually included the code for that, though.

Another python extension but without cython. The source code is available at the bottom of this answer or from this gist. On Mac with Python3:
string
   0.38166s total,
   1310067.4 strings per second
   84.4% increase over baseline
seqpy
   0.20524s total,
   2436182.7 strings per second
   91.6% increase over baseline

On Linux with Python2 (seqpy is the first):
seqpy
   0.10960s total,
   4561842.5 strings per second
   93.3% increase over baseline
string
   0.17811s total,
   2807330.9 strings per second
   89.1% increase over baseline

This is file setup.py:
from setuptools import setup, Extension

seqpy_module = Extension('seqpy', sources = ["seqpy.c"])

setup(
    name = "seqpy",
    description = "seqpy",
    ext_modules = [seqpy_module])

This is file seqpy.c:
#include <stdlib.h>
#include <stdint.h>
#include <Python.h>

static unsigned char seq_comp_table[256] = {
      0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,  13,  14,  15,
     16,  17,  18,  19,  20,  21,  22,  23,  24,  25,  26,  27,  28,  29,  30,  31,
     32,  33,  34,  35,  36,  37,  38,  39,  40,  41,  42,  43,  44,  45,  46,  47,
     48,  49,  50,  51,  52,  53,  54,  55,  56,  57,  58,  59,  60,  61,  62,  63,
     64, 'T', 'V', 'G', 'H', 'E', 'F', 'C', 'D', 'I', 'J', 'M', 'L', 'K', 'N', 'O',
    'P', 'Q', 'Y', 'S', 'A', 'A', 'B', 'W', 'X', 'R', 'Z',  91,  92,  93,  94,  95,
     64, 't', 'v', 'g', 'h', 'e', 'f', 'c', 'd', 'i', 'j', 'm', 'l', 'k', 'n', 'o',
    'p', 'q', 'y', 's', 'a', 'a', 'b', 'w', 'x', 'r', 'z', 123, 124, 125, 126, 127,
    128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143,
    144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159,
    160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175,
    176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191,
    192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207,
    208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223,
    224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239,
    240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255
};

static PyObject *seqpy_revcomp(PyObject *self, PyObject *args)
{
    PyObject *r;
    char *seq, *rev;
    int i, len;
    PyArg_ParseTuple(args, "s#", &seq, &len);
    rev = (char*)malloc(len);
    for (i = 0; i < len; ++i)
        rev[len - i - 1] = seq_comp_table[(uint8_t)seq[i]];
    r = Py_BuildValue("s#", rev, len);
    free(rev);
    return r;
}

static PyMethodDef seqpy_methods[] = {
    {"revcomp", seqpy_revcomp, METH_VARARGS, "Reverse complement a DNA sequence"},
    {NULL, NULL, 0, NULL}
};

#if PY_MAJOR_VERSION >= 3
static struct PyModuleDef seqpy_module = { PyModuleDef_HEAD_INIT, "seqpy", NULL, -1, seqpy_methods };
PyMODINIT_FUNC PyInit_seqpy(void) { return PyModule_Create(&seqpy_module); }
#else
PyMODINIT_FUNC initseqpy(void) { Py_InitModule3("seqpy", seqpy_methods, NULL); }
#endif

To compile and test:
python setup.py build_ext -i
python -c 'import seqpy; print(seqpy.revcomp("GGGTT"))'

Use a bytearray instead of a string and then employ maketrans to translate
You do not need the more advanced string encoding capabilities of string to store a string of bases, but you're still paying for it in performance. Devon Ryan's suggestion of maketrans is the huge improvement, 10x faster than your naive implementation. Using the same approach, but swapping everything out for bytes allows a further 40% speed improvement, however:
# DNA sequence generation also needs to be swapped to bytes
tab_b = bytes.maketrans(b"ACTG", b"TGAC")

def reverse_complement_bytes(seq):
    return seq.translate(tab_b)[::-1]

In my comparison:
Table implementation
       0.24339s total,
       2054315.7 strings per second
       82.5% increase over baseline
Bytes table implementation
       0.14985s total,
       3336755.9 strings per second
       89.2% increase over table

If one needs to convert back to string to interface with the rest of the code, what is the impact on speed ?

@bli It is still about 10% faster if you can work with bytes all the way through and then transfer to a string at the end. If you need to go string->bytes->string then it is about 25-30% slower than staying with strings

Here is a revision of my original Cython answer which incorporates my suggestion to use a char lookup array:
from libc.stdlib cimport malloc

cdef int seq_len = 17
cdef char *seq_dest = <char *>malloc(seq_len + 1)
seq_dest[seq_len] = '\0'

cdef char *basemap = [ '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0',   '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0',  'T', '\0',  'G', '\0', '\0', '\0',  'C', '\0', '\0', '\0', '\0', '\0', '\0',  'N', '\0', '\0', '\0', '\0', '\0',  'A',  'A', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0',  't', '\0',  'g', '\0', '\0', '\0',  'c', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0',  'a',  'a' ]

def reverse_complement_c_v2(str seq):
    cdef bytes py_bytes = seq.encode('UTF-8')
    cdef char *seq_src = py_bytes
    cdef int i = 0
    for i in range(seq_len):
        seq_dest[seq_len - i - 1] = basemap[<int>seq_src[i]]
    return seq_dest[:seq_len].decode('UTF-8')

In testing under OS X, Python 3.6.3:
$ python ./revcomp.py
['ACTGCAATAGACCACGC', 'CAGCTTGAGCCATTAAC', 'GGCCCAAGAGTTCGAAC', 'CGACTGTCTCGAATTGT', 'ATCCGCTATGCGCCGTC']
...
Cython implementation (v1)
       0.20664s total,
       2419626.4 strings per second
       88.3% increase over baseline
Cython implementation (v2)
       0.14170s total,
       3528668.5 strings per second
       92.0% increase over baseline

Using my lookup array approach ("v2") adds a very decent performance bump over using if blocks ("v1"), and you can keep everything as a Python string.

Since at least version 1.71 of biopython you can use Bio.Seq.reverse_complement, which also works on plain strings natively (no conversion to Seq objects).  On my mac I get 800k strings converted with that implementation ("biopython just rc") when using the benchmark.
10000 strings and 250 repetitions
                         percent increase over baseline  seconds total strings per second
name
alexreynolds Cython (v2)                          94.0%       0.002491          4015188.7
user172818 seqpy.c                                93.2%       0.002833          3529944.4
alexreynolds Cython (v1)                          90.7%       0.003885          2573863.7
devonryan string                                  86.0%       0.005804          1722983.5
jackaidley bytes                                  85.7%       0.005959          1678085.3
jackaidley bytesstring                            82.1%       0.007459          1340693.2
biopython just rc                                 70.5%       0.012275           814638.6
revcomp_translateSO                               49.9%       0.020822           480254.6
revcom from SO                                    41.7%       0.024253           412320.2
string_replace                                    35.2%       0.026946           371107.3
lambda from SO                                    10.7%       0.037111           269464.0
global dict                                        0.0%       0.041556           240641.7
biopython seq then rc                            -30.9%       0.054430           183721.9
naive (baseline)                                  -0.9%       0.041962           238313.6

def revcomplement(s):
    s = s.upper()
    c = { "A":"T","C":"G","T":"A","G":"C"}
    co = str()
    for i in s:
        co = c[i] + co
    return co

I don't think this piece of code actually "reverts" the sequence but just changes the bases with their complementary bases.

No it reverts the sequence and gives reverse complement.

My bad, co = c[i] + co does the trick. Could be useful if you run a benchmark along with the other offered solutions.

Thanks for joining the community Amardeep. I don't doubt that your code works, but I am a bit sceptical if it answers the original question (seeking for fastest solution). As a matter of fact, your solution is sort of included in the question already (reverse_complement_naive). There are plenty of questions that need to be answered, just make sure that you are really addressing what the person have asked :-)

Sorry.  I implement what you said. Thanks

authors: Костя Дряхлов

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

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


text:
Post ipyrad filtering of SNP loci from GBS, based on SNP quality
I want to use SNPs produced by Ipyrad, which is a python script for RADseq, using maize genome data via RAxML to examine the monophyly of a highly-variable focal species and its phylogenetic relationships with about 40 congeners.
Essentially phylogenomics using maximum likelihood trees.
An initial analysis suggested that samples from two outgroup species might share few SNPs with other samples, so the threshold for shared loci was set to two samples.  As a result, the .u.snps.phy alignment has over 600,000 SNPs.
For subsequent use in RAxML, I would like to implement two filters to construct alignments with subsets of SNPs:

SNPs that occur in at least one outgroup sample and in some (potentially maximal) number of other samples,
SNPs constituting a random sample from an arbitrary subset of samples, i.e. potentially excluding the outgroup samples or other samples,
merge the alignments into one

How can I do this efficiently, maximizing use of existing functionality in R or Python? Either code or a (likely successful) strategy would be helpful. I have access to all output formats from ipyrad as potential starting points.

Sorrey but I was lost at the first sentence, could you add a link to what ipyrad, GBS and RAxML is? Also, could you explain what have you tried ? (It is better to show the effort done before asking random people in internet) You seem to have already a plan, so you might have already coded something...

My bad.  Links are now available. The initial analysis wasn't appropriate and I didn't code it. I have no code yet, or I would have posted it. I do have some idea about what I want the code to accomplish. Since I haven't post-processed an alignment before, I aimed my question towards a non-random someone with knowledge of ways to post-process alignments (and specifically from ipyrad).  I naively tried a tibble (in R), with SNPs in columns, but it hung--must not be an appropriate approach with this much data. If I could work on rows and cols, it would be easy, but another approach must exist.

I don't know the answer but have a look at Bioconductor packages (R) and Biopython (python).

Its easy, just Biopython's AlignIO, all manipulations can be done therein. RADseq is restriction digest genome sequencing, bits of random contigs. I think the question is an example of maxium obstrufication, its quite simple really..

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


text:
Difference between 1D and 1D^2 data
I am working on ONT data. Initially, I have worked on data from 1D using Minimap2 aligner. I came across 1D^2 from ONT website. I was wondering how different is the data from two techniques? What is the difference between 1D and 2D? and do we get 2  fastq files (one per strand) for each sample?
TIA
I would greatly appreciate any help regarding my question.

1D² data that is called by Albacore generates two overlapping sets of files. Here's what Devon O'Rourke mentioned on the ONT community forums and on Twitter:

Following Albacore basecalls on a 1d2 library I get two sets of .fq files, summary stats, etc.: one for the 1D basecalling script, and one for the 1D2 script. I'd like to use the reads generated from these scripts in a genome assembly.
After a little bit of grep searching, it seems like there are overlapping reads among the 1D and 1D2 directories - namely, those 1d reads that generated 1d2 reads!

Link here, for those who have ONT community access.
I've recommended generating a list of the 1D^2 reads, then excluding those reads from the 1D fastq files. I have a couple of scripts that can help for this:
$ cat 1d2_reads.fq | ~/scripts/fastx-length.pl | perl -lane '$F[1] =~ s/(.{36})/$1\n/; print $F[1];' > readNames.txt
$ cat 1d_reads.fq | ~/scripts/fastx-fetch.pl -v -i readNames.txt > filtered_1d_reads.fq

More information on this can also be found in this StackExchange answer.


text:
Link to UCSC gtexGene track description page for selected gene symbol
I have this link that comes up when I am in the browser window and click on the gene name, which leads me to description page about the gene (in this example it is BOK gene).
How can I generate that link? Below are some of my failed attempts. There is a UCSC FAQ page about the links, didn't help much, maybe I am reading it wrong.
Link that works, generated when clicked on gene on the website:
http://genome-euro.ucsc.edu/cgi-bin/hgc?hgsid=226968178_C6XmUjsjtt3Mm9fOmePo6VqRN6TV&c=chr2&l=241667282&r=242920780&o=242498135&t=242513546&g=gtexGene&i=BOK
Which gives me this page:

My failed attempts to generated custom link:
http://genome-euro.ucsc.edu/cgi-bin/hgTracks?db=hg19&position=chr2%3A241667283-242920780&gtexGene&i=BOK
http://genome.ucsc.edu/cgi-bin/org=human&db=hg19&gtexGene&i=BOK
http://genome.ucsc.edu/cgi-bin/gtexGene?org=human&db=hg19&i=BOK
http://genome.ucsc.edu/cgi-bin/gtexGene?org=human&db=hg19&hgg_gene=BOK
http://genome.ucsc.edu/cgi-bin/gtexGene?org=human&db=hg19&i=BOK
http://genome.ucsc.edu/cgi-bin/gtexGene?org=human&db=hg19&g=gtexGene&i=BOK
http://genome.ucsc.edu/cgi-bin/org=human&db=hg19&g=gtexGene&i=BOK
http://genome.ucsc.edu/cgi-bin/hgTracks?org=human&db=hg19&g=gtexGene&i=BOK

The first one of your failed attempts work for me, but the returned page is different to the reference link. Could you clarify what is the expected result?

@Llopis see screenshot, I just again tested all the links, no they don't work for me :(

I can't reproduce with the first and the last link (those with hgTracks). But the page is not like the screenshot. Maybe you (or your IP) hit an API limit ?

@Llopis the first long link that was generated from website works for you?

Yes, the link with "hgsid=226968..." works for me too.

Great, that is the page I want to be able to access by pasting together some bits, programmatically.

Let us continue this discussion in chat.

You never want to include the hgsid parameter in your url because that is the link to your current browsing session and all its specific contents. So if you send that link to someone else and then they start browsing around that link may not work correctly again in the future if you or the person you sent the link to modify your browsing session.
As to how to create links to hgc pages, the usually follow a somewhat general format:
http://genome-euro.ucsc.edu/cgi-bin/hgc?db=hg19&c=chr2&o=242498135&t=242513546&g=gtexGene&i=BOK
You need the db and g parameters to specify the assembly build and track, and then you need to provide the position in the form of the c,o and t parameters. The last thing you want is the i paramter to specify the item name. The reason the i parameter alone is not enough is because we do not require uniqueness on item names (there may be many copies of an item across the genome for instance).
You can obtain the c, o, t and i parameters by selecting the chrom, chromStart, chromEnd and name fields from the gtexGene table:
$ mysql --user=genome --host=genome-mysql.soe.ucsc.edu -A -P 3306 -e "select chrom, chromStart, chromEnd, name from gtexGene limit 1" hg38
+-------+------------+----------+---------+
| chrom | chromStart | chromEnd | name    |
+-------+------------+----------+---------+
| chr1  |      11868 |    13052 | DDX11L1 |
+-------+------------+----------+---------+

Great, thank you. Is the link matching on the positions from c,o,t parameters or by the name i=BOK? In other words can that c,o,t region be wider than the gene size, meaning look in this wider region for gene named BOK?

By the way, what does c, o, t stand for, c=chromosome?

OK, just tested, c,o,t must match with gene positions. :(


text:
How to make chromosome color maps for bed ranges
I have genomic .bed file data of 4 different types; type A,B,C,D. These are some examples-
Type A:
1   101380000   101710000   A
1   110085000   110320000   A

Type B:
1   100930000   101335000   B
1   10430000    10560000    B

Type C:
1   101775000   101825000   C
1   103570000   104070000   C

Type D:
1   153325000   153475000   D
1   156355000   156415000   D

This data is for chr1, but I such data for all the chromosomes. I want to draw a chromosome map (possibly using the Gviz package?) which can assign different colors to the 4 types of genomic ranges and make me an idiogram like plot for the chromosomes I wish to see, so I know where my ranges physically are, on the chromosome.
Can this be done using the Gviz package in R? If not, what is the best way and code to do this?

On my (retina) screen some of these ranges would be less than a pixel in width. Are you sure this is the best approach?

@heathobrien Aah then I guess perhaps not :/ Can you show me how to do it anyway, then I can think about it. My bed files are ~1000-4000 lines in number so maybe I can see something.

@heathobrien on smaller chromosomes perhaps

I haven't done that specifically, but it should be easy with ggbio. If you want to see the distribution of your ranges across a chromosome, I suggest plotting the centres of the ranges as points. There are a bunch of examples here

You could also use karyoploteR too instead of Gviz, see the vignette in Bioconductor to see the resulting plots

Take a look at Circos highlights, zoomed into a subregion of chr1, e.g. http://circos.ca/documentation/tutorials/highlights/wedges/images

You can use karyoploteR to plot your regions on an ideogram quite easily.
Disclosure: I'm one of the authors of the tool.
For this example I'll start creating 4 GRanges with random regions, but with your data you should load your bed into R and split it into four GRanges.
library(karyoploteR)

typeA <- createRandomRegions(nregions = 100, length.mean = 5e6, length.sd = 3e6)
typeB <- createRandomRegions(nregions = 20, length.mean = 15e6, length.sd = 10e6)
typeC <- createRandomRegions(nregions = 1000, length.mean = 5e5, length.sd = 2e5)
typeD <- createRandomRegions(nregions = 100, length.mean = 2e6, length.sd = 5e6)

colors <- c(A="#FFD700", B="#00BFFF", C="#FF7F50", D="#90EE90")

And we can plot them. We have to start creating an empty ideogram with plotKaryotype (we are working with human here, but you can use any genome you want) and then plot the regions in the GRanges with kpPlotRegions. As a bonus, if two regions overlap it will plot them one above the other.
kp <- plotKaryotype(genome = "hg19")

kpPlotRegions(kp, data=typeA, r0=0, r1=0.2, col=colors["A"])
kpPlotRegions(kp, data=typeB, r0=0.25, r1=0.45, col=colors["B"])
kpPlotRegions(kp, data=typeC, r0=0.5, r1=0.7, col=colors["C"])
kpPlotRegions(kp, data=typeD, r0=0.75, r1=0.95, col=colors["D"])
legend(x = "right", legend = names(colors), fill = colors, border=NA)

And you'll get something similar to this

In this example we have fairly large regions and not a lot of them. As @heathobrien pointed out, plotting the actual regions in the whole genome plot might not be the best visualization option. You can either use zoom to create different plots of smaller parts of the genome or change your visualization to show the density of regions along the genome with kpPlotDensity.
Plotting the density and tweaking a bit the visualization (and changing the example to have 10 times more regions) with this codei
typeA <- createRandomRegions(nregions = 1000, length.mean = 5e6, length.sd = 3e6)
typeB <- createRandomRegions(nregions = 200, length.mean = 15e6, length.sd = 10e6)
typeC <- createRandomRegions(nregions = 10000, length.mean = 5e5, length.sd = 2e5)
typeD <- createRandomRegions(nregions = 1000, length.mean = 2e6, length.sd = 5e6)

kp <- plotKaryotype(genome = "hg19", plot.type = 4,
                    ideogram.plotter = NULL, labels.plotter = NULL, main="Density")
kpAddCytobandsAsLine(kp)
kpAddChromosomeNames(kp, srt=45)

kpPlotDensity(kp, data=typeA, r0=0, r1=0.2, col=colors["A"], window.size     = 10e6)
kpAddLabels(kp, labels = "Type A", r0=0, r1=0.2)
kpPlotDensity(kp, data=typeB, r0=0.25, r1=0.45, col=colors["B"], window.size = 10e6)
kpAddLabels(kp, labels = "Type B", r0=0.25, r1=0.45)
kpPlotDensity(kp, data=typeC, r0=0.5, r1=0.7, col=colors["C"], window.size = 10e6)
kpAddLabels(kp, labels = "Type C", r0=0.5, r1=0.7)
kpPlotDensity(kp, data=typeD, r0=0.75, r1=0.95, col=colors["D"], window.size = 10e6)
kpAddLabels(kp, labels = "Type D", r0=0.75, r1=0.95)

it would look something like this

The plots could be tweaked a lot changing various parameters and you could add other data on top of it if needed. You can find more information on the Bioconductor landing page or at the karyoploteR's Tutorial and Examples page.

@bernatgel thank you so much for the elaborate answer! I used karyoploteR as the comments suggested, but your code gave me a better answer, which I can now tweak about. Thanks!

I moved the discussion to chat so as not keep this excellent answer clean and not to bother bernatgel with the discussion.

The downloaded source packages are in
‘/tmp/RtmpKhPald/downloaded_packages’
Warning message:
package ‘karyoploteR’ is not available (for R version 3.4.3)
Can you give me a solution for this? Do I upgrade my R?

@rishi karyoploteR is available from R3.4 onward, so it should be available for 3.4.3. Are you installing it with biocLite?
source("https://bioconductor.org/biocLite.R")
biocLite("karyoploteR")

@bernatgel yes, that's exactly what I did. I even updated my biocLite to be sure

@rishi Can you send me the exact code you are using to install karyoploteR and the output of sessionInfo()? maybe by email would be best: bgel at igtp.cat   and I'll try to understand what's going on

@bernatgel This does not work for me


text:
What are the differences between BNGL (BioNetGen Language) and SBML (System Biology Markup Langague) formats?
Per request from meta comment.
I am self-learning about whole cell modeling, specifically An introduction to whole-cell modeling and Fundamentals of Systems Biology: From Synthetic Circuits to Whole-cell Models. While my background is decades of programming, for cellular biology it is only skimming through a few books.
In the process of learning I soon discovered software/apps that use SBML (System Biology Markup language) but in referring back to "An introduction to whole-cell modeling" see that they use BioNetGen (GitHub) which has its own language BioNetGen language (BNGL)
What are the differences between BNGL (BioNetGen Language) and SBML (System Biology Markup Langague) formats?

Both formats have formal grammars for describing similar concepts, however they are different and to my knowledge one is not a proper subset of the other. While one could do a comparison of a syntactic nature, a simple Google search of just both acronyms reveals BioNetGen 2.2: Advances in Rule-Based Modeling
In here can be found:

2.4 SBML-to-BNGL translation
  SBML is a widely-used model exchange format in systems biology (Hucka et al., 2003). Models encoded in SBML
  are ﬂat, i.e., their species do not have internal structure, which
  limits their utility for rule-based modeling. BioNetGen 2.2 includes
  an SBML-to-BNGL translator, called Atomizer (also available as a web
  tool at ratomizer.appspot.com), that can extract implicit molecular
  structure from ﬂat species (Tapia and Faeder, 2013). A full report on
  Atomizer and its application to the BioModels database (Li et al.,
  2010) is currently in preparation. However, Tapia and Faeder (2013)
  reported that an early version of the tool could recover implicit
  structure for about 60% of species in models within the database that
  contain ≥20 species. Thus, Atomizer makes available a large set of
  pre-existing models in a rule-based format, facilitating their
  visualization (Wenskovitch et al., 2014) and extension (Chylek et al.,
  2015).

While not directly related to the question this also helped in understanding more about BioNetGen
Wikipedia: Multi-state modeling of biomolecules

Biological signaling systems often rely on complexes of biological
  macromolecules that can undergo several functionally significant
  modifications that are mutually compatible. Thus, they can exist in a
  very large number of functionally different states. Modeling such
  multi-state systems poses two problems: The problem of how to describe
  and specify a multi-state system (the "specification problem") and the
  problem of how to use a computer to simulate the progress of the
  system over time (the "computation problem"). To address the
  specification problem, modelers have in recent years moved away from
  explicit specification of all possible states, and towards rule-based
  formalisms that allow for implicit model specification, including the
  κ-calculus, BioNetGen, the Allosteric Network Compiler and others. To
  tackle the computation problem, they have turned to particle-based
  methods that have in many cases proved more computationally efficient
  than population-based methods based on ordinary differential
  equations, partial differential equations, or the Gillespie stochastic
  simulation algorithm. Given current computing technology,
  particle-based methods are sometimes the only possible option.
  Particle-based simulators further fall into two categories:
  Non-spatial simulators such as StochSim, DYNSTOC, RuleMonkey, and
  NFSim and spatial simulators, including Meredys, SRSim and MCell.
  Modelers can thus choose from a variety of tools; the best choice
  depending on the particular problem. Development of faster and more
  powerful methods is ongoing, promising the ability to simulate ever
  more complex signaling processes in the future.

Answering your own question immediately is not frowned upon. you have the option since you ask! But is usually recommended to not accept an answer between minutes of asking it.

It would help tremendously if this were written and formatted as a proper answer rather than a meta-discussion.

@KonradRudolph I edited both Q and A, but you can improve them further :D in case I missed something

I have still to read the sources, but I hope so. Thanks (+1) for sharing your knowledge. One more thing, in case you cite a paper, usually is better to cite the peer-reviewed and accepted paper if possible.


text:
Can I upload only co-assemblies to ENA, using existing ENA samples?
Lets say I am using reads downloaded from ENA, assembling them and doing further downstream analyses. Would it be possible to then upload the results of the co-assemblies and other analyses back to ENA, and just linking them up with the reads I originally downloaded, without having to re-upload the original reads?


text:
Can I use my computer while MinION is sequencing without negatively affecting the run?
I have a Mac Book Pro and I am about to start a sequencing run on the MinION. MinKNOW is going to be running for 2 days.
Can I use my computer during sequencing? Can I browse the web and use Excel, etc? I am not speaking about running a parallelised genome assembly on all my cores using my entire RAM but more mundane activities, like checking my emails or taking notes. Furthermore, can I lock my screen using Control + Shift + OFF without negatively affecting the sequencing?

By locking do you include suspending the machine or not?

Yes, that should be fine, especially if basecalling is disabled. Locking a screen should be fine as well.
I prefer to do basecalling after a run anyway, because there's a chance that it could chew up too much processor time and result in lots of skipped reads.

Since biological samples and reagents cost time and money, I wouldn't push it too much.  Checking your email would be fine, but I wouldn't be browsing the web (besides gmail/calendar) in case a plug-in (e.g. flash), bad javascript, etc. causes any trouble.
I'd disable any auto-updates too.
Locking the screen would be fine.
Personally, I'd buy a cheap chromebook and let the Mac work.


text:
Set directory where MinKNOW writes FAST5 files
I would like to set the folder in which MinKNOW writes the raw data. How can I do this? I do not currently know where MinKNOW will output my data. What is the default directory where MinKNOW outputs the data? My OS is Mac OS X.

Important to mention your operating system. At least for Linux, there is a user.conf file somewhere in /opt/ONT/ (don't remember the exact path, but find will help you find it).

I found a page explaining where this is set up in a user_conf file here:
https://community.nanoporetech.com/support/faq/test1/minknow/minknow/how-do-i-change-the-directory-for-the-reads-folder?search_term=user_conf
Windows
Open Notepad with administrative rights (search for notepad, right click and then select run as administrator). Select open file and navigate to C:\ProgramFiles\OxfordNanopore\MinKNOW\conf and then open user_conf to view the script. Under output_dirs change the two instances of "C:\data" (default read location) to the directory required. Please note the "\" is important.
Mac OS
Go to /Applications/MinKNOW.app/Contents/Resources/conf/user_conf and open user_conf to view the script. Under output_dirs replace the default path "/Library/MinKNOW/data" with your path of choice.
Ubuntu
Go to /opt/ONT/MinKNOW/conf and open user_conf to view the script. Under output_dirs replace the default path "/var/lib/MinKNOW/data/reads" with your path of choice.


text:
Basecalling process carried out by full_1dsq_basecaller.py/.exe
I am unable to understand what exactly the full_1dsq_basecaller.py/.exe script does. Is it giving the consensus sequence for linear and complementary strand, or it is just giving the reads that came from the complimentary strand?
If we are getting two fastq files as output which one we should use for alignment?
Thank you in advance.
I greatly appreciate any help

Possible duplicate of Removing repeated reads from nanopore 1D² reads

My understanding is that the 1D$^2$ process is carried out in three stages:

Carry out a 1D calling of reads (i.e. what the standard Albacore does at the moment)
Link reads that are similar to each other (in reverse-complement orientation), from the same channel, and have a very short pore transition time
Create a 2D consensus base call, using information from the two linked reads

The two outputs from the basecaller are the 1D calls from step 1, and the consensus calls from step 3. Ideally just the 2D consensus calls should be used for assembly (assuming it represents a random distribution of reads), but  if that gives too low coverage (or is not random enough), then the non-converted reads from the 1D calls could be included as well.
Because the linking is not a physical link, it's not a good idea to use 1D$^2$ for amplicon reads; the likelihood of two complementary reads from different templates appearing one after another is too high.


text:
efficient counting of dinucleotides/trinucleotides on fastq reads?
What's an efficient way of counting of dinucleotides/trinucleotide pattern on fastq.gz file reads?
I know there are tools like seqtk that will be very efficient at reading through the .fastq.gz files, and tools like cutadapt that are good at finding long(ish) adapter patterns in the fastq.gz PE of files, but, what would be the best tool to count the number of reads where a di-/trinucleotide is present N-number of times in the read?
E.g.
countingtool --pattern CGC file_R1.fastq.gz

Result: column two is number of times the pattern is seen (exclusive numbers, no overlap)
CGC    0     51341234
CGC    1     13948130
CGC    2      1234344
CGC    3        51344
CGC    4         1343

Are those numbers inclusive? I mean, are the 1343 reads that have 4 CGC included in the 13948130 that have 1? Are the reads with >N included in the count of reads with N?

should CGCGC as one occurrence or two?

number of times the pattern is seen (exclusive numbers, no overlap)

Are reads sequenced from a "stranded" library? If not, you need to consider the reverse complement. I would choose the strand, for each read, that gives the most count.

Well, I'm sure you can find a more efficient approach, but here's a simplistic solution using standard UNIX tools:
zcat file.fastq.gz | grep CGC |
    awk -F'CGC' 'NR%4==2{a[NF-1]++}END{for(i in a){print i,a[i]}}' |
         sort -rnk2

I tested this on a 2.6G fastq file (exome data) and it took:
real    1m18.987s
user    1m41.392s
sys     0m10.684s

On a 66G WGS sequence file, it took:
real    34m6.771s
user    40m16.452s
sys     4m10.176s

The extra grep CGC isn't needed, you can do the match in awk, but it does speed things up a little bit. You can always get the reads with 0 matches by counting the number of reads and subtracting the ones with 1. Alternatively, remove the grep:
zcat file.fq.gz | awk -F'CGC' 'NR%4==2{a[NF-1]++}END{for(i in a){print i,a[i]}}' |
    sort -rnk2

Following the OP's suggestion, I also tested with pigz and, it turns out, pigz can fly!
time unpigz -c exome.fastq.gz | grep CGC |
    awk -F'CGC' 'NR%4==2{a[NF-1]++}END{for(i in a){print i,a[i]}}' |
         sort -rnk2

real    0m44.012s
user    1m19.544s
sys     0m14.900s

And for the 66G WGS:
real    17m29.935s
user    29m23.816s
sys     5m17.856s

Brilliant, sometimes awk is all you need! Would pigz speed-up the process?

@719016 good idea! And yes, it certainly does speed it up, see updated answer.

Check out bioawk https://github.com/lh3/bioawk . This prevents counting CGC located in other lines. I wonder if it is faster than your current solutions?

@conchoecia I doubt it because I am taking advantage of awk's internal splitting by setting the field separator (FS) to CGC and using the % function to ignore irrelevant lines. Using bioawk with -c fastx adds a whole level of needless complexity since that will also look at the name, quality and comment lines as well. It also means I can't just use FS but would need to do something much splower like split(). And, indeed, I just tested on the exome file and stopped it after 7 minutes since it was taking too long already.


text:
Albacore wheel not supported on this platform (Mac OS X El Capitan)
I download the Albacore wheel from Oxford Nanopore.
I activated my conda environment that has Python 3 as the default python version.
I tried to install Albacore with
pip3 install ont_albacore-2.1.10-cp36-cp36m-macosx_10_11_x86_64.whl

but got the following error:
ont_albacore-2.1.10-cp36-cp36m-macosx_10_11_x86_64.whl is not a supported wheel on this platform.
You are using pip version 8.1.2, however version 9.0.1 is available.
You should consider upgrading via the 'pip install --upgrade pip' command.`

I then upgraded my pip with:
pip install --upgrade pip
Requirement already up-to-date: pip in /Users/.../anaconda3/envs/python3env/lib/python3.6/site-packages`

So it seems that my pip was already up to date. I do not understand why it complains about my pip if I am using pip3 anyways. I tried to install Albacore again but I got the same message as before.
Can anyone help me install Albacore on my Mac El Capitan laptop, please?
Update based on Wouter's suggestion:
python -m pip install ont_albacore-2.1.10-cp36-cp36m-macosx_10_11_x86_64.whl
Processing ./ont_albacore-2.1.10-cp36-cp36m-macosx_10_11_x86_64.whl
Requirement already satisfied: python-dateutil in /Users/userName/anaconda3/lib/python3.6/site-packages (from ont-albacore==2.1.10)
Requirement already satisfied: setuptools in /Users/userName/anaconda3/lib/python3.6/site-packages (from ont-albacore==2.1.10)
Requirement already satisfied: h5py in /Users/userName/anaconda3/lib/python3.6/site-packages (from ont-albacore==2.1.10)
Collecting progressbar33 (from ont-albacore==2.1.10)
Requirement already satisfied: numpy>=1.12.1 in /Users/userName/anaconda3/lib/python3.6/site-packages (from ont-albacore==2.1.10)
Collecting ont-fast5-api>=0.4.1 (from ont-albacore==2.1.10)
  Using cached ont_fast5_api-0.4.1-py2.py3-none-any.whl
Requirement already satisfied: six>=1.5 in /Users/userName/anaconda3/lib/python3.6/site-packages (from python-dateutil->ont-albacore==2.1.10)
Installing collected packages: progressbar33, ont-fast5-api, ont-albacore
Successfully installed ont-albacore-2.1.10 ont-fast5-api-0.4.1 progressbar33-2.4

At first sight, it looks like Albacore was successfully installed but if I run: read_fast5_basecaller.py --help
I get:
Traceback (most recent call last):
  File "/Users/userName/anaconda3/bin/read_fast5_basecaller.py", line 69, in <module>
    from albacore.input_utils import list_input_files
  File "/Users/userName/anaconda3/lib/python3.6/site-packages/albacore/input_utils.py", line 7, in <module>
    from albacore import log_utils
  File "/Users/userName/anaconda3/lib/python3.6/site-packages/albacore/log_utils.py", line 5, in <module>
    from albacore.pipeline_core import get_debug_level
ImportError: dlopen(/Users/userName/anaconda3/lib/python3.6/site-packages/albacore/pipeline_core.cpython-36m-darwin.so, 2): Library not loaded: /usr/local/opt/python3/Frameworks/Python.framework/Versions/3.6/Python
  Referenced from: /Users/userName/anaconda3/lib/python3.6/site-packages/albacore/pipeline_core.cpython-36m-darwin.so
  Reason: image not found

Is there a library missing?

Are you using a 32bit kernel by any chance?

Which python version?
python --version
You are using a wheel for python3.6

RELEASE_X86_64 x86_64, so 64 bit kernel

and Python version 3.6 python --version
Python 3.6.4 :: Anaconda, Inc.

could you try python -m pip install ont_albacore-2.1.10-cp36-cp36m-macosx_10_11_x86_64.whl

Not sure about conda as I don't use it, but I can highly recommend pipenv which is the python.org recommended package management system for python.
I have tested the following on my Mac and it worked like a dream
mkdir albacore && cd albacore
mv ../ont_albacore-2.1.10-cp36-cp36m-macosx_10_11_x86_64.whl .

# install pipenv if you don't have it
brew install pipenv  # or pip3 install pipenv

# the following will create a virtual environment and install albacore into it
pipenv install ont_albacore-2.1.10-cp36-cp36m-macosx_10_11_x86_64.whl

# activate the virtual environment
pipenv shell

read_fast5_basecaller.py -h

And that should bring up the albacore help menu and you're ready to go.


text:
What does it mean to over sample?
One of my studies (A) involve sequencing the microbiome. After selecting the variable region and the primers for the targeted region of the 16S sequence. The samples were sent to a platform.
For another study (B) we used the same strategy (same primers, region, machine...) but another platform. However we obtained from this platform a lot more sequences:
In study A the mean sequences per sample are 10k, while in the study B the mean sequences per sample are 80k
I consulted with others and it seems that the study A we under-sampled and in study B we over-sampled.
I can understand that by under-sampling we didn't obtain enough reads to be sure about the microbiome present in the sample. But what does it mean to over-sampling?
That expert I consulted suggested normalize the data to fewer sequences, discarding information we already have. I can understand that the sequencing has errors and can create chimeras but the resulting consensus sequences of the OTUs have found a match with databases of microorganisms. What or how should I take into account that I over-sampled?

As you sequence a library you are sampling from the pool of all possible oligonucleotides found in it. This is a stochastic process, meaning that you are more likely to sequence a piece of DNA that is more prevalent.
As you increase the number of reads that you sequence for a library, you decrease the chances of seeing a new read that you have not yet seen before. If you plot the number of unique reads per total reads sequenced, it will look something like the cumulative density function of the exponential distribution, where you have diminishing returns on unique sequences as you continue to sequence.
If you sequence a library and quickly stop seeing unique reads, we say that this library has low complexity. Therefore, you have "oversampled" it if you sequence too deeply. I am guessing that your metagenomics library actually has a lot fewer species/unique sequences than they anticipated, and that you sequenced deeply anticipating more.
Edit: In response to your comments, the bigger problem in 16S analyses is not sequencing too deeply, it is over-amplifying your samples such that they no longer represent what is actually in your gut microbiome samples. In this case if you used read normalization you would not accurately represent what was in the library that you sequenced since the amplified molecules (16S) are so small that you only have len(16S molecule)-(read len) + 1 possible reads in each direction. You might consider removing PCR duplicates with GATK but this will still have the same effect of not representing what was in the sequencing pool. It is very possible that you have many PCR duplicates after 16S amplification. I generally don't trust amplicon sequencing for quantitive analyses, and don't really trust for low coverage "is-this-species-there-or-not" either, especially when dealing with the 'dark matter' of uncharacterized species. You might be well off just asking "is this species present in high coverage at this time point?"

Thanks for your answer. The problem raise when the same samples in both platforms yield different results (but take into account that I am using OTUs not unique reads). I can upload an image if needed. Anyway, how to take oversampling into account (If it is something to take care of)? I am worried that I can't then compare the results between the studies

BTW, wouldn't that be considered then deep-sequencing?

Yes, this is "deep-sequencing", however I really think that term is overused and overhyped. As you can see in your case, deep sequencing actually causes its own problems with normalization on platforms like Illumina or with PCR-based libraries. :)
Whether or not the sample coverage requires special handling depends on the downstream analysis that you would like to perform. Is this whole-genome shotgun metagenomics, something like 16S amplicons, or something else?

As described in the question it is 16S done in MiSeq. I want to perform integration with RNA-seq data and analyze them on its own via comparing different times, response, localization (it is gut sequencing), sex ...

I don't think there's anything inherently wrong with sequencing to high depth (though naive approaches to dealing with sequencing errors such as singleton removal may fail if the depth is so high that the same errors may be encountered repeatedly by chance). The bigger issue is that they are over sampled relative to the undersampled microbiomes. This does complicate comparisons, particularly for things like OTU richness. The obvious solution is to down-sample the oversampled microbiomes so you have the same number of reads for each. This is probably unnecessary if your analyzing the oversampled microbiomes on their own.

I think like this answer specially because the tool I use to transform the fastq files to OTUs is specially designed to remove singletons and chimeras.


text:
How to calculate the memory usage of storing kmers in RAM
I want to write a program in C++ that stores kmers in a hash or in a trie. How can I calculate how much RAM I would need for each type of data structure? For this application my kmers are strand-specific, so I cannot reduce the size complexity by considering canonical kmers.
For a hash, for example, I know that there are $4^k$ possible kmers of length $k$. I also am assuming that DNA character with an alphabet $\in \{G,A,T,C\}$ can be stored in two bits ($2^2 = 4$). Using this logic the kmer hash should occupy $(4^k) \cdot 2 \cdot k$ bits. Plugging this into Wolfram Alpha ((4^13) * 2 * 13 bits to gigabytes), I calculated the amount of RAM that some small kmer sizes would require to store the entire hash in memory in GB.
k = 13 requires    0.2181 GB
k = 15 requires    4.027  GB
k = 17 requires   73.01  GB
k = 19 requires 1306   GB

This escalated rather quickly, and I think that my calculations must be wrong. Otherwise, how could programs count kmers in RAM?
The kmers that I want to store are sometimes sparse compared to the entire space (for k=13, only 29 million or so out of 67 million possibilities).
Constraints

I do not need to store the number of occurences of a kmer, just its existence
I cannot risk false positives. Bloom filters won't work in this case. :^/

To follow up on Devon Ryan's answer, I thought it would be a little fun to write a Python script that demonstrates using a bit array to maintain a presence/absence table.
Note: I wrote a C++ port that includes a custom bitset implementation that can be sized at runtime. This and the Python script are available on Github: https://github.com/alexpreynolds/kmer-boolean
#!/usr/bin/env python

import sys
import os
import bitarray

# read FASTA
def read_sequences():
    global seqs
    seqs = []
    seq = ""
    for line in sys.stdin:
        if line.startswith('>'):
            if len(seq) > 0:
                seqs.append(seq)
                seq = ""
        else:
            seq += line.strip()
    seqs.append(seq)

# build and initialize bit array
def initialize_bitarray():
    global ba
    ba = bitarray.bitarray(4**k)
    ba.setall(False)
    sys.stderr.write("Memory usage of bitarray.bitarray instance is [%ld] bytes\n" % (ba.buffer_info()[1]))

# process sequences
def process_sequences():
    global observed_kmers
    observed_kmers = {}
    for seq in seqs:
        for i in range(0, len(seq)):
            kmer = seq[i:i+k]
            if len(kmer) == k:
                observed_kmers[kmer] = None
                idx = 0
                for j in range(k-1, -1, -1):
                    idx += 4**(k-j-1) * bm[kmer[j]]
                ba[idx] = True

def test_bitarray():
    test_idx = 0
    for j in range(k-1, -1, -1):
        test_idx += 4**(k-j-1) * bm[test_kmer[j]]
    test_result = ba[test_idx]
    if test_result:
        sys.stdout.write("%s found\n" % (test_kmer))
        sys.exit(os.EX_OK)
    else:
        sys.stdout.write("%s not found\n" % (test_kmer))
        sys.exit(os.EX_DATAERR)

def main():
    global k
    k = int(sys.argv[1])

    global bm
    bm = { 'A' : 0, 'C' : 1, 'T' : 2, 'G' : 3 }

    read_sequences()
    initialize_bitarray()
    process_sequences()

    try:
        global test_kmer
        test_kmer = sys.argv[2]
        if len(test_kmer) == k:
            test_bitarray()
        else:
            raise ValueError("test kmer (%s) should be of length k (%d)" % (test_kmer, k))
    except IndexError as err:
        keys = list(observed_kmers.keys())
        for i in range(0, len(keys)):
            sys.stdout.write("%s found\n" % (keys[i]))

    sys.exit(os.EX_OK)

if __name__== "__main__":
    main()

Note that this doesn't look at canonical kmers, e.g., AG is considered a distinct 2mer from its reverse complement CT.
To use this script, you pipe in your FASTA, specify the k, and an optional kmer that you want to test for presence/absence, e.g.:
$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-bool.py 2 CG
Memory usage of bitarray.bitarray instance is [2] bytes
CG found

Or:
$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-bool.py 3 AAA
Memory usage of bitarray.bitarray instance is [8] bytes
AAA not found

Or if the optional test kmer is left out:
$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-bool.py 5
Memory usage of bitarray.bitarray instance is [128] bytes
CATTC found
ATTCT found
TTCTC found
TCTCG found
CTCGG found
TCGGG found
CGGGA found
GGGAC found
TTATA found
TATAT found
TTTAT found
TTATT found
TATTA found
ATTAG found
TTAGA found
TAGAC found
AGACC found
GACCT found
ACCTC found
CCTCT found

Or for the ~67M kmers in a 13-mer set, for which a roughly 8.4MB bit array is reserved:
$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-bool.py 13
Memory usage of bitarray.bitarray instance is [8388608] bytes
TTTATTAGACCTC found
TTATTAGACCTCT found

I think that my calculations must be wrong. Otherwise, how could programs count kmers in RAM?

Hash table based k-mer counters only keep k-mers seen in data. For $16<k\le32$, you need 64-bit integers to store a k-mer. Given $n$ distinct k-mers, a naive implementation with open addressing hash table would roughly take $2n\cdot 64/8=16n$ bytes. We assume the table is half full – this is where factor 2 comes from. There are various tricks to reduce the memory footprint. You'd better read papers. On a side note, don't use std::unordered_map. It wastes extra memory.
A bitmap described by other answers takes $4^k/8=2^{2k-3}$ bytes, which is a better choice sometimes, but rarely used by k-mer counters as it does not work with long k, either. Trie or full-text indices are overkilling for fixed k-mers and are usually slower.
By the way, this question reminds me of your question on efficient reverse complement. Note that for k-mer counting, you traverse both strands at the same time. You should not reverse complement the whole k-mer.

If you simply need a set of present/absent k-mers and your kmer length is short, then convert the sequence length to an integer (a 13-mer would require 26 bits or 4 bytes) and use that as an index into a bit array. Initialize the array to 0 and use 1 for "k-mer present". If you have ~67 million k-mers then your array is 67 million bits, which is only about 8.5 megs.
This will scale "OK", but eventually balloon to a huge size if you use larger k-mers. At that point you simply have to use other data structures. The reality is that most applications are fairly tolerant to small amounts of error.


text:
How many reads has my sequencing run produced on minion?
I am running a 48 h sequencing protocol on a FLO-MIN107, on MinION, with DNA library-prepped with SQK-RAD004. MinKNOW is installed on my Mac and controlling the MinION. I have found the following histogram in MinKNOW. What do the numbers mean? In particular, how many reads have been produced? I understood that the events are consecutive electrical signals and are then transformed into reads? Do I have 4029 reads? The run has started 24 hours ago. 4029 reads surely is really low, isn't it? What does 4,120,722 mean?

The Platform QC of the flow cell was as follows:

group 1: 480
group 2: 374
group 3: 201
group 4: 59
total = 1114

That's indeed the number of reads and that's quite low. How was your pore occupancy (number of pores sequencing) and your flow cell QC (number of pores good enough for sequencing)? How was your library concentration?
Events is an approximation of nucleotides: it's a guess based on the current signal. As far as I know if you multiply this by 1.8 you get approximately your throughput in nucleotides.
The 4,120,722 is the number of events in the currently selected (mouse hovering?) bin of the histogram.

I added my platform QC in my question. Regarding the library concentration, I did not do one. I only used 400 ng input material based on Qubit Broad Range measurement. In the SQK-RAD004, there is no library concentration measurement. Should I measure the concentration of my library before loading it into the MinION?

Oh I've only done ligation preps. The FC looks quite good. How many pores were actually sequencing?

How do I know how many pores were actually sequencing? Right now (28h after starting), there are only 2 single pores but I remember having hundreds at the beginning.

The "Physical Layout" tells you how many pores are sequencing, you need to look at the number in "Strand". The "pore occupancy" is the ratio between Strand / (Strand + single pore), meaning the total number of "good" pores which you are actually using. This is a parameter which tells you someting about how good your library is.
See for example this screenshot: https://pbs.twimg.com/media/DU32ZAaXkAEojB9.jpg:large
Showing 445 pores actually sequencing out of 445+39 = 92% pores occupied.

very useful to know! So the single pores are the pores available and the strand pores are the single pores actually sequencing? My run is finished now and I don't remember how many 'strand' pores I had but I have the feeling that there were not many... Do you know if these parameters are stored somewhere after the run?

That's right, but no, I don't think those are stored.


text:
Why would someone use a CRAM instead of a BAM?
I had this question from a graduate student yesterday, and I was stuck.
What should I say? Why use a CRAM instead of a BAM?

When is it a good idea to use a CRAM instead of a BAM?
When is it a bad idea?

Whenever you want to save space (this can be a substantial savings). Until quite recently (samtools/htslib 1.7), only CRAM supported long CIGAR strings.
If you need to guarantee that any random obscure downstream program will be able to handle it.

Uptake of CRAM has been pretty slow. Java programs using htsjdk (e.g., picard, IGV and GATK) have only relatively recently added support for CRAM. If you need to use an old version of those for some very odd reason then CRAM may not be supported.
There are a lot of programs written in python that use pysam to open BAM files and these should, theoretically, support CRAM. The issue is that some of the functions may fail and one can't assume that authors will have always written the code needed to handle this. I'll use deepTools as an example, since I'm one of its developers. One of the things about CRAM files is that they (by default) are made such that you require a reference genome in order to construct the sequence field in each alignment. This works fine if you're using a standard genome (htslib, via pysam, can fetch many standard genomes from the web automatically), but if you're not, then you need to specify a fasta file to use for decompression. Every tool, then, needs to add an option for this. With pysam 0.14 and htslib 1.7 this can be circumvented by not decompressing the sequence, but behavior has to be explicitly requested.
Another issue is that many tools will use features from the file index, such as the .mapped accessor, to get the number of mapped reads in a file. CRAM files contain very very little information, so this then fails. Consequently, tool authors need to check for CRAM files and both derive and propagate this information through their functions if it's needed. This can be a time-consuming task (e.g., it took me a couple days to get this implemented in deepTools). Relatedly, samtools idxstats is useless on CRAM files, since there are no statistics stored in the index.
That having been said, it's likely that CRAMs slowly gaining acceptance will eventually make it the standard. It's already a convenient archival format, it's just a matter of time before users can assume that most analysis programs are written to handle it.

@EB2127 If this answer solved your issue, please take a moment and accept it by clicking on the check mark to the left.

@terdon Thank you for the reminder. I normally try to leave the question open a few days for other answers :)


text:
How to fix 5'UTR annotation with RNA-seq data?
I am studying the transcriptome of Arabidopsis. Interestingly, the 5' UTR of the latest annotation is usually too long.
Here is an example (this gene is right to left).

You can see both EST (orange bars) and the RNA-seq reads cover a much smaller region of the annotated 5' UTR. Unfortunately, many 5'UTR of the genes have the same problem. Is there a way to fix that? I was wondering if any assembly tools can help me to get the correct 5'UTR ranges. Thanks!!!

What evidence do you have that the gene stops earlier? This observation could just be due to selective amplification of the 3' end in the RNASeq reads.

I think it starts late (not stop early) as I mentioned 5'UTR, not 3'UTR.

Do you really think standard RNAseq is perfect enough to judge the start of a transcript with single-base resolution? There's a reason that people have developed methods explicitly for that.

It is not only my data. It is also from many people's data. There is no read at 5'UTR in the old EST data, and in >100 set of RNA seq data. There is just nothing there. I do not understand why those upstream sequences are annotated as 5'UTR, and it could mislead many sequence analysis.

@Devon Ryan, Could you please enlighten me the reasoning (why did they give a long 5'UTR with no evidence?). I am asking this because there are lots of genes behaving like this. I am afraid it will cause a problem in motif finding. Thanks!
I updated the image to include the EST data.

To be frank, you need to ask the people who made the annotation.

So is there any tool that can help me to reannotate the 5'UTR based on RNAseq data? I am asking this because, as I checked, most cases RNAseq data are quite consistent with the EST data. Thanks!

@l0110 Could you [edit] your comments in the question. I think I would get a better idea of what do you want.

@Llopis, I cannot edit my comments after 5 minutes. Is there any part unclear? Please let me know. Thanks!

In the case of light-grown seedling it looks like there actually are reads (very few) that are actually the entire length of the 5' UTR. I think what you're after is "what is the absolute 5' transcription start site of my genes?" Right? If that is the case then you're better off doing a few SMRT IsoSEQ or ONT cDNA sequencing flowcells to get the starts and stops. Keep in mind that when you select for polyA you are going to get RNAs that have undergone more processing and may have lost some 5' UTR.

Genes, particularly transcribed genes in Eukaryotes, can have alternate start and stop sites (and splicing).  Alternate transcriptional start sites may be found in different tissues or at different developmental periods due to alternate transcription factor binding sites being used.  The only reliable way to determine the transcriptional start and stop sites is by read coverage.
Genome annotations are generally the product of bioinformatic pipelines and are mostly predictions based upon computational models.  Different gene-finders may call different boundaries; for example, in eukaryotic genomes, an exon may be assigned to the wrong up/down-stream gene.  Arabadopsis, as model organism, has had it's gene boundaries considerably refined/corrected/validated, but even so, there could still exist errors in the annotation or, more likely, the boundaries indicated are the most common start/stop sites (see point 1 above).

In short, your UTRs are whatever your experimental data (reads) show it to be.  Don't worry too much about it not matching up to the standard annotation boundaries.


text:
Negative binomial modeling of RNA-Seq data
A common way to model RNA-seq data is using a negative binomial distribution, where each sample-gene pair is modeled by a different negative binomial distribution with mean $\mu_{ij}$ where $i$ and $j$ are indices for genes and samples, respectively (see article). My question is, given we only have a single observation ($X_{ij}$) from each sample-gene pair, how come we can learn the $\mu_{ij}$ from that single observation? We need multiple observations to estimate a mean, right? Maybe I am missing something critical.

The negative binomial distribution is not fitted for a single gene-sample, but along a sample

“how come we can learn the $\mu_{ij}$ from that single observation?” — we can’t. This fundamentally doesn’t work.

@Llopis thanks, but if they are fitting for each gene separately, where is the sample-specific information coming from such that they can learn a different $\mu_{ij}$ for each pair?

I wasn't sure of my first comments and I'm not with this one but, the  value for each gene can't be estimated for a single sample. But it is estimated from all the samples. See the first section of the methods in the linked paper.

I don’t think this is common. I’ve seen Poisson distributions suggested for modelling discrete read counts and normal distribution for continuous log-transformed data. Negative binomial would be a relatively uncommon approach. With such a wide dynamic range (even for old arrays) it’s safe to assume the central limit theorem: the most common way to simulate log-expression values is with a multivariate-normal distribution.

@TomKelly It is, in fact, extremely common. The two most commonly used packages for differential expression use it (edgeR, DESeq[2]).

Is the model wrong ? I thought the NB has as parameters the mean and dispersion (and not the variance).

If reads were independently sampled from a population with given, fixed fractions of genes, the read counts would follow a multinomial distribution, which can be approximated by the Poisson distribution. However, the Poisson distribution assumes the mean equals the variance, which is usually not satisfied by the data. And for the RNA sequencing data, the variance is usually larger than the mean. So we can assume the data follows a $Y|\lambda\sim Poisson(\lambda), \lambda\sim Gamma(r_0,b_0)$ hierarchy.
$$\begin{align}
P(Y=y)&=\int_{0}^{\infty}P(Y=y|\lambda)f(\lambda)d\lambda\\
&=\int_{0}^{\infty}\frac{e^{-\lambda}\lambda^y}{y!}\frac{b_0^{r_0}}{\Gamma(r_0)}\lambda^{r_0-1}e^{-b_0\lambda}d\lambda\\
&=\frac{b_0^{r_0}}{\Gamma(r_0)y!}\int_{0}^{\infty}e^{-\lambda}\lambda^y\lambda^{r_0-1}e^{-b_0\lambda}d\lambda\\
&=\frac{b_0^{r_0}}{\Gamma(r_0)y!}\int_{0}^{\infty}e^{-\lambda(b_0+1)}\lambda^{y+r_0-1}d\lambda\\
&=\frac{\Gamma(r_0+y)}{\Gamma(r_0)y!}\frac{b_0^{r_0}}{(b_0+1)^{r_0+y}}\int_{0}^{\infty}\frac{1}{\Gamma(r_0+y)}e^{-\lambda(b_0+1)}\lambda^{y+r_0-1}(b_0+1)^{r_0+y}d\lambda\\
&=\frac{\Gamma(r_0+y)}{\Gamma(r_0)y!}\frac{b_0^{r_0}}{(b_0+1)^{r_0+y}}\int_{0}^{\infty}\frac{1}{\Gamma(r_0+y)}e^{-\lambda(b_0+1)}[\lambda(b_0+1)]^{r_0+y}\frac{1}{\lambda}d\lambda\\
&=\frac{\Gamma(r_0+y)}{\Gamma(r_0)y!}\frac{b_0^{r_0}}{(b_0+1)^{r_0+y}}\\
&=\frac{\Gamma(r_0+y)}{\Gamma(r_0)y!}\left(\frac{1}{b_0+1}\right)^y\left(\frac{b_0}{b_0+1}\right)^{r_0}\\
\end{align}$$ which is negative binomial distribution.
Negative binomial distribution depicts the number of failures $k$ before $r$th success of a series of independent Bernoulli experiments, with success probability $p$:
$$f(k;r,p)\equiv Pr(X=k)=\frac{\Gamma(k+r)}{k!\Gamma(r)}p^r(1-p)^k\quad k=0,1,2,\cdots$$
$$\mu=\mathbb E(X)=\frac{r(1-p)}{p}$$
$$\sigma^2=\mathbb V(X)=\frac{r(1-p)}{p^2}$$
Then
$$p=\frac{\mu}{\sigma^2}$$
$$r=\frac{\mu^2}{\sigma^2-\mu}$$
The probability mass function (PMF) is:
$$f(k;r,p)\equiv Pr(X=k)={k+\frac{\mu^2}{\sigma^2-\mu}-1\choose{\frac{\mu^2}{\sigma^2-\mu}-1}}\left(\frac{\sigma^2-\mu}{\sigma^2}\right)^k\left(\frac{\mu}{\sigma^2}\right)^{\frac{\mu^2}{\sigma^2-\mu}}$$
The $\mu$ and $\sigma^2$ of each gene can be estimated from data.
Assume that the number of reads in sample $j$ that are assigned to gene $i$ can be modeled by a negative binomial (NB) distribution,
$$K_{ij}\sim NB(\mu_{ij},\sigma^2_{ij})\tag1$$ which has two parameters, the mean $\mu_{ij}$ and the variance $\sigma^2_{ij}$. The read counts $K_{ij}$ are non-negative integers. $q_{i, \rho(j)}$ is proportional to the expectation value of the true (but unknown) concentration of fragments from gene $i$ under condition $\rho(j)$. The size factor $s_j$ represents the coverage, or sampling depth, of library $j$.

First, the mean parameter $\mu_{ij}$ , that is, the expectation value of the observed counts for gene $i$ in sample $j$, is the product of a condition-dependent per-gene value $q_{i, \rho(j)}$(where $\rho(j)$ is the experimental condition of sample $j$) and a size factor $s_j$,

$$\mu_{ij} = q_{i, \rho(j)}s_j\tag2$$

Second, the variance $\sigma^2_{ij}$ is the sum of a shot noise term and a raw variance term,

$$\sigma^2_{ij}=\mu_{ij}+s_j^2v_{i,\rho(j)}\tag3$$

Third, we assume that the per-gene raw variance parameter $v_{i,\rho(j)}$ is a smooth function of $q_{i, \rho(j)}$,

$$v_{i,ρ(j)}=v_ρ(q_{i, \rho(j)}) \tag4$$
This assumption is needed because the number of replicates is typically too low to get a precise estimate of the variance for gene $i$ from just the data available for this gene. This assumption allows us to pool the data from genes with similar expression strength for the purpose of variance estimation (see article).

The $\mu_{ij}$ in that section of the DESeq manuscript is the expected value of sample $j$ in gene $i$ given its group association (with expected value $\mu_{ij}$). This is why $K_{ij}$ (this is presumably what you meant by $X_{ij}$) is a function of $\mu_{ij}$, rather than $\mu_{ij}$ itself. I think you're just swapping $K_{ij}$ and $\mu_{ij}$ in your head.

"given its group association" -- what is the group here? a gene?

@user5054 An experimental group (e.g., WT or treatment).

Where exactly is this group modeled?

See the context surrounding equation 2 in the methods.

Then you cannot estimate $\mu_{ij}$ if you don't have any information about the grouping of samples?

@user5054 you can’t do differential expression without groups or covariates.

You can however estimate the parameters of a negative binomial model by treating all samples as if they come from a single group.


text:
Why does cutadapt remove low quality bases from the ends of reads only?
I use cutadapt to remove low quality bases from my Illumina reads. The algorithm only removes low quality bases from the end until it reaches a good quality base. If there is a bad quality base beyond that, it is not trimmed. Why? Why doesn't the algorithm remove all low quality bases? It could replace low quality bases in the middle of a read by an N, for example. Or is it that the low quality base is still most probably correct and hence one does not want to lose the base information for mapping?

I am commenting on this part:

The algorithm only removes low quality bases from the end until it reaches a good quality base. If there is a bad quality base beyond that, it is not trimmed.

According to its user guide, cutadap is designed this way: it trims off bases from the 3'-end until it sees a base with quality higher than a threshold. This is not a good algorithm. For example, when we see a quality string 30,30,30,3,3,3,3,3,20,3,3, we would prefer to trim it down to 30,30,30 because a 3,3,3,3,3,20 tail is still not useful.
Some tools such as trimmomatic uses a sliding window to avoid this issue. In my view, a better algorithm is the so-called the "modified Mott algorithm" used by phred 20 years ago. seqtk among others implements this algorithm. Note that the Mott algorithm always trims from both ends, which is sometimes not preferred. BWA implements a variant of this algorithm, trimming from the 3'-end only.
In practice, though, different quality trimming algorithms probably work equally well because it is relatively rare to see a high-quality base in the middle of a low-quality tail. In addition, modern aligners can well handle low-quality tails; assemblers and error correctors can correct through such tails, too. It is usually not necessary to apply quality trimming.

It could replace low quality bases in the middle of a read by an N, for example. Or is it that the low quality base is still most probably correct and hence one does not want to lose the base information for mapping?

With a low-quality base, you might have an error/mismatch, but with an N, you always have a mismatch. Masking low-quality bases to N is not as good.
EDIT: responding to the following comment from OP:

What if the base is low-quality and then it matches the reference by mistake? Wouldn't it be better if it did not match and was penalized? I guess this should be modeled mathematically.

If a read could be mismapped due to a sequencing error, its true location is often in the genome. In this case, a competent mapper will give the alignment a low mapping quality (mapQ). A quality-aware mapper such as novoalign will penalize mapQ even more. In comparison, if you hard mask this sequencing error to "N", you will get a mapQ=0. You can see the difference between the two approaches mostly comes from mapQ.
In modern Illumina data, it is quite often to see a Q8 base in the middle of high-quality bases. >80% of them (in theory) are still correct. My hunch is that hard masking all of them would lead to considerable data loss.

Is it possible to expand on the last part? 'With a low-quality base, you might have an error/mismatch, but with an N, you always have a mismatch. Masking low-quality bases to N is not as good.' What if the base is low-quality and then it matches the reference by mistake? Wouldn't it be better if it did not match and was penalised? I guess this should be modelled mathematically.

If you check the read QC statistics of an Illumina run in e.g. fastQC, you will see that at the end of the read the quality decreases. This is because of exhaustion of chemicals at the end of the run. This is a general trend seen in all runs, therefore you can remove these low quality bases from the end of your run. If you have incidentally a bad quality base in the middle of a read this is not the same general trend but is sporadic. Which is hard to automatically remove with tools such as cutadapt. If you really want an N for each low quality bases within the read, you can make a simple program for that, but why would you need that? Most aligners can handle a few bad quality bases within a read, and for variant calling you'll need many reads (so the incidental error base will be overruled by the good ones).

It's also worth noting that many aligners penalize an N by default, so it's not like replacing super-low quality bases with an N is going to greatly boost alignment scores.

For Illumina (and 454) reads, the quality decreases with read length.  It's not linear and is run/library-dependent.  It has less to do with exhaustion of reagents and more to do with the strands on a spot being out of phase due to incomplete/missed base incorporation during sequencing.
It is common practice to trim off the low quality 3' ends as the first step in QC.  Most analyses and most can deal with a few low errors in the middle of the read and what you do about such sequencing errors may vary depending on what type and use of the data.  But in general, nobody wants the garbage at the end of the reads.
P.S. Even for the same library prep, the fwd/rev runs will have different quality curves.
P.P.S. Sometimes the first couple bases are also poor too.

"in general, nobody wants the garbage at the end of the reads." For RNA-seq, quality trimming may be ok, but for de novo assembly, trimming may reduce N50.

If you don't trim prior to assembly then a DBG assembler will have a lot of extra kmers, blowing up RAM.  It's been my experience that trimming with increase N50, not decrease it.  Although I don't think N50 is a great metric to evaluate assembly quality.  I suppose as an alternative to end-trimming, one could do error correction, e.g. with bfc.

Modern assemblers (e.g. allpaths, sga and spades) often come with an error correction step. They only trim an end when they can't correct it. Only a small fraction of rubbish k-mers can go into the graph.

Thanks, user172818, you're probably right.  But we have found it beneficial to use bfc prior to metaspades.


text:
Can I index a compressed FASTA file using STAR?
I am using STAR to align RNA-seq reads to a reference genome. Before the alignment, I need to generate an index of the reference genome. I use the following code to generate the index successfully:
STAR --runThreadN 8 --runMode genomeGenerate --genomeDir output/index/star --genomeFastaFiles ref.fa --sjdbGTFfile ref.gtf --sjdbOverhang 100
This works fine. However, I would like to keep my reference genome compressed to save disk space. So I am trying the following command:
STAR --runThreadN 8 --runMode genomeGenerate --genomeDir output/index/star --genomeFastaFiles ref.fa.gz --readFilesCommand "gunzip -c" --sjdbGTFfile ref.gtf --sjdbOverhang 100
but I get the following error:
EXITING because of INPUT ERROR: the file format of the genomeFastaFile: ref.fa.gz is not fasta: the first character is '' (31), not '>'.
 Solution: check formatting of the fasta file. Make sure the file is uncompressed (unzipped).
I am using the readFilesCommand successfully with compressed RNA-seq fastq files. Does anybody know if there is a similar option to use compressed references? Is there a workaround using Unix commands (maybe piping?) or do I need to decompress the reference, index it and then compress it again?

I can think of two possible workarounds. The simplest would be to try using zcat instead of gunzip -c, in case that works. Assuming it doesn't, you can use process substitution, if your shell supports it. Something like this:
STAR --runThreadN 8 --runMode genomeGenerate --genomeDir output/index/star \
    --genomeFastaFiles <(zcat ref.fa.gz) --sjdbGTFfile ref.gtf --sjdbOverhang 100

The <() creates a file descriptor and passes that to the command. It is very useful for passing streams to commands that expect files.
If your shell doesn't support this construct, assuming you are on a system that supports named pipes (Linux does), you can do the same thing a bit more manually:

Create a named pipe (FIFO)
mkfifo foo

Set it to read the uncompressed genome
zcat ref.fa.gz > foo &

The & is important since the command above will hang until something starts reading from foo. Alternatively, you can just run this in one terminal and the next command in another.
Read the input from the pipe
STAR --runThreadN 8 --runMode genomeGenerate --genomeDir output/index/star \
    --genomeFastaFiles foo --sjdbGTFfile ref.gtf --sjdbOverhang 100


text:
How to you find read statistics for bam files with terminal mismatch?
I have a bam file aligned to a reference sequence. I am interested in finding the number of reads for read length 18 to 34, where the last base is mis-matched. I got this bam file from bowtie alignment allowing 3 mismatches (-v 0) option, so I am expecting some of the reads to have terminal bases mis-matched. How can we count the number of sequences of a specific length with terminal mismatch in bam file?
**UPDATE:
I have created a bowtie aligned .sam file with this command:
bowtie -S -p 18 /media/owner/ref.fa test1.fastq test1.sam
which I then converted to .bam file like this:
samtools view -bS -o test1.bam test1.sam
From this test1.bam, I was able to get the total number of plus and minus strand of reads aligned to the reference loci using bedtools coverage maping (for positives reads and negative reads) like this:
bedtools coverage -a my_region.bed -b test1.bam -bed -d -s | gzip > test1_region_pos.tsv.gz
Since I used default bowtie option which allows for three mismatches, I want to know how many reads are aligned with terminal mismatches (i.e,last three, last two and last one base mismatched at 3' region of the aligned reads at their loci position).

Have you looked into using the R Bioconductor Rsamtools::scanBam function? It can import a BAM file into R and then you can easily get the CIGAR string and check the CIGAR string for mismatches?

Yes I tried using pileup package, but couldn't figure out how to get cigar string.

Could you provide some examples of lines that should be counted? I'm not sure how mismatches are encoded in the version of Bowtie you're using. As far as I remember, Bowtie doesn't distinguish between = and X in CIGARs (and uses M for both), therefore, mismatches must be detected from the MD tag.
You can use SAMsift for filtering the alignments. The following filtering criterion should keep only reads with  read length from 18 to 34 that have an MD tag with 0 as the last character, which must be preceded by a non-number (this should signalise a mismatch at the last position).
samsift \
  -i test1.bam \
  -f '18 <= len(SEQ) <= 34 and MD[-1]=="0" and not MD[-2].isdigit()' \
  -m nonstop-remove

The obtained alignments can be counted by adding | samtools view | wc -l.
Explanation:

-f ... – filtering criterion
18 <= len(SEQ) <= 34 – the allowed read lengths
MD[-1]=="0" – the last character of the MD tag is 0
not MD[-2].isdigit() – the penultimate character of the MD tag is not a digit
-m nonstop-remove – remove alignments without any MD tag (probably unaligned, they would cause an error).

Thanks for your answer. I have updated the question with how I created BAM file.

Could you explain this a bit please MD[-1]=="0" and not MD[-2].isdigit()?

Does this also consider for + and - strands?

Please, could you show at least one line that should pass the filter? I need it to verify the MD filtering step. Yes, it's possible to filter using +/- via FLAG.

The bowtie -v 1 option allows one mismatch, -v 2, alows 2  mismatches anywhere in the sequences. The default option allows maximum of 3 mismatches.

I mean 1 line from the SAM file.

Did you mean: samtools view -H R1_001.sam        @HD VN:1.0 SO:unsorted @SQ SN:NC_013116 LN:2166 @PG ID:Bowtie VN:1.1.2 CL:"bowtie --wrapper basic-0 -S -p 18 /m

Why am I getting permission denied while using samsift? I did chmod u+x.

/usr/bin/python: can't find '__main__' module in 'samsift' Getting this error now?

If you have problems with installation, could you create an issue on Github, please? It might be related to python2 vs. python3. https://github.com/karel-brinda/samsift/issues