# Question
Title: Detecting portions of human proteins with high degree of microbial similarity
I'm a newcomer to the world of bioinformatics, and in need of help solving a problem.

My goal is to take a list of human proteins, and identify segments (13-17aa in length) with a high degree of similarity to microbial sequences. Ideally, I would like to start with list of FASTA sequences, and have an easy way to generate an output of the corresponding high similarity segments of each protein.

Are there existing tools or software that I should be aware of that will make my life easier?

Thanks in advance.

# Answer
> 4 votes
Sounds like precisely the job BLAST was developed for. Now, which flavor will depend on what you want to do and what data you have available. Some options:

1. **PSI-BLAST**: this is usually the best choice if you are trying to find protein homologs. It works by building a hidden markov model describing your query sequence and using that model to query a database of proteins. The advantage is that it is run in multiple iterations, giving you the chance to add or remove results (so you add the ones that are true positives and remove false ones), eventually building a pretty good model of your protein. This is far moer powerful than a simple homology-based approach since proteins work via protein domains and simple homology is not as important as specific conserved functional residues.

For this, go to the NCBI protein blast page and select PSI-BLAST:
2. **BLASTp**: Simple protein-protein blast. It will identify homologous proteins based on sequence similarity. Whether or not that also implies functional homology is not that simple and will depend on each case you investigate.

As above, go to the NCBI protein blast page, but this time use the defaults.
3. **tBLASTn**: this is a tool that takes protein sequences as input and compares them to a database of DNA which is dynamically translated in all 6 possible reading frames. Very good for finding homologous sequences when you don't have well annotated protein information for the target species. It has the benefit of being more sensitive and able to find more distant homologies than basic nucleotide BLASTn and the go-to approach when your target species is distant and not well annotated.

NCBI's tBLASTn page.

All of these can be run online through the NCBI's BLAST page. If you want to investigate hundreds of proteins, I suggest you install blast locally. You can then either download the relevant target sequences from NCBI and rebuild the blast database locally (if so, I suggest you ask a new question about how to do that) or, use NCBI's remote blast client which lets you use a locally stored query file and will run blast on the NCBI's servers.

Now, these programs will return what are known as High Scoring Pairs (HSPs), the regions of your query sequence(s) that align well to the target. There are various options you can play with to improve the sensitivity or the specificity but a discussion of those would require far more details about what you're doing and would be best in a new question as well.

Once you have your HSPs, you can relatively easily parse them to select regions with a given range of sequence similarity values and of a specific length. Once again, that would be better discussed in a separate question once you have your results and can show an example.

---
Tags: fasta, sequence-homology
---


# Question
Title: How to extract RNA sequence and secondary structure restrains from a PDB file
I'm trying to find a programmatic way to automatically extract the following information from a PDB file:

1. RNA sequence
2. Secondary structure restraints in bracket format, e.g. `. (( . ( . ) . ))`

Does software exist that can take a PDB file as input and generate these two pieces of information?

e.g. file. 3NDB\_ba.pdb

The nucleotide sequence is:

`gUCUCGUCCCGUGGGGCUCGGCGGUGGGGGAGCAUCUCCUGUAGGGGAGAUGUAACCCCCUUUACCUGCCGAACCCCGCCAGGCCCGGAAGGGAGCAACGGUAGGCAGGACGUCGGCGCUCACGGGGGUGCGGGAC`.

And the secondary structure :

`.(((.(..(((((((((.(((((..(((((.(((((((((....)))))))))..)))))....((((((.....(((.....(((....))).....)))..)))))).))))))).))))))).....).))).`

# Answer
> 8 votes
I suggest you take a look at `rna-pdb-tools` we do way more than you need! :-) The tools can get you a sequence, secondary structure and much more using various algorithms, and all is well documented http://rna-pdb-tools.readthedocs.io/en/latest/

To get sequence http://rna-pdb-tools.readthedocs.io/en/latest/main.html#get-sequence

```
$ rna_pdb_tools.py --get_seq 5_solution_1.pdb
> 5_solution_1.pdb A:1-576
CAUCCGGUAUCCCAAGACAAUCUCGGGUUGGGUUGGGAAGUAUCAUGGCUAAUCACCAUGAUGCAAUCGGGUUGAACACUUAAUUGGGUUAAAACGGUGGGGGACGAUCCCGUAACAUCCGUCCUAACGGCGACAGACUGCACGGCCCUGCCUCAGGUGUGUCCAAUGAACAGUCGUUCCGAAAGGAAG

```

