text:
probeset to probeset mappings between Affymetrix arrays
I am interested in identifying mappings between different types of Affymetrix arrays. I am aware that mappings between gene and probeset can be extracted using Ensembl's Biomart database.
Ensembe gene id ENSG00000181019 maps to
1. AFFY HG-U133_Plus_2's 210519_s_at
2. AFFY HuGene-1_0-st-v1's 8002303

Is there any way to extract mappings of probeset ids between two arrays(eg: HG-U133_Plus_2, HuGene-1_0-st-v1)?
eg: 210519_s_at and 8002303

Is your question how to find overlapping probe sets between arrays?  Or how to identify probe sets in general?  It is not clear.

I have added an example. Does this make sense now?

So basically you can retrieve two gene-id mappings from Ensembl's Biomart. Correct? What is preventing you from joining these two tables based on the Ensembl gene id? You can use merge function in R or JOIN statement in SQL.

@IakovDavydov or a Unix join would work too

The problem is: All mappings are not exactly one-to-one. Either probesets of HG-U133_Plus_2 or HuGene-1_0-st-v1 mapping to one ensembl gene id.

@Prradep if you explain what the underlying problem is that you're trying to solve, we might be able to offer better advice

You should not necessarily expect to see one-to-one mappings, see my answer below.

If your question is: can probeset IDs from different platforms be mapped to one another in a similar way as mapping probesets to genes, then the answer is: Yes. BioMart allows you to map almost anything that has an ID to anything else that has an ID.
You can use BioMart either via the web interface or programatically. A brief guide to using the web interface, for mapping HG U133 Plus 2 to HuGene 1.0:

Go to the start page
Select H. sapiens Ensembl genes for your database; Human genes (GRCh38.p10) for your dataset
Click Filters in the left-hand column
Expand "REGION", scroll down to GENE and select "Input microarray probes/probesets ID list [Max 500 advised]"
Select AFFY HG U133 PLUS 2 probe ID(s) and either copy/paste or upload a list, one per line
Click Attributes in the left-hand column
Scroll through, uncheck what you don't want to see and choose what you do, for example Gene Stable ID, AFFY HuGene 1 0 st v1 probe and AFFY HG U133 Plus 2 probe
Finally click "Results" in the menu at the top of the left-hand column

And you should see a result like this (you'll need to click the "Results" button).
You should not expect that there be a one-to-one mapping, for two reasons:

Genes have multiple transcripts and probesets are mapped to each transcript
Some transcripts have more than one probeset: in this case, the HuGene IDs 8002301 and 8002303 map to transcripts for this gene

Finally: here's a programmatic example using R/BioMart:
library(biomaRt)
mart.hs <- useMart("ensembl", "hsapiens_gene_ensembl")
results <- getBM(attributes = c("ensembl_gene_id", "affy_hugene_1_0_st_v1", "affy_hg_u133_plus_2"),
                 filters = "affy_hg_u133_plus_2",
                 values = c("210519_s_at"),
                 mart = mart.hs)

results
  ensembl_gene_id affy_hugene_1_0_st_v1 affy_hg_u133_plus_2
1 ENSG00000181019               8002301         210519_s_at
2 ENSG00000181019               8002303         210519_s_at

Instead of biomaRt, it is also possible to use the mapping databases built into Bioconductor itself, and map from probe to gene, and then from gene to probe in the second. Some R code to convert between hgu133 and hgu95 using the same probe ID provided in another:
library(hgu133plus2.db)
library(hgu95av2.db)

query_probe <- "210519_s_at"

hgu133_ensembl <- select(hgu133plus2.db, keys = query_probe, columns = "ENSEMBL")

ensembl_hgu95 <- select(hgu95av2.db, keys = hgu133_ensembl$ENSEMBL, keytype = "ENSEMBL", columns = "PROBEID")

dplyr::inner_join(hgu133_ensembl, ensembl_hgu95, by = "ENSEMBL", suffix = c(".133", ".95"))

Other answers explain why there might not be one to one mapping between the probes.
The AbsID database does conversion based on mapping the probe sequences to a genome build, and then determines mappings based on overlapping genome alignment coordinates. This is really useful if you want to be sure that two probes are actually likely measuring the same transcript.
It is dependent on both probes aligning to the genome, however.
Disclaimer: I worked in the group that provides the AbsID mapping tool.

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


text:
How can the cell line contribution be estimated from RNASeq data?
Using a laser-capture microdissection of cells a group of cells stained with the marker of interest was sequenced. In another cohort of patients (this is all  human liver tissue) the whole tissue was sequenced (RNA-seq in both cases)
Can I estimate the contribution of the cells marked in the whole liver ("weight of these cells" in the liver in words of my PI)?
My gut feeling is that it can't be done this way, it would require both single cell sequencing and whole tissue sequencing to estimate the contribution of each cell line. But perhaps there is a tool that given the cell lines or the main expression of other cells lines it can be compared to using GSVA or some similar tool.

Have you looked at tools for admixture estimation (that's what you're doing)? How precise do you need the estimate to be?

I haven't heard of admixture estimation, so I haven't look for tools with that keyword. I don't have a requirement of precision, the more, the better :D. But I suspect my data is not too much good (I have just 6 technical replicates of the laser microdissection) so I can't expect much.

There are couple of computational methods which try to do this (I never used them, so no experience):

CellMix, based on sets of marker gene lists
Subset Prediction from Enrichment Correlation, which is based on correlations with subset-specific genes across a set of samples.
Cell type enrichment, which uses our highly expressed, cell specific gene database
Cell type-specific significance analysis using differential gene expression for each cell type

You might have to get some reference expression levels from public databases or papers for some of the methods.
One thing to keep in mind: you cannot really compute the cells proportion, only RNA proportion. If you have a good reason to assume that RNA quantity
per cell is very similar, this is a good proxy for the cells proportion in a tissue.

Nice references. I'll have that limitation in mind. Thanks

Cell deconvolution is mentioned in this Biostars post, which mentions CIBERSORT for immune cell mixes, and the Bioconductor package DeconRNASeq.
As far as I'm aware, it is only possible at best to get proportional representation for transcript expression from standard high-throughput sequencing results, because the sequencers and sample preparation workflow are designed in such a way that the same number of reads are output regardless of the input amount.

CIBERSORT sounds like a nice tool, worth a try.

You can download public single-cell RNA-seq data or bulk RNA-seq for purified primary cells for reference and then use CIBERSORT or MuSiC to deconvolve.
One thing to notice: to make bulk RNA-seq and single-cell reference compatible, TPM from RNA-seq can be compared with CPM in 3' (10X,Drop-seq etc.) single-cell data and TPM in full-length single-cell data (SMARTseq etc.)

Actually a month ago, I was searching for the same question and I found a few packages that can be used. First of all the functionality of them can be cell proportion estimation, cell specific gene expression estimation, or both of these functions. Secondly, their primary input, as expected, is bulk RNA-seq data and some also want a single cell RNA-seq data as secondary input.
The most popular one is the child of famous cibersort, CIBERSORTX. I went through their website and read the tutorials. They have several example datasets that can be used for learning their platform. However, the two biggest flaws were that they haven't shared their source code which means we can't reproduce the result on our own systems and we have to upload all of our data to Stanford servers (security concerns!). The second problem is that there is no API access to their server and everything should be done online which means very low compatibility with other tools.
The idea, using deconvolution and deep learning to infer subpopulation of cells from bulk RNA-seq, is great. So, I searched to see if others have also done something similar and maybe even better. I found three other tools:

MuSiC:   https://www.nature.com/articles/s41467-018-08023-x
Published in nature communication and a few months older than cibersortx. I haven't studied their algorithm yet and so I can't say if their method is better than cibersortx or not. However, the part that makes the difference is that all of their works, including their source code, is open-source. So, I can access their code and we don't need to upload our unpublished data into some not safe online server (unlike cibersortx). The second good thing is that it is implemented in R language which I use for most of my RNA-seq data analysis, so it is easily compatible with other tools.

Scaden : https://advances.sciencemag.org/content/6/30/eaba2619.full
This one is also open-source and accessible easily. It is written in Python and runs in the shell. So, it is somewhat compatible with other tools (not as compatible as MuSic). We again can process our data locally and don't need to upload our data to any servers (which is good). They claim that their method has a better correlation with the real scRNA-seq data than both cibersortx and MuSic (not always!!). Running their tools seems easy based on their tutorials (I haven't run it on my system yet!). However, their tutorials seem to be a little confusing and it seems like it is still in pre-release version and they have bugs!!

CDSeq:
https://github.com/kkang7/CDSeq_R_Package
Finally, I found the CDSeq package which doesn't need any scRNA-seq input to find single cell proportions and cell-specific gene expression in bulk data. It just gets your buld counts as input and info about very few parameters and does the whole thing automatically, and suprisingly good, I checked! If you could find gene expression pattern of your specific cells and also use it as an input, it is optional, it will give you even a better estimation.

source: Stack Exchange

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

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

sort: votes


text:
Normalizing RNAseq for PCA and CCA
Usually the expression data is transformed to log space using either RPKM, FPKM or CPM, this is required when looking for differential expression because the data is tested against the normal distribution(limma) or the negative bionimal distribution (DESeq2).
In weighted gene co-expression analysis (WGCNA) there is also the recommendation to normalize the data (see FAQ 4) but says "RPKM, FPKM, or simply normalized counts doesn't make a whole lot of difference for WGCNA" because it is based in correlations.
For PCAs the data should be on the same scale to avoid missing the scale as a factor. See this for more references. But should it be normalized too to look like a normal distribution using logarithms?
How does this extend to other methods like CCA? Should the data look like a normal distribution for CCA? As it is based on correlation I don't expect much changes, so I am using raw counts (RSEM output), but the examples I've seen use normalized data.
This question comes after reading a comment

I've given an answer about using how I normalise for PCA here. PCA works best if the data are normal-ish, but as models go it is fairly robust to fairly large deviations from normal. Transforming is not essential, but if data can be easily transformed into something that looks a bit more normal, then conclusions derived from looking at the PCA should be a bit more robust.
I'm not familiar with canonical-correlation analysis (CCA), but because I see sums and square roots in the formulas, I would assume it has similar data cleaning requirements to PCA. Unfortunately the R documentation for cancor doesn't have any information about any assumptions made on the data.


text:
Accuracy of the original human DNA datasets sequenced by Human Genome Project?
The Human Genome Project was the project of 'determining the sequence of nucleotide base pairs that make up human DNA, and of identifying and mapping all of the genes of the human genome'. It was declared complete in 2003, i.e. 99% of the euchromatic human genome completed with 99.99% accuracy.
Are the datasets provided by HGP still accurate, or as accurate as was claimed in 2003?
Given the technology in the past (such as using old techniques), or any other reason (newer research studies), is it possible that the datasets are not as accurate as originally expected?

It should be noted that the reason why a technology is replaced by a newer one is not always because the old one was not accurate. Sometimes the newer technology can be a bit less accurate on some aspects, but so much cheaper that it is deemed worth switching. The sequencing and assembly techniques used to generate the data at the basis of the first release of the human genome may have been quite accurate. That said, using newer techniques to improve on this basis will of course make newer versions better.

