# Question
Title: Given a Genomic Ranges of SNPs, how to inject these SNPs in genome via BSGenome?
Let's say I have the genome `hg19` loaded into R via BSGenome

```
library("BSgenome")
hg19genome = getBSgenome('BSgenome.Hsapiens.UCSC.hg19', masked=FALSE)

```

I then have a list of SNPs loaded as a GRanges object, `gr`

```
library(GenomicRanges)

 > gr
 GRanges object with 212 ranges and 3 metadata columns:
        seqnames               ranges strand |     width     REF    ALT
           <Rle>          v <IRanges>  <Rle> | <numeric>
    [1]        1 [86099032, 86099032]      * |         1     C      T
    [2]        1 [86099033, 86099033]      * |         1     C      A
    [3]        1 [86099199, 86099199]      * |         1     G      A
  ....

```

Is there a straightforward way to inject these SNPs into `hg19genome`?

# Answer
This seems relatively complicated given the structure of a BSGenome object.

The creator of the package answered this question previously on the Bioconductor support forums: https://support.bioconductor.org/p/86665/#86757

> We don't provide an easy way to inject arbitrary SNPs in an arbitrary BSgenome at the moment. However, it should not be too hard to forge the BSgenome ... package. First you would need to compute the sequences of the mutated chromosomes (you can use replaceLetterAt for this), then write them to a 2bit file (put them in a DNAStringSet object and call rtracklayer::export on it), then use that 2bit file to forge the BSgenome .... package (see the BSgenomeForge vignette in the BSgenome package for how to do this). \- Herve Pages

Here is another similar question with more info about replaceLetterAt: https://support.bioconductor.org/p/26199/
> 3 votes

---
Tags: r, bioconductor, reference-genome, snp, biostrings
---


title: Human body sites knowledge base

text:
# Question
Title: Human body sites knowledge base
I'm looking for a lightweight knowledgebase describing the human body to annotate disease sites. I do not need a great level of detail, I just need a kind of basic organs/sub-organs taxonomy. I checked out the *Foundational Model of Anatomy ontology*, but it's over-kill for my purposes and really cumbersome (Protege fails when it tries to load FMA.owl).

Is there a simple online knowlegebase available for this?

# Answer
> 6 votes
I am the developer of Uberon and I would be happy to help you with what you need, either from Uberon, or from the FMA.

You mentioned you need something simpler than FMA. There are a variety of tools for creating custom subsets. Additionally, some ontologies provide ready-made subsets for particular purposes.

For example, we make a subset called 'basic.obo' from Uberon: http://uberon.github.io/downloads.html#subsets

You can examine this directly in either oboedit or Protege5 (note Protege5 parses obo so you don't need to convert).

From your question, you didn't indicate you needed programmatic access, but if you do there are a variety of options, as well as more programmer-friendly JSON exports.

Most ontologies have a number of different options for web-browsing, you can see some of the options for Uberon here: http://obofoundry.org/ontology/uberon.html

You mentioned annotating disease sites. There may be a tools (either desktop or web-based) for helping with this - e.g. providing autocomplete over a desired subset of an ontology.

You may also want to ask your question on a mailing list like obo-discuss. Unfortunately many ontology developers don't read stackexchange, but most are always keen to hear of community requirements, to help them provide better and more useful products.

# Answer
> 3 votes
Search in bioportal: For example: **SNOMED**, **UBERON**. Please, note that **FMA** is also available as RDF, csv , etc...

You could also try to convert FMA to obo and try to open it with OBO-edit

---
Tags: database, ontology
---


title: Merging bed records based on name

text:
# Question
Title: Merging bed records based on name
I generated a file starting with the following bed lines:

```
$ head -6 /tmp/bed_with_gene_ids.bed
I   3746    3909    "WBGene00023193"    .   -
I   3746    3909    "WBGene00023193"    .   -
I   4118    4220    "WBGene00022277"    .   -
I   4118    4358    "WBGene00022277"    .   -
I   4118    10230   "WBGene00022277"    .   -
I   4220    4223    "WBGene00022277"    .   -

```

I would like to merge them based on the name field (the 4-th column), taking the min for the start and the max for the end. Other fields are expected to be the same for all records having the same name.

Expected result:

```
I   3746    3909    "WBGene00023193"    .   -
I   4118    10230   "WBGene00022277"    .   -

```

I found a potential solution based on `bedtools groupby` here: https://www.biostars.org/p/145751/#145775

---

> Sample data:

```
cat genes.bed
chr14    49894259    49895806    ENSMUST00000053290    0.000000    ...
chr14    49894873    49894876    ENSMUST00000053290    0.000000    ...
chr14    49894876    49895800    ENSMUST00000053291    0.000000    ...
chr14    49895797    49895800    ENSMUST00000053291    0.000000    ...
chr14    49901908    49901941    ENSMUST00000053291    0.000000    ...

```

> Example output:

```
sort -k4,4 genes.bed \
| groupBy -g 1,4 -c 4,2,3 -o count,min,max \
| awk -v OFS='\t' '{print $1, $4, $5, $2, $3}'

chr14    49894259    49895806    ENSMUST00000053290    2
chr14    49894876    49901941    ENSMUST00000053291    3

```

---

However:

1. I don't understand the groupBy behaviour (Why `-g 1,4` and not just `-g 4`?, Why `-c 4,2,3` in this order and then rearrange things using `awk`?)
2. This code doesn't work for me.

Here is what happens when I try the solution given above:

```
$ head -3 /tmp/bed_with_gene_ids.bed | bedtools groupby -g 1,4 -c 4,2,3 -o count,min,max | awk -v OFS='\t' '{print $1, $4, $5, $2, $3}'
3           3746    4220

```

Here are attempt based on what I thought could work according to the documentation:

```
$ head -6 /tmp/bed_with_gene_ids.bed | bedtools groupby -g 4 -c 1,2,3,4,5,6 -o first,min,max,distinct,first,first
I   3746    10230   "WBGene00022277","WBGene00023193"   .   -

$ head -6 /tmp/bed_with_gene_ids.bed | bedtools groupby -g 4 -c 1,2,3,4,5,6 -o first,min,max,last,first,first
I   3746    10230   "WBGene00022277"    .   -

$ head -6 /tmp/bed_with_gene_ids.bed | bedtools groupby -g 4 -c 1,2,3,5,6 -o first,min,max,first,first
I   3746    10230   .   -

```

I don't get why when I group based on the 4-th column, for which I have two distinct values, I cannot obtain two lines in the resulting output.

I understand based on the comments on the documentation page that the documentation is not up-to-date. In particular, there is a `-full` option that is needed if one wants all fields to be outputted. Re-reading the solution mentioned above, I think I now understand the reason for the multiple columns for the `-g option` and for the `awk` rearrangement. Hence the following attempt.

```
$ head -6 /tmp/bed_with_gene_ids.bed | bedtools groupby -g 1,4,5,6 -c 2,3 -o min,max -full
    I   3746    3909    "WBGene00023193"    .   -   3746    10230

```

But this still doesn't give me two lines.

Are there other tools that could do what I want efficiently?

---

### Edit: Solution

According to this answer, the problem with bedtools is that there is a bug in the latest release (2.26.0 as of august 2017). In order to have a functional `bedtools groupby`, one needs to get the development version from github.

With the github version of bedtools, I can now get the expected result as follows:

```
$ head -6 /tmp/bed_with_gene_ids.bed | bedtools groupby -g 1,4,5,6 -c 2,3 -o min,max | awk -v OFS="\t" '{print $1,$5,$6,$2,$3,$4}'
I   3746    3909    "WBGene00023193"    .   -
I   4118    10230   "WBGene00022277"    .   -

```

I include fields 1, 5 and 6 in `-g` (besides field 4) in order to have them printed out. In my bed file, they should be the same for a given value of field 4. The `awk` part is needed because one has apparently not total control on the output order: the `-g` fields come before the `-c` fields.

Note (added 23/12/2020): As specified in `bedtools groupby` documentation, the input bed file should be sorted according to (chromosome (1) and) name (4) field(s) for this to work properly. Otherwise, duplicate (i.e. non-merged) names may appear in the output.

# Answer
Although you don't mention it, I'm guessing you're using bedtools v2.26.0. Version 2.26.0 of groupBy has a bug in it, which you've encountered (it was fixed shortly after release, so you'll either have to use a version before the bug was introduced, or compile the current source code yourself from https://github.com/arq5x/bedtools2)

v2.26.0:

```
local10:~/Documents/tmp$ cat asdf.bed
I 3746 3909 WBGene00023193 . -
I 3746 3909 WBGene00023193 . -
I 4118 4220 WBGene00022277 . -
I 4118 4358 WBGene00022277 . -
I 4118 10230 WBGene00022277 . -
I 4220 4223 WBGene00022277 . -
local10:~/Documents/tmp$ groupBy -i asdf.bed -g 4 -c 2,3 -o min,max
3746    10230

```

v2.26.0-125-g52db654 (I.E. compiling the source code from github):

```
local10:~/Documents/tmp$ bedtools2/bin/groupBy -i asdf.bed -g 4 -c 2,3 -o min,max
WBGene00023193  3746    3909
WBGene00022277  4118    10230

```

---

To answer your questions:

1) You might notice that my output above gives the grouped columns first; you'll have to reorder the output via awk in order to get it back in order. As for why they chose to group on both columns 1 and 4: if you have the same name on multiple chromosomes, you may want to treat them as separate features.

2) Version differences, as stated in the first part of my answer.

---

To actually merge the file:

Make sure to run this with a version other than v2.26.0 (As Devon Ryan writes in the comments, you may want to add column 6 to `-g` to make it strand-specific):

```
./bedtools2/bin/groupBy -i asdf.bed -g 1,4 -c 2,3,5,6 -o min,max,first,first \
     | awk -v OFS='\t' '{print $1, $3, $4, $2, $5, $6}'
I   3746    3909    WBGene00023193  .   -
I   4118    10230   WBGene00022277  .   -

```
> 6 votes

# Answer
You could do this with the `CGAT` toolkit:

`cgat bed2bed --method=merge --merge-by-name -I bed_with_gene_ids.bed`

Installing such a massive package might be overkill for this task though.
> 4 votes

# Answer
You can do this easily with Hail. Hail primarily uses BED files to annotate genetic datasets (see the last annotate\_variants\_table example), but you can manipulate BED files using Hail's general facilities for manipulating delimited text files. For example:

```
$ cat genes.bed
I   3746    3909    "WBGene00023193"    .   -
I   3746    3909    "WBGene00023193"    .   -
I   4118    4220    "WBGene00022277"    .   -
I   4118    4358    "WBGene00022277"    .   -
I   4118    10230   "WBGene00022277"    .   -
I   4220    4223    "WBGene00022277"    .   -

```

The Hail script (python code):

```
from hail import *
hc = HailContext()
(hc
 .import_table('genes.bed', impute=True, no_header=True)
 .aggregate_by_key('f0 = f0, f3 = f3',
    'f1 = f1.min(), f2 = f2.max(), f4 = ".", f5 = "-"')
 .select(['f0', 'f1', 'f2', 'f3', 'f4', 'f5'])
 .export('genes_merged.bed', header=False))

```

