text:
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:

RNA sequence
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 :
.(((.(..(((((((((.(((((..(((((.(((((((((....)))))))))..)))))....((((((.....(((.....(((....))).....)))..)))))).))))))).))))))).....).))).

I suggest you have a look at Biopython. They have a PDB module: http://biopython.org/wiki/The_Biopython_Structural_Bioinformatics_FAQ
I'm not familiar with PDB, so I don't know if this does what you want.

Hi Peter, thanks for posting your question, and welcome to Bioinformatics Stack Exchange. I have fixed up the grammar and spelling of your question to make it easier for other people to read and answer. I am not familiar with the PDB file format, so it would be great if you could add in a bit more context / story around your question. For example, where are these PDB files from (i.e. public, or private data source)? Is there a reason why taking transcript sequences from NCBI is not appropriate?

Could you give an example of a specific PDB file we can test on?

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)

I added in the example from that page for getting transcripts, but couldn't [quickly] find anything about bracket notation

@gringer yeah, me neither (and that documentation has some pretty bad suggestions like for f in \ls *txt`; do . . .`(!) so I would take it with a grain of salt anyway. Thanks for the edit though!

what's wrong with for f in ls *txt; do . . . ? I can fix it but why?

@MarcinMagnus two things: first, it will fail on weird file names which is why it is generally bad practice to parse the output of ls. ls is designed to produce output for humans, not machines and can't deal with odd filenames gracefully (although recent versions of GNU ls have improved on that). Second, it's just inelegant since it isn't needed here. A simple for f in *txt will do the same thing and is far more robust since it can deal with any file name, even those containing newlines.

Also, I apologize for my tone. The parsing ls thing is a classic meme in the *nix world and one of the more common pitfalls for novice shell programmers. I spend a lot of time in the *nix communities (I am a mod on the [linux.se] and [ubuntu.se] sites) so I am probably more sensitive about this than I need to be. In any case, I didn't mean to disparage your work and my apologies if it came across that way.

OK, I see it now ;-) I knew about abuse of cat, and echo. Now I see ls on the list http://porkmail.org/era/unix/award.html I didn't know that before. Thanks for pointing this out. I will fix the docs in all places ;-) thanks! I'm using Linux for more than 10y but still some things remind to be learned ;-)

@MarcinMagnus thanks :) While you're at it, it would be great if you could also fix the missing quotes (for i in *; do command "$i"; done) since those are also essential when dealing with arbitrary file names. For more than you even wanted to know about the importance of quoting variables in *nix shells, see Security implications of forgetting to quote a variable in bash/POSIX shells and Why does my shell script choke on whitespace or other special characters?.

I'll throw another suggestion into the mix regarding backticks, which can be replaced in sh scripts with $() (which is nestable), e.g. for i in $(ls *pdb) [but as per @terdon's comment the ls is redundant, so for i in *pdb would work better]

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

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.


text:
samtools depth print out all positions
I am trying to use samtools depth (v1.4) with the -a option and a bed file listing the human chromosomes chr1-chr22, chrX, chrY, and chrM to print out the coverage at every position:
cat GRCh38.karyo.bed | awk '{print $3}' | datamash sum 1
3088286401

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

I tried it for a few bam files that were aligned the same way, and I get differing number of entries:
3087003274
3087005666
3087007158
3087009435
3087009439
3087009621
3087009818
3087010065
3087010408
3087010477
3087010481
3087012115
3087013147
3087013186
3087013500
3087149616

Is there a special flag in samtools depth so that it reports all entries from the bed file?
If samtools depth is not the best tool for this, what would be the equivalent with sambamba depth base?
sambamba depth base --min-coverage=0 --regions $bedfile $inputfile

Any other options?

Out of interest, does -aa at least give the same number of bases each time? -a -a (or -aa) output absolutely all positions, including unused ref. sequences
There is also an old but relevant discussion about -a and BED behaviour on github/samtools/#374

I tried with -a and with -aa, and given a bed file, it seems to produce the same number of entries.

Two samtools commits from May 2016; Added 30ish new tests for depth/mpileup -a and -aa and Depth/Mpileup -a fixes. Fixes #374 are somewhat enlightening. My gut feeling is that although possible, using -a/-aa with a BED file is not entirely supported with samtools depth (i.e. has potentially undefined behaviour) and perhaps you would want to investigate another tool as you've suggested?

I dont work with human genomes, so excuse my naive question...why do you need to provide a bed file if you want every position?

If every reference defined in the BAM file is touched by at least one read, then -a and -aa will produce the same output. The second option, -aa is useful if you want to produce zero-depth information for sequences that are completely uncovered.

@gringer Yes, the reason why I use GRCh38 bed file is to have the depth values for only the karyotypic chromosomes, rather than the scaffolds and decoys.

... except it appears that that there's a bug in samtools around this (i.e. I just noticed Sam Nicholls' comment).

You might try using bedtools genomecov instead. If you provide the -d option, it reports the coverage at every position in the BAM file.
bedtools genomecov -d -ibam $inputfile > "${inputfile}.genomecov"

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

Just replace the -a option with -aa:
samtools depth -b $bedfile -aa $inputfile

I see that you're using the GRCh38 human reference genome build, which includes alternate scaffolds that represent a wider variety of genomic variation in the human genome. This was quite a substantial change from previous reference genome builds, and many bioinformatics tools still haven't been updated to accept the GRCh38 reference.
The discrepancy in entry size is probably because there are alternate contigs that are uncovered by the reads that have been mapped. This is expected, but can cause a bit of confusion.

