text:
How to extract / convert gff3 CDS sequences to multifasta
I would like to extract all the CDS from a batch of genomes. I have found a perl script from BioStars but this does not seem to work for me. I would preferably like to have a script/ method which will use the locus tag as a header.
Example
My gff files are from prokka output
ATCC0000    Prodigal:2.6    CDS 243 434 .   +   0   ID=FKKLIMLP_00001;Parent=FKKLIMLP_00001_gene;inference=ab initio prediction:Prodigal:2.6;locus_tag=FKKLIMLP_00001;product=hypothetical protein
ATCC0000    prokka  gene    243 434 .   +   .   ID=FKKLIMLP_00001_gene;locus_tag=FKKLIMLP_00001
ATCC0000    Prodigal:2.6    CDS 1727    2131    .   -   0   ID=FKKLIMLP_00002;Parent=FKKLIMLP_00002_gene;inference=ab initio prediction:Prodigal:2.6;locus_tag=FKKLIMLP_00002;product=hypothetical protein

If you are comfortable with python/biopython you can look into this GFF parser: https://github.com/chapmanb/bcbb/tree/master/gff

I think gffread is the best choice.

I like bedtools getfasta.  My typical option set is bedtools getfasta -fi <reference> -bed <gff_file> -name -s.  Be aware of the -s to make sure you are pulling the correct strand.  I like bedtools because it is a versatile tool overall for handling bed, gff and vcf file manipulations.
# bedtools getfasta

Tool:    bedtools getfasta (aka fastaFromBed)
Version: v2.26.0
Summary: Extract DNA sequences from a fasta file based on feature coordinates.

Usage:   bedtools getfasta [OPTIONS] -fi <fasta> -bed <bed/gff/vcf>

Options:
    -fi Input FASTA file
    -bed    BED/GFF/VCF file of ranges to extract from -fi
    -name   Use the name field for the FASTA header
    -split  given BED12 fmt., extract and concatenate the sequencesfrom the BED "blocks" (e.g., exons)
    -tab    Write output in TAB delimited format.
        - Default is FASTA format.

    -s  Force strandedness. If the feature occupies the antisense,
        strand, the sequence will be reverse complemented.
        - By default, strand information is ignored.

    -fullHeader Use full fasta header.
        - By default, only the word before the first space or tab is used.

The gffread utility in Cufflinks package might be interesting for you. To generate a multi-fasta file with nucleotide sequences from your GFF file, then you can try:
gffread -w output_transcripts.fasta -g reference_genome.fa input_transcripts.gff

Using python and this GFF parser that mimics Biopython's SeqIO parsers:
from BCBio import GFF

# Read the gff
for seq in GFF.parse('my_file.gff'):
    # only focus on the CDSs
    for feat in filter(lambda x: x.type == 'CDS',
                       seq.features):
        # extract the locus tag
        locus_tag = feat.qualifiers.get('locus_tag',
                                        ['unspecified'])[0]
        # extract the sequence
        dna_seq = seq = str(feat.extract(seq).seq)
        # simply print the sequence in fasta format
        print('>%s\n%s' % (locus_tag, dna_seq))

AFAIK GFF3 is a version of GFF, and the page that advertises the BCBio parser indicates that it is compatible in reading/writing GFF3: http://biopython.org/wiki/GFF_Parsing

The xtractore program from the AEGeAn Toolkit was designed for this type of use case. Just set --type=CDS.
$ xtractore -h

xtractore: extract sequences corresponding to annotated features from the
           given sequence file

Usage: xtractore [options] features.gff3 sequences.fasta
  Options:
    -d|--debug            print debugging output
    -h|--help             print this help message and exit
    -i|--idfile: FILE     file containing a list of feature IDs (1 per line
                          with no spaces); if provided, only features with
                          IDs in this file will be extracted
    -o|--outfile: FILE    file to which output sequences will be written;
                          default is terminal (stdout)
    -t|--type: STRING     feature type to extract; can be used multiple
                          times to extract features of multiple types
    -v|--version          print version number and exit
    -V|--verbose          print verbose warning and error messages
    -w|--width: INT       width of each line of sequence in the Fasta
                          output; default is 80; set to 0 for no
                          formatting

It looks like you're processing prokaryotic genomes, but this program works also on eukaryotic genomes where CDSs often have introns interspersed.

You can use agat_sp_extract_sequences.pl from AGAT.
See here for all option: https://agat.readthedocs.io/en/latest/tools/agat_sp_extract_sequences.html.
Below example of use:


text:
How to run enrichment analysis of protein functional annotation?
I have a lot of protein clusters. I want to perform an enrichment analysis of their functional annotations, against reference datasets or list of genes I select.
More precisely: a method yields cluster of proteins. I want to decide, for each cluster (which is a set of proteins identifiers), if it holds meaning regarding the proteins it contains, by comparing the annotations found in the cluster and the annotations found in common or specific datasets.
Initially, I used DAVID, which compute the GO annotations from the protein list, then perform the enrichment analysis against common datasets.
However, DAVID suffer of the following drawbacks:

since I didn't find any other way to use it, I'm limited by the web interface and its limitations (number of genes, url size, number of request/day).
automation over the web interface is a hacky script.

The team behind seems to offer a way to run DAVID in-house (allowing, hopefully, to relax the limitations), but I can't find any way to download the database from their website.
What are the alternative ways to get enrichment analysis over proteins in a reliable and automated way ?

This paper on the obsolescence of some enrichment analysis tools might be a good starting point to pick one (note that DAVID has been updated since): https://www.nature.com/nmeth/journal/v13/n9/full/nmeth.3963.html

What is the biological question? Why do you want to compare the functional annotations of your proteins against reference datasets or your list of genes? (I have usually found the opposite: I have a list of genes and I need to know what functions/pathways/process are they involved)

The question is do a cluster expose an enrichment in functional annotations. The method studied yields many clusters, and i want to decide if these clusters hold meaning regarding the proteins it contains, by comparing the annotations found in the cluster and the annotations found in common or specific datasets.

To make it even clearer, the clusters you mean you have a list of interacting proteins (in a complex or in a pathway) or that some tools grouped your proteins in a cluster? Could you [edit] the question to include the information in the body of your question?

A cluster is a set of proteins. There is no assumption about their interaction. I updated the question accordingly.

And the annotations you refer, what are they about? Sorry by these questions, but as you explain more is easier for me to think about an answer that fits you

The proteins -> annotations mapping was initially performed by DAVID. I think each database would have its own way to perform that mapping. I do not want to perform that step myself, but, as for the enrichment analysis itself, i will do it myself with some scripting only if its necessary. AFAIK, considered annotations was GO annotations. I updated the question. There is no problem ; i understand the value of your questions.

DAVID depend on a couple of databases from the Gene Ontology Consortium, Reactome, KEGG,... most of them are accessible via Bioconductor. To perform an enrichment analysis you can have a look at the tutorial of the several packages in Bioconductor that do this.
Some of the most important for analyzing enrichment in GO terms are topGO, goseq, GOstats. I would also recommend GOSemSim if you want to compare between GO  to focus on a specific GO terms.
Other important packages are the fgsea to test any kind of gene set (which is similar to the one hosted by the Broad Institute), gsva for enrichment analysis by sample, limma has some functions for functional enrichment too. Piano, GSCA, SPIA are also worth mentioning.
Bioconductor has "standard" data sets of the expressions of some cells, like airway and ALL frequently used in vignettes. They are not reference data sets because there isn't a reference expression for a cell of an organisms. It depends on the type of cell, the experiment, the conditions...

I would recommend R instead of sticking to websites.
There are many tools for enrichment analysis. I like to use GOseq, which is initially made for RNAseq data, but can also be used for proteins (I have used it for proteomics). You can use the length of the protein instead of transcript length (to correct for length bias), or you can exclude the length information by just using the hypergeometric method.
See manual for code.
Subsequently you can use gogadget, a tool I made for visualization of GOseq results with a heatmap or cytoscape network.
Additionally, you might wanna take a look at clusterProfiler. With this package you can do enrichment analysis for many clusters in parallel.

Under what circumstances is it appropriate to correct for bias using the length of the protein?

What I understand from the correction of bias in GOseq, is that if there is no bias a correction would not harm the analysis (because there is nothing to correct for). See previous discussion on Bioconductor. However, I do the hypergeometric test often, because the researchers (biologists) I work with prefer a test that they understand. The hypergeometric (Fisher exact) test is something many biologists understand.

@cmungall This would be a good question on its own. Please consider asking it

source: Stack Exchange

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

sort: votes


text:
Correct for gene length or read counts in GO enrichment analysis
It is a well reported fact that GO analysis of RNAseq results is affected by a number of biases, including length bias and expression level bias.
The bioconductor package goseq allows you to correct for these biases.
By default it corrects for length bias, but you can also get it to do read count bias. Using read counts to do the correction is attractive because in theory it should account for both sources of bias ($read counts\approx expression \times length$).
I'm doing an enrichment analysis were I have tried both options (length and read counts) and get very different answers. If I run a binomial regression on expression and length vs probability of being differential, I can see that both are independently important.

> model <- glm(sig ~  expression + log(length), data=retained_genes,  family=binomial(link="logit"))
> print(anova(model, test="Chisq"))

Analysis of Deviance Table

Model: binomial, link: logit

Response: sig

Terms added sequentially (first to last)

                       Df Deviance Resid. Df Resid. Dev  Pr(>Chi)
NULL                                    6676     4507.1
expression              1  114.998      6675     4392.1 < 2.2e-16 ***
log(length)             1  102.553      6674     4289.5 < 2.2e-16 ***
expression:log(length)  1   34.094      6673     4255.4 5.252e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

So I'm know unsure what to do, should I use the analysis corrected for length or read count. Or perhaps take only terms significant in both? Or only in one?

Did you test for the interaction of them (expression*log(length))?

I had not. Now I have. See updated question.

My understanding of the subject is that the bias of the gene length (and other bias) should be taken care when analyzing the expression and before enrichment analysis. The enrichment analysis should be done once the corrections has been performed. Because as the abstract of the GOseq paper states:

GO analysis is widely used ... but standard methods give biased results on RNA-seq data due to over-detection of differential expression for long and highly expressed transcripts.

So first, take care of the differential expression bias by length and then use the GO to reduce complexity. How you take care of the bias in the RNA-seq data is another question. But the cqn package of Bioconductor can correct the expression by gene length and GC content. However, this correction might hurt the differential tool used (see this discussion in Bioconductor), so it might be better in some cases to use GOSeq.
Now, to the question itself:

So I'm know unsure what to do, should I use the analysis corrected for length or read count. Or perhaps take only terms significant in both? Or only in one?

