text:
How to correct alpha, and not p-values themselves, for visualization purposes
I have a set of differentially methylated/expressed/whatever entities with p-values attached (example below).
entity_name    p-value    magnitude
entity1        0.04459    0.68
entity2        0.02283    0.99
...
entity_n       0.78       0.025

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

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

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

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

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

Aren’t alpha thresholds usually exclusive? I.e. it should be < 0.05, not <= 0.05.

If it matters, you're not using p-values correctly. See the ASA position statement.

@KonradRudolph Yes, good catch, I've fixed that. I had >= and min() in my mind apparently.

When you're doing a p-value adjustment, the same unadjusted p-value in different genes can be given different adjusted p-values depending on other factors. That means that you can't directly draw a line associated with the FDR on a plot of unadjusted p-value.
One possibility would be to take a range of values that are close to the FDR threshold (e.g. the 20 values closest to threshold), and draw a p-value greyzone within that region:
#!/usr/bin/Rscript
values <- c(rnorm(10000),rnorm(100, mean=1.5));
val.mean <- median(values);
val.diffs <- abs(values - median(values));
val.reldiffs <- (values - median(values));

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

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

Why the dummy <- assignment?

It's my default workaround to stop reporting null device / 1 to standard out. I've changed it to invisible(dev.off()), which does the same thing (and is probably more descriptive about its purpose).

Ah. The writeout never bothered me here. It’s not like it would be written to the output in an actual script/document.

Sometimes I send stdout to files and parse it. I've just got into the habit of getting rid of things like this so that I don't have weird errors cropping up for unknown reasons. I also have an issue with the scan function, which I silence when I use it.


text:
How to perform functional analysis on a gene list in R?
From an RNA-seq experiment I have about 17000 gene ids for 2 sample conditions arranged according to their log2 fold changes when compared to a control. I need to annotate these, but I've never done annotation before and am wondering how to do this in R.
I'm interested in up and down regulation of the pathways of the involved genes. I'd love to get enrichment and p-value analysis too. I'm interested in whether or not a package similar to PANTHER exists in R for this analysis. I've been reading how different annotation programs can be better or worse; it seems PANTHER is much better than DAVID, so I was wondering if there is a package that provides PANTHER-like analyses in R.
There seems to be multiple packages available, are there any that stand out as being the best?
I'm primarily interested in human samples and annotated pathways.

I've updated your question with that information. What gene IDs do you currently have (e.g., Ensembl, refseq, etc.)?

How do you define best? What do you want to know, the genes that are coordinated, the pathways that are down-regulated or up-regulated? You want a GSEA single-sample or whole sample, competitive or non-competitive? For me this question is like I have a list of number I want to do some operation to get the most important ones, which method is available?

@J0HN_TIT0R to answer a commenter use @ nickname. Please [edit] the question to include this new information.

As it stands, this question is not answerable. @J0HN_TIT0R you don't provide enough background information as to which function you are interested in and why.: this very much decides which enrichment and enrichment tool we can recommend. Also, questions about "the best tool" are discouraged on SE because they are primarily opinion-based.

N.B., I'm avoiding discussion of "best", since that's more or less impossible to answer.
Your question can actually be divided into two:

What's a good tool for pathways analysis (ideally in R)?
What are good sources of pathway information?

