text:
How to find hairpin loop structure in a large set of long non coding RNA transcripts
I have done de novo assembly of pair end raw read sequences and resulted transcripts sequences were separated based on coding potential into two categories: long non-coding RNA transcripts and coding RNA transcripts.
I want to find how many of long non coding can make hairpin loop structure and how can I annotate these miRNA which are generated from long non coding RNA in plants.

Welcome to Bioinformatics! Is there anything you tried before posting? Is there any reason why you tagged the question with Bioconductor (Maybe you want to restrict to use a Bioconductor package, or have used one and you can't manage to make it classify your lnRNA)? Also, I think there are two questions here: how to find if a lncRNA can make hairpins and annotate RNA to miRNA. It would be better to focus in just one question (perhaps the first one and later on the second one).

For the first question: "how many of long non coding can make hairpin loop structure?"
The easiest thing to would be to run your sequences through an RNA secondary structure prediction tool. There's a few tools for doing this but the most commonly used are RNAfold from the ViennaRNA package and mfold.
This will give you an output that looks like this, where the parentheses indicate which bases are paired with each other.

GGGCUAUUAGCUCAGUUGGUUAGAGCGCACCCCUGAUAAGGGUGAGGUCGCUGAUUCGAAUUCAGCAUAGCCCA
(((((((..((((.........)))).(((((.......))))).....(((((.......)))))))))))).

To see if there's hairpins, you'll need to annotate the elements within it. Take a look at this question on biostars for an example of how to do this. In short, you can use the forgi library to produce an annotation like this:

(((((((..((((.........)))).(((((.......))))).....(((((.......)))))))))))).
sssssssmmsssshhhhhhhhhssssmssssshhhhhhhsssssmmmmmssssshhhhhhhsssssssssssse

If there's an h in the output annotation, you have a hairpin. The remainder of your question about annotating miRNAs will be left to somebody else.

hairpins are quite easy to identify using the bracket notation; just look for brackets that are in a left-right orientation with nothing in between (e.g. (.), (..), (...)).


text:
For what bioinformatics tasks is Biopython more adapted than Bioperl?
Are there any advantages to learning Biopython instead of learning Bioperl?
Ideally, we would learn both, but someone starting out in bioinformatics may have to choose what to learn first depending on the kind of problems actually encountered.
Are there problems for which Biopython is better than Bioperl (or vice-versa)?

This question is almost entirely opinion based, and therefore generally discouraged on Stack Exchange sites. This is easier to see if you replace 'Biopython' with 'python' and 'Bioperl' with 'perl'.

The only possible objective answer is - advantage of biopython is that it is module in python, on the other hand bioperl has an advantage of being module of perl. So you can reduce your question to "python or perl?" and that's a holy war, not question.

If you can, learn both, and even learn R.
But prefer Python if you are to write code that others will need to read, understand and eventually modify. In my (admittedly biased) experience, Python code tends to be easier to understand.

Currently you could use either but a major question is which platform will others be using in the future. AFAIK Perl is only superior to Python for regex.
Based on the trend I see for new programmers and new software being released: Perl is on the way out and Python is still growing.
https://trends.google.com/trends/explore?date=all&q=bioperl,biopython

“AFAIK Perl is only superior to Python for regex.” — Andy use strict;, which has absolutely no equivalent in Python. There are linters, but nothing substitutes the requirement — provided in Perl by use strict — to explicitly declare identifiers (via my). I like Python but it gives incredibly weak correctness guarantees; even JavaScript surpasses it in this regard.

Regading the perl vs python discussion, there is no final answer which language is better, but I have some advice for you:
Learn the language your colleagues or your advisor use. This way you are able to discuss your code with them and also get help if you run into problems.

I think that the question is more about the libraries than about the language itself.

This usually comes down to religious issues, so let me try and steer it back to more objective grounds:

What language do you know (better)? Use the library for that one.
If you know neither and will be learning a language to use the library, the majority opinion would be that Python is easier to learn.
However, some people say that they "click" with Perl better
Python is what the majority of bioinformaticians are using at the moment (about 60%). And there is a virtue and aid in using what everyone else is.
What are your colleagues / collaborators going to be using?
Broadly, newer libraries tend to be written in Python.
Conversely, many sing the praises of regular expressions and string handling in Perl, citing speed and ease of use.
I think the numerics support in Python may be markedly better than that in Perl.

Agree with agapow's breakdown (even though I'll probably always prefer Ruby over Python). Just as a comment you do find that alot of legacy projects at places like the Sanger use Perl. Depending on what you're aiming to to do initially I wouldn't discount R especially if your project is stats heavy. Later languages to pick up could include Java and other languages that run on the JVM like Clojure / Scala / Groovy which can use Biojava (via Java interop), or something like Rust with Rust-Bio if you want performance. Just make sure you don't start coding in Python 2.7...

I think your 60% number is about right, but I was wondering where you got it from. Do you have a source or is that just your impression?

The 60% is from a survey of published software from a few years a go. No, I don't have the reference.

I agree with @story and @agapow here. As it stands right now BioPerl is still much more mature than BioPython. If there's a tool or method that only exists in the BioPerl library then you should probably use BioPerl to accomplish what you need to accomplish. But it's undeniable that BioPython is progressing very quickly and will soon overtake BioPerl. That being said there are still a lot of plus side regarding Perl in the Bioinformatics realm and it will definitely help to be somewhat competent in it.
So in short my answer is to use both! I find myself mostly using Python and BioPython but if I need to accomplish something using a tool that isn't available in the BioPython library or can be accomplished more quickly using Perl I will look towards BioPerl/Perl.

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 are all the reference files produced by bwa index, and are these dependent upon whether the reference is zipped?
I have indexed a gzipped reference with bwa: bwa index reference.fa.gz, which produces a series of other files reference.fa.gz.{amb,ann,bwt,pac,sa}. These are working fine with bwa alignment.
I have discovered that samtools does not take a gzipped reference, so I am planning to use an unzipped version of the reference for my workflow instead of dealing with two separate representations of the reference.
Are the files reference.fa.gz.{amb,ann,bwt,pac,sa} dependent upon the reference being gzipped? Do I need to reindex the unzipped reference.fa with bwa, or can I just rename the current files to remove the .gz portion of the filename?

You'll get the exact same index (the amb, ann, bwt, pac and sa files) whether the reference is gzipped or not. BWA also makes its own packed reference sequence (the .pac file) so you don't even need the genome around after you index.


text:
Where to download JASPAR TFBS motif bed file?
I am interested in determining if any transcription factor binding site motifs are enriched in some BED files from a DNA methylation experiment.
I am looking for a database that has BED Files containing regions enriched for specific transcription factor binding motifs in homo sapiens.
I know that JASPAR contains TFBS motifs, but I can't find a download for regions of hg19 that are enriched with those motifs.
EDIT:  I suppose I could determine TF motif enriched regions using HOMER if there are no curated databases with this information, but I would rather not.  Your help is appreciated!

You could make your own such file with FIMO, the JASPAR MEME files available via the MEME download site, and your genome build of interest (per-chromosome FASTA files from UCSC goldenpath, say), e.g., for assembly hg19 and UCSC chromosome naming scheme:
$ for chr in `seq 1 22` X Y; do echo ${chr}; wget -qO- http://hgdownload.cse.ucsc.edu/goldenpath/hg19/chromosomes/chr${chr}.fa.gz | gunzip -c - > hg19.chr${chr}.fa; done

For each chromosome (omitting the for loop code for brevity), you can use the fasta-get-markov in the MEME suite to generate a 5th-order Markov model of kmer frequencies (why 5th order?):
$ fasta-get-markov -m 5 hg19.chr${chr}.fa > hg19.chr${chr}.5th_order_background_model.txt

Here is one way to download a non-redundant, vertebrate JASPAR database (MEME-formatted):
$ wget http://jaspar.genereg.net/download/CORE/JASPAR2020_CORE_vertebrates_non-redundant_pfms_meme.txt

To run a FIMO query with the 5th-order background model and JASPAR database, using a p-value threshold of 1e-4:
$ fimo --verbosity 1 --bgfile hg19.chr${chr}.5th_order_background_model.txt --thresh 1e-4 --text JASPAR2018_CORE_vertebrates_non-redundant_pfms_meme.txt hg19.chr${chr}.fa \
  | tail -n+2 \
  | awk -vOFS="\t" -vFS="\t" '{ print $3, $4, $5, $1":"$2">"$10, $8, $6 }' \
  | sort-bed - \
  | starch --omit-signature - \
  > hg19.chr${chr}.bed.starch

The output of a FIMO search is a text file which can be modified into a sorted BED file via sort-bed, an example of which is shown above.
To get a faster answer, each per-chromosome query can be put onto its own node on a computational cluster.
Using a p-value threshold of 1e-4 will result in very large files (when uncompressed). Using tools to compress BED files is useful here, as is using a more stringent threshold (1e-5, etc.).
Use of BEDOPS starch to compress BED files (as shown above, for instance) will give better compression performance than bgzip with per-chromosome random access and integration with BEDOPS bedops and bedmap CLI tools. The bgzip option offers integration with tabix and offers finer-grained random access, but it will use more disk space and doesn't integrate with BEDOPS tools.

This is a decent answer, so I am accepting it.  If you ever hear of these files being added to a database or of a database of this type of data please let me know :)