Use whatever correction method that yields better differentially expressed genes (DEG. If you find that the correction for length improves the accuracy of the predictions of DEG better than correcting by length and GC, then use that one.
Another option to obtain accurate GO terms, then you could use other testing procedures which don't rely uniquely in the Fisher test, such as those that take into account the structure of the GO graph. TopGO use this approach (note that it is a bit difficult to work with this package),this will reduce the role of the gene length bias (and probably other bias) in the resulting significant GO.

These bias' work through power. You cannot magically generate power out of thin air, so there are two ways you could remove the bias pre or at the DE stage: 1) reducing variance between samples 2) down-weighting more highly powered genes. As Mike Love notes in that link, modeling in the row-wise read counts is not compatible with DESeq2s methods (it were possible to adjust for extra power of highly expressed genes at the DE stage, presumably DESeq would already do it). Also who does one measure the accuracy of DEG prediction with no ground truth?

Taking into account the length of the genes won't reduce the variance between samples, so your only option is to down-weight highly powered genes. Maybe there are other methods to take it into account (or you could develop them). To find the true DEG for a given comparison: You can explore the FDR values of the genes in the comparison. The FDR values are assumed to follow an uniform distribution, because each gene has the same probability to be DEG. By exploring how uniform is that distribution you can find if the model captures the biases or not. You could also use the SeqGSA method

But I'm not interested in finding a better set of DEGs, rather, in better assessing whether their division in categories is what you would expect by chance. GOSeq allows you to do this while correcting for bias, but only one bias at once.

The problem in this approach is that the DEG you test, are affected by this length and GC bias. If you want to to assess their division in categories use topGO, because it is the only method which does an ORA analysis taking into account the structure of the GO graph.  GOSeq is biased by not taking into account the relationships between GO terms.

I am not sure if it makes sense to use read counts as bias instead of gene length (and I certainly wouldn't expect the same results).
Do you use total read counts of all your samples (library size)?
The correction for gene length is a pure technical one, the longer a gene the more reads will align (and higher read count genes are easier significant, since they are way above the noise threshold). If you use read counts, you also have a biological factor (expression) in there, which (I think) is the stuff you test with statistics (e.g., with edgeR) and thus not the bias you want to correct for.

Imagine you have two genes: Gene A has 10,000 counts in Condition A and 5,000 in condition B. Gene B has 1,000 counts in Condition A and 500 in condition B. The lfcs are the same, but Gene A is more likely to be call differential than gene B. Now consider two GO categories. Category X is made up of 10,000/5,000 genes and Category Y is made up of 1,000/500 genes. Cat X is more likely to be enriched than Cat Y, despite the real LFC in the genes being the same.

I don't know if the second example is true (that category X is more likely enriched). Maybe you should ask real statisticians on bioconductor?

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


text:
Identity plus Similarity Interpretation
I am helping a colleague to interpret protein multiple sequence alignment results. My colleague used a third-party proprietary solution (I know...) to calculate the alignment, and one of the resulting tables is a pairwise distance matrix with identity % in the top triangle and identity plus similarity in the bottom triangle.
I am struggling to interpret identity plus similarity, despite the fact I found the metric in several papers.

A guess; Identity: the amino acids that are equal. Similarity: The amino acids of the same group.

I guess you already know what identity means. But for completeness, it is the % of pairwise alignment position which have an identical amino acid residue. For the identity+similarity it is important to know that amino acids can be grouped according to their physicochemical properties e.g. charge, size etc. While there are standard groupings for amino acids it may be that the proprietary algorithm uses different groupings. Either way, % identity+similarity would imply the number of pairwise alignment positions which contain either the same amino acid or amino acids from the same group (however they define groups).


text:
Is there any difference between SNPs 'AG' and 'GA' in association analyses?
I am using plink to check the association between a phenotype and the SNPs of a gene. plink says the phenotype is significantly associated with a SNP on that gene, and when I check the SNP alleles, for 6 out of 49 samples, the nucleotides are 'GA', and for the rest, they are 'AG', so only the order is different. I am not a biologist, but my basic genetics knowledge says that those two are semantically the same. However, obviously, plink does not agree with me and it associates this SNP with a small p-value since 'GA' samples have a higher value for the phenotype than 'AG' samples. I am very new to SNP analysis, any help is appreciated.

Could you please show us the context in which this appears, as you seem to be interpreting this differently to Devon.
If it's appearing as you say, GA and AG, then Devon is right, this usually means that the sequence goes ###AG### or ###GA###, which are two very different sequences.
If, however, as you're implying, the sequences are actually ###G### or ###A###, then semantically, the two ways of writing it have the same outcome. Normally, however, this is written as A/G or G/A.
Alternatively, you're looking at genotypes, which would usually appear as A|G and G|A. In this case, both individuals are heterozygotes, and in terms of the individual SNP they are the same. However the side of the pipe is relevant as this tells you which homologous chromosome it is on. This allows you to calculate LD with respect to neighbouring SNPs, and to determine the actual sequences of the proteins in that individual.

Thanks for the reply. I am looking at genotypes. I realized that I had a bug while retrieving the genotypes, so what plink thinks as being different was not in fact GA vs. AG (each of which would really refer to heterozygous individuals as I replied to Devon), but it was GA and AA, where one is homozygous and the other is heterozygous. So, it is totally normal that those alleles are different for SNP association analysis.

"However the side of the pipe is relevant as this tells you which homologous chromosome it is on. This allows you to calculate LD with respect to neighbouring SNPs, and to determine the actual sequences of the proteins in that individual." -- Could you please explain this a little further?

If you have two variants, separated by, say 50, bp and they both fall in one gene. The genotype in one individual is G|A for the first variant and C|A for the second. That means that one of their genes contains a G and a C, the other gene contains an A and an A.
This might be relevant in various ways. In a population study, you might just be generally interested to know whether G...C and A...A are more common than G...A and A...C. This would also matter in a phenotype association study, as only one may be causative, but if they always turn up together then how would you know which?

(ran out of space!) It could also be that the phenotype only occurs when both variants occur, or is perhaps alleviated by the pairing. If both those variants affect amino acids, and those two amino acids are in contact, mutating one only could disrupt the connection, but mutating both the connection is maintained. Or if they're functionally redundant amino acids, one is able to compensate for the other, so when only one occurs, there's no phenotype, but when both do there is a phenotype.

There's a pair of courses here on introduction to genetic variation, which explains a lot about things like LD, effect on proteins etc: http://www.ebi.ac.uk/training/online/course/human-genetic-variation-i-introduction http://www.ebi.ac.uk/training/online/course/human-genetic-variation-ii-exploring-publicly-available-data

Is "user5054" the same as "user5504"? No? Exactly. Not only does order matter, it's incredibly vitally important. AG and GA are completely and totally different from each other. If this is in a coding region, then the resulting amino acid is undoubtedly changed (fun fact, the only exception is TAG and TGA). If it's at a splice site then splicing is likely altered (AG is part of a splice acceptor site). If it's at the binding site of something then it wouldn't be surprised if it didn't bind any more.

I think "user5054" vs. "user5504" example is misleading. My example would be that the sets {1,2} or {2,1} mean exactly the same thing, while if those represent permutations, they are different. So it depends on the case. In the SNP case, each of 'G' or 'A' in 'GA' or 'AG' are the nucleotides in the one of the two strands of DNA, right? So, both 'GA' and 'AG' are heterozygous alleles for the same SNP.Then what I understand is that the ordering of the strands matter,because maybe the transcribed strand's nucleotide is written first while representing the SNP?I think your reply needs more detail.

Is this actually a SNP or is it an InDel? If it's a SNP then you're correct that it doesn't matter (assuming no phasing information). My presumption is that this is not a SNP, but rather an InDel.

Yeah, this is a SNP analysis. I should have been clearer when I asked the question. Please see my comment to Emily_Ensembl above. The statement in my previous comment to your reply saying "...the ordering of the strands matter, because maybe the transcribed strand's nucleotide is written first while representing the SNP?" was completely out of confusion due to being new to genetic analyses. Apparently, Each of G or A in "GA" or "AG" do not refer to different strands of DNA, but refer to the different chromosomes of an individual.

Most association analyses are carried out at a single SNP level, so AG and GA are likely to indicate a heterozygous genotype at a particular location. However, the precise notation matters.
As @Emily_Ensembl has alluded to, for VCF files, A/G indicates unphased SNPs (order unknown, and shouldn't be considered in analyses), whereas A|G indicates phased SNPs (order matters).
It is possible to carry out association analyses by combining multiple variants into a single haplotype, in which case the order would also matter. It's unlikely that this is the case, because such haplotype variant analyses are substantially more complicated to carry out, and will typically be paired with information about SNP linkage (statistic D´ or r²).


text:
Cellranger aggr and Seurat merge difference?
I have multiple libraries of 10x Chromium single-cell RNA-seq data, which I'd like to combine. One option is using cellranger aggr which by default does a depth normalization:
mapped: (default) Subsample reads from higher-depth libraries until they all have an equal number of confidently mapped reads per cell. raw: Subsample reads from higher-depth libraries until they all have an equal number of total (i.e. raw, mapping-independent) reads per cell.
none: Do not normalize at all.
https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/aggregate
The new version of Seurat (2.0) provides another way via the MergeSeurat() (or AddSamples()) functions.
http://satijalab.org/seurat/merge_vignette.html
These have the arguments:
do.logNormalize whether to normalize the expression data per cell and transform to log space.
total.expr      scale factor in the log normalization
do.scale        In object@scale.data, perform row-scaling (gene-based z-score)
do.center       In object@scale.data, perform row-centering (gene-based centering)

I'm wondering how do these methods compare? What are the pros/cons etc? What do we need to consider when we combine samples?

Cell Ranger aggregate subsamples reads (unless you select none), so you will end up with less total reads in samples that have more initially. The output is still raw counts, but you will have more or less per cell.
Seurat just merges the raw counts matrices and normalizes those.


text:
Introduce errors in reference transcripts according to external dataset error model
I would like to modify some reference transcripts from Ensembl (D. melanogaster) to introduce a controlled rate of random errors in the sequences. The idea would be to introduce random base substitutions in these sequences, no indels for now, because I would like to keep the transcript sequence length as it is in the reference.
The rate of error per transcript will be determined according to an error profile computed from an external set of RNA-seq reads (e.g., generated with ONT MinION)
The aim of this modification would be to establish a rough benchmark of the performances of aligners to use over transcripts from spliced reads (rna-to-genome), aka with more than one exon.
Any idea of which software would be best for this purpose?

Do any of the answers for this question help? Karel Brinda has mentioned a few read simulators in the answer to that question, and has a thesis with more information.
Excluding INDEL errors doesn't sound like a good idea; length can still be preserved even if doing that, it just needs an adjustment at the end of the sequence. Note that if you're trying to model nanopore reads, what you're really modelling is the base-caller, rather than the sequencer. I mention this in more detail in my answer.
In most cases where errors are modelled, I find it better to use publicly-available data instead. Especially for nanopore data, there are unmodelled systematic errors in the base-callers and sequencer that can't be simulated using any programs (because they are unmodelled). The following paper would be a good place to start for cDNA sequences, which looks at single-cell data from mouse (C57Bl/6) B1a cells:
http://www.biorxiv.org/content/early/2017/04/13/126847
Illumina and ONT reads for that study can be found in SRA under accession number SRP082530.
I don't know of any recent D. melanogaster studies that have been done using nanopore. There's always the option of spending $1000 on a purchase of a MinION with an RNA starter kit to do the study yourself. Here's an older targeted gene study, but bear in mind that it was using an R7.3 flow cell, so errors rates will be much higher than what is currently available:
https://genomebiology.biomedcentral.com/articles/10.1186/s13059-015-0777-z

This preprint uses pbsim to simulate ONT RNA-seq reads for fruit fly. It is probably worth reading if you want to do the same thing.
You should include INDEL errors. Those are what make RNA-seq alignment challenging. For the benchmark purpose, adding INDELs does not increase the complexity at all. You can parse splice junctions on the reference from CIGAR and compare them to the annotation. You don't need to worry about the base-level alignment.
In addition, there are public real ONT data (AC:SRP082530) for SIRV spike-in control and mouse B cells. You don't actually need simulation.
PS: just noticed that you are an author of the first preprint I cited. I would use real data for evaluation.

It sounds like what you're really looking for is a read simulator. A cursory search turns up NanoSim, which is designed to simulate reads from a MinION. This has the benefit of at least having been used in some of the published literature, which is always a nice sign.
You may also find this review article on read simulators useful. It doesn't specifically mention NanoSim, but it should prove to be a useful review of the general concepts anyway if you need to read up on them.

Hum, not really what I am looking for because the simulation tool uses the model built in the previous step to produce in silico reads for a given reference genome - and the thing I need is to first, compute the error rate from the experimental reads (can be done parsing an alignment file in sam), and second, substitute in the reference as many nucleotides as needed to reach the average error rate of real ONT reads. for NanoSim I have the impression they are generating totally denovo 'reads' from the genome directly.

@AminaEchchiki If you don't want to use one of their pre-trained models, then you can have it train on real data and then it'll output reads with the appropriate error profile. Of course, since it uses LAST then whatever is most similar will have the best performance in your benchmarks.

The executable fastq-sim in DNemulator package is able to modify a set of input sequences in fasta format according to an external set of quality scores reported in a fastq file.


text:
Microbial diversity analysis using whole-genome metagenomic data
I have data, obtained from a single metagenomic DNA sample, that consists of two MiSeq FASTQ files (R1 and R2) that I merged using PEAR.
Now I want to estimate the abundances of the bacteria taxa to generate a figure like this one:

Figure from: Panosyan, Hovik, and Nils‐Kåre Birkeland. "Microbial diversity in an Armenian geothermal spring assessed by molecular and culture‐based methods." Journal of basic microbiology 54.11 (2014): 1240-1250.
The problem is that there wasn't a step of amplification of the 16S region as the goal of the sequencing was to discover new genes. I've already isolated 16S reads from my sample using SortMeRNA, but it seems like softwares that do OTU picking, taxonomic assignment and diversity analyses (such as mothur and QIIME) require that all the reads come from the same region of the 16S gene.
Is there a way of using these 16S reads that I've filtered using SortMeRNA in a diversity analysis using mothur/QIIME?

Cross-posted on biostars

New methods have been published in the time since this question was asked. Now, the relative abundance of organisms across different taxonomic levels can be estimated directly from metagenomic shotgun sequencing data without the need to identify and pull rRNA genes. Here are a few examples of popular tools:

Similar to the approach described in the question, mOTUs employs a database of single-copy phylogenetic marker genes to create taxonomic profiles from short-read metagenomic data.
Kraken associates short genomic substrings (k-mers) with the lowest common ancestor in its database. The output can then be processed with Bracken to get estimates of relative abundance at any taxonomic level, and the Bracken output table can then be passed to a tool like bracken_plot to visualize differences in relative abundance between samples.
Like mOTUs, MetaPhlAn uses a large marker gene database to assign reads to bacterial, archaeal, viral, and eukaryotic clades. This allows for relative abundance estimation in terms of number of cells. If their paper is to be believed, MetaPhlAn3 outperforms both mOTUs2 and Kraken2/Bracken in terms of abundance profile accuracy.

In addition to more accurate species detection, MetaPhlAn3 also quantified taxonomic abundance profiles more accurately compared to MetaPhlAn 2, mOTUs2, and Bracken based on Bray-Curtis similarity in most datasets ... While it was slightly outperformed by mOTUs2 (only in high-recall mode) on the synthetic mouse gut dataset, even in this case, correlation-based measures (Pearson Correlation Coefficient between estimated and expected relative abundances) found MetaPhlAn3 to be more accurate (r = 0.73) than the other considered profilers (MetaPhlAn2 r = 0.63, mOTUs2 precision r = 0.60, mOTUs2 recall r = 0.71, Bracken r = 0.43).

Although it is now dated, How to Choose Your Metagenomics Classification Tool from the Salzberg lab lists more options and links the associated publications, including Centrifuge and Pavian, the latter of which is geared specifically toward visualization of metagenomic classification data.

Its a thorough answer. Machine learning has enabled classification of bacteria, presumably on 16S. I'm aware of early papers in bioinformatics.

This is maybe a very simplified approach but have you thought about simply blasting them against the rRNA/ITS database of NCBI/ SILVA with ie. diamond?
First download the database files from either of the database:
SILVA: https://www.arb-silva.de/no_cache/download/archive/current/Exports/
NCBI: https://ftp.ncbi.nlm.nih.gov/blast/db/
Second: Use DIAMOND to align fastas of interest (in your case these will be the 16S rRNA Sequences. These will give you hits per fasta which you can save into tsv or csv format and parse in excel.
To get the abundances, take your fastas of interest and align them back to your raw files (paired reads). You can use BOWTIE + SAMTOOLs for this.
Concatenate the results from the DIAMOND output with the abundances and voila, you should get a nice table with taxa and their abundances. If you have more samples and you would like to compare them, ie. compare the taxa composition. I would strongly suggest that you do some further normalisation ie. normalise the abundances with sequencing depths.
Maybe it's a long way around it but in theory, it should do the job.

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


text:
Is there a way to retrieve several SAM fields faster than `samtools view | cut -f`?
I am constructing a bit of software which pipes the outputs of the bam file via samtools view into a script for parsing. My goal is to (somehow) make this process more efficient, and faster than samtools view.
I am only using 3-4 fields in the bam. Naturally, if I only wanted those fields, I could do something like samtools view file1.bam | cut -f ##.
However, I'm trying to be faster than this approach. More specifically:
(1) Given the bam file format, is it theoretically possible to output several requested SAM fields faster than samtools view|cut -f?
(2) Is there software available to achieve (1), or would one need to hack samtools?

Which fields do you need exactly? Using htslib is relatively easy if you know C, but if you can get by with something like sambamba then that'd be vastly easier still.

@DevonRyan Outputting the bam could be faster with sambambda? These are "mandatory" fields, so they're present in each bam "row"

samtools view, especially in the way it is used here, produces a text representation of a BAM file. Using a library (or tool that taps into a library) allows faster access at the binary level by avoiding full text conversion of each alignment prior to filtering.

@gringer I'm not following this exactly. Could you provide more details/an explanation?

Sambamba is generally faster than samtools, or at least that used to be the case. Samtools is a reference implementation, it doesn't always have to be the fastest.

The BAM file format is not a text-based format. It has a specific binary structure, specified in reasonable detail in the SAM file format specification. Whenever this information is displayed on a screen as text, it needs to be converted from the binary format to a text format, which takes a bit of time and processing power.
As this question suggests, if only a specific field from alignments is needed (or fields), then it will probably be better to extract just those fields and do any necessary conversion only on those fields. While this can be done by writing a BAM parser from scratch, many developers have already written software libraries to process BAM files in this way.
Devon Ryan has suggested htslib, which is a C library written by a group that includes the people who wrote the SAM/BAM/CRAM file format specifications. There's also pysam, which is a python wrapper around htslib.
The particular tool that is used will depend on your familiarity with programming and the specific thing that you want to do. If you just want to "output a BAM file", then cat is the quickest:
cat file1.bam

... but you probably don't want to do that, because it seems like you want to process a text-based representation with a script. Because you haven't specified which fields you are interested in, it's not possible to suggest the best thing to use. In the end, I expect that an "efficient" solution to your problem would involve htslib in some form. This is not really hacking samtools, it's using the backend of samtools to process BAM data.
However, it's useful to also think about the time cost of coding. What application do you have which means that the text-based processing method is not sufficient? It takes time to write code, and a lot more time to debug that code to make sure it's doing the right thing in all situations. If this is for a one-off thing, then samtools view output fed into your script may be the quickest solution.

I'm a bit familiar with pysam. It seems like this would be very inefficient, as I would need to input the bam into a python data structure, e.g.  pysam.AlignmentFile(). As for the script, yes, I currently am parsing the text-format of the outputted bam file. I have done nothing with compression/binary yet.

Whether or not it is inefficient depends on the way the library is implemented. It's possible for data structures to contain pointers to files (or places within files) and only retrieve information from the files when requested. There's no need to load the entire file into memory to explore just one aspect of it.

I modified your original question: as you are extracting 4 fields, you are not outputting BAM. The answer to the modified question is: yes, you can write a C program with htslib (or with bamtools, bioD, bioGo or rust-bio). Formatting an entire SAM is fairly expensive. You can see this by comparing samtools view aln.bam > /dev/null and samtools view -u aln.bam > /dev/null. With a C program, you can select fields to output. This will give you a noticeable performance boost, depending on the fields you extract.
In general, if you really care about performance, don't use script. The time you save from optimizing samtools view|cut -f will be quickly neutralized by the inefficiency of most scripting languages.

With a C program, you can select fields to output.

Can you detail a bit more how this would be done? The BAM format is in binary, so one would have to convert to text first and then extract the fields desired?

@ShanZhengYang The library parses the file record for you and places it in memory. From there, you print (or otherwise process) from the in-memory data only your desired fields in whatever format you want. There is no text parsing required.

Another library you can use for this purpose is the htsjdk, which is written in java. Using htsjdk with java is analogous to using htslib with C; the BAM format is already handled by the library and you can manipulate fields in your own code. The same basic analysis applies to java as C; you do not need to convert the BAM file to text and then parse it. Generally, a program implemented in java will be slower than C but still significantly faster than scripting. You do get language features of java: portability, array bounds checking, garbage collection, etc.


text:
Error rate setting in Canu error correction
I want to use Canu to correct my nanopore long read (version: MinION R9.5), but I am not quite sure how to set the correctErrorRate. Should I follow the Canu manual (Nanopore R7 2D and Nanopore R9 1D Increase the maximum allowed difference in overlaps from the default of 14.4% to 22.5% with correctedErrorRate=0.225), or you guys have a better option?

It's generally a good idea to trust the "official" suggestions. You can also adjust the error rate based on coverage according to parameter reference:

For low coverage datasets (less than 30X), we recommend increasing
  correctedErrorRate slightly, by 1% or so.
For high-coverage datasets (more than 60X), we recommend decreasing
  correctedErrorRate slighly, by 1% or so.
Raising the correctedErrorRate will increase run time. Likewise,
  decreasing correctedErrorRate will decrease run time, at the risk of
  missing overlaps and fracturing the assembly.

That being said, the sequencing quality can vary from between libraries or flow cells. The error rate is not a constant.
According to the developer (full thread):

if you want to assemble 1D data or have a bad run, it is possible you
  would need to increase the error rate. You could run with a high error
  rate (0.1) and look at the distribution of overlap error rates in the
  unitig step to look for a peak in the distribution. If you have a near
  neighbor you could also map the corrected reads to it and estimate the
  residual error that way.

There seems to be some confusion here.
For the first stage of the Canu pipeline (read correction), i think the parameter to use is rawErrorRate.
From the docs:
rawErrorRate  The allowed difference in an overlap between two uncorrected reads, expressed as fraction error. Sets corOvlErrorRate and corErrorRate. The rawErrorRate typically does not need to be modified. It might need to be increased if very early reads are being assembled. The default is 0.300 For PacBio reads, and 0.500 for Nanopore reads.
I think that correctedErrorRate is what you use at the next stages of the pipeline (trimming and assemble).


text:
Is it wise to use RepeatMasker on prokaryotes?
I'm looking for a way to identify low complexity regions and other repeats in the genome of Escherichia coli. I found that RepeatMasker may be used for example when drafting genomes of prokaryotes (E. coli example). But RepeatMasker works on a limited dataset of species, neither of them being prokaryotes. By default, when running RepeatMasker, if no species is specified, it will compare with homo sapiens data.
This seems rather inadequate, but the most relevent alternative, PRAP, requires a "dead" tool (VisCoSe, by Michael Spitzer).

Is it still wise to to use RepeatMasker on Escherichia coli?
If yes, which settings would maximise relevance ?

RepeatMasker isn't designed for use with prokaryotic genomes.  Nonetheless, it performs a contamination check for E. coli so you could play around with the -is_only flag to try and detect bacterial repeats, probably better to find an alternative tool or repeat library though

"Repeats" have different meanings. What types of repeats to mask are highly dependent on the downstream analyses – in fact, repeat masking is often discouraged. Explaining why you want to mask repeats would give you a more accurate answer.

@user172818 In fact, I am not interested in masking, but really finding any type of repeats, from low complexity regions to small repeats, as can be obtained using RepeatMasker.
I would then use these repeats as (sort of) an explanatory variable.

If I understood correctly your question, you want to mask those regions in a (FASTA?) genome. I think you could identify those regions using mummer and mask them using bedtools.
# align genome against itself
nucmer --maxmatch --nosimplify genome.fasta genome.fasta

# select repeats and convert the corrdinates to bed format
show-coords -r -T -H out.delta | awk '{if ($1 != $3 && $2 != $4) print $0}' | awk '{print $8"\t"$1"\t"$2}' > repeats.bed

# mask those bases with bedtools
bedtools maskfasta -fi genome.fasta -bed repeats.bed -fo masked.fasta

Have a look at nucmer and bedtools maskfasta options to fine-tune your analysis.

This approach may work, but it seems quite an un-orthodox way to mask repeats, have you used this or seen people use it for bacterial genomes?

I have used it, following the advice found on the mummer manual: http://mummer.sourceforge.net/manual/#identifyingrepeats

thanks, you used it for prokaryotes? they says in the docs it's not really designed for that and only identifies a limited number of repeat types, so I assumed this was not really advised?

yes, I did use it on E. coli K-12, finding 1324 repeats. I never tried other methods so I'm unsure how it compares with those. It would definitely be interesting. I wouldn't say the authors of mummer discourage you to use nucmer for finding repeats, just that it wasn't specifically made for that purpose.

From your comment, it seems that masking regions are not your priority, but you would rather find them (correct me if I am wrong):

not interested in masking, but really finding any type of repeats, from low complexity regions to small repeats

To find these regions, you can try RepeatFinder. From their paper, it seems that is suitable for bacterial genomes too. It also seems to be faster because based on suffix tree data structure instead of working with string-matching algorithm (as in RepeatMasker).


text:
Is it possible to create a DESeqDataSet with a user-provided design matrix?
I'm trying to run a differential gene expression analysis using DESeq2, with counts coming from kallisto. I have imported them using tximport and I'm creating the DESeqDataSet (dds) using the DESeqDataSetFromMatrix function.
> dds <- DESeqDataSetFromMatrix(counts,
                                samples,
                                ~ strain + batch)

And I get the following error, expected given my experimental design:
Error in checkFullRank(modelMatrix) :
  the model matrix is not full rank, so the model cannot be fit as specified.
  One or more variables or interaction terms in the design formula are linear
  combinations of the others and must be removed.

Now, I know that I can just remove one column from the design matrix to make it work, but is there a way to supply my own design matrix to DESeq2? The following code raises an error:
> design <- model.matrix(~strain+batch, s2c)
> design = design[, -9] # removing a column to avoid full-rank
> dds <- DESeqDataSetFromMatrix(counts, s2c, design=design)
converting counts to integer mode
Error: $ operator is invalid for atomic vectors

Is there a way to provide my own model.matrix?
p.s. the modified model works in sleuth, but I would like to use DESeq2 for this particular analysis.

Provide rank sufficient design to DESeqDataSetFromMatrix and then use your custom model matrix in DESeq. In essence:
dds = DESeqDataSetFromMatric(counts, s2c, design=~batch)
design <- model.matrix(~strain+batch, s2c)
design = design[, -9]
DESeq(dds, full=design)

See this thread on the bioconductor site for details.

If your strain and batch are confounded, (like you did all of strain one on one day, and all of strain two on another), then there is no clever way to fix that.  Just use one column, and inform whomever designed the experiment that you can't separate batch effect from strain difference.


text:
What is indel calling and what is its purpose?
I'm having a difficulty in grasping the general purpose and concept of indel calling.
What exactly is this process?

Other than finding InDels, which tend to be highly disruptive?

Insertions and deletions (indels) are one type among many different types of genetic variation, such as single nucleotide variants (SNVs), copy number variants (CNVs), and structural variants (SVs). I'll assume here that you know how indels are defined, but are simple trying to understand the importance of discovering and analyzing them.
The goal of indel calling, like the goal of any variant calling, is to identify genetic variants that can subsequently be associated with important phenotypes, esp. disease. For example, if 60% of patients with disease XYZ have an indel in the promoter region of gene 123, then that is information of extreme interest and value in research and in clinical care.
Genome-wide association studies (GWAS) have been trying to correlate SNVs to disease and other phenotypes for years. Much less work has been done with indels, but their discovery and analysis remains an area of intense interest.
As far as the process of indel calling, large indels can usually be found  by mapping paired reads and looking for large discrepancies in the expected distance between pairs and the observed distance.

Huh, my average insert size is 400bp, but the aligned read pairs flanking this area are 1200bp apart. Must be an 800bp deletion in there!

Smaller indels are much more difficult to detect using this strategy, since they are harder to distinguish from noise (i.e. the variation in length of sequenced fragments). However, as another answer mentions, many short indels are simply reflected as short gaps in the alignment.

There are two types of INDELs: short indels and long indels. Some put the threshold at 50bp; others choose 1000bp. Short and long indels are called differently.
<50bp short indels are called from read alignment directly. Modern indel callers essentially build a multi-alignment of reads and call an indel with enough read supports. Short indels may break open reading frames.
For Illumina reads, long indels are usually called with split alignment or read pair alignment because mainstream aligners are unable to directly align through a long indel. Long indels are a type of structural variation. They may affect gene structures in a more dramatic way (e.g. delete a whole exon; screw up transcription by transposon insertion; create pseudogenes).
Ultimately, both short and long indels are types of genetic variants. Calling them helps to understand how genetics shapes phenotypes.


text:
RIP-seq analysis?
Given an experiment consisting of an input (baseline RNA) and IP (pulldown to find RNAs bound to certain protein of interest)... Is a DE analysis performed over the RNA-seq data from the samples (lets say with EdgeR or DESEQ2) suitable to reveal the preferentially bound RNAs? What other software tools would you recommend?

You need to be careful of terminology. To me, a RIP-seq experiment involves a pull down, followed by a RNAseq library prep. Thus, the whole transcript is captured, not just the binding site of the protein (as in CLIP-Seq, HITS-CLIP, PAR-CLIP, iCLIP or eCLIP). Thus "peak-callers" whether they be designed for calling protein-DNA or protein-RNA binding sites are not suitable as the signal will not be punctuate or peaky. Of the methods mentioned by @Devon Ryan, only RIPSeeker seems setup to deal with this sort of data.
In some of the papers that talk about calling RNA-protein interactions they lump RIP-seq in with the CLIP techniques (see piranha and ASPeak papers referenced by @Devon Ryan). They appear to be talking about a technique where the RNA is fragmented before it is pulled down. Thus you would capture peaky binding sites (actually, the protocol would be remarkably similar to CLIP).
You should work out which of these you have.
As for using DESeq/EdgeR etc, in the RIPSeeker paper they say:

Furthermore, programs for de novo transcript assembly followed by
  differential expression (DE) analysis, such as the Cufflinks/Cuffdiff
  suite (15,16), and for DE on a set of known transcripts, such as DESeq
  (17), may appear applicable to RIP-seq analysis. Unlike peak-calling
  strategy, however, the transcript-based methods assume the full
  transcriptome being sequenced at a fairly deep coverage (as usually
  the case in RNA-seq) and thus may be sensitive to background noise
  typical to the IP-based protocols, which is due to both the
  non-specific RNA interactions with a protein of interest and the
  non-specific RNA input from the pull-down of the (mutant) control
  (Supplementary Figures S3 and S4).

I don't know to what extent RIPSeeker performs better than the naive DESeq approach as its not included as one of the comparators in the paper.

How important is peakiness though? Macs2 seems to deal well with broad peaks. Couldn’t it be used here?

There are number of reasons not to use a protein-DNA peak caller, many detailed in the RIPSeeker paper. Many of these are to do with splicing: fragment sizes estimated from mapped RIPseq will not reflect genuine fragment sizes, using either extended reads or fragment midpoints will likely leave you in the middle of an intron.
MACS2 broadpeak works by merging narrow peaks connected by sub-significant sequence. Its not clear this is a good description of the data you'd get in RIP. The RIPSkeeker paper compares to Macs, and looks like it out-performs by quite a way to me.

You can try doing standard differential expression, but I worry that the between-sample normalization will work poorly. Personally, I would do peak calling instead, followed by diffBind. You have a few tools to choose from when it comes to this. In the past, I've rolled my own methods for this using MACS2 and genomic alignments (I then converted those to bedGraph files where entries are transcripts, so peak calling is then in transcript-coordinates).
It's probably easier, though, to use premade software. There are a number of packages out there for RIP-seq. These include piranha, RIPSeeker and ASPeak. Particularly if you're just starting out, you'd be well served with using one of those tools.

I've got the RIP-seq analysis outputs. We used an isotype control antibody VS antibody against an RBP to pull down corresponding RNAs!
I did the following pipeline: Hisat2 (mapping the reads after QC) -> htseq-count (getting read counts) -> DESeq2 (getting enriched RNAs in test VS isotype control) -> volcano plot shows up-regulated and down-regulated genes in the test sample.
Logically, I should have seen just enriched RNAs (up-regulated genes), but surprisingly there are some RNAs down-regulated in the test sample (or up-regulated in the isotype control sample)!!!!!
How can we explain and interpret this???? it seems some genes are pull down in isotype control sample as well. Are these background? Is my computational pipeline right?
BTW, it's worth to mention that all enriched RNA (called up-regulated genes in DESeq2 pipeline) are significantly relevant and meaningful to my experiment. However, I am not sure how much confident I can be on this output? And how should I explain down-regulated RNAs? because normally we just expect enriched RNAs in RIP-seq experiment.
Any thoughts?


text:
Ensembl id to GeneSymbol with biomart
I have 3224 Ensembl id's as rownames in a dataframe "G". To convert Ensembl ids into Genesymbols I used biomart like following.
library('biomaRt')
mart <- useDataset("hsapiens_gene_ensembl", useMart("ensembl"))
genes <- rownames(G)
G <-G[,-6]
G_list <- getBM(filters= "ensembl_gene_id", attributes= c("ensembl_gene_id","hgnc_symbol"),values=genes,mart= mart)

Now in G_list I can see only 3200 ensembl ids showing Genesymbols / No Gene_symbols. Why the other 24 ensembl ids are not seen in G_list? If there are no gene_symbol for those 24 ensembl ids it should atleast show "-"
Examples of problematic IDs are: ENSG00000257061, ENSG00000255778, ENSG00000267268. These are not at all shown in G_list (biomaRt). So, I gave them in biodbnet, which seems to handle them.
what is the problem here?

I used gProfileR where there are no problematic IDs with that. It converted ENSEMBL IDs to Gene_symbols and made my work easier with GO analysis too. There is an R package and also API.

could you post a code snippet? or is it trivial?

Here is the reference manual for gProfileR [https://cran.r-project.org/web/packages/gProfileR/gProfileR.pdf]

It looks like you were using an old annotation. The problematic IDs you posted existed in the GRCh37 annotations, but don't in the most recent GRCh38 annotation. For that reason they were excluded. The IDs that have - as symbols don't have associated symbols, but are present in the database.
To use an archived version in biomart:
mart = useDataset("hsapiens_gene_ensembl", useEnsembl(biomart="ensembl", version=84))

That's an example for release 84.

biomaRt will return NA for genes that are in the database, but don't have symbols.

@IanSudbery Yeah, I'm assuming that OP replaced NA with - at some point.

As you said, I thought the same that the problematic ID's are existed in GRCh37. And tried with that too. And among 3224, I found only 3089 in the results.

Those IDs existed in the GRCh38 release 78 annotation (just checked a local GTF), so perhaps they're still GRCh38, but just an older release of it.

I actually did it in this way. install.packages("devtools")
devtools::install_github("stephenturner/annotables")
library(dplyr)
library(annotables)
grch38

So ENSG00000257061 was last seen in ENSEMBL 84 (you can look this up by searching on the ENSEMBL website). You need to search against the exact same version as you got your gene list from. Replace useMart with useEnsembl(version=X) where X is the version of ensembl you used to generate the gene list.

So...can anyone of you please tell how can I use the latest version in the above code? or Is there any R package for biodbnet?

@user3351523 I've updated my answer. The biodbnet question should be a separate post.

Just by-the-by, I'm going to guess that any ID that has been retired is unlikely to map to a HGNC symbol.

@DevonRyan With the above mentioned version 84 I got 3220 now. Still 4 ids missing.

Try version 83 then.

@user3351523, you have to use the same ensembl version as was used for generating your genes object, or you'll never get 100% conversion.

I used "gProfileR". There are no problematic IDs with that. It converted ENSEMBL IDs to Gene_symbols and made my work easier with GO analysis too. There is an R package and also API. Thank you all for the help !!


text:
Greedy Motif Search Using Probability Matrices
In my greedy search with pseudocounts algorithm in my bioinformatics course, I did not follow the pseudocode since I wanted to solve the problem in my own way. Unfortunately, although my answer is correct, the answer we are supposed to give is the first correct answer, not the last.
The problem works as such: Given t strings, you make a profile matrix from the first k-mer in the first string, and then you look at the next string and using that profile matrix, you find the best k-mer match, and then you combine that k-mer with the one from the first string to make a new profile matrix, etc. You eventually do this for all k-mer windows in your first string of DNA. The way I did it was by making a probability profile matrix from the best k-mer groups in each case, and then multiplying the greatest probabilities in each column. This gives me a score, and if the next score is higher (thus, more probable) that group of k-mers is used. Below is my code:
 def highest_probability(string, k, matrix):
    score = 1
    temp_string = ""
    c = 0
    best_string = ""
    temp = 0
    while c+k < len(string):
        for i in range(c, c+k):
            if string[i] == 'A':
                score *= matrix[i-c]['A']
                temp_string += 'A'
            elif string[i] == 'C':
                score *= matrix[i-c]['C']
                temp_string += 'C'
            elif string[i] == 'G':
                score *= matrix[i-c]['G']
                temp_string += 'G'
            elif string[i] == 'T':
                score *= matrix[i-c]['T']
                temp_string += 'T'
            if i == c+k-1:
                if score >= temp:
                    temp = score
                    best_string = temp_string
        temp_string=""
        score=1
        c+=1
    print(best_string)
    return best_string

def GreedyMotifSearch(DNA, k, t):
    best_motifs = [i[0:k] for i in DNA]
    score=0
    c=0
    while c+k<len(DNA[0]):
        k_mer = DNA[0][c:c+k]
        motif=[k_mer]
        matrix = make_matrix(k_mer, matrix=[], i=0)
        for i in range(1, t):
            compare_string = DNA[i]
            good_motif = highest_probability(compare_string, k, matrix)
            motif.append(good_motif)
            matrix = make_matrix(good_motif, matrix, i)
        temp_score=1
        for dictionary in matrix:
            temp_score *= max(dictionary.values())
        print(temp_score)
        if temp_score > score:
            score = temp_score
            best_motifs = motif
        c+=1
    return best_motifs

def make_matrix(string, matrix, i):
    if len(matrix)==0:
        for ch in string:
            new_dict = {'A':0, 'C':0, 'G':0, 'T':0}
            new_dict[ch] = 1
            matrix += [new_dict]
        print(matrix)
        return matrix
    elif len(matrix)==3:
        for j in range(len(string)):
            matrix[j][string[j]] += 1
        Sum = [sum(matrix[0].values()),sum(matrix[1].values()),sum(matrix[2].values())]
        for columns in range(len(matrix)):
            for keys in matrix[columns]:
                matrix[columns][keys] /= Sum[columns]
        print(matrix)
        return matrix

DNA=['GGCGTTCAGGCA',
'AAGAATCAGTCA',
'CAAGGAGTTCGC',
'CACGTCAATCAC',
'CAATAATATTCG']
k=3
t=5
print(GreedyMotifSearch(DNA, k, t))

The output gives me ['TTC', 'GTC', 'TTC', 'GTC', 'TTC'] which is equally correct to the answer given in the problem output which is TTC,ATC,TTC,ATC,TTC. Unfortunately, my answer is the last correct occurrence. I tried adjusting my highest_probability function with if score is greater than temp instead of greater or equal, but that gives me the answer ['GGC', '', '', '', ''] which I don't understand at all. My question is 1) Why is this the case? and 2) How can I fix this? And does anyone have a simpler way of doing this?

Assuming you're talking about Implement GreedyMotifSearch from Rosalind, while I haven't actually finished this particular one, I've got an idea why TTC,ATC,TTC,ATC,TTC doesn't work. Here's my explanation from the "Questions" section of that problem:

The profile is created from the first kmers in each string, and the most likely string from those is CAG. The answers will generally be similar to this.

I've been trying to do Rosalind by only using the information provided in the questions; this is one of the many problems where there wasn't enough information in the Question for me to arrive at a correct solution when I tried it. It's possible that the debug datasets were added to this question since I last looked at it, as they seem to describe this situation:
TEST DATASET 1:
Input:
3 4
GCCCAA
GGCCTG
AACCTA
TTCCTT

Output:
GCC
GCC
AAC
TTC

This dataset checks that your code always picks the first-occurring Profile-most Probable k-mer in a given sequence of Dna. In the first sequence (“GCCCAA”), “GCC” and “CCA” are both Profile-most Probable k-mers. However, you must return “GCC” since it occurs earlier than “CCA”. Thus, if the first sequence of your output is “CCA”, this test case fails your code.


text:
Is there any R package for biodbnet?
I want to convert Ensembl ids to Gene_symbols. I used biomaRt package but all Ensembl ids are not seen in the results.
Examples of problematic IDs are: ENSG00000255414, ENSG00000274358, ENSG00000235565, ENSG00000278309. I gave them in bioDBnet, which seems to handle them.
I would like to know whether there is any R package for bioDBnet or not?
If so, could you please tell how to use it?

Welcome to Bioinformatics! This question seems to need more information to get a good answer. Could you please [edit] adding what is the G_list, and how did you create it? What do you mean by "I gave them in biodbnet"? Also some ensemble genes might not have an official name from the HGNC

Sorry. I edited the question now. Please check the link I gave !!

You might consider ignoring or removing those genes. The reason you dont find these genes are that they have been removed from the current data bases meaning the Ensembl people thinks it's a false postive (see http://www.ensembl.org/Homo_sapiens/Gene/Idhistory?g=ENSG00000255414 for example)

As @KristofferVitting-Seerup says, the real issue here is that you should not expect every ENSG to map to a HGNC symbol, for a variety of reasons.

Ok. "gProfileR" is the best one for converting into gene_symbols and GO analysis too. There are no problematic IDs with that. Anyways thankyou all for the help.

I couldn't find a package neither in Github nor in Bioconductor. Neither could I find any paper about a package linking to bioDBnet. So I don't believe there is any (R) package at the moment which links to biodbnet in order to retrieve the IDS.
However it has a couple of APIs so you could build one or use them to retrieve the information you want.

I used "gProfileR". There are no problematic IDs with that. It converted ENSEMBL IDs to Gene_symbols and made my work easier with GO analysis too. There is an R package and also API. I recommend everyone to use this. Thank you all for the help !!

@user3351523 You are commenting to something different to what you asked. In this question you want to access to biodbnet, while in this comment you want to convert some IDs from ENSEMBL to Gene Symbols (Which can be done with several packages).

library(EnsDb.Hsapiens.v86)
library(ensembldb)
genes(EnsDb.Hsapiens.v86)
genes(EnsDb.Hsapiens.v86, filter=GeneIdFilter('ENSG00000235857'))
transcripts(EnsDb.Hsapiens.v86, filter = GeneIdFilter("ENSG00000184895"),
              order.by = "tx_seq_start")

(The ids you listed are older than v86). See the vignette for the ensembldb Bioconductor package.

Although helpful, I think it won't help the original poster (and other readers) if you do not explain how did you identify that the ids provided were older than a specific version. Also explaining that this package is able to convert these ids would help to the original poster


text:
Where are .motif files from homer knownResults?
I have been using homer's findMotifsGenome.pl, but with my new version (v4.9.1) of homer I don't get .motif files in the knownResults folder. I do get them in the homerResults (de novo) folder, though. With my previous version of homer I did get them and used them to find which peaks are in which motif. Which is what I would like to do
Does anyone know why there are no .motif files? I have only .svg files. I have seen that for the homerResults the .motif files are concatenated already into a file called homerMotifs.all.motifs, but no such file for knownResults.
PS. I tried to call Chuck for help, but he was not answering his phone...
PS2. I will ask this question on Biostars to hopefully get more response than here.

You seem to have (relative easy) access to the code, have you checked it?

What do you mean Llopis?

who is Chuck :')

Norris, it's a running gag thru out the homer user guide ;P

@b.nota findMotifsGenome.pl is a perl script, so you could read the code and see what happens.

Sorry Llopis, that is not helpful. Thanks anyway.

I have provided a solution in Biostars:
If you edit the file findKnownMotifs.pl and move the call to printMotif to come before testing $svgFlag you will find that .motif files are printed again.
In the following I have commented out the call from where it need not be, and introduced it as the first line, right before the if block:
   printMotif($_, $directory . "/$knownDir/$matrixFile"); # this should be done regardless of $svgFlag
    if ($svgFlag) {
        my $str = HomerSVGLogo::motif2svg($_,$bitFlag,$rnaFlag2);
            open SVG, ">$directory/$knownDir/$logoFile.svg";
            print SVG $str;
            close SVG;
            print MAIN "$str\n";
    } else {
  #printMotif($_, $directory . "/$knownDir/$matrixFile"); # this should be done regardless of $svgFlag
            `profile2seq.pl "$directory/$knownDir/$matrixFile" 100 > "$tmpFile"`;
            my $width = $widthFactor * $_->{'len'};
            `seqlogo -a $bits -f "$tmpFile" -F PNG -S -c -o "$directory/$knownDir/$logoFile" -h $height -w $width`;
  `seqlogo -a $bits -f "$tmpFile" -F PDF -S -c -o "$directory/$knownDir/$logoFile" -h $height -w $width 2> /dev/null`;
  print MAIN "<IMG src=\"$knownDir/$logoFile.png\"/>";
}

I expect this change was mistakenly introduced while refactoring the code to support svg as alternate output format for logos.
I have reported this to the author of this code.

Please consider adding the relevant parts of your answer here too, in case Biostars isn't accessible and to prevent rotten links.

Great job, malcook! Upvote from me!

The relevant help page of homer seems to be this one.
On section 13 it says:

An HTML page is created in the output directory named homerResults.html along with a directory named "homerResults/" that contains all of the image and other support files to create the page.  These pages are explicitly created by running a subprogram called "compareMotifs.pl".

You seeem to have this images on the expected folder.
Also on section 12 it says:

They [The *.motifs files] are reported in the output directories from findMotifs.pl and findMotifsGenome.pl

So you'll need to specify the output directory for findMotifsGenome.pl.

Thanks again, I appreciate it. Section 13 is about homerResults, which is the de novo results folder. This seems fine, they contain .motif files. However, my question is about the knownResults folder, which is for known motifs, not de novo. As mentioned in my question they do not contain .motifs files, but according to the help page you sent me they should contain them.


text:
What format is this? Pretty sure it's not a BED file
Time for guess the bioinformatics file format... :)
I found this '.bed' file on GEO (GSE84660), but it's clearly not a BED file. Anyone know what it might be? And what can view it? It's from a HiC expt.
#column headers: bait1_chr, bait1_start, bait1_end, bait2_chr, bait2_start, bait2_end, contact_ID, NA, NA, B1_D0_counts, B1_D3_counts, B1_D6_counts, B2_D0_counts, B2_D3_counts, B2_D6_counts

chr1    831895  848168  chr1    850619  874081  b2b_1   .       .       114     125     127     75      97      92
chr1    831895  848168  chr1    889424  903640  b2b_2   .       .       15      12      16      4       15      20
chr1    831895  848168  chr1    903641  927394  b2b_3   .       .       13      13      10      18      15      8
...

Could you explain what you need it for? I mean, that's just a simple text file so any text editor (or even excel) can view it. What information do you need to extract from it?

Are you trying to visualize the contact matrix? This data comes from this paper.

@DevonRyan ideally yes, but first I wanted to check that this wasn't some specific format I wasn't aware of. Have any suggestions to visualise a contact matrix?

@ithinkiam I'd suggest HiCExplorer, but I'm biased. To your original question, you'll still need to munge this file, but I'll see if I can come up with some code for that.

@ithinkiam I take that back, it's not 100% clear what the appropriate way is to go from large bait regions to single positions. Hopefully someone else has a better idea.

Perhaps it's been called .bed to suggest that the start/end coordinates are to be interpreted in a BED-style 0-based half-open (or equivalently 0-based start, 1-based end) way?  The pair of genomic intervals per record is reminiscent of BEDPE format, but I'm not familiar with anything HiC so wouldn't comment on whether this is mere coincidence…

It looks like a weird way to represent data. You're right, it does not at all look like a BED.
Digging in GEO made me find this information:

Supplementary_files_format_and_content: Tab-delimited file reports one
  contact per row. Interacting HindIII fragments are represented in bed
  format (chr/start/stop) with bait (or upstream bait, in the case of
  b2b contacts) listed first. The following fields correspond to a
  unique contact ID, fragment IDs for first and second fragments, and
  raw read counts supporting contacts in each replicates.

So I think the authors merged a BED3 format plus additional BED-unrelated fields generating a tab-separated file - not clearly stated in their methods though. The easiest would be to contact the authors maybe?
Depending on your aims, if you only need the first 3 fields (but I doubt it) I would reformat this "bed" file (to make it processable as a BED file) in two sub-files:
cat original.bed | cut -f 1,2,3 > bait1.txt
cat bait1.txt | head -3
chr1    831895  848168
chr1    831895  848168
chr1    831895  848168

cat original.bed | cut -f 4,5,6 > bait2.txt
cat bait2.txt | head -3
chr1    850619  874081
chr1    889424  903640
chr1    903641  927394

Not sure if it would fit your aims though.
Also, for visualization of contact matrices, maybe you'll find this paper useful.

You could simplify your command to cut -f 1-3 original.bed > bait1.txt if you like. The cat isn't needed, cut can take a file path as input and lists of fields can be defined with x-y.

This is one of the umpteen ways one can think of representing a contact matrix and is likley particular to the supplemental data in this paper (since it's unusual to store 6 contact matrices in a single file). If OP want to visualize this, then it'll likely need to be munge it into a very different format.

Thanks for the help, but it's not really what I'm after.

ok, so maybe would be wise to edit your question and make your aims clear, so users can give more meaningful answers. or you just wanted to verify is not some kind of documented format you might have skipped?


text:
Problems to Extract uniquely mapping reads from BWA MEM alignment
I did a mapping of genomic paired-end reads to a reference assembly using bwa mem. I need to extract the reads that mapped only once to my reference.
For that, I have tried to follow the method proposed in another post Obtaining uniquely mapped reads from BWA mem alignment and unfortunately in my case, after using:
$ samtools view -h output.bam | grep -v -e 'XA:Z:' -e 'SA:Z:' | samtools view -b > output_filtered.bam
$ samtools view output_filtered.bam | awk '{print $1}' | sort | uniq -d

I still obtain some outputs. Here is an example:
NB501461:7:HH3GGAFXX:4:21612:26866:20400    99  ctg7180000256968    441 0   147M4S  =   560 230 AAGCTTTAGTCGCATGGCCATTCACAAGACCAGTGATTTTTTATAGAAGTCGTTCAAGAGATACTGATGAGGAGTACAGCCATTCTTTAAGCTATACAAAAAGCACCACCCACACCACACAATCCTCAAAAAGATCAAGGCAGTACAGATT AAAAAEEAEEAEAAEEEAAEEEAEEEEE/EEEEEEAA/EEEEEEEEEEEEEEAEE/EEEE///E6EAEEEEAEA<EE/EA/</EEAE/EEEEE//EEEEEE6E/<EAE/EEE/EE</E/E//AE<6EA//A<A///<A/6A/6//A///// NM:i:7  MD:Z:34A56T11A4A3A2A22T8    AS:i:112    XS:i:117
NB501461:7:HH3GGAFXX:4:21612:26866:20400    147 ctg7180000256968    560 0   32S111M8S   =   441 -230    CCAGCTCCACCACCCGACCCCCCCACCCACCCCAATCCTCAAACAGATCCATGCAGTACATTTTCATTTGTTCTGGCCCTGCGGCTCTTAGTTTAGGAAATATTCAATGACTTGCTACCTAGTAAGAACATTATTGTGCACAACGAATAGA ///E/////A////A/6A/EE/<//////A//AA/EEE/A/<E///A/A<//</////A/////<///EE////E/E/E/E</E//EA//E//A/</E//E//AE<//A<///EE///////E//E//EA/A//EA/AE/A/E/E//6/A/ NM:i:7  MD:Z:11A5A22A4A4T27A7A24    AS:i:76 XS:i:79
NB501461:7:HH3GGAFXX:4:21612:6399:20400 113 ctg7180000282074    357 8   73S50M28S   ctg7180000432624    28  0   CTTATCCAATTTTTAAACATAGGATTAGGGTACCTTCCAAATGTTATTGAAACTTAATCTGAAGATGACCTCAAATGGAAGATGACCTCCAATGGAAGACAACCTTCCATGGAAGACGACCTCAAAAATGAAGAAGAAGAAAACATCTAAG /A//<//<<//E/66A//<A////A/<</E/EA/EE//<E<E//<AE//A/EEAEEEE/E//A/A/</EEE//E<EAEEEEE/EEEEEE<EEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE6EEEEEEEEEEEEEEAAAAA NM:i:2  MD:Z:32C1A15    AS:i:40 XS:i:35
NB501461:7:HH3GGAFXX:4:21612:6399:20400 177 ctg7180000432624    28  26  76S75M  ctg7180000282074    357 0   CTCCTTACTGTTGTTTCTTCTTTTTTCTTCTCTTTCATTTTTGAGGTTGTCTTCCATGTTTTGTTGTCTTCCATTTGAGGTCGTCTTCATATTAAGTCTCTATAACATTTGGCTTTGTCCCCTAAGCCTAAGTTTCAAATTTGGAAAAGAT /////</////6/////A////EAEE<//////<A////E/E<//A//<///E//<E/////E/A/E//E/EAEEAE/EA///E//E/AE/<EAEE/6E6EEA/EE/AEEAE6EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAAAAA NM:i:4  MD:Z:5T15T17G29G5   AS:i:55 XS:i:44

As I do not see any XA or SA flags there anymore, I wonder what those can be?
I have tried to continue the steps I want to pursue until the "ValidateSamFile" from Picard tools, and there I do not know if it is linked or not, but I get the message:
"Error: mate not found" for a few thousands of reads.
Any idea what could be wrong, and how to fix it?
I hope I was clear enough, otherwise please feel free to ask me for specific details.
Any help is welcome,
Thanks!

The issue is that you are filtering reads without regard to the pairing so you might filter out one of the mates if one aligns uniquely and the other does not.  You will need to filter out the un-mated reads from the file to have a valid BAM.

As I do not see any XA or SA flags there anymore, I wonder what those can be?

I am not sure what this means.  The grep should remove those flags and only leave records without them.
As a side note if you don't want split alignments I would recommend using a different aligner such as BWA aln or BowTie2

my goal with using grep -v  was to remove all the lines corresponding to the reads flagged with XA and SA, which mean all the reads that are mapped more than once. I am surprised that I still get duplicates of read names afterwards because there should be none left.

You have paired reads so one of those is the read 1
and the other is the read2.  The 99 maps to read 1, among other things and the 147 maps to read 2.  See https://broadinstitute.github.io/picard/explain-flags.html

extract the reads that mapped only once to my reference

Every read can be mapped to multiple places if you allow enough mismatches. You heard "unique mapping" from time to time. Strictly speaking, however, this concept is either ill defined or practically useless given 150bp reads. For a proper aligner, the best way to measure the reliability of an alignment is simply mapping quality. All the other approaches, including the answer in the other thread, have problems.
By the way, many tools allow you to set a threshold on mapping quality. You can apply the threshold on the fly. It is a common and good practice to avoid writing filtered alignments to another BAM.


text:
How to convert DNAbin to FASTA in R?
I am trying to convert DNAbin files to fasta format. Rationale: I want to use the fasta file to calculate non-synonymous/synonymous mutation rate.
my_dnabin1 is a DNAbin file of 55 samples and I am using the following code to convert it into a fasta file.
dnabin_to_fasta <- lapply(my_dnabin1, function(x) as.character(x[1:length(x)]))
This generates a list of 55 samples which looks like:
$SS.11.01
[1] "t" "t" "a" "c" "c" "t" "a" "a" "a" "a" "a" "g" "c" "c" "g" "c" "t" "t" "c" "c" "c" "t" "c" "c" "a" "a"
 [27] "c" "c" "c" "t" "a" "g" "a" "a" "g" "c" "a" "a" "a" "c" "c" "t" "t" "t" "c" "a" "a" "c" "c" "c" "c" "a"
$SS.11.02
[1] "t" "t" "a" "c" "c" "t" "a" "a" "a" "a" "a" "g" "c" "c" "g" "c" "t" "t" "c" "c" "c" "t" "c" "c" "a" "a"
 [27] "c" "c" "c" "t" "a" "g" "a" "a" "g" "c" "a" "a" "a" "c" "c" "t" "t" "t" "c" "a" "a" "c" "c" "c" "c" "a"
and so on...
However, I want a fasta formatted file as the output that may look something like:
>SS.11.01 ttacctga
>SS.11.02 ttacctga

The following should work:
cat(file="out.fa", paste(paste0(">",names(dnabin_to_fasta)),
          sapply(dnabin_to_fasta, paste, collapse=""), sep="\n"), sep="\n");

A bit of explanation, starting from the middle outwards:

Collapse the letter arrays inside dnabin_to_fasta into single strings
Append ">" to the front of the names of each sequence
Concatenate the results of steps 1 & 2 together, separating with a line break
Write each concatenated string to a file, separating by a line break

Thank you so much! It solved a major part of my problem. When I print it to the screen, I get the exact output I want. However, it does not produce any output files. Could you please help me with this, too?

Note that R will save files in the "working directory", which is not necessarily the directory that the last file was loaded from. If you're using cat with the file option, it's probably saving somewhere. You can find out where by using the getwd() command.

Thank you so much! It solved my problem! I was stuck in it since a long time.


text:
Extracting expression data from GSE dataset downloaded from GEO
I have downloaded GSE16146 dataset from GEO using GEOquery R package.
I would like to extract "Data table" from downloaded GSE16146.
>library("GEOquery")
>GSE16146 <- getGEO("GSE16146")
>Table(GSE16146)[1:5,]

This returns the following error:
> Table(GSE16146)[1:5,]
Error in (function (classes, fdef, mtable)  :
unable to find an inherited method for function ‘Table’ for signature ‘"list"’

Can anyone help?

Have you tried the function getGSEDataTables to access GSE Data Tables from GEO, as reported in this tutorial? Also, what version of R/GEOquery are you using?

GEOquery version 2.42.0, R version 3.4.0  > getGSEDataTables(gset)
Error in as.character.default(<S4 object of class "GSE">) :
no method for coercing this S4 class to a vector

The error is telling you that you can't calculate a Table for a list. Also you might have tried to use table instead of Table.

According to the manual, all you need to do is:
library('GEOquery')
gseGSE16146 <- getGEO('GSE16146', GSEMatrix=FALSE)

As explanation, getGEO() outputs by default to GSEMatrix=TRUE and returns a list of ExpressionSet objects. You should get what you were looking for with:
Table(GSMList(gseGSE16146)[[1]])[1:5,]

The manual has also a paragraph about this:

The GSE has  a  metadata section, just like the other classes.
  However, it doesn’t have a GEODataTable.  Instead, it contains two
  lists, accessible using GPLList and GSMList, that are each lists of
  GPL and GSM objects.

I tested the code I posted on R 3.3.2 and I get access to the Data Table just fine.
update:
Answering how to extract the actual expression data (not what asked in the original question):
library(GEOquery)
data = getGEO("GSE16146")
datExpr = exprs (data[[1]])

Thanks. It works. Would you suggest how can I extract all the expression matrix for 80 samples (list) from the gseGSE16146?

Hi thanks for your answer! What I didn't get is how do I access the actual expression data?

@Tapper I updated my answer above with how to extract the expression data from the experiment. hope it helps?

@aechchiki Thank you! The snippet does not contain the data but a list in my case though: GSM4446615 GSM4446616 GSM4446617 GSM4446618. The file getGEO() got was anyway to small to contain the dataset imho. I finally got the data by downloading the big data file myself and processing it with: gse = read.table("GSE123456.txt", sep = "\t"); data=CreateSeuratObject(raw.data = gse, ..) -- What do you think?


text:
Database of copy number alterations in different cancer types
Is there a database of somatic copy number alterations (CNAs) in different cancer types? That is, annotations of which CNAs are more commonly seen in different cancer types.

I don't know if there's a database that does exactly what you want, but there are some places that might help you figure this out, especially if you already have a list of CNAs/genes/regions in mind.
ICGC is the first that comes to mind, as it has samples from many types of cancer and copy number data for most of them. Depending on what you want though, you might have to apply for access to download the data you want. It's rather annoying, and their data download tools are atrocious to use, but they have a lot of data.
DGV might be another option, though their data browser is rather dated feeling. They publish a lot of gold-standard germline CNVs and such.
COSMIC is a large database for cancer-associated mutations, including CNVs, so it might be worth a look for you.


text:
Convert local alignments to spliced alignments in SAM file
I mapped RNA reads to reference genome, using LAST in split mode, and converted the MAF alignment to SAM with maf-convert.
My problem is that the transcripts are not reported in a spliced manner, meaning that a transcript_ID is reported several times in the alignment SAM file with an identical bitwise flag in $2. From what I understand, this is due to the fact that only the exons are mapped (one exon per row) and reported as local alignments, as the software cannot deal with splicing patterns combination exon-intron (yet), which behaviour is clear from the CIGAR string.
For a more concrete and visual example let's consider the mapping of transcript FBtr0344900 to the reference genome as given by LAST:
$ cat last_aln.sam | grep FBtr0344900
FBtr0344900 0   4   42774   100 384=144H    *   0   0   TGCGACATTGTTCTACGATGACTACAAAAAATGACCAATAACTTCTATAAACCAATACGATATGTCAGGAGTTTCGGTCCCATACGAAGTCGCCGACTTAAGTATTTTATttttattttgatATGTGTTTGCTATTTTACCTTGTCGAATGCTTCCACACGCTATGAGAATACCATCGTGAGCGTAGCTTACTACTAGAATTTTGTTGAAGTTATTGACAAGCGATGTCTCAATATCTTCCGGACAGCCTCCAGCGTGACATTGCGGGGAATCATGTAACGGCCCAGTAACAGCCTCGGCCAGCACTCGAAGGTTTTCGTTAAGTTTAAGTATTTTATTTGTAGCACCCGCAAACAAAACATTGTGCATAAAGTCGAAGCTCAT    *   NM:i:0  AS:i:2304
FBtr0344900 0   4   43231   100 384H144=    *   0   0   CTGGAAGCTGTTGATTGAACTGGTATTGATGGCAAGTTAAACTGGGCGACTATGTCATTTAAGGGAGATAACGCCTGAGCCGGCAGTTCTTCAATGCAGTTAACGCAATAATGCTGAGAACCGAGTATGATAATAATACACAGT    *   NM:i:0  AS:i:864

And here is the mapping of the same transcript FBtr0344900 as given by STAR - software which reports the alignment right the way I need it:
cat star_aln.sam | grep FBtr0344900
FBtr0344900    0    4    42774    255    384M73N144M    *    0    0TGCGACATTGTTCTACGATGACTACAAAAAATGACCAATAACTTCTATAAACCAATACGATATGTCAGGAGTTTCGGTCCCATACGAAGTCGCCGACTTAAGTATTTTATTTTTATTTTGATATGTGTTTGCTATTTTACCTTGTCGAATGCTTCCACACGCTATGAGAATACCATCGTGAGCGTAGCTTACTACTAGAATTTTGTTGAAGTTATTGACAAGCGATGTCTCAATATCTTCCGGACAGCCTCCAGCGTGACATTGCGGGGAATCATGTAACGGCCCAGTAACAGCCTCGGCCAGCACTCGAAGGTTTTCGTTAAGTTTAAGTATTTTATTTGTAGCACCCGCAAACAAAACATTGTGCATAAAGTCGAAGCTCATCTGGAAGCTGTTGATTGAACTGGTATTGATGGCAAGTTAAACTGGGCGACTATGTCATTTAAGGGAGATAACGCCTGAGCCGGCAGTTCTTCAATGCAGTTAACGCAATAATGCTGAGAACCGAGTATGATAATAATACACAGT    *    NH:i:1    HI:i:1    NM:i:0    MD:Z:528

From discussion with the author, it seems that I cannot get what I need straight from the current LAST version, and is not a technical problem. So, I have to modify the output myself. The aim would be to get at least a CIGAR line representing a complete transcript.
My question is, do you know any software to do this? What I would need is a file containing one line per uniquely mapped transcripts, featuring three fields: transcript_ID, start position and CIGAR string.
From my side, I proceeded like this :
1) extract the fields that I need for the  interested in from the SAM file:
$ cut -f1,4,6
FBtr0344900  42774  384=144H
FBtr0344900  43231  384H144=

2) split the CIGAR line to remove items I am not interested in - I simplify the command here, just assuming I only have perfect matches (which I am interested in) and hard-clipping (which I am not interested in):
$ cut -f3 | sed 's/H/_ /g'| sed 's/=/= /g' | sed 's/\w*_\s*//' | sed 's/=//g'
384
144

3) paste the modified CIGAR in the original file, with paste, resulting in:
$ paste (1) (2) | cut -f1,2,4
FBtr0344900  42774  384
FBtr0344900  43231  144

4) merge rows which start ba the same transcript_id:
$ awk -F'\t' -v OFS='\t' '{ x=$1 ; $1="" ; a[x]=a[x]$0 } END { for(x in a) print x,a[x] }' |
FBtr0344900    42774    384    43231   144

5) compute the new cigar, calculating the intron length as the arithmetic formula intron_length= (next_exon_start_coordinate - exon_length - previous_exon_start_coordinate), in this simple case above: intron_length= 43231-384-42774
$ awk '{ printf("%s", $1) }; { for (i=4; i<=NF; i+=2) { printf("%s%d%s%d", OFS, $(i-1), OFS, $i - $(i-1) - $(i-2)) } }; { printf("%s%d%s%s", OFS, $NF, OFS, RS) }'
FBtr0344900 384 73 144

6) ideally, with some method that I will figure out, I will modify the CIGAR string (adding alternate M,N at the end of each field except first one), this is how the final file should look like:
FBtr0344900 42774  384M73N144M