or you can get sequence and secondary structure via x3dna (http://x3dna.org/)

```
[mm] py3dna$ git:(master) ✗ ./rna_x3dna.py test_data/1xjr.pdb
test_data/1xjr.pdb
 >1xjr nts=47 [1xjr] -- secondary structure derived by DSSR
 gGAGUUCACCGAGGCCACGCGGAGUACGAUCGAGGGUACAGUGAAUU
 ..(((((((...((((.((((.....))..))..))).).)))))))

```

The problem is not trivial. I hope rna-pdb-tools has stuff that can solve this issue for you (http://rna-pdb-tools.readthedocs.io/en/latest/want.html#module-rna\_pdb\_tools.utils.rna\_x3dna.rna\_x3dna)

# Answer
> 4 votes
## Context

The **PDB file format** is a fixed-column file format designed in 1970s for storing structural models of macromolecules. The format has been around for long time, has many uses, and although it has official spec the files in circulation may not strictly conform to it. It always has a list of atoms with coordinates (the first two lines are added to stress that it's a fixed-column format, they are not part of the file):

```
         1         2         3         4         5         6         7         8
12345678901234567890123456789012345678901234567890123456789012345678901234567890
ATOM      1  N   VAL A   1       3.320  14.780   4.844  1.00 35.53           N
ATOM      2  CA  VAL A   1       3.577  16.239   4.984  1.00 35.39           C
ATOM      3  C   VAL A   1       4.896  16.398   5.727  1.00 30.43           C
ATOM      4  O   VAL A   1       5.143  15.702   6.732  1.00 30.51           O
ATOM      5  CB  VAL A   1       2.343  16.975   5.494  1.00 45.39           C
ATOM      6  CG1 VAL A   1       2.586  18.497   5.590  1.00 60.06           C
ATOM      7  CG2 VAL A   1       1.103  16.811   4.634  1.00 53.93           C
ATOM      8  N   LEU A   2       5.748  17.241   5.158  1.00 28.04           N
ATOM      9  CA  LEU A   2       7.116  17.471   5.661  1.00 24.31           C
ATOM     10  C   LEU A   2       7.166  18.490   6.792  1.00 24.00           C
...

```

or, to show some RNA:

```
ATOM     42  N3    G B   2       9.252  12.871  -1.168  1.00 36.98           N
ATOM     43  C4    G B   2       8.424  12.964  -2.233  1.00 34.09           C
ATOM     44  P     A B   3      10.376   8.321  -4.834  1.00 40.53           P
ATOM     45  OP1   A B   3      11.773   8.279  -5.364  1.00 38.97           O
ATOM     46  OP2   A B   3       9.396   7.283  -5.218  1.00 39.32           O
ATOM     47  O5'   A B   3      10.429   8.473  -3.211  1.00 36.55           O
ATOM     48  C5'   A B   3      11.698   8.232  -2.554  1.00 35.13           C

```

The **Protein Data Bank** -- international institution that archives and annotates structural models of biological molecules that anyone can deposit -- is now using **mmCIF** as the primary format. The fixed-column PDB format had inherent limitations (max. 99,999 atoms, but also it was hard to include additional info). They still generate PDB files if possible (i.e. except the largest structures), but the only format accepted now for depositions is mmCIF.

mmCIF has a rather not intuitive CIF syntax (think JSON or XML, but designed before XML), plus ontology defined by PDBx/mmCIF dictionary (think XML Schema). The content of mmCIF is divided into table-like categories with relations between tables (it was designed at the peak of RDMBS popularity) and is a bit harder to work with, so the old PDB format is more popular.

## Sequence from PDB

The PDB format has a record called `SEQRES` that explicitly lists the sequence, for example (5NEO):

```
SEQRES   1 A   18    G   G   U   G   G   G   G   A   C   G   A   C   C
SEQRES   2 A   18    C   C   A CBV   C

```

The mmCIF format has more details:

```
_entity_poly.pdbx_seq_one_letter_code       'GGUGGGGACGACCCCA(CBV)C'
_entity_poly.pdbx_seq_one_letter_code_can   GGUGGGGACGACCCCACC
_entity_poly.pdbx_strand_id                 A

loop_
_entity_poly_seq.entity_id
_entity_poly_seq.num
_entity_poly_seq.mon_id
_entity_poly_seq.hetero
1 1  G   n
1 2  G   n
1 3  U   n
1 4  G   n
1 5  G   n
1 6  G   n
1 7  G   n
1 8  A   n
1 9  C   n
1 10 G   n
1 11 A   n
1 12 C   n
1 13 C   n
1 14 C   n
1 15 C   n
1 16 A   n
1 17 CBV n
1 18 C   n

```

The sequence is conveniently extracted for you.

Additionally, `_entity_poly_seq` in mmCIF includes information about microheterogeneity, and the SEQRES record doesn't. It is relevant when the model has two alternative residues in the same place, because part of the sample had a point mutation.

PDB files that didn't come from PDB usually don't have the SEQRES record. You may get the sequence from the list of atoms (residue names in columns 18-20), but some residues may be missing in the atom list if atomic positions could not be determined.

So it should be easy to extract the sequence either way. For example, BioPython has module BioSeqIO with two PDB pseudo-formats:

* pdb-seqres - Reads a Protein Data Bank (PDB) file to determine the complete protein sequence as it appears in the header (no dependencies).
* pdb-atom - Uses Bio.PDB to determine the (partial) protein sequence as it appears in the structure based on the atom coordinate section of the file (requires NumPy for Bio.PDB).

*I just tried BioPython 1.66 and while `SeqIO.parse()` can extract protein sequence, it fails with RNA, such as the 5NEO above. Ooops.*

## Secondary structure restraints

We may be thinking about different restraints here. I'll write about restraints used in refinement. The experimental data alone is normally not sufficient to refine a model, so one needs restraints that represent prior knowledge - lengths of atomic bonds, angles, planarity restraints, restraints based on local similarity to other structures that were determined from higher resolution data, etc, whatever can help to make a sensible model. Refinement programs (such as Refmac, BUSTER, Phenix.refine) try to fit the model to the data and satisfy geometrical restraints at the same time.

Macromolecular crystallographic software is quite fragmented and one normally uses several different programs in the process. It happens that the programs to generate secondary structure restraints are separate from the actual refinement programs. (Secondary structure restraints are less essential than restraints for covalent bonds). CCP4 has a program called LibG that makes DNA/RNA restraints for Refmac. Phenix has a program called phenix.secondary\_structure\_restraints, or you can use a server from UCSC.

These restraints don't have a "bracket format", but they explicitly specify expected distances and angles between atoms. For example (for phenix.refine):

```
  bond {
    action = *add
    atom_selection_1 = chain  A and resid    1  and name  O6
    atom_selection_2 = chain  A and resid   18  and name  N4
    distance_ideal = 2.91
    sigma = 0.1
  }

```

or (for Refmac):

```
exte dist first chain A resi 16 ins  . atom N6 second chain A resi 3 ins . atom O4 value 2.94 sigma 0.15 type 1
exte dist first chain A resi 16 ins  . atom N1 second chain A resi 3 ins . atom N3 value 2.84 sigma 0.1 type 1
exte torsion first chain A resi 16 ins  . atom C2 second chain A resi 16 ins . atom N1 third chain A resi 3 ins . atom N3 fourth chain A resi 3 ins . atom C4 value 180 sigma 15 type 1

```

# Answer
> 1 votes
If you know how to use python:

1 - Download the ModeRNA module from here and install

2 - from a python IDLE execute:

```
from moderna import *
m = load_model('file.pdb', 'A') #A is the chain
seq = m.get_sequence()
sec = m.get_secstruc()

```

the variable `sec` store the secondary structure in dot-bracket notation.

---
Tags: public-databases, pdb, rna, rna-structure
---


# Question
Title: How to read and interpret a gene expression quantification file?
I have a gene expression quantification file from TCGA that contains the following lines:

```
ENSG00000242268.2   591.041000514
ENSG00000270112.3   0.0
ENSG00000167578.15  62780.6543066
ENSG00000273842.1   0.0
ENSG00000078237.5   36230.832883
ENSG00000146083.10  189653.152706
ENSG00000225275.4   0.0
ENSG00000158486.12  420.761140072
ENSG00000198242.12  2914738.3675
ENSG00000259883.1   1632.83700531
ENSG00000231981.3   0.0
ENSG00000269475.2   0.0
ENSG00000201788.1   0.0
ENSG00000134108.11  925529.547944
ENSG00000263089.1   2646.63769677
ENSG00000172137.17  23162.6989867
ENSG00000167700.7   291192.25157

```
1. What is the `.number` that are added to the Gene e.g. the `.2` in `ENSG00000242268.2`
2. Why isn't the quantified value an integer; what does `591.041000514` mean?

In case someone one wants more info about the gene LINC02082.

# Answer
> 9 votes
1. The first column contains Ensembl gene identifiers, and the suffix is a version number that can be used to track changes to the gene annotations over time. From the **Ensembl Stable IDs** documentation:

> Ensembl annotation uses a system of stable IDs that have prefixes based on the species scientific name plus the feature type, followed by a series of digits and a version e.g. ENSG00000139618.1. The version may be omitted.
2. Following the first link you provided leads to a page with details for the file `2edcaaa7-63b4-40b4-abbe-5d7a84012e60.FPKM-UQ.txt.gz`. The first thing that caught my eye about this filename was `FPKM`, or "fragments per kilobase of exon per million reads", which is a commonly used unit of RNA expression. Since these are not raw read counts, there's no expectation that these values should be integers.

The best explanation I've seen of FPKM comes from a blog post written by Harold Pimentel of kallisto and sleuth fame. From the blog post:

> The interpretation of FPKM is as follows: if you were to sequence this pool of RNA again, you expect to see `FPKM_i` fragments for each thousand bases in the feature for every `N/10^6` fragments you’ve sequenced. It’s basically just the rate of fragments per base multiplied by a big number (proportional to the number of fragments you sequenced) to make it more convenient.

---

More generally though, even when FPKM is not the unit of expression abundance used, most quantification methods and associated units of expression will not produce integer estimates.

---
Tags: gene, normalization, identifiers, quantification
---


title: Working with old genome builds

text:
# Question
Title: Working with old genome builds
Is working with and relying on old genome builds still valid?

For example NCBI36/hg18. Would results from papers based on old builds require LiftOver and re-analysis to be useful?

A bit of context, this is related to other post, where we have aCGH results based on old build: How do I validate a single sample ArrayCGH result?

# Answer
> 8 votes
In my opinion, it is not very reliable. LiftOver is very limited in terms of transformations it can support. The LiftOver Chain format can capture only matching regions in the same order. It means that it can account for indels, but even simple structural variations become problematic.

For instance, when a newer assembly is available, it is usually a recommended practice to remap all the reads rather than transform the existing alignments.

# Answer
> 4 votes
I think that right now, the only human builds that are worth considering are hg19/GRCh37 as many data bases such as gnomAD still exclusively use this release. On the other hand, hg38/GRCh8 has many important fixes and the useful (but yet underused) feature of alternative loci.

Anything from older releases should be remapped to a more recent one.

# Answer
> 2 votes
You could use liftOver which isn't always great.

Whenever I encounter this (especially NGS data readily available on the SRA), I often just get the raw files (e.g. fastqs) and re-align/re-map.

In your case (arrays) it may be a bit tough. Not impossible though, as I recently took some old yeast DNA/RNA microarray data and updated it to the newest genome. Just requires the right data (like DNA for normalization) and a good understanding of the entire process.

A last resort/alternative is to align your **new** data to the old genome to be able to make comparisons. This isn't ideal but works in cases where upgrading one source isn't possible or is a HUGE amount of time/effort. I've done this for a few fly experiments where all the available/previous data was done in dm3. All the old genomes can generally be found on http://archive.ensembl.org .

# Answer
> 0 votes
For mouse, I still see people using mm9/NCBI37 in high profile publications even though mm10/GRCm38 was released more than 5 years ago (2011). I personally don't think that's a great idea, but it's certainly valid according to the peer reviewers.

It also depends on your application. If you are working with coding regions (likely well known for a long time) or extracting genome-wide stats (enrichment at TSS, for example), the differences should be negligible.

---
Tags: human-genome, array-cgh, liftover
---


# Question
Title: How can I create my own GO association file (gaf)?
This question is based on a question on BioStars posted \>2 years ago by user jack.

It describes a very frequent problem of generating GO annotations for non-model organisms. While it is based on some specific format and single application (Ontologizer), it would be useful to have a general description of the pathway to getting to a GAF file.

Note, that the input format is lacking a bit of essential information, like how it was obtained. Therefore, it is har to assign evidence code. Therefore, lets assume that the assignments of GO terms were done automagically.

---

> I want to do the Gene enrichment using Ontologizer without a predefined association file(it's not model organism).
>
> I have parsed a file with two columns for that organism like this :
>
> ```
> geneA  GO:0006950,GO:0005737
> geneB  GO:0016020,GO:0005524,GO:0006468,GO:0005737,GO:0004674,GO:0006914,GO:0016021,GO:0015031
> geneC  GO:0003779,GO:0006941,GO:0005524,GO:0003774,GO:0005516,GO:0005737,GO:0005863
> geneD  GO:0005634,GO:0003677,GO:0030154,GO:0006350,GO:0006355,GO:0007275,GO:0030528
>
> ```
>
> I have downloaded the .ob file from Gene ontology file which contain this information (from here) :
>
> ```
> !
> ! GO IDs (primary only) and name text strings
> ! GO:0000000 [tab] text string [tab] F|P|C
> ! where F = molecular function, P = biological process, C = cellular component
> !
> GO:0000001  mitochondrion inheritance   P
> GO:0000002  mitochondrial genome maintenance    P
> GO:0000003  reproduction    P
> GO:0000005  ribosomal chaperone activity    F
> GO:0000006  high affinity zinc uptake transmembrane transporter activity    F
> GO:0000007  low-affinity zinc ion transmembrane transporter activity    F
> GO:0000008  thioredoxin F
> GO:0000009  alpha-1,6-mannosyltransferase activity  F
> GO:0000010  trans-hexaprenyltranstransferase activity   F
> GO:0000011  vacuole inheritance P
>
> ```
>
> What I need as output is .gaf file in the following format (in the format of the files here):
>
> ```
> !gaf-version: 2.0
>
> !Project_name: Leishmania major GeneDB
>
> !URL: http://www.genedb.org/leish
>
> !Contact Email: mb4@sanger.ac.uk
>
>  GeneDB_Lmajor    LmjF.36.4770    LmjF.36.4770        GO:0003723    PMID:22396527    ISO    GeneDB:Tb927.10.10130    F    mitochondrial RNA binding complex 1 subunit, putative    LmjF36.4770    gene    taxon:347515    20120910    GeneDB_Lmajor
>  GeneDB_Lmajor    LmjF.36.4770    LmjF.36.4770        GO:0044429    PMID:20660476    ISS        C    mitochondrial RNA binding complex 1 subunit, putative    LmjF36.4770    gene    taxon:347515    20100803 GeneDB_Lmajor             GeneDB_Lmajor    LmjF.36.4770    LmjF.36.4770        GO:0016554    PMID:22396527    ISO    GeneDB:Tb927.10.10130    P    mitochondrial RNA binding complex 1 subunit, putative    LmjF36.4770    gene   taxon:347515    20120910    GeneDB_Lmajor
>  GeneDB_Lmajor    LmjF.36.4770    LmjF.36.4770        GO:0048255    PMID:22396527    ISO    GeneDB:Tb927.10.10130    P    mitochondrial RNA binding complex 1 subunit, putative    LmjF36.4770    gene    taxon:347515    20120910    GeneDB_Lmajor
>
> ```
>
> How to create your own GO association file (gaf)?

# Answer
> 8 votes
Here's a Perl script that can do this:

```
#!/usr/bin/env perl
use strict;
use warnings;

## Change this to whatever taxon you are working with
my $taxon = 'taxon:1000';
chomp(my $date = `date +%Y%M%d`);

my (%aspect, %gos);
## Read the GO.terms_and_ids file to get the aspect (sub ontology)
## of each GO term.
open(my $fh, $ARGV[0]) or die "Need a GO.terms_and_ids file as 1st arg: $!\n";
while (<$fh>) {
    next if /^!/;
    chomp;
    my @fields = split(/\t/);
    ## $aspect{GO:0000001} = 'P'
    $aspect{$fields[0]} = $fields[2];
}
close($fh);

## Read the list of gene annotations
open($fh, $ARGV[1]) or die "Need a list of gene annotattions as 2nd arg: $!\n";
while (<$fh>) {
    chomp;
    my ($gene, @terms) = split(/[\s,]+/);
    ## $gos{geneA} = (go1, go2 ... goN)
    $gos{$gene} = [ @terms ];
}
close($fh);

foreach my $gene (keys(%gos)) {
    foreach my $term (@{$gos{$gene}}) {
        ## Warn and skip if there is no aspect for this term
        if (!$aspect{$term}) {
            print STDERR "Unknown GO term ($term) for gene $gene\n";
            next;
        }
        ## Build a pseudo GAF line
        my @out = ('DB', $gene, $gene, ' ', $term, 'PMID:foo', 'TAS', ' ', $aspect{$term},
                             $gene, ' ', 'protein', $taxon, $date, 'DB', ' ', ' ');
        print join("\t", @out). "\n";
    }
}

```

Make it executable and run it with the `GO.terms_and_ids` file as the 1st argument and the list of gene annotations as the second. Using the current `GO.terms_and_ids` and the example annotations in the question, I get:

```
$ foo.pl GO.terms_and_ids file.gos
DB  geneD   geneD       GO:0005634  PMID:foo    TAS     C   geneD       protein taxon:1000  20170308    DB
DB  geneD   geneD       GO:0003677  PMID:foo    TAS     F   geneD       protein taxon:1000  20170308    DB
DB  geneD   geneD       GO:0030154  PMID:foo    TAS     P   geneD       protein taxon:1000  20170308    DB
Unknown GO term (GO:0006350) for gene geneD
DB  geneD   geneD       GO:0006355  PMID:foo    TAS     P   geneD       protein taxon:1000  20170308    DB
DB  geneD   geneD       GO:0007275  PMID:foo    TAS     P   geneD       protein taxon:1000  20170308    DB
DB  geneD   geneD       GO:0030528  PMID:foo    TAS     F   geneD       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0016020  PMID:foo    TAS     C   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0005524  PMID:foo    TAS     F   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0006468  PMID:foo    TAS     P   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0005737  PMID:foo    TAS     C   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0004674  PMID:foo    TAS     F   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0006914  PMID:foo    TAS     P   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0016021  PMID:foo    TAS     C   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0015031  PMID:foo    TAS     P   geneB       protein taxon:1000  20170308    DB
DB  geneA   geneA       GO:0006950  PMID:foo    TAS     P   geneA       protein taxon:1000  20170308    DB
DB  geneA   geneA       GO:0005737  PMID:foo    TAS     C   geneA       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0003779  PMID:foo    TAS     F   geneC       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0006941  PMID:foo    TAS     P   geneC       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0005524  PMID:foo    TAS     F   geneC       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0003774  PMID:foo    TAS     F   geneC       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0005516  PMID:foo    TAS     F   geneC       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0005737  PMID:foo    TAS     C   geneC       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0005863  PMID:foo    TAS     C   geneC       protein taxon:1000  20170308    DB

```

Note that this is very much a pseudo-GAF file since most of the fields apart from the gene name, GO term and sub-ontology are fake. It should still work for what you need, however.

---
Tags: sequence-annotation, gene-ontology
---


title: What to use to edit RNA alignments?

text:
# Question
Title: What to use to edit RNA alignments?
I have many alignments from Rfam Database, and I would like to edit them. I saw that many tools are used for Protein sequence alignments, but there is something specific to edit RNA alignments ?

e.g. Stockholm Alignment of Pistol (just few entries).

```
FP929053.1/1669026-1668956   AGUGGUCACAGCCACUAUAAACA-GGGCUU-UAAGCUGUG-AGCGUUGACCGUC----------ACAA-----CGGCGGUCAGGUAGUC
AFOX01000025.1/1981-1912     ACUCGUCUGAGCGAGUAUAAACA-GGUCAU-UAAGCUCAG-AGCGUUCACCGGG----------AUCA------UUCGGUGAGGUUGGC
HE577054.1/3246821-3246752   ACUCGUCUGAGCGAGUAUAAACA-GGUCAU-UAAGCUCAG-AGCGUUCACCGGG----------AUCA------UGCGGUGAGGUUGGC
CP000154.1/3364237-3364168   GUUCGUCUGAGCGAACGCAAACA-GGCCAU-UAAGCUCAG-AGCGUUCACUGGA----------UUCG------UCCAGUGAGAUUGGC`
`#=GC SS_cons                 <<<<__AAAAA_>>>>-------..<<<<-.----aaaaa.----<<<<<<<<<..........____....._>>>>>>>>>-->>>>`
`#=GC RF                      acUCGUCuggGCGAguAUAAAuA..cgCaU.UAgGCccaG.AGCGUcccggcgg..........uUau.....uccgccgggGGUuGcg
//

```

# Answer
I would suggest use RALEE—RNALignment Editor in Emacs. It can get for you the consensus secondary structure, you can move left/right sequences and their secondary structures (you can't do it in JalView!), and more.

It's an Emacs mode, so could be a bit hard to start off, but just try, you don't have to use all Emacs features to edit your alignments!

> The RALEE (RNA ALignment Editor in Emacs) tool provides a simple environment for RNA multiple sequence alignment editing, including structure-specific colour schemes, utilizing helper applications for structure prediction and many more conventional editing functions.

Sam Griffiths-Jones Bioinformatics (2005) 21 (2): 257-259.

Fig. You can move left/right sequences and their secondary structures (you can't do it in JalView!)
> 11 votes

# Answer
There's nothing really special about RNA alignments, you can use any alignment editor, including whichever one you use for protein. That said, a classic and very useful tool for this sort of thing is JalView. It can be installed locally or run as a Java webapp from your browser.

> Jalview has built in DNA, RNA and protein sequence and structure visualisation and analysis capabilities. It uses Jmol to view 3D structures, and VARNA to display RNA secondary structure. The Jalview Desktop can also connect with databases and analysis services, and provides a graphical interface to the alignment and analysis services provided by the JavA Bioinformatics Analysis Web Services framework.
> 7 votes

---
Tags: rna, rna-alignment
---


# Question
Title: What are the ways to process a list of differentially expressed genes?
We are studying six different human macrophage/dendritic cell types isolated from healthy skin. They all differ from each other in a few cell surface markers.

We are interested in the characteristics of each cell type ("marker" genes or biological processes), the differences between them (especially cell surface proteins), and their relations (e.g. "which are the closest relatives?") as can be inferred from the transcriptome, from an immunological point of view. The wider context is HIV-infection, thus HIV infection-related differences (and similarities) are of particular interest.

One naive approach is to contrast two of the most similar cell types (as inferred from previous knowledge, the sorting strategy and the PCA plot), and produce a list of differentially expressed (DE) genes. I have now a list of ~500 DE genes between two cell types, together with fold changes, expression means etc., produced with DESeq2 on bulk RNA-seq data.

What are the best practices to process this list to answer some of the above?

I've found PPI analysis on NetworkAnalyst, and GO and pathway analysis on InnateDB useful.

What other options are available?

I'm particularly interested in approaches using `R`.

# Answer
You originally had asked a very broad question, so I'll try to demonstrate why that is such a hard question to answer. I've done two fairly large differential analysis studies (and a few smaller ones) covering very different areas of research, and the approaches that other researchers used subsequent to my differential expression calculations were unsurprisingly also very different.

In the first study, I did a *de-novo* transcriptome assembly of *Schmidtea mediterranea (Smed)* from Illumina RNASeq reads, then carried out a DE analysis (DESeq) on a *β-catenin* knockout mutant population compared to wildtype. The other researchers had some very impressive tools available at their disposal, allowing them to subsequently visualise the positional expression of differentially-expressed genes in the knockout and wildtype populations after amputation. They found a few candidate genes that had graded expression along the anterior-posterior axis in the wildtype population, and this expression disappeared in the knockout mutants. One of these genes was a *teashirt (tsh)* family member, and a good candidate for further discovery based on the researchers' existing knowledge of the gene pathways associated with cellular regeneration. To cut a very long story short, they experimentally validated the *teashirt* link through further RNAi studies, generating a double-headed *Smed* as a result. This was followed up in Zebrafish *(Danio rerio)*, because there happened to be another building about 100m away that specialised in Zebrafish development.

The second study was/is a bit more similar to your current situation. We looked at sorted dendritic cell expression in mice after the application of two Th2-inducing treatment conditions (injected *Nippostrongylus brasiliensis*; applied DBP+FITC) and one control condition (injected PBS). I did a differential expression analysis with DESeq2, and developed a Shiny App for result exploration after I got tired of getting asked to do a whole bunch of mechanical follow-up tasks (mostly generating graphs, lists, and heatmaps for different gene subsets). We ended up with thousands of differentially expressed genes, so it wasn't feasible to look individually at every gene for research evidence and pathways -- one scientist spent a couple of weeks looking at one interesting-looking candidate gene. The differentially-expressed genes were fed into the oh-so-expensive IPA (which gives results, but not necessarily good ones) and DAVID, among a few others. These helped a little bit, and led to what was the central idea of our paper (i.e. that there were at least two different Th2 responses), but we still had the problem of too many "significant" pathways, and lots of differentially-expressed genes that were members of a whole bunch of other pathways. We tried WGCNA; we did a bit of GSEA; and we compared Limma and DESeq2, but still couldn't find anything concrete or obvious to demonstrate what was triggering the Th2 response. There were experimental studies as well: blocking and mutant studies to investigate the two types of Th2 response (and experimentally confirm the differential expression). Our experiments and exploration are still ongoing, but we felt it was necessary to get one paper out at the intermediate stage, so just ended up taking the top 25 differentially-expressed genes in each experiment (by fold change), and finding some way of grouping some of them together (via expert knowledge and literature searches). It kind-of-sort-of worked, but we're left with a bit of frustration in that we (and the tools we used) were basically stumped by having too many biologically-relevant results.
> 7 votes

---
Tags: rna-seq, r, deseq2, networks
---


# Question
Title: How to select a power for a scale-free topology network
In a weighted gene co-expression network analysis (using WGCNA), the soft-threshold power is recommended as a noise filtering. It consists on raising the correlation to a certain number. To decide this power the scale-free topology is estimated for some powers. The function to estimate this prints:

```
   Power SFT.R.sq  slope truncated.R.sq mean.k. median.k. max.k.
1      1   0.9300  3.110          0.996  2960.0    3060.0   3970
2      2   0.7510  1.010          0.964  1750.0    1780.0   2900
3      3   0.1730  0.258          0.806  1170.0    1150.0   2280
4      4   0.0942 -0.183          0.713   833.0     782.0   1870
5      5   0.3800 -0.463          0.777   623.0     559.0   1580
6      6   0.5350 -0.656          0.834   481.0     412.0   1360
7      7   0.6270 -0.797          0.872   381.0     312.0   1190
8      8   0.6870 -0.910          0.900   307.0     241.0   1050
9      9   0.7270 -1.000          0.918   252.0     189.0    936
10    10   0.7490 -1.080          0.928   210.0     150.0    841
11    12   0.7850 -1.190          0.948   150.0      98.0    693
12    14   0.8090 -1.280          0.958   111.0      65.9    582
13    16   0.8290 -1.360          0.968    84.0      45.6    497
14    18   0.8410 -1.410          0.973    65.2      32.1    429
15    20   0.8490 -1.450          0.977    51.6      23.0    375

```

The recommendations of the FAQ indicate that a SFT.R.sq value should be above

> 0.8 for reasonable powers (less than 15 for unsigned or signed hybrid networks, and less than 30 for signed networks)

and the mean connectivity below the hundreds.

Others have used the power just as noise filtering without caring much about the scale-free topology fit. I would pick the first power even if the mean connectivity is in the order of thousands because the scale-free topology fit is pretty high, however the slope is puzzling me.

How should the soft-threshold power be selected?

Based on an example question

# Answer
I would be concerned about the input data that produced these numbers. This dataset has a really odd shape with a dip at the start:

```
data.df <- read.table("soft_threshold.txt", header=TRUE);
png("scale-free-194.png");
plot(data.df$SFT.R.sq);
dummy <- dev.off();

```

Here's what I expect this to look like (from the supplementary information in this paper):

It's easier to select the power/threshold when the plot looks like the second situation. With the data as presented, it doesn't look like any value would be suitable. but if I were forced to pick, then I'd choose 12 because it is the first number above the magic threshold (ignoring the weirdness at the start). However, I'll admit that I don't understand the concepts behind the power, and am just following the recommended instructions, so can't explain why it is what it is.
> 2 votes

---
Tags: statistics, networks
---


# Question
Title: Protein folding problem model for quantum computers?
It's supposed that quantum computers could help to find the native state of all human DNA proteins.

Does anyone know how to formulate the protein folding problem into a decision problem solvable either by quantum annealers or universal quantum computers?

It would also be interesting to know, how to model a decision problem for classical computers, like Folding@Home.

# Answer
> 5 votes
Protein folding problem can be viewed as a minimization problem. One can use quantum annealing to perform the minimization. Running this on quantum computers would improve the performance since they can perform the tunneling directly.

In fact, quantum annealing was used (blog post) for lattice protein folding on the the D-Wave quantum computer (128 qubits). But this was an extremely simplified model, with only 40 discrete states.

# Answer
> 1 votes
I would expect that without any further insights into protein folding, a quantum method would look fairly similar to existing methods, but might have a bit more capability of exploring multiple states at the same time (without multitasking).

Here's the approach used by Rosetta, which is used in both the Rosetta@home distributed computing project, and the foldit distributed gaming project:

> Rosetta's strategy for finding low energy shapes looks like this:
>
> 1. Start with a fully unfolded chain (like a metal chain with its ends pulled).
> 2. Move a part of the chain to create a new shape.
> 3. Calculate the energy of the new shape.
> 4. Accept or reject the move depending on the change in energy.
> 5. Repeat 2 through 4 until every part of the chain has been moved a lot of times.
>
> We call this a trajectory. The end result of a trajectory is a predicted structure.

---
Tags: proteins, computation
---


title: How to use Python to count k-mers?

text:
# Question
Title: How to use Python to count k-mers?
I have some FASTQ sequence files and a FASTA file for some regions I'm interested in.

I would like:

1. Build an index for the FASTA file
2. Use the index to count number of k-mers occurred in my sequence files

I know how to do this in many k-mer counting tools, such as Kallisto, Jellyfish, DSK, Salmon. However, I'm more interested in Python (khmer is a Python package, all the other k-mer counting tools are C++).

Unfortunately, I can't figure out how to do that in khmer from the documentation:

> http://khmer.readthedocs.io/en/v2.1.1/introduction.html

It talks about RNA-Seq, metagenomics, de Bruijn graph etc... All I need are just counts. Khmer is a popular k-mer tool, so it should be able to do it, but how?

**Q:** How to build an index and count k-mers using khmer?

**EDIT:**

I don't need number of unique k-mers, I'm looking for a solution like Jellyfish where we build an index and count k-mers to the index.

Something like:

```
MY_SEQ_1  100 counts
MY_SEQ_2  134 counts

```

# Answer
If you want to count the number of unique k-mers that occur in your data set, you should use the unique-kmers.py script, which implements a HyperLogLog-based cardinality estimator.

If the number of unique k-mers is not what you're after, please clarify.

**EDIT**:

khmer uses probabilistic data structures internally, which store k-mers as hashed values that cannot be un-hashed back into k-mer sequences uniquely. Querying a k-mer's abundance requires you to know which k-mer(s) you're looking for, which in this case would require a second pass over the reads.

The closest thing to what you're asking provided in khmer's command-line scripts is abundance-dist.py (or the alternative abundance-dist-single.py), which will produce a k-mer abundance histogram but not per-kmer abundance.

The load-into-counting.py scripts computes the k-mer abundances and stores them in a probabilistic data structure, which is written to disk and can subsequently be re-loaded into memory quickly. Many of the scripts in khmer expect the reads to have been pre-processed by `load-into-counting.py`, while others will invoke the counting routines directly themselve. It's also fairly easy to do this with the Python API. In every case, the accuracy of the k-mer counts is dependant on memory considerations.

```
>>> import khmer
>>> counts = khmer.Counttable(31, 1e7, 4)
>>> counts.consume_seqfile('reads.fq.gz')
(25000, 1150000)
>>> counts.get('TGACTTTCTTCGCTTCCTGACGGCTTATGCC')
117

```

To check the false positive rate:

```
>>> fpr = khmer.calc_expected_collisions(counts)
>>> fpr

```

If you *must* have the counts of each k-mer, it's not much more work...although reporting each k-mer's abundance only once would consume a lot of memory with a naive approach.

```
>>> outfile = open('outfile.txt', 'w')
>>> seenkmers = set()  # Consumes a lot of memory for large input!!!
>>> for read in khmer.ReadParser('reads.fq.gz):
...     for kmer in counts.get_kmers(read.sequence):
...         if kmer not in seenkmers:
...             print(kmer, counttable.get(kmer), file=outfile)
...             seenkmers.add(kmer)

```
> 5 votes

# Answer
I wrote a command-line k-mer counter called `kmer-counter` that will output results in a form that your Python script can consume: https://github.com/alexpreynolds/kmer-counter

You can grab, build and install it like so:

```
$ git clone https://github.com/alexpreynolds/kmer-counter.git
$ cd kmer-counter
$ make
$ cp kmer-counter /usr/local/bin

```

Once the binary is in your path, you might use it in Python like so:

```
k = 6
fastaFile = '/path/to/some/seqs.fa'
kmerCmd = 'kmer-counter --fasta --k=%d %s' % (k, fastaFile)
try:
    output = subprocess.check_output(kmerCmd, shell=True)
    result = {}
    for line in output.splitlines():
        (header, counts) = line.strip().split('\t')
        header = header[1:]
        kmers = dict((k,int(v)) for (k,v) in [d.split(':') for d in counts.split(' ')])
        result[header] = kmers
    sys.stdout.write("%s" % (str(result)))
except subprocess.CalledProcessError as error:
    sys.stderr.write("%s" % (str(error)))

```

Given example FASTA like this:

```
>foo
TTAACG
>bar
GTGGAAGTTCTTAGGGCATGGCAAAGAGTCAGAATTTGAC

```

For k=6, you would get an iterable Python dictionary like this:

```
{'foo': {'TTAACG': 1, 'CGTTAA': 1}, 'bar': {'GTTCTT': 1, 'AGAACT': 1, 'GAGTCA': 1, 'ATGGCA': 1, 'GAACTT': 1, 'ATTCTG': 1, 'CTAAGA': 1, 'CTTCCA': 1, 'ATTTGA': 1, 'GGAAGT': 1, 'AGGGCA': 1, 'CCTAAG': 1, 'CTCTTT': 1, 'AATTTG': 1, 'TCTGAC': 1, 'TTTGCC': 1, 'CTTAGG': 1, 'TTTGAC': 1, 'GAAGTT': 1, 'CCCTAA': 1, 'AGAATT': 1, 'AGTCAG': 1, 'CTGACT': 1, 'TCTTAG': 1, 'CGTTAA': 1, 'GTGGAA': 1, 'TGCCAT': 1, 'ACTCTT': 1, 'GGGCAT': 1, 'TTAGGG': 1, 'CTTTGC': 1, 'TGGAAG': 1, 'GACTCT': 1, 'CATGCC': 1, 'GCAAAG': 1, 'AAATTC': 1, 'GTCAAA': 1, 'TGACTC': 1, 'TAGGGC': 1, 'AAGTTC': 1, 'ATGCCC': 1, 'TCAAAT': 1, 'CAAAGA': 1, 'AACTTC': 1, 'GTCAGA': 1, 'CAAATT': 1, 'TAAGAA': 1, 'CATGGC': 1, 'AAGAAC': 1, 'AAGAGT': 1, 'TCTTTG': 1, 'TTCCAC': 1, 'TGGCAA': 1, 'GGCAAA': 1, 'AGTTCT': 1, 'AGAGTC': 1, 'TCAGAA': 1, 'GAATTT': 1, 'AAAGAG': 1, 'TGCCCT': 1, 'CCATGC': 1, 'GGCATG': 1, 'TTGCCA': 1, 'CAGAAT': 1, 'AATTCT': 1, 'GCATGG': 1, 'ACTTCC': 1, 'TTCTTA': 1, 'GCCATG': 1, 'GCCCTA': 1, 'TTCTGA': 1}}

```

You can use standard Python calls to manipulate this dictionary object and get sums of counts per record, for sequence, etc. which seems to answer your question. Please feel free to clarify what you're looking for if this object representation is not clear.

For a fully-fleshed out demonstration, see: https://github.com/alexpreynolds/kmer-counter/blob/master/test/kmer-test.py
> 7 votes

---
Tags: genome, k-mer, metagenome, khmer
---


title: Annotation format design

text:
# Question
Title: Annotation format design
Bashing file formats is a favorite pastime in bioinformatics, and annotation file formats such as GFF and BED seem to get special attention. A lot of this frustration stems from community's shockingly inconsistent adherence to specifications and conventions, but there are also some (dare I say objectively) problematic design choices in each of these formats.

* GFF (and its more common derivatives GTF and GFF3) use 1-based closed interval notation, which optimizes for human comprehension but is far inferior to 0-based half-open interval notation (such as used by BED) for computations involving interval arithmetic.
* Although BED and GTF were designed for very specific use cases (visualization and gene prediction, respectively), they have been used and abused in a much wider set of contexts. For example, the BED fields related to the **thick part** are irrelevant if you're not plotting them in a genome browser.
* BED supports a single level of feature decomposition (a feature can be broken up into blocks). GTF supports two levels (exons grouped by transcript\_id, transcripts grouped by gene\_id). In contrast, GFF3 supports an arbitrary number of levels, and uses parent/child relationships defined by `ID` and `Parent` attributes to declare a directed acyclic graph of features.
* Data that does not fit into mandatory pre-defined fields must be relegated to optional fields or free-form attribute key/value pairs. While this flexibility is powerful, a common complaint is that "all the action" happens in these optional/free-form fields.
* There is a dearth of validation tools, and those that do exist focus primarily on validating syntax and not semantics. To use an aging analogy, it's one thing to say an XML file is valid, but it's completely different to validate it against a schema. There are essentially no widely used tools that do the latter for annotation files.

If we were tasked with creating a new annotation format, and if we were guaranteed the resources needed to develop it, and interest and wide adoption from the wider community (one can dream!), what design criteria should be considered in the development of this new format? **What, if anything, makes an objectively good annotation data format?**

# Answer
Presuming we consider "human readable", "easily parsable", and "quickly queryable" to be objectively good qualities (and if not, I worry for the future):

1. Text-based: It's absurdly common to want to use `grep` or `awk` on annotations. Sure, one could make variants of these that are binary-format aware, but why reinvent the wheel. Of course text files don't innately allow region-based queries of their contents, so on to point 2...
    * This should further explicitly be line-based. Fields would need to be tab-separated (ideally it'd use the ascii record separator character instead, but I fear that ship has long ago sailed).
2. A strictly defined line order: This is one of my personal pet-peeves about BED/GFF/GTF files. If you're going to make a text-based annotation format, everyone would end up ahead if said format made explicit that it should be sorted. This then allows things like tabix to be used so the "query a region" problem is then solved. But I would go further than that. The issue with things like GFF is that there are multiple inter-dependent lines and there's no strict rule about whether a parent line absolutely must come before a child line. This just makes implementing things a nightmare and more often than not a randomly ordered file will just break tools. Since GTF or GFF-style interdependent lines likely be how any annotation format works, this ordering between lines with the same start position should be made explicit in the format.
> 4 votes

# Answer
The problem is that GFF, fundamentally, is a relational format: it provides tags that relate individual rows via one-to-many relationships (e.g. gene–exon). This indirectly highlights the second complication: individual rows have different types, and therefore store different attributes in the 9th column.

Over the last few decades (!), we have accumulated a wealth of theory and tools to work with this kind of data. And the usual solution is to create a database schema for a relational database, and to use database drivers and database query languages (e.g. SQL, but increasingly also data-relational mappers such as LINQ and dplyr).

Using a text-based format is attractive for many reasons that Devon mentions, but it is fundamentally at odds with a lot of the theory and tools for relational data. This creates an impedance mismatch.

I’m convinced that, in the long run, the solution is going to be to revert to using relational databases for complex annotations. I say “revert” even though these databases already exist, they are simply often ignored in bioinformaticians’ day-to-day work (*I* never use them). Because this is objectively the best technical solution, and we have the research to back this up.
> 4 votes

# Answer
I quite like BED and GFF3 (I don't like GTF/GFF2, though). As text-based formats, I don't think they leave us much room for improvement. Anyway, if you want a new format, here is one. The following is a hybrid between GFF3 and BED. It is a TAB-delimited text-based format with the following fields:

1. chr, required
2. start (0-based), required
3. end, required
4. strand
5. ID
6. type (see below)

Like BED, only the first 3 columns are required, the rest are optional. What is different starts with the "type" field as in GFF. This "type", when present, defines the columns following it. For example, if type==coding, we could have cdsStart, cdsEnd, blockCount, blockSizes and blockStarts like BED; if type==exon, column 7 could be a "phase", indicating the phase of the first base. This way, "type" makes this format highly extensible while still relatively easy to parse in comparison to using optional tags all the way. In addition, we may have semi-colon-separated "key"="value" pairs at the end of each line as in GFF3. Example:

```
chr1   10000  50000  -  x1 coding   10100   40800  2  1000,2000   10000,48000
chr1   10000  11000  -  *  exon     *       foo1=bar1;foo2=bar2
chr1   48000  50000  -  *  exon     2
chr2   10000  50000

```

The above only gives a bare bone of the format. There are subtle questions such as: 1) the uniqueness/scope of ID; 2) whether to have parent ID as a fixed field or a type-specific field; 3) where to put display name; 4) sorting order. These are also important in practice.

PS: I like SQL a lot and think it should get used more often in bioinformatics, but I don't think it replaces text formats completely. Formats are useful for serialization and data transfer. They require less skill to work with and less software/hardware resource to deploy. Carefully engineered binary representations of formats can be much more efficient in a lot of use cases.
> 4 votes

# Answer
Just throwing this out there, as all the best stuff has already been said, but why does one need to specify a strict specification for a bioinformatic data format? As you say, what typically happens is all the action ends up in the optional fields.

There's a lot to be said for the libertarian values of just letting people figure it out on their own. Take optional fields in the BAM tag spec. for example. Here you can have your own tags, but they must start with an "X", a "Y" or a "Z" and be followed by "\[A-Za-z0-9\]". Why? Why can't people use their own tag names like "read is unique" or "edit distance to genome" or whatever they want. Are we not to be trusted with the power of naming things and recalling things arbitrarily? The result is one must look up the tag name's actual meaning in the publication of the tool that produced it, or ask on certain forums, etc. And this is assuming one knows the tool that produced the reads - if not then who knows wtf 'XT' stands for.

Essentially, downstream readers of optional fields already trust that 'XT' is what it thinks it is. So if we're happy to use optional fields in any capacity, why stop at fields?

Extend this logic to the columns of the data format. Let users determine column names. You may once in a blue moon come across a file with the chromosome column as "chr" rather than "chromosome", but most of the time you'll be looking for "read is quantum-radioactive", and detecting and fixing naming-errors is not going to be as much of a problem as using a dataformat that does not store what it is you want to store. Or worse, it can store it, but in a really illogical way that results in everyone who uses the dataformat to have to ask at least 3 questions here or on other forums before they understand whats really going on.

So you end up with a general solution to the problem of storing tabular information in a compatible way across many different types of software, also known as SQL (or any other kind of general-purpose database such as REDIS, Neo4J, Numpy, etc). In fact, does it even matter what the datastore is, so long as it's tabular and has the "chromosome" and "position" values for each item?

TL;DR - We will not think of tomorrow's best data format today, because the nature of tomorrow's data is not yet known. Less policing in this area would most likely result in more robust software, where nothing is taken for granted and no assumptions about the data can be made until the schema has been parsed.

In truth, the only reason we think we need to specify our dataformats is because we know if we didn't, other people will do something we didn't think of - and we, i think falsely, assume that will be bad.
> 0 votes

---
Tags: file-formats, bed, gff3, sequence-annotation
---


title: How to calculate Fst from AMOVA

text:
# Question
Title: How to calculate Fst from AMOVA
I calculated an AMOVA from a genind object, with one hierarchical factor. In the table I obtain there are SSD values (for my grouping factor,"Error" and total) and sigma2 values (for my grouping factor and "Error")

I have two questions:

* What does "error" stand for?
* How do I calculate Fst? Which values do I have to use?

A genind object in R is an object which contains allelic information about a set of individuals.
This is the table that I get:

```
Analysis of Molecular Variance

Call: amova(formula = plantula.dist ~ g, nperm = 100)

           SSD         MSD  df
g     354992.9 354992.9297   1
Error 310814.9    693.7834 448
Total 665807.9   1482.8683 449

Variance components:
       sigma2 P.value
g     1574.69       0
Error  693.78

Variance coefficients:
       a
224.9956

```

# Answer
> 1 votes
Error is the estimation that can't be fitted by your parameters (in this case g). The amova function in pegas package does the following:

> This function performs a hierarchical analysis of molecular variance as described in Excoffier et al. (1992). This implementation accepts any number of hierarchical levels

And as per the definition of Fixation index, you need the variance in the frequency of the allele between different subpopulations, weighted by the sizes of the subpopulations, the variance of the allelic state in the total population. Which is not provided by amova. Thus you cannot compute the fixation index from amova results.

---
Tags: r
---


# Question
Title: Is it possible to calculate p-values for fold changes of single replicate RNA-seq samples?
I'm thinking that this isn't actually possible, but I'd like to check before I write it off. I have RNA-seq data from cells under three different conditions, there are no replicates for any of the conditions. The data has been processed with RSEM, and log2 fold changes have been calculated for each control-test pairing using the normalized expected read counts using EBseq.

If possible, I'd like to also calculate the p-value for each of these fold-changes, however, because there are no replicates I don't think that this is possible. Is there any way to do this?

# Answer
> 3 votes
DESeq2 is able to adjust p-values for single-replicate samples by estimating shared dispersion across all samples. DESeq2 will give a warning, but try its best to carry out the analysis. More information on doing that can be found in the DESeq2 vignette and the DESeq2 manual page.

The results will not be anywhere near the sensitivity that can be achieved when replicate information is available, but it might be good enough if only very obvious large-scale changes are of importance.

---
Tags: ebseq, rsem
---


# Question
Title: What are the de facto required fields in a SAM/BAM read group?
The SAM specification indicates that each read group must have a unique ID field, but does not mark any other field as required.

I have also discovered that htsjdk throws exceptions if the sample (SM) field is empty, though there is no indication in the specification that this is required.

Are there other read group fields that I should expect to be required by common tools?

# Answer
> 9 votes
The sample tag (i.e. SM) was a mandatory tag in the initial SAM spec (see the `.pages` file; you need a mac to open it). When transitioned to Latex, this requirement was mysteriously dropped. Picard is conforming to the initial spec. Anyway, the sample tag is important to quite a few tools. I would encourage you to add it.

# Answer
> 3 votes
The read group identifier needs to be specified in both the header lines of the BAM/SAM file and the alignment line. No other fields are required, but note that because the additional information is only stored once (i.e. in the header), it won't add much to file sizes or processing time if additional fields are included. If any of the other optional fields cause problems (either by inclusion or exclusion), it would be helpful to report an issue about it.

Older tools required you to specify both a header read group ID and an alignment read group ID, but most tools of the recent tools I've used seem to be intelligent enough that they will add headers for any read groups without that information.

---
Tags: file-formats, sam
---


# Question
Title: Do variant calls change when you call from CRAM?
We're considering switching our storage format from BAM to CRAM. We work with human cancer samples, which may have very low prevalence variants (i.e. not diploid frequency).

If we use lossy CRAM to save more space, how much will variants called from those CRAM files change? Which compression strategy has the lowest impact?

Are there any other impacts on downstream tools that we're not considering?

# Answer
> 6 votes
By default, a CRAM you create with samtools is lossless. It typically halves the input BAM in terms of file size. If you want to compress more, you can let samtools convert most read names to integers. You won't be able to tell optical duplicates from read names, but this is a minor concern. You can also drop useless tags depending on your mapper and the downstream caller in use. For cancer data, I wouldn't reduce the resolution of base quality without comprehensive benchmarks. Unfortunately, base quality takes most of space in CRAM. Discarding the original read names and some tags probably won't save you much space.

# Answer
> 3 votes
The main concern has always been the "binning" of quality scores that occurs via CRAM compression (and is also standard on the HiSeqX, HiSeq4000, and NovaSeq platforms). Anecdotally, I can report very little difference between 4-bin quality scores and full quality scores on cancer samples, though I don't know if I've seen a direct head-to-head comparison.

---
Tags: sam, file-formats, variant-calling, storage
---


# Question
Title: What unit do I get on the y-axis of a metagene profile plot?
I start with a sorted and indexed bam file ("mapped.bam") representing the mapping of small reads on a reference genome, and a bed file ("genes.bed") containing the coordinates of a set of features of interest (let's say they are genes), for which I want to compute an average profile using programs from deeptools. **I would like to understand the steps involved to be sure of what the vertical axis of the final profile represents.**

### First step: making a bigwig file

I create a bigwig file ("mapped.bw") from the bam file using `bamCoverage` as follows:

```
bamCoverage -b mapped.bam -bs 10 -of=bigwig -o mapped.bw

```

The help of `bamCoverage` says:

> The coverage is calculated as the number of reads per bin, where bins are short consecutive counting windows of a defined size.

In my case, the bins are 10 bp long. My reads are longer than that.

For a given bin, a given read can:

1. completely overlap the bin
2. overlap the bin on n bp, n \< 10
3. not overlap the bin at all

Please correct me if I'm wrong: My guess it that the read is counted as 1 in cases 1. and 2., and 0 otherwise, and I also suppose that a read can be counted for several successive bins if it is long enough.

### Second step: averaging over genes and plotting

I compute a "meta profile matrix" ("mapped\_on\_genes.gz") using `computeMatrix scale-regions` as follows:

```
computeMatrix scale-regions \
    -S mapped.bw \
    -R genes.bed \
    --upstream 300 \
    --unscaled5prime 500 \
    --regionBodyLength 2000 \
    --unscaled3prime 500 \
    --downstream 300 \
    -out mapped_on_genes.gz

```

(There is a `-bs` parameter which default value is 10 according to the help of the command.)

I use this to plot a profile using `plotProfile`:

```
plotProfile -m mapped_on_genes.gz \
    -out mapped_on_genes_meta_profile.pdf

```

I obtain a profile in with values on the y axis. In what units are these values?

My guess is the following:

For the upstream (300 bp) and internal 5-prime (500 bp), since the bin size was the same in `bamCoverage` and `computeMatrix`, each point on the x axis probably represents a 10 bp window, and its y coordinate is the average over the regions present in the bed file of the corresponding bins in the bigwig file, so it is an average number of reads overlapping a 10 bp bin.

Same thing at the 3-prime and downstream side.

For the central 100 bp portion, before averaging over regions some shrinking or spreading of the bins must have been performed, I guess by averaging between neighbouring bins. So the final unit is still a **number of reads overlapping a 10 bp bin**.

And if I use larger bins, I should end up with proportionally higher values.

**Am I correct?**

# Answer
Feel free to @ me in deepTools questions, since I'm the primary developer.

For a given bin, the count assigned to it is the number of reads overlapping it, regardless of whether they overlap by 1 or 10 bases. So a read overlapping only partially and one overlapping completely are treated the same.

Since your bigWig file is in units of "alignments" (i.e., it's not 1x normalized), the resulting profile will also be in units of "alignments" (i.e., profiles and heatmaps are in whatever units the input files are in).

Upstream/downstream regions and unscaled regions are also 10 base bins. Note that these are then the average of the per-base value, since the bins here may not perfectly correspond to the bins in the bigWig files. The line in the profile plot is indeed the average (by default, you can choose median, max, min, etc.) of the underlying regions for each bin.

Regarding the scaled section in the middle, the the number of genomic bases per bin is changed such that the region will have `"length"/(regionBodyLength/binSize)` bases each. As above, the per-base value is then averaged (or whatever you specify) to derive the per-bin value. The `length` here is decreased if you have unscaled regions, since otherwise bases would get counted twice.
> 4 votes

---
Tags: hts, deeptools, bigwig
---


# Question
Title: How do I efficiently subset a very large line-based file?
This has come up repeatedly recently: I have a very large text file (in the order of several GiB) and I need to perform line-based subsetting for around 10,000 lines. There exist solutions for specific scenarios (e.g. `samtools view -s` for randomly sampling BAM files) but sometimes my use-case doesn’t fit into these categories.

Unfortunately a naïve `sed`-based solution is *extremely* slow:

```
time sed -n -f <(awk -vOFS='' '{print $0, "p"}' line_numbers.txt) input_file > selected_lines.txt

```

Where `line_numbers.txt` is a file containing one line number per line.

Forget running this for 10,000 lines; it’s already grinding to a halt for a mere 1000.

How can I speed this up, ideally so that it scales only with the size of the input file, and has more or less constant runtime n the number of lines that I subset?

# Answer
> 7 votes
Turns out, simply keeping track of the next candidate line (after sorting the sample line numbers) fixes the performance issue, and most of the remaining slowness seems to be due to the overhead of actually reading the file so there’s not very much to improve.

Since I don’t know how how to do this in `sed`, and it’s not trivial in `awk` either, here’s a Perl script:

```
#!/usr/bin/env perl

use strict;
use warnings;

my $file = $ARGV[0];
my $lines_file = $ARGV[1];

open my $lines_fh, '<', $lines_file or die "Cannot read file $lines_file";
chomp (my @lines = <$lines_fh>);
close $lines_fh;

@lines = sort {$a <=> $b} @lines;

open my $fh, '<', $file or die "Cannot read file $file";
my $line = 1;
my $next_line = 0;
while (<$fh>) {
    last if $next_line == scalar @lines;
    if ($line++ == $lines[$next_line]) {
        $next_line++;
        print;
    }
}
close $fh;

```

I’ve implemented a similar function in C++ for an R package, that's only slightly longer than the Perl script. It is ~3 times faster than the Perl script on my test file.

# Answer
> 5 votes
Perl should be fairly fast with this when using a hash set to store the list of lines. A structure like this also works for subsetting based on a field value, where the comparison would be with the field rather than "$.":

```
#!/usr/bin/perl

use strict;
use warnings;

my $lines_file = $ARGV[0];
my %include_lines = ();

open my $lines_fh, '<', $lines_file or die "Cannot read file $lines_file";
while(<$lines_fh>){
  chomp;
  $include_lines{$_} = 1;
}
close $lines_fh;

while(<>){
  if($include_lines{$.}){ # "$." -- line number of current file
    print;
  }
}

```

Note that according to this SO answer, the "$." operator is not strictly the current line number, and can be influenced by different file operations or other settings.

Edit: just saw your comment about speed in your answer, comparing hash sets to a sorted list. The `$lines[$next_line]` bit feels a bit odd to me. Have you tried out using `shift` or `pop` on a sorted list to fetch the next line:

```
#!/usr/bin/perl

use strict;
use warnings;

my $lines_file = $ARGV[0];

open my $lines_fh, '<', $lines_file or die "Cannot read file $lines_file";
chomp (my @lines = <$lines_fh>);
close $lines_fh;

@lines = sort {$a <=> $b} @lines;
my $next_line = shift(@lines);

while (<>) {
    if ($. == $next_line) {
        $next_line = shift(@lines);
        print;
        last if (!@lines);
    }
}

```

I changed the `shift` to a `pop` (reversing the sort order to match), and got times of 4.8s, 5.8s and 4.1s for Konrad's original code, my hash code, and my pop code respectively, fetching 10,000 lines from `/usr/share/dict/british-english-insane` 25 times over (after copying input files to `/tmp`). They're all the same from my perspective: quick enough that it would take me longer to type out the command than it would to run it. Using `shift` instead of `pop` doesn't seem to change the time noticeably.

# Answer
> 4 votes
Some related questions appear in other sites, with potentially interesting solutions, which I report here:

To sample approximately 1% of the non-empty lines:

```
awk 'BEGIN {srand()} !/^$/ { if (rand() <= .01) print $0}' input_file

```

(from https://stackoverflow.com/a/692321/1878788)

To select 1000 random lines:

```
shuf -n 1000 input_file

```

(from https://stackoverflow.com/a/15065490/1878788, and https://unix.stackexchange.com/a/108604/55127)

### Edit: Python solutions using a list of lines

Using a set of line indices and selecting lines by testing set membership:

```
#!/usr/bin/env python3

import sys

with open(sys.argv[2], "r") as line_numbers_file:
    line_indices = set(int(line) - 1 for line in line_numbers_file)

with open(sys.argv[1], "r") as input_file:
    print(*(line.strip() for (idx, line) in enumerate(input_file)
            if idx in line_indices), sep="\n")

```

Using a numpy boolean array together with `itertools.compress`:

```
#!/usr/bin/env python3

import sys
from itertools import compress
from numpy import zeros

with open(sys.argv[2], "r") as line_numbers_file:
    line_indices = [int(line) - 1 for line in line_numbers_file]

selector = zeros(max(line_indices) + 1, dtype=bool)
selector[line_indices] = 1

with open(sys.argv[1], "r") as input_file:
    print(*(line.strip() for line in compress(input_file, selector)), sep="\n")

```

I did some tests on a file containing 15774756 sam records and a list of 10000 pre-generated random line numbers.

The perl script proposed by Konrad Rudolph (https://bioinformatics.stackexchange.com/a/454/292) runs in about 5.3 seconds.

The set membership testing python solution runs in about 4.45 seconds.

The compress based solution runs in about 3.4 seconds. I suspect this may vary a lot depending on the highest line number we want, since the number of iterations will depend on the length of the boolean array. Here the highest line number was 15773768, so pretty high compared with the total number of lines.

I tried with python 3.6. I suspect that python 2.7 could be slightly faster, but haven't tested.

# Answer
> 1 votes
I wrote a command-line (C++14) tool called `subset` which is up on Github: https://github.com/alexpreynolds/subset

This should be reasonably memory efficient and fast. The `subset` tool does not store input lines in a table, but instead streams through the file once, storing a 4 or 8k buffer chunk of the input file (depending on OS).

It stores line numbers in an array, but eight bytes per integer * 100k is 800kB, for that use case — not very much memory there.

There's an O(nlogn) sort penalty on the line number array, but again this list will be much smaller than the query file, and integer sorting is fairly optimized, so the hit should be small.

If your line number list is already sorted, I could add an option to skip sorting; let me know if that would be useful.

The filtering step walks through the line number array and input file linearly, printing lines where there are index matches, and skipping over the rest.

Indeed, `subset` will quit early in parsing the input file if there are no more line numbers to query. So this feature is especially useful for speeding up filtering of very large query files. (If your query file has 1M rows, say, and your last line number of interest is 12345, there's no reason to read through the rest of the file.)

You can grab, build and install it like so:

```
$ git clone https://github.com/alexpreynolds/subset.git
$ cd subset
$ make
$ cp subset /usr/local/bin

```

Once the binary is in your path, there are a couple ways to use it.

For example, you can specify a start index and length value. The following grabs seven lines starting with the 33rd line (32 as a 0-indexed value):

```
$ subset --prefix-with-indices -s 32 -n 7 -i query.txt > answer.txt

```

Or you can specify a text file containing line numbers, each on a separate line. The following reads in a file called `line-numbers.txt` and uses that to filter `query.txt`:

```
$ subset --prefix-with-indices -l line-numbers.txt -i query.txt > answer.txt

```

The indices in `line-numbers.txt` should be positive, 0-indexed integers. The list of numbers does not need to be sorted, as `subset` will sort the list of numbers for you. This is so that an efficient single pass through the input/query file can be done.

You can leave out `--prefix-with-indices` to leave out the debug prefix. This is there so that you can do a sanity check on the result.

The `test/makefile` tests demonstrate options and usages for the two types of filtering.

---
Tags: shell, benchmarking, subset, text
---


# Question
Title: What is the best method to estimate a phylogenetic tree from a large dataset of >1000 loci and >100 species
I have a large phylogenomic alignment of \>1000 loci (each locus is ~1000bp), and \>100 species. I have relatively little missing data (\<10%).

I want to estimate a maximum-likelihood phylogenetic tree from this data, with measures of statistical support on each node.

There are many phylogenetics programs that claim to be able to analyse datasets like this (e.g. RAxML, ExaML, IQtree, FastTree, PhyML etc). Given that I have access to a big server (512GB RAM, 56 cores), what are the pros and cons of each program.

Which programme likely to give me the most accurate estimate of the ML tree for a dataset of this size within a reasonable timeframe?

# Answer
This paper claims that FastTree is almost as accurate as RAxML, while being much faster. You just have to be careful, however, that the support values output by FastTree are not bootstrap values, they are based on the Shimodaira-Hasegawa test. (Also, see this comment for the case you have very short branch lengths). \[**update:** However, according to the recent comparison paper mentioned below FastTree performed quite poorly in comparison to RAxML or IQ-tree.\]

From what I understand, you should use ExaML only if your data is too large to be handled by RAxML in a single node. ExaML should perform like RAxML but with some parallelization overhead. For all effects I treat them as the same. I don't know of relevant advantages of phyML over RAxML (for me, it's easier to use but I am very used to phyML).

I am not familiar with IQ-tree, but its authors claim that even given the same time as RAxML or phyML, IQ-tree already finds better likelihoods more often than not (although by default it takes a bit longer to converge). A recent comparison between all these programs favoured IQ-TREE for both single-gene and concatenation analyses (with RAxML very close). It may also estimate branch support through a SH-like test only, but I'm not sure. \[**update:** IQ-tree offers 3 measures of support, standard bootstrap, aLRT, and ultrafast bootstrap. See OP's comment below for details.\]

However, since you have few missing data, you might also want to try a single-locus tree inference followed by gene tree clustering (using treescape or treeCL) to see how spread your data is, or to see the effect of removal of outliers, or to use ideas similar to statistical binning.
> 11 votes

---
Tags: phylogenetics, phylogeny
---


# Question
Title: How can I combine a kinship matrix with subset individuals when using rvtests?
I've been processing data for a consortium project that uses the rvtests toolkit. Our data analysis process uses a kinship matrix and the calculation of inverse-normalised transformed statistics after covariate adjustment (for 1kG-imputed data from a SNPchip). One of the operations we have been asked to do is to analyse males and females separately for a number of different phenotypes.

Here's an example of the command line that I'm running:

```
rvtest --inVcf chr${ch}.imputed.poly.vcf.gz \
 --pheno ../scripts/phenotypes.raw.2017-May-01.ped \
 --pheno-name ${traitName} \
 --out rv_results/${traitName}_FEMALE_chr${ch} \
 --kinship kinship_matrix.kinship \
 --meta score,cov[windowSize=500000] \
 --peopleIncludeFile ../scripts/female.iids.txt \
 --covar ../scripts/covar.2017-May-11.ped \
 --covar-name ${covs} --useResidualAsPhenotype --dosage DS

```

Unfortunately, the kinship matrix doesn't interact well with the `--peopleinclude` parameter. When I run this, the program complains about unexpected columns in the kinship matrix (starting at one column greater than the number of females in the dataset). I suspect what is happening is that rvtests is assuming that the kinship matrix applies to the sample set *after* filtering, rather than applying to the total sample set. This sets off a few warning bells in my head, because it means that rvtest is probably ignoring the column/row labels in the kinship matrix that it has generated.

How can I work around this problem?

# Answer
> 2 votes
Adam Locke (a collaborator of this project) suggests that removing covariate information for the unselected individuals (i.e. setting it to NA) works around this problem:

> I believe the problem is that he is using a pre-computed kinship matrix including both males and females, and when using the “—peopleIncludeFile” it can’t properly select the right people from the kinship matrix. I don’t know that others have tried to do analysis this way. An alternative, that I think has worked for others, is to generate phenotype variables that have values only for males or females (so, e.g., BMI, BMI\_females, BMI\_males). For whatever reason, this appears to correctly only use the and select the samples desired without the kinship matrix.

Presumably a kinship matrix that is generated for just the subset individuals will work as well, although I don't like the idea that rvtests isn't able to do a label lookup on the kinship matrix to find the correct individuals.

---
Tags: microarray
---


# Question
Title: How do I compare measured ratios to calculated ratios in peptide mass spectrometry?
I have two different sets of internal standards in a sample (1), and a sample (2) with these internal standards as well as about 30 peptides. Each peptide has one of each of these kind of standards. From sample 1, I can get the true ratio between standard 1 (S1) and standard 2 (S2), and from sample 2 I can get the ratio between the analyte (L) to each of the standards:

Sample 1: S1/S2

Sample 2: L/S1, L/S2

I can thus calculate the ratio S1/S2 from the measured ratios in sample 2. Each sample is run in 5 technical replicates (that is; 5 injections into the LC-mass spectrometer), so S1/S2 ratios for one particular analyte could look like this (mock data):

```
measured    calculated
0.967       0.987
1.007       0.967
1.044       1.012
1.041       1.025
1.048       1.046

```

I want to compare the measured ratios to the calculated, to see if there are any bias in my acquisition of data. As the calculated and the measured ratio are from different samples, the replicates aren't "linked" so I could not use for example Bland Altman analysis. Should I use some kind of ANOVA? Or do you have any other recommendation?

I would also like to make an overall assessment of the ratios of all 30 peptides to see if the bias is different for different peptides. How could I do that?

Any input would be greatly appreciated.

# Answer
From a statistical perspective ratios are icky because the variance of the sampling distribution of a ratio decreases with the size of the numerator and denominator being sampled to create the ratio. If you can get at the numerator and denominator separately, then you should absolutely do that and look at count based tests like the Chi-square test, Fisher's exact test, or Poisson regression.

That said, a t-test would be appropriate in the two-sample case if its assumptions are met. You need to check that both samples are normally distributed, ideally via a qq-plot. You also need to check the variance of both samples and choose a variant of the t-test that does or does not assume equal variances.

If the samples are not normally distributed, then the next approach would be a Mann-Whitney U test.
> 2 votes

---
Tags: proteins
---


# Question
Title: Can a customized GRCh38 .gtf file be used with any of the GRCh38 released patches?
I got a customized GRCh38.79 .gtf file (modified to have no MT genes) and I need to create a reference genome out of it (for 10xGenomics CellRanger pipeline). I suspect that the .79 part is the Ensembl number, which according this ensembl archive list is paired with the GRCh38.p2 patch.

Should I use this patch's fasta file, or it would be fine to use any of the GRCh38 patches?

# Answer
You can see the exact assembly by looking in the README file in Ensembl's genome download page. ftp://ftp.ensembl.org/pub/release-89/fasta/homo\_sapiens/dna/README

As you can see the current assembly is GCA\_000001405.25. Ensembl 79 uses version GCA\_000001405.17. You can expect some differences between assembly versions.

If you take a look at https://www.ncbi.nlm.nih.gov/assembly?term=GRCh38&cmd=DetailsSearch

You can see that 20,466,394 bp have been added since 17 -\> 25. This is only about 0.6% change so I would guess there is not that much difference between these versions. However, if I were you I would use the corresponding version of Ensembl. Even better would be using the latest version of Ensembl, as not only the genomes primary sequence but also the positions of genes can change between versions (https://doi.org/10.1093/bib/bbw017).
> 2 votes

---
Tags: human-genome, reference-genome, gtf, 10x-genomics
---


title: Why Bioconductor?

text:
# Question
Title: Why Bioconductor?
What are the advantages of having Bioconductor, for the bioinformatics community?

I've read the 'About' section and skimmed the paper, but still cannot really answer this.

I understand Bioconductor is released twice a year (unlike `R`), but if I want to use the latest version of a package, I'll have to use the dev version anyway. A stamp of *approval* could be achieved much easier with a tag or something, so it sounds just like an extra (and unnecessary) layer to maintain.

Related to this, what are the advantages as a developer to have your package accepted into Bioconductor?

# Answer
> 18 votes
**Benefits of central repository for Community**

Having a central repository for packages is very useful. For couple of reasons:

* It makes very easy to resolve **dependencies**. Installing all the dependencies manually would be exhausting but also dangerous (point 2).
* Package **compatibility**! If I install package with dependencies, I would like to be sure that I install correct versions of all the dependencies.
* **Reliability** thanks to unified and integrated testing. `Bioconductor` is trying really hard to force developers to write good test, they also have people manually testing submitted packages. They also remove packages that are not maintained. Packages in `Bioconductor` are (reasonably) reliable.

In the end, installing dev versions of R packages is in my opinion **very bad practise** for reproducible science. If developers delete GitHub repo, commit hash you have used won't be enough to get the code.

**Benefits of central repository for developers**

I forgot about the advantages for you as developer to submit your package to `Bioconductor`:

1. Your package will be more visible
2. users will have a guarantee that your code was checked by third person
3. Your package will be for users easier to install
4. Your package will be forced to use standardized vignettes, version tags and tests -\> will be more accessible by community to build on your code

**Bioconductor specific advantages over CRAN**

I see the big advantage in the community support page, provided by `Bioconductor`. @Llopis' comprehensive elaboration.

# Answer
> 10 votes
Here is a list of the advantages of having Bioconductor for the bioinformatic community:

* **Outreach**: You have a **repository for the field**, in that language.

Some packages related to bioinformatics (in R) are distributed through personal repositories, CRAN, github, bitbucket, sourceforge, but they are less used and harder to find.
There are such efforts in other languages too: Biopython, Bioperl, Biojava, ...
Also is harder to find the repositories related to a subject in CRAN, you don't have the BiocViews, the equivalent is optional and not usually filled, which is quite useful when looking for a method.

* **Quality**: In Bioconductor each package is tested in Linux, Windows, and iOS, to make sure it works in all major operative systems (with all the dependencies).

In some rare cases like this one, a packages is not supported for certain platform, but you can known it by checking the build report.
You are required to provide a vignette and examples in every exported element (and the vignette, examples and tests should pass). You are required to be able to install the package with stricter quality than CRAN, because there is a manual review (They pointed out a *comment* in one of my functions!).
They also provide docker images of the base packages. You don't need to install the latest R version to use the Bioconductor! But developers do so (at least when checked by Bioconductor servers) to ensure that the package will keep working in next R release.

* **Reusing**: Bioconductor provides the basic elements to a big number of applications.

For example the summarizedExperiment class is provided so that any package that needs a similar object can (should) use it. Or GSEABase is the base package to deal with GSEA enrichment analysis, providing functions, methods and classes for gene sets, and collections, making easier for anyone to create their own GSE analysis.
It is easier to build upon the work of others if you know you are following the same quality standards.

* **Support**: To support Q&A is mandatory, the package maintainer must be registered in the webpage.

While in CRAN usually the support is given by each package in its own way, in Bioconductor you can directly reach the maintainer and the users by posting in the same central place.

# Answer
> 3 votes
Regarding what the advantage is to you as a developer of having a bioconductor package rather than using CRAN:

1. There's a hierarchy of package quality, with bioconductor on the top (followed by CRAN and then "random github repos"). While there are many excellent CRAN packages, the average bioconductor package is better tested and documented. So if I as a user have two different packages that I can use and one is on bioconductor and the other CRAN, then I use the bioconductor package.
2. Higher visibility. Since bioconductor packages are held in higher regard, they also become more visible. Further, the different "views" (e.g., "Alternative Splicing" and "Transcription") make it convenient to find relevant packages. You can always search CRAN, but allowing tagging like this aids in discovery.

# Answer
> 2 votes
The most important reason is that Bioconductor has growing set of common data structures and base packages. If package X and package Y needs to work with the same type of data, having a common data structure in the core Bioconductor package Z makes our lives so much easier. I could do something in package X, take out the results and keep working on my data with package Y. This works because I'm using versions of package X and Y that are compatible with the data structure defined in package Z. The Bioconductor team makes sure that all packages use the common data structures and already existing packages where possible, so that we're all on the same page and so that people don't reinvent the wheel again and again.

Also, there's a sophisticated process around getting your package accepted to Bioconductor. This ensures that packages use data structures and functions already available in other Bioconductor packages (where reasonable), it ensures that packages are well written, that they have good documentation, and that they are well tested.

---
Tags: r, bioconductor
---


# Question
Title: What kind of "gff" format does bioawk parse?
I was wondering if I could use the gff parsing capability of bioawk to facilitate the parsing of gtf files, and I looked at the following help message:

```
$ bioawk -c help
bed:
    1:chrom 2:start 3:end 4:name 5:score 6:strand 7:thickstart 8:thickend 9:rgb 10:blockcount 11:blocksizes 12:blockstarts
sam:
    1:qname 2:flag 3:rname 4:pos 5:mapq 6:cigar 7:rnext 8:pnext 9:tlen 10:seq 11:qual
vcf:
    1:chrom 2:pos 3:id 4:ref 5:alt 6:qual 7:filter 8:info
gff:
    1:seqname 2:source 3:feature 4:start 5:end 6:score 7:filter 8:strand 9:group 10:attribute
fastx:
    1:name 2:seq 3:qual 4:comment

```

I see that there are 10 fields defined for gff parsing.

However, when I look at the pages for gtf and gff in the ensembl website (gff/gtf and gff3), all have 9 fields.

I'm curious about these "filter" and "group" fields. What are they meant for? Are they extracted from some of the columns mentioned in the above pages?

# Answer
I would consider the description there a bug. The `filter` is actually the strand, `strand` is the frame, `group` is the attribute, and `attribute` does nothing. These are really meant to be the 9 columns.

**Edit**: There's a bug report related to this.

**Edit 2**: I've made a pull request to clarify this and fix the aforementioned bug report.

**Edit 3**: I realized that I never directly answered the title of your question (*mea culpa*). bioawk itself will work with gff, gff3, or gtf files. It really is just treating them as tab-separated files with named columns (this is surprisingly convenient, since it's a PITA to remember what column does what).

**Edit 4**: The PR has been merged. If you install from github then you'll see corrected field names.
> 5 votes

---
Tags: file-formats, gff3, gtf, bioawk
---


# Question
Title: How to correct alpha, and not p-values themselves, for visualization purposes
I have a set of differentially methylated/expressed/whatever entities with p-values attached (example below).

```
entity_name    p-value    magnitude
entity1        0.04459    0.68
entity2        0.02283    0.99
...
entity_n       0.78       0.025

```

Typically, I apply the p.adjust function in R with the "fdr" (Benjamini-Hochberg) approach to leave me with p-values adjusted to control the FDR.

```
adjusted <- p.adjust(mydata,"fdr")

```

However, I am interested in showing a volcano plot with the unadjusted p-values, and two alpha levels: 0.05 and one that corresponds to the correction. What is the best way to get this alpha? Is it appropriate to set the "corrected alpha" to the lowest original p-value that doesn't pass FDR correction?

# Answer
The only way to get the alpha levels is to determine what they will be with `p.adjust()`, since they will depend on the distribution of your unadjusted p values. The general steps you should follow will be:

1. Add a column of adjusted p-values to your dataframe (`mydata$padj = p.adjust(mydata, method="BH")`, which is the same as FDR and saves a character).
2. Use `which` and `max` to determine your two alpha threshold (e.g., `max(mydata$pvalue[mydata$padj < 0.05])`

Then you can adjust your plots however you like (presumably with some horizontal lines at the various alphas). Whether you take the smallest non-significant value or the largest significant value is up to you, just describe what "dots on the line" represent.
> 7 votes

# Answer
When you're doing a p-value adjustment, the same unadjusted p-value in different genes can be given different adjusted p-values depending on other factors. That means that you can't directly draw a line associated with the FDR on a plot of unadjusted p-value.

One possibility would be to take a range of values that are close to the FDR threshold (e.g. the 20 values closest to threshold), and draw a p-value greyzone within that region:

```
#!/usr/bin/Rscript
values <- c(rnorm(10000),rnorm(100, mean=1.5));
val.mean <- median(values);
val.diffs <- abs(values - median(values));
val.reldiffs <- (values - median(values));

val.pval <- pnorm(val.diffs, mean = mean(val.diffs),
                   sd=sd(val.diffs), lower.tail=FALSE);
val.padj <- p.adjust(val.pvals, method="BH");
fdr.threshold <- 0.1;
close.bh <- order(abs(val.padj - fdr.threshold))[1:20];

png("SE.663.png");
plot(val.reldiffs, -log10(val.pval),
     col=ifelse(1:10100 <= 10000,"darkblue","darkgreen"));
abline(h=-log10(0.05), col="red");
text(0,-log10(0.05),"p=0.05", pos=1);
abline(h=range(-log10(val.pval[close.bh])), col="#00000040", lty="dashed");
rect(xleft=min(val.reldiffs)*2, xright=max(val.reldiffs)*2,
     ytop=max(-log10(val.pval[close.bh])),
     ybottom=min(-log10(val.pval[close.bh])), col="#00000020", border=NA);
text(0,min(-log10(val.pval[close.bh])),"FDR=0.1", pos=1);
invisible(dev.off());

```
> 3 votes

---
Tags: r
---


# Question
Title: Finding residues physically interacting in a PDB structure and sequence related to that structure
I am going to describle the problem completely. First I should say that there are two kinds of interactions: transient and stable interactions. We are supposed to work on stable interaction like interactions between two monomers in a heterodimer. In a heterodimer there are two chains in which there are some residues between monomers which take part in the interaction of monomers and building the structure of heterodimer and also there are other residues interacting inside each monomer for building the structure of that monomer. There are some PDB IDs related to heterodimers in PDB. Now imagine that we are trying to show which residues are physically interacting between two monomers in a heterodimer and then we are trying two show the site of these physically interacting residues in the sequence related to that special structure. For that, we need to find the sequence related to that structure and we can download it from PDB. But if we want to make a multiple sequence alignment of that sequence, first we need to find the refseq sequence of that and then we must blast the sequence and find homologues of that sequence. Here there are some challenges: first, we must find residues interacting in the structure of heterodimer between two monomers. Second, we must find the sequences of monomer chains of the heterodimer and map the interacting residues in the structure to the sequence. Third, we know that when a structure is solved and its monomers are sequenced, may be they are not able to completely sequence the heterodimer and residue numbers are maybe different from the related refseq sequence. I would like to know how can we find the related refseq sequence (or sequences, in the case of heterodimer) and then how can we map the physically interacting residues in the structure of PDB ID to refseq sequence (finding the sites of these interacting residues in that refseq sequence) regarding three challenges mentioned above. Imagine that I have 100 PDB IDs of 100 nonredundant heterodimers and I would like to find and download the structures, sequence and finally refseq sequence. What should I do?

# Answer
My reputation isn't high enough to make a comment, but have you considered a modelling software such as PyMol or Chimera? You can execute a command line script with softwares such as these and these have the functionality to predict H-bonds.

In the PDB files there is information regarding the GenBank and UNIPROT info of the protein, and you can find more information here:

http://www.wwpdb.org/documentation/file-format-content/format33/sect3.html.

If you want RefSeq ID, you can use GenBank info you extract and convert by piping into a biomaRt (Bioconducor) script, although maybe there is a more direct way.

As for finding which residues in the sequence are interacting, I'm fairly certain that the software packages like PyMol or Chimera will give you information about the interacting residues.

From what I understand, you will want to filter out only the residues that interact and are on different chains (e.g, Chain A vs. Chain B). This process will require you to work with the syntax of the different software packages when you are scripting. For example, Chimera has a basic primer here:

https://www.cgl.ucsf.edu/chimera/current/docs/ProgrammersGuide/basicPrimer.html

Note that I only gave a reference to Chimera because I am more comfortable with that then other packages and there is likely a more efficient way to do this. With Chimera, you can find the position within the chain of residues participating in the interaction and then you'd store these in a file to do whatever you want later.

I hope this helps.
> 2 votes

---
Tags: protein-structure
---