The HGP developed the first "reference" human genome - a genome that other genomes could be compared to, and was actually a composite of multiple human genome sequences.
The standard human reference genome is actually continually updated with major and minor revisions, a bit like software. The latest major version is called GRCh38, was released in 2013, and has since had a number of minor updates.

Are the datasets provided by HGP still accurate?

Yes, in a sense, but we certainly have better information now. One way to measure the quality of the assembly is that the initial release from the HGP had hundreds of thousands of gaps - sequences that could not be resolved (this often occurs because of repetitive sequences). The newest reference genome has less than 500 gaps.

While the quality of the reference human assembly keeps improving, there are still misassemblies in it. A common problem is recent segmental duplications are occasionally collapsed into one sequence in the reference. Another issue is that the centromeric sequences in the reference are computationally generated, which are probably different from real sequences. Issues like these often complicate data analyses.
You should also beware that each human has a different genome. A large region having one copy in the reference genome may have two copies in a specific sample. While this is not really the problem with the reference, such copy-number changes will have the same effect as reference errors and mess up your pipelines.
There is still room for improvement to the human reference genome. In some regions, the CHM1 and CHM13 PacBio assemblies are better than the current reference genome at the larger scale. Illumina population data can produce better consensus at the single base level. GRC is continuously releasing new patches to the latest assembly.


text:
Plotting gene models - alternatives to ggbio
I am curious how people make genome track plots such as what can be done with the Bioconductor package ggbio.
For example this plot was made w/ ggbio:

I think ggbio is great for making fast plots, but it is kind of a pain to make publication quality figures because it is hard to manipulate the plots after they are constructed. I have noticed that ggbio constructed plots do not return a ggplot or gtable object. I can access the grobs from a ggbio tracks object, but they won't be aligned as they would be in the final plot. I like to make shared legends and include them in the plot, and without having a gtable or ggplot object it is hard to add the shared legends in, and I end up having to save them as a separate file which is tedious. FYI I mainly use ggbio for gene tracks + x-axis alignment.
Normally I would use ggplot2 to customize the plot extensively, but it doesn't work very well after plotting with ggbio.
How can I use ggbio for these types of plots? Is there a better tool instead?
EDIT:
I would get the legends for each individual plot and combine them before adding back to the main plot using something like this:
leg1.grob = cowplot::get_legend(plot1)
leg2.grob = cowplot::get_legend(plot2)
sharedleg.grob = gridExtra::gtable_rbind(leg1,leg2)

mainplot.grob = ggplotGrob(main_plot)
mainplot.grob = gtable_add_cols(mainplot.grob,
                                widths = sum(sharedleg.grob$widths))
mainplot.grob = gtable_add_grobs(mainplot.grob,
                                 grobs = sharedleg.grob,
                                 t=t,b=b,l=l,r=r)
cowplot::ggdraw(mainplot.grob)

I've used gviz before, but long ago. Do you know if you can extract the grid plot grobs of the gene track plots with Gviz?

It seems like the ggbio object already contains the corresponding ggplot object.
You should be able to extract the ggplot from the ggbio_obj@ggplot slot and extract your gtables.
Small example:
library(ggbio)
library(ggplot2)
library(patchwork)

fl.bam <- system.file("extdata", "wg-brca1.sorted.bam", package = "biovizBase")
wh <- as(c("chr17:41239394-41319151:+"), "GRanges")
p_cov <- autoplot(fl.bam, which = wh)
#> reading in as Bamfile
#> Parsing raw coverage...
#> Read GAlignments from BamFile...
#> extracting information...

df <-  data.frame(x= sample(1:100,size = 100,replace = TRUE))
p_hist <- ggplot(df) +
    geom_histogram(aes(x))

p_hist / p_cov
#> Error: Can't add `e2` to a ggplot object.
#> Backtrace:
#>     █
#>  1. └─patchwork:::`/.ggplot`(p_hist, p_cov)
#>  2.   └─GGally:::`+.gg`(e1, e2)
#>  3.     └─e1 %+% e2
#>  4.       └─ggplot2:::add_ggplot(e1, e2, e2name)
#>  5.         ├─ggplot2::ggplot_add(object, p, objectname)
#>  6.         └─ggplot2:::ggplot_add.default(object, p, objectname)

p_hist / p_cov@ggplot
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Created on 2021-02-05 by the reprex package (v1.0.0)
Session info

sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 3.6.0 (2019-04-26)
#>  os       Red Hat Enterprise Linux
#>  system   x86_64, linux-gnu
#>  ui       X11
#>  language (EN)
#>  collate  en_US.UTF-8
#>  ctype    en_US.UTF-8
#>  tz       Europe/Berlin
#>  date     2021-02-05
#>
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package              * version  date       lib source
#>  acepack                1.4.1    2016-10-29 [1] CRAN (R 3.6.0)
#>  AnnotationDbi          1.48.0   2019-10-29 [1] Bioconductor
#>  AnnotationFilter       1.10.0   2019-10-29 [1] Bioconductor
#>  askpass                1.1      2019-01-13 [1] CRAN (R 3.6.0)
#>  assertthat             0.2.1    2019-03-21 [1] CRAN (R 3.6.0)
#>  backports              1.1.8    2020-06-17 [1] CRAN (R 3.6.0)
#>  base64enc              0.1-3    2015-07-28 [1] CRAN (R 3.6.0)
#>  Biobase                2.46.0   2019-10-29 [1] Bioconductor
#>  BiocFileCache          1.10.2   2019-11-08 [1] Bioconductor
#>  BiocGenerics         * 0.32.0   2019-10-29 [1] Bioconductor
#>  BiocManager            1.30.10  2019-11-16 [1] CRAN (R 3.6.0)
#>  BiocParallel           1.20.1   2019-12-21 [1] Bioconductor
#>  biomaRt                2.42.1   2020-03-26 [1] Bioconductor
#>  Biostrings             2.54.0   2019-10-29 [1] Bioconductor
#>  biovizBase             1.34.1   2019-12-04 [1] Bioconductor
#>  bit                    1.1-15.2 2020-02-10 [1] CRAN (R 3.6.0)
#>  bit64                  0.9-7    2017-05-08 [1] CRAN (R 3.6.0)
#>  bitops                 1.0-6    2013-08-17 [1] CRAN (R 3.6.0)
#>  blob                   1.2.1    2020-01-20 [1] CRAN (R 3.6.0)
#>  BSgenome               1.54.0   2019-10-29 [1] Bioconductor
#>  checkmate              2.0.0    2020-02-06 [1] CRAN (R 3.6.0)
#>  cli                    2.0.2    2020-02-28 [1] CRAN (R 3.6.0)
#>  cluster                2.0.8    2019-04-05 [2] CRAN (R 3.6.0)
#>  colorspace             1.4-1    2019-03-18 [1] CRAN (R 3.6.0)
#>  crayon                 1.3.4    2017-09-16 [1] CRAN (R 3.6.0)
#>  curl                   4.3      2019-12-02 [1] CRAN (R 3.6.0)
#>  data.table             1.12.8   2019-12-09 [1] CRAN (R 3.6.0)
#>  DBI                    1.1.0    2019-12-15 [1] CRAN (R 3.6.0)
#>  dbplyr                 1.4.2    2019-06-17 [1] CRAN (R 3.6.0)
#>  DelayedArray           0.12.2   2020-01-06 [1] Bioconductor
#>  dichromat              2.0-0    2013-01-24 [1] CRAN (R 3.6.0)
#>  digest                 0.6.25   2020-02-23 [1] CRAN (R 3.6.0)
#>  dplyr                  0.8.5    2020-03-07 [1] CRAN (R 3.6.0)
#>  ellipsis               0.3.1    2020-05-15 [1] CRAN (R 3.6.0)
#>  ensembldb              2.10.2   2019-11-20 [1] Bioconductor
#>  evaluate               0.14     2019-05-28 [1] CRAN (R 3.6.0)
#>  fansi                  0.4.1    2020-01-08 [1] CRAN (R 3.6.0)
#>  farver                 2.0.3    2020-01-16 [1] CRAN (R 3.6.0)
#>  foreign                0.8-71   2018-07-20 [2] CRAN (R 3.6.0)
#>  Formula                1.2-3    2018-05-03 [1] CRAN (R 3.6.0)
#>  fs                     1.4.0    2020-03-31 [1] CRAN (R 3.6.0)
#>  GenomeInfoDb           1.22.1   2020-03-27 [1] Bioconductor
#>  GenomeInfoDbData       1.2.2    2020-04-01 [1] Bioconductor
#>  GenomicAlignments      1.22.1   2019-11-12 [1] Bioconductor
#>  GenomicFeatures        1.38.2   2020-02-15 [1] Bioconductor
#>  GenomicRanges          1.38.0   2019-10-29 [1] Bioconductor
#>  GGally                 1.5.0    2020-03-25 [1] CRAN (R 3.6.0)
#>  ggbio                * 1.34.0   2019-10-29 [1] Bioconductor
#>  ggplot2              * 3.3.2    2020-06-19 [1] CRAN (R 3.6.0)
#>  glue                   1.4.1    2020-05-13 [1] CRAN (R 3.6.0)
#>  graph                  1.64.0   2019-10-29 [1] Bioconductor
#>  gridExtra              2.3      2017-09-09 [1] CRAN (R 3.6.0)
#>  gtable                 0.3.0    2019-03-25 [1] CRAN (R 3.6.0)
#>  highr                  0.8      2019-03-20 [1] CRAN (R 3.6.0)
#>  Hmisc                  4.4-0    2020-03-23 [1] CRAN (R 3.6.0)
#>  hms                    0.5.3    2020-01-08 [1] CRAN (R 3.6.0)
#>  htmlTable              1.13.3   2019-12-04 [1] CRAN (R 3.6.0)
#>  htmltools              0.5.0    2020-06-16 [1] CRAN (R 3.6.0)
#>  htmlwidgets            1.5.1    2019-10-08 [1] CRAN (R 3.6.0)
#>  httr                   1.4.1    2019-08-05 [1] CRAN (R 3.6.0)
#>  IRanges                2.20.2   2020-01-13 [1] Bioconductor
#>  jpeg                   0.1-8.1  2019-10-24 [1] CRAN (R 3.6.0)
#>  knitr                  1.29     2020-06-23 [1] CRAN (R 3.6.0)
#>  labeling               0.3      2014-08-23 [1] CRAN (R 3.6.0)
#>  lattice                0.20-38  2018-11-04 [2] CRAN (R 3.6.0)
#>  latticeExtra           0.6-29   2019-12-19 [1] CRAN (R 3.6.0)
#>  lazyeval               0.2.2    2019-03-15 [1] CRAN (R 3.6.0)
#>  lifecycle              0.2.0    2020-03-06 [1] CRAN (R 3.6.0)
#>  magrittr               1.5      2014-11-22 [1] CRAN (R 3.6.0)
#>  Matrix                 1.2-17   2019-03-22 [2] CRAN (R 3.6.0)
#>  matrixStats            0.56.0   2020-03-13 [1] CRAN (R 3.6.0)
#>  memoise                1.1.0    2017-04-21 [1] CRAN (R 3.6.0)
#>  mime                   0.9      2020-02-04 [1] CRAN (R 3.6.0)
#>  munsell                0.5.0    2018-06-12 [1] CRAN (R 3.6.0)
#>  nnet                   7.3-12   2016-02-02 [2] CRAN (R 3.6.0)
#>  openssl                1.4.1    2019-07-18 [1] CRAN (R 3.6.0)
#>  OrganismDbi            1.28.0   2019-10-29 [1] Bioconductor
#>  patchwork            * 1.0.0    2019-12-01 [1] CRAN (R 3.6.0)
#>  pillar                 1.4.3    2019-12-20 [1] CRAN (R 3.6.0)
#>  pkgconfig              2.0.3    2019-09-22 [1] CRAN (R 3.6.0)
#>  plyr                   1.8.6    2020-03-03 [1] CRAN (R 3.6.0)
#>  png                    0.1-7    2013-12-03 [1] CRAN (R 3.6.0)
#>  prettyunits            1.1.1    2020-01-24 [1] CRAN (R 3.6.0)
#>  progress               1.2.2    2019-05-16 [1] CRAN (R 3.6.0)
#>  ProtGenerics           1.18.0   2019-10-29 [1] Bioconductor
#>  purrr                  0.3.4    2020-04-17 [1] CRAN (R 3.6.0)
#>  R6                     2.4.1    2019-11-12 [1] CRAN (R 3.6.0)
#>  rappdirs               0.3.1    2016-03-28 [1] CRAN (R 3.6.0)
#>  RBGL                   1.62.1   2019-10-30 [1] Bioconductor
#>  RColorBrewer           1.1-2    2014-12-07 [1] CRAN (R 3.6.0)
#>  Rcpp                   1.0.4.6  2020-04-09 [1] CRAN (R 3.6.0)
#>  RCurl                  1.98-1.1 2020-01-19 [1] CRAN (R 3.6.0)
#>  reprex                 1.0.0    2021-01-27 [1] CRAN (R 3.6.0)
#>  reshape                0.8.8    2018-10-23 [1] CRAN (R 3.6.0)
#>  reshape2               1.4.4    2020-04-09 [1] CRAN (R 3.6.0)
#>  rlang                  0.4.6    2020-05-02 [1] CRAN (R 3.6.0)
#>  rmarkdown              2.1      2020-01-20 [1] CRAN (R 3.6.0)
#>  rpart                  4.1-15   2019-04-12 [2] CRAN (R 3.6.0)
#>  Rsamtools              2.2.3    2020-02-23 [1] Bioconductor
#>  RSQLite                2.2.0    2020-01-07 [1] CRAN (R 3.6.0)
#>  rstudioapi             0.11     2020-02-07 [1] CRAN (R 3.6.0)
#>  rtracklayer            1.46.0   2019-10-29 [1] Bioconductor
#>  S4Vectors              0.24.3   2020-01-18 [1] Bioconductor
#>  scales                 1.1.1    2020-05-11 [1] CRAN (R 3.6.0)
#>  sessioninfo            1.1.1    2018-11-05 [1] CRAN (R 3.6.0)
#>  stringi                1.4.6    2020-02-17 [1] CRAN (R 3.6.0)
#>  stringr                1.4.0    2019-02-10 [1] CRAN (R 3.6.0)
#>  styler                 1.3.2    2020-02-23 [1] CRAN (R 3.6.0)
#>  SummarizedExperiment   1.16.1   2019-12-19 [1] Bioconductor
#>  survival               3.1-11   2020-03-07 [1] CRAN (R 3.6.0)
#>  tibble                 3.0.0    2020-03-30 [1] CRAN (R 3.6.0)
#>  tidyselect             1.0.0    2020-01-27 [1] CRAN (R 3.6.0)
#>  VariantAnnotation      1.32.0   2019-10-29 [1] Bioconductor
#>  vctrs                  0.2.4    2020-03-10 [1] CRAN (R 3.6.0)
#>  withr                  2.4.1    2021-01-26 [1] CRAN (R 3.6.0)
#>  xfun                   0.15     2020-06-21 [1] CRAN (R 3.6.0)
#>  XML                    3.99-0.3 2020-01-20 [1] CRAN (R 3.6.0)
#>  xml2                   1.3.0    2020-04-01 [1] CRAN (R 3.6.0)
#>  XVector                0.26.0   2019-10-29 [1] Bioconductor
#>  yaml                   2.2.1    2020-02-01 [1] CRAN (R 3.6.0)
#>  zlibbioc               1.32.0   2019-10-29 [1] Bioconductor