The issue of my basic approach are that:

I am not sure how to account for the 1-based SAM : should I add +1 at every exon_start_coordinate? Does not look like, as the STAR output has the exact same cigar string that I calculated on the STAR output.
BIG issue: This will only work for the reads mapped in forward strand: how to make it feasible with reads mapped in the reverse strand? If I keep my current approach, I will have negative intron sizes...

Any suggestion is welcome!

I strongly encourage you to not do this with normal unix tools and to instead code something in python (or whatever language you prefer).

I would use an aligner that outputs proper SAM, such as star, gmap, spaln or minimap2.

yes, I already run the same alignment with GMAP, STAR. we have suspect that LAST might work quite well, but as for now the output is not directly helpful. I will try minimap2 (nice to know you can use it for RNA also)

@user172818 minimap2 works amazingly good... many thanks for redirecting me to it

Just use minimap2 in split alignment mode to realign the reads.
If that is not an option, then you could try using pysam to modify the CIGAR strings. I do not recommend this, as there are many opportunities for subtle bugs because the SAM spec is complex. You would need to:

Sort the BAM on read ID so that you can efficiently retrieve reads that you want to merge
Iterate through the BAM with pysam while loading all reads of the same read ID
For each set of reads with the same read ID:

Order the reads of the same read ID by the number of hard-clipped bases at the start of the CIGAR string.
Build a new CIGAR string from the ordered CIGAR strings
Merge the meta tags from the ordered reads
Write out a single, merged read

