title: Mapping drug names to ATC codes

text:
# Question
Title: Mapping drug names to ATC codes
I'm interested working with the medication information provided by the UK Biobank. In order to get these into a usable form I would like to map them to ATC codes. Since many of the drugs listed in the data showcase include dosage information, doing an exact string match between drug names is not very effective. I've considered using something like fuzzywuzzy to do string matching between the medications in the data showcase and the ATC drug names but validating the matches could still be a laborious process. Does anyone know of a tool that can match drug names to ATC codes or some other drug ontology? If not, maybe there's a better way to do it that I haven't thought of.

# Answer
> 5 votes
The CART tool let's you upload a set of names and map them (optionally in a fuzzy way) to STITCH 4 identifiers, and then use those to map to ATC codes (using the chemicals sources download file). It's a bit indirect, and I'm not sure what CART will do with the dosage info you mention.

# Answer
> 5 votes
DrugBank seems to have a tool to map ATC codes to drug names and DrugBank IDs.

A quick look in the XSD schema on the release page suggests the complete database includes ATC codes for drugs, you could then do a fuzzy match of the BioBank names against all of DrugBank's synonyms, or match on some other data (e.g. canonicalised SMILES) if available.

The downside is that there may not be complete overlap between UK BioBank and DrugBank. Additionally, DrugBank is under licence for commercial use.

---
Tags: drugs, ontology
---


# Question
Title: Accuracy of the original human DNA datasets sequenced by Human Genome Project?
The Human Genome Project was the project of 'determining the sequence of nucleotide base pairs that make up human DNA, and of identifying and mapping all of the genes of the human genome'. It was declared complete in 2003, i.e. 99% of the euchromatic human genome completed with 99.99% accuracy.

Are the datasets provided by HGP still accurate, or as accurate as was claimed in 2003?

Given the technology in the past (such as using old techniques), or any other reason (newer research studies), is it possible that the datasets are not as accurate as originally expected?

# Answer
> 18 votes
The HGP developed the first "reference" human genome - a genome that other genomes could be compared to, and was actually a composite of multiple human genome sequences.

The standard human reference genome is actually continually updated with major and minor revisions, a bit like software. The latest major version is called GRCh38, was released in 2013, and has since had a number of minor updates.

> Are the datasets provided by HGP still accurate?

Yes, in a sense, but we certainly have better information now. One way to measure the quality of the assembly is that the initial release from the HGP had hundreds of thousands of gaps - sequences that could not be resolved (this often occurs because of repetitive sequences). The newest reference genome has less than 500 gaps.

# Answer
> 8 votes
While the quality of the reference human assembly keeps improving, there are still misassemblies in it. A common problem is recent segmental duplications are occasionally collapsed into one sequence in the reference. Another issue is that the centromeric sequences in the reference are computationally generated, which are probably different from real sequences. Issues like these often complicate data analyses.

You should also beware that each human has a different genome. A large region having one copy in the reference genome may have two copies in a specific sample. While this is not really the problem with the reference, such copy-number changes will have the same effect as reference errors and mess up your pipelines.

There is still room for improvement to the human reference genome. In some regions, the CHM1 and CHM13 PacBio assemblies are better than the current reference genome at the larger scale. Illumina population data can produce better consensus at the single base level. GRC is continuously releasing new patches to the latest assembly.

---
Tags: hgp, phylogenetics
---


# Question
Title: Compare alignment quality of multiple sequencing runs aligned against the same reference genome
I have run Oxford Nanopore Technologies' MinION sequencing on the same DNA sample using three flowcells, each aligned against the same reference genome (E.coli K12 MG1655) using both BWA MEM and GraphMap and stored as BAM files.

How can I quantitatively and efficiently analyse the quality of alignment (percentage identity, insertion rate, deletion rate) of each of these files?

# Answer
> 6 votes
Qualimap will do this for you.

1. Go to qualimap.bioinfo.cipf.es
2. Run qualimap (default params are fine) on each BAM file
3. Open up the HTML output, and you can read off the %identity (they measure the opposite, i.e. mismatch rate, but 100% - mismatch rate is %identity of course), indel rate, etc.