@Alex Does this solution work when aligning multiple genome tracks using ggbio as in my original example?  Either way, useful info for others. My original question was actually aiming to get feedback on how others make these types of plots rather than seeking a solution to a specific problem in my workflow, so if you have any thoughts on that I'd be interested

@Reilstein this solution requires more work to align multiple genome tracks based on coordinates, so actually does not really answer your initial question.
I was using Gviz before, but I like ggbio more due to the ease of use and how it integrates with ggplot2.

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


text:
Is there an established method for comparing expression of groups of genes (gene sets)?
I have a figure in a paper where I show the log2(fold-change) (obtained with DESeq2) between two groups (based on genotype) for genes within a specific hallmark gene set, and a reviewer is asking whether the differences are statistically significant.
I have not been able to find any methods for comparing expression of groups of genes like this, and i'm sure there is a method or a reason why it isn't done that I'm not aware of.
To be extra clear, what i have is two sets of log2FC values for a set of genes (one for one genotype (drug-treated vs. untreated) and the same for a second genotype), and I'm looking for a statistical comparison of these two groups of Log2FC measurements.
Thanks

They want you to find genes whose response to the drug is different based on the genotype? If so, you can add an interaction term to your DESeq2 regression formula ~ condition + genotype + condition:genotype. Refer to the DESeq2 documentation and help("results", package="DESeq2") for more information on multifactor design and interaction terms.

Okay yeah that makes sense for how to obtain differences in the response, but then how might I statistically compare these differences in response for a group of genes (I.e a gene set)? Would you just do something simple like a t test / Wilcox test?

You could apply the GSVA transformation to later compare them. You would need to give a step back and start from the expression values not log2FC, but I think it would be better than using the log2FC directly.
If you need to use them you could compare the gene set enrichment value of these two sets on the same comparison. I would use fgsea with the log2FC and the two sets of genes. If the normalized enrichment value (NES) is different then you could explain that these sets expression are indeed different.


text:
Annotation format design
Bashing file formats is a favorite pastime in bioinformatics, and annotation file formats such as GFF and BED seem to get special attention. A lot of this frustration stems from community's shockingly inconsistent adherence to specifications and conventions, but there are also some (dare I say objectively) problematic design choices in each of these formats.

GFF (and its more common derivatives GTF and GFF3) use 1-based closed interval notation, which optimizes for human comprehension but is far inferior to 0-based half-open interval notation (such as used by BED) for computations involving interval arithmetic.
Although BED and GTF were designed for very specific use cases (visualization and gene prediction, respectively), they have been used and abused in a much wider set of contexts. For example, the BED fields related to the thick part are irrelevant if you're not plotting them in a genome browser.
BED supports a single level of feature decomposition (a feature can be broken up into blocks). GTF supports two levels (exons grouped by transcript_id, transcripts grouped by gene_id). In contrast, GFF3 supports an arbitrary number of levels, and uses parent/child relationships defined by ID and Parent attributes to declare a directed acyclic graph of features.
Data that does not fit into mandatory pre-defined fields must be relegated to optional fields or free-form attribute key/value pairs. While this flexibility is powerful, a common complaint is that "all the action" happens in these optional/free-form fields.
There is a dearth of validation tools, and those that do exist focus primarily on validating syntax and not semantics. To use an aging analogy, it's one thing to say an XML file is valid, but it's completely different to validate it against a schema. There are essentially no widely used tools that do the latter for annotation files.

If we were tasked with creating a new annotation format, and if we were guaranteed the resources needed to develop it, and interest and wide adoption from the wider community (one can dream!), what design criteria should be considered in the development of this new format? What, if anything, makes an objectively good annotation data format?

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

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

I quite like BED and GFF3 (I don't like GTF/GFF2, though). As text-based formats, I don't think they leave us much room for improvement. Anyway, if you want a new format, here is one. The following is a hybrid between GFF3 and BED. It is a TAB-delimited text-based format with the following fields:

chr, required
start (0-based), required
end, required
strand
ID
type (see below)

Like BED, only the first 3 columns are required, the rest are optional. What is different starts with the "type" field as in GFF. This "type", when present, defines the columns following it. For example, if type==coding, we could have cdsStart, cdsEnd, blockCount, blockSizes and blockStarts like BED; if type==exon, column 7 could be a "phase", indicating the phase of the first base. This way, "type" makes this format highly extensible while still relatively easy to parse in comparison to using optional tags all the way. In addition, we may have semi-colon-separated "key"="value" pairs at the end of each line as in GFF3. Example:
chr1   10000  50000  -  x1 coding   10100   40800  2  1000,2000   10000,48000
chr1   10000  11000  -  *  exon     *       foo1=bar1;foo2=bar2
chr1   48000  50000  -  *  exon     2
chr2   10000  50000

The above only gives a bare bone of the format. There are subtle questions such as: 1) the uniqueness/scope of ID; 2) whether to have parent ID as a fixed field or a type-specific field; 3) where to put display name; 4) sorting order. These are also important in practice.
PS: I like SQL a lot and think it should get used more often in bioinformatics, but I don't think it replaces text formats completely. Formats are useful for serialization and data transfer. They require less skill to work with and less software/hardware resource to deploy. Carefully engineered binary representations of formats can be much more efficient in a lot of use cases.

The problem is that GFF, fundamentally, is a relational format: it provides tags that relate individual rows via one-to-many relationships (e.g. gene–exon). This indirectly highlights the second complication: individual rows have different types, and therefore store different attributes in the 9th column.
Over the last few decades (!), we have accumulated a wealth of theory and tools to work with this kind of data. And the usual solution is to create a database schema for a relational database, and to use database drivers and database query languages (e.g. SQL, but increasingly also data-relational mappers such as LINQ and dplyr).
Using a text-based format is attractive for many reasons that Devon mentions, but it is fundamentally at odds with a lot of the theory and tools for relational data. This creates an impedance mismatch.
I’m convinced that, in the long run, the solution is going to be to revert to using relational databases for complex annotations. I say “revert” even though these databases already exist, they are simply often ignored in bioinformaticians’ day-to-day work (I never use them). Because this is objectively the best technical solution, and we have the research to back this up.

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

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

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

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

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

Text-based: It's absurdly common to want to use grep or awk on annotations. Sure, one could make variants of these that are binary-format aware, but why reinvent the wheel. Of course text files don't innately allow region-based queries of their contents, so on to point 2...