I did that at the end. At the time, I just wanted to reformat the non-spliced alignments from LAST (for example) into spliced alignments. Why? Benchmark the most tools possible. But it's ok, I gave up at the end and just used tools that already gave that spliced alignment directly.

That's great. However, the point of SE is to also help others that might have similar issues in the future by providing answers to questions so that users can upvote good answers. If you feel that my answer is a good one to the question you posed, then please consider up voting and accepting it.

If you feel that this answer does not answer the question you really had, then please consider rewriting your question for clarity.


text:
JSmol - hide "JSmol" logo
I'm using JSmol for my SimRNAweb server. I would like to hide "JSmol" logo from all small visualizations that I have there. I would like to just put information about JSmol under the panel in text.

I wanted to hide this JSmol little logo. Now I'm done :-)

OK, I got an answer here. Use:
frank off
to get:


text:
Meaning of BWA-MEM MAPQ scores
Does anyone know what the MAPQ values produced by BWA-MEM mean?
I'm looking for something similar to what Keith Bradnam discovered for Tophat v 1.4.1, where he realized that:

0 = maps to 5 or more locations
1 = maps to 3-4 locations
3 = maps to
  2 locations
255 = unique mapping

I'm familiar with the notion that MAPQ should be theoretically be related to the probability of an "incorrect" alignment (10^(-MAPQ/10)), although this is vaguely-specified enough that aligners in practice tend to actually just use something like the above.
Note that this is different than the BWA MAPQ scoring interpretation, because BWA-MEM gives MAPQ scores in the range $[0,60]$, rather than $[0,37]$ as has been established for BWA.