Yes, the reason why I use GRCh38 bed file is to have the depth values for only the karyotypic chromosomes, rather than the scaffolds and decoys. The -aa option gave the same output as the -a option.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/

sort: votes


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

dna_sm - Repeats soft-masked (converts repeat nucleotides to lowercase)
dna_rm - Repeats masked (converts repeats to to N's)
dna - No masking

Second part options:

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

Right now I usually use a non-masked primary assembly for analysis, so in the case of humans:
Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
Does this make sense for standard RNA-Seq, ChIP-Seq, ATAC-Seq, CLIP-Seq, scRNA-Seq, etc... ?
In what cases would I prefer other genomes? Which tools/aligners take into account softmasked repeat regions?

Relevant blog post: http://genomespot.blogspot.ch/2015/06/mapping-ngs-data-which-genome-version.html

What kind of "alignments"? Protein tblastn? Whole genome alignments? NGS read alignments? Gene-level alignments?

Generally, you should use the soft-masked or unmasked primary assembly. Cross-species whole-genome aligners, especially older ones, do need to know soft-masked regions; otherwise they can be impractically slow for mammalian genomes. Modern read aligners are designed to work with repeats efficiently and therefore they don't need to see the soft mask.
For GRCh38, though, I would recommend to use the official build at GRC FTP. Most people will probably choose "no_alt_analysis_set". Using the Ensembl version is discouraged due to its chromosome naming. We more often use "chr1" instead of "1" for GRCh38. At one point, Ensembl actually agreed to use "chr1" as well, but didn't make that happen due to technical issues, I guess.
As to alternate haplotypes, most aligners can't work with them; no variant callers can take the advantage of these sequences, either. When you align to a reference genome containing haplotypes with an aligner not supporting these extra sequences, you will get poor mapping results.

There's rarely a good reason to use a hard-masked genome (sometimes for blast, but that's it). For that reason, we use soft-masked genomes, which only have the benefit of showing roughly where repeats are (we never make use of this for our *-seq experiments, but it's there in case we ever want to).
For primary vs. toplevel, very few aligners can properly handle additional haplotypes. If you happen to be using BWA, then the toplevel assembly would benefit you, but only if you use a dedicated wrapper to handle the ALT information, see bwakit. If you use BWA (bwa-mem) right from the command line without this wrapper then do not use the toplevel assembly. For STAR/hisat2/bowtie2/BBmap/etc. the haplotypes will just cause you problems due to increasing multimapper rates incorrectly. Note that none of these actually use soft-masking.

Are there any de bruijn graph type read aligners that use haplotype info?

Not that I know of at least.

Which tools/aligners take into account softmasked repeat regions?

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

TOPLEVEL
These files contains all sequence regions flagged as toplevel in an Ensembl
schema. This includes chromsomes, regions not assembled into chromosomes and
N padded haplotype/patch regions.
E.g: I used the soft masked assemblies for genome annotation pipelines like MAKER, also toplevel unmasked ones for RNA-seq, ChipSeq analysis
PRIMARY ASSEMBLY
Primary assembly contains all toplevel sequence regions excluding haplotypes
and patches. This file is best used for performing sequence similarity searches
where patch and haplotype sequences would confuse analysis.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Publicly available genome sequence database for viruses?
As a small introductory project, I want to compare genome sequences of  different strains of influenza virus.
What are the publicly available databases of influenza virus gene/genome sequences?

There area few different influenza virus database resources:

The Influenza Research Database (IRD) (a.k.a FluDB - based upon URL)

A NIAID Bioinformatics Resource Center or BRC which highly curates the data brought in and integrates it with numerous other relevant data types

The NCBI Influenza Virus Resource

A sub-project of the NCBI with data curated over and above the GenBank data that is part of the NCBI

The GISAID EpiFlu Database

A database of sequences from the Global Initiative on Sharing All Influenza Data. Has unique data from many countries but requires user agree to a data sharing policy.

The OpenFluDB

Former GISAID database that contains some sequence data that GenBank does not have.

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

Virus Pathogen Resource (VIPR)

A companion portal to the IRD, which hosts curated and integrated data for most other NIAID A-C virus pathogens including (but not limited to) Ebola, Zika, Dengue, Enterovirus, and Hepatitis C

LANL HIV database

Los Alamos National Laboratory HIV database with HIV data and many useful tools for all virus bioinformatics

PaVE: Papilloma virus genome database (from quintik comment)

NIAID developed and maintained Papilloma virus bioinformatics portal

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

Glad that my upvote allowed you to expand your answer with more links :) I hear ViralZone http://viralzone.expasy.org/ is decent too (although I have not used it myself)

Thanks Chris! I agree Viral Zone is GREAT! I think it is one of the best virus education site around!

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

Hi quantik, thanks for your answer and welcome to Bioinformatics Stack Exchange. As a new user with reputation under 50 it's not possible yet for you to comment on other questions or answers (as @zx8754 has recommended), so in this case your answer was the most appropriate thing to do. If you want to improve on someone else's answer in the interim, feel free to edit your answer and incorporate a summary of the other answer(s) that you think are appropriate, making sure to include attribution where appropriate.

