text:
How do you generate read-length vs read-quality plot for long-read sequencing data (e.g., MinION)?
How do you generate read-length vs read-quality plot (heat map with histograms in the margin) for long-read sequencing data from the Oxford Nanopore Technologies (ONT) MinION? The MinKNOW software from ONT provides a plot like this during base calling.
This would also be very helpful for PacBio data.

I also wrote a package to create various plots from Oxford Nanopore sequencing data and alignments: NanoPlot. It can be installed through pip (see also the README on Github). In addition to multiple plots also a limited NanoStats output is created (see also NanoStat). Data can be presented using:

A fastq file (optionally compressed)
A bam file
The sequencing_summary.txt file generated by albacore

Using optional flags you can:

Log transform the read lengths
Use aligned reads rather than sequenced reads
Downsample the reads
Set a maximum read length

I've added an example below, plotting log transformed read length vs average read quality (using a kernel density estimate). More examples can be found in the gallery on my blog.
I welcome all feedback and suggestions!

It's important to always consider read length and quality jointly with high-error read data, and current long-read technologies (e.g., MinION and PacBio) have high error rates. Considering read length and quality jointly will help you determine how successful the run was, how many reads were 'high quality', whether the longer reads are 'real' (or just pore noise), etc.
I've had a recent spike of interest in similar plots and came across a project called pauvre (french for 'poor', play on 'pore') through the Oxford Nanopore Technologies (ONT) community that I think is even better than MinKNOW's base calling plot.  Plus, you can generate these plots from a fastq file when ever you want to, unlike with MinKNOW.
[Note: I am not the original author, but I am now contributing because I liked (and needed) it.]

Pauvre will also report useful statistics:
fastq stats for fastq_runid_bb8b8ddedb22bdd6802b2bfa2b4e424c92c30d28_0.fastq
numReads: 2164829
numBasepairs: 4970615217
meanLen: 2296.077527139557
medianLen: 1495.0
minLen: 5
maxLen: 392031
N50: 3450
L50: 402786

                               Basepairs >= bin by mean PHRED and length
minLen          Q0          Q5         Q10         Q15       Q17.5         Q20      Q21.5  Q25  Q25.5  Q30
     0  4970615217  4970611559  4835461771  3889995868  2900103275  1087779109  165637656  429      0    0
 50000    11531044    11531044      270324      160128       50729       50729          0    0      0    0
100000     6260554     6260554           0           0           0           0          0    0      0    0
150000     3504240     3504240           0           0           0           0          0    0      0    0
200000     2501101     2501101           0           0           0           0          0    0      0    0
250000     1609592     1609592           0           0           0           0          0    0      0    0
300000     1033423     1033423           0           0           0           0          0    0      0    0
350000      392031      392031           0           0           0           0          0    0      0    0

                    Number of reads >= bin by mean Phred+Len
minLen       Q0       Q5      Q10      Q15    Q17.5     Q20  Q21.5  Q25  Q25.5  Q30
     0  2164829  2164605  2083436  1626706  1183812  435687  77341    1      0    0
 50000      109      109        5        3        1       1      0    0      0    0
100000       36       36        0        0        0       0      0    0      0    0
150000       15       15        0        0        0       0      0    0      0    0
200000        9        9        0        0        0       0      0    0      0    0
250000        5        5        0        0        0       0      0    0      0    0
300000        3        3        0        0        0       0      0    0      0    0
350000        1        1        0        0        0       0      0    0      0    0

These plots and stats would be equally useful with PacBio, but that's not super easy (though it is possible) with current raw output from the Sequel sequencer: Which quality score encoding does PacBio use?
Pauvre currently uses Biopython to parse the fastq and matplotlib for the actual plot, and will let you choose the output image format (e.g., .png, .pdf, etc.). You can also choose whether the background is transparent or white (for .png output).
The parser is currently super slow because it's using SeqIO.parse, but we're changing parsers to speed that up. We're also adding some extra features (e.g., choose whether to include y-axes in margin histograms, print some stats directly to the plot for documentation, etc.)
Purple is currently the only color choice (which I personally love), but adding options to change that will be super easy.


text:
What is the standard way to work with a diploid reference genome? Complementary strands?
At the moment, the standard reference genomes (e.g. hg19, hg38) are haploid genomes. We know that the human genome is diploid. Naturally, the latter would be the respectively correct representation of the human genome.
More and more biologists are using emerging technologies into order to capture the diploid nature of genetic information, e.g. phasing SNPs between mother and father chromosomes.
What is the standard way that bioinformaticians generated a standard diploid reference genome?
Actually, reference genomes are not truly haploid (by my understanding). Given that reference genomes are 5'-3' coordinated, in order to create a complementary strand, one would need to take the 3'-5' complement. In order to have a diploid genome, you need two reference genomes and two 3'-5' complements.
More importantly, how have large-scale genomic studies dealt with the fact that the haploid reference genome is a consensus-based "half" of a human genome?

This  BioStar post might clarify bit the confusion about orientation.

When you say "work with", are you only thinking of variant calling?

This post is asking two questions. Perhaps it would be worth splitting the post into two?

For calling small variants, the standard way is to simply call diploid genotypes. You can already do a variety of research with unphased genotypes. You may further phase genotypes with imputation, pedigree or with long reads/linked reads, but not many are doing this because phasing is more difficult, may add cost and may not always give you new insight into your data. For these analyses, we use a haploid genome. For human samples, the vast majority of "large-scale genomic studies" are done this way.
A diploid reference actually doesn't help much with reference-based analysis; it only complicates algorithms. What could help a lot is a population reference, which may be represented by a graph or a compressed full-text index or both. In theory, if you have a comprehensive population reference and a capable mapping algorithm, you may call extra variants that would not be callable with short reads. In practice, however, there are quite a few technical challenges. Handling population references is a research topic. There are no "standards" yet.
If the goal is to assemble a new reference genome from a diploid sample, we almost always prefer to produce a diploid assembly. Unfortunately, I believe there are no "standard" procedures, either. SuperNova from 10x genomics builds the diploid information into a graph. Falcon from PacBio uses "unzip". I don't think they have got widely used and evaluated so far.
PS: saw your edit while writing the above. The fact that the genome only represents one strand does not mean we have to create the complement strand explicitly in analyses. We do most of reverse complement on the fly in algorithms as well as in mind.

reference genomes are not truly haploid

That depends on how the reference is assembled. If you sequence a haploid sample (e.g. bacteria), your assembly will be haploid. If you sequence an inbreed lab strain that is almost homozygous (e.g. mouse and fruit fly), your assembly will be nearly haploid. If you sequence a diploid sample, your assembly is very likely to be a mosaic of the two haplotypes. In case of the human reference genome, it is more complicated. It is largely a mosaic of several humans by stitching ~150kb haplotypes from these samples.

Sorry, but lots of groups phase their samples. That is because phasing of genotypes is a prerequisite for state of the art genotype imputation from a haplotype reference panel.  The paper of one such panel has 3437 citations according to the publisher's website (https://www.nature.com/articles/nature11632)

1000g just had to do phasing. It didn't really have a choice. However, not many projects follow the 1000g design nowadays. I was referring to the new genome sequencing projects.

I am not talking about 1000g. I am talking about the large subset of the 3437 studies that used the 1000g (phase1) reference panel. And then there are the 224 citations of the 2015 paper on 1000g phase3 https://www.nature.com/articles/nature15394, and the 63 citations of the 2016 Haplotype Rereferenc Consortium paper https://www.nature.com/articles/ng.3643. I think it really depends on what you are trying to achieve. For large GWAS, SNP chips + phasing + imputation are still the way to go. That may of course be different in your field.  Full disclosure: I am an author on the latter two papers.

Note the context of my sentence: "for calling small variants". It is just an example that you can work with genotypes directly. I of course know GWAS uses imputation all the time.

At the moment, the standard reference genomes (e.g. hg19, hg38) are haploid genomes. We know that the human genome is diploid. Naturally, the latter would be the respectively correct representation of the human genome.

The premise of the OP's question is false.  The natural reference representation of the human genome is not diploid.
Think of a reference genome as a map, and not as a specific example of a human being's DNA.
Not only is the human genome reference haploid, but it is also a composite genome. This means that the human genome reference sequence is composed of sequences from multiple individuals.  In other words, the human reference does not correspond to any one human sequence.
Any particular read from a DNA sequencer will be a read from a human genome that diverges from the reference genome. So an algorithm that tries to match the read to the reference genome will always need to handle potential discrepancies. Adding a second map against which to match a read would not change that fact.  Therefore, there is little value in providing a second haploid reference genome.
Side note:
There are parts of the human genome that "are too complex to be represented by a single path", and the Genome Reference Consortium provides "alternate loci" for such regions of the genome.
Answer to the OP's second question

More importantly, how have large-scale genomic studies dealt with the fact that the haploid reference genome is a consensus-based "half" of a human genome?

I interpret this question as "How do large genomic studies represent genome diversity that cannot be represented by a haploid reference genome?"
One standard approach that has been very popular in genome-wide association studies is the use of a haplotype reference panel such as for example from the 1000 Genomes Project.  Modern phasing programs exploit the shared ancestry of the samples in a study and reference haplotypes to phase a study's samples. One popular phasing program is Impute2

