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

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

# Answer
> 29 votes
> What does this soft masking actually mean?

A lot of the sequence in genomes are repetitive. Human genome, for example, has (at least) two-third repetitive elements.\[1\].

These repetitive elements are soft-masked by converting the upper case letters to lower case. An important use-case of these soft-masked bases will be in homology searches: An `atatatatatat` will tend to appear both in human and mouse genomes but is likely non-homologous.

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

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

> What does a lowercase n represent?

UCSC uses Tandom Repeat Finder and RepeatMasker for soft-masking potential repeats. NCBI most likely uses TANTAN. 'N's represents no sequence information is available for that base. It being replaced by 'n' is likely an artifact of the repeat-masking software where it soft-masks an 'N' by an 'n' to indicate that portion of the genome is likely a repeat too.

\[1\] http://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1002384

# Answer
> 13 votes
The use of lower/upper case letters and `N`/`n` letters in genomes sequences is not completely standardised and you should always check the specification of the resource you are using.

Lower case letters are most commonly used to represent “soft-masked sequences”, a convention popularised by RepeatMasker, where interspersed repeats (which covers transposons, retrotransposons and processed pseudogenes) and low complexity sequences are marked with lower case letters. Note that larger repeats, such as sizable tandem repeats, segmental duplications, and whole gene duplications are not generally masked.

However, there are other uses for lower/upper case letters, for example, Ensembl have used upper/lower case letters to represent exonic and intronic sequences respectively.

`N` and `n` nucleotides may represent “hard masked sequences”, where interspersed repeats and low complexity sequences are replaced by `N`s. But `N`/`n`s may alternatively represent ambiguous nucleotides, indeed this is the IUPAC specification.

Also note occasionally (although fortunately rarely) `X`/`x` is used to represent ambiguous nucleotides or “hard-masked sequences” too.

# Answer
> 4 votes
1. Lower case nucleotides commonly denotes a **soft** masked sequences. How exactly the genome was masked you can find in FAQ of NCBI:

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

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

**An example of usage of the soft mask**

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

---
Tags: fasta, genome
---


# Question
Title: samtools mpileup empty when filtering out flags
I produced a bam file by aligning reads to a small set of synthetic sequences using bwa-mem.

I am heavily filtering reads that are not paired and of a certain orientation. Applying the filtering, I get a few thousands of reads:

```
samtools view -h $myfilebam | \
samtools view -h -F4 - | \
samtools view -h -F8 - | \
samtools view -h -F256 - | \
samtools view -h -F512 - | \
samtools view -h -F1024 - | \
samtools view -h -F2048 - | \
samtools view -h -f16 - | \
samtools view -h -f32 -  | wc -l

```

Gives me `89502` reads.

If I then pipe this into `samtools mpileup`, I get no results:

```
samtools view -h $myfilebam | \
samtools view -h -F4 - | \
samtools view -h -F8 - | \
samtools view -h -F256 - | \
samtools view -h -F512 - | \
samtools view -h -F1024 - | \
samtools view -h -F2048 - | \
samtools view -h -f16 - | \
samtools view -h -f32 -  | \
samtools mpileup --excl-flags 0 -Q0 -B -d 999999 - | wc -l

```

Returns 0.

I tried different combinations of filtering, and when I do both `-f 16` and `-f 32` returns empty, but if I do either of those, then it works:

```
samtools view -h $myfilebam | \
samtools view -h -F4 - | \
samtools view -h -F8 - | \
samtools view -h -F256 - | \
samtools view -h -F512 - | \
samtools view -h -F1024 - | \
samtools view -h -F2048 - | \
samtools view -h -f16 - | \
samtools mpileup --excl-flags 0 -Q0 -B -d 999999 - | wc -l

```

Returns `1056`.

Any ideas why? My thinking was that it would work with `--excl-flags 0`.

EDIT: substituting `mpileup` for `depth` does work, and prints out each position and the depth as expected.

EDIT2: adding `-q 0` to `mpileup` gives the same empty result.

Thanks in advance

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

```
...
samtools view -h -f32 -  | \
samtools view | wc -l

```

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

```
...
samtools mpileup -A -Q0 -B -d 999999 - | wc -l

```

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

---
Tags: sam, samtools, mpileup
---


# Question
Title: Classifying samples based on marker gene expression
I have a few sets of marker genes that I can classify RNA-seq samples using semi-supervised clustering. I would like to automate the process, however, I am struggling to find the ideal algorithm that could generate some kind of score for marker gene set from a given sample.

I presume that this is a standard analyses in many groups but I am not sure which method(s) are yielding good results in practice.

# Answer
> 7 votes
I would consider using gene expression signatures to classify samples (especially cancer subtypes but the same principles apply to other problems of this type) one of the classic problems of bioinformatics. Quite a lot of work has been done on methods to derive gene sets that provide good classification performance. This is slightly different from your problem since you already have a gene signature but it may still prove useful.

These methods will typically fit a model that selects a (small) number of genes from genome-wide expression data that distinguish between the cell types/conditions in question, i.e. they derive a gene signature. The resulting model then allows for the classification of new samples. I've had success using GeneRave for this purpose (but note that this was designed for microarray data, I haven't used it with RNA-seq data and don't know how well it holds up there). A more recent paper relating to this issue can be found here.

So how does that help you? One option would be to fit one of these classifiers to gene expression data for the genes you already know to obtain a model that can then be applied to new samples automatically.

---
Tags: rna-seq
---


# Question
Title: What is the difference between ENCODE Tier 1, 2 and 3 cell types?
The ENCODE Experiment Matrix at UCSC lists the different available cell types under the categories "Tier 1", "Tier 2" and "Tier 3". What is the difference between these classifications?

What, for example, makes GM12878 a Tier 1 cell type and A549 a Tier 2 cell type?

# Answer
The tiers denote the original priority of sequencing and processing the samples. This is understandable given the number of marks and experiment types that ENCODE tried to sequence. Further details are available on ENCODE's website.

**Edit**: Just to expand a bit, tier one cells were supposed to get processed first. Then tier 2 came second. Of course tier 2 had a LOT more cells, so it got split into 2 and 2.5, since there were already apparently tier 3 cells.
> 6 votes

# Answer
From Encode's website:

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

---
Tags: encode
---


# Question
Title: Large-scale gRNA design for a CRISPR screen
What are the best tools to design gRNAs in a high-throughput way for a CRISPR screen, e.g. targeting all protein-coding genes in a genome? I would like to take into account possible off-target effects, as well as to allow for flexibility in the PAM sequence, to have the possibility of using non-canonical PAMs.

I'm aware of CRISPRseek, but I'm afraid it might be quite slow to run for a genome-scale search.

# Answer
> 2 votes
Perhaps GenomeCRISPR, I've personally never used it, but it has a well-documented API and appears you could automate it to go through the whole genome.

---
Tags: crispr
---


# Question
Title: Does GISTIC (v 2.0) estimate amplified/deleted probabilities on a single sample basis?
Does GISTIC 2.0 estimate the background model:

`G = -log(Probability | Background)`

by permuting within the sample or across all samples in the set?

The paper describes the probabilistic scoring method based on permutations, but I could not understand if this permutation is performed within the sample only. The documentation page seems to suggest across samples, but this would mean that different sized sets might lead to different outcomes on the sample sample.

Basically, does it matter if the set is constituted by, for example, 10 samples of the same tissue of origin or ~1000 from multiple tissues?

Thanks,

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

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

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

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

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

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

However, the relevant section (of the supplementary materials) seems to be "Null Hypothesis Generation: An Analytic Derivation of the Null Distribution", where it describes the semiexact approximation used.

In a supplementary file it describes it as "Generate all permutations of SNP labels within each sample to simulate datasets with random aberrations"

**Conclusion**:

It seems that the background probabilities are calculated within sample.

I can't say if it matters how the sets are constituted but I would say that the broader the sets are, the better estimation of the structural variations it will perform but the same estimation of the background will be done.

Ultimately you can check the code of the program, or test with 10 samples and replace one of them to see if the results are changing accordingly.
> 3 votes

---
Tags: cnv
---


# Question
Title: variant calling on ChIP-seq style data: samtools mpileup with minimal filters
I am running samtools mpileup (v1.4) on a bam file with very choppy coverage (ChIP-seq style data). I want to get a first-pass list of positions with SNVs and their frequency as reported by the read counts, but no matter what I do, I keep getting all SNVs filtered out as not passing QC.

What's the magic parameter set for an initial list of SNVs and frequencies?

EDIT: this is a question I posted on "the other" website, but didn't get a reply there.

# Answer
I used this in the past for ChIP-seq data and it generated SNVs:

```
samtools mpileup \
--uncompressed --max-depth 10000 --min-MQ 20 --ignore-RG --skip-indels \
--fasta-ref ref.fa file.bam \
| bcftools call --consensus-caller \
> out.vcf

```

This was samtools 1.3 in case that makes a difference.
> 7 votes

# Answer
Another approach is htsbox. You can get a candidate list with:

```
htsbox pileup -Cvcf ref.fa -q20 -Q20 -s5 file.bam > out.vcf

```

Here, `-q` sets min mapping quality, `-Q` sets min base quality, `-v` outputs variants only `-c` outputs VCF, `-C` gives you base counts on both strands and finally `-s5` requires at least 5 high-quality bases to call out an allele. It is useful when your data are failing the assumptions made by typical variant callers.