Thanks quantik...nice addition. A woeful oversight on my part as I am part of the group that develops PaVE. :-)

Influenza virus resource at NCBI or FluDB.

These are both included above. FluDB is actually the NIAID Influenza Research Database (IRD) Bioinformatics Resource Center (BRC) that has fludb.org as a URL.

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


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

Loads :) Could you please [edit] your question and tell us i) how many sequences you are talking about; ii) whether or not you already have the human protein sequences or just their names; iii) if you are just looking for simple sequence homology or if you will use the homology to infer a functional homology (in which case protein domains should be taken into account); iv) why 13-17aa specifically? v) what microbial species? Any? Specific ones?

Hi bluescholar1212, thanks for your question and welcome to the Bioinformatics Stack Exchange. Bioinformatics is a large area, and can include a number of different programs that could make your life easier, almost all of which are unrelated to microbial homology. Being more specific in the questions that you ask can be a great help to answerers, as it allows them to answer without going down the wrong path. What sort of output are you interested in? You've mentioned a problem exists, but haven't stated what that problem is. What is your story around the problem that you are trying to solve?

Welcome to Bioinformatics.SE! Just to split some hairs :)…you'll want to be careful about using terms like high homology. Sequences are either homologous (share ancestry) or not. Sequence similarity is a commonly-used proxy for homology, & it is appropriate to say high similarity. But there is no sliding scale of homology. Some pairs of sequences are so similar that the only feasible explanation is shared ancestry. Some are so different that it's clear there is no homology. There is some grey area in the middle, but this refers only to our uncertainty. Hope this makes sense!

I agree with @daniel-standage, there is no high degree of homology, or high homology segments. The segments are homologous or not. It is like saying high degree of pregnancy. Maybe you can replace the terms by something like "high probability of homology"?

@DanielStandage Since we are splitting hairs, all sequences have a shared ancestry if you go back far enough, if you accept that there was a single origin of life

@Chris_Rands True enough. :-)

@DanielStandage Duly noted! That is a good and important point. I am interested in regions with a high percentage of sequence similarity, which does not necessarily entail homology as you pointed out.

@bluescholar1212 I have then edited the question. Please feel free to revert back or modify as you wish if you think I misrepresented your question.

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:

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:

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

thank you for your helpful answer! BLAST is one of the tools that I was considering. I hadn't looked into the PSI-BLAST algorithm as an option and it sounds like exactly what I am looking.
In addition to running the search, which I will plan on doing with blast locally, I would like to automate the curation of the results. Would my best bet for doing this be to write a basic script around blast specifying which results to keep and filtering the rest?

@bluescholar1212 yes, probably. That would also be on topic here, by the way. Once you have your results, ask another question and we'll be happy to help you parse them. You might also want to ask a new question once you decide which of these tools you use, explaining what output you want to eventually keep and asking about what output format top use. Blast can return various output formats depending on the options you use it with and some are easier to parse than others.


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

Welcome to Bioinformatics.SE! I'm afraid in its current state this question is very broad and may not be a good fit for this site in its current state.

@DanielStandage thanks! ooo! ok understand! hopefully someone helps, since there is a big curiosity about this :D

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.

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:

Start with a fully unfolded chain (like a metal chain with its ends pulled).
Move a part of the chain to create a new shape.
Calculate the energy of the new shape.
Accept or reject the move depending on the change in energy.
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.

wow that's surprising. I would have expected the protein to start folding as soon as one end of it is able to, which is probably the first end to leave the ribosome. The last amino acid probably hasn't even been attached yet, let alone available for moving and calculating the energy :)

That does happen to some degree, but the Rosetta [computer] model doesn't do this. A whole bunch of odd things happen with protein folding. For example, there are chaperone proteins that encourage folding in a particular fashion to get through low-energy humps to reach a higher energy state. The foldit researchers have found a way to model chaperones using a technique called "humans".


text:
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

What is the .number that are added to the Gene e.g. the .2 in ENSG00000242268.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.

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.

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.

Tools such as DESeq2 are able to generate other expression statistics that are similar to, but not the same as, FPKM. I'm not sure if it makes sense to answer this for the specific question (i.e. because of FPKM normalisation), or in the more general case, where I would say that the values are floating point due to various different normalisation procedures (e.g. correction for transcript length, sequencing depth, shot noise, confounding covariates).


text:
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:

Build an index for the FASTA file
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

Does this accomplish what you're looking for: https://github.com/dib-lab/khmer/blob/master/examples/python-api/exact-counting.py

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)

Not exactly what I'm looking for. EDITED.

See this question if you run into problems not counting kmers above 255

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


text:
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?

Are you only asking about a format describing genomic features? "Annotation" is a very broad term but it looks like you're only considering genomic regions here or, at least, things that have i) a defined "region" and ii) a defined "function". That would still exclude phenotype annotations for proteins or GI annotations for genes etc. Could you [edit] and clarify what kind of "annotations" you are considering?

The BED concept of autoSql is a pretty nice feature of an annotation format and allows for a lot of extensibility. The concept of feature hierarchy is still basically single level though

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

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

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.

Regarding point 1, I would add a "line-based" further requirement, not just "text-based". It is also easier to manipulate if all the information is organized in tab-separated columns, but this may reduce the flexibility compared to what is currently allowed within the attributes column of gff format.