I can't presume to speak for groups like UCSC or Ensembl, but I'd guess that each analysis of this kind is short-lived, based on the fact that it really comes down to what motif models are used (which get updated), what statistical parameters get used (which are down to the experimenter), and the reference genomes involved. It's perhaps not something unique enough that you'd generate once in a blue moon, that it could be efficiently made available in a genome browser or via a biomart-like interface. It is basically something you make to meet your particular experimental needs, I think.

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


text:
DESeq2 complicated design - effect of replicated samples
I have RNAseq data from a relatively complicated experimental design with variables = genotype, treatment, time, and batch. I have 2 biological replicates for each genotype/condition, however unfortunately in the first iteration of the experiment one sample had poor RNA quality and is unusable. To get around this lack of a replicate, we repeated the experiment for the two biological replicates in the condition of the lost sample, and did a new library prep/sequencing from this repeated experiment.  The repeated experiment obviously lead to a batch effect, and although it isn't the main driver of variance, it is a significant factor.  I would like to model or remove the batch effect from the experiment using DESeq2 if possible, but after reading this similar post: can I model technical replicates in DESeq2? I am concerned about inflating the power of the test for those conditions that were replicated.
Here is a diagram of the experimental design to make it more easy to understand:

My "colData" for the DESeq2 object looks like this:

My previous thought was to simply collapse the genotype+time+treatment variables into a single combined factor "condition" and then model the batch effect using a design of "~batch + condition".  However this doesn't really address the concern of inflating the power of the tests due to technical replication.
Should I Just collapse the counts of these replicated conditions even though the batch affect will affect the combined count values? Or should I use a different tool than DESeq2 ?  Or is there a way for me to remove the batch effect prior to collapsing the counts for replicates?
What would  you guys recommend?

Were the libraries in the second batch simply resequenced or were they actually recreated (presumably from the same samples, though I wonder if this is cell culture so these represent different passages).

@DevonRyan The libraries in the second batch were actually recreated from the same samples --- this is indeed a cell culture experiment, so different passage number is a concern, however these repeated samples were within ~1 passage of the samples in the original experiment (they may have actually been from the same batch of cryopreserved aliquots).  So, the experiment was repeated for the two KO cell lines and then libraries were prepared using the same techniques/kits and sequenced on the same machine.

Different passage numbers are actually a good thing in this case, since it makes them not technical replicates. I would suggest proceeding with your current design matrix.

Good point.  It is likely that there were small variations in timing of samples that also would have affected the replication.  Are you aware of any biases that could arise from this type of experimental design? In the post I linked you and others mentioned that technical replicates could inflate the power of certain comparisons. I'm specifically concerned about comparing the # of DEG's between the wild-type and KO conditions for drug vs. mock treated conditions -- if power is inflated I may see more DEG's for the KO comparisons due to more replicates right?

I agree with @devon-ryan that the design matrix looks reasonable, i.e. that it's okay to treat the new samples as biological replicates of the original batch, and that the batch number should be included in that matrix (assuming DESeq2 doesn't complain about it) to filter out most of the batch effects.
DESeq2 tends to be conservative in its differential expression calculations. While it's true that more replicates in one sample will increase the number of differentially-expressed genes that are statistically significant (because the biological variation can be better modelled), I wouldn't expect this to change the overall picture.
However, whole-genome differential expression testing should not be used on its own to demonstrate biological relevance. I recommend as a first-pass sanity check that people compare any differential expression results to what is seen on a genome browser to identify any quirky things going on (e.g. expression across only a part of the gene). Ideally, results should be confirmed with targeted analyses on the same (or similar) biological samples, either at the transcript level (e.g. qPCR) or at the protein level (e.g. Western Blot). When aggregating evidence from multiple analyses, the specific way in which one of those analyses was done is less important.
Aside: It'd obviously be better [in the future] to make sure that your sample replicates are both robust to drop-out, and allow for outliers to be noticed, reducing the batch-effect problem. In practise this means a minimum of four replicates per informative condition, but because funders / researchers are well known for only choosing the minimum, I'd just state that six replicates per sample as what is needed to produce good results.

Thanks for your comments.  I agree that it is important to validate any important findings by other methods like qPCR/Western.  And yes, I learned an important lesson about experimental design here when I had a sample drop out.  Its tough though when you want to do a time-series experiment comparing drugs and genotype because the amount of samples is x10 for each additional replicate.  Either way, I agree that in the future I will go for at least 4 replicates.  Thanks!

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


text:
Publicly available, free, complete database for antibiotics names and classes?
This is a tough one I think: is there a publicly available, up-to-date, free, complete database for antibiotics names and classes?
I am specifically looking for information like, e.g., cefpirome (is a) cephalosporin.
I've looked in a couple of different places, but nothing really comes close. E.g.:

CARD comes closest, but is not free (companies need to pay a license) and not complete/up-to-date (e.g., cefpirome is missing)
@Karel mentioned KEGG in his answer, but this also needs a license (companies for sure, academics if they download data via FTP).
Novel Antiobiotics Database has no class ... and stops in 2003
The Antimicrobial Index is web only, and terms of use apparently prohibit any kind of web-scraping.
A suggestion by @Pierre on Twitter was to try Antimicrobial stubs on Wikipedia. Good idea, but the class is not stored in the ChemBox and must be free-text parsed. Or writing a script which evaluates the Wikipedia category tags.
I drew blank at different government agencies (EMA, FDA, EFSA, NIH, etc.). The closest I found was from the WHO a list of critically important antibiotics. But not complete, and in a PDF.

Any pointer appreciated.

The FDA's NDC may be a good resource. It is freely available and was updated: 6/22/2017. You'll have to do some grepping to get the info you need but I think it is mostly all there. Cefpirome is not in the db but many cephalasporins are. There is also an unfinished drug list which I did not check.

Good one, didn't know it. But unfortunately error prone. Some entries contain multiple ingredients, e.g., "Erythromycin and Benzoyl Peroxide", which fall into different categories and cannot be re-assigned unambiguously with the table data.

Can you define what classes you want? For example, are efflux pumps a class? Classification schemes are not universally; I suggested CARD on Twitter because it's a full ontology rather than just 2 columns

Well, e.g., I'd like to have the association that cefmenoxime is a cephalosporin (bonus points for 3rd, 4th or 5th generation, but not a must). Twitter didn't show me a notification of your response until I specifically looked for it, sorry. But as I wrote in examples: CARD is not free, companies need a license.

ChEBI
It should be possible to get this information from the ChEBI database, see the exported tables. You could download the ontologies (in OWL / OBO), parse them using some ad-hoc parser or using a dedicated library (e.g., Pronto), build a directed acyclic graph based on the is_a relationship, and extract the antibiotic subtree.
KEGG
You could also use KEGG: Kyoto Encyclopedia of Genes and Genomes, but a subscription is probably required. The exported drug table should be easy to parse (e.g., in Python).
For a quick verification of the completeness, you can download the table, filter out everything except antibacterial agents and check the number of records:
curl "ftp://ftp.genome.jp/pub/kegg/medicus/drug/drug" > keg_all.txt

cat keg_all.txt \
   | perl -pe 's/\n/BREAK/g' \
   | perl -pe 's@///BREAK@///\n@g' \
   | grep 'CLASS       Antibacterial' \
   | perl -pe 's/BREAK/\n/g' \
   > keg_antibacterial.txt

cat keg_antibacterial.txt | grep "///" | wc -l

There are 646 records. However, some of them are almost identical (e.g., only 1 atom difference).
If you add | grep -i "antibiotic" \ to the pipeline, then you get only 445 records.

Excellent resource, but does not pass the "freely available" filter: KEGG requires a yearly license. I will update my question to include that.

Would ChEBI satisfy your definition of "freely available"?

ChEBI would, indeed! And it contains approximately what I need. OBO or OWL are not necessarily formats from which it is easy to extract the wanted information as one needs to build and traverse graphs to piece together required data, but it is doable. Can you please amend your answer to include ChEBI before KEGG? I'll accept the answer then.


text:
Have DNA motifs 6-12bp long, trying to get conservation scores
I have about 200 short nucleotide motifs (6-12 bp in length) from the human genome, and I'm trying to see how conserved they are across vertebrates.
I was thinking that I'd need to make a bed file for each motif that lists all of its occurrences in the human genome. From there, I could map the beds to a bigwig files of PhastCons scores (essentially doing the reverse of what the PhastCons software was designed to do). Does that sound like the best approach?
I'm getting stuck at the step of going from motifs to bed files. I've tried using BLAST to find all occurrences of motifs, but their short length is causing issues.
I've tried messing with the e-value threshold, word size, and filter parameters, but I still don't get any hits.
Is there a work-around for this issue, or should I just rethink my entire approach?

I did something similar, but with a list of well defined n-mers instead of "motifs", all with the same value of n. In case this can be of any help, the code I used to make the bed file is the following: https://bitbucket.org/blaiseli/conservation_mirna_targets/src/20b7a0894b957e69912f637ef2d0493779036029/trie.py?at=master&fileviewer=file-view-default
The code might not be super memory efficient.

In case you have only ACGT in your motifs
The short motifs make it sound as if you are in the business of looking for a kmer counter. You can choose to use existing software or build your own.