First of all, if you want to understand mapping quality (mapQ), ignore RNA-seq mappers. They often produce misleading mapQ because mapQ is not important to RNA-seq anyway.
Strictly speaking, you have two questions, one in the title: the meaning of mapQ; and the other in a comment: how mapQ is computed. On the meaning, mapQ is nearly the same as baseQ – the phred scaled probability of the alignment/base being wrong. It often amuses me that we question mapQ but take baseQ for granted. BaseQ is also scaled and discretized differently; even fewer people know how Illumina/pacbio/nanopore/historical sequencers estimates baseQ.
On the second question, Section 2 of the MAQ supplementary explains the theoretical aspects of mapQ, which is still correct today. Briefly, mapping quality consists of three components: 1) the probability of contamination; 2) the effect of mapping heuristics and 3) the error due to the repetitiveness of the reference. Only 3) can be modeled theoretically.
In case of bwa-mem, if we assume the matching score is 1, type-3 error is estimated with:
$$10/\log10\cdot[\log4\cdot(S_1-S_2)-\log n_{\rm sub}]$$
where $S_1$ is the best alignment score, $S_2$ is the second best and $n_{\rm sub}$ is the number of suboptimal alignments. Factor $\log 4$ comes from the scoring matrix. Factor $10/\log10$ is the Phred scale. This equation assumes gap-free alignment and is very close to Section 2.5.2 in the MAQ supplementary. It is ok-ish for short reads, but often overestimates for long reads. I am not aware of a practical approach in general cases. In addition to this method, you can estimate mapQ by read simulation: just try to find a function such that it fits empirical mapQ. Some have tried machine learning, too.