The result:

```
$ cat genes_merged.bed
I   3746    3909    WBGene00023193  .   -
I   4118    10230   WBGene00022277  .   -

```

I aggregate over chrom and name so this solution won't merge entries on different chromosomes. The `select` is necessary to reorder the fields because `aggregate_by_key` places the keys being aggregated over first.

Disclosure: I work on Hail.
> 3 votes

# Answer
```
$ cut -f4-6 in.bed | sed 's/\t/_/g' | sort | uniq | awk -F'_' '{ system("grep "$1" in.bed | bedops --merge - "); print $0; }' | paste -d "\t" - - | sed 's/_/\t/g' | sort-bed - > answer.bed

```

Given your sample input:

```
$ more in.bed
I   3746    3909    "WBGene00023193"    .   -
I   3746    3909    "WBGene00023193"    .   -
I   4118    4220    "WBGene00022277"    .   -
I   4118    4358    "WBGene00022277"    .   -
I   4118    10230   "WBGene00022277"    .   -
I   4220    4223    "WBGene00022277"    .   -

```

The `answer.bed` file:

```
$ more answer.bed
I   3746    3909    "WBGene00023193"    .   -
I   4118    10230   "WBGene00022277"    .   -

```

Sorting with `sort-bed` is useful at the end, so that you can pipe it or work with it with other BEDOPS tools, or other tools that now accept sorted BED input.

Streaming is a pretty efficient way to do things, generally.

---

**How this works**

Here's the pipeline again:

```
$ cut -f4-6 in.bed | sed 's/\t/_/g' | sort | uniq | awk -F'_' '{ system("grep "$1" in.bed | bedops --merge - "); print $0; }' | paste -d "\t" - - | sed 's/_/\t/g' | sort-bed - > answer.bed

```

We start by cutting columns 4 through 6 (id, score and strand), replacing tabs with underscores, sorting and removing duplicates:

```
cut -f4-6 in.bed | sed 's/\t/_/g' | sort | uniq

```

What we get out of this is a sorted list of "needles" — one for each ID-score-strand combination: an *ID-needle* — that we can use to `grep` or filter the original BED file.

This list is piped to `awk` which, for each ID-needle, runs `grep` against the original BED file and pipes the subset to `bedops --merge -`, which merges overlapping intervals.

Note that merging only works for overlapping intervals. *Merging is not necessarily the same as returning a min-max pair, and this pipeline will break if there are intervals that do not overlap.* But you could modify the `awk` statement to process the input intervals and return the minimum and maximum interval coordinates, if that is really what you want, by tracking the min and max values over all intervals that come into `awk`, and printing a final interval with an `END` block.

The `system` command prints the merged interval on one line. The following `print $0` statement prints the needle on the next line:

```
awk -F'_' '{ system("grep "$1" in.bed | bedops --merge - "); print $0; }'

```

We take each pair of alternating lines and re-linearize them with `paste`. This result now contains four columns: the three columns of each merged interval, and the ID-needle.

We then use `sed` to replace underscores with tabs, so that we turn the ID-needle back into three, tab-separated ID-score-strand columns:

```
paste -d "\t" - - | sed 's/_/\t/g'

```

The output is now a six-column BED file, but it is ordered by the sort order we applied on ID-needles further up the pipeline, which we don't want. What we really want is BED that is sorted per BEDOPS `sort-bed`, so that we can do more set operations and get a correct result. So we pipe this to `sort-bed -` to write a sorted file to `answer.bed`:

```
sort-bed - > answer.bed

```
> 2 votes

# Answer
If you're 100% sure that everything but the start and end positions will be the same for all lines sharing a name, you could just do it yourself. For example, in Perl:

```
$ perl -lane '$start{$F[3]}||=$F[1];
              if( $F[1] < $start{$F[3]} ){
                 $start{$F[3]} = $F[1]
              }
              if( $F[2] > $end{$F[3]} ){
                 $end{$F[3]} = $F[2]
              }
              $chr{$F[3]} = $F[0];
              $rest{$F[3]} = join "\t", @F[4,$#F];
              END{
                 foreach $n (keys %chr){
                  print "$chr{$n}\t$start{$n}\t$end{$n}\t$n\t$rest{$n}"
                 }
              }' file.bed
I   3746    3909    "WBGene00023193"    .   -
I   4118    10230   "WBGene00022277"    .   -

```
> 0 votes

---
Tags: bed, bedtools
---


# Question
Title: Tools to create annotated table of variants from VCF
The problem: I have a VCF file, a reference genome, and a bunch of annotations for the reference (genes, repeat regions, etc.) as GFF or BED files.

What I would like is a tool that takes all of this as input and outputs a tab- or comma-delimited table containing as much information as possible. Potential columns in the output include:

1. Type of variant (SNV, structural, etc)
2. Details of variant (e.g. reference base, variant base, coverage, position, etc)
3. Annotations overlapping variant
4. Annotations near variant (e.g. is it just upstream of a gene)
5. If it appears in a coding region, does it change the amino acid.

There are many tools that do *something like this*. But to the newbie (like me) it is not clear which tools are worth starting with. Since the majority of tools take a little effort to get working in the first place, my question is:

What tool comes closest to doing what I have described, is known to reliable, and is likely to be maintained for the next few years.

# Answer
> 5 votes
I have used FEATnotator and I think it can provide all of the columns you would like to see. There are many output files generated, but the consolidated output has the following columns:

* Chromosome
* Position
* Column\_3
* Consensus\_Allele
* Annotation\_Signature
* Reference\_Base
* Alternate\_Base
* Transition/Transversion SNP\_Type
* Premature\_STOP\_Gained STOP\_Lost
* START\_CODON
* STOP\_CODON
* SPLICE\_SITE
* InterGenic
* Gene\_Body
* Intron
* Exon
* Coding
* UTR
* Transcription\_Start\_Site
* Nearest\_gene
* Distance

here is some example output records from using a VCF, reference genome and GFF annotation file (sorry about the crappy formatting...there are a lot of fields!):

```
Chr01   28177   T       G       NA      T     G TRANSVERSION    NA      NO      NO      NO      NO      NO      YES     NO      NO      NO      NO      NO      NO      Eucgr.A00211.v2.0       16054
Chr01   44876   C       G       gene:Eucgr.A00211.v2.0   mRNA:Eucgr.A00211.1.v2.0[intron]       C       G       TRANSVERSION    NA      NO      NO      NO      NO      NO      NO      YES     YES     NO      NO      NO      NO      NA  NA
Chr01   46819   A       G       gene:Eucgr.A00211.v2.0   mRNA:Eucgr.A00211.1.v2.0   CDS:Eucgr.A00211.1.v2.0.CDS.4[AAA - K => GAA - E (MISSENSE)]   exon:Eucgr.A00211.1.v2.0.exon.4      A       G       TRANSITION      MISSENSE        NO  NO       NO      NO      NO      NO      YES   NO        YES     YES     NO      NO      NA      NA

```

# Answer
> 9 votes
Hail might be an option for you.

It is actively developed by a growing team at the Broad. It is rigorously tested (continuous integration, continuous deployment, bug reports get regression tests, blah blah blah).