@bli Excellent suggestion! I've made a note of this. Flexibility is a double-edged sword; once you can do almost anything you can't reliably depend on a format :(

I hadn't realized there existed a record separator ASCII character. I think tab has the merit of making files more human readable (with current text-based computing tools).

I’m honestly not sure whether human readable is an objective criterion for goodness. In fact, it seems to get in the way so there’s at least an argument to be made for a trade-off.

@KonradRudolph Sure, all of the tabix work-arounds can be considered bandages over the "human readable" part.

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.

I guess that as long as one only does very simple things with the annotation information, one (=me) is happy with line-based format that can be processed with common command-line tools. However, your point of view is interesting and you seem to have a much better understanding of these issues than me. I'll try to get at least an idea of what "impedance mismatch" is.

Totally agree, flat files are not the way forwards for complex data structures, although I am guilty of this. Any thoughts on SQL vs NoSQL?

@Chris_Rands I have no experience working with NoSQL but since the data here is specifically relational, I don’t think NoSQL would offer any specific benefits.

Seems more like "nosql" rather than relational in some case. Just look at the mess that Chado is, with hundreds of "linking tables"

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:

chr, required
start (0-based), required
end, required
strand
ID
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.

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.


text:
How can I build a protein network pathway from a gene expression quantification file?
Assume I have found the top 0.01% most frequent genes from a gene expression file.
Let's say, these are 10 genes and I want to study the protein protein interactions, the protein network and pathway.
I thought to use string-db or interactome, but I am not sure what would be a plausible way to approach this problem and to build the protein network etc. Are there other more suitable databases?
How can I build a mathematical graph or network for these data?

There are multiple ways to do this, and multiple protein interaction databases besides the ones you mentioned, such as BioGRID or IntAct. Interaction databases are different in how interactions are defined, sometimes it can be experimental evidence of interaction, sometimes coexpression, orthology-based predictions, etc.
There is no single solution to your problem. For String-DB you can use their R package STRINGdb.
# string-db
library(STRINGdb)

# retrieve full graph for H. sapiens
string_db <- STRINGdb$new(version="10", species=9606,
                          score_threshold=400, input_directory="" )

# define genes of interest
genes.of.interest <- ...

# get their neighbors
neighbors <- string_db$get_neighbors(genes.of.interest)

# get the subgraph of interest (your genes + their neighbors)
my.subgraph <- string_db$get_subnetwork(c(genes.of.interest, neighbors))

# look how many genes and interactions you have
my.subgraph

If you prefer GUI, you might consider using Cytoscape. They have multiple addons for interaction databases, e.g. for STRINGdb. There you can just provide your gene list and get a network in a few clicks.

Note that STRING is very dirty. if you want to use it toi build a network it is essential to filter the results by detection method (it has many, many inferred interactions).

Many interaction databases now work with PSI format files. Most of the main databases can do this and the EBI has set up PSICQUIC View, a very useful page where you can query multiple databases at once.
Note that it is very important to limit the results according to the detection method. There is a lot of noise in protein interaction databases. Depending on what you want to do you could limit to only experimentally verified interactions or to only direct, binary interactions
(so exclude the results of, for example, ChIP analyses which can also find complexes) etc.
That said, here's a simple example script that will query the APID database using its PSICQUIC service:
#!/usr/bin/perl
use strict;
use warnings;
use LWP::Simple;

my @proteins;
## Read a list of target proteins, one per line (this expects UniProt names)
open(my $fh, "$ARGV[0]") or die "Need a list of proteins as the 1st argument: $!\n";
while (<$fh>) {
    chomp;
    push @proteins, $_;
}
close($fh);
## Get human interactions only
my $species="9606";

## Get the interactions for each target protein
foreach my $protein (@proteins) {
    my $queryUrl= "http://cicblade.dep.usal.es/psicquic-ws/webservices/current/search/query/$protein";
    $queryUrl  .= "?taxidA:$species%20ANDtaxidB:$species";

    my $tries=1;
    my $content = get $queryUrl;
    while ($tries<=10) {
        if (defined($content)) {
            $tries=11;
        } else {
            print STDERR "Could not retrieve $queryUrl, retrying($tries)...\n";
            $content = get $queryUrl;
        }
        $tries++;
    }

    # Now list all interactions
    my @lines = split(/\n/, $content);
    my $LINES= @lines;
    my $count = 0;
    for my $line (@lines) {
        $count++;
        my @flds = split(/\t/, $line); # split tab delimited lines
        # split fields of a PSIMITAB 2.5 line
        my ($idA, $idB, $altIdA, $altIdB, $aliasA, $aliasB, $detMethod, $author, $pub, $orgA, $orgB, $intType, $sourceDb, $intID, $conf) = @flds;
        ## Here you can add logic to limit the interactions by specific detection method codes ($detMethod)
        ## or database of origin etc.

        ## Print
        print "$line\n"
    }
}

I tested it by using a file with the UniProt ID for human TP53 (P04637) and it returns a list of 3988 interactions (showing 10 randomly selected results below):
$ foo.pl names.txt | shuf -n 10
uniprotkb:Q9BWC9    uniprotkb:P04637    -   -   uniprotkb:CCDC106(gene_name)    uniprotkb:TP53(gene_name)   psi-mi:"MI:0007"(anti tag coimmunoprecipitation)    Zhou, J. et al.(2010)   pubmed:20159018 taxid:9606(Homo sapiens)    taxid:9606(Homo sapiens)    -   psi-mi:"MI:0469"(intact)    intact:EBI-7812926  -
uniprotkb:P04637    uniprotkb:O15126    -   -   uniprotkb:TP53(gene_name)   uniprotkb:SCAMP1(gene_name) psi-mi:"MI:0018"(two hybrid)    Lim, J. et al.(2006)    pubmed:16713569 taxid:9606(Homo sapiens)    taxid:9606(Homo sapiens)    -   psi-mi:"MI:0463"(biogrid)   biogrid:720622  -
uniprotkb:P63165    uniprotkb:P04637    -   -   uniprotkb:SUMO1(gene_name)  uniprotkb:TP53(gene_name)   psi-mi:"MI:0018"(two hybrid)    Minty, A. et al.(2000)  pubmed:10961991 taxid:9606(Homo sapiens)    taxid:9606(Homo sapiens)    -   psi-mi:"MI:0463"(biogrid)   biogrid:262339  -
uniprotkb:P04637    uniprotkb:P31350    -   -   uniprotkb:TP53(gene_name)   uniprotkb:RRM2(gene_name)   psi-mi:"MI:0416"(fluorescence microscopy)   Xue, L. et al.(2003)    pubmed:12615712 taxid:9606(Homo sapiens)    taxid:9606(Homo sapiens)    -   psi-mi:"MI:0465"(dip)   dip:DIP-40167E  -
uniprotkb:P04637    uniprotkb:Q00987    -   -   uniprotkb:TP53(gene_name)   uniprotkb:MDM2(gene_name)   psi-mi:"MI:0004"(affinity chromatography technology)    Dai, MS. et al.(2004)   pubmed:15308643 taxid:9606(Homo sapiens)    taxid:9606(Homo sapiens)    -   psi-mi:"MI:0463"(biogrid)   biogrid:478073  -
uniprotkb:Q99576    uniprotkb:P04637    -   -   uniprotkb:TSC22D3(gene_name)    uniprotkb:TP53(gene_name)   psi-mi:"MI:0428"(imaging technique) Ayroldi, E. et al.(2015)    pubmed:25168242 taxid:9606(Homo sapiens)    taxid:9606(Homo sapiens)    -   psi-mi:"MI:0463"(biogrid)   biogrid:1255896 -
uniprotkb:P04637    uniprotkb:Q00987    -   -   uniprotkb:TP53(gene_name)   uniprotkb:MDM2(gene_name)   psi-mi:"MI:0415"(enzymatic study)   Lui, K. et al.(2013)    pubmed:23572512 taxid:9606(Homo sapiens)    taxid:9606(Homo sapiens)    -   psi-mi:"MI:0463"(biogrid)   biogrid:859223  -
uniprotkb:P25685    uniprotkb:P04637    -   -   uniprotkb:DNAJB1(gene_name) uniprotkb:TP53(gene_name)   psi-mi:"MI:0004"(affinity chromatography technology)    Qi, M. et al.(2014) pubmed:24361594 taxid:9606(Homo sapiens)    taxid:9606(Homo sapiens)    -   psi-mi:"MI:0463"(biogrid)   biogrid:938952  -
uniprotkb:P04637    uniprotkb:Q8IW41    -   -   uniprotkb:TP53(gene_name)   uniprotkb:MAPKAPK5(gene_name)   psi-mi:"MI:0424"(protein kinase assay)  Sun, P. et al.(2007)    pubmed:17254968 taxid:9606(Homo sapiens)    taxid:9606(Homo sapiens)    -   psi-mi:"MI:0469"(intact)    intact:EBI-1202077  -
uniprotkb:Q13526    uniprotkb:P04637    -   -   uniprotkb:PIN1(gene_name)   uniprotkb:TP53(gene_name)   psi-mi:"MI:0096"(pull down) Mantovani, F. et al.(2007)  pubmed:17906639 taxid:9606(Homo sapiens)    taxid:9606(Homo sapiens)    -   psi-mi:"MI:0469"(intact)    intact:EBI-6112688  -

Thank you, I am not sure what's EBI nor how to convert/map gene names to list of proteins.

@0x90 EBI is just the European Bioinformatics Institute. Mapping is always complicated, I am afraid, and one of the main practical problems for this sort of analysis is that both genes and proteins have many, many different identifiers in the different databases. However, there are many tools that can help you map. For example: 1, 2.

If you run this Rscript using the gene name as an argument, you'll get a file with the pathway written to the directory.
#!/usr/bin/Rscript
args = commandArgs(trailingOnly=TRUE)
library(paxtoolsr)
id <- args[1]
write.table(graphPc(source=id,kind='neighborhood',
            format='BINARY_SIF',verbose=TRUE),
            file=paste(id,'pathway',sep='_'), sep='\t',quote=FALSE)

Why not just use string-db's online tool? There you can adjust all parameters of interest such as interaction confidence as well as pull in extra protein interactions.
If you make a user you can also upload a custom background gene-set set which means the analysis of go-terms and interconnectedness that sting provides is usable.


text:
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)?

