text:
How can I call structural variants (SVs) from pair-end short read resequencing data?
I have a reference genome and now I would like to call structural variants from Illumina pair-end whole genome resequencing data (insert size 700bp).
There are many tools for SV calls (I made an incomplete list of tools bellow). There is also a tool for merging SV calls from multiple methods / samples - SURVIVOR. Is there a combination of methods for SV detection with optimal balance between sensitivity and specificity?
There is a benchmarking paper, evaluating sensitivity and specificity of SV calls of individual methods using simulated pair-end reads. However, there is no elaboration on the combination of methods.
List of tools for calling structural variants:

Lumpy
BreakDancer
Manta
Delly
GRIDSS
Meerkat
Pindel
Softsv
Prism

I will just add a comment, as it's not a complete answer. Check the Genome in a Bottle consortium. There are discussions now on how to determine the best caller(s) and definition on a standard set of calls for benchmarking and testing new approaches. In my work I had good results with Socrates, now replaced with GRIDSS.

@nuin - I didn't know about Genome in a Bottle consortium, looks interesting, but I could not find any public record of a discussion. Do you have a link?

I think the best method or combination of methods will depend on aspects of the data that might vary from one dataset to another. E.g. the type, size, and frequency of structural variants, the number SNVs, the quality of the reference, contaminants or other issues (e.g. read quality, sequencing errors) etc.
For that reason, I'd take two approaches:

Try a lot of methods, and look at their overlap
Validate a subset of calls from different methods by wet lab experiments - in the end this is the only real way of knowing the accuracy for a particular case.

What do you mean by "Validate a subset of calls from different methods"?

