text:
What motif finding software is available for multiple sequences ~10Kb?
I have around ~3,000 short sequences of approximately ~10Kb long. What are the best ways to find the motifs among all of these sequences? Is there a certain software/method recommended?
There are several ways to do this. My goal would be to:
(1) Check for motifs repeated within individual sequences
(2) Check for motifs shared among all sequences
(3) Check for the presence of "expected" or known motifs
With respect to #3, I'm also curious if I find e.g. trinucleotide sequences, how does one check the context around these regions?
Thank you for the recommendations/help!

Are you looking for motifs shared by all the sequences? Motifs repeated in each sequence? Are you looking for the presence of specific, known motifs?

"Are you looking for motifs shared by all the sequences? Motifs repeated in each sequence?" I was looking for both motifs shared by all sequences and repeated in sequences, but actually it would be interesting to check each sequence for a repeated sequence. "Are you looking for the presence of specific, known motifs?" No, but it would be interesting to look into this in retrospect

OK, please [edit] your question and ass this information. Comments are easy to miss, hard to read and can be deleted without warning. The three problems (de-novo shared motif identification, de-novo repeated motif identification and detecting known motifs) are different and each requires its own approach.

@terdon Of course

Do you mean 10kb in total, or 3000 sequences, each of length around 10kb?

@CharlesE.Grant 3000 sequences, each of length around 10kb

For (3), this page has a lot of links to pattern/motif finding tools. Following through the YMF link on that page, I came across the University of Washington Motif Discovery section. Of these projection seemed to be the only downloadable tool. I find it interesting how old all these tools are; maybe the introduction of microarrays and NGS has made them all redundant.
Your sub-problem (2) seems similar to the problem I'm having with Nippostrongylus brasiliensis genome sequences, where I'd like to find regions of very high homology (length 500bp to 20kb or more, 95-99% similar) that are repeated throughout the genome. These sequences are killing the assembly.
The main way I can find these regions is by looking at a coverage plot of long nanopore reads mapped to the assembled genome (using GraphMap or BWA). Any regions with substantially higher than median coverage are likely to be shared repeats.
I've played around in the past with chopping up the reads to smaller sizes, which works better for hitting smaller repeated regions that are such a small proportion of most reads that they are never mapped to all the repeated locations. I wrote my own script a while back to chop up reads (for a different purpose), which produces a FASTA/FASTQ file where all reads are exactly the same length. For some unknown reason I took the time to document that script "properly" using POD, so here's a short summary:

Converts all sequences in the input FASTA file to the same length.
     Sequences shorter than the target length are dropped, and sequences longer
     than the target length are split into overlapping subsequences covering
     the entire range. This prepares the sequences for use in an
     overlap-consensus assembler requiring constant-length sequences (such as
     edena).

And here's the syntax:
$ ./normalise_seqlengths.pl -h
Usage:
    ./normalise_seqlengths.pl <reads.fa> [options]

  Options:
    -help
      Only display this help message

    -fraglength
      Target fragment length (in base-pairs, default 2000)

    -overlap
      Minimum overlap length (in base-pairs, default 200)

    -short
      Keep short sequences (shorter than fraglength)

