# Question
Title: How can I systematically detect unknown barcode/adapter sequences within a set of samples?
I have often downloaded datasets from the SRA where the authors failed to mention which adapters were trimmed during the processing.

Local alignments tend to overcome this obstacle, but it feels a bit barbaric.

fastQC works occasionally to pick them up, but sometimes fails to find the actual adapter sequences.

Usually, I ended up looking up the kits they used and trying to grep for all the possible barcodes.

Is there a more robust/efficient way to do this?

# Answer
> 1 votes
The `minion` utility from the kraken/reaper toolkit may be helpful for this: http://wwwdev.ebi.ac.uk/enright-dev/kraken/reaper/src/reaper-latest/doc/minion.html

# Answer
> 4 votes
You mention that FastQC *"fails to find the actual adapter sequences"* \- I guess you mean in the Adapter Sequence Contamination plot. However, the kmer and Sequence Content Plots are often useful even when the former fails. I've used these in the past - you can sometimes just read off the adapter sequence from the start of the Sequence Content Plot (or at least see how many bases to trim).

# Answer
> 2 votes
I'm not aware of any existing methods to do this, but here are a couple of ideas about how it might be done:

Canu has a method of adapter trimming which involves looking for the absence of overlap for reads. If there are no other reads which share sequence across a particular region, then the read is broken up at the point of low coverage, and small pieces are discarded. It would be possible to use a method like this to hunt for possible adapter/barcode sequences by preserving the short reads.

Another option is to do a kmer search at the start of reads, and see if any of the high-abundance kmers can be assembled together and/or matched to existing known adapters or barcodes.

# Answer
> 1 votes
If you happen to know a sequence that should be highly abundant in the library, you can grep its beginning or end (with pattern match highlighting) and see if the same sequence systematically comes just before or just after respectively. This kind of visual inspection can help you finding the adaptor.

For instance, in a previous lab, we were working on *D. melanogaster* small RNA sequencing data and my colleague knew from previous experience with this kind of data that the following small RNA was likely to be abundant: http://flybase.org/reports/FBgn0065042.html

We just had to grep it in the fastq file to see many lines with this sequence, next to another sequence that happened to be always the same: the unknown adapter.

---
Tags: barcode, adapter, trimming
---


# Question
Title: What is the difference between FASTA, FASTQ, and SAM file formats?
I'd like to learn the differences between 3 common formats such as FASTA, FASTQ and SAM. How they are different? Are there any benefits of using one over another?

Based on Wikipedia pages, I can't tell the differences between them.

# Answer
> 62 votes
Let’s start with what they have in common: All three formats store

1. sequence data, and
2. sequence metadata.

Furthermore, all three formats are text-based.

However, beyond that all three formats are different and serve different purposes.

Let’s start with the simplest format:

### FASTA

FASTA stores a variable number of sequence records, and for each record it stores the sequence itself, and a sequence ID. Each record starts with a header line whose first character is `>`, followed by the sequence ID. The next lines of a record contain the actual sequence.

The Wikipedia artice gives several examples for peptide sequences, but since FASTQ and SAM are used exclusively (?) for nucleotide sequences, here’s a nucleotide example:

```
>Mus_musculus_tRNA-Ala-AGC-1-1 (chr13.trna34-AlaAGC)
GGGGGTGTAGCTCAGTGGTAGAGCGCGTGCTTAGCATGCACGAGGcCCTGGGTTCGATCC
CCAGCACCTCCA
>Mus_musculus_tRNA-Ala-AGC-10-1 (chr13.trna457-AlaAGC)
GGGGGATTAGCTCAAATGGTAGAGCGCTCGCTTAGCATGCAAGAGGtAGTGGGATCGATG
CCCACATCCTCCA

```

The ID can be in any arbitrary format, although several conventions exist.

In the context of nucleotide sequences, FASTA is mostly used to store reference data; that is, data extracted from a curated database; the above is adapted from GtRNAdb (a database of tRNA sequences).

### FASTQ

FASTQ was conceived to solve a specific problem arising during sequencing: Due to how different sequencing technologies work, the confidence in each base call (that is, the estimated probability of having correctly identified a given nucleotide) varies. This is expressed in the Phred quality score. FASTA had no standardised way of encoding this. By contrast, a FASTQ record contains a sequence of quality scores for each nucleotide.

A FASTQ record has the following format:

1. A line starting with `@`, containing the sequence ID.
2. One or more lines that contain the sequence.
3. A new line starting with the character `+`, and being either empty or repeating the sequence ID.
4. One or more lines that contain the quality scores.

Here’s an example of a FASTQ file with two records:

```
@071112_SLXA-EAS1_s_7:5:1:817:345
GGGTGATGGCCGCTGCCGATGGCGTC
AAATCCCACC
+
IIIIIIIIIIIIIIIIIIIIIIIIII
IIII9IG9IC
@071112_SLXA-EAS1_s_7:5:1:801:338
GTTCAGGGATACGACGTTTGTATTTTAAGAATCTGA
+
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII6IBI

```

FASTQ files are mostly used to store short-read data from high-throughput sequencing experiments. The sequence and quality scores are usually put into a single line each, and indeed many tools assume that each record in a FASTQ file is exactly four lines long, even though this isn’t guaranteed.

As with FASTA, the format of the sequence ID isn’t standardised, but different producers of FASTQ use fixed notations that follow strict conventions.

### SAM

SAM files are so complex that a complete description <sup>\[PDF\]</sup> takes 15 pages. So here’s the short version.

The original purpose of SAM files is to store mapping information for sequences from high-throughput sequencing. As a consequence, a SAM record needs to store more than just the sequence and its quality, it also needs to store information about where and how a sequence maps into the reference.

Unlike the previous formats, SAM is tab-based, and each record, consisting of either 11 or 12 fields, fills exactly one line. Here’s an example (tabs replaced by fixed-width spacing):

```
r001  99  chr1  7 30  17M         =  37  39  TTAGATAAAGGATACTG   IIIIIIIIIIIIIIIII
r002  0   chrX  9 30  3S6M1P1I4M  *  0   0   AAAAGATAAGGATA      IIIIIIIIII6IBI    NM:i:1

```

For a description of the individual fields, refer to the documentation. The relevant bit is this: SAM can express exactly the same information as FASTQ, plus, as mentioned, the mapping information. However, SAM is also used to store read data *without* mapping information.

In addition to sequence records, SAM files can also contain a *header*, which stores information about the reference that the sequences were mapped to, and the tool used to create the SAM file. Header information precede the sequence records, and consist of lines starting with `@`.

SAM itself is almost never used as a storage format; instead, files are stored in BAM format, which is a compact, gzipped, binary representation of SAM. It stores the same information, just more efficiently. And, in conjunction with a search index, allows fast retrieval of individual records from the middle of the file (= fast random access). BAM files are also much more compact than compressed FASTQ or FASTA files.

---

The above implies a *hierarchy* in what the formats can store: FASTA ⊂ FASTQ ⊂ SAM.

In a typical high-throughput analysis workflow, you will encounter all three file types:

1. FASTA to store the reference genome/transcriptome that the sequence fragments will be mapped to.
2. FASTQ to store the sequence fragments before mapping.
3. SAM/BAM to store the sequence fragments after mapping.

# Answer
> 19 votes
In a nutshell,