Funny, the review paper (sim data) kind of agree with you - all methods have lot of false positives - it is better to take overlap.
However, the paper "An integrated map of structural variation" (Germain's answer) have done opposite - tweaking every software to reduce false positives and than taking them all.

By validate, I mean go in and resequence a bunch of predicted variants. All of the bioninformatic methods will give lots of false positives (and lots of false negatives, though it's harder to find the negatives of course). So the only genuine way to know the accuracy of your calls is to get independent information e.g. from resequencing.

In case if you are really dedicated to obtain perfect results, you can use strategy described there, in 1000GP 3rd Phase SV detection paper - use these tools, validate your calls with IRS test, merge calls into one callset.
If you do not wanna spend thousands human-hours as was spent during this paper preparation, from my experience, it is better to use 1 paired-end insert distance method and one read-depth based method. Each of them cover "different" regions in the genome. (even thou they have huge overlap, paired-end detection requires both SV breakpoints to be located within the regions with good mappability which is not always the case, but resolution of read-depth methods is lower in general, paired-ends works well for deletions/tandem duplications/inversions, but have troubles with non-tandem duplications).
Hope it helps.

I went though the supplement of the paper, it is really crazy - it seems that there is very little overlap of different computational methods, but the False discovery rate is for non inversions types of SVs pretty good (2 - 10%). However, they have used a huge population dataset, I have sequencing of 9 individuals including the reference.

Yeap, the overlap is kinda small, but this is mainly because different tools look for different type of CNVs (as I explained above comparing Read Depth vs Insert Distance-based methods). This is the best accuracy you can get. Having only 9 individuals, you'd better exclude CNVs that happen with frequency >1% in the population (see the paper) and then use eg DELLY and cn.mops (just an example, not sure if this is an ideal choice, but both of them are quite good).


text:
What is the best way to account for GC-content shift while constructing nucleotide-based phylogenetic tree?
Let's say I want to construct a phylogenetic tree based on orthologous nucleotide sequences; I do not want to use protein sequences to have a better resolution. These species have different GC-content.
If we use a straightforward approach like maximum likelihood with JC69 or any other classical nucleotide model, conserved protein coding sequences of distant species with similar GC-content will artificially cluster together. This will happen because GC-content will mainly affect wobbling codon positions, and they will look similar on the nucleotide level.
What are possible ways to overcome this? I considered the following options so far:

Using protein sequence. This is possible of course, but we lose a lot of information on the short distance. Not applicable to non-coding sequences.
Recoding. In this approach C and T can be combined into a single pyrimidine state Y (G and A could be also combined in some implementations). This sounds interesting, but, first, we also lose some information here. Mathematical properties of the resulting process are not clear. As a result, this approach is not widely used.
Excluding third codon position from the analysis. Losing some short-distance information again. Also, not all synonymous substitution are specific to the third codon positions, so we still expect to have some bias. Not applicable to non-coding sequence.

It should be possible in theory to have a model which allows shifts in GC-content. This will be a non time-reversible Markov process. As far as I understand there are some computational difficulties estimating likelihood for such models.

I would just add that I think there's a key assumption in the setup here: "I do  not want to use protein sequences to have a better resolution". We can decompose 'better' here - it's likely to be more precise but also more biased, the latter for all the reasons you outline.

In case you might be interested, I tested some of the approaches you mention, plus a few other recoding schemes (http://dx.doi.org/10.6084/m9.figshare.732758) in the following papers: http://arxiv.org/abs/1307.1586 and http://dx.doi.org/10.1093/molbev/msu105

There are models that take into account compositional heterogeneity both under the maximum likelihood and Bayesian frameworks. Although the substitution process is not time-reversible, the computations are simplified by assuming that the instantaneous rate matrix can be decomposed into an "equilibrium frequency vector" (non-homogeneous) and a symmetric, constant exchange rate matrix.
I guess all your suggestions are also valid, and I remember recoding being used successfully to reduce the GC-content bias (examples in the references above and here).

The following 2004 paper describes a way to model compositional changes across the tree, in a Bayesian framework: https://doi.org/10.1080/10635150490445779
A python package implementing this ("p4"), and improvements added along the years, is available here: https://github.com/pgfoster/p4-phylogenetics
To get started, you may find useful examples here: http://p4.nhm.ac.uk/scripts.html
This has been used in a few large-scale phylogenetic analyses.

The answer is the logDet algorithm was constructed to overcome GC% clustering.
At the time it was devised only a distance method was available/implemented, so it wasn't very powerful. The posts here imply that a Bayesian or ML approach are available and these do tightly stick to the model.
Original publication here

Do you have a publication or webpage in mind? Can you link it?

Link provided above. Goes back a long way ... to 1996

source: Stack Exchange

license: 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/

sort: votes


text:
How to filter out cross alignments from a BED file?
I have a BAM file:
@SQ SN:chr1 LN:248956422
@SQ SN:chrx LN:248956423
ST-E00110:348:HGVKKALXX:1:1201:5822:48670   323 chr1    9999    0   67H66M16H   chrx    1000    0   GATAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAAC  JJJJJJJJJJJJJJJJAJJJJJJJJJJJJFJJJJJJFJFJJJJJJFJJJJJJJJJJJA77FJFJJJ  NM:i:0  MD:Z:66 AS:i:66 XS:i:65 SA:Z:chr5,18606834,-,73S76M,34,0;   RG:Z:g1

There is a read aligned to chr1, and it's mate aligned to chrx.
I have a BED file:
chr1    0   100000  TestOnly

I would like to filter out everything that falls outside my BED region, that includes cross-alignments. In my example, although my read aligned to chr1 but it's mate is not. I don't want this read.
When I do:

samtools view -L test.bed test.bam

the command gives me the read because it doesn't check cross-alignments.
My solution:

samtools view -L test.bed test.bam | grep -v chrx

but this is very slow and clumsy. In my production pipeline I would have to do something like:

samtools view -L test.bed test.bam | grep -v chrx | grep -v ... | grep -v ... | grep -v ... | grep -v ...

Q: Is there a better solution?

According to the SAM specification, the 3rd field of a SAM line (RNAME) is:

RNAME: Reference sequence NAME of the alignment. If @SQ header lines
  are present, RNAME (if not ‘*’) must be present in one of the SQ-SN
  tag. An unmapped segment without coordinate has a ‘*’ at this field.
  However, an unmapped segment may also have an ordinary coordinate such
  that it can be placed at a desired position after sorting. If RNAME is
  ‘*’, no assumptions can be made about POS and CIGAR.

And the 7th field is (emphasis mine, missing "to" theirs):

RNEXT: Reference sequence name of the primary alignment of the NEXT
  read in the template. For the last read, the next read is the first
  read in the template. If @SQ header lines are present, RNEXT (if not
  ‘*’ or ‘=’) must be present in one of the SQ-SN tag. This field is set
  as ‘*’ when the information is unavailable, and set as ‘=’ if RNEXT is
  identical RNAME. If not ‘=’ and the next read in the template has one
  primary mapping (see also bit 0x100 in FLAG), this field is identical
  to RNAME at the primary line of the next read. If RNEXT is ‘*’, no
  assumptions can be made on PNEXT and bit 0x20

So, you want to remove those lines whose 7th field isn't = and, just in case, those lines whose 7th field isn't = and isn't the same as the 3rd field. You can therefore use something like this:
samtools view -L test.bed test.bam | awk '$7=="=" || $3==$7

And, to save as a bam file again:
samtools view -L test.bed test.bam | awk '$7=="=" && $3==$7 |
    samtolls view -b > fixed.bam

On a separate note, there's very rarely a need to chain multiple grep commands like that. You can just use \| (or | with the -E or -P options) to separate them. Something like:
samtools view -L test.bed test.bam | grep -v 'chrx\|chr2\|chr10\|chrN'

Or
samtools view -L test.bed test.bam | grep -Ev 'chrx|chr2|chr10|chrN'

If you do it like this, the fixed.bam file is missing the header, which in my experience generates a lot of problems. I recommend always adding the header back in; either by specifying -h when reading the original BAM, or by adding it separately: (samtools view -H infile.bam; samtools view …) > samtools view -b > outfile.bam.

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


text:
Merge hundreds of small BAM files into a single BAM file
I am working with over a million (long) reads, and aligning them to a large genome. I am considering running my alignment jobs in parallel, distributing horizontally across hundreds of nodes rather than trying to run a single job with dozens of cores.
I would like to merge the sorted BAM files together for further downstream analysis. What is the most efficient way to do so while maintaining a valid file header and taking advantage of the fact that the input bam files are already sorted?

samtools merge merged.bam *.bam is efficient enough since the input files are sorted. You can get a bit faster with sambamba and/or biobambam, but they're not typically already installed and IO quickly becomes a bottleneck anyway.

Any idea of how this compares to bamtools merge?

Samtools tends to be a bit slower than other tools since it's the reference implementation, so I'd expect bamtools to be a bit faster.

Bamtools is usually slower because it is not optimized and now unmaintained. Sambamba etc are developed with a specific goal to be faster.

Merging sorted files is a linear operation, so any well-implemented tools that do it will do it with approximately the same efficiency.  So samtools merge (use the most up-to-date version, as there have been improvements in merge header handling in the 1.3.x and 1.4.x versions), picard MergeSamFiles, etc.
These tools need to hold all the input BAM files open simultaneously, so depending on how many hundred of input files you have you may run into the system limit on open file descriptors.  Use ulimit to maximise this first; if there are still too many, you may need to merge the first 500 files, then merge the next 500 into that, etc.  Samtools does not do this internally; I'm not sure whether any of the other merge implementations do.

Maybe merging stepwise as you mention in your second paragraph could speed up things as well as one doesn't hit IO limitations?

No, some of the data is now being written and read more than once, so it must be slower.  With crafty use of pipes and process substitution (samtools merge … <(samtools merge …) …) this can be minimised and avoid hitting disk, so the slowdown can be trivial (at the cost of using a few more cores).


text:
How can I use annotations to remove variants not relevant to cancer risk?
I currently have ~180 whole germlines and around 10M SNPs/indels. I would like to build a predictive model using Machine Learning (ML) techniques to predict cancer risk according to these germline variants. The thing is, most of these 10M variants are not relevant to cancer, so first I have to annotate them and remove non-relevant SNPs/indels.
Which tools/measurements do you think are important to use in order to filter out irrelevant variants and keep only, let's say, around 10,000 SNPs/indels?
I guess the most obvious is CADD score, so I could keep those SNPs/indels with PHRED >= 10.
Which one would you suggest taking into account that downstream analysis would be ML algorithms that take as features variants that really contribute to cancer?

While your question is specific to cancerous germline mutations, I'd suggest you look at the COSMIC database of somatic mutations to include in your analysis.
There are other factors to include in this kind of analysis you're suggesting, such as predictive deleterious effects (PolyPhen for example can perform such predictions).
If you have 10M variants/indels, for the same form of cancer, then look for common variants, or maybe the frequency of variants identified across exons of a gene.
What you've asked is highly broad and we'd need more specific details about your cohort to suggest more specific approaches.

I don't think it will be possible to do what you ask, right now with current knowledge. Selecting variants relevant to cancer risk is still an open problem and usually requires quite a lot of human intervention.
You can use different measures of population frequency to filter common variants with the assumption that frequent variants won't be pathogenic.
You can use multiple in-silico predictors of pathogenicity to try to remove "benign" or "tolerated" variants, but predictors are far from perfect and you have to take into account that many of them are essentially measuring the same, so when combining their outputs you may be biasing your analysis towards nucleotide conservation, etc...
You can use in-silico predictors of splicing effects. These are even farther from perfect and return very noisy data.
You can take into account where the variants lie in the gene. In many genes, variants in the last few nucleotides won't be pathogenic, but may be important in others. In some genes only variants in some domains have an effect, while not in others.
You have to take into account if only one or both alleles of certain genes are affected. For that you may need phasing information. And obviously the list of genes.
Some genes have an impact on cancer risk when they are "deactivated" (or broken somehow) (tumor supressors) while others have an effect when their function is improved (oncogenes).
And for many many genes you'll simply not know if they have any effect in cancer risk at all.
You can get data from different databases and annotate your genomes to remove many variants but in the end you won't have very high quality data to enter to your ML algorithms. Maybe this will work for your project, but take into account that this "first step" of data cleaning is not easy or obvious with current biological knowledge.


text:
How to make a distinction between the "classical" de Bruijn graph and the one described in NGS papers?
In Computer Science a De Bruijn graph has (1) m^n vertices representing all possible  sequences of length n over m symbols, and (2) directed edges connecting nodes that differ by a shift of n-1 elements (the successor having the new element at the right).
However in Bioinformatics while condition (2) is preserved, what is called a De Bruijn graph doesn't seem to respect condition (1). In some cases the graph doesn't look anything like a de Bruijn graph at all (e.g. http://genome.cshlp.org/content/18/5/821.full).
So my question is, if I want to make it explicit that I am using the Bioinformatics interpretation of a de Bruijn graph, is there a term for it? Something like "simplified de Bruijn graph", "projection of a de Bruijn graph", or "graph of neighbouring k-mers"? Are there any papers making this distinction, or did I get it all wrong?

Basically the condition 1 means that even edge-less vertices should be present in the graph, right?

I mean, I wonder if any non-bioinformatics implementation of De Bruijn graph actually store them, since they do not carry any useful information.

There is one more difference in De Brujin graphs used for genome assembly - edges are weighted.

Hi @Slim re. Q1, I believe de Bruijn graphs are connected (one component). You can build them just by providing m and n (http://mathworld.wolfram.com/deBruijnGraph.html).  Q2: yes, implementations don't need all nodes; de Bruijn graph is an abstract entity, a combinatorial structure, like a "complete graph".  But if my very important graph miss some edges (b/c useless) I can't call it "complete". It doesn't make it less important BTW! Q3: that's true! Thanks for editing the question.

Several papers have made this distinction, and a few indeed use different terms to distinguish between them. For example, Kazaux et al. (2016) acknowledge that:

These constraints favour the use of a version of the de Bruijn Graph (dBG) dedicated to genome assembly – a version which differs from the combinatorial structure invented by N.G. de Bruijn.

Kingsford et al. (2010) also recognise the distinction:

Note that this definition of a de Bruijn graph differs from the traditional definition described in the mathematical literature in the 1940s that requires the graph to contain all length-k strings that can be formed from an alphabet (rather than just those strings present in the genome).

The oldest reference I found for a specific term to refer to the assembly-related structure is Skiena and Sundaram (1995), where they call it a subgraph of the de Bruijn digraph. Later, in 2002, Błażewicz et al. will refer to it as a de Bruijn induced subgraph.  The term de Bruijn subgraph is also formally defined in Quitzau’s thesis (2009). There, and also in the article (Quitzau and Stoye, 2008) the authors describe the sequence graph as a modification of the sparse de Bruijn subgraph (commonly used in assembly problems), where non-branching paths are replaced by a single vertex. The term sparse de Bruijn graph is also used by Chauve et al. (2013).
Another term that I found was word graph, described by both Malde et al. (2005) and by Heath and Pati (2007) as a subgraph or as a generalization of a de Bruijn graph. Rødland (2013) summarises some of the terms used for this data structure:

The data structure is best understood in terms of the de Bruijn subgraph representation of S[k]. (...) Some authors may refer to this as a word graph, or even just a de Bruijn graph.

Although we can recognise that the distinction is not very relevant, the question is asking specifically for the situation where one wants to make such a distinction.

As many papers and myself said, assembly de Bruijn graph is just a subgraph of the full de Bruijn graph. Anyone saying differently fails to acknowledge this simple relationship. "Sequence graph" is too general and used in other context (e.g. sequence assembly graph). "Sparse de Bruijn graph" is more appropriate for a graph constructed by skipping some k-mers in reads (e.g. in sparse assembler). Directed acyclic word graph (DAWG) is a pre-existing concept, at least dated back to the 80's, which makes "word graph" ambiguous, too. People should stop inventing new names for a subgraph.

Pevzner did seminal work in using de Bruijn graphs in assembly (http://www.pnas.org/content/98/17/9748.full) and alternative splicing (https://www.ncbi.nlm.nih.gov/pubmed/12169546)

In addition to the regular De Bruijn graph as depicted on the wikipedia, some implementations in bioinformatics feature additional processing. I guess the main reason figure 1 in the paper you linked (concerning the Velvet genome assembler) is slightly different is that a node represents a series of overlapping k-mers. In order to visualize this as a more classic De Bruin graph you would have to connect the k-mers depicted above the nodes. The caption next to figure one describes the processing quite clearly.
As per your last question: I don't think there is a 'Bioinformatic interpretation of a De Bruijn graph'. There are different implementations, which all have there specifics. Thus it would be best to refer to the actual implementation.
As an example: this is a nice paper on how to construct a pan-genome De Bruijn graph of multiple genomes simultaneously.

But an "implementation" of a de Bruijn graph that doesn't include all k-mers is not a de Bruijn graph (in the original sense) any more, right? If the implementation doesn't satisfy condition (1) above, I wonder if there is another name (or a qualifier) being used.

I am quite sure that all original k-mers are present in some form.

Let's first assume DNA only has one strand. An assembly de Bruijn graph is a subgraph of a complete de Bruijn graph. It contains a vertex u if u is a k-mer in reads; it contains an edge u->v, if u and v are adjacent k-mers on a read. Alternatively, we note that an edge u->v is represented by a (k+1)-mer. An assembly de Bruijn graph can be considered a subgraph edge induced from all (k+1)-mers in reads – in fact, some assemblers take the list of (k+1)-mer as a succinct representation of de Bruijn graphs.
DNA has two strands. We just need to induce an assembly de Bruijn graph from all (k+1)-mers and their reverse complement. It is still a subgraph of a complete de Bruijn graph.
Because an assembly de Bruijn graph is just a subgraph. It is not necessary to give it a new name.
PS: I deleted my old answer as that was not what you are asking for based on your comments. I was confused by your mentioning velvet. Velvet uses an equivalent but uncommon representation of de Bruijn graphs, which complicates your question.


text:
Are fgsea and Broad Institute GSEA equivalent?
Several gene set enrichment methods are available, the most famous/popular is the Broad Institute tool. Many other tools are available (See for example the biocView of GSE which list 82 different packages). There are several parameters in consideration :

the statistic used to order the genes,
if it competitive or self-contained,
if it is supervised or not,
and how is the enrichment score calculated.

I am using the fgsea - Fast Gene Set Enrichment Analysis package to calculate the enrichment scores and someone told me that the numbers are different from the ones on the Broad Institute despite all the other parameters being equivalent.
Are these two methods (fgsea and Broad Institute GSEA) equivalent to calculate the enrichment score?
I looked to the algorithms of both papers, and they seem fairly similar, but I don't know if in real datasets they are equivalent or not.
Is there any article reviewing and comparing how does the enrichment score method affect to the result?

According to the FGSEA preprint:

We ran reference GSEA with default parameters. The permutation number
  was set to 1000, which means that for each input gene set 1000
  independent samples were generated. The run took 100 seconds and
  resulted in 79 gene sets with GSEA-adjusted FDR q-value of less than
  10−2. All significant gene sets were in a positive mode. First, to get
  a similar nominal p-values accuracy we ran FGSEA algorithm on 1000
  permutations. This took 2 seconds, but resulted in no significant hits
  due after multiple testing correction (with FRD ≤ 1%).

Thus, FGSEA and GSEA are not identical.
And again in the conclusion:

Consequently, gene sets can be ranked more precisely in the results
  and, which is even more important, standard multiple testing
  correction methods can be applied instead of approximate ones as in
  [GSEA].

The author argues that FGSEA is more accurate, so it can't be equivalent.
If you are interested specifically in the enrichment score, that was addressed by the author in the preprint comments:

Values of enrichment scores and normalized enrichment scores are the
  same for both broad version and fgsea.

So that part seems to be the same.

Yes, the preprint talk much about how the p-values are differently calculated, but I was asking about if the enrichment scores is the same or not, I was aware of the different ways of calculating the p-values. I find the answer by the author himself

Sorry. I didn't realize that you were asking specifically about the enrichment score. You may want to edit the question to make that more obvious.

I found that the author responded that in the discussion in the preprint,

I wonder is the enrichment score calculated the same way as in broad gsea?
Values of enrichment scores and normalized enrichment scores are the same for both broad version and fgsea.


text:
Parallel processing of scripts that use obitools
I have an obitools script (de Barba et al. 2016) that I would like to run faster. How would you run it in parallel to cut down on time?
illuminapairedend -r rawdata_scandinavia_R2.fastq rawdata_scandinavia_R1.fastq | tee rawdata_scandinavia.fastq | obiannotate -S goodAli:'"Alignement" if score>40.00 else "Bad"' | obisplit -t goodAli -p rawdata_scandinavia.
touch rawdata_scandinavia.Bad.fastq
touch rawdata_scandinavia.Alignement.fastq
touch rawdata_scandinavia.fastq
ngsfilter -t rawdata_scandinavia.ngsfilter -u rawdata_scandinavia.unidentified.fastq rawdata_scandinavia.Alignement.fastq > rawdata_scandinavia.filtered.fastq
obisplit -p MICROSAT.PCR_ -t experiment rawdata_scandinavia.filtered.fastq

Can someone describe what tools and steps one needs to take to run this job in parallel, given that the functions are not designed to be run in parallel out of the box?
De Barba, M., Miquel, C., Lobréaux, S., Quenette, P. Y., Swenson, J. E., & Taberlet, P. (2016). High-throughput microsatellite genotyping in ecology: improved accuracy, efficiency, standardization and success with low-quantity and degraded DNA. Molecular Ecology Resources, 1–16. https://doi.org/10.1111/1755-0998.12594

I think this is a more general programming question and so would be more suited for stack overflow

Divide & conquer!

Bioinformatics encompasses biology, statistics, and computer science. Parallel processing is a common need in bioinformatics, and can be a challenging task to implement depending on the particular programs that are used. The programs used here in the example script are specific to bioinformatics, and answers in stack overflow might be lower quality due to the inexperience of the community in working with bioinformatics programs.

@gringer this question was downvoted and closed because it's a programming question. Fast forward a few weeks later, they're suddenly fine.

Pure programming questions (as was the impression from reading only the question title for this) are better placed in SO. However, this question is about a specific type of parallel processing that is associated with a particular bioinformatic tool. I can imagine that this question might be closed in SO for offtopic reasons.

@gringer I hope we can get one more vote. I've solved the problem in the mean time and would like to offer my fruits of labor to anybody who might come up with this kind of a problem in the future.

@RomanLuštrik the FASTA/FASTQ question has now been updated with more context.

Your only option here is to first split the fastq files into small chunks, process each and then merge them back together. Alternatively, find something else that runs faster.

Here is a bash script and GNU parallel along with some other tools to split-apply-combine the data. It may or may not work out of the box, but it's not far off and it should give the general idea of how one way of how to approach this.
#!/bin/bash

# ====== USER PROVIDED PARAMETERS ========
N=20 # number of cores to use
ngsfilter="rawdata_scandinavia.ngsfilter"
unidentified="unidentified_samples.fastq"

R1="rawdata_scandinavia_R1.fastq"
R2="rawdata_scandinavia_R2.fastq"
# ==== END USER PROVIDED PARAMETERS =====

# fastqutils available @ http://ngsutils.org/modules/fastqutils
fastqutils split $R1 to_illumina_R1_ $N &
fastqutils split $R2 to_illumina_R2_ $N &
wait

R1=$( ls | grep -P "to_illumina_R1_" )
R2=$( ls | grep -P "to_illumina_R2_" )

# Do alignment of two reads.
for i in $( seq $N )
do
    partR1=$( ls | grep -P "^to_illumina_R1_\\.$i\\.fastq$" )
    partR2=$( ls | grep -P "^to_illumina_R2_\\.$i\\.fastq$" )
    illuminapairedend -r $partR1 $partR2 | tee to_ngsfilter_$i.fastq | obiannotate -S goodAli:'"Align" if score>40.00 else "Bad"' | obisplit -t goodAli -p to_ngsfilter_$i. &
done

# Remove intermediate aligned files.
echo $R1 | xargs rm -r
echo $R2 | xargs rm -r
ls | grep -P "\\.seq$|\\.err$" | xargs rm -r

input=$( ls | grep -P "to_ngsfilter_[0-9]+\\.Align\\.fastq" )

# For some reasons, files need to be "touched" in order to work. Perhaps
# it's not closing them...
touch $unidentified

touchedbyangel=$( ls | grep -P "to_ngsfilter_"  )
for i in $touchedbyangel
do
    touch $i
done

parallel -j$N --result filtered_{} ngsfilter -t $ngsfilter -u $unidentified {} ::: $input

# Remove intermediate results.
ls | grep -P "\\.seq$|\\.err$" | xargs rm -r

ls | grep -P "filtered_.*\\.fastq$" | xargs cat > filtered_data.fastq

# remove intermediate results to conserve space
rm to_ngsfilter*

# Split reads by locus.
filtered=$( ls | grep -P "^filtered_"  )
parallel -j$N --result to_split_{} obisplit -p MS.PCR_ -t experiment {} ::: $filtered

ls | grep -P "\\.err$|\\.seq$" | xargs rm -r
ls | grep -P "^to_split_" | xargs rm -r

# Wasn't able to parallelize this, so it's just pushing tasks into the background.
obiuniq -m sample MS.PCR_UA_MxRout1_03.fastq > MS.PCR_UA_MxRout1_03.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_06.fastq > MS.PCR_UA_MxRout1_06.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_14.fastq > MS.PCR_UA_MxRout1_14.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_16.fastq > MS.PCR_UA_MxRout1_16.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_17.fastq > MS.PCR_UA_MxRout1_17.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_25.fastq > MS.PCR_UA_MxRout1_25.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_51.fastq > MS.PCR_UA_MxRout1_51.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_57.fastq > MS.PCR_UA_MxRout1_57.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_63.fastq > MS.PCR_UA_MxRout1_63.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_64.fastq > MS.PCR_UA_MxRout1_64.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_65.fastq > MS.PCR_UA_MxRout1_65.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_67.fastq > MS.PCR_UA_MxRout1_67.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_68.fastq > MS.PCR_UA_MxRout1_68.uniq.fasta &
obiuniq -m sample MS.PCR_UA_MxRout1_ZF.fastq > MS.PCR_UA_MxRout1_ZF.uniq.fasta &
wait

uniqfiles=$( ls | grep -P ".{2}+\\.uniq\\.fasta$" )

for i in $uniqfiles
do
    # this section removes the file suffix
    # https://stackoverflow.com/questions/125281/how-do-i-remove-the-file-suffix-and-path-portion-from-a-path-string-in-bash
    tab=${i%.fasta}
    tab=${tab##*/}
    obigrep -p 'count>1' $i | obiannotate -k merged_sample -k count | obiannotate --length | obisort -r -k seq_length | obitab --no-definition --output-seq > $tab.tab
done

ls | grep -P "\\.uniq\\.fasta" | xargs rm -r

echo "Done. I think."


text:
Understanding DESeq2 design, contrast and results
I have a set of high-troughput experiments with 2 genotypes ("WT" and "prg1") and 3 treatments ("RT", "HS30" and "HS30RT120"), and there are 2 replicates for each of the genotype x treatment combinations.
The read counts for the genes are summarized in a file that I load as follows in R:
> counts_data <- read.table("path/to/my/file", header=TRUE, row.names="gene")
> colnames(counts_data)
 [1] "WT_RT_1"          "WT_HS30_1"        "WT_HS30RT120_1"   "prg1_RT_1"
 [5] "prg1_HS30_1"      "prg1_HS30RT120_1" "WT_RT_2"          "WT_HS30_2"
 [9] "WT_HS30RT120_2"   "prg1_RT_2"        "prg1_HS30_2"      "prg1_HS30RT120_2"

I describe the experiments as follows:
> col_data <- DataFrame(
    geno = c(rep("WT", times=3), rep("prg1", times=3), rep("WT", times=3), rep("prg1", times=3)),
    treat = rep(c("RT", "HS30", "HS30RT120"), times=4),
    rep = c(rep("1", times=6), rep("2", times=6)),
    row.names = colnames(counts_data))
> col_data
DataFrame with 12 rows and 3 columns
                        geno       treat         rep
                 <character> <character> <character>
WT_RT_1                   WT          RT           1
WT_HS30_1                 WT        HS30           1
WT_HS30RT120_1            WT   HS30RT120           1
prg1_RT_1               prg1          RT           1
prg1_HS30_1             prg1        HS30           1
...                      ...         ...         ...
WT_HS30_2                 WT        HS30           2
WT_HS30RT120_2            WT   HS30RT120           2
prg1_RT_2               prg1          RT           2
prg1_HS30_2             prg1        HS30           2
prg1_HS30RT120_2        prg1   HS30RT120           2

I want to build a DESeq2 object that I could use to either:

find differentially expressed genes when the treatment varies for a given fixed genotype

or:

find differentially expressed genes when the genotype varies for a given fixed treatment

In the bioconductor help forum I think I've found a somewhat similar situation, and I read the following:

Try a design of ~ genotype + genotype:condition
Then you will have a condition effect for each level of genotype, including the reference level.
You can constrast pairs of them using the list style of the 'contrast' argument.

However, this doesn't explain how to apply this "list style" to the "contrast" argument. And the above situation seems to be asymmetrical. By that I mean that genotype and condition do not seem to have an interchangeable role.
So I tried the following more symmetric formula:
> dds <- DESeqDataSetFromMatrix(
                                countData = counts_data,
                                colData = col_data,
                                design = ~ geno + treat + geno:treat)
> dds <- DESeq(dds)

Now, can I for instance get the differential expression results when comparing treatment "HS30" against "RT" as a reference, in genotype "prg1"?
And how?
If I understand correctly, the above-mentioned "list style" uses names given by the resultsNames function. In my case, I have the following:
> resultsNames(dds)
[1] "Intercept"               "geno_WT_vs_prg1"
[3] "treat_HS30RT120_vs_HS30" "treat_RT_vs_HS30"
[5] "genoWT.treatHS30RT120"   "genoWT.treatRT"

I guess I would need a contrast between "genoprg1.treatRT" and a "genoprg1.treatHS30", but these are not in the above results names.
I'm lost.

this is a more similar example (I think): https://support.bioconductor.org/p/63201/. Conclusion is similar to Devon's answer - use LRT.

It seems that the "combining factors" trick described in part 3.3 of DESeq2 current "vignette" (as of may 2017) under the title "Interaction" is a way to access to the desired contrasts.
It seems possible do do it directly when building the colData and when calling DESeqDataSetFromMatrix:
Let's add a combined "geno" and "treat" factors to the future colData parameter:
> col_data$geno_treat <- as.factor(paste(col_data$geno, col_data$treat, sep="_"))
> col_data
DataFrame with 12 rows and 4 columns
                        geno       treat         rep     geno_treat
                 <character> <character> <character>       <factor>
WT_RT_1                   WT          RT           1          WT_RT
WT_HS30_1                 WT        HS30           1        WT_HS30
WT_HS30RT120_1            WT   HS30RT120           1   WT_HS30RT120
prg1_RT_1               prg1          RT           1        prg1_RT
prg1_HS30_1             prg1        HS30           1      prg1_HS30
...                      ...         ...         ...            ...
WT_HS30_2                 WT        HS30           2        WT_HS30
WT_HS30RT120_2            WT   HS30RT120           2   WT_HS30RT120
prg1_RT_2               prg1          RT           2        prg1_RT
prg1_HS30_2             prg1        HS30           2      prg1_HS30
prg1_HS30RT120_2        prg1   HS30RT120           2 prg1_HS30RT120

We can now use a design where differential expression will be explained by these combined factors:
> dds <- DESeqDataSetFromMatrix(
    countData = counts_data,
    colData = col_data,
    design = ~ geno_treat)

We run the analysis:
> dds <- DESeq(dds)

Then we can query results for a particular contrast between such factor combinations. For instance, to have the results for the effect of treatment "HS30" against the reference state "RT" in genotype "prg1":
res <- results(dds, contrast=c("geno_treat", "prg1_HS30", "prg1_RT"))

The simplest manner is to not use a wald test, but rather an LRT with a reduced model lacking the factor of interest:
dds = DESeq(dds, test="LRT" reduced=~geno+geno:Treatment)

The above would give you results for Treatment regardless of level while still accounting for a possible interaction (i.e., a "main effect of treatment, regardless of the type of treatment").
As an aside, this is probably a case where the edgeR-preferred way of creating groups of genotype-treatment combinations and then using a model of ~0+group might make your life a bit easier. You'll get the same results (more or less) regardless, but it'll probably be easier for you to think in those terms rather than remembering that the base level will be treatment HS30 and geno prg1.

I tried the following:
dds <- DESeqDataSetFromMatrix(countData = counts_data, colData = col_data, design = ~ geno + treat + geno:treat) then dds = DESeq(dds, test="LRT", reduced=~geno+geno:treat) but this fails with Error in nbinomLRT(object, full = full, reduced = reduced, betaPrior = betaPrior,  : less than one degree of freedom, perhaps full and reduced models are not in the correct order. Besides, what I want is not "regardless of the type of treatment". I want results for a given fixed treatment or for a given fixed genotype. Maybe my question was not clear.

I'd have to go over the model matrix to see why it's producing the error. Anyway, I would suggest going the route the edgeR authors prefer and making the treatment-genotype combinations into individual groups. Then the contrasts are more obvious.

"The simplest manner is to not use a wald test, but rather an LRT" -- could you explain why?

Because treatment has multiple levels, so with a wald test one would need to handle the contrasts appropriately, whereas with an LRT one can simultaneously omit all the levels at once.


text:
Why do some assemblers require an odd-length kmer for the construction of de Bruijn graphs?
Why do some assemblers like SOAPdenovo2 or Velvet require an odd-length k-mer size for the construction of de Bruijn graph, while some other assemblers like ABySS are fine with even-length k-mers?

From the manual of Velvet:

it must be an odd number, to avoid palindromes. If you put in an even
  number, Velvet will just decrement it and proceed.

the palindromes in biology are defined as reverse complementary sequences. The problem of palindromes is explained in this review:

Palindromes induce paths that fold back on themselves. At least one
  assembler avoids these elegantly; Velvet requires K, the length of a
  K-mer, to be odd. An odd-size K-mer cannot match its reverse
  complement.

It is possible to construct graph with palindromes, but then the interpretation will be harder. Allowing only graphs of odd k-mers is just an elegant way to avoid writing a code for interpretation of a more complicated graph.

Lest someone in the future misinterpret this, it should be noted that palindrome in this context has a slightly more specific meaning than it normally would in English.

To expand on the answer above, in case it isn't clear, we show:

Why palindromic sequences must be even in length
Why palindromic sequences induce self-loops in a de Bruijn graph
Why self loops in a de Bruijn graph are problematic

1. Palindromic sequence ⇒ sequence is of even length

Idea: in an odd-length k-mer, its middle nucleotide is 'flipped' in its reverse complement, so the two can never be equal.

Suppose you have a palindromic sequence $X$. Then $X$ is identical to its reverse complement, which we will label $\bar{X}$.
Suppose $X$ is of odd length. Then it is of the form $AbC$, where $len(A) = len(C) = \frac{len(X)-1}{2}$, and $len(b) = 1$.
Then
$X = \bar{X} \implies AbC = \overline{AbC} = \bar{C}\bar{b}\bar{A}$
And hence:
$b = \bar{b}$
(since $len(A) = len(C) = len(\bar{C}))$. But this is a contradiction, since $b$ is a single nucleotide, and cannot be equal to its complement. Therefore k-mers of odd length cannot form palindromes.
Hence the length of a k-mer that forms a palindrome must be even.

2. Why palindromic k-mers induce self loops
Each node in a traditional de Bruijn graph is a unique string, but in most bio-informatics implementations each pair of reverse-complementary k-1-mers are identified as a single node, e.g for $k=6$:

A palindromic k-mer (of $k \geq 2$) is of the form:
$xAy$
where $len(A) = k-2$, $x=\bar{y}$ and $A=\bar{A}$ (possibly the empty string).
Hence it will contribute two nodes in the de Bruijn graph:

its left k-1-mer $xA$
its right k-1-mer $Ay$

And an edge going from 1 to 2.
But since this k-mer is palindromic, $xA = \overline{Ay}$ and hence these two nodes are reverse-complementary, and thus the 'same' node, and so this edge is a self loop at this node.

3. Why are self loops problematic?
Self-loops (if they occur in a node with $in\_degree \geq 2$ and $out\_degree \geq 1$) increase the number of possible Eulerian paths in a de Bruijn graph (or more specifically, in the connected component containing this node, which represents a contig, of which their may be multiple), as you have an additional possible Eulerian path for each time you traverse this node.
This increases ambiguity in reading the graph, as each possible Eulerian path is an extra possible reconstruction of the full sequence.
Consider the example:

There is only one possible Eulerian path:

$A B C D B E$

However, if we include a self loop at $B$, which is visited twice above, this doubles to two possible Eulerian paths:

$A B B C D B E$
$A B C D B B E$

Depending on whether we traverse the self loop during the first time we reach $B$, or the second.

https://homolog.us/Tutorials/book4/p2.4.html "Genome assembly programs also avoid even k, because with even k, many k-mers become reverse complements of their own sequences. That causes ambiguities in the strand-specificness of the graph. Therefore, odd k-values are preferred."

Nice answer @ukemi. It took me a while to understand the conclusion of the point 1, so I have added there a sentence that would have helped me. I you don't like it you can reverse the change, but I would say that a small clarification would be good there.

@KamilSJaron no worries, the clearer the better - yeah, technically I should have also shown existence for the implication that they must be even to follow (as opposed to just not-odd), but showing existence is trivial via example (e.g. AT, ATAT etc).

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


text:
Are soft-clipped bases used for variant calling in samtools + bcftools?
If there are soft clipped base pairs specified in the CIGAR string for a read in a SAM/BAM file, will these be used for variant calling in a samtools + bcftools workflow?
The GATK HaplotypeCaller, for example, has an explicit option --dontUseSoftClippedBases for whether to use soft clipped bases. The samtools documentation does not mention clipped bases.

No, samtools (and therefore bcftools) does not use soft-clipped bases. You can quickly confirm this by using either samtools depth or samtools mpileup to look at a region with a soft-clipped alignment. You'll note that the soft-clipped region isn't used in the depth/pileup (both tools use the same underlying code, so it doesn't matter which you use). If you're curious, samtools ignores soft-clipped bases because it's based on making a per-base stack of alignments covering each position. In the BAM format, alignments are sorted and assigned to bins according to their start/end positions, which won't include soft-clipping. Consequently, when samtools is making the pileup it won't even see the alignments that would overlap a given base if soft-clipped bases were included.
This then sort of begs the question of what GATK's HaplotypeCaller is doing differently. There, regions in the genome are essentially assembled in a small de Bruijn graph, which allows for soft-clipped bases around indels to then be resolved, given that the graph would start/end a little-way on past each side of indels. This is also why you don't need to do indel realignment with the HaplotypeCaller (this was needed in the old UnifiedGenotyper).
Edit: For more details regarding the HaplotypeCaller, see this nice page on GATK's website, which goes into much more detail than I did here.


text:
Phyre2 vs ITasser, completely different models generated
Does anyone have experience generating pdb structures with Phyre and ITasser online tools. The results generated from each given the same amino acid sequence input are very different and I am wondering whether or not this is a usual experience. I know ITasser was the number 1 rated in the CASP trails, but should the results really be this disparate?

You should always pass your models through a model quality assessment program to evaluate reliability, and select the best one.

I'm less familiar with Phyre, but I-TASSER is a really sophisticated system that takes the results of a search using multiple threaders and plugs them into an ab initio simulation which tries to minimize the energy of the models by sampling many possible 3D conformations, which I don't think Phyre does.
https://en.wikipedia.org/wiki/I-TASSER#/media/File:I-TASSER-pipeline.jpg
Compare with a similar workflow schematic for Phyre:
http://www.nature.com/nprot/journal/v10/n6/images/nprot.2015.053-F1.jpg
http://www.nature.com/nprot/journal/v10/n6/fig_tab/nprot.2015.053_F2.html
Structure prediction still has a long way to go, and you'll always get better results if there are close homologues available in the PDB, but given the consistent high performance of I-TASSER in CASP I would treat those results as more significant. That said, it can't hurt to consider multiple answers.
Edit: included blmoore's link to second half of Phyre protocol schematic

You're only linking to the first half of the phyre2 workflow: http://www.nature.com/nprot/journal/v10/n6/fig_tab/nprot.2015.053_F2.html

This may not be the whole answer to the question, but one key difference between Phyre and ITasser (I think I'm correct in saying), is that Phyre only returns regions of the model it considers to be sufficiently well modelled, so it may not be a full length model. ITasser on the other hand, always generates full size models.
As one of the other answers alludes to, it is very dependent on your protein of interest. If it isn't well represented in the databases, the initial threading stages will not yield good models, and the subsequent steps all suffer as a result.


text:
How are MACS2's narrow peak and broad peak algorithms different?
The peak calling tool MACS2 can call peaks in either narrow peak mode (for focused signals like transcription factor ChIPseq) or broad peak mode (for more defuse signals, like certain histone modifications).
The algorithm for narrow peak calling is well described in the MACS publication. But I don't find much documentation for how peak calling is different in broad peak mode. The manual only contains the following:

--broad
When this flag is on, MACS will try to composite broad regions in
  BED12 ( a gene-model-like format ) by putting nearby highly enriched
  regions into a broad region with loose cutoff. The broad region is
  controlled by another cutoff through --broad-cutoff. The maximum
  length of broad region length is 4 times of d from MACS

But this doesn't really describe exactly how this is performed.
So what is the algorithm MACS uses for calling broad peaks?

I found the algorithm, but honestly I could not really get the try/except construction on line 629...

The key function is call_broadpeaks:
The description attached to the function says:

This function try to find enriched regions within which, scores are
  continuously higher than a given cutoff for level 1, and link them
  using the gap above level 2 cutoff with a maximum length of
  lvl2_max_gap.
scoring_function_s: symbols of functions to calculate
  score. 'p' for pscore, 'q' for qscore, 'f' for fold change, 's' for
  subtraction. for example: ['p', 'q']
lvl1_cutoff_s:  list of cutoffs
  at highly enriched regions, corresponding to scoring functions.
lvl2_cutoff_s:  list of cutoffs at less enriched regions,
  corresponding to scoring functions.
min_length :minimum peak length,
  default 200.
lvl1_max_gap   :  maximum gap to merge nearby enriched
  peaks, default 50.
lvl2_max_gap   :  maximum length of linkage
  regions, default 400.
Return both general PeakIO object for highly
  enriched regions and gapped broad regions in BroadPeakIO.

To give some basic explanation, the algorithm (briefly) appears to be as follows:

Two separate levels of peaks are called, level 1 (a higher pval, ie more significant) and level 2 (a lower pval). Level 1 is controlled by -p and level 2 is controlled by --broad-cutoff. When each peakset is called, they are immediately linked by the max gap parameter for each set.
Then, assuming that all level 1 peaks should be inside level 2 peaks (this is an explicit assumption by MACS2), the algorithm groups level 1 peaks inside level 2 peaks to output a broad peak.

...
This has a few implications:

The broad peak calls really come from the level 2 peaks alone (+ linking). The level 1 peak calls allow you to distinguish sub peaks (so that you can have gapped peaks).
Aside from the linking, the broad peak calls would be the same as narrow peak calls, if you called both with the same pval threshold (for example, if you set --broad-cutoff 0.1 in broad peak mode, and the -p 0.1 for narrow peak mode)

Could you clarify in your answer that "higher pval" is actually "lower" or "more significant" p value and that "lower pval" is actually "higher" or "less significant" p value? I foresee someone getting confused by the wording.

Added a clarification - thanks for the catch!


text:
How to retrieve logical expressions (KO based) for reactions from KEGG?
The completeness of a module can easily be checked by looking at the Definition entry associated with the module. For example, in module M00010, it is given as:
Definition  K01647 (K01681,K01682) (K00031,K00030)

which can be translated to:
K01647 AND (K01681 OR K01682) AND (K00031 OR K00030)

If this expression evaluates to TRUE, the module is considered complete.
I am wondering whether analogous information exists for a single reaction. For example, when looking at reaction R00352, one finds the following information about its Orthology:

But I can't easily find the logical relation that correspond to the KOs for this reaction. It could be:
K01648 AND K15230 AND K15231

or:
K01648 OR K15230 OR K15231

or:
K01648 OR (K15230 AND K15231)

and so on.
In the above example, I know that the correct expression is:
K01648 OR (K15230 AND K15231)

i.e. one either needs either K01648 or both of the other two subunits.
Can this information be retrieved from KEGG for each reaction and if so, how?

Have you tried to asked KEGG directly via http://www.kegg.jp/feedback/ ? If anything, that might give them the information that there's something missing which people need.

You might need to look into the FTP dump; I looked into the KEGG REST API (http://www.kegg.jp/kegg/rest/keggapi.html) and couldn't find anything related to your question...

Some information seems to be available on the module level utilising the API.
Reaction R00352 points to the module M00173. While some information is in the Definition, it's not trivial to link it to the particular reaction. This linkage is however easier accessible via API:

Thanks. There are two issues with this as far as I can see: i) the third KO (K01648) is missing in the list and ii) it is not clear which logical relationship exists between these KOs; it seems it is always an AND but it is not entirely clear (e.g. from the link you provide: K00239,K00240,K00241,K00242,K00244,K00245,K00246,K00247 - do I need all of these or not?.

This table seems to point where the indicated members should all be required as one. In a couple of cases I looked it is indeed the case. Hence ad i) K01648 is separate so should not be noted ad ii) a very detailed decomposition to subunits, but yes.


text:
Given a VCF of a human genome, how do I assess the quality against known SNVs?
I'm looking for tools to check the quality of a VCF I have of a human genome. I would like to check the VCF against publicly known variants across other human genomes, e.g. how many SNPs are already in public databases, whether insertions/deletions are at known positions, insertion/deletion length distribution, other SNVs/SVs, etc.? I suspect that there are resources from previous projects to check for known SNPs and InDels by human subpopulations.
What resources exist for this, and how do I do it?

To achieve (at least some of) your goals, I would recommend the Variant Effect Predictor (VEP). It is a flexible tool that provides several types of annotations on an input .vcf file.  I agree that ExAC is the de facto gold standard catalog for human genetic variation in coding regions.  To see the frequency distribution of variants by global subpopulation make sure "ExAC allele frequencies" is checked in addition to the 1000 genomes.
Output in the web-browser:

If you download the annotated .vcf, frequencies will be in the INFO field:
##INFO=<ID=CSQ,Number=.,Type=String,Description="Consequence annotations from Ensembl VEP. Format: Allele|Consequence|IMPACT|SYMBOL|Gene|Feature_type|Feature|BIOTYPE|EXON|INTRON|HGVSc|HGVSp|cDNA_position|CDS_position|Protein_position|Amino_acids|Codons|Existing_variation|DISTANCE|STRAND|FLAGS|SYMBOL_SOURCE|HGNC_ID|TSL|SIFT|PolyPhen|AF|AFR_AF|AMR_AF|EAS_AF|EUR_AF|SAS_AF|AA_AF|EA_AF|ExAC_AF|ExAC_Adj_AF|ExAC_AFR_AF|ExAC_AMR_AF|ExAC_EAS_AF|ExAC_FIN_AF|ExAC_NFE_AF|ExAC_OTH_AF|ExAC_SAS_AF|CLIN_SIG|SOMATIC|PHENO|MOTIF_NAME|MOTIF_POS|HIGH_INF_POS|MOTIF_SCORE_CHANGE

The previously mentioned Annovar can also annotate with ExAC allele frequencies.  Finally, should mention the newest whole-genome resource, gnomAD.

Very good suggestion. My understanding is that ExAC holds approx. 123K exomes and  gnomAD now has around 15K genomes

These are all great answers, but as this seems to be the community favorite, I will mark is as "the answer". But future readers: please look at the others!

The greatest protein coding variant catalogue is definitely ExAC (>65k individuals). They also published a blogpost where they describe how to reproduce figures in the paper (it is a good start how to get familiar with the dataset).
For the whole-genome variants I would look at the data created by 1000 genomes project (the latest release has more than 3k individuals). The integrated variant call sets can be downloaded though the portal and catalogue of SVs can be found here.
In this paper (also 1000 genomes project) they speak about non-precise placement of SVs by SV callers. I would keep this in mind for the comparison of your genome to the known variants.

Thanks for the help. I'm a bit confused by the hyperlink at guide which describes the perl script needed. Using this script, this will create a VCF from an aggregated ~2.5K individuals?

The script is part of VCFtools, they link it to manual of the package instead of the source code, you are expected to install the package, not just to download the script. I also found thay have started a portal on the welpage so I edited the answer.

Your best bet is to use programs that provide you an complete annotation of variants present in your VCF. Two examples are snpEff and Annovar. These programs work on known variants deem different sources and provide you with information on each item in your file, that you can filter after to try to understand the effects of each variant.

There is now access to gnomAD with ANNOVAR, which is pretty cool.


text:
How to deal with heterozygosity during polishing of genome assembly based on long reads?
All the long-read sequencing platforms are based on single-molecule sequencing which causes higher per-base error rates. For this reason a polishing step was added to genome assembly pipelines - mapping raw reads back to assembly and correcting details of the assembly.
I have decent PacBio RSII dataset of single individual genome of heavily heterozygous non-model species. Assembly went well, but when I tried to polish the assembly using quiver it could not converge over a couple of iterations and I bet it is because of too great divergence of haplotypes.
Is there any other way to polish a genome with such properties?
For instance, is there a way to separate long reads by haplotype, so I could polish using one haplotype only?

A few possibilities:
Falcon
Try falcon and falcon-unzip. These are designed exactly for your problem and your data: https://github.com/PacificBiosciences/FALCON
Not Falcon
If you think you have assembled haplotypes (which seems reasonable to expect given enough coverage), you should be able to see the two haplotypes by just doing all pairwise alignments of your contigs. Haplotypes should show up as pairs of contigs that are MUCH more similar (even with a lot of between-haplotype divergence) than other pairs. Once you have all such pairs, you can simply select one of each pair to polish.

I indeed do have both haplotype sequences. I got them using tool called haplomerger2. But this tool produces a chimeric haploid assembly, therefore they are not really correctly phased haplotypes.
Falcon-unzip is indeed software that could work. It was too young to try at the time, but I could try to give it another shot now.

You could also have a go at Canu. It's designed for long-read assembly (both PacBio and Nanopore), although not specifically for complex population sequencing. It tries to strip a genome down into its unique components, and generates paths from those components that are well-supported from the reads.
With regards to polishing, it seems to be the case that polishing doesn't converge, and there will be a lot of variants that just oscillate between two possibilities. For me and at least one other person at London Calling this year, there was basically no gain in accuracy for polishing past the third iteration. I used my own error correction algorithm, but they used the more "standard" polishing with Pilon. For what it's worth, the nanopore WGS consortium used Racon for polishing their Canu assemblies.

I actually have assembled the genome using Canu, I got ~2x haploid size of the genome, which I collapsed to haplotypes the using HaploMerger2.
I know that globally the assembly is good. It just need to be polished.

Oh, yes. Sorry, I looked at the first answer and had assumed that this was just about assembly. I realise now that the question was discussing polishing, rather than assembly.

@gringer I was also trying to polish a highly heterozygous genome assembly (generated by canu), using Racon (Quiver would collapse haplotypes), but couldn't get a satisfying output (basically, no statistics has changed). any advice?

My general recommendation at the moment would be to use nanopolish in methylation mode to correct, then Pilon with Illumina reads to only correct the homopolymer fragments (i.e. no SNP correction, and no long-range scaffolding). Based on this:
https://github.com/rrwick/Basecalling-comparison#methylation

There is a new tool called Hapo-G for diploid-aware polishing. Unfortunately, it seems it's designed for short reads only. It seems to be a more practical choice with respect to inferring the phase information on individual reads.


text:
Expected allele frequency distribution of SNVs in real NGS data
I have a huge amount of ~20x human WGS samples, aligned, and all SNVs that were called with GATK under standard germline parameters set.
What I need to do is to model SNVs Allele Frequency (AF) for different underlying Copy Numbers. I'd better provide a toy example. For particular genomic region X:
If X is presented by 2 copies for the particular samples, we expect AF to be super-close to 1 or to 0.5.
If X is presented by 4 copies, I expect any particular AF to be close to 0.25, 0.5, 0.75 or 1.
Of course, I can use Binomial Distribution for these purposes. However, as we know, the distribution is not exactly Binomial due to alignment/sequencing biases and the median AF for all heterozygous SNVs is more close to 0.48 but not to 0.5 as we would expect. Another thing: for high copy numbers we expect higher coverages. And GATK use several filters so I suppose that we will not see SNVs with AF like 0.125 (in case if the segment has ploidy 8) - despite the super high coverage there GATK may reject this "weird" AF.
I have read several papers that model SNVs AFs (and I agree that Beta Binomial Distribution may be quite accurate), however, I was not convinced enough that I should use the particular modelling. From your experience (in case if you do SNVs calling), which probabilistic distribution should I use? How should I estimate parameters for each of them (should I expect for CN4 AF=0.5 more frequent than AF=0.75 or vice versa, how to estimate this from data)?
UPD: For simplicity we can say that we have a lot of previously identified regions with ploidy different from CN2, and I can take these coordinates from here. So I can use more or less "supervised" learning for parameters' estimation.

I don't have enough experience to answer which probabilistic distribution should be used.
However, this questions also also asks how to estimate parameters of the distributions. If a binomial distribution is chosen, then Heng Li's paper titled "A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data" 1 is probably the definitive one. Section 2.3.1 of that paper describes an EM algorithm for estimating allele frequencies from multiple samples under the assumption of Hardy-Weinberg equilibrium for arbitrary but constant ploidy.

Ideally even most popular tools like GATK also uses EM algorithm for the record to estimate the likelihood while samtools uses both EM and Brent's method. The real problems with the setting up of ploidy apriori but am not aware as of now about tools that estimates ploidy and then uses that information for SNV calls and finds out the AF. I know tools like ABSOLUTE have the power to estimate the ploidy but then you have to use them as an input for your later SNP calls. But I reckon most tools usually works with the assumption with ploidy set apriori. This is what is coming to my mind as of now.


text:
How can I extract normalized read count values from DESeq2 results?
The results obtained by running the results command from DESeq2 contain a "baseMean" column, which I assume is the mean across samples of the normalized counts for a given gene.
How can I access the normalized counts proper?
I tried the following (continuing with the example used here):
> dds <- DESeqDataSetFromMatrix(countData = counts_data, colData = col_data, design = ~ geno_treat)
> dds <- DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
> res <- results(dds, contrast=c("geno_treat", "prg1_HS30", "prg1_RT"))

Here is what I have for the first gene:
> res["WBGene00000001",]$baseMean
[1] 181.7862
> mean(assays(dds)$mu["WBGene00000001",])
[1] 231.4634
> mean(assays(dds)$counts["WBGene00000001",])
[1] 232.0833

assays(dds)$counts corresponds to the raw counts. assays(dds)$mu seems to be a transformation of these counts approximately preserving their mean, but this mean is very different from the "baseMean" value, so these are likely not the normalized values.

The normalized counts themselves can be accessed with counts(dds, normalized=T).
Now as to what the baseMean actually means, that will depend upon whether an "expanded model matrix" is in use or not. Given your previous question, we can see that geno_treat has a bunch of levels, which means that expanded models are not in use. In such cases, the baseMean should be the mean of the base factor in geno_treat.

mean(assays(dds)$mu["WBGene00000001",] / sizeFactors(dds)) gives something very close to res["WBGene00000001",]$baseMean (181.9139). The vignette says that assays(dds)$mu correspond to s_j * q_ij where, according to the DESeq2 paper, s_j is the size factor for library j and q_ij is, if I understand correctly, an estimation of a quantity "proportional to the concentration of cDNA fragments from the gene in the sample". Is this q_ij what you call "base factor" ?

Not really. With an "expanded model matrix" the intercept (aka, "baseMean") should be approximately the average of the normalized counts. For non-"expanded model matrices" then one level of the group factor (the "base level") becomes the intercept, so the baseMean is then just the average of its counts. The base level is the alphanumerically first level in R by default, if you have grp = factor(c("WT", "MUT")) then the MUT group is used to calculate the baseMean (assuming no expanded model matrices).

Here, I find that the mean of all normalized values (mean(counts(dds, normalized=T)["WBGene00000001",])), is actually equal to the baseMean (res["WBGene00000001",]$baseMean). This is not the case if I compute the mean only for a given level of the geno_treat factor.

I guess it using an "expanded model matrix" then. I know Michael used to code in exactly the cases when that would happen, but perhaps he's since generalized the code, which is nice :)

It depends what you mean by “normalised”. As Devon said, the normalized = TRUE argument to the count function gives you normalised counts. However, these are “only” library-size normalised (i.e. divided by the sizeFactors(dds)).
However, as the vignette explains, downstream processing generally requires more advanced normalisation, to account for the heteroscedasticity of the counts. This is often done by simply logging the counts but this has obvious drawbacks (most trivially, what do we do with 0 counts? A workaround is to add a pseudocount but that’s problematic too).
DESeq2 offers two different methods to perform a more rigorous analysis:

rlog — a regularised log, and
vst — a variance stabilising transformation.

You’d generally use either of these for downstream analysis, not count(dds, normalized = TRUE).

to make my confusion clear let's say im going for any parametric or non parametric test between any comparison for genes ,can i take the rlog or vst ? As rlog and vst i take it for clustering ,pca  ,correlation ,wgcna .So is it possible ?

@krushnachChandra It depends on the test and the comparison. DESeq for instance does not use either rlog or vst for its test of differential expression, as these would be inappropriate.

yes that i read in the paper , for example i got a cluster and i wanted to see what if the genes are having significant difference in their expression so if i have to run say a non parametric test can i use those rlog ? values or i have to use the count normalised ones ?

@krushnachChandra You cannot use count normalised values in this scenario, nor simple rlog values, since they are not normalised for transcript length (and thus longer transcripts will show higher expression on average). You could use something like size-factor normalised TPM instead. You can also log those TPMs but for a nonparametric test it would make no difference — that’s kind of the point of non-parametric tests.

So for PCA ,clustering i can use rlog i hope that part im doing correct let me know, i read your post , so can i use the deseq2 value and convert them into TPM ,but there is one part where even deseq2 also calculate size factor "dds <- estimateSizeFactors(dds)" this i guess . So can i use your TPM conversion method on deseq2 values? or just counts .

im curious about this paper https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4986529/ Figure 6 last box plot "(Right) Mean expression of the two gene sets in each population. Error bars represent SD. CPM, counts per million" where they use cpm values i  guess these are edge R output ,so here they used that value to make comparison ,so i was wondering if it can be done with rlog values ?sorry for mundane question my concept about these are very fuzzy as im not a statistician ..


text:
How to select a power for a scale-free topology network
In a weighted gene co-expression network analysis (using WGCNA), the soft-threshold power is recommended as a noise filtering. It consists on raising the correlation to a certain number. To decide this power the scale-free topology is estimated for some powers. The function to estimate this prints:
   Power SFT.R.sq  slope truncated.R.sq mean.k. median.k. max.k.
1      1   0.9300  3.110          0.996  2960.0    3060.0   3970
2      2   0.7510  1.010          0.964  1750.0    1780.0   2900
3      3   0.1730  0.258          0.806  1170.0    1150.0   2280
4      4   0.0942 -0.183          0.713   833.0     782.0   1870
5      5   0.3800 -0.463          0.777   623.0     559.0   1580
6      6   0.5350 -0.656          0.834   481.0     412.0   1360
7      7   0.6270 -0.797          0.872   381.0     312.0   1190
8      8   0.6870 -0.910          0.900   307.0     241.0   1050
9      9   0.7270 -1.000          0.918   252.0     189.0    936
10    10   0.7490 -1.080          0.928   210.0     150.0    841
11    12   0.7850 -1.190          0.948   150.0      98.0    693
12    14   0.8090 -1.280          0.958   111.0      65.9    582
13    16   0.8290 -1.360          0.968    84.0      45.6    497
14    18   0.8410 -1.410          0.973    65.2      32.1    429
15    20   0.8490 -1.450          0.977    51.6      23.0    375

The recommendations of the FAQ indicate that a SFT.R.sq value should be above

0.8 for reasonable powers (less than 15 for unsigned or signed hybrid networks, and less than 30 for signed networks)

and the mean connectivity below the hundreds.
Others have used the power just as noise filtering without caring much about the scale-free topology fit. I would pick the first power even if the mean connectivity is in the order of thousands because the scale-free topology fit is pretty high, however the slope is puzzling me.
How should the soft-threshold power be selected?
Based on an example question

I would be concerned about the input data that produced these numbers. This dataset has a really odd shape with a dip at the start:
data.df <- read.table("soft_threshold.txt", header=TRUE);
png("scale-free-194.png");
plot(data.df$SFT.R.sq);
dummy <- dev.off();

Here's what I expect this to look like (from the supplementary information in this paper):

It's easier to select the power/threshold when the plot looks like the second situation. With the data as presented, it doesn't look like any value would be suitable. but if I were forced to pick, then I'd choose 12 because it is the first number above the magic threshold (ignoring the weirdness at the start). However, I'll admit that I don't understand the concepts behind the power, and am just following the recommended instructions, so can't explain why it is what it is.

Yes, I am aware of this particularity, I think it is due to a bug in the WGCNA package, when pickSoftThreshold doesn't use expression data. I probably should have picked a more typical data.  But simply cutting at certain threshold leaves outside the consideration about the mean connectivity part. (Also it seems to treat the same that the slope of the power =1  is positive and the others become negative)


text:
The state, limitations and comparisons of large variant stores
Background: We're increasingly needing some way of storing lots of variant data associated with lots of subjects: think clinical trials and hospital patients, looking for disease-causing or relevant genes. A thousand subjects is where we'd start, there's talk of millions on the horizon. With various genomic medicine initiatives, this is likely a wider need.
The problem: While there's plenty of platforms out there, it's a rapidly evolving field. It's difficult to get a feel for how (and if) they perform and how they line up against each other:

What's scalable and can handle a lot of data? What sort of limits?
What's robust and not a teetering pile of hacked-together components?
What has a large community behind it and is actually used widely?
What makes for easy access and search from another service? (Commandline, REST or software APIs)
What sort of variants they handle?
What sort of parameters can be used in searching?

Solutions I've seen so far:

BigQ: used with i2b2, but its wider use is unclear
OpenCGA: looks the most developed, but I've heard complaints about the size of data it spits out
Using BigQuery over a Google Genomics db: doesn't seem to be a general solution
Gemini: recommended but is it really scalable and accessible from other services?
SciDb: a commercial general db
Quince
LOVD
Adam
Whatever platform DIVAS & RVD run on: which may not be freely available
Several graphical / graph genome solutions: We (and most other people) are probably not dealing with graph genome data at the moment, but is this a possible solution?
Roll your own: Frequently recommended but I'm sceptical this is a plausible solution for a large dataset.

Anyone with experience give a review or high-level guide to this platform space?

My two cents: use MongoDB wrapped in a simple REST framework.  Allows flexible model and queries and should scale into the billions of records on a single node.  Working on a FLOSS project for this at the moment, but it is not production ready yet.

@woemler How is it compared to other approaches? Someone I know tried MongoDB ~5 years ago on 1000g genotypes. He said MongoDB was over 10x slower than bcf2 on parallel queries while having a much larger disk/memory footprint. That said, he was new to MongoDB back then and might not be doing it the optimal way.

@user172818: The newer versions of MongoDB (3.2+) are significantly faster than the versions from several years ago.  I have benchmarked it against other free RDBMS's, and it typically performs as-well-as or better, especially for complex data representations, like variant calls.

Is storing the data more important here, or is processing statistics (using Python, R, etc..) about the data more important?

@macgyver : good observation. The data - it's supposedly people will want to be mining and querying the data, rather than looking at summary stats and analyses.

An epic question. Unfortunately, the short answer is: no, there are no widely used solutions.
For several thousand samples, BCF2, the binary representation of VCF, should work well. I don't see the need of new tools at this scale. For a larger sample size, ExAC people are using spark-based hail. It keeps all per-sample annotations (like GL, GQ and DP) in addition to genotypes. Hail is at least something heavily used in practice, although mostly by a few groups so far.
A simpler problem is to store genotypes only. This is sufficient to the majority of end users. There are better approaches to store and query genotypes. GQT, developed by the Gemini team, enables fast query of samples. It allows you to quickly pull samples under certain genotype configurations. As I remember, GQT is orders of magnitude faster than google genomics API to do PCA. Another tool is BGT. It produces a much smaller file and provides fast and convenient queries over sites. Its paper talks about ~32k whole-genome samples. I am in the camp who believe specialized binary formats like GQT and BGT are faster than solutions built on top of generic databases. I would encourage you to have a look if you only want to query genotypes.
Intel's GenomicDB approaches the problem in a different angle. It does not actually keep a "squared" multi-sample VCF internally. It instead keeps per-sample genotypes/annotations and generates merged VCF on the fly (this is my understanding, which could be wrong). I don't have first-hand experience with GenomicDB, but I think something in this line should be the ultimate solution in the era of 1M samples. I know GATK4 is using it at some step.
As to others in your list, Gemini might not scale that well, I guess. It is partly the reason why they work on GQT. Last time I checked, BigQuery did not query individual genotypes. It only queries over site statistics. Google genomics APIs access individual genotypes, but I doubt it can be performant. Adam is worth trying. I have not tried, though.

+1 for hail, clearly The Right Answer at this point

You can query individual genotypes using BigQuery.  The biggest challenge at this point is having to write your own queries to do analysis.


text:
Designing a lab NGS file database schema
I am the resident Bioinfo Geek in a hospital academic lab that routinely employs NGS as well as CyTOF and other large volume data producing technologies. I am sick of our current "protocol" for metadata collection and association with the final products (miriad excel sheets and a couple poorly designed RedCap DBs).
I want to implement a central structured, controlled datastore that will take care of this. I know that the interface to the technicians how will be inputing the data is crucial to its adoption, but this is not the focus of THIS particular question: Does there exist a schema or schema guidelines for this type of database?
I would rather use a model that has been developed by people who know how to do this well. I know of BioSQL but it seems more geared towards full protein/nucleotide records like those found in uniprot or genbank. That is not what we have here. What I want is something similar to the system touched on in this preprint: http://biorxiv.org/content/early/2017/05/10/136358
Alternatively, can anyone provide links to where I might find relevant guidelines or supply personal advice?

Are you looking to store processed or unprocessed data?  What would an example file format be that you would be trying to capture?

This is mostly for primary data organization: we get 800 BAMs of WES and I want the file location of each BAM associated with metadata like: PROJECT, READ_LENGTH, SAMPLE_NAME, FAMILY_ID, DATA_TYPE, DIAGNOSIS, etc.

Hey Gus, we are also doing the same research and created this question in Biostars, let us know if you found something! https://www.biostars.org/p/350514/

For metadata, I would use a SQL schema something like the following:
CREATE TABLE Project (
    ac TEXT, -- project/Study accession
    PRIMARY KEY (ac)
);
CREATE TABLE Sample ( -- biological sample/biopsy
    ac TEXT,
    PRIMARY KEY (ac)
);
CREATE TABLE AnalysisSample (
    prj_ac TEXT, -- project acccession (Project.ac)
    symbol TEXT, -- a short name unique in the project
    sample_ac TEXT, -- sample accession (Sample.ac)
    PRIMARY KEY (prj_ac, symbol)
);
CREATE TABLE Collection ( -- a BAM file
    ac TEXT, -- collection/alignment file accession
    prj_ac TEXT, -- project accession (Project.ac)
    PRIMARY KEY (ac)
);
CREATE TABLE ReadGroup (
    cl_ac TEXT, -- collection accession (Collection.ac)
    rg_id TEXT, -- @RG-ID
    sample_sym TEXT, -- @RG-SM; matching AnalysisSample.symbol
    PRIMARY KEY (cl_ac, rg_id)
);
CREATE TABLE VariantSet ( -- a VCF file
    ac TEXT, -- VCF file accession
    prj_ac TEXT, -- project accession (Project.ac)
    PRIMARY KEY (ac)
);
CREATE TABLE VariantSample (
    vs_ac TEXT, -- VCF file accession (VariantSet.ac)
    sample_sym TEXT, -- sample symbol in the VCF file; matching AnalysisSample.symbol
    PRIMARY KEY (vs_ac, sample_sym)
);

In the schema, you have Project and biological Sample tables, which are independent of each other at the high level. An AnalysisSample describes a sample used in BAM or VCF and connects Project and biological Sample. Importantly, each AnalysisSample has a symbol unique in a project (see the primary index). This is the symbol on a BAM read group line or on a VCF sample line. A Collection is in effect a BAM/CRAM file. In theory, a BAM file may contain more than one samples (though rare in practice), which is addressed by a separate ReadGroup table. Finally, a VariantSet is a VCF file. VariantSample tells you which samples are included in each VCF file.
This is the skeleton of a full schema. You can add extra fields to appropriate tables (e.g. file path and hg19/hg38/etc to Collection, read length to ReadGroup and family ID to Sample). You also need indices for efficient table joining and perhaps more tables for complex structures (e.g. pedigree).
For the projects I have participated, this schema should work most of time. It is inspired by GA4GH's JSON schema, but my version is in SQL, is simpler and also has a slightly different structure which I think is better.

The Global Alliance for Genomics and Health has been working on the issue of representing sequencing data and metadata for storage and sharing for quite some time, though with mixed results.  They do offer a model and API for storing NGS data in their GitHub repository, but it can be a bit of a pain to get a high-level view.  I am not sure if any better representation of this exists elsewhere.
I can say from personal experience (having built over a dozen genomic databases), there is no ideal data model and storage best practices.  Genomic data comes in many shapes and sizes, and your needs are going to vary from every other organization, so what works for one bioinformatics group won't necessarily work for you.  The best thing to do is design and implement a model that will cover all of the data types in your workflow and downstream analyses you might do with the data and metadata.

I agree that there is no ideal data model that is going to be stable for very long in a quick-moving field like genome informatics. Perhaps a schema-less (NoSQL or some other document-based system, such as MongoDB) database approach would work better? This gives you ultimate flexibility to attach whatever information is relevant to database entries you're adding to your database now, without the need to rebuild the database later if you want to attach more/different information to subsequent database entries.


text:
Tools for simulating Oxford Nanopore reads
Are there any free open source software tools available for simulating Oxford Nanopore reads?

What are you planning on simulating data for? In my opinion, given the wealth of publicly available real datasets, there is little need for simulation except in very specific use cases.

It's often useful to simulate a variant or some other type of genomic artifact to evaluate whether a tool or set of tools can correctly identify the artifact. This is much easier to do on simulated data where you know the exact answer, and is extremely valuable even if the simulated data isn't a perfect reflection of real data.

Simulators designed specifically for Oxford Nanopore:

NanoSim
NanoSim-H
SiLiCO
ReadSim
DeepSimulator

General long read simulators:

Loresim
Loresim 2
FASTQsim
LongISLND

For an exhaustive list of existing read simulators, see page 15 of my thesis, Novel computational techniques for mapping and
classifying Next-Generation Sequencing data.

By chance, just today I've heard of a nanopore read simulator, NanoSim. It is  released under a GPL license. I have never used it, though...

In addition to the already mentioned NanoSim, there is also SiLiCO and ReadSim (although it hasn't been updated in over 2 years, so I am not sure how relevant it is at this point considering how fast the technology is progressing).

The best nanopore read simulators would be associated with the best base-callers. For a base-caller to effectively model the DNA strand, it needs to take into account the expected underlying electrical model together with the associated signal noise (both in the time dimension as well as the amplitude dimension). In theory, it should be possible to reverse the algorithm and generate an electrical signal given an underlying sequence.
Unfortunately, I'm not aware of any tools that attempt to simulate nanopore reads at the electrical level. Any "nanopore read simulator" that concentrates only on base sequence would need to encompass all the possible base-calling software models that exist, which is an impossible task (particularly given how quickly ONT updates their own base callers).


text:
Which reference to use for read mapping for popular model organisms
What is the "best" assembly for the popular model organisms:

human (GRCh37 and GRCh38 are obvious, I'd pick whatever bwakit uses)
mouse (GRCm37/GRCm38, OK)

but what about non-human/mouse ones?

fruit fly
zebrafish
E. coli
any other idea?

I will answer one of the points – E.coli.
TL;DR Bacteria, and in particular E.coli, are highly variable and there is usually no single best assembly. Large scale WGS studies should come with multiple assemblies for individual monophyletic clusters.
Long answer: Whereas a single reference sequence can make sense for human or mice (to a certain extent), bacteria are not the case. The variability between them even within a single species can be enormous and various lineages can contain very different genes. This is the reason why biologists usually distinguish a core genome consisting of the genes shared among a vast majority (typically 95%) of individuals within a given phylogenetic clade, (i.e., of those genes which are essential for this clade – imagine a CPU and RAM in a computer), and an accessory genome consisting of the other, unnecessary, genes – imagine a joystick in our analogy). The size of the accessory genome can vary a lot for different bacteria. For instance, Escherichia coli is said to have an open accessory genome (i.e., much bigger than the core genome; imagine a Raspberry-PI and all the extension kits) whereas Chlamydia trachomatis has a closed accessory genome (i.e., small compared to the core genome; imagine an ATM).
For species having a small accessory genome, a single reference might be sufficient. However, if a core genome is large, we can use either pan-genomes, or multiple references. The word pan-genome tends to be used for different things: the computational biologists usually use it as a representation of all genomic content in a certain phylogenetic clade (which can be mathematically modeled using finite automata, hence as graphs) whereas some biologists define it as a concatenation of all (core and accessory) genes from the clade (without considering SNPs, etc.). Unfortunately, methods for graph references (especially for read mapping) are still not sufficiently developed, therefore, researchers still have to resort to reference sequence-based methods and use either concatenated genes or multiple reference sequences.
When particular bacterial species are studied using massive whole-genome sequencing of many isolates, researchers usually cluster the isolates to monophyletic sequence clusters and infer a special reference sequence for each of these clusters (see, e.g., this example for the pneumococcus). All of the inferred references should be then published together with the paper or in a separate data paper (see, e.g., the sequences from the previous example).
Disclaimer: I am not a biologist. If anything in this post is too simplified or incorrect, please, leave me a comment and I will try to fix it. I tried to make this text accessible for my younger self.

I would add that even for eukaryotes, the "reference" genome is sometimes not ideal. For instance, the lab where I was working was studying transposons in D. melanogaster, and their lab strain had a transposon that was absent from the "official" reference genome.
The reference choice also depends on the available annotations. In our case, we had to use an older genome because some of the annotations that were published in a paper could not be converted to the most recent version. This likely indicated that some of these annotations were artefactual, and arguably we might have ignored them.

If the programs you are using allow for it, take the most recent available genome, as it will be most likely to have the fewest errors.
Many of the well-known model organisms have an official release site (e.g. flybase for drosophila, wormbase for nematodes) from which genomes are fairly promptly propagated through to larger public databases, of which the most well-known are probably Ensembl and RefSeq Genome. There's usually enough cross-talk between the large data repositories that it doesn't matter all that much where the source is, as long as it is at least a few months since the last genome release.

The "best" assembly depends on the species. For more common species, Ensembl/UCSC/NCBI has their own version that corresponds to the most popular assembly. Usually the actual genetic sequence is identical across all of them, but the chromosome names and gene annotations will vary. If those three sources agree for a species, there is a "standard" reference.
For fruit fly specifically, FlyBase Consortium/Berkeley Drosophila Genome Project is the gold standard. Other species have their own organizations.


text:
QC measures for NGS sequencing
What are good means for performing quality control (QC) or NGS reads?
I'm aware of:

FastQC
NGS Screen
Kraken (e.g., for screening against contaminants)

What are other popular means for such QC?

This question is too broad. There are many available QC programs which all have their own advantages/disadvantages. Do you have a specific use case that you would like the QC to be applied to?

agree with @gringer. Right now the answer to this question is, roughly, whatever you get by googling 'NGS QC'. If you can provide a detailed use case - e.g. maybe you have 200GB of paired-end illumina data from 12 individuals of a non-model species, and you have a distant reference available - it will be a much more useful question and one that google couldn't solve.

Agreed. I'd be happy to rescind my downvote if the question text were updated.

MultiQC can merge all your different reports into a single one.
Which could be useful once you manage to know which QC tools to use.

FWIW, Multiqc is also really pretty...

We routinely run both FastQC and FastQ Screen on all of our raw sequencing reads. FastQ Screen is a tool for detecting cross-species contamination. MGA is another similar tool.
There are then lots of QC tools specific to different types of data, most of which run after alignment. For example RSeQC (RNA data), Qualimap and many many others. Without specifying what kind of data you have this is a bit difficult to make recommendations for though.
Phil

The quality control of ngs reads is heavily dependent on type of the project.
For genome assembly projects based on short reads, beside already covered checking quality of sequencing, you would like to look at the kmer spectra to find out, if your reads are going to make sense when they will be translated to De Brujin graph. You will get also a clues about genome coverage, genome size, repetitive content and small clue about heterozygocity. A lot of useful info about interpretation you can find in the README of GenomeScope.
A list tools I used:

jellyfish - to count k-mer frequencies
GenomeScope - a package for analysis and visualisation of k-mer frequencies (they recommend to use jellyfish for counting k-mer frequencies)
kmergenie - for prediction of optimal kmer for assembly

Using these tools can save you a lot of frustration if you accidentally sequence a contaminated sample or if your colleagues / a sequencing company have sent you a wrong file!

The kind of QC you do routinely depends on what your lab's focus is. We do a lot of low-quality, multiplexed DNA and RNA. If you routinely do fresh frozen whole genomes, your QC will be different.
Weighing in from the resequencing side of things (i.e. sequencing an organism that has a good reference), there are several things we are testing with quality control:

Did the library preparation succeed? Was the capture successful? Is our input material of high enough quality?
Did the sequencing work? For Illumina, did enough clusters pass filter? What's our Q30? Did we sequence into adapters?
Can we proceed with analysis? Are the indices on multiplexed samples correct?
Do we have the correct reference sequence? Do we have the right targets?

You may use a combination of tools:

bcl2fastq : if you multiplexed, check the undetermined indices file after converting to fastq to make sure you caught all of the indexed reads
FastQC: General sequencing quality, sequencing into adapters, handy flagging for obvious problems
Illumina interop/SAV: Cluster pass filter, Q30
Alignment (with BWA, bowtie, etc), bedtools, samtools, Picard: check for contamination, make sure we're using the right reference, overlap with target regions, depth of coverage, soft-clipping for low-quality sequence

If you are looking for NGS QC for your fastq, bam, bed and vcf files I would suggest a commercial tool called omnomicsQ.
It automates coverage analysis and sample pass, warn fail according according to your defined SOPs. Unlike fastqc it also comes with a database so that you can compare protocols, samples, and your performance to peer organisations. It also charts performance over time, exceptions,and  correlations between metrics


text:
Stable download URLs
One big problem that I'm regularly facing is that URLs for downloading Bioinformatics data (e.g., RefSeq releases or NCBI genome releases) disappear.
Does anyone have any good solution for this?

Do you have examples of specific URLs that have stopped working? Some projects try very hard to provide permalinks.

Genome releases are versioned, so the URL may move but the data is still there.

Could you specify which platform or which kind of data you are looking for a stable URL? Many sites provide a ftp server for big downloads

I think this is a good question, maybe with minor modifications: URIs for big databases rarely disappear but stable identifiers aren’t always obviously marked (e.g. Ensembl defaults to “current” release rather than a stable URI).

The question should be flavoured with some examples, but it rises a very good point.

I was referring to NCBI stuff, e.g., dbSNP where the latest release for GRCh37 has disappeared and has to be interpolated from an UCSC table.

The Persistent uniform resource locator or PURL is one such solution, these are designed to be a bit more robust than permalinks in so much as they are supposed to survive the change of domain name.  The bio ontology community already use them http://purl.bioontology.org/docs/index.html


text:
Structural variant calling for low-coverage PacBio data
PacBio is selling ~10x PacBio SEQUEL long reads as an upgrade to Illumina data for SV discovery.
In a clinical setting, the main requirements are proper sensitivity and specificity but also the processing of cohorts, at least families. This requires a genotyping step, such that it can be identified whether a given variant is shared by two or more individuals or whether it is not.
What are the tools of the trade for this task?
As having 50-60x PacBio reads is not an option from a economic point of view, one has to make do with 10x coverage.

There is an evaluation of PB Honey and Sniffles algorithms for low coverage PacBio datasets in this preprint and another evaluation is shown on this poster. Both reports agree that optimal is (surprisingly) combination of PB Honey and Sniffles.
Author of Sniffles have benchmarked Sniffles against PB Honey, where he shown that Sniffles performs significantly better. Take a look on this presentation (slide 15).
Another option is SMRT-SV, but I am not aware of any benchmarking.

SMRT-SV uses assembly-quality coverage of PB reads and sadly of low use. PB Honey and Sniffles are essentially useless for clinical germline studies as neither supports multi-sample calling and when comparing calls from both tools on patients where we also had Illumina X Ten genomes showed promise in PB data but severe problems with existing tools.

I can imagine that if the SV calls would take into account all the family data in one go, the results would be more reliable. However is it really a problem to do variant calls individually and then compare / merge them? I would guess that it would be still better than Illumina...

Has anyone successfully used either of these programs on MinION generated data?

@MatthewBashton for this you should open a new question...

I'm aware of the following (very few and suboptimal) options:

Sniffles -- sadly not very reliable in my experience, also no genotyping step or multi-sample support
PB Honey -- no genotyping step or multi-sample support