Please add some tags for GO, non-model organism, gaf, overrepresentation analysis etc. I just used a random tag because there are not many sensible tags in the list yet.

I modified a bit the question, added tags and I centered the question in creating the file in that format, if you disagree, just revert my changes

I added [tag:gene-ontology] and [tag:gaf]. I'm not sure whether a [tag:over-representation-analysis] or a [tag:enrichment-analysis] would be more appropriate.

@Llopis: Hopefully I didn't just accidentally overwrite your changes!

No luck, my edit was rejected (or at least doesn't appear) :\

If you want to bring questions from another site, please don't just copy/paste them verbatim. Instead, ask them again in your own words. And if you do copy, it would be better to put the quoted text in a quote block to clearly indicate that it is a quote.

@terdon, sorry need to adapt to the new rules ;) I accepted the edits, however I don't have time to rewrite the question in my own words, and I don't think it is strictly necessary. Even if that might seem a contradiction with my previous post, it is not imo, because the question is pretty aged.

@Michael no, I didn't expect you to rewrite this now! But please do so next time. Going around just pinching questions verbatim from other sites doesn't seem right. Not unless you're in a position to offer clarifications or this is an issue that affects you too.

Indeed, I can offer ample clarification for this question, we are doing exactly the same thing, except that we are using InterproScan to assign GO terms not Blast2GO as was likely used here. I have a similar script to convert IPS output to GAF.