Using existing software might be your easiest path. An older post from 2014 will probably give you a first idea what's out there: http://homolog.us/blogs/blog/2014/04/07/kmer-counting-a-2014-recap/ . Note that a couple of algorithms mentioned there got successors, so it is worthwhile digging a bit around. The small kmer size will make most of them usable for your needs.
As the maximum size of your kmers is comparatively small (12 nt need 24 bits, i.e., max 16.7 million entries in your kmer table), you should be able to easily roll your own kmer counting in about any language you like and on about any of nowadays computer. The pseudocode section on the Wikipedia entry for kmers will give you first pointers for that. Might be a bit more work, but maybe more flexible depending on your needs.

In case you have IUPAC bases (N, W, etc.) in your motifs
I don't know any any pre-existing software doing what you need. I could imagine that the short motifs make using regular expressions doable for this kind of search, but I may be wrong. Testing this should be easy in a simple script as all major programming languages have modules or libraries for REs. Even if it should take a couple of hours to run on your data set, that would be good enough for a one-off calculation.

To scan motifs in a genome (or database) I would use FIMO which will give you the exact locations of these motifs in your genome.
Once you have the locations, you can use a phastCons bigiwig from UCSC to calculate the basewise conservation scores. However, please remember that phastCons scores are smoothed across windows and it might not be the best metric if you are trying to compare the conservation levels at your motif matching sites as compared to the sequences flanking them.
I wrote a package a while back to do this, including doing de-novo motif discovery. However, it might be an overkill for your use case.

I've used FIMO before actually, but for comparing motifs against PWMs (from JASPAR). Might be a dumb question, but is it pretty easy to compare against a genome rather than a database of PWM's?

@EricBrenner You can simply upload a sequence file in the web version.

Where transcription factor binding sites are concerned, in our papers we have used per-base phyloP conservation data instead of smoothed phastCons scores.
We use BEDOPS bedmap to map scores over multiple, padded binding sites (BED-formatted) for a given motif model. The phyloP scores are WIG files obtained from UCSC goldenpath and converted to BED via wig2bed.
The resulting matrix of binding site intervals and their per-base scores can be turned into a ranked heatmap or aggregated to determine average per-base conservation for a motif model.
While TF binding sites innately show high information content and therefore high conservation, matrices can be further sorted by score-maps of ChIP-seq or DNaseI-seq tag density.
This can help filter for interesting low-noise, high-signal patterns and is useful when the matrix is windowed and the window includes low-information regions with a few high-information residues outside of the main motif — as is the case with CTCF, for instance.


text:
Detecting portions of human proteins with high degree of microbial similarity
I'm a newcomer to the world of bioinformatics, and in need of help solving a problem.
My goal is to take a list of human proteins, and identify segments (13-17aa in length) with a high degree of similarity to microbial sequences. Ideally, I would like to start with list of FASTA sequences, and have an easy way to generate an output of the corresponding high similarity segments of each protein.
Are there existing tools or software that I should be aware of that will make my life easier?
Thanks in advance.

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

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

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

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

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

@Chris_Rands True enough. :-)

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

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

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

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

BLASTp: Simple protein-protein blast. It will identify homologous proteins based on sequence similarity. Whether or not that also implies functional homology is not that simple and will depend on each case you investigate.
As above, go to the NCBI protein blast page, but this time use the defaults.
tBLASTn: this is a tool that takes protein sequences as input and compares them to a database of DNA which is dynamically translated in all 6 possible reading frames. Very good for finding homologous sequences when you don't have well annotated protein information for the target species. It has the benefit of being more sensitive and able to find more distant homologies than basic nucleotide BLASTn and the go-to approach when your target species is distant and not well annotated.
NCBI's tBLASTn page.

All of these can be run online through the NCBI's BLAST page. If you want to investigate hundreds of proteins, I suggest you install blast locally. You can then either download the relevant target sequences from NCBI and rebuild the blast database locally (if so, I suggest you ask a new question about how to do that) or, use NCBI's remote blast client which lets you use a locally stored query file and will run blast on the NCBI's servers.
Now, these programs will return what are known as High Scoring Pairs (HSPs), the regions of your query sequence(s) that align well to the target. There are various options you can play with to improve the sensitivity or the specificity but a discussion of those would require far more details about what you're doing and would be best in a new question as well.
Once you have your HSPs, you can relatively easily parse them to select regions with a given range of sequence similarity values and of a specific length. Once again, that would be better discussed in a separate question once you have your results and can show an example.

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

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


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

One big downside of KEGG is the licensing issue. One big advantage of Reactome are various crosslinks to other databases and data.
ad 1, This depends on which pathway, they are both primary databases. Sometimes other databases that for instance combine data of primary databases have better annotation of pathways (there is an example in the review paper bellow)
ad 3, There is very extensive relatively new (2015) review on this topic focused on human pathways: Comparison of human cell signaling pathway databases—evolution, drawbacks and challenges. However I could not find there which one is more complete ...

For you the main point would be whether an enrichment analysis is going to give you an informative answer. That's what will make a particular database better. And there's all sorts of subjective decisions made in their construction as to what what a pathway is, what to include, where to draw boundaries, etc. so different databases will give different answers and it will not be clear which is more correct.
So browse over the references / pointers people have given, select a service and use it and no others. Don't jump around databases until you get an answer you like, that's just fishing.

+1 for the last sentence, is so easy to say I will check with X once one have a "good" results and if both agree then go on.

One big advantage of Reactome, in my opinion, is its visualization using the web interface.
Many pathways (in Reactome and KEGG) consist of genes / proteins that are up- and down-regulated through the respective pathway. If you do a simple overrepresentation analysis this is not taken into consideration. Therefore, you might end up seeing a pathway as "overexpressed" although only the down-regulated genes were observed more frequently.
In Reactome, you can zoom in on the different pathways very conveniently and then pick up these inconsistencies. I haven't really found a public database and tool that can take these different regulations into consideration. Therefore, you'll probably always will need some manual investigation of your data. In my opinion, this is easier with Reactome than with KEGG.