One thing to watch out for (you don't mention it in your question, but just in case) is that you cannot directly compare Q scores - these are a bit of a mess and calculated very differently in each piece of software.

Unsolicited suggestion: you might also try NGM-LR for mapping MinION data. We've found it beats the others for our data (though we map to a distant reference).

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


# Question
Title: What are the optimal parameters for docking a large ligand using Hex?
I'm looking to dock a large ligand (~90kDa) to a receptor slightly larger receptor (~125kDa) using Hex. If anyone is familiar with docking large structures, are there any recommended parameters for finding the best docking solution?

Parameters in particular:

* Number of Solutions
* N order of correlation for initial and final searches
* Receptor Range
* Ligand Range

# Answer
> 1 votes
While I have no experience with this *specific* question you have, according to Feinstein & Brylinski (2015), a "fully automated procedure", i.e. a Perl script, can be used to optimize the box size itself, and it can be found here.

---
Tags: proteins, docking
---


# Question
Title: What is the difference between a transcriptome and a genome?
I have a computer engineering background, not biology.

I started working on a bioinformatics project recently, which involves *de-novo* assembly. I came to know the terms `Transcriptome` and `Genome`, but I cannot identify the difference between these two.

I know a transcriptome is the set of all messenger RNA molecules in a cell, but am not sure how this is different from a genome.

# Answer
> 18 votes
In brief, the “genome” is the collection of all DNA present in the nucleus and the mitochondria of a somatic cell. The initial product of genome expression is the “transcriptome”, a collection of RNA molecules derived from those genes.

# Answer
> 5 votes
They are two very different things. Your genome is a large section of about 3 billion DNA nucleotide bases. It has no concept of exon and introns.

Transcriptome is a study of transcriptions. You have introns and exons. We can now talk about alternative splicing and gene expression.

You can think your genome is like a cooking recipe. While it's good to have a good recipe, you can't do much if you don't use it for cooking.

---
Tags: genome, transcriptome
---


# Question
Title: Alignment based vs reference-free (transcriptome analysis)?
I want to focus on transcriptome analysis. We know it's possible to analyze RNA-Seq experiment based on alignment or k-mers.

Possible alignment workflow:

Possible reference-free workflow:

* Quantify sequence reads with Kallisto reference-free index

Both strategy generate gene expression table.

**Q:** What are pros and cons for each of the approach? Can you give guideline?

# Answer
First of all, I would emphasize that "alignment-free" quantification tools like Salmon and Kallisto are *not* reference-free. The basic difference between them and more traditional aligners is that they do not report a specific position (either in a genome or transcriptome) to which a read maps. However, their overall purpose is still to quantify the expression levels (or differences) of a known set of transcripts; hence, they require a reference (which could be arbitrarily defined).

The most important criterion for deciding which approach to use (and this is true of almost everything in genomics) is exactly what question you would like to answer. If you are primarily interested in quantifying and comparing expression of mature mRNA from known transcripts, then a transcriptome-based alignment may be fastest and best. However, you may miss potentially interesting features outside of those known transcripts, such as new isoforms, non-coding RNAs, or information about pre-mRNA levels, which can often be gleaned from intronic reads (see the EISA method).

This paper also has some good considerations about which tools may work best depending on the question you want to answer.

Finally, another fast and flexible aligner (which can be used with or without a reference transcriptome) is STAR.
> 11 votes

# Answer
I wouldn't say Kallisto (or Salmon) are reference-free. They use a transcriptome as reference anda concept called *pseudo-alignment* which greatly speed up the process of assigning your reads to a transcript.

That said, both approaches of (i) mapping against a reference genome (what you called *alignment workflow* ) and (ii) mapping against a reference transcriptome will serve different purposes

Transcriptome mapping using pseudoalignent is becoming the method of choice for gene/transcript quantification and differential expression analysis. The drawback is that you only focus on known transcripts

Two typical workflows are:

* Kallisto followed by sleuth
* Salmon, followed by tximport and DESeq2/EdgeR

Genome mapping is useful for, per example discovery of new isoforms. You shouldn't use TopHat anymore as it has been discontinued by the author.

A typical workflow would be:

* Hisat2 (alignment)
* StringTie (transcript assembly and abundance estimation)
* Ballgown (differential expression)
> 7 votes

---
Tags: transcriptome, rna-seq
---


# Question
Title: How to run HYPHY on multiple files non-interactively?
I would like to test for positive selection in a large set of genes.

I want to have a yes/no answer to the question if gene was evolving under positive selection. Therefore I chose model BUSTED for my analysis, which is implemented in HYPHY.

When I launch HYPHY in the command-line, I get a series of question I have to answer in the interactive form.

Is there way to perform a batch analysis involving multiple datasets, i.e. many sequence alignments?

# Answer
**TL;DR**

I made the following bash script to generate an input file for HYPHY.

```
#!/bin/bash
cat << EOF
inputRedirect = {};
inputRedirect["01"]="Universal"; // genetic code
inputRedirect["02"]="$(readlink -f $1)"; // codon data
inputRedirect["03"]="$(readlink -f $2)"; // tree
inputRedirect["04"]="${3:-All}"; // Test for selection on a branch
inputRedirect["05"]=""; // complete selection

ExecuteAFile (HYPHY_LIB_DIRECTORY+"TemplateBatchFiles/BUSTED.bf", inputRedirect);
EOF

```

If you want to test all the branches jointly run:

```
./gen_busted.sh alignment.phy tree.nwk > script.bf

```

Or for a particular branch:

```
./gen_busted.sh alignment.phy tree.nwk branchLabel > script.bf

```

And then run:

```
HYPHYMP script.bf

```

Now you can use the script to generate a `.bf` file for every sequence alignment you have, and run HYPHY for every file.

**Longer version**

I will describe how to make a similar script for any model, such as MEME, RELAX or PARRIS.

First, find the `.bf` script which performs the analysis. On GNU/Linux it should be at `HYPHY_INSTALLATION_PATH/lib/hyphy/TemplateBatchFiles` (it might be different between different OS). In case of BUSTED your script has a name `BUSTED.bf`.

Now run this model in the interactive mode in order to record all the inputs required by the model. Make sure to specify **full paths** to all the files.

```
HYPHYMP HYPHY_INSTALLATION_PATH/lib/hyphy/TemplateBatchFiles/BUSTED.bf

```

In case of BUSTED the inputs are:

1. Genetic code.
2. Codon data (sequence alignment in phylip format).
3. Tree (in newick format).
4. Which branches to test.

Now for every input alignment you need to generate a batch file which looks like this:

```
inputRedirect = {};
inputRedirect["01"]="Universal"; // genetic code
inputRedirect["02"]="/path/to/alignment.phy"; // codon data
inputRedirect["03"]="/path/to/tree.nwk"; // tree
inputRedirect["04"]="All"; // Test for selection on all branches
inputRedirect["05"]="BRANCH1"; // Test for selection on branch1
inputRedirect["06"]="BRANCH2"; // Test for selection on branch2
inputRedirect["07"]=""; // complete selection

ExecuteAFile (HYPHY_LIB_DIRECTORY+"TemplateBatchFiles/BUSTED.bf", inputRedirect);

```

Now you can create a script which generates the file. The most important thing, don't forget to use full paths in the `.bf` file.
> 5 votes

---
Tags: codon, positive-selection
---


title: Single-sample vs. joint genotyping

text:
# Question
Title: Single-sample vs. joint genotyping
I am trying to understand the benefits of joint genotyping and would be grateful if someone could provide an argument (ideally mathematically) that would clearly demonstrate the benefit of joint vs. single-sample genotyping.

This is what I've gathered from other resources (Biostars, GATK forums, etc.)

* Joint-genotyping helps control FDR because errors from individually genotyped samples are added up, and amplified when merging call-sets (by Heng Li on https://www.biostars.org/p/10926/)

If someone understands this, can you please clarify what is the difference on the overall FDR rate between the two scenarios (again, with an example ideally)

I don't understand how the presence of a confidently called variant at the same locus in another individual can affect the genotyping of an individual with low coverage. Is there some valid argument that allows one to consider reads from another person as evidence of a particular variant in a third person? What are the assumptions for such an argument? What if that person is from a different population with entirely different allele frequencies for that variant?

Having read several of the papers (or method descriptions) that describe the latest haplotype-aware SNP calling methods (HaplotypeCaller, freebayes, Platypus) the overall framework seems to be:

* 1. Establish a prior on the allele frequency distribution at a site of interest using one (or combination) of: non-informative prior, population genetics model-based prior like Wright Fisher, prior based on established variation patterns like dbSNP, ExAC, or gnomAD.
* 2. Build a list of plausible haplotypes in a region around the locus of interest using local assembly.
* 3. Select haplotype with highest likelihood based on prior and reads data and infer the locus genotype accordingly.

At which point(s) in the above procedure can information between samples be shared or pooled? Should one not trust the AFS from a large-scale resource like gnomAD much more than the distribution obtained from other samples that are nominally party of the same "cohort" but may have little to do with each other because of different ancestry, for example?

I really want to understand the justifications and benefits offered by multi-sample genotyping and would appreciate your insights.

# Answer
> 10 votes
Say you are sequencing to 2X coverage. Suppose at a site, sample S has one reference base and one alternate base. It is hard to tell if this is a sequencing error or a heterozygote. Now suppose you have 1000 other samples, all at 2X read depth. One of them has two ALT bases; 10 of them have one REF and one ALT. It is usually improbable that all these samples have the same sequencing error. Then you can assert sample S has a het. Multi-sample calling helps to increase the sensitivity of not so rare SNPs. Note that what matters here is the assumption of error independency. Ancestry only has a tiny indirect effect.

Multi-sample calling penalizes very rare SNPs, in particular singletons. When you care about variants only, this is for good. Naively combining single-sample calls yields a higher error rate. Multi-sample calling also helps variant filtering at a later stage. For example, for a sample sequenced to 30X coverage, you would not know if a site at 45X depth is caused by a potential CNV/mismapping or by statistical fluctuation. When you see 1000 30X samples at 45X depth, you can easily know you are looking at a CNV/systematic mismapping. Multiple samples enhance most statistical signals.

Older methods pool all BAMs when calling variants. This is necessary because a single low-coverage sample does not have enough data to recover hidden INDELs. However, this strategy is not that easy to massively parallelized; adding a new sample triggers re-calling, which is very expensive as well. As we are mostly doing high-coverage sequencing these days, the old problem with INDEL calling does not matter now. GATK has this new single-sample calling pipeline where you combine per-sample gVCFs at a later stage. Such sample combining strategy is perhaps the only sensible solution when you are dealing with 100k samples.

The so-called haplotype based variant calling is a separate question. This type of approach helps to call INDELs, but is not of much relevance to multi-sample calling. Also, of the three variant callers in your question, only GATK (and Scalpel which you have not mentioned) use assembly at large. Freebayes does not. Platypus does but only to a limited extent and does not work well in practice.

I guess what you really want to talk about is imputation based calling. This approach further improves sensitivity with LD. With enough samples, you can measure the LD between two positions. Suppose at position 1000, you see one REF read and no ALT reads; at position 1500, you see one REF read and two ALT reads. You would not call any SNPs at position 1000 even given multiple samples. However, when you know the two positions are strongly linked and the dominant haplotypes are REF-REF and ALT-ALT, you know the sample under investigation is likely to have a missing ALT allele. LD transfers signals across sites and enhances the power to make correct genotyping calls. Nonetheless, as we are mostly doing high-coverage sequencing nowadays, imputation based methods only have a minor effect and are rarely applied.

# Answer
> 1 votes
The benefit to additional samples is seen in your point 1. The likelihood of making a variant call is a function of (1) the depth of coverage supporting a given variant (ignoring mapping/base quality considerations) and (2) the likelihood of that variant existing given background knowledge. With low depth and no background knowledge, poorly covered variants will be assumed to be sequencing errors. Adding more samples can just serve then to increase the background knowledge on a position.

---
Tags: genotyping, gatk
---


# Question
Title: How to calculate overlapping genes between two genome annotation versions
I have two annotations of the same genome generated with different annotation pipelines. I want to identify overlapping gene models.

*An important feature of this genome is that there are many 'genes within genes', i.e. a genemodel in the intron of another genemodel*. Therefore, I only want to count two genemodels as overlapping when their coding sequence exon annotations overlap.

Using something like `bedtools intersect` it is straightforward to calculate overlap between the gene-level annotations.

However: I am not sure how to select genes as overlapping when only their coding sequence exons (CDS features) overlap.

# Answer
> 6 votes
**Short Answer:** In my opinion, my approach would be to pull out the CDS exons and run bedtools on those.

**A Few More Details:** When you pull out the exons, make sure that you assign them all IDs if the don't already have them assigned and record which IDs "belong" to which genes. Now when you get exons that overlap, you know that they are coding and you can tie them back to which genes they originate from.

---
Tags: sequence-annotation
---


title: How to compute RPKM in R?

text:
# Question
Title: How to compute RPKM in R?
I have the following data of fragment counts for each gene in 16 samples:

```
> str(expression)
'data.frame':   42412 obs. of  16 variables:
 $ sample1 : int  4555 49 122 351 53 27 1 0 0 2513 ...
 $ sample2 : int  2991 51 55 94 49 10 55 0 0 978 ...
 $ sample3 : int  3762 28 136 321 94 12 15 0 0 2181 ...
 $ sample4 : int  4845 43 193 361 81 48 9 0 0 2883 ...
 $ sample5 : int  2920 24 104 151 50 20 32 0 0 1743 ...
 $ sample6 : int  4157 11 135 324 58 26 4 0 0 2364 ...
 $ sample7 : int  3000 19 155 242 57 12 18 2 0 1946 ...
 $ sample8 : int  5644 30 227 504 91 37 11 0 0 2988 ...
 $ sample9 : int  2808 65 247 93 272 38 1108 1 0 1430 ...
 $ sample10: int  2458 37 163 64 150 29 729 2 1 1049 ...
 $ sample11: int  2064 30 123 51 142 23 637 0 0 1169 ...
 $ sample12: int  1945 63 209 40 171 41 688 3 2 749 ...
 $ sample13: int  2015 57 432 82 104 47 948 4 0 1171 ...
 $ sample14: int  2550 54 177 59 201 36 730 0 0 1474 ...
 $ sample15: int  2425 90 279 73 358 34 1052 3 3 1027 ...
 $ sample16: int  2343 56 365 67 161 43 877 3 1 1333 ...

```

How do I compute RPKM values from these?

# Answer
> 31 votes
First off,

**Don’t use RPKMs**.

They are truly deprecated because they’re confusing once it comes to paired-end reads. If anything, use *FPKM*s, which are mathematically the same but use a more correct name (do we count paired reads separately? No, we count *fragments*).

Even better, use TPM (= transcripts per million), or an appropriate cross-library normalisation method. TMP is defined as:

$$ \text{TPM}\_\color{orchid}i = {\color{dodgerblue}{\frac{x\_\color{orchid}i}{{l\_\text{eff}}\_\color{orchid}i}}} \cdot \frac{1}{\sum\_\color{tomato}j \color{dodgerblue}{\frac{x\_\color{tomato}j}{{l\_\text{eff}}\_\color{tomato}j}}} \cdot \color{darkcyan}{10^6} $$

where

* $\color{orchid}i$: transcript index,
* $x\_i$: transcript raw count,
* $\color{tomato}j$ iterates over all (known) transcripts,
* $\color{dodgerblue}{\frac{x\_k}{{l\_\text{eff}}\_k}}$: rate of fragment coverage per nucleobase ($l\_\text{eff}$ being the effective length),
* $\color{darkcyan}{10^6}$: scaling factor (= “per millions”).

That said, FPKM an be calculated in R as follows. Note that most of the calculation happens in log transformed number space, to avoid numerical instability:

```
fpkm = function (counts, effective_lengths) {
    exp(log(counts) - log(effective_lengths) - log(sum(counts)) + log(1E9))
}

```

Here, the *effective length* is the transcript length minus the mean fragment length plus 1; that is, all the possible positions of an average fragment inside the transcript, which equals the number of all distinct fragments that can be sampled from a transcript.

This function handles *one* library at a time. I (and others) argue that this is the way functions should be written. If you want to apply the code to multiple libraries, nothing is easier using ‹dplyr›:

```
tidy_expression = tidy_expression %>%
    group_by(Sample) %>%
    mutate(FPKM = fpkm(Count, col_data$Lengths))

```

However, the data in the question isn’t in tidy data format, so we first need to transform it accordingly using ‹tidyr›:

```
tidy_expression = expression %>%
    gather(Sample, Count)

```

---

This equation fails if all your counts are zero; instead of zeros you will get a vector of NaNs. You might want to account for that.

---

And I mentioned that TPMs are superior, so here’s their function as well:

```
tpm = function (counts, effective_lengths) {
    rate = log(counts) - log(effective_lengths)
    exp(rate - log(sum(exp(rate))) + log(1E6))
}

```

# Answer
> 8 votes
RPKM is defined as:

> RPKM = numberOfReads / ( geneLength/1000 * totalNumReads/1,000,000 )

As you can see, you need to have gene lengths for every gene.

Let's say `geneLength` is a vector which have the same number of rows as your `data.frame`, and every value of the vector corresponds to a gene (row) in `expression`.

```
expression.rpkm <- data.frame(sapply(expression, function(column) 10^9 * column / geneLength / sum(column)))

```

**Regarding numerical stability**

It is suggested in one of the answers, that all the computations should be performed in a log-transformed scale. In my opinion there is absolutely no reason for doing that. IEEE binary64 stores a number as binary number 1.b\_{51}b{50}...b\_0 times 2^{e-1023}. The *relative* precision doesn't depend on the exponent value given that a number is in the range \[~10^-308; 10^308\].

In case of RPKM we can only get out of the range if total number of reads is around 10^300, which is not realistic at all.

On the bright site there is not much harm in doing computations in the log-scale either. Worst case you loose a bit of precision.

# Answer
> 1 votes
If you are planning to do a differential expression analysis, you will probably don't need the RPKM calculation.

RPK= No.of Mapped reads/ length of transcript in kb (transcript length/1000)

RPKM = RPK/total no.of reads in million (total no of reads/ 1000000)

**The whole formula together:**

RPKM = (10^9 * C)/(N * L) Where,

C = Number of reads mapped to a gene

N = Total mapped reads in the experiment

L = exon length in base-pairs for a gene

# Answer
> 0 votes
If you are looking for a more visual solution (in addition to the other answers), NCI Genomic Data Commons (TCGA repository) offers a nice formula:

> Where:
>
> RCg: Number of reads mapped to the gene
>
> RCpc: Number of reads mapped to all protein-coding genes
>
> RCg75: The 75th percentile read count value for genes in the sample
>
> L: Length of the gene in base pairs

As others have pointed out, FPKMs have some problems. GDC also calculates FPKM-UQ values that are upper quartile normalized. Those are recommended for cross-sample comparison and differential expression analysis.

---
Tags: rna-seq, r, transcriptome, fpkm
---


title: Linear models of complex diseases

text:
# Question
Title: Linear models of complex diseases
A popular framework to analyze differences between groups, either experiments or diseases, in transcriptomics is using linear models (limma is a popular choice).

For instance we have a disease D with three stages as defined by clinicians, A, B and C. 10 samples each stage and the healthy H to compare with is RNA-sequenced. A typical linear model would be to observe the three stages`~A+B+C` independently. The data of each stage is not from the same person. (but for the question assume it isn't)

My understanding is that such a model would not take into account that stage C appears only on 30% of patients in stage B. And that a healthy patient upon external factors can jump to stage B.

If we want to find the role of a gene in the disease we should include somehow this information in the model. Which makes me think about mixing linear models and hidden Markov chains.

How can such a disease be described in terms of linear models with such data and information?

# Answer
> 4 votes
There are two potential sources of bias in this design.

1. We cannot distinguish correlation from causation.

Imagine two cases. In the first, the disease progression is inducing immune response. Later stages will be associated with the higher gene expression levels. In the second scenario, the disease is caused by overexpression of a gene. Later stages will be also associated with the higher expression.

This is typical for observational studies. But I just want to mention that a special care should be taken during interpretation of the results.

2. If we are not following our individuals, we cannot distinguish correlation and avoidance.

Let's say the disease is lethal in certain cases; the survival is negatively correlated with the disease stage. Now imagine there is a gene, which causes severe symptoms when highly expressed. On the later stages you will only observe those patients in which the gene was not highly expressed. From that you would conclude that gene expression is decreasing with the disease progression. In reality this gene is very important and causal, you just do not have patients which are alive to see this.

This is similar to Wald's studies of aircrafts (Survival bias).

> Researchers from the Center for Naval Analyses had conducted a study of the damage done to aircraft that had returned from missions, and had recommended that armor be added to the areas that showed the most damage. Wald proposed that the Navy instead reinforce the areas where the returning aircraft were unscathed, since those were the areas that, if hit, would cause the plane to be lost.

I think the second point is crucial, and can and will lead to false conclusions. I do not see immediately how to avoid this bias if you study different individuals.

Therefore I suggest that the same individuals are followed for a long time.

There are different approaches you can use later. For example you can have two-stage procedure:

1. Identify genes which are differentially expressed between healthy (H) and sick (A, B or C).
2. Build a linear model of disease stage stage ~ gene1 + gene2 + ..., using genes identified at step 1.
3. Similarly build a linear model of survival as a function gene expression.
4. It is possible you can use logistic regression to infer probability of transition from stage B to stage C.

This is not the only possible approach of course. As you suggested, Markov model is also applicable.

While following the same individuals it is possible to use continuous-time Markov model for such a study. In this case discrete or continuous measurements are recorded at certain time points, and the model parameters are inferred using maximum likelihood.

I'm not an expert in this field, but I think this paper describing the msm package for R will be useful. It also supports hidden Markov models.

---
Tags: transcriptome
---


# Question
Title: Understanding the significance of BLAT score calculations for read fragment alignment
I'm attempting to reconstruct the read fragments that form a hit to a target, excluding/skipping gaps. *BLAST*'s `outfmt 6`, and alternatives such as Diamond's `M8` format only give the number of gap opens on the alignment. One can't just use the `start` and `end` positions of the alignment as indices to slice the read DNA sequence, as this will include the nucleotides that were gapped by the alignment.

`PSL` appeared to be the only format I could find that gave the actual start positions (and lengths) of each of the pieces that formed a hit (that is, each gap forces a new block).

The regular *BLAST* output is human readable and is a nightmare to parse. The `XML` format mentioned has a structure for each HSP, but an HSP can still include gaps:

```
  [...]
  <Hit_hsps>
    <Hsp>
      <Hsp_num>1</Hsp_num>
      <Hsp_bit-score>51.6</Hsp_bit-score>
      <Hsp_score>122</Hsp_score>
      <Hsp_evalue>9.8e-06</Hsp_evalue>
      <Hsp_query-from>2</Hsp_query-from>
      <Hsp_query-to>98</Hsp_query-to>
      <Hsp_hit-from>35</Hsp_hit-from>
      <Hsp_hit-to>65</Hsp_hit-to>
      <Hsp_query-frame>2</Hsp_query-frame>
      <Hsp_hit-frame>0</Hsp_hit-frame>
      <Hsp_identity>24</Hsp_identity>
      <Hsp_positive>27</Hsp_positive>
      <Hsp_gaps>1</Hsp_gaps>
      <Hsp_align-len>32</Hsp_align-len>
         <Hsp_qseq>ITAIGAGLQGPAGCEVIDAGGLLVMPGGIDTH</Hsp_qseq>
         <Hsp_hseq>IAAVGTGLE-PAGAEIIDAGGLLVMPGGIDVH</Hsp_hseq>
      <Hsp_midline>I A+G GL+ PAG E+IDAGGLLVMPGGID H</Hsp_midline>
    </Hsp>
  </Hit_hsps>
</Hit>

```

I'm trying to be lazy, I don't want to have to parse the `Hsp_qseq` and `Hsp_hseq` to determine exactly where my gaps are. I just want co-ordinates of where my read hits a protein. The `PSL` format (below) tells me this (`tStarts`):

```
match   mis-    rep.    N's     Q gap   Q gap   T gap   T gap   strand  Q               Q       Q       Q       T               T       T       T       block   blockSizes      qStarts  tStarts
        match   match           count   bases   count   bases           name            size    start   end     name            size    start   end     count
---------------------------------------------------------------------------------------------------------------------------------------------------------------
[...]
21      4       0       0       1       8       1       24      ++      <PROTEIN>   662     321     354     <READ>      101     2       101     2       8,17,   321,337,        2,50,

```

As an alternative, I've executed the `BLAT` aligner as a means to find alignments of translated DNA reads against a small protein database. I selected `BLAT` as it provided the easiest means to access information on the actual alignment blocks (with the smallest amount of parsing). That is, the PSL has a comma delimited field that enumerates each of the starting positions of a query to a target.

However, `BLAT` outputs a significantly greater number of alignments than alternative tools. This is likely due to its identity percentage threshold defaulting to 25% when emulating a `blastx` style query. I can cope with this (apologies to my disks for the I/O), but for filtering purposes, the output `PSL` format appears to leave me nothing but raw match counts to work with.

The `PSL` gives me what I want, except now I don't have information on its significance. I know that the BLAT FAQ describes a Perl script for adding "web-based identity and score calculations". I can actually get identity percentages and scores from this script, but I cannot tell what the score actually represents. Is it comparable to bit-score? How can I recover something like an e-value from a `PSL` output file?

# Answer
It would appear that `BLAT` and `exonerate` (as suggested by @terdon) are the only pieces of software that provide enough information for "HSPFragments" to be parsed from their outputs, at least according to BioPython's documentation, anyway:

> Most search programs only have HSPs with one HSPFragment in them, making these two objects inseparable. However, there are programs (e.g. BLAT and Exonerate) which may have more than one HSPFragment objects in any given HSP. If you are not using these programs, you can safely consider HSP and HSPFragment as a single union.

It would seem the options are sticking with `BLAT` and generating scores afterwards with the script described in their FAQ, or using `exonerate` which provides a raw score as part of its default output. What these scores actually mean I will have to investigate, but for now I think the case is closed!
> 2 votes

# Answer
I would suggest using blastx and parsing the output XML. As an example, have a look at this: http://www.perlmonks.org/?node\_id=1006367.

The relevant part of the hit would be the following:

```
<Hit>
  <Hit_num>1</Hit_num>
  <Hit_id>gnl|BL_ORD_ID|1515029</Hit_id>
  <Hit_def>43989.cce_0262 (Cyanothece ATCC 51142)</Hit_def>
  <Hit_accession>1515029</Hit_accession>
  <Hit_len>65</Hit_len>
  <Hit_hsps>
    <Hsp>
      <Hsp_num>1</Hsp_num>
      <Hsp_bit-score>40.0466</Hsp_bit-score>
      <Hsp_score>92</Hsp_score>
      <Hsp_evalue>0.00664016</Hsp_evalue>
      <Hsp_query-from>155</Hsp_query-from>
      <Hsp_query-to>253</Hsp_query-to>
      <Hsp_hit-from>12</Hsp_hit-from>
      <Hsp_hit-to>44</Hsp_hit-to>
      <Hsp_query-frame>-1</Hsp_query-frame>
      <Hsp_hit-frame>0</Hsp_hit-frame>
      <Hsp_identity>17</Hsp_identity>
      <Hsp_positive>27</Hsp_positive>
      <Hsp_gaps>0</Hsp_gaps>
      <Hsp_align-len>33</Hsp_align-len>
      <Hsp_qseq>LRGAICSMEHIEEALGKLKDWARKLIELLLGPR</Hsp_qseq>
      <Hsp_hseq>ITGAVCLMDYLEKVLEKLRELAQKLIETLLGPQ</Hsp_hseq>
      <Hsp_midline>+ GA+C M+++E+ L KL++ A+KLIE LLGP+</Hsp_midline>
    </Hsp>
  </Hit_hsps>
</Hit>

```

You should be able to get the position information from the following tags: `<Hsp_query-from>`,`<Hsp_query-to>`,`<Hsp_hit-from>`,`<Hsp_hit-to>`.

To translate back to you original DNA sequence I can imagine you need the `<Hsp_query-frame>` tag.

And then of course you get the evalue for free: `<Hsp_evalue>`
> 1 votes

---
Tags: sequence-alignment, blat, blast
---


# Question
Title: probeset to probeset mappings between Affymetrix arrays
I am interested in identifying mappings between different types of Affymetrix arrays. I am aware that mappings between gene and probeset can be extracted using Ensembl's Biomart database.

```
Ensembe gene id ENSG00000181019 maps to
1. AFFY HG-U133_Plus_2's 210519_s_at
2. AFFY HuGene-1_0-st-v1's 8002303

```

Is there any way to extract mappings of probeset ids between two arrays(eg: HG-U133\_Plus\_2, HuGene-1\_0-st-v1)?

```
eg: 210519_s_at and 8002303

```

# Answer
> 9 votes
If your question is: can probeset IDs from different platforms be mapped to one another in a similar way as mapping probesets to genes, then the answer is: Yes. BioMart allows you to map almost anything that has an ID to anything else that has an ID.

You can use BioMart either via the web interface or programatically. A brief guide to using the web interface, for mapping HG U133 Plus 2 to HuGene 1.0:

1. Go to the start page
2. Select H. sapiens Ensembl genes for your database; Human genes (GRCh38.p10) for your dataset
3. Click Filters in the left-hand column
4. Expand "REGION", scroll down to GENE and select "Input microarray probes/probesets ID list \[Max 500 advised\]"
5. Select AFFY HG U133 PLUS 2 probe ID(s) and either copy/paste or upload a list, one per line
6. Click Attributes in the left-hand column
7. Scroll through, uncheck what you don't want to see and choose what you do, for example Gene Stable ID, AFFY HuGene 1 0 st v1 probe and AFFY HG U133 Plus 2 probe
8. Finally click "Results" in the menu at the top of the left-hand column

And you should see a result like this (you'll need to click the "Results" button).

You should not expect that there be a one-to-one mapping, for two reasons:

* Genes have multiple transcripts and probesets are mapped to each transcript
* Some transcripts have more than one probeset: in this case, the HuGene IDs 8002301 and 8002303 map to transcripts for this gene

Finally: here's a programmatic example using R/BioMart:

```
library(biomaRt)
mart.hs <- useMart("ensembl", "hsapiens_gene_ensembl")
results <- getBM(attributes = c("ensembl_gene_id", "affy_hugene_1_0_st_v1", "affy_hg_u133_plus_2"),
                 filters = "affy_hg_u133_plus_2",
                 values = c("210519_s_at"),
                 mart = mart.hs)

results
  ensembl_gene_id affy_hugene_1_0_st_v1 affy_hg_u133_plus_2
1 ENSG00000181019               8002301         210519_s_at
2 ENSG00000181019               8002303         210519_s_at

```

# Answer
> 7 votes
Instead of biomaRt, it is also possible to use the mapping databases built into Bioconductor itself, and map from probe to gene, and then from gene to probe in the second. Some R code to convert between hgu133 and hgu95 using the same probe ID provided in another:

```
library(hgu133plus2.db)
library(hgu95av2.db)

query_probe <- "210519_s_at"

hgu133_ensembl <- select(hgu133plus2.db, keys = query_probe, columns = "ENSEMBL")

ensembl_hgu95 <- select(hgu95av2.db, keys = hgu133_ensembl$ENSEMBL, keytype = "ENSEMBL", columns = "PROBEID")

dplyr::inner_join(hgu133_ensembl, ensembl_hgu95, by = "ENSEMBL", suffix = c(".133", ".95"))

```

# Answer
> 2 votes
Other answers explain why there might not be one to one mapping between the probes.

The AbsID database does conversion based on mapping the probe sequences to a genome build, and then determines mappings based on overlapping genome alignment coordinates. This is really useful if you want to be sure that two probes are *actually likely* measuring the same transcript.

It is dependent on both probes aligning to the genome, however.

Disclaimer: I worked in the group that provides the AbsID mapping tool.

---
Tags: ensembl, biomart, microarray
---


title: How to convert FASTA to BED

text:
# Question
Title: How to convert FASTA to BED
I have a FASTA file:

```
> Sequence_1
GCAATGCAAGGAAGTGATGGCGGAAATAGCGTTAGATGTATGTGTAGCGGTCCC...
> Sequence_2
GCAATGCAAGGAAGTGATGGCGGAAATAGCGTTAGATGTATGTGTAGCGGTCCC....
....

```

I want to generate a BED file for each sequence like:

```
Sequence_1 0 1500
Sequence_2 0 1700

```

The BED regions will simply be the size of the sequences.

**Q:** I did that before with a one-line command. I don't remember what that is, it was on Biostars. I can't find the post now. What's the simplest way to do the conversion?

# Answer
We have many excellent answers! This will be an excellent reference for future users.

I found what exactly what I was asking in my question:

> https://www.biostars.org/p/191052/

```
$ pip install pyfaidx
$ faidx --transform bed test.fasta > test.bed

```

This is the one-line command I was asking. The other answers also work, but I want to accept my own answer.
> 7 votes

# Answer
You can do this easily with bioawk, which is a version of awk with added features facilitating bioinformatics:

```
bioawk -c fastx '{print $name"\t0\t"length($seq)}' test.fa

```

`-c fastx` tells the program that the data should be parsed as fasta or fastq format. This makes the `$name` and `$seq` variables available in the awk commands.
> 15 votes

# Answer
It's good practice to have your `FASTA` indexed, so you can leverage the `.fai` you are likely to already have. If not, you can just generate the index with `samtools` and use some `awk` to make your `BED`:

```
samtools faidx $fasta
awk 'BEGIN {FS="\t"}; {print $1 FS "0" FS $2}' $fasta.fai > $fasta.bed

```

This will maintain tab separation but you can drop the `BEGIN` statement to use spaces. The BED spec only requires "whitespace" for the simple `BED` format.
> 13 votes

# Answer
You could adapt this awk one-liner. Note that it assumes that sequence IDs are not longer than 100 characters and that there is no description following the sequence ID on the header line.

```
cat myseqs.fasta | awk '$0 ~ ">" {print c; c=0;printf substr($0,2,100) "\t0\t"; } $0 !~ ">" {c+=length($0);} END { print c; }'

```

Otherwise, any Bio\* library (Perl, Python, Ruby) provides FASTA format parsers which will extract sequence IDs and lengths.

I'd point out that whilst this resembles BED it is not, strictly-speaking, since BED maps to coordinates on a chromosome or some longer sequence object.
> 6 votes

# Answer
Inspired by this answer to a related question on read length distributions, you could do this with Biopython:

```
from Bio.SeqIO import parse
with open("regions.bed", "w") as bed:
    for record in parse("regions.fasta", "fasta"):
        print(record.id, 0, len(record.seq), sep="\t", file=bed)

```
> 4 votes

# Answer
Here is an approach with `BioPython`. The `with` statement ensures both the input and output file handles are closed and a *lazy* approach is taken so that only a single fasta record is held in memory at a time, rather than reading the whole file into memory, which is a bad idea for large input files. The solution makes no assumptions about the sequence ID lengths or the number of lines that the sequences are spread across:

```
from Bio import SeqIO

with open('sequences.fasta') as in_f, open('sequences.bed','w') as out_f:
    for record in SeqIO.parse(in_f, 'fasta'):
        out_f.write('{}\t0\t{}\n'.format(record.id, len(record)))

```
> 4 votes

# Answer
If the FASTA sequences are all on a single line, then the following perl one-liner should work:

```
cat myseqs.fasta | perl -ne 'if(/^>([^ ]+)/){print $1} else {print " 0 ",length,"\n"}'

```

Explanation:

* If the line begins with a '\>', then print everything up to the first space (but don't put a line break at the end)
* Otherwise, print " 0 ", followed by the length of the line, followed by a line break
> 0 votes

---
Tags: fasta, bed, file-formats, format-conversion
---


# Question
Title: Merge hundreds of small BAM files into a single BAM file
I am working with over a million (long) reads, and aligning them to a large genome. I am considering running my alignment jobs in parallel, distributing horizontally across hundreds of nodes rather than trying to run a single job with dozens of cores.

I would like to merge the sorted BAM files together for further downstream analysis. What is the most efficient way to do so while maintaining a valid file header and taking advantage of the fact that the input bam files are already sorted?

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

# Answer
> 8 votes
Merging sorted files is a linear operation, so any well-implemented tools that do it will do it with approximately the same efficiency. So `samtools merge` (use the most up-to-date version, as there have been improvements in merge header handling in the 1.3.x and 1.4.x versions), `picard MergeSamFiles`, etc.

These tools need to hold all the input BAM files open simultaneously, so depending on how many hundred of input files you have you may run into the system limit on open file descriptors. Use `ulimit` to maximise this first; if there are still too many, you may need to merge the first 500 files, then merge the next 500 into that, etc. Samtools does not do this internally; I'm not sure whether any of the other merge implementations do.

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


# Question
Title: How to filter out cross alignments from a BED file?
I have a BAM file:

```
@SQ SN:chr1 LN:248956422
@SQ SN:chrx LN:248956423
ST-E00110:348:HGVKKALXX:1:1201:5822:48670   323 chr1    9999    0   67H66M16H   chrx    1000    0   GATAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAAC  JJJJJJJJJJJJJJJJAJJJJJJJJJJJJFJJJJJJFJFJJJJJJFJJJJJJJJJJJA77FJFJJJ  NM:i:0  MD:Z:66 AS:i:66 XS:i:65 SA:Z:chr5,18606834,-,73S76M,34,0;   RG:Z:g1

```

There is a read aligned to `chr1`, and it's mate aligned to `chrx`.

I have a BED file:

```
chr1    0   100000  TestOnly

```

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

When I do:

> samtools view -L test.bed test.bam

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

My solution:

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

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

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

**Q:** Is there a better solution?

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

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

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

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

So, you want to remove those lines whose 7th field isn't `=` and, just in case, those lines whose 7th field isn't `=` *and* isn't the same as the 3rd field. You can therefore use something like this:

```
samtools view -L test.bed test.bam | awk '$7=="=" || $3==$7

```

And, to save as a bam file again:

```
samtools view -L test.bed test.bam | awk '$7=="=" && $3==$7 |
    samtolls view -b > fixed.bam

```

---

On a separate note, there's very rarely a need to chain multiple grep commands like that. You can just use `\|` (or `|` with the `-E` or `-P` options) to separate them. Something like:

```
samtools view -L test.bed test.bam | grep -v 'chrx\|chr2\|chr10\|chrN'

```

Or

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

```
> 6 votes

---
Tags: sam, bed, genomics
---


# Question
Title: Are soft-clipped bases used for variant calling in samtools + bcftools?
If there are soft clipped base pairs specified in the CIGAR string for a read in a SAM/BAM file, will these be used for variant calling in a samtools \+ bcftools workflow?

The GATK HaplotypeCaller, for example, has an explicit option `--dontUseSoftClippedBases` for whether to use soft clipped bases. The samtools documentation does not mention clipped bases.

# Answer
No, samtools (and therefore bcftools) does not use soft-clipped bases. You can quickly confirm this by using either `samtools depth` or `samtools mpileup` to look at a region with a soft-clipped alignment. You'll note that the soft-clipped region isn't used in the depth/pileup (both tools use the same underlying code, so it doesn't matter which you use). If you're curious, samtools ignores soft-clipped bases because it's based on making a per-base stack of alignments covering each position. In the BAM format, alignments are sorted and assigned to bins according to their start/end positions, which won't include soft-clipping. Consequently, when samtools is making the pileup it won't even see the alignments that would overlap a given base if soft-clipped bases were included.

This then sort of begs the question of what GATK's HaplotypeCaller is doing differently. There, regions in the genome are essentially assembled in a small de Bruijn graph, which allows for soft-clipped bases around indels to then be resolved, given that the graph would start/end a little-way on past each side of indels. This is also why you don't need to do indel realignment with the HaplotypeCaller (this was needed in the old UnifiedGenotyper).

**Edit**: For more details regarding the HaplotypeCaller, see this nice page on GATK's website, which goes into much more detail than I did here.
> 15 votes

---
Tags: sam, samtools
---


# Question
Title: How can I call structural variants (SVs) from pair-end short read resequencing data?
I have a reference genome and now I would like to call structural variants from Illumina pair-end whole genome resequencing data (insert size 700bp).

There are many tools for SV calls (I made an incomplete list of tools bellow). There is also a tool for merging SV calls from multiple methods / samples - SURVIVOR. Is there a combination of methods for SV detection with optimal balance between sensitivity and specificity?

There is a benchmarking paper, evaluating sensitivity and specificity of SV calls of individual methods using simulated pair-end reads. However, there is no elaboration on the combination of methods.

List of tools for calling structural variants:

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

For that reason, I'd take two approaches:

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

# Answer
In case if you are really dedicated to obtain perfect results, you can use strategy described there, in 1000GP 3rd Phase SV detection paper - use these tools, validate your calls with IRS test, merge calls into one callset.

If you do not wanna spend thousands human-hours as was spent during this paper preparation, from my experience, it is better to use 1 paired-end insert distance method and one read-depth based method. Each of them cover "different" regions in the genome. (even thou they have huge overlap, paired-end detection requires both SV breakpoints to be located within the regions with good mappability which is not always the case, but resolution of read-depth methods is lower in general, paired-ends works well for deletions/tandem duplications/inversions, but have troubles with non-tandem duplications).

Hope it helps.
> 3 votes

---
Tags: ngs, structural-variation, variant-calling
---


# Question
Title: How are MACS2's narrow peak and broad peak algorithms different?
The peak calling tool MACS2 can call peaks in either narrow peak mode (for focused signals like transcription factor ChIPseq) or broad peak mode (for more defuse signals, like certain histone modifications).

The algorithm for narrow peak calling is well described in the MACS publication. But I don't find much documentation for how peak calling is different in broad peak mode. The manual only contains the following:

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

But this doesn't really describe exactly how this is performed.

So what is the algorithm MACS uses for calling broad peaks?

# Answer
> 10 votes
The key function is `call_broadpeaks`:

The description attached to the function says:

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

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

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

...

This has a few implications:

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

---
Tags: algorithms, chip-seq
---


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

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

I am using the fgsea - Fast Gene Set Enrichment Analysis package to calculate the enrichment scores and someone told me that the numbers are different from the ones on the Broad Institute despite all the other parameters being equivalent.

Are these two methods (fgsea and Broad Institute GSEA) equivalent to calculate the enrichment score?

I looked to the algorithms of both papers, and they seem fairly similar, but I don't know if in real datasets they are equivalent or not.

Is there any article reviewing and comparing how does the enrichment score method affect to the result?

# Answer
> 7 votes
According to the FGSEA preprint:

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

Thus, FGSEA and GSEA are not identical.

And again in the conclusion:

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

The author argues that FGSEA is more accurate, so it can't be equivalent.

If you are interested specifically in the enrichment score, that was addressed by the author in the preprint comments:

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

So that part seems to be the same.

# Answer
> 4 votes
I found that the author responded that in the discussion in the preprint,

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

---
Tags: r, gse, bioconductor
---


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

What resources exist for this, and how do I do it?

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

Output in the web-browser:

If you download the annotated .vcf, frequencies will be in the `INFO` field:

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

```

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

# Answer
> 5 votes
The greatest **protein** coding variant catalogue is definitely ExAC (\>65k individuals). They also published a blogpost where they describe how to reproduce figures in the paper (it is a good start how to get familiar with the dataset).

For the **whole-genome** variants I would look at the data created by 1000 genomes project (the latest release has more than 3k individuals). The integrated variant call sets can be downloaded though the portal and catalogue of SVs can be found here.

In this paper (also 1000 genomes project) they speak about non-precise placement of SVs by SV callers. I would keep this in mind for the comparison of your genome to the known variants.

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

---
Tags: vcf, public-databases, variation, structural-variation, snv
---


# Question
Title: How to quickly determine mutations in a read of a sam file?
After DNA sequencing, I generated a sam file through alignment of a fastq file. Before using well known variant calling programs (eg. Annovar etc.), I want to pick some reads and know what kinds of mutations are there.

Is there any quick way of doing it?

# Answer
For qualitative analysis, you're probably better off using something less granular like IGV or IGB. However, if you really want to look at a couple of reads:

If you're willing to ignore sequencing errors, you can inspect the CIGAR string or the MD tag, both of which give information on the alignment of a single read.

The CIGAR string gives details on insertions, deletions, clippings, matches and mismatches. From Genome Analysis Wiki,

> The sequence being aligned to a reference may have additional bases that are not in the reference or may be missing bases that are in the reference. The CIGAR string is a sequence of of base lengths and the associated operation. They are used to indicate things like which bases align (either a match/mismatch) with the reference, are deleted from the reference, and are insertions that are not in the reference. For example:

```
RefPos:     1  2  3  4  5  6  7     8  9 10 11 12 13 14 15 16 17 18 19
Reference:  C  C  A  T  A  C  T     G  A  A  C  T  G  A  C  T  A  A  C
Read:                   A  C  T  A  G  A  A     T  G  G  C  T
With the alignment above, you get:
POS: 5
CIGAR: 3M1I3M1D5M

```

Most common usage of the CIGAR string uses M (match/mismatch), I (insertion), D (deletion), S (soft clipping), and H (hard clipping). Note that = (match) and X (mismatch) are available as alternatives to the less informative M, but they are less widely used.

The MD tag gives specific details on mismatches and deletions. From the SAMtools tags specification,

> The MD field aims to achieve SNP/indel calling without looking at the reference. For example, a string ‘10A5^AC6’ means from the leftmost reference base in the alignment, there are 10 matches followed by an A on the reference which is different from the aligned read base; the next 5 reference bases are matches followed by a 2bp deletion from the reference; the deleted sequence is AC; the last 6 bases are matches. The MD field ought to match the CIGAR string.

Note that neither of these will give you any idea of structural variants in short reads, and neither will be particularly readable (or helpful, due to the higher error rate) in long reads.
> 8 votes

# Answer
I wrote a program, ASCIIGenome, that I find handy in cases where you want to have a quick look at genomic data. It's a genome browser for the command line.

To view only reads containing mismatches you can use the internal function `awk`. To filter for reads where the NM tag (number of mismatches) is \>0:

```
ASCIIGenome -fa genome.fa aln.bam
...

[h] for help: awk 'getSamTag("NM") > 0'

```

The view on terminal screen may look something like this:

Similarly, to get only reads containing indels you can use `awk '$6 ~ "D|I"'`

Hope this helps and feel free to report bugs & issues.
> 7 votes

# Answer
`samtools mpileup` can do this quickly:

```
samtools mpileup -f reference.fasta -uv input.sam > variants.vcf

```

This will produce a VCF-formatted file containing information about what variants have been seen in the SAM file, aggregated for all the mapped reads.
> 4 votes

---
Tags: variant-calling
---


# Question
Title: Understanding DESeq2 design, contrast and results
I have a set of high-troughput experiments with 2 genotypes ("WT" and "prg1") and 3 treatments ("RT", "HS30" and "HS30RT120"), and there are 2 replicates for each of the genotype x treatment combinations.

The read counts for the genes are summarized in a file that I load as follows in R:

```
> counts_data <- read.table("path/to/my/file", header=TRUE, row.names="gene")
> colnames(counts_data)
 [1] "WT_RT_1"          "WT_HS30_1"        "WT_HS30RT120_1"   "prg1_RT_1"
 [5] "prg1_HS30_1"      "prg1_HS30RT120_1" "WT_RT_2"          "WT_HS30_2"
 [9] "WT_HS30RT120_2"   "prg1_RT_2"        "prg1_HS30_2"      "prg1_HS30RT120_2"

```

I describe the experiments as follows:

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

```

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

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

or:

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

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

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

However, this doesn't explain how to apply this "list style" to the "contrast" argument. And the above situation seems to be asymmetrical. By that I mean that genotype and condition do not seem to have an interchangeable role.

So I tried the following more symmetric formula:

```
> dds <- DESeqDataSetFromMatrix(
                                countData = counts_data,
                                colData = col_data,
                                design = ~ geno + treat + geno:treat)
> dds <- DESeq(dds)

```

Now, can I for instance **get the differential expression results when comparing treatment "HS30" against "RT" as a reference, in genotype "prg1"**?

And how?

If I understand correctly, the above-mentioned "list style" uses names given by the `resultsNames` function. In my case, I have the following:

```
> resultsNames(dds)
[1] "Intercept"               "geno_WT_vs_prg1"
[3] "treat_HS30RT120_vs_HS30" "treat_RT_vs_HS30"
[5] "genoWT.treatHS30RT120"   "genoWT.treatRT"

```

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

I'm lost.

# Answer
> 10 votes
It seems that the "combining factors" trick described in part 3.3 of DESeq2 current "vignette" (as of may 2017) under the title "Interaction" is a way to access to the desired contrasts.

It seems possible do do it directly when building the `colData` and when calling `DESeqDataSetFromMatrix`:

Let's add a combined "geno" and "treat" factors to the future `colData` parameter:

```
> col_data$geno_treat <- as.factor(paste(col_data$geno, col_data$treat, sep="_"))
> col_data
DataFrame with 12 rows and 4 columns
                        geno       treat         rep     geno_treat
                 <character> <character> <character>       <factor>
WT_RT_1                   WT          RT           1          WT_RT
WT_HS30_1                 WT        HS30           1        WT_HS30
WT_HS30RT120_1            WT   HS30RT120           1   WT_HS30RT120
prg1_RT_1               prg1          RT           1        prg1_RT
prg1_HS30_1             prg1        HS30           1      prg1_HS30
...                      ...         ...         ...            ...
WT_HS30_2                 WT        HS30           2        WT_HS30
WT_HS30RT120_2            WT   HS30RT120           2   WT_HS30RT120
prg1_RT_2               prg1          RT           2        prg1_RT
prg1_HS30_2             prg1        HS30           2      prg1_HS30
prg1_HS30RT120_2        prg1   HS30RT120           2 prg1_HS30RT120

```

We can now use a design where differential expression will be explained by these combined factors:

```
> dds <- DESeqDataSetFromMatrix(
    countData = counts_data,
    colData = col_data,
    design = ~ geno_treat)

```

We run the analysis:

```
> dds <- DESeq(dds)

```

Then we can query results for a particular contrast between such factor combinations. For instance, to have the results for the effect of treatment "HS30" against the reference state "RT" in genotype "prg1":

```
res <- results(dds, contrast=c("geno_treat", "prg1_HS30", "prg1_RT"))

```

# Answer
> 16 votes
The simplest manner is to not use a wald test, but rather an LRT with a reduced model lacking the factor of interest:

```
dds = DESeq(dds, test="LRT" reduced=~geno+geno:Treatment)

```

The above would give you results for Treatment regardless of level while still accounting for a possible interaction (i.e., a "main effect of treatment, regardless of the type of treatment").

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

---
Tags: r, bioconductor, deseq2, differential-expression
---


# Question
Title: Confirm success or failure of RNA-Seq normalization
I am working with a set of (bulk) RNA-Seq data collected across multiple runs, run at different times of the year. I have normalized my data using library size / quantile / RUV normalization, and would like to check (quantitatively and/or qualitatively) whether or not normalization has succeeded in removing the batch effects.

It is important to note that by "normalization has succeeded", I simply mean that unwanted variation has been removed - further analysis is required to ensure that biological variation has not been removed. What are some plots / statistical tests / software packages which provide a first-step QC for normalization?

# Answer
> 17 votes
You should use box plots and PCA plot. Let's take a look at the RUV paper:

> https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4404308/

**Before normalization and after UQ normalization**:

> Libraries do not cluster as expected according to treatment. ... for UQ-normalized counts. UQ normalization does not lead to better clustering of the samples...

Before normalization, the medians in the box-plot obviously look very different among replicates.

After UQ normalization, the medians look closer but Trt.11 look like an outlier. Furthermore, the treatments aren't clustered on the PCA plot. Since they are replicates, you'd like them be close on the plot.

**After RUV normalization**

> ... RUVg shrinks the expression measures for Library 11 towards the median across libraries, suggesting robustness against outliers. ... Libraries cluster as expected by treatment. ...

The RUV has made the distribution more robust and the samples closer on the PCA plot. However, it's still not perfect as one of the treatments is not close to the other two on the first PC.

The vignettes for Bioconductor `RUVSeq` describes the two functions: `plotRLE` and `plotPCA`.

# Answer
> 8 votes
Visual inspection with histograms, boxplots, or some other distribution visualization is the way to go. Prior to normalization, your abundances may look something like this.

Post-normalization, they should look something like this.

See this blog post for example code.

---
Tags: rna-seq, normalization
---


# Question
Title: Designing a lab NGS file database schema
I am the resident Bioinfo Geek in a hospital academic lab that routinely employs NGS as well as CyTOF and other large volume data producing technologies. I am sick of our current "protocol" for metadata collection and association with the final products (miriad excel sheets and a couple poorly designed RedCap DBs).

I want to implement a central structured, controlled datastore that will take care of this. I know that the interface to the technicians how will be inputing the data is crucial to its adoption, but this is not the focus of THIS particular question: **Does there exist a schema or schema guidelines for this type of database?**

I would rather use a model that has been developed by people who know how to do this well. I know of BioSQL but it seems more geared towards full protein/nucleotide records like those found in uniprot or genbank. That is not what we have here. What I want is something similar to the system touched on in this preprint: http://biorxiv.org/content/early/2017/05/10/136358

Alternatively, can anyone provide links to where I might find relevant guidelines or supply personal advice?

# Answer
> 5 votes
For *metadata*, I would use a SQL schema something like the following:

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

```

In the schema, you have `Project` and biological `Sample` tables, which are independent of each other at the high level. An `AnalysisSample` describes a sample used in BAM or VCF and connects `Project` and biological `Sample`. Importantly, each `AnalysisSample` has a symbol unique in a project (see the primary index). This is the symbol on a BAM read group line or on a VCF sample line. A `Collection` is in effect a BAM/CRAM file. In theory, a BAM file may contain more than one samples (though rare in practice), which is addressed by a separate `ReadGroup` table. Finally, a `VariantSet` is a VCF file. `VariantSample` tells you which samples are included in each VCF file.

This is the skeleton of a full schema. You can add extra fields to appropriate tables (e.g. file path and hg19/hg38/etc to `Collection`, read length to `ReadGroup` and family ID to `Sample`). You also need indices for efficient table joining and perhaps more tables for complex structures (e.g. pedigree).

For the projects I have participated, this schema should work most of time. It is inspired by GA4GH's JSON schema, but my version is in SQL, is simpler and also has a slightly different structure which I think is better.

# Answer
> 9 votes
The Global Alliance for Genomics and Health has been working on the issue of representing sequencing data and metadata for storage and sharing for quite some time, though with mixed results. They do offer a model and API for storing NGS data in their GitHub repository, but it can be a bit of a pain to get a high-level view. I am not sure if any better representation of this exists elsewhere.

I can say from personal experience (having built over a dozen genomic databases), there is no ideal data model and storage best practices. Genomic data comes in many shapes and sizes, and your needs are going to vary from every other organization, so what works for one bioinformatics group won't necessarily work for you. The best thing to do is design and implement a model that will cover all of the data types in your workflow and downstream analyses you might do with the data and metadata.

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

---
Tags: ngs, database
---


# Question
Title: How to make a distinction between the "classical" de Bruijn graph and the one described in NGS papers?
In Computer Science a De Bruijn graph has (1) `m^n` vertices representing all possible sequences of length `n` over `m` symbols, and (2) directed edges connecting nodes that differ by a shift of `n-1` elements (the successor having the new element at the right).

However in Bioinformatics while condition (2) is preserved, what is called a De Bruijn graph doesn't seem to respect condition (1). In some cases the graph doesn't look anything like a de Bruijn graph at all (e.g. http://genome.cshlp.org/content/18/5/821.full).

So my question is, if I want to make it explicit that I am using the Bioinformatics interpretation of a de Bruijn graph, is there a term for it? Something like "simplified de Bruijn graph", "projection of a de Bruijn graph", or "graph of neighbouring k-mers"? Are there any papers making this distinction, or did I get it all wrong?

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

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

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

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

The oldest reference I found for a specific term to refer to the assembly-related structure is Skiena and Sundaram (1995), where they call it a **subgraph of the de Bruijn digraph**. Later, in 2002, Błażewicz et al. will refer to it as a **de Bruijn induced subgraph**. The term **de Bruijn subgraph** is also formally defined in Quitzau’s thesis (2009). There, and also in the article (Quitzau and Stoye, 2008) the authors describe the **sequence graph** as a modification of the sparse de Bruijn subgraph (commonly used in assembly problems), where non-branching paths are replaced by a single vertex. The term **sparse de Bruijn graph** is also used by Chauve et al. (2013).

Another term that I found was **word graph**, described by both Malde et al. (2005) and by Heath and Pati (2007) as a *subgraph* or as a *generalization* of a de Bruijn graph. Rødland (2013) summarises some of the terms used for this data structure:

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

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

# Answer
> 4 votes
In addition to the regular De Bruijn graph as depicted on the wikipedia, some implementations in bioinformatics feature additional processing. I guess the main reason figure 1 in the paper you linked (concerning the Velvet genome assembler) is slightly different is that a node represents *a series of overlapping k-mers*. In order to visualize this as a more classic De Bruin graph you would have to connect the k-mers depicted *above* the nodes. The caption next to figure one describes the processing quite clearly.

As per your last question: I don't think there is a 'Bioinformatic interpretation of a De Bruijn graph'. There are different implementations, which all have there specifics. Thus it would be best to refer to the actual implementation.

As an example: this is a nice paper on how to construct a pan-genome De Bruijn graph of multiple genomes simultaneously.

# Answer
> 3 votes
Let's first assume DNA only has one strand. An assembly de Bruijn graph is a subgraph of a complete de Bruijn graph. It contains a vertex u if u is a k-mer in reads; it contains an edge u-\>v, if u and v are adjacent k-mers on a read. Alternatively, we note that an edge u-\>v is represented by a (k+1)-mer. An assembly de Bruijn graph can be considered a subgraph edge induced from all (k+1)-mers in reads – in fact, some assemblers take the list of (k+1)-mer as a succinct representation of de Bruijn graphs.

DNA has two strands. We just need to induce an assembly de Bruijn graph from all (k+1)-mers and their reverse complement. It is still a subgraph of a complete de Bruijn graph.

Because an assembly de Bruijn graph is just a subgraph. It is not necessary to give it a new name.

PS: I deleted my old answer as that was not what you are asking for based on your comments. I was confused by your mentioning velvet. Velvet uses an equivalent but uncommon representation of de Bruijn graphs, which complicates your question.

---
Tags: ngs, assembly, k-mer, de-bruijn-graphs
---


# Question
Title: The state, limitations and comparisons of large variant stores
**Background:** We're increasingly needing some way of storing lots of variant data associated with lots of subjects: think clinical trials and hospital patients, looking for disease-causing or relevant genes. A thousand subjects is where we'd start, there's talk of millions on the horizon. With various genomic medicine initiatives, this is likely a wider need.

**The problem:** While there's plenty of platforms out there, it's a rapidly evolving field. It's difficult to get a feel for how (and if) they perform and how they line up against each other:

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

**Solutions I've seen so far:**

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

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

# Answer
> 14 votes
An epic question. Unfortunately, the short answer is: no, there are no widely used solutions.

For several thousand samples, BCF2, the binary representation of VCF, should work well. I don't see the need of new tools at this scale. For a larger sample size, ExAC people are using spark-based hail. It keeps all per-sample annotations (like GL, GQ and DP) in addition to genotypes. Hail is at least something heavily used in practice, although mostly by a few groups so far.

A simpler problem is to store genotypes only. This is sufficient to the majority of end users. There are better approaches to store and query genotypes. GQT, developed by the Gemini team, enables fast query of samples. It allows you to quickly pull samples under certain genotype configurations. As I remember, GQT is orders of magnitude faster than google genomics API to do PCA. Another tool is BGT. It produces a much smaller file and provides fast and convenient queries over sites. Its paper talks about ~32k whole-genome samples. I am in the camp who believe specialized binary formats like GQT and BGT are faster than solutions built on top of generic databases. I would encourage you to have a look if you only want to query genotypes.

Intel's GenomicDB approaches the problem in a different angle. It does not actually keep a "squared" multi-sample VCF internally. It instead keeps per-sample genotypes/annotations and generates merged VCF on the fly (this is my understanding, which could be wrong). I don't have first-hand experience with GenomicDB, but I think something in this line should be the ultimate solution in the era of 1M samples. I know GATK4 is using it at some step.

As to others in your list, Gemini might not scale that well, I guess. It is partly the reason why they work on GQT. Last time I checked, BigQuery did not query individual genotypes. It only queries over site statistics. Google genomics APIs access individual genotypes, but I doubt it can be performant. Adam is worth trying. I have not tried, though.

---
Tags: genomics, variation, human-genome, data-management
---


title: Influential papers in bioinformatics

text:
# Question
Title: Influential papers in bioinformatics
As a newcomer, what are some influential research papers published in the past few years that I should read to catch up in the field?

Edit: @Devon Ryan answered it perfectly. I tried a google search but it turned up a bunch of opinionated articles that varied. The biostar discussions are perfect and was the link I was searching for. I hope someone else new to the field finds this useful!

# Answer
> 4 votes
Over on biostars there's a thread like this every year or so. I'll link to the 2016 edition and the (much shorter) 2015 edition.

My personal picks from those would be:

---
Tags: literature-search
---


# Question
Title: The effects of incomplete bisulfite conversion upon mapping efficiency
*This question has also been posted on Biostars*

I have sequenced numerous multiplexed pools of BS amplicon-seq libraries derived from human samples on a MiSeq over the past few weeks. I have been utilising trim-galore and Bismark for alignment and am finding the mapping efficiency to be really low for two pools (55% & 30% respectively) both of which had a cytosine per base sequence of around 10-20% throughout the entire read in their fastqc files.

The high C content visible in the fastqc files makes me think the poor mapping efficiency is due to poor bisulfite conversion, as this would be expected to be close to zero is bisulfite conversion had actually took place.

I am new to the field so any help would be greatly appreciated.

# Answer
> 7 votes
Bisulfite conversion efficiency has no effect on the mapping rate in bismark and similar tools. The reason is that the reads are fully bisulfite converted *in silico* before alignment to minimize mapping bias.

I would suggest that you play around with the settings handed to bowtie2, such as using local alignment and modifying the `--score-min` option to allow more mismatches. Alternatively, you might try a different aligner like bwameth, which will always do local alignment.

---
Tags: sequence-alignment, fastq
---


# Question
Title: How can I extract normalized read count values from DESeq2 results?
The results obtained by running the `results` command from DESeq2 contain a "baseMean" column, which I assume is the mean across samples of the normalized counts for a given gene.

**How can I access the normalized counts proper?**

I tried the following (continuing with the example used here):

```
> dds <- DESeqDataSetFromMatrix(countData = counts_data, colData = col_data, design = ~ geno_treat)
> dds <- DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
> res <- results(dds, contrast=c("geno_treat", "prg1_HS30", "prg1_RT"))

```

Here is what I have for the first gene:

```
> res["WBGene00000001",]$baseMean
[1] 181.7862
> mean(assays(dds)$mu["WBGene00000001",])
[1] 231.4634
> mean(assays(dds)$counts["WBGene00000001",])
[1] 232.0833

```

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

# Answer
> 13 votes
The normalized counts themselves can be accessed with `counts(dds, normalized=T)`.

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

# Answer
> 6 votes
It depends what you mean by “normalised”. As Devon said, the `normalized = TRUE` argument to the `count` function gives you normalised counts. However, these are “only” library-size normalised (i.e. divided by the `sizeFactors(dds)`).

However, as the vignette explains, downstream processing generally requires more advanced normalisation, to account for the heteroscedasticity of the counts. This is often done by simply `log`ging the counts but this has obvious drawbacks (most trivially, what do we do with 0 counts? A workaround is to add a pseudocount but that’s problematic too).

DESeq2 offers two different methods to perform a more rigorous analysis:

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

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

---
Tags: r, bioconductor, deseq2, normalization
---


# Question
Title: Can I model technical replicates in DESeq2?
I’d normally use `collapseReplicates` (or do the collapsing upstream) to handle technical replicates.

However, in my current RNA-seq experimental design, samples were sequenced twice using different library preparation protocols, leading to marked differences in the resulting count estimates (in fact, the choice of different protocol explains most of the variance in the bulk data according to PCA). I would therefore like to model these differences by adding them as a covariate to the DESeq2 model.

I *thought* I could just add another factor to the design formula, equivalent to the “pasilla” example in the DESeq2 vignette, which adds the factor `type` for the library type (single-end vs paired-end). However, the pasilla vignette states that these different libraries are actually independent *biological* replicates, not technical replicates as I had previously assumed.

I’m now concerned that having technical replicates in the design might artificially inflate the power. Can the design be salvaged?

It’s worth noting that we *also* have biological replicates, and a full rank design; i.e. every combination of treatment and library prep, with biological replicates.

# Answer
> 6 votes
If they're truly technical replicates, then there's no way to model them using DESeq2\*, as you've alluded to with the `collapseReplicates` function. DESeq2/ Mike Love's general recommendation with `collapseReplicates` is to just add the reads together for technical replicates.

If you want to model them instead of collapsing them down, you can `voom` transform your data and follow an example similar to section 11.3 in the Limma users guide.

**(Code example from 11.3):**

*Design*

```
| FileName  | Cy3   | Cy5  |
| --------- |:-----:| ----:|
| File1     | wt1   | mu1  |
| File2     | wt1   | mu1  |
| File3     | wt2   | mu2  |
| File4     | wt2   | mu2  |

```

*Example*

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

```

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

# Answer
> 5 votes
Practically speaking, there's no way to include the technical replicates in that design (in DESeq2 at least). Your concern regarding inflating the power is exactly correct and the only way to combat that would be to add a pairing factor like one might do with case-control or tumor-normal studies. That is, something like:

```
  group libraryPrep sample
1    WT           A      1
2    WT           B      1
3   MUT           A      2
4   MUT           B      2
5    WT           A      3
6    WT           B      3
7   MUT           A      4
8   MUT           B      4

```

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

---
Tags: rna-seq, deseq2
---