`FASTA` file format is a DNA sequence format for specifying or representing DNA sequences and was first described by Pearson `(Pearson,W.R. and Lipman,D.J. (1988) Improved tools for biological sequence comparison. Proc. Natl Acad. Sci. USA, 85, 2444–2448)`

`FASTQ` is another DNA sequence file format that extends the FASTA format with the ability to store the sequence quality. The quality scores are often represented in ASCII characters which correspond to a phred score)

Both FASTA and FASTQ are common sequence representation formats and have emerged as key data interchange formats for molecular biology and bioinformatics.

`SAM` is format for representing sequence alignment information from a read aligner. It represents sequence information in respect to a given reference sequence. The information is stored in a series of tab delimited ascii columns. The full SAM format specification is available at http://samtools.sourceforge.net/SAM1.pdf

# Answer
> 7 votes
FASTA (officially) just stores the name of a sequence and the sequence, unofficially people also add comment fields after the name of the sequence. FASTQ was invented to store both sequence and associated quality values (e.g. from sequencing instruments). SAM was invented to store alignments of (small) sequences (e.g. generated from sequencing) with associated quality values and some further data onto a larger sequences, called reference sequences, the latter being anything from a tiny virus sequence to ultra-large plant sequences.

# Answer
> 3 votes
FASTA and FASTQ formats are both file formats that contain sequencing reads while SAM files are these reads aligned to a reference sequence. In other words, FASTA and FASTQ are the "raw data" of sequencing while SAM is the product of aligning the sequencing reads to a refseq.

A FASTA file contains a read name followed by the sequence. An example of one of these reads for RNASeq might be:

```
>Flow cell number: lane number: chip coordinates etc.
ATTGGCTAATTGGCTAATTGGCTAATTGGCTAATTGGCTAATTGGCTAATTGGCTAATTGGCTA

```

The FASTQ version of this read will have two more lines, one + as a place holder and then a line of quality scores for the base calls. The qualities are given as characters with '!' being the lowest and '~' being the highest, in increasing ASCII value. It would look something like this

```
@Flow cell number: lane number: chip coordinates etc.
ATTGGCTAATTGGCTAATTGGCTAATTGGCTAATTGGCTAATTGGCTAATTGGCTAATTGGCTA
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65

```

A SAM file has many fields for each alignment, the header begins with the @ character. The alignment contains 11 mandatory fields and various optional ones. You can find the spec file here: https://samtools.github.io/hts-specs/SAMv1.pdf .

Often you'll see BAM files which are just compressed binary versions of SAM files. You can view these alignment files using various tools, such as SAMtools, IGV or USCS Genome browser.

As to the benefits, FASTA/FASTQ vs. SAM/BAM is comparing apples and oranges. I do a lot of RNASeq work so generally we take the FASTQ files and align them the a refseq using an aligner such as STAR which outputs SAM/BAM files. There's a lot you can do with just these alignment files, looking at expression, but usually I'll use a tool such as RSEM to "count" the reads from various genes to create an expression matrix, samples as columns and genes as rows. Whether you get FASTQ or FASTA files just depends on your sequencing platform. I've never heard of anybody really using the quality scores.

---
Tags: fasta, fastq, file-formats, sam
---


# Question
Title: Why is SMALT better for microbial genomics than other mappers?
SMALT seems to be one of the most used read mappers for bacterial data, see, e.g., this query. I do not say that it is not a great mapper, but I cannot easily see what are its main strengths compared to mappers such as BWA-MEM, Bowtie2, NovoAlign or GEM. Moreover, it is not even published.

Could you name some of its distinguishing features (e.g., user support by Sanger Pathogens)?

So far I have heard only arguments like "We use SMALT because everyone does it.", but this is not convincing enough for me.

# Answer
I have used smalt a couple of times. Smalt is more flexible than bwa/bowtie and can be tuned to be more sensitive to divergent hits, which is useful to certain applications. I heard Sanger did evaluate a few mappers for some non-typical applications (not sure what). They found smalt to perform better. Also smalt was developed at Sanger. I guess they could modify smalt based on their needs, though I don't know the details.
> 3 votes

---
Tags: read-mapping
---


# Question
Title: What are all the reference files produced by bwa index, and are these dependent upon whether the reference is zipped?
I have indexed a gzipped reference with bwa: `bwa index reference.fa.gz`, which produces a series of other files `reference.fa.gz.{amb,ann,bwt,pac,sa}`. These are working fine with bwa alignment.

I have discovered that samtools does not take a gzipped reference, so I am planning to use an unzipped version of the reference for my workflow instead of dealing with two separate representations of the reference.

Are the files `reference.fa.gz.{amb,ann,bwt,pac,sa}` dependent upon the reference being gzipped? Do I need to reindex the unzipped reference.fa with bwa, or can I just rename the current files to remove the `.gz` portion of the filename?

# Answer
You'll get the exact same index (the `amb`, `ann`, `bwt`, `pac` and `sa` files) whether the reference is gzipped or not. BWA also makes its own packed reference sequence (the .pac file) so you don't even need the genome around after you index.
> 7 votes

---
Tags: bwa, samtools, reference-genome
---


# Question
Title: 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?

# Answer
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:

1. 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.
2. 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
3. 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.
4. 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."
5. 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.
> 16 votes

---
Tags: sequence-alignment, sam
---


# Question
Title: 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?

# Answer
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.
> 7 votes

# Answer
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.
> 4 votes

---
Tags: sam
---


title: Do you use SNAP for short-read mapping?

text:
# Question
Title: 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!

# Answer
> 5 votes
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.

---
Tags: bwa
---


title: Stable download URLs

text:
# Question
Title: Stable download URLs
One big problem that I'm regularly facing is that URLs for downloading Bioinformatics data (e.g., RefSeq releases or NCBI genome releases) disappear.

Does anyone have any good solution for this?

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

---
Tags: data-download
---


# Question
Title: 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"?

# Answer
> 10 votes
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.

---
Tags: bwa, sam
---


# Question
Title: Is there any publicly available multi-omics dataset?
I am searching for a multi-omics dataset which may include genomics, transcriptomics and proteomics (eg. snyderome). I need such a dataset for an introductory data-exploration purpose. So it can be from any host organism and any experimental methods/conditions.

Is there any database/repository where I can find it?

# Answer
> 10 votes
There is a database called OmicsDI, where one can search for multi-omics datasets.

Here's a link of the associated publication (Perez-Riverol, Yasset, et al. "Omics Discovery Index-Discovering and Linking Public Omics Datasets." bioRxiv (2016): 049205.) for more details.

# Answer
> 3 votes
There are several datasets available on GEO, though you do have to search for them. For example, here are three data sets that have both Illumina methylation and gene expression microarray profiling:

# Answer
> 2 votes
Next to OmicsDI the EBI has a special repository for multi-omics datasets: https://www.ebi.ac.uk/biosamples

It links the different datasets between repositories, ie. PRIDE for MS/MS based data and ArrayExpress for RNASeq data.

---
Tags: multi-omics
---


# Question
Title: Expected allele frequency distribution of SNVs in real NGS data
I have a huge amount of ~20x human WGS samples, aligned, and all SNVs that were called with GATK under standard germline parameters set.

What I need to do is to model SNVs Allele Frequency (AF) for different underlying Copy Numbers. I'd better provide a toy example. For particular genomic region X:

*If X is presented by **2** copies for the particular samples, we expect AF to be super-close to 1 or to 0.5.*

*If X is presented by **4** copies, I expect any particular AF to be close to 0.25, 0.5, 0.75 or 1.*

Of course, I can use *Binomial Distribution* for these purposes. However, as we know, the distribution is not exactly Binomial due to alignment/sequencing biases and the median AF for all heterozygous SNVs is more close to 0.48 but not to 0.5 as we would expect. Another thing: for high copy numbers we expect higher coverages. And GATK use several filters so I suppose that we will not see SNVs with AF like 0.125 (in case if the segment has ploidy 8) - despite the super high coverage there GATK may reject this "weird" AF.

I have read several papers that model SNVs AFs (and I agree that Beta Binomial Distribution may be quite accurate), however, I was not convinced enough that I should use the particular modelling. From your experience (in case if you do SNVs calling), which probabilistic distribution should I use? How should I estimate parameters for each of them (should I expect for CN4 AF=0.5 more frequent than AF=0.75 or vice versa, how to estimate this from data)?

**UPD:** For simplicity we can say that we have a lot of previously identified regions with ploidy different from CN2, and I can take these coordinates from here. So I can use more or less "supervised" learning for parameters' estimation.

# Answer
> 4 votes
I don't have enough experience to answer which probabilistic distribution should be used.

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

---
Tags: snv, modelling, statistics
---


# Question
Title: 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.

# Answer
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).
> 2 votes

---
Tags: sam, linux, bedtools
---


# Question
Title: 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.

