text:
Minion channel ID's from Albacore
The latest version of Albacore from Oxford Nanopore Technologies calls bases from raw fast5 files. A useful piece of output is the sequence_summary.txt, which is a big tab-delimited file with information on each read.
One of the columns in there is the channel, which references the channel number from 1 to 512 on the flowcell.
These channels obviously have a physical layout on the flowcell. What is the layout? I.e. how many channels are there per row and per column, and how do the channel numbers map to the rows and columns?
My best guess is that it's a 16*32 array, but I don't know for sure, and even if I did I wouldn't know whether the channels map across the shorter or the longer axis in order (or maybe they are not in order at all, but that would be silly).

thanks to the selected answer, we have it figured out. In case anyone needs an implementation in R, here's one: https://gist.github.com/roblanf/df47b9748c3aae00809cc675aca79989

There are actually 2048 usable sequencing wells, hexagonally packed with four wells connected to the same sequencing sensor/channel via a multiplex (mux) selector. The combination of the mux and the channel number determines the physical location of the well. Unfortunately, the association between channel number and physical location is not obvious, and within each channel the muxes are also not in an obvious order, i.e. [3,4,1,2,2,1,4,3] for two adjacent channels.
Here's an "official" text description of the layout as it was for R7; it hasn't changed for the most recent flow cells (R9.5):

Channels 1-64 occur at the top of the chip (the other low channel numbers are at the bottom)
Channels order down the chip: 1-64, 449-512, 385-448, 321-384, 257-320, 193-256, 129-192, 65-128
Muxes run from left to right in the order: 3, 4, 1, 2, 2, 1, 4, 3

To confirm this layout, I pulled the channel/position lookup table definitions out of MinKNOW, and it seems to be the same as this:

Here's the code to produce this image (I haven't adjusted it to have a hexagonal structure):
#!/usr/bin/Rscript
library(rjson);

layout.fname <- "/opt/ONT/MinKNOW/Client/resources/app/resources/channels.json";

data.list <- fromJSON(file=layout.fname);

png("MinION_FC_Layout.png", width=3072, height=1536, pointsize=24);
par(mar=c(4,4,0.5,0.5), lwd=3);
data.mat <- matrix(unlist(data.list), ncol=2, byrow=TRUE);
plot(data.mat, pch=21,
     bg=(colorRampPalette(c("white","grey"))(2048))[1:2048],
     col=c("red","yellow","green","blue"),
     xlab="locX", ylab="locY", cex=2,
     xlim=c(7,76), ylim=c(0,33));
text(x=data.mat, labels=rep(1:512, each=4), cex=0.4);
legend("top",legend=c(1,2,3,4), fill=c("red","yellow","green","blue"),
       inset=0.05, title="MUX");
text(x=6, y=data.mat[data.mat[,1] == 8,2],
     labels = ceiling(which(data.mat[,1] == 8) / 4));
text(x=77, y=data.mat[data.mat[,1] == 75,2],
     labels = ceiling(which(data.mat[,1] == 75) / 4));
arrows(x0=c(8,75), x1=c(39,44), y0=33, length=0.2);
invisible(dev.off());

an interesting thing to note is that this answer is different from @wouter-de-coster answer that comes from poRe. As far as I can tell, this is the correct answer.

I think the layout in poRe was based on the numbering in the original MinKNOW GUI (which was also incorrect).

for those interested I built an annoying by-hand implementation of the mapping to use for plotting in R here: https://gist.github.com/roblanf/df47b9748c3aae00809cc675aca79989

poRe has a show.layout() function which shows you the 32*16 grid on which the channels are arranged. EDIT: THIS APPEARS NOT TO BE THE CORRECT PHYSICAL LAYOUT.
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
 [1,]  125  121  117  113  109  105  101   97   93    89    85    81    77
 [2,]  126  122  118  114  110  106  102   98   94    90    86    82    78
 [3,]  127  123  119  115  111  107  103   99   95    91    87    83    79
 [4,]  128  124  120  116  112  108  104  100   96    92    88    84    80
 [5,]  253  249  245  241  237  233  229  225  221   217   213   209   205
 [6,]  254  250  246  242  238  234  230  226  222   218   214   210   206
 [7,]  255  251  247  243  239  235  231  227  223   219   215   211   207
 [8,]  256  252  248  244  240  236  232  228  224   220   216   212   208
 [9,]  381  377  373  369  365  361  357  353  349   345   341   337   333
[10,]  382  378  374  370  366  362  358  354  350   346   342   338   334
[11,]  383  379  375  371  367  363  359  355  351   347   343   339   335
[12,]  384  380  376  372  368  364  360  356  352   348   344   340   336
[13,]  509  505  501  497  493  489  485  481  477   473   469   465   461
[14,]  510  506  502  498  494  490  486  482  478   474   470   466   462
[15,]  511  507  503  499  495  491  487  483  479   475   471   467   463
[16,]  512  508  504  500  496  492  488  484  480   476   472   468   464
      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
 [1,]    73    69    65    61    57    53    49    45    41    37    33    29
 [2,]    74    70    66    62    58    54    50    46    42    38    34    30
 [3,]    75    71    67    63    59    55    51    47    43    39    35    31
 [4,]    76    72    68    64    60    56    52    48    44    40    36    32
 [5,]   201   197   193   189   185   181   177   173   169   165   161   157
 [6,]   202   198   194   190   186   182   178   174   170   166   162   158
 [7,]   203   199   195   191   187   183   179   175   171   167   163   159
 [8,]   204   200   196   192   188   184   180   176   172   168   164   160
 [9,]   329   325   321   317   313   309   305   301   297   293   289   285
[10,]   330   326   322   318   314   310   306   302   298   294   290   286
[11,]   331   327   323   319   315   311   307   303   299   295   291   287
[12,]   332   328   324   320   316   312   308   304   300   296   292   288
[13,]   457   453   449   445   441   437   433   429   425   421   417   413
[14,]   458   454   450   446   442   438   434   430   426   422   418   414
[15,]   459   455   451   447   443   439   435   431   427   423   419   415
[16,]   460   456   452   448   444   440   436   432   428   424   420   416
      [,26] [,27] [,28] [,29] [,30] [,31] [,32]
 [1,]    25    21    17    13     9     5     1
 [2,]    26    22    18    14    10     6     2
 [3,]    27    23    19    15    11     7     3
 [4,]    28    24    20    16    12     8     4
 [5,]   153   149   145   141   137   133   129
 [6,]   154   150   146   142   138   134   130
 [7,]   155   151   147   143   139   135   131
 [8,]   156   152   148   144   140   136   132
 [9,]   281   277   273   269   265   261   257
[10,]   282   278   274   270   266   262   258
[11,]   283   279   275   271   267   263   259
[12,]   284   280   276   272   268   264   260
[13,]   409   405   401   397   393   389   385
[14,]   410   406   402   398   394   390   386
[15,]   411   407   403   399   395   391   387
[16,]   412   408   404   400   396   392   388

this solution (while I can see that it should be thought to be correct) appears to be at odds with the current information (in the selected answer). I wonder if poRe uses a layout from an older flowcell?

Hmm good catch.

for those interested, I wrote a quite compact Python function to generate the (correct!) layout in a numpy array: https://gist.github.com/wdecoster/1a86d4ba1039afebb250352fc4e502d1


text:
What is the etymology of "Entrez ID"?
Since I have seen NCBI gene names were called "Entrez ID" for the first time, I am wondering where that comes from. Such a weird name!
Does anybody know where that originates?
My hypothesis is: in French, "entrez" can be translated to "please enter" (it is to enter in the imperative mood), so I'm guessing that might come from a past website/database in French where Entrez ID as in "Please enter the ID" was written in the search bar...
But I could not find any reference... I cannot be the only one to ask himself this!

A Brief History of NCBI’s Formation and Growth only mentions Entrez starting as a search and retrieval system in CD form, and the Entrez wikipedia page doesn't elaborate on the name, so given this information I'd assume it's just a name given to the original database interface. Someone who was involved back then might know more though!

The Entrez Wikipedia page elaborates somewhat on the name. "The name "Entrez" (a greeting meaning "Come in!" in French) was chosen to reflect the spirit of welcoming the public to search the content available from the NLM." This is mentioned also in a Google Books search result. But "citation needed".

Update
As pointed out by neilfws:
According to Wikipedia and Practical bioinformatics by Michael Agostino: The name "Entrez" means "Come in!" in French (pronounced “on-tray”); it was chosen to reflect the spirit of welcoming the public to search the content.
Old (probably incorrect) theory
The first version of Entrez database was distributed by NCBI in 1991 (on CD-ROM). At that time, it consisted of nucleotide sequences, protein sequences and associated citations and abstracts from MEDLINE.
So even at the early stage it contained objects (entries) of very different types (sequences and abstracts). It is hard to find a word which will accurately and briefly describe all the types of data.
While I don't know the correct answer and this is purely a speculation, my guess is that Entrez is an intentionally misspelled word "entries"; which accurately describes all the datatypes in the database.
Using a misspelled word as a name is not uncommon, e.g. Google comes from misspelled "googol" and Flickr is a misspelled version of "flicker".