Although technically composite, 70% of GRCh38 is from RP11, a single male donor (ref: http://genome.cshlp.org/content/27/5/849.full).

To add to all the other great answers, I would mention that the question is somewhat misleading. If the reference genome is for a single individual, then it should be diploid. However, it's a reference for all humans. It should really contain billions of copies to fully account for all the diversity. Since that is not realistic, the reference serves as a simple approximation.
This point was addressed by the recent Korean genome paper:

Human genomes are routinely compared against a universal reference.
  However, this strategy could miss population-specific and personal
  genomic variations, which may be detected more efficiently using an
  ethnically relevant or personal reference. ... Systematic comparison
  of human assemblies shows the importance of assembly quality,
  suggesting the necessity of new technologies to comprehensively map
  ethnic and personal genomic structure variations.

There are some assemblers that produce assembly graphs that attempt to describe all the possible haploid paths within a set of reads. Such an assembly attempts to capture all the diploid variation (and/or population variation) in a sample at the expense of not having full-length chromosomes.
Canu (for example) will produce contigs that are extended as long as consensus is maintained across different reads, but when there is a reliable break in coverage (i.e. an area where chromosomes are heterozygous) then the contigs will be broken up. Canu provides as output a GFA file (assembly graph) that can be used to determine which paths might combine together into a single chromosome.

Since nobody has tackled you second question, I'll give it a go:
The complementary strand to a genome sequence is the  5'-3' reverse complement, not the 3'-5' complement. This sequence is not captured in standard reference genomes, but the information is. Given a sequence, it's trivial to calculate it's reverse complement, and any tool that's designed to work with sequence data takes this feature into account.

source: Stack Exchange

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

sort: votes


text:
How to merge sam files together with adding read groups
I have three sequencing libraries of single individual mapped to a reference using bwa-mem. I would like to merge the three unsorted .sam files I have so, I can call variants and heterozygosity estimates using atlas. Atlas requires one input mapping file (bam) with defined read groups because it estimates the error profiles of different libraries separately.
How can I merge multiple sam files? Preferably avoiding java (Picard tools).
I tried to figure out a solution using samtools 1.3. I sorted individual files using samtools sort, then I used samtools merge -r merged.bam s1.sort.sam s2.sort.sam s3.sort.sam to merge the sorted files. However, the read group didn't make it to the header (and the variant caller I use is complaining about it), also the read group is stupidly the file name.
I tried to define meaningful readgroup names using procedure described on BioStars, but I found that this will just change the header, it does not adjust the names of read groups defined by samtools merge (the file names).
Following this related thread on SeqAnswers, I tried to define the correct header with read groups corresponding to merged file names:
samtools -rh rg.txt merged.bam s1.sort.sam s2.sort.sam s3.sort.sam

where rg.txt is a file with header
@RG s1.sort
@RG s2.sort
@RG s3.sort
...
output of samtools view -H s1.sort

However, the header still had not the read group, I guess because sam header accepts only tagged items to be specified (something like @RG    XY:s1.sort). So, I looked into the merged bam and I found that the tag of RG is Z:. So I tried to just rename header of the merged file using samtools reheader, but then samtools complain about the fact that a tag needs to be of length 2:
Malformed key:value pair at line 123: "@RG    Z:s1.sort"
Segmentation fault (core dumped)

I have opened an issue to report this strange incompatibility of read groups generated by samtools merge with samtools reheader.
I would like solution to :

create a bit more standardized read group names (SM:Sample\tLB:library format)
avoid pointless writing to disk like in sam -> sorted.sam -> merged.bam case (can be probably achieved through "pipes and tees", thanks @bli)

I also know that I can specify RG to bwa, so the sam files will have read groups defined in the first place. But I do not like the idea of remapping three libraries just to create correct formatting of read groups.

I have used https://github.com/ekg/bamaddrg to add read groups quickly to multiple sam files.  And then you can do a samtools merge of the tagged files.

This looks like elegant and faster version of awk. I have to take a look!

Not so elegant but working solution
I found a solution that satisfy several of my conditions, basically I just have to assign read groups to individual mapping files, which can be just added to sorting process and then merge the sorted sam files with read groups.
function sort_and_assign_RG {
    # $1 input file
    # $2 read group id
    # $3 library

    OFILE=$(basename $1 .sam).sort.RG.sam
    HEADER=$(basename $1 .sam).header.sam

    # since I take header from unsorted file, so I need to add this information to header as well as the readgroup
    echo -e "@HD\tVN:1.3\tSO:coordinate" > $HEADER
    # add read group to header
    echo -e "@RG\tID:"$2"\tLB:""$3" >> $HEADER
    # and the rest of the header
    samtools view -H $1 >> $HEADER

    # now sort the input sam, remove header,
    # attach to each alignment read group ID (RGID)
    # and cat header and alignment together
    samtools sort -- $1 | samtools view - | \
      awk -v RGID="$2" '{ printf "%s\tRG:Z:%s\n", $0, RGID; }' | \
      cat $HEADER - > $OFILE

    # remove temp header
    rm $HEADER
}

sort_and_assign_RG s1.sam s1 is180
sort_and_assign_RG s2.sam s2 is350
sort_and_assign_RG s3.sam s3 is550

samtools merge merged.bam s{1,2,3}.sort.RG.sam

It would be nice to avoid creation of three sam files on the way, but at least I get standardized read groups in a merged bam tile.
--- edit ---
I also found where my confusion is coming from. The read group in header is in format
@RG    ID:foo    LB:lib_foo    ...`
@RG    ID:bar    LB:lib_bar    ...`

where ID is read group ID (not sample ID as I thought) and other tags are just specifications of the read group. Then individual mapping have read group assigned in different format
{aliment1}    RG:Z:foo
{aliment2}    RG:Z:bar

where RG is just a tag for read group and Z is just a mark that says that this tag is just a "printable string". Therefore I think that another solution would just to merge these sorted files and then just adding three correct lines to the header.
I have to say, .sam is a unbelievable intuitive format.

You may try to use named pipes instead of files for your sam files: mkfifo $OFILE before the samtools sort. Then you call your sort_and_assign_RG in the background with & and do rm -f s{1,2,3}.sort.RG.sam after the samtools merge. I'm not sure this will really work in parallel, though: I suspect that the named pipes can be written only once the samtools sort are completed.

If the files already have defined read groups inside them, then doing a merge without the -r option should work:
samtools merge merged.bam s1.sort.sam s2.sort.sam s3.sort.sam

If there are multiple input files that share the same read group, then by default they will have random strings appended to make the read groups unique. This can be stopped by using the -c option, as mentioned in man samtools merge:
   -c      When  several  input files contain @RG headers with the same ID, emit only
           one of them (namely, the header line from the first file we find  that  ID
           in) to the merged output file.  Combining these similar headers is usually
           the right thing to do when the files being merged originated from the same
           file.

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:
Running differential expression analyses on count matrices with many zeroes
Note: I also posted this issue (with less context) in the bioconductor support site: https://support.bioconductor.org/p/97424/
I'm working on a snakemake workflow that identifies various small RNA species in C. elegans small RNA-seq libraries. Some are endogenous siRNAs supposedly generated from RNA templates (through RNA-dependant RNA polymerases (RdRP)) that can be variously classified (protein-coding genes, transposons, and other repeat types).
I count such small RNA reads in a set of libraries and try to identify differentially producing sources. For this, I use DESeq2 (that I run using rpy2 from within snakemake).
I'm not sure DESeq2 is always appropriate for this kind of data, but so far, the analyses would at least complete.
However, I recently added new potential types of small RNA source (simple repeats and satellites), and these happen to have low counts. I'm not 100% sure, but I suspect these low counts are the reason for failures during DESeq2 analyses (for debugging purposes, I ran this manually in R):
> dds <- DESeq(dds, betaPrior=T)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.
Error in lfproc(x, y, weights = weights, cens = cens, base = base, geth = geth,  :
  newsplit: out of vertex space
In addition: There were 12 warnings (use warnings() to see them)
> warnings()
Warning messages:
1: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
2: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
3: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
4: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
5: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
6: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
7: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
8: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
9: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
10: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
11: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
12: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight

The count matrix indeed contains a high proportion of zeroes, the rest are mostly ones, but there are some sources that seem a little more significantly producing small RNAs:
> mean(counts_data == 0)
[1] 0.7488889
> mean(counts_data == 1)
[1] 0.1507407
> max(counts_data)
[1] 34

How would you recommend handling such kind of data?
Should I for instance discard rows without enough counts?
If so, what would be a reasonable threshold?
Edits: trying to get a dispersion plot:
The same error occurs when trying to estimate dispersion:
> dds <- estimateSizeFactors(dds)
> dds <- estimateDispersions(dds)
gene-wise dispersion estimates
mean-dispersion relationship
-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.
Error in lfproc(x, y, weights = weights, cens = cens, base = base, geth = geth,  :
  newsplit: out of vertex space
In addition: There were 12 warnings (use warnings() to see them)

This seems to prevent the generation of dispersion plots:
> plotDispEsts(dds)
Error in min(py[py > 0], na.rm = TRUE) :
  invalid 'type' (list) of argument
In addition: Warning message:
In structure(x, class = unique(c("AsIs", oldClass(x)))) :
  Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
  Consider 'structure(list(), *)' instead.

Using fitType="local" also fails:
> dds <- DESeq(dds, betaPrior=T, fitType="local")
using pre-existing size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
Error in lfproc(x, y, weights = weights, cens = cens, base = base, geth = geth,  :
  newsplit: out of vertex space
In addition: There were 12 warnings (use warnings() to see them)

But fitType="mean" works:
> dds <- DESeq(dds, betaPrior=T, fitType="mean")
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

The dispersion plot then looks as follows:

What do these fitType options mean, and are there good reasons why they would have to be changed in order to accomodate for low counts?

I don't think that the issue is the low counts, but rather the number of features without any real variance (the black dots at the bottom).
So what the heck is the dispersion plot and why does one need to fit it anyway?
In a typical RNAseq experiment, one measures many thousands of genes with only a few replicates per biological group. This then leads to issues when performing testing, since your variance measurement for any given gene will be larger or smaller than reality, simply due to the lack of a large number of replicates. Consequently, at least since the limma package came out, most RNAseq packages have tried to fit the mean-dispersion with some sort of line and then used that as a background distribution toward which the variance assigned to each gene should be "shrunken" (i.e., moved toward). Getting this at least approximately correct is important, since if you have few replicates and the final "shrunken variance" is significantly wrong, then your p-values will also be wrong.
What are the various fitType options and what are they doing?
The three options that one can use to fit this relationship are parametric, local, and mean. parametric is the default and basically results in fitting to a bent line. This usually works quite well, but you can see that you have two clouds of data, one shaped like a bent line (at the top) and another (the black dots at the bottom) with no variance. That's causing this fit type to fail. Likewise, local is doing a local fit and failing for the same reason. mean is just taking the mean (average) of all of the estimates. This is fool-proof, but you're then regressing to something that really doesn't fit your data. mean is really a last choice if nothing else works.
For your particular dataset, my main concerns would be figuring out (1) why the maximum mean counts are around 20 (that's really low) and (2) why so many of these features have incredibly low variance.

Remove low-count features in advance. This is standard for most tools including DESeq2 and edgeR (see section 2.6).
This will keep you from testing a lot of features that cannot be differentially expressed (using NB model you need at least four reads - below that the uncertainty is too high to ever be DE) and you will make the false discovery correction unnecessarily hard.


text:
Tool for predicting interactions in the cell
What tools are available to predict, based on the structure of a certain protein, its interactions within a cell?
For example, I am considering the split GFP protein and I am trying to predict if there will be any interactions which might hinder its binding efficiency to the other half of the GFP protein.
Note that I am not only looking for protein interactions, but other molecules as well.

You could try one of these tools to predict protein-protein interactions:

Struct2Net

Given two protein sequences, the structure-based interaction prediction technique threads these two sequences to all the protein complexes in the PDB and then chooses the best potential match. Based on this match, the method generates alignment scores, z-scores, and an interfacial energy for the sequence pair. Logistic regression is then used to evaluate whether a set of scores corresponds to an interaction or not. The algorithm is also extended to find all potential partners given a single protein sequence. Further details about the method are described here.

PRED

Pred_PPI is a web-based system that serves for predicting PPIs from different organisms. This server is freely available to any researcher wishing to use it for non-commercial purposes. Based on auto covariance (AC) and support vector machine (SVM), this tool is capable of predicting PPIs for any target protein pair only using their primary sequences, and assigning an interaction probability to each SVM prediction as well. So the user can use this tool to predict novel PPIs with high confidence.

For protein-RNA interactions, you can try:

catRapid

Through the calculation of secondary structure, hydrogen bonding and van der Waals contributions, catRAPID is able predict protein-RNA interaction propensities with great accuracy (up to 89% on the ncRNA-protein interaction database, NPinter).

See updated question, the OP ist not only looking for Protein-Protein interactions...

@Llopis ah, damn. That will be harder.

Thank you for the original answer this does help though

@jaslibra catRapid should be good too. I haven't used it but have used other tools by the same author and they tend to be very good. What else do you need though? Protein-protein and protein-RNA, and what else?

IBIS reports protein interaction with other biomolecules.  Your protein or homolog must exist in the Protein Data Bank.
IBIS @ NCBI: https://www.ncbi.nlm.nih.gov/Structure/ibis/ibis.cgi
"For a given protein sequence or structure query, IBIS reports protein-protein, protein-small molecule, protein nucleic acids and protein-ion interactions observed in experimentally-determined structural biological assemblies. IBIS also infers/predicts interacting partners and binding sites by homology, by inspecting the protein complexes formed by close homologs of a given query. To ensure biological relevance of inferred binding sites, the IBIS algorithm clusters binding sites formed by homologs based on binding site sequence and structure conservation."


text:
How do I generate a color-coded tanglegram?
I want to compare two phylogenies and colour the association lines based on some metadata I have. I have been using ape cophyloplot but I have not had any success in getting the lines to colour accurately according to my data (see previous question).
Note that in my actual work flow I define the colour scheme using a palette to control the colour outcome.
I want a means to make a tangle using phylogenies which I can format. Preferably in R. I like to get an output like this:

Thanks everyone for the suggestions @b.nota answer is useful to colour according to clade groups but does not address my actual question on using metadata to colour the lines. For this answer see below:
dendextend states for color_lines
a vector of colors for the lines connected the labels. If the colors are shorter than the number of labels, they are recycled (and a warning is issued). The colors in the vector are applied on the lines from the bottom up.
# Extract labels from dendrogram on the left
labels <- dendA %>% set("labels_to_char") %>% labels

#Using a metadata table with colours create a vector of colours
labels <- as.data.frame(labels)
labels2 <- merge(labels, metadata, by.x="labels", by.y="Sample.name", sort=F)
cols <- as.character(labels2$Colours)

# Make tanglegram
tanglegram(dendA, dendC, color_lines = cols)

Side note if anyone with a phylogeny needs to convert it to a dendrogram to try this on their data see my other question.

Hi AudileF, the vector of colors for color_lines takes the order of the labels in the tree. If you want to have it follow the labels in another order, you'll need to put the colors in the initial dataframe with your metadata first, and direct color_lines to the right column of the data frame. It is hard to explain without any example.

Hi b.nota I have a workable example in my previous thread linked above, if that is of any use? Just swap rtree for rcoal, like you mentioned earlier.

I think you can try dendextend, in this manual there is an example of coloring the branches. I don't think it is exactly like your coloring, but with a little tweaking you might get your colorscheme in there.
The manual mentions an argument called color_lines for the function  tanglegram():
# The `which` parameter allows us to pick the elements in the list to compare
iris_dendlist %>% dendlist(which = c(1,4)) %>% ladderize %>%
   # untangle(method = "step1side", k_seq = 3:20) %>%
   set("rank_branches") %>%
   tanglegram(common_subtrees_color_branches = TRUE)

According to the manual, the code above produces an image like this:

Try to see if you can fit your color palette in there.

To get your colors in the original order, you can make a dataframe with the labels in the right order with your metadata. You can use merge for this.
So from your example code it would be something like this:
site <- structure(list(name = structure(c(1L, 3L, 4L, 5L, 6L, 7L, 8L,9L, 10L, 2L), .Label = c("t1", "t10", "t2", "t3", "t4", "t5","t6", "t7", "t8", "t9"), class = "factor"), site = c(1L, 1L,1L, 2L, 2L, 3L, 1L, 3L, 2L, 2L)), .Names = c("name", "site"), row.names = c(NA,10L), class = "data.frame")

library(ape)
library(dendextend)

t1 <- rcoal(10)
t2 <- rcoal(10)

str(site)
# The name are factors now, so make chr first
site$name <- as.character(site$name)

# make a data.frame of your labels
labels_df <- data.frame(t1$tip.label)

#merge the 2 data.frames in the right order (hence sort=F)
colors <- merge(labels_df,site,by.x="t1.tip.label", by.y="name",all.x=T, all.y=F,sort=F)

tanglegram(t1, t2, color_lines=colors$site)

Thanks b.nota. this is something akin to a workflow I had before its good but still inaccurate to get tip labels in correct order. See my answer from that thread.  I did try this just in case and it still doesnt fix my issue. Thanks though.

As an alternative (though I realise this doesn't really answer the OPs question directly), Dendroscope and Treemap can also produce these types of plots. I don't know if they'll do coloured connections off the top of my head, but I think they can differentially colour the tree branches.
Just in case that is of use!


text:
How to convert a phylogeny to a dendrogram in R
I want to convert my phylogeny into a dendrogram so I can use it with dendextend in R to produce a tanglegram. I have made some progress but I keep encountering errors, see below:
library(ape)
library(dendextend)

Tree <- rtree(10, rooted=F)

Tree <- read.tree(file="clipboard", text=NULL)
# test <- as.dendrogram(Tree)
## Error in ape::as.hclust.phylo(object) : the tree is not ultrametric

Tree.ultra <- chronos(Tree)
# test <- as.dendrogram(Tree.ultra)
## Error in ape::as.hclust.phylo(object) : the tree is not rooted

Tree.ultra$root.edge <- 0
# test <- as.dendrogram(Tree.ultra)
## Error in ape::as.hclust.phylo(object) : the tree is not binary

This thread seems to be the closest to answering my issues but the workflow is different and Im not sure how to address this binary issue. I have tried using as.hclust, as.hclust.phylo and hclust.

Please [edit] your question and give us a minimal reproducible example. Something that lets us test your code (and ours) and make sure our answers would work for you. A simple version of the tree you are using, for example, or even a way to generate random data that are similar to what you're working with.

Its not really clear to me (and maybe it is just me) what you mean? A dendrogram and a phylogeny are synonyms. Unless theres some suble difference I'm not aware of, so I'll let others enlighten me if thats the case.

Ive added a small example which incurs the same errors. @JoeHealey. Depends whom you ask about dendrogram vs phylogeny. Many say they are synonyms but classically a dendrogram is seen as a simple graphical representation of hierarchy vs phylogeny which is based on more sound and substantial, evidence and data.

If dendrogram is the same as cladogram than it contains just a tree topology, phylogeny is more general term, it frequently contains branch lengths corresponding to mutation rate * time. See this post.

I was googling and it seems that dendrogram is a term used for visualisation of hierarchical clustering, which could be analogous to NJ-tree, but then it is not just a topology. If you mean cladogram, can you use this word?

After searching high and low I have found an answer from this thread
Workflow goes like:
library(DECIPHER)
dend1 <- ReadDendrogram(file="clipboard")

# Alternatively to load from existing environment
Tree$node.label <- NULL # Need to remove node labels
dend1 <- ReadDendrogram(textConnection(write.tree(Tree)))

This object can then be used in the dendextend package and tanglegram function.

Does that give different results than the chronos function?

Contrary to what you say, you can read a character string from within R by passing a textConnection to ReadDendrogram.

Thanks @KonradRudolph I was not aware of that. Im only a novice to be honest. Ill change my answer above to include this.

Hi @KonradRudolph I have been trying the text connect as so: dend1 <- ReadDendrogram(textConnection(Tree)) but I get an error. Error in textConnection(Tree) : invalid 'text' argument

@AudileF Did you read the documentation? What’s Tree? A representation of the dendrogram in text format?

@AudileF not sure if you've figured it out already but what you should have done is dend1 <- ReadDendrogram(textConnection(write.tree(Tree)))

I hadn't actually @august thanks a mill I'll include it above now.

For your example to work you should replace the rtree function with the rcoal function, see example below.
library(ape)
library(dendextend)
Tree <- rtree(10, rooted=F)
is.ultrametric(Tree)
[1] FALSE
is.binary.tree(Tree)
[1] TRUE

Tree_rcoal <- rcoal(10, rooted=F)
is.ultrametric(Tree_rcoal)
[1] TRUE
is.binary.tree(Tree_rcoal)
[1] TRUE

Tree_rcoal2 <- rcoal(10, rooted=F)

tanglegram(Tree_rcoal, Tree_rcoal2)

To make use of a file in newick format, you can try to convert with the chronos function.
s <- "owls(((Strix_aluco:4.2,Asio_otus:4.2):3.1,Athene_noctua:7.3):6.3,Tyto_alba:13.5);"
cat(s, file = "ex.tre", sep = "\n")
tree.owls <- read.tree("ex.tre")
is.ultrametric(tree.owls)
[1] FALSE
dendrogram <- chronos(tree.owls)

Setting initial dates...
Fitting in progress... get a first set of estimates
         Penalised log-lik = -24.04992
Optimising rates... dates... -24.04992
Optimising rates... dates... -24.04908

Done.
is.ultrametric(dendrogram)
[1] TRUE

Hi b.nota. Thats all well and good but I have an actual newick file I need convert to a dendrogram. rtree was to give an example for troubleshooting. Thanks though.

Yes, I understand. But you didn't mention any newick file in your original post. Check [link] https://stackoverflow.com/questions/7445684/how-to-convert-a-tree-to-a-dendrogram-in-r this post from years ago, I think the chronos function might be helpful in your case.


text:
How can I dock a protein to a nucleic acid?
I have a protein of interest and I would like to now how it interacts with RNA. I have structures of both molecules.
What tool can I use?

You mean a specific RNA molecule, right? Not any RNA in general?

Autodock Vina is a common tool for ligand binding studies, if the stretch of RNA isn't too large I suspect it may work with Vina.

I'm not sure what you meant but you can take a look at NPDock (disclaimer we wrote that tool). If you have a structure of your protein of interest, you can dock it to the structure of your DNA/RNA of interest. Mind that this is a rigid body docking which means that the structure will not change upon binding.


text:
What are the labels in my rstats distance matrix?
I have a distance matrix generated by hierfstat, thusly, with a link to the fasta file here:
library(adegenet)
library(hierfstat)
snps <- fasta2DNAbin('test.fa', chunkSize = 50)
gi <- DNAbin2genind(snps)
# manually define the populations for passing to pairwise.fst
p <- c('Botswana', 'Botswana', 'Botswana', 'Botswana', 'France', 'France', 'Vietnam', 'Vietnam', 'Uganda', 'Uganda', 'Uganda', 'Uganda', 'Vietnam', 'Vietnam', 'Laos', 'Laos', 'Laos', 'Vietnam', 'Vietnam', 'Vietnam', 'Vietnam', 'Vietnam')
f <- pairwise.fst(gi, p, res.type = c('dist', 'matrix'))
as.matrix(f)

which produces
    1   2   3   4   5
1   0   0.2189008   0.225567409 0.1821518   0.259409722
2   0.2189008   0   0.130736953 0.1648034   0.191050772
3   0.2255674   0.130737    0   0.1669077   0.006396789
4   0.1821518   0.1648034   0.166907683 0   0.203931457
5   0.2594097   0.1910508   0.006396789 0.2039315   0

Is it safe to assume that 1 = Botswana, 2 = France, 3 = Vietnam, 4 = Uganda and 5 = Laos? i.e. that the distance matrix rows/columns follow the order that each population first appeared in p?
Is there a way to determine for certain, which row/column in the distance matrix corresponds to which population?

Welcome to Bioinformatics! You can change the order of your columns in the input and see how is the output, if the columns in the output follow the columns in the input, then  it would be safer. But this is from adegenet package, isn't it? Could you provide some test data to try myself? Otherwise you could check the source code of pairwise.fst

added a test file and the libraries I was using :-) i guess i could go into the source code or permute the samples to check, but seems sensible that there would be a option just to print them somehow...

Disclaimer: I don't have adagenet and hierfstat installed, I just looked at the source code. It seems like pairwise.fst (the deprecated code from the adagenet package, which you seem to be using) already returns a matrix. Set res.type="matrix" and try calling rownames(f) and colnames(f) afterwards.

Worked a treat, thanks! I'm using pairwise.fst because it is the method mentioned in this tutorial from 2015 http://adegenet.r-forge.r-project.org/files/tutorial-basics.pdf. Where did you read that it is deprecated?

I was searching for the source code. There are pairwise.fst functions in adagenet and in hierfstat. If you look at the function in the adagenet repo, you will see that it is commented out (R/fstat.R). It is, to my (limited) understanding, one of the functions that have been outsourced/exported to the hierfstat package.


text:
Is there public RESTful api for Gnomad?
I currently find Harvard's RESTful API for ExAC extremely useful and I was hoping that a similar resource is available for Gnomad?
Does anyone know of a public access API for Gnomad or possibly any plans to integrate Gnomad into the Harvard API?

Just to add a comment that Cellbase annotation now includes Gnomad data (Exomes and Genomes) and can be accessed via RESTful API -
http://bioinfo.hpc.cam.ac.uk/cellbase/webservices

As far as I know, no but the vcf.gz files are behind a http server that supports Byte-Range, so you can use tabix or any related API:
$ tabix "https://storage.googleapis.com/gnomad-public/release-170228/vcf/exomes/gnomad.exomes.r2.0.1.sites.vcf.gz" "22:17265182-17265182"
22  17265182    .   A   T   762.04  PASS    AC=1;AF=4.78057e-06;AN=209180;BaseQRankSum=-4.59400e+00;ClippingRankSum=2.18000e+00;DP=4906893;FS=1.00270e+01;InbreedingCoeff=4.40000e-03;MQ=3.15200e+01;MQRankSum=1.40000e+00;QD=1.31400e+01;ReadPosRankSum=2.23000e-01;SOR=9.90000e-02;VQSLOD=-5.12800e+00;VQSR_culprit=MQ;GQ_HIST_ALT=0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|1;DP_HIST_ALT=0|0|0|0|0|0|0|0|0|0|0|1|0|0|0|0|0|0|0|0;AB_HIST_ALT=0|0|0|0|0|0|0|0|0|0|0|0|1|0|0|0|0|0|0|0;GQ_HIST_ALL=1591|589|120|301|650|589|1854|2745|1815|4297|5061|2921|10164|1008|6489|1560|7017|457|6143|52950;DP_HIST_ALL=2249|1418|6081|11707|16538|9514|28624|23829|7391|853|95|19|1|0|0|1|0|1|0|0;AB_HIST_ALL=0|0|0|0|0|0|0|0|0|0|0|0|1|0|0|0|0|0|0|0;AC_AFR=0;AC_AMR=0;AC_ASJ=0;AC_EAS=0;AC_FIN=1;AC_NFE=0;AC_OTH=0;AC_SAS=0;AC_Male=1;AC_Female=0;AN_AFR=11994;AN_AMR=31324;AN_ASJ=7806;AN_EAS=13112;AN_FIN=20076;AN_NFE=94516;AN_OTH=4656;AN_SAS=25696;AN_Male=114366;AN_Female=94814;AF_AFR=0.00000e+00;AF_AMR=0.00000e+00;AF_ASJ=0.00000e+00;AF_EAS=0.00000e+00;AF_FIN=4.98107e-05;AF_NFE=0.00000e+00;AF_OTH=0.00000e+00;AF_SAS=0.00000e+00;AF_Male=8.74386e-06;AF_Female=0.00000e+00;GC_AFR=5997,0,0;GC_AMR=15662,0,0;GC_ASJ=3903,0,0;GC_EAS=6556,0,0;GC_FIN=10037,1,0;GC_NFE=47258,0,0;GC_OTH=2328,0,0;GC_SAS=12848,0,0;GC_Male=57182,1,0;GC_Female=47407,0,0;AC_raw=1;AN_raw=216642;AF_raw=4.61591e-06;GC_raw=108320,1,0;GC=104589,1,0;Hom_AFR=0;Hom_AMR=0;Hom_ASJ=0;Hom_EAS=0;Hom_FIN=0;Hom_NFE=0;Hom_OTH=0;Hom_SAS=0;Hom_Male=0;Hom_Female=0;Hom_raw=0;Hom=0;POPMAX=FIN;AC_POPMAX=1;AN_POPMAX=20076;AF_POPMAX=4.98107e-05;DP_MEDIAN=58;DREF_MEDIAN=5.01187e-84;GQ_MEDIAN=99;AB_MEDIAN=6.03448e-01;AS_RF=9.18451e-01;AS_FilterStatus=PASS;CSQ=T|missense_variant|MODERATE|XKR3|ENSG00000172967|Transcript|ENST00000331428|protein_coding|4/4||ENST00000331428.5:c.707T>A|ENSP00000331704.5:p.Phe236Tyr|810|707|236|F/Y|tTc/tAc||1||-1||SNV|1|HGNC|28778|YES|||CCDS42975.1|ENSP00000331704|Q5GH77||UPI000013EFAE||deleterious(0)|benign(0.055)|hmmpanther:PTHR14297&hmmpanther:PTHR14297:SF7&Pfam_domain:PF09815||||||||||||||||||||||||||||||,T|regulatory_region_variant|MODIFIER|||RegulatoryFeature|ENSR00000672806|TF_binding_site|||||||||||1||||SNV|1||||||||||||||||||||||||||||||||||||||||||||,T|regulatory_region_variant|MODIFIER|||RegulatoryFeature|ENSR00001729562|CTCF_binding_site|||||||||||1||||SNV|1||||||||||||||||||||||||||||||||||||||||||||

UPDATE: 2019: the current server for gnomad doesn't support Byte-Range requests.

As of 15 Jan 2020, cellbase (bioinfo.hpc.cam.ac.uk/cellbase/webservices) does. See Pasted's comment on the question. Byte-range requests are supported, and results come from GnomAD, ExAC, and other data sets. A pretty awesome resource actually, even if byte-range requests were supported at GnomAD or UCSC.

The new gnomAD site (as of August 2019) says no, no API yet:

How do I query a batch of variants? Do you have an API?
We currently do not have a way to submit batch queries on the browser, but we are actively working on developing an API for ExAC/gnomAD. If you would like to learn about GraphQL, which we will use to work with the API, an overview can be found at https://graphql.org.  You can also obtain information on all variants from the VCFs and Hail Tables available on our downloads page.

But, the web interface itself already makes POST requests to https://gnomad.broadinstitute.org/api to send and receive JSON/GraphQL.  So, you can make those same queries programmatically right now, even if it's not officially a public API.
Here's an example in Python to get some basic info on variants for a particular gene.  This way you get simple nested Python objects to work with:
  { 'consequence': 'intron_variant',
    'pos': 6928442,
    'rsid': 'rs782435448',
    'variant_id': '12-6928442-C-A'},
  { 'consequence': 'splice_region_variant',
    'pos': 6928462,
    'rsid': None,
    'variant_id': '12-6928462-C-A'},
  { 'consequence': 'splice_acceptor_variant',
    'pos': 6928464,
    'rsid': 'rs782577109',
    'variant_id': '12-6928464-G-A'},
  { 'consequence': 'missense_variant',
    'pos': 6928466,
    'rsid': 'rs782208003',
    'variant_id': '12-6928466-C-T'},

(I found it useful to go this route because then the full metadata visible in the gnomAD web interface is then available, including the per-variant details like allele counts by population.  I couldn't find this information in the other APIs described here.)

can you expand on your example a bit. the graphAPI is unfamiliar to me. How would you turn the json above into an actual query?

You can browse gnomAD variants with ClinGen Allele Registry (there is API spec available).

I faced same issue recently, I found those link and python script:
gnomAD GraphQL api
https://gnomad.broadinstitute.org/api
It works great but it is a kind of different query language.
Please check here for the docs: https://graphql.org/learn/queries/
gnomAD Python Api
https://github.com/furkanmtorun/gnomad_python_api

I found Jesse's code quite usefull ! For those who try to reproduce it, you should now add the reference genome ID, such as :
#!/usr/bin/env python

import requests
# import pprint
# prettyprint = pprint.PrettyPrinter(indent=2).pprint

def fetch(jsondata, url="https://gnomad.broadinstitute.org/api"):
    # The server gives a generic error message if the content type isn't
    # explicitly set
    headers = {"Content-Type": "application/json"}
    response = requests.post(url, json=jsondata, headers=headers)
    json = response.json()
    if "errors" in json:
        raise Exception(str(json["errors"]))
    return json

def get_variant_list(gene_id, dataset="gnomad_r2_1"):
    # Note that this is GraphQL, not JSON.
    fmt_graphql = """
    {
        gene(gene_id: "%s", reference_genome: GRCh38) {
          variants(dataset: %s) {
            consequence
            pos
            rsid
            variant_id: variantId
          }
        }
      }
    """
    # This part will be JSON encoded, but with the GraphQL part left as a
    # glob of text.
    req_variantlist = {
        "query": fmt_graphql % (gene_id, dataset),
        "variables": {"withFriends": False}
        }
    response = fetch(req_variantlist)
    return response["data"]["gene"]["variants"]

prettyprint(get_variant_list("ENSG00000010610"))

Thanks, I updated my code to include the reference genome as well (you're right, it no longer works without that).

I created a python package based on SQLite databases, where you can easily query all gnomAD variants for GRCh37/38. https://github.com/KalinNonchev/gnomAD_DB I have precomputed SQLite databases for gnomAD WGS for GRCh37/38 in the description of the package. Please take a look there."

https://github.com/KalinNonchev/gnomAD_MAF would solve your problem I think.
It is in python and you can annotate fast variants with their allele frequencies.


text:
Fast way to count number of reads and number of bases in a fastq file?
I am looking for a tool, preferably written in C or C++, that can quickly and efficiently count the number of reads and the number of bases in a compressed fastq file. I am currently doing this using zgrep and awk:
zgrep . foo.fasq.gz |
     awk 'NR%4==2{c++; l+=length($0)}
          END{
                print "Number of reads: "c;
                print "Number of bases in reads: "l
              }'

The zgrep . will print non-blank lines from the input file and the awk 'NR%4==2 will process every 4th line starting with the second (the sequence).
This works fine, but can take a very long time when dealing with large files such as WGS data. Is there a tool I can use (on Linux) that will give me these values? Or, if not, I'm also open to suggestions for speeding up the above command.

I don’t think zgrep . fulfils any tangible purpose. You should be able to leave it off entirely (replaced with zcat).

It might be worth noting that the FASTQ specification (such as it is) allows for line breaks in the sequence and qual strings, so simply taking the second of every group of 4 lines is not guaranteed to work. (see https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2847217/#__sec7title)

@KonradRudolph I use that for the unlikely case that a fastq file contains blank lines. The zgrep . will only print non-blank lines (granted, that will still count lines that have nothing but whitespace since those aren't technically empty, but it's better than nothing).

I don’t think there’s a risk of that. The bigger risk is outlined by @sjcockwell (but depending on the provenance of the files even that can be ignored).

@sjcockell yeah, I know and that's another reason why I'd rather use a more sophisticated, dedicated tool. However, I have yet to actually encounter a file with >4 lines per record in the wild, so I am not too worried about that. I should though, you're quite right.

@KonradRudolph agreed, but it is a simple enough check to add and replacing it with zcat wouldn't make all that much of a difference, I think. My main concern here is to find a more sophisticated tool that can deal with such issues robustly.

@terdon you're right that this is so rare it can be effectively ignored. I thought it was important to point out though, as so few people realise it's a possibility.

@sjcockell yes, absolutely. Very valid point.

Do you have access to a block-compressed file? It might be possible to use a parallel approach to this work to speed things up considerably.

@AlexReynolds I don't know. I'm afraid I don't know much about how compression works. The files I need to process are regular ASCII text files (fastq) compressed using gzip (usually GNU gzip, possibly BSD gzip sometimes if clients are using macs).

Use pigz https://zlib.net/pigz/ to do the decompression you'll get a good speed up from that alone.

@MatthewBashton I doubt that will make much difference. According to the pigz manual, it will compress in parallel but not decompress: Decompression can’t be parallelized, at least not without specially prepared deflate streams for
that purpose. As a result, pigz uses a single thread (the main thread) for decompression, but will create three other threads for reading, writing, and check calculation, which can speed up decompression under some circumstances. So maybe a little faster (testing it now) but I don't expect much difference for decompression.

@MatthewBashton it turns out (I tested on two different sets of reads) that unpigz is actually slower than zgrep ..

@terdon I disagree, yes the deflate process is only single threaded, I knew this already, BUT the other 3 threads for read, write and checksum offer a noticeable speedup (especially for larger multi GB files which are the norm) over a single thread in gzip doing all the work.
In my testing gzip is over 10 seconds faster than zgrep BUT pigz is faster again ~30 seconds, it's also faster than the kseq using the latest klib.h. I'm not sure how you're getting that result, perhaps this only manifests on larger gzipped files.  I've included the timings and test file in my answer.

Checkout Heng's bencharks https://github.com/lh3/biofast

It's difficult to get this to go massively quicker I think - as with this question working with large gzipped FASTQ files is mostly IO-bound. We could instead focus on making sure we are getting the right answer.
People deride them too often, but this is where a well-written parser is worth it's weight in gold. Heng Li gives us this FASTQ Parser in C.
I downloaded the example tarball and modified the example code (excuse my C...):
#include <zlib.h>
#include <stdio.h>
#include "kseq.h"
KSEQ_INIT(gzFile, gzread)

int main(int argc, char *argv[])
{
    gzFile fp;
    kseq_t *seq;
    int l;
    if (argc == 1) {
        fprintf(stderr, "Usage: %s <in.seq>\n", argv[0]);
        return 1;
    }
    fp = gzopen(argv[1], "r");
    seq = kseq_init(fp);
    int seqcount = 0;
    long seqlen = 0;
    while ((l = kseq_read(seq)) >= 0) {
        seqcount = seqcount + 1;
        seqlen = seqlen + (long)strlen(seq->seq.s);
    }
    kseq_destroy(seq);
    gzclose(fp);
    printf("Number of sequences: %d\n", seqcount);
    printf("Number of bases in sequences: %ld\n", seqlen);
    return 0;
}

Then make and kseq_test foo.fastq.gz.
For my example file (~35m reads of ~75bp) this took:
real    0m49.670s
user    0m49.364s
sys     0m0.304s

Compared with your example:
real    0m43.616s
user    1m35.060s
sys     0m5.240s

Konrad's solution (in my hands):
real    0m39.682s
user    1m11.900s
sys     0m5.112s

(By the way, just zcat-ing the data file to /dev/null):
real    0m38.736s
user    0m38.356s
sys     0m0.308s

So, I get pretty close in speed, but am likely to be more standards compliant. Also this solution gives you more flexibility with what you can do with the data.
And my horrible C can almost certainly be optimised.

Same test, with kseq.h from Github, as suggested in the comments:
My machine is under different load this morning, so I've retested. Wall clock times:
OP: 0m44.813s
Konrad: 0m40.061s
zcat > /dev/null: 0m34.508s
kseq.h (Github): 0m32.909s
So most recent version of kseq.h is faster than simply zcat-ing the file (consistently in my tests...).

Wow. So the C approach is actually marginally slower than the zgrep/awk I used and significantly slower than Konrad's trick of switching to wc.  That's surprising. I guess there must be some fancy parsing going on inside kseq_read.

The kseq.h from the tarball is outdated. The one from github is faster. That said, decompression is the bottleneck as you have observed.

@user172818 will that always be the case, irrespective of the size of the file? I would have thought that the decompression would be less of a bottleneck than the processing of each line when we're dealing with millions of lines. Is that wrong?

I think with respect to speed of core utility base solutions, people often overlook just how optimised these programs are (which of course are themselves C binaries), they have been around for a long time and are not going away.

Decompression will always be the bottleneck, it's the same situation we saw with https://bioinformatics.stackexchange.com/questions/361/what-is-the-fastest-way-to-calculate-the-number-of-unknown-nucleotides-in-fasta

pigz | awk | wc is the fastest method
First off for benchmarks with FASTQ it's best to use a specific real-world example with a known answer.  I've chosen this file:
ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG01815/sequence_read/ERR047740_1.filt.fastq.gz
as my test file, the correct answers being:
Number of reads: 67051220
Number of bases in reads: 6034609800

Next we want to find the fastest way possible to count these, all timings are the average wall-clock time (real) of 10 runs collected with the bash time on an otherwise unloaded system:
zgrep
zgrep . ERR047740_1.filt.fastq.gz |
     awk 'NR%4==2{c++; l+=length($0)}
          END{
                print "Number of reads: "c;
                print "Number of bases in reads: "l
              }'

This is the slowest method with an average run-time of 125.35 seconds
gzip awk
Using gzip we gain about another 10 seconds:
gzip -dc ERR047740_1.filt.fastq.gz |
     awk 'NR%4==2{c++; l+=length($0)}
          END{
                print "Number of reads: "c;
                print "Number of bases in reads: "l
              }'

Average run-time is 116.69 seconds
Konrad's gzip awk wc variant
fix_base_count() {
    local counts=($(cat))
    echo "${counts[0]} $((${counts[1]} - ${counts[0]}))"
}

gzip -dc ERR047740_1.filt.fastq.gz \
    | awk 'NR % 4 == 2' \
    | wc -cl \
    | fix_base_count

This runs slower on this test file than the gzip awk variant of the solution, average run-time is 122.28 seconds.
kseq_test using latest kseq.h from klib
Code compiled with: gcc -O2 -o kseq_test kseq_test.c -lz where kseq_test.c is Simon's adaptation of Heng Li's FASTQ parser.
kseq_test ERR047740_1.filt.fastq.gz
Average run-time is 99.14 seconds, which is better than the gzip core utilities based solution so far, but we can do better!
piz awk
Using Mark Adler's pigz as a drop-in replacement for gzip, note that pigz gives us a speed gain as on top of gzip as in addition to the main deflate thread it uses another 3 threads for reading, writing and checksum calculations, see the man page for details.
pigz -dc ERR047740_1.filt.fastq.gz |
     awk 'NR%4==2{c++; l+=length($0)}
          END{
                print "Number of reads: "c;
                print "Number of bases in reads: "l
              }'

Average run-time is now 93.86 seconds, this is ~5 seconds faster than the kseq based C code but we can further improve the benchmark.
pigz awk wc
Next we use pigz as a drop in replacment for Konrad's wc variant of the awk based solution.
fix_base_count() {
    local counts=($(cat))
    echo "${counts[0]} $((${counts[1]} - ${counts[0]}))"
}

gzip -dc ERR047740_1.filt.fastq.gz \
    | awk 'NR % 4 == 2' \
    | wc -cl \
    | fix_base_count

Average run-time is now down to 83.03 seconds, this is ~16 seconds faster than the kseq based solution and ~42 seconds faster than the OPs zgrep based solution.
Next as a baseline lets see just how much of this run-time is due to decompression of the input fastq.gz file.
gzip alone
gzip -dc ERR047740_1.filt.fastq.gz > /dev/null
Average run-time: 105.95 seconds, so the gzip based solutions (which also includes zcat and zgrep as these are provided by gzip) are never going to be faster than kseq_test.
pigz alone
pigz -dc ERR047740_1.filt.fastq.gz > /dev/null
Average run-time: 77.66 seconds, so quite clearly the additional three threads for read, write and checksum calculation offer a useful advantage.  What's more this speed-up is greater when leveraging the awk | wc based solution, it's not clear why, but I expect this is due to the extra write thread.
Interestingly average CPU usage across all threads is quite revealing for the various answers, I've collated these stats using GNU time /usr/bin/time --verbose
zgrep based solution 133% - must be more than one thread somehow
gzip | awk based solution 99% - all gzip based solutions run single-threaded at 99% CPU usage
pigz | awk 147%
gzip | awk | wc 99% as with gzip
pgiz | awk | wc 155%
kseq_test 99%
gzip > dev/null 99%
pigz > dev/null 155%
Whilst the main deflate thread in pigz will run at 100% CPU load the extra 3 don't quite fully occupy additional cores to 100% (as is evidenced by average CPU usage of ~150%) they do however clearly result in reduced run-time.
I'm using Ubuntu 16.04.2 LTS**, my gzip, zcat, zgrep versions are all gzip 1.6 and pigz is version 2.3.1.  gcc is version 5.4.0
** I think my patch level is actually 16.04.4 but I've not rebooted for 170 days :p

Nice and thorough! Thank you. I don't know why your tests showed pigz to be faster when it was slower in mine. I suspect it might be because I was testing using files stored on an NFS volume so my times could well have depended more on the network latency than anything else. I'll repeat your tests using a local drive and report back.

OK, confirmed. On your test file and on a local drive, pigz was indeed faster than gzip or zcat | grep.

wow excellent. thanks

The following is more than twice as fast; however, wc counts newline characters as well. We thus need to subtract the line count from the base count (using Bash):
fix_base_count() {
    local counts=($(cat))
    echo "${counts[0]} $((counts[1] - counts[0]))"
}

gunzip -c "$file" \
| awk 'NR % 4 == 2' \
| wc -cl \
| fix_base_count

All the caveats from Simon’s comment apply: this assumes the “simple” FASTQ format, where each record consists of exactly four lines. I think this is true for all files produced by Illumina sequencers and downstream tools.

What did you test this on? I used a 2.6G exome fastq file and this approach (without the bash function, so even faster) took  1m6.754s while my version using awk to count everything took 1m14.228s. What did you measure to see a two-fold increase in speed?

@terdon I’ve used the time command for benchmarking, and working on a 2.2GiB FASTQ file. The times on my machine were 9m5.628s and 3m5.825s real time, respectively. The user time in both cases was slightly elevated (11m vs 5m) to make the  difference a factor 2.

@terdon But given Simon’s answer and common sense, this is a fishy result indeed. I don’t know what my awk version is doing here.

I get fairly quick results with my fastx-length.pl script, with the added bonus of being able to handle multi-line FASTQ files and displaying additional read-length QC statistics:
time zcat albacored_all.fastq.gz | /bioinf/scripts/fastx-length.pl > /dev/null
Total sequences: 301135
Total length: 283.902419 Mb
Longest sequence: 5.601 kb
Shortest sequence: 6 b
Mean Length: 942 b
Median Length: 999 b
N50: 111835 sequences; L50: 1.103 kb
N90: 245243 sequences; L90: 608 b

real    0m8,802s
user    0m16,584s
sys 0m0,260s

Versus the script you have provided:
zcat albacored_all.fastq.gz | awk 'NR%4==2{c++; l+=length($0)}
          END{
                print "Number of reads: "c;
                print "Number of bases in reads: "l
              }'
Number of reads: 301135
Number of bases in reads: 283902419

real    0m8,382s
user    0m10,216s
sys 0m0,332s

Cat to /dev/null for comparison:
time zcat albacored_all.fastq.gz > /dev/null

real    0m7,877s
user    0m7,856s
sys 0m0,020s

I suspect that something using bioawk might be a bit faster (and similarly FASTQ-compliant).

bioawk is not that efficient, at least in my tests: https://bioinformatics.stackexchange.com/a/961/292

I hava implemented seqtk_counts using kseq.h from  klib
Just a few line of Codes:
#include <stdio.h>
#include <zlib.h>
#include "kseq.h"

KSEQ_INIT(gzFile, gzread)

int main(int argc, char *argv[]){

gzFile fp;
kseq_t *seq;

int l = 0;

int64_t total = 0;
int64_t lines = 0;

if (argc == 1) {
    fprintf(stderr, "Usage: %s <fastq> <sample>\n", argv[0]);

    return 1;
}

fp = strcmp(argv[1], "-")? gzopen(argv[1], "r") : gzdopen(fileno(stdin), "r");
seq = kseq_init(fp);

while ((l = kseq_read(seq)) >= 0){
    total += seq->seq.l;
    lines += 1;
}

printf("%s\t%lld\t%lld\n", argv[2] ,(long long)lines, (long long)total);
kseq_destroy(seq);
gzclose(fp);
return 0;

}
Compile it:
gcc  -O2  seqtk_counts.c  -o  seqtk_counts  -Iklib  -lz

Usage:
seqtk_counts foo.fasq.gz foo
or
cat foo.fasq.gz | seqtk_counts  - foo

Thanks for the edit! :) But isn't this basically the same program as shown in  sjcockell's answer?

Using pyGATB
(I use the same file as in https://bioinformatics.stackexchange.com/a/400/292, same workstation as in https://bioinformatics.stackexchange.com/a/380/292)
$ time python3 -c "from gatb import Bank; seq_lens = [len(seq) for seq in Bank('SRR077487_2.filt.fastq.gz')]; print('Number of reads: %d' % len(seq_lens), 'Number of bases in reads: %d' % sum(seq_lens), sep='\n')"
Number of reads: 23861612
Number of bases in reads: 2386161200

real    0m41.122s
user    0m40.788s
sys     0m0.312s

It is quite faster than bioawk:
$ time bioawk -c fastx '{nb_seq+=1; nb_char+=length($seq)} END {print "Number of reads: "nb_seq"\nNumber of bases in reads: "nb_char}' SRR077487_2.filt.fastq.gz
Number of reads: 23861612
Number of bases in reads: 2386161200

real    1m3.182s
user    1m2.916s
sys     0m0.268s

But not so much than the OP example:
$ time zgrep . SRR077487_2.filt.fastq.gz | awk 'NR%4==2{c++; l+=length($0)} END{print "Number of reads: "c; print "Number of bases in reads: "l}'
Number of reads: 23861612
Number of bases in reads: 2386161200

real    0m47.127s
user    1m36.292s
sys     0m6.796s

Or than the wc based solution:
$ fix_base_count() {
>     local counts=($(cat))
>     echo "${counts[0]} $((${counts[1]} - ${counts[0]}))"
> }
$ time gunzip -c SRR077487_2.filt.fastq.gz | awk 'NR % 4 == 2' | wc -cl | fix_base_count
23861612 2386161200

real    0m44.915s
user    1m12.000s
sys     0m6.972s

I didn't compare with C-based solutions.
The zcat to /dev/null reference is the following:
$ time zcat SRR077487_2.filt.fastq.gz > /dev/null

real    0m39.745s
user    0m39.464s
sys     0m0.252s

I'm still impressed by pyGATB speed

if the data is in SRA, there is sra-stat utility that returns reads,bases and quality distribution. these are stored in the SRA file. use --quick to get the stored stats or --statistics to calculate additional values broken down per readgroup/barcode.
sra-stat --quick SRR077487

Use assembly-stats!
$ assembly-stats barcode13_filtered.fastq
stats for barcode13_filtered.fastq
sum = 2080834976, n = 656192, ave = 3171.08, largest = 15321
N50 = 3598, n = 225863
N60 = 3263, n = 286569
N70 = 2920, n = 353904
N80 = 2548, n = 429982
N90 = 2081, n = 519671
N100 = 500, n = 656192
N_count = 0
Gaps = 0

real    0m7.418s
user    0m5.976s
sys     0m1.275s

Works on any fasta/fastq files, even those that aren't assemblies. And it's pretty fast!
Download from here: https://anaconda.org/bioconda/assembly-stats

Thanks, that looks interesting. How does the time compare to the existing answers here? Is this faster? And what are the numbers reported? Is that 656192 reads? And 2080834976 bases?

Yes, yes, and yes! 656192 reads, 2080834976 bases, with the respective N50 etc values. And compared to the ones above, it is lightning fast.
real 0m7.418s
user 0m5.976s

Could you show an example of that? You are running this on a single, decompressed file here, while the question is asking about compressed files and the other answers have shown that the rate limiting step is usually the decompression. So how fast is this approach when run on a compressed file and taking the decompression time into account? Also, in order for the reported times to be meaningful, please also include the time taken by one of the other methods (you can't compare times across different machines and files).

This should too good to be true, is it because the fastq is uncompressed? @sjcockell benchmark is on a comparable ~2.6e9 bases file, but it's about 5x slower... But it's a compressed file that was used.


text:
Trouble using biomaRt to retrieve hgnc symbols from Ensembl transcript ids
I have a matrix of gene counts which I'm going to use as input for DESeq.  Right now, each gene is labeled by its Ensemble transcript ID, but I'd like to convert these to their HGNC symbols before I input them into DESeq for analysis. I'm attempting to do this conversion using biomaRt, however, this is my first time using the program, and I'm running into issues.  Here is my current code:
library("biomaRt")
mart <- useMart(biomart = "ENSEMBL_MART_ENSEMBL", dataset = "mmusculus_gene_ensembl")
transcript.ids <- rownames(txi.kallisto$counts)
hgnc_symbols <- getBM(attributes = "hgnc_symbol", filters = "ensembl_transcript_id", values = transcript.ids, mart = mart)

Printing hgnc_symbols shows that it isn't what I want:
>ghnc_symbols
[1] hgnc_symbol
<0 rows> (or 0-length row.names)

Lastly, here is what the input transcript IDs look like:
> head(rownames(txi.kallisto$counts))
[1] "ENSMUST00000178862.1" "ENSMUST00000178537.1" "ENSMUST00000196221.1"
[4] "ENSMUST00000179664.1" "ENSMUST00000177564.1" "ENSMUST00000179520.1"

I may be missing something here, but why would you be looking for HGNC (Human Genome Organisation  Gene Nomenclature Committee) names for mouse genes?

You need to specify the number without the version. Instead of "ENSMUST00000178862.1" just "ENSMUST00000178862":
You can do this with one more line:
g <- gsub("\\..*", "", rownames(txi.kallisto$counts))
(hgnc_symbols <- getBM(attributes = c("hgnc_symbol", "chromosome_name", "ensembl_transcript_id"), filters = "ensembl_transcript_id", values = g, mart = mart))
##  hgnc_symbol chromosome_name ensembl_transcript_id
##1          NA              14    ENSMUST00000177564
##2          NA               6    ENSMUST00000178537
##3          NA               6    ENSMUST00000178862
##4          NA              12    ENSMUST00000179520
##5          NA              14    ENSMUST00000179664
##6          NA              14    ENSMUST00000196221

With these transcripts it doesn't seem that there is an equivalent hgnc_symbol. But is recognized, otherwise they are not recognized:
(hgnc_symbols <- getBM(attributes = c("hgnc_symbol", "chromosome_name", "ensembl_transcript_id"), filters = "ensembl_transcript_id", values = rownames(txi.kallisto$counts), mart = mart))
## [1] hgnc_symbol           chromosome_name       ensembl_transcript_id
##<0 rows> (or 0-length row.names)

However you might be interested in doing something like this:
human = useMart("ensembl", dataset = "hsapiens_gene_ensembl")
mouse = useMart("ensembl", dataset = "mmusculus_gene_ensembl")
getLDS(attributes = "ensembl_transcript_id",
       filters = "ensembl_transcript_id", values = g,mart = mouse,
      attributesL = "hgnc_symbol", martL = human)
## The query to the BioMart webservice returned an invalid result: the number of
## columns in the result table does not equal the number of
## attributes in the query. Please report this to the mailing list.


text:
How to find hairpin loop structure in a large set of long non coding RNA transcripts
I have done de novo assembly of pair end raw read sequences and resulted transcripts sequences were separated based on coding potential into two categories: long non-coding RNA transcripts and coding RNA transcripts.
I want to find how many of long non coding can make hairpin loop structure and how can I annotate these miRNA which are generated from long non coding RNA in plants.

Welcome to Bioinformatics! Is there anything you tried before posting? Is there any reason why you tagged the question with Bioconductor (Maybe you want to restrict to use a Bioconductor package, or have used one and you can't manage to make it classify your lnRNA)? Also, I think there are two questions here: how to find if a lncRNA can make hairpins and annotate RNA to miRNA. It would be better to focus in just one question (perhaps the first one and later on the second one).

For the first question: "how many of long non coding can make hairpin loop structure?"
The easiest thing to would be to run your sequences through an RNA secondary structure prediction tool. There's a few tools for doing this but the most commonly used are RNAfold from the ViennaRNA package and mfold.
This will give you an output that looks like this, where the parentheses indicate which bases are paired with each other.

GGGCUAUUAGCUCAGUUGGUUAGAGCGCACCCCUGAUAAGGGUGAGGUCGCUGAUUCGAAUUCAGCAUAGCCCA
(((((((..((((.........)))).(((((.......))))).....(((((.......)))))))))))).

To see if there's hairpins, you'll need to annotate the elements within it. Take a look at this question on biostars for an example of how to do this. In short, you can use the forgi library to produce an annotation like this:

(((((((..((((.........)))).(((((.......))))).....(((((.......)))))))))))).
sssssssmmsssshhhhhhhhhssssmssssshhhhhhhsssssmmmmmssssshhhhhhhsssssssssssse

If there's an h in the output annotation, you have a hairpin. The remainder of your question about annotating miRNAs will be left to somebody else.

hairpins are quite easy to identify using the bracket notation; just look for brackets that are in a left-right orientation with nothing in between (e.g. (.), (..), (...)).


text:
Is there any way of using biopython to write Swissprot files?
I have the following mwe for filtering a Swissprot file based on a certain feature, in this case, transmembrane regions.
from Bio import SeqIO
records=[]
for record in SeqIO.parse("Input.txt", "swiss"):
    transmembrane_protein=False
    print record.id
    for i, feature in enumerate(record.features):
        if feature.type == "TRANSMEM":
            transmembrane_protein=True
    if transmembrane_protein==True:
        records.append(record)
SeqIO.write(records, "Output.txt", "swiss")

The script works when SeqIO.write(records, "Output.txt", "swiss") becomes SeqIO.write(records, "Output.txt", "fasta")
However, such a method is not yet supported.
ValueError: Reading format 'swiss' is supported, but not writing

From the docs, I see that writing is not supported for swiss:

Note that while Bio.SeqIO can read all the above file formats, it cannot write to all of them.

Is there any unofficial way of using biopython/python to write swissprot files from parsed swissprot files?

I read that "SeqIO cannot output the file format (e.g. the plain text SwissProt format, "swiss") ". But you could paste the proteins without using Biopython once you know the id of the features you want

It must be in some of the docs somewhere I'm sure, I just haven't seen it explicitly stated, but it is pointed out in the error handling so that's good enough for me to know there is no easy way. I'm dealing with several proteomes. Manual pasting wouldn't be fun ;) I could really do with a programmatic solution.

You could use grep, or awk to select the records and output them in another file using a pipe.

@Llopis That would certainly be one way forward. However I'd like to keep the script in a single python file because it's cross platform and transportable.

Using SeqIO.index rather than SeqIO.parse lets you read all the records into a dict, from which you can then extract the raw entry:
from Bio import SeqIO

record_dict = SeqIO.index('Input.txt', 'swiss')
for key in record_dict:
    print(record_dict.get_raw(key).decode())

Now you should be able to apply your test for a transmembrane protein to each entry, and write out only the ones you want to keep.
from Bio import SeqIO

record_dict = SeqIO.index('Input.txt', 'swiss')
out = open('Output.txt', 'wb')
for key in record_dict:
    record = record_dict[key]
    for i, feature in enumerate(record.features):
    if feature.type == "TRANSMEM":
        out.write(record_dict.get_raw(key).decode())

Perfect solution, worked like a charm.

I don't know how to do this with biopython, but it's simple enough in Perl using Swisskinfe:
#!/usr/bin/env perl
use strict;
use warnings;
use SWISS::Entry;

my $flat=$ARGV[0]||die("Need a flat file\n");

## Change the line termination string so we read an
## entire entry at a time
local $/ = "\n//\n";

my $fh;
## Deal with compressed or uncompressed files
if($flat=~/\.gz$/){
    open($fh,"zcat $flat |")|| die("cannot open flat file $flat : $! \n");
}
else{
    open($fh,"$flat")|| die("cannot open flat file $flat : $!: $@\n");
}
## Parse the file.
entry: while (my $current_entry = <$fh>) {
    ## Read the current entry as an "entry" object
    my $entry = SWISS::Entry->fromText($current_entry);
    ## Iterate over the entry's features
    foreach my $feature ($entry->FTs->elements()) {
        ## If any of the features are 'TRANSMEM', print the entry
        if ($feature->[0] eq 'TRANSMEM') {
            print "$current_entry\n";
            ## Skip to the next entry
            next entry;
        }
    }
}

You need to install the Swissknife module first:

Download Swissknife_1.73.tar.gz from here (SourceForge link).
Install
tar xvzf Swissknife.tar.gz
cd Swissknife_1.73/
perl Makefile.PL
sudo make install

Now, make the script above executable and run with:
foo.pl file.flat


text:
How to convert the .vcf (imputed) file with GT:GP format to GT:DS?
I have the genotyped data from impute2 output in .gen format (imputed to 1000G P3). The file has genotype posterior probabilities (GP:3 values per variant). I have converted .gen to .vcf using qctools and the .vcf file has GT:GP format. I need to convert the .vcf file with GT:GP format to GT:DS. Genotype dosages are recommended for use in qtltools/fastqtl analysis. However, I cannot find a tool that would keep the .vcf format and convert GP to DS. Any help much appreciated!

You can do this in Hail.
Here's the rough code to do it (0.1 versions).
Setup:
from hail import *
hc = HailContext()

Import the .gen file. VCF works too:
dataset = hc.import_gen(
    'src/test/resources/example.gen',
    'src/test/resources/example.sample')

Remap the genotype schema and export to VCF:
dataset.annotate_genotypes_expr('g = {GT: g.call(), DS: g.dosage()}')\
    .export_vcf('/tmp/out.vcf.bgz')

Take a look at the getting started page if you want to try it out!
I should note that you may be able to do QTL analyses in Hail, depending on the method you want to use. See blog post here.

Hm, I didn't know that the plugin existed so I wrote my own script to convert GP to minor allele dosage on github. Maybe someone else will find it useful :) https://github.com/7methylg/VCF-GP-to-DS

There's the dosage plugin for bcftools, but it only outputs tab separated values. It would not be too hard to extend the plugin to output a VCF with the DS tag instead, but it has not been done yet.  There's a good chance the bcftools devs would respond to a feature request...
In any case, this code:
curl https://raw.githubusercontent.com/samtools/bcftools/develop/test/convert.vcf > convert.vcf
bcftools +dosage convert.vcf > output.tsv
head -2 output.tsv

has the output:
#[1]CHROM       [2]POS  [3]REF  [4]ALT  [5]NA00001      [6]NA00002      [7]NA00003      [8]NA00004      [9]NA00005      [10]NA00006     [11]NA00007     [12]NA00008    [13]NA00009     [14]NA00010
X       2698560 G       A       0.1     0.0     0.1     0.2     0.3     0.2     0.2     0.2     0.2     0.1

This is using bcftools version 1.3.1.
Here is an excerpt from the bcftools manual for the dosage plugin:

dosage
print genotype dosage. By default the plugin searches for PL, GL and GT, in that order.


text:
Minimizing particle-flow grid mapping time of a protein surface using spherical-coordinates and differentiation encoding
This is an interesting problem - I was wondering if anyone has a creative solution.
So I have a vector of vertices representing atoms in a protein, as well as 6 variables containing the absolute minimum/maximum bound of the set at each direction.
I need to build a grid that surrounds the protein/vertices in a homogeneous layer.
Right now I just loop through every square of the grid bound by the minimum-maximum values and then loop through every vertex to calculate the minimum distance between my point and every atom in the protein, and see if that point is at an acceptable distance. This process is quite expensive.
Do you guys have any creative answers as to how minimize the time?
Thanks!
P.S. the list is ordered by phyiical linkage of amino acids in a  protein, which means it does not always have a relationship to the location of the amino acid in protein space.

The question is clear to me apart from the sentence "build a grid that surrounds...". Usually you define a grid with a spacing of your choice and then only assign numbers to the nodes in the grid. Maybe you could edit the question and elaborate.

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


text:
Pooling data in metagenome assembly
I have 12 human gut microbiome WGS Nextseq reads (151 bp paired end). What will be an effective strategy to assemble a metagenome?
Let us say I have already filtered the fastq for quality, adapter sequence and host contamination (human, in this case).
1) Should I concatenate all the R1 reads as one single R1 read and one Single R2 read?
cat Sample[1..12].R1 > Single_R1.fastq
cat Sample[1..12].R2 > Single_R2.fastq

and then use Diginorm to normalize the Single_R1.fastq and Single_R2.fastq. Subsequently, feed these fastq files into any metagenome assembler such as Megahit, MetaSPAdes?
Normalize the output by using CD-HIT or similar tool to remove duplicates and filter by contig length.
OR
2) Perform metagenome assembly for each of the samples individually after applying filtering, removing adapters and host contamination.
R1=(*_R1_001.filtered.fastq)
R2=(*_R2_001.filtered.fastq)
for ((i=0;i<=${#R1[@]};i++)); do
  /bin/metagenome-assembler -1 "${R1[i]}" -2 "${R2[i]}" -o ${R1[i]%.*}.contigs.fa;
done

Followed by combining all the contigs.fa into one mega_contigs.fa
cat *.contigs.fa > Mega_contigs.fa

and use CD-HIT or similar tool to remove duplicates.

I think assembling data separately is better idea, since there will more noise in combined signals of different individuals and also resource wise. However, I guess you should not just cat assemblies together, you want to avoid redundancy (discard what is shared between individuals), not entirely sure how to achieve that though.

With metagenomes, there's several strategies and no one-size-fits-all.
 You can spend a lot of time playing with the data to try to get the "best assembly."  But it depends on what you're trying to do, e.g. whether or not you are mostly interested in the highest abundant species or are interested in the less abundant ones too.
But to propose a strategy that seems reasonable and hasn't been mentioned yet, my initial thoughts are that since the gut microbiome isn't complex, you can likely get decent assemblies from each run separately.  Then you can run a contig clustering algorithm (e.g. MetaBAT) to bin the contigs into genomes.  Next compare the bins from the various assemblies (i.e. cluster the clusters).
It may or may not be worth the trouble to try to combine the related bins' contigs into a consensus (pan)genome sequence.  If you want to, methods that have been used successfully include shredding the contigs into overlapping chunks and assembling with an OLC assembler.  But if you're interested in KO:read-counts (for example), a consensus genome sequence isn't really required.  Cross-mapping the unassembled reads can be useful too (i.e. a sample may have a low-abundance OTU that doesn't produce good contigs, but it can be associated with an OTU that is better assembled in another sample).
If you provide more detail on your goal, I can comment further.

If your goal is to bin the resulting contigs into genomes, then you should do option #1, pooling the reads and assembling into one set of contigs.

I also think pooling is the better option, followed by partitioning by coverage/taxonomy of contigs.
Maybe check out BlobTools, which helps with filtering read-pairs by taxonomy of contigs they contribute to and also does nice visualisations of assemblies.
Workflow B seems to be what you want.
Disclaimer: I am the developer of this tool.


text:
Run kallisto iteratively across many samples
I am on a Mac using UNIX. I am trying to use the kallisto quant command on all files in a directory (instead of manually entering them). Because I'm running the analysis against the same index file, I first enter the following:
./kallisto index -i --index --make-unique index.fa

This successfully creates an index file. Then, I tried this:
for file in *.fasta; do kallisto quant --single -l 200 -s 0.1 -o $file-aligned
"$file"; done

While the $file-aligned folders are created, they're empty. I get an error that states that the index file is missing. So I assumed I would need to specify the index file in the command line. I then tried this:
for file in *.fasta; do kallisto quant -i index.fa --single -l 200 -s 0.1 -o
$file-aligned  "$file"; done

I get an error that says "Error: incompatible indices. Found version 0, expected version 10.
I have now also tried the following:
for file in *.fasta; do ./kallisto quant -i index --single - 200 -s 0.01 -o
"${file}-aligned" "${file}";

but I get an error saying the index file is not found.

I am now running
or file in *.fasta; do ./kallisto quant -i index --single -l 200 -s 0.1 -o $file-aligned

"$file"; done
This sort of works. The index file error messages are gone, but it runs quant on the same file 10 times instead of each file once. If I run echo $file it spits back the file that my command is running repeatedly.
If I close terminal, start a new session and run the same command, I get
Error: file not found ""

Why are you using fasta files for reads? Kallisto expects FASTQ files as input.

Does kallisto use quality scores of reads? I had an impression that it does not, but I don't have a reference to backup it....

I down-voted because this post does not pose a question.

Running the jobs in parallel
Iteratively is one solution, but since you're not using bootstrapping, which can be run multithread (and is not needed if you're taking the popular tximport route) you can do better using GNU parallel, this will enable you to run as many jobs as you have execution threads on your system simultaneously.  Note that Kalliso only runs multithreaded with -t if you're using bootstrapping.

You'll need to make the index first:
kallisto index --make-unique -I index.fa.idx index.fa

Next lets define the input FASTQ files and the sample IDs you want to use using a tab-delimited input file like so:
sample_1    sample_1_R1.fastq.gz
sample_2    sample_2_R1.fastq.gz
sample_3    sample_3_R1.fastq.gz

Here the first column is the sample ID which will be used for output and the second separated by a tab character is the input file.
You can produce this using a text editor but as a short-cut you can feed the list of fastq.gz you have into a list via:
ls -1 *.fastq.gz > editme

Next we create a list of jobs to run on our system using a utility script Make_job_list_kallisto.sh which consumes our tab delimited input file list above.  I've used your -l and -s parameters below, note that with paired end data normally these are estimated for you.
#!/bin/bash -eu

[ $# -ne 2 ] && { echo -en "\n*** This script generates jobs for GNU parallel. *** \n\n Error Nothing to do, usage: < input tab delimited list > < output run list file >\n\n" ; exit 1; }
set -o pipefail

# Get command-line args
INPUT_LIST=$1
OUTPUT=$2

# Set counter
COUNT=1
END=$(wc -l $INPUT_LIST | awk '{print $1}')

echo " "
echo " * Input file is: $INPUT_LIST"
echo " * Number of runs: $END"
echo " * Output job list for GNU parallel saved to: $OUTPUT"
echo " "

# Main bit of command-line for job
CMD="kallisto quant -i index.fa.idx --single -l 200 -s 0.1"

# Main Loop
[ -e $OUTPUT ] && rm $OUTPUT
while [ $COUNT -le $END ];
do
    LINE=( $(awk "NR==$COUNT" $INPUT_LIST) )
    # Make file list
    echo "Working on $COUNT of $END Sample ID: ${LINE[0]}, Files ${LINE[@]:1}"
    echo "$CMD -o ${LINE[0]} ${LINE[@]:1}" >> $OUTPUT
    ((COUNT++))
done

We use the script like so:
 Make_job_list_kallisto.sh input_file_list.txt job_list.txt

Now we have a job list for GNU parallel we can run our jobs in parallel like so:
parallel --progress --jobs 4 --joblog kallisto_joblog.txt < job_list.txt

This will run 4 of the kallisto jobs simultaneously (Increase this if you have more threads - don't mind consuming all the IO/CPU threads on your system) producing a nice progress report of those completed, additionally a job log will be written to the file kallisto_joblog.txt.

For paired end data I already have these scripts + similar instructions on GitHub which people might find useful.

@ZincFingers, so -I upper case for specifying the index in my run generation scrip should in fact be lower case -i.  I've tested this and it works for me.  This might be your issue or you might have some other error.  If it's still not working try pasting in a few lines of your jobs list here.

The error I was getting was my fault, I fixed it. Now, when I run the parallel command, I get an error from kallisto saying
Error: missing read files
Error: file srr3884878.gz exists and is not a directory
It says this for all 17 of my files

I suspect you still have an issue with you input file to Make_job_list_kallisto.sh or the jobs it's generating, without seeing the input sample ID file list or the job list it's very difficult to debug.

Check the diagnostic output of Make_job_list_kallisto.sh is sane, I get:
* Input file is: test_for_BS.txt
* Number of runs: 4
* Output for GNU parallel saved to: test_job_list.txt
Working on 1 of 4 Sample ID: test1, Files ATCACG_L004_R1_001.fastq.gz
Working on 2 of 4 Sample ID: test2, Files GTCCGC_L004_R1_001.trimmed.fastq.gz
etc.

I think you might have forgotten to give an output dir or sample ID in the first column of your tab delimited input list you're giving to Make_job_list_kallisto.sh because it looks like kallisto is missing the variable which would be supplied to -o and treating the input fastq.gz file as the output dir, it's then not finding the gz.fastq file as it's trying to use that as the output dir.

You were right! I fixed the input file and it seems to have resolved the issue. Many thanks for your fantastic solution and bearing with me while I trudged through this.

Glad you got it working!

For generating the index, you should be using one of -i or --index, but not both:
./kallisto index --make-unique -i index.fa.idx index.fa

You need to point kallisto to the generated index file (usually has a .idx extension), rather than the fasta file:
for file in *.fasta;
  do kallisto quant -i index.fa.idx --single -l 200 -s 0.1 -o "${file}-aligned" "${file}";
done

Why run it sequentially? Shouldn't this be parallelized, i.e. multiprocessed

You're asking the wrong person; try adding that as a comment to the question. Kallisto quant allows multi-threaded operation with -t <threads>.

I get the same error when I point to the index file
Error: incompatible indices. Found version 0, expected version 10

maybe it's the index command that needs a bit of attention... I'll update the commands to include index generation

I managed to get this to work after fiddling with the index generation:
' for f in *.fasta' do kallisto quant -i index --single -l 200 -s 0.1 -o "${file}-aligned" "${file)"; '
But it appears that it's repeatedly analyzing one of the 10 .fasta files I placed in the directory (also I'm not sure how to get that line to appear as a code)

@gringer - just an additional comment on multithreading: as of version 0.43.1, the -t <threads> parameter doesn't work should you wish to look at the alignments (which kallisto quant can output upon supplying --pseudobam. I usually get around it by running several individual kallisto quant commands concurrently.

@gringer Using -t is fine but it’s actually almost certainly much more efficient to process the files in parallel rather than in a sequential loop.

-t is not going to help, it's only useful if you've enabled bootstrapping, otherwise you can set it as high as you want but kallisto will still run single core.  Multi-threaded execution is used only for concurrently running the bootstrapping.

If you're using tximport for downstream analysis in DESeq2 rather than Sleuth the bootstrapping won't be used as it's not taking into account quantification uncertainty via the bootstrapping procedure, so you need to consider if this is needed.