This should further explicitly be line-based. Fields would need to be tab-separated (ideally it'd use the ascii record separator character instead, but I fear that ship has long ago sailed).

A strictly defined line order: This is one of my personal pet-peeves about BED/GFF/GTF files. If you're going to make a text-based annotation format, everyone would end up ahead if said format made explicit that it should be sorted. This then allows things like tabix to be used so the "query a region" problem is then solved. But I would go further than that. The issue with things like GFF is that there are multiple inter-dependent lines and there's no strict rule about whether a parent line absolutely must come before a child line. This just makes implementing things a nightmare and more often than not a randomly ordered file will just break tools. Since GTF or GFF-style interdependent lines likely be how any annotation format works, this ordering between lines with the same start position should be made explicit in the format.

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

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

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

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

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

Just throwing this out there, as all the best stuff has already been said, but why does one need to specify a strict specification for a bioinformatic data format? As you say, what typically happens is all the action ends up in the optional fields.
There's a lot to be said for the libertarian values of just letting people figure it out on their own. Take optional fields in the BAM tag spec. for example. Here you can have your own tags, but they must start with an "X", a "Y" or a "Z" and be followed by "[A-Za-z0-9]". Why? Why can't people use their own tag names like "read is unique" or "edit distance to genome" or whatever they want. Are we not to be trusted with the power of naming things and recalling things arbitrarily? The result is one must look up the tag name's actual meaning in the publication of the tool that produced it, or ask on certain forums, etc. And this is assuming one knows the tool that produced the reads - if not then who knows wtf 'XT' stands for.
Essentially, downstream readers of optional fields already trust that 'XT' is what it thinks it is. So if we're happy to use optional fields in any capacity, why stop at fields?
Extend this logic to the columns of the data format. Let users determine column names. You may once in a blue moon come across a file with the chromosome column as "chr" rather than "chromosome", but most of the time you'll be looking for "read is quantum-radioactive", and detecting and fixing naming-errors is not going to be as much of a problem as using a dataformat that does not store what it is you want to store. Or worse, it can store it, but in a really illogical way that results in everyone who uses the dataformat to have to ask at least 3 questions here or on other forums before they understand whats really going on.
So you end up with a general solution to the problem of storing tabular information in a compatible way across many different types of software, also known as SQL (or any other kind of general-purpose database such as REDIS, Neo4J, Numpy, etc). In fact, does it even matter what the datastore is, so long as it's tabular and has the "chromosome" and "position" values for each item?
TL;DR - We will not think of tomorrow's best data format today, because the nature of tomorrow's data is not yet known. Less policing in this area would most likely result in more robust software, where nothing is taken for granted and no assumptions about the data can be made until the schema has been parsed.
In truth, the only reason we think we need to specify our dataformats is because we know if we didn't, other people will do something we didn't think of - and we, i think falsely, assume that will be bad.


text:
Tools for simulating Oxford Nanopore reads
Are there any free open source software tools available for simulating Oxford Nanopore reads?

What are you planning on simulating data for? In my opinion, given the wealth of publicly available real datasets, there is little need for simulation except in very specific use cases.

It's often useful to simulate a variant or some other type of genomic artifact to evaluate whether a tool or set of tools can correctly identify the artifact. This is much easier to do on simulated data where you know the exact answer, and is extremely valuable even if the simulated data isn't a perfect reflection of real data.

Simulators designed specifically for Oxford Nanopore:

NanoSim
NanoSim-H
SiLiCO
ReadSim
DeepSimulator

General long read simulators:

Loresim
Loresim 2
FASTQsim
LongISLND

For an exhaustive list of existing read simulators, see page 15 of my thesis, Novel computational techniques for mapping and
classifying Next-Generation Sequencing data.

By chance, just today I've heard of a nanopore read simulator, NanoSim. It is  released under a GPL license. I have never used it, though...

In addition to the already mentioned NanoSim, there is also SiLiCO and ReadSim (although it hasn't been updated in over 2 years, so I am not sure how relevant it is at this point considering how fast the technology is progressing).

The best nanopore read simulators would be associated with the best base-callers. For a base-caller to effectively model the DNA strand, it needs to take into account the expected underlying electrical model together with the associated signal noise (both in the time dimension as well as the amplitude dimension). In theory, it should be possible to reverse the algorithm and generate an electrical signal given an underlying sequence.
Unfortunately, I'm not aware of any tools that attempt to simulate nanopore reads at the electrical level. Any "nanopore read simulator" that concentrates only on base sequence would need to encompass all the possible base-calling software models that exist, which is an impossible task (particularly given how quickly ONT updates their own base callers).


text:
Determine reference for reference-compressed SRA file
Note: this question was also asked on Github.
I have 241 SRA files that appear to be reference compressed. I didn't even know this was a thing until I tried to convert them to Fastq files without an internet connection. I got the "name not found while resolving..." error described at the bottom of this page when I attempted to do the conversion.
The instructions say:

If you have elected to prevent the toolkit from contacting NCBI, you will need to manually acquire the reference(s) here

We have many NCBI databases mirrored locally, including RefSeq, so it's likely that we already have all the references. I have two questions:

How do I determine the reference used to compress a reference-compressed SRA file? (Note this needs to be scripted for 241 SRA files, so manual inspection is out of the question.)
Once I've determined the correct reference and ensured it is available locally, how to I tell fastq-dump where to find it?

Run align-info  to see the references used.
Run prefetch </path/to/sra/file> to download missed references.


text:
How are snakemake's --cluster and --drmaa options implemented?
I'm fairly new to snakemake and I'm trying to understand the difference between the --cluster and --drmaa flags, both of which are used to submit jobs to compute clusters/nodes.
The docs give a few hints about the advantages of using --drmaa here:

If your cluster system supports DRMAA, Snakemake can make use of that
  to increase the control over jobs. E.g. jobs can be cancelled upon
  pressing Ctrl+C, which is not possible with the generic --cluster
  support.

And here:

If available, DRMAA is preferable over the generic cluster modes
  because it provides better control and error handling.

So I have a conceptual understanding of the advantages of using --drmaa. However, I don't consider the above a very complete explanation, and I don't know how these flags are implemented in snakemake under the hood, can someone elaborate?
Note, that although this could be considered a more general programming question, snakemake is predominantly used in bioinformatics; I was persuaded that this question would be considered on-topic because of this meta-post and this answer.

Hi Chris! What's the actual question here? Seems that DRMAA is an open API describing submission and control of cluster/grid jobs. So if you have a cluster or submission system that has support for it (such as SLURM), snakemake will afford you the additional controls as supported by that API -- including as per your question, the ability to pass on a Ctrl+C at your console to the grid -- a scenario that would require use of the qdel command on Sun Grid Engine, for example. So if you've infrastructure with a scheduler that has DRMAA support, I can see why you might as well use --drmaa.

Maybe wait until we've actually poached the snakemake developers :)

The code is open source and fairly easy to follow, I found it quite easy to work out how DRMAA jobs were submitted by nosing around the source for 10 minutes. I'm not entirely sure what the question is looking for, though.

I'd always kind of wondered how this worked too, so I took this as an excuse to look into the snakemake code. At the end of the day this becomes a question of (1) how are jobs actually submitted and (2) how is it determined if jobs are done (and then whether they failed)?
For DRMAA, python has a module (appropriately named "drmaa") that wraps around the libdrmaa library that comes with most schedulers. This is a very popular route, for example the Galaxy project uses this for dealing with most clusters (e.g., I use it to connect our internal Galaxy instance to our slurm cluster). The huge benefit here is that DRMAA does magic to allow you to simply submit commands to your cluster without needing to know whether you should execute qsub or srun or something else. Further, it then provides methods to simply query whether a job is running or not and what its exit status was.
Using the --cluster command requires a LOT more magic on the side of snakemake. At the end of the day it creates a shell script that it then submits using the command you provided. Importantly, it includes some secret files in that script that it can then watch for (ever noticed that .snakemake directory where you execute a script? That appears to be among the things it's used for.). These are named {jobid}.jobfinished and {jobid}.jobfailed and one of them will get touched, depending on the exit status of your command/script. Once one of those is there, then snakemake can move on in its DAG (or not, if there's a failure). This is obviously a LOT more to keep track of and it's then not possible for snakemake to cancel running jobs, which is something one can do easily with DRMAA.

Having written a package for submitting jobs to Sun Grid Engine without DRMAA, can definitely confirm the magic must be strong here. Though instead of monitoring whether files have been touched, I chose to fire a POST request at the monitoring server. I can see that DRMAA must expose some nice API to get job status without having to do some odd workarounds.

Yup, it's more convenient than what I would have come up with too :)

DRMAA (Distributed Resource Management Application API) appears to be an open API that describes a specification for submission and management of work submitted to some grid/cluster. If your scheduler is DRMAA compliant, my guess would be that using snakemake's --drmaa flag will afford you the additional controls exposed by that API.
As mentioned in your question with DRMAA support, a Ctrl+C on one's console would be passed on to the grid to kill or stop submitted jobs. A task that would otherwise require use of the qdel command on Sun Grid Engine, for example.
I could see why the suggestion would be to enable --drmaa where possible as it's likely to expose functionality (such as control of jobs, and availability of logs and errors) that makes the submission and management of jobs a little easier for you and snakemake. Perhaps someone with snakemake experience can explain a little more, but I hope this helps?
For what it's worth, I use Sun Grid Engine without DRMAA and it's quite a pain to get good information out of qstat, and its qacct log format is quite possibly the worst file format I have ever encountered. I suspect DRMAA provides a nice (or at least more reasonable) API for snakemake (and others) to get at this information more readily.
Appended:
Under the hood, the flags are parsed by Python's argparse module in the __init__.py. A mutually exclusive group of options forces a selection of --cluster or --drmaa (or --cluster-sync).
--cluster or -c requires you to specify an argument; the command for which to submit a job on your cluster (the example names qsub). The --drmaa flag just seems to indicate to snakemake that DRMAA is to be used, which doesn't change a lot in terms of how the command you run from your console is processed in the __init__.py or workflow.py.
However, when it is time to interface with your scheduler (as specified in scheduler.py), an elif statement checks whether you're in some form of cluster mode and the else catches the case where you've raised the --drmaa flag instead.
Here's where the magic happens, as now your job is submitted with the DRMAAExecutor as specified in executors.py, instead of the default executor.
Without investigating much more, I can see that the DRMAAExecutor features some more class attributes and exposes additional functions when compared to the others.
Hope this is of more use? Again, I've only just looked at the code myself but in lieu of a snakemake user chiming in, I thought I'd have a nose.

Thanks, but this tells me what I already know; i have a conceptual understanding but i'd like the details on how the 2 different approaches are implemented in snakemake, i'll try to clarify the question a bit

Sorry I pre-empted what you might have been asking for as I didn't have much more room to write in the comments. I'll try and expand to address your edited question if I can help!

@Chris_Rands I've expanded, is this relevant?

Thanks, upvoted, sorry if I wasn't very clear with my question, I just thought this might be useful question for others,i was planning on answering myself if i didn't get any answers, but i've not delved into the source yet

@Chris_Rands No problem (and thanks)! I enjoy digging around code bases to see how they work. Between Devon and I, have we found what you were looking for? Happy to chat on this with you if not.


text:
How do I validate a single sample ArrayCGH result?
We have arrayCGH (aCGH) results for one sample. There is a 0.5 Mb terminal duplication on chromosome 19 (62995490-63407936, according to NCBI36/hg18). The duplication is rare: a literature review suggests there are only 3-4 samples with clinical information.
What are the steps to validate the results? How do we ascertain that this duplication is the cause of the clinical symptoms?
I have some ideas:

aCGH the parents. Not sure how this would help.
whole genome exome sequencing. Worried this might make it more difficult to pinpoint genetic cause.
whole genome sequencing?
other ideas?

Note: I am new to aCGH and high-throughput sequencing, any advice is welcome.

Is the duplication in the telomere or before it?

@DevonRyan I think, it is before, on the last band of chr19. Does it matter, sorry not being ignorant, just no clue about it.

Telomeric sequences are highly repetitive, so yes it matters :)

Telomere starts at the end of the last band, or is it part of chromosome?

If you look at the UCSC genome browser, the telomere is the GGGTTA simple repeat at the end. In others it's just a stretch of NNNN in the reference (no clue how this looks in hg18).

Just as a little suggestion, especially in the beta you should take the recommendations of the SE interface very seriously and always leave a comment when down-voting, or you demonstrate just the prejudice of hostility! Also be a little open minded, questions about omics data and how to validate omics results are typically very relevant. Otherwise you will notice that nothing can be asked anymore. Just saying, from and old bs mod. P.S.: I have the feeling that people down-vote because they do not understand, that's not always the fault of the Op.

Another observation, this question might be ethically questionable! You are putting forward info on a real case to diagnosis, possibly identifying information is revealed.

Why not just good old qPCR? That's (A) quick, (B) cheap and (C) easy to analyze. If you care about the exact location of the break point (I'm guessing from the context that you don't), then targeted sequencing with a custom capture kit would work.
Regarding validating the biological relevance of this, there are multiple (parallel) routes one can take. Firstly, screen unaffected family members for this same alteration. If you find this in unaffected individuals then it's obviously not the causative alteration. Ideally, one would also do either a cell-line experiment or a mouse (or other model organism) experiment to see if you can reconstitute at least some component of the clinical phenotype. This may not always be possible, of course.


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:
Extract nanopore read ID & start times from fastq file
I have a fastq file from minION (albacore) that contains information on the read ID and the start time of the read. I want to extract these two bits of information into a single csv file.
I've been trying to figure out a grep/awk/sed solution, but without success.
E.g.
@93a12f52-95e5-40c7-8c3e-70bf94ed0720 runid=17838b1d08f30a031bf60afabb146a8b0fba7486 read=12217 ch=492 start_time=2017-07-04T06:42:43Z
CTATTGTCCCCTGCCGCTGCCCCTCCTGCTACGCCCCACTGCCTCACCAGCCGTTACGGTCGCCCCCCATCGCATGCCTTTACACACACACTTCTTTACACATGCTATCTTCCC
+
""*)&$-.,-(#"&'$%''+16#"$##&)%%/"+*(*(&#"&'%1"+)#)""%$$#&&'1%"'8>MJ<#'&%'.2'.$(&#'()'&&%'$('%"%%%..$#"&"#+&,*$%"#"
@ff37e422-a25f-404c-8314-ef1733f9c30c runid=17838b1d08f30a031bf60afabb146a8b0fba7486 read=8432 ch=200 start_time=2017-07-04T06:56:41Z
CGATGGCCGTATGCTTTTGTTATGAAGCGAAAAGCTGCTCGCTTCTCTAGATAATAATGATGTGGCGAAAACGCTATGCGATTCGTTGACATACATGATGGCGGATTTATCTACCACTTTGTGGCATGCTTTTCTCGCCAGATAATGGAATGTTTCTCTGGCGGTAATGGATAGTATCAAATCTCACTAGCCCATTCTATAAGCGCATCCGCATGCACTAGTTCTTGATTCGATCGTCCTCTAGCATGTTCGAAGAAGATAGCATTCACTATCATCATCGCCTCAGGTAAGTTTATTCGGTTGGGGCGTGTGAAGGCAAACACCTTGTTGTCCAGTAAGTTTTCAGTTACTATAACTTAAAGTCGCACATGAATCTAGTCTCCTATTCCCCACCCATGATCCCACTCACACATTTCTACAGAGATGTGGTTAGAAATTTTCATATTAGGTCAGCTTTGACTCAATAAGACATAATTCTTCACTGAATGACTTTTTAAGAACCACCAGGACCAGAGAGAACCAAGAGAGTGGTACCTCTTAAAACACAATAAAGTGATTCAGCCTTAGCCATTGGATTCTGGAGGACCTTGAACCATGTGGGAAGCAGCTCAGGGTGGCCATGTACTATACTGGCGGGTAAGCTTCTGGAGTGCTAGGTTCTTTTTGTCTTTTCTTAAGCATTGCCGCCAGTTGATTGGGTTTTGAACATAAAATAATGCGCCACCAGCAATTCCAGATTTGTTCCTACGGGATAGATTTGTTCAGTTCTAGCATTATGCTTCACTAACCAGATGCGGGCCCTAAGTCCTTCACTTGGAATATTGGATTGGATCATGAGAATATTCTGTCTGAAGCTCGTCATTAATTTTGTTACAAAATAGAGCTTTTTGACTGGAAGTACCACCATACGTGTTCTCAAACTTCAGCATTTTTAGAACTTCCCACGGCATCTTGACCCTTTTCACAGCATGGATAGTCAGGCAGCAGTGAACTTTGTGACTCTTTAATGCCTTCACTTTTCTCTCAGTTTCCCCGCCTTGCGTTATCTTTACTCGTCTTGGGACTTTTATCCCAATGCCAGCCTTCTACCCTGAGACCTCAGTGGGTCATCATCCCAGCCCGGGACATCTCATCCCATCATTTATGGGCTGTTGTGTTTTTTTCAAAACCTAGCCCTCTCAGGAGGAGGAGGAGTGGGAGTCAGTTCAGTGAGGAGGATTAGGATGATCTGAAATGTAAGCACATATAAGCGAAGCACTTATTTTGGGTTGGGTCCTCACGGTGGACATAAGATCGCCTTATGTGTTTAGTAAGCCATTCCTAGCTCTCAATGGCGTGATTACATAGAAGCGTGAGGGATCAGTCCTATGGAAGACTAGGAAGTAAATGAACAAAATATATTAACCATAGAAGTCTCATGGGTCGCTGTAGCCAAAAGATTAACACTTTTGACTACATTGTGGTTTTAGGCATTGAAACAAAAACTTTGAGTCTCCTAAACAAATGAATGGAAAATAGTAGCGAACTTCGATTCCTAACATTAAATCTAGAAATAGCAAGTTAGTTTAAAGACTTTATTTAGCTTTGCTTGCTATAATGAAAACCTTGCCTCCCGGTCGGGGCCATTGTGCCTGAAGCTAGCTTATTGTCTCCTCGAGCTCCCAGCTTCAGCAACTCCTTTTGAAGCGTTTGTCTCAGCTTGGATCTTCAGCAGCTCTTGGTGGCTCTTTTGAGCTAGCTCCTCTGAGATCTTGTATTTGGTAGGTCGCTTAGTCATAGTACTTTTCTTTTAACACCCTTCAGCTCTACGATTACATTTGGTTTTGTGGATATCATAATGATTGATGTGAAGATACATTGTACATGTG
+
"#""#"##"'"&%&?JP7+80)'&*+&&(,.>3&*(#%(')*&'*#(-$)"&,63;?844&&#$'3+.>@9;-/259...&:"$"/*"#*#%(.&%)+,/76/+'4B93:;70.1)'1:)*#&()(+'03589/++&)'$-:;@,B5(9+.JAU;7-+&,1212.+329NPSDHRHH6),&)%)+).-*$,<1(.+&$$2&,-0+.&'//"4(%"+#,/)<56%*%*$<4*/;DI+6>((-%*).+--;2-'04-062<@-<>6CP)#.0/+,#(*2+#"$$-2'5(%%&'5+)$$%8+0;>.@9*(,((,*-2/49G9/1/3./.0CQ93E+/-D-'67J*)($-2+'*/&*)+%&'092<C9/5()$&%(.,+-'**)'2/.(//6<59-5-)"**40:A4+834)-(*+;K0-,&)+%)(#-&+&&))+$>9.6<?A;3'/85$/--..)/3DP99HJ9."2;I?@C/5%.4<0+0+&$-05C46208AFQE80"%356+02*9I73//(-*/9706192*/)%)(&'($1749C*KKA*/#"/*)5($'#+7@=$%)/&+%##*&&;.),#(4.,2/'-G03%""#,(2)./,%%1;7#+%$'(03'5)+&.%1)&(%'#*,$3#''-3,),$"),'+,,*$%'%$33ILADPPY@JW\O2-8'9[Y@--*.,""$%$"%#"65',0/4(6.*43B1.#&)+),(6BI729=9&#$.)#,(,(/-286D&),'C;P=)'-'#$)*91*$%(.3LH31,($/(+)'$$*#+'((('*)00:7C1E24(%+$&-,$'++)<)'**"&''"(*(.14++$5&*,%?.F1<10"&%$,)&&*(&$%$/.""+3"#"#"*%5*->87FQ>AKY<9855,&(*@>HAZRQ)'%$"--&/#%.+,/A9F8-1'+"6--*34P-9;=997;7<SUU8((-+.:784/%).&)#*'(*/%7H(1JI-+&+)&9.,/N)"%"&))22/&,)%)'+*/17:8$''(*.*..(/.(%+*&22(+)7GM99A@9)0*1%/470((&((%(937'15084#,($"&'%"'4.,+5<7($*0+*291$*&$%%('+/)'.1-(-.7=K/'%"'*)03<:/<=1(>.((-9-'+>&&&%';,)69V<'1:>FNXA+01:7*#"$'$.%&4BQT7/&03++)'&'*3E.+0(*6..*<.)9*)))'%#(*.&&'"*-IU=&,-))-.,->'G--6..90067%(,0/*%9)&#'""#$,5H/$/#-/&2(.0.+;>9--7.3@+2(/4.*0&,-&'/*%0-**-#,/2&-$*)'+#"*)-*##&$3):C:D6H+',./2?<0>)%#&&3::9;8$#%#('$#$',,,*-346@=.%&&(&%('*+&*&1)'&-..#'+2*)/6499=-2*,0,82,$)($:2%*##%$,.5=-5/((#)2./('($@=4,)%1+176:AHQT,(,)0(()@2'+(18G?($)(*'.)&&+*12SD7,18I2)()*+C:01ETZIPDUOB+.55K7B:A5)*/('')/-1=7E<IGR=43@A68;K/:8.'--)BGEE5-0/5+,2''DXBLW7)#DHTC5AL82+--,&#);+).,'3)/.336/&'&(,+8E346-/%1-)1*$+),@-&-%'$+)5.&&%&#(,2&03&/((,(-5%%)?7I4*').+/7+&***845(*/.)(.0))<?YK3=GHH2&*)%$%#%KLY<9&&+,8;/2*./;+))-)/03$&*$+0'36<J+5/+&+(HGFFM57,59)-,0)/:**"%$&#$'&$%#)*%33VX=@FOL00(''4IGIYOG0=1,,#(*$)-0CFM?8+31E8)(+)(%'"'/""%+3/,1;+7-*''&-%98.$$/'&)&+&(&&"

Should produce

93a12f52-95e5-40c7-8c3e-70bf94ed0720, 2017-07-04T06:42:43Z
ff37e422-a25f-404c-8314-ef1733f9c30c, 2017-07-04T06:56:41Z

awk '{if(NR%4==1) print $1, $5}' file.fastq | sed -e "s/ start_time=/, /" -e "s/^@//"

The awk command gets the first of every 4 lines, printing the first and fifth "word". sed is then used to strip the initial @ and replace start_time= with ,. The output on your example file is:
93a12f52-95e5-40c7-8c3e-70bf94ed0720, 2017-07-04T06:42:43Z
ff37e422-a25f-404c-8314-ef1733f9c30c, 2017-07-04T06:56:41Z

You could simplify your awk to awk 'NR%4==1{print $1, $5}', if you like.

Since the string start_time will only appear on the header line, or else you don't have a valid fastq file, you can simply do:
$ perl -ne '/^@(\S+).*start_time=(.*)/ && print "$1, $2\n"' file.fastq
93a12f52-95e5-40c7-8c3e-70bf94ed0720,2017-07-04T06:42:43Z
ff37e422-a25f-404c-8314-ef1733f9c30c,2017-07-04T06:56:41Z

Alternatively, since you mentioned awk and sed:
$ awk -v OFS=", " '/start_time/{print $1,$NF}' file.fastq | sed 's/start_time=//'
@93a12f52-95e5-40c7-8c3e-70bf94ed0720, 2017-07-04T06:42:43Z
@ff37e422-a25f-404c-8314-ef1733f9c30c, 2017-07-04T06:56:41Z

Or, doing the whole thing in awk:
$ awk 'sub(/start_time=/,""){print $1", "$NF}' file.fastq
@93a12f52-95e5-40c7-8c3e-70bf94ed0720, 2017-07-04T06:42:43Z
@ff37e422-a25f-404c-8314-ef1733f9c30c, 2017-07-04T06:56:41Z

And if the @ annoy you:
$ awk 'sub(/^@/,"") && sub(/ .*start_time=/,", ")' file.fastq
93a12f52-95e5-40c7-8c3e-70bf94ed0720, 2017-07-04T06:42:43Z
ff37e422-a25f-404c-8314-ef1733f9c30c, 2017-07-04T06:56:41Z

And in sed:
$ sed -n 's/^@\([^ ]*\).*start_time=\(.*\)/\1, \2/p' file.fastq
93a12f52-95e5-40c7-8c3e-70bf94ed0720, 2017-07-04T06:42:43Z
ff37e422-a25f-404c-8314-ef1733f9c30c, 2017-07-04T06:56:41Z

Or, if your sed supports it:
$ sed -En 's/^@(\S+).*start_time=(.*)/\1, \2/p' file.fastq
93a12f52-95e5-40c7-8c3e-70bf94ed0720, 2017-07-04T06:42:43Z
ff37e422-a25f-404c-8314-ef1733f9c30c, 2017-07-04T06:56:41Z

Finally, since grep can't do replacements, to use it you would have to do something like:
$ grep -oP '^@\K\S+|start_time=\K.*' file.fastq | paste - -
93a12f52-95e5-40c7-8c3e-70bf94ed0720    2017-07-04T06:42:43Z
ff37e422-a25f-404c-8314-ef1733f9c30c    2017-07-04T06:56:41Z

And to get the commas:
$ grep -oP '^@\K\S+|start_time=\K.*' file.fastq | paste - - | sed 's/\t/, /'
93a12f52-95e5-40c7-8c3e-70bf94ed0720, 2017-07-04T06:42:43Z
ff37e422-a25f-404c-8314-ef1733f9c30c, 2017-07-04T06:56:41Z

Albacore produces a sequencing_summary.txt file (actually TSV, not CSV) in the same directory as the workspace folder that might have the data that you want in it. Here are the fields present in that file:

filename
read_id
run_id
channel
start_time
duration
num_events
template_start
num_events_template
template_duration
num_called_template
sequence_length_template
mean_qscore_template
strand_score_template
complement_start
num_events_complement
complement_duration
num_called_complement
sequence_length_complement
mean_qscore_complement
strand_score_complement
sequence_length_2d
mean_qscore_2d

We'd thought the same. The issue is that the start_time that albacore calculates is actually the time since the start of a run. This is fine if you keep your minION switched on for a whole run. But if you reset it occasionally (which we do, because it really seems to increase yields) then the start_time column is a lot less useful.


text:
Building STAR Genome Index for nanopore RNA sequencing
I am aligning a dataset of 1,000,000 reads oh human mRNA sequenced on Oxford Nanopore Technologies' MinION, and would like to use the STAR aligner, using the parameters recommended by Pacific Biosciences for long reads.
According to this Google Groups thread, in setting up the genome index for short reads, the parameter sjdbOverhang should be set to 1 less than the read length. Obviously, with long (mean 1.7Kb, max >50Kb) reads, this doesn't make sense.
Running STARlong --runMode genomeGenerate without setting sjdbOverhang sets the parameter to 0. Does anyone know what this means, how it might affect mapping, and what I should set it to for my long reads?

Not an answer, but for full-length transcripts, it is worth considering traditional cDNA mappers such as gmap and blat. How well they perform highly depends on the error rate. They are probably ok at ~5% error rate. If you are using the older chemistry, they may not work, though.

A very good point. We're using GMAP as our baseline aligner, but I wanted to try out STAR to confirm the rather surprising results in this bioRxiv preprint: http://biorxiv.org/content/early/2017/04/11/126656

The parameter is used to determine how much sequence STAR indexes on each side of a splice junction to improve its alignment accuracy.  For very long reads, this may not be ideal. I am not sure if STAR is capable of including multiple splice junctions since a long read is more than likely to span more than one.
It may be worthwhile to consider aligning your reads directly to transcripts and taking the reads that align poorly and using them to look for novel transcripts.

I've found the STAR manual to be incredibly helpful with trying to navigate my way round all the STAR command line parameters. Here's the section on the sjdbOverhang parameter:

--sjdbOverhang specifies the length of the genomic sequence around the annotated junction to be used in constructing the splice junctions database. Ideally, this length should be equal to the ReadLength-1, where ReadLength is the length of the reads. For instance, for Illumina 2x100b paired-end reads, the ideal value is 100-1=99. In case of reads of varying length, the ideal value is max(ReadLength)-1. In most cases, the default value of 100 will work as well as the ideal value.

[emphasis preserved from the original text]
I would recommend trying out running STAR with this value unset, or set to 100, and seeing how it goes. If you want to be particularly adventurous, you could set it to 1000 and see if anything changes.


text:
Read length distribution from FASTA file
I have a single ~10GB FASTA file generated from an Oxford Nanopore Technologies' MinION run, with >1M reads of mean length ~8Kb. How can I quickly and efficiently calculate the distribution of read lengths?
A naive approach would be to read the FASTA file in Biopython, check the length of each sequence, store the lengths in a numpy array and plot the results using matplotlib, but this seems like reinventing the wheel.
Many solutions that work for short reads are inadequate for long reads. If I'm hey output a single (text) line per 1/10 bases, which would lead to a text output of upwards of 10,000 lines (and potentially more than 10x that) for a long read fasta.

If you want something quick and dirty you could rapidly index the FASTA with samtools faidx and then put the lengths column through R (other languages are available) on the command line.
samtools faidx $fasta
cut -f2 $fasta.fai | Rscript -e 'data <- as.numeric (readLines ("stdin")); summary(data); hist(data)'

This outputs a statistical summary, and creates a PDF in the current directory called Rplots.pdf, containing a histogram.

You can go even simpler with something like: faidx --transform chromsizes | cut -f2 | Rscript -e 'data <- as.numeric (readLines ("stdin")); summary(data); hist(data)'. This requires pyfaidx: pip install pyfaidx.

How can the Rscript be modified so the pdf file has a name based on the sample name?

Statistics for nanopore reads are tricky because of the huge range of read lengths that can be present in a single run. I have found that the best way to display lengths is by using a log scale on both the x axis (length) and the y axis (sequenced bases, or counts, depending on preference).
I have written my own scripts for doing this: one for generating the read lengths, and another for plotting the length distribution in various ways. The script that generates read lengths also spits out basic length summary statistics to standard error:
$ ~/scripts/fastx-length.pl > lengths_mtDNA_called.txt
Total sequences: 2110
Total length: 5.106649 Mb
Longest sequence: 107.414 kb
Shortest sequence: 219 b
Mean Length: 2.42 kb
Median Length: 1.504 kb
N50: 336 sequences; L50: 3.644 kb
N90: 1359 sequences; L90: 1.103 kb

$ ~/scripts/length_plot.r lengths_mtDNA_called.txt
lengths_mtDNA_called.txt ... done
Number of sequences: 2110
Length quantiles:
      0%      10%      20%      30%      40%      50%      60%      70%
   219.0    506.9    724.4    953.0   1196.2   1503.0   1859.2   2347.3
     80%      90%     100%
  3128.2   4804.7 107414.0

Here are a couple of of the produced graphs:

The scripts to generate these can be found here:

https://github.com/gringer/bioinfscripts/blob/master/fastx-length.pl
https://github.com/gringer/bioinfscripts/blob/master/length_plot.r

Using Biopython and matplotlib would seem like the way to go, indeed.
It really just boils down to three lines of code to get that graph:
import Bio, pandas
lengths = map(len, Bio.SeqIO.parse('/path/to/the/seqs.fasta', 'fasta'))
pandas.Series(lengths).hist(color='gray', bins=1000)

Of course you might want to make a longer script that's callable from the command line, with a couple options. You are welcome to use mine:
#!/usr/bin/env python2

"""
A custom made script to plot the distribution of lengths
in a fasta file.

Written by Lucas Sinclair.
Kopimi.

You can use this script from the shell like this:
$ ./fastq_length_hist --input seqs.fasta --out seqs.pdf
"""

###############################################################################
# Modules #
import argparse, sys, time, getpass, locale
from argparse import RawTextHelpFormatter
from Bio import SeqIO
import pandas

# Matplotlib #
import matplotlib
matplotlib.use('Agg', warn=False)
from matplotlib import pyplot

################################################################################
desc = "fasta_length_hist v1.0"
parser = argparse.ArgumentParser(description=desc, formatter_class=RawTextHelpFormatter)

# All the required arguments #
parser.add_argument("--input", help="The fasta file to process", type=str)
parser.add_argument("--out", type=str)

# All the optional arguments #
parser.add_argument("--x_log", default=True, type=bool)
parser.add_argument("--y_log", default=True, type=bool)

# Parse it #
args        = parser.parse_args()
input_path  = args.input
output_path = args.out
x_log       = bool(args.x_log)
y_log       = bool(args.y_log)

################################################################################
# Read #
lengths = map(len, SeqIO.parse(input_path, 'fasta'))

# Report #
sys.stderr.write("Read all lengths (%i sequences)\n" % len(lengths))
sys.stderr.write("Longest sequence: %i bp\n" % max(lengths))
sys.stderr.write("Shortest sequence: %i bp\n" % min(lengths))
sys.stderr.write("Making graph...\n")

# Data #
values = pandas.Series(lengths)

# Plot #
fig   = pyplot.figure()
axes  = values.hist(color='gray', bins=1000)
fig   = pyplot.gcf()
title = 'Distribution of sequence lengths'
axes.set_title(title)
axes.set_xlabel('Number of nucleotides in sequence')
axes.set_ylabel('Number of sequences with this length')
axes.xaxis.grid(False)

# Log #
if x_log: axes.set_yscale('symlog')
if y_log: axes.set_xscale('symlog')

# Adjust #
width=18.0; height=10.0; bottom=0.1; top=0.93; left=0.07; right=0.98
fig.set_figwidth(width)
fig.set_figheight(height)
fig.subplots_adjust(hspace=0.0, bottom=bottom, top=top, left=left, right=right)

# Data and source #
fig.text(0.99, 0.98, time.asctime(), horizontalalignment='right')
fig.text(0.01, 0.98, 'user: ' + getpass.getuser(), horizontalalignment='left')

# Nice digit grouping #
sep = ('x','y')
if 'x' in sep:
    locale.setlocale(locale.LC_ALL, '')
    seperate = lambda x,pos: locale.format("%d", x, grouping=True)
    axes.xaxis.set_major_formatter(matplotlib.ticker.FuncFormatter(seperate))
if 'y' in sep:
    locale.setlocale(locale.LC_ALL, '')
    seperate = lambda x,pos: locale.format("%d", x, grouping=True)
    axes.yaxis.set_major_formatter(matplotlib.ticker.FuncFormatter(seperate))

# Save it #
fig.savefig(output_path, format='pdf')

EDIT - an example output:

There are several potential approaches. For example:

histogram of sequence lengths in the Biopython tutorial
plot_distribution from the Ruby-based biopieces framework
various solutions to get sequence length including bioawk and EMBOSS infoseq

As to which of these are "quick and efficient" using a 10 GB file...it's hard to say in advance. You may have to try and benchmark a few of them.

You specifically asked about FASTA files, but it's important to always consider read length and quality jointly when assessing high-error long-read data. FASTA files do not provide the quality. This information will help you determine how successful the run was, how many reads were 'high quality', etc.
I originally posted a full answer here, suggesting pauvre, but I decided it was a little off topic since you seem to only have the FASTA files.
I recommend generating FASTQ files, but I'm not sure whether you have the original base-called fast5 files. If so, generate FASTQs using poretools as follows (poretools doc for generating FASTQ files):
poretools fastq fast5/

Then I recommend generating a heat map and histogram margin plot using pauvre with both read length and quality like shown below.

[Note: I am not the original author for pauvre, but I am now contributing to this project]

I would appreciate any feedback on why my answer may be disagreeable.

Dear Mark, I am having problems with pauvre. Can you help me?

bioawk could be reasonably efficient for this kind of task.
$ bioawk -c fastx '{histo[length($seq)]++} END {for (l in histo) print l,histo[l]}' \
    | sort -n
0   33270
1   1542
2   1132
3   3397
4   8776
5   11884
6   12474
7   14341
8   13165
9   15467
10  21089
11  30469
12  45204
13  62311
14  88744
15  115767
16  140770
17  191810
18  313088
19  518111
20  1097867
21  4729715
22  6575557
23  2734062
24  1015476
25  493323
26  323827
27  164419
28  107120
29  72487
30  40120
31  24538
32  22295
33  13121
34  9382
35  4847
36  3858
37  3161
38  2852
39  2388
40  1639
41  961
42  686
43  377
44  199
45  114
46  78
47  59
48  50
49  52
50  48
51  42
52  39
53  28
54  49
55  59
56  55
57  51
58  55
59  43
60  52
61  56
62  48
63  67
64  95
65  488

The -c fastx tells the program to parse the data as fastq or fasta. This gives access to the different parts of the records as $name, $seq (and $qual in the case of fastq format) in the awk code (bioawk is based on awk, so you can use whatever language features you want from awk).
Between the single quotes come a series of <condition> {<action>} blocks.
The first one has no <condition> part, which mean it is executed for each record. Here, it updates the lengths counts in a table which I named "histo". length is a predefined function in awk.
In the second block, the END condition means we want it to be executed after all the input has been processed. The action part consists in looping over the recorded length values and print them together with the associated count.
The output is piped to sort -n in order to sort the results numerically.
On my workstation, the above code took 20 seconds to execute for a 1.2G fasta file.

I realize that the output will not be convenient when dealing with sparse length values, because there is no binning (or, equivalently, the bins have width 1).

It is not exactly what you asked, but you can generate a histogram of read length distribution of your nanopore data directly from the HDF5 files using poretools

Due to disk space constraints, default behaviour for nanopore basecallers no longer creates the FAST5 (HDF5) files as output. While I actually have them, most users will not, and additionally this method would not generalize to other sequencers like PacBio.


text:
Improve a reference genome with sequencing data
I have a DNA sample which I know doesn't quite match my reference genome - my culture comes from a subpopulation which has undergone significant mutation since the reference was created.
The example I have in mind is E.coli. We've tried assembly using a couple of different tools and the de-novo assembly isn't as high quality as we would like, despite having tonnes of data. Approaching this from a Bayesian point of view, the reference genome provides a very good prior if we could use it wisely
From visual inspection with IGV, a significant number of both SNPs and SVs appear to be present, but an assembly built entirely from my own sequencing data is not high enough quality for my purposes.
How can I modify this reference genome to match my sample with new sequencing data (preferably with Oxford Nanopore Technologies long reads, but I can also use these to scaffold short reads if necessary), taking advantage of my knowledge that the existing reference is mostly very good, without having to access the reads which were originally used to construct the reference genome?
The goal of the project isn't to determine where the SVs are, I just need a reference that accurately represents my sample in order to use the data for downstream analysis (as the training set for machine learning.) So by a high quality reference, I mean one which represents as well as possible the sample that was sequenced. To make matters worse, this may not be the one which has the highest alignment identity if there are systematic sequencing errors, as in nanopore sequencing!

Depending on the coverage of your data and the complexity of the genome, you could either reassemble the genome de novo or run a reference guided (or reference assisted) assembly. It sounds like you're leaning more towards the latter.
There are a couple of reference-guided assembly tools available: AlignGraph and Ragout. These may or may not be appropriate depending on the organism of interest and your data types. For example, these tools are very unlikely to work well on Oxford Nanopore reads which have not been error corrected using Nanopolish or Canu -correct.

One approach to this is to use whatever data you have to iteratively update the reference genome. You can keep chain files along the way so you can convert coordinates (e.g. in gff files) from the original reference to your new pseudoreference.
A simple approach might be:

Align new data to existing reference
Call variants (e.g. samtools mpileup, GATK, or whatever is best for you)
Create new reference incorporating variants from 2
Rinse and repeat (i.e. go to 1)

You can track some simple stats as you do this - e.g. the number of new variants should decrease, the number of reads mapped should increase, and the mismatch rate should decrease, with every iteration of the above loop. Once the pseudoreference stabilises, you know you can't do much more.

If it is a short-read draft assembly and you have long-reads (ONT or Pacbio) run links to scaffold the genome and then run Pilon iteratively to try to polish and fill gaps using the short-reads.

You can use nanopolish using the illumina reads. Also have a look at pilon.

Kindel (which I wrote) can infer consensus from low quality alignments of short reads to viral references, and extending it to work with single molecule reads and larger genomes is on my to-do list, though I imagine this will require some redesign.
Presumably you're dealing with a bacterial or fungal genome in this case ? I do also have a basic C++ version, but it's a long way away from being user friendly. Anyway, it may be worth a look – feel free to get in touch with any issues you encounter. I'd use Pilon as mentioned above


text:
Samtools sort: most efficient approach to sort a single sample after aligning split .fastq files
Related to my other question (Samtools sort: most efficient memory and thread settings for many samples on a cluster), we need to optimize samtools sort as we prepare to analyze thousands of .bam files, beginning with re-alignment, sorting, etc.
Ultimate question: what is the fastest way to sort a single sample after aligning the sample across many separate cluster jobs? i.e., after splitting .fastqs for a single sample and aligning as distinct jobs. I'm looking for recent and clear time/thread/memory comparisons.
Unless I'm missing something, samtools merge was designed exactly for the purpose of merging pre-sorted bams. An excerpt from the manpage (http://www.htslib.org/doc/samtools-merge.html): "Merge multiple sorted alignment files, producing a single sorted output file that contains all the input records and maintains the existing sort order."
Is it fastest to:

samtools sort each 'mini-bam' and then samtools merge?
samtools cat the 'mini-bams' together and then samtools sort?
Another option I haven't considered?

Note: I have seen the following similar questions on other forums, but they are outdated (~8 years ago) and do not answer this question directly:

Best Way To Merge A Many Thousand Small Bam Files Into One Big Bam File?
Merging BAM files before sorting

As part of its operations, samtools sort will split the input up into multiple temporary BAM files (which are each individually sorted) if the reads can't fit entirely within memory. Those files are merged for the final output. From the manpage:

The  sorted output is written to standard output by default, or to the specified file (out.bam) when -o is used. This command will also create temporary files tmpprefix.%d.bam as needed when the entire alignment data cannot fit into  memory (as controlled via the -m option).

Given this, I would expect that a concatenated, unsorted BAM will take more time to process than merging sorted BAM files, because in addition to creating the concatenated BAM file, the splitting would need to be done again (assuming the reads don't fit in memory) prior to the final merge carried out by samtools sort.
I appreciate that this is not a quantified answer; it's just my guess based on what I have observed to happen when doing my own analyses.

It's funny. Now that you mention that, it seems so obvious! Great point. I may still quantify it just to see the difference.


text:
Nextflow: how to create many-to-many tuple for process
Using Nextflow, I need to submit x * y jobs where x is the number of input .bam files and y is the number of genome intervals (e.g., <chrom>:<start>-<end>). i.e., for every .bam file, submit a single job for every defined interval.
I have tried a few variations with no luck (see below).
Ultimate question: For each .bam file, how do I submit a single job for every defined genome interval using NextFlow?
Attempts
Attempt 1 (submit both as individual Channels):
Only submits one .bam for all intervals.
sample_bams = Channel.fromPath(sample_input_path + "*.bam")
intervals = Channel.from(['1:10000-20000', '5:55555-77777'])
PROC(sample_bams, intervals)

Attempt 2 (submit as tuple of mapped Channels):
Submits only 10 jobs and fails because the interval Channel is passed as DataflowBroadcast around DataflowStream[?] object. The PROC was changed to receive a tuple rather than individual arguments.
sample_interval_tuples = Channel.fromPath(sample_input_path + "*.bam")
                         .map { sample_file -> tuple(sample_file, align_to_ref, DRF_jar, Channel.from(['1:10000-20000', '5:55555-77777'])) }

Attempt 3 (submit as tuple of mapped .bams to non-Channel of intervals):
Submits 10 jobs and fails because the intervals are passed as a single list. The PROC was changed to receive a tuple rather than individual arguments.
sample_interval_tuples = Channel.fromPath(sample_input_path + "*.bam")
                         .map { sample_file -> tuple(sample_file, align_to_ref, DRF_jar, ['1:10000-20000', '5:55555-77777']) }

Really appreciate your help!
related query: for DSL v2, we wouldn't include from combined_inputs in the input, yeah?

The combine operator can be used to produce the Cartesian product:
sample_bams = Channel.fromPath( './path/to/bams/*.bam' )
intervals = Channel.of( '1:10000-20000', '5:55555-77777' )

sample_bams
    .combine(intervals)
    .set { combined_inputs }

process test {

    echo true

    input:
    tuple path(bam), val(interval) from combined_inputs

    """
    echo -n "${interval} ${bam}"
    """
}

Regarding the related query: Correct. The from and into 'bind' declarations can be omitted

we should crown Steve as King Of Nextflow Answers

I think this cheat sheet might help:
https://github.com/danrlu/Nextflow_cheatsheet/blob/main/nextflow_cheatsheet.pdf

Love the cheat sheet!


text:
Samtools sort: most efficient memory and thread settings for many samples on a cluster
We're preparing to analyze thousands of .bam files beginning with re-alignment, sorting, etc.
Complicated question: Has anyone investigated the optimal thread and memory settings to sort hundreds to thousands of .bam files using samtools sort on a cluster?
By 'optimal', I'm looking to sort as many samples as possible in the shortest time frame. I'm mostly interested in optimizing samtools sort because the other steps (e.g., alignment) are much more obvious to fine-tune.
For a single sample, throwing large amounts of memory + threads is trivial, but that isn't scalable for thousands of samples. There is also likely a diminishing return at a certain point, even for a single sample.
Ideally, I'd like to see plots identifying the 'sweet spot' for threads + memory for a single sample (e.g., x == threads, y == memory, z == time), and then plots maximizing throughput for more samples. I would expect the data to follow a -log function, where we could identify the optimal settings at the knee in the curve.
There are obviously many other factors (e.g., cluster size, administrative scheduler settings, etc.), but I suspect findings on a given cluster will generalize to other clusters.
Really appreciate any feedback.

Update 2
Ok, for what I expect will be my final update for this answer, I compared the following:

samtools 1.14 + zlib
samtools 1.15 + libdeflate (this is a different version, but shouldn't have a major effect. They released a new version just before I asked cluster admins to compile against libdeflate and didn't notice until later.)
sambamba v0.8.2

Note: This time I used the following settings:

Only two samples instead of three
Max Nextflow queue size of 30 to avoid too many threads reading from the same two files.
I tested CPU requests from 1 to 7 (step size 1) and then from 9 to 17 (step size 2).
For jobs where I allotted few CPUs and low memory, I provided a 200% buffer to prevent the jobs from failing with OUT OF MEMORY errors. i.e., I gave the SLURM job 2x the amount of memory than I told the tools they could use. sambamba was the worst offender for exceeding memory in these conditions, but I gave them all 200% buffer to keep comparisons equal.

Seems pretty clear to me that sambamba is the fastest and uses requested CPU resources most efficiently (i.e., is able to use all of the CPUs requested up until ~8 or 9 CPUs). For reference, the difference between sambamba and samtools + libdeflate with 9 CPUs is ~10 minutes on average.
We will certainly be using sambamba going forward, probably with ~9 CPUs and ~9GB per thread since I had to provide a 200% buffer in lower CPU/Mem jobs to avoid job failure, anyway.
Again, all code and results are available on my GitHub repository (Samtools sort optimization test).
(Angle 1) Comparison between samtools + zlib, samtools + libdeflate (Lsamtools), and sambamba
Here, I'm showing the general trend that real time decreases with increased CPUs for all three tools, but Mem Per CPU doesn't have a huge effect. sambamba out performs samtools + zlib regardless of settings. samtools + libdeflate out performs samtools + zlib until ~11 CPUs, where they perform the same (plot 3 below explains why).
Interesting finding: samtools + libdeflate performs better than sambamba with a single thread, but sambamba takes over from there as CPUs increase (again, plot #3 explains why).

(Angle 2) Comparison between samtools + zlib, samtools + libdeflate (Lsamtools), and sambamba
Here, I'm showing that samtools + zlib and sambamba performance flattens out at ~11 CPUs (again, plot 3 below explans why).

Requested CPUs vs. CPU utilization for samtools + zlib, samtools + libdeflate (Lsamtools), and sambamba
Here, I'm showing how well each tool is able to utilize the CPUs allotted to it. It was fascinating to see that samtools + libdeflate flattens out very early in it's CPU utilization (i.e., ~3 CPUS) and struggles to cross the 400% CPU utilization with even >10 CPUs. This is why samtools + libdeflate performance flattens out quickly.
samtools + zlib does a much better job utilizing allotted CPUs, but still lags behind sambamba throughout (and is slower, anyway).

Update 1

My first plots (below) were showing duration on the z-axis, which is the time from job submission to completion rather than job start to completion. That's why the larger memory requests were much longer in the original plots. Duration is still a useful measure, but I'm now including real time, which is from start to completion.
Judging solely by real time (above), there is a clear gain between 5 & 7 CPUs (~10 minutes, on average when compared to the same mem-per-thread).
More data coming, including samtools compiled against libdeflate and sambamba.

Original post
I made time to semi-formally look at this. I'm glad I did because the results surprised me. I don't answer all of my questions directly, but these results are good enough to get me going in the right direction.
My code and results are in my GitHub repository (Samtools_sort_optimization_test) where you will also find an interactive plot you can download to see the results more clearly (see SAMTOOLS_SORT/samtools_sort_CPU_and_memory_comparison.tar.gz in the repository). It's hard to truly appreciate the 3D plotly plot until you start rotating it.
NOTES:

Even though this is from only three samples, I think my results are accurate. I welcome any feedback in case there's a bug.
In my experience, samtools often exceeds its allotted memory (at least on a SLURM scheduler). i.e., the job gets killed by SLURM. To avoid this, I requested an additional 20% of memory to provide a buffer for samtools. I don't know if this discrepancy is because of SLURM or samtools, though.
These results are using samtools compiled against zlib. We are preparing a version of samtools compiled against libdeflate, per @user172818's suggestion. I will add those results once completed. I may also compare against sambamba, per @acvill's suggestion.
These results are related to my other question about the best sorting/merging method. I will post those results in that question and provide a brief summary here.

tl;dr (too long, didn't read)
Primary conclusions based on the plots below

Optimal CPUs (results as expected): There's a significant gain with some additional CPUs, but it levels off around 7 CPUs (total, not 'additional'; i.e., samtools sort -@ 6). I thought we might still see some nice gains above 7, but this is probably what most would expect (e.g., similar to what @user172818 suggested).
Optimal memory per thread (results NOT as expected): There is very little to gain with more than 1GB of memory (per thread); that surprised me. There are some performance gains up until ~7GB per thread (~2-4 minutes faster). What really surprised me, however, is that performance dramatically decreases as memory per thread goes >7GB. The updated plot (above) shows there is some gain with >1GB mem per thread. When running lots of samples at once (hundreds to 1000), mem per thread seems to play a much larger role (data coming).
Effect from sorting status (results as expected): Largely in an attempt to do this as quickly as possible (laziness?), I used three samples, where one was from an aligned and unsorted .bam file, and the other two were from .bam files that had been aligned and sorted, but I then collated them (samtools collate), rather than re-aligning them. Unsurprisingly, the collated samples (not the average use case) sorted slightly faster than the
unsorted sample (average use case).

Initial opinion

"Small" sample set: For a small sample set, I would probably use 7 total CPUs and 7GB per thread.
"Large" sample set: I will probably use 7 CPUs (maybe 5) with only 1GB memory per thread. That will allow the scheduler to keep as many jobs going as possible.

samtools sort: duration vs. CPUs and mem per thread
As mentioned above, an interactive version of this plot is available on my GitHub.
Angle 1

Angle 2

Additional details
Analyses were run on University of Kentucky's MCC Cluster, which uses AMD 7702P processors, using samtools v.1.14

I usually use samtools sort -@4 -m4g. Samtools will be faster with more threads and more RAM, but I feel the gain is not significant. I haven't done a proper evaluation, but I have run samtools sort with different settings from time to time on clusters when I have to consider speed-resource tradeoffs.
When you use samtools, make sure you compile it with libdeflate, which is about twice as fast as zlib on compression. I don't know if the bioconda version built with libdeflate. You may consider the precompiled binary here.
Samtools sort, especially with libdeflate, is fairly fast these days. Performance metrics in the other answer are outdated.

Super helpful info about libdeflate. I was not aware of that. I really want to do a comparison, but it's difficult to make the time. I may do it soon, and will try to include a comparison between zlib and libdeflate.

Side note: if you want to see if the samtools binary you're using is using zlib or libdeflate, you can see which libraries it was compiled against using ldd /path/to/<samtools_binary>. zlib will show as something similar to libz.so.1.


text:
efficient counting of dinucleotides/trinucleotides on fastq reads?
What's an efficient way of counting of dinucleotides/trinucleotide pattern on fastq.gz file reads?
I know there are tools like seqtk that will be very efficient at reading through the .fastq.gz files, and tools like cutadapt that are good at finding long(ish) adapter patterns in the fastq.gz PE of files, but, what would be the best tool to count the number of reads where a di-/trinucleotide is present N-number of times in the read?
E.g.
countingtool --pattern CGC file_R1.fastq.gz

Result: column two is number of times the pattern is seen (exclusive numbers, no overlap)
CGC    0     51341234
CGC    1     13948130
CGC    2      1234344
CGC    3        51344
CGC    4         1343

Are those numbers inclusive? I mean, are the 1343 reads that have 4 CGC included in the 13948130 that have 1? Are the reads with >N included in the count of reads with N?

should CGCGC as one occurrence or two?

number of times the pattern is seen (exclusive numbers, no overlap)

Are reads sequenced from a "stranded" library? If not, you need to consider the reverse complement. I would choose the strand, for each read, that gives the most count.

Well, I'm sure you can find a more efficient approach, but here's a simplistic solution using standard UNIX tools:
zcat file.fastq.gz | grep CGC |
    awk -F'CGC' 'NR%4==2{a[NF-1]++}END{for(i in a){print i,a[i]}}' |
         sort -rnk2

I tested this on a 2.6G fastq file (exome data) and it took:
real    1m18.987s
user    1m41.392s
sys     0m10.684s

On a 66G WGS sequence file, it took:
real    34m6.771s
user    40m16.452s
sys     4m10.176s

The extra grep CGC isn't needed, you can do the match in awk, but it does speed things up a little bit. You can always get the reads with 0 matches by counting the number of reads and subtracting the ones with 1. Alternatively, remove the grep:
zcat file.fq.gz | awk -F'CGC' 'NR%4==2{a[NF-1]++}END{for(i in a){print i,a[i]}}' |
    sort -rnk2

Following the OP's suggestion, I also tested with pigz and, it turns out, pigz can fly!
time unpigz -c exome.fastq.gz | grep CGC |
    awk -F'CGC' 'NR%4==2{a[NF-1]++}END{for(i in a){print i,a[i]}}' |
         sort -rnk2

real    0m44.012s
user    1m19.544s
sys     0m14.900s

And for the 66G WGS:
real    17m29.935s
user    29m23.816s
sys     5m17.856s

Brilliant, sometimes awk is all you need! Would pigz speed-up the process?

@719016 good idea! And yes, it certainly does speed it up, see updated answer.

Check out bioawk https://github.com/lh3/bioawk . This prevents counting CGC located in other lines. I wonder if it is faster than your current solutions?

@conchoecia I doubt it because I am taking advantage of awk's internal splitting by setting the field separator (FS) to CGC and using the % function to ignore irrelevant lines. Using bioawk with -c fastx adds a whole level of needless complexity since that will also look at the name, quality and comment lines as well. It also means I can't just use FS but would need to do something much splower like split(). And, indeed, I just tested on the exome file and stopped it after 7 minutes since it was taking too long already.