Haa, so my hypothesis was close to being correct! (French is my native language, this is why it surprised me and why I'm never sure about how to pronounce it, French or English way). Thanks a lot @neilfws and Iakov Davydov

source: Stack Exchange

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

authors: АЛЕКСЕЙ ТРОФИМОВ

sort: votes


text:
How to access the patient-derived xenografts (PDXs) repository?
I saw this nature news item, it sounds that Cellminer is obsolete, is it right?
What are the new tools to analyze the "new cell lines"? Where's the PDX repository?
What about the European initiative that the article refers to?
I would be happy to see correlation between drugs and genes and similar capabilities as we can see on Cellminer.

I already see for links when I try to search. Take a look at them below. Also if you can be precise about kind of data are you trying to find, it would be better
https://pdmr.cancer.gov/
http://cdt.northwestern.edu/news/patient-derived-xenograft-repository-now-available-researchers
http://www.epo-berlin.com/epo-tumor-models-xenografts.html
http://www.proxe.org/
http://data-analysis.charite.de/care/
http://www.oasis-genomics.org/
http://www.cbioportal.org/
Hope one of these above will serve the purpose for you.
Update: added some more links that might serve. As for the Cellminer the last update was 2012.

not sure if you will get something exactly like cellminer but take a look at those link and download the files and see to what extent you can use them. Cellminer is still not retired but probably it is not being updated. I could not find any log in their database so probably thinking they are not updating it now.

also take a look at this http://data-analysis.charite.de/care/ . to be honest cellminer update was not done since 2012, that is what I understand. Also take a look at OASIS and cbioportal as well

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


text:
What are the ways to keep track of branches in the analysis?
I'm going through an RNA-seq pipeline in R/Bioconductor and want to try multiple parameters at subsequent steps, for example, running clustering with different settings, running RegressOut or not on unwanted effects etc. That's a lot of "versions", even if I don't do combinations of these steps.
How can I keep track of this, and my conclusions?
Not necessarily want to save the results.

Save the different scripts with git (seems overkill)
Make notes in the script itself

The main purpose of git is to version code, which usually means sequential improvement of the codebase. While it is possible to use branches for multiple variants of the software, permanent branches are traditionally used for gradual integration of new features (i.e. dev/testing/master branches). Supporting multiple independent branches requires some investment, i.e. distributing common changes among branches via merge or cherry-pick. This is hard to manage when you have more than two-three branches.
If you compare different methods of analysis, you probably want to compare the results between methods. Having them on separate branches makes it hard.
In my opinion you should integrate all methods of analysis into the master branch. To avoid copy & paste, it is better to put common code in a library or an independent script. You can also specify a method as a run-time parameter of your pipeline, and create a meta-script which will execute all methods of interests.
Once you performed benchmarking, you shouldn't remove unused methods from you master branch. Having them is important for reproducible research, and your scripts could be used in the future for new datasets.

"You can also specify a method as a run-time parameter of your pipeline, and create a meta-script which will execute all methods of interests." This is the approach I use. Allows maximum flexibility for code reuse and keeps a straightforward record of everything that I tried.

Don't forget to store the output of the code (with the sessionInfo) if you run the meta-script.R

Save the different scripts with git (seems overkill)

Whoa. I did an actual double take when reading this:1 it’s the opposite of overkill.
Version controlling your scripts (using Git or something similar) is the absolute minimum, and should become completely automatic. For every new project I begin, one of the very first steps is to issue the git init command, and to set up a remote repository (on Github).
To keep track of different analyses, I use a combination of the following approaches:

Write reusable functions/scripts and parametrise. The parameters are kept either inside the script itself (that then calls the relevant function repeatedly), or in a declarative workflow.2
Document the alternative analysis approaches; once again, this could be a declarative workflow with different rules for alternative analyses, or a set of notebooks (via R Markdown).
Have different Git branches for mutually exclusive approaches. At the end of the analysis one of these branches gets merged into the main branch, and published. If I want to publish several analysis approaches, I merge all these branches into the main branch, and use approaches (1) or (2) enable them simultaneously.

In fact, I recommend creating a declarative workflow script for every analysis; I have found that this is the most practical, self-documenting way to run an analysis. It most closely resembles a wet-lab lab notebook. The advantage over a single R Markdown document is that rerunning just parts of the analysis can be completely automated, and dependencies in the workflow are apparent from the dependencies of the workflow rules. This is much harder in R Markdown.
Some time ago I create an example analysis workflow to show how this can be structured.
Regarding your other point:

Make notes in the script itself

“Notes” are a dangerous beasts: documentation is important, but experience shows that it’s sometimes very hard to keep documentation synchronised with the code. There is no mechanism that ensures that documentation and code actually agree. Differences between presumed analysis and actually executed analysis can become very problematic.
People therefore prefer using self-documenting code as much as possible; that is: writing code so that its meaning becomes immediately clear, without comments, even to somebody who hasn’t worked on the code before. Doing this well is hard and takes practice, but improves overall code quality. Once again, using a declarative workflow helps here because the dependencies between the rules are self-documenting the kind of analysis that was performed.
Jeff Atwood has written two seminal essays on this subject:

Coding without comments, and
Code tells you how, comments tell you why.

They are two of the best pieces of advice on programming that I could give.

1 To emphasise: take a look at the edit history of this answer.
2 I recommend ‘targets’ (for R), Nextflow or Snakemake. I used to use (and recommend) GNU make, but I no longer believe this to be an adequate tool.

I appreciate the edit, but it does distort the meaning of my answer: The outrage implicit in the “What????!” is entirely intentional. I’ll change it to something else since people seem to disapprove.

+1 for the emphasis on makefiles. Makefile and version control make an awesome combo for reproducibility!

Maybe I should have stated that I meant branching. That seemed to be a lot of overhead for running an R function with a different parameter setting. I'd like to accept your answer too.. but can accept only one.

The question was clearly about storing scripts in separate branches, not just general git usage.

@BradyGilg Yes, and so is my answer — in particular, the middle part.

I quite agree with this answer by @Konrad Rudolph.
I would like to emphasis that using parametrization for your scripts is what will help you avoid multiplying branches in git. So yes, use git, but you don't necessarily have to go "overkill" creating lots of branches.
Then, I would command these scripts from a workflow managment tool that will somehow take care of generating the branches of your analyses. If you use Snakemake, the various options taken along the path from your data to your results will be represented by the wildcards system, and this will be visible in the structure of your folders and file names, due to the fact that Snakemake works by inferring what should be done to produce a file based on its name.
This of course is not an excuse for not using other documentation approaches: Comments in the Snakefile and in the scripts, README files explaining how the workflow was run, using what configuration file. Put your scripts, the Snakefile, its configuration files and the README files under version control and document again using commit messages.

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


text:
GRCh38 vcf file with common cancer mutations
Is there a vcf file on the GRCh38 assembly with common cancer mutations I can download somewhere? Maybe from one of the big international cancer genomics consortia?
By common, I mean whichever mutations have been found recurrent in different types of cancer.

How would you define this? Any variant found in >1 cancer type? In how many patients? Do you need to exclude variants found in only one type of cancer then? Should these only be causative variants? What if a variant has been found in three different cancer patients but has absolutely no connection to their cancer? And yes, that is not only possible, but almost certain to happen since many variants are very common in the population. Would you want to limit this by variant frequency? Please [edit] your question and narrow it down.

@terdon It sounds like you are assuming he is doing population genetics, which may or may not be the case. In cancer genetics, they usually only care about somatic variants, so that automatically excludes "variants are very common in the population". I guess what the question really needs is a clarification regarding germline versus somatic mutations.

Your question doesn't give enough information for a specific answer but this should do for a start.

Get the VCF file describing all variants in Clinvar from NCBI:
wget ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/clinvar.vcf.gz

Extract any variants whose VCF line contains the word "cancer":
zgrep -iE '^#|CLNDBN=[^;]*cancer' clinvar.vcf.gz > cancer.vcf

Now, you should probably filter that further based on variant frequency, type of cancer, etc etc but that should be a good start.

If you are looking for cancer mutations, the primary resource is COSMIC and they provide GRCh38 VCFs. The download page is here: http://cancer.sanger.ac.uk/cosmic/download
It'll be up to you how you define "common", but the VCF includes a lot of information you can use.

The trouble with looking for cancer-associated variants is that it can be difficult to tease out spurious effects (e.g. ethnicity) from causative variants. If you're interested in what genes are implicated in most types of cancer, the annotations for the NanoString panels are quite good:

cancer pathways
cancer progression

The "Support documents" section of those pages has a link to an Excel file with gene lists and probe sequences. While these lists don't give information about variants that disrupt function, I'm struggling a bit to think through why it would be useful in advance to know that, unless the aim is to induce cancer through CRISPR or something like that.


text:
How to find GEO data sets using Drug Bank ID and bioDBnet
I want to find all experiments in GEO that are associated with a drug (for example tolvaptan). Is there any quick and scalable way to to this? I want to query more than 100 drugs.
I tried to use bioDBnet to map Drug Bank ID to look up data sets, but I don't know what output format should I select.

I'd say the quickest way is simply to search GEO datasets at the NCBI with the drug name. Did you have some kind of programmatic pipeline in mind?

there are more than 100 drugs, so a faster method would help.

Then I'd suggest programmatic search of GEO using the EUtils API. Happy to post an answer based on that if it would help? (will have to wait a few hours though)

@hhoomn Welcome to our community. I tried to add the scalability requirement to the question. Feel free to improve it. Generally we are trying to keep all questions precise, therefore do not hesitate to add further specifications directly into question instead of comments.

It would also be useful to see some example input (presumably a list of DrugBank IDs?) and some idea of the desired output (GEO dataset IDs, or summaries?)

You could take a look at the rentrez package for R. You can use the entrez_search function to get the GEO IDs for the studies related to your query, and then the entrez_summary function for the results. So, for example:
x <- rentrez::entrez_search(db="gds", term=paste0("tolvaptan AND gse[ETYP]"),
                            use_history = T)
res <- rentrez::entrez_summary(db="gds", id=x$ids)
str(res)

Where res will be a list containing objects of the esummary class. Then, you can use the following line to get the GSE accession numbers related to tolvaptan:
lapply(res, `[[`, "gse")

So you could put this function in a loop replacing "tolvaptan" with your drugs:
drugs <- c("drugA", "drugB", "drugC")
for(drug in drugs){
    x <- rentrez::entrez_search(db="gds", term=paste0(drug, " AND gse[ETYP]"))
    res <- rentrez::entrez_summary(db="gds", id=x$ids)
    lapply(res, `[[`, "gse")
}

Take a look at the rentrez documentation for more details.

One of the problems I've found with GEO and even ArrayExpress is that it seems like there are tons of false positives that come up in the search. Additionally, there may be a ton of results that you miss because you didn't include the synonym of the drug.
In order to make sure I have all the information relating to name, etc, you can write a script which downloads the result of a search in PubChem (this can be done with drugbank as well I think) and then parse the search result, extracting important information such as synonyms of drug name.
Finally, take those synonyms and use them in an Entrez (EUtils) search for the experiments. Then try to filter the results to keep only what you're actually looking for. One bad way to do this is just by making sure the experiment summary contains the term your looking for. Another way to potentially do this, for example if there is a large experiment which tests serveral different factors is to programmatically search through associated files to ensure they contain a desired keyword.
You'll get different experiment types and so the files associated with the results will depend on the type of experiment in the search result.
When you get have a list of IDs of the experiments you want, you can feed them through to GEOquery in Bioconductor.

I would recommend using GEOmetadb from Bioconductor. GEOmetadb is a simple wrapper that will download an up-to-date SQLite database with the metadata for all entities in GEO. You can then use your favorite SQL client to query the database quickly, since it's offline, refining your search as you go.

Welcome to the site. Although the recommendation is useful how does GEOmetadb could be used to download GEO data associated with a drug?

GEOmetadb will allow you to do fast full text search searches so you could search all GSEs for specific drug names or drug ids.


text:
Techniques for analyzing and quantifying sample bleed through in genotyping with Illumina
I am looking at the presence of viral genotypes within individual samples within an assay. Often times there is a sample whose read counts are firing off the charts and this sample tends to "bleed through" to the other samples.
I have recently tried confronting this problem with chi square analysis of conditional probability. For instance let's say I observe with 10 sample viral genotypes A and B. I then calculate the expected number of samples having both genotypes A and B from the number of A positive samples and number of B positive samples.
So far this is all I have. I do not know if anyone else has any ideas on how to determine and quantify bleed through in viral genotyping.
Slightly more information about this situation: these samples were run on Illumina HiSeq, we are planning on doing another HiSeq run with even more samples shortly so would like to take proper precautions with that, additionally the current pipeline utilizes Novobarcode as a demultiplexer so if anyone has suggestions on what I can implement during processing to minimize the effect of cross contamination that'd be helpful.

Are you using dual-index barcodes with different barcodes at each end? There is a known phenomena of "index switching" that occurs in Illumina reads. One way to control for this is to add in unamplified water samples that have unused barcode combinations where each of the indexes is shared with another sample.
More information can be found on this SeqAnswers thread. Here are some useful associated links:

Signal-spreading preprint -- detailed investigation of the phenomenom, including methods, graphs, and associated data
Crossblock -- tool by Brian Bushnell that attempts to remove cross contamination
Summary by James Hadfield -- mentions a few fixes suggested by Illumina, many centred around not using barcode combinations that could conflict
Illumina's quick-fire response -- they've known about it for a long time, and point out that it shouldn't be a problem for most purposes

I have suspicions that at least some of the index switching is being caused by ligation of two separate barcoded fragments during the sample prep (which definitely happens during sample prep for nanopore), but I haven't done any experimentation to work out if this is also the case for Illumina (due to money/time/goal constraints). The experimental design would be to take a prepared Illumina library just prior to sequencing, ligate on ONT adapters, then sequence on both the ONT and Illumina machines to see if chimeric reads discovered via nanopore sufficiently matches the index switching observed via short-read sequencing.

Yes, these are dual indexed samples with unique barcodes. I have not heard of doing this with water-samples but I'll look into it, and perhaps suggest it to my PI.


text:
Can you tell me where these RNAi’s cleave on the Drosophila melanogaster CG2316 mRNA?
How can I work out at which sites on the Drosophila melanogaster CG2316 mRNA the RNAi`s:

12170/FBti0089992,
12168/FBtp0030589,
107343/FBtp0042163 and
41984/FBst0034349

(stock numbers/flybase ID) cleave?
I attempted to work this out myself some time ago but with limited success. I found the Fasta sequence of CG2316 and searched for the two complementary sites of the RNAi hairpin using the word doc find function. I feel this rudimentary method was flawed on a few levels, not least because I only got a single hit. As such I have posted the problem here.

From FlyBase, FBti indicates a transposable element insertion site. Searching for these in combination with your stock numbers gives nothing. Are these commercial RNAis? If so, how are we to know what their sequences are without further information?

Firstly, I would download some sort of editor that is built to handle nucleotide sequences. There are a handful out there, but APE is a fairly popular software designed for plasmid editing. It is much more efficient at searching for nucleotide strings as it can search for reverse complements and to a certain extent, mismatches. It's also free!
http://biologylabs.utah.edu/jorgensen/wayned/ape/
For your question, I'm not exactly sure what you are asking, but I would use the BLAST alignment tool to try to see similarity between your RNAi sequences and the target mRNA. You input fasta files and the program searches for any possible overlap between the two inputs. The higher the "max-score" the more likely the two strands (or subsections of strands) have high similarity. This program also considers the complement
https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastn&PAGE_TYPE=BlastSearch&LINK_LOC=blasthome
You can use the results from BLAST to annotate in APE using different colors which is extremely helpful. You can also use BLAST in the user interface of APE.
Also keep in mind RNAi's are not always effective, even though they are predicted to bind at a given location, that may not be a biological reality
Good luck!


text:
How can I extract only insertions from a VCF file?
I'm looking to subset a standard VCF file to generate one which only includes insertions (i.e. not indels).
I can get part of the way there with:
bcftools view -v indels <vcf> | awk '{if(length($4) == 1) print}'

However this wouldn't catch an insertion that was part of a multi-allelic record with an insertion and deletion, where the reference length would be greater than 1 bp. Potentially then one way to go is a chain of decomposition, normalisation and then the same reference length filtering — surely there's a better way in one of the many VCF manipulation utils?

What does your VCF look like? Many tools produce VCFs that include the type of variant.

Safe to assume there's not a type annotation you can just grep out, though running an annotation tool first could be one answer.

Oh well, things are never that easy. Still, could you please [edit] your question and give an example of the kind of multi-allelic record you're thinking of? Ideally, give us a minimal example of your vcf file and include a couple of entries you want to keep and a couple you want to skip and show the desired output. That way we can use your example to test our solutions and we can be sure we're giving you what you need.

A one-liner:
zcat my.vcf.gz |
  perl -ane '$x=0;for $y (split(",",$F[4])){$x=1 if length($y)>length($F[3])}print if /^#/||$x'

or equivalently
zcat my.vcf.gz |
  perl -ane '$x=0;map{$x=1 if length>length($F[3])}split(",",$F[4]);print if /^#/||$x'

For simple VCF operations, I generally recommend to write a script. This may be faster than those using heavy libraries. With a script, you only parse fields you care about; most libraries unnecessarily parse every field.

On a related note, I recommend not to decompose multi-allelic sites unless necessary. Decomposing is tricky, makes VCF harder to parse and to understand and may be a potential source of errors. Here is an example:
#CHROM POS ID REF  ALT               QUAL FILTER  INFO  FORMAT  S1  S2  S3  S4
11     101 .  GCGT G,GCGA,GTGA,CCGT  199  PASS    .     GT      0/1 1/2 2/3 2/4

vt decompose+normalize produces the following VCF:
#CHROM POS ID REF  ALT QUAL FILTER INFO FORMAT  S1   S2   S3   S4
11     101 .  GCGT G   199  PASS   .    GT      0/1  1/.  ./.  ./.
11     101 .  G    C   199  PASS   .    GT      0/.  ./.  ./.  ./1
11     102 .  CGT  TGA 199  PASS   .    GT      0/.  ./.  ./1  ./.
11     104 .  T    A   199  PASS   .    GT      0/.  ./1  1/.  1/.

In theory, you can reconstruct the original VCF from this output. However, it is very challenging for a program to do that. When you compute allele frequency line-by-line, this VCF will give you wrong results. bcftools norm -m- replaces "." with "0". You can get a correct ALT allele frequency from the bcftools output, but a wrong REF allele frequency. Furthermore, vt is also imperfect in that "CGT=>TGA" is not decomposed.
My preferred output is:
#CHROM POS    ID     REF    ALT    QUAL   FILTER INFO   FORMAT S1     S2     S3     S4
11     101    .      GCGT   G,<M>  0      .      .      GT     0/1    1/2    2/2    2/2
11     101    .      G      C,<M>  0      .      .      GT     0/2    2/0    0/0    0/1
11     102    .      C      T,<M>  0      .      .      GT     0/2    2/0    0/1    0/0
11     104    .      T      A,<M>  0      .      .      GT     0/2    2/1    1/1    1/0

Here we use a symbolic allele <M> to represent "another ALT allele". You can calculate the allele frequency by looking at one line, and won't confuse other ALT alleles with REF. bgt can produce such a VCF indirectly. However, it discards all INFO, so is not a practical solution, either.
In summary, it is very difficult to decompose multi-allelic sites. When you get decomposing wrong, your downstream analyses may be inaccurate. Decomposition should be used with caution.

Just to highlight that all the steps can be done within bcftools capabilities, and since I can't just comment on @blmoore 's answer:
bcftools view --types indels <vcf> |
  bcftools norm -m - |
  bcftools filter --include 'strlen(REF)<strlen(ALT)' |
  bcftools view -H

more bult-in functions for 'bcftools filter' here

Very nice, I had only ever used bcftools view to uncompress and then parse the output into either the header, or the data records.

One method is to decompose multi-allelic records so that they're represented as one-allele, one-record using vt or similar:
bcftools view -v indels <vcf> |
  vt decompose - |
  bcftools view -H |
  awk '{if(length($5)>length($4)) print}'

Some decomposed alleles will come with excess reference padding. To left-shift and trim these, add a normalize step (NB matching these back to your input VCF becomes non-trivial):
bcftools view -v indels <vcf> |
  vt decompose - |
  vt normalize -r <reference.fasta - |
  awk '{if(length($4)==1) print}'

edit: As gringer suggests, this can also be done without vt:
bcftools view -Ou -v indels <vcf> |
  bcftools norm -Ou -Nm - |
  bcftools view -H |
  awk '{if(length($5)>length($4)) print}'

To also include complex alleles, use view -V snps (!snvs) instead of -v indels

bcftools norm will also split multi-allelic records with the -m - option

Using vcffilterjs

get the length of the REF;
loop over the ALT, ignore the symbolic
accept the variant if it's an insertion , eq:  len(ALT)>len(REF)

.
java -jar dist/vcffilterjs.jar -e 'function accept(vc){var a=vc.getAlleles();var lenRef=a.get(0).length();for(i=1;i<a.size();++i) {var alt=a.get(i);if(alt.isSymbolic()) continue;var lenAlt=alt.length(); if(lenRef<lenAlt) return true; } return false; }accept(variant);' input.vcf

.

If you're doing set operations, you could use vcf2bed:
$ vcf2bed --insertions < in.vcf > out.bed

Based on the VCF v4.2 specification, --snvs, --insertions, and --deletions are options available to filter input. In each case, the length of the reference and alternate alleles is used to determine which type of variant is being handled.


text:
10x Genomics Chromium single-cell RNA-seq data analysis options?
Provide an overview of 10x data analysis packages.
10x provides Cell Ranger which prepares a count matrix from the bcl sequencer output files and other files (see bottom of page https://support.10xgenomics.com/docs/license for the programs it uses).
What can we do with the output files?

Data preparation
Cell Ranger uses the Illumina sequencing output (.bcl) files

Make fastq files:
cellranger mkfastq ==> .fastq
Prepare count matrix: cellranger count ==> matrix.mtx, web_summary.html, cloupe.cloupe
Optional: combine multiple matrix.mtx files (libraries): cellranger aggr

Data analysis
Loupe Cell Browser visualization of cloupe.cloupe files
Count table matrix.mtx analysis options:

Python
R Cell Ranger R Kit: cellrangerRkit::load_cellranger_matrix() ==> ExpressionSet
R Scater: scater::read10XResults() ==> SCESet object
R Seurat: Seurat::Read10X() ==> Seurat object


text:
Analysis of differential transcript usage (DTU)
Recent breakthroughs in bioinformatics tools for quantification (e.g.
Cufflinks/Kallisto/Salmon etc.) and tools which can identify differential transcript usage (DTU) (e.g. DRIMSeq, Cufflinks etc.) mean that from RNA-seq data we can now relatively easy obtain a genome wide analysis of transcripts that are differentially used betwen conditions.
What can you use these results for? What systematic analysis does this transcript level data enable?

I would recommend to remove the software-recommendation tag. Here you don't ask for a tool to analyze this type of data, but what can be done analyzing this kind of data.

But I'm also asking for tools that can do it...?

I will take the liberty of giving one possible answers to my own question – but I’m very interested in other answers.
One analysis type that such data enables is the analysis of transcript switches with predicted potential consequences.
I myself have recently developed such a tool called IsoformSwitchAnalyzeR. IsoformSwitchAnalyzeR enables statistical identification (via DRIMSeq) of isoform switches with predicted functional consequences. The consequences analyzed can be chosen from a long list which includes gain/loss of protein domains, signal peptides changes in NMD sensitivity etc. The R package also enables easy visualization of isoform switches along with their consequences and it directly supports the output of Cufflinks/Cuffdiff, RSEM, Salmon and Kallisto.
Apart from enabling identification of interesting examples switch identification also enables systematic analysis of what genes are affected. For inspiration I will recommend one of my own articles (describing results obtained with IsoformSwitchAnalyzeR) as well Hector et al’s recent paper (which is not using IsoformSwitchAnalyzeR). Of particular interest and finesse is Hector’s analysis of how isoform switches can disrupt protein-protein interactions.
Looking forward to hear more ideas.

I am not sure if this has been done, is common, or are research lines, but here is what I think can be done with transcripts' differences (aside from comparing the change of the expression of each transcript):
Having the different transcripts, enables us to study also the relation between transcripts, which are co-expressed together or differently co-expressed. Maybe at the individual level the switch between transcripts of a gene is not related to a condition, but together with other transcripts it is.
Transcripts usage could be also linked to structural or genomic changes. Either short/small changes or big. At genomic level SNP, inversions or methylation..., or structural changes such as compression of a region, different distances between chromosomes/regions...

Could you elaborate on what you mean by "differently co-expressed" ?

The way a transcript (or a gene) is expressed with other transcripts can change between conditions. If in a condition transcript A is co-expressed with B, C and D, in another condition it can be co-expressed with D, E and F. See this paper

Additional insights can frequently be made in research by studying things at a lower level. By peeling back the higher layers, effects can be seen that are hidden (or reversed) when looking at aggregate data. This is generally referred to as Simpson's Paradox.
We know that RNA is functional at an isoform/transcript level, and in some cases the expression of different transcripts from the same gene can have different biologically-relevant outcomes. By looking at transcript-level information, even if the transcripts for a gene average out to be a similar level between different treatment conditions, variation at an individual transcript level can still be discovered. This is important when looking at differential expression of populations of cells, because it can clue researchers into finding out that a particular cell sub-population might be of interest; single cells typically will only produce single isoforms for each gene at a specific point in time.

I understand that this post is rather old but there has been some development in the field and this 2019 paper might be a good resource for interested individuals.


text:
Searching for gene expression data by cell line
I have two cancer cell lines (OCI-Ly18 & riva) that I want to find gene expression data for, but I'm not aware of many gene expression databases that allow searching by cell-line without searching by gene.
I tried Genevestigator on the recommendation of this thread
but found no data on either cell line.
What would be the best way to find gene expression data for these cell lines?

Try the Gene Expression Omnibus - it looks like they have some datasets.

@JSneathThompson Consider marking the answer as accepted if it seems to be the best working solution for you.

Try in this section of Cancer Cell Line Encyclopedia. It has expression data for a lot of cancer cell lines. I have tried few times and found all cell lines I wanted, except HeLa. You just have to select mRNA expression data and download the file you prefer (raw data, gene-centric, etc.)


text:
Handling sample identity in aggregated 10x libraries?
cellranger aggr can combine multiple libraries (samples), and appends each barcode with an integer (e.g. AGACCATTGAGACTTA-1). The sample identity is not recorded in the combined matrix.mtx file.
https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/aggregate#gem_groups
How can we keep and assign sample information to each cell after loading the data into R? (e.g. Seurat::Read10X())

An example solution for Seurat:
Retrieve sample IDs from the .csv used with cellranger:
samples <- read.csv(file.path("/path/to/csv", "nameof.csv"), stringsAsFactors=F)

Load the 10x dataset and initialize the Seurat object:
cells.data <- Read10X("path/to/filtered_gene_bc_matrices")

cells <- new("seurat", raw.data=cells.data)

Get barcodes and suffix:
cellcodes <- as.data.frame(cells@raw.data@Dimnames[[2]])
colnames(cellcodes) <- "barcodes"
rownames(cellcodes) <- cellcodes$barcodes

cellcodes$libcodes <- as.factor(gsub(pattern=".+-", replacement="", cellcodes$barcodes))
cellcodes$samples <- as.vector(samples$library_id[cellcodes$libcodes])

Create dataframe for meta.data argument and set up object:
sampleidentity <- cellcodes["samples"]

cells <- Setup(cells,
               meta.data=sampleidentity,
               min.cells=3,
               min.genes=200,
               do.logNormalize=T, total.expr=1e4, project="projectname")

Use the group.by argument of various functions.

The new Seurat 2.0 has the MergeSeurat() and AddSamples() functions. See this question and seurat vignette

You have to keep track of the file order you used for aggr.  The suffix number represents which represents the aggregated sample. You could store this information in a text file and load it into R independently of the Read10X function and combine it with the cell names to get their sample names. From the page you link:

Gem Groups
This is an integer that is appended to each barcode in the
  gene-barcode matrix. For example, AGACCATTGAGACTTA-1 and
  AGACCATTGAGACTTA-2 are distinct cell barcodes from different
  libraries, despite having the same nucleotide sequence.
The numbering of the GEM groups will reflect the order that the
  libraries were provided in the Aggregation CSV.


text:
Getting a "system is computationally singular" error in sleuth
I am analysing 142 samples belonging to 6 batches. Additionally, those samples belong to 72 strains, which means that for most of the strains there are two samples.
I could fit simple models (for strain and batches for instance), but when I get to the "full" model (~batch+strain), I get the following error:
so <- sleuth_fit(so, ~strain+batch, 'full')
Error in solve.default(t(X) %*% X) :
  system is computationally singular: reciprocal condition number = 5.2412e-19

I should point out that of the 72 strains, only 15 have samples in distinct batches. This means that most strains (57) have both samples in the same batch.
Is the error due to an unknown bug or rather to the experimental design? Does it mean that the information on batches cannot be used?
Thanks
EDIT I've posted the experimental design in a gist
batch   strain  replica
batch_1     strain_41   1
batch_4     strain_41   2
batch_1     strain_28   1
batch_4     strain_28   2
batch_1     strain_26   1
[...]

You have a strain completely in a batch or something along those lines. You'll have to post the entire design (the samples and their strain and batch association) somewhere.

In fact, I did mention that: "I should point out that of the 72 strains, only 15 have samples in distinct batches.". Is it correct to say that that is the root of the problem?

It depends on what you really mean by that. Post the data frame on pastebin or something.

Hi, I just edited the question with a link to the experimental design matrix

You should be able to remove any one of the following strains to end up with a rank-sufficient model matrix: 5, 10, 12, 13, 14, 15, 19, 26, 28, 3, 30, 32, 36, 39, 41, 45, 46, 49, 5, 50, 52, 53, 58, 59, 60, 69, 8. As an aside you can figure this sort of thing out as follows (I read your dataframe in a the d object):
> m = model.matrix(~batch+strain, d)
> dim(m) # 142 row, 77 columns, so minimum rank is 77
> qr(m)$rank # 76, so just barely rank insufficient
> #see if we can remove a single column and still get rank 76
> colnames(m)[which(sapply(1:77, function(x) qr(m[,-x])$rank) == 76)]

You obviously don't want to remove the batch columns or the intercept. The normal tricks that you can sometimes use to get around this issue with case-control studies don't appear to help here, which is why I would just drop a strain and call it done. Keep in mind that your power is still likely terrible. I generally recommend at least 6 replicates per group (scale down the number of groups to fit your budget).
EDIT:
Once the desired strain is removed from the model, it can be fit directly into the sleuth_fit function to obtain the full model:
> m = m[, -9] # whatever column to drop to get the appropriate rank
> so <- sleuth_fit(so, m, 'full')

That's just briliant; I removed a strain and I can now run the full model to take into account batch effects. I will update your answer to include a snippet on how to feed the "chopped" model into sleuth

This happen when the variables (strain +batch) create a design matrix like this:
batch strain
1 1 #
1 1 #
1 2
2 2
3 3
4 3
...
16 72

Which means that some of the covariates are not linearly independent (ie batch 1 and strain 1), all the strain 1 is in batch 1.
You can correct for batch effects, but not in this design of the linear model (if you want to take into account the strain). You could do one batch more with those 15 strains that are in a single batch (if they are in different batch between them) that way you would get an independent design.
Many strains in a single batch are in the same batch. You need to increase the number of samples (recommended anyway due to the low number of samples per strain) to avoid this problem.

There are a lot of related question in Bioconductor support forum, from where I expanded an answer.

Yes, this is how my design matrix looks like. I guess that means that I will have no way of correcting for batch effects?

You can correct for batch effects, but not in this design of the linear model (if you want to take into account the strain).
I am unsure (at the moment) if you could use combat or other batch effect removal tool to incorporate the batch effect, other option could be to try to use surrogate variables (which will detect the batch effect AFAIK).Or you could do one batch more with those 15 strains that are in a single batch (if they are in different batch between them) that way you would get an independent design

@mgalardini As I said above, please post the entire design dataframe somewhere. Then we can actually look an tell you what your options are rather than speculating.


text:
How can we find the distance between all residues in a PDB file?
If we have a PDB structrure, how can we find residues physically interacting with each other in space? I know that we must find the distance between residues and if the distance is less than 5-6 Angstrom, we say that residues are physically interacting. But how can we find the distance between all residues and how can we finally determine the distances between all residues? Is there a software or webserver for that?

These interaction may be already annotated in the PDB, see Connectivity Section and Connectivity Annotation Section in the PDB spec.

One thing you're going to have to define for yourself is what you mean by "distance between residues". Do you mean closest atom pair? Closest heavyatom pair? Distance between C-alphas? The exact mechanism is probably going to be similar, regardless, but the details and the results will differ based on the definition.

I mean the distance between C-alphas

Is the Ring 2.0 webserver (http://protein.bio.unipd.it/ring/) good enough for this purpose?

If anyone is interested in speed benchmarks of distance calculations using various software, there are some at https://github.com/jgreener64/pdb-benchmarks. There are also example scripts for each software in that repository.

If you need to process multiple files, you could use Biopython to parse a PDB structure.
from Bio.PDB import PDBParser

# create parser
parser = PDBParser()

# read structure from file
structure = parser.get_structure('PHA-L', '1fat.pdb')

model = structure[0]
chain = model['A']

# this example uses only the first residue of a single chain.
# it is easy to extend this to multiple chains and residues.
for residue1 in chain:
    for residue2 in chain:
        if residue1 != residue2:
            # compute distance between CA atoms
            try:
                distance = residue1['CA'] - residue2['CA']
            except KeyError:
                ## no CA atom, e.g. for H_NAG
                continue
            if distance < 6:
                print(residue1, residue2, distance)
        # stop after first residue
        break

If you need to look at one structure, using a viewer perhaps would be easier. You could try PyMOL: (how to measure distance). There are other PDB viewers, some of which can work even through a browser.

Dear lakov, I installed Anaconda and Biopython packages and run the above code in PyCharm 5.0 and got the following error:

Traceback (most recent call last):
File "C:/Users/VN7-592/PycharmProjects/untitled/sara.py", line 7, in
structure = parser.get_structure('PHA-L', '1FAT.pdb')

File "C:\Users\VN7-592\Desktop\Anaconda\lib\site-packages\Bio\PDB\PDBParser.py", line 81, in get_structure
with as_handle(file, mode='rU') as handle:
File "C:\Users\VN7-592\Desktop\Anaconda\lib\contextlib.py", line 17, in enter
return self.gen.next()
File "C:\Users\VN7-592\Desktop\Anaconda\lib\site-packages\Bio\File.py", line 88, in as_handle
with open(handleish, mode, **kwargs) as fp:
IOError: [Errno 2] No such file or directory: '1FAT.pdb'

You don't have the 1FAT.pdb file in the working directory, or your operating system is case sensitive.

But this gives only intermolecular distances? BioPython ignores crystal symmetry? In general, to check for contacts between residues X and Y one would need to calculate distance between X and the nearest symmetrical image of Y.

Could you use CCP4's NCONT program? There's a GUI and a command line interface, whatever suits. You can specify which chains you want to target and interact with and set a cut off for distance. The bonus here is once you're in you have a nice suite of other structural tools to use.
If you're just doing it once, the GUI is friendly enough to work things out, if you're doing a batch then you can run it across several files via the command line.
Download

I have downloaded the program but I haven't been able to work with it yet. I'm trying to run it

@Sara what is it you're having a problem with?

It doesn't run with an error. The program says choose a project name and even when I have chosen a project name, again says choose a project name and stops running

As part of a project me and some teammates did a script that outputs visual maps of distances between residues. It uses Biopython.
The module contact_map.py does what you are looking for. As an example, if you want to find the residues whose CA are below 5 you can run the following command:
python3 contact_map.py pdb1cd8.ent -a CA -CA 5

This will produce three files:
distance_map_pdb1cd8_CA.png # Heatmap of the distance between the residues
contact_map_pdb1cd8_CA.png # Black/White heatmap: If it is at that min distance
contact_map.log  # The actions taken

If you don't have downloaded already the pdb structue you can use the main module:
python3 cozmic.py real 1cd8 -a CA -CA 5

How it calculates distances? Inside the asymmetric unit only, or as the shortest distance between the two residues in the crystal (or in biological assembly)?

@marcin It calculates the distance between the C alpha (or beta if you chose so) between each residue in the file.

Usually, the file == asymmetric unit. If you'd like to find all contacts in biological assembly or in crystal you'd need to take the symmetry into account.

Thanks for the tips, maybe you want to contribute to the modules.

I'm working now on a project that in the future may be competitive to Bio.PDB, but probably not this year yet :-)

I think Bio.PDB is being actively modified, but good luck with your work

You can use MDtraj. The package is easy to install using Anaconda.
You can get the interacting residues using the following snippet (taken from http://mdtraj.org/1.6.2/examples/native-contact.html)
heavy_pairs = np.array(
    [(i,j) for (i,j) in combinations(heavy, 2)
        if abs(native.topology.atom(i).residue.index - \
               native.topology.atom(j).residue.index) > 3])

# compute the distances between these pairs in the native state
heavy_pairs_distances = md.compute_distances(native[0], heavy_pairs)[0]
# and get the pairs s.t. the distance is less than NATIVE_CUTOFF
native_contacts = heavy_pairs[heavy_pairs_distances < NATIVE_CUTOFF]
print("Number of native contacts", len(native_contacts))

Distance between C-alphas can also be found using the BioStructures.jl package in Julia:
using BioStructures
struc = read("1AKE.pdb", PDB)
calphas = collectatoms(struc['A'], calphaselector)
dm = DistanceMap(calphas)

Then dm[5, 10] gives 12.19677604943208 etc.
Other variations also work, see the docs.


text:
Dataset: Locations of regulatory sequences in the human genome?
I am looking for the positions of annotated regulatory sequences (promoters, enhancers and suppressors) in the human genome.
I looked at Ensembl regulatory Build and PAZAR but I am not used to look for datasets and I failed to find what I was looking for.
Can you help to find annotated regulatory sequences (start and end positions) in the human genome?

Nice to see you here! Could you clarify what type of regulatory sequences you'd want? All of them? Would you include promoters in that or only things like splicing enhancers/supressors?

Yes, both promoters, enhancers and suppressors (TBH, I did not think of promoter when asking my question but I'd like to consider them as well). Annotation on what type of sequence (promoter, enhancer/supressor) it is would be great. Post edited to add this info. Thanks @terdon

This can be done quite easily using Ensebl's BioMart.

Choose the Ensembl Regulation database:

Select the "Human Regulatory Features" dataset:

That's basically it right there, just click on "Results":

Export to file and click "Go":

This will download a file called mart_export.txt which looks like this (I chose TSV for tab separated values):
$ head mart_export.txt
Chromosome/scaffold name    Start (bp)  End (bp)    Feature type
18  76429380    76430144    Open chromatin
8   66405962    66406502    Open chromatin
4   61184401    61184600    CTCF Binding Site
X   40733600    40737000    Promoter
5   97407001    97407200    CTCF Binding Site
X   73946201    73946600    Promoter Flanking Region
15  19948201    19949200    CTCF Binding Site
5   11302601    11303143    Open chromatin
2   208407801   208408000   CTCF Binding Site

As you can see the format is pretty straightforward. The first column has the chromosome or scaffold name (it should always be a chromosome for human), the second and third are the star and end positions and the last field is the type of region.
You can also choose to limit by region type by choosing a filter (click the "Filters" link) in BioMart before downloading or by simply parsing the file once you've downloaded it.

Amazing! The Images link are broken but your step by steps explanations were very easy to follow even without images. Thanks a lot @terdon! +1

@Remi.b they are? They work fine for me, that sounds like a bug. Anyway, glad I could help :)

Promoters and enhancers based on CAGE-seq from FANTOM5:
http://fantom.gsc.riken.jp/5/data/
Just choose the category and download bed-file, e.g. CAGE-peaks would be promoters. There are READMEs in each folder that explain the data.


text:
How to estimate whether a long-read is meaningful sequence?
The setup
Imagine that I work on an organism without a reference genome, and that the closest reference genome I can get is quite diverged. E.g. ~10% diverged in terms of SNVs when measured with short reads, and also has a lot of structural variants too.
Now imagine I get a 1 million base-pair long-read (e.g. from Nanopore data) for my organism. The question is this:
How can I estimate the proportion of the read that is meaningful sequence vs garbage?
Some things that probably won't work
Most standard approaches won't work here. E.g. I could try mapping the read to the reference, but even if the read was perfectly good I wouldn't expect most of it to map thanks to true structural variations between the read and the reference. The same goes for standard alignment or BLAST.
Some things that might work
The best naive method here seems to be to cut the read up into smaller pieces (either overlapping or not) and use standard approaches to map/align each of these.
So, what have people tried for this? And what tools have you used and why?

As a first pass, you could check if the read is chimeric. Porechop searches for known nanopore adaptors both on the ends and through the middle of the read. This won't resolve issues around blocked or empty pores, but it will at least check if you have found two long-ish reads lumped into the same file.
By default, Porechop splits chimeric reads into two (or more, I suppose!) non-chimeric reads, but the option --discard-middle would be a quick and easy way to check - run it on a fasta file containing only your long read, and if the output is empty, the read is chimeric.

If you're looking at a single organism, in the absence of a reference genome you can map other reads to suspect reads and look at coverage. Looking at the actual sequence can also be useful: real DNA usually doesn't have an abundance of two different bases.
Nanopore reads also give another way to see if the read looks weird by having a look at the raw signal. Our chimeric reads paper gives a few examples of what DNA sequence should look like under normal circumstances. Here's the first raw signal figure from that paper:

If there's lots of contiguous signal with very similar current levels (i.e. it looks like the "stall" region, then it's not a good read and any basecalls should be ignored.


text:
Finding homologs of a protein sequence
I have a refseq ID of a protein from E.coli and I want to find homologs of this protein. I ran Blast against refseq database but I got a lot of sequences most of which were from Ecoli again. I decided to run PSI-Blast to get more divergent species, but I do not exactly know if my result are real homologs or false positives. what is your idea for finding homologous protein sequences from more divergent species? And what can I do for selecting the real (but not false positive) hits?

False positive can be judged by the E and P-values so remember to use them.

I would suggest that if you can robustly obtain the same topology using a variety of phylogenetic methods and models on your proteins, there is a chance that the homology is true. Unfortunately, this might be difficult to obtain when dealing with quite distant species, even if the homology is true.

It depends on what you're looking for. If you're just looking for sequence homology, then you can simply pick the best hits from a blast search. If, however, you are referring to functional homology, if you are looking for the protein which has the same functions as your query, then it's more complicated.
Sequence homology is not enough to infer functional homology. For example, you can have cases of gene duplication and subsequent functional divergence. Such paralogs are still homologs (paralogs are a subset of homologs), but they don't necessarily have the same function. It is also often the case that the homolog (be it orthologous or paralogous) of a protein in species B has a completely different function than its homolog in species A despite a high level of sequence similarity. This is usually very hard to determine in silico.
To find the functionally homologous protein(s), you would ideally need to identify the essential residues that allow your protein to perform its function. This could be done using something like PFam which will identify protein domains. You can then check whether the homologs you find also have this domain.
This is essentially what PSI-blast does. Although it doesn't take domains into account, each successive iteration is used to build a model of your proptein. The model is built under the reasonable assumption that highly conserved residues are important. So it will consider more diverged sequences as homologous if those residues are conserved.
If you know how your protein works and what residues are important, you can use that knowledge to refine the results of your PSI-blast. If you don't, you'll have to use only "good" hits to make the model. One way to do this, for well studied proteins, is to only add proteins that are already annotated as homologs of what you are looking for to build your model, then use that model to search in un-annotated species.
If you don't know, you could try looking for recognizable protein domains in your query protein (use PFam) and then use the HMM (hidden markov model) of the domain to identify important residues. For example, this is the HMM logo for the zf-A20 zinc finger domain:

The huge cysteine (C) residues are shown that size because they are very conserved across proteins carrying this domain and, presumably, are functionally important for the domain. So, if you pass your protein through PFam and identify domains, find the important residues and make sure all your homologs have those conserved. If using PSI-blast, only include sequences where those residues are conserved in the results you keep.
Finally, another useful tool that works in the same way is HMMER. This takes a protein alignment as input, like PSI-blast builds an HMM model from it and then can use this model to query a protein database for more hits. Methods like HMMER and PSI-blast are far better than simple sequence similarity approaches when looking for homologs.

Great tool HMMER, is as fast as BLAST and works pretty well

Dear terdon, your answer was the best answer someone could have given to this question and it was exactly the answer I was looking for. Thank you very much for your excellent answer

@Sara you're very welcome :)

It sounds like Smart BLAST might do what you want. Here's the description of it's goal:

SmartBLAST is a new and experimental NCBI tool that makes it easier to
  complete common sequence analysis tasks, such as finding a candidate
  protein name for a sequence, locating regions of high sequence
  conservation, or identifying regions covered by database sequences but
  missing from the query.
To do this, SmartBLAST performs the following tasks in much less time
  than it takes to run a typical BLASTp search:

a BLASTp comparison of the query with the closest matching sequences available;
a parallel BLASTp search to find the closest matches to high quality sequences from model organisms;
a multiple alignment between the query and five of the closest matching sequences (usually including two high quality sequences);
an analysis that produces a phylogenetic tree from the multiple sequence alignment.

[from NCBI Insights]

SmartBLAST is a very good idea as it gives so much more divergent organism, but is there any way to give its result to PSI-BLAST for more iterations?

@Sara just note the accession numbers of the sequences found by Smart BLAST and make sure to add them to your alignment in the psi-blast. They should be found as hits there too.


text:
What command/invocation is used to generate NCBI 16SMicrobial blastdb
I'm looking for the exact invocation used to generate the 16SMicrobial database that you can download from here:
https://ftp.ncbi.nlm.nih.gov/blast/db/
I'm hoping to create the same type of blastdb with the same type of metadata with custom sequences.
Platform isn't an issue, but let's say on ubuntu 14.04 or 16.04.
I would like to replicate the creation of the database as closely as possible.
 The most important feature is the taxonomic information as can be seen here:

The databases on the FTP site contain taxonomic information for each
  sequence, include the identifier indices for lookups, and can be up to
  four times smaller than the FASTA. The original FASTA can be generated
  from the BLAST database using blastdbcmd

Creation of a blastdb using makeblastdb from a set of a fasta sequences is not an issue and can be achieved via:
makeblastdb -in <your_file.fasta> -dbtype nucl -out <database_name>

My question is specifically about the invocation NCBI uses to add the metadata that is present in the NCBI's 16SMicrobial blast database as I am keen to make sure I have replicated the process as closely as possible.

What properties, specifically?

Most crucially the taxonomic information, but I would like to replicate this creation of this database as closely as possible.

Running blastdbcmd -db foo -info provides a little information but I haven't seen anything which will report exactly how a blastdb was created.
A blastdb readme suggests that only the parse_seqids option has been added to the standard parameters.

For those from NCBI, the following makeblastdb commands are recommended:

For nucleotide fasta file:   makeblastdb -in input_db -dbtype nucl -parse_seqids
For protein fasta file:      makeblastdb -in input_db -dbtype prot -parse_seqids

The -parse_seqids requires NCBI fasta headers in the correct format (spec can be found here).

I don't know what specific properties you are considering. As far as I know, that's just a normal blast database, like any other. These are produced using the formatdb command which is part of the ncbi-blast software. You haven't specified what system you are using so I can't help you find and install it, but the command will be something like:
formatdb -p F -i yourSeqs.fasta -n yourDBName

formatdb : the command
-p F : protein: false. This is needed to make databases of nucleotide sequences.
-i yourSeqs.fasta : your input file with all your sequences in fasta format.
-n yourDBName : the name (choose whatever you like) for your database.

You can then use yourDBName as a database to blast against. The above assumes yuou are

I would recommend using the newest ncbi-blast+ suite instead of the (now deprecated?) legacy-blast. In that cases formatdb would become makeblastdb, with slightly different arguments.

@mgalardini that sounds very reasonable. I'm afraid I'm showing my age. This is the only one I've ever used but I last did so more than 10 years ago now. Please post an answer with this new tool you mention. I'm sure it's the Right Way© to go today.

@terdon - whatever works I say ;).  Thanks for the help, I've updated the question to make it (hopefully) more clear.  I am specifically looking to add metadata to a blastdb and looking for a reference or similar for what metadata is added to the 16SMicrobial database.

I am not sure whether you are looking for how to produce taxonomy file or rep_set file (qiime reference file types). But you can produce fasta file (rep_Set) like that:
$blastdbcmd -db 16SMicrobial/16SMicrobial -outfmt '%f' -entry 'all' > 16SMicrobial.fasta

Then you can create QIIME-compatible id_to_tax maps with this code:
https://github.com/mtruglio/QIIME_utilities/blob/master/id_to_tax_mapmaker.sh
Hope this helps.


text:
Random access on a FASTQ file
I would like to select a random record from a large set of n unaligned sequencing reads in log(n) time complexity (big O notation) or less. A record is defined as the equivalent of four lines in FASTQ format. The records do not fit in RAM and would need to be stored on disk. Ideally, I would like to store the reads in a compressed format.
I would prefer a solution that does not require any extra files such as for example a reference genome.
The title of this question mentions a FASTQ only because FASTQ is a common format for storing unaligned reads on disk. I am happy with answers that require a single limited transformation of the data to another file format in time complexity order n.
Update
A clarification: I want the random record to be selected with probability 1/n.

Related question: https://bioinformatics.stackexchange.com/q/453/292

As wkretzsch suggested this was worthy of an actual answer, I feel the obvious solution is missing here; index the FASTQ.
Index it
As much as I typically hesitate to jump to a solution that requires a script or framework (as opposed to just unix command line tools), there is sadly no samtools fqidx (perhaps there should be), and existing answers suggest a lot of munging. Whilst they probably work, some appear cumbersome and have many steps in which you could make a mistake.
So, to keep things simple - a quick and dirty alternative approach might be to use biopython, seeing as it already has this functionality implemented to do this, and if installed, is trivial to use:
from Bio import SeqIO
fq = SeqIO.index("myfastq.fq", "fastq")

Once you've acquired an index, you'll get fast random access for any read:
# Random access by read name (I like owls)
record = fq["HOOT"]
record
#> SeqRecord(seq=Seq('ACGTACGT', SingleLetterAlphabet()), id='HOOT', name='HOOT', description='HOOT', dbxrefs=[])

# We can get the sequence
record.seq
#> Seq('ACGTACGT', SingleLetterAlphabet())

# and qualities
record.letter_annotations
#> {'phred_quality': [41, 41, 41, 41, 41, 41, 41, 41]}

If you want to select arbitrary random records, you could use something like randrange to select between 0 and the length of the references list.
from random import randint

# Coerce keys to a list, as `dictionary-keyiterator` has no
#   __getitem__ attribute that would allow integer indices
# Note also that this doesn't necessarily guarantee a sorted order
#  but I guess that doesn't matter if you just want random records
key_list = list(fq.keys())

# Select a random key
# Note we use len(key_list)-1 as randint endpoints are inclusive
random_readname = key_list[ randint(0, len(key_list)-1) ]

# Get your record
rand_record = fq.get( random_readname )

If you want multiple records, you'll probably want sample (to avoid replacement) instead:
from random import sample
N = 100
random_indices = sample(xrange(len(key_list)), N)

for key_i in random_indices:
    random_readname = key_list[ key_i ]
    rand_record = fq.get( random_readname )
    # ...

For what it's worth, I think biopython holds this index in RAM, so if your file is absolutely massive, you might need to be more clever. If that's the case, you could iterate through the FASTQ once, and output the readname, file offset and length - akin to a FASTA fai.
Update
See samtools fqidx

Bio.SeqIO.index does hold the index in RAM as a Python dict, but Bio.SeqIO.index_db puts the index in a reusable SQLite3 database file.

@peterjc Aha! That's good to know.

Arbitrary record access in constant time
To get a random record in constant time, it is sufficient to get an arbitrary record in constant time.
I have two solutions here: One with tabix and one with grabix. I think the grabix solution is more elegant, but I am keeping the tabix solution below because tabix is a more mature tool than grabix.
Thanks to user172818 for suggesting grabix.
Update
This answer previously stated that tabix and grabix perform lookups in log(n) time. After taking a closer look at the grabix source code and the tabix paper I am now convinced that lookups are independent of n in complexity. However, both tools use an index that scales in size proportionally to n. So, the loading of the index is order n. However, if we consider the loading of the index as "...a single limited transformation of the data to another file format...", then I think this answer is still a valid one. If more than one record is to be retrieved, then the index needs to be stored in memory, perhaps with a framework such as pysam or htslib.
Using grabix

Compress with bgzip.
Index the file and perform lookups with grabix

In bash:
gzip -dc input.fastq.gz | bgzip -c > output.fastq.gz

grabix index output.fastq.gz

# retrieve 5-th record (1-based) in log(n) time
# requires some math to convert indices (4*4 + 1, 4*4 + 4) = (17, 20)
grabix grab output.fastq.gz 17 20

# Count the number of records for part two of this question
export N_LINES=$(gzip -dc output.fastq.gz | wc -l)

Using tabix
The tabix code is more complicated and relies on the iffy assumption that \t is an acceptable character for replacement of \n in a FASTQ record.  If you are happy with a file format that is close to but not exactly FASTQ, then you could do the following:

Paste each record into a single line.
Add a dummy chromosome and line number as the first and second column.
Compress with bgzip.
Index the file and perform lookups with tabix

Note that we need to remove leading spaces introduced by nl and we need to introduce a dummy chromosome column to keep tabix happy:
gzip -dc input.fastq.gz | paste - - - - | nl | sed 's/^ *//' | sed 's/^/dummy\t/' | bgzip -c > output.fastq.gz
tabix -s 1 -b 2 -e 2 output.fastq.gz

# now retrieve the 5th record (1-based) in log(n) time
tabix output.fastq.gz dummy:5-5

# This command will retrieve the 5th record and convert it record back into FASTQ format
tabix output.fastq.gz dummy:5-5 | perl -pe 's/^dummy\t\d+\t//' | tr '\t' '\n'

# Count the number of records for part two of this question
export N_RECORDS=$(gzip -dc output.fastq.gz | wc -l)

Random record in constant time
Now that we have a way of retrieving an arbitrary record in log(n) time, retrieving a random record is simply a matter of getting a good random number generator and sampling. Here is some example code to do this in python:
Using grabix
# random_read.py
import os
import random

n_records = int(os.environ["N_LINES"]) // 4
rand_record_start = random.randrange(0, n_records) * 4 + 1
rand_record_end = rand_record_start + 3
os.system("grabix grab output.fastq.gz {0} {1}".format(rand_record_start, rand_record_end))

Using tabix
# random_read.py
import os
import random

n_records = int(os.environ["N_RECORDS"])
rand_record_index = random.randrange(0, n_records) + 1
# super ugly, but works...
os.system(
    "tabix output.fastq.gz dummy:{0}-{0} | perl -pe 's/^dummy\t\d+\t//' | tr '\t' '\n'".format(
        rand_record_index)
)

And this works for me:
python3.5 random_read.py

Disclaimer
Please note that os.system calls a system shell and is vulnerable to shell injection vulnerabilities. If you are writing production code, then you probably want to take extra precautions.
Thanks to Chris_Rands for raising this issue.

You could shuffle the FASTQ once and then read sequences off the top of the file as you need them:
gzip -dc input.fastq.gz | paste - - - - | shuf | tr '\t' '\n'| gzip -c > output.fastq.gz

I would recommend pigz as a replacement for gzip in the compression step if you have it available.
The downside of this approach is that you only get n reads before you need to run the shuffle again, and apparently shuf holds all data in RAM, so it would die with an out of memory error if the FASTQ file does not fit into RAM as is specified in the question.
Using sort -R is complexity n log(n) and uses temporary files, so it should not run out of memory:
gzip -dc input.fastq.gz | paste - - - - | nl | sort -R | perl -pe 's/\s*\d+\t//' | tr '\t' '\n'| gzip -c > output.fastq.gz

The nl and perl commands are necessary to make sure that identical records are not sorted next to each other.

One possibility is to:

reformat the data such that each record is a single line containing the read description, bases, and quality scores
pad out each record to a maximum length in each field such that every record in the file is the same number of bytes
the total number of records can now be calculated as file size / record size
choose a random record number between 0 and the total number of records
binary search over the reformatted file until you obtain your read

This would get you the log(n) lookup time you want.
Of course, the data wouldn't be compressed.  You could 2-bit encode the bases and quantize the quals to save some space, but that'd be lossy and is perhaps not what you're looking for.  Alternatively, you could block-gzip the reformatted data and keep a record of how many blocks are in the file and how many reads are in each block (since the filesize will no longer reflect the number of records in the file).  Then to obtain a specific read, you'd calculate the block number it'll appear in, decompress the block, and return the appropriate read.

I wrote a tool called sample that you can use to do random sampling without reading the entire file into memory.
It can be used where GNU shuf fails for lack of sufficient memory.
It requires two passes through the file to do a random sample, but the second pass is generally fast(er) as it uses mmap routines to do cached reads.
If you do repeated samples, the repeated samples are also mmap-ed (cached) and will run quickly.
You might use it on a FASTQ file like so:
$ sample -k 1234 -l 4 in.fq > out.fq

It parses the input file into records by every four newline characters (such as the format of a FASTQ file), reading line offset positions into memory. So the memory overhead is relatively very low.
It then applies reservoir sampling on those line offsets to write out a random sample (say, 1234 records in this example) to standard output.

Your tool sample uses reservoir sampling, which is order n if I am not mistaken.  However, the question asks for an order log(n) algorithm.  This may be an XY problem, since the tool does only need one pass to retrieve k records (1234 in your case). However, In my particular case I do not know k ahead of time, which is the reason for the log(n) complexity requirement.

You can omit -k to shuffle the entire file, or specify it as a shell variable if derived at runtime (i.e. -k ${K}).

One of the most thorough treatments of this question (or a similar question: grabbing a random subset of reads) was given by Jared Simpson in a blog post a few years ago. http://simpsonlab.github.io/2015/05/19/io-performance/
If you just want to grab a single random read, Jared's benchmarks suggest that seeking to a random position in the file and then retrieving the next complete read should be the most performant option.
However, more generally, if you want to take a random subset of reads, there are many factors that affect performance.

Caching can dramatically improve query times. It would be interesting if they test out a mmap-based approach, over fseek. Especially with SSD hardware.

To emphasise the issue (as outlined in the 1/n update), consider an input file with one record that is 5 million bases long, and one records that is 100 bases long. You want an equal probability of selecting either of these two records. Any random seek methods will overwhelmingly pick out the long record.
I expect that indexing the locations of record starts is really the only workable option here, particularly if multi-line records are possible. Create an index file containing the start locations of each record (as identically sized integers, e.g. 64-bit), then sample from the index file (which is an identical length for each record) to fetch the start location. I'd envisage that this file would only contain the start locations; any additional metadata (including sequence name) would require seeking in the original file.
Once the indexing is done, the file can be compressed with bgzip, with specific offsets retrieved using the -b and -s options. However, I expect that compression would not be particularly efficient if multiple random records were desired.

In reference to your third paragraph: I am probably missing something, but would indexing not be order n?  And I believe that's acceptable as a pre-processing step in the question?

n is defined as the number of "unaligned sequencing reads" right in the first sentence.  I'm not sure how to be more clear. But please do suggest an edit now that you know what is meant.

I am conflating "record" with "read" in my question.  I shall make an edit to be more specific. I was adopting the wording of biopython, which talks of a sequencing record as the combination of a sequence, a list of qualities and a description.

"Anything that needs to read through all the lines will need at least n, so anything that does more than that has complexity greater than n" - That is not true of complexity in big O notation, which is what I am using here.  I shall make that clearer as well.

I really like the random seek method with block compression.  The problem I have with it is that it does not guarantee uniform sampling of records if the records are not all of identical size.

Here's another approach that doesn't require any indexing, using BEDOPS bedextract to do a log(n) sample on a sorted BED file. Your sample contains random records with equal probability 1/n.
This approach requires a single O(n) pass through the file to transform it to a BED file:
$ cat records.fastq | paste - - - - | awk '{ print "chrZ\t"s"\t"(s+1)"$0 }' > records.bed

Store the intervals in a separate file:
$ cut -f1-3 records.bed > intervals.bed

To do a random sample of k elements, shuffle the intervals file and preserve the order of shuffled elements.
You can do this with the sample tool I outlined earlier:
$ sample -k ${K} -s intervals.bed > intervals-sample.bed

Or you can shuf and sort-bed to do the same thing:
$ shuf -n ${K} intervals.bed | sort-bed - > intervals-sample.bed

There's an O(klog(k)) cost here, but if k <<< n, i.e., you're working with whole-genome scale input, this cost is amortized over the log(n) search performance.
Next, use bedextract to do a binary search on the records, and delinearize to get back to FASTQ:
$ bedextract records.bed intervals-sample.bed | cut -f4 | tr '\t' '\n' > sample.fq

With Unix I/O streams, this can be done in one pass:
$ sample -k ${K} -s intervals.bed | bedextract records.bed - | cut -f4 | tr '\t' '\n' > sample.fq

By baking the sort order into records.bed, you're guaranteed the ability to do a binary search, which is log(n).
Note: Further, by linearizing the FASTQ input to a BED file and querying on BED intervals, you have equal probability of picking any one interval (interval == FQ record). You can draw an unbiased sample without the hassle of creating and storing a separate index.

I wrote a tool called strandex that matches entire FASTQ records starting at a random offset inside a decompressed file (since non-block compressed gzip files are not seek-able). It was a response to this comment about the screed FASTQ parser in khmer, and was initially just a proof of concept but it's been useful for me and a few other people. The general idea is to:

seek to a random offset in the file
read a chunk of the file
test if the regex pattern @.+[\n\r]+.+[\n\r]+\+.*?[\n\r].+[\n\r] matches
if so, extract the matching FASTQ record using regex capture groups
if no match, goto step 2

The CLI script and python library can be installed with pip install strandex, and by default the script samples -n reads using file offsets starting at a random number (reproducable by setting -s seed), then moves a number of bytes to uniformly sample the entire file, determined by the file size. In this way the process is not truly random, but random enough, and avoids sampling any part of the file too much.
If you specify -n reads less than the total number in the file you get downsampling, and if -n is greater than the total number of reads you get upsampling.
The advantage of this method over others is - it's a one-pass method, and since it's using the python C internals regex engine it's quite fast.

Cool! Thanks for this answer. Does this tool have a way of guaranteeing uniform sampling of records with uneven sizes?

I suppose you could perform rejection sampling of the record with probability 1/L that you keep the record, where L is the number of chars in a record...

@wkretzsch unfortunately it does not. The tool even fails if you have paired FASTQ files that have been trimmed to different lengths.

I believe some of these answers misinterpret the original question.  The OP (I believe) doesn't wish to find particular sequences, so indexing is unnecessary.  I believe the OP just wants a random sampling of the dataset.

If the number of reads you need in small, sure you can seek on the file and return the next complete record.  Remember the read IDs to sample without replacement.  Jumping around on disk is slow, however, and you'll be reading 4k or whatever your disk buffer size is every time.
If the number of reads you need is large, read the entire file and decide for each read whether or not to sample it.
The way I do it (for Illumina): estimate the number of reads using the file size and read length and header length, decide how many reads I want, skip the first bunch of reads (due to edge effect, e.g. 20-100k; you can use seek or just read from start and discard), and read the next batch sequentially.  I suppose it depends on what you're using the reads for; i.e. if you want random reads from your sample for biological analysis, or you want random reads from your run for QC analysis.  This method is more appropriate for the former, but works well enough for the latter too.

I suppose I can provide code for this if needed, but keep in mind that it'd probably have to be done in C.
One could make this easier by using a bgzf compressed fastq file. Yes, BAM uses that already, but since Illumina's software now defaults to producing that there should be a reasonable amount of it already available.

Count the number of blocks in the file.
Randomly select one of these blocks.
Perform reservoir selection on the entries in that block.

That won't be strictly O(log N), but it's lower on memory and IO requirements than parsing through the whole file and doing shuffling or otherwise completely munging the data into an otherwise not-terribly-useful format. This also has the benefit of allowing selecting more than one random entry a bit faster, since you can memoize things.

A quick and dirty solution could be to convert the FASTQ file to two FASTAs, one of the them storing the bases and the other one storing the qualities, and then use standard methods for random access for FASTA.
The only complication is that base qualities can contain >. However, this is a problem only when it is the first character of a line and we can fix it by prepending, e.g., I to each quality sequence.
If also other people think that this could be a sufficient solution, I will expand the answer and add all the commands.

I think this would be an acceptable solution if it is possible to reconstruct the entire randomly selected FASTQ record (including description) from this scheme?

Yes. Each record will be identical in length (assuming the 'I' hack included a corresponding '-' hack for sequence), so if it is known where in file A the start of a sequence record is, then the same location will be the start of a sequence record in file B.

@gringer, What about the read descriptions?

@wkretzsch Do you have any specific read in mind? If yes, could you create a gist for it? Basically all FASTQ I have ever seen had the + and @ lines identical.

The headers for the quality and sequence would be identical in such a paired situation.

This read is from the Genome in a Bottle Ashkenazi Jew trio. It has an @ line with several characters, but only a single character on the + line.  Furthermore, the assumption that all @ lines are of identical length seems evil to me.

It's not that the header lines would be the same length; it's that the header line for each read would be the same as its quality counterpart. Regardless, this also doesn't deal with the updated request for an equal probability for each sequence.


text:
How to extract fasta from a blastdb
How to extract the sequence used to create a blast database.  This is useful when you download a blastdb from somewhere else e.g. one of the databases provided by NCBI including the 16SMicrobial database. Or alternatively, when you want to double check which version of a sequence you have included in a blastdb.

You can extract fasta sequence from a blastdb constructed from a fasta file using blastdbcmd which should be installed when you install blast/makeblastdb.
blastdbcmd -entry all -db <database label> -out <outfile>

If you had a database called my_database which contained the files:

my_database.nhr
my_database.nsq
my_database.nin

and you wanted your fasta output file to be called reference.fasta you would run the following:
blastdbcmd -entry all -db my_database -out reference.fasta


text:
How does blastn sort output?
I used blastn to search one a genome database for the sequences in a file.
blastn -out output.txt -outfmt 6 -query sequence_list -db genome_database -perc_identity 100

These sequences are all from one contig of the genome so the output was all hits in one contig. I then appended more sequences to the sequence file and used blastn again. However the weird thing is the output did not start with the same hits as the last blast. I expected if blastn iterated through sequences in a query file that it would output the original sequences first then go through the appended sequences. Since this was not the case, how does blast order its output?
I also checked if it was sorted by which contig the hit was in and this is not the case.

You mean which query's results are shown first? I've always seen them appear in the order they were found in the original file. Is that the actual file you used? That would be considered as one long sequence, not many. Is that what happened?

Edit I stand corrected. @Daniel Standage correctly pointed out that the below is only true for a single query, or for the sorting of hits within those found by a given query. Otherwise they're sorted in the order they were given.
BLAST results are sorted in descending order by the "Bit score" (column 12). Sorting by E. value would not be quite the same.
i.e. column 11 (E value) will be presented in ascending order (an E. value of 0 is as good as it gets), and as it gets higher, your hit is getting worse. Likewise, bitscores are highest at the top and decrease down the list. A hit with E. value = 0 should have a high score. How high exactly depends on the hit length, so it's common to see many E values with 0.0, but the scores differ. Here's a snippet of one of my blast results:
PVCcif_ATCC43949  PAU_06042014  100.000  24180  0     0   1      24180  2233012  2257191  0.0  43606
PVCcif_ATCC43949  PAU_06042014  78.098   5593   1164  16  2765   8323   3219956  3225521  0.0  4547
PVCcif_ATCC43949  PAU_06042014  77.549   5394   1139  16  2948   8323   2376163  2370824  0.0  4260
PVCcif_ATCC43949  PAU_06042014  74.140   5669   1342  34  8490   14100  3225688  3231290  0.0  3573
PVCcif_ATCC43949  PAU_06042014  76.835   3665   825   9   16507  20150  3234191  3237852  0.0  2764
PVCcif_ATCC43949  PAU_06042014  76.871   2646   571   11  16625  19248  2362597  2359971  0.0  2001
PVCcif_ATCC43949  PAU_06042014  70.733   4131   1005  51  9047   13063  2368869  2364829  0.0  1954
PVCcif_ATCC43949  PAU_06042014  75.040   1891   448   7   8500   10383  2370771  2368898  0.0  1272
PVCcif_ATCC43949  PAU_06042014  78.986   1361   269   7   22766  24119  890051   891401   0.0  1150
PVCcif_ATCC43949  PAU_06042014  77.819   1082   217   5   536    1610   3217983  3219048  0.0  868

As you can see, they all have E values of 0.0, but they differ quite a lot in score, this is because of the length of the alignment, and the different percentage IDs.

This is true only if there is a single query sequence. If there are multiple queries, the output is reported in the order of the query sequences, with bitscore used to sort the hits for each query.

BLASTN should definitely report alignments in the order that the query sequences were provided.
Based on your description and my personal experience, I think the most likely explanation is some kind of mix up. It happens to all of us, especially when we're in a rush to get an answer, but even sometimes when we're being careful and disciplined.
If you can go back and recreate the result, I would be very interested to see the two pairs of input and output files. If not, I wouldn't spend too much more time worrying about it. :-)
UPDATE Alignments are reported in the order that the query sequences are provided, even in tabular mode. The order of each query's hits is sorted by bitscore or e-value, but that is a secondary sorting. I've included an example below that has two queries, seq1 and seq2 with -outfmt 6.
seq1    scaffold2       100.000 280     0       0       1       280     11131   11410   7.73e-146       518
seq1    scaffold294     88.693  283     25      4       1       279     12160   12439   6.46e-92        339
seq1    scaffold1525    85.315  286     29      9       1       277     27142   26861   3.07e-75        283
seq1    scaffold1874    83.505  291     32      10      1       280     8558    8273    1.86e-67        257
seq1    scaffold147     82.653  294     34      9       1       279     10797   10506   1.45e-63        244
seq1    scaffold2478    82.374  278     38      8       11      280     7622    7896    1.13e-59        231
seq1    scaffold3405    87.940  199     17      5       1       196     14390   14584   1.46e-58        228
seq1    scaffold395     83.060  183     23      4       106     280     12021   11839   5.40e-38        159
seq1    scaffold2853    85.965  114     13      3       1       113     20864   20753   9.16e-26        119
seq2    scaffold9       100.000 140     0       0       1       140     1121    1260    2.34e-68        259
seq2    scaffold2950    90.909  110     6       3       1       106     110     1       6.84e-34        145
seq2    scaffold3416    95.385  65      3       0       1       65      11882   11818   1.16e-21        104
seq2    scaffold3103    95.385  65      3       0       1       65      9736    9672    1.16e-21        104
seq2    scaffold5297    88.608  79      7       1       64      140     150     228     6.99e-19        95.3

"BLASTN should definitely report alignments in the order that the query sequences were provided."
This is wrong. The OP is using BLAST tabular output so it's sorted by Bitscore.This might be true when using the alignment/pairwise output format, but I can't remember exactly off the top of my head.

@JoeHealey If that's the reason you downvoted, please read the update to my answer and reconsider.

Ah right you are indeed. Rectified the vote

I have analyzed the situation more and found that blastn automatically sorts by e-value. http://ix.io/xLh (column 11)

If your input file has 2 queries ("contig1" and "contig2"), the output will contain all of the "contig1" hits sorted by e-value, followed by all of the "contig2" hits sorted by e-value. Unless I'm mistaken...