text:
Why do BAM files created by different tools have different file sizes?
I have a BAM created by Picard. I want to filter alignments by flags with samtools view. However, I noticed that even if I apply no filters, the output BAM is different from my input BAM. Are BAMs produced by different tools also different in size? How can I check if two BAMs are the same?

Hi medbe, thanks for asking your question, which fits well with this bioinformatics site. To improve the answers that you will get, it can be helpful to edit your question and build a bit of a story around it. What led you to want to filter alignments by flags? Why does it matter that the BAM files have a different format?

Could you elaborate on exactly how samtools/picard output is different? Is it just the file size or are you missing data from the optional fields of the alignment section?

It's worth bearing in mind that when outputting compressed BAM, as most tools do by default, they may well be using different levels of compression and/or different libraries, or versions of said, libraries for doing (de)compression which will result in different file sizes.  Additionally coordinate sorted BAM will compress more than unsorted BAM. The current version of Picard uses HTSJDK which in turn uses java.util.zip.Deflater/Inflater, current versions of samtools should be using the HTSlib library which in turn depends on the standard zlib library. You can see the effect of different implementations of zlib have on file size and execution time in benchmarking done by the samtools team.
However, in your case the best way to see if there is any difference between the BAM files is to rule out the effect of different levels of compression or libraries used for compression and save both BAM files as uncompressed.  Both samtools and Picard have options for disabling or changing the levels of compression, since the BAM compression standard BGZF is implemented on top of the gzip format it has inherited the ability, just like with gzip, to change the compression level from 0 to 9.
samtools view -bu will allow you to produce uncompressed BAM output (which is also handy for piping into other programs as it saves time wasted compressing decompressing what is essentially a stream).  Also note that samtools sort has a -l INT setting where INT can be set between 0 (compression off, as with -u) 1 (for fastest compression, but increased file size) or -9 (for maximal compression, with increased run time).  Some of the effects of increased runtime for higher compression settings might be ameliorated using the -@ argument which allows you to set the number of extra threads used for BAM compression, by default samtools won't use any.
Picard tools has a general setting COMPRESSION_LEVEL which is applicable to most of its tools setting this to 0, COMPRESSION_LEVEL=0 should disable compression.
So re-running whatever Picard tool you used in the first instance with COMPRESSION_LEVEL=0 will then enable you to check that the file is not further modified by samtools view -bu. The assumption here being that if both files have exactly the same content they should be the same size uncompressed, of course if they have trivial differences with regard to white space formatting things may still differ.

It's unlikely that two different mapping tools will give exactly the same alignment, scores, and match strings for the same sequence mapped to the same reference. For some sequence/reference alignments, it's impossible to determine which is the "best" alignment, and slight differences in code can have large effects on the chosen alignment.
However, even if the actual mapping location and match string are exactly the same (for example when using a tool like Picard to filter BAM/SAM files), different tools will incorporate different metadata with each mapping. This is allowed in the SAM file format specification by the addition of optional fields beyond the 11th column. There are a few standard optional tags that can be used in these fields, and additional custom non-standard tags can be used as well. It is very likely that Picard is adding additional metadata to the alignments in the BAM/SAM file.
There's an additional complication in that the underlying SAM alignment (and metadata) could be identical, but the BAM file can still have different file sizes. One reason for this is that BAM file compression methods can be changed. For example, alignment tools might choose a quick compression method, while filtering tools might choose a method that results in better compression.
Checking for alignment similarity is more difficult than just comparing files at a binary level, and your particular application (or context, or story) will change what the best method of comparison is. It would be helpful to know why you want to compare BAM files in order to provide a better answer to your question.


text:
How can I downsample a BAM file while keeping both reads in pairs?
I know how to downsample a BAM file to lower coverage. I know I can randomly select lines in SAM, but this procedure can't guarantee two reads in a pair are always sampled the same time. Is there a way to downsample BAM while keeping pairing information intact?

Could you please give a bit more explanation as to why it would be useful to downsample reads? There are a few different ways to do it, and the answers will depend on your goal. For example, you can do a digital normalisation to reduce the input read count for genome assembly, but that method would not be appropriate for a quantitative analysis of reads (e.g. for RNASeq).

samtools has a subsampling option:

-s FLOAT:
  Integer part is used to seed the random number generator [0]. Part after the decimal point sets the fraction of templates/pairs to subsample [no subsampling]

samtools view -bs 42.1 in.bam > subsampled.bam

will subsample 10 percent mapped reads with 42 as the seed for the random number generator.

To avoid possible confusion in the future, it should be noted that this works by hashing the read name. If for some terrible reason the mates don't have the exact same name (e.g., /1 and /2 weren't stripped) then this will break.

And also due to the way this works, the number of subsampled entries is not guaranteed to be the exact required percentage: some transformation of the hash residing in the [0,1] interval is compared against the sampling proportion. Depending on the read names present in the file, the number of sam entries effectively sampled will vary.

With this function, you can subsample a given BAM file to a user-specified number of reads using SamBamba. The function automatically calculates the scaling factor. Ideas in part taken from here
function SubSample {

## Calculate the sampling factor based on the intended number of reads:
FACTOR=$(samtools idxstats $1 | cut -f3 | awk -v COUNT=$2 'BEGIN {total=0} {total += $1} END {print COUNT/total}')

if [[ $FACTOR > 1 ]]
  then
  echo '[ERROR]: Requested number of reads exceeds total read count in' $1 '-- exiting' && exit 1
fi

sambamba view -s $FACTOR -f bam -l 5 $1

}

## Usage example, selecting 100.000 reads:
SubSample in.bam 100000 > subsampled.bam

Note that $2 is the number of reads. Say you choose 10000 reads but have paired-end data, you'll end up with 5000 fragments.

A totally different tactic...it helps if you aren't super picky about how many reads you have at the end.
Illumina read names contain the lane, tile, and xy coordinates of the read in them.  If you filtered by these, you would be guaranteed to get either both members of a pair, or neither.  I would recommend trying to avoid the edges of the flowcell.

It have used both sambamba and samtools. Either one works great to downsample to a desired read depth after calculating the downsample factor.
Sambamba can use multiple threads by using -t INT, where INT is the number of threads. It can also set the compression level for the BAM file using -l.
BAM="Sample.bam"
READS="10000000"
fraction=$(samtools idxstats $BAM | cut -f3 | awk -v ct=$myReads 'BEGIN {total=0} {total += $1} END {print ct/total}')

# With samtools
samtools view -b -s  ${fraction} ${BAM} > output.bam

# With sambamba
sambamba view -f bam -l 5 -s  ${fraction} ${BAM} > output.bam

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:
What is the difference between samtools, bamtools, picard, sambamba and biobambam?
After some google searches, I found multiple tools with overlapping functionality for viewing, merging, pileuping, etc. I have not got time to try these tools, so will just see if anyone already know the answer: what is the difference between them? Performance? Features? Or something else? Which one is generally preferred? Samtools?

This question is likely to inspire an opinionated debate about the merits (or otherwise) of different tools, which is discouraged on StackExchange. It would be helpful if a single question were asked, preferably one that tries to stay away from people's preferences.

In the "something else" you could add the quality of the documentation, which may be an important factor coming into play when deciding which tool to adopt.

The obvious answer is that different people wrote them. It's fairly common in bioinformatics for people with a computer science background to get frustrated with existing tools and create their own alternative tool (rather than improving an existing tool). Over time, tools with similar initial aims will have popular functionality implemented in them (and eventually have bugs fixed), such that it matters less which particular tool is used for common methods.
Here's my impression of the tools:

samtools -- originally written by Heng Li (who also wrote BWA). The people who now work on samtools also maintain the alignment file format specification for SAM, BAM, and CRAM, so any new file format features are likely to be implemented in samtools first.
bamtools -- this looks like it was written by Derek Barnett, Erik Garrison, Gabor Marth, Michael Stromberg to mirror the samtools toolkit, but using C++ instead of C
picard -- Java tools written by the Broad Institute for manipulating BAM/SAM files. Being written in Java makes it easier to port to other operating systems, so it may work better on Windows systems. I'm more familiar with picard being used at a filtering level (e.g. removing PCR duplicates), and for statistical analysis, but it links in with the Java HTS library from samtools, so probably shares a lot of the functionality.
sambamba -- a GPL2-licensed toolkit written in the D programming language (presumably by Artem Tarasov and Pjotr Prins). I haven't used it (and don't know people who have used it), but the github page suggests "For almost 5 years the main advantage over samtools was parallelized BAM reading. Finally in March 2017 samtools 1.4 was released, reaching parity on this."
biobambam -- written by German Tischler in C++. I also have no experience with this toolkit. This seems to have some multithreading capability, but is otherwise similar to other toolkits.

A comparison of sorting speed between SAMtools (version 1.2) and sambamba (version 0.6.3) here.

In light of older versions of Samtools being slower than Sambamba, sometimes you also have to consider the needs of the overall pipeline. For example, some older software requires old versions of Samtools to run, which may make it hard to gain the speed advantages from newer Samtools and lead to choosing to use a different program entirely instead of having to support different versions of the same tool in your pipeline.

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


text:
Why most aligners do not output the "X" CIGAR operation?
As I read the SAM spec, the "X" CIGAR operator represents a mismatch. This seems useful as we can know where are the mismatches without looking at the reference genome. However, many popular aligners such as BWA do not output "X". Why do they omit "X"?

Strangely enough, the mismatch information is stored in the MD tag (at least for BWA), so it appears to have been a conscious decision rather than a matter of simplicity.

They also generally don't output =.

@DevonRyan I think it will either output M or X/= so if they omit X they will also automatically omit =