For (1), there are a few different possibilities, but I prefer either roast or camera from limma or goseq (a stand-alone bioconductor package). Of those, I expect what roast or goseq are doing are the most similar to what panther is doing in the sense that they're not competitive. However, since genes do correlate with each other, I prefer competitive tests like that offered by camera.
For (2), your options are basically KEGG or reactome. KEGG is problematic, since you need a license for anything remotely recent (there are some R packages to get around this, but I haven't a clue how kosher they are legally). Given that, reactome might be your best bet.
Having written all of that, the best database out there is what IPA uses. This is a commercial product (I have no affiliation), but you can get a demo license that will work for a couple weeks at least. If you're doing analyses relatively often then it makes sense to spend a bit of cash on a license (maybe share with neighboring groups).
Regardless of what you use, you might also want to use pathview, which can produce easily interpretable plots related to pathway enrichment.

I kindly disagree with IPA being best (I thought you wanted to avoid 'best' in your answer). IPA is commercial and is not transparent, but a black box. You have to trust that the company did a good job making this database behind scenes, which is not really how science should be. Just my thoughts.

You requested a tool similar to PANTHER but in R. First, the PANTHER (Protein ANalysis THrough Evolutionary Relationships) tool does a classification  based on evolutionarily related proteins, gene ontologies ( molecular function, and biological process) and pathways.
AFAIK there isn't a tool in R that integrates all these into one, but several packages do pieces of it. However, my observations after using it a couple of times is that is based a lot in gene ontologies. Gene Ontologies (GO) are not pathways, but the best tool I found for GO are topGO (A bit hard to use but useful to get insights on the biology) which test over representations of a term in a given list.

I would go for the topGO package. The logical steps in your pipeline would be:
1) Build a matrix with the following variables:
|gene_ID|GO_ID|GO_term_description|
You can enter a list of uniprot IDs to uniprotkb and retrieve the GO annotations for each gene in your list/file
2) filter that list for the genes that were captured in your differential expression analysis even tho they are not differentially expressed. This filter is defined by DESeq/Sleuth by "at least 4 read counts in at least 49% of the libraries"
3)Define this gene list as gene_universe
4) define your gene set and get the gene_IDs and if possible, also a ranked score (p-val/FDR) of the genes you are interested in
5) plot the graphs and get tables of enriched GO terms (these would be ranked by p-value without need of performing correction analysis
Note: you are going to need some data manipulation functions; I would read the "tidyverse" documentation

source: Stack Exchange

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

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

sort: votes


text:
Why are there missing calls in a VCF file from exome sequencing?
My data is a VCF file generated from an exome sequencing variant call pipeline. I'm not very familiar with the sequencing and variant calling process. I noticed that there are some missing genotypes, which are recorded as "./." at the GT field. From googling I learned that they're not homozygous reference genotype ("0/0"), but missing calls due to some sequencing failure.
In order to make my data "cleaner", I thought it would be better to filter out the loci with missing calls if there're not too many of them. I also checked the corresponding "DP" of the loci with and without missing calls. For example:
chr1:123 GT:DP 0/0:2 1/1:4 ./.
chr1:234 GT:DP 0/0:10 1/1:11 1/1:20
chr1:345 GT:DP 0/1:40 1/1:37 0/0:78
chr1:456 GT:DP 0/1:7 0/0:23 ./.
chr1:567 GT:DP 0/1:34 1/1:39 0/0:58

In the above toy example, there are 3 people and 5 loci. I checked the mean DP of all 5 loci and I found that the mean DP of loci with missing calls (1st and 4th loci) to be significantly lower than the mean DP of loci without missing calls (2nd, 3rd and 5th loci). This happens to my real data with a lot more loci and samples than this toy example. Is that a coincidence? Or are there any specific reason that the loci with missing calls have lower coverage than the "normal" loci? Thanks!

Hi Yan, thanks for your question and welcome to Bioinformatics Stack Exchange. I have modified your question slightly to fix up grammar and make it clearer what is being asked; please feel free to change your question back if you don't think the changes are appropriate.

Missing variant calls due to lack of coverage shouldn't happen in the targeted capture region and I'd think most of these would come from off-target regions where some samples had reads mapped. I'd filter out the VCF to only include on-target loci before proceeding with further analyses.
If you did your variant calling on samples separately and then merged them together, you'll have missing calls at loci called variant in one sample and homozygous reference in the other.
You can get over this by asking the variant caller to emit all sites within the target region, which would lead to a large VCF or do multi-sample calling, which is more computationally intensive.
One other way is to do variant calling on all samples separately, generating a BED file of all variant regions in the cohort and then asking the variant caller to genotype those loci across all samples before merging them together.

I've just been generating data like this, so can tell you about why/how missing calls are created in my dataset. There are two main reasons:
1. Sequencing failure
When reads don't map across the variant region, then it's impossible to accurately determine a genotype for that region. This will commonly happen just outside the borders of the selected regions for exome sequencing, but can also happen through natural random sampling of the genome, or through repeated sequences coupled with systematic error in either the assembled genome or the sequencer. In this case, it would be expected that called variant regions including missing data would be lower coverage. It may be helpful to plot the coverage within a particular region to work out if this is likely to be happening.
2. Merged non-variant data
Variant calling can be processed in two different modes, one which takes into account all samples when doing calling, and another which does the calling one sample at a time, followed by subsequent merging of the results. Multi-sample variant calling is more accurate, but computational limits (or experimental design) might mean that it is more appropriate for samples to be called one-by-one. It's fairly common (for the purpose of preserving space and computational power) for calling algorithms to only report variant information so if a single sample is identical to the reference, and that sample is the only one being called, then the variant and coverage information for that sample will be lost. When samples are merged after being separately called, the declared call for the "same as reference" samples will be set to missing. In this case, it would be generally expected that the other called samples would have high-coverage within the variant region.

Thanks gringer, for the detailed explanation and the grammer editing of my post. Both answers are very good so I just chose the first one and I hope you don't mind! :) I don't have enough reputation to up vote the answer right now but I will when I do.

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


text:
Given a transcription factor, what genes does it regulate?
I have a list of transcription factors and I am interested in finding out which which genes might be transcribed as a result of the formation of transcription factor complexes with that transcription factor.
Any ideas on good databases? I've seen that ENCODE hosts data from CHIPseq experiments, but I'm not sure if I can find conclusions regarding interactions from this site. I'm looking for a database with known or putative transcription/gene interactions.
Thanks in advance.

It's challenging to know all the true targets of a TF. See this answer.

MSigDB has a collection (C3:TFT) of gene sets corresponding to transcription factor targets.
Harmonizome has functional terms for genes extracted from over a hundred publicly available resources.

iRegulon takes a sequence-based approach to finding transcription factor targets. There's a Cytoscape app that you can use to find the regulators of a given gene list, or the targets of a particular transcription factor.
Transcription factor binding sites are predicted using a collection of position weight matrices (PWMs) from a number of sources, including Jaspar and TRANSFAC.
iRegulon: From a Gene List to a Gene Regulatory Network Using Large Motif and Track Collections PLoS Comput Biol. 2014 Jul 24;10(7):e1003731. doi: 10.1371/journal.pcbi.1003731. eCollection 2014 Jul.

In general, there is two options to identify targets for transcription factors: experimental (ChIP-seq) and sequence-based predictions.
TF binding from experimental data
The are multiple projects that produce binding data of transcription factors and quantify their peaks across the genome. The advantage here is that you know that binding actually occurs, as opposed to motif predictions. But you need the corresponding experiments.
You already mentioned ENCODE, which is probably the biggest producer of such data. NIH's Roadmap Epigenomics is another one, but that focusses less on TFs.
If you have a gene list and would like to know which transcription factor was likely involved, you can do an enrichment test of your gene set in known TF targets. The ChEA (ChIP Enrichment Analysis) database does this.
Prediction using motifs
Another possibility is to look at binding motifs and see whether your gene list is enriched in those. These will, however, be inactive in a given tissue or cell type if the chromatin is packed or the methylation state of the promotor is unfavorable.
Examples for those are the JASPAR and TRANSFAC databases, or the MSigDB motif gene set (as burger mentioned). You can also query those features using the Ensembl genome database (BioMart, REST).
Calculating enrichment in gene sets
Most likely, an analysis you will want to perform is to calculate the enrichment in your gene list, e.g. that you got from differential expression in two conditions.
The most convenient way is to use the Enrichr platform, which is a web page that accepts a gene list and will compute enrichment in ChEA, JASPAR, TRANSFAC, etc. You can also download their gene sets.

TRUSST is a manually curated database with around 790 TFs and their target genes. TRUSST Website


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.

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.

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.

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.


text:
how to find the bound form of an enzyme structure?
For my undergraduate research I'm looking for a database that gives the bound form of a particular protein structure. Is there any database that provide us with such data? So far I've found following proteins related to my work

Ascorbate peroxidase (APX)
Beta-glucosidase
Calcineurin b-like protein-interacting protein kinase.

I've searched in PDB and UNIPROT.

Hi Lock23, thanks for your question and welcome to Bioinformatics Stack Exchange. After reading your question, I did a quick check on the Biology and Chemistry Stack Exchanges, and was surprised at how little "database" stuff there was.

For experimentally determined structures PDB is the database and it's unlikely you'll find structures that are publicly available but not in the PDB. There is also protein structure prediction software (servers).

Tried looking for an explicit database? i.e.

ComSin: database of protein structures in bound (complex) and unbound (single) states in relation to their intrinsic disorder: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2808974/ or
LigASite: a database of biologically relevant binding sites in proteins with known apo-structures https://www.ncbi.nlm.nih.gov/pubmed/17933762

Else the PDB itself has a very powerful advanced search that should allow you to do this: http://www.rcsb.org/pdb/search/advSearch.do

Are you aware of BRENDA? I was just introduced to it today for a completely separate reason (looking at carbohydrate enzyme families in the Nippostrongylus brasiliensis proteome), and it seems to be a fairly comprehensive database. There is at least a literature link there for the ascorbate-complexed crystalisation of Ascorbate peroxidase on that site. Following through the paper, I see that the complex has been uploaded to PDB.


text:
Get the mapping statistics of a single read from a BAM file
I have a BAM file, and I have a read ID. What is the simplest way to get mapping statistics of that read in human-readable format?
E.g. I might want: % identity of aligned bases; number of insertions and deletions; number of bases aligned to the reference; etc.
Extra points for one-liners that give detailed summary statistics, and extra-extra points if you can extend the answer to output summaries for all reads longer than some predetermined length.

samtools stats seems to be able to do most of this, excluding the CIGAR-string parsing stuff (i.e. INDELs):
$ samtools view -h mapped.bam | grep -e '^@' -e 'readName' |
    samtools stats | grep '^SN' | cut -f 2-

raw total sequences:    2
filtered sequences: 0
sequences:  2
is sorted:  1
1st fragments:  2
last fragments: 0
reads mapped:   2
reads mapped and paired:    0   # paired-end technology bit set + both mates mapped
reads unmapped: 0
reads properly paired:  0   # proper-pair bit set
reads paired:   0   # paired-end technology bit set
reads duplicated:   0   # PCR or optical duplicate bit set
reads MQ0:  0   # mapped and MQ=0
reads QC failed:    0
non-primary alignments: 0
total length:   13632   # ignores clipping
bases mapped:   13632   # ignores clipping
bases mapped (cigar):   13502   # more accurate
bases trimmed:  0
bases duplicated:   0
mismatches: 3670    # from NM fields
error rate: 2.718116e-01    # mismatches / bases mapped (cigar)
average length: 6816
maximum length: 6816
average quality:    8.9
insert size average:    0.0
insert size standard deviation: 0.0
inward oriented pairs:  0
outward oriented pairs: 0
pairs with other orientation:   0
pairs on different chromosomes: 0

For a length filter, an awk length check can be slotted in instead of the grep:
$ samtools view -h mapped.bam |
    awk -F'\t' '{if((/^@/) || (length($10)>500)){print $0}}' |
    samtools stats | grep '^SN' | cut -f 2-

raw total sequences:    131
filtered sequences: 0
sequences:  131
is sorted:  1
1st fragments:  131
last fragments: 0
reads mapped:   131
reads mapped and paired:    0   # paired-end technology bit set + both mates mapped
reads unmapped: 0
reads properly paired:  0   # proper-pair bit set
reads paired:   0   # paired-end technology bit set
reads duplicated:   0   # PCR or optical duplicate bit set
reads MQ0:  0   # mapped and MQ=0
reads QC failed:    0
non-primary alignments: 0
total length:   569425  # ignores clipping
bases mapped:   569425  # ignores clipping
bases mapped (cigar):   453884  # more accurate
bases trimmed:  0
bases duplicated:   0
mismatches: 113755  # from NM fields
error rate: 2.506257e-01    # mismatches / bases mapped (cigar)
average length: 4346
maximum length: 16909
average quality:    9.1
insert size average:    0.0
insert size standard deviation: 0.0
inward oriented pairs:  0
outward oriented pairs: 0
pairs with other orientation:   0
pairs on different chromosomes: 0

all extra points awarded. This is (hopefully) super useful for folks working on long reads, and wanting to summarise the longer bits of their data...

@terdon please don't edit this to make it less readable than it already is. I'm not interested in playing awk golf; I want to show people how it works and make the logic explicit.

Fair enough. I tend to favor brevity in one-liners, but you're right in that it's less clear to the non-expert. Sorry. I would urge you to remove the \  since they are not needed and, IMO, harm readability though. I would also remove the needless parentheses from the awk call for the same reason: {if(/^@/ || length($10)>500){print}}

I've removed the end of line escapes, but not the brackets. I'm not completely comfortable doing either, because I've had my code do unexpected things for both reasons in the past (although obviously not these exact situations).

Use htsbox:
htsbox samview -p in.bam | less -S
htsbox samview -pS in.sam | less -S

It outputs mapping positions in the PAF format, which looks something like:
read1  4983  774  4982  +  chr18  80373285  26911072  26915544  3835  4631  60 \
       mm:i:214  io:i:119  in:i:159  do:i:339  dn:i:423

This is one long line in terminal. I folded the line for display purposes. For the first 12 fixed fields, see the table in PAF page. They give mapping positions and identity. The optional fields tell you #mismatches (mm), #insOpens (io), #insertions (in), #delOpens (do) and #deletions (dn). Note that the alignment MUST have an "NM" tag; otherwise the number of matching bases (col. 11) is overestimated and "mm" will always be zero.
You can easily compute summary statistics from PAF if you are familiar with command lines. For example:
htsbox samview -p in.bam | awk '{x+=$10;y+=$11}END{print x/y}'  # identity
htsbox samview -p in.bam | awk '$2>1000{x+=$10;y+=$11}END{print x/y}'  # identity for >1000bp
htsbox samview -p in.bam \
  | perl -ane '{$g+=$1 while /[id]n:i:(\d+)/g;$y+=$F[10]}END{print $g/$y,"\n"}' # gap rate

If you don't like long one-liners, feel free to write a proper script to collect all statistics in one go.


text:
What methods are available to find a cutoff value for non-expressed genes in RNA-seq?
I have a gene expression count matrix produced from bulk RNA-seq data. I'd like to find genes that were not expressed in a group of samples and were expressed in another group.
The problem of course is that not all effectively non-expressed genes will have 0 counts due to sequencing errors, or because they were expressed in a small subset of cells.
I'm interested in solutions using R.

I'd like to find genes that were not expressed in a group of samples and were expressed in another group.

This is, fundamentally, a differential expression analysis, with a twist. To solve this, you’d first use a differential expression library of your choice (e.g. DESeq2) and perform a one-tailed test of differential expression.
Briefly, you’d perform the normal setup and then use
results(dds, altHypothesis = 'greater')

To perform a one-tailed test. This will give you only those genes that are significantly upregulated in one group. Check chapter 3.9 of the vignette for details.
Of course this won’t tell you that the genes are unexpressed in the other group. Unfortunately I don’t know of a good value to threshold the results; I would start by plotting a histogram of the (variance stabilised) expression values in your first group, and then visually choose an expression threshold that cleanly separates genes that are clearly expressed from zeros:
vst_counts = assay(vst(dds))
dens = density(vst_counts[, replicate])
plot(dens, log = 'y')

(This merges the replicates in the group, which should be fine.)
Counts follow a multimodal distribution, with one mode for unexpressed and one or more for expressed genes. The expression threshold can be set somewhere between the clearly unexpressed and expressed peaks:

Here I used identify(dens) to identify the threshold interactively but you could also use an analytical method:
threshold = identify(dens)
quantile = sum(dens$x < dens$x[threshold]) / length(dens$x)

# Using just one replicate here; more robust would be to use a mean value.
nonzero_counts = counts(dds, normalized = TRUE)[, replicates[1]]
nonzero_counts = nonzero_counts[nonzero_counts > 0]

(expression_threshold = quantile(nonzero_counts, probs = quantile))

26.5625%
4.112033

I like this approach but I have a few concerns. In the DE step, how does the differential expression test perform with zero count values? If one group has zero values and the other group has near zero, or clearly expressed, but with high variability etc.
In the cutoff step, (assuming that the blue line is the threshold) it seems you have a very low cutoff value (~1), so you will basically accept all DE genes (or accept all genes above count 0)?

Also the graph is not bimodal, but has 4 local maxima. The second peak corresponds to roughly 1-2 counts, and when you have 20 million reads, they are probably not not expressed. Is it not better to have the cutoff after the 3rd peak, around rlog=3?

@Peter Your observation about the different modes is a good one that I glossed over. In fact, the distinction should be between a mode for unexpressed genes, and “all the rest”. In order to determine which mode corresponds to truly unexpressed genes, it might be necessary to plot the density with different bandwidths and see how the modes change. My plot isn’t perfect in this regard since I used simulated data based on averaged counts. Real counts would behave better.

@Peter Your comments made many good points so I redid the example with actual read count data, and used VST instead of rlog, which appears to work somewhat better here.

A common method is to use zFPKMs, which you can find implemented in R here.
Having said that, there's an inherent problem in declaring a difference between two things on either side of a given threshold. Given that, I would encourage you to use more than "expressed in one and not in another", likely adding at least a "with a minimal difference of X between them" metric. You may also find the tau metric useful, which is implemented in R here. This is meant to measure tissue-specificity, which is more akin to what you're probably interested in doing.

Depending on how much effort you wish to put into this, here is one suggestion I have see used before (uses more than just R). Steps 1-4 come from this paper (see supplimentary material section "Calculation of per gene Local FDR").

Take your gene models and collapse the introns down to maybe 100bp
Thus if we have a gene with two exons and a 1kb intron eg. exons (1000, 1100) and (2000,2100), we reduce the intron size so that the exons are (1000,1100) and (1200,1300). This is because we need to find sufficient gene free space to fit the null gene in.
Shift each gene into the nearest genomic space at least 5 kp from a genomic annotation and free of ESTs.
Calculate the FPKM distribution of this shifted null set
From these two distributions it should be possible to calculate a local FDR for expression.
You will need to use a "mixing proportion" to do this. The above reference used Qvality to do this.
Carry out differential expression analysis as described in Konrad's answer
Subset the differentially expressed genes to only consider those that are not expressed in the control condition
You might want to also recalculate the FDRs for differential expression as you are only considering a subset of the tests you might have.

This sounds like a promising approach for calculating a threshold. However,  I don’t entirely understand your step 1, could you expand on that?

Edited to include more detail on step 1.

Yes the calculation of local FDR is something that has come up lately and been supported both by Gordon Smith and Mike Love. In any case if you do not want to do that one can always start from a list of genes that are above a certain threshold of expression levels to restrict the universe of expressed genes and then perform your DE with DESeq2.

Unfortunatly that is not doable - there is now way of distinguishing between features not expressed and features not expressed for technical reasons. Furthermore a lot of genes might seem lowly expressed to to technical noise.
This also means that you will have to choose a cutoff and simply use that. For choosing a cutoff looking at the replicate distribution of log10(counts) and log10(FPKM) are probably the best options.

I think you’ve posted your answer on the wrong question: did you mean to post it here instead? https://bioinformatics.stackexchange.com/q/693/29

@KonradRudolph Actually, the linked question did not exist when this answer was posted.


text:
How can I convert SCF trace files to ABI files?
I have a large collection of sequence trace files which are in SCF file format (binary files that in addition to the string of DNA bases also contain the electropherogram of the sample and quality information about the base calls).
I have an existing BioPython pipeline that accepts files in AB1 format (pdf), so I need to convert my SCF files to AB1 files.
I found a question on Biostars that asks for the opposite: converting from SCF to AB1.
The accepted answer there was to use convert_trace from the Staden package.
./convert_trace -out_format scf < trace.ab1 > trace.scf

So I tried the opposite:
./convert_trace -out_format abi < trace.scf > trace.ab1

This did indeed produce a file with *.ab1 extension, but upon importing it to my BioPython pipeline, I got the error
File should start ABIF, not '\xaeZTR'

So apparently, convert_trace converted my files to ZTR format, even altough I specified -out_format abi.
Closer reading of the manual states that [convert_exe] can write CTF, EXP, PLN, SCF and ZTR formats. (so indeed no ABI), but it confuses me why ABI is still a valid output format.
So, my question is, how do I convert SCF files into ABI files?

Apparently convert_trace does not do a good parameter checking and silently sets ZTR as default if it does not recognise a valid output format.
I don't know of any freely available command-line converters in that specific direction. Even at the height of Sanger sequencing (late 90s, early 2000s), pipelines usually tried to reduce complexity and footprint by converting AB1 to anything else, often SCF (ZTR came too late and never got enough traction).
That being said, I dug out an SCF I had from 1997, loaded it into Geneious, and re-exported it as .ab1. I was a bit astonished to see that this worked. So Geneious (or similar packages from other companies) may be your best bet for a quick, one-off thing.
Depending on your needs, you might want to think about other ways: either expand BioPython to read SCFs, or maybe to expand the Staden io_lib to be able to write .ab1 (James still maintains that package and is extremely helpful).
Good luck.

I downloaded Geneious and requested a trial license. I'm indeed able to import SCF files, but I don't see 'ab1' as an export format? So, I'm guessing that I will need to go with one of your other two options...

I checked again and my version of Geneious (10.1.3) definitely offers an "ABI chromatogram (*.ab1)" as type in the "File->Export selcted document ..." dialog. Was your SCF selected?

DNA Baser Assembler will convert ABI to SCF (and other formats).
Just double click the ABI file you want to convert then go to the "Samples" menu and choose "Save current samples as...". Plus, it can automatically trim the low quality ends.

There is another tool in that package what will convert all ABI files at batch: Chromatogram Explorer. There is a Lite version that is free, but the conversion to SCF might only work for the first 100 files (if I remember correctly).
Related to the comment posted by BioGeek about the Genious license: DNA Baser Assembler does not require a license key. It will work right away.

The OP was asking for SCF->ABI I think your response is only SCF->ABI.  Do those tools work in both directions?


text:
How can we distinguish between true zero and dropout-zero counts in single-cell RNA-seq?
In single-cell RNA-seq data we have an inflated number of 0 (or near-zero) counts due to low mRNA capture rate and other inefficiencies.
How can we decide which genes are 0 due to gene dropout (lack of measurement sensitivity), and which are genuinely not expressed in the cell?
Deeper sequencing does not solve this problem as shown on the below saturation curve of 10x Chromium data:

Also see Hicks et al. (2017) for a discussion of the problem:

Zero can arise in two ways:

the gene was not expressing any RNA (referred to as structural zeros) or
the RNA in the cell was not detected due to limitations of current experimental protocols (referred to as dropouts)

Possible duplicate of What methods are available to find a cutoff value for non-expressed genes in RNA-seq?

I don't think this is a duplicate: this question is about distinguishing between zero counts and the linked question is about distinguishing between zero and non-zero counts.

Definitely not a duplicate, this is a very different question.

Sorry, I missed the "single cell" mention. Despite that, I expect that some answers are likely to overlap, even if the problem is different. See, for example, Kristoffer's answer.... er, which may have been posted in the wrong question.

hey @Peter if the answers here are not satisfactory you could maybe add more info to the question or comment on the answers

Actually this is one of the main problems you have when analyzing scRNA-seq data, and there is no established method for dealing with this. Different (dedicated) algorithms deal with it in different ways, but mostly you rely on how good the error modelling of your software is (a great read is the review by Wagner, Regev & Yosef, esp. the section on "False negatives and overamplification"). There are a couple of options:

You can impute values, i.e. fill in the gaps on technical zeros. CIDR and scImpute do it directly. MAGIC and ZIFA project cells into a lower-dimensional space and use their similarity there to decide how to fill in the blanks.
Some people straight up exclude genes that are expressed in very low numbers. I can't give you citations off the top of my head, but many trajectory inference algorithms like monocle2 and SLICER have heuristics to choose informative genes for their analysis.
If the method you use for analysis doesn't model gene expression explicitly but uses some other distance method to quantify similarity between cells (like cosine distance, euclidean distance, correlation), then the noise introduced by dropout can be covered by the signal of genes that are highly expressed. Note that this is dangerous, as genes that are highly expressed are not necessarily informative.
ERCC spike ins can help you reduce technical noise, but I am not familiar with the Chromium protocol so maybe it doesn't apply there (?)

since we are speaking about noise, you might consider using a protocol with unique molecular identifiers. They remove the amplification errors almost completely, at least for the transcripts that you capture...
EDIT: Also, I would highly recommend using something more advanced than PCA to do the analysis. Software like the above-mentioned Monocle or destiny is easy to operate and increases the power of your analysis considerably.

Some people use imputation to differentiate between true zeros and dropout in single-cell data. Some approaches you can look into:

CIDR which was published in Genome Biology, here
scImpute , published in Nature, here
Seurat (addImputedScore) , i.e. the famous R package, also here
SAVER which was reviewed in Frontiers of Genetics, here


text:
Converter between PDB or mmCIF and MMTF
I'd like to test MMTF, a new format for storing biomolecular structures which is promoted by RCSB as a more compact alternative to mmCIF and PDB.
From MMTF FAQ:

How do I convert a PDBx/mmCIF file to an MMTF file?
The BioJava library contains methods to read and write PDBx/mmCIF files and MMTF files.

Can I do such a conversion, ideally from command line, but without writing my own Java program?

I wrote a very quick and dirty script to handle conversion between file types using BioJava.
https://github.com/eedlund/Utils/tree/master/BioUtils
Download the jar file here
To run: java -jar BioUtils.jar $FILE $TYPE
where \$FILE is a PDB or mmCIF file you'd like to convert and \$TYPE is the format of the output file [PDB, CIF, MMTF].

Does this take into account biological assemblies, or simply convert e.g. MMTF to a PDB of the asymmetric unit?

You looked at this https://github.com/rcsb/mmtf-python
"The python implementation of the MMTF API, decoder and encoder."

but how to use it to convert between MMTF and something else?

There's also a C++ encoder/decoder available here: https://github.com/rcsb/mmtf-cpp

You can do this with BioStructures.jl in Julia. All the 6 transformations between PDB/mmCIF/MMTF are possible.
For example, PDB to MMTF:
using BioStructures
struc = read(in_filepath, PDB)
writemmtf(out_filepath, struc)

mmCIF to MMTF:
using BioStructures
struc = read(in_filepath, MMCIF)
writemmtf(out_filepath, struc)

You can use PyMol to change the PDB to mmtf or other format. First open the PDB file and then choose file->export molecule-> choose which format you want to save as

This python script uses Biopython and does the job well even on large structures. For example:
python cif2pdb.py 4ckh-assembly-1.cif 4ckh-assembly-1.pdb
This generates the structure below.

The question was about converting to/from MMTF.
For PDB <-> mmCIF there are too many tools to list them all. For example, I'm maintaining gemmi which should be at least an order of magnitude faster than BioPython and it also converts more fields (but I'm obviously biased).