The MEME Suite web site contains a collection of tools for motif analysis (I'm one of the maintainers). It contains two de novo motif discovery tools: MEME and DREME. Public web applications are provided, but you can also download and build command line tools for a local installation.
For your first goal you could use MEME and select the "Any number of repetitions model" (ANR). For your second goal, you'd use MEME with the "Zero or One Occurrences Per Sequence" (ZOOPS) model. For your third goal you could use FIMO (Find Individual Motif Occurrences), and one or more of the motif databases provided on the software and database download page.
It sounds like your sequence data is about 30Mb. The MEME web application is limited to 60kb of sequence data, so you'd have to install a local copy of the MEME Suite. MEME would take a long time to analyze a 30Mb sequence database unless you have MPI configured, and lots of cores available. You might want to consider analyzing a randomly selected subset of your sequences. The running time of MEME grows as the cube of the number of sequences.
For short motifs, you may want to use DREME rather than MEME. DREME is better than MEME at identifying short motifs, but is limited to motifs <= 8 positions wide.

Check out HOMER. "Software for motif discovery and next generation sequencing analysis", it's what my lab uses currently for finding eRNA motifs.
Edit: For @ShanZhengYang "HOMER was designed as a de novo motif discovery algorithm..."
HOMER De Novo Motif

Can one do de novo motif discovery? If not, I'm not entirely sure how I would create a background with this model.

@ShanZhengYang See my edit.

Thanks. I'm still not sure given my set-up how to choose a background for Homer however...

I've never used it for De Novo, so take this with a grain of salt, but try the defaults first. Perhaps shoot the developers an email.

To expand on @EMiller 's comment, for many model organisms, HOMER can figure out a background model for you. You just have to install the supporting files for your model organism using  configureHomer.pl

Most tools I know of looks for enrichment of specific motifs - but that requires that you have a set of sequences which are of special interest and a background set to test against.
Is that your case?
Update after comments 12th June 2017.
You could try the meme suite more specifically the motif finder

Hi Kristoffer Vitting-Seerup, thanks for your comment and welcome to Bioinformatics Stack Exchange. We appreciate any comments like this that help people ask better questions. I have flagged your response as "not an answer" because I think it should be added as a comment to the question (which you can't do yourself until you reach a reputation of 50).

Hi there, yes this would be a better comment to the original question. This is the SO style. Regarding your question, I don't have a background set

it is under development, but maybe BaMMmotif! is something for you? Its main selling point is that it can look for motifs enriched in a set of sequences of equal length de novo. If you can't/don't want to supply a negative set it learns one from the positive sequences. There is a wealth of options to choose from if you have more information about your sequences: there are different models for "zero or one", "one" and "multiple" occurrences of the motif.
You can also use it to look for known motifs, if you encode them as an XXmotif PWM. If you have a file with motifs (like binding sites) you can use that as initialization as well.
While I have not used the software myself, the authors are very responsive on git and the installation instructions seem pretty straightforward.
EDIT: Apparently the software was developed with ChIP experiments as the usual use case scenario, it might misbehave for larger sequences or take super long to run.

Hi galicae, thanks for your answer and welcome to Bioinformatics Stack Exchange. Your answer is a bit short, and could do with a bit more explanation. Could you please add a bit more information into your answer about BaMMmotif? Is it a tool that you have used before? Would this be an appropriate program to find novel shared motifs, repeated motifs within the same sequence, or known motifs?

@gringer fair enough, edited.


text:
How can longest isoforms (per gene) be extracted from a FASTA file?
Is there a convenient way to extract the longest isoforms from a transcriptome fasta file? I had found some scripts on biostars but none are functional and I'm having difficulty getting them to work.
I'm aware that the longest isoforms aren't necessarily 'the best' but it will suit my purposes.
The fasta was generated via Augustus. Here is what the fasta file looks like currently (sequence shortened to save space)
>Doug_NoIndex_L005_R1_001_contig_2.g7.t1
atggggcataacatagagactggtgaacgtgctgaaattctacttcaaagtctacctgattcgtatgatcaactcatca
ttaatataaccaaaaacctagaaattctagccttcgatgatgttgcagctgcggttcttgaagaagaaagtcggcgcaagaacaaagaagatagaccg
>Doug_NoIndex_L005_R1_001_contig_2.g7.t2
atggggcataacatagagactggtgaacgtgctgaaattctacttcaaagtctacctgattcgtatgatcaactcatca

The format is as such:
Gene 1 isoform 1
Gene 1 isoform 2
Gene 2 isoform 1
Gene 2 isoform 2

and so forth. There are several genes that have more than one pair of isoforms (up to 3 or 4). There are roughly 80,000 total transcripts, probably 25,000 genes. I would like to extract the single longest isoform for each gene.

It depends on the exact format of your fasta. Where did you get it? What do the fasta header lines look like? Can you paste a small example?

Doug_NoIndex_L005_R1_001_contig_4.g13.t1

(sequence data here)

Doug_NoIndex_L005_R1_001_contig_4.g13.t2

(sequence data here)

What exactly have you tried? 2. What does them not being functional mean? 3. Do you just want the single longest one or the top N longest?

I've tried all of the solutions here: https://www.biostars.org/p/107759/

Non functional in that I couldn't get it work

I would like the single longest read for each set of isoforms

@ZincFingers: Please put all of this in your original post and include whether isoforms are ever separated from each other by unrelated genes.

You don't have the same sequence IDs there. Should we assume that the part after the last . in the ID line needs to be discarded?

I'm not sure what you mean. The ".t1" or ".t2" indicates which isoform the read is from, I would like to discard the shorterst isoforms and retain the longest isoform.

OK, so they do need to be removed.

While the solution from https://bioinformatics.stackexchange.com/users/96/daniel-standage should work (after adjusting for possible python3 incompatibility), the following is a shorter and less memory demanding method that uses biopython:
#!/usr/bin/env python
from Bio import SeqIO
import sys

lastGene = None
longest = (None, None)
for rec in SeqIO.parse(sys.argv[1], "fasta"):
    gene = ".".join(rec.id.split(".")[:-1])
    l = len(rec)
    if lastGene is not None:
        if gene == lastGene:
            if longest[0] < l:
                longest = (l, rec)
        else:
            lastGene = gene
            SeqIO.write(longest[1], sys.stdout, "fasta")
            longest = (l, rec)
    else:
        lastGene = gene
        longest = (l, rec)
SeqIO.write(longest[1], sys.stdout, "fasta")

If you saved this as filter.py, then filter.py original.fa > subset.fa would be the command to use.

Here is a solution in R. Could get really slow with big files. Works for the example you posted.
library(Biostrings)

## read your fasta in as Biostrings object
fasta.s <- readDNAStringSet("sample.fa")

## get the read names (in your case it has the isoform info)
names.fasta <- names(fasta.s)

## extract only the relevant gene and isoform id (split name by the period symbol)
gene.iso <- sapply(names.fasta,function(j) cbind(unlist(strsplit(j,'\\.'))[2:3]))

## convert to good data.frame = transpose result from previous step and add relevant column names
gene.iso.df <- data.frame(t(gene.iso))
colnames(gene.iso.df) <- c('gene','isoform')

## and length of isoforms
gene.iso.df$width <- width(fasta.s)

## split data.frame into list with entry for each gene
gene.iso.df.split <- split(gene.iso.df,gene.iso.df$gene)

## optional to keep all the information but really you just need indices
##gene.iso.df.split.best <- lapply(gene.iso.df.split,function(x) x[order(x$width)[1],])

## pull out the longest isoform ID for each gene (in case of a tie just take the first one)
best.id <- sapply(gene.iso.df.split,function(x) row.names(x)[order(x$width)[1]])

## subset your original reads with the subset
fasta.s.best <- fasta.s[best.id]

## export new fastafile containing longest isoform per gene
writeXStringSet(fasta.s, filepath='sample_best_isoform.fasta')

Late to the party here, but I like to try to avoid writing scripts when some command line magic will do. It's good practice to index your FASTA so use it.
samtools faidx <myfasta.fa>

A little bit of awk and sort can determine the largest isoform of each gene (and they don't even have to be sorted by name in the source FASTA).
awk -F'[\t.]' '{print $1,$2,$3,$4}' <myfasta.fa>.fai | sort -k4nr,4 | sort -uk1,2 | cut -f1-3 -d' '| tr ' ' '.' > selection.ls

The -F[\t.] specifies both the fai tab delimiters, and the full stops in your contig names (Doug_NoIndex_L005_R1_001_contig_2.g7.t1) as delimiters for awk.
For each line we print four columns. $1 to $3 are the three components of your contig name (after splitting on the .), and column $4 is the column in the index for sequence length.
We use sort to organise the lines by $4 (the contig length). Reverse sort with -r so the longest genes appear first.
We then subsort by columns $1 and $2 (the contig and gene name (g1...)). -u then drops any repeat of column $1 and $2 pairs (that is, all lines after the first pair -- the longest length).
Finally, we reassemble $1 to $3 with cut and tr to give you a list of sequence names.

A loop with our old friend samtools faidx can put your chosen sequences out to a new FASTA.
while read contig;
    do samtools faidx <myfasta.fa> $contig >> selection.fa;
done < selection.ls

Don't forget to index it, they are useful!
samtools faidx selection.fa

Oh and if you wanted to remove the .tN part, just switch out that cut -f1-3 -d' ' for cut -f1-2 -d' '

An ex-coworker (Josep Avril) has written a couple of very useful little scripts that convert fasta to tbl (seqID<TAB>Sequence) and back again. These are extremely handy for this sort of thing (and are included at the end of this answer). Using them, you can convert your fasta to a one sequence per line format, keep the longest sequence with a simple awk script and convert back to fasta.
I am assuming that the part after the last . in the ID line should be removed (so that Doug_NoIndex_L005_R1_001_contig_2.g7.t1 and Doug_NoIndex_L005_R1_001_contig_2.g7.t2 both map to Doug_NoIndex_L005_R1_001_contig_2.g7). If that is a correct assumption, this should work for you:
$ FastaToTbl file.fa | sed 's/\.[^.]*\t/\t/' |
    awk -F"\t" '{
                    if(length($2)>length(a[$1])){
                        a[$1]=$0
                    }
                }
                END{
                    for(i in a){
                        print a[i]
                    }
                }' | TblToFasta
>Doug_NoIndex_L005_R1_001_contig_2.g7
atggggcataacatagagactggtgaacgtgctgaaattctacttcaaagtctacctgat
tcgtatgatcaactcatcattaatataaccaaaaacctagaaattctagccttcgatgat
gttgcagctgcggttcttgaagaagaaagtcggcgcaagaacaaagaagatagaccg

One possible issue with that approach is that it needs to keep one sequence per ID in memory. If you have a huge fasta file and not much memory, that might be a problem. If that's the case, you can first sort the file to ensure that similar IDs always appear together and then print each line as soon as you find the next ID:
FastaToTbl file.fa | sed 's/\.[^.]*\t/\t/' | LC_ALL=C sort -t$'\t' -k1,1 |
    awk -F"\t" '{
                    if(NR>1 && prev!=$1){
                        print a[$1];
                        prev=$1
                    }
                    if(length($2)>length(a[$1])){
                        a[$1]=$0
                    }
                }
                END{
                    print a[$1]
                }' | TblToFasta

FastaToTbl
#!/usr/bin/awk -f
{
        if (substr($1,1,1)==">")
        if (NR>1)
                    printf "\n%s\t", substr($0,2,length($0)-1)
        else
            printf "%s\t", substr($0,2,length($0)-1)
        else
                printf "%s", $0
}END{printf "\n"}

TblToFasta
#! /usr/bin/awk -f
{
  sequence=$NF

  ls = length(sequence)
  is = 1
  fld  = 1
  while (fld < NF)
  {
     if (fld == 1){printf ">"}
     printf "%s " , $fld

     if (fld == NF-1)
      {
        printf "\n"
      }
      fld = fld+1
  }
  while (is <= ls)
  {
    printf "%s\n", substr(sequence,is,60)
    is=is+60
  }
}

This solution is will work for your example and probably all Augustus-derived Fastas, but mileage will vary beyond that.
#!/usr/bin/env python
from __future__ import print_function
import sys

def parse_fasta(data):
    """Stolen shamelessly from http://stackoverflow.com/a/7655072/459780."""
    name, seq = None, []
    for line in data:
        line = line.rstrip()
        if line.startswith('>'):
            if name:
                yield (name, ''.join(seq))
            name, seq = line, []
        else:
            seq.append(line)
    if name:
        yield (name, ''.join(seq))

isoforms = dict()
for defline, sequence in parse_fasta(sys.stdin):
    geneid = '.'.join(defline[1:].split('.')[:-1])
    if geneid in isoforms:
        otherdefline, othersequence = isoforms[geneid]
        if len(sequence) > len(othersequence):
            isoforms[geneid] = (defline, sequence)
    else:
        isoforms[geneid] = (defline, sequence)

for defline, sequence in isoforms.values():
    print(defline, sequence, sep='\n')

Save as longest.py, and invoke on the command line like so.
python longest.py < intput.fasta > output.fasta

The line geneid = '.'.join(defline[1:].split('.')[:-1]) is key. Let me break it down.

defline[1:]: ignore the first character of the defline (the > symbol)
defline[1:].split('.'): split the string on . symbol
defline[1:].split('.')[:-1]: ignore the last value after the split (the isoform name)
'.'.join(defline[1:].split('.')[:-1]): join the split values again

This value will be the same for all isoforms for the same gene, so we use it in the script to keep track of the longest sequence associated with each gene ID.

This could save a lot of memory and a bit of time if you write out and reset isoforms whenever geneid not in isoforms. OP says isoforms are always next to each other.

Hi Daniel, I tried running this script and received the following error message:
File "longest.py", line 29
     print(defline, sequence, sep='\n')
                            ^
SyntaxError: invalid syntax

@ZincFingers: Change it to print(defline) and print(sequence) on separate lines.

Python 2 / 3 incompatibility. Will update.

source: Stack Exchange

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

sort: votes


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

Have you considered emailing 10xGenomics and asking? You're already paying them money for their products, so they'll reply.

I know how to make the reference genome using cellranger, in which case they could help. My question is more on how the GRCh## versioning works, for which I think here is a better and definitely faster place to ask rather than 10xGenomics.

The different Ensembl releases are associated to genome patches with different patch contigs, so it becomes a question whether their software can handle missing contigs.

Okay, right! This is a question I should ask them. Thanks @DevonRyan . However, if their answer is positive, would it still be significantly better to use both parts of the same release, or it doesn't make that much of a difference?

If it doesn't break their software then it shouldn't matter either way :) In reality, it's trivial to strip a whole chromosome out of a GTF, so just use GTF/fasta from the most recent Ensembl release and you'll get slightly more reliable results.

BTW, feel encouraged to post their reply, since then it'll be indexed on google.