Ha! I only now realized you were the person who first raised the whole copying from other sites issue as a problem :)

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.

I think the evidence code should be changed to IEA, because it is unlikely one would have manually annotated GO's for a non-model organism.

@Michael, it doesn't really make any difference so I chose TAS since these are annotations provided by the OP. But sure, IEA would also work. The field will be ignored in any case when doing an enrichment analysis so you can put whatever you prefer.

Sure, IEA might be even discarded by some tools.


text:
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
//

Is there any reason you can't use your preferred tool for editing protein alignments? The format is the same, I've never come across an alignment editor that can deal with protein but not nucleotide alignments.

Yes, because I do not have a preferred tool yet, and this question would broad my horizon of possibilities. Since it is possible that someone that use to edit RNAs has a preference between software because, maybe, one is more optimized than another. Since the bioinformatic is so vast and new tools are born every day, I don't see how this could be a stupid question.

Oh, hey, nobody said it was a stupid question! I was just asking why you specified RNA in case you needed a specific feature regarding it. In general, a sequence alignment is a sequence alignment and it doesn't make much difference if it is protein, DAN or RNA. That said, RALEE does have some RNA_specific features.

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!)

Ha! +1 for pure geekiness and for using the world's best editor! :P Are you sure about jalview not supporting such moving though? I haven't used it in a while but as I recall you could just select multiple sequences (and the structure line is just another "sequence") and any edit applied to one would be applied to all. Is that wrong?

Since I use emacs, I think this will be really good for me ! Thanks.

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.


text:
How to convert species names into common names?
I’m trying to find common names from a list of scientific names (not all will have them though).
I was attempting to use taxize in R but it aborts if it doesn’t find an entry in EOL and I don’t know a way around this other than manually editing the list - in which case I might as well check all this way.
Does anyone have a better way of doing this?

The answer to your question might be taxon-specific. For instance, if you deal with fish, the fishbase database has quite important common-name information. Are you dealing with a particular kind of organism?

Commons names in any language? Or specifically in English? Or another specific language?

to answer the above two, it's a mix of all types of prokaryotes. The returns would need to be in English, it seems EOL gives other languages as well (but returns nested lists of diff lengths- i gave up trying to use those).

As Pierre mentioned, NCBI is a good resource for this kind of transformation.
You can still use taxize to perform the conversion:
library("taxize")

species <- c('Helianthus annuus', 'Mycobacterium bovis', 'Rattus rattus', 'XX', 'Mus musculus')

uids <- get_uid(species)

# keep only uids which you have in the database
uids.found <- as.uid(uids[!is.na(uids)])
# keep only species names  corresponding to your ids
species.found <- species[!is.na(uids)]

common.names <- sci2comm(uids.found, db = 'ncbi')
names(common.names) <- species.found

common.names

## output:
## $`Helianthus annuus`
## [1] "common sunflower"
##
## $`Mycobacterium bovis`
## character(0)
##
## $`Rattus rattus`
## [1] "black rat"
##
## $`Mus musculus`
## [1] "house mouse"

This works correctly if no common name is available, or if there is no species in the database. You can also try a different database, e.g. db='itis'.
The advantage of this approach is that you do not have to parse a file yourself. The downside is that it can be slow for large list, since for every species a database request is performed.

great stuff, thanks folks- this works perfectly for the ncbi data. one small addition is that due to the complexities of the eol database being searchable in a different way you can't use that option without some other wizzy functions (i have a chap working on this but no joy so far)

use the NCBI taxon dump under ftp://ftp.ncbi.nih.gov/pub/taxonomy
in taxdmp.zip you'll find all the names for a given NCBI taxon
$ grep -w ^9606 names.dmp
9606    |   Homo sapiens    |       |   scientific name |
9606    |   Homo sapiens Linnaeus, 1758 |       |   authority   |
9606    |   human   |       |   genbank common name |
9606    |   man |       |   common name |

Getting back to the original question about using taxize to search specifically in EOL, it looks like the most recent version of the package does not produce this error.
The current version of R I'm using (3.4.2) produces the following output:

sci2comm( c("Cancer", "Lophopanopeus"), db="eol" )
$Cancer
   [1] "En art taskekrabbe"              "Felsenkrabbe"                    "Atlantic rock crab"
   [4] "Jaiba de roca amarilla"          "Tourteau poïnclos"               "Granciporro atlantico giallo"
   [7] "Atlantische rotskrab"            "Sapateira de rocha do Atlântico" "Bogenkrabbe"
  [10] "red rock crab"                   "Pacific rock crab"               "Jaiba del Pacífico"
  [13] "Tourteau du Pacifique"           "Granciporro del Pacifico"        "Pacifische rotskrab"
  [16] "Sapateira de rocha do Pacífico"  "Jonahkrabbe"                     "Jonah crab"
  [19] "Jaiba de roca Jonás"             "Tourteau jona"                   "Grancipporo atlantico rosso"
  [22] "Jonakrab"                        "Sapateira boreal"                "Bou"
  [25] "Almindelig taskekrabbe"          "Taskekrabbe"                     "Taschenkrebs"
  [28] "Edible crab"                     "rock crab"                       "Buey"
  [31] "Buey de mar"                     "Jaiba de roca masera"            "Masera"
  [34] "Paguro"                          "Pato"                            "Sabago"
  [37] "Isotaskurapu"                    "Crabe dormeur"                   "Dormeur"
  [40] "Tourteau"                        "crabe-dormeur"                   "Granchio di mare"
  [43] "Granciporro"                     "Noordzeekrab"                    "grote krab"
  [46] "grote noordzeekrab"              "hoofdkrab"                       "slagkrab"
  [49] "steenkrab"                       "toogkrab"                        "zeekrab"
  [52] "Sapateira"
$Lophopanopeus
  [1] "blackclaw crestleg crab" "knobkneed crestleg crab" "knobknee crestleg crab"  "molarless crestleg crab"
$'Pugettia foliata'
  character(0)

A slightly older version of R (3.3.2) produces an error when one of the taxa was not found:

sci2comm( c("Cancer", "Lophopanopeus", "Pugettia foliata"), db="eol" )
Error in eol_pages(taxonconceptID = x, common_names = TRUE, ...) :
    Not Found (HTTP 404).

I am very thankful for the package updates!

As Pierre suggested, you can get a dump of such names from NCBI. Then, to query it for the common name for a species using its Latin name, you can do:
$ awk -F'\t' -vname="Mus musculus" '($7=="scientific name" && $3==name){
                                      a[$1]=$3
                                    }
                                    ($1 in a && /genbank common name/){
                                      print $3
                                    }' names.dmp
house mouse

please note : the tabulation is not the delimiter, it's a mix of spaces and pipes

@Pierre is it? I admit I've never worked with this file before but I downloaded it and saw that it is tab-separated. The | are their own fields and are actually \t|\t in the file, so using tab as the separator seems the best approach. Am I wrong?

By the way, @Pierre, this is only making your answer a bit more specific. It might be better to just add it into your answer instead of cluttering up the page with multiple answers. I would have done so but I don't have enough rep to do so directly yet. So feel free to just copy this (or your improved version of it) into your answer. If you do, let me know and I'll delete this one.

ah yes, with your awk indexes, it should be ok. And no, leave your answer here please.

Here is a code on github that can generate a csv file as well as a python dictionary.
https://github.com/rahulnutron/scientific_name_to_common-name_converter


text:
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

Welcome to Bioinformatics! Could you please clarify what is a geneind object? And could you provide the exact output of the table? (If you provide how you get that table that would probably be also helpful to answer you)

I tried to insert the table, though I can't quite format it right... sorry, I'm still figuring this out!

By Fst are you refering to Fixation index?

Please also indicate the non-standard packages you are using. I’m assuming amova comes from ‹pegas›?

Yes @Llopis I mean fixation index, and yes Konard Rudolph, the amova I am using is from the package pegas.
Sorry for not specifying!

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.

That's what I thought. Thank you for confirming it. And thanks for the explanation for error!

@BeatriceBaldi if this answers your question, you can mark as solved by checking the tick at the left of the answer. This way you let others know that this problem is solved


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

The question you should be asking isn't "what can we do with this?", but rather "what biological question do we want to answer?" You sound like you're aimlessly fishing for anything, which is likely not worth your time.

What do you mean by "process"? Process the list in order to do what? What question are you asking? Please [edit] your post and clarify what you are trying to do. As it stands, a valid answer would be that you can use wc to count the number of lines in the file. That is a type of "processing" but, obviously, not very helpful for you. The more specific you make your question, the higher the chances of your getting a useful answer.

Hi Peter, thanks for your question and welcome to Bioinformatics Stack Exchange. Bioinformatics is a very large subject area, frequently with many different solutions to the same problem. It would be helpful if you could give more story / context surrounding your problem to guide answers. For example, what organism are you working with? Are these conditions treatment vs control, or separate treatments (without a control)? What kind of results are you interested in getting out of this analysis?

Thank you all for the feedback. I understand that each dataset has to be studied according to the questions of interest, just wanted to get a few ideas. Also tried to make the question general so that it will be useful for more people, but I now added the specifics.

I think the question is still too broad. What knowledge do you expect to learn? Could you narrow more to one of your interests (Process, pathways, or relations between cells lines)?  Which is your "immunological point of view" ?

It's best to stick to questions that are as specific as possible, and let other people decide on their own generalisations about the answers. There are so many different things that can be done in bioinformatics even for a single task that it's no useful to speculate about what other similar things could be done.

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.