The SAM format originally had only M, I, D, N, S, H, and P CIGAR operators.  See the original SAM specification (if you can view Apple Pages documents) and Table 1 in The Sequence Alignment/Map format and SAMtools (Li et al, 2009).  This was in line with previous tools using CIGAR strings, notably exonerate which introduced them with just the M, I, and D operators.
BWA-backtrack was written contemporaneously with the SAM format in 2008 and 2009 (and published in May 2009).  Its ChangeLog shows that it was outputting mismatch information in an MD tag from January 2009 and that MD was defined in the SAM specification of the time (and that the tag value's syntax was somewhat in flux in February 2009).  The MD tag is also described in that early v1.0 Pages-formatted SAM specification.
The = and X CIGAR operators were introduced later in SAM v1.3 as a result of this lengthy samtools-devel mailing list thread.  The characters used for the operators and the initial implementations were essentially in place by November 2009.
Since then, the =/X operators have not really taken over from M, as you've seen.  There's probably a number of factors contributing to this:

Their later introduction, long after the basic M/I/D operators and MD tag were well-established;
Conceivably their being specific to SAM and unavailable in other CIGAR flavours;
The MD tag still provides more information — X doesn't tell you what the mismatched reference bases were.


text:
Do you use SNAP for short-read mapping?
I am calling SNPs from WGS samples produced at my lab. I am currently using bwa-mem for mapping Illumina reads as it is recommended by GATK best practice. However, bwa is a bit slow. I heard from my colleague that SNAP is much faster than bwa. I tried it on a small set of reads and it is indeed faster. However, I am not sure how it works with downstream SNP callers, so here are my questions: have you used SNAP for short-read mapping? What is your experience? Does SNAP work well with SNP callers like GATK and freebayes? Thanks!

If you have the option of using a paid solution, I would strongly urge you to try Sentieon's implementation of the bwa algorithm. I've been using it for a few months and it is several times faster than vanilla bwa.

GATK best practices are explicably meant to consume BWA MEM generated BAM.  Whilst SNAP may be faster, the Broad will not have tested it for compatibility with GATK as such you can't guaranty using it won't have unexpected consequences.
As such you'd be better off using BWA MEM because I assume accurately called variation is always better than fast and incorrectly called variation.  The main issue you'll have is ensuring shorter split hits and mapping quality are reported in the same way as bwa MEM -M which GATK/Picard is expecting.  Ultimately however you'd be better off posting this question on the GATK forum.
It's also worth noting that the soon to be released GATK 4 will utilise bwaspark which can distribute it's alignment processes across Apache Spark for increase performance.  Consequently I can't see SNAP being adopted anytime soon.


text:
What's the best way to download data from the SRA? Is it really this slow?
I'm trying to download three WGS datasets from the SRA that are each between 60 and 100GB in size.  So far I've tried:

Fetching the .sra files directly from NCBI's ftp site
Fetching the .sra files directly using the aspera command line (ascp)
Using the SRA toolkit's fastqdump and samdump tools

It's excruciatingly slow.  I've had three fastqdump processes running in parallel now for approximately 18 hours.  They're running on a large AWS instance in the US east (Virginia) region, which I figure is about as close to NCBI as I can get.  In 18 hours they've downloaded a total of 33GB of data.  By my calculation that's ~500kb/s.  They do appear to still be running - the fastq files continue to grow and their timestamps continue to update.
At this rate it's going to take me days or weeks just to download the datasets.  Surely the SRA must be capable of moving data at higher rates that this?  I've also looked, and unfortunately the datasets I'm interested have not been mirrored out to ENA or the Japanese archive, so it looks like I'm stuck working with the SRA.
Is there a better way to fetch this data that wouldn't take multiple days?

The SRA is definitely capable of moving this data at higher rates. The throttling happens not at their end but somewhere else in the network connection. In my test just now, using fastq-dump, I get a throughput of ~11.5 MiB/s, more than 20 times the rate you observe.

I don't know if this can help answering your question, but providing the commands you used might be at least useful for other readers.

In general I find that its the conversion from SRA to fastq that is slow, not the downloading. You can try this by running NCBI's prefetch tool first to download the file (through using ascp if possible) and then fastq-dumping it. Also note that ascp works better when you limit the band width its allowed. We limit it to 100m/sec.

Proximity to NCBI may not necessarily give you the fastest transfer speed. AWS may be deliberately throttling the Internet connection to limit the likelihood that people will use it for undesirable things. There's a chance that a home network might be faster, but you're likely to get the fastest connection to NCBI by using an academic system that is linked to NCBI via a research network.
Another possibility is using Aspera for downloads. This is unlikely to help if bandwidth is being throttled, but it might help if there's a bit of congestion through the regular methods:
https://www.ncbi.nlm.nih.gov/public/
NCBI also has an online book about best practises for downloading data from their servers.
On a related note, just in case someone sees this and EBI/ENA is an option, there's a great guide for how to do file transfer using Aspera on the EBI web site:
https://www.ebi.ac.uk/ena/browse/read-download#downloading_files_aspera

Your command should look similar to this on Unix:

ascp -QT -l 300m -i <aspera connect installation directory>/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:<file or files to download> <download location>

In my case, I've just started downloading some files from a MinION sequencing run. The estimated completion time via standard FTP was 12 hours for about 32GB of data; ascp has reduced that estimated download time to about an hour. Here's the command I used for downloading:
ascp -QT -l 300m -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:/vol1/ERA932/ERA932268/oxfordnanopore_native/20160804_Mock.tar.gz .

In general, the best way to download SRA data is: don't download from SRA. However, as ENA has not be sync'd yet, I would recommend to download from SRA ftp and then convert to fastq locally. You can find files in the SRA format here. Downloading and then converting locally is much faster than direct retrieval from NCBI for some mysterious reasons.

Agree. I always download via ENA, it's reliably fast and you get the fastq files natively. I usually search with SRA as it's terrible on ENA.

By far the fastest method in my experience has been to use the SRAdb library in R. For most entries, you can download fastq files directly. Some older experiments don't have them, but I've still found it much faster to download SRA files via getSRAfile() and then to convert them using fastqdump than to use fastqdump directly.

+1 for providing an R solution. That said, I don't think NCBI keeps FASTQ. I speculate that SRAdb just downloads from ENA. Its source code seems to confirm my guess.

NCBI is moving slowly to cloud (AWS/Google) for hosting. As such the tradiational FTP urls do not always work. My solution involves using pysradb.
An entire project can be downloaded by:
pysradb download -p <SRP_ID>

A Google Colab notebook with metadata retrieval commands is here.

Thanks thats good to know.

Another possibility in 2020 is to use fasterq-dump which is promoted to be faster than previous fastq-dump.

I recommend using fastq-dl: downloads from either ENA or SRA, multithreaded, automatic retries.


text:
Merging sequencing data for ChIP-seq experiments
I need to merge sequencing data from different sequencing runs but for the same ChiP-seq library (HiSeq 2000).
Are there any potential advantages or disadvantages when merging files at .fastq or .BAM stage (alignment with Bowite/1.1.2)?

I don’t think it matters. Both are easy to merge (BAM via samtools merge, and (gzipped) FASTQ via cat), and neither method has specific disadvantages, unless your FASTQ files are sorted for some reason (but they generally shouldn’t be).
One advantage of keeping the FASTQ files separate is that it makes it slightly easier to parallelise the mapping step: just run the mapper in parallel on the separate FASTQ files. Although bowtie has an option (-p) for this, throughput from that is slightly worse than running the mapping on split files.

For ChIP-seq it shouldn't really matter. But do be aware that by default, samtools merge retains read group information (the @RG field in the header) from each input file. This could pose a problem for some downstream analyses (e.g. for the GATK HaplotypeCaller) if you want the merged data to be considered as all part of the same sample. You can change this behaviour using the -c option.

I do not think one need @RG information for a ChIP-Seq , I mean it is very unlikely that someone would like to do a variants calling with ChIP-Seq. So in any case it it would hardly matter. I would just not mention about the @RG here since people might get confused.

Fair point. But since olga did ask about advantages or disadvantages, I thought it would be worth mentioning. It might be useful in the future, or for other users, when merging files for other applications besides ChIP-seq, since samtools merge is widely used.

Agree with the others that it doesn't really matter. One thing to note though - if you're deduplicating your BAM files (you probably should for ChIP-seq data), make sure that you do this after merging.. :)

Note that if instead of one library, you have multiple technical replicates (ie multiple libraries from the same sample), you should do aligning for each technical replicate separate, then deduplicate, and THEN merge. This allows you to keep unique fragments from each library that might look like PCR duplicates (because they map to the same place) but are actually unique.


text:
Why is bwa-mem the standard algorithm when using bwa?
The industry standard for aligning short reads seems to be bwa-mem. However, in my tests I have seen that using bwa backtrack (bwa-aln + bwa-sampe + bwa-samse) performs better. It is slightly slower, but gives significantly better results in terms of both sensitivity and specificity. I have tested it using the genome in a bottle data and public samples (NA12878 and NA12877 among others) and found that backtrack consistently outperformed bwa-mem.
So why is bwa-mem the standard? Am I wrong and other tests have shown the opposite? I don't really see how since I tested using the most common datasets and validation data. Is it that the slight increase in efficiency outweighs the decrease in performance?
The only other explanation I can see is that bwa backtrack is designed specifically for Illumina reads and all my tests have been on Illumina data. Is it just that bwa-mem is "sequencer agnostic"? So that we can use the same algorithm irrespective of what sequencing platform is used? In that case, it makes sense to use backtrack for if we only deal with Illumina data and mem if we can have different sequencers. But, if so, seeing as Illumina is so widespread, why isn't backtrack used more often on Illumina data? I feel I must be missing something.

Heng Li would be the best person to respond to this....
Also it might be worth citing existing benchmarks or giving details of those which you have run.

There is always the bwa mem paper too https://arxiv.org/pdf/1303.3997.pdf

@MatthewBashton yes, I know. Unfortunately, I don't have them in a form that's easy to show and I ran them more than a year ago now so they might not even be relevant anymore. I am hoping that others have seen similar issues and know why one would be preferred over the other. Or that a review article has been published comparing them and my pubmed-fu just failed me.  If not, I guess I'll rerun them myself, look into it more deeply and post a new question.

What's your read length? http://bio-bwa.sourceforge.net recommends mem as best for >70 bp. It would be interesting to see a reproducible benchmark if you've found otherwise

@blmoore it's mostly ~150nt on average. Which is precisely why I was surprised to see this. However, I realize that unless I show my data (which I will have to regenerate, I am afraid) it is hard for you to just take my word for it. I was kind of hoping this was a well known issue and I was just ignorant.

Hi @terdon; our sequencing core team has looked at BWA-MEM vs backtrack for several data sets and have consistently found that MEM is both faster and more accurate than backtrack. We've also found that if you have significant quality drop-offs at the tail backtrack's performance suffers heavily... and that's another advantage of MEM: you don't need quality trimming, where backtrack needs reads to be mapped in full-length, which we have observed independently of, and in addition to Heng Li (see this thread).