# Answer
> 7 votes
`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.

---
Tags: ngs, read-mapping, bwa
---


# Question
Title: 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?

# Answer
> 9 votes
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.

---
Tags: terminology, blast, read-mapping
---


# Question
Title: Finding peaks and estimating cell population sizes in multi-dimensional flow cytometry data
Our research institute processes a lot of flow cytometry data, but the produced data is under-utilised due to the effort required to process it. A typical run will produce 5 million events (ideally one event per cell), with up to 14 dimensions of \[ideally\] log-normal fluorescence values for particular groups of cells ("populations"). Due to various systematic errors, negative values, scatter, and a non-zero "zero" value can happen, but I'm going to ignore those for the purpose of this question and assume that the data are well-distributed.

Researchers will typically probe these data using manual filters (set up specifically for each experiment) to find populations of interest. I suppose a picture might help. This one shows three identifiable cell populations, in order of size one at about (X/Ly6C:2,Y/CD86:2), one at about (4.5,3), and one small population at (2,4).

Here's another plot that has two cell populations that are close to each other, such that the population "humps" overlap substantially.

Manual filters are typically used because it can be very difficult to distinguish between a noisy data point from a large cell population and a less noisy data point from a small cell population, particularly when considering populations that make up about 0.01% of the total cells.

As an additional complication, counting these populations can be difficult when populations overlap (as in the second image). A filter/slice through the plot that separates populations could count many cells as being members of the wrong population.

Can these populations be detected in an automated way? If it is assumed that the populations are spread in a gaussian fashion at some point in each dimension, is there some method that can be used to approximate the number of cells in each population, even when populations are close by?

# Answer
There are quite a few algorithms developed for the automatic classification of multidimensional flow cytometry data, you can see a (not so recent?) review here: http://www.nature.com/nmeth/journal/v10/n3/full/nmeth.2365.html. In your case, you are interested in the **unsupervised** methods, since you do not have data from single populations (which would provide you "labels" associated to each sample), as I understood.

In other words, if you have your data set as a matrix of 5 million samples (cells) with 14 dimensions (FC values) each, then you are looking basically for algorithms that will do a *clustering* of the samples into groups. By the way your description of each sample belonging to a multivariate Normal distribution leads to a **Gaussian Mixture Model**, which can be estimated from your data using e.g. the scikit-learn library for python: http://scikit-learn.org/stable/modules/mixture.html
> 2 votes

# Answer
Automated gating methods are gaining popularity over manual and time consuming analysis in say FlowJo, you should take a look at the openCyto bioconductor package this is a framework which builds on top of existing Bioconductor infrastructure for flow cytometry (PLOS compbio paper) with the aim of making various analysis routes more accessible to the researcher by assisting getting data into and out various gating algorithms. You can find a nice intro here which starts with importing a manual gating scheme from FlowJo and a how-to on 1D and 2D gating can also be found here.

An excellent YouTube intro to analysing flow data in R by Ryan Duggan Cytometry on Air is worth watching especially for those which have no R experience and are curious as to what can be achieved. Additionally it's worth having a look at the growing list of Bioconductor packages and datasets that exists for flow data: https://www.bioconductor.org/help/search/index.html?q=cytometry/
> 2 votes

---
Tags: flow-cytometry, peak-calling
---


title: Generic HMM solvers in bioinformatics?

text:
# Question
Title: Generic HMM solvers in bioinformatics?
Hidden Markov models (HMMs) are used extensively in bioinformatics, and have been adapted for gene prediction, protein family classification, and a variety of other problems. Indeed, the treatise by Durbin, Eddy and colleagues is one of the defining volumes in this field.

Although the details of each of these different applications of HMMs differ, the core mathematical model remains unchanged, and there are efficient algorithms for computing the probability of the observed sequence given the model, or (perhaps more useful) the most likely hidden sequence given the sequence of observed states.

Accordingly, it seems plausible that there could be a generic software library for solving HMMs. As far as I can tell that's not the case, and most bioinformaticians end up writing HMMs from scratch. Perhaps there's a good reason for this? (Aside from the obvious fact that it's already difficult, nigh impossible, to get funding to build and provide long-term support for open source science software. Academic pressures incentivize building a new tool that you can publish a paper on much more than building on and extending existing tools.)

Do any generic HMM solver libraries exist? If so, would this be tempting enough for bioinformaticians to use rather than writing their own from scratch?

# Answer
> 7 votes
I would also recommend to take a look at pomegranate, a nice Python package for probabilistic graphical models. It includes solvers for HMMs and much more. Under the hood it uses Cythonised code, so it's also quite fast.

# Answer
> 5 votes
There are certainly software libraries for working with HMMs. For a general-purpose implementation in C++, take a look at the SeqAn HMM algorithms.

For your purposes, i.e. “computing … the most likely hidden sequence given the sequence of observed states”, you’d invoke `viterbiAlgorithm` with your observed sequence and the HMM graph.

More fundamentally I think that most existing, mature implementations are probably found in the domain of signal processing, which has been using them longer than biology, and where most of the underlying theory was developed.

# Answer
> 4 votes
If I remember correctly Ewan Birney's Dynamite (a compiler-compiler) as presented at ISMB 1997 had this functionality, there is also some code here on GitHub https://github.com/birney/wise3 which at least mentions Dynamite. Suspect Ewan is too busy these days to work on this, although he has tweeted about blowing dust of his old Dynamite, sorry Dynamite code: https://twitter.com/ewanbirney/status/788121636973142016

---
Tags: hidden-markov-models, sequence-analysis
---


# Question
Title: Tumor purity/contamination/admixture estimation
Can anyone recommend a good tool for estimating the tumor content given a matched tumor and normal file for DNA NGS whole genome sequencing data or whole exome data?

Is it possible to estimate this without a normal sample as well?

# Answer
I have previously estimated tumour purity with the EXPANDS an inferred tumour heterogeneity program which is designed to calculate the number of clonal subpopulations in matched tumour/normal samples. The purity is essentially the size of the largest subpopulation identified in that sample - this is discussed in the programs FAQ. In addition to a matched tumour / normal exome or genome sequencing you'll also need match somatic copy number as a seg file as input. Other programs also exist for this sort of inferred heterogeneity analysis - some of which may also give you a measure of purity are: Absolute, ThetA, SciClone, CHAT, PyClone and Canopy. A more complete list looks to be here.

The only other thing I'd suggest is have an estimate via some orthogonal measure of the purity to try at least judge how the different methods are performing with you data, *e*.*g*. you might have clonality indications from cytogenetics/FISH or histology work, or perhaps FACS. Picking samples known to be pure / impure via one or more of these might help you get a handle on how well various methods are performing.

As regards estimation without a normal sample the only program which looks like it might help is QuantumClone, but that requires more than one tumour sample form said patient either spatially or temporally distinct to perform the analysis.

Also If you have low coverage DNA sequencing CNAnorm looks like it could be useful with just a single sample.
> 3 votes

# Answer
It's usually CNV callers that make use of Tumour/Normal WGS pairs to estimate purity. It can also be done with WES (exome) Tumour/Normal pairs.

There are several tools out there, I have some experience with the one written by Illumina (public on Github):

https://github.com/Illumina/canvas

It requires realigning things with bowtie2, so I don't think it can take existing data aligned with other aligners. Not sure why.
> 2 votes

# Answer
Celluloid emits copy number profiles and tumour purity / ploidy information. There's a nice tutorial https://github.com/mathieu-lemire/celluloid\_0.11\_tutorial. One thing to keep in mind is that celluloid will find multiple solutions and taking advantage of its plotting tools is essential to determining which solution is correct. Typically the first solution is correct, but occasionally it may not be.

Another tool similar is TITAN, which may also require hand annotated the ideal solution.
> 1 votes

---
Tags: cancer, wgs
---


# Question
Title: Detecting structural variants with MinION data
Working on various cancers I have an interest in detecting structural variation (SV) in human, we've successfully used various tools like Pindel, SVDetect, Manta, and LUMPY, to name a few for detecting SVs in illumina short-read sequencing. I curious if anyone has successfully used ONTs MinION sequencer for detecting SV, as there are many cases where longer reads would be beneficial for SV detection especially where events are long and/or occur in repetitive regions. Has anyone tried or had success with tools previously designed for PacBio data such as Sniffles?

Supplementary question, obviously the MinION has not quite got the throughput of say the PacBio SEQUEL so can't generate low-pass coverage for the whole human genome as easily, has anyone got experience trying to generate reads targeting specific areas of the genome in areas that are known for translations, long inversions, duplications, *etc*?

# Answer
There was a Structural Variant breakout session at the London Calling conference this year. Unfortunately I didn't attend that session, but MinION community members have access to Constance Donnell's summary of that:

https://community.nanoporetech.com/posts/breakout-structural-varia

Here are my attempts at grabbing non-creative chunks from those notes:

* Professor Wigard Kloosterman has been developing a bioinformatics pipeline called NanoSV for mapping genomic structural variants in patients with congenital abnormalities
* Tomas Sesani at the University of Utah studied the precise copy number of duplicated genes in individual virus genes over time, and tracked the SNP within copy number variable regions of the genome \[no public tool mentioned\]
* Dr. Sudha Rao from Genotypic Technology in India used the MinION sequencer with Sniffles for structural variation calling

ONT have said that the sessions/talks from London Calling 2017 will all be made publicly available at some point in the future.
> 5 votes

---
Tags: nanopore, structural-variation, long-reads, cancer
---


# Question
Title: Filtering imputed GWAS SNPs based on a MAF difference of 10%
There are many posts on the web regarding QC steps pre and post-imputation.

Does applying below (new?) 10% MAF difference rule make sense, pitfalls?

Here is the process:

1. Get MAF for imputed set, using SNPTEST with flag `-summary_stats_only`
2. Convert imputed set to hard-calls using gtools with flag `--threshold 0.9`
3. If the MAF from step 1 and step 2 differs more than 10% than exclude the variant.

More info:

* GWAS is 50K vs 50K case control samples.
* This step is applied after `info > 0.4` filter.

# Answer
I cannot think of any principled rationale for choosing this filtering strategy.

However, I am going to take a guess that this filtering strategy is supposed to filter out SNPs for which imputation did not work well? In that case the appropriate statistic to filter on is the INFO score as described here.

You might consider picking a higher threshold than 0.4 if you feel that you aren't filtering out enough variants.
> 3 votes

---
Tags: gwas, imputation, qc
---


# Question
Title: How can I speed up INDEL calling/correction on BAM files?
The `samtools mpileup` command has quite a neat feature that it is able to correct mapping errors associated with misalignment of INDELs. By default, the `mpileup` command will not work for reads that have more than 250X coverage of the reference genome. While this limit can be increased, very high coverage causes the mpileup program to grind to a halt, so it'd be nice to know if there's some easy way to make that faster.

To add a bit more context, I've been doing this with mitochondrial genome reads that were extracted from both Illumina whole-genome sequencing (coverage ~1000X), and from targeted amplicon sequencing done on the IonTorrent (coverage up to ~4000X).

I see that @rightskewed has mentioned the downsampling ability of samtools with `samtools view -s <float>` (see here), which seems like it might work as a solution for this if used prior to the mpileup operation.

# Answer
> 4 votes
I wasn't aware of the samtools subsampling when I had this problem a couple of years ago, so ended up writing my own digital normalisation method to deal with *mapped* reads. This method reduces the genome coverage, but preserves reads where coverage is low.

Because I was working with IonTorrent reads (which have variable length), I came up with the idea of selecting the longest read that mapped to each location in the genome (assuming such a read existed). This meant that the highly variable coverage for different samples (sometimes as low as 200X, sometimes as high as 4000X) was flattened out to a much more consistent coverage of about 100-200X. Here's the core of the Perl code that I wrote:

```
  if(($F[2] ne $seqName) || ($F[3] != $pos) || (length($bestSeq) <= length($F[9]))){
    if(length($bestSeq) == length($F[9])){
      ## reservoir sampling with a reservoir size of 1
      ## See https://en.wikipedia.org/wiki/Reservoir_sampling
      ## * with probability 1/i, keep the new item instead of the current item
      $seenCount++;
      if(!rand($seenCount)){
        ## i.e. if rand($seenCount) == 0, then continue with replacement
        next;
      }
    } else {
      $seenCount = 1;
    }
    if(($F[2] ne $seqName) || ($F[3] != $pos)){
      if($output eq "fastq"){
        printSeq($bestID, $bestSeq, $bestQual);
      } elsif($output eq "sam"){
        print($bestLine);
      }
    }
    $seqName = $F[2];
    $pos = $F[3];
    $bestLine = $line;
    $bestID = $F[0];
    $bestFlags = $F[1];
    $bestSeq = $F[9];
    $bestQual = $F[10];

```

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

---
Tags: sam, samtools, variant-calling, indel
---


# Question
Title: What causes the difference in total length of assembled contigs and scaffolds in SOAPdenovo2?
I use SOAPdenovo2 to assemble a large genome (4.8G) using ~20X paired-end reads. The total length of contig sizes is 6.3G while total length of scaffolds is 2.7G. Note that this is a haploid genome, so there is no issue of heterozygosity for scaffolding.

I am wondering what happened during the scaffolding procedure? How comes that sum of scaffolds is so much smaller than sum of contigs?

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

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

---
Tags: assembly, scaffold
---


# Question
Title: Pseudo-temporal ordering in heterogeneous populations
What's the most widely accepted tool for doing pseudo-temporal ordering from scRNAseq data? Also is there away to separate differential expression that occurs based on "cell identity" or maybe more accurately cell type fate from that which arises from cells being in different stages differentiation.

To be more concrete, lets say there's a population of cells, some of which were born at time 1, time 2, and time 3. The progression along the temporal trajectory can be described via a set of genes that are fluctuating as the cell matures. So you might have the same cell type which was born at time 1 be transcriptionally distinct from a younger one born at time 3. On the other hand, within this population there are subpopulations which will have different cell fates and are transcriptionally distinct. Is there away to reliably separate the temporal axis from the "cell fate axes". If not is this something people are working on or is flawed logic to think this kind of thing is possible?

# Answer
Regarding the first part of the question: scRNA-seq is a rapidly developing field so it may be hard to talk about "widely accepted tool for doing pseudo-temporal ordering from scRNAseq data". Few of the tools aiming to do just that include Monocole, Waterfall or Sincell (see this paper for references https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4122333/)

The second part is a bit more complex. Many aspects have come together here, e.g. cell-cycles, sub-populaitons identification and pseudo-temporal ordering, to get a true reflection of the biological processes. There are efforts on all (and more) these fronts (see again e.g. the above paper, though not the latest) and there are most likely people working on their integration. I'm not aware of any studies published yet at this depth though
> 3 votes

---
Tags: scrnaseq
---


title: How to read structural variant VCF?

text:
# Question
Title: 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?

# Answer
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.
> 8 votes

# Answer
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.
> 4 votes

---
Tags: vcf, variation, structural-variation
---


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

# Answer
> 2 votes
I can see that there would be a three step process to doing this:

1. Merge counts from all samples in the group and then resample pseudo-replicates from this. If `x` is a matrix with samples in the group being columns and genes being rows

```
s <- sample(row.names(x), n = mean(colSums(x)), probs=rowSums(x)/sum(rowSums(x))
stab <- table(s)
s <- as.vector(stab)
names(s) <- names(stab)

```

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

---
Tags: r, edger
---


title: How to compare groups using WGS data?

text:
# Question
Title: How to compare groups using WGS data?
We have whole genome sequencing data for patients (not-cancer) (n=60) and for healthy controls (n=20). The sequencing centre has provided us with the best practice bioinformatics analyses including reads mapping (.BAM) and variant calling using GATK (.vcf) as well as annotation (annotated .vcf and .gVCF).

What should be our next steps? We are interested to see if there are any differences (global and/or specific) between the groups.

# Answer
> 2 votes
I'm not familiar with the program, but apparently Hail is setting itself up as a swiss-army chainsaw project for doing downstream analysis on variant-called datasets.

An overview of Hail can be found here:

http://blog.cloudera.com/blog/2017/05/hail-scalable-genomics-analysis-with-spark/

A tutorial on association testing can be found here:

https://hail.is/hail/tutorial.html#Association-testing

# Answer
> 2 votes
If you have gVCFs, the first thing you should try is joint variant calling. According to GATK, joint variant calling "empowers variant discovery by providing the ability to leverage population-wide information from a cohort of multiple sample\[sic\], allowing us to detect variants with great sensitivity and genotype samples as accurately as possible." source.

Once you have two sets of high-quality variants, what you do next depends on your research question. Are you looking for druggable mutations? Underlying causes? Biomarkers? Patient prognosis? You can look at the most frequently mutated positions in the patient cohort and compare them to your reference cohort, check for co-occurring mutations, cluster them, do principal component analysis, do some machine learning to stratify them, etc.

---
Tags: wgs
---


# Question
Title: 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)?

# Answer
> 7 votes
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.

# Answer
> 3 votes
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.

# Answer
> 2 votes
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.. :)

---
Tags: chip-seq
---


# Question
Title: What methods exist to calculate RNA expression profile similarity
Some of the work in our lab requires a comparison of a strain across several experimental conditions. We are looking to identify most similar experimental conditions based on the gene transcription response similarity from the cell.

While we could easily invent and create home-grown methods to do it, their implementation and testing are a laborious project in themselves and are outside the scope of our current work.

Are there any methods for RNA expression profile similarity computation that has been already published? If yes, what is your experience using them?

# Answer
There's a new paper that's just appeared on the subject of "proportionality", including a method by which RNA expression might be compared.

This is a new concept for me, and the article is not easy enough for me to read that I can write up a quick summary; the authors don't seem to devote a section in the paper to defining "proportionality". However, here's an interesting chunk from the article:

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

# Answer
The suggestion to use proportionality is probably the correct one if you are interested in similar *patterns* between samples. But not if you are interested in *absolute* differences.

For example: the following two samples are similar in pattern, but not similar on absolute levels:

```
        Sample 1    Sample 2    Sample 3
Gene A    10          100           80
Gene B     8           80          100
Gene C    12          120          120

```

Samples 1 and 2 have a perfect proportionality (phi is 0) and also a perfect correlation (as a side note any pair with a perfect proportionality will always have a perfect correlation). However, in terms of logfold changes, samples 2 and 3 are more like each other.

Of course in real life you would never see a comparison like the sample 1 - sample 2 one because normalisation would have removed the scale difference. This was exactly the point brought up by the proportionality paper. But normalisation methods don't normally guarantee that the sum of expression for each sample is identical, and such differences can still occur.

An alternative that might suit more in the second case is either the euclidean distance between the samples, or the euclidean distance on the first two components of a principle component or multi-dimensional scaling. The latter is effectively using mean logFC between samples.

Assuming that `x` is a matrix containing normalised, log transformed expression values, you could use R and limma to calculate distance in multi-dimensional scaled space as follows:

```
library(limma)
mds <- plotMDS(x, plot=FALSE)
mds <- data.frame(mds$x, mds$y)
distances <- dist(mds)

```
> 5 votes

# Answer
If you have a lot of samples with very different environmental conditions, then a Weighted Gene Correlation Network Analysis (WGCNA) might be appropriate.

This type of analysis looks for genes that trace similar (or opposing) expression patterns throughout different conditions (e.g. high-medium-medium-low-absent-high would be highly negatively correlated with low-medium-medium-high-very high-low).

That particular paper introduces the concept of "modules", which are groups of genes that share similar expression patterns. Functions are available for plotting how the expression of canonical module members changes throughout the different conditions, and for identifying which module (or modules) a particular gene is likely to be a member of.

WGCNA works best when there's a lot of different expression changes in the different conditions, which sounds like it would fit well with your project. However, it concentrates more on the genes, rather than the conditions (which seems like it would be less useful for you).
> 3 votes

# Answer
The euclidean distance is probably the simplest, both conceptually and in terms of implementation. It's by no means an elegant solution, and it may not perform well in certain circumstances.

Euclidean distance is easiest to conceptualize as the distance between two points in a two dimensional space.

```
Y
^
|
|
|
|            * p = (3, 3)
|
|
|
|        * q = (2, 1)
|
----------------------------> X

```

In this example, the distance between the two points is

```
d(p, q) = sqrt( (p_x-q_x)^2 + (p_y-q_y)^2 )
        = sqrt( (2-3)^2 + (1-3)^2 )
        = sqrt(5)
        ≈ 2.24

```

For a gene expression profile with two genes, this is exactly how the euclidean distance would be calculated, using expression values from one gene as the X axis and expression values from the other gene as the Y axis. Realistically, though, gene expression profiles typically contain thousands or tens of thousands of genes, so instead we use the generalization of the distance calculation for N dimensions.

```
d(p, q) = sqrt( (p_1-q_1)^2 + (p_2-q_2)^2 + ... + (p_N-q_N)^2 )

```

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

---
Tags: rna-seq
---


# Question
Title: How are snakemake's --cluster and --drmaa options implemented?
I'm fairly new to snakemake and I'm trying to understand the difference between the `--cluster` and `--drmaa` flags, both of which are used to submit jobs to compute clusters/nodes.

The docs give a few hints about the advantages of using `--drmaa` here:

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

And here:

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

So I have a conceptual understanding of the advantages of using `--drmaa`. However, I don't consider the above a very complete explanation, and I don't know how these flags are implemented in snakemake under the hood, can someone elaborate?

Note, that although this could be considered a more general programming question, snakemake is predominantly used in bioinformatics; I was persuaded that this question would be considered on-topic because of this meta-post and this answer.

# Answer
> 5 votes
I'd always kind of wondered how this worked too, so I took this as an excuse to look into the snakemake code. At the end of the day this becomes a question of (1) how are jobs actually submitted and (2) how is it determined if jobs are done (and then whether they failed)?

For DRMAA, python has a module (appropriately named "drmaa") that wraps around the libdrmaa library that comes with most schedulers. This is a very popular route, for example the Galaxy project uses this for dealing with most clusters (e.g., I use it to connect our internal Galaxy instance to our slurm cluster). The huge benefit here is that DRMAA does magic to allow you to simply submit commands to your cluster without needing to know whether you should execute `qsub` or `srun` or something else. Further, it then provides methods to simply query whether a job is running or not and what its exit status was.

Using the `--cluster` command requires a LOT more magic on the side of snakemake. At the end of the day it creates a shell script that it then submits using the command you provided. Importantly, it includes some secret files in that script that it can then watch for (ever noticed that `.snakemake` directory where you execute a script? That appears to be among the things it's used for.). These are named `{jobid}.jobfinished` and `{jobid}.jobfailed` and one of them will get touched, depending on the exit status of your command/script. Once one of those is there, then snakemake can move on in its DAG (or not, if there's a failure). This is obviously a LOT more to keep track of and it's then not possible for snakemake to cancel running jobs, which is something one can do easily with DRMAA.

# Answer
> 4 votes
**DRMAA** (Distributed Resource Management Application API) appears to be an open API that describes a specification for submission and management of work submitted to some grid/cluster. If your scheduler is DRMAA compliant, my guess would be that using `snakemake`'s `--drmaa` flag will afford you the additional controls exposed by that API.

As mentioned in your question with DRMAA support, a Ctrl+C on one's console would be passed on to the grid to kill or stop submitted jobs. A task that would otherwise require use of the `qdel` command on `Sun Grid Engine`, for example.

I could see why the suggestion would be to enable `--drmaa` where possible as it's likely to expose functionality (such as control of jobs, and availability of logs and errors) that makes the submission and management of jobs a little easier for you and `snakemake`. Perhaps someone with `snakemake` experience can explain a little more, but I hope this helps?

For what it's worth, I use Sun Grid Engine without DRMAA and it's quite a pain to get good information out of `qstat`, and its `qacct` log format is quite possibly the worst file format I have ever encountered. I suspect DRMAA provides a nice (or at least more reasonable) API for `snakemake` (and others) to get at this information more readily.

**Appended:**

Under the hood, the flags are parsed by Python's `argparse` module in the `__init__.py`. A mutually exclusive group of options forces a selection of `--cluster` or `--drmaa` (or `--cluster-sync`).

`--cluster` or `-c` requires you to specify an argument; the command for which to submit a job on your cluster (the example names `qsub`). The `--drmaa` flag just seems to indicate to `snakemake` that DRMAA is to be used, which doesn't change a lot in terms of how the command you run from your console is processed in the `__init__.py` or `workflow.py`.

However, when it is time to interface with your scheduler (as specified in `scheduler.py`), an `elif` statement checks whether you're in some form of cluster mode and the `else` catches the case where you've raised the `--drmaa` flag instead.

Here's where the magic happens, as now your job is submitted with the `DRMAAExecutor` as specified in executors.py, instead of the default executor.

Without investigating much more, I can see that the `DRMAAExecutor` features some more class attributes and exposes additional functions when compared to the others.

Hope this is of more use? Again, I've only just looked at the code myself but in lieu of a `snakemake` user chiming in, I thought I'd have a nose.

---
Tags: snakemake
---


title: Is spark widely used in bioinformatics?

text:
# Question
Title: 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.

# Answer
> 8 votes
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.

# Answer
> 6 votes
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.

# Answer
> 6 votes
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).

---
Tags: gatk
---


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

# Answer
EPIC data can be processed in the same manner as the previous iteration of methylation array data from Illumina (450k). This means that starting with .idat files, normalization should be performed (for example, via the minfi package). A recent paper from the creators of minfi is particularly helpful because it makes clear that normalized EPIC data from their package can be immediately compared against, for example, level 3 TCGA data.

After that, I suggest using the manifest to attach genomic coordinates to your probes and segregating them into functional regions. By testing differential methylation in only regulatory regions, for example, you can increase the statistical power by reducing the overall number of tests to the ones you expect to yield major differences.

There are existing packages out there for differential methylation analysis, but without knowing your replicate structure or aims, it is difficult to point you in the right direction.
> 6 votes

---
Tags: r, normalization, methylation, microarray
---


# Question
Title: State of the art in predicting Translation Initiation Sites
I'm working on a university project of predicting Translation Initiation Sites in human DNA. I searched the net for papers and documentation to get guidelines and inspiration, but I feel uncertain that I was able to find the state of the art in solving this problem, so I thought of asking the community here: What's the state of the art method of predicting Translation Initiation Sites?

(I'm also posting the question to help this nascent Q&A site by participating, please close it if it completely misses the mark of being a good question.)

# Answer
Ribosome Profiling experiments profile the positions of ribosome-protected mRNA. a cycloheximide treatment fixes the translating ribosomes. A naive approach to define Translation Initiation Sites (TIS) would be to look at positions where you see distinct 'peaks' in the distribution of these fragments since the ribosome density will be higher at the start of the CDS (at the TIS).

TISDB is a database comprised of GTI-Seq global mapping of TIS codons at nearly single-nucleotide resolution. GTI-Seq uses two 'fixers' that fix ribosomes where one of them specifically fixes the initiating ribosomes.
> 3 votes

---
Tags: human-genome, phylogenetics
---


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

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

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

# Answer
> 3 votes
You can use centrifuge with the NT library to profile the taxonomy, and then remove the reads from eukaryotes.

# Answer
> 1 votes
*Answer from /u/Romanticon on reddit. The original answer can be found here*

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

---
Tags: read-mapping, metagenome
---


# Question
Title: Building STAR Genome Index for nanopore RNA sequencing
I am aligning a dataset of 1,000,000 reads oh human mRNA sequenced on Oxford Nanopore Technologies' MinION, and would like to use the STAR aligner, using the parameters recommended by Pacific Biosciences for long reads.

According to this Google Groups thread, in setting up the genome index for short reads, the parameter `sjdbOverhang` should be set to 1 less than the read length. Obviously, with long (mean 1.7Kb, max \>50Kb) reads, this doesn't make sense.

Running `STARlong --runMode genomeGenerate` without setting `sjdbOverhang` sets the parameter to `0`. Does anyone know what this means, how it might affect mapping, and what I should set it to for my long reads?

# Answer
> 8 votes
The parameter is used to determine how much sequence STAR indexes on each side of a splice junction to improve its alignment accuracy. For very long reads, this may not be ideal. I am not sure if STAR is capable of including multiple splice junctions since a long read is more than likely to span more than one.

It may be worthwhile to consider aligning your reads directly to transcripts and taking the reads that align poorly and using them to look for novel transcripts.

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

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

*\[emphasis preserved from the original text\]*

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

---
Tags: rna-seq, sequence-alignment, nanopore, star
---


title: What's the scaling for HOMER metagenes?

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

What am I missing?

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

---
Tags: ngs, bed, metagenome, homer
---


title: Entrez or Ensembl gene IDs?

text:
# Question
Title: Entrez or Ensembl gene IDs?
I'm using several datasets that are encoded using either Entrez or Ensembl IDs to specify genes, and need to decide on which to standardise on.

Are there any major reasons to use one over the other?

Which tools and conversion tables are best to use for this?

# Answer
I think there is no general answer to your question.

Which type of ID you will use depends on how practical it will be for your downstream processing. This depends on what other source of information you have, how much conversion work you will need to do in order to match your data to these information if you use one or the other.

I think the choice may emerge by trial and error, as you realize that a given choice is actually not convenient in your particular use case. Try to set up easy to use conversion procedures that you can easily insert in some steps of your workflow, so that you can easily switch from one choice to the other.

Concerning the tools to do the conversion, see the following question and its answers: Converting gene names from one public database format to another
> 2 votes

---
Tags: conversion, gene, public-databases
---


# Question
Title: R packages for data analyses of pooled CRISPR screens
We are designing a CRISPR/Cas9 experiment and thinking of the down-stream data analyses.

Are there any R packages for analysing raw NGS read count data from pooled genetic screens using CRIPSR/Cas9 to disrupt gene expression in a population of cells? I guess we will need to start with basics, i.e. sequences processing, data exploration, visualisation etc.

# Answer
> 5 votes
Assuming you mean CRISPR screens targeting many loci (for example using the GeCKO library) there is an R package here.

It uses a linear mixed effect model to compare guide counts before and after a selection step, allowing for multiple guides/gene and multiple replicates. It can also do the initial read alignment AFAIK and the author should be quite responsive in case you have questions.

Another option could be the *python* program MAGeCK

# Answer
> 1 votes
From my experience, there is few about CRISPR pooled screenings in R. One of them is ScreenBEAM. It uses a linear mixture model to determine significance at gene level by considering all the guides together at the same step in the analysis. However, it does not work very well for me.

There are other ways of analyze these screenings outside of R that are easy to use such the suggested MAGeCK or others such HiTSelect, BAGEL or CRISPR-Analyzer.

CRISPR-Analyzer is an online tool that could be useful and that recapitulates some of the mentioned methods in a single analysis. However, if you want to work from command line I would recommend you BAGEL for gene essentiality experiments and MAGeCK for other experimental setup (i.e, control vs treatment)

Sorry for not putting all the links but I can only put two at the moment!

---
Tags: r, crispr
---


# Question
Title: How can HISAT2/StringTie report decimal coverage values
I have performed RNA-seq analysis using HISAT2 & StringTie workflow suggested in: Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown.

Some of the transcripts/exons have decimal coverage values (eg., 1.1 or 2.59).

**My question is:** How can these tools report decimal coverage? If only half of the read overlaps exon will this exon have 0.5 coverage?

# Answer
> 6 votes
I suspect this is average coverage across the transcript or exon.

Formula for average coverage:

Where:

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

---
Tags: rna-seq
---


# Question
Title: What Ensembl genome version should I use for alignments? (e.g. toplevel.fa vs. primary_assembly.fa)
When you look at all the genome files available from Ensembl. You are presented with a bunch of options. Which one is the best to use/download?

You have a combination of choices.

First part options:

* **dna\_sm** \- Repeats soft-masked (converts repeat nucleotides to lowercase)
* **dna\_rm** \- Repeats masked (converts repeats to to N's)
* **dna** \- No masking

Second part options:

* **.toplevel** \- Includes haplotype information (not sure how aligners deal with this)
* **.primary\_assembly** \- Single reference base per position

Right now I usually use a non-masked primary assembly for analysis, so in the case of humans: *Homo\_sapiens.GRCh38.dna.primary\_assembly.fa.gz*

Does this make sense for standard RNA-Seq, ChIP-Seq, ATAC-Seq, CLIP-Seq, scRNA-Seq, etc... ?

In what cases would I prefer other genomes? Which tools/aligners take into account softmasked repeat regions?

# Answer
> 9 votes
Generally, you should use the soft-masked or unmasked primary assembly. Cross-species whole-genome aligners, especially older ones, do need to know soft-masked regions; otherwise they can be impractically slow for mammalian genomes. Modern read aligners are designed to work with repeats efficiently and therefore they don't need to see the soft mask.

For GRCh38, though, I would recommend to use the official build at GRC FTP. Most people will probably choose "no\_alt\_analysis\_set". Using the Ensembl version is discouraged due to its chromosome naming. We more often use "chr1" instead of "1" for GRCh38. At one point, Ensembl actually agreed to use "chr1" as well, but didn't make that happen due to technical issues, I guess.

As to alternate haplotypes, most aligners can't work with them; no variant callers can take the advantage of these sequences, either. When you align to a reference genome containing haplotypes with an aligner not supporting these extra sequences, you will get poor mapping results.

# Answer
> 13 votes
There's rarely a good reason to use a hard-masked genome (sometimes for blast, but that's it). For that reason, we use soft-masked genomes, which only have the benefit of showing roughly where repeats are (we never make use of this for our \*-seq experiments, but it's there in case we ever want to).

For primary vs. toplevel, very few aligners can properly handle additional haplotypes. If you happen to be using BWA, then the toplevel assembly would benefit you, but only if you use a dedicated wrapper to handle the ALT information, see bwakit. If you use BWA (bwa-mem) right from the command line without this wrapper then **do not use** the toplevel assembly. For STAR/hisat2/bowtie2/BBmap/etc. the haplotypes will just cause you problems due to increasing multimapper rates incorrectly. Note that none of these actually use soft-masking.

# Answer
> 4 votes
> Which tools/aligners take into account softmasked repeat regions?

If you're doing whole genome - whole genome alignment (rather than read alignment) then using the softmasked genome is definitely best. Tools suitable for such large scale alignments task tend to skip marked repeats completely in their initial steps to prevent the build up of bogus short alignments that can have a massive performance impact in terms of time and memory usage. For example, LASTZ skips lower case letters during the seeding stage.

# Answer
> 2 votes
TOPLEVEL

These files contains all sequence regions flagged as toplevel in an Ensembl schema. This includes chromsomes, regions not assembled into chromosomes and N padded haplotype/patch regions.

E.g: I used the soft masked assemblies for genome annotation pipelines like MAKER, also toplevel unmasked ones for RNA-seq, ChipSeq analysis

PRIMARY ASSEMBLY

Primary assembly contains all toplevel sequence regions excluding haplotypes and patches. This file is best used for performing sequence similarity searches where patch and haplotype sequences would confuse analysis.

---
Tags: fasta, genome, ensembl, biomart
---


# Question
Title: Difference between BWA-backtrack and BWA-MEM
Many of my colleagues recommend I use BWA-MEM instead of regular old BWA. The problem is I don't understand why and reading the BWA man page doesn't seem to help the matter.

What is the difference between BWA and BWA-MEM? And, in which instances would you employ one over the other?

# Answer
> 12 votes
**TL;DR:**

BWA-backtrack is based on backtracking. This approach is appropriate only when the dissimilarity between the reads and the reference is low, or when you want to find all best hits or enumerate all possible alignments up to a specified number of errors.

In all other situations, BWA-MEM is preferable as it can, thanks to its sophisticated strategy based on maximum exact matches, deal with errors better and also automatically switch between local and global alignment modes.

**Long description:**

I would like to provide some algorithmic insight since I believe that it can be very useful in this case. Both BWA-backtrack and BWA-MEM use the same indexing strategy (heavily relying on BWT-index), but the actual algorithms are quite different.

BWT-index (and also other full-text indexes such as suffix arrays or suffix trees) can easily find exact matches (imagine Ctrl+F-like search in a text editor), but any differences between the read and the reference, such as sequencing errors or genomic variants, make the situation complicated. One then needs to somehow transform inexact matching to exact matching, and all three BWA mappers (note that there exists also BWA-SW, but it is deprecated) use quite different strategies.

**BWA-backtrack** looks for substrings of the reference, which would be similar to the entire read (end-to-end) using an algorithm called backtracking. First, it searches occurences of the read without any "corrections". If nothing found, it consideres all possible single edits; then two edits, etc. To make mapping efficient, one usually wants to stop with the first found alignment as this would be the best one. When required, it is also possible to find the other equally good alignments or to enumerate all alignments up to some edit distance or withing some divergance rate (see the -N option of BWA-backtrack).

It turns out that the time required for finding an alignment can be exponential in the number of errors, which is probably the main problem of backtracking-based approaches. To prevent huge overheads due to dissimilar reads, one needs to limit the number of allowed errors to some reasonable number (see maxDiff in the BWA man page) and consider the other reads unaligned. In the case of BWA-backtrack, the minimum required identity level is ~97% with the default options (see the -n option).

In fact, the algorithm is more complicated and uses various heuristics such as seed-and-extend or Z-dropoff in order to make the computation fast enough (at the price of lower accuracy). If you are interested in more details, all these tricks are well described in the paper.

**BWA-MEM** uses quite a different strategy. It detects long exact matches between the read and the reference, and then chains them into local or global alignments, based on what is more appropriate in that specific case. Such an automatic local-global switching can be very powerful and BWA-MEM works well with various types of data (short reads, long reads, low error rates, high error rates, etc.).

# Answer
> 9 votes
To quote the Introduction to `BWA` on sourceforge:

> 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.

In short, for anything where you have read lengths over 70bp BWA-MEM is faster, and more accurate.

---
Tags: sequence-alignment, bwa, read-mapping, reads
---


# Question
Title: Publicly available genome sequence database for viruses?
As a small introductory project, I want to compare genome sequences of different strains of influenza virus.

What are the publicly available databases of influenza virus gene/genome sequences?

# Answer
> 17 votes
There area few different influenza virus database resources:

* The Influenza Research Database (IRD) (a.k.a FluDB - based upon URL)
    + A NIAID Bioinformatics Resource Center or BRC which highly curates the data brought in and integrates it with numerous other relevant data types
* The NCBI Influenza Virus Resource
    + A sub-project of the NCBI with data curated over and above the GenBank data that is part of the NCBI
* The GISAID EpiFlu Database
    + A database of sequences from the Global Initiative on Sharing All Influenza Data. Has unique data from many countries but requires user agree to a data sharing policy.
* The OpenFluDB
    + Former GISAID database that contains some sequence data that GenBank does not have.

For those who also may be interested in other virus databases, there are:

* Virus Pathogen Resource (VIPR)
    + A companion portal to the IRD, which hosts curated and integrated data for most other NIAID A-C virus pathogens including (but not limited to) Ebola, Zika, Dengue, Enterovirus, and Hepatitis C
* LANL HIV database
    + Los Alamos National Laboratory HIV database with HIV data and many useful tools for all virus bioinformatics
* PaVE: Papilloma virus genome database (from quintik comment)
    + NIAID developed and maintained Papilloma virus bioinformatics portal

Disclaimer: I used to work for the IRD / VIPR and currently work for NIAID.

# Answer
> 4 votes
In addition to what others have suggested, I would also recommend PaVE as a resource. This is a curated database maintained by the NIAID and current holds over 300 papilloma virus genomes.

# Answer
> 2 votes
Influenza virus resource at NCBI or FluDB.

---
Tags: database, public-databases, covid-19, genome-sequencing, sars-cov-2
---


title: samtools depth print out all positions

text:
# Question
Title: samtools depth print out all positions
I am trying to use `samtools depth` (v1.4) with the -a option and a bed file listing the human chromosomes chr1-chr22, chrX, chrY, and chrM to print out the coverage at every position:

```
cat GRCh38.karyo.bed | awk '{print $3}' | datamash sum 1
3088286401

```

I would like to know how to run `samtools depth` so that it produces 3,088,286,401 entries when run against a GRCh38 bam file:

```
samtools depth -b $bedfile -a $inputfile

```

I tried it for a few bam files that were aligned the same way, and I get differing number of entries:

```
3087003274
3087005666
3087007158
3087009435
3087009439
3087009621
3087009818
3087010065
3087010408
3087010477
3087010481
3087012115
3087013147
3087013186
3087013500
3087149616

```

Is there a special flag in `samtools depth` so that it reports all entries from the bed file?

If `samtools depth` is not the best tool for this, what would be the equivalent with `sambamba depth base`?

```
sambamba depth base --min-coverage=0 --regions $bedfile $inputfile

```

Any other options?

# Answer
You might try using bedtools genomecov instead. If you provide the -d option, it reports the coverage at every position in the BAM file.

```
bedtools genomecov -d -ibam $inputfile > "${inputfile}.genomecov"

```

You can also provide a BED file if you just want to calculate in the target region.
> 6 votes

# Answer
Just replace the `-a` option with `-aa`:

```
samtools depth -b $bedfile -aa $inputfile

```

I see that you're using the GRCh38 human reference genome build, which includes alternate scaffolds that represent a wider variety of genomic variation in the human genome. This was quite a substantial change from previous reference genome builds, and many bioinformatics tools still haven't been updated to accept the GRCh38 reference.

The discrepancy in entry size is probably because there are alternate contigs that are uncovered by the reads that have been mapped. This is expected, but can cause a bit of confusion.
> 2 votes

---
Tags: sam, samtools, sambamba
---