text:
Deciding which samples go in which batch
I have 370 samples to sequence, we probably will end up using only 96 samples per run (due to the barcode with the primers we'll use). This means running 4 batches. To minimize the batch effect I need to decide which samples are sequenced in which run.
I have several variables from these samples, some of them are categoric (sex, disease status, treatment stage, ...) and some of them are numeric (age, clinical score...).
To minimize the batch effect I wouldn't group all samples from the same category in one run, and I would distribute as uniformly as possible the samples according to the numeric variables. However, I could find many papers about adjusting, taking into account, but I couldn't find a paper about software to help designing the runs.
Is there any software to help designing each batch?
If there are several I would like to point out that I would define best as the one with easier user interface

I found the OSAT package which does almost what I want.

However it is specific for Illumina and for a limited number of plates: "plates that hold 2, 4, 8 IlluminaBeadChip chips and have 24, 48, 96 wells, respectively."
It doesn't take in consideration the numeric values:
library(OSAT)
inPath <- system.file('extdata', package='OSAT')
pheno <- read.table(file.path(inPath, 'samples.txt'),
                    header=TRUE, sep="\t", colClasses="factor")
pheno$AgeGrp <- runif(576, 0, 100)
gs <- setup.sample(pheno, optimal=c("SampleType", "Race", "AgeGrp"))
gc <- setup.container(IlluminaBeadChip96Plate, 6, batch='plates')
gSetup <- create.optimized.setup(sample=gs, container=gc, nSim=1000)

The minDiff package reported in Stats.SE.

It handles both numeric and categorical data, but it can only optimize for two nominal criteria.

I wrote experDesign which does everything I wanted

It can handle both numeric and categorical data for an arbitrary number of variables. It is not restricted to an even number of samples per batch.

I've had to go through this pain, and found the best way to do it was to simulate the model fit to start assessing the balance of covariates. I've got two notebooks available on Github, which should help you deal with this issue, although I've not thought about how to accurately assess continuous variables.
While it's not wrapped in a package, you can use these notebooks as a template to fit to your problem. I may wrap this in a package at some point, but I haven't had the chance to do so yet.
This Notebook will allow you to set up all your variables, max size of a batch (may vary due to multiplexing, library set up, plate size, etc). Bin them into appropriate batches, then create a design matrix, model fit, and test.
Apologies in advance that it's a bit of a mess, but if anything isn't clear, open an issue on the repo.

Interesting, thanks for sharing! However I don't understand why do you perform a test? If I understood the script it is for a mock data, so what is the purpose of this? To estimate what would be a completely random output?

That's a fair point, and you're correct, there's no reason for it. Generally, it was just for completeness so that nothing was widely unexpected when performing the test.

After trying it. I realized that this assumes even number of sample per each category (except for batch), which is not my case

Yeah, unfortunately you'll need to modify it to deal with uneven cases.


text:
Sequencing center returns data in several files
We send some samples to sequence and we got several (fastq.gz) files for each sample. The files are distributed at two or three folders with different dates (more than a week apart). The dates of the folders where the data is in:
2019-08-12
2019-08-21
2019-08-21
2019-09-04
2019-09-14

I asked the facility about why there are the same file in different folders, and they told us that:

we have split sequencing of that project on different flowcells, so all
  files are correct and you have to merge the individual (duplicate) files
  to get the requested total read numbers.

How can different flow cells be from the same experiment and be run with more than a week of difference? (I thought that all the flow cells are run at the same time in each sequencing)
Can be there some batch effect? How can I check if there is batch effect if we merge the files before trimming, mapping and counting?

Usually, a single sequencing library is created from each sample. Assuming you are using barcoded indexes to identify multiple samples, the index is attached to each library, and the libraries are pooled to make a library pool.  A library pool can be frozen, and is very stable.
A pool can be sequenced multiple times by pipetting a portion of the pool and supplying it to the sequencing process.  Each time, a sample of the molecules in the pool will be sequenced.
Generally, it is true that sequencing runs are very repeatable, and that you will not suffer any batch effects.
If you want to check for batch effects you would need to analyse all of your sequence files separately, and check this first, or in parallel, to your main analysis.
One thing you might want to check is whether the same library pool was sequenced multiple times, or whether separate pools or libraries were made for each sequencing run.  In the second case, there may be more likely to be batch effects.

Oh, thanks. I think they submitted a portion of the library pool because the index are the same. However, why they didn't add more sample from the pool from the beginning?

Regarding your sample pooling question, probably this is not happening in your case but here is our workflow, which generates multiple output files per sample like you get: We use 10x and the most expensive part is the sequencing run. As single cells are delicate and not all "droplet formation" steps go smoothly, first we do "shallow sequencing" on NextSeq to see if there are enough number of good quality cells in a given sample (after mapping, read counting, single cell QC) and then continue with deeper sequencing (one or more runs) with NovaSeq until we reach saturation with number of reads.

The number of sequences has a limit per sequencing run, even if you add more of the pool to the run - it is limited by the flow cell.  For this reason, it is normal to only sequence a part of the pool in any one run - enough to fill the flow cell - and if you need more reads, then carry out more runs with additional flow cells.

Another alternative would be to allocate some samples to one run, some samples to another run, and to fill up the flowcells that way.  However, this would mean that, if there are any batch effects it would be much more difficult to separate them from sample effects.  At least if all samples are included on every flow cell, the batch effects will not affect any between-samples comparison.

@Jonathan Well, in our case not all samples where included in all the runs, that's partially why I ask about the batch effect.


text:
How to select the most representative pathways from a gene enrichment analysis?
I have perform an enrichment analysis to a cluster of genes. The output is a list of pathways and their p-value (the pathways are selected because p-value < 0.05). The list is still quite long, so I want to reduce it. For that purpose I have a calculated the Dice coefficient of the pathways in a matrix $p$x$p$ where $p$ is the number of pathways in the list. I want both the ones that are more different (they overlap less, their Dice coefficient is lower) and the pathways more representative of the most similar pathways (So if a there is a group of 5 pathways that overlap over 0.8 take just one).
How can I select the most representatives pathways?
There is a similar tool for GO but it relays on discarding not significant GO, while here all the initial pathways are already significant.
If I do a clustering of the pathways using the Dice coefficient matrix I don't know where (or how) to cut.

I tried using the height to select the pathways. But I am unsure of the interpretation of height.
Some other tools I have seen use a multidimensional scaling plot, but I am not sure if performing it and cutting at certain point of the first dimension would help.

Are you using p-value < 0.05 or p-value < 0.05 / (# of terms tested)? You'll want to use the former to avoid false positives due to multiple testing

I'm guessing you'll have to settle with a "reasonably good" heuristic solution, since there are a few different variables in the problem you're trying to solve (number of pathways in group, desired overlap/distance between pathways, etc.) I can imagine many possible answers.

@CloudyGloudy correcting for multiple testing is already done, sorry I didn't mention it previously. Yes, I was playing with the idea of selecting an overlap/distance of 0.5 and keep those above that, but that would leave outside the pathways that are completely different from the others, but maybe I could select those above 0.75 and those below 0.25

One way to group similar significant pathways, is to quantify how many genes overlap between the pathways, and then use this in clustering (heatmap). I have made a tool in R which calculates the overlap index between GO terms and subsequently clusters them in a heatmap. Overlap index is the fraction of genes that overlap (number between 0-1). Also Pearson correlation can be used for clustering instead of the overlap index. My package (gogadget) works only with  goseq analysis, but you can use goseq also for reactome or kegg data.
I have used this overlap clustering approach for different data sets now, and usually we can reduce 200-300 GO terms into 10-20 functional groups.

Which overlap index does it use? There are several overlap indices. I would advice against using this approach in GO, there are specific methods to measure how similar two GO terms are. See GOSemSim.

The overlap index is defined by the number of overlapping genes divided by the number of genes in the smaller of the two gene sets. It is described in Bioconductor Case studies, chapter 13.3. This book is written by big names in bioinformatics (such as R. Gentleman and W. Huber).

Thanks for the new overlap index, I didn't knew it. BTW the way these similarity between genes is calculated is using other overlap indexes so this won't help.

Okay, but my point is not about which index to use but the clustering approach (same as gringer is suggesting). If you look in the userguide of my package on page 26-28, you'll see a heatmap. Usually I get good results with Ward.D and Euclidean. The tree can be cut in R, but you first might want to see and evaluate which gene sets are clustered together in which branches before you can find the right cutoff.

Sorry, my last comment was thinking in another question. Yes, this approach could work

This sounds like something that might be amenable to a clustered heatmap plot, or a correlation matrix plot, or something similar. Have you looked at a correlation matrix of the dice coefficient matrix (or maybe just a heatmap plot of that matrix without the correlation matrix)?
The corrplot package looks like it might be useful, in particular the hclust / drawing rectangles presentation.
I can't vouch for this package though; it's just something I found by a search for "R plot correlation matrix".

With a heatmap or a correlation plot I could observe the similarities between the pathways, as I currently do with the dendrogram and the MDS plot. However the question is how to select those more representative pathways. Thanks for your corrections and comments btw.

If you're happy with a more confident ranking of the most representative gene sets, rather than necessarily cutting down the list, you might try EGSEA. It uses an ensemble approach to give a ranking of the most relevant gene sets, and also produces an interactive HTML output with statistics, heatmaps, pathway maps, summary plots and GO graphs which allows you to examine the output at varying levels of granularity.
You can read the paper on bioRxiv or download the package from Bioconductor.

The ranking is already done, through the p-value (one could argue that a more elaborated test like the proposed by EGSEA would be better) but this doesn't seem to answer how to select relevant pathways from the output of EGSEA or other software/tools/methods

You shouldn't be ranking by p-value. The p-value is only an indication of whether the observed change is statistically significant, not an indication of the magnitude of the observed change.


text:
What does it mean to over sample?
One of my studies (A) involve sequencing the microbiome. After selecting the variable region and the primers for the targeted region of the 16S sequence. The samples were sent to a platform.
For another study (B) we used the same strategy (same primers, region, machine...) but another platform. However we obtained from this platform a lot more sequences:
In study A the mean sequences per sample are 10k, while in the study B the mean sequences per sample are 80k
I consulted with others and it seems that the study A we under-sampled and in study B we over-sampled.
I can understand that by under-sampling we didn't obtain enough reads to be sure about the microbiome present in the sample. But what does it mean to over-sampling?
That expert I consulted suggested normalize the data to fewer sequences, discarding information we already have. I can understand that the sequencing has errors and can create chimeras but the resulting consensus sequences of the OTUs have found a match with databases of microorganisms. What or how should I take into account that I over-sampled?

As you sequence a library you are sampling from the pool of all possible oligonucleotides found in it. This is a stochastic process, meaning that you are more likely to sequence a piece of DNA that is more prevalent.
As you increase the number of reads that you sequence for a library, you decrease the chances of seeing a new read that you have not yet seen before. If you plot the number of unique reads per total reads sequenced, it will look something like the cumulative density function of the exponential distribution, where you have diminishing returns on unique sequences as you continue to sequence.
If you sequence a library and quickly stop seeing unique reads, we say that this library has low complexity. Therefore, you have "oversampled" it if you sequence too deeply. I am guessing that your metagenomics library actually has a lot fewer species/unique sequences than they anticipated, and that you sequenced deeply anticipating more.
Edit: In response to your comments, the bigger problem in 16S analyses is not sequencing too deeply, it is over-amplifying your samples such that they no longer represent what is actually in your gut microbiome samples. In this case if you used read normalization you would not accurately represent what was in the library that you sequenced since the amplified molecules (16S) are so small that you only have len(16S molecule)-(read len) + 1 possible reads in each direction. You might consider removing PCR duplicates with GATK but this will still have the same effect of not representing what was in the sequencing pool. It is very possible that you have many PCR duplicates after 16S amplification. I generally don't trust amplicon sequencing for quantitive analyses, and don't really trust for low coverage "is-this-species-there-or-not" either, especially when dealing with the 'dark matter' of uncharacterized species. You might be well off just asking "is this species present in high coverage at this time point?"

Thanks for your answer. The problem raise when the same samples in both platforms yield different results (but take into account that I am using OTUs not unique reads). I can upload an image if needed. Anyway, how to take oversampling into account (If it is something to take care of)? I am worried that I can't then compare the results between the studies

BTW, wouldn't that be considered then deep-sequencing?

Yes, this is "deep-sequencing", however I really think that term is overused and overhyped. As you can see in your case, deep sequencing actually causes its own problems with normalization on platforms like Illumina or with PCR-based libraries. :)
Whether or not the sample coverage requires special handling depends on the downstream analysis that you would like to perform. Is this whole-genome shotgun metagenomics, something like 16S amplicons, or something else?

As described in the question it is 16S done in MiSeq. I want to perform integration with RNA-seq data and analyze them on its own via comparing different times, response, localization (it is gut sequencing), sex ...

I don't think there's anything inherently wrong with sequencing to high depth (though naive approaches to dealing with sequencing errors such as singleton removal may fail if the depth is so high that the same errors may be encountered repeatedly by chance). The bigger issue is that they are over sampled relative to the undersampled microbiomes. This does complicate comparisons, particularly for things like OTU richness. The obvious solution is to down-sample the oversampled microbiomes so you have the same number of reads for each. This is probably unnecessary if your analyzing the oversampled microbiomes on their own.

I think like this answer specially because the tool I use to transform the fastq files to OTUs is specially designed to remove singletons and chimeras.


text:
Linear models of complex diseases
A popular framework to analyze differences between groups, either experiments or diseases, in transcriptomics is using linear models (limma is a popular choice).
For instance we have a disease D with three stages as defined by clinicians, A, B and C. 10 samples each stage and the healthy H to compare with is RNA-sequenced. A typical linear model would be to observe the three stages~A+B+C independently. The data of each stage is not from the same person. (but for the question assume it isn't)
My understanding is that such a model would not take into account that stage C appears only on 30% of patients in stage B. And that a healthy patient upon external factors can jump to stage B.
If we want to find the role of a gene in the disease we should include somehow this information in the model. Which makes me think about mixing linear models and hidden Markov chains.
How can such a disease be described in terms of linear models with such data and information?

Of course, your model is not taking into account that observations (i.e 30% of patients from B evolve to C). Here, limma would look for genes with most significant differences among states based on F-statistic. If you are more interested in disease progression, you should read more about time-course patterns (states follow a progression, so positions are not transposable).

But time-course patterns wouldn't handle the jump from healthy to B. A quick look for disease time course patterns show the top hit with an article from 2001. Could you point me to a review or an introductory title

Could you clarify (probably in the question) if stages A, B & C come from the same individuals?

@IakovDavydov Done, the stages are from the disease, the data might come from the same individuals or not. But we can assume the data from each stage doesn't come from the same individuals

Why do you think it's useful to include the progression information in the model ("we should include somehow this information in the model"). You're only able to discern when/if the gene is changing and no statistical model will change that.

What about saying that if you progressed to B from C, then you are positive for both those variables, rather than that you just have B?

@DevonRyan My rationale behind is the following: imagine that a gene between B and C is found to change the expression. Is the change in the gene expression a consequence or a reason of changing the stage?

@IanSudbery I thought about that possibility but clinicians don't think that way. But I would actually like to test this model with the dataset I have at hands.

@Llopis You're not going to get causality from this sort of data.

I don't understand your description of the linear model as ~A+B+C. If these are indeed stages of a disease - then the model would rightly just be ~Disease. And the levels of the data column or vector would be the factors (H,A,B,C) - where H is reset to the reference level.

@StephenHenderson Doing that way the model I couldn't compare between stages.

@Llopis I believe you are confused and have that the wrong way round. I would refer you to section 9.3 page 42 of the limma manual on dealing with several groups. Note the actual congrats of interest are specified separately.. https://www.bioconductor.org/packages/devel/bioc/vignettes/limma/inst/doc/usersguide.pdf

@StephenHenderson AFAIK, section 9.3 use ~0+f (where f is the factor with the groups) which is the same as ~A+B+C (that's why later the columns are renamed to RNA1, RNA2, RNA3 in the tutorial, if It had one Disease/Group column how could it rename three columns?)

@Llopis that's the design matrix - 3 columns for 3 groups in the single covariate (f). ~A+B+C is like a model with 3 covariates Disease+Sex+Age... see for instance section 9.42 on the next page 43.

@Stephen The same design matrix can be build of different ways. I think we are talking about the same model but we build the design matrix differently. But that would be another question IMHO

Like others, I am confused by your model. You say A, B, C are not all observed in a given sample - some individuals never make it to C. Whether C is binary or a quantitative trait, if it never happens for some individuals, then no data is observed that would enter the model.  These sample’s estimates would be functions of their observations of A, B.  If you wish to more strongly impact the model estimation process, you can look at like a general likelihood, where you add priors to model coefficients, as well as encode other prior information.

There are two potential sources of bias in this design.

We cannot distinguish correlation from causation.

Imagine two cases. In the first, the disease progression is inducing immune response. Later stages will be associated with the higher gene expression levels. In the second scenario, the disease is caused by overexpression of a gene. Later stages will be also associated with the higher expression.
This is typical for observational studies. But I just want to mention that a special care should be taken during interpretation of the results.

If we are not following our individuals, we cannot distinguish correlation and avoidance.

Let's say the disease is lethal in certain cases; the survival is negatively correlated with the disease stage. Now imagine there is a gene, which causes severe symptoms when highly expressed. On the later stages you will only observe those patients in which the gene was not highly expressed. From that you would conclude that gene expression is decreasing with the disease progression. In reality this gene is very important and causal, you just do not have patients which are alive to see this.
This is similar to Wald's studies of aircrafts (Survival bias).

Researchers from the Center for Naval Analyses had conducted a study
  of the damage done to aircraft that had returned from missions, and
  had recommended that armor be added to the areas that showed the most
  damage.
  Wald proposed that the Navy instead reinforce the areas where the
  returning aircraft were unscathed, since those were the areas that, if
  hit, would cause the plane to be lost.

I think the second point is crucial, and can and will lead to false conclusions. I do not see immediately how to avoid this bias if you study different individuals.
Therefore I suggest that the same individuals are followed for a long time.
There are different approaches you can use later. For example you can have two-stage procedure:

Identify genes which are differentially expressed between healthy (H) and sick (A, B or C).
Build a linear model of disease stage stage ~ gene1 + gene2 + ..., using genes identified at step 1.
Similarly build a linear model of survival as a function gene expression.
It is possible you can use logistic regression to infer probability of transition from stage B to stage C.

This is not the only possible approach of course. As you suggested, Markov model is also applicable.
While following the same individuals it is possible to use continuous-time Markov model for such a study. In this case discrete or continuous measurements are recorded at certain time points, and the model parameters are inferred using maximum likelihood.
I'm not an expert in this field, but I think this paper describing the msm package for R will be useful. It also supports hidden Markov models.

You explained better than myself why this model is not good to find causation :). However in your proposed solutions I have two problems. 1) In certain diseases is not easy to follow the same individuals for a long time (for instance you might have to operate to extract the data, so he/she can't be operated every year). 2) Your two models +survival analysis doesn't take into account that stage C appears on only 30% of patients in stage B.

I agree that it is not always possible, although I do not immediately see a way how to overcome this bias. 2) I expanded question on this (see point #4 and msm).

Ok, I guess there isn't a way to mix the data in linear models :(


text:
PCA of methylation-values normalized by coverage
I have started analyzing methylation (EM-seq) data for the first time ~0.8M positions.
In this data set I have 27 samples of 20 patients. I want to perform a PCA of the dataset to check for possible batch effects or problems.
Some publications and guides use the raw M-values to plot it without taking into account the coverage of each position. If I do that, I can detect a previously unknown source of variation, which is already great.
Usually there is also the information of the coverage of each position (not the number of reads per sample). Example of the data, which some packages store in the same object:
Cov:
      A      B      C      D      E
     19     10     18     14     20
     13      8     20     12     12
     12      6      9     16      7
     10      7      9     16      8
      8      5      6     12      8

M-values:
      A      B      C      D      E
   9.88   0.00   0.00   0.00   0.00
   7.93   0.00  10.00   3.00   0.00
   9.00   1.98   5.94  13.92   4.97
   9.00   2.94   3.96   8.96   4.96
   2.96   4.00   1.98   4.92   4.96

As you can see the coverage of each position is not the same for each sample.
However, I haven't found any specific method for calculating a PCA taking into account the coverage: there doesn't seem to be a method that normalizes the M-values by the coverage of each position.
I checked the distribution of coverage and 91% is covered by less than 14 reads, but some positions are covered by more than 300 reads.
What is the best way to normalize the position by coverage and plot the PCA ?

This is because of the statistical issue of Eigen value calculations - transformation are generally not welcome. In the changing age of unsupervised learning the concept of PCA is changing admittedly. You can superimpose coverage post PCA to assess with it is affecting the output, i.e. if it starts forming clusters there's a problem and your approach is justified.

I have experience in the analysis of EM-seq sequencing datasets but I honestly never really used M-values instead of classic fractions of methylation. One common approach is too remove highly covered positions (outliers) as it is often the results of sequencing or mapping errors. If you have a average coverage of ~ 14X, I think a threshold at a 100X is good. On the detection of methylation variation accross sample conditions, you won't need to deal with coverage variations a thes DSS soft (for DMC and DMR detection) works with raw counts. Not sure if it helps!

@PaulEndymion I also used the beta-values for a PCA (it revealed two other sample outliers). Interesting point about the mapping and the errors, I haven't actually mapped it so I don't know how well it was done. Are you suggesting that I don't need to perform PCA because they are already taken into account? I will look if removing those highly covered positions changes anything. Thanks!

Yes mapping errors on high repetition/GC% regions is a common bias, especially with short reads. These regions have the potential to attract reads with sequencing errors. We would do a PCA after the preprocessing of raw counts from "stranded_CpG_reports" files. Preprocessing included normalization of counts and removing of known SNP falling on "G" or "C" positions in "CG" motifs as well as higly covered positions. However as @M__ said if you found clusters of samples with PCA you can't explain with your samples experiment conditions, you might have a problem with these samples.

Just to be clear on the PCA. You can perform PCA but more as part of an exploratory/descriptive analysis, I don't remember reading about differential methylation studies using only PCA for DMR or DMC detection (even though DNA methylation studies are tedious so there is room for new approaches).


text:
How can I read multiple different regions from a BAM file in R?
I’m trying to process BAM records (of a coordinate-sorted, indexed BAM file) successively in fixed-size, non-overlapping, sliding genome coordinate windows. Unfortunately I am unable to read records via scanBam in different windows (which parameter): the first call to scanBam succeeds but subsequent calls return zero-length results.
By contrast, if I don’t specify a coordinate window and instead use yieldSize when opening the BAM file, I can successfully read records in a loop:
library(Rsamtools)
bamfile = system.file('extdata', 'ex1.bam', package = 'Rsamtools')
process = function (reads) message(length(reads$seq))

bam = open(BamFile(bamfile, yieldSize = 1000L))
param = ScanBamParam(what = 'seq')

repeat {
    reads = scanBam(bam, param = param)[[1L]]
    if (length(reads$seq) == 0L) break
    process(reads)
}

… where process performs the read processing (unimportant for the example).
Now I’m trying to use the following code to work in genomic coordinate windows instead, and it fails as explained above:
windowsize = 1000L

bam = open(BamFile(bamfile))
chr = seqlevels(bam)[1L] # Only do one chromosome for now
chrlengths = seqlengths(bam)

result = lapply(seq(1L, chrlengths[chr] - windowsize, by = windowsize), function (pos) {
    which = setNames(RangesList(IRanges(pos, pos + windowsize)), chr)
    param = ScanBamParam(what = 'seq', which = which)
    reads = scanBam(bam, param = param)[[1L]]
    process(reads)
})

This only returns data from the first iteration, not any subsequent iterations.
Conversely, it works if I close and re-open the BAM file in each iteration (i.e. before each scanBam call) — but that’s undesirable: it carries a substantial overhead (running on an empty BAM file with 16 chromosomes in its header takes several hours).
Since the process function is nontrivial and written in R, I’m stuck with R for now. And there’s no reason why this shouldn’t work in R.
Am I overlooking something? Is this a bug in {Rsamtools}? Is it “by design”?

@gringer I’ve noticed that you’ve added the tag into the question title — why? On Stack Overflow there’s a broad consensus against redundant mentions of the tag in the question title, and in fact I frequently (several times a day, on average) edit to remove such mentions from titles on Stack Overflow (most useless taggings are more egregious, but removing “… in X” is a fair share of the edits I make).

Following the same logic, why not remove BAM from the title? I don't think that "in R" is redundant here. In this case I considered that there were multiple valid Bioinformatics questions that could begin with "How can I read multiple different regions from a BAM file", even with the "R" tag. In this case it seems like you specifically want an R solution, which is why I added it in.

Are you working under linux? and do you absolutely have to do it in R?

@user324810 I mean, this question is three years old. And at any rate, regardless of whether somebody “absolutely” has to do this in R, it should work in R. But for posterity: yes, I absolutely had to do this in R since the process placeholder in my case was a piece of nontrivial logic developed in R over the span of months (in fact, the final implementation was in C++ anyway, the R code was my research prototype for faster iteration). I ended up processing fixed-sized chunks rather than scanning over genomic coordinate windows.


text:
Can I model technical replicates in DESeq2?
I’d normally use collapseReplicates (or do the collapsing upstream) to handle technical replicates.
However, in my current RNA-seq experimental design, samples were sequenced twice using different library preparation protocols, leading to marked differences in the resulting count estimates (in fact, the choice of different protocol explains most of the variance in the bulk data according to PCA). I would therefore like to model these differences by adding them as a covariate to the DESeq2 model.
I thought I could just add another factor to the design formula, equivalent to the “pasilla” example in the DESeq2 vignette, which adds the factor type for the library type (single-end vs paired-end). However, the pasilla vignette states that these different libraries are actually independent biological replicates, not technical replicates as I had previously assumed.
I’m now concerned that having technical replicates in the design might artificially inflate the power. Can the design be salvaged?
It’s worth noting that we also have biological replicates, and a full rank design; i.e. every combination of treatment and library prep, with biological replicates.

If they're truly technical replicates, then there's no way to model them using DESeq2*, as you've alluded to with the collapseReplicates function. DESeq2/ Mike Love's general recommendation with collapseReplicates is to just add the reads together for technical replicates.
If you want to model them instead of collapsing them down, you can voom transform your data and follow an example similar to section 11.3 in the Limma users guide.
(Code example from 11.3):
Design
| FileName  | Cy3   | Cy5  |
| --------- |:-----:| ----:|
| File1     | wt1   | mu1  |
| File2     | wt1   | mu1  |
| File3     | wt2   | mu2  |
| File4     | wt2   | mu2  |

Example
> biolrep <- c(1, 1, 2, 2)
> corfit <- duplicateCorrelation(MA, ndups = 1, block = biolrep)
> fit <- lmFit(MA, block = biolrep, cor = corfit$consensus)
> fit <- eBayes(fit)
> topTable(fit, adjust = "BH")

*Edit: As @Devon Ryan mentions, you can design a paired model in DESeq2, but beware of making this full rank.

This answer makes me think that we need table layout support on this website: tables of features/counts/… will be extremely common here.

@KonradRudolph - Agreed, It was a shame when I found out that markdown tables weren't supported!

Practically speaking, there's no way to include the technical replicates in that design (in DESeq2 at least). Your concern regarding inflating the power is exactly correct and the only way to combat that would be to add a pairing factor like one might do with case-control or tumor-normal studies. That is, something like:
  group libraryPrep sample
1    WT           A      1
2    WT           B      1
3   MUT           A      2
4   MUT           B      2
5    WT           A      3
6    WT           B      3
7   MUT           A      4
8   MUT           B      4

Here sample pairs the technical replicates together. However this ends up either being rank deficient or in the end not more informative.


text:
Are there any RepBase alternatives for genome-wide repeat element annotations?
I’m using the RepBase libraries in conjunction with RepeatMasker to get genome-wide repeat element annotations, in particular for transposable elements.
This works well enough, and seems to be the de facto standard in the field.
However, there are two issues with the use of RepBase, which is why I (and others) have been looking for alternatives (so far without success):

RepBase isn’t open data. Their academic license agreement includes a clause that explicitly forbids dissemination of data derived from RepBase. It’s unclear to what extent this is binding/enforceable, but it effectively prevents publishing at least some of the data I’m using and generating. This is unacceptable for open science.

Subordinate to this, the subscription model of RepBase also makes it impossible to integrate RepBase into fully automated pipelines, because user interaction is required to subscribe to RepBase, and to provide the login credentials.

RepBase is heavily manually curated. This is both good and bad. Good, because manual curation of sequence data is often the most reliable form of curation. On the flip side, manual curation is inherently biased; and worse, it’s hard to quantify this bias — this is acknowledged by the RepBase maintainers.

Were you only asking about defined repeat libraries? I interpreted it slightly more broadly as about tools used to build the libraries also (which becomes relevant when genomes from new taxa are sequenced)

@Chris_Rands Both (libraries and tools). Your answer is spot-on.

The goal is to build the annotated library of repetitions, or to mask repetitive parts of a genome?

@KamilSJaron I’m working with TEs, so I need the annotated library, not (merely) a repeat masked sequence.

Ouch, but for TEs just a subset of repetitive regions, there are specialised tools to annotate them (like DNApipeTE and REPET). Maybe you could specify it in the question.

@KamilSJaron Nice, and this might be worth an answer. I’ll also update teh question. That said, I am also asking for repetitive elements beyond TEs.

Also are you sure that RepBase is against publications that derive information based on their data? Could it be that they just don't want you sharing the raw data and files they provide. Since they have a good number of citations in different fields: https://scholar.google.de/scholar?um=1&ie=UTF-8&lr&cites=11574259945967474319

@story They literally say so in the academic user agreement that I link to. Here’s the relevant quote: “You agree NOT to make the Repbase (or any part thereof, including Repbase Reports, Repeat Maps and other derived materials, modified or not) available to anyone outside your research group.” Emphasis mine. In fact, another clause in the agreement technically even forbids me from signing it because my institute requires public data deposition, so I’m probably not allowed to sign such agreements.”

Ya that seems to agree with my previous statement. I guess my point is what exactly did you need to share (based on your original post) that would be considered from their database? I feel like this wouldn't include counts of features but sequences might be an issue.

@story I need to potentially share all data that was used/generated in my analysis. This particularly includes the specific repeat annotation I used, which are derived from RepBase, as well as potentially sequence data from these repeats.

This might be an old question, but someone is trying to set-up a new, open, alternative to repBase (which i snow going full commercial), or at least that is how I perceive it: https://twitter.com/TransposableMan/status/1060519887897067521

Dfam has recently launched a sister resource, Dfam_consensus, whose stated aim is to replace RepBase. From the annoucement:

Dfam_consensus provides an open framework for the community to store both seed alignments (multiple alignments of instances for a given family) and the corresponding consensus sequence model.

Both RepeatMasker and RepeatModeler have been updated to support Dfam_consensus.
I haven’t tried it yet but it looks promising.

For pre-existing reliabe TE libraries it is a bit of a mess, because not everybody deposits the species-specific TE libraries to a database like RepBase. And as far as I know DFAM contains only human resources, or am I wrong?
As for de novo generation of species-specific TE libraries (which should be done for any species not already present in eg. RepBase):
There is no "gold-standard" how to tackle this best.
In principle one has to think about two main parts
-repeat detection
-annotation
For repeat detection I would recommend using a combination of two things (which is necessary, because TE copies might miss in the assemblies as repetitive regions tend to be difficult to assemble and thrown away in the final assembly).
I) Repeat detection from raw reads (as with e.g. DNApipeTE or tedna or RepeatExplorer). For me, DNAPipeTE worked quite nicely, but everything has pros and cons.
II) Repeat detection from assemblies (as with e.g. REPET or as mentioned before RepeatModeler)
Then the annotation of these repeats is tricky too, because most methods are relying on homology between the de novo TEs and the TEs from some (probably distantly) related species. But some programs take also structure into account (like REPCLASS). REPET can do both detection and annotation, but is a pain to get to run.
I would recommend using some programs to do de novo repeat detection on your species of interest on both the raw reads and assembly, clustering these libraries together (with e.g. uclust and 95% identity) and then run an annotation with homology and structural identification.
Probably the programs will not give you complete, full-length TEs but rather consensus sequences of several copies from TE families.
If you want, you could search all copies of one family, extract them from the contigs plus boundaries and align them manually and curate boundaries manually. Then extend boundaries if not hitting the surrounding (non-alignable) regions or landmarks of TEs like LTRs or TIRs or so. But this is very time consuming if you only want to compare TE abundance between species for example, I would not do this and rather compare the abundance using read coverage (as in Bast et al. 2016). Depends all on the questions you want to ask.