This has an ambiguous definition, at least because it doesn't define what happens when there are two equally similar and optimal matches in different parts of the reference index. If "second best" is allowed to have the same alignment score as "best", then this score reduces to $$10/\log10\cdot[-\log n_{\rm sub}]$$ in such a case.

Beautiful. Do you have a reference for the type-3 error estimation expression for bwa-mem?

@gringer Read the MAQ supplementary. It is an approximate on the condition that $S_1\gg S_2$. Furthermore, you do need to consider type-2 errors and pairing. The real estimate is full of ad hoc. ijoseph, that would the MAQ paper.

Those numbers are not arbitrarily picked (well... maybe 255/60/40 is arbitrarily picked).
To convert from log10 Q values like these (also used for error rates in FASTQ files) to probabilities, divide the number by 10, negate it, then raise 10 to the power of the result.
Another way of looking at it is to consider the decade to mean the number of 9s in the [mapping] accuracy, e.g.:
Q10: 90%
Q20: 99%
Q30: 99.9%
Q40: 99.99%
With the formula adjusted to make all the numbers in between a smooth transition.
With this transformation, a MAPQ score of 3 corresponds to a probability of an incorrect mapping of almost exactly 50%, which is as confident as you can be about a specific mapping if you know that it maps equally well to two locations:
> 10^-0.3
[1] 0.5011872

A MAPQ score of 2 is an incorrect mapping probability of 63% (i.e. close to 3 different locations), and a MAPQ of 1 is an incorrect mapping probability of 79% (close to 5 different locations). These fit with the numbers that you have described.
However, as you point out, BWA doesn't do this; it seems to ignore the number of mappable locations when considering the MAPQ score. I explore the reason why this is a problem in my question about multi-mapped reads.
I notice that there's a bit of an attempt to answer this in a similar previous question:

$C * (s_1 - s_2) / s_1,$
where $s_1$ and $s_2$ denoted the alignment scores of two best alignments and C was some constant.

See the answer by user172818 for more insight into the details of BWA-MEM MAPQ scores.

"Similar" question is generous, that's basically a duplicate of my question that I didn't do due diligence to find!

So it basically comes down to understanding this uncommented C code. Maybe I'll try contacting Heng Li directly. https://github.com/lh3/bwa/blob/master/bwamem.c#L952-L976

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