Why not samtools+bcftools or varscan? Transparency and speed. This command line simply counts based on the parameters you use. It does not apply additional operations. And because of this, it is over an order of magnitude faster than samtools mpileup or varscan. It is worth noting that samtools uses BAQ by default, which reduces FPs occasionally. However, BAQ is not quite necessary for longer Illumina reads and it hurts sensitivity at the same time.
> 5 votes

---
Tags: variant-calling, samtools, chip-seq, mpileup
---


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

# Answer
**Simulators designed specifically for Oxford Nanopore:**

**General long read simulators:**

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

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

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

# Answer
The best nanopore read simulators would be associated with the best base-callers. For a base-caller to effectively model the DNA strand, it needs to take into account the expected underlying electrical model together with the associated signal noise (both in the time dimension as well as the amplitude dimension). In theory, it should be possible to reverse the algorithm and generate an electrical signal given an underlying sequence.

Unfortunately, I'm not aware of any tools that attempt to simulate nanopore reads at the electrical level. Any "nanopore read simulator" that concentrates only on base sequence would need to encompass all the possible base-calling software models that exist, which is an impossible task (particularly given how quickly ONT updates their own base callers).
> 3 votes

---
Tags: nanopore, genome-sequencing, software-recommendation
---


# Question
Title: How can I improve a long-read assembly with a repetitive genome?
I'm currently trying to assembly a genome from a rodent parasite, *Nippostrongylus brasiliensis*. This genome does have an existing reference genome, but it is highly fragmented. Here are some continuity statistics for the scaffolds of the current Nippo reference genome (assembled from Illumina reads):

```
Total sequences: 29375
Total length: 294.400206 Mb
Longest sequence: 394.171 kb
Shortest sequence: 500 b
Mean Length: 10.022 kb
Median Length: 2.682 kb
N50: 2024 sequences; L50: 33.527 kb
N90: 11638 sequences; L90: 4.263 kb

```

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

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

Canu seems to deal quite well with the first two types of repeats, despite the abundance of repetitive structure in the genome. Here's the unitigging summary produced by Canu on one of the assemblies I've attempted. Notice that about 30% of the reads either span or contain a long repeat:

```
category            reads     %          read length        feature size or coverage  analysis
----------------  -------  -------  ----------------------  ------------------------  --------------------
middle-missing        694    0.07     7470.92 +- 5552.00        953.06 +- 1339.13    (bad trimming)
middle-hump           549    0.05     3770.05 +- 3346.10         74.23 +- 209.86     (bad trimming)
no-5-prime           3422    0.33     6711.32 +- 5411.26         70.92 +- 272.99     (bad trimming)
no-3-prime           3161    0.30     6701.35 +- 5739.86         87.41 +- 329.42     (bad trimming)

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

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

```

However, the third type of repeat is giving me a bit of grief. Using the above assembly, here are the continuity parameters from the assembled contigs:

```
Total sequences: 3505
Total length: 322.867456 Mb
Longest sequence: 1.762243 Mb
Shortest sequence: 2.606 kb
Mean Length: 92.116 kb
Median Length: 42.667 kb
N50: 417 sequences; L50: 194.126 kb
N90: 1996 sequences; L90: 35.634 kb

```

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

The repetitive regions are typically over 500bp in length, average about 3kb, and I've seen at least one case which seems to be a 20kb sequence duplicated in multiple regions.

The Canu defaults seem to give the best assembly results for the few parameters that I've tried, with one exception: trimming. I've tried playing around a little bit with the trimming parameters, and curiously a trimming coverage of 5X (with overlap of 500bp) seems to give a more contiguous assembly than with a trimming coverage of 2X (with the same overlap).

If anyone is interested in having a look at these data themselves, called FASTQ files from Nippo sequencing runs can be found here. Raw nanopore signal files are available within ENA project PRJEB20824. There's also a Zenodo archive here that contains the GFA and assembly contigs.

I can use Illumina data to correct the Canu assembly, but that doesn't help with resolving the "type 3" repeats. The regions are sufficiently similar that illumina reads get mapped to multiple points in the genome. The Illumina contigs are high quality (i.e. they have good BUSCO scores, indicating few variant errors), but quite short. Any sniff of a repeat and the contig ends. I've got more than a few examples of regions that would make an Illumina read (even 10x linked reads) cower in fear.

Does anyone have any other suggestions on how I could resolve these complex repeats?

Computational solutions would be preferred, but resequencing is not out of the question.

# Answer
"A few" 100kb reads won't help much. You need to apply the ultra-long protocol, which is different from the standard protocol.

You can't resolve 20kb near identical repeats/segdups with 10kb reads. All you can do is to bet your luck on a few excessively long reads spanning some units by chance. For divergent copies, it is worth looking at this paper. It uses Illumina reads to identify k-mers in unique regions and ignores non-unique k-mers at the overlapping stage. The paper said that this strategy is better than using standard overlappers, which I buy, but probably it can't resolve a 20kb segdup with a handful of mismatches, either.

Such mismatch-based approaches always have limitations and may not work for recent segdups/repeats. The ultimate solution is to get long reads, longer than your repeat/segdup units. The ~100kb reads in the recent preprint will be a game changer for you. If your ~20kb repeats are not tandem, 10X's ~100kb linked reads may help, too.
> 7 votes

---
Tags: assembly, nanopore, long-reads, canu, repeat
---


# Question
Title: Quantifying reads mapping to multiple loci
I have been using STAR for our RNA-Seq samples. The `final.out` log file reports percentage of uniquely mapped reads along with percentage of reads that map to `multiple loci` (less than or equal to 10) and percentage of reads mapping to `too many loci` (greater than 10). However, I want to break down the `multiple loci` part to individual counts: Reads mapping to 2 locations, 3 locations, 4 locations .. 10 locations.

The `NH` tag seems to be used by `STAR`. However a naive read counting approach results in it reporting more number of reads than total reads.

For example, my `final.out` looks like this:

```
   Mapping speed, Million of reads per hour |       1403.36

                      Number of input reads |       53015978
                  Average input read length |       26
                                UNIQUE READS:
               Uniquely mapped reads number |       368916
                    Uniquely mapped reads % |       0.70%
                      Average mapped length |       26.45
                   Number of splices: Total |       1057
        Number of splices: Annotated (sjdb) |       0
                   Number of splices: GT/AG |       802
                   Number of splices: GC/AG |       1
                   Number of splices: AT/AC |       0
           Number of splices: Non-canonical |       254
                  Mismatch rate per base, % |       0.31%
                     Deletion rate per base |       0.00%
                    Deletion average length |       1.45
                    Insertion rate per base |       0.00%
                   Insertion average length |       1.00
                         MULTI-MAPPING READS:
    Number of reads mapped to multiple loci |       45766732
         % of reads mapped to multiple loci |       86.33%
    Number of reads mapped to too many loci |       3757890
         % of reads mapped to too many loci |       7.09%
                              UNMAPPED READS:
   % of reads unmapped: too many mismatches |       0.00%
             % of reads unmapped: too short |       5.89%
                 % of reads unmapped: other |       0.00%

```

Counting histogram of number of positions a read maps to using `pysam`:

```
def get_reads_hist(bam):
    bam = pysam.AlignmentFile(bam, 'rb')
    counts = Counter()
    for query in bam.fetch():
        nh_count = Counter(dict(query.get_tags())['NH'])
        counts += nh_count
    return counts

```

results in

```
Counter({1: 330606,
     2: 86772164,
     3: 329,
     4: 38083,
     5: 31,
     6: 1094,
     7: 129,
     8: 50,
     10: 50})

```

The count `1` reads are fine even though they do not match the counts in `final.out` file since I am counting a certain category of reads (say those mapping to `tRNA` only), but the reads mapping to 2 locations are highly overestimated. Why is that?

# Answer
You almost had the correct python code already, you just need to filter out secondary alignments:

```
def get_reads_hist(bam):
    bam = pysam.AlignmentFile(bam, 'rb')
    counts = Counter()
    for query in bam.fetch():
        if query.is_secondary:
            continue
        nh_count = Counter(dict(query.get_tags())['NH'])
        counts += nh_count
    return counts

```
> 9 votes

---
Tags: rna-seq, sam, samtools, star, pysam
---


# Question
Title: How to correlate two zero inflated bedgraph-like signals?
This question pertains to iCLIP, but it could just as easily be ChIP-seq or ATAC-seq or mutation frequencies.

I have iCLIP read counts across the transcriptome and I wish to know if the signals are correlated - that is, where one of them is high, is the other likely to be high.

Often when dealing with such data (e.g. iCLIP data) we know that the data is generally sparse - that is at most positions both signals are zero and this is correct, and also zero-inflated - that is many bases that "should" have a signal are missing that data. So just calculating the Spearman's correlation is likely to give an artificially low value.

**What might be a way to asses the association? I should add that the aim is the assess association of binding patterns within genes, rather than (or as well as) between genes.**

Things I have thought of:

* Apply some sort of smoothing to the data (eg a rolling mean). Remove any bases with 0 in both samples. Compute the spearmans.
* Calculate the average pairwise distance between every read in sample one and every read in sample two. Compare this to data where the reads have been randomised within genes.

In the first case removing all bases with 0 in both samples seems wrong. But if 99.99% of all bases have zero in both samples, then this seems like its necessary for Spearman.

In the second case, the result seems like it would be non-intuitive to interpret. And also calculating this would be massively computationally intensive.