@marcin Thanks for the heads up. Since this is probably a more common conversion, I've opened a new question.


text:
Count genomic ranges
I have a set of genomic ranges that are potentially overlapping. I want to count the amount of ranges at certain positions using R.
I'm Pretty sure there are good solutions, but I seem to be unable to find them.
Solutions like cut or findIntervals don't achieve what I want as they only count on one vector or accumulate by all values <= break.
Also countMatches {GenomicRanges} doesn't seem to cover it.
Probably one could use Bedtools, but I don't want to leave R.
I could only come up with a hilariously slow solution
# generate test data
testdata <- data.frame(chrom = rep(seq(1,10),10),
                       starts = abs(rnorm(100, mean = 1, sd = 1)) * 1000,
                       ends = abs(rnorm(100, mean = 2, sd = 1)) * 2000)

# make sure that all end coordinates are bigger than start
# this is a requirement of the original data
testdata <- testdata[testdata$ends - testdata$starts > 0,]

# count overlapping ranges on certain positions
count.data <- lapply(unique(testdata$chrom), function(chromosome){
    tmp.inner <- lapply(seq(1,10000, by = 120), function(i){
        sum(testdata$chrom == chromosome & testdata$starts <= i & testdata$ends >= i)
    })
    return(unlist(tmp.inner))
})