You can see the exact assembly by looking in the README file in Ensembl's genome download page. ftp://ftp.ensembl.org/pub/release-89/fasta/homo_sapiens/dna/README
As you can see the current assembly is GCA_000001405.25. Ensembl 79 uses version GCA_000001405.17. You can expect some differences between assembly versions.
If you take a look at https://www.ncbi.nlm.nih.gov/assembly?term=GRCh38&cmd=DetailsSearch
You can see that 20,466,394 bp have been added since 17 -> 25. This is only about 0.6% change so I would guess there is not that much difference between these versions. However, if I were you I would use the corresponding version of Ensembl. Even better would be using the latest version of Ensembl, as not only the genomes primary sequence but also the positions of genes can change between versions (https://doi.org/10.1093/bib/bbw017).

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


text:
What unit do I get on the y-axis of a metagene profile plot?
I start with a sorted and indexed bam file ("mapped.bam") representing the mapping of small reads on a reference genome, and a bed file ("genes.bed") containing the coordinates of a set of features of interest (let's say they are genes), for which I want to compute an average profile using programs from deeptools. I would like to understand the steps involved to be sure of what the vertical axis of the final profile represents.
First step: making a bigwig file
I create a bigwig file ("mapped.bw") from the bam file using bamCoverage as follows:
bamCoverage -b mapped.bam -bs 10 -of=bigwig -o mapped.bw

The help of bamCoverage says:

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

In my case, the bins are 10 bp long. My reads are longer than that.
For a given bin, a given read can:

completely overlap the bin
overlap the bin on n bp, n < 10
not overlap the bin at all

Please correct me if I'm wrong: My guess it that the read is counted as 1 in cases 1. and 2., and 0 otherwise, and I also suppose that a read can be counted for several successive bins if it is long enough.
Second step: averaging over genes and plotting
I compute a "meta profile matrix" ("mapped_on_genes.gz") using computeMatrix scale-regions as follows:
computeMatrix scale-regions \
    -S mapped.bw \
    -R genes.bed \
    --upstream 300 \
    --unscaled5prime 500 \
    --regionBodyLength 2000 \
    --unscaled3prime 500 \
    --downstream 300 \
    -out mapped_on_genes.gz

(There is a -bs parameter which default value is 10 according to the help of the command.)
I use this to plot a profile using plotProfile:
plotProfile -m mapped_on_genes.gz \
    -out mapped_on_genes_meta_profile.pdf

I obtain a profile in with values on the y axis. In what units are these values?
My guess is the following:
For the upstream (300 bp) and internal 5-prime (500 bp), since the bin size was the same in bamCoverage and computeMatrix, each point on the x axis probably represents a 10 bp window, and its y coordinate is the average over the regions present in the bed file of the corresponding bins in the bigwig file, so it is an average number of reads overlapping a 10 bp bin.
Same thing at the 3-prime and downstream side.
For the central 100 bp portion, before averaging over regions some shrinking or spreading of the bins must have been performed, I guess by averaging between neighbouring bins. So the final unit is still a number of reads overlapping a 10 bp bin.
And if I use larger bins, I should end up with proportionally higher values.
Am I correct?

Feel free to @ me in deepTools questions, since I'm the primary developer.
For a given bin, the count assigned to it is the number of reads overlapping it, regardless of whether they overlap by 1 or 10 bases. So a read overlapping only partially and one overlapping completely are treated the same.
Since your bigWig file is in units of "alignments" (i.e., it's not 1x normalized), the resulting profile will also be in units of "alignments" (i.e., profiles and heatmaps are in whatever units the input files are in).
Upstream/downstream regions and unscaled regions are also 10 base bins. Note that these are then the average of the per-base value, since the bins here may not perfectly correspond to the bins in the bigWig files. The line in the profile plot is indeed the average (by default, you can choose median, max, min, etc.) of the underlying regions for each bin.
Regarding the scaled section in the middle, the the number of genomic bases per bin is changed such that the region will have "length"/(regionBodyLength/binSize) bases each. As above, the per-base value is then averaged (or whatever you specify) to derive the per-bin value. The length here is decreased if you have unscaled regions, since otherwise bases would get counted twice.

When you write "the average of the per-base value" you mean that if a region 10 bp bin overlaps with 3 bp from genomic bin i (having m reads) and 7 bp from genomic bin i+1 (having n reads), then the value for this region bin will be (3*m + 7*n) / 10?

Correct, the average will be weighted by the amount of overlap exactly as you showed.


text:
How do I compare measured ratios to calculated ratios in peptide mass spectrometry?
I have two different sets of internal standards in a sample (1), and a sample (2) with these internal standards as well as about 30 peptides. Each peptide has one of each of these kind of standards. From sample 1, I can get the true ratio between standard 1 (S1) and standard 2 (S2), and from sample 2 I can get the ratio between the analyte (L) to each of the standards:
Sample 1: S1/S2
Sample 2: L/S1, L/S2
I can thus calculate the ratio S1/S2 from the measured ratios in sample 2. Each sample is run in 5 technical replicates (that is; 5 injections into the LC-mass spectrometer), so S1/S2 ratios for one particular analyte could look like this (mock data):
measured    calculated
0.967       0.987
1.007       0.967
1.044       1.012
1.041       1.025
1.048       1.046

I want to compare the measured ratios to the calculated, to see if there are any bias in my acquisition of data. As the calculated and the measured ratio are from different samples, the replicates aren't "linked" so I could not use for example Bland Altman analysis. Should I use some kind of ANOVA? Or do you have any other recommendation?
I would also like to make an overall assessment of the ratios of all 30 peptides to see if the bias is different for different peptides. How could I do that?
Any input would be greatly appreciated.

Hi Max Jonatan Karlsson, thanks for your question and welcome to Bioinformatics Stack Exchange. My first glance at this question suggests it might be a pure statistics question, and therefore would be more appropriate on the stats Stack Exchange, Cross Validated. If you can give more context/story as to why this question would be be helped by computer science and/or biological knowledge (e.g. you're calculating thousands of these statistics), feel free to update your question.

The reason why I posted it here is because it is proteomics data, but I stripped down the context as I deemed it not to be important for the question. I see that you are right that this belongs in statistics however, so I will post my question there, as well as update the question here with some context. Thanks!

Context is always helpful for bioinformatics questions. It's frequently the case that there are multiple tools that serve a similar purpose, and the best tool (or tools) can change depending on the nature of the problem. Context also helps to identify when the core problem doesn't match the question (i.e. the XY problem); that doesn't seem to be the case here, but it's worth bearing in mind when asking questions in the future.

From a statistical perspective ratios are icky because the variance of the sampling distribution of a ratio decreases with the size of the numerator and denominator being sampled to create the ratio. If you can get at the numerator and denominator separately, then you should absolutely do that and look at count based tests like the Chi-square test, Fisher's exact test, or Poisson regression.
That said, a t-test would be appropriate in the two-sample case if its assumptions are met.  You need to check that both samples are normally distributed, ideally via a qq-plot. You also need to check the variance of both samples and choose a variant of the t-test that does or does not assume equal variances.
If the samples are not normally distributed, then the next approach would be a Mann-Whitney U test.


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

There are some helpful hints here!

DESeq2 is able to adjust p-values for single-replicate samples by estimating shared dispersion across all samples. DESeq2 will give a warning, but try its best to carry out the analysis. More information on doing that can be found in the DESeq2 vignette and the DESeq2 manual page.
The results will not be anywhere near the sensitivity that can be achieved when replicate information is available, but it might be good enough if only very obvious large-scale changes are of importance.

Agreed. DESeq2 uses an empirical Bayes model that incorporates knowledge about the variation of all genes to estimate the true expression value of any single gene in any single sample. Unfortunately due to its complexity I find it incredibly hard to figure out what's going on under the hood when I use it.

I just accept that I can never have the depth of knowledge of analysis tools that the tool authors have, and treat data analysis as similar to a biological experiment -- you don't have to know everything about how things work in order to gain understanding from them.


text:
What are the de facto required fields in a SAM/BAM read group?
The SAM specification indicates that each read group must have a unique ID field, but does not mark any other field as required.
I have also discovered that htsjdk throws exceptions if the sample (SM) field is empty, though there is no indication in the specification that this is required.
Are there other read group fields that I should expect to be required by common tools?

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

The mystery was explained in a samtools-devel message at the time: Sometimes, we simply want to group reads without sample information. See also hts-specs issue #286.

The read group identifier needs to be specified in both the header lines of the BAM/SAM file and the alignment line. No other fields are required, but note that because the additional information is only stored once (i.e. in the header), it won't add much to file sizes or processing time if additional fields are included. If any of the other optional fields cause problems (either by inclusion or exclusion), it would be helpful to report an issue about it.
Older tools required you to specify both a header read group ID and an alignment read group ID, but most tools of the recent tools I've used seem to be intelligent enough that they will add headers for any read groups without that information.


text:
Spliced vs. unspliced ratios for transcripts in RNA-seq data
Is there a computational tool for measuring what percentage of RNA is spliced in an RNAseq experiment?
I'm not particularly interested in complicated analyses that give ratios for all possible alternative splicing variations. I'd rather have a binary classification. For example: (RNA1 = 35% spliced, 65% unspliced)
Edit: (Additional information)
We're very early in the exploratory stage, so a specific hypothesis is still forming. But we're interested in the sequence content of lncRNAs, many of which are poorly spliced. So we want to get a sense for how prevalent unspliced version of these transcript are.
Because this is currently exploratory, estimations and noise are fairly acceptable. A very basic approach might be to calculate the ratio of reads that map to exon-to-exon junctions, verses those that map from exon-to-intro junctions.
A more advanced approach I'm also considering is Salmon.

Did you ever settle on something? I have the same question as you.

We ended up using aligning with STAR then quantifying with salmon. We gave salmon all known splice variants, including the unspliced annotation. Then we summed all spliced isoforms to get "spliced vs. unspliced".
This pipeline was easy, accurate, and we were really satisfied with the result.

I don't know about tools, but i've used the following python code to calculate the ratio of reads that overlap the 5' or 3' ends of introns or that are spliced. We sum these across all introns in a gene set (we actaully use this for iCLIP analysis to see if RNA binding proteins bind pre-mRNA or spliced RNA).
import pysam
from collections import Counter
from CGAT import GTF, IOTools

def calculateSplicingIndex(bamfile, gtffile, outfile):

    bamfile = pysam.AlignmentFile(bamfile)

    counts = Counter()

    for transcript in GTF.transcript_iterator(
            GTF.iterator(IOTools.openFile(gtffile))):

        introns = GTF.toIntronIntervals(transcript)

        for intron in introns:
            reads = bamfile.fetch(
                reference=transcript[0].contig,
                start=intron[0], end=intron[1])

            for read in reads:
                if 'N' in read.cigarstring:
                    blocks = read.get_blocks()
                    starts, ends = zip(*blocks)
                    if intron[0] in ends and intron[1] in starts:
                        counts["Exon_Exon"] += 1
                    else:
                        counts["spliced_uncounted"] += 1
                elif (read.reference_start <= intron[0] - 3
                      and read.reference_end >= intron[0] + 3):
                    if transcript[0].strand == "+":
                        counts["Exon_Intron"] += 1
                    else:
                        counts["Intron_Exon"] += 1
                elif (read.reference_start <= intron[1] - 3
                      and read.reference_end >= intron[1] + 3):
                    if transcript[0].strand == "+":
                        counts["Intron_Exon"] += 1
                    else:
                        counts["Exon_Intron"] += 1
                else:
                    counts["unspliced_uncounted"] += 1

    header = ["Exon_Exon",
              "Exon_Intron",
              "Intron_Exon",
              "spliced_uncounted",
              "unspliced_uncounted"]

    with IOTools.openFile(outfile, "w") as outf:

        outf.write("\t".join(header)+"\n")
        outf.write("\t".join(map(str, [counts[col] for col in header]))
                   + "\n")

Unfortunately this uses a bunch of libraries you may or may not have, including CGAT (for the GTF parser and IOTools package) and pysam.
Once you've got these statistics you can calculate the "splicing index" as
the log2 ratio of 2 times the number of spliced reads divided by the number of reads overlapping the 3' and 5' ends of introns.

Take this first comment with a grain of salt, since this isn't an area I've worked in much, but: is binary classification possible? If a gene has 3 introns, and 2 are spliced out but 1 is retained, is this "spliced" or "unspliced". My first impression is that an analysis would be a bit more nuanced than a binary classification.
That said, I'm aware of a tool called Keep Me Around for intron retention analysis (preprint, code). I've never used this software, but it looks like it's actively maintained and was created by a research group with a lot of experience and influence in this area, for whatever that's worth. :-)
UPDATE: It looks like the original post was modified since I submitted my answer, but I'll leave the text about binary classification around for posterity. :-)


text:
What methods should I use from PythonCyc API to query metabolites in BioCyc database?
I am using PythonCyc API in order to write a query for metabolites in BioCyc. The purpose of this API is to communicate with the database software of BioCyc- Pathway Tools. Pathway Tools is in lisp therefore, PythonCyc creates a bridge between python and common lisp. To use this this API one must first create a PGDB object with a specified organism ID (orgid). In the example below I create a PGDB with orgid "meta". After this, I am able to call methods from PythonCyc with this object:
import sys
import pythoncyc
#this creates PGDB object associated with meta(MetaCyc)
meta = pythoncyc.select_organism('meta')
#lists pathways of specified compound
pathways = meta.pathways_of_compound('sucrose')
print pathways

The query for the metabolite above, 'sucrose', provides a list a pathways:
[u'|PWY-7347|', u'|SUCSYN-PWY|', u'|PWY-7238|', u'|PWY-6527|', u'|PWY-5343|', u'|PWY-3801|', u'|PWY-7345|', u'|PWY-862|', u'|SUCUTIL-PWY|', u'|PWY66-373|', u'|PWY-621|', u'|PWY-822|', u'|SUCROSEUTIL2-PWY|', u'|PWY-6525|', u'|PWY-6524|', u'|PWY-5337|', u'|PWY-5384|']

However if I switch this metabolite name to a common amino acid name, such as 'valine':
import sys
import pythoncyc
#this creates PGDB object associated with meta(MetaCyc)
meta = pythoncyc.select_organism('meta')
#lists pathways of specified compound
pathways = meta.pathways_of_compound('valine')
print pathways

The query provides an error message that states "non coercible frame", which means it cannot find this ID. This is a common metabolite, just as sucrose, which should very well have pathway entries in database, however there are none found with this method. I have also tried synonyms of valine and other methods such as "all_pathways(cpds)" stated in API, which give me the same error message.
What is the reason for this? Does BioCyc not list any pathways for amino acids? Or am I using an incorrect method to access amino acid information?
Link to PythonCyc API

Hi astrika, thanks for your question and welcome to Bioinformatics Stack Exchange. I've had a glance over this question, and it looks like it has the potential to get good answers: it is well-researched, with a good story, and a specific question. I look forward to seeing more questions from you in the future.

UPD: As suggested by Llopis, maybe the correct identifier is VAL?
I think there is no metabolite with such name in the database. I tried their search, the correct name seems to be L-valine. I expect that they use similar naming for other amino acids. If you cannot find something, you can use their compound search through the website.

I have tried 'L-valine' in the script above and it has also gave me the same error. I am aware it is possible to search in the site, however for my purpose I must write a query. My purpose is to obtain possible pathway networks for a given metabolics dataset

@astrika maybe the name is Protein-L-valine, or VAL. One thing is how they store things an another is how they show them. Sorry I didn't test how to search using python

Pathway-tools has a unique identifier (called a Frame-ID) for each metabolite in its database.  You can get the list of all compounds with this api call:
all_cpds = meta.get_class_all_instances('|Compounds|')
for each compound, you can also get its common name, or all its synonyms:

cpd_names = {}
for cpd_frame_id in all_cpds:
      cpd_names[cpd_frame_id] = meta.get_name_string(cpd_frame_id)


text:
What is the best method to estimate a phylogenetic tree from a large dataset of >1000 loci and >100 species
I have a large phylogenomic alignment of >1000 loci (each locus is ~1000bp), and >100 species. I have relatively little missing data (<10%).
I want to estimate a maximum-likelihood phylogenetic tree from this data, with measures of statistical support on each node.
There are many phylogenetics programs that claim to be able to analyse datasets like this (e.g. RAxML, ExaML, IQtree, FastTree, PhyML etc). Given that I have access to a big server (512GB RAM, 56 cores), what are the pros and cons of each program.
Which programme likely to give me the most accurate estimate of the ML tree for a dataset of this size within a reasonable timeframe?

You want us to review all the phylogenetic programs or just the ones you listed: RAxML, ExaML, IQtree, FastTree, PhyML (Which still I think is quite broad). How do you measure the accuracy of the ML tree?

I'm interested in any and all opinions, evidence, and links to comparisons between any software capable of estimating ML trees from large datasets like this. There would be many ways of measuring accuracy, including: (i) evidence from simulation; (ii) comparisons of Likelihood scores of estimated trees in a common framework.

RAxML is the goto program. The problem with FastTree is that its accuracy is limited to the datasets used, i.e. yours might be an outlier.

This paper claims that FastTree is almost as accurate as RAxML, while being much faster. You just have to be careful, however, that the support values output by FastTree are not bootstrap values, they are based on the Shimodaira-Hasegawa test. (Also, see this comment for the case you have very short branch lengths).   [update: However, according to the recent comparison paper mentioned below FastTree performed quite poorly in comparison to RAxML or IQ-tree.]
From what I understand, you should use ExaML only if your data is too large to be handled by RAxML in a single node. ExaML should perform like RAxML but with some parallelization overhead. For all effects I treat them as the same. I don't know of relevant advantages of phyML over RAxML (for me, it's easier to use but I am very used to phyML).
I am not familiar with IQ-tree, but its authors claim that even given the same time as RAxML or phyML, IQ-tree already finds better likelihoods more often than not (although by default it takes a bit longer to converge). A recent comparison between all these programs favoured IQ-TREE for both single-gene and concatenation analyses (with RAxML very close). It may also estimate branch support through a SH-like test only, but I'm not sure.   [update: IQ-tree offers 3 measures of support, standard bootstrap, aLRT, and ultrafast bootstrap. See OP's comment below for details.]
However, since you have few missing data, you might also want to try a single-locus tree inference followed by gene tree clustering (using treescape or treeCL) to see how spread your data is, or to see the effect of removal of outliers, or to use ideas similar to statistical binning.

In the recent comparison you mention, I note that RAxML only performed better than IQtree when they did 10 independent searches per replicate. So one might just favour IQtree here too, noting that it's sensible to do a bunch of independent searches. Also FastTree performed very poorly in all comparisons of that ms.

Measures of support in IQtree: (i) standard bootstrap; (ii) aLRT (approximate likelihood ratio test, which is something like asking whether a given branch length is significantly >0); (iii) Ultrafast bootstrap (not the same interpretation as a standard bootstrap, more like a posterior probability if I've understood it right).

Note that likelihood comparisons between program outcomes also depend on the choice of model and data partitioning. I haven't looked at the details, but some programs may implement models that are not available in others. This can be a criterion of choice.

Thank you for the comments @roblanf, I took the liberty of updating my answer with this info.

note also that there is a newly released re-write of RAxML available here. It combines the best parts of RAxML and ExaML, while being faster in general. However not all current RAxML features have been implemented yet.

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


text:
How can I combine a kinship matrix with subset individuals when using rvtests?
I've been processing data for a consortium project that uses the rvtests toolkit. Our data analysis process uses a kinship matrix and the calculation of inverse-normalised transformed statistics after covariate adjustment (for 1kG-imputed data from a SNPchip). One of the operations we have been asked to do is to analyse males and females separately for a number of different phenotypes.
Here's an example of the command line that I'm running:
rvtest --inVcf chr${ch}.imputed.poly.vcf.gz \
 --pheno ../scripts/phenotypes.raw.2017-May-01.ped \
 --pheno-name ${traitName} \
 --out rv_results/${traitName}_FEMALE_chr${ch} \
 --kinship kinship_matrix.kinship \
 --meta score,cov[windowSize=500000] \
 --peopleIncludeFile ../scripts/female.iids.txt \
 --covar ../scripts/covar.2017-May-11.ped \
 --covar-name ${covs} --useResidualAsPhenotype --dosage DS

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

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

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

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


text:
Rule Extraction from nnet results
I used a script in R language that uses nnet library to predict promoter bacteria and i would like to know how to extract rules from this neural network results.
As my input of the neural network i have n positive examples of promoter and n false examples. I use as input these examples and the FASTA file with the genome of the bacteria. As a result, i have a value of 0 to 1 for each network tested corresponding to its learning. I want to discover roles that can improve my network in future experiments in the same bacteria.
Which algorithms or softwares could i use?

Welcome to Bioinformatics! Could you post an example? What kind of rules do you expect? What have you found/tried (to avoid proposing things you already did/search) ?

Please [edit] it into the question. What have you looked for to solve your problem?

I am starting to search for possible solutions to rule extraction.

That would be the first thing before other people start to freely invest their time

Short answer: you can't.
Neural networks use positive and negative examples to add weights to the neural network architecture that is provided to it. Trying to deconvolve the meaning behind the weights is nigh-on impossible except for the simplest perceptron.
This is common to many machine learning algorithms: they work very much like black boxes. One exception is decision trees. They will report what features are used to classify the positive vs negative datasets.
However, as with all ML methods, you need to be very careful of your training datasets. This is especially true of motif searching and especially for negative datasets. It's quite easy to find a negative dataset which is completely inappropriate for learning and give you misleading accuracies.
Like @user172818, I would try traditional methods as they can work well, again, if given appropriate data. The MEME-suite would be a good start.

You may have a look at the paper under a link in my answer. For CNN, the weights in the first convolutional layer is interpretable.

Agree. Which is what I said about "the simplest perceptron".

R's nnet package only supports fully connected neural networks with one hidden layer. This is the most primitive type of network. I doubt it will work well for promotors finding. In addition, such network won't give you useful interpretations.
If you want to explore neural networks, you should use a one-dimension convolution layer as is described in this paper. This layer effectively represents a positional weight matrix (PWM). You can know which motifs are used. With backtracking, you can also identify the coordinates of motifs. To deploy such a model, you need to be fairly familiar with deep learning and to learn one deep learning framework. For framework, you may start with keras. Once you get used to keras, you can implement your full promotor finder in less than 100 lines of python code. Alternatively, you may try dragonn. It is supposed to simplify deploying models for DNA sequences. I have no experience with it, though.
For promotor finding, it is also worth trying traditional methods. It is interesting that few/no neuralNet-based publications have evaluated traditional methods, probably because many ML people know little about classical motif finding. These methods could work well, too.


text:
How does Cellminer's "Cross-correlations of transcripts, drugs, and microRNAs" work
I have seen one of Cellminer tools. I am not sure how do they calculate the cross correlation of the genes, what does it actually mean? Based on what databases?
For example if I take their example input (HUGO format):
abcb1
BRCA2
CNBP

I get the following cross correlation matrix (based on what signals).
Identifier      abcb1   BRCA2   CNBP
abcb1           1      -0.142   0.069
BRCA2          -0.142   1       0.176
CNBP           0.069    0.176   1

You're probably best asking the authors

In the article describing Cellminer says about the cross-correlations:

... a cross-correlation table of the resultant z-scores can be
  generated ...

On page 3503
And the z-scores are defined previously as:

The tool output includes relative transcript intensity presented
  as z-scores...


text:
How to detect a mutation and predict its consequence?
I read lecture notes about mutations, and wondered what kind of algorithms are there to detect mutations?
How do we know if the gene is mutated or whether it is a sequencing error?
I saw this thread which is related, but I am not sure how to use the CIGAR information, e.g., is it 100% accurate? What is the underline mechanism to detect mutation? Is there a way to predict what the mutation will cause?

I'll follow up to the great answer from Kamil S Jaron:
Regarding predicting what the variant ("mutation" is a very loaded term) will do, there are a variety of tools. Chief among these are annovar and VEP. The general idea behind these is to classify the variants according to their overlap with genes, which codons they change (if any), how big that change is (e.g., changes in charge are more likely detrimental) and so on. One could also consider conservation, since if a position is highly conserved then changes in it are more likely to be detrimental.
If you really want to predict how a variant will change a protein's function then that usually requires prior knowledge about the proteins in question. Eventually someone will use machine learning to cull the literature and provide good predictions, but I haven't seen that yet.

Nice summary, to add, PhastCons is the standard for conservation scores, although not the best as the Sipel lab has many newer tools alternatives and GERP++ is also better IMO. Also if you're looking at just a few individual loci, nothing beats a stare at a genome browser like UCSC

Part 1 : how to detect mutations
The keywords you are searching for are "variant calling". Basically you have to map sequencing reads to a reference genome (or gene) and then estimate for each position of the genome if the observed difference of mapped reads and the reference is more likely a sequencing error or a mutation (in genomic glossary - variant).
Popular tools for variant calling are GATK, FreeBayes or bcftools (previously part as Samtools package).
The question you linked asks for quick alternatives to simply looking for variants. Indeed, you can just visualise the mapped reads to the reference sequence and see if the variant is there or not. CIGAR is just a notation of read alignment used in sam files (files with mapped reads), you can find good explanation of CIGAR strings here.
The follow-up about "How to estimate an effect of mutation" is in @DevonRyan 's awesome answer.

Part 3

I thought I'd add a part 3 . Having found via VEP a truncation (missense or frameshift) or a missense, how do you find what are the consequences at the protein level?

The effect a mutation has on a protein can be diverse and how much imbalance in function can the cell's equilibrium absorb depends on the component. As a result you get recessive, dominant and embryonically lethal mutations.
In this figure are shown some of the effects, divided into core and surface mutations:

Core variants
Core mutations are the most black-and-white to analyse, but require a protein structure available. A severe truncation and a highly deleterious missense will result in a protein that will likely be degraded. The difference in folding Gibbs free energy (∆∆G, where $\Delta \Delta G = \Delta G_{\rm{mutant}} – \Delta G_{\rm{WT}}$)) can be calculated for a missense, but how large this must be to result in a certain percentage being degraded depends on factors that cannot be calculated analytically. Alternatively, the Gibbs free energy of ligand binding is affected resulting in loss of catalysis.
To calculate the ∆∆G effect (stability), there are many tools:

Online SDM uses ML to estimate the ∆∆G quickly (do note that positive kcal/mol normally is a loss of stability, but here the sign is inverted so positive is a gain in stability
Rosetta toolkit. A highly flexible panel of movers  (samplers) can calculate ∆∆G
FoldX
MD runs
Missense3D assesses a variant based on several structural criteria

In the case of ligand binding modelling as above but using the ligand bound structure and calculating the binding energy as the difference between bound and unbound (far far away). For altered catalysis, modelling as above but using the transition state and assessing the difference in the transition hill can be done. However, an active site variant is very likely to affect either ligand binding, transition state binding (catalysis) or product binding that it's surely deleterious.
This is all (mostly) straightforward. If the domain is a regulatory domain however, there is loss of inhibition, which results in a gain of function (dominant variant).
Surface variants
Then stuff gets messier when you get to surface variants and these can cause dominant effects. However, if you are close to a phosphorylation, ubiquitination or other PTM site, as identified in phosphosite plus, you likely have a case of deregulation. For linear motifs, ELM is useful, but requires a lot of literature scanning to confirm any find. Alternative a loss of an interface may break the protein function, however, what binds where can only be ascertained by available structures, cross-linking mass spec or in some cases well defined Western experiments, so this data is rare.
For this type of variant there are two pipelines:

Miscast collates information on the protein
VENUS (beta) calculates ∆∆G and finds any relevant close-by PTM or gnomAD variants [disclaimer: this is mine].

Footnote
As the variants can have a variety of effects, none of these score (unlike the genetics parts), but instead try to inform as to what may be going on.

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


text:
Why Bioconductor?
What are the advantages of having Bioconductor, for the bioinformatics community?
I've read the 'About' section and skimmed the paper, but still cannot really answer this.
I understand Bioconductor is released twice a year (unlike R), but if I want to use the latest version of a package, I'll have to use the dev version anyway. A stamp of approval could be achieved much easier with a tag or something, so it sounds just like an extra (and unnecessary) layer to maintain.
Related to this, what are the advantages as a developer to have your package accepted into Bioconductor?

+1 … the idea of having centralised repositories is becoming more and more out of fashion because there’s no clear evidence for its supposed benefits.

Is the question "What are the advantages of having Bioconductor over CRAN?" or "What are the advantages of having any kind of central repository? - lets take Bioconductor as example"?

@KamilSJaron My question is neither -- or both: What are the advantages of having Bioconductor, or any similar central repository, over CRAN? If anything, the focus is on having any R repo for bioinformatics, as opposed to a website, where the approved packages are listed.

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

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

In the end, installing dev versions of R packages is in my opinion very bad practise for reproducible science. If developers delete GitHub repo, commit hash you have used won't be enough to get the code.
Benefits of central repository for developers
I forgot about the advantages for you as developer to submit your package to Bioconductor:

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

Bioconductor specific advantages over CRAN
I see the big advantage in the community support page, provided by Bioconductor. @Llopis' comprehensive elaboration.

I think you're missing: data packages. CRAN has quite strict file size limits (5 MB) and biology needs e.g. genomes. Your points 1 and 3 are debatable (1: bio yes, otherwise no; 3: I'd argue it's harder, needs an external installer).

You do not want write package with genomic data stored inside anyway - we rather write packages for pulling data from databases. But point taken.
Now I see I made a smart move, when I numbered both overall advantages for community and advantages for developers.

Ok, I turned upper numbers to bullet points and I made clear that I am speaking about central repositories in general.

@KamilSJaron Regarding data packages, they’re an established best practice in Bioconductor (see BSGenome) though I tend to agree: I dislike how the Bioconductor data packages clutter my disk.

Here is a list of the advantages of having Bioconductor for the bioinformatic community:

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

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

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

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

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

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

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

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

Did you try submitting to CRAN? As far as I understand they also run tests, and perform a package review; some people even complain that CRAN people are too pedantic.

Regarding what the advantage is to you as a developer of having a bioconductor package rather than using CRAN:

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

I have been using a lot of packages from both CRAN and Bioconductor, I haven't noticed that CRAN packages are worse in any sense. I even noticed that documentation of  some packages from Bioconductor was worse (at least that was the state several years ago).

You've had nice luck with CRAN packages then.

Pretty sure the quality hierarchy is wrong: if a package is maintained on Github, odds are its quality is in fact superior than that of the average CRAN or Bioconductor package. This might soon change once Github becomes more widespread amongst beginning R programmers but I’m pretty certain it currently still holds.

@KonradRudolph the "if" in your first sentence is a big one.

@DevonRyan My mistake, I meant "is hosted on". There's a risk it's unmaintained (but the same is true for centralised repositories) but the odds of quality given location = Github are still high. I'll grant that it's possible that the variance of quality is also higher on Github, but even for that I'm not sure.

@KonradRudolph the most of Bioconductor packages are maintained on GitHub (at least the most of I have used). I have GitHub R package as well, but it is definitely not passing any tests, it is working on my computer in the way that I use it. No being on GitHub is not any sign of quality... but true, if the package is only there, nobody will know about it, so there won't be any bug reports.

The most important reason is that Bioconductor has growing set of common data structures and base packages. If package X and package Y needs to work with the same type of data, having a common data structure in the core Bioconductor package Z makes our lives so much easier. I could do something in package X, take out the results and keep working on my data with package Y. This works because I'm using versions of package X and Y that are compatible with the data structure defined in package Z. The Bioconductor team makes sure that all packages use the common data structures and already existing packages where possible, so that we're all on the same page and so that people don't reinvent the wheel again and again.
Also, there's a sophisticated process around getting your package accepted to Bioconductor. This ensures that packages use data structures and functions already available in other Bioconductor packages (where reasonable), it ensures that packages are well written, that they have good documentation, and that they are well tested.


text:
Getting structure and sequence related to PDB IDs
I have two kinds of interactions: transient and stable. We are supposed to work on stable interaction, like interactions between two monomers in a heterodimer.
In a heterodimer there are two chains in which there are some residues between monomers which take part in the interaction of monomers and building the structure of heterodimer, and there are other residues interacting inside each monomer for building the structure of that monomer. There are some PDB IDs related to heterodimers in PDB.
Now imagine that we are trying to show which residues are physically interacting between two monomers in a heterodimer, and then we are trying two show the site of these physically interacting residues in the sequence related to that special structure. For that, we need to find the sequence related to that structure and we can download it from PDB.
But if we want to make a multiple sequence alignment of that sequence, first we need to find the refseq sequence of that and then we must blast the sequence and find homologues of that sequence.
Here there are some challenges:

we must find residues interacting in the structure of heterodimer between two monomers.
we must find the sequences of monomer chains of the heterodimer and map the interacting residues in the structure to the sequence.
we know that when a structure is solved and its monomers are sequenced, may be they are not able to completely sequence the heterodimer and residue numbers are maybe different from the related refseq sequence.

I would like to know how can we find the related refseq sequence (or sequences, in the case of heterodimer) and then how can we map the physically interacting residues in the structure of PDB ID to refseq sequence (finding the sites of these interacting residues in that refseq sequence) regarding three challenges mentioned above.
Imagine that I have 100 PDB IDs of 100 nonredundant heterodimers and I would like to find and download the structures, sequence and finally refseq sequence. What should I do?

So you want to map IDs to RefSeq IDs? Do you want protein or gene IDs?   Have you tried anything? If so, please [edit] and tell us what you tried and how it failed so we don't just give you the same solutions. What species are the PDBs from? What to PPIs have to do with it? PDB isn't a database of interactions. Could you [edit] your question and add a bit more detail so we can understand what you need, please?

You can find mappings PDB <-> Uniprot and Uniprot <-> RefSeq, but if you just want to use Blast, can't you use the sequence from PDB?

That works but I do not just want to do blast. I want to do multiple sequence alignment after that (after finding homologous sequences by running blast). The main challenge is mapping residue-residue interaction in PDB structure to refseq sequences and finding residues sites which participate in interaction. what can I do for that?

I think you could use the sequence from PDB also in MSA, and the alignment will show what the residues of interest correspond to in other sequences (but take what I write with a grain of salt, I don't have first-hand experience here)

You can download seqs and structures based on a list of PDB ids using http://www.rcsb.org/pdb/download/download.do#FASTA

The post might seem old, but I wish to share my thaughts based on how I solved similar problem. Here it goes:
I will use an example to illustrate. Other posts have already explained how you could get the protein ID from PDB and so on. I used the database PDBSUM (www.ebi.ac.uk/pdbsum/) and with protein PDB ID: 1KJF.
By following these steps I got to the results below:
Go to --> www.ebi.ac.uk/pdbsum/ -->(put in PDB ID) or (others e.g UniProt ID) --> Enter --> Click on "Ligands" -->(a LIGPLOT appears at the botoom of the page).
You can as well search just by pasting in the sequence.
Results:
You can see that 1KJF is made up of 3 chains A, B & P and both chains have clearly anntated and interacting residues. You can now see clearly read inidividual types of interactions.
Another good option is http://protein.bio.unipd.it/ring. It shows you a graph plot of all inteacting residues and how they are linked to each other. You can adjust the default result settings to meet your needs.
Hope it helps. :)

Yes.  'pdbsum' and 'ring' are definitely two useful tools for fast analytical work.


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

Dear jaslibra. Thank you for your nice comment

My reputation isn't high enough to make a comment, but have you considered a modelling software such as PyMol or Chimera? You can execute a command line script with softwares such as these and these have the functionality to predict H-bonds.
In the PDB files there is information regarding the GenBank and UNIPROT info of the protein, and you can find more information here:
http://www.wwpdb.org/documentation/file-format-content/format33/sect3.html.
If you want RefSeq ID, you can use GenBank info you extract and convert by piping into a biomaRt (Bioconducor) script, although maybe there is a more direct way.
As for finding which residues in the sequence are interacting, I'm fairly certain that the software packages like PyMol or Chimera will give you information about the interacting residues.
From what I understand, you will want to filter out only the residues that interact and are on different chains (e.g, Chain A vs. Chain B). This process will require you to work with the syntax of the different software packages when you are scripting. For example, Chimera has a basic primer here:
https://www.cgl.ucsf.edu/chimera/current/docs/ProgrammersGuide/basicPrimer.html
Note that I only gave a reference to Chimera because I am more comfortable with that then other packages and there is likely a more efficient way to do this. With Chimera, you can find the position within the chain of residues participating in the interaction and then you'd store these in a file to do whatever you want later.
I hope this helps.

Hi jaslibra, thanks for your comment, and welcome to Bioinformatics Stack Exchange. This comment is well-structured, and looks long enough and verbose enough to be an answer for this question, so I haven't flagged it as a comment. I look forward to reading more comments and answers such as this in the future.


text:
What are the standard ways to visualize protein-protein or gene-gene interactions?
I would like to visualize this interaction list.  Is there an online/web based way to do it? Is there a way to analyze the data?
Assuming I have exported the list into lines such as:
2  A2M 348 APOE

What tools have you found (To avoid suggesting tools you already tried) ?
What do you want to analyze from it (what is your biological question)?

since I am not yet allowed to comment, I will have to pots this as an answer.
I guess Cytoscape (http://cytoscape.org/) would be a nice way to analyse such interaction lists/networks (still depends a bit on what you plan to do).
However, for this you will have to preprocess the downloadable excel sheet a bit.

save as tab-separated file (tsv)
run the below python3 script to format into an easily importable format for cytoscape. You can run the script as python3 script.py <path-to-tsv-file> <delimeter for cytoscape> . The new format is written to std out on the console, therefore maybe you need to write the output to a new file. A sample call could look like python3 testbench.py human_interactome_2017-06-12.tsv.csv ";". Make sure to add the delimeter in quotes.
you can then import the new file into cytoscape (import network as table) where the first column is your source, the third column your target and the middle column your interaction type.
Within cytoscape you can then adapt your edge style to vary the line type in a discrete mapping depending on the value of the interaction.

Screenshot from cytoscape:

Code below:
import sys

HIIIidx = 4
HII5idx = 5
venkatesanidx = 6
yuidx = 7
Litidx = 8

print(sys.argv)

if len(sys.argv) <= 1:
    print("Usage: script.py file delimeter")

delim = ';'
inFile = sys.argv[1]

if len(sys.argv) > 2:
    delim = sys.argv[2]

def printInteraction(syma, symb, inter):
    print(syma + delim + inter + delim + symb)

with open( inFile , 'r') as file:

    file.readline()

    for line in file:

        aline = line.strip('\n').split('\t')

        if aline[HIIIidx] == '1':
            printInteraction(aline[1], aline[3], 'HIII')

        if aline[HII5idx] == '1':
            printInteraction(aline[1], aline[3], 'HII5')

        if aline[venkatesanidx] == '1':
            printInteraction(aline[1], aline[3], 'VEN')

        if aline[yuidx] == '1':
            printInteraction(aline[1], aline[3], 'YU')

        if aline[Litidx] == '1':
            printInteraction(aline[1], aline[3], 'LIT')

Instead of the lineCount counter, if you want to skip the first line you can do: with open(inFile, "r") as my_file: my_file.readline(); for line in my_file:. The readline will advance to the file to its next line. It is considered best practice to open a file inside a with context manager.

Also, you can do directly aline = line.strip("\n").split("\t"), since you don't re-use the modified line (and I would call this variable fields rather than aline).

Indeed the code was written down quite fast, however, I added your comments to the solution.

Your list has two identifiers for the same node per line. In order to use it, you will need to change that. If you want to use the gene name (2nd and 4th fields, in your example), just run:
awk 'print $2,$4' netw.txt > netw.gr

If you want to use the Entrez geneIDs instead, run:
awk 'print $1,$3' netw.txt > netw.gr

Then, as others already mentioned, install Cytoscape, launch it and import your gene list:

In the welcome screen, select "From Network File...":

Then, select your network file (netw.gr in the example above) and choose "Advanced Options":

Set the delimiter to SPACE and uncheck "Use first line as column names":

Finally, click on the header of each column and set the first to "source node" and the second to "target node":

Now, click "OK" and you will have successfully imported your network into Cytoscape.

Usually interactions are represented in a graph, with edges representing the interactions (colors, thickness or values on the edges represent different properties of those edges).
As far as I know, there aren't on line programs to represent this kind of information. For an off line representation, you can use the igraph package in R, or Cytoscape program.

Several have mentioned the excellent Cytoscape desktop package for visualising interactions as a graph.  If you want to build your own web-based visualisations (rather than using a pre-built website to visualise them) then the Cytoscape javascript library might fulfil your needs.


text:
How do I generate a variant list (i.e. VCF file) using Illumina reads from a human genome?
This is a problem I have to solve frequently, and I'd be interested in knowing what other methods people use to solve the same problem.
About twice a year, I get asked to determine variants from Illumina reads, usually from either mouse or human. These are things that have good reference genomes (e.g. human genome from ensembl, mouse genome from ensembl), and are frequently tested out on various bioinformatics tools.
I have at my disposal three computers, in order of usage preference: a laptop with 4 processing threads and 8GB of memory, a desktop with 12 processing threads and 64GB of memory, and a server with 24 processing threads and 256GB of memory. Occasionally I get access to better computers when I need to do lots of processing, or to do things quickly, but these comprise my basic bioinformatics team.
I've had whole-genome projects that involved processing data from 4 individuals, up to projects with about 120 individuals, with each individual having about 10-40M paired-end reads (100-125bp each). Beyond this level, I suspect other researchers will typically have their own bioinformatics team and/or capabilities, which is why I don't get any of the larger-scale jobs that I hear about from salaried bioinformaticians.
What is the standard approach that you would use to process these reads into reference-mapped variant data (i.e. gzipped VCF) for use in downstream analysis tools (e.g. VEP)?

1. Adapter Trimming
One of the first things I do after encountering a set of reads is to remove the adapter sequences from the start and end of reads. Most basecalling software includes some amount of built-in adapter trimming, but it is almost always the case that some adapter sequence will remain. Removing adapters is helpful for mapping because it reduces the effort the mapper needs to go through to make a match (which may help some borderline reads to be mapped properly).
My preferred adapter trimmer is Trimmomatic. It has the ability to search for palindromic adapter read through, and includes an adaptive sliding window option. Trimmomatic is threaded, and seems to work reasonably fast when run through files one at a time in threaded mode. Here's an example script that I might run to trim reads from some samples. This puts trimmed FASTQ files into a trimmed subdirectory:
TRIMPATH=/data/all/programs/trimmomatic/Trimmomatic-0.36;
JARPATH="${TRIMPATH}/trimmomatic-0.36.jar";
ADAPTfile="${TRIMPATH}/adapters/TruSeq3-PE-2.fa";
CMDS=" ILLUMINACLIP:${ADAPTfile}:2:30:10:5:true LEADING:3 TRAILING:3  SLIDINGWINDOW:10:20 MINLEN:40 ";

mkdir -p trimmed;

for INFILE1 in *_1.fastq.gz; do
  base=$(basename "$INFILE1" _1.fastq.gz);
  echo ${base};
  INFILE2="${base}_2.fastq.gz";
  OUTFILE_P1="trimmed/${base}_P1.fastq.gz";
  OUTFILE_P2="trimmed/${base}_P2.fastq.gz";
  OUTFILE_U1="trimmed/${base}_U1.fastq.gz";
  OUTFILE_U2="trimmed/${base}_U2.fastq.gz";
  java -jar "${JARPATH}" PE -threads 20 -phred33 \
        "${INFILE1}"  "${INFILE2}" \
        "${OUTFILE_P1}" "${OUTFILE_U1}" "${OUTFILE_P2}"  "${OUTFILE_U2}" \
        "${CMDS}";
done;

Trimmomatic is somewhat sensitive to bad input data, and will loudly fail (i.e. stop running) if it sees quality strings of different length to sequence strings:
Exception in thread "Thread-1" java.lang.RuntimeException: Sequence and quality length don't match: 'TACATGGCCCTGAAATGACTTTCACCCAGGCAACCAGTGCCCCCTGTATAGACACATGCCTTGGGCGCTCCCCACCCTTCCTCGCGTGGCCACACCTCTGT' vs '-AAFFJFJJ7A-FJJJFJFJFJJFFA-FFAF-<A-<FFFJA-<-A7-F<<FFJJAJJJJJJJJJ--<--7-7AJ7<AAJA--J7<ACTGCTGTGGGGCACCCAGCCCCCCAGATAGCCTGGCAGAAGGATGGGGGCACAGACTTCCCAGCTGCACGGGAGAGAC'
        at org.usadellab.trimmomatic.fastq.FastqRecord.<init>(FastqRecord.java:25)
        at org.usadellab.trimmomatic.fastq.FastqParser.parseOne(FastqParser.java:89)
        at org.usadellab.trimmomatic.fastq.FastqParser.next(FastqParser.java:179)
        at org.usadellab.trimmomatic.threading.ParserWorker.run(ParserWorker.java:42)
        at java.lang.Thread.run(Thread.java:745)

2. Read Mapping
Read mapping finds the most likely location for a read within a target genome.
There are extremely fast approximate mappers available, but these don't yet work for variant calling, and an exact mapping approach is necessary. My current preferred mapper is HISAT2. I use this instead of BWA due to the double-read issue, and because of the local variant-aware mapping. HISAT2 is made by the same computing group as Bowtie2 and Tophat2 (JHUCCB), and is their recommended tool for replacing those other programs.

HISAT2 is a successor to both HISAT and TopHat2. We recommend that HISAT and TopHat2 users switch to HISAT2.

The HISAT2 page includes a link to a genomic index file for the human genome, including SNP variants, which I use when mapping human reads. HISAT2 is also threaded, so I run it on the files one at a time and pipe through samtools to create a sorted BAM file:
mkdir mapped
for r1 in trimmed/*_P1.fastq.gz;
  do base=$(basename "${x}" _P1.fastq.gz);
  r2="trimmed/${base}_P2.fastq.gz";
  echo ${base};
  hisat2 -p 20 -t -x /data/all/genomes/hsap/hisat_grch38_snp/genome_snp -1 \
    "${r1}" -2 "${r2}" 2>"mapped/hisat2_${y}_vs_grch38_stderr.txt" | \
    samtools sort > "mapped/hisat2_${y}_vs_grch38.bam";
done

3. Variant Calling
I currently use samtools/bcftools to do this, but would be interested in other opinions as we recently had a grant reviewer response that samtools was a dinosaur program and better approaches were available. Samtools doesn't currently work with threads for variant calling, so I save the commands to a text file and then run them through GNU Parallel. This requires the genome files to be first downloaded and indexed via samtools faidx. This produces a set of gzipped VCF files in the variants directory:
mkdir -p variants;
(for x in mapped/*.bam;
   do echo samtools mpileup -v -f /data/all/genomes/hsap/hisat_grch38_snp/Homo_sapiens.GRCh38.dna.primary_assembly.fa "${x}" \| \
   bcftools call --ploidy GRCh38 -v -m -O z -o variants/$(basename "${x}" .bam).vcf.gz; done) > call_jobs.txt
cat call_jobs.txt | parallel -u;

As a slow, but more accurate alternative, variants can be called for all samples at the same time:
samtools mpileup -v -f /data/all/genomes/hsap/hisat_grch38_snp/Homo_sapiens.GRCh38.dna.primary_assembly.fa mapped/*.bam | \
  bcftools call --ploidy GRCh38 -v -m -O z -o variants/hisat2_allCalled_vs_grch38.vcf.gz

If a faster processing is desired, then the variant calling can be parallelised across chromosomes using parallel and merged afterwards using bcftools norm. The implementation of this is left as an exercise to the reader.
Aside: I see that GATK's best practise guidelines call per-sample with HaplotypeCaller, then process those VCF files as aggregated data. @Geraldine-vdauwera says this two-step joint calling workflow was designed to be scalable for any size cohort, so you can apply it pretty much the same way whether you have two samples or twenty thousand. Once you go higher than that there are some additional tweaks to keep scaling up.

The GATK is probably the most widely used program for this. There are several workflows you can use depending on what types of variants you're interested in (germline or somatic, short variants, copy number, structural). The most common use case is Germline short variant discovery (SNPs and Indels), for which the overall workflow is described here:
https://gatk.broadinstitute.org/hc/en-us/articles/360035535932
The GATK team at the Broad Institute now provides workflow scripts written in WDL (Worfklow Description Language) that can be run on pretty much any platform with the Cromwell engine. The workflows include data pre-processing, including mapping (alignment) and some important cleanup operations, plus variant calling and filtering.
One of the advantages of using those workflows is that the Broad uses them in production for all the data they sequence, which nowadays amounts to abut a new human whole genome every five minutes. That means there's a lot of pressure to make sure the tools are reliable and efficient because any interruption of the pipelining operations would be a big problem. They're also highly standardized so that the outputs will be compatible with other datasets and tools.
Disclosure: I work at the Broad Institute and wrote a book about using GATK (https://oreil.ly/genomics-cloud).


text:
What are the available cloud computing services for bioinformatics?
I'm looking for cloud computing services that can be used for doing bioinformatics. An example I found is InsideDNA and there is Amazon of course.
A little description of these would be appreciated.

Can you mention Google compute as well and further specify what exactly you're interested in knowing? As is I expect this will end up getting closed otherwise.

What do you want to use the computer services for? Do you have access to a university, research institute, or other academic facility? Is there a particular reason why your desktop or laptop is not appropriate?

What kind of computation are you interested in? There are specialised companies for bioinformatics that are outsourcing the computations on google cloud or amazon. I think for different types of computations you should contact a different company...

At the moment this is a purely theoretical question as the computing facilities at my university are more than adequate for RNA-seq data analysis. I'm just looking for a simple list of companies that provide this service, because a simple search does not return many. In my ideal environment I do not depend on pre-installed software (i.e. I can install software), command line (bash) interaction is fast and I can trust the company for backups etc -- so I'm also interested in the location of the provider, how long has it been around and the computing resources they have.

I do not think that this is a good question since there is no one good answer (or at least the answer should a list, not a one provider), so I will just link the providers I know : DNAnexus (private) and Vital-it (academic)

This post is too broad.

Related post: https://www.biostars.org/p/86463/

What kind of bioinformatics analysis you are planning to do?

@arupgsh See my comment above

I put this on hold since asking for lists like this is really not a very good fit for the Q&A format and better suited to something like a forum where a discussion can be had. See this old blog post for an explanation of the general rationale behind considering this type of question off topic.

It really depends on what you are trying to do, but here are a few services that I know of.

GATK on Google Genomics Cloud: Google and the Broad offer a cloud instance tailored to GATK pipelines.
Genomics on Amazon Web Services: I don't think there is anything that makes this unique, but Amazon offers some resources to help get started with genomics/life sciences-centered cloud solutions.
Illumina Bioinformatics: Illumina is working on a whole suite of bioinformatics software for the cloud.
Cancer Genomics Cloud: This is specific to cancer genomics, but I believe Seven Bridges allows you to push all sorts of data into the tool and analyze it.

I have trialed the free version of InsideDNA, and these were my notes:

Cost: \$225/month for a team of 5 with 50TB storage or \$45/month with 10TB storage for individuals (assuming 6 month package: https://insidedna.me/pricing).
Software installed: Around 600 bioinformatic tools available and standard command line tools; some popular tools missing (like CD-HIT), but should be possible to install on request.
Jobs: Maximum of 32 CPUs and 208 RAM per job submission. Test jobs generally worked, although a larger test job failed.
Other points: Command line was sometimes slow, wget queries were slow, and scp was blocked. However, these may be resolvable issues.

Overall, I felt InsideDNA could be useful for groups without their own computational infrastructure and could be used for easily sharing resources between groups. The packages on offer seem not expensive, but I had a few issues, and I don't know how good their sys admin support would be.
I have not used the Amazon service, so can't comment beyond the details on their website. Also there are a few alternative companies, such as Genestack and DNAnexus, but I haven't directly tested them either.

I'm not sure what kinds of bioinformatics tasks you would like to perform, therefore it is difficult to give a good recommendation.
If you're specifically working on statistical genetics, I can recommend Hail [1]. Hail is an open-source tool for analyzing genetics data at the tens of terabyte scale. Most of Hail's users do their science in Jupyter notebooks that are backed by Google Cloud Platform Dataproc clusters. Hail permits you to perform a variety of statistical genetics tasks including:

filtering and aggregation for quality control
subsetting, linear regression, linear mixed model regression, and linear burden testing
utilities for computing various measures of relatedness
principal components analysis
variant splitting
import/export from a variety of formats including PLINK, VCF, and BGEN, and
a python API which enables the use of libraries like matplotlib for plotting analysis results

To learn specifically about using Hail with the Google Cloud Platform and Jupyter notebooks, I strongly recommend Liam's Hail forum post about his cloud-tools repository.
Here's an example, from the Hail tutorial, of using Hail to perform some quality control and display a scatter plot of the first two principal components of the individuals:
from hail import *
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

hc = HailContext()

table = hc.import_table('data/1kg_annotations.txt', impute=True).key_by('Sample')
common_vds = (hc.read('data/1kg.vds')
              .annotate_samples_table(table, root='sa')
              .sample_qc()
              .filter_samples_expr('sa.qc.dpMean >= 4 && sa.qc.callRate >= 0.97')
              .filter_genotypes('''let ab = g.ad[1] / g.ad.sum() in
                         ((g.isHomRef && ab <= 0.1) ||
                          (g.isHet && ab >= 0.25 && ab <= 0.75) ||
                          (g.isHomVar && ab >= 0.9))''')
              .variant_qc()
              .filter_variants_expr('va.qc.AF > 0.01')
              .ld_prune(memory_per_core=512, num_cores=4))

pca = common_vds.pca('sa.pca', k=5, eigenvalues='global.eigen')
pca_table = pca.samples_table().to_pandas()

colors = {'AFR': 'green', 'AMR': 'red', 'EAS': 'black', 'EUR': 'blue', 'SAS': 'cyan'}
plt.scatter(pca_table["sa.pca.PC1"], pca_table["sa.pca.PC2"],
            c = pca_table["sa.SuperPopulation"].map(colors),
            alpha = .5)
plt.xlim(-0.6, 0.6)
plt.xlabel("PC1")
plt.ylabel("PC2")
legend_entries = [mpatches.Patch(color=c, label=pheno) for pheno, c in colors.items()]
plt.legend(handles=legend_entries, loc=2)
plt.show()

[1] Disclaimer: I work on Hail

Depending on your applications and uses, you might be interested in checking out CyVerse. It is an NSF funded initiative that provides you with data storage, high performance computing resources, and easy access to commonly used tools. As far as I know, it is free to use once you have an account. I also usually encounter it being used with plant and microbial genomics, so not sure how it will work with something like human genomics projects. But might be worth checking out at least. :)
More information: http://www.cyverse.org/about

Google Genomics
Google has an API called Google Genomics.
SNPedia
"SNPedia is a wiki investigating human genetics."
snpedia.com
Promethease
"Promethease is a literature retrieval system that builds a personal DNA report"
promethease.com
DNA Land
"Compare DNA with reference data from different populations"
dna.land
The CyDAS Project
And, there's the CyDAS project which has an API that can analyze ISCN formulae. Per their web site: their API "lets you analyze a Karyotype for virtually all information which can be extracted from karyotypes and the rearrangements therein: gains and losses of chromosomal material, break points, junctions..." It's a free service, but I don't know how up to date it is.

DNAnexus -- http://dnanexus.com
BaseSpace -- http://basespace.illumina.com
Seven Bridges Genomics -- http://www.sbgenomics.com
Curoverse http://curoverse.com
InsideDNA http://insidedna.me/

There are tonnes of them. On top of the excellent ones everybody mentioned

iRods
Arvados
Galaxy

Neither Galaxy nor iRods are cloud computing services. I suspect that Arvados isn't either, but I've never used it.


text:
What kind of "gff" format does bioawk parse?
I was wondering if I could use the gff parsing capability of bioawk to facilitate the parsing of gtf files, and I looked at the following help message:
$ bioawk -c help
bed:
    1:chrom 2:start 3:end 4:name 5:score 6:strand 7:thickstart 8:thickend 9:rgb 10:blockcount 11:blocksizes 12:blockstarts
sam:
    1:qname 2:flag 3:rname 4:pos 5:mapq 6:cigar 7:rnext 8:pnext 9:tlen 10:seq 11:qual
vcf:
    1:chrom 2:pos 3:id 4:ref 5:alt 6:qual 7:filter 8:info
gff:
    1:seqname 2:source 3:feature 4:start 5:end 6:score 7:filter 8:strand 9:group 10:attribute
fastx:
    1:name 2:seq 3:qual 4:comment

I see that there are 10 fields defined for gff parsing.
However, when I look at the pages for gtf and gff in the ensembl website (gff/gtf and gff3), all have 9 fields.
I'm curious about these "filter" and "group" fields. What are they meant for? Are they extracted from some of the columns mentioned in the above pages?

I would consider the description there a bug. The filter is actually the strand, strand is the frame, group is the attribute, and attribute does nothing. These are really meant to be the 9 columns.
Edit: There's a bug report related to this.
Edit 2: I've made a pull request to clarify this and fix the aforementioned bug report.
Edit 3: I realized that I never directly answered the title of your question (mea culpa). bioawk itself will work with gff, gff3, or gtf files. It really is just treating them as tab-separated files with named columns (this is surprisingly convenient, since it's a PITA to remember what column does what).
Edit 4: The PR has been merged. If you install from github then you'll see corrected field names.