You could use RepeatScout, which has defined repeat libraries for a limited number of species (including human, mouse, and rat). If your taxon is not represented, you can also do de novo repeat prediction with RepeatScout to build your own library to feed to RepeatMasker. The RepeatScout publication includes some comparisons with RepBase. Another related tool is RepeatModeler, which wraps RepeatScout with RECON and some other programs, and shares authors with the RepeatMasker team.
On the plus side RepeatScout/RepeatModeler are open source and do not use manual curation, meeting your criteria. On the negative, I'm not sure exactly how RepeatModeler and the component tools are maintained. The RepeatScout web and github pages have not been updated for several years, although the RepeatModeler page shows its latest release was in 2017. Anyway, I know that some combination of RepeatScout/RepeatModeler have been used to annotate repeats for some fairly recent newly sequenced genomes, e.g. for cichlids, coelacanth, and Darwin's finch, so I think it's fair to say this kind of approach is accepted in the field, at least for vertebrate genome projects.

AFAIK Dfam and Repbase are currently the two best sources of (a variety of) TE sequences.
In my genome annotations I have used RepeatModeler+RepeatMasker and then later used Repbase+tblastx and Dfam+nhmmer to classify them.
The classification process in my pipeline PhyLTR (https://github.com/mcsimenc/PhyLTR) is based on Dfam and Repbase. The process I used for LTR identification is

Putative ID with LTRHarvest (based on structural sequence characteristics)
Classification-by-homology to Repbase and Dfam
Removal of elements without homology to sequences in Repbase or Dfam.

This results in a set of LTR-Rs which are full-length and have evidence that they are LTR-Rs.

+1 for taking issue with RepBase.
I use the annotations from the Hammell Lab GTFs that they put out with TEtoolkit. It is similar to what you described to be using so this may be a redundant and useless answer but from the digging I’ve done they do seem to be comprehensive and well curated (for Drosophila, at least).

I know this question is a bit old, but this is still an issue for a lot of researchers not being able to access RepBase. It seems now that the most recent version of RepeatMasker is dependent upon RepBase for full functionality if masking anything other than human (currently DFAM only has human models). I recently discovered a de novo repeat masking approach called REpeat Detector (Red).  This might be a solution for some looking to mask repeats on a genome assembly for annotation.  The paper is here. I also then wrote a wrapper around Red to make it a bit more easy to soft-mask a genome, which you can find here.
One of the limitations with Red is that the repeats are not classified, so they are only identified. You would have to use some of the other tools mentioned above to try to classify them.


text:
How are positions visually indicated in the UCSC Genome Browser?
I'm a little confused about how the position indices work in the browser. Based on the picture below, is base number 755440 an A or a G?
The general question is how do the indexing guides work? Is the base right under the | the one corresponding to the index or is it the one just before the |?

It's an A. The number refers to the position just before the |. This is clearer on different zoom levels, but I double checked by going to the very beginning of a chromosome:

As you can see in the image above, the | is always after the base it refers to. I also confirmed by checking a few random positions using another genome browser which places the positions unambiguously above each base.
For more details, see this UCSC blog post.

Yep and the reason for this is because we store everything internally as zero based coordinates but then display everything in a 1-based format. For more information see the following blog post: http://genome.ucsc.edu/blog/the-ucsc-genome-browser-coordinate-counting-systems/

@ChristopherLee thank you! I had actually written out an answer stating that it must be the position after since bed files start at 0 and then thought to check, just to make sure, and was very much surprised that it was actually the position before!


text:
Generating the reconstructed alignment from BAM
I have a (small) BAM file with CIGAR and MD fields.
Question 1: What tools exists in Python and/or R to reconstruct the alignment between the reference and the read in a BAM? Given that this is a very standard analysis, I feel that there should be R packages or Python libraries which have functions which do this task given a BAM input...
My goal was to create a *tsv file with two columns, one with the reconstructed alignment and the other column with the corresponding sequence. This would be useful in my analysis. (I would either use R data.table, or pandas in Python.)
Here is read and reference:
CGGGCCGGTCCCCCCCCGCCGGGTCCGCCCCCGGC
|||||||||||||||| |||||||||||||||||
CGGGCCGGTCCCCCCC-GCCGGGTCCGCCCCCGGG

Here is the desired output dataframe (with only the one example row). We could use an R data.table/data.frame, or a pandas DataFrame. Column1 are the alignments (i.e. the reconstruction of the alignment between the reference and the read based on the CIGAR and MD strings), Column 2 are the corresponding sequences. :
ReadAligned    Reference
CGGGCCGGTCCCCCCC-GCCGGGTCCGCCCCCGGG    CGGGCCGGTCCCCCCCCGCCGGGTCCGCCCCCGGC
...    ....

Question 2:
Using the CIGAR/MD tags, one could then create columns with the number of insertions and the number of deletions. This problem feels much more straightforward, e.g. the python library cigar or the cigar parsing functions of pysam

Could you clarify your question: You want to store the alignments in a matrix? Or the CIGAR format in a matrix?

@Llopis A tsv/csv file with two columns: the first column has the alignments, the second column has the sequence. Does this make sense?

Do you want sam2pairwise to generate the alignments for you, or do you want to parse the cigar strings yourself in R/python? Writing a cigar parser is tedious, but do-able. It might be easier if you just use sam2pairwise to output the alignments for each read, then you parse that into a table.

Now your question is reduced to "how to count the number of insertions and deletions" in cigar, which is much easier. With perl, it is just an one-liner echo 78M2I10M5D30M|perl -ne 's/(\d+)([ID])/$h{$2}+=$1/eg; print "$h{I}\t$h{D}\n"'. Someone will probably post python answers.

@user172818 I've tried to edit this to be more specific. Is this more organized? The second problem is easier; I'm primarily motivated by the first problem, which is related.

If all you want are the number of insertions/deletions, you can use this to extract that info:
from cigar import Cigar
c = c=Cigar('20M5I10M3I5M')
num_insertions = len(list(filter(lambda x: x[1] == 'I', c.items())))
insertion_length = sum([x[0] for x in list(filter(lambda x: x[1] == 'I', c.items()))])

Deletions would be the same thing with 'D' in place of 'I'

Is there a standard method to grab the string from the BAM within python? It looks like this is parsed automatically with pysam

I guess pysam would be the standard way

I am answering question 1. However, as I am unfamiliar with python, I am providing a javascript solution only to show the basic logic. To try the code at the end, you need to install node.js and run it with node this-script.js.

Given that this is a very standard analysis

I tend to believe this is a very rare analysis. I have only seen it done twice. Converting CIGAR+MD to padded alignment is complicated in coding, as you have to keep sequence, cigar and MD in sync all the time. The code below generates two strings in one go. It may be logically easier to generate the reference and the query strings separately.
In general, MD is ill defined. It is too difficult to work with. If possible, avoid this tag. You can do a lot of things with the combination of CIGAR and NM.
function md2pad(CIGAR, MD, SEQ) {
    var m, re_cigar = /(\d+)([MIDSHNX=])/g, re_MD = /(\d+)|(\^[A-Za-z]+)|([A-Za-z])/g;
    var cigar = [];
    while ((m = re_cigar.exec(CIGAR)) != null) // parse CIGAR into an array
        if (m[2] != 'H') // do nothing for hard clipping, as it has no effect
            cigar.push([parseInt(m[1]), m[2]]);
    var k = 0, sx = '', sy = ''; // k: k-th cigar operator; x: ref; y: query
    // cx/cy: start coordinate of cigar[k] on ref/query; mx/my: start of current MD op
    var cx = 0, cy = 0, mx = 0, my = 0;
    while ((m = re_MD.exec(MD)) != null) {
        if (m[2] != null) { // deletion from the reference
            var len = m[2].length - 1; // length of deletion
            sx += m[2].substr(1), sy += Array(len+1).join("-");
            mx += len, cx += len, ++k; // consume a deletion D from ref
        } else { // copy or mismatch
            var ml = m[1] != null? parseInt(m[1]) : 1; // length of this MD operation
            while (k < cigar.length && cigar[k][1] != 'D') {
                var cl = cigar[k][0], op = cigar[k][1];
                if (op == 'M' || op == 'X' || op == '=') {
                    if (my + ml < cy + cl) { // an MD ends in the middle of cigar M
                        if (ml > 0) {
                            if (m[3] != null) sx += m[3], sy += SEQ[my];
                            else sx += SEQ.substr(my, ml), sy += SEQ.substr(my, ml);
                        }
                        mx += ml, my += ml, ml = 0;
                        break;
                    } else {
                        var dl = cy + cl - my;
                        sx += SEQ.substr(my, dl), sy += SEQ.substr(my, dl);
                        cx += cl, cy += cl, ++k; // consume an M operator
                        mx += dl, my += dl, ml -= dl;
                    }
                } else if (op == 'I') {
                    sy += SEQ.substr(cy, cl), sx += Array(cl+1).join("-");
                    cy += cl, my += cl, ++k; // consume an insertion I to ref
                } else if (op == 'S') {
                    cy += cl, my += cl, ++k; // consume a soft-clipping S
                } else throw Error("inconsistent MD"); // not working with N
            }
            if (ml != 0) throw Error("inconsistent MD");
        }
    }
    if (cx != mx || cy != my) throw Error("inconsistent MD");
    return [sx, sy];
}

var cigar = '16025H3M2I55M1D108M306H';
var MD = '7T15C24T9^C0G2T6T16T3A1A7G21T17G6C1T6T5A4';
var seq = 'GATCACAGGCCTATCACCCTATTAATCACTCACGGGAGCTCTCCATGCATCTGGTATTTTTTCGGGGGGGGATGCACGCGATAGCATCGCGGGCCGCTGGAACCGGAGCACCCTATGTCGCAGGATCTGTCTTTGATTCCTACCTCATGCCATTATTAATCGCGCCTA';
var s = md2pad(cigar, MD, seq);
console.log(s[0]);
console.log(s[1]);

great reference code with the caveat that MD is indeed really annoying to work with. oftentimes easier to just recreate the mismatches from the sequence itself

This program in C++ is able to perform the task https://github.com/mlafave/sam2pairwise


text:
How to compute RPKM in R?
I have the following data of fragment counts for each gene in 16 samples:
> str(expression)
'data.frame':   42412 obs. of  16 variables:
 $ sample1 : int  4555 49 122 351 53 27 1 0 0 2513 ...
 $ sample2 : int  2991 51 55 94 49 10 55 0 0 978 ...
 $ sample3 : int  3762 28 136 321 94 12 15 0 0 2181 ...
 $ sample4 : int  4845 43 193 361 81 48 9 0 0 2883 ...
 $ sample5 : int  2920 24 104 151 50 20 32 0 0 1743 ...
 $ sample6 : int  4157 11 135 324 58 26 4 0 0 2364 ...
 $ sample7 : int  3000 19 155 242 57 12 18 2 0 1946 ...
 $ sample8 : int  5644 30 227 504 91 37 11 0 0 2988 ...
 $ sample9 : int  2808 65 247 93 272 38 1108 1 0 1430 ...
 $ sample10: int  2458 37 163 64 150 29 729 2 1 1049 ...
 $ sample11: int  2064 30 123 51 142 23 637 0 0 1169 ...
 $ sample12: int  1945 63 209 40 171 41 688 3 2 749 ...
 $ sample13: int  2015 57 432 82 104 47 948 4 0 1171 ...
 $ sample14: int  2550 54 177 59 201 36 730 0 0 1474 ...
 $ sample15: int  2425 90 279 73 358 34 1052 3 3 1027 ...
 $ sample16: int  2343 56 365 67 161 43 877 3 1 1333 ...

How do I compute RPKM values from these?

What did you tried to solve this question? Did you visit Bioconductor?

First off,
Don’t use RPKMs.
They are truly deprecated because they’re confusing once it comes to paired-end reads. If anything, use FPKMs, which are mathematically the same but use a more correct name (do we count paired reads separately? No, we count fragments).
Even better, use TPM (= transcripts per million), or an appropriate cross-library normalisation method. TMP is defined as:
$$
\text{TPM}_\color{orchid}i =
    {\color{dodgerblue}{\frac{x_\color{orchid}i}{{l_\text{eff}}_\color{orchid}i}}}
    \cdot
    \frac{1}{\sum_\color{tomato}j \color{dodgerblue}{\frac{x_\color{tomato}j}{{l_\text{eff}}_\color{tomato}j}}}
    \cdot
    \color{darkcyan}{10^6}
$$
where

$\color{orchid}i$: transcript index,
$x_i$: transcript raw count,
$\color{tomato}j$ iterates over all (known) transcripts,
$\color{dodgerblue}{\frac{x_k}{{l_\text{eff}}_k}}$: rate of fragment coverage per nucleobase ($l_\text{eff}$ being the effective length),
$\color{darkcyan}{10^6}$: scaling factor (= “per millions”).

That said, FPKM an be calculated in R as follows. Note that most of the calculation happens in log transformed number space, to avoid numerical instability:
fpkm = function (counts, effective_lengths) {
    exp(log(counts) - log(effective_lengths) - log(sum(counts)) + log(1E9))
}

Here, the effective length is the transcript length minus the mean fragment length plus 1; that is, all the possible positions of an average fragment inside the transcript, which equals the number of all distinct fragments that can be sampled from a transcript.
This function handles one library at a time. I (and others) argue that this is the way functions should be written. If you want to apply the code to multiple libraries, nothing is easier using ‹dplyr›:
tidy_expression = tidy_expression %>%
    group_by(Sample) %>%
    mutate(FPKM = fpkm(Count, col_data$Lengths))

However, the data in the question isn’t in tidy data format, so we first need to transform it accordingly using ‹tidyr›:
tidy_expression = expression %>%
    gather(Sample, Count)

This equation fails if all your counts are zero; instead of zeros you will get a vector of NaNs. You might want to account for that.

And I mentioned that TPMs are superior, so here’s their function as well:
tpm = function (counts, effective_lengths) {
    rate = log(counts) - log(effective_lengths)
    exp(rate - log(sum(exp(rate))) + log(1E6))
}

Can I ask sort of a meta question?  I've started seeing '%>%' in R code recently and had never noticed it before.  What does that do exactly?

@Greg It’s a pipe. It has been around for quite a long time (though different libraries used different operator symbols; I myself used to use %|%, similar to the shell pipe) but only recently gained mainstream traction, mainly through the dplyr library.

RPKM is defined as:

RPKM =   numberOfReads / ( geneLength/1000 * totalNumReads/1,000,000 )

As you can see, you need to have gene lengths for every gene.
Let's say geneLength is a vector which have the same number of rows as your data.frame, and every value of the vector corresponds to a gene (row) in expression.
expression.rpkm <- data.frame(sapply(expression, function(column) 10^9 * column / geneLength / sum(column)))

Regarding numerical stability
It is suggested in one of the answers, that all the computations should be performed in a log-transformed scale. In my opinion there is absolutely no reason for doing that. IEEE binary64 stores a number as binary number 1.b_{51}b{50}...b_0 times 2^{e-1023}. The relative precision doesn't depend on the exponent value given that a number is in the range [~10^-308; 10^308].
In case of RPKM we can only get out of the range if total number of reads is around 10^300, which is not realistic at all.
On the bright site there is not much harm in doing computations in the log-scale either. Worst case you loose a bit of precision.

and then picking the right gene length itself can be a nightmare depending on the organism...   So maybe you should also develop what could be the strategies to pick the correct one.

@Mitra I think this is slightly out of scope. But if you have an experience in doing that, could you add an answer here? Would be really great!

@Mitra given it a second thought, this deserves a separate question. If there will be one, I'll link it from this answer.

I don't think it is out of scope: you need the gene length to compute the RPKM. Other answers bypassed this problem by using transcript or exon length which are easily retrieved from the annotation. That is not the case with gene length which you are using in your formula.

If you are planning to do a differential expression analysis, you will probably don't need the RPKM calculation.
RPK= No.of Mapped reads/ length of transcript in kb (transcript length/1000)
RPKM = RPK/total no.of reads in million (total no of reads/ 1000000)
The whole formula together:
RPKM = (10^9 * C)/(N * L)
Where,
C = Number of reads mapped to a gene
N = Total mapped reads in the experiment
L = exon length in base-pairs for a gene

If you are looking for a more visual solution (in addition to the other answers), NCI Genomic Data Commons (TCGA repository) offers a nice formula:

Where:
RCg: Number of reads mapped to the gene
RCpc: Number of reads mapped to all protein-coding genes
RCg75: The 75th percentile read count value for genes in the sample
L: Length of the gene in base pairs

As others have pointed out, FPKMs have some problems. GDC also calculates FPKM-UQ values that are upper quartile normalized. Those are recommended for cross-sample comparison and differential expression analysis.

The denominator for FPKM-UQ is still sample specific, so this will not be an appropriate normalisation for differential expression analysis (it's within sample normalisation, not between).

GDC seems to imply otherwise: https://gdc.cancer.gov/about-data/data-harmonization-and-generation/genomic-data-harmonization/high-level-data-generation/rna-seq-quantification