We should move to a chat,... but study 2 looks like a study I did: how did you make your enrichment analysis (outside IPA), specially of GO? How did you try WGCNA? Did you do a GSEA of the genes up-regulated in the comparison of a treatment vs control in the comparison of treatment1 vs treatment2?

WCGNA: We had about 80 different biological replicates comprising about 25 different treatment conditions, so there was enough expression variation to be able to pull out a few shared gene sets

GO enrichment analysis wasn't done; I don't know of any GOEA tools that use the hierarchical structure of the ontology map, and without that it's very prone to selection bias. For GSEA, we were mostly looking at treatment vs control (for various treatments and cell populations), and I think using the gene sets in the database that seemed to be immune-related. I'm not completely sure on that; it was other scientists that did the GSEA.

For the WGCNA, I first did a VST transform of the expression data in DESeq, normalised by the length of the longest transcript for that gene, then determined the 5,000 genes that were most variable (by MAD) across all our samples. This matrix of "lowFiltered" genes was fed into WGCNA by another scientist, and we spent a fair amount of time trying to understand the process and get the results looking right.

topGO package of Bioconductor uses the DAG structure of GO to find the enriched terms. WGCNA requires homogeneous samples, so if you mixed the 25 treatment conditions in those 80 samples to build a network, it would not detect the network for each treatment. You would need to do a consensus analysis (with multiset.* functions) to retrieve the network shared between each treatment, but WGCNA requires above 15 samples per condition to effectively work.

Thanks, I might give topGO a try once I've sorted out the proportionality calculations


text:
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?

CRAM doesn't need to be lossy, is there a reason you need it to be?

Saving disk space. We pay by the GB and need to keep the data around for 10 years.

Can't argue that budget isn't a good reason :)

Interesting question. I think this is the kind of thing that makes a nice side project. Take a bam file call the variants, transform it to cram and run the variant caller. Measure the difference and the variant concordance between the two approaches using a number of different files.

If one is concerned about absolute integrity/reproducibility of the data then budget consideration is not a good reason.

Unless the answer is: some types of lossy compression don't have any impact on variant calls.

BAM files are zipped with standard gzip compression. Unzip them to "naked BAM" - not my terminology - and re-zip them with something stronger like 7zip/LZMA. You can always re-zip them again with the bgzip tool when you need them back in true BAM format again. This gets you most of the way there filesize-wise without really changing the format, which could be good if you've got things set up how you like. Not an answer because it doesn't answer your question but it might solve your problem.

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.

These are all great suggestions for reducing file size without losing information, but don't address the main question: the effect of lossiness on variant calls.

@DanielS If you don't touch bases and qualities and name pairing, you won't change the variant calls.

Yes, but then that's not really lossy, is it? Doesn't lossy compression conventionally involve changing the sequence and/or quality values for greater compression efficiency?

That depends on the definition of "lossy" :) To me, losing read names and tags is lossy.

¯\(ツ)/¯ You yourself said that calls shouldn't change if the sequence and quality are unchanged. So everything else is ancillary. Don't get me wrong, I think it's valuable to point out that it's possible to reduce file size without changing sequence or quality, but it seemed pretty clear to me that the OP was talking about lossy compression of the sequence and/or quality scores.

Then again, this answer is a reasonable response to the question "Which compression strategy has the lowest impact?" Ok, I take it all back! :-)

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.

+1. My colleagues have done some benchmarks to show 4-bin has little effect on germline samples. I have seen the similar. Cancer samples always make me wary, though. It would be great if someone do a systematic evaluation on cancer samples. I haven't seen one so far.


text:
Is there a software or web server to calculate thermodynamic parameter enthalpy change of cellulase?
I want to calculate thermodynamic paremeter enthalpy change of cellulase, but I don't know what software can do this work. Please tell me, thank you very much !

Hi and welcome to the site! Please [edit] your question and clarify what exactly you mean by "thermodynamic paremeter". A thermodynamic parameter for what? And what parameter? Also, you have tagged this with 5 tags, are they really relevant to your question? What is the connection with structural variation for instance?

Enthalpy, how to calculate this by a software or a web server ? Thank you!

Again, the enthalpy of what? A folded RNA structure? The enthalpy of a chemical reaction? Please [edit] your question and clarify what you need.

This is the first question I asked, thank you for your advice.

This looks more like computational chemistry. Not sure if we have experts here to answer this question. 2) Perhaps it will help if you clarify what kind of change; e.g. you can find enthalpy of formation and enthalpy of  combustion in wikipedia

Hi user239932, this looks like a pure biochemistry question, and would probably be better asked in SE Biology or SE Chemistry, depending on what type of enthalpy you're looking for.

I'm against closing this question as off topic. Yes, it's borderline, but it involves computation and enzyme structure, and I'd hate to give people working at the interface of biochemistry/biophysics and structural protein/RNA informatics the impression that their questions are not welcome here. That said, the question is unclear. If the question is closed, this should be the reason.

I think question is fine, but @user239932 should edit it to make it answerable. Please, specify what exactly you want to compute. If you do not know, try to at least write why you want to do that. You can also skip "Please tell me, thank you very much !", we all believe that you are a nice person, kind behaviour is implicit here.

You could use PISA assuming you have a PDB file of cellulase? Sorry it's a short answer and should be a comment, I don't have the reputation, however.

Thank you for your answer, I will try it. What about YASARA＋FoldX ?