It was designed to solve this problem (among others). It can import a variety of formats, including VCF, TSV, UCSC BED, JSON and interval files. (We don't have explicit support for GFF, but we can probably handle them with general facilities. If not, get in touch and we'll add support.) It can call out to VEP (and soon Nirvana, Illumina's VEP rewrite). It has general facilities to transform, filter, clean and query data.

What's more, we've curated a large collection of annotation resources (currently ~22 databases of annotations + VEP), hosted in a public bucket in Google cloud, and built an interactive query builder to select which resources you want to use. Get in touch if you'd like us to add additional resources. You can read more about it here.

Disclaimer: I work on Hail.

# Answer
> 5 votes
snpEff is a great tool for annotating VCF files and you can add custom reference sequences.

http://snpeff.sourceforge.net/

Guide to add custom annotation files in snpEff

https://gatkforums.broadinstitute.org/gatk/discussion/50/adding-genomic-annotations-using-snpeff-and-variantannotator

There are a bunch of pre curated annotation datasets available in their database.

---
Tags: sequence-annotation, vcf
---


# Question
Title: How can I do an overlapping sequence count in Biopython?
Biopython's `.count()` methods, like Python's `str.count()`, perform a non-overlapping count, how can I do an overlapping one?

For example, these code snippets `return 2`, but I want the answer `3`:

```
>>> from Bio.Seq import Seq
>>> Seq('AAAA').count('AA')
2
>>> 'AAAA'.count('AA')
2

```

# Answer
For Biopython 1.70, there is a new `Seq.count_overlap()` method, which includes optional `start` and `end` arguments:

```
>>> from Bio.Seq import Seq
>>> Seq('AAAA').count_overlap('AA')
3
>>> Seq('AAAA').count_overlap('AA', 1, 4)
2

```

This method is also implemented for the `MutableSeq` and `UnknownSeq` classes:

```
>>> from Bio.Seq import MutableSeq, UnknownSeq
>>> MutableSeq('AAAA').count_overlap('AA')
3
>>> UnknownSeq(4, character='A').count_overlap('AA')
3

```

---

Disclaimer: I co-contributed the `.count_overlap()` methods with Peter Cock, see `97709cc`
> 18 votes

# Answer
I've encountered this problem before, and used python `re` module to solve this problem.

```
import re
all = re.findall(r'(?=(AA))','AAAA')
counts = len(all)

```

You can get more details in this thread
> 9 votes

# Answer
You can use `finditer` from python's `re` module. The advantage of this approach is it allows for getting the indices of those matches, which could be handy down the track.

```
>>> import re
>>> matches = re.finditer(r'(?=(AA))', 'AAAA')
>>> indices = [match.span(1) for match in matches]
>>> indices
[(0, 2), (1, 3), (2, 4)]
>>> num_matches = len(indices)
>>> num_matches
3

```
> 4 votes

---
Tags: python
---


# Question
Title: How to filter intervals (reads or genomic coordinates) that have the exact same 5' or 3' ends?
Say I have reads that overlap some genes that produce small RNAs, but I want only those reads that start at exactly the TSS of the loci. In other words, reads whose 5' end match the 5' end of a genomic feature.

`5'....3'
...+++++++... Gene
...0000000... R1
...00000000.. R2
..000000000.. R3
..00000000... R4`

The output should be R1 and R2.

I believe this is a fairly common bioinformatics operation but somehow it doesn't seem to be an option in the general use tools I looked at (py-bedtools, samtools, bedops).

A solution I thought of would be:

1. reduce the gene coordinates to its 5' end, and
2. using bedtools `bedClosest` to annotate the distance between reads and genes, and finally
3. select reads that overlap (d=0) and write them to a (bam) file.

This involves quite a bit of wrangling data with `pybedtools` or `awk/bash`, and I wonder is there is a more elegant solution to this. Solutions in `R/Bioconductor` are also welcome.

# Answer
> 1 votes
Presuming you have deepTools installed:

```
#/!usr/bin/env python
from deeptoolsintervals.parse import GTF
import pysam

anno = GTF(["some/gtf/or/bed/file", "or/more/than/one/if you like/"])
bam = pysam.AlignmentFile("something.bam")
for read in bam:
    s = read.reference_start
    if read.is_reverse:
        s = read.reference_end
    e = s + 1
    overlaps = anno.findOverlaps(read.reference_name s, e, matchType=4)
    # For a GTF, parse the results to see if this is the correct type.
    # For a BED3/BED6 file there will be no exons, so anything other than
    # None indicates an appropriate overlap

    do something

```

`deeptoolsintervals` in the GTF/BED parsing and custom interval tree library used by deepTools. `anno` holds the interval tree and the associated information (gene/transcript name or ID, exons, strand, etc.). the `matchType=4` is the pivotal part since that specifies that only overlaps having the exact same start position should be included. Assuming you have a stranded library, you might instead do something like switch between `matchType=4` (the same start position) and `matchType=5` (the same end position) depending on the orientation of the alignment and the the `strand=` of the gene.

Note that you'll need to make a `set()` to hold read names to find the mates (that or jump around the file to get them, which is probably slower).

---
Tags: bedtools
---


# Question
Title: Difference between CPM and TPM and which one for downstream analysis?
What the difference between TPM and CPM when dealing with RNA seq data?

What metrics would you use if you have to perform some down stream analysis other than Differential expression for eg.

Clustering analysis using Hclust function and then plotting heat map to find differences in terms of expression levels, correlation and pca

Is it wrong to use TPM for such analysis, if yes then when does one use TPM versus CPM.

# Answer
> 13 votes
You can find the various equations in this oft-cited blog post from Harold Pimentel. CPM is basically depth-normalized counts, whereas TPM is length-normalized (and then normalized by the length-normalized values of the other genes).

If one has to choose between those two choices one typically chooses TPM for most things, since generally the length normalization is handy. Realistically, you probably want `log(TPM)` since otherwise noise in your most highly expressed genes dominates over small expression signals.

# Answer
> 3 votes
Neither CPM nor TPM are well suited here, because neither performs robust cross-sample normalisation (see the blog post Devon linked to).

DESeq2 provides two robust log-space normalisation methods for downstream analysis, the *regularised log* (`rlog`), and the *variance stabilising transformation* (`vst`). The DESeq2 vignette explains how to use these for things like hclust.

---

On a more general note, CPM does not account for transcript length differences, while TPM does. If the choice is between TPM and CPM I would therefore use TPM. However, if you are only comparing the same transcripts across experiments, the transcript length is actually invariant so it doesn’t matter (but CPM is still not a good cross-experiment normalisation).

---
Tags: rna-seq
---


title: How to convert fastq to fast5

text:
# Question
Title: How to convert fastq to fast5
`fast5` is a variant of `HDF5` the native format in which raw data from Oxford Nanopore MinION are provided. You can easily extract the reads in fast5 format into a standard fastq format, using for example `poretools`.

Say I have aligned these reads in `fastq` format to an external reference genome, resulting in a `SAM` file. Say I have then taken a subset of the `SAM` file, according to the bitwise flag, to include only the reads that map to the reference. With the read ID, I can then grep them out from the file containing the reads in `fastq` format, generating a subset file in `fastq` format containing only the IDs that have mapped to the reference.

Now my question is, can we subset reads from the `fast5` archive according to the list of mapping reads as taken from the file with reads in `fastq` format? This is for educational purposes, so that we have a smaller starting archive, and the `fast5` -\> `fastq` extraction takes less cpu time.

# Answer
> 6 votes
**NOTICE**: I have altered my answer slightly from the original as I have turned the original script into a `pip` installable program (with tests) and have updated the links and code snippets accordingly. The essence of the answer is still exactly the same.

---

This is something I have been meaning to get around to for a while, so thanks for the prompt.

I have created a python program called `fast5seek` to do what (I think) you're after.

As you have mentioned this is for educational purposes I have added a tonne of comments to the code too so I think you shouldn't have any issues following it.

The docs on the GitHub repo have all the info, but for those reading along at home

```
pip3 install fast5seek
fast5seek -i /path/to/fast5s -r in.fastq in.bam in.sam -o out.txt

```

What it does is read in `<in.fastq|in.bam|in.sam>` and extract the read id from each header. It then goes through all the fast5 files under `/path/to/fast5s` and checks whether their read id is in the set of read ids from `<in.fastq|in.bam|in.sam>`. If it is, the path to the file is written to it's own line in `out.txt`.

If no output (`-o`) is given, it will write the output to stdout.

So if you wanted to pipe these paths into another program, you could do something like

```
mkdir subset_dir/
fast5seek -i /path/to/fast5s/ -r in.fastq | xargs cp -t subset_dir/

```

The above example would copy the fast5 files that are found in your fastq/BAM/SAM to `subset_dir/`.

# Answer
> 3 votes
I do this very frequently, using the read file name to identify FAST5 files associated with particular reads. If a FASTQ record includes the channel number and read number (and preferably the runID as well), then I use that information to find the associated FAST5 read.

If reads are called with Albacore, then the `sequencing_summary.txt` file has additional information. The first column is the exact name of the file that was called, the second column is the read ID of the FASTQ sequence, and the third column is the run ID associated with that read ID. This is much easier to work with, but does require calling reads with Albacore (which, admittedly, is what seems to be producing the best results at the moment).

---
Tags: nanopore, software-recommendation, fastq, format-conversion, fast5
---


# Question
Title: Can gene co-expression networks be used to help identify differentially expressed genes?
In my RNAseq dataset of differentiated stem-cell lines, some samples have far fewer significantly differentially expressed genes than others. QC shows that this is because there are way fewer reads for one sample than the others.

Can gene co-expression networks be used to infer differential expression? For example, if certain genes just miss the cutoff q-value cutoff If so, what papers have done so?

# Answer
> 2 votes
While you can use networks to find differentially expressed genes (see the WGCNA package, which does this) in my experience this ends up largely matching what you'd get using a traditional package with a looser threshold for significance. Given the time savings of traditional packages, there's rarely any gain to using networks (it won't hurt to try, just note that it'll take some time). If some very interesting genes are just above your significance threshold then change your threshold. You need to validate your findings in some way anyway, so your p-value threshold is partly just a way to protect you from wasting your time (but if background knowledge suggests that your time wouldn't be wasted...).

# Answer
> 2 votes
Co-expression network will give you an idea about the genes having similar expression patterns and the nodes will be decided on the basis of correlation scores and nothing much you will get from differential gene expression perspective.

I would suggest you to try some other models for differential expression analysis like baySeq , DESeq, edgeR, NOIseq. Depending upon the model you use for analysis you might get some extra significant DEG. Also, try different significance cutoff and see how it's changing the analysis result to get the threshold that suits the data best.

---
Tags: rna-seq, differential-expression, networks
---


# Question
Title: Why are some of the gene peptides returned by biomaRt missing an asterisk in the end?
*Note: this question has also been asked on BioStars*

I am getting peptides using biomaRt library in R for the gene 'BRCA1'. 27 different aminoacid sequences are returned, and 12 of those have an asterisk (\*) in the end, while 15 do not. What does that mean? I know that asterisk is used to indicate stop codon. So, does that mean the peptides without an asterisk do not have a stop codon? Is that possible at all? Are all of those valid peptides or not?

Below is the R code I used:

```
library(biomaRt);
mart = useMart('ensembl', dataset='hsapiens_gene_ensembl');
seq = getSequence(id='BRCA1', type='hgnc_symbol', seqType='peptide', mart = mart);

```

`MDLSALRVEEVQNVINAMQKILECPICLELIKEPVSTKCDHIFCKVLLCCPSWSTVVRS*` is an example sequence with an asterisk in the end. `MDLSALRVEEVQNVINAMQKILECPICLELIKEPVSTKCDHIFCKSLQESTRFSQLVEELLKIICAFQLDTGLEYANSYNFAKKENNSPEHLKDEVSI` is an example without an asterisk.

# Answer
If it's 3' incomplete that means the evidence used to create it was a fragment. Here's the evidence used to construct BRCA1-214 ENST00000477152.5, a 3' incomplete.

You can see that there's a full length cDNA from EMBL, AK307553.1, which was used to create this model.

The sequence was mapped against the genome sequence to create the transcript model. When searching for the translation, it was found that the open reading frame was started, but there was no stop codon. This suggests that the cDNA AK307553.1 is actually a fragment, that the mRNA was broken or cleaved before it was reverse-transcribed sequenced. We display it in Ensembl in the hope that this will lead someone to identify the full length transcript that it represents.

Many people choose not to work with these incomplete transcripts. You can do this by filtering by Gencode Basic in biomaRt.
> 8 votes

# Answer
If you look at the original data at Ensembl you'll notice that most of these are labeled "CDS 3' incomplete" ~~and have a TSL (transcript support level) of 1, which is as low as it goes~~. It seems likely that that is simply an incomplete annotation. I'm not surprised that there are a bunch of these for BRCA1, since there was a long time when there were a LOT of legal issues surrounding working on it.
> 4 votes

---
Tags: gene, proteins, biomart
---


# Question
Title: Split FASTQ and matching BAM into matching chunks
I am running a slow downstream analysis on a large set of nanopore reads (approx 3 million), and would like to split them into smaller chunks, run the analysis in massively parallel, and then recombine. Originally I just split the FASTQ into chunks, re-aligned each chunk, and then merged the output, but here I would like to use an existing alignment so I can compare results with existing analyses (i.e. the alignments must be the same).

How can I efficiently split a FASTQ file and a BAM file giving the alignment of the FASTA file into chunks, ensuring that all of the reads in FASTQ chunk 1 are in BAM chunk 1, vice versa and so on?

My FASTQ is approximately 45GB and my BAM is 33GB, so I would prefer to avoid storing one of the two files in memory if possible.

EDIT: Here's some pseudocode of exactly what I'm trying to do:

```
# input: in.bam, in.fastq, chunk_size
i <- 0
for fastq_read in in.fastq:
    bam_read <- extract fastq_read.read_name from in.bam
    n <- i modulo chunk_size
    write fastq_read to out.n.fastq
    write bam_read to out.n.bam
    i <- i + 1

```

I could swap the above to iterate through the bam file and fetch from the fastq if that's easier.

# Answer
> 5 votes
Hmm, it's hard to think of a super efficient way of doing this (assuming the files aren't ordered the same - if they are then this whole answer is basically redundant). And also assuming the read ids for both files aren't a perfect intersection.

Off the top of my head you probably want to build a set of read ids for the fastq file and another for the bam. Some python code to get started on that:

```
import pysam
import itertools

def get_read_id_fastq(ref_path):
    """Extracts the read ids from a fastq file."""
    read_ids = set()
    with open(ref_path, 'r') as ref:
        for line in ref:
            if line.startswith('@'):  # i.e if line is header
                # split the line on spaces, take the first element, remove @
                read_id = line.split(' ')[0].replace('@', '')
                read_ids.add(read_id)

    return read_ids

def get_read_id_bam(ref_path):
    """Extract the read ids from a BAM file."""

    read_ids = set()
    with pysam.AlignmentFile(ref_path, 'r', ignore_truncation=True) as ref:
        for read in ref:
            # query_name is the query template name
            read_ids.add(read.query_name)

    return read_ids

fastq_ids = get_read_id_fastq(fastq_path)
bam_ids = get_read_id_bam(bam_path)

```

Then take the intersection of these two sets and that's your common read ids.

```
common_ids = fastq_ids & bam_ids

```

The next part will be a bit more involved. You will have to iterate through each file one at a time. I would suggest that for the first one you iterate through, create a running dictionary with the read id that was written as the key and the 'chuck number' it was written to as the value. You could create a cycle for your chunk size to manage this effectively

```
chunk_cycle = itertools.chain(*zip(range(chunk_size)*len(common_ids)))
write_idx = {}

```

The next part will probably require you to put in some hard-fought times getting it to work. I'll put in some rough pseudo-code to give you an idea.

```
for line in file:
    read_id = # get line read_id
    if read_id in common_ids:
        chunk_num = chunk_cycle.next()
        write_idx[read_id] = chunk_num
        file_to_write_to = 'out.{}.bam'.format(chunk_num)
        # open this file or write to it if already open

```

When you go to do the next file when you find a read id is in the common set, you look up it's value in `write_idx` and this will give you the chunk number to write to.

The reads wont be in the exact same order between the two files, but you could sort later if you needed it (not sure it would matter?).

Hope this helps. Sorry I couldn't give you more, but this should give you a head start hopefully.

---
Tags: sam, fasta, nanopore, fastq, subset
---


# Question
Title: How to count reads in bam per bed interval with bedtools
I recently installed Ubuntu 16.04 (because I was still using 12.04). But it seems my bedtools scripts don't work properly anymore. I can't figure out how to use the new bedtools for my old ways. What I want to do is get the number of reads from a bam file, per interval from a bed file. It was very simple with my old version of Ubuntu and bedtools:

```
bedtools coverage -abam file.bam -b All_peaks.bed > file.cov.txt

```

But this doesn't seem to work anymore. I used to get small files as results (~3MB), but now huge files are created (~4GB) including read names??

I am using bedtools v2.25 now.

# Answer
The order of `-a` and `-b` switched at some point. You want:

```
bedtools coverage -a All_peaks.bed -b file.bam > file.cov.txt

```

For reference, this is the end of the help output in version 2.25:

```
Default Output:
     After each entry in A, reports:
       1) The number of features in B that overlapped the A interval.
       2) The number of bases in A that had non-zero coverage.
       3) The length of the entry in A.
       4) The fraction of bases in A that had non-zero coverage.

```

And this is the equivalent output from version 2.19:

```
Default Output:
     After each entry in B, reports:
       1) The number of features in A that overlapped the B interval.
       2) The number of bases in B that had non-zero coverage.
       3) The length of the entry in B.
       4) The fraction of bases in B that had non-zero coverage.

```

**Update**: The change in behavior happened in version 2.24:

> We have changed the behavior of the coverage tool such that it is consistent with the other tools. Specifically, coverage is now computed for the intervals in the A file based on the overlaps with the B file, rather than vice versa.
> 10 votes

---
Tags: sam, bedtools
---


# Question
Title: Mathematical notation for formulating a rank score
This is a cross-post from Mathematics forum. As no one has answered it yet, I fIgured post here as well.

I would like to describe a transform I used to rank my data points. I have recoded my variable columns from $10$ measurements with different range to $n$ (for each measurement) where $n \in \\{0,0.5,1\\}$.

For columns where measurements were expressed as categorical values I simply recoded the categorical value to either $0$, $0.5$ or $1$ and for continuous variables I recoded the lower quartile as $0$, interquartile as $0.5$ and upper quartile as $1$ and finally I summed the recoded values to produce a single score $\sum n\_{i1,..i10} $ for each row.

I am trying to write my method and I'd like to have it in mathematical notation. I was wondering if anyone could help me with this!

# Answer
Here is my solution to the problem. I am posting it here in case anyone else come up with the same idea but did not know how to formulate it in mathematical notation!

> $$ \psi\_i = \sum\_{i=1}^9 S\_i $$
>
> where $S\_i$ is the score function and defined as:
>
> $$ S\_i = \cases{ 0 & \text{if $\theta(x\_i) \< Q\_1$ or $\theta(x\_i) = $ benign/neutral} \\ 0.5 & \text{if $Q\_1 \leq \theta(x\_i) \< Q\_3$ or $\theta(x\_i) = $ possibly damaging/uncertain} \\ 1 & \text{if $\theta(x\_i) \geq Q\_3$ or $\theta(x\_i) = $ damaging} } $$
>
> The $\theta(x\_i)$ is the pathogenicity or conservation score for variant $x$ as defined by model $i$ and $Q$ denotes the quartile range for scores from M-CAP, CADD, GERP and Phylop models.
> 1 votes

---
Tags: modelling
---


title: GRCh38 centromeres mappability

text:
# Question
Title: GRCh38 centromeres mappability
Is the mappability of the centromeres in the GRCh38 genome reference similar to each other?

As far as I can remember when GRCh38 came out, the sequence of the centromeres was determined by a combination of sequencing data and software prediction.

Given the way the sequence of the centromeres has been determined, should we expect Illumina 2x150bp (or shorter, 2x75bp) reads to map relatively equally to all centromere sequences?

# Answer
> should we expect Illumina 2x150bp (or shorter, 2x75bp) reads to map relatively equally to all centromere sequences?

No. It has long been established that different chromosomes are associated with different centromeric sequences. It is sometimes possible to tell which chr a read is originated from based on its sequence.

The GRCh38 centromeres are trickier. As I remember, the centromeric sequences were computationally generated with a Markov chain (or something alike) modeled after Venter's genome. GRC can distinguish most chromosomes, but not all. Some alpha arrays are placed onto 2 or 4 chromosomes. The original GRCh38 keeps all 4 copies. When you download GRCh38 for mapping purposes, only one copy is retained; additional copies are hard masked along with PARs on chrY.

If you want to know more, see this paper and this.
> 3 votes

# Answer
I doubt it, unless you are asking if mapping would be similarly bad for all centromeres. Here are some repetitive structures (probably not centromeric) that I've found in the nanopore reads for "human" sample NA12878, produced by the nanopore-WGS consortium:

These structures are consistent in that they repeat lots of times, but the internal patterns can be quite different. Here are a few more:

Given that centromeres need to be uniquely linkable to a single chromosome, it would make sense to me if the centromere internal structure is unique to each chromosome.

It's possible to have highly-repetitive structures which don't map to each other. While I haven't dug too deeply into the human reads, I've looked at the 5 most-compressible regions in an assembled rodent parasite (*Nippostrongylus brasiliensis*) genome, and found no internal similarities between them:

One of the issues with assembly from Illumina reads is that these highly-repetitive regions get collapsed into a single repeat (or at best a region up to twice the fragment length). With internal repeat units that have identity of over 98%, assembling the true sequence is very difficult, even with an exact knowledge of paired read separation. Even if this were possible, it can be impossible to correctly place a read because multiple internal units could be identical (or similarly different) to the sequenced read.
> 3 votes

---
Tags: illumina, centromere, mappability
---


# Question
Title: How to calculate p-values for fold changes?
I'm currently working with data from a Luminex multiplex assay. In this assay, the concentrations of 17 different analyte proteins were identified in 12 groups in triplicate. One of these 17 groups was used as the control, and the log2 fold changes were calculated for the analyte concentration of each sample in each group using the average control concentration for that analyte.

However, now I would like to calculate a p-value for the identified fold changes if possible.

My current preliminary idea is to perform the t test for each group compared with the control group for each analyte, but I don't think that my group sizes are large enough for this. My statistical knowledge is lacking, so if there is a better way of calculating these p-values I don't know about it.

Is there a better way of obtaining these p-values?

# Answer
The key point here is whether or not the values are approximately normally distributed and whether any transformations can be applied to make them so.

For a t-test, the most important thing is that there is no relationship between the mean and the variance.

Its also correct that you don't want to divide by the mean of the controls. You want to retain all the information.

So here is how i'd start, first calculate the mean and standard deviation of each of your 18 groups (17 treatments + control), and plot this on a scatter plot (i.e. mean on X, SD on Y).

Is there a relationship between the two?

Repeat the exercise with log transformed values. Is the relationship more or less.

Pick the set of values that has the weakest mean-SD relationship. I'm guessing the log transformed might well be better.

If $X = (x\_1, x\_2, x\_3)$ is your controls and $Y\_i = (y\_{i,1}, y\_{i,2}, y\_{i,3})$ are the treated samples from treatment i, to test the difference between treatment and control, do a t-test on X and Y (not `T/mean(c)`)

for example in R

```
X = c(x_1, x_2, x_3)
Y_1 = c(y_11, y_12, y_13)
t.test (X, Y_1)

```

or in excel is if the three values for X are in column A and the three values for $Y\_1$ in column B

```
=T.TEST(A1:A3,B1:B3,2,3)

```

If you are going to do 17 of these, remember to do a multiple testing correction by multiplying the resulting p-values by 17.

If you want to compare treatment 1 ($Y\_1$) vs treatment 2 ($Y\_2$) do the comparison directly, not taking the control into account:

```
t.test(Y_1, Y_2)

```

Remember that if you use log transformed values, you calculate log fold changes using `mean(Y)-mean(X)` rather than `mean(Y)/mean(X)`
> 4 votes

# Answer
Your null hypothesis would be that the fold-changes are 0, so you can either do the T-test accordingly or simply do away with the fold changes and perform the T-test between the raw values in the group (this is preferable to performing a T-test of one group vs. 0, since it allows you to assess the expected variability around 0). Note, however, that 3 samples is pretty much the bare minimum needed for any statistics, so your power will be terrible and your estimation of variance will likely be pretty inaccurate. Consequently, take any p-values with an appropriate grain of salt.

To my knowledge there aren't any good alternatives for your situation without having a good expected background distribution.
> 3 votes

# Answer
You can't calculate a p-value on the fold-change values, you need to use the concentrations in triplicate thus giving a measure of the variance for the t-test to use.

t-test assumes your data are normally distributed, if they aren't you're going to get spurious p-values. If you aren't sure a non-parametric test like Wilcoxon is better. It will be less sensitive although with only 3 replicates your experiment is low powered anyway.
> 1 votes

# Answer
So to get a p value you must have one null hypothesis like the fold changes are significant or not. Then you test the hypothesis against the data and calculate the p value (likeliness or unlikeliness) for a given threshold.

In your case, t test is good enough to calculate the p-value.

http://willett.ece.wisc.edu/wp-uploads/2016/01/05b-TandP.pdf
> 0 votes

---
Tags: statistics
---


# Question
Title: Which measure should be used in a PCA or RNA-seq data? TPM or counts?
I'm trying to understand the magnitude of batch effects in my RNA-seq samples, and I was wondering which expression units are more suitable to draw a PCA. I'm thinking of either `counts` or `TPM`, but things like `rlog` or `vst` could work too.

Additionally, I'm wondering whether either of this units should be log-transformed first, to avoid high-abundance transcripts driving the PCA.

# Answer
tldr: log transform counts and TPMs, but rlog/vst are preferred

TPM should be log transformed to get more useful results. If you're using DESeq2 already (given the reference to `rlog` and `vst`, this seems likely), then please go ahead and use `rlog` or `vst`. That will give you more reasonable results than raw counts. If you're stuck with counts for some reason, then first use normalized counts so they're at least a bit more comparable and then log transform them so your highly expressed genes aren't driving everything.

**Edit**: As an aside, if you know what the batch effect is (e.g., library prep. date), it's sometimes convenient to include that in your model. You can then assess the genes that are actually changed due to that, which is sometimes handy to know (e.g., which genes might be more/less prone to degradation).
> 5 votes

# Answer
PCA works best when the input data is approximately normally distributed on each dimension. It would be a good idea to do some initial data quality checks to verify that this is the case (and transform the data appropriately if not), or at least verify that the data is approximately normally distributed in the aggregate.

For looking at Illumina RNASeq data, what worked best for me (i.e. produced the most normal-looking data) was the following steps:

1. Removing genes that had low raw counts in all samples
2. Using DESeq's variance-stabilized transform (which transforms counts into a log-like distribution)
3. Further normalising the VST values by dividing by the longest transcript length within each gene (which I call VSTPk)

These steps are stated in a bit more detail in our Th2 paper that was published at the end of last year:

http://jem.rupress.org/content/early/2016/12/01/jem.20160470#materials-methods
> 4 votes

---
Tags: rna-seq
---


title: Analyzing Illumina Counts

text:
# Question
Title: Analyzing Illumina Counts
I'm pretty new to all of this--forgive me if this is a simple question.

When I download illumina counts from GEO (like the supplementary file in GSE89225). Can I do comparisons directly on that file? Is there some normalization procedure I should go through?

Thanks!

# Answer
The counts files for GSE89225 is the output of HTSeq-count as a large matrix. Unless you are developing a differential expression package yourself you should not attempt to directly use this. Rather, you should load it into R and use packages such as DESeq2, edgeR, or limma (those are the most popular ones).

For convenience, in DESeq2 you would want the `DESeqDatasetFromMatrix()` function after loading this.
> 5 votes

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


# Question
Title: Is there a way to quickly verify the presence of some SNPs in Fastq files?
Is there a tool that can scan fastq files without assembling them for a custom list of user defined snps?

# Answer
Try this pipeline:

1. Clean raw fastq files (trimmomatic)
2. Align clean fastq files to the reference genome (dna-to-dna aligner of your choice)
3. Convert the alignment sam file to bam (samtools)
4. Call SNPs from bam file (samtools)

Here's a start tutorial on this: http://ged.msu.edu/angus/tutorials-2013/snp\_tutorial.html
> 4 votes

# Answer
You could use ARIBA, which will map you reads to the regions of your interest (more on that in a second) and then assemble/call variants.

The tool is mostly used to find known variants related to AMR (antimicrobial resistance) from a set of databases, but according to the project's wiki you can build your own database with your region of interest.

```
Additionally, reference sequences can be either of:

    Presence/absence sequences. ARIBA will look for these sequences in the input reads and report any that it finds, and also any variants between these sequences and the input reads.

    Variants only sequences. These should have known variant details specified in the metadata file (see below). ARIBA reports only when it finds at least one of the given variants in each of these these sequences. If a sample has one of these sequences, but does not have one of the supplied variants, then it is not reported. If you supply a variants only sequence, but no variant, then the sequence will be removed during sanity checks.

```
> 1 votes

---
Tags: fastq, vcf, snp
---


# Question
Title: Correct for gene length or read counts in GO enrichment analysis
It is a well reported fact that GO analysis of RNAseq results is affected by a number of biases, including length bias and expression level bias.

The `bioconductor` package `goseq` allows you to correct for these biases.

By default it corrects for length bias, but you can also get it to do read count bias. Using read counts to do the correction is attractive because in theory it should account for both sources of bias ($read counts\approx expression \times length$).

I'm doing an enrichment analysis were I have tried both options (length and read counts) and get very different answers. If I run a binomial regression on expression and length vs probability of being differential, I can see that both are *independently* important.

```
> model <- glm(sig ~  expression + log(length), data=retained_genes,  family=binomial(link="logit"))
> print(anova(model, test="Chisq"))

Analysis of Deviance Table

Model: binomial, link: logit

Response: sig

Terms added sequentially (first to last)

                       Df Deviance Resid. Df Resid. Dev  Pr(>Chi)
NULL                                    6676     4507.1
expression              1  114.998      6675     4392.1 < 2.2e-16 ***
log(length)             1  102.553      6674     4289.5 < 2.2e-16 ***
expression:log(length)  1   34.094      6673     4255.4 5.252e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

```

So I'm know unsure what to do, should I use the analysis corrected for length or read count. Or perhaps take only terms significant in both? Or only in one?

# Answer
My understanding of the subject is that the bias of the gene length (and other bias) should be taken care when analyzing the expression and **before** enrichment analysis. The enrichment analysis should be done once the corrections has been performed. Because as the abstract of the GOseq paper states:

> GO analysis is widely used ... but standard methods give biased results on RNA-seq data due to over-detection of differential expression for long and highly expressed transcripts.

So first, take care of the differential expression bias by length and then use the GO to reduce complexity. How you take care of the bias in the RNA-seq data is another question. But the cqn package of Bioconductor can correct the expression by gene length and GC content. However, this correction might hurt the differential tool used (see this discussion in Bioconductor), so it might be better in some cases to use GOSeq.

Now, to the question itself:

> So I'm know unsure what to do, should I use the analysis corrected for length or read count. Or perhaps take only terms significant in both? Or only in one?

Use whatever correction method that yields better differentially expressed genes (DEG. If you find that the correction for length improves the accuracy of the predictions of DEG better than correcting by length and GC, then use that one.

Another option to obtain accurate GO terms, then you could use other testing procedures which don't rely uniquely in the Fisher test, such as those that take into account the structure of the GO graph. TopGO use this approach (note that it is a bit difficult to work with this package),this will reduce the role of the gene length bias (and probably other bias) in the resulting significant GO.
> 2 votes

# Answer
I am not sure if it makes sense to use read counts as bias instead of gene length (and I certainly wouldn't expect the same results).

Do you use total read counts of all your samples (library size)?

The correction for gene length is a pure technical one, the longer a gene the more reads will align (and higher read count genes are easier significant, since they are way above the noise threshold). If you use read counts, you also have a biological factor (expression) in there, which (I think) is the stuff you test with statistics (e.g., with edgeR) and thus not the bias you want to correct for.
> 1 votes

---
Tags: rna-seq, differential-expression, go
---


title: Identity plus Similarity Interpretation

text:
# Question
Title: Identity plus Similarity Interpretation
I am helping a colleague to interpret protein multiple sequence alignment results. My colleague used a third-party proprietary solution (I know...) to calculate the alignment, and one of the resulting tables is a pairwise distance matrix with identity % in the top triangle and **identity plus similarity** in the bottom triangle.

I am struggling to interpret identity plus similarity, despite the fact I found the metric in several papers.

# Answer
> 3 votes
I guess you already know what identity means. But for completeness, it is the % of pairwise alignment position which have an identical amino acid residue. For the identity+similarity it is important to know that amino acids can be grouped according to their physicochemical properties e.g. charge, size etc. While there are standard groupings for amino acids it may be that the proprietary algorithm uses different groupings. Either way, % identity+similarity would imply the number of pairwise alignment positions which contain either the same amino acid or amino acids from the same group (however they define groups).

---
Tags: proteins, sequence-alignment
---


# Question
Title: Is there any difference between SNPs 'AG' and 'GA' in association analyses?
I am using plink to check the association between a phenotype and the SNPs of a gene. plink says the phenotype is significantly associated with a SNP on that gene, and when I check the SNP alleles, for 6 out of 49 samples, the nucleotides are 'GA', and for the rest, they are 'AG', so only the order is different. I am not a biologist, but my basic genetics knowledge says that those two are semantically the same. However, obviously, plink does not agree with me and it associates this SNP with a small p-value since 'GA' samples have a higher value for the phenotype than 'AG' samples. I am very new to SNP analysis, any help is appreciated.

# Answer
Could you please show us the context in which this appears, as you seem to be interpreting this differently to Devon.

If it's appearing as you say, GA and AG, then Devon is right, this usually means that the sequence goes ###AG### or ###GA###, which are two very different sequences.

If, however, as you're implying, the sequences are actually ###G### or ###A###, then semantically, the two ways of writing it have the same outcome. Normally, however, this is written as A/G or G/A.

Alternatively, you're looking at genotypes, which would usually appear as A|G and G|A. In this case, both individuals are heterozygotes, and in terms of the individual SNP they are the same. However the side of the pipe is relevant as this tells you which homologous chromosome it is on. This allows you to calculate LD with respect to neighbouring SNPs, and to determine the actual sequences of the proteins in that individual.
> 8 votes

# Answer
Is "user5054" the same as "user5504"? No? Exactly. Not only does order matter, it's incredibly vitally important. `AG` and `GA` are completely and totally different from each other. If this is in a coding region, then the resulting amino acid is undoubtedly changed (fun fact, the only exception is `TAG` and `TGA`). If it's at a splice site then splicing is likely altered (AG is part of a splice acceptor site). If it's at the binding site of something then it wouldn't be surprised if it didn't bind any more.
> 4 votes

# Answer
Most association analyses are carried out at a single SNP level, so AG and GA are likely to indicate a heterozygous genotype at a particular location. However, the precise notation matters.

As @Emily\_Ensembl has alluded to, for VCF files, A/G indicates unphased SNPs (order unknown, and shouldn't be considered in analyses), whereas A|G indicates phased SNPs (order matters).

It is possible to carry out association analyses by combining multiple variants into a single haplotype, in which case the order would also matter. It's unlikely that this is the case, because such haplotype variant analyses are substantially more complicated to carry out, and will typically be paired with information about SNP linkage (statistic D´ or r²).
> 3 votes

---
Tags: snp, gwas, plink
---


# Question
Title: How to adapt the fgseaL function to perform rapidGSEA computation of gene ranks across 9 different phenotype labels?
I wish to adapt the r language function fgseaL, https://github.com/ctlab/fgsea , to perform rapidGSEA, https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-016-1244-x , computation of inter-class deviation per gene and the subsequent gene rank sorting operation on 9 different phenotype labels as illustrated in the diagram immediately below:

I thought of applying the R-language rank() function on the Expression Data Matrix D. If that is not correct, what sequence of R language commands should we apply to the Expression Data Matrix D to calculate a key value sorted deviation measure across 8 labeled human leukemia groups and a healthy labeled normal control group prior to running fgseaL?

I show below how fgseaL finds the correlation matrix between the R language variable , mat, which corresponds to the Expression Data Matrix D and the R language variable , labels , which is a vector of gene phenotype labels

```
    tmatSc <- scale(t(mat))
    labelsSc <- scale(labels)[, 1]

    minSize <- max(minSize, 1)

    pathwaysFiltered <- lapply(pathways, function(p) { as.vector(na.omit(fmatch(p, rownames(mat)))) })
    pathwaysSizes <- sapply(pathwaysFiltered, length)

    toKeep <- which(minSize <= pathwaysSizes & pathwaysSizes <= maxSize)
    m <- length(toKeep)

    if (m == 0) {
        return(data.table(pathway=character(),
                          pval=numeric(),
                          padj=numeric(),
                          ES=numeric(),
                          NES=numeric(),
                          nMoreExtreme=numeric(),
                          size=integer(),
                          leadingEdge=list()))
    }

    pathwaysFiltered <- pathwaysFiltered[toKeep]
    pathwaysSizes <- pathwaysSizes[toKeep]

    corRanks <- var(tmatSc, labelsSc)[,1]
    ranksOrder <- order(corRanks, decreasing=T)
    ranksOrderInv <- invPerm(ranksOrder)
    stats <- corRanks[ranksOrder]

    pathwaysReordered <- lapply(pathwaysFiltered, function(x) ranksOrderInv[x])

    gseaStatRes <- do.call(rbind,
                           lapply(pathwaysReordered, calcGseaStat,
                                  stats=stats,
                                  returnLeadingEdge=TRUE))

I found a problem with the algorithm shown immediately below.
correcttest <- data.frame(names = row.names(normal))
correcttest <- cbind(correcttest3, normal)
correcttest <- cbind(correcttest3, ALL3m)
rownames(correcttest) <- correcttest$names
correcttest$names <- NULL
correctlabelnormal <- rep(0:0, 73)
correctlabelALL3m <- rep(1:1, 122)
correctlabel <- as.vector(c(correctlabelnormal,correctlabelALL3m))
s <- apply(correcttest, 1, function(x) coef(lm(x~correctlabel))[2])
o <- rank(s)
o <- max(o) - o + 1
res <- fgseaL(df,o,correctlabel,nperm = 2000,minSize = 1, maxSize=50000)
empty data table (0 rows) of 8 columns:   pathway,pval,padj,ES,NES,nMoreExtreme,size

I found the binary phenotype labeled group fgseaL test results below looked satisfactory.
correcttest <- data.frame(names = row.names(normal))
correcttest <- cbind(correcttest3, normal)
correcttest <- cbind(correcttest3, ALL3m)
rownames(correcttest) <- correcttest$names
correcttest$names <- NULL
correctlabelnormal <- rep(0:0, 73)
correctlabelALL3m <- rep(1:1, 122)
correctlabel <- as.vector(c(correctlabelnormal,correctlabelALL3m))
fgseaL(df,correcttest,correctlabel,nperm = 2000,minSize = 1, maxSize=50000)
       pathway        pval        padj         ES       NES nMoreExtreme  size
1: Gene.Symbol 0.003940887 0.003940887 -0.2460126 -1.180009            3 45714
                                 leadingEdge
1: AKIRIN2,LRRC20,HSPA5,HSPA5,DTWD2,ZFYVE28,

```

Thank you for your consideration.

# Answer
> 2 votes
You rank the fit coefficient rather than the original score matrix. So, given a score matrix, `D`:

```
D = matrix(c(22,20,9,8,46,22,18,10,3,18,3,29,2,1,5,45,43,47,17,5,14,44,21,36), byrow=T, ncol=6)
cl = c(0,0,0,1,1,1)
s = apply(m, 1, function(x) coef(lm(x~cl))[2]) # [1]
o = rank(s)
o = max(o) - o + 1 # [2]

```

`o` is then the rank of each row.

\[1\] This fits each row as a linear model of `cl` and extracts the `cl` coefficient.

\[2\] This converts the ranking to be the same as shown in the figure. I don't know if this is important, but I would assume so given how GSEA methods tend to work.

# Answer
> 3 votes
I finally found another answer to my question. Please read this great article in May 12 2017 BioMed Central (BMC) Bioinformatics article titled Ranking metrics in gene set enrichment analysis: do they matter?.

Also, please read this blog , Diving into Genetics and Genomics: Gene Set Enrichment Analysis (GSEA) explained.

After reading these two articles, my choice for the best rapidGSEA local ranking measure is Minimum Significant Difference (i.e., MSD), because it has the best overall false positive rate (i.e., FPR) for larger sample sizes.

Finally, it is important to realize that fgseaL's phenotype labeling can hypothetically be emulated by GSEA with either one of it's sixteen possible ranking metrics or a custom ranking metric.

# Answer
> 0 votes
This weekend, I installed the R language packages rapidGSEA and the Broad Institute GSEA , the Nvidia CUDA Toolkit 6.5-- Custom installation option and MINGW makefiles on my son's Windows 8.1 notebook computer. This step was very tricky to accomplish because I had to maneuver around the "NVIDIA Installer cannot continue. The NVIDIA graphics driver is not compatible with this version of Windows" installation error message which prevents the NVIDIA compilers from being installed correctly. My objective is to run rapidGSEA without the CUDA enhancements since I am not using a computer with a NVIDIA GPU. I did this step because rapidGSEA and Broad Institute GSEA have a significantly richer set of gene set enrichment analytics than fgseaL as shown below:

The GSEA method takes five default arguments and three optional arguments GSEA \<- function(exprsData, labelList, geneSets, numPermutations, metricString, dumpFileName="", checkInput=TRUE, doublePrecision=FALSE) {...}

exprsData, labelList and geneSets refer to the data obtained in the previous section. numPermutations denotes the number of permutations in the resampling test, metricString denotes the local ranking measure (one of the following):

```
naive_diff_of_classes
naive_ratio_of_classes
naive_log2_ratio_of_classes
stable_diff_of_classes
stable_ratio_of_classes
stable_log2_ratio_of_classes
onepass_signal2noise
onepass_t_test
twopass_signal2noise
twopass_t_test
stable_signal2noise
stable_t_test
overkill_signal2noise
overkill_t_test

```

---
Tags: r, genomics, gsea
---


# Question
Title: How to run enrichment analysis of protein functional annotation?
I have a lot of protein clusters. I want to perform an enrichment analysis of their functional annotations, against reference datasets or list of genes I select.

More precisely: a method yields cluster of proteins. I want to decide, for each cluster (which is a set of proteins identifiers), if it holds meaning regarding the proteins it contains, by comparing the annotations found in the cluster and the annotations found in common or specific datasets.

Initially, I used DAVID, which compute the GO annotations from the protein list, then perform the enrichment analysis against common datasets.

However, DAVID suffer of the following drawbacks:

* since I didn't find any other way to use it, I'm limited by the web interface and its limitations (number of genes, url size, number of request/day).
* automation over the web interface is a hacky script.

The team behind seems to offer a way to run DAVID in-house (allowing, hopefully, to relax the limitations), but I can't find any way to download the database from their website.

What are the alternative ways to get enrichment analysis over proteins in a reliable and automated way ?

# Answer
> 3 votes
DAVID depend on a couple of databases from the Gene Ontology Consortium, Reactome, KEGG,... most of them are accessible via Bioconductor. To perform an enrichment analysis you can have a look at the tutorial of the several packages in Bioconductor that do this.

Some of the most important for analyzing enrichment in GO terms are topGO, goseq, GOstats. I would also recommend GOSemSim if you want to compare between GO to focus on a specific GO terms.

Other important packages are the fgsea to test any kind of gene set (which is similar to the one hosted by the Broad Institute), gsva for enrichment analysis by sample, limma has some functions for functional enrichment too. Piano, GSCA, SPIA are also worth mentioning.

Bioconductor has "standard" data sets of the expressions of some cells, like airway and ALL frequently used in vignettes. They are not reference data sets because there isn't a reference expression for a cell of an organisms. It depends on the type of cell, the experiment, the conditions...

# Answer
> 3 votes
I would recommend R instead of sticking to websites.

There are many tools for enrichment analysis. I like to use GOseq, which is initially made for RNAseq data, but can also be used for proteins (I have used it for proteomics). You can use the length of the protein instead of transcript length (to correct for length bias), or you can exclude the length information by just using the hypergeometric method.

See manual for code.

Subsequently you can use gogadget, a tool I made for visualization of GOseq results with a heatmap or cytoscape network.

Additionally, you might wanna take a look at clusterProfiler. With this package you can do enrichment analysis for many clusters in parallel.

---
Tags: go, go-enrichment
---


# Question
Title: Cellranger aggr and Seurat merge difference?
I have multiple libraries of **10x Chromium single-cell RNA-seq** data, which I'd like to combine. One option is using `cellranger aggr` which by default does a depth normalization:

`mapped`: (default) Subsample reads from higher-depth libraries until they all have an equal number of confidently mapped reads per cell. `raw`: Subsample reads from higher-depth libraries until they all have an equal number of total (i.e. raw, mapping-independent) reads per cell. `none`: Do not normalize at all.

https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/aggregate

The new version of `Seurat` (2.0) provides another way via the `MergeSeurat()` (or `AddSamples()`) functions.

http://satijalab.org/seurat/merge\_vignette.html

These have the arguments:

```
do.logNormalize whether to normalize the expression data per cell and transform to log space.
total.expr      scale factor in the log normalization
do.scale        In object@scale.data, perform row-scaling (gene-based z-score)
do.center       In object@scale.data, perform row-centering (gene-based centering)

```

I'm wondering how do these methods compare? What are the pros/cons etc? What do we need to consider when we combine samples?

# Answer
> 3 votes
Cell Ranger aggregate subsamples reads (unless you select `none`), so you will end up with less total reads in samples that have more initially. The output is still raw counts, but you will have more or less per cell.

Seurat just merges the raw counts matrices and normalizes those.

---
Tags: rna-seq, single-cell
---


# Question
Title: What is indel calling and what is its purpose?
I'm having a difficulty in grasping the general purpose and concept of indel calling.

What exactly is this process?

# Answer
> 9 votes
Insertions and deletions (indels) are one type among many different types of *genetic variation*, such as single nucleotide variants (SNVs), copy number variants (CNVs), and structural variants (SVs). I'll assume here that you know how indels are defined, but are simple trying to understand the importance of discovering and analyzing them.

The goal of indel calling, like the goal of any variant calling, is to identify genetic variants that can subsequently be associated with important phenotypes, esp. disease. For example, if 60% of patients with disease XYZ have an indel in the promoter region of gene 123, then that is information of extreme interest and value in research and in clinical care.

Genome-wide association studies (GWAS) have been trying to correlate SNVs to disease and other phenotypes for years. Much less work has been done with indels, but their discovery and analysis remains an area of intense interest.

As far as the *process* of indel calling, large indels can usually be found by mapping paired reads and looking for large discrepancies in the expected distance between pairs and the observed distance.

> Huh, my average insert size is 400bp, but the aligned read pairs flanking this area are 1200bp apart. Must be an 800bp deletion in there!

Smaller indels are much more difficult to detect using this strategy, since they are harder to distinguish from noise (i.e. the variation in length of sequenced fragments). However, as another answer mentions, many short indels are simply reflected as short gaps in the alignment.

# Answer
> 10 votes
There are two types of INDELs: short indels and long indels. Some put the threshold at 50bp; others choose 1000bp. Short and long indels are called differently.

\<50bp short indels are called from read alignment directly. Modern indel callers essentially build a multi-alignment of reads and call an indel with enough read supports. Short indels may break open reading frames.

For Illumina reads, long indels are usually called with split alignment or read pair alignment because mainstream aligners are unable to directly align through a long indel. Long indels are a type of structural variation. They may affect gene structures in a more dramatic way (e.g. delete a whole exon; screw up transcription by transposon insertion; create pseudogenes).

Ultimately, both short and long indels are types of genetic variants. Calling them helps to understand how genetics shapes phenotypes.

---
Tags: genome, variant-calling, indel
---


# Question
Title: Introduce errors in reference transcripts according to external dataset error model
I would like to modify some reference transcripts from Ensembl (*D. melanogaster*) to introduce a controlled rate of random errors in the sequences. The idea would be to introduce random base substitutions in these sequences, no indels for now, because I would like to keep the transcript sequence length as it is in the reference.

The rate of error per transcript will be determined according to an error profile computed from an external set of RNA-seq reads (e.g., generated with ONT MinION)

The aim of this modification would be to establish a rough benchmark of the performances of aligners to use over transcripts from spliced reads (rna-to-genome), aka with more than one exon.

Any idea of which software would be best for this purpose?

# Answer
Do any of the answers for this question help? Karel Brinda has mentioned a few read simulators in the answer to that question, and has a thesis with more information.

Excluding INDEL errors doesn't sound like a good idea; length can still be preserved even if doing that, it just needs an adjustment at the end of the sequence. Note that if you're trying to model nanopore reads, what you're really modelling is the base-caller, rather than the sequencer. I mention this in more detail in my answer.

In most cases where errors are modelled, I find it better to use publicly-available data instead. Especially for nanopore data, there are unmodelled systematic errors in the base-callers and sequencer that can't be simulated using any programs (because they are unmodelled). The following paper would be a good place to start for cDNA sequences, which looks at single-cell data from mouse (C57Bl/6) B1a cells:

http://www.biorxiv.org/content/early/2017/04/13/126847

Illumina and ONT reads for that study can be found in SRA under accession number SRP082530.

I don't know of any recent *D. melanogaster* studies that have been done using nanopore. There's always the option of spending $1000 on a purchase of a MinION with an RNA starter kit to do the study yourself. Here's an older targeted gene study, but bear in mind that it was using an R7.3 flow cell, so errors rates will be much higher than what is currently available:

https://genomebiology.biomedcentral.com/articles/10.1186/s13059-015-0777-z
> 3 votes

# Answer
This preprint uses pbsim to simulate ONT RNA-seq reads for fruit fly. It is probably worth reading if you want to do the same thing.

You should include INDEL errors. Those are what make RNA-seq alignment challenging. For the benchmark purpose, adding INDELs does not increase the complexity at all. You can parse splice junctions on the reference from CIGAR and compare them to the annotation. You don't need to worry about the base-level alignment.

In addition, there are public real ONT data (AC:SRP082530) for SIRV spike-in control and mouse B cells. You don't actually need simulation.

PS: just noticed that you are an author of the first preprint I cited. I would use real data for evaluation.
> 3 votes

# Answer
It sounds like what you're really looking for is a read simulator. A cursory search turns up NanoSim, which is designed to simulate reads from a MinION. This has the benefit of at least having been used in some of the published literature, which is always a nice sign.

You may also find this review article on read simulators useful. It doesn't specifically mention NanoSim, but it should prove to be a useful review of the general concepts anyway if you need to read up on them.
> 1 votes

# Answer
The executable `fastq-sim` in DNemulator package is able to modify a set of input sequences in `fasta` format according to an external set of quality scores reported in a `fastq` file.
> 1 votes

---
Tags: rna-seq, software-recommendation, nanopore, simulated-data
---


# Question
Title: Is there a way to retrieve several SAM fields faster than `samtools view | cut -f`?
I am constructing a bit of software which pipes the outputs of the bam file via `samtools view` into a script for parsing. My goal is to (somehow) make this process more efficient, and faster than `samtools view`.

I am only using 3-4 fields in the bam. Naturally, if I only wanted those fields, I could do something like `samtools view file1.bam | cut -f ##`.

However, I'm trying to be faster than this approach. More specifically:

(1) Given the bam file format, is it theoretically possible to output several requested SAM fields faster than `samtools view|cut -f`?

(2) Is there software available to achieve (1), or would one need to hack samtools?

# Answer
> 6 votes
The BAM file format is not a text-based format. It has a specific binary structure, specified in reasonable detail in the SAM file format specification. Whenever this information is displayed on a screen as text, it needs to be converted from the binary format to a text format, which takes a bit of time and processing power.

As this question suggests, if only a specific field from alignments is needed (or fields), then it will probably be better to extract just those fields and do any necessary conversion only on those fields. While this *can* be done by writing a BAM parser from scratch, many developers have already written software libraries to process BAM files in this way.

Devon Ryan has suggested htslib, which is a C library written by a group that includes the people who wrote the SAM/BAM/CRAM file format specifications. There's also pysam, which is a python wrapper around htslib.

The particular tool that is used will depend on your familiarity with programming and the specific thing that you want to do. If you just want to "output a BAM file", then cat is the quickest:

```
cat file1.bam

```

... but you probably don't want to do that, because it seems like you want to process a text-based representation with a script. Because you haven't specified which fields you are interested in, it's not possible to suggest the best thing to use. In the end, I expect that an "efficient" solution to your problem would involve `htslib` in some form. This is not really hacking *samtools*, it's using the backend of samtools to process BAM data.

However, it's useful to also think about the time cost of coding. What application do you have which means that the text-based processing method is not sufficient? It takes time to write code, and a lot more time to debug that code to make sure it's doing the right thing in all situations. If this is for a one-off thing, then `samtools view` output fed into your script may be the quickest solution.

# Answer
> 6 votes
I modified your original question: as you are extracting 4 fields, you are not outputting BAM. The answer to the modified question is: yes, you can write a C program with htslib (or with bamtools, bioD, bioGo or rust-bio). Formatting an entire SAM is fairly expensive. You can see this by comparing `samtools view aln.bam > /dev/null` and `samtools view -u aln.bam > /dev/null`. With a C program, you can select fields to output. This will give you a noticeable performance boost, depending on the fields you extract.

In general, if you really care about performance, don't use script. The time you save from optimizing `samtools view|cut -f` will be quickly neutralized by the inefficiency of most scripting languages.

# Answer
> 3 votes
Another library you can use for this purpose is the htsjdk, which is written in java. Using htsjdk with java is analogous to using htslib with C; the BAM format is already handled by the library and you can manipulate fields in your own code. The same basic analysis applies to java as C; you do not need to convert the BAM file to text and then parse it. Generally, a program implemented in java will be slower than C but still significantly faster than scripting. You do get language features of java: portability, array bounds checking, garbage collection, etc.

---
Tags: sam, samtools
---


title: VerifyBamID freemix score

text:
# Question
Title: VerifyBamID freemix score
We all have used freemix score at some stage to check contamination or swaps in our sequencing experiment, but can anyone once for all explain how the score is calculated?

# Answer
Here's my current understanding of the freemix metric output by verifyBamID. The freemix score is a measure of contamination that is calculated without knowing the genotypes of the individuals. In the paper it is equation (2), and the authors refer to it as the "sequence-only" method. However, in addition to the sequencing reads, it requires some knowledge of SNPs and their allele frequency from the population. They demonstrate in the paper that using the SNPs and allele frequency estimates from a different population actually reduces the estimated level of contamination.

As for intuition for what equation (2) is measuring, imagine what would happen to the number of heterozygous sites detected if DNA from two different individuals was sequenced together. The total number of heterozygous sites would increase because it would count all sites where the two individuals are homozygous for opposite alleles or at least one of the two individuals is heterozygous. This would lead to an increased fraction of heterozygous sites over that expected from Hardy-Weinberg Equilibrium (i.e. 2pq).

Here's the explanation of freemix directly from the author, Hyun Min Kang:

> The key idea of FREEMIX estimate is to use excessive heterozygosity to estimate the level of contamination. Especially for common SNPs, you will observe higher fraction of heterozygous alleles than 2\*p\*(1-p), and it turns out that you can quantify the contamination very well if you know the population allele frequency already. If you do not have accurate population allele frequency information, than it would be harder to estimate FREEMIX parameters using verifyBamID.

Additional sources:
> 4 votes

---
Tags: sam, ngs
---


# Question
Title: Greedy Motif Search Using Probability Matrices
In my greedy search with pseudocounts algorithm in my bioinformatics course, I did not follow the pseudocode since I wanted to solve the problem in my own way. Unfortunately, although my answer is correct, the answer we are supposed to give is the first correct answer, not the last.

The problem works as such: Given t strings, you make a profile matrix from the first k-mer in the first string, and then you look at the next string and using that profile matrix, you find the best k-mer match, and then you combine that k-mer with the one from the first string to make a new profile matrix, etc. You eventually do this for all k-mer windows in your first string of DNA. The way I did it was by making a probability profile matrix from the best k-mer groups in each case, and then multiplying the greatest probabilities in each column. This gives me a score, and if the next score is higher (thus, more probable) that group of k-mers is used. Below is my code:

```
 def highest_probability(string, k, matrix):
    score = 1
    temp_string = ""
    c = 0
    best_string = ""
    temp = 0
    while c+k < len(string):
        for i in range(c, c+k):
            if string[i] == 'A':
                score *= matrix[i-c]['A']
                temp_string += 'A'
            elif string[i] == 'C':
                score *= matrix[i-c]['C']
                temp_string += 'C'
            elif string[i] == 'G':
                score *= matrix[i-c]['G']
                temp_string += 'G'
            elif string[i] == 'T':
                score *= matrix[i-c]['T']
                temp_string += 'T'
            if i == c+k-1:
                if score >= temp:
                    temp = score
                    best_string = temp_string
        temp_string=""
        score=1
        c+=1
    print(best_string)
    return best_string

def GreedyMotifSearch(DNA, k, t):
    best_motifs = [i[0:k] for i in DNA]
    score=0
    c=0
    while c+k<len(DNA[0]):
        k_mer = DNA[0][c:c+k]
        motif=[k_mer]
        matrix = make_matrix(k_mer, matrix=[], i=0)
        for i in range(1, t):
            compare_string = DNA[i]
            good_motif = highest_probability(compare_string, k, matrix)
            motif.append(good_motif)
            matrix = make_matrix(good_motif, matrix, i)
        temp_score=1
        for dictionary in matrix:
            temp_score *= max(dictionary.values())
        print(temp_score)
        if temp_score > score:
            score = temp_score
            best_motifs = motif
        c+=1
    return best_motifs

def make_matrix(string, matrix, i):
    if len(matrix)==0:
        for ch in string:
            new_dict = {'A':0, 'C':0, 'G':0, 'T':0}
            new_dict[ch] = 1
            matrix += [new_dict]
        print(matrix)
        return matrix
    elif len(matrix)==3:
        for j in range(len(string)):
            matrix[j][string[j]] += 1
        Sum = [sum(matrix[0].values()),sum(matrix[1].values()),sum(matrix[2].values())]
        for columns in range(len(matrix)):
            for keys in matrix[columns]:
                matrix[columns][keys] /= Sum[columns]
        print(matrix)
        return matrix

DNA=['GGCGTTCAGGCA',
'AAGAATCAGTCA',
'CAAGGAGTTCGC',
'CACGTCAATCAC',
'CAATAATATTCG']
k=3
t=5
print(GreedyMotifSearch(DNA, k, t))

```

The output gives me \['TTC', 'GTC', 'TTC', 'GTC', 'TTC'\] which is equally correct to the answer given in the problem output which is TTC,ATC,TTC,ATC,TTC. Unfortunately, my answer is the last correct occurrence. I tried adjusting my highest\_probability function with if score is greater than temp instead of greater or equal, but that gives me the answer \['GGC', '', '', '', ''\] which I don't understand at all. My question is 1) Why is this the case? and 2) How can I fix this? And does anyone have a simpler way of doing this?

# Answer
Assuming you're talking about Implement GreedyMotifSearch from Rosalind, while I haven't actually finished this particular one, I've got an idea why `TTC,ATC,TTC,ATC,TTC` doesn't work. Here's my explanation from the "Questions" section of that problem:

> The profile is created from the first kmers in each string, and the most likely string from those is CAG. The answers will generally be similar to this.

I've been trying to do Rosalind by only using the information provided in the questions; this is one of the many problems where there wasn't enough information in the Question for me to arrive at a correct solution when I tried it. It's possible that the debug datasets were added to this question since I last looked at it, as they seem to describe this situation:

## TEST DATASET 1:

## Input:

```
3 4
GCCCAA
GGCCTG
AACCTA
TTCCTT

```

## Output:

```
GCC
GCC
AAC
TTC

```

> This dataset checks that your code always picks the first-occurring Profile-most Probable k-mer in a given sequence of Dna. In the first sequence (“GCCCAA”), “GCC” and “CCA” are both Profile-most Probable k-mers. However, you must return “GCC” since it occurs earlier than “CCA”. Thus, if the first sequence of your output is “CCA”, this test case fails your code.
> 1 votes

---
Tags: python
---


# Question
Title: How to simulate "base error rate" in art_illumina?
I'd like to simulate 10% sequencing error using `art_illumina`. The simulator doesn't have a parameter that I can just give the 10%, but it has this:

```
-qs  --qShift   the amount to shift every first-read quality score by
-qs2 --qShift2  the amount to shift every second-read quality score by
                NOTE: For -qs/-qs2 option, a positive number will shift up
                quality scores (the max is 93) that reduce substitution sequencing
                errors and a negative number will shift down quality scores that
                increase sequencing errors. If shifting scores by x, the error
                rate will be 1/(10^(x/10)) of the default profile.

```

Heng Li's `wgsim` has an option for base error rate. Can I do the same for `art_illumina`?

**Q:** Does the "error rate" in the description mean sequencing error? If I set both `qShift` and `qShift2` to 10, `1/(10^(x/10)) == 0.1, where x is 10`, does that mean 10 is the value I should give to the simulator?

# Answer
> 4 votes
You need to specify an error profile to do this. However, you can shift quality scores of an existing profiles to get 10% error, which is easier than defining a complete profile.

For example, Illumina HiSeq 1000 has the error rate of 0.1% and you want 100 time more errors, so you should shift down the quality scores by 20 by setting both `qShift` and `qShift2` to -20. The following code will simulate paired end reads with 10% error.

```
art_illumina -ss HS10 -i reference.fa -qs -20 -qs2 -20 -sam -p -l 100 -f 20 -m 200 -s 10 -o output

```

---
Tags: ngs, simulated-data
---


title: Ensembl id to GeneSymbol with biomart

text:
# Question
Title: Ensembl id to GeneSymbol with biomart
I have 3224 Ensembl id's as rownames in a dataframe "G". To convert Ensembl ids into Genesymbols I used biomart like following.

```
library('biomaRt')
mart <- useDataset("hsapiens_gene_ensembl", useMart("ensembl"))
genes <- rownames(G)
G <-G[,-6]
G_list <- getBM(filters= "ensembl_gene_id", attributes= c("ensembl_gene_id","hgnc_symbol"),values=genes,mart= mart)

```

Now in G\_list I can see only 3200 ensembl ids showing Genesymbols / No Gene\_symbols. Why the other 24 ensembl ids are not seen in G\_list? If there are no gene\_symbol for those 24 ensembl ids it should atleast show "-"

Examples of problematic IDs are: ENSG00000257061, ENSG00000255778, ENSG00000267268. These are not at all shown in G\_list (biomaRt). So, I gave them in biodbnet, which seems to handle them.

what is the problem here?

# Answer
> 2 votes
I used gProfileR where there are no problematic IDs with that. It converted ENSEMBL IDs to Gene\_symbols and made my work easier with GO analysis too. There is an R package and also API.

# Answer
> 6 votes
It looks like you were using an old annotation. The problematic IDs you posted existed in the GRCh37 annotations, but don't in the most recent GRCh38 annotation. For that reason they were excluded. The IDs that have `-` as symbols don't have associated symbols, but are present in the database.

To use an archived version in biomart:

```
mart = useDataset("hsapiens_gene_ensembl", useEnsembl(biomart="ensembl", version=84))

```

That's an example for release 84.

---
Tags: r, bioconductor, biomart, ensembl
---


# Question
Title: What format is this? Pretty sure it's not a BED file
Time for guess the bioinformatics file format... :)

I found this '.bed' file on GEO (GSE84660), but it's clearly not a BED file. Anyone know what it *might* be? And what can view it? It's from a HiC expt.

```
#column headers: bait1_chr, bait1_start, bait1_end, bait2_chr, bait2_start, bait2_end, contact_ID, NA, NA, B1_D0_counts, B1_D3_counts, B1_D6_counts, B2_D0_counts, B2_D3_counts, B2_D6_counts

chr1    831895  848168  chr1    850619  874081  b2b_1   .       .       114     125     127     75      97      92
chr1    831895  848168  chr1    889424  903640  b2b_2   .       .       15      12      16      4       15      20
chr1    831895  848168  chr1    903641  927394  b2b_3   .       .       13      13      10      18      15      8
...

```

# Answer
It looks like a weird way to represent data. You're right, it does not at all look like a BED.

Digging in GEO made me find this information:

> Supplementary\_files\_format\_and\_content: Tab-delimited file reports one contact per row. Interacting HindIII fragments are represented in bed format (chr/start/stop) with bait (or upstream bait, in the case of b2b contacts) listed first. The following fields correspond to a unique contact ID, fragment IDs for first and second fragments, and raw read counts supporting contacts in each replicates.

So I think the authors merged a BED3 format plus additional BED-unrelated fields generating a tab-separated file - not clearly stated in their methods though. The easiest would be to contact the authors maybe?

Depending on your aims, if you only need the first 3 fields (but I doubt it) I would reformat this "bed" file (to make it processable as a BED file) in two sub-files:

```
cat original.bed | cut -f 1,2,3 > bait1.txt
cat bait1.txt | head -3
chr1    831895  848168
chr1    831895  848168
chr1    831895  848168

cat original.bed | cut -f 4,5,6 > bait2.txt
cat bait2.txt | head -3
chr1    850619  874081
chr1    889424  903640
chr1    903641  927394

```

Not sure if it would fit your aims though.

Also, for visualization of contact matrices, maybe you'll find this paper useful.
> 6 votes

---
Tags: file-formats
---


# Question
Title: Is it wise to use RepeatMasker on prokaryotes?
I'm looking for a way to identify low complexity regions and other repeats in the genome of *Escherichia coli*. I found that RepeatMasker may be used for example when drafting genomes of prokaryotes (*E. coli* example). But RepeatMasker works on a limited dataset of species, neither of them being prokaryotes. By default, when running RepeatMasker, if no species is specified, it will compare with homo sapiens data.

This seems rather inadequate, but the most relevent alternative, PRAP, requires a "dead" tool (VisCoSe, by Michael Spitzer).

1. Is it still wise to to use RepeatMasker on *Escherichia coli*?
2. If yes, which settings would maximise relevance ?

# Answer
> 5 votes
If I understood correctly your question, you want to mask those regions in a (FASTA?) genome. I think you could identify those regions using mummer and mask them using bedtools.

```
# align genome against itself
nucmer --maxmatch --nosimplify genome.fasta genome.fasta

# select repeats and convert the corrdinates to bed format
show-coords -r -T -H out.delta | awk '{if ($1 != $3 && $2 != $4) print $0}' | awk '{print $8"\t"$1"\t"$2}' > repeats.bed

# mask those bases with bedtools
bedtools maskfasta -fi genome.fasta -bed repeats.bed -fo masked.fasta

```

Have a look at nucmer and bedtools maskfasta options to fine-tune your analysis.

# Answer
> 1 votes
From your comment, it seems that masking regions are not your priority, but you would rather find them (correct me if I am wrong):

> not interested in masking, but really finding any type of repeats, from low complexity regions to small repeats

To find these regions, you can try RepeatFinder. From their paper, it seems that is suitable for bacterial genomes too. It also seems to be faster because based on suffix tree data structure instead of working with string-matching algorithm (as in RepeatMasker).

---
Tags: genome, repeat-elements, sequence-analysis
---


title: How to convert DNAbin to FASTA in R?

text:
# Question
Title: How to convert DNAbin to FASTA in R?
I am trying to convert DNAbin files to fasta format. Rationale: I want to use the fasta file to calculate non-synonymous/synonymous mutation rate.

`my_dnabin1` is a DNAbin file of 55 samples and I am using the following code to convert it into a fasta file.

`dnabin_to_fasta <- lapply(my_dnabin1, function(x) as.character(x[1:length(x)]))` This generates a list of 55 samples which looks like:

$SS.11.01

\[1\] "t" "t" "a" "c" "c" "t" "a" "a" "a" "a" "a" "g" "c" "c" "g" "c" "t" "t" "c" "c" "c" "t" "c" "c" "a" "a" \[27\] "c" "c" "c" "t" "a" "g" "a" "a" "g" "c" "a" "a" "a" "c" "c" "t" "t" "t" "c" "a" "a" "c" "c" "c" "c" "a"

$SS.11.02

\[1\] "t" "t" "a" "c" "c" "t" "a" "a" "a" "a" "a" "g" "c" "c" "g" "c" "t" "t" "c" "c" "c" "t" "c" "c" "a" "a" \[27\] "c" "c" "c" "t" "a" "g" "a" "a" "g" "c" "a" "a" "a" "c" "c" "t" "t" "t" "c" "a" "a" "c" "c" "c" "c" "a"

and so on...

However, I want a fasta formatted file as the output that may look something like:

`>SS.11.01 ttacctga`

`>SS.11.02 ttacctga`

# Answer
The following should work:

```
cat(file="out.fa", paste(paste0(">",names(dnabin_to_fasta)),
          sapply(dnabin_to_fasta, paste, collapse=""), sep="\n"), sep="\n");

```

A bit of explanation, starting from the middle outwards:

1. Collapse the letter arrays inside dnabin\_to\_fasta into single strings
2. Append "\>" to the front of the names of each sequence
3. Concatenate the results of steps 1 & 2 together, separating with a line break
4. Write each concatenated string to a file, separating by a line break
> 5 votes

---
Tags: r, fasta
---