# generate a data.frame containing all data
df.count.data <- ldply(count.data, rbind)

# ideally the chromosome will be columns and not rows
t(df.count.data)

So you only want overlaps where the start position is not contained within your interval range? What's the actual biological context to this?

You could start by optimizing the code by using data.table (and probably using vapply instead of lapply). Also if one nucleotide overlaps between ranges then you count that too?

@DevonRyan no what i want is that the respective counting position is with in the range, which is a genomic deletion (hence the biological context)

@Llopis yes i want that counted as well. The list contains overlapping areas and i want to know how many of them are overlapping a certain position. Which is then used later for plotting.

Hi sargas, thanks for your question and welcome to Bioinformatics Stack Exchange. If you have any additional story/context associated with your question (such as a desire to count genomic deletions), it can be helpful to include that in your question. This makes it easier to understand the question, and helps answerers to solve the problem you have rather than the question you're asking (see the XY problem for more information).

GenomicRanges::countOverlaps seems to be what you’re after:
position_range = GRanges(position$chrom, IRanges(position, position, width = 1))
ranges_at_position = countOverlaps(position_ranges, granges)


text:
Preparing binary matrix data for Scikit classification algorithms
I made this post in regular stack overflow but I was told about this awesome feature by @nbryans.
I am a researcher (my programming knowledge is small) conducting analysis on a set of antibiotic (methicillin) resistant and a set of antibiotic (methicillin) susceptible. The sole difference between these two sets is assumed to be the resistance (R) vs susceptibility (S). I want to find if a genetic element in the resistant or susceptible genomes is correlated to resistance (positively or negatively).
I have done a complete fisher exact test already and will detail how I did so below. However I would like to figure out how to prepare my data for a more complex learning algorithm such as those offered by sci-kit.
(many of my techniques were learned from this article)
First I used KMC to split each fasta file for both R and S genomes into k-mers (I don't have enough "reputation" so I can't use >2 links but k-mers are essentially genetic sequences of k length). The output for KMC was a dump of k-mers present in a genome.
Example: kmc dump
AAAAAAAAAACGAATGTACACAATCGAA    1
AAAAAAAAACGAATGTACACAATCGAAC    1
AAAAAAAAATCAAATCCTGACTATTTAG    1
AAAAAAAAATCGGTCAATTCATTAAAAG    1
AAAAAAAACAAACATGAAGACCTTGTTA    1
AAAAAAAACGAATGTACACAATCGAACA    1
AAAAAAAACGTGTTAAAGTGAATCACAC    1

I then merged the data into one binary matrix which had columns as genomes and rows as k-mers.
Example: 1 resistant genome and 1 susceptible genome
kmer    R_1    S_1
AAAAAAAAAACGAATGTACACAATCGAA    1   0
AAAAAAAAACGAATGTACACAATCGAAC    1   0
AAAAAAAAATCAAATCCTGACTATTTAG    1   0
AAAAAAAAATCGGTCAATTCATTAAAAG    1   0
AAAAAAAAATTCCCTTCTAATCTTGAAT    0   1
AAAAAAAACAAAAATTATATAAAGCGAA    0   1
AAAAAAAACAAACATGAAGACCTTGTTA    1   1
AAAAAAAACAACCACCCATACATTGAGT    0   1
AAAAAAAACCCTTACAACAAATATGTAA    0   1