@WyattMcDonnell thank you, that looks interesting although I have trouble reading that kind of chaotic thread. Do you think you could whip up an answer with these points?

our department uses bwa backtrack for illumina sequencing data as well. as a side note, minimap2 (written by heng li) outperforms bwa mem for long reads, so I am having trouble figuring out where bwa mem falls into the picture.

Closely related (arguably a duplicate): https://bioinformatics.stackexchange.com/q/15/3144. Still no empirical comparison of the effectiveness of the two approaches of the sort that the OP here was asking for, though.

Possible duplicate of Difference between BWA-backtrack and BWA-MEM

bwa mem is newer, faster, and [should be] more accurate, particularly for longer reads.
From the bwa man page (presumably in Heng Li's own words):

BWA is a software package for mapping low-divergent sequences against a large reference genome, such as the human genome. It consists of three algorithms:  BWA-backtrack, BWA-SW and BWA-MEM. The first algorithm  is designed for Illumina sequence reads up to 100bp, while the rest two for longer sequences ranged from 70bp to 1Mbp. BWA-MEM and BWA-SW share similar  features such as long-read support and split alignment, but BWA-MEM, which is the latest, is generally recommended for high-quality queries as it is faster and  more  accurate. BWA-MEM also has better performance than BWA-backtrack for 70-100bp Illumina reads.

Yes, I know that's what he claims. That's just not what I've seen. Do you have any benchmarks or a review article or anything demonstrating that bwa-mem is indeed more accurate?

Heng Li's paper on BWA is here. It compares the accuracy of BWA-MEM with BWA-SW, and demonstrates that BWA-MEM has a higher number of mapped reads for the same number of mismapped reads (see Fig. 1).

Yes, but he doesn't compare it to bwa backtrack which is what I am asking about.

Brian Bushnell seems to think that BWA-MEM is more accurate than BWA-backtrack.

Yes, many people do. But I haven't seen anyone give evidence supporting this and my tests have shown it isn't.

To be fair, you also haven't given any evidence supporting this. It would be great if you could add a graph, or other results into your question (preferably with a downloadable read/reference set).

Yes indeed. I should show my data. To be honest, I did this more than a year ago and just don't have them anymore. It is absolutely possible that I'm wrong and there is no reason whatsoever for anyone to believe me just because I say so. I was hoping someone could answer with benchmarks of their own or, even better, with a peer-reviewed article showing comparison data.


text:
What is the difference between NGS, 2GS, SBS and HTS?
I've come across a bit of confusion about the initialism NGS, so think it would be a good idea to clarify this term (and similar terms like 2GS, SBS, and HTS) for this site with a bit of discussion.
What are the most common initialisms used to describe different sequencing technologies?

I think your answer is good, but also feel these terms are mainly used for marketing; much more informative to name the sequencing company and platform. NGS is the worst, some 'NGS' like SOLID and 454 are basically retired

NGS isn't as bad as SMRT, which, despite the expansion of "single-molecule real-time sequencing", is considered specific to PacBio and frowned upon as a tag for nanopore sequencing (which is also single-molecule and real-time).

I am not sure how appropriate it is at this point to still refer to sequencing as next-generation sequencing. The leading NGS technology is Illumina/Solexa that has been around for over 10 years at this point. 454 was around even earlier. It's not really "next" at this point.
Opinions aside, I would also add "third-generation sequencing" to that list, referring to long read technologies like PacBio and Oxford Nanopore. See this question on Biology SE for more details.

I have a bit of trouble classifying PacBio in a sequencing generation. It is a model-based sequencing technology (dependent on a polymerase and pre-mixed nucleotides for sequencing), but produces single-molecule reads sequentially over the course of a run. I therefore consider it to be somewhere between second and third generation sequencing technologies.

Here are my attempts at definitions:
Sanger: A method of sequencing that depends on chain-terminatiing dideoxynucleotides. This sequencing uses the differential flow of DNA sequences of different lengths through a gel to determine the original DNA sequence, producing a single sequence per reaction container.
NGS: Next-generation sequencing, also referred to as 2GS (second-generation sequencing). This term is used to describe the first wave of sequencing technologies that followed Sanger sequencing technology. The use of NGS has become more confusing with the advent of long-read sequencing, because it's a common assumption that "next-generation" refers to the most recent technology (which is incorrect in this case).
HTS: High-throughput sequencing. This term describes any type of sequencing technology that produces large amounts of data, usually in the form of millions of different sequences produced from the same sequencing run.
SBS: Sequencing by synthesis. This term describes a method of sequencing that depends on the synthesis of [DNA] bases in order for sequencing to be carried out. This definition can extend into long-read sequencing (e.g. PacBio sequencers depend on synthesis during sequencing), but is more typically associated with only the second-generation sequencing technology.

For me ngs and hts are the same thing - any massively parallel sequencing (anything after sanger). That is further divided to 2nd generations sequencing (short reads) and 3rd generation sequencing (long reads). Sequencing by synthesis can by whatever sequencing that uses template and synthesis of free nucleotides (454, Illumina, PacBio).

If you have to divide "next-generation" into multiple generations, it's probably not a good definition.

Some people will distinguish the ABI SOLiD and Nanostring Hyb+Seq sequencing technologies as "sequencing by ligation" rather than "sequencing by synthesis". I currently lump these all together as SBS because I haven't found a better descriptive name for the collective term.

I did bit of digging about "Sequencing by synthesis" : it is defined on an Illumina flyer, page 15, it is defined though florescent died nucleotides... PacBio write on their web : "Also, while it is a form of ‘sequencing by synthesis’, it operates in real-time rather than incorporating a single base (or type of base) at a time."

It's not quite "real-time" sequencing, as the base-calling is done by decoding fluorescence "movies" for each photodetector (ZMWG)... but then again nothing is truly real-time unless it's detecting and calling single bases at the same time.

Also, look at answer on very similar question on biology.SE. What you write here is really not what the most of scientist see behind these terms...

The terms "second-generation sequencing", "third-generation sequencing", "short-read sequencing", and "long-read sequencing" are used commonly but (as far as I know) lack commonly used acronyms/initialisms.
In the forensic genetics community, the acronym MPS is used commonly to refer to "massively parallel sequencing." This definitely includes 2G/short-read sequencing (such as Illumina). It's less clear whether MPS includes 3G/long-read sequencing such as PacBio or Oxford Nanopore, but since these are rarely used in a forensics context it's a moot point.
Of course, it's important to acknowledge (as @burger has done) that there is a lot of imprecision and opinion regarding the appropriate scope of each term. Some folks seem to want a term to refer to everything post-Sanger, which includes an incredible diversity of platforms with respect to sequencing technology, read length, throughput, and error profiles. Some very reasonable questions should make us think seriously about our chosen terminology: How long do (for example) Illumina reads have to get before they're no longer classified as short reads? How much throughput is required of a sequencer to be classified as "high" throughput? When do these distinctions really matter?

authors: ลูกชิ้นเด้ง เตี๋ยวหนับน้ำตกตับ


text:
What is the difference between SAM mapping quality and Blast E-value?
Blast reports E-values, but short-read mappers report mapping qualities. Are they the same thing? Can they be converted to each other? If not, why blast doesn't report mapping quality while short-read mappers do not report E-values?

The E-value and the mapping qualities are two very different things.
The E-value is "a parameter that describes the number of hits one can 'expect' to see by chance when searching a database of a particular size". More details can be found here: https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=FAQ#expect
The mapping quality is an attempt to estimate the probability that a given base from a given read is mapped correctly to a particular place in the reference genome. Different aligners calculate mapping qualities very differently, so there is no simple way to compare them across aligners. But in general the mapping quality will consider things like: (i) the quality of the base call; (ii) the repeat structure of the reference; (iii) the alignment algorithm; (iv) whether the read has a mapped pair; (v) anything else the author of the mapping software thought might help.
E-values and mapping qualities are, therefore, measuring two fundamentally different properties with very different uses.


text:
Is spark widely used in bioinformatics?
I learned that GATK 4 is using Spark for parallelization. I googled around, though I am still not quite sure how spark really works and how to use it in practice. Besides GATK 4, are any other bioinformatics tool using spark? Generally, is spark widely used? Is it a necessary skill to learn? Thanks in advance.

The only bioinformatics tool other than GATK4 that I am aware of that uses spark is Hail (a Spark based replacement for Plink). Hail is also supported by researchers at the Broad.
Most places I have worked at have not switched over to Spark.  As such, I don't think it is widely used generally. Therefore, I don't think knowing Spark qualifies as a necessary skill to learn for bioinformatics at this time.
However, I think Spark is a superior paradigm for handling the large amount of data that bioinformaticians routinely deal with, and I think we will soon see the field move towards using it more.  Any bioinformatician would do well to acquaint themselves with Spark and to play around with Hail.

From reading into that a bit more, it looks like calling Hail a Spark-based replacement for Plink is selling it short. It's mostly true, but a bit like calling a car a wheel-based replacement for a horse. The efficiency of variant storage, and ability to distribute the workload, seems like quite a disruptive change in how variant analysis, filtering and QC is currently done.

I agree with you.  On the other hand, I think there is value in thinking of a car as a wheel-based replacement for a horse O,o

This is my impression as well: few tools are using Spark as yet, but there is a lot of interest in Spark as an idiom for concurrency. It's more suited to biological problems than a many other choices and is less complicated. So we'll probably see a lot more use of it in the future.

ADAM and avocado are Spark-based alignment and variant calling tools under active development by a collaboration (http://bdgenomics.org) that also includes the Broad, but I do not believe they have wide adoption.

I think others have provided fine answers to your question already. I just thought it might be relevant to mention about a tool that allows using Spark to distribute computations that re-use existing commandline bioinformatics software ... so, avoiding the need to re-implement proven algorithms in a language that Spark supports, EasyMapReduce.
This perceived need to re-implement algorithms, has otherwise been one of the biggest road blocks for wider adoption of Spark in the bioinformatics community in my understanding.
EasyMapReduce was presented SparkSummit East this year (video link here, as shown on the repository home page).


text:
How to extract exome on-target reads in batch?
I was given a list of target regions in BED and many exome alignments in BAM. I was asked to extract on-target alignments from these BAMs to save disk space. I know I can use bedtools to extract sub-BAMs. I am thinking to write a script to apply bedtools to all BAMs at my hand, but I speculate there may be some more convenient command lines to achieve my goal. How would you do? Thanks.

It sounds like bedtools intersect will work for you:
bedtools intersect -wa -a <alignment.bam> -b <region.bed> > <intersect_alignment.bam>

When a BAM file is used for the A file, the alignment is retained if overlaps exist, and exlcuded if an overlap cannot be found.  If multiple overlaps exist, they are not reported, as we are only testing for one or more overlaps.

It's easy enough to wrap this up into a for loop for batch processing (line breaks can be removed if desired):
 for alnBam in alignmentFiles*.bam;
   do echo ${alnBam};
   bedtools intersect -wa -a ${alnBam} -b <region.bed> > intersect_${alnBam};
 done

There are a lot of other options which should cover most related uses. For more details, just run "bedtools intersect" (with no arguments).

In your example shell code, ${x} should be ${alnBam}.


text:
Influential papers in bioinformatics
As a newcomer, what are some influential research papers published in the past few years that I should read to catch up in the field?
Edit: @Devon Ryan answered it perfectly. I tried a google search but it turned up a bunch of opinionated articles that varied. The biostar discussions are perfect and was the link I was searching for. I hope someone else new to the field finds this useful!

I think this is far too broad and opinion based for this site

I was pretty sure someone would point that out. It depends on how broad of a range we want to cover. Are we just here to answer specific technical questions? Or is this going to be more of a community? I personally agree and think this question might be better for the bioinformatics subreddit or Biostar.  I was just testing the waters while looking for an answer.

EMiller, The concern with this question is there's no problem statement... so the question is to generate an unbounded list of whatever folks care to post. An occastional, well-regarded "big list" question isn't generally a problem unless it becomes excessive, but please be careful of the inference that, if folks don't agree, this "isn't a community." I don't believe that was your intention here, but folks sometimes throw that out there as a backhanded way to discourage further discussion.

@EMiller there are some rules that apply to all Stack Exchange sites. The SE model is strictly Q&A so broad, discussion-based questions or requests whose answer will perforce be a long list (such as this one) will always be too broad.

@terdon "there are some rules that apply to all Stack Exchange sites. The SE model is strictly Q&A". I thought it was up to the moderators to decide when a question is unfit. That's a pity because on biostars I find questions like this often very interesting and still the volume of "bad" questions isn't too annoying, IMHO. Personally I don't see much point in having a site strictly for Q&A and another for general discussion. As I said, biostars is both and works well for me. Just my 2p...

@Robert I meant that a new commer with a specific question could just pop in and ask a question, without watching the community first like you might have to somewhere like biostar.

Over on biostars there's a thread like this every year or so. I'll link to the 2016 edition and the (much shorter) 2015 edition.
My personal picks from those would be:

ExAC
salmon, which is now published
kallisto, which is also now published


text:
Uppercase vs lowercase letters in reference genome
I am using a reference genome for mm10 mouse downloaded from NCBI, and would like to understand in greater detail the difference between lowercase and uppercase letters, which make up roughly equal parts of the genome. I understand that N is used for 'hard masking' (areas in the genome that could not be assembled) and lowercase letters for 'soft masking' in repeat regions.

What does this soft masking actually mean?
How confident can I be about the sequence in these regions?
What does a lowercase n represent?

What does this soft masking actually mean?

A lot of the sequence in genomes are repetitive. Human genome, for example, has (at least) two-third repetitive elements.[1].
These repetitive elements are soft-masked by converting the upper case letters to lower case. An important use-case of these soft-masked bases will be in homology searches: An atatatatatat will tend to appear both in human and mouse genomes but is likely non-homologous.

How confident can I be about the sequence in these regions?

As you can be about in non soft-masked based positions. Soft-masking is done after determining portions in the genome that are likely repetitive. There is no uncertainty whether a particular base is 'A' or 'G', just that it is part of a repeat and hence should be represented as an 'a'.

What does a lowercase n represent?

UCSC uses Tandom Repeat Finder and RepeatMasker for soft-masking potential repeats. NCBI most likely uses TANTAN. 'N's represents no sequence information is available for that base. It being replaced by 'n' is likely an artifact of the repeat-masking software where it soft-masks an 'N' by an 'n' to indicate that portion of the genome is likely a repeat too.
[1] http://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1002384

Informative answer, but I think it's controversial to say the human genome is "(at least) two-third repetitive elements"; the P-clouds method you cite is quite permissive and half is a more commonly accepted figure. And soft-masking doesn't involve masking all repeats generally, just interspersed repeats and low complexity sequences. Also there is always uncertainty around base calling and assembly building, and more so for repetitive sequences, although mm10 is one of the best assemblies of course...

The use of lower/upper case letters and N/n letters in genomes sequences is not completely standardised and you should always check the specification of the resource you are using.
Lower case letters are most commonly used to represent “soft-masked sequences”, a convention popularised by RepeatMasker, where interspersed repeats (which covers transposons, retrotransposons and processed pseudogenes) and low complexity sequences are marked with lower case letters. Note that larger repeats, such as sizable tandem repeats, segmental duplications, and whole gene duplications are not generally masked.
However, there are other uses for lower/upper case letters, for example, Ensembl have used upper/lower case letters to represent exonic and intronic sequences respectively.
N and n nucleotides may represent “hard masked sequences”, where interspersed repeats and low complexity sequences are replaced by Ns. But N/ns may alternatively represent ambiguous nucleotides, indeed this is the IUPAC specification.
Also note occasionally (although fortunately rarely) X/x is used to represent ambiguous nucleotides or “hard-masked sequences” too.

Lower case nucleotides commonly denotes a soft masked sequences. How exactly the genome was masked you can find in FAQ of NCBI:

Are repetitive sequences in eukaryotic genomes masked?
Repetitive sequences in eukaryotic genome assembly sequence files, as
  identified by WindowMasker, have been masked to lower-case.
The location and identity of repeats found by RepeatMasker are also
  provided in a separate file. These spans could be used to mask the
  genomic sequences if desired. Be aware, however, that many less
  studied organisms do not have good repeat libraries available for
  RepeatMasker to use.

IMHO, low complexity regions are always more likely to be missassembled than high complexity sequences. However, this will be problem for non-model organisms. I would guess that the reliability of the softmasked regions of mouse genome will very high.
No idea, looks like an artefact.

An example of usage of the soft mask
Mapping of sequence to reference usually starts with perfect matches of seeds (substrings) of the mapped reads and the reference sequence. Soft masked (low complexity) regions are not used for matches of seeds, but they are used only for the extension of the alignment if there was a seed in a neighbouring region. This application of softmasking applied to problem of long read assembly is described on this blog.


text:
The effects of incomplete bisulfite conversion upon mapping efficiency
This question has also been posted on Biostars
I have sequenced numerous multiplexed pools of BS amplicon-seq libraries derived from human samples on a MiSeq over the past few weeks. I have been utilising trim-galore and Bismark for alignment and am finding the mapping efficiency to be really low for two pools (55% & 30% respectively) both of which had a cytosine per base sequence of around 10-20% throughout the entire read in their fastqc files.
The high C content visible in the fastqc files makes me think the poor mapping efficiency is due to poor bisulfite conversion, as this would be expected to be close to zero is bisulfite conversion had actually took place.
I am new to the field so any help would be greatly appreciated.

It may be good to familiarize with the steps behind mapping bisulfite converted reads. In this post on biostars I give a quick breakdown of the steps involved if you want to have a look (but of course there is a lot of material available on it).

Bisulfite conversion efficiency has no effect on the mapping rate in bismark and similar tools. The reason is that the reads are fully bisulfite converted in silico before alignment to minimize mapping bias.
I would suggest that you play around with the settings handed to bowtie2, such as using local alignment and modifying the --score-min option to allow more mismatches. Alternatively, you might try a different aligner like bwameth, which will always do local alignment.

Thank you. --local flag in the bowtie2 docs looks promising, I'll give this a shot when I get into the office. Will Bismark parse any bowtie2 options/flags to bowtie2 (I ask as --local is not present in bisamark docs but only within bowtie2 docs)? What can I conclude about the high cytosine per base sequence in the fastqc files? I've generally found this is at 0% in previous analyses or 10% at the ends, at which point i would trim the ends.

I don't think that bismark will handle the --local flag correctly, since it expects reads to either fully align or be discarded, without any soft clipping.

If after running bismark the MBIAS plots don't look relatively flat, you can pass a filter parameter to bismark_methylation_extractor, so that e.g. you ignore the last 10bp of each read.

@719016: You might be correct, it's been a few years since I've looked at the bismark code, I generally try to avoid using it since it's so incredibly slow and used to give poor results.

As Devon Ryan also proposes, you may want to try bwameth, which is a simple wrapper around bwa-mem, and does keep the soft-clipped alignments if they exist. You can then apply another script for counting the CpGs from the bam, e.g. this one: https://github.com/dariober/bioinformatics-cafe/tree/master/bam2methylation

I generally recommend MethylDackel for methylation extraction and QC, but I'm biased :P


text:
What are the advantages and disadvantages between using KEGG or Reactome?
As enrichment analysis a usual step is to infer the pathways enriched in a list of genes. However I can't find a discussion about which database is better. Two of the most popular (in my particular environment) are Reactome and KEGG (Maybe because there are tools using them in Bioconductor).
KEGG requires a subscription for ftp access, and for my research I would need to download huge amounts of KGML files I am now leaning towards Reactome
Which is the one with more genes associated to pathways ?
Which is more completely annotated ?
Is there any paper comparing them ?

Can you rephrase this as "what are advantages and disadvantages" instead of "which one is better"? SE doesn't like the latter for the QA format, but I think the former is fine.

Done @MichaelSchubert you could change that yourself. But I think both are fine if the better is posted not in general (better for X and Y) YMMV

This BioStars post may be useful, although it's several years old https://www.biostars.org/p/3432/

Thanks @Chris but that question is 6.5 years old and things have changed in that time (including the fee for ftp access to KEGG).

@Llopis yes, I did edit in that detail just before you responded. I'm not saying this makes your question redundant, just pointing towards a possibly useful link

You should take a look at http://www.pathwaycommons.org/ for downloading pathway data.  Maybe you could go with both (or more), rather than just choosing one.

@woemier Yes, I actually I am using both Reactome and Kegg, but the problem is comparing if the pathways of one database and the other are equivalent regardless of the name

One big downside of KEGG is the licensing issue. One big advantage of Reacome are various crosslinks to other databases and data.
ad 1, This depends on which pathway, they are both primary databases. Sometimes other databases that for instance combine data of primary databases have better annotation of pathways (there is an example in the review paper bellow)
ad 3, There is very extensive relatively new (2015) review on this topic focused on human pathways: Comparison of human cell signaling pathway databases—evolution, drawbacks and challenges. However I could not find there which one is more complete ...

For you the main point would be whether an enrichment analysis is going to give you an informative answer. That's what will make a particular database better. And there's all sorts of subjective decisions made in their construction as to what what a pathway is, what to include, where to draw boundaries, etc. so different databases will give different answers and it will not be clear which is more correct.
So browse over the references / pointers people have given, select a service and use it and no others. Don't jump around databases until you get an answer you like, that's just fishing.

+1 for the last sentence, is so easy to say I will check with X once one have a "good" results and if both agree then go on.

One big advantage of Reactome, in my opinion, is its visualization using the web interface.
Many pathways (in Reactome and KEGG) consist of genes / proteins that are up- and down-regulated through the respective pathway. If you do a simple overrepresentation analysis this is not taken into consideration. Therefore, you might end up seeing a pathway as "overexpressed" although only the down-regulated genes were observed more frequently.
In Reactome, you can zoom in on the different pathways very conveniently and then pick up these inconsistencies. I haven't really found a public database and tool that can take these different regulations into consideration. Therefore, you'll probably always will need some manual investigation of your data. In my opinion, this is easier with Reactome than with KEGG.


text:
Can I model technical replicates in DESeq2?
I’d normally use collapseReplicates (or do the collapsing upstream) to handle technical replicates.
However, in my current RNA-seq experimental design, samples were sequenced twice using different library preparation protocols, leading to marked differences in the resulting count estimates (in fact, the choice of different protocol explains most of the variance in the bulk data according to PCA). I would therefore like to model these differences by adding them as a covariate to the DESeq2 model.
I thought I could just add another factor to the design formula, equivalent to the “pasilla” example in the DESeq2 vignette, which adds the factor type for the library type (single-end vs paired-end). However, the pasilla vignette states that these different libraries are actually independent biological replicates, not technical replicates as I had previously assumed.
I’m now concerned that having technical replicates in the design might artificially inflate the power. Can the design be salvaged?
It’s worth noting that we also have biological replicates, and a full rank design; i.e. every combination of treatment and library prep, with biological replicates.

If they're truly technical replicates, then there's no way to model them using DESeq2*, as you've alluded to with the collapseReplicates function. DESeq2/ Mike Love's general recommendation with collapseReplicates is to just add the reads together for technical replicates.
If you want to model them instead of collapsing them down, you can voom transform your data and follow an example similar to section 11.3 in the Limma users guide.
(Code example from 11.3):
Design
| FileName  | Cy3   | Cy5  |
| --------- |:-----:| ----:|
| File1     | wt1   | mu1  |
| File2     | wt1   | mu1  |
| File3     | wt2   | mu2  |
| File4     | wt2   | mu2  |

Example
> biolrep <- c(1, 1, 2, 2)
> corfit <- duplicateCorrelation(MA, ndups = 1, block = biolrep)
> fit <- lmFit(MA, block = biolrep, cor = corfit$consensus)
> fit <- eBayes(fit)
> topTable(fit, adjust = "BH")

*Edit: As @Devon Ryan mentions, you can design a paired model in DESeq2, but beware of making this full rank.

This answer makes me think that we need table layout support on this website: tables of features/counts/… will be extremely common here.

@KonradRudolph - Agreed, It was a shame when I found out that markdown tables weren't supported!

Practically speaking, there's no way to include the technical replicates in that design (in DESeq2 at least). Your concern regarding inflating the power is exactly correct and the only way to combat that would be to add a pairing factor like one might do with case-control or tumor-normal studies. That is, something like:
  group libraryPrep sample
1    WT           A      1
2    WT           B      1
3   MUT           A      2
4   MUT           B      2
5    WT           A      3
6    WT           B      3
7   MUT           A      4
8   MUT           B      4

Here sample pairs the technical replicates together. However this ends up either being rank deficient or in the end not more informative.


text:
samtools mpileup empty when filtering out flags
I produced a bam file by aligning reads to a small set of synthetic sequences using bwa-mem.
I am heavily filtering reads that are not paired and of a certain orientation.
Applying the filtering, I get a few thousands of reads:
samtools view -h $myfilebam | \
samtools view -h -F4 - | \
samtools view -h -F8 - | \
samtools view -h -F256 - | \
samtools view -h -F512 - | \
samtools view -h -F1024 - | \
samtools view -h -F2048 - | \
samtools view -h -f16 - | \
samtools view -h -f32 -  | wc -l

Gives me 89502 reads.
If I then pipe this into samtools mpileup, I get no results:
samtools view -h $myfilebam | \
samtools view -h -F4 - | \
samtools view -h -F8 - | \
samtools view -h -F256 - | \
samtools view -h -F512 - | \
samtools view -h -F1024 - | \
samtools view -h -F2048 - | \
samtools view -h -f16 - | \
samtools view -h -f32 -  | \
samtools mpileup --excl-flags 0 -Q0 -B -d 999999 - | wc -l

Returns 0.
I tried different combinations of filtering, and when I do both -f 16 and -f 32 returns empty, but if I do either of those, then it works:
samtools view -h $myfilebam | \
samtools view -h -F4 - | \
samtools view -h -F8 - | \
samtools view -h -F256 - | \
samtools view -h -F512 - | \
samtools view -h -F1024 - | \
samtools view -h -F2048 - | \
samtools view -h -f16 - | \
samtools mpileup --excl-flags 0 -Q0 -B -d 999999 - | wc -l

Returns 1056.
Any ideas why? My thinking was that it would work with --excl-flags 0.
EDIT: substituting mpileup for depth does work, and prints out each position and the depth as expected.
EDIT2: adding -q 0 to mpileup gives the same empty result.
Thanks in advance

Hi, could you add some details of how you got the synthetic reads and how you have aligned the to the reference?

Are you aware that -f 16 combined with -f 32 (btw, you can just -f 48 and skip the piping) are filtering for discordant reads? What happens if you add -q 0 to mpileup?

Since depth and mpileup share a fair bit of code, I have to wonder if you simply have no variants there. What happens if you specify -a and filter out the positions with no coverage?

BTW, I suspect this will be a case where you'll need to post the BAM file somewhere. There are a few (currently commented out) printf()s in the samtools code that could probably be used to track down what's actually happening here.

Devon already mentioned it in passing but it’s too important to get lost: Samtools flags are bit masks that you need to combine! Don’t call samtools repeatedly; call it once, with the combined bit flag. You can use the BROAD online tool to compute the combined flag if you don’t want to do it manually.

By using -h in the samtools view command, you're including all the header lines in your word count. If you happen to have about 89500 reference sequences, then the lengths of those would all appear in the header and inflate the -h word count, but not the mpileup count. Try piping it through an additional samtools view (i.e. without -h) and see if the counts change:
...
samtools view -h -f32 -  | \
samtools view | wc -l

Also, samtools mpileup by default only considers high-quality bases and concordant reads. Try adding a -A to your mpileup line (which stops anomalous read pairs from being discarded):
...
samtools mpileup -A -Q0 -B -d 999999 - | wc -l

Whether or not this is actually a good idea will be dependent on what you want to get out of the analysis, and what the downstream programs / analyses are expecting.

Bingo!  Adding -A, --count-orphans     do not discard anomalous read pairs did the trick!

i do have a follow up question. Is using the -A option a good idea? I have same problem as OP (i.e. empty mpileup), but can't find any useful info. of course with -A it seems to work but...how reliable will all downstream analysis be?

As with most questions of this nature, there's not a universal answer. If you want a useful answer, please post this as a separate question with enough context to explain what the eventual end-goal is, and what that downstream analysis is.


text:
Classifying samples based on marker gene expression
I have a few sets of marker genes that I can classify RNA-seq samples using semi-supervised clustering.  I would like to automate the process, however, I am struggling to find the ideal algorithm that could generate some kind of score for marker gene set from a given sample.
I presume that this is a standard analyses in many groups but I am not sure which method(s) are yielding good results in practice.

There was recently a similar question on Biostars that yielded no answers: https://www.biostars.org/p/239228/

I am surprised. It seems like such an important problem. Especially with scRNA-seq gaining popularity.

Since you mentioned scRNA-seq data, you might be interested in Buettner & al.: “Computational analysis of cell-to-cell heterogeneity in single-cell RNA-sequencing data reveals hidden subpopulations of cells”. It doesn’t quite address your problem but it shows some of the issues associated with identifying cell populations in scRNA-seq in particular, which are largely smoothed out in bulk RNA-seq.

I would consider using gene expression signatures to classify samples (especially cancer subtypes but the same principles apply to other problems of this type) one of the classic problems of bioinformatics. Quite a lot of work has been done on methods to derive gene sets that provide good classification performance. This is slightly different from your problem since you already have a gene signature but it may still prove useful.
These methods will typically fit a model that selects a (small) number of genes from genome-wide expression data that distinguish between the cell types/conditions in question, i.e. they derive a gene signature. The resulting model then allows for the classification of new samples. I've had success using GeneRave for this purpose (but note that this was designed for microarray data, I haven't used it with RNA-seq data and don't know how well it holds up there). A more recent paper relating to this issue can be found here.
So how does that help you? One option would be to fit one of these classifiers to gene expression data for the genes you already know to obtain a model that can then be applied to new samples automatically.

That's really helpful thank you so much.  I will give those a try or at least see how I can adapt their methods.

Tacking on to @Peter Humberg's caveat of GeneRave being designed for microarray data, you could voom transform your counts using limma to make them microarray-like.

When I need to compare cDNASeq expression with microarray, I use a transcript length normalisation applied to DESeq's VST transformation (which I call 'VSTPk'). More details of that can be found in the methods section of our Th2 RNASeq paper:
http://dx.doi.org/10.1084/jem.20160470


text:
How can I improve a long-read assembly with a repetitive genome?
I'm currently trying to assembly a genome from a rodent parasite, Nippostrongylus brasiliensis. This genome does have an existing reference genome, but it is highly fragmented. Here are some continuity statistics for the scaffolds of the current Nippo reference genome (assembled from Illumina reads):
Total sequences: 29375
Total length: 294.400206 Mb
Longest sequence: 394.171 kb
Shortest sequence: 500 b
Mean Length: 10.022 kb
Median Length: 2.682 kb
N50: 2024 sequences; L50: 33.527 kb
N90: 11638 sequences; L90: 4.263 kb

This genome is most likely difficult to assemble because of the highly repetitive nature of the genomic sequences. These repetitive sequences come in (at least) three classes:

Tandem repeats with a repeat unit length greater than the read length of Illumina sequencers (e.g. 171bp)
Tandem repeats with a cumulative length greater than the fragment length of Illumina sequencers, or the template length for linked reads (e.g. 20kb)
Complex (i.e. non-repetitive) sequence that appears at multiple places throughout the genome

Canu seems to deal quite well with the first two types of repeats, despite the abundance of repetitive structure in the genome. Here's the unitigging summary produced by Canu on one of the assemblies I've attempted. Notice that about 30% of the reads either span or contain a long repeat:
category            reads     %          read length        feature size or coverage  analysis
----------------  -------  -------  ----------------------  ------------------------  --------------------
middle-missing        694    0.07     7470.92 +- 5552.00        953.06 +- 1339.13    (bad trimming)
middle-hump           549    0.05     3770.05 +- 3346.10         74.23 +- 209.86     (bad trimming)
no-5-prime           3422    0.33     6711.32 +- 5411.26         70.92 +- 272.99     (bad trimming)
no-3-prime           3161    0.30     6701.35 +- 5739.86         87.41 +- 329.42     (bad trimming)

low-coverage        27158    2.59     3222.51 +- 1936.79          4.99 +- 1.79       (easy to assemble, potential for lower quality consensus)
unique             636875   60.76     6240.20 +- 3908.44         25.22 +- 8.49       (easy to assemble, perfect, yay)
repeat-cont         48398    4.62     4099.55 +- 3002.72        335.54 +- 451.43     (potential for consensus errors, no impact on assembly)
repeat-dove           135    0.01    16996.33 +- 6860.08        397.37 +- 319.52     (hard to assemble, likely won't assemble correctly or even at all)

span-repeat        137927   13.16     9329.94 +- 6906.27       2630.06 +- 3539.53    (read spans a large repeat, usually easy to assemble)
uniq-repeat-cont   155725   14.86     6529.83 +- 3463.16                             (should be uniquely placed, low potential for consensus errors, no impact on assembly)
uniq-repeat-dove    28248    2.70    12499.99 +- 8446.95                             (will end contigs, potential to misassemble)
uniq-anchor          5721    0.55     8379.86 +- 4575.71       3166.22 +- 3858.35    (repeat read, with unique section, probable bad read)

However, the third type of repeat is giving me a bit of grief. Using the above assembly, here are the continuity parameters from the assembled contigs:
Total sequences: 3505
Total length: 322.867456 Mb
Longest sequence: 1.762243 Mb
Shortest sequence: 2.606 kb
Mean Length: 92.116 kb
Median Length: 42.667 kb
N50: 417 sequences; L50: 194.126 kb
N90: 1996 sequences; L90: 35.634 kb

It's not a bad assembly, particularly given the complexity of the genome, but I feel like it could be improved by tackling the complex genomic repeats in some fashion. About 60Mb of the contigs in this assembly are linked with each other in a huge web (based on the GFA output from Canu):

The repetitive regions are typically over 500bp in length, average about 3kb, and I've seen at least one case which seems to be a 20kb sequence duplicated in multiple regions.
The Canu defaults seem to give the best assembly results for the few parameters that I've tried, with one exception: trimming. I've tried playing around a little bit with the trimming parameters, and curiously a trimming coverage of 5X (with overlap of 500bp) seems to give a more contiguous assembly than with a trimming coverage of 2X (with the same overlap).
If anyone is interested in having a look at these data themselves, called FASTQ files from Nippo sequencing runs can be found here. Raw nanopore signal files are available within ENA project PRJEB20824. There's also a Zenodo archive here that contains the GFA and assembly contigs.
I can use Illumina data to correct the Canu assembly, but that doesn't help with resolving the "type 3" repeats. The regions are sufficiently similar that illumina reads get mapped to multiple points in the genome. The Illumina contigs are high quality (i.e. they have good BUSCO scores, indicating few variant errors), but quite short. Any sniff of a repeat and the contig ends. I've got more than a few examples of regions that would make an Illumina read (even 10x linked reads) cower in fear.
Does anyone have any other suggestions on how I could resolve these complex repeats?
Computational solutions would be preferred, but resequencing is not out of the question.

"A few" 100kb reads won't help much. You need to apply the ultra-long protocol, which is different from the standard protocol.
You can't resolve 20kb near identical repeats/segdups with 10kb reads. All you can do is to bet your luck on a few excessively long reads spanning some units by chance. For divergent copies, it is worth looking at this paper. It uses Illumina reads to identify k-mers in unique regions and ignores non-unique k-mers at the overlapping stage. The paper said that this strategy is better than using standard overlappers, which I buy, but probably it can't resolve a 20kb segdup with a handful of mismatches, either.
Such mismatch-based approaches always have limitations and may not work for recent segdups/repeats. The ultimate solution is to get long reads, longer than your repeat/segdup units. The ~100kb reads in the recent preprint will be a game changer for you. If your ~20kb repeats are not tandem, 10X's ~100kb linked reads may help, too.

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


text:
What is the difference between ENCODE Tier 1, 2 and 3 cell types?
The ENCODE Experiment Matrix at UCSC lists the different available cell types under the categories "Tier 1", "Tier 2" and "Tier 3". What is the difference between these classifications?
What, for example, makes GM12878 a Tier 1 cell type and A549 a Tier 2 cell type?

You should read A User's Guide to the Encyclopedia of DNA Elements (ENCODE)

The tiers denote the original priority of sequencing and processing the samples. This is understandable given the number of marks and experiment types that ENCODE tried to sequence. Further details are available on ENCODE's website.
Edit: Just to expand a bit, tier one cells were supposed to get processed first. Then tier 2 came second. Of course tier 2 had a LOT more cells, so it got split into 2 and 2.5, since there were already apparently tier 3 cells.

From Encode's website:

... To aid in the integration and
  comparison of data produced using different technologies and
  platforms, the ENCODE Consortium has designated cell types that will
  be used by all investigators. These common cell types include both
  cell lines and primary cell types, and plans are being made to explore
  the use of primary tissues and embryonic stem (ES) cells.
Cell types were selected largely for practical reasons, including
  their wide availability, the ability to grow them easily, and their
  capacity to produce sufficient numbers of cells for use in all
  technologies being used by ENCODE investigators. Secondary
  considerations were the diversity in tissue source of the cells, germ
  layer lineage representation, the availability of existing data
  generated using the cell type, and coordination with other ongoing
  projects. Effort was also made to select at least some cell types that
  have a relatively normal karyotype...


text:
Does GISTIC (v 2.0) estimate amplified/deleted probabilities on a single sample basis?
Does GISTIC 2.0 estimate the background model:
G = -log(Probability | Background)
by permuting within the sample or across all samples in the set?
The paper describes the probabilistic scoring method based on permutations, but I could not understand if this permutation is performed within the sample only. The documentation page seems to suggest across samples, but this would mean that different sized sets might lead to different outcomes on the sample sample.
Basically, does it matter if the set is constituted by, for example, 10 samples of the same tissue of origin or ~1000 from multiple tissues?
Thanks,

I expanded a bit my question.

There doesn't seem to be much differences between GISTIC 1.0 and 2.0 as it says:

As with GISTIC 1.0, we obtain P-values for each marker by comparing the score at each locus to a background score distribution generated by random permutation of the marker locations in each sample

But on the supplementary material of the GISTIC 1.0 there is a more detailed explanation of the method. See the "Stage 2" section:

Second, we compare these G-scores to the distribution of scores expected if only random aberrations were observed. This distribution can be determined by rescoring the genome after permuting marker locations within each sample; we instead derive a semiexact estimate.

Furthermore in another section ("Stage 2: Aggregation of Data from Different Tumors to Differentiate Between Driver and Passenger Aberrations") it says:

To determine which of the aberrations identified in Stage 1 are likely to represent driver events, we aggregate the data from all tumors used in the analysis to generate summary scores for amplifications, deletions, and LOH. The statistical significance of each score is determined by comparison to the distribution of scores obtained by all permutations of the data (using a semiexact approximation), with correction for multiple hypothesis testing.

However, the relevant section (of the supplementary materials) seems to be "Null Hypothesis Generation: An Analytic Derivation of the Null Distribution", where it describes the semiexact approximation used.
In a supplementary file it describes it as "Generate all permutations of SNP labels within each sample to simulate datasets with random aberrations"
Conclusion:
It seems that the background probabilities are calculated within sample.
I can't say if it matters how the sets are constituted but I would say that the broader the sets are, the better estimation of the structural variations it will perform but the same estimation of the background will be done.
Ultimately you can check the code of the program, or test with 10 samples and replace one of them to see if the results are changing accordingly.

Thanks, yes it seems to be within sample. I'll do a test of the same sample with different data-sets and report here later.


text:
variant calling on ChIP-seq style data: samtools mpileup with minimal filters
I am running samtools mpileup (v1.4) on a bam file with very choppy coverage (ChIP-seq style data). I want to get a first-pass list of positions with SNVs and their frequency as reported by the read counts, but no matter what I do, I keep getting all SNVs filtered out as not passing QC.
What's the magic parameter set for an initial list of SNVs and frequencies?
EDIT: this is a question I posted on "the other" website, but didn't get a reply there.

Would you be able to try another caller just to check? Something like varscan?

@nuin I may give it a try, does varscan allow for switching off the filtering?

Could you add some examples of commands you tried?

@719016 Try this command, and change options like _--min-coverage _, --min-reads2 and --min-avg-qual. This will also help you check if your problem is samtools related or BAM file related

am just a bit worried that if the depth is not high and also ChIP-Seq data has biases then it is not high confidence SNPs. The whole point of HC SNP are the read depth as well which is not so much in traditional ChIP-Seq

Yes, maybe the calling should be restricted to near peak summits, then compare sample to sample on the subset of such SNVs

I used this in the past for ChIP-seq data and it generated SNVs:
samtools mpileup \
--uncompressed --max-depth 10000 --min-MQ 20 --ignore-RG --skip-indels \
--fasta-ref ref.fa file.bam \
| bcftools call --consensus-caller \
> out.vcf

This was samtools 1.3 in case that makes a difference.

I got it to work with these parameters. Thanks!

Another approach is htsbox. You can get a candidate list with:
htsbox pileup -Cvcf ref.fa -q20 -Q20 -s5 file.bam > out.vcf

Here, -q sets min mapping quality, -Q sets min base quality, -v outputs variants only -c outputs VCF, -C gives you base counts on both strands and finally -s5 requires at least 5 high-quality bases to call out an allele. It is useful when your data are failing the assumptions made by typical variant callers.
Why not samtools+bcftools or varscan? Transparency and speed. This command line simply counts based on the parameters you use. It does not apply additional operations. And because of this, it is over an order of magnitude faster than samtools mpileup or varscan. It is worth noting that samtools uses BAQ by default, which reduces FPs occasionally. However, BAQ is not quite necessary for longer Illumina reads and it hurts sensitivity at the same time.

I am definitely going to try htsbox, as you say it works by simply counting based on the parameters.


text:
How to select the most representative pathways from a gene enrichment analysis?
I have perform an enrichment analysis to a cluster of genes. The output is a list of pathways and their p-value (the pathways are selected because p-value < 0.05). The list is still quite long, so I want to reduce it. For that purpose I have a calculated the Dice coefficient of the pathways in a matrix $p$x$p$ where $p$ is the number of pathways in the list. I want both the ones that are more different (they overlap less, their Dice coefficient is lower) and the pathways more representative of the most similar pathways (So if a there is a group of 5 pathways that overlap over 0.8 take just one).
How can I select the most representatives pathways?
There is a similar tool for GO but it relays on discarding not significant GO, while here all the initial pathways are already significant.
If I do a clustering of the pathways using the Dice coefficient matrix I don't know where (or how) to cut.

I tried using the height to select the pathways. But I am unsure of the interpretation of height.
Some other tools I have seen use a multidimensional scaling plot, but I am not sure if performing it and cutting at certain point of the first dimension would help.

Are you using p-value < 0.05 or p-value < 0.05 / (# of terms tested)? You'll want to use the former to avoid false positives due to multiple testing

I'm guessing you'll have to settle with a "reasonably good" heuristic solution, since there are a few different variables in the problem you're trying to solve (number of pathways in group, desired overlap/distance between pathways, etc.) I can imagine many possible answers.

@CloudyGloudy correcting for multiple testing is already done, sorry I didn't mention it previously. Yes, I was playing with the idea of selecting an overlap/distance of 0.5 and keep those above that, but that would leave outside the pathways that are completely different from the others, but maybe I could select those above 0.75 and those below 0.25

One way to group similar significant pathways, is to quantify how many genes overlap between the pathways, and then use this in clustering (heatmap). I have made a tool in R which calculates the overlap index between GO terms and subsequently clusters them in a heatmap. Overlap index is the fraction of genes that overlap (number between 0-1). Also Pearson correlation can be used for clustering instead of the overlap index. My package (gogadget) works only with  goseq analysis, but you can use goseq also for reactome or kegg data.
I have used this overlap clustering approach for different data sets now, and usually we can reduce 200-300 GO terms into 10-20 functional groups.

Which overlap index does it use? There are several overlap indices. I would advice against using this approach in GO, there are specific methods to measure how similar two GO terms are. See GOSemSim.

The overlap index is defined by the number of overlapping genes divided by the number of genes in the smaller of the two gene sets. It is described in Bioconductor Case studies, chapter 13.3. This book is written by big names in bioinformatics (such as R. Gentleman and W. Huber).

Thanks for the new overlap index, I didn't knew it. BTW the way these similarity between genes is calculated is using other overlap indexes so this won't help.

Okay, but my point is not about which index to use but the clustering approach (same as gringer is suggesting). If you look in the userguide of my package on page 26-28, you'll see a heatmap. Usually I get good results with Ward.D and Euclidean. The tree can be cut in R, but you first might want to see and evaluate which gene sets are clustered together in which branches before you can find the right cutoff.

Sorry, my last comment was thinking in another question. Yes, this approach could work

This sounds like something that might be amenable to a clustered heatmap plot, or a correlation matrix plot, or something similar. Have you looked at a correlation matrix of the dice coefficient matrix (or maybe just a heatmap plot of that matrix without the correlation matrix)?
The corrplot package looks like it might be useful, in particular the hclust / drawing rectangles presentation.
I can't vouch for this package though; it's just something I found by a search for "R plot correlation matrix".

With a heatmap or a correlation plot I could observe the similarities between the pathways, as I currently do with the dendrogram and the MDS plot. However the question is how to select those more representative pathways. Thanks for your corrections and comments btw.

If you're happy with a more confident ranking of the most representative gene sets, rather than necessarily cutting down the list, you might try EGSEA. It uses an ensemble approach to give a ranking of the most relevant gene sets, and also produces an interactive HTML output with statistics, heatmaps, pathway maps, summary plots and GO graphs which allows you to examine the output at varying levels of granularity.
You can read the paper on bioRxiv or download the package from Bioconductor.

The ranking is already done, through the p-value (one could argue that a more elaborated test like the proposed by EGSEA would be better) but this doesn't seem to answer how to select relevant pathways from the output of EGSEA or other software/tools/methods

You shouldn't be ranking by p-value. The p-value is only an indication of whether the observed change is statistically significant, not an indication of the magnitude of the observed change.


text:
Is there an easy way to create a summary of a VCF file (v4.1) with structural variations?
I got a bunch of vcf files (v4.1) with structural variations of bunch of non-model organisms (i.e. there are no known variants). I found there are quite a some tools to manipulate vcf files like VCFtools, R package vcfR or python library PyVCF. However none of them seems to provide a quick summary, something like (preferably categorised by size as well):
type    count
DEL     x
INS     y
INV     z
....

Is there any tool or a function I overlooked that produces summaries of this style?
I know that vcf file is just a plain text file and if I will dissect REF and ALT columns I should be able to write a script that will do the job, but I hoped that I could avoid to write my own parser.
--- edit ---
So far it seems that only tool that aims to do summaries (@gringer answer) is not working on vcf v4.1. Other tools would provide just partial solution by filtering certain variant type. Therefore I accept my own parser perl/R solutions, till there will be a working tool for stats of vcf with structural variants.

What sort of summary are you looking for? Just raw counts (of number of insertions/deletions)? I'm just curious as I think personally I'd be interested to know how that variation is actually spread along a genome/contig/sequence spatially, would that be of interest to you?

I am comparing couple of relative species. However, I have not identified homologous regions yet (the annotation is so far under construction), therefore positioning of SVs is right now not very relevant. However, what can I do already is take a look if there is at least a difference in counts / sizes of SVs between species.

Maybe even better than summary would be to get type of SV and its size for every SV call. Then I could a summary in R would be trivial.

You've probably seen this already, but vcftools has a sibling called bcftools, which has a query function, that allows users to query a VCF/BCF to pull out fields and information and output their own format. It might not do exactly what you want, but might get you very close (enough to maybe just need a little post-processing in R?).

Great, that sounds promising. I will take a look.

The "my own parser" solutions. The information I was searching for in part of column INFO, namely in variables SVLEN and SVTYPE.
very fast SV types + counts (by @user172818 in comment) :
zcat var.vcf.gz | perl -ne 'print "$1\n" if /[;\t]SVTYPE=([^;\t]+)/' | sort | uniq -c

quite slow SV types + counts + sizes :
SV_colnames <- c('CHROM', 'POS', 'ID', 'REF', 'ALT', 'QUAL', 'FILTER', 'INFO', 'FORMAT', 'SAMPLE1')

ssplit <- function(s, split = '='){
    unlist(strsplit(s, split = split))
}

# note, capital letters just respect original naming conventions of the VCF file
getSVTYPE <- function(info){
    ssplit(grep("SVTYPE", info, value = T))[2]
}

getSVLEN <- function(info){
    SVLEN <- ssplit(grep("SVLEN", info, value = T))
    ifelse(length(SVLEN) == 0, NA, as.numeric(SVLEN[2]))
}

load_sv <- function(file){
    vcf_sv_table <- read.table(file, stringsAsFactors = F)
    colnames(vcf_sv_table) <- SV_colnames
    # possible filtering
    # vcf_sv_table <- vcf_sv_table[vcf_sv_table$FILTER == 'PASS',]
    vcf_sv_table_info <- strsplit(vcf_sv_table$INFO, ';')
    vcf_sv_table$SVTYPE <- unlist(lapply(vcf_sv_table_info, getSVTYPE))
    vcf_sv_table$SVLEN <- unlist(lapply(vcf_sv_table_info, getSVLEN))
    return(vcf_sv_table)
}

sv_df <- load_sv('my_sv_calls.vcf')
table(sv_df$SVTYPE)

According to the bcftools man page, it is able to produce statistics using the command bcftools stats. Running this myself, the statistics look like what you're asking for:
# This file was produced by bcftools stats (1.2-187-g1a55e45+htslib-1.2.1-256-ga356746) and can be plotted using plot-vcfstats.
# The command line was: bcftools stats  OVLNormalised_STARout_KCCG_called.vcf.gz
#
# Definition of sets:
# ID    [2]id   [3]tab-separated file names
ID  0   OVLNormalised_STARout_KCCG_called.vcf.gz
# SN, Summary numbers:
# SN    [2]id   [3]key  [4]value
SN  0   number of samples:  108
SN  0   number of records:  333
SN  0   number of no-ALTs:  0
SN  0   number of SNPs: 313
SN  0   number of MNPs: 0
SN  0   number of indels:   20
SN  0   number of others:   0
SN  0   number of multiallelic sites:   0
SN  0   number of multiallelic SNP sites:   0
# TSTV, transitions/transversions:
# TSTV  [2]id   [3]ts   [4]tv   [5]ts/tv    [6]ts (1st ALT) [7]tv (1st ALT) [8]ts/tv (1st ALT)
TSTV    0   302 11  27.45   302 11  27.45
# SiS, Singleton stats:
...
# IDD, InDel distribution:
# IDD   [2]id   [3]length (deletions negative)  [4]count
IDD 0   -9  1
IDD 0   -2  4
IDD 0   -1  6
IDD 0   1   4
IDD 0   2   1
IDD 0   3   3
IDD 0   4   1
# ST, Substitution types:
# ST    [2]id   [3]type [4]count
ST  0   A>C 2
ST  0   A>G 78
ST  0   A>T 2
ST  0   C>A 5
ST  0   C>G 0
ST  0   C>T 66
ST  0   G>A 67
ST  0   G>C 0
ST  0   G>T 1
ST  0   T>A 1
ST  0   T>C 91
ST  0   T>G 0
# DP, Depth distribution
# DP    [2]id   [3]bin  [4]number of genotypes  [5]fraction of genotypes (%)    [6]number of sites  [7]fraction of sites (%)
DP  0   >500    0   0.000000    333 100.000000

This is annotating what is explicitly in the VCF file, and that's SNVs (and INDELs). If you want a structural variant analysis (i.e. on a larger scale than single nucleotides), then you'll need to use something that does more than a summary of the VCF file.
While Inversions, large-scale deletions, and breakpoints are part of the 4.1 specifications, they are unfortunately not currently supported by bcftools.

You might want to try Bio-VCF. From the authors description

Bio-vcf is a new generation VCF parser, filter and converter. Bio-vcf
  is not only very fast for genome-wide (WGS) data, it also comes with a
  really nice filtering, evaluation and rewrite language and it can
  output any type of textual data, including VCF header and contents in
  RDF and JSON.

In addition it is a very fast parser. The rewrite DSL might help you customise your filtering and needs.

Has Bio-vcf got any functionality for structural variation as the OP requested? A quick scan of the GitHub page reveals nothing.  Whilst many of these tools are useful for filtering SNVs and short indels most fail to do anything of worth for SVs.

For structural variants  I think you can use SURVIVOR like SURVIVOR stats which was design specifically for such purpose (statistics on SV file).

Let's say you have an INFO field called SVTYPE to indicate the structure variant type.
Here is how you get the stats easily:

Install vcfstats: pip install vcfstats
Define a macro to extract the SVTYPE information, say in svtype.py:

from vcfstats.macros import cat
@cat
def SVTYPE(variant):
    return variant.INFO['SVTYPE']

Generate the statistics:

vcfstats --vcf <your vcf file> \
   --outdir <the output directory> \
   --macro svtype.py \
   --formula 'COUNT(1) ~ SVTYPE[DEL,INS,INV...]' \ # desired svtypes
   --title 'Number of SV types'

Then in the output directory you will have the stats file named Number_of_SV_Types.txt and a plot for it as well.
Check the details at https://github.com/pwwang/vcfstats

The advantage of this solution is that it using a standard parser. The disadvantage is that it requires several steps (unlike zcat var.vcf.gz | perl -ne 'print "$1\n" if /[;\t]SVTYPE=([^;\t]+)/' | sort | uniq -c  that will make the summary in one go).


text:
How can I output an identity matrix in progressiveMauve?
I'm just getting started with the Mauve aligner and I'm finding the documentation a bit lacking. I'm using the progressiveMauve tool from the command line and would like to output an identity matrix file. The set of output files it generates by default for progressiveMauve --output=outfile are outfile, outfile.backbone and output.bbcols.
The documentation for the tool lists an --input-id-matrix option but this doesn't seem to have any effect. The documentation page for output files does describe an identity matrix file but doesn't say anything about how to generate it.

mauve documentation have a similar command but the usage is --id-matrix=<file>. So very similar to the progressivemauve command (just doesn't include --input). Try it and see if it works.
I think it's a bug for sure as several other people on biostar have complained.

It indeed seems that the --input-id-matrix option doesn't work as expected, as Amar has pointed out.
I tried the following command:
./progressiveMauve genome1.fasta genome2.fasta genome3.fasta --output=test > test.log

And saw the following bits in test.log (truncated):
[...]
Genome conservation distance matrix:
0       0.0767545       0.187454
0.0767545       0       0.193385
0.187454        0.193385        0

[...]
genome content distance matrix:
0       0.0767545       0.187454
0.0767545       0       0.193385
0.187454        0.193385        0

bp distance matrix:
0.9     0.0125  0.075
0.0125  0.9     0.0375
0.075   0.0375  0.9
[...]

Could any of those three tables contain the information you are looking for?
p.s. as an aside, if you are looking into a way to obtain a whole genome similarity matrix you might want to look into mash or sourmash.


text:
How to read structural variant VCF?
The IGSR has a sample for encoding structural variants in the VCF 4.0 format.
An example from the site (the first record):
#CHROM  POS   ID  REF ALT   QUAL  FILTER  INFO  FORMAT  NA00001
1 2827693   . CCGTGGATGCGGGGACCCGCATCCCCTCTCCCTTCACAGCTGAGTGACCCACATCCCCTCTCCCCTCGCA  C . PASS  SVTYPE=DEL;END=2827680;BKPTID=Pindel_LCS_D1099159;HOMLEN=1;HOMSEQ=C;SVLEN=-66 GT:GQ 1/1:13.9

How to read it? From what I can see:

This is a deletion (SVTYPE=DEL)
The end position of the variant comes before the starting position (reverse strand?)
The reference starts from 2827693 to 2827680 (13 bases on the reverse strand)
The difference between reference and alternative is 66 bases (SVLEN=-66)

This doesn't sound right to me. For instance, I don't see where exactly the deletion starts. The SVLEN field says 66 bases deleted, but where? 2827693 to 2827680 only has 13 bases between.
Q: How to read the deletion correctly from this structural VCF record? Where is the missing 66-13=53 bases?

I'm pretty sure the END value is simply wrong in that example. Maybe send the info@ email address a note regarding this so they can fix it.

I've taken the liberty of asking the 1000 genomes project about this, I'll post the reply should it still be relevant when I receive it.

This line is wrong in several aspects. VCF4.0 is something like 5 years old and has long been deprecated. Its spec has also been moved to github. You should read the latest spec instead. The SV example there is correct.

This was fixed in the maintained VCF specifications in 2015 upon being reported as https://github.com/samtools/hts-specs/issues/84.

I just received a reply from 1000Genomes regarding this. I'll post it in its entirety below:

Looking at the example you mention, I find it difficult to come up with an
  interpretation of the information whereby the stated end seems to be correct,
  so believe that this may indeed be an error.
Since the v4.0 was created, however, new versions of VCF have been introduced,
  improving and correcting the specification. The current version is v4.3
  (http://samtools.github.io/hts-specs/). I believe the first record shown on
  page 11 provides an accurate example of this type of deletion.
I will update the web page to include this information.

So we can take this as official confirmation that we were all correct in suspecting the example was just wrong.

that was an amazing follow up! :)

So, first off, as others have pointed out, I'm pretty sure that example is just wrong. At least, the numbers don't match as you've pointed out.
That said, it is impossible to be sure without showing us the header of the VCF file as well. The INFO field (the 5th field of a VCF file) is very, very variable and depends entirely on the header lines. Each program (or human) implementing a VCF is free to choose to have whatever they feel like in the INFO field. However, each IDENTIFIER= needs to have an associated INFO line at the beginning of the file.
So, the SVTYPE, SVLEN, HOMLEN etc will have commented (start with a #) lines at the beginning of the file explaining what these values are. So check those, even though they're relatively standard, you never know, the obvious reading you used might be wrong despite its seeming so reasonable.
Here's a newer example of a VCF line for an SV taken from the current VCF specification:
##fileformat=VCFv4.1
##fileDate=20100501
##reference=1000GenomesPilot-NCBI36
##assembly=ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/sv/breakpoint_assemblies.fasta
##INFO=<ID=BKPTID,Number=.,Type=String,Description="ID of the assembled alternate allele in the assembly file">
##INFO=<ID=CIEND,Number=2,Type=Integer,Description="Confidence interval around END for imprecise variants">
##INFO=<ID=CIPOS,Number=2,Type=Integer,Description="Confidence interval around POS for imprecise variants">
##INFO=<ID=END,Number=1,Type=Integer,Description="End position of the variant described in this record">
##INFO=<ID=HOMLEN,Number=.,Type=Integer,Description="Length of base pair identical micro-homology at event breakpoints">
##INFO=<ID=HOMSEQ,Number=.,Type=String,Description="Sequence of base pair identical micro-homology at event breakpoints">
##INFO=<ID=SVLEN,Number=.,Type=Integer,Description="Difference in length between REF and ALT alleles">
##INFO=<ID=SVTYPE,Number=1,Type=String,Description="Type of structural variant">
##ALT=<ID=DEL,Description="Deletion">
##ALT=<ID=DEL:ME:ALU,Description="Deletion of ALU element">
##ALT=<ID=DEL:ME:L1,Description="Deletion of L1 element">
##ALT=<ID=DUP,Description="Duplication">
##ALT=<ID=DUP:TANDEM,Description="Tandem Duplication">
##ALT=<ID=INS,Description="Insertion of novel sequence">
##ALT=<ID=INS:ME:ALU,Description="Insertion of ALU element">
##ALT=<ID=INS:ME:L1,Description="Insertion of L1 element">
##ALT=<ID=INV,Description="Inversion">
##ALT=<ID=CNV,Description="Copy number variable region">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GQ,Number=1,Type=Float,Description="Genotype quality">
##FORMAT=<ID=CN,Number=1,Type=Integer,Description="Copy number genotype for imprecise events">
##FORMAT=<ID=CNQ,Number=1,Type=Float,Description="Copy number genotype quality for imprecise events">
#CHROM POS     ID        REF             ALT QUAL FILTER INFO                                               FORMAT NA00001
1      2827694 rs2376870 CGTGGATGCGGGGAC C   .    PASS   SVTYPE=DEL;END=2827708;HOMLEN=1;HOMSEQ=G;SVLEN=-14 GT:GQ  1/1:13.9

Note how the numbers do match and also note how each of the subfields in the INFO field is explained with an ##INFO line.