# Answer
> 4 votes
Honestly I'd just use `multiBigwigSummary` and then `plotCorrelation` from deepTools for this, but I'm a bit biased. There, the idea would be to consider each gene as a unit (you could instead use bins, but I don't think that would as nicely do what you want), namely by giving the tools a BED or GTF file input. It would then calculate the average signal in each gene/transcript and you could do your spearman's correlation. Features with 0 in all samples could optionally be removed (`plotCorrelation --skipZeros`).

While you certainly could go the whole 9 yards and use per-base comparisons, that seems a bit overkill and I suspect that it won't really yield appreciably more information (especially when one considers the additional time overhead).

# Answer
> 3 votes
Rather than working on the base level, you could probably work on say gene level counts. Kendall's tau, an ordinal association metric, can then be used as an appropriate correlation measure.

If $X$ and $Y$ are your iCLIP replicates, $i$ represents gene index and $(x\_i, y\_i)$ represents the number of RBP binding sites in $X$ and $Y$ respectively for the $i^{th}$ gene, Kendall's tau is defined as :

$$ \frac{\text{#(concordant pairs)} - \text{#(discordant pairs)}}{n(n-1)/2} $$

Where any two pairs $(x\_i, y\_i)$ and $(x\_j, y\_j)$ are concordant if:

* $x\_i \> x\_j$ AND $y\_i \> y\_j$

OR

* $x\_i \< x\_j$ AND $y\_i \< y\_j$

Correspondingly they are discordant if:

* $x\_i \< x\_j$ AND $y\_i \> y\_j$

OR

* $x\_i \> x\_j$ AND $y\_i \< y\_j$

# Answer
> 3 votes
It depends whether you want to treat the peak intensities as binary (comparing presence/absence of peaks in the sets) or continuous (comparing the relative magnitudes of the peaks).

**Binary**

For starting out, a simple binary comparison may be appropriate. You can use a peak caller of your choosing to identify peaks in each sample according to your desired criteria. Then you can use a similarity metric such as the Jaccard index to quantify the level of agreement among the peaks in the two samples.

One potential obstacle is that defining the boundaries of your peaks won't be totally straightforward. For example, a peak in one sample might have 2 overlapping peaks in the other sample, one on each end. A rough solution for this is to divide the genome into bins (maybe around 100-1000 bp, depending on your desired resolution). You can treat a peak as present in a bin if more than half of the peak lies in the bin. That way, bins in one sample can be directly compared to the corresponding bins in the other sample. Obviously, this isn't the only way to do this; other appropriate methods exist too.

**Continuous**

If you want to treat the peak intensities as continuous, you could apply a similar binning method, taking the "score" of a bin to be the average peak intensity at positions within that bin. You could then throw away all bins with no peaks or only low-intensity peaks throughout the genome. Then you could compute the Spearman's correlation for the remaining bins. I'm guessing it will be harder to find a strong correlation for continuous intensities, because of the amount of experimental variability that is inherently present.

If, after following these steps, the Spearman's correlation is still "artificially low" as you suggested, then this is likely a problem with the underlying data, not the overall analysis; maybe your two datasets actually don't agree that well.

# Answer
> 2 votes
It is one of my favorite stories.

Drop a glance to StereoGene software, it for genomic track correlation, it described in a preprint.

You also can run MACS or another peak caller and estimate the correlation of two interval sets using the GenomtriCorr package.

---
Tags: statistics
---


# Question
Title: Are there any rolling hash functions that can hash a DNA sequence and its reverse complement to the same value?
A common bioinformatics task is to decompose a DNA sequence into its constituent k-mers and compute a hash value for each k-mer. Rolling hash functions are an appealing solution for this task, since they can be computed very quickly. A rolling hash does not compute each the hash value from scratch with each k-mer: rather it updates a running hash value using an update strategy and a sliding window over the data.

It's also very useful for many applications to have a k-mer hash to the same value as its reverse complement. Unless the data were generated using a strand-specific sample prep, it's impossible to distinguish a k-mer from its reverse complement, and they should be treated as the same sequence.

Are there any rolling hashes that will map reverse complements to the same value? If not, how would we develop such an algorithm?

**UPDATE**: Ideally the hash function would be able to support k \> 32, which would be lossy unless using something larger than a 64-bit integer.

**ANOTHER UPDATE**: I don't think it's necessary to *store* both the running k-mer and its reverse complement in a single value. If storing two k-mer strings and/or two hash values makes this easier, I'm totally cool with that.

# Answer
A rolling hash function for DNA sequences called ntHash has recently been published in Bioinformatics and the authors dealt with reverse complements:

> Using this table, we can easily compute the hash value for the reverse-complement (as well as the canonical form) of a sequence efficiently, without actually reverse- complementing the input sequence, as follows: ...

EDIT (by @user172818): I will add more details about how ntHash works. The notations used in its paper are somewhat uncommon. The source code is more informative.

Let's first define rotation functions for 64-bit integers:

```
rol(x,k) := x << k | x >> (64-k)
ror(x,k) := x >> k | x << (64-k)

```

We then define a hash function `h()` for each base. In the implementation, the authors are using:

```
h(A) = 0x3c8bfbb395c60474
h(C) = 0x3193c18562a02b4c
h(G) = 0x20323ed082572324
h(T) = 0x295549f54be24456
h(N) = 0

```

The rolling hash function of a forward k-mer `s[i,i+k-1]` is:

```
f(s[i,i+k-1]) := rol(h(s[i]),k-1) ^ rol(h(s[i+1]),k-2) ^ ... ^ h(s[i+k-1])

```

where `^` is the XOR operator. The hash function of its reverse complement is:

```
r(s[i,i+k-1]) := f(~s[i,i+k-1])
               = rol(h(~s[i+k-1]),k-1) ^ rol(h(~s[i+k-2]),k-2) ^ ... ^ h(~s[i])

```

where `~` gives the reverse complement of a DNA sequence. Knowing `f(s[i,i+k-1])` and `r(s[i,i+k-1])`, we can compute their values for the next k-mer:

```
f(s[i+1,i+k]) = rol(f(s[i,i+k-1]),1) ^ rol(h(s[i]),k)  ^ h(s[i+k])
r(s[i+1,i+k]) = ror(r(s[i,i+k-1]),1) ^ ror(h(~s[i]),1) ^ rol(h(~s[i+k]),k-1)

```

This works because `rol`, `ror` and `^` can all be switched in order. Finally, for a k-mer `s`, the hash function considering both strands is the smaller between `f(s)` and `r(s)`:

```
h(s) = min(f(s),r(s))

```

This is a linear algorithm regardless of the k-mer length. It only uses simple arithmetic operations, so should be fairly fast. I have briefly tested its randomness. It seems comparable to murmur. ntHash is probably the best algorithm so far if you want to hash an arbitrarily long k-mer into 64 bits.
> 8 votes

# Answer
Let `f` and `r` be two integers. They always keep the k-mer on the forward and reverse strand, respectively. At a new base `c` in a proper 2-bit encoding, we update the two integers as follows:

```
f = (f<<2|c) & ((1ULL<<2*k) - 1)
r = r>>2 | (3ULL-c)<<2*(k-1)

```

With this updating rule, `f` keeps the forward strand k-mer ending at `c` and `r` is `f`'s reverse complement. The hash of the k-mer can be `min(f,r)`. If integer operations take constant time, computing the hashes of all k-mers of a sequence of length L takes O(L) time.

The above only works if k-mer fits a word, which means `k<=32` on x86. For long k-mers, you can overload bit operators in C++. In case of `32<k<=64`, a simpler solution is to split into lower and higher bits:

```
c1 = c&1, c2 = c>>1
f1 = (f1<<1|c1) & ((1ULL<<k) - 1)
f2 = (f2<<1|c2) & ((1ULL<<k) - 1)
r1 = r1>>1 | (1ULL-c1)<<(k-1)
r2 = r2>>1 | (1ULL-c2)<<(k-1)

```

This is twice as slow as the `k<=32` version. It is possible to implement these lines with SSE2 intrinsics, but that is overkilling.

When you want to hash a long k-mer into fewer bits, there are a few options:

```
xor:      f^r
min:      min(f,r)
min+hash: hash(min(f,r))

```

where `hash()` is a generic randomization hash function, which could be Murmur3, Wang's integer hash function, etc.

The following is a microbenchmark. Here, we squeeze all k-mers into a 32-bit hash table. `collision_rate = 1 - #distinct_hash/#distinct_kmer`. This is not the best metric, but it somehow measures the randomness, which is better than nothing.

```
====================================================
 Algorithm    data  max-k  k   %collision  CPU time
----------------------------------------------------
 xor          chr11  32    31     8.5%       0.9s
 xor+Wang     chr11  32    31     9.6%       1.3s
 min+Wang     chr11  32    31     1.4%       1.3s
 min+Wang     chr11  64    31     1.4%       1.9s
 min+murmur3  chr11  inf   31     1.4%       5.7s
 min+Wang     chr11  64    51     1.5%       2.1s
 min+murmur3  chr11  inf   51     1.5%       6.8s
 min+Wang     chr11  64    63     1.5%       2.1s
 min+murmur3  chr11  inf   63     1.5%       7.5s
 radix sort   chr11  -     -      -        5.6-7.3s
----------------------------------------------------
 min+Wang     hg38   64    63    26.2%        37s
 min+murmur3  hg38   inf   63    26.2%       192s
 radix sort   hg38   -     -      -         ~138s
====================================================

```

Observations:

* Don't use XOR. In fact, XOR hashes all palindromes to `0xffff...`. This is already worrying enough.
* Unless you want to work with very long k-mers, don't use generic string hash functions like FNV and Murmur. Murmur is even slower than radix sorting all hashes.
* For `k<=64`, min+Wang is the best here. It is fast and simple to compute and has randomness nearly as good as murmur.
> 8 votes

# Answer
If your goal is to minimise storage by just having one hash per kmer and its reverse complement, there is a simple solution for non-rolling hashes. For any sequence S, you compute and store the hash of the smaller of S and its reverse complement. A simple lexicographical comparison for "smaller" is enough. In programming terms:

```
hash=computeHash(min(S,rev(S));

```

If your goal is to minimise computing time via rolling hashes AND store only one kmer value, this will be hard. Hashes are normally designed to minimise collisions and you are asking it not only to collide, but also to collide in very, very specific circumstances.

Question is: is the computation of a hash a bottleneck for an application? Maybe, maybe not.

This post on StackExchange seems to imply that, in 2012, MurmurHash was able to compute de-novo a hash for a UUID (similar to a standard kmer in size and complexity) in ~250ns, i.e., ~4 million hashes per second. That's already quite fast.

My guess is that most applications will use a hash to look up things in memory or disk. And here, even with look-ups in memory, real life data will let you run into CPU cache miss speed penalties very quickly and I expect the impact of this to be higher than the speed of the hash computation itself.
> 6 votes

# Answer
At its easiest, you just store the *forward* ($F$) and *reverse* ($R$) hash value.

You update the forward hash value by conventional means, e.g. bit-shifting the base value into its lower bits:

$$ F\_{n + 1} = ((F\_n \ll B \mid x) \mathop\\& M, $$

$B$ is the bit size of the base encoding, $M$ is the word mask for a word of length $W$ bits, and can be omitted if the width of the hash is the width of the machine word; and the reverse hash value by doing the mathematical inverse, e.g. bit-shifting the value of the reverse complement of the base into the upper bits:

$$ R\_{n + 1} = ((R\_n \gg B) \mid (\operatorname{compl}x \ll (W - B))) \mathop\\& M. $$

And then you compute a *combined hash* value by xoring the two hashes, $H = F \oplus R$.
> 4 votes

---
Tags: hashing, k-mer, algorithms
---


# Question
Title: Is it possible to perform MinION sequencing offline?
I vaguely remember, that the original plan of Oxford Nanopore was to provide cheap sequencers (MinION), but charge for base-calling. For that reason the base-calling was performed in the cloud, and the plan was to make it commercial once the technology is established.

Of course, this limits the potential uses of MinION in the field, since huge areas do not have decent internet connection. Also, not all the data can be legally transferred to a third-party company in the clinical studies.

For example for the Ebola paper, they had to create a special version of their software:

> An offline-capable version of MinKNOW, with internet ‘ping’ disabled and online updates disabled was made available to us by Oxford Nanopore Technologies specifically for the project

There are couple of third-party base-callers available today. I am aware of Nanocall and DeepNano, but since they are not official, it can be hard for them to keep-up with the latest versions of sequencers and cells.

1. Is it possible as of today to sequence offline without a special arrangement (like the Ebola one).
2. If not, what's the policy of Oxford Nanopore toward third-party base-callers? Are they going to help them, or to sue them eventually?

# Answer
> 8 votes
Short answer: yes, but you need to get permission (and modified software) from ONT before doing that.

... but that doesn't tell the whole story. This question has the potential to be very confusing, and that's through no fault of the questioner. The issue is that for the MinION, sequencing (or more specifically, generating the raw data in the form of an electrical signal trace) is distinct and separable from base calling. Many other sequencers also have distinct raw data and base-calling phases, but they're not democratised to the degree they are on the MinION.

The "sequencing" part of MinION sequencing is carried out by ONT software, namely MinKNOW. As explained to me during PoreCampAU 2017, when the MinION is initially plugged into a computer it is missing the firmware necessary to carry out the sequencing. The most recent version of this firmware is usually downloaded at the start of a sequencing run by sending a request to ONT servers. In the usual case, you can't do sequencing without being able to access those servers, and you can't do sequencing without ONT knowing about it. However, ONT acknowledge that there are people out there who won't have Internet access when sequencing (e.g. sequencing Ebola in Africa, or metagenomic sequencing in the middle of the ocean), and an email to `<support@nanoporetech.com>` with reasons is likely to result in a quick software fix to the local sequencing problem.

Once the raw signals are acquired, the "base-calling" part of MinION sequencing can be done anywhere. The ONT-maintained basecaller is Albacore, and this will get the first model updates whenever the sequencing technology is changed (which happens a lot). Albacore is a local basecaller which can be obtained from ONT by browsing through their community pages (available to anyone who has a MinION); ONT switched to only allowing people to do basecalling locally in about April 2017, after establishing that using AWS servers was just too expensive. Albacore is open source and free-as-in-beer, but has a restrictive licensing agreement which limits the distribution (and modification) of the program. However, Albacore is not the only available basecaller. ONT provide a FOSS basecaller called nanonet. It's a little bit behind Albacore on technology, but ONT have said that all useful Albacore changes will eventually propagate through to nanonet. There is another non-ONT basecaller that I'm aware of which uses a neural network for basecalling: deepnano. Other basecallers exist, each varying distances away technology-wise, and I expect that more will appear in the future as the technology stabilises and more change-resistant computer scientists get in on the act.

Edit: ONT has just pulled back the curtain on their basecalling software; all the repositories that I've looked at so far (except for the Cliveome) have been released under the Mozilla Public License (free and open source, with some conditions and limitations). Included in that software repository is Scrappie, which is their testing / bleeding-edge version of Albacore.

# Answer
> 5 votes
In short, yes offline basecalling is supported.

Oxford nanopore has made available an official offline basecaller, called albacore. If you have an account to the nanpore community website you can download it from here.

Concerning running it in the field, last time I checked live basecalling was still requiring a lot of computing power, so you might want to run the MinION in the field, but basecall once back a your lab / computing facility.

I'm not aware of any plans to sue the third-party basecallers, I think Oxford Nanopore is pretty positive about them :-)

---
Tags: ngs, minion
---


# Question
Title: How to manipulate protein interaction network from String database in R?
How can I manipulate protein-interaction network graph from the String database using `STRINGdb` bioconductor package and R?

I have downloaded the whole graph for *Homo sapiens* from STRING, which has about 20.000 proteins.

1. How do I read the file using that package?
2. How do I filter things I don't need? Supposing that I want to keep tumor data, as an example.

# Answer
> 4 votes
I think the easiest way is to download the graph using `STRINGdb`.

```
library(STRINGdb)
string_db <- STRINGdb$new(version="10", species=9606,
                          score_threshold=400, input_directory="" )
full.graph <- string_db$get_graph()

```

Now you can use `igraph`, to manipulate the graph. Let's assume you want to take 200 proteins with the highest degree, i.e. number of edges they have.

```
library(igraph)

# see how many proteins do you have
vcount(full.graph)

# find top 200 proteins with the highest degree
top.degree.verticies <- names(tail(sort(degree(full.graph)), 200))

# extract the relevant subgraph
top.subgraph <- induced_subgraph(full.graph, top.degree.verticies)

# count the number of proteins in it
vcount(top.subgraph)

```

**How to get disease specific genes?**

There's no GO annotation for cancer or Alzheimer's disease. It is out of scope of the GO consortium.

What you can do, you can either take KEGG Pathways annotation, or manually select list of relevant GO-terms. Or acquire the list from one of the papers. For example annotation term `05200` corresponds to the cancer KEGG pathway. You can easily retrieve proteins associated with the annotation:

```
cancer.pathway.proteins <-
    string_db$get_term_proteins('05200')$STRING_id

```

And then perform subgraphing as described above.

Alternatively you can try to get an enrichment score for an every gene given it's neighbors (the way enrichment is shown on the string-db website). Then you can keep only those having top enrichment scores. Probably `get_ppi_enrichment_full` or `get_ppi_enrichment` functions will help you to do that.

---
Tags: r, bioconductor, database, networks
---


# Question
Title: Changing the record id in a FASTA file using BioPython
I have the following FASTA file, `original.fasta`:

```
>foo
GCTCACACATAGTTGATGCAGATGTTGAATTCACTATGAGGTGGGAGGATGTAGGGCCA

```

I need to change the record id from `foo` to `bar`, so I wrote the following code:

```
from Bio import SeqIO

original_file = r"path\to\original.fasta"
corrected_file = r"path\to\corrected.fasta"

with open(original_file) as original, open(corrected_file, 'w') as corrected:
    records = SeqIO.parse(original_file, 'fasta')
    for record in records:
        print record.id             # prints 'foo'
        if record.id == 'foo':
            record.id = 'bar'
        print record.id             # prints 'bar' as expected
        SeqIO.write(record, corrected, 'fasta')

```

We printed the record id before and after the change, and get the expected result. We can even doublecheck by reading in the corrected file again with BioPython and printing out the record id:

```
with open(corrected_file) as corrected:
    for record in SeqIO.parse(corrected, 'fasta'):
        print record.id                  # prints 'bar', as expected

```

However, if we open the corrected file in a text editor, we see that the record id is not `bar` but `bar foo`:

```
>bar foo
GCTCACACATAGTTGATGCAGATGTTGAATTCACTATGAGGTGGGAGGATGTAGGGCCA

```

We can confirm that this is what is written to the file if we read the file using plain Python:

```
with open(corrected_file) as corrected:
    print corrected.readlines()[0][1:] # prints 'bar foo'

```

Is this a bug in BioPython? And if not, what did I do wrong and how do I change the record id in a FASTA file using BioPython?

# Answer
If I use `SeqIO.parse(filehandle, 'fasta')` to parse a FASTA file, then it will return a `SeqRecord` object where the `id` and `name` are the first word (everything before the first whitespace) of the line beginning with `>` and the `description` is the complete line (all not including the initial `>`). (This behaviour can overruled by providing a custom `title2ids` function).

So, if I for example have a `original.fasta` file like:

```
>Lorem|ipsum>dolor sit amet
GCTCACACATAGTTGATGCAGATGTTGAATTCACTATGAGGTGGGAGGATGTAGGGCCA

```

Then I will get:

```
>>> from Bio import SeqIO
>>> path = r'C:\path\to\original.fasta'
>>> records = SeqIO.parse(open(path), 'fasta')
>>> record = next(records)
>>> record.id
'Lorem|ipsum>dolor'
>>> record.name
'Lorem|ipsum>dolor'
>>> record.description
'Lorem|ipsum>dolor sit amet'

```

So, what happens in my case when you do

```
    if record.id == 'foo':
        record.id = 'bar'

```

is that the `record.id` is successfully changed from `foo` to `bar`, but that the `record.description` is not changed and stays `foo`. That's why, when the FASTA file is printed out, I see the described behaviour of

```
>bar foo
GCTCACACATAGTTGATGCAGATGTTGAATTCACTATGAGGTGGGAGGATGTAGGGCCA

```

So, the solution to my problem is to both change the `id` AND the `description`:

```
from Bio import SeqIO

original_file = r"path\to\original.fasta"
corrected_file = r"path\to\corrected.fasta"

with open(original_file) as original, open(corrected_file, 'w') as corrected:
    records = SeqIO.parse(original_file, 'fasta')
    for record in records:
        if record.id == 'foo':
            record.id = 'bar'
            record.description = 'bar'
        SeqIO.write(record, corrected, 'fasta')

with open(corrected_file) as corrected:
    records = SeqIO.parse(corrected, 'fasta')
    for record in records:
        print record.id             # prints bar
        print record.name           # prints bar
        print record.description    # prints bar

```
> 11 votes

---
Tags: fasta, python
---


# Question
Title: How to transfer gff annotations in genome with extensive duplications?
Microbial genomes can contain extensive duplications. Often we'd like to transfer annotations from an annotated species to one that is newly sequenced.

Existing tools (e.g. RATT, LiftOver, Kraken) either make specific assumptions about how closely related the species are or fail to transfer when multiple matches are found in the new genome, especially if the sequences are highly similar.

Specifically, I have a synthetic biology application where genes can duplicate extensively. They are identical in sequence but duplicated many times and be relocated (i.e., not just adjacent to each other). None of the above mentioned tools are able to transfer coordinates of annotations to genomes with multiple copies of features.

Are there any pre-existing tools or software that transfer annotations in this scenario? Ideas for ways to do this robustly?

# Answer
There is one very simplistic way I use which *might* work for what you are doing, it is similar to what terdon proposed.

Take a de-novo microbial genome annotation tool (I have my own, but you could use/modify prokka). Tools like these often first predict gene boundaries (with other tools like prodigal or glimmer) and then try to assign a function to found genes. This function assignment is often done with BLAST and other tools ... and that is where you can go in and modify to do what you need.

I use a "knowledge" protein database of genes I want to have very strictly annotated as a first line of annotation (e.g. in your case: the annotated genomes). For that I loop through very strict identity/similarity parameters which get gradually relaxed.

E.g.: Loop 0: only transfer annotations at 100% DNA identity, same length. Loop 1: only transfer annotations at 100% similarity, same length. Loop 2: only transfer annotations at 99% similarity, length +/- 1%. ... Loop n: only transfer annotations at 100-(n-1)% similarity, length +/- (n-1)%.

In each loop, obviously only annotate what has not been annotated in previous loops.

After that, use "normal" annotation pipeline of the tool to annotate the rest.
> 5 votes

# Answer
I think you will have to first identify the regions homologous to the ones defined in your GFF and then transfer the annotations. Of course, the assumption there is that the homolog will also have the same annotation which is often not true. However, I don't see how you can do it in any other way since you cannot use genomic coordinates (and you would still be making the same assumption even if you could, anyway) when the genomes are so different.

For a very simplistic approach (which might be enough if, as you say, your sequences are almost identical), you can do something like:

1. Collect the sequences of interest from your already annotated species.
2. Use a tool like `genewise` or `exonerate` to map these into the target genome. Both tools can return gff-formatted output and both can find multiple hits in the target genome. For what you want, I would suggest using a very high threshold of sequence similarity and query coverage (where the target sequence found covers all or most of the query sequence used).

Since these are microbial genomes and therefore splicing isn't a problem, you could do the same thing with even a simple BLASTn or tBLASTn if you start from protein sequences.
3. At this point, you should have a list of homologs (some of which will be orthologs and others paralogs) and you can transfer the annotations of the query sequence over to the target.

Again, I stress that this is making a whopping huge assumption: that homologous sequences have the same function and can automatically be annotated as whatever you had in the query genome. This is going to be true for many cases but it will also be false for others. Especially if you are looking at paralogs (genes whose duplication occurred after the speciation event and are therefore likely to have diverged in function).

However, as I said before, this problem would be exactly the same even if you did manage to transfer annotations just by identifying the syntenic regions of the genomes1, so there's not much difference there.

---

1 As I said in the comments, I don't see how this could be possible. By definition, if you have extensive duplications, the genomic coordinates will be completely different and it is impossible to map from one genome into the other.
> 3 votes

---
Tags: sequence-annotation, ngs
---


# Question
Title: Rapid metagenomics classifiers on long read data
I recently used the minION (Nanopore, 9.4 flow cell, RAD001 kit) to generate a metagenome out of environmental samples. Passed reads weren't brilliant (196, average 1,594bp lenght), but working with centrifuge the classification outputs turned out to have quite low hitLength to queryLength scores (average 2%, max 14%). This plus score values don't give me a lot of confidence towards the results I got.

Has anyone else used *centrifuge* and experienced the same?

# Answer
You should also keep in mind that single read accuracy of ONT still is a bit lacking, a 2D or 1D^2 accuracy of 95% still means that on average there's one error every 20 bp, and due to the nature of the data some stretches may be more junky than others. Maybe centrifuge doesn't like that.

Getting a feel of what the data can tell you (and what not) is pretty important.

You may want to try other classification tools to get an idea how well centrifuge works, Kraken comes to mind. But also: what about BLASTing some of these reads at NCBI and simply look at the results to get a first impression?
> 1 votes

---
Tags: nanopore, long-reads, minion, metagenome
---


# Question
Title: Which tools can detect chimeric RNA (fusion genes) from WGS or RNA-Seq data?
Given WGS data or RNA-seq data, which tools can I use to detect gene fusions?

# Answer
Most of these use RNA-seq data, some use WGS data, and some use both. They are listed alphabetically. I will add to the list when I discover more.

1. Barnacle: http://bmcgenomics.biomedcentral.com/articles/10.1186/1471-2164-14-550
2. Bellerophontes: http://bioinformatics.oxfordjournals.org/content/28/16/2114.long
3. BreakDancer: http://www.nature.com/nmeth/journal/v6/n9/abs/nmeth.1363.html
4. BreakFusion: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3389765/
5. BreakPointer: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3561864/
6. ChimeraScan: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3187648/
7. Comrad: http://bioinformatics.oxfordjournals.org/content/27/11/1481.long
8. CRAC: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC4053775/
9. deFuse: http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1001138
10. Dissect: http://bioinformatics.oxfordjournals.org/content/28/12/i179.abstract
11. EBARDenovo: http://bioinformatics.oxfordjournals.org/content/early/2013/03/01/bioinformatics.btt092
12. EricScript: http://bioinformatics.oxfordjournals.org/content/28/24/3232
13. FusionAnalyser: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3439881/
14. FusionCatcher: http://biorxiv.org/content/early/2014/11/19/011650.full-text.pdf+html
15. FusionFinder: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3384600/
16. FusionHunter: http://bioinformatics.oxfordjournals.org/content/27/12/1708.long
17. FusionMap: http://bioinformatics.oxfordjournals.org/content/27/14/1922
18. FusionQ: http://www.biomedcentral.com/1471-2105/14/193
19. FusionSeq: http://www.genomebiology.com/2010/11/10/R104
20. IDP-fusion: http://nar.oxfordjournals.org/content/early/2015/06/03/nar.gkv562.full
21. iFUSE: http://bioinformatics.oxfordjournals.org/content/29/13/1700.long
22. InFusion: https://bitbucket.org/kokonech/infusion/wiki/Home
23. INTEGRATE: http://www.ncbi.nlm.nih.gov/pubmed/26556708
24. JAFFA: http://www.genomemedicine.com/content/7/1/43
25. LifeScope: http://www.thermofisher.com/no/en/home/life-science/sequencing/next-generation-sequencing/solid-next-generation-sequencing/solid-next-generation-sequencing-data-analysis-solutions/lifescope-data-analysis-solid-next-generation-sequencing/lifescope-genomic-analysis-software-solid-next-generation-sequencing.html
26. MapSplice: http://www.ncbi.nlm.nih.gov/pubmed/20802226
27. MOJO https://github.com/cband/MOJO
28. nFuse: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3483554/
29. Pegasus: http://www.biomedcentral.com/1752-0509/8/97
30. PRADA: http://www.ncbi.nlm.nih.gov/pubmed/24695405
31. ShortFuse: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3072550/
32. SnowShoes-FTD: http://nar.oxfordjournals.org/content/39/15/e100
33. SOAPFuse: http://www.genomebiology.com/2013/14/2/R12
34. SOAPFusion: http://www.ncbi.nlm.nih.gov/pubmed/24123671
35. STAR: http://bioinformatics.oxfordjournals.org/content/29/1/15
36. STAR-Fusion: https://github.com/STAR-Fusion/STAR-Fusion/wiki
37. TopHat-Fusion: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3245612/
38. TRUP: http://www.genomebiology.com/2015/16/1/7
39. ViralFusionSeq: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3582262/

**Other useful programs:**

Chimeraviz (visualization tools for gene fusions): https://www.ncbi.nlm.nih.gov/pubmed/28525538 (disclaimer: I created this)

Chimera: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC4253834/
OncoFuse: http://bioinformatics.oxfordjournals.org/content/29/20/2539.long
FuMa: http://bioinformatics.oxfordjournals.org/content/early/2015/12/09/bioinformatics.btv721.abstract
> 9 votes

# Answer
1- Chimerascan + Star

2- Tophat Fusion

Two pipelines I have used extensively and recommend. The others listed by @L42 I can't speak to but Star Fusion sounds promising. Star is fast.
> 0 votes

---
Tags: rna-seq, wgs
---


# Question
Title: Gene set enrichment analysis on differential phosphorylation sites
I have:

* A list of differentially phosphorylated sites in a knockout condition. Some genes contain as many as 70 possible phosphorylation sites; others contain only one.
* A list of genes belonging to a specific gene set annotation.

**How can I test the differentially phosphorylated proteins for enrichment of this annotation?**

A few ideas I’ve considered:

1. Ignore the number of phosphorylation events detected within a gene and simply count the gene as differentially phosphorylated if it contains at least a single site that is differentially phosphorylated. Compare the enrichment score for this selected set to the enrichment score for the set genes containing at least one site that is not differentially phosphorylated. The problem here is that genes with a large number of phosphorylation sites have almost no influence in the enrichment score, since it’s almost certain that they’ll have at least one site that is differentially phosphorylated and at least one that is not.
2. Mark each candidate phosphorylation site as either “in the set” or “not in the set” based on the protein in which it’s found. Then perform the enrichment analysis using the set of annotated phosphorylation sites instead of using the traditional enrichment analysis performed at the gene level. The potential problem with this approach is that it may place too much influence on genes with many potential phosphorylation sites.
3. Aggregate all candidate phosphorylation sites within a gene and use some numerical threshold to determine whether the gene is differentially phosphorylated or not. (There are various ways that this could be done.) Then perform enrichment analysis using the resulting set of differentially phosphorylated genes. A possible problem here is that some of the phosphorylation sites may be more functionally important than others, so it’s not clear how to weight the relative importance of individual phosphorylation sites within a gene.

I realize the goal here isn’t well-defined mathematically; I’m mainly curious what approach makes the most sense given the biological context. Of the above approaches, I’m currently leaning towards approach 2 because it’s straightforward to implement and it at least *attempts* to account for the variable number of phosphorylation events within a gene.

NOTE: I also have normalized phosphoprotein and protein abundances for all of these sites, obtained from mass spec. So if the solution requires an alternative method of computing differential phosphorylation, that’s fine.

# Answer
> 1 votes
Your dataset looks perfect for the SUMSTAT enrichment test.

1. You need to come up with a statistic representing your gene. The simplest ideas here are the number of sites or, probably better, a proportion of phosphorylated sites.
2. Now you can compute a statistic for every gene set, for example just sum (SUMSTAT). You can also have average, sum of squares or something else.
3. Geta a null distribution of your gene set statistics by permutations. You can either randomize statistics for genes, or just randomly assign phosphorylated sites across the genome keeping the total number of sites.

Now you can compute p-value just by comparing your value with the null-distribution.

You have to be cautious of two things:

1. You are testing multiple gene sets, i.e. you perform a statistical test for every gene set. Therefore correct your significance for multiple testing. I suggest using FDR for this.
2. Beware of potential biases when doing permutations. The obvious one is gene length, i.e. longer genes have higher chance of getting at least one site. But there can be other ones like GC-content, chromosomes, etc. You can overcome this by using more realistic permutations, or by controlling for the potential correlations. You can get a few ideas on fighting biases from this paper (sorry for the self-advertisement).

---
Tags: gse, differential-expression
---


# Question
Title: How do I validate a single sample ArrayCGH result?
We have arrayCGH (aCGH) results for one sample. There is a 0.5 Mb terminal duplication on chromosome 19 (62995490-63407936, according to NCBI36/hg18). The duplication is rare: a literature review suggests there are only 3-4 samples with clinical information.

What are the steps to validate the results? How do we ascertain that this duplication is the cause of the clinical symptoms?

I have some ideas:

* aCGH the parents. Not sure how this would help.
* whole genome exome sequencing. Worried this might make it more difficult to pinpoint genetic cause.
* whole genome sequencing?
* other ideas?

**Note:** I am new to aCGH and high-throughput sequencing, any advice is welcome.

# Answer
Why not just good old qPCR? That's (A) quick, (B) cheap and (C) easy to analyze. If you care about the exact location of the break point (I'm guessing from the context that you don't), then targeted sequencing with a custom capture kit would work.

Regarding validating the biological relevance of this, there are multiple (parallel) routes one can take. Firstly, screen unaffected family members for this same alteration. If you find this in unaffected individuals then it's obviously not the causative alteration. Ideally, one would also do either a cell-line experiment or a mouse (or other model organism) experiment to see if you can reconstitute at least some component of the clinical phenotype. This may not always be possible, of course.
> 2 votes

---
Tags: ngs, array-cgh, exome
---


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

Which is the one with more genes associated to pathways ? Which is more completely annotated ? Is there any paper comparing them ?

# Answer
One big downside of KEGG is the licensing issue. One big advantage of Reacome are various crosslinks to other databases and data.

**ad 1,** This depends on which pathway, they are both primary databases. Sometimes other databases that for instance combine data of primary databases have better annotation of pathways (there is an example in the review paper bellow)

**ad 3,** There is very extensive relatively new (2015) review on this topic focused on human pathways: Comparison of human cell signaling pathway databases—evolution, drawbacks and challenges. However I could not find there which one is more complete ...
> 4 votes

# Answer
For you the main point would be whether an enrichment analysis is going to give you an *informative* answer. That's what will make a particular database better. And there's all sorts of subjective decisions made in their construction as to what what a pathway is, what to include, where to draw boundaries, etc. so different databases will give different answers and it will not be clear which is more correct.

So browse over the references / pointers people have given, select a service and use it and no others. Don't jump around databases until you get an answer you like, that's just fishing.
> 4 votes

# Answer
One big advantage of Reactome, in my opinion, is its visualization using the web interface.

Many pathways (in Reactome and KEGG) consist of genes / proteins that are up- and down-regulated through the respective pathway. If you do a simple overrepresentation analysis this is not taken into consideration. Therefore, you might end up seeing a pathway as "overexpressed" although only the down-regulated genes were observed more frequently.

In Reactome, you can zoom in on the different pathways very conveniently and then pick up these inconsistencies. I haven't really found a public database and tool that can take these different regulations into consideration. Therefore, you'll probably always will need some manual investigation of your data. In my opinion, this is easier with Reactome than with KEGG.
> 1 votes

---
Tags: database
---


# Question
Title: Downloading a reference Genome for Bowtie2
How do I download a reference genome that I can use with bowtie2? Specifically HG19. On UCSC there are a lot of file options.

# Answer
> 9 votes
tl;dr: Just use the either the downloads on the Bowtie2 homepage or the Illumina iGenomes. Or just uncompress and concatenate the FASTA files found on UCSC goldenpath and then build the index.

A bit longer answer:

There are two components to "genome for a read mapper" such as Bowtie or BWA.

First, you need to choose the actual sequence (genome release such as GRCh37/hg19 or GRCh38/hg38). There are patch releases such as GRCh37.p3 where some bases might be exchanged and depending on the release, some "unmapped" loci contigs might be added, but generally GRCh37.p1 is roughly the same as GRCh37.p2, for example. Usually, people have agreed on some specific patch version for each read and use this for read mapping.

Generally, there is the UCSC flavour hg19/hg38 etc. and the NCBI/GRC flavour GRCh37, GRCh38 etc. (similar with mouse). UCSC has no versioning besides the genome release and (to the best of my knowledge) does not update the genome sequence after releasing a hg19 FASTA file.

Second, you have to build the index files for each genome. Depending on the read mapper you use, you might or might not need the original FASTA files for the alignment. For Bowtie and Bowtie 2, you don't need the original FASTA files after building the index as Bowtie 1/2 can reconstruct the sequence "on the fly" from the index files.

HTH

# Answer
> 11 votes
It’s a matter of preference I guess but I recommend the **Ensembl** builds. Decide whether you want the toplevel or primary assembly, and whether you want soft-masked, repeat-masked or unmasked files. The naming schema is very straightforward; the combinations are described in the `README` file, and all files reside in one directory.

For example, if you want the unmasked primary assembly, the file to download would be `Homo_sapiens.GRCh37.75.dna.primary_assembly.fa.gz`.

As for **GoldenPath/UCSC**, there’s no need to download and concatenate separate chromosomes (contrary to what the other answer said); you can download the whole (toplevel) reference from the `bigZips` directory; from the `README`:

> This directory contains the Feb. 2009 assembly of the human genome (hg19, GRCh37 Genome Reference Consortium Human Reference 37 (GCA\_000001405.1)), as well as repeat annotations and GenBank sequences.

There are essentially three options here:

1. `chromFa.tar.gz`, which contains the whole genome in one chromosome per file;
2. `chromFaMasked.tar.gz`, the same with repeats masked by `N`;
3. `hg19.2bit`, which is the whole genome in one file, but needs to be extracted using the utility program `twoBitToFa`, which needs to be downloaded separately.

In any case, I always download the reference and build my own index for mapping, since this allows me more control; not everybody might need this much control, but then building the index once is fairly fast anyway.

---
Tags: human-genome, reference-genome, genome-sequencing
---


# Question
Title: What's the difference between VCF spec versions 4.1 and 4.2?
What are the key differences between VCF versions 4.1 and 4.2? It looks like v4.3 contains a changelog (specs available here) but earlier specifications do not.

This biostar post points out one difference: the introduction of `Number=R` for fields with one value per allele including REF — can anyone enumerate the other changes between these two versions?

# Answer
> 16 votes
This is easy to check, you can download both specs in .tex format and do `diff`.

Changes to the v4.2 compared to v4.1:

1. Information field format: adding source and version as recommended fields.
2. INFO field can have one value for each possible allele (code `R`).
3. For all of the ##INFO, ##FORMAT, ##FILTER, and ##ALT metainformation, extra fields can be included after the default fields.
4. Alternate base (ALT) can include `*`: missing due to a upstream deletion.
5. Quality scores, a sentence removed: *High QUAL scores indicate high confidence calls. Although traditionally people use integer phred scores, this field is permitted to be a floating point to enable higher resolution for low confidence calls if desired.*
6. Examples changed a bit.

---
Tags: vcf, htslib, file-formats
---


# Question
Title: Structural variant calling for low-coverage PacBio data
PacBio is selling ~10x PacBio SEQUEL long reads as an upgrade to Illumina data for SV discovery.

In a clinical setting, the main requirements are proper sensitivity and specificity but also the processing of cohorts, at least families. This requires a genotyping step, such that it can be identified whether a given variant is shared by two or more individuals or whether it is not.

What are the tools of the trade for this task?

As having 50-60x PacBio reads is not an option from a economic point of view, one has to make do with 10x coverage.

# Answer
> 2 votes
There is an evaluation of PB Honey and Sniffles algorithms for low coverage PacBio datasets in this preprint and another evaluation is shown on this poster. Both reports agree that optimal is (surprisingly) combination of PB Honey and Sniffles.

Author of Sniffles have benchmarked Sniffles against PB Honey, where he shown that Sniffles performs significantly better. Take a look on this presentation (slide 15).

Another option is SMRT-SV, but I am not aware of any benchmarking.

# Answer
> 1 votes
I'm aware of the following (very few and suboptimal) options:

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

---
Tags: structural-variation, long-reads, pacbio
---


# Question
Title: Imputing missing genotypes from separate genotyping panels
What is the current standard for imputing missing genotypes between two genotyping panels? I have two populations genotyped using two different panels (A & B), and I would like to impute all the genotypes in population B for those positions on used in panel A.

I've read the examples for impute2, and I think the closest thing to what I am looking for is this example, "Imputation with one unphased reference panel".

Simply put, I want to provide a list of SNPs, some variant file for population B, and haplotype information from 1,000 Genomes and get imputed genotypes for each SNP in the list. Is impute2 the state of the art for this?

# Answer
Given that you mention wanting to use 1000 Genomes as a reference panel for imputing genotypes into your two SNP chip panels, I am going to assume that you are working with human data.

In that case there are several options you can go with:

* If your two panels are of European descent, then you are probably best off using the HRC reference panel together with a fast genotype imputation tool such as Beagle 4.1 to impute genotypes in each of your two SNP chip panels separately.
* If your panels are not of European descent, then you will likely want to use the 1000 Genomes phase 3 reference panel with Beagle 4.1, Impute2, or Minimac3.

In either case, there are two phasing services available that will do much of the heavy lifting for you ,.

The second Wellcome Trust Case-Control Consortium paper performed a cross-imputation analysis as you describe. I don't see many studies using multiple SNP chip panels. You will need to take care in your analysis that you are not hit by batch effects from using two different SNP chip panels.

Also, none of these methods will work if the region you are imputing into has too few variants. I'm not sure what the minimum number of variants is, but if you are using a whole genome genotyping panel of at least 500k SNPs, then you should be ok if you impute a whole chromosome at a time.
> 4 votes

---
Tags: imputation, impute2
---


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

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

but what about non-human/mouse ones?

* fruit fly
* zebrafish
* E. coli
* any other idea?

# Answer
I will answer one of the points – E.coli.

**TL;DR** Bacteria, and in particular E.coli, are highly variable and there is usually no single best assembly. Large scale WGS studies should come with multiple assemblies for individual monophyletic clusters.

**Long answer:** Whereas a single reference sequence can make sense for human or mice (to a certain extent), bacteria are not the case. The variability between them even within a single species can be enormous and various lineages can contain very different genes. This is the reason why biologists usually distinguish a *core genome* consisting of the genes shared among a vast majority (typically 95%) of individuals within a given phylogenetic clade, (i.e., of those genes which are essential for this clade – imagine a CPU and RAM in a computer), and an *accessory genome* consisting of the other, unnecessary, genes – imagine a joystick in our analogy). The size of the accessory genome can vary a lot for different bacteria. For instance, Escherichia coli is said to have an *open* accessory genome (i.e., much bigger than the core genome; imagine a Raspberry-PI and all the extension kits) whereas Chlamydia trachomatis has a *closed* accessory genome (i.e., small compared to the core genome; imagine an ATM).

For species having a small accessory genome, a single reference might be sufficient. However, if a core genome is large, we can use either pan-genomes, or multiple references. The word pan-genome tends to be used for different things: the computational biologists usually use it as a representation of all genomic content in a certain phylogenetic clade (which can be mathematically modeled using finite automata, hence as graphs) whereas some biologists define it as a concatenation of all (core and accessory) genes from the clade (without considering SNPs, etc.). Unfortunately, methods for graph references (especially for read mapping) are still not sufficiently developed, therefore, researchers still have to resort to reference sequence-based methods and use either concatenated genes or multiple reference sequences.

When particular bacterial species are studied using massive whole-genome sequencing of many isolates, researchers usually cluster the isolates to monophyletic sequence clusters and infer a special reference sequence for each of these clusters (see, e.g., this example for the pneumococcus). All of the inferred references should be then published together with the paper or in a separate data paper (see, e.g., the sequences from the previous example).

**Disclaimer:** I am not a biologist. If anything in this post is too simplified or incorrect, please, leave me a comment and I will try to fix it. I tried to make this text accessible for my younger self.
> 5 votes

# Answer
If the programs you are using allow for it, take the most recent available genome, as it will be most likely to have the fewest errors.

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

# Answer
The "best" assembly depends on the species. For more common species, Ensembl/UCSC/NCBI has their own version that corresponds to the most popular assembly. Usually the actual genetic sequence is identical across all of them, but the chromosome names and gene annotations will vary. If those three sources agree for a species, there is a "standard" reference.

For fruit fly specifically, FlyBase Consortium/Berkeley Drosophila Genome Project is the gold standard. Other species have their own organizations.
> 0 votes

---
Tags: read-mapping, reference-genome
---


# Question
Title: Missing genes and normalisation of RSEM output using EBSeq
Without going into too much background, I just joined up with a lab as a bioinformatics intern while I'm completing my masters degree in the field. The lab has data from an RNA-seq they outsourced, but the only problem is that the only data they have is preprocessed from the company that did the sequencing: filtering the reads, aligning them, and putting the aligned reads through RSEM. I currently have output from RSEM for each of the four samples consisting of: gene id, transcript id(s), length, expected count, and FPKM. I am attempting to get the FASTQ files from the sequencing, but for now, this is what I have, and I'm trying to get something out of it if possible.

I found this article that talks about how expected read counts can be better than raw read counts when analyzing differential expression using EBSeq; it's just one guy's opinion, and it's from 2014, so it may be wrong or outdated, but I thought I'd give it a try since I have the expected counts.

However, I have just a couple of questions about running EBSeq that I can't find the answers to:

1: In the output RSEM files I have, not all genes are represented in each, about 80% of them are, but for the ones that aren't, should I remove them before analysis with EBSeq? It runs when I do, but I'm not sure if it is correct.

2: How do I know which normalization factor to use when running EBSeq? This is more of a conceptual question rather than a technical question.

Thanks!

# Answer
> 6 votes
Yes, that blog post *does* represent just one guy's opinion (hi!) and it does date *all the way back to 2014*, which is, like, decades in genomics years. :-) By the way, there is quite a bit of literature discussing the improvements that expected read counts derived from an Expectation Maximization algorithm provide over raw read counts. I'd suggest reading the RSEM papers for a start\[\]\[\].

But your main question is about the mechanics of running RSEM and EBSeq. First, RSEM was written explicitly to be compatible with EBSeq, so I'd be very surprised if it does not work correctly out-of-the-box. Second, EBSeq's `MedianNorm` function worked very well in my experience for normalizing the library counts. Along those lines, the blog you mentioned above has another post that you may find useful.

But all joking aside, these tools *are indeed* dated. Alignment-free RNA-Seq tools provide orders-of-magnitude improvements in runtime over the older alignment-based alternatives, with comparable accuracy. Sailfish was the first in a growing list of tools that now includes Salmon and Kallisto. When starting a new analysis from scratch (i.e. if you ever get the original FASTQ files), there's really no good reason not to estimate expression using these much faster tools, followed by a differential expression analysis with DESeq2, edgeR, or sleuth.

---

1**Li B, Ruotti V, Stewart RM, Thomson JA, Dewey CN** (2010) RNA-Seq gene expression estimation with read mapping uncertainty. *Bioinformatics*, 26(4):493–500, doi:10.1093/bioinformatics/btp692.

2**Li B, Dewey C** (2011) RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. *BMC Bioinformatics*, 12:323, doi:10.1186/1471-2105-12-323.

# Answer
> 3 votes
1. Include all genes/transcripts in your analysis.

A transcript that is not detected could be undetected through sampling error (i.e. the sequencer / library prep just happened to miss that transcript), or it could be because the transcript isn't generated in a particular sample. It's not uncommon for genes to be switched off in response to different biological factors, so zero-count genes shouldn't be ignored. I can't speak from experience with EBSeq, but as long as the analysis package treats a zero count as "unobserved" rather than "absent" (and makes relevant corrections), it's a good idea to keep them in.

---
Tags: rna-seq, differential-expression, ebseq, rsem
---


title: Using shells other than bash

text:
# Question
Title: Using shells other than bash
As someone who's beginning to delve into bioinformatics, I'm noticing that like biology there are industry standards here, similar to Illumina in genomics and bowtie for alignment, many people use bash as shell.

Is using a shell besides bash going to cause issues for me?

# Answer
> 18 votes
Bioinformatics tools written in shell and other shell scripts generally specify the shell they want to use (via `#!/bin/sh` or e.g. `#!/bin/bash` if it matters), so won't be affected by your choice of user shell.

If you are writing significant shell scripts yourself, there are reasons to do it in a Bourne-style shell. See *Csh Programming Considered Harmful* and other essays/polemics.

A Bourne-style shell is pretty much the industry standard, and if you choose a substantially different shell you'll have to translate some of the documentation of your bioinformatics tools. It's not uncommon to have things like

> Set some variables pointing at reference data and add the script to your PATH to run it:
>
> ```
> export FOO_REF=/path/to/stuff
> export PATH=/path/to/foo-x.y:$PATH
> foo blah blah
>
> ```

These will typically be shown in Bourne-shell syntax. By using a different shell you have to translate the `export` commands to your local syntax, and especially PATH munging is somewhat shell-dependent.

If you're experienced in Unix, this will be only a minor niggle. If you're a beginner, IMHO this will add a non-negligible amount of friction on top of all the other things you're learning.

# Answer
> 11 votes
**TL;DR**

SH adheres to an official industry standard, but it is not suitable for scientific computing. Bash is considered an informal standard (e.g., by Google). Bash 3 is preferable in most of situations in the world of bioinformatics.

**Long answer**

As already described in other answers, SH (`/bin/sh`, plain Bourne shell, the original UNIX shell) should fully adhere to POSIX which is a real industry standard. However, SH is too limited for scientific computing since many key features were incorporated later in SH successors, especially in Bash (`/bin/bash`, Bourne Again Shell): `set -o pipefail`, `[[ ... ]]`, or process substitutions `< ()` to name at least few.

In practice, it is much more difficult to write "safe" scripts in pure SH and only shell experts are usually capable to prevent unexpected behavior. For example, it may be hard to ensure that no command in a pipeline failed in the middle of computation. For Bash, various easy-to-follow defensive-programming recommendations have been developed and they should prevent all these problems. From this reason, many computer scientists, software engineers and companies use Bash as a kind of a standard. For instance, Google internal policy allows only Bash for writing any shell scripts.

Even though we cannot expect that Bash is present on completely every Unix machine (e.g., on mobile devices as @terdon pointed out), a vast majority of \*nix machines used for scientific computation should have it. We also should be aware of the fact that Bash can be slower than SH and that it has recently suffered from major security issues. Moreover, various Bash versions exist and scripts working on modern Linux machines with Bash 4 might not work on OS X, which is still based on Bash 3.

To sum up, Bash 3 is probably be the most reasonable choice for scientific computing.

**Edit**

I addressed the comments from @terdon and @John Marshall. In particular, I added an explanation why Bash is more suitable for scientific computing than SH (in my opinion).

# Answer
> 7 votes
The Open Group Base Specifications Issue 7 IEEE Std 1003.1™-2008, 2016 Edition, or "The POSIX Standard" for short, is the standard that defines the interfaces and utilities provided by a Unix system. Among these is the command line shell language and tools (see "Shell & Utilities" in the main index on the page linked above).

As far as I know, there is no shell that implement *exactly* what's specified by the standard, but both `bash` and `ksh93` does a pretty good job of adhering to the standard along with their own, sometimes conflicting, extensions. The `ksh93` shell in particular has had a big impact on the past development of the POSIX shell specification, but future POSIX specifications *may* borrow more from `bash` due to its wide use on Linux.

The `bash` shell is pretty much ubiquitous on Linux systems, and may be installed on all other Unices too. `ksh93` is also available for most Unices but is usually not installed by default on Linux. `ksh93` is available by default on at least macOS (as `ksh`) and Solaris.

If you are concerned about portability when writing a shell script (which is IMHO a good thing to be concerned about), you should make sure that you use only the POSIX utilities and their POSIX command line flags, as well as only use POSIX shell syntax. You should then ensure that you script is executed by `/bin/sh` which is supposed to be a shell that understands the POSIX specification. `/bin/sh` is often implemented by `bash` running in "POSIX mode", but it may also be `dash`, `ash` or `pdksh` (or something else) depending on what Unix you are using.

For a Linux user, the most difficult bit in writing a portable script is often not the shell per se, but the multitude of non-standard command line flags provided by the GNU implementation of the many shell utilities. The GNU coreutils (basic shell utilities) may, like `bash`, be installed on all Unices though.

Also note that `bash`, when running in POSIX mode (either when invoked as `/bin/sh` or with its `--posix` command line flag), is not *strict* about its POSIX conformance and may accept some syntax extensions to the POSIX standard.

# Answer
> 5 votes
I would not say bash as a "standard", but it is indeed likely to be the most widely used unix shell and available by default on most modern unix/linux distros. There are a few other more convenient shells like zsh that are broadly compatible with `/bin/sh`, but they are not as widely available. There is also C-shell and in particular its open-source implementation tcsh. C-shell is quite different from bash. Over ten years ago, I saw it was used from time to time, but nowadays, I rarely see its use, except by programmers from older generations.

# Answer
> 5 votes
The generic command `sh` is quite literally an industry standard, a POSIX standard, to be precise (IEEE 1003.2 and 1003.2a, available for purchase for hundreds of dollars at various websites). In theory, any script that starts with `#!/bin/sh` should conform to this standard. In practise, most Linux systems have a shell that is close to this standard, but has a few quirks and extensions.

Problems crop up when these quirks and extensions become standard practise in shell scripts. The Debian operating system changed to `dash` as their `sh` shell to encourage people to stop using "bashisms" in shell scripts that didn't specify a particular shell, i.e. those that began with `#!/bin/sh`. The `dash` shell tries to be as standards-compliant as possible:

> dash is the standard command interpreter for the system. The current version of dash is in the process of being changed to conform with the POSIX 1003.2 and 1003.2a specifications for the shell. This version has many features which make it appear similar in some respects to the Korn shell, but it is not a Korn shell clone (see ksh(1)). Only features designated by POSIX, plus a few Berkeley extensions, are being incorporated into this shell. This man page is not intended to be a tutorial or a complete specification of the shell.

I'm not familiar with the differences, and generally try to stick to the `sh` man pages to instruct me with regards to correct standards-compliant shell scripts.

---
Tags: linux, shell
---


# Question
Title: Parallel processing of scripts that use obitools
I have an obitools script (de Barba et al. 2016) that I would like to run faster. How would you run it in parallel to cut down on time?

```
illuminapairedend -r rawdata_scandinavia_R2.fastq rawdata_scandinavia_R1.fastq | tee rawdata_scandinavia.fastq | obiannotate -S goodAli:'"Alignement" if score>40.00 else "Bad"' | obisplit -t goodAli -p rawdata_scandinavia.
touch rawdata_scandinavia.Bad.fastq
touch rawdata_scandinavia.Alignement.fastq
touch rawdata_scandinavia.fastq
ngsfilter -t rawdata_scandinavia.ngsfilter -u rawdata_scandinavia.unidentified.fastq rawdata_scandinavia.Alignement.fastq > rawdata_scandinavia.filtered.fastq
obisplit -p MICROSAT.PCR_ -t experiment rawdata_scandinavia.filtered.fastq

```

Can someone describe what tools and steps one needs to take to run this job in parallel, given that the functions are not designed to be run in parallel out of the box?

De Barba, M., Miquel, C., Lobréaux, S., Quenette, P. Y., Swenson, J. E., & Taberlet, P. (2016). High-throughput microsatellite genotyping in ecology: improved accuracy, efficiency, standardization and success with low-quantity and degraded DNA. Molecular Ecology Resources, 1–16. https://doi.org/10.1111/1755-0998.12594

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

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

```
#!/bin/bash

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

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

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

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

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

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

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

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

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

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

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

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

# remove intermediate results to conserve space
rm to_ngsfilter*

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

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

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

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

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

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

echo "Done. I think."

```

---
Tags: ngs, linux
---