This is the data I have been working with and would like to analyze further with sci kit or other classification algorithms. Essentially each k-mer is a feature whose correlation with resistance is based on its presence in R vs S genomes. This is where I would like help transitioning this data into a dataset format suitable for sci kit.
I conducted a Fisher Exact test on each k-mer using R. For each k-mer a 2x2 matrix was created with the first column as resistant and the second column as susceptible. The two rows were present or not. So four numbers (# times k-mer was present in R, # in susceptible, # not in R, # not in S).
      R    S
yes   #    #
no    #    #

I used the following code in R:
phenotype = as.numeric(grepl('^R_', colnames(raw_kmer_table)[2:ncol(raw_kmer_table)]))
fe_results = apply(kmer_pres_abs_matrix, 1,
               FUN = function(row) {
                 fe_mat = matrix(0, ncol = 2, nrow = 2)
                 fe_mat[1,1] = sum(row == 1 & phenotype == 1)
                 fe_mat[1,2] = sum(row == 1 & phenotype == 0)
                 fe_mat[2,1] = sum(row == 0 & phenotype == 1)
                 fe_mat[2,2] = sum(row == 0 & phenotype == 0)

                 fe = fisher.test(fe_mat)
                 return(fe$p.value)
               }
)

Now that you have read how my data is currently formatted I would like to hear suggestions on how I can fit this data into a more complex test using sci kit or another resource. Algorithms I am interested in using are any you see fit and previously field-tested algorithms such as adaboost and forest learning algorithms. Also please ask me any clarifications or things I may have left out of this post!

Since you are using R, you probably don't want to use scikit-learn, which is for Python. However, there is a similar R library mlr ("R package to make machine learning in R easy") that provides a unified interface to all popular machine learning methods. Check their tutorial on how to get started.


text:
Should the cell sorting marker genes be excluded during clustering?
We sort different populations of blood cells using a number of fluorescent flow cytometry markers and then sequence RNA. We want to see what the transcriptome tells us about the similarity and relation between these cells. In my experience on bulk RNA-seq data, there is a very good agreement between flow cytometry and mRNA expression for the markers. It's good to remember we sort using a few markers only, while there are hundreds if not thousands of cell surface proteins.
Should we exclude genes of these sorting (CD) marker proteins when we perform PCA or other types of clustering?
The argument is that after sorting, these genes may dominate clustering results, even if the rest of the transcriptome would tell otherwise, and thus falsely confirm similarity relations inferred from sorting.
Should we at least check the influence of these genes on clustering results?

I do not think there is a simple "yes" or "no" answer here.
A good starting point would be, as you suggest, use all the genes and assess the results in the light of the marker genes and expected results. This could both serve as as good quality control as well as give you overview of all the processes happening in the cells.
Depending on the marker genes effect and biological question you may then want to remove the marker genes with potential ordering effect, or even other genes, e.g. by GO terms or pathways.
You most likely want to account for the cell cycle phases as well. And check for other technical factors. I recommend Bioconductor pipeline to get inspired https://www.bioconductor.org/help/workflows/simpleSingleCell/ when it comes to scRNA-seq analyses.


text:
How should the SAM MD tag match the CIGAR string?
I am trying to understand how the MD:Z tag is used. The following is from the SAM Optional Fields Specification, which gives an example but is not thorough.

The MD field aims to achieve SNP/indel calling without looking at the
  reference. For example, a string ‘10A5^AC6’ means from the leftmost
  reference base in the alignment, there are 10 matches followed by an A
  on the reference which is different from the aligned read base; the
  next 5 reference bases are matches followed by a 2bp deletion from the
  reference; the deleted sequence is AC; the last 6 bases are matches.
  The MD field ought to match the CIGAR string.

Suppose I have a read that I want to soft clip at both ends. If the read starts with CIGAR 100M and MD 50G49, and I want to change the CIGAR to 7S86M7S, what should the MD field become?
Is there a more complete explanation of what appears in the MD tag field?

The MD string doesn't apply to soft or hard clipped regions,. so your example read becomes 43G42. Since variant calling using something simplistic like this is only ever going to use the aligned portion, I guess there was never much reason to bother with adding in soft-clipping information. This also means that you need to parse the CIGAR string if you're using the MD string.


text:
Customizing bigWig file
I generate bigWig files using bamCoverage from deeptools, in part for my colleagues to visualize their mapped libraries in the IGV viewer.
A problem is that the displayed track name is apparently the file name, which is not convenient because some files for different libraries have the same name. The reason is that the files are organized with one directory per library. I would prefer to avoid changing the way I name and organize the files.
Is there a way to edit my bigWig files in order to have control of the track name independently from the file name?
As far as I understand, bigWig format derives from wiggle format, which have a track definition line that can contain a name and a definition attribute: https://genome.ucsc.edu/goldenpath/help/customTrack.html#TRACK.
Maybe there is a way to convert a bigWig file back to wiggle, edit the track definition, and then rebuild an updated bigWig file.

There's no equivalent to the wiggle header in bigWig (or bigBed) files, which is why UCSC uses the file name. This is actually the reason for the track line stuff that you linked to, since you can then specify a name and just point to where the bigWig (or other format) file is on the internet.
BTW, you can certainly convert your bigWig to wiggle, add the track line, and convert it back with UCSC tools (e.g., bigWigToWig and wigToBigWig). It shouldn't change the track name when you display it in UCSC (though if it does, then let me know, since I'll modify pyBigWig and libBigWig to support that then).

Using bigWigToWig, I can indeed obtain a wig file from the bigwig and add a custom track name that is displayed in IGV. Then, rebuilding a bigwig using wigToBigWig indeed looses the track name and I get the file name instead in IGV. This is for biwig files obtained directly with bamCoverage. I failed at reconverting wig to bigwig when the initial bigwig file was "manually" generated using pyBigWig. wigToBigWig complained as follows: line 1507430 of tests/tmp.wig: chromosome I has 15072434 bases, but item ends at 15072440. I guess I should post another question regarding this issue.

Thanks for the update, at least my understanding of bigWig headers was correct :)


text:
Why do BLASTn and prokka not seem to be searching the whole fasta file?
When I use blastn and prokka (I will detail exactly how I did so below) on a 2.8 million bp fasta file I get output start/end numbers that do not seem to cover the entire genome.
Starting with a .fna genome such as genome.fna I ...
1 - blastn
Searched the genome for 3135 different 28-mers using BLASTn.
makeblastdb -in genome.fna -dbtype nucl -parse_seqids -out ./output/genome

next command in python
blast_tsv = NcbiblastnCommandline(query=Q, db=DB, perc_identity=100, outfmt=6, out=(OUT))
    stdout, stderr = blast_tsv()

Q is this list of k-mers. DB is the database created
This outputs a list of search results but the important thing is that none of the start/end (columns 7 and 8) integers were greater than 100,000 yet the entire genome is 2.8 million base pairs long. I can provide the link to this file in the comments. Does this have to do with blastn stopping searching after finding one match? And if so how can I tell blastn to search for every match for each k-mer? (I'm open to using other blastn programs other than biopython)
2 - prokka
I used prokka to create gff files for the genome.
prokka-1.12/bin/prokka —setupdb
prokka-1.12/bin/prokka -kingdom Bacteria -rfam -outdir ./prokka_database/genome -force -prefix "${genome/.fna/}" genome.fna

Th gff file produced by this command only described genes with start and end positions less than 200k. I can provide the link to this file as well.

With a k-mer size of 28 it shouldn't be finding that many matches. And the prokka results are suspicious as well. Maybe you have multiple contigs (none larger than 100kb) in that file? What is the result of

grep ^'>' fasta_file | wc -l

? This would show how many contigs you have in the file.

Over at http://github.com/dib-lab/khmer we'd be happy to help you work up a script that does what you need. Just post an issue.

I searched 3135 different 28-mers and about half of them actually matched. What I'm concerned with is why they only happen in the first 100 kpb of the genome. Using grep ^'>' fasta_file | wc -l , I found 32 contigs in the original fasta file.

Link to blast results (searched 3135 unique k-mers on one genome): http://ix.io/xqN

Hi Daniel, so it looks like you have multiple contigs. 32 * 100kb is about 3.2 Mbp so many of those contigs could be less than 100kb. It looks like maybe all the accession numbers are the same, from the BLAST results. Weird. Not sure what's going on, would have to take a look at the genome file to be sure.

I see! Here is the genome file I am using. ftp://ftp.patricbrc.org/patric2/current_release/AMR_genome_sets/Staphylococcus/methicillin/Resistant/1413323.3.fna

I have tested a bunch of these genomes from that ftp directory though with the same blast command and they all search 1 contig only but it varies which contig.


text:
Time and memory efficient processing of many intervals for a bigWig file
I am trying to calculate the mappability adjusted length of introns as described by Boutz et al.
Briefly, for each intron I wish to calculate the length minus the number of bases that are non-uniquely mappable. Mappability tracks can be downloaded as bigWigs from here. The the score at each base is 1/number of mapping positions, so 1 indicates a read can be uniquely mapped at this location.
I am struggling to do this computation in either a reasonable amount of time, or a reasonable amount of memory.
First I tried importing the whole bigwig:
library(rtracklayer)
mappability <- import(BigWigFile("wgEncodeCrgMapabilityAlign50mer.bigWig"),
                      as="NumericList",
                      selection=BigWigSelection(intron_ranges))

where intron_ranges is a GRanges object with the introns (about 800,000 of them).
This used way too much memory, and soon caused by machine to fall over.
Second I tried processing one intron at a time:
mappability_file = BigWigFile("wgEncodeCrgMapabilityAlign50mer.bigWig")
effective_length <- function (gr) {
  intron_selection <- BigWigSelection(gr)
  scores <- import(mappability_file, as="NumericList", selection=intron_selection)
  non_unique <- sum(scores[[1]] < 1.0)
  eff_len = width(gr)[1] - non_unique
  return(eff_len)
}

mappability <- sapply(intron_ranges, effective_length)

This has been running for hours and shows no sign of finishing.
Is there a way to do this that uses less than, say 4GB of RAM, but finishes in say less than 30 minutes? It feels like there should be.
I'm not tied to R, just what I've tried so far. Happy with answers using binary packages, common shell tools, python or R.

4 GiB is not a lot of memory, especially when working with genomic data. In fact, 800.000 GRanges isn’t all that much (I’m working with repeat elements, I have a lot more). … What I’m trying to say: time to upgrade your hardware.

I have 32GB on my machine. Its not the Granges that takes the space, its the passing them through import.bw. I siad 4GB as I don't want to have to close everything else down to do this sort of thing (processing BigWigs must be a pretty common task). I have access to a machine with 256GB on the local cluster so it would be possible to use that, but it seems unnecessary. Genome aligners that store indexes in memory use less that than. Even STAR which store all sorts fits in 24GB.

Assuming you can make a BED file of introns, you can then use the pyBigWig module in python:
import pyBigWig
import numpy

bw = pyBigWig.open("some file.bw")
bed = open("introns.bed")
for line in bed:
    cols = line.strip().split("\t")
    vals = bw.values(cols[0], int(cols[1]), int(cols[2]), numpy=True)
    effLen = cols[2] - cols[1] - (vals < 1.0).sum()
    # Do something with this.

One can make that parallel (with the deeptools API) should it prove too slow. This will only use an appreciable amount of memory for very large introns, since there you need to store the value at each base. One could make this much more memory efficient by using bw.intervals() instead, but that'd require quite a bit more code (you'd need to determine the number of bases overlap per interval and then assume that either each base you want to query has an overlapping interval or, if not, keep track of that fact).

Well, this looks very promising! I've only come across bx-python before for accessing BigWig files before.

It should also be a bit faster than bx-python (I wrote it largely to speed up deepTools).

This solution takes about 2 hours to process the file on our system in case you are interested. Although I wasn't able to use the numpy version. Uses very little memory as far as I can tell.

@IanSudbery: Thanks for the update. I don't expect the numpy version to be appreciably faster. If you need this faster then I can make a multithreaded version using the deepTools API.

Out of interest, why is BigWig access so slow. If this were a BAM file or a tabix indexed bedgraph, the process would have taken seconds with psyam or Rsamtools or similar.

The format is optimized for using zoom levels, which you're not using in this case.


text:
Is the optional SAM NM field strictly computable from the MD and CIGAR?
From SAM Optional Fields Specification the NM field is

Edit distance to the reference, including ambiguous bases but excluding clipping

Assuming both the MD and CIGAR are present, is the edit distance simply the number of characters [A-Z] appearing in the MD field plus the number of bases inserted (xI, if any) from the CIGAR string? Are there any other complications?

Assuming both the MD and CIGAR are present and correct, then yes, you can parse both to get the edit distance (NM auxiliary tag). One big caveat to this is that there's a reason that the samtools calmd command exists, since it's historically been the case that not all aligners have output correct MD strings. It's rare for the CIGAR string to be wrong and that'd be more of a catastrophic error on the part of an aligner. For what it's worth, if the NM auxiliary is absent on a given alignment but present on others produced by the same aligner then it's fair to assume NM:i:0 for a given alignment by default (many aligners only produce NM:i:XXX if the edit distance is at least 1).


text:
Visualisation of long read RNA-Seq splicing
I have a dataset of Oxford Nanopore cDNA reads. Many of my reads are full-length or close to full-length transcripts, and I and am interested in examining alternative splicing. For this, I would like to begin by visualising my reads and comparing variants qualitatively.
I have tried visualising my data in both IGV and SeqMonk, but neither has given a satisfying result.
IGV

IGV shows some links between aligned sections, shown as fine blue lines and thick black lines. I cannot work out what the difference between these lines is, and even more confusingly, we have numerous alignments where two aligned exons come from the same sequence but are not joined by any line. This means that, to confirm or deny a spliced variant, I need to manually examine the read ID of each exon.
SeqMonk

SeqMonk gives a much cleaner visualisation, but unfortunately shows no links between aligned exons, and worse, I cannot find the read IDs. This means there is functionally (as far as I can see) no way to tell which exons come from the same read.
EDIT: Sashimi Plot
As recommended in this answer, I have tried using IGV's Sashimi Plot. Unfortunately, I believe the high error rate and relatively low coverage is causing this to create somewhat confusing output. The Sashimi Plot shows mostly nearly every junction to have just a single read supporting it.

Is there a visualisation tool which allows simple examination of splicing of full-length transcripts?

Are you sure you are mapping the reads correctly? Do they have a lot of soft clips, perhaps suggesting that they aren't fully mapped?

I couldn't say for sure, but it looks okay to me.. They were mapped with GMAP, and with read lengths in the range of 1500-4000bp we have soft clipping sitting around 20-50 for most of those reads. About half of the reads have mapping quality 40 (the others are very low, about a quarter have MAPQ 0.)

IGB gives the intended result without much hassle. It doesn't have the bells and whistles of IGV, but presents a clean and intuitive view of the individual reads.
In IGB, each read is presented as a single, contiguous line.

If you zoom in you can see the read ID for each read.

For a more holistic view, IGV's Sashimi plot can be modified to exclude spurious junctions supported by few reads. Here is a section of the above gene with default settings:

The minimum junction coverage can be set by right-clicking on the Sashimi plot itself, or application-wide under View -> Preferences -> Alignments -> Splice Junction Track Options. With minimum junction coverage 3, the noisy plot above now looks like this:

The numbers are still hard to read due to numerous overlapping almost-identical junctions, but at least the plot is now readable.

DNASTAR's software is for purchase, but high quality. GenVision Pro does genomic visualization, including Sashimi plots.
Edit: not sure why this answer is being downvoted, unless it's because the software isn't free. OP has tried IGV and SeqMonk, I mentioned an alternative he might not have heard of.
Here is a video demonstrating the use of Sashimi plots in GenVision Pro:
http://www.dnastar.com/t-support-videos.aspx?video=YJvcERoSIsg

Please downvoters, explain the reason for downvoting.

Maybe add example plot outputs?

The Sashmi Plot feature built into IGV. It gives a nice summary of the spliced transcripts and the coverage of each exon.
For example:


text:
Assembling sequence data generated by RADseq
I'm working on a project where I have to assemble sequences generated by RADseq. At the end I hope to compare two species of woodpeckers in Sri Lanka by using SNPs.
I tried to assemble it using Stacks pipeline, but it crushed on the way with:

*** Error in `ustacks': malloc(): memory corruption: 0x000055db58024be0 ***

Are there are any other tools or pipelines that I can use for assembling my data instead?

There are plenty, see a list here.
If you want a better answer, write a specific question, write what you tried, why Stacks doesn't work and what is the desired output.

@ShanikaLasanthi I tried to make your question clearer, it would be good to add the details about how you run stacks and what kind of error message you got before it crushed.

I don't think you should be assembling RADseq, you should be aligning it.

I am getting this error saying "*** Error in `ustacks': malloc(): memory corruption: 0x000055db58024be0 ***"
what am I supposed to do now?

@ShanikaLasanthi Which command did you use?

I used the code denovo_map.pl -s trim_ln7_brds_KD03S01.fastq -o ./stacks -b 1 -B trial1_radtags -t -T 15 -m 3 -M 2 -n 1.

I figured out what happened. My sequences were not in the same length,therefore it interfered with the algorithm giving this error. once I trimmed the sequences into the same length, pipeline was implemented successfully without giving any error.
Thank you all for the suggestions and the advice.

@ShanikaLasanthi it would be nice if you write here an actual answer for your question. It seems that the solution is not very intuitive and it could be very helpful for others to see the answer.

Well, first my sequences were in different lengths between 50-100bp. I got the error when I used these sequences in the pipeline. When I checked the denovo log file for the run ( you can find it in the output folder you have mentioned for the run) there was a warning saying "different sequence lengths detected, this will interfere with Stacks algorithms."
Therefore I trimmed the samples into 80bp and filtered out sequences shorter than 80bp.
After that I did not get the error message and the pipeline was implemented successfully.
And I hope this helps.


text:
How to select a cutoff for interaction confidence in STRINGdb?
I have a list of 100 genes that are called as hits in a genetic screening. I want to have a network of the interactions between the proteins of these 100 genes. I am using both STRINGdb web and its R API.
In both situations you can select the threshold of confidence of interaction (default 0.4) for displaying the network. Changing this parameter drastically changes the network that I am generating. I have read the FAQ section of STRINGdb and they recommend to choose some arbitrary number based on the number of interactions you need for you analysis. If I use the default threshold, I have few interactions and I don't know if lowering the threshold would be correct.
Therefore, my question is, is there any established threshold based on habitual user's experience used for that, beyond the recommendation of setting an arbitrary one?

I have used STRING pretty heavily, and have compared it to various other databases of protein interactions and signaling pathways.  I do feel like it has a lot of quality interaction annotations, but you have to sift through a lot of noise to get to them.  The simplest method I have found for doing this is to look at the individual scores for each interaction, and accept it if it passes one of the following tests:

Experiment Score > 0.4
Database Score > 0.9

Anything that passes one of these thresholds we consider at least an interaction of acceptable quality.  Those interactions with Experiment Scores > 0.9 are high-quality, and have been experimentally validated.  The other scores represent inaccurate methods for determining signaling events, and should be ignored.
You are not going to catch every actual protein signaling event this way, but you will at least be spared a lot of false positives.  The best way to construct an actual network of signaling events is to combine interaction records from multiple databases.

This depends on what you are trying to do and whether you value specificity over sensitivity. We can't tell you since it is entirely dependent on the biological question you want to answer.
However, I would recommend two things:

Don't use STRING. The creators of STRING made the choice to value sensitivity over all else, so they include any interaction they can get their hands on. That is sometimes fine, depending on what you want to do, but is more often a problem. Personally, I tend to avoid STRING as much as possible and instead use more curated databases like APID or IntAct. You can find a list of many PPI databases in the EBI's Psiqcuick View page. You might also be interested in my answer here which gives an example script using APID which can easily be modified to query other DBs.

Whatever you use, I recommend you filter by detection method. There are many different interaction detection methods ranging from those that only identify direct, binary interactions between two proteins (e.g. Yeast two-hybrid), those that also find proteins that are in the same complex even if there is no direct interaction between them (e.g. ChIP), to various non-experimental methods which are used to infer interactions.
Luckily, all of these have been collected into a controlled vocabulary which you can see in the Interaction Detection Method Ontology Lookup Service page. Go there, collect the interaction detection methods you consider good enough for whatever you want to do and filter your interactions using those.

For example, for a high quality human interactome I built for a project I used to work on, I was only interested in direct binary interactions, so I only kept interactions detected by the following methods:
    "MI:0008" => "array technology",
    "MI:0009" => "bacterial display",
    "MI:0010" => "beta galactosidase complementation",
    "MI:0011" => "beta lactamase complementation",
    "MI:0012" => "bioluminescence resonance energy transfer",
    "MI:0013" => "biophysical",
    "MI:0014" => "adenylate cyclase complementation",
    "MI:0016" => "circular dichroism",
    "MI:0017" => "classical fluorescence spectroscopy",
    "MI:0018" => "two hybrid",
    "MI:0020" => "transmission electron microscopy",
    "MI:0030" => "cross-linking study",
    "MI:0031" => "protein cross-linking with a bifunctional reagent",
    "MI:0034" => "display technology",
    "MI:0040" => "electron microscopy",
    "MI:0041" => "electron nuclear double resonance",
    "MI:0042" => "electron paramagnetic resonance",
    "MI:0043" => "electron resonance",
    "MI:0047" => "far western blotting",
    "MI:0048" => "filamentous phage display",
    "MI:0049" => "filter binding",
    "MI:0051" => "fluorescence technology",
    "MI:0052" => "fluorescence correlation spectroscopy",
    "MI:0053" => "fluorescence polarization spectroscopy",
    "MI:0055" => "fluorescent resonance energy transfer",
    "MI:0065" => "isothermal titration calorimetry",
    "MI:0066" => "lambda phage display",
    "MI:0073" => "mrna display",
    "MI:0081" => "peptide array",
    "MI:0084" => "phage display",
    "MI:0089" => "protein array",
    "MI:0090" => "protein complementation assay",
    "MI:0091" => "chromatography technology",
    "MI:0092" => "protein in situ array",
    "MI:0095" => "proteinchip(r) on a surface-enhanced laser desorption/ionization",
    "MI:0097" => "reverse ras recruitment system",
    "MI:0098" => "ribosome display",
    "MI:0099" => "scintillation proximity assay",
    "MI:0107" => "surface plasmon resonance",
    "MI:0108" => "t7 phage display",
    "MI:0111" => "dihydrofolate reductase reconstruction",
    "MI:0112" => "ubiquitin reconstruction",
    "MI:0114" => "x-ray crystallography",
    "MI:0115" => "yeast display",
    "MI:0226" => "ion exchange chromatography",
    "MI:0227" => "reverse phase chromatography",
    "MI:0231" => "mammalian protein protein interaction trap",
    "MI:0232" => "transcriptional complementation assay",
    "MI:0255" => "post transcriptional interference",
    "MI:0369" => "lex-a dimerization assay",
    "MI:0370" => "tox-r dimerization assay",
    "MI:0397" => "two hybrid array",
    "MI:0398" => "two hybrid pooling approach",
    "MI:0399" => "two hybrid fragment pooling approach",
    "MI:0400" => "affinity technology",
    "MI:0401" => "biochemical",
    "MI:0405" => "competition binding",
    "MI:0406" => "deacetylase assay",
    "MI:0410" => "electron tomography",
    "MI:0411" => "enzyme linked immunosorbent assay",
    "MI:0415" => "enzymatic study",
    "MI:0416" => "fluorescence microscopy",
    "MI:0419" => "gtpase assay",
    "MI:0420" => "kinase homogeneous time resolved fluorescence",
    "MI:0423" => "in-gel kinase assay",
    "MI:0424" => "protein kinase assay",
    "MI:0425" => "kinase scintillation proximity assay",
    "MI:0426" => "light microscopy",
    "MI:0428" => "imaging technique",
    "MI:0432" => "one hybrid",
    "MI:0434" => "phosphatase assay",
    "MI:0435" => "protease assay",
    "MI:0437" => "protein three hybrid",
    "MI:0440" => "saturation binding",
    "MI:0508" => "deacetylase radiometric assay",
    "MI:0509" => "phosphatase homogeneous time resolved fluorescence",
    "MI:0510" => "homogeneous time resolved fluorescence",
    "MI:0511" => "protease homogeneous time resolved fluorescence",
    "MI:0512" => "zymography",
    "MI:0513" => "collagen film assay",
    "MI:0514" => "in gel phosphatase assay",
    "MI:0515" => "methyltransferase assay",
    "MI:0516" => "methyltransferase radiometric assay",
    "MI:0655" => "lambda repressor two hybrid",
    "MI:0657" => "systematic evolution of ligands by exponential enrichment",
    "MI:0678" => "antibody array",
    "MI:0695" => "sandwich immunoassay",
    "MI:0696" => "polymerase assay",
    "MI:0726" => "reverse two hybrid",
    "MI:0727" => "lexa b52 complementation",
    "MI:0728" => "gal4 vp16 complementation",
    "MI:0809" => "bimolecular fluorescence complementation",
    "MI:0813" => "proximity enzyme linked immunosorbent assay",
    "MI:0824" => "x-ray powder diffraction",
    "MI:0825" => "x-ray fiber diffraction",
    "MI:0827" => "x-ray tomography",
    "MI:0841" => "phosphotransferase assay",
    "MI:0870" => "demethylase assay",
    "MI:0872" => "atomic force microscopy",
    "MI:0879" => "nucleoside triphosphatase assay",
    "MI:0880" => "atpase assay",
    "MI:0887" => "histone acetylase assay",
    "MI:0889" => "acetylase assay",
    "MI:0892" => "solid phase assay",
    "MI:0894" => "electron diffraction",
    "MI:0895" => "protein kinase A complementation",
    "MI:0899" => "p3 filamentous phage display",
    "MI:0900" => "p8 filamentous phage display",
    "MI:0905" => "amplified luminescent proximity homogeneous assay",
    "MI:0916" => "lexa vp16 complementation",
    "MI:0920" => "ribonuclease assay",
    "MI:0921" => "surface plasmon resonance array",
    "MI:0946" => "ping",
    "MI:0947" => "bead aggregation assay",
    "MI:0949" => "gdp/gtp exchange assay",
    "MI:0953" => "polymerization",
    "MI:0968" => "biosensor",
    "MI:0969" => "bio-layer interferometry",
    "MI:0972" => "phosphopantetheinylase assay",
    "MI:0976" => "total internal reflection fluorescence spectroscopy",
    "MI:0979" => "oxidoreductase assay",
    "MI:0984" => "deaminase assay",
    "MI:0989" => "amidase assay",
    "MI:0990" => "cleavage assay",
    "MI:0991" => "lipid cleavage assay",
    "MI:0992" => "defarnesylase assay",
    "MI:0993" => "degeranylase assay",
    "MI:0994" => "demyristoylase assay",
    "MI:0995" => "depalmitoylase assay",
    "MI:0996" => "deformylase assay",
    "MI:0997" => "ubiquitinase assay",
    "MI:0998" => "deubiquitinase assay",
    "MI:0999" => "formylase assay",
    "MI:1000" => "hydroxylase assay",
    "MI:1001" => "lipidase assay",
    "MI:1002" => "myristoylase assay",
    "MI:1003" => "geranylgeranylase assay",
    "MI:1004" => "palmitoylase assay",
    "MI:1005" => "adp ribosylase assay",
    "MI:1006" => "deglycosylase assay",
    "MI:1007" => "glycosylase assay",
    "MI:1008" => "sumoylase assay",
    "MI:1009" => "desumoylase assay",
    "MI:1010" => "neddylase assay",
    "MI:1011" => "deneddylase assay",
    "MI:1016" => "fluorescence recovery after photobleaching",
    "MI:1019" => "protein phosphatase assay",
    "MI:1024" => "scanning electron microscopy",
    "MI:1026" => "diphtamidase assay",
    "MI:1030" => "excimer fluorescence",
    "MI:1031" => "protein folding/unfolding",
    "MI:1036" => "nucleotide exchange assay",
    "MI:1037" => "Split renilla luciferase complementation",
    "MI:1038" => "silicon nanowire field-effect transistor",
    "MI:1087" => "monoclonal antibody blockade",
    "MI:1088" => "phenotype-based detection assay",
    "MI:1089" => "nuclear translocation assay",
    "MI:1111" => "two hybrid bait or prey pooling approach",
    "MI:1112" => "two hybrid prey pooling approach",
    "MI:1113" => "two hybrid bait and prey pooling approach",
    "MI:1137" => "carboxylation assay",
    "MI:1138" => "decarboxylation assay",
    "MI:1142" => "aminoacylation assay",
    "MI:1145" => "phospholipase assay",
    "MI:1147" => "ampylation assay"

If you also want those that might not be direct binary but simply in the same complex, include:
      "MI:0019" => "coimmunoprecipitation",
      "MI:0006" => "anti bait coimmunoprecipitation",
      "MI:0007" => "anti tag coimmunoprecipitation",
      "MI:0858" => "immunodepleted coimmunoprecipitation",
      "MI:0096" => "pull down",
      "MI:0963" => "interactome parallel affinity capture",
      "MI:0676" => "tandem affinity purification"

Thank you for the answer, Terdon. Could you tell me how to retrieve MI's from String for a particular network?

@SashkoLykhenko please ask a new question about that. It is hard to answer in the comments.


text:
How do I find the most similar sequences from a large set of short sequences?
I've been asked to evaluate whether the MinION will be sufficient to distinguish between 20bp CRISPR guide RNAs from the GeCKO v2 set. I know that this set has some exactly identical sequences, which wouldn't be distinguishable regardless of the sequencing method used, so I would like to know whether any two non-indentical sequences differ by less than 15% of the 20bp (i.e. only 3 bases different).
There are around 130k 20bp sequences in the dataset, and I'm pretty sure that an all-vs-all approach is not feasible, but perhaps there's some trick that I'm missing.
What I thought of doing was looking at dimer and trimer/codon counts to hunt down sequences that were most similar. I found about 100 sequences that shared the same trimer count signature (but were not identical), and within this group found [only] one pair with 85% identity; all the remainder were 70% or less.
... but I'm having trouble ignoring the doubt I have that these two sequences are the most similar. It might be possible that there are sequences in the dataset that don't have identical trimer count signatures, but differ by less than three bases. I'll give a very obvious example (which isn't in the dataset):
1: AAAAA AAAAA AAAAA AAAAA
2: AAAAC AAAAA AAAAA AAAAA

The trimer counts for these two sequences are different (Sequence 1 -- AAA: 18; Sequence 2 -- AAA: 15; AAC: 1; ACA: 1; CAA: 1), and yet the sequences are only 1 base different (which might be a problem for the MinION). How can I work out if any such sequence pairs exist in the dataset?

130k * 20bp is a small data set. At this scale, SSE2 Smith-Waterman may work well:
git clone https://github.com/attractivechaos/klib
cd klib && gcc -O2 -D_KSW_MAIN ksw.c -o ksw -lz
./ksw -a1 -b1 -q1 -r1 -t14 20bp.fa 20bp.fa > out.tsv

I simulated 130k 20bp reads from E. coli K-12. The last command line takes about ~2 CPU hours based on the partial output.
It is possible to make this much faster by using 5bp seeds. You can also try ssearch and fasta from the fasta aligner package. They are slower, probably because they are not optimized for such input. Computing Hamming distance with SSE2 will be much faster, but is probably of little use as it does not allow gaps. There are also Gene Myers' edit-distance based algorithms, which can be faster than SSE2-SW.
Bowtie2 (as well as most short-read mappers) won't work well at the default setting. Bowtie2 uses long exact seeds. It will miss many 1-mismatch hits, let alone 3-mismatch ones. The number you get from bowtie2 is an underestimate. You might be able to tune bowtie2, but to find 3-mismatch hits, you have to make it a lot slower.

For this application, you could probably also do something like calculate the Hamming distances between all of the strings in an all-vs-all approach (it should not take too long or too much overhead). You could use something like the Hamming distance tools in Julia. Here is an example of what I mean (using Julia):
using StringDistances

k = ["AATTGGCC", "AATTGGCA", "AATTCCGG"]

for s in k
    for y in k
        print(s, " compared to ", y, ":  ", compare(Hamming(), s, y), "\n")
    end
end

The output looks like this:
AATTGGCC compared to AATTGGCC:  1.0
AATTGGCC compared to AATTGGCA:  0.875
AATTGGCC compared to AATTCCGG:  0.5
AATTGGCA compared to AATTGGCC:  0.875
AATTGGCA compared to AATTGGCA:  1.0
AATTGGCA compared to AATTCCGG:  0.5
AATTCCGG compared to AATTGGCC:  0.5
AATTCCGG compared to AATTGGCA:  0.5
AATTCCGG compared to AATTCCGG:  1.0

This is just a simple small example but the core is there. Should not be too bad scaling it up. Hope this helps and provides another helpful option. :)

I'm pretty sure that an all-vs-all approach is not feasible, but perhaps there's some trick that I'm missing.

The "trick" is that I hadn't actually tried this, and Bowtie2 happens to be quite good at doing this. When there are other similar matches, the MAPQ score that Bowtie2 produces are reduced, so all that is needed is to identify those sequences with less than a MAPQ of 42. I've found 2810 such sequences that share some similarity to others (2416 if I exclude reverse-complement matches):
$ bowtie2-build Mouse_GeCKOv2_Library.fasta Mouse_GeCKOv2_Library.fasta
$ bowtie2 -f -x Mouse_GeCKOv2_Library.fasta -U Mouse_GeCKOv2_Library.fasta |
   samtools view | awk '{if($5 != 42){print ">"$1"\n"$10}}' > similar_seqs.fasta
125795 reads; of these:
  125795 (100.00%) were unpaired; of these:
    1 (0.00%) aligned 0 times
    122982 (97.76%) aligned exactly 1 time
    2812 (2.24%) aligned >1 times
100.00% overall alignment rate


text:
Aligning many long sequences
I'm faced with having to align many (some 100s) bacterial genomes, where the genome length is in the millions. Obviously, this is beyond normal alignment techniques and it's unclear to me what the best practice is for such circumstances:

conventional alignment on a very powerful computer with lots of memory
break up the genome into smaller fragments and align them individually
some exotic different procedure

What possible avenues of attack are there?
(I've attempted to use Mafft and Clustal with little success)

A standard MSA tool like Clustal will assume that there are no rearrangements between sequences. Unless all of your sequences are completely colinear, you'll need a dedicated genome alignment method like MUMmer.

Adding some clarifications: there are no large scale rearrangements (otherwise I wouldn't try this). Most solutions are pairwise or to a reference, which isn't what I want to do here.

I wouldn't expect Clustal to work but what went wrong with MAFFT? I'm pretty sure that's been used for whole genome alignments on vertebrate data, so I would be surprised if it failed for bacteria.

There are some tips here to run mafft on a relatively small number of long sequences.

Also, could you explain why you are doing this? I ask on the off chance that this is an XY problem and whole genome alignments aren't actually the best way to achieve whatever your end goal is.

Answers to above questions: The "quick" Mafft alignment was visibly poor, while the higher quality one never completed. I'm trying to understand variation across a large set of bacteria gathered from a variety of sites and circumstances and see what may influence phenotype.

Please [edit] your answer to provide clarifications. Comments are hard to read, easy to miss and can be deleted without warning. That said, if you are looking for variation and things that influence phenotype, why would you go for whole genome alignments? That doesn't seem likely to be of much help. You might want to post a new question explaining what you are trying to do and ask for help designing your approach. For example, I would most likely first identify the genes of your species (either using homology based approaches or de-novo gene prediction) and align those.

As for the quick mafft, yes, I'm not surprised it was bad, it isn't designed for WG alignments. The full one will of course take a long time. This is a big job. How long did you run it for and on what hardware? Multiple WG alignments aren't the sort of thing you can do on your laptop but, luckily, they are very rarely necessary.

Have you considered using LASTZ (http://www.bx.psu.edu/~rsharris/lastz/ )? It is meant to align chromosome sized entities?

Whole genome aliment can be done using Progressive Mauve, LAST or Mummer. For bacteria I used Mauve since it has also very nice visualisation engine. A very new tool is Minimap2, a super fast mapper that supposed beside read mapping be able to handle reference vs reference. However, I do not know how performance of it compares to the tools mentioned above.
If you are interested in a rough idea of the shared genome regions, you can use bevel. Bevel is not really an aligner, it is more like a dot-plot, but it is super fast (even for mammalian sized genomes).


text:
Get RefSeq accession numbers with versions
Google searching for NM_002084 gives the following result:
NM_002084.4
This, I assume, is the latest version v4, hence the .4 suffix.
Searching for previous versions I get the following results, along with notes saying it was updated or removed.
NM_002084.3

This sequence has been updated. See current version.

NM_002084.2
NM_002084.1

Record removed. This record was replaced or removed.

Using biomaRt I can get the latest(?) version as follows:
library("biomaRt")

# define db
ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")

# get refseqs
getBM(attributes = c('refseq_mrna',
                     'chromosome_name',
                     'transcript_start',
                     'transcript_end',
                     'strand'),
      filters = c('refseq_mrna'),
      values = list(refseq_mrna = "NM_002084"),
      mart = ensembl)
#  refseq_mrna chromosome_name transcript_start transcript_end strand
#1   NM_002084               5        151020438      151028992      1

But querying for specific versions gives nothing:
getBM(attributes = c('refseq_mrna',
                     'chromosome_name',
                     'transcript_start',
                     'transcript_end',
                     'strand'),
      filters = c('refseq_mrna'),
      values = list(refseq_mrna = c("NM_002084.1", "NM_002084.2", "NM_002084.3", "NM_002084.4")),
      mart = ensembl)
# [1] refseq_mrna      chromosome_name  transcript_start transcript_end   strand
# <0 rows> (or 0-length row.names)

Question: How can I get all versions (preferably using R)?

Can you please explain why it's important to get all versions? The most recent version is typically an update of the sequence to correct errors in previous versions.

@gringer True having the latest version should suffice, but I find it strange that it is not possible to get all versions automatically, when data is available online through manual search.

I don't believe this is possible using biomaRt, nor using AnnotationHub.
I have two suggestions, neither of them very satisfactory. First, you can specify an Ensembl archive for biomaRt, for example:
mart72.hs <- useMart("ENSEMBL_MART_ENSEMBL", "hsapiens_gene_ensembl",
                      host = "jun2013.archive.ensembl.org")

Of course, that requires that you have some idea of the date for each accession version and that the archives span that date - so not especially useful.
The other option is to access EUtils using e.g. rentrez, which does allow search by version number:
library(rentrez)
es <- entrez_search("nuccore", "NM_002084.1 NM_002084.2")
es$ids

[1] "4504104" "6006000"

So knowing the accession, you could simply append 1, 2, 3... to it, run the search and see if UIDs come back, then get them using entrez_fetch.
EDIT: here's a quick and dirty function which takes an accession as query, appends version = 1, fetches the ID, then increments the version and repeats until no more results are returned. It is not well-tested!
getVersions <- function(accession) {
  require(rentrez)
  ids <- character()
  version <- 1
  repeat({
    es <- entrez_search("nuccore", paste0(accession, ".", version))
    if(length(es$ids) == 0) {
      break
    }
    ids[version] <- es$ids
    version <- version + 1
  })
  ids
}

Example:
getVersions("NM_002084")
[1] "4504104"    "6006000"    "89903006"   "1048339180"

rentrez library looks like a good solution, could you please expand with code the last sentence?

Edited with some extra code.

Are the ID's returned GI numbers? If yes, then they are not going to be very useful.

Well, that depends on what you want to do with the versions. If you want chromosome and coordinates as in the original question: that's always going to be difficult with old versions from old builds. If you want to retrieve the RefSeq records to compare versions, then UIDs are just fine.