text:
How do I find identical sequences in a FASTA file?
I want to create a database for a proteomics study. Therefore, the mapping from a given sequence to a protein ID has to be unique. I am wondering whether there is already a built-in function in Biopython for that, but I could not find any. The duplicated sequences should be merged into a single entry and there respective IDs should be just concatenated (see expected outcome below).
Using the following example input
>prot_id_1
MAWIGLDISKLFVENSRDAAA
>prot_id_2
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS
>prot_id_3
MRTRPSRDFHYIVGRITRESPEEKGIFYVFH
>prot_id_4
MEMVLSSANPSTTFADSYVV
>prot_id_5
MAWIGLDISKLFVENSRDAAA
>prot_id_6
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS
>prot_id_7
MAWIGLDISKLFVENSRDAAA

I currently do it as follows:
from Bio import SeqIO
import numpy as np
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import IUPAC

# Read entries into dictionary, gi: sequence
record_dict = SeqIO.to_dict(SeqIO.parse("test.fasta", "fasta"))

# Extract sequence as string
record_dict = {gi: str(record_dict[gi].seq) for gi in record_dict}

# Check for duplicated sequences
seq_id_dict = {}
for k, v in record_dict.items():
    seq_id_dict.setdefault(v, []).append(k)

# Filter sequences that appear more than once
filter_dict = {seqi: li for seqi, li in seq_id_dict.items() if len(li) > 1}

# Merge identifiers to a single one
merged_ids = {"_".join(sorted(li)): seqi for seqi, li in filter_dict.items()}

# Remove all entries that are associated with duplicate IDs
for gi in set(np.concatenate(filter_dict.values())):
    record_dict.pop(gi, None)

# Add the merged IDs again
record_dict.update(merged_ids)

# Create list of Seq objects
final_seq = [SeqRecord(Seq(seqi, IUPAC.protein), id=gi, name='', description='') for gi, seqi in record_dict.items()]

# Write file
SeqIO.write(final_seq, 'test_output.fasta', 'fasta')

which gives me the desired output:
>prot_id_4
MEMVLSSANPSTTFADSYVV
>prot_id_1_prot_id_5_prot_id_7
MAWIGLDISKLFVENSRDAAA
>prot_id_3
MRTRPSRDFHYIVGRITRESPEEKGIFYVFH
>prot_id_2_prot_id_6
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS

Is there a smarter way of accomplishing this?

The trick would be to swap the key in the dictionary to be the sequence itself.
Also I would recommend using a different separator that "_" since that is what the current ids have so that you can easily separate the individual ids from the concatenated id. I used a pipe "|" in this example. Also I just manually wrote the FASTA output instead of using Biopython to do it.
from Bio import SeqIO
from collections import defaultdict

dedup_records = defaultdict(list)
for record in SeqIO.parse("test.fasta", "fasta"):
    # Use the sequence as the key and then have a list of id's as the value
    dedup_records[str(record.seq)].append(record.id)
with open("Output.fasta", 'w') as output:
    for seq, ids in dedup_records.items():
        # Join the ids and write them out as the fasta
        output.write(">{}\n".format('|'.join(ids)))
        output.write(seq + "\n")

Yes, that works fine and is more condensed than my solution (your dedup_records is the same as seq_id_dict). But true, all the additional filtering is not needed; Friday afternoon... :)

(I wait with acceptance for a while just to see if there is already some kind of built-in for that)

There is no built in function that I know of that does the id concatenation.  If it did exist it would likely be doing the same logic.

Seems like a good approach if fitting the file into memory is ok, but why not use biopython to write also? This asymmetry seems odd, plus that way you won’t dump all the sequence on a single line

Mainly because I find the BioPython interface cumbersome to write out sequences in a fasta format when the data is not a Seq object already.  @cleb can use Biopython if he wants the main changes are in the use of the defaultdict with the swapped key.

I see, well still consider using record.description not record.id to ensure you get the whole headers. Also you could save memory by storing the sequences as a hash rather than the full sequences

OP referred to the protein ID so thats what I used.  I agree that hashes could be useful to handle large input amounts.

@Chris_Rands: Mind adding the hash solution!? In my case I only have the IDs, no further description but if you could add a more detailed answer I am more than happy to upvote it.

@Cleb I meant to say, the reason I didn't add a hash solution is because it would be tricky to resolve hash collisions (although they'd be very rare)

@Chris_Rands: No worries, the solutions posted are currently sufficient; would have just been a 'nice-to-see' type of thing...

Just as a small variation to @Bioathlete's answer in case you want to write the fasta using Biopython (e.g. to add names and description):
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import IUPAC
from collections import defaultdict

dedup_records = defaultdict(list)
for record in SeqIO.parse("test.fasta", "fasta"):
    # Use the sequence as the key and then have a list of id's as the value
    dedup_records[str(record.seq)].append(record.id)

# this creates a generator; if you need a physical list, replace the outer "(", ")" by "[" and "]", respectively
final_seq = (SeqRecord(Seq(seqi, IUPAC.protein), id="|".join(gi), name='', description='') for seqi, gi in dedup_records.items())

# write file
SeqIO.write(final_seq, 'test_output_final.fasta', 'fasta')

This will give the desired output:
>prot_id_1|prot_id_5|prot_id_7
MAWIGLDISKLFVENSRDAAA
>prot_id_2|prot_id_6
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS
>prot_id_4
MEMVLSSANPSTTFADSYVV
>prot_id_3
MRTRPSRDFHYIVGRITRESPEEKGIFYVFH

I think you mean it creates a generator not an iterator.

@Bioathlete: Of course, fixed, thanks for spotting it.

using datamash ans some awk
$ cat input.fa | paste - -  |\
cut -c2- | sort -t $'\t' -k2,2 |\
datamash  collapse 1 -g 2  | tr "," "_" |\
awk '{printf(">%s\n%s\n",$2,$1);}'

>prot_id_1_prot_id_5_prot_id_7
MAWIGLDISKLFVENSRDAAA
>prot_id_4
MEMVLSSANPSTTFADSYVV
>prot_id_2_prot_id_6
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS
>prot_id_3
MRTRPSRDFHYIVGRITRESPEEKGIFYVFH

This can also be done quite easily with awk:
awk '{if(NR%2){printf("%s\t",$1)}else{print}}' filename | sort -k2,2| awk '{if($2==SEQ){gsub(">","",$1);ID=ID"|"$1} else{if(SEQ!=""){printf("%s\n%s\n", ID,SEQ);}SEQ=$2;ID=$1;} }END{printf("%s\n%s\n", ID,SEQ)}'
The first awk converts the fasta file to a tab separated file with format ID\tSequence, which is then sorted by sequence by sort. The last awk goes through the sorted file looking at the sequences: if the sequence in the current line is the same as that in the previous line, it appends its id to the ID variable. Otherwise (i.e. we have encountered a new sequence), it prints the previous concatenated ids and the corresponding sequence.

Here is my free program on Github  Sequence database curator
(https://github.com/Eslam-Samir-Ragab/Sequence-database-curator)
It is a very fast program and it can deal with:

Nucleotide sequences
Protein sequences

It can work under Operating systems:

Windows
Mac
Linux

It also works for:

Fasta format
Fastq format

Best Regards

Hello Eslam, thanks for your contribution. It would be nice if you could also include how to use your program to achieve what the original poster has asked for... Here we try to be as concrete as possible both when asking and answering :-).

I agree with @KamilSJaron: Please show how to use your package to solve the issue from above (ideally with the same input). That will help others to get started. Thanks!

If you have a very big fasta file(60 GB) why use SeqIO.parse and write the output file using biopython?
Why not use SimpleFastaParser instead (you get a tuple(name, seq)) and you dont need to wait so long as using the parser and using Seqrecords.
What do you think about it guys?
Thanks
Paulo
Kind of:
from Bio.SeqIO.FastaIO import SimpleFastaParser
from collections import defaultdict

dedup_records = defaultdict(list)
for (name, seq) in SimpleFastaParser(filename):
    # Use the sequence as the key and then have a list of id's as the value
    dedup_records[seq)].append(name)
with open(filename, 'w') as output:
    for seq, ids in dedup_records.items():
        # Join the ids and write them out as the fasta
        output.write(">{}\n".format('|'.join(ids)))
        output.write(seq + "\n")


text:
Highlight CpG sites with Texshade
I am using Texshade to produce figures for MSA. Upon reading the documentation I could not find anyway to highlight how conserved CpG sites within the alignment. Does anyone know if this is possible? If not is there another MSA visualizer I can use that can produce a figure like this. Thanks!

The best solution I could find was to use the \shaderegion function on all the sequences so it ended up looking something like this:
\shaderegion{1}{CG}{Red}{OrangeRed}

\shaderegion{2}{CG}{Red}{OrangeRed}

...

Then I cast the rest of shading similarity to grey so the CG pairs were more easily distinguishable.
\shadingcolors{grays}

Two examples one with \fingerprints and one regular.

source: Stack Exchange

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

sort: votes


text:
What is local realignment and what is the problem it solves?
I am trying to understand local realignment but I could not get a clear idea of what is the problem solved by it.
For example, reading Homer and Nelson (2010):

Because alignment algorithms map each read independently to the reference genome, alignment artifacts could result, such that SNPs, insertions, and deletions are improperly placed relative to their true location. This leads to local alignment errors due to a combination of sequencing error, equivalent positions of the variant being equally likely, and adjacent variants or nearby errors driving misalignment of the local sequence. These local misalignments lead to false positive variant detection, especially at apparent heterozygous positions. For example, insertions and deletions towards the ends of reads are difficult to anchor and resolve without the use of multiple reads. In some cases, strict quality and filtering thresholds are used to overcome the false detection of variants, at the cost of reducing power [13]. Since each read represents an independent observation of only one of two possible haplotypes (assuming a diploid genome), multiple read observations could significantly reduce false-positive detection of variants.

I don't understand how the different artifacts (such as SNPs) can result from the alignment.

Homer, Nils, and Stanley F. Nelson. "Improved variant discovery through local re-alignment of short-read next-generation sequencing data using SRMA." Genome biology 11.10 (2010): R99.

SNPs are likely to be created and InDels are likely to be missed. Suppose you have a read, ACTGACTGACTGTAC and you align it to a reference sequence ACTGACTGACTGTTAAGAACGACTACGAC. If you aligned that, you would either get:
ACTGACTGACTGTac    (lower case denotes soft-clipping)
ACTGACTGACTGTTAAGAACGACTACGAC

or
ACTGACTGACTGTAC    (N.B., you've created some possible SNPS)
ACTGACTGACTGTTAAGAACGACTACGAC

But what happens if you look at a bunch of reads? See below:
ACTGACTGACTGTAC
  TGACTGACTGTACTA
    ACTGACTGTACTAAG
ACTGACTGACTGT--TAAGAACGACTACGAC

I've made the 2 base insertion obvious. As you can see, using many alignments, we can "realign" our original read to see that it actually ends in a 2 base insertion, something that's completely impossible to find with a single read. This realignment is part of what the GATK haplotype caller does (it makes a local de Bruijn graph). That's also what my MethIndelRealigner does. I have another real-world example in its documentation.
Raw alignments in IGV:

Alignments after realignment:

One can find plenty of similar scenarios that result in SNPs rather than simply losing coverage for InDels.

Awesome.. I was not considering the problem of mapping each read independently. Now it's clear, thanks!

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


text:
How to query a network for a bacterium, specifically Streptomyces caatingaensis or Streptomyces thioluteus?
I have a list of gene/protein IDs and I could not query a network using String, BioGrid, Intact and IIS, for Streptomyces thioluteus.
The bacteria I am working with is the Streptomyces caatingaensis. Since it is newer than others, I thought I would find an interaction network looking for Streptomyces thioluteus (which shares one pathway I am analyzing in common, others do not share). As a result, I have this list of proteins, but I could not query any network from the IDs. Sometimes I could not find the bacterium in the tools and when it was found, as in IIS the result was zero interactions.
Is there anything I could use to have some network that could be used as insight for the analysis?

I realized that some tools work better with Gene Names instead of some type of ID.
Using IIS I was capable of creating a network using this list of proteins (Streptomyces caatingaensis) by querying proteins by Genes Names having as base the E. Coli.
The tool has 3 main tabs: proteins, annotation, and interactome.
Search the Gene Names in proteins section, go to Interactome, set the parameters and create the network.
P.S. If you are not used to the tool you should read the research article available on the home page before using the network for research purposes.


text:
Matlab scatter diagram
I have two columns of data, x and y, how can I make a ‘thermodynamic diagram’ with matlab? The z axis is the density value of dots, like below figure.

The x-axis is angle, y-axis is distance.
I want to construct a heatmap like this figure:

Hi, how is this related to bioinformatics? What have you tried? I don't know if you want your data like in the second plot or like in the first one. Could you provide a some data to test solutions?

A reverse image search on google leads to the exact code used to produce the second image:
[Xcoarse, Ycoarse] = meshgrid([0 1 2 3], [0 1 2 3]);
[Xfine, Yfine] = meshgrid(linspace(0,3,3000), linspace(0,3,3000));
DataCoarse = [ 1 2 4 1; ...
           6 3 5 2; ...
           4 2 1 5; ...
           5 4 2 3];
DataBicubicFine  = interp2(Xcoarse, Ycoarse, DataCoarse, Xfine, Yfine, 'bicubic');
figure
surf(Xfine, Yfine, DataBicubicFine); shading flat; colormap(jet); view(0, 90)
 hold on; plot3(Xcoarse, Ycoarse, 10*ones(size(Xcoarse)), 'k.', 'MarkerSize', 20)
colorbar
print -dpng -r300 BicubicInterpolationExample.png


text:
Converting a VCF into a FASTA given a reference with Python, R
I am interested in converting a VCF file into a FASTA file given a reference sequence with Python or R.
Samtools/BCFtools (Heng Li) provides a Perl script vcfutils.pl which does this, the function vcf2fq (lines 469-528)
This script has been modified by others to convert InDels as well, e.g. this by David Eccles
./vcf2fq.pl -f <input.fasta> <all-site.vcf> > <output.fastq>

Given that converting a VCF into a FASTA given some reference is a somewhat standard procedure in bioinformatics, are there standard R or Python implementations used by the community?

I use the pybedtools package which is a wrapper around bedtools. Bedtools will work on VCF files in addition to BED and GTF. Why do you want a Python or R implementation instead of using Heng's solution?

"Why do you want a Python or R implementation instead of using Heng's solution?" Sometimes it's easier to integrate R code with R code, and Python code with Python

Based on the name and a quick look at the code, the function you mention produces fastq and not fasta. Do you want fasta or fastq?

@terdon Either works

Typically I find using the best tool is more important than what language it was written in.

You can whip up something quite easily in Python using pyfaidx, which allows you to access a FASTA file by genomic coordinates. So you'd just have to pull the coordinates from the VCF to get the sequence, then output in FASTA format with whatever unique identifier for each sequence you want to use.
Edit: The FastaVariant class is new in pyfaidx, I hadn't seen it before. It makes this even easier, as you won't have to manually mess with the VCF to get the positions. Something like the following should work:
consensus = FastaVariant('yourgenome.fasta', 'variants.vcf.gz', het=True, hom=True)

out = open("variants.fasta", "w")

for chrom in consensus.keys:
    for var in consensus[chrom].variants_sites:
        record = consensus[chrom][var-1:var]
        print(record, file=out)

out.close()

Hi Jared, could you expand your answer, with the code to use to perform this? It is my the first time with pyfaidx and I have problems finding how to use it. It seems something like consensus = FastaVariant('tests/data/chr22.fasta', 'tests/data/chr22.vcf.gz', het=True, hom=True) should work?

@Llopis yes, it should be something like that. The FastaVariant class is actually new, I wasn't aware it existed, so thanks! I've updated my answer with code I believe should do as intended.

It seems like FastaVariant ignores indels.

I tried the pyfaidx approach reported by Jared, but I got unhandled exceptions that I could not handle easily within the framework. I had bcftools lying around, so I tried bcftools consensus and it worked like a charm. However, you do need an indexed VCF.
bcftools consensus all-site.vcf.gz < input.fasta > output.fasta

NB: bcftools consensus has a few options specified with the --haplotype argument for choosing which alleles should be incorporated into the FASTA file.

Assuming you have samtools-indexed FASTA files for your reference genome, and have BEDOPS and htslib installed, you could use bed2faidxsta.pl to convert BED to FASTA:
$ vcf2bed < variants.vcf | bed2faidxsta.pl --fastaDir=/path/to/fasta > variants.fa

To index your FASTA files, you could do the following:
$ cd /path/to/fasta
$ bgzip -c my_genome.fa > my_genome.fa.bgz
$ samtools faidx my_genome.fa.bgz

Thank you. This wouldn't be too difficult to convert into Python/R, but I'm asking for scripts which have this already completed.

If you have Python available, you very probably already have Perl. In any case, you can use system in R to call external scripts (written in any language). See ?system for more details.

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


text:
Verify a predicted protein in one genome in a different genome of the same species
I have two genome assemblies of the same non-model species, call them Assembly 1 (generated from Illumina data) and Assembly 2 (generated from PacBio data).
For Assembly 1, I also have predicted proteome data, generated with EVM. Say there is a protein, call it Protein X, for which I have the fasta sequence from predictions, but suspect an error in Assembly 1 (because the Protein X has a long segment of auto-identity). For this reason, I would like to verify the assembly of the locus coding for Protein X in Assembly 2 and/or verify the presence of Protein X in the Assembly 2.
Data in my possession:

both genome assemblies (fasta)
raw rna-seq data
reference-based transcriptome assembly for each genome assembly (fasta/gtf)
set of predicted protein sequence based on Assembly 1 (fasta)

Things I was thinking of:

build the EVM predictions on the Assembly 2 and verify the identity of Protein X by simple fasta identity check
look for the protein sequence in the Assembly 2. I could use something similar to Reverse Translate to get the most likely nucleotide sequence from the input amino-acid sequence of Protein X, then align it to the Assembly 2 to check the coordinates of the locus it originates from - but I am not very sure of the relevance and/or correctness of this procedure
something allowing 'blast' of the protein sequence in the Assembly 2 (but I am not aware of any software that could do that)

Any ideas on how I could do that?

This should be very easy to do. Here are some options:

Use a tool like Exonerate or GeneWise that can align protein sequences to genomic DNA while attempting to model splice sites etc.
As you said, blast. There's nothing special in what you describe, really, this is what tBLASTn is for. Just use your protein as the query sequence and tBLASTn as the blast flavor which will take a protein query and match it against all 6 possible reading frames of a DNA database. Set up your assembly as the blast DB, and then blast as you would normally.

Of the two, using the first approach is much better since it will try to build a valid gene model for you instead of simply finding the regions of high homology as blast will do.
Finally, just for completeness, you can go old school and run a full de-novo gene prediction tool on your assembly. Something like GeneID for instance. But that's certainly overkill if you're just looking for a single gene.


text:
What happens if a major bug is discovered in a bioinformatic package that has been used in published literature?
Yesterday I was debugging some things in R trying to get a popular Flow Cytometry tool to work on our data. After a few hours of digging into the package I discovered that our data was hitting an edge case, and it seems like the algorithm wouldn't work correctly under certain circumstances.
This bug is not neccessarily the most significant issue, but I'm confident that it has occurred for other users before me, perhaps in a more insidious and less obvious way.
Given that since this tool was released, literature has been published that utiltised it, how does the profession handle the discovery of bugs in these widely used packages?
I'm sure that given the piecemeal nature of a lot of the libraries out there, this is going to happen in a much more significant way at some point (i.e. invalidating a large number of published results)
For some context I'm a programmer who's been working across various other parts of software development, and is quite new to Bioinformatics.

Have you reported the bug to the developers?  Is the tool still maintained?  Really it depends on the developers how the bug will be handled.  Sometime it will go out in the next point release and sometime it will be determined to be a low risk and low probability occurrence and the potential fix would introduce great risk so it will stay as a known issue.

I think this question might benefit from a more general tag. This might come to other packages/libraries in C++, python, perl, S, or whatever language. Maybe a good tag would be [tag:software-quality]?

Might also be a relevant question over on Academics.SE

@WayneWerner there is a relevant question there

I'm not sure I understood the question. "it seems like the algorithm wouldn't work correctly under certain circumstances" - do you mean that the algorithm is incorrect from a biologist's point of view, but was implemented correctly in the package, or that the general algorithm would compute everything properly, but its R implementation in the specific package does not meet the specification?

@Llopis I wanted to create a new general tag, unfortunately I didn't have the reputation at the time. I'll have a think about it and retag if I can.

@Bioathlete Yep, I probably should have been more clear - I've reported and also fixed the bug, but I was musing about the larger implications on public literature in a more general sense.

@rumtscho What I mean by that statement is that for some cases it works as expected, and for other cases it doesn't. I.e. for certain data sets the results would be erroneous, but for other data sets it works fine.

@NicBarker You seem to have created issues on the repository package but you haven't arised these concerns in the support site of Bioconductor. I would first report there your findings. You can always fork the most updated version and modify the code yourself.

Thanks for the tip @Llopis, I'm still finding my feet and figuring out best practise in the field.

I prefer to treat software tools and computers in a similar fashion to laboratory equipment, and in some sense biology in general. Biologists are used to unexpected things happening in their experiments, and it's not uncommon for a new discovery to change the way that people look at something. Things break down, cells die off quicker on a Wednesday afternoon, results are inconsistent, and that third reviewer keeps on about doing that thing that's been done a hundred times before without anything surprising happening (just not this time). It's a good idea to record as much as can be thought of that might influence an experiment, and for software that includes any input data or command line options, and especially software version numbers.
In this sense, a discovered software bug can be treated as a new discovery of how the world works. If the discovery is made public, and other people consider that it's important enough, then some people might revisit old research to see if it changes things.
Of course, the nice thing about software is that bugs can be reported back to the creators of programs, and possibly fixed, resulting in an improved version of the software at a later date. If the bug itself doesn't spark interest and the program gets fixed anyway, people unknowingly use newer versions, and there might be a bit more confusion and discussion about why results don't match similar studies carried out before the software change.
If you want a bit of an idea of the biological equivalent of a major software bug, have a look at the barcode index switching issue, or the cell line contamination issue.

Great answer (+1), except for introducing the IMHO unnecessary and confusing distinction between 'corner' and 'non-corner' case bugs. A code is not like a lab equipment in that you do not have a general understanding over the state-space of how the tool is used (e.g. you do not need to label all lab equipment with 'this tool might not work in outer space').  Another way to put this: all bugs are necessarily corner cases to the person that wrote the code.

One person's bug could be another person's (feature)[https://xkcd.com/1172/]. A person doesn't need to know (and often can't know) everything about software in order to use it, and developers don't need to fix every known bug in order for a program to be useful.

Thanks for those two links @gringer, definitely an interesting read.

"people unknowingly use newer versions, and there might be a bit more confusion and discussion about why results don't match similar studies carried out before the software change." Which is why you should note the versions of the software you used in your publications. However small and however subtle you do it, make sure it's in there.

While I agree with gringer's answer I would first report the bug to the maintainers with a test case (if possible) in order to  have the bug corrected.
The bug might be an erroneous implementation in a corner case or in a more common use case of the software. But it is a good opportunity to collaborate to increase the quality of the software, either by introducing tests (or expanding them), by re-evaluating the state of the art with other tools, or to communicate with users of other software that depends on the tool.
Depending on the bug and the time one is able to spend, do a collaboration to study the effect of such a bug in the published studies using the software.
Sometimes the bug doesn't affect much, but sometimes it can have quite an impact which can be difficult to estimate. For instance, I found that the reactome.db package didn't have the same annotations that the available online database (in version 1.58.0), due to an incorrect mapping. However, it was impossible for me to evaluate the damage of such an error on the published literature, because gene set enrichment tools (usually) don't report the data they are based on and the versions of the databases used are often under reported in the literature.

Don't worry about it! That's the beauty of a collaboratively edited site! And I cheat: I'm a native speaker :)

If you have a software development background you probably know about filing bug reports. Your question seems to be more about how to correct the scientific literature which has been affected by the bug.
The traditional way that this has been done is by writing an paper that explains the bug and gives examples of how it may impact published results (see gringer's examples of index switching and HELA contamination). This way the community will be aware that older results using the implicated method need to be verified before building upon them. Unfortunately, awareness of these issues doesn't always diffuse through the community as quickly as it should, but it's probably still the best tool we have to ensure that science is self-correcting.

I wonder if you could publish enough papers purely from bug-hunting to earn your way to a normal career in academia ;)

In extreme cases, it can destroy careers:

One of the most spectacular flameouts in science happened last year.
  In a short letter (barely over 300 words long) published in Science in
  the very last issue of 2006, Geoffrey Chang, a crystallographer,
  retracted 3 Science articles, a Nature article, a PNAS article and a
  JMB article. The sum of 5 years of work was destroyed, apparently,
  over a single sign error in a data-processing program.
  – A Sign, a Flipped Structure, and a Scientific Flameout of Epic Proportions

Most cases are not that extreme. In fact, realistically, I would guess that most papers have some data processing errors of some kinds, yet still we muddle on towards truth. And this isn't unique to software – the whole idea of independent replication in science is that we're always probably messing things up in ways we don't realize, so we need to re-do things and hope that the errors are different this time. But software is unusual in how tiny errors can cause very large effects, and often the effects are hard to detect. If you drop a test tube on the floor, you'll probably notice. If a piece of new software gives incorrect results... well, we didn't know what the results were supposed to look like, that's what makes it science! There are things you can do to mitigate this somewhat, but to some extent it's just inherent to the scientific process.
I still like to show that blog post to new grad students to scare them into writing tests. I actually think it's a good habit to, whenever you start using a new software package, start by writing some tests for it – it's a good way to check your assumptions and understanding about what it actually does, and even if it's a popular package by eminent researchers, you'll still sometimes be shocked by what you find.

It depends also where the package is available from. I mean if it is from bioconductor, authors and maintainers are responsible to respond to bug reports. It could be reported at support sites, so others are also aware of the bug. If the maintainer doesn't fix it, it might lose its spot in bioconductor. However, other sites, such as CRAN or Sourgeforge might be less stringent.

I always wonder what the "Response to bug reports" means, but that would be for the bioconductor mailing list. I read it as they don't need to be corrected, but just have a (reasonable) response. And AFAIK it is not enforced, I reported several bugs and I didn't get any response from the maintainers :(

In my opinion an adequate response would be a bug fix. Not just a "thanks, you'll hear from us" response on the support site. It's a pity to hear that you don't even got a response, I thought bioconductor was stringent is these matters.

authors: ə̷̶̸͇̘̜́̍͗̂̄︣͟


text:
Why does a very strong BLAST hit get lost when I change num_alignments, num_descriptions or max_target_seqs parameter?
Disclaimer: This is a self answered question for documentation purpose and I adapted this from the following github gist. Especially from users terrycojones and peterjc as well as sujaikumar who raised the issue.

I have a strange situation. I have two BLAST databases (20170821 and 20170824) that both contain a certain subject sequence LC074724. When I query them both with a given query, Q, in one case (20170824) there is a very strong match and in the other (20170821) no match at all.
Why is there not even a lower quality match found although the databases are nearly the same?
BLAST version

$ blastn -version
blastn: 2.6.0+
 Package: blast 2.6.0, build May 11 2017 22:22:40

Both databases contain an identical subject sequence LC074724:
$ blastdbcmd -entry LC074724 -db 20170824 | md5
fcfbc167c7f81dfd75aad1bf3db6220b

$ blastdbcmd -entry LC074724 -db 20170821 | md5
fcfbc167c7f81dfd75aad1bf3db6220b

Only difference for the databases is 20170821 is just 20170824 with an extra 643 sequences. Rest is completely identical.

I'm trying to figure out why. This will be a longer read, Tl;dr at the end:
Doing a match
There is a good match of Q against LC074724 in 20170824:
$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170824 -query Q.fasta \
  -task blastn -max_hsps 1 | grep LC074724
Q   LC074724.1  2581

Note that the bitscore of the LC074724 match is 2581.
But there is no match of LC074724 at all in 20170821:
$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170821 -query Q.fasta \
-task blastn -max_hsps 1 | grep LC074724

The option -outfmt '6 qaccver saccver bitscore' tells BLAST to print the query accession number and version, the subject accession and version, and the bitscore of the match.
What are the best hits for 20170821
So what are the bitscores of the Q matches in 20170821?
$ blastn -outfmt '6 bitscore' -db 20170821 -query Q.fasta -task blastn -max_hsps 1 | sort -nr | head
3301
2616
2590
2590
2563
2554
2547
2544
2540
2540

The bitscore of Q against LC074724 in 20170824 (2581) would have placed it in 5th position in the top 10 matches of Q, but LC074724 is not matched at all in 20170821! The lowest bitscores matched in 20170821 is 2354.
What's the difference between the databases?
Get the ids of the subjects in each database
$ blastdbcmd -entry all -db 20170824 | egrep '^>' | cut -c2- | sort > 20170824.ids
$ blastdbcmd -entry all -db 20170821 | egrep '^>' | cut -c2- | sort > 20170821.ids

$ wc -l 2017082?.ids
    9754 20170821.ids
    9111 20170824.ids
   18865 total

20170821 has 643 sequences that are not in 20170824:
$ comm -23 20170821.ids 20170824.ids | wc -l
     643

20170824 has no sequences that are not in 20170821:
$ comm -13 20170821.ids 20170824.ids | wc -l
       0

and the two have 9111 sequences in common:
$ comm -12 20170821.ids 20170824.ids | wc -l
    9111

These numbers match with the numbers from wc -l above.
So 20170821 is just 20170824 with an extra 643 sequences.
Maybe LC074724 matches, but with a low bit score?
BLAST only shows the top 500 matched subjects by default. Let's ask it to show more:
$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170821 -query Q.fasta \
  -task blastn -max_hsps 1 -max_target_seqs 10000 > 20170821-10000-matches.txt

And....
$ head 20170821-10000-matches.txt
Q   Q   3301
Q   X234A   2616
Q   CS388973.1  2590
Q   DM059402.1  2590
Q   LC074724.1  2581
Q   KJ843188.1  2576
Q   AB697496.1  2576
Q   AB697499.1  2576
Q   AB697503.1  2576
Q   AB697508.1  2576

There it is! And it's in 5th position, just as it should be (based on the bitscores discussed above).
Just to confirm, taking away the -max_target_seqs option gets no match:
$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170821 -query Q.fasta \
  -task blastn -max_hsps 1 | grep LC074724
$

Three confusingly named BLAST options
I've always found the names and descriptions of the following three BLAST command-line options very confusing.
From the Formatting options section of blastn -help:

-num_descriptions =0>   Number of database sequences to show one-line descriptions for   Not applicable for outfmt > 4
  Default = `500'    * Incompatible with:  max_target_seqs
-num_alignments =0>   Number of database sequences to show alignments for   Default = `250'    * Incompatible with:
  max_target_seqs

And from the Restrict search or results section:

-max_target_seqs =1>   Maximum number of aligned sequences to keep   Not applicable for outfmt <= 4   Default = `500'    *
  Incompatible with:  num_descriptions, num_alignments

From the section names, it seems like only the latter option (-max_target_seqs) would have any effect on the search, and that the former two are just about what is displayed.
But... using either just -num_descriptions or -num_alignments with a big value also gets a match:
$ blastn -db 20170821 -query Q.fasta -task blastn -max_hsps 1 -num_descriptions 10000 | grep LC074724
LC074724.1  Hepatitis B virus DNA, complete genome, isolate: p621     2581    0.0
>LC074724.1 Hepatitis B virus DNA, complete genome, isolate: p621

$ blastn -db 20170821 -query Q.fasta -task blastn -max_hsps 1 -num_alignments 10000 | grep LC074724
LC074724.1  Hepatitis B virus DNA, complete genome, isolate: p621     2581    0.0
>LC074724.1 Hepatitis B virus DNA, complete genome, isolate: p621

So BLAST is apparently finding the LC074724 match even without the high -max_target_seqs option, but it's not displaying it unless there is a high -num_alignments or -num_descriptions value. (Keep in mind, this is the 5th best hit! Not somewhere >500)
Just to confirm, when none of the 3 options are given and no -outfmt option is either, the match is not found:
$ blastn -db 20170821 -query Q.fasta -task blastn -max_hsps 1 | grep LC074724
$

Tl;dr:
In summary, the options
 -evalue
 -max_target_seqs
 -num_alignments
 -num_descriptions

are NOT just filters applied to the results, but instead directly influence the results found.
This means you can End up with a bacterial sequence as BEST hit when using a value x, and with no bacterial sequences at all in your top20 when using value y.
To quote the NCBI user service:

We don't consider this a bug, but I agree that we should document this
  possibility better. This can happen because limits, including max
  target sequences, are applied in an early ungapped phase of the
  algorithm, as well as later.

And

The cutoffs are applied during an ungapped extension phase, and some
  alignments may improve in score during gapped extension. The safest
  approach is to relax (increase) the threshold and deal with a larger
  result set.

This was totally unexpected behaviour for me and is as far as I know is nowhere documented in the BLAST help or similar.


text:
Quantification of EST/cDNA based retained intron transcripts
Note: this question has also been asked on Biostars
We are seeing a retained intron transcript event for some RNAseq samples, and we want to assess at the sequence level which intron retention events are actually occurring so that we can do protein domain predictions. For example, we have transcript level expression support for DROSHA-203, which according to Ensembl is a retained intron transcript that potential gives rise to multiple splicing variants via the cDNA sequence information:

In theory, the orange sites indicate where a mutation might occur which could obliterate a splicing recognition site, thereby causing an retained intron event. So in total there are 6 potential mutational events which might give rise to 4 possible retention events. Given that 4 intronic regions might be retained, the number of unique retention events is a combinatorics question equal to 2^4, e.g., intron1 retained only, vs intron1 + 2 retained, all possible unique combination, etc, correct?
Is there a better way to more granularly quantify these possible transcripts? Or do we have to PCR validate all of them to figure out which one's are actually occurring? Is there a way that I can look at the reads mapping to the corresponding intronic regions to make even a probabilistic assessment of which transcripts are present and to what extent they are expressed? The libraries are prepped for total RNA with a ribozero depletion, so I should still be able to quantify any enrichment to these intronic regions, correct?

Did you look up if there is any kind of LD between these positions? Maybe whenever a mutation in 359 occurs, there's also a mutation in the the previous splicing event.

You can do full-length whole-transcriptome sequencing using either the direct RNA kit or a cDNA kit (with or without PCR) on a MinION. This will allow you to quantify isoforms directly, rather than estimating by looking at read pairs.
However, it'd be a bit expensive to do it that way if you're looking at a single transcript (and aren't barcoding to multiplex samples on the same run).
If you've already done Illumina whole-transcriptome sequencing (e.g. strand-specific paired-end sequencing), tools like Kallisto, RSEM, or Salmon can be used to estimate isoform expression based on the mapped reads. These tools use the information about transcript placement and/or similarity to alter the observed read counts in an attempt to better represent reality.
Unfortunately, if you want these tools to quantify intronic sequence (i.e. the $2^4$ possibilities you have suggested), then those sequences need to be included in the sequences provided. They won't estimate frequencies for non-existent sequence.
As an alternative to including every possible gene model, something like a SuperTranscript model might work. I'm not familiar with Alicia Oshlack's method, but it at least gives the appearance that it might reduce the complexity of searching for multiple unannotated isoforms within a gene.

So the current transcript level quantification is performed with Salmon. Essentially salmon says that the retained intronic transcript DROSHA-203 (ie., ENST00000504133.5) is present in the samples. However when you look at the transcriptome fasta or the actual annotation outlined by Ensembl, it only lists the exonic portions of the transcript. My issue is quantifying which introns are being retained in the resulting transcript. As it stand it looks like there are multiple possibilities based upon the cDNA sequence provided through the annotation.


text:
How can I colour boxes in Gviz AnnotationTrack in R?
I'm learning the Gviz bioconductor package, I generate a plot as follows:
library(Gviz)
track <- AnnotationTrack(start=c(1,5,7), end=c(2,6,10), strand=c('*','*','*'), stacking="dense", showFeatureId=TRUE, id=c('red','blue', 'red'))
pdf(file=paste("test","pdf", sep="."))
plotTracks(track)
dev.off()

test.pdf looks like this:
However, what I want is for the boxes to be coloured rather than labelled with the colour word, how to do this?
Note that my real example is much more complex than this.

Did you try the fill argument?
Something like this:
track <- AnnotationTrack(start=c(1,5,7), end=c(2,6,10), strand=c('*','*','*'), stacking="dense", showFeatureId=TRUE, id=c('red','blue', 'red'), fill=c('red','blue', 'red'))

Thanks, seems to work, can't believe I missed that! What colours are accepted as arguments for fill?

Probably all colors, here a list of color names.

@Chris_Rands see ?plotTracks for the options accepted by plotTracks which will then refer you to ?settings for the various things you can tweak (long, long documentation). Finally, availableDisplayPars(track) will show you what you can set for track.


text:
Least present k-mers in the human genome
What are the least present k-mers in the human genome at different sizes?
Starting with k=4 and going up in size until k=10, what are the k-mers least seen (or not at all) in the human genome? I am only interested in the reference human genome, so I am discounting SNPs/Indels in the population.
If this is not pre-calculated somewhere, what tool is recommended to use that starts with the GRCh38 reference as input?

If you are interested in least present k-mers, you might be interested also in the concept of minimum absent words (these two things are related). There are many papers on this topic, see, e.g.: http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0029344. Another phrase to google is "word avoidance".

You can use the Jellyfish software to calculate the k-mer profiles up to length 31.
From the instructions in the user guide:
The basic command to count all k-mers is as follows:
jellyfish count -m 21 -s 100M -t 10 -C reads.fasta

To compute the histogram of the k-mer occurrences, use the histo subcommand (see section 3.1):
jellyfish histo mer_counts.jf

To query the counts of a particular k-mer, use the query subcommand (see section 3.3):
jellyfish query mer_counts.jf AACGTTG

To output all the counts for all the k-mers in the le, use the dump subcommand (see section 3.2):
jellyfish dump mer_counts.jf > mer_counts_dumps.fa

I am trying jellyfish now.

I tried jellyfish, and got the list of kmers. But I couldn't get the list of kmers of count=0.

You can also use R. I give you an example of only chr1 and only kmer=4.
library(BSgenome.Hsapiens.UCSC.hg38)
library(Biostrings)

genome <- BSgenome.Hsapiens.UCSC.hg38
kmers <- oligonucleotideFrequency(genome$chr1, 4)
kmers

m <- as.matrix(kmers)

m[order(m),]

Grabbing chromosomes for hg38:
$ wget ftp://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/*.fa.gz
$ for fn in `ls *.fa.gz`; do gunzip $fn; done

Via kmer-counter and Python, here's how to search for kmers of length 7 from chromosome chrY:
#!/usr/bin/env python

import sys
import subprocess
import itertools

k = 7
chr = 'chrY'
fastaFile = '%s.fa' % (chr)
kmerCmd = 'kmer-counter --fasta --no-rc --k=%d %s' % (k, fastaFile)

try:
    output = subprocess.check_output(kmerCmd, shell=True)
    result = {}
    for line in output.splitlines():
        (header, counts) = line.strip().split('\t')
        header = header[1:]
        kmers = dict((key,int(val)) for (key,val) in [d.split(':') for d in counts.split(' ')])
        result[header] = kmers
except subprocess.CalledProcessError as error:
    sys.stderr.write("%s\n" % (str(error)))

kmers = result[chr]
comp = {'A':'T', 'C':'G', 'G':'C', 'T':'A'}
for kmerList in itertools.product('ACGT', repeat=k):
    kmerKey = ''.join(kmerList)
    kmerCompKey = ''.join(reversed([comp.get(b,b) for b in kmerList]))
    if kmerKey not in kmers and kmerCompKey not in kmers:
        kmers[kmerKey] = 0

for key, val in sorted(kmers.iteritems(), key=lambda (key,val):(val,key)):
    sys.stdout.write("%s\t%s\n" % (key, val))

This script will print a two-column tab-delimited text file to standard output, where the first column is the 7mer and the second column is the count of that 7mer and its reverse complement over chrY (including zero-counts):
CGACGCG 20
CGTCGCG 20
CGCGATA 23
TACGCGC 25
...
AAATAAA 33521
TTTCTTT 34014
GAATGGA 35361
AATGGAA 36906
TTTTTTT 103093

With a couple tweaks, this script could run over a range of k and chromosomes for your reference genome.

authors: Hášíl Páůďýál


text:
How to find the sum of all sequences in a multifasta file
I have a multifasta file containing several genes from a genomic island. I want to know the total length of this island. How can I do this?

What have you tried? How many sequences do you have? How are the position of each sequence stored in your multifasta or elsewhere?

Most tools that calculate genome assembly summary stats will report total length. I use assembly-stats

grep -v ">" multi.fasta | tr -d '\n' | wc -c

grep -v ">" find all non header lines (i.e. without ">")
tr -d '\n' delete all newlines
wc -c count all characters

Here’s a function that I keep around for calculating the individual lengths of fasta sequences
falens(){
awk '/^>/ {if (seqlen){print seqlen}; print ;seqlen=0;next; } { seqlen += length($0)}END{print seqlen}' $1
}

And one I use for summing a column of numbers
sum(){
awk '{ sum += $1 } END { print sum}'
}

You can combine these by piping the output to do the sum for you:
falens myseqs.fa | grep -v “>” | awk '{ sum += $1 } END { print sum}'

(Can’t test this myself currently as I’m not at my laptop, but it’s based on 2 functions that work normally)


text:
How to colour multiple residues in Pymol?
We can use
color yellow, resi 3
color yellow, resi 3-8

to colour one residue or a range of residues. But what if I want to colour residue 3-8 and 10-20, how to write it?
Also, my protein has two chains, light chain (L) with 100 residues and heavy chain (H) with 100 residues. How to colour the residue 10-20, 30-40 of light chain, and the residue 5-10, 20-30 of heavy chain in one go?

What you’ve written isn’t valid Python code. Is this in Pymol or similar?

@KonradRudolph it is Pymol code (you can use python with pymol but in the GUI you can type pymol code)

I don't think that python is an appropriate tag for this post.  This really isn't a bioinformatics question.

there is no tag called "pymol", and where else would you recommend?

@Bioathlete why wouldn't this be a bioinformatics question? Where would you have such question otherwise?

It is a general question for the use of a tool.  It belongs in a mailing or user list for that tool.

of course each software may have its own mailing list, like I use Gromacs as well. But we should still be able to ask questions here @Mitra thank you for changing it to tag pymol. I thought I do not have the authority as a newbie.

Note: as I don't have any working version of pymol at the moment,
I am not able to test the solution properly.
However, an easy way is to first create a selection, and then to color it:
select toBecolored, resi 3-8
color yellow, toBecolored

Regarding the L/H chain part, if they are effectively annotated as L and H, something like the following should work:
select toBeColored, ((i. 10-20 or i. 30-40) and c. L ) or ((i. 5-10+20-30) and c. H)
color yellow, toBeColored

Combining both commands is probably possible as:
color yellow, ((i. 10-20 or i. 30-40) and c. L ) or ((i. 5-10+20-30) and c. H)

Pymol allows synonyms (eg. resi and i.  or chain and c.).
The following pages may be also helpful:

Pymol wiki: help on Selection
SourceForge

And the two following links, if you want to get a bit deeper in proper programming with python for pymol (as they can provide you with inspiration):

Pymol wiki: plugins
Pymol scripts

@Llopis: Thank you for the edit! It's a lot clearer now :)


text:
How to determine a protein's cellular location based on its sequence?
I am wondering about the appropriate workflow to determine a protein's cellular location based on its sequence.
Let's say I have a sequence like this from a fasta file
MTIKEHKVVYEAHNVKALKAPQHFYNSQPGKGYVTDMQHYQEMYQQSINEPEKFFDKMAK
EYLHWDAPYTKVQSGSLNNGDVAWFLNGKLNASYNCVDRHAFANPDKPALIYEADDESDN
KIITFGELLRKVSQIAGVLKSWGVKKGDTVAIYLPMIPEAVIAMLAVARIGAIHSVVFAG
FSAGSLKDRVVDANSKVVITCDEGKRGGKTINTKKIVDEGLNGVDLVSRILVFQRTGTEG
IPMKAGRDYWWHEEAAKQRTYLPPVSCDAEDPLFLLYTSGSTGSPKGVVHTTGGYLLGAA
LTTRYVFDIHPEDVLFTAGDVGWITGHTYALYGPLTLGTASIIFESTPAYPDYGRYWRII
QRHKATHFYVAPTALRLIKRVGEAEIAKYDTSSLRVLGSVGEPISPDLWEWYHEKVGNKN
CVICDTMWQTESGSHLIAPLAGAVPTKPGSATVPFFGINACIIDPVTGVELEGNDVEGVL
AVKSPWPSMARSVWNHHDRYMDTYLKPYPGHYFTGDGAGRDHDGYYWIRGRVDDVVNVSG
HRLSTSEIEASISNHENVSEAAVVGIPDELTGQTVVAYVSLKDGYLQNNATEGDAEHITP
DNLRRELILQVRGEIGPFASPKTIILVRDLPRTRSGKIMRRVLRKVASNEAEQLGDLTTL
ANPEVVPAIISAVENQFFSQKKK

I can feed it to NCBI's blastp which then gives me a list of genes. I can check e.g. the first entry which also provides a link to SwissProt. There I finally see that this protein is located in the nucleus and also the cytosol.
My questions are:
1) Is there a more straightforward way of achieving this information?
2) Let's say someone annotated this protein to be a mitochondrial one. How could I check that i.e. (dis)prove that?

Any update on this?

You could use Blast2go to find the predicted relevant gene ontologies of your sequences. If you look for the cellular compartment sub ontology you should find the location of such protein. I have never used this tool, but here is the paper describing how it works.
I wouldn't restrict to the first predicted match of blastp see if some more genes with a good fit are also located at the same sites. But otherwise your workflow seems right to me.
But to prove this is in X or Y you, or your team or collaborators, should always check with histochemistry if it really allocates on the predicted place. All the softwares predictions must be validated with some biological experiments to prove them.

Thanks for the suggestion, I will give it a try once I find time.

@Cleb I forgot to write that your worklow is how I would do it too.

You can run the sequence through tools like SignalP to look for cell surface signals and so on, which will give you a clue about where I is bound for.
None of these tools will give you an absolutely definitive answer though, often because we simply don’t have the experimental data.

Thanks, was not aware of this tool, will give it a try.


text:
IMPUTE2 missing filter option
I went through the manuals and can't find the meaning of the below flag:
-filt_rules_l 'filter==0'

Here is the link to manual for this flag: filt_rules_l
Can anyone explain?
IMPUTE version:
impute_v2.3.1_x86_64_dynamic

Edit: The full IMPUTE2 code is irrelevant, I am not using the command to impute, I received an imputed data and the script used to impute. This one flag was a mystery to me. I will update the post with impute software version when I find out.

I have never used this program, but reading the manual in the link it has to mean the following.
The manual says:

Writing -filt_rules_l 'eur.maf<0.05' on the command line would tell the program to remove any variants with eur.maf values less than 0.05 from the reference panel.

So, -filt_rules_l 'filter==0' would mean that the program removes any variants with filter values of 0 from the reference panel. In the legend file there should be a column named filter, which are filtered out when the value is 0.
Here is the example legend file, in your case there must be a custom additional column named "filter":
$ head ALL_1000G_phase1integrated_v3_chr1_impute.legend

id position a0 a1 afr.aaf amr.aaf asn.aaf eur.aaf afr.maf amr.maf asn.maf eur.maf
rs58108140 10583 G A 0.0407 0.1713 0.1311 0.2071 0.0407 0.1713 0.1311 0.2071
rs189107123 10611 C G 0.0142 0.0276 0.0140 0.0211 0.0142 0.0276 0.0140 0.0211

Thank you, are you suggesting that reference panel has column named "filter"? I will have to check. But I think "eur.maf" column is calculated on the fly from the reference file, need to check.

If you want to use your command it should. If not it won't work, and you need to find the column on which you want to filter.

Yes, sorry I have missed the _l bit, meaning this filter works from legend file, which in my case apparently had a modified version of Phase3 legend file with extra "filter" column. All makes sense now. Thank you.


text:
How do PCR duplicates arise and why is it important to remove them for NGS analysis?
I am trying to understand PCR duplicates in NGS analyses (actually whole-genome). I searched, and the best answer I found is in this blog.
However I don't understand if I understood how PCR duplicates arise  correctly because I cannot see the problem of having them in the downstream analysis - aside from computational problems, i.e. unnecessary redundancy.
If I understood correctly, PCR duplicates arise during library preparation when PCR amplifies the fragments with adapters. In this case, if you have duplicate fragments, you'll end up amplifying some of them twice or more.
However, when you are doing the mapping, they should fall on the same region, most likely decreasing the quality of the mapping (since they increase the consensus on a specific sequence, which could have been subject to sequencing errors). But aside from that, you should have the same mapping with respect to the mapping where duplicates are removed, although with decreased quality.
Is the quality problem the actual reason for removing PCR duplicates or is there anything that I'm missing?

In any scenario where depth of coverage is an important factor, PCR duplicates erroneously inflate the coverage and, if not removed, can give the illusion of high confidence when it is not really there.
For example, consider the following hypothetical scenario.
                                         *
TTTCATACTAACTAGCCTGCGGTCTGTGTTTCCCGACTTCTGAGTCATGGGGTTTCAATGCCTATAGATTC
               ..........................C.
                ............................
                      .............T..............
                       ..................C.........
                           ............................
                               ............................
                                   ............................
                                    .....C......................
                                       ..C.........................

The marked position has a coverage of 9 (9 overlapping reads), 4 of which suggest the presence of an alternate allele. Since these 4 reads map to different positions, they are independent observations supporting the same putative single nucleotide variant (SNV). The T in the third read can safely be ignored as a sequencing error, since it only appears once.
Now consider the following hypothetical scenario.
                                         *
TTTCATACTAACTAGCCTGCGGTCTGTGTTTCCCGACTTCTGAGTCATGGGGTTTCAATGCCTATAGATTC
                ............................
                      .............T..............
                       ..................C.........
                       ..................C.........
                       ..................C.........
                       ..................C.........
                               ............................
                               ............................
                               ............................

In this case, the marked position also has a coverage of 9, but 7 of those reads appear to be PCR duplicates. If we remove the duplicates, the marked position has only a coverage of 4, and only one C at the marked position, which is not enough confidence to call the SNV.

Best explanation I've seen so far

I agree with BCArg, very simple and clear explanation. Thank you.

PCR polymerases introduce errors. When an error arises in the first few cycles of amplifications, it will appear in a reasonably high fraction of DNA fragments in the library. After sequencing, you may see the same error occur to multiple reads. If you remove PCR duplicates when calling variants, all errors are reduced down to one read. For high-coverage data, you won't call the error as a false variant. However, if you keep PCR duplicates, the multiple occurrences of the error may make the error look like a real variant. This is particularly bad for calling SNVs from impure tumor samples.


text:
Java error when launching Pilon
I am trying to use Pilon to improve a reference based on some Illumina data I have got. So, I aligned the Illumina reads to the reference using bwa. Then, I want to use Pilon to improve the reference using the alignment file in BAM format.
Information about Pilon is found here. I installed Pilon using conda.
The output of 'pilon' is:
Pilon version 1.22 Wed Mar 15 16:38:30 2017 -0400

    Usage: pilon --genome genome.fasta [--frags frags.bam] [--jumps jumps.bam] [--unpaired unpaired.bam]
                 [...other options...]
           pilon --help for option details

The output of:
pilon --genome ref.fa --frags readsSorted.bam --output pilo --outdir output/scaffold

is:
 Exception in thread "main" java.lang.reflect.InvocationTargetException
    at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
    at sun.reflect.NativeMethodAccessorImpl.invoke(NativeMethodAccessorImpl.java:62)
    at sun.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43)
    at java.lang.reflect.Method.invoke(Method.java:498)
    at com.simontuffs.onejar.Boot.run(Boot.java:340)
    at com.simontuffs.onejar.Boot.main(Boot.java:166)
Caused by: java.lang.OutOfMemoryError: Java heap space
    at org.broadinstitute.pilon.GenomeRegion.<init>(GenomeRegion.scala:55)
    at org.broadinstitute.pilon.GenomeFile$$anonfun$contigRegions$1.apply(GenomeFile.scala:72)
 at org.broadinstitute.pilon.GenomeFile$$anonfun$contigRegions$1.apply(GenomeFile.scala:72)
    at scala.collection.TraversableLike$$anonfun$map$1.apply(TraversableLike.scala:234)
 at scala.collection.TraversableLike$$anonfun$map$1.apply(TraversableLike.scala:234)
    at scala.collection.immutable.Range.foreach(Range.scala:160)
    at scala.collection.TraversableLike$class.map(TraversableLike.scala:234)
 at scala.collection.AbstractTraversable.map(Traversable.scala:104)
 at org.broadinstitute.pilon.GenomeFile.contigRegions(GenomeFile.scala:72)
 at org.broadinstitute.pilon.GenomeFile$$anonfun$2.apply(GenomeFile.scala:52)
    at org.broadinstitute.pilon.GenomeFile$$anonfun$2.apply(GenomeFile.scala:52)
    at scala.collection.immutable.List.map(List.scala:277)
    at org.broadinstitute.pilon.GenomeFile.<init>(GenomeFile.scala:52)
    at org.broadinstitute.pilon.Pilon$.main(Pilon.scala:98)
    at org.broadinstitute.pilon.Pilon.main(Pilon.scala)
    ... 6 more

Note that the error is a memory issue:
java.lang.OutOfMemoryError: Java heap space

Like gringer already commented, the cause of your error is nicely given by the java program. Your default max heap settings are probably too low for the program. If you have more memory on your machine you can give max heap a higher value. See this post and answers on Stack overflow.
You'll probably need to use the -Xmxn argument in java.

"If pilon was installed by condo, edit the file "/Users/user/anaconda/bin/pilon" by changing the maximum memory limit from 1G to higher numbers like "14G". It works. "
check the thread: https://github.com/broadinstitute/pilon/issues/26


text:
TruSeq Exome Kit bed file old version
I have to run an analysis using Whole Exome Sequencing (WXS) data prepared with the TruSeq Exome Kit (2012).
I need the .bed file for the kit and I found it here.
However, this file was posted on Nov. 16, 2015, while the kit used in my case is from 2012. I searched but I didn't find the bedfile for that specific version online.
Will I run into problems if I use this version, and in this case, where can I retrieve the bedfile for 2012's version?

Care to explain why the downvote?

I haven't down voted you, but you don't show what have you done to understand what problems can you have when using an older file, nor the places you have searched for the bed files of 2012. Or maybe it was from biology.se, who knows ?

I usually search before writing a question, so also in this case I searched online to understand what could be the difference. Since online I couldn't find anything I wrote this question. Also, the fact that I put the reference to the Illumina kit page (the only page in which I could find reference to the bedfile of the kit) disproves that I didn't put the places I have searched. So I don't see the point of downvoting, since I: clearly stated the problem, provided details about what I've done, provided links to what I've found.

Well, those are imaginary internet points, don't mind them too much!

Actually I don't mind them :) but since I put some effort to search and write I couldn't understand why the downvote.

Usually it is more problematic to work with older versions of data with newer protocols than the other way around.
In this case, I think that the features annotated in the (newer) bed file should be more accurate than the older one. Unless you specifically want to compare your experiment with a previous experiment performed with a bed file from 2012 you'll be fine using the newer one.
However make sure that you are correctly using the right file and that you report (or record) each version of protocol and data used in case some doubts arise.


text:
TCGA bam file duplicates (aligned reads) and related last part of TCGA barcode
I am running an analysis on TCGA bam files (GCRh37/hg19) of whole exome sequencing data.
I usually find one bam file for each sample. However, at least one sample has two bam files, in which the last part of the name (before the extension) has the suffix .1 or .2 or .3...
For example, C500.TCGA-BL-A3JM-01A-12D-A21A-08.1.bam and C500.TCGA-BL-A3JM-01A-12D-A21A-08.2.bam.
I didn't find an explanation for this last suffix. Are they duplicates, or.. ? Is there some wiki about the last part of the filename?

If you check the bam list of TCGA, you'll see that they use different numbers for the bam files (maybe different settings). What the difference is, I don't know, but if you want the file that they used for their paper, it is C500.TCGA-BL-A3JM-01A-12D-A21A-08.2.bam, the other one is not in the list.


text:
How to align genomic sequence with corresponding amino acid sequence
Does anyone know of a program that can align a genomic sequence with introns with the corresponding amino acid sequence?
I have both the genomic sequence and the correct amino acid sequence but no information on the genemodel in e.g. gff or genbank format.
This is not what one would call a 'normal' alignment, but I would like it to look like this:
Hypothetical example:
nucleotide sequence with intron in italic:
ATGCACGATACGACTGACGTACGTACGTACGTACGTACGTACGTACGGAGACGTAGACTC
corresponding amino acid sequence:
MHDTTDVRTYVRRRRL
The intron does not encode for amino acid which should thus result in a gap in the 'alignment'.
Desired 'alignment':
ATGCACGATACGACTGACGTACGTACGTACGTACGTACGTACGTACGGAGACGTAGACTC
  M  H  D  T  T  D  V  R              T  Y  V  R  R  R  R L

Have a look at exonerate's protein2genome!
From the documentation:

Aligning a protein to genomic sequence:
Similarly, it is possible to align a protein sequence to the genome, (similar to GeneWise, but with
  heuristics).
exonerate --model protein2genome query.fasta target.fasta
This model will allow introns in the alignment, but also allow frameshifts, and
  exon phase changes when a codon is split by an intron.

The classic tools for this are genewise and the already mentioned exonerate. Both are capable fo aligning protein sequences to genomic DNA (or cDNA) while modelling splice sites to give you a gene model.
For example, this is the output1 of genewise on your sequences:
EMBOSS_001         1 MHDTTDVRTYVR----RRRL
                     MHDTTDVRTYVR    RRRL
                     MHDTTDVRTYVRTYVRRRRL
EMBOSS_001         1 acgaaggcatgcatgcacac
                     taaccatgcatgcatggggt
                     gctgtcatgcatgcagatac

The output also includes a GFF of your sequence:
EMBOSS_001  GeneWise    match   1   60  22.88   +   .   EMBOSS_001-genewise-prediction-1
EMBOSS_001  GeneWise    cds 1   60  0.00    +   0   EMBOSS_001-genewise-prediction-1

Back when I was using these tools daily to align proteins to genomes, I found that exonerate was much faster but sometimes slightly less accurate. This was a few years ago now though, so that might have changed. I would always suggest trying exonerate first and only going for genewise if it fails.

1This link will only be valid for a short while, but I just used the sequences from the OP as input for the GeneWise webserver.

Thanks, The GeneWise seems to work as well, but the Resulting GFF is quite useless. It would be great if it would give an annotation for the intron, but that is unfortunately not the case.

@user1817 what do you mean? Of course it annotates the intron. It won't in your example since that's too small to be a real intron, but try it on a real sequence and the GFF will include the introns. I just ran p54 against its gene, for example, and the introns are there.


text:
What's a good ontology for drug names?
I have ...

A database with patient phenotypes in, stored as HPO terms
Genetic data in whatever format I need

I want ...

To store drug names in a way that won't make my life difficult
If it's feasible, to make them so that you could draw links between HPO terms, medications, genetic changes, and outcomes.  Note that I'm not actively trying to do that, but I don't want to exclude the possibility of it in the future

Ideally, I'd like those drug names in some kind of ontology, and even more ideally an ontology that's available on http://bioportal.bioontology.org
So - what's the best drug ontology to use in this case?
Update

What relationships would this drug ontology describe? Structural similarity of the compounds? Affected pathways? How would they be linked? Do you want all drugs that can cause the same side effect linked, for instance? Or all drugs that can be used for the same disease?
  Drugs by disease would be the most useful.  The information that would be most useful would be something like

Atenolol

Atenolol 25mg
Atenolol 50mg
Atenolol 100mg

Propanolol

Propanolol 10mg
Propanolol 40mg
Propanolol 80mg
...

In short, being able to get a list of dosages for a given drug.  It'd be extra-useful if I could also retrieve a set of drugs for a given condition, so in the above example to retrieve a list of beta blockers (Atenolol, Metoprolol, Nebivolol, ...)

Have you looked at any of the various drug ontologies listed in bioportal.bioontology.org? Which ones have you discarded so we don't suggest the same ones? What were they missing that you need?

MDDB doesn't seem to be structured as above
Likewise The Drug Ontology
RxNorm seems to be the best available that I've seen

What relationships would this drug ontology describe? Structural similarity of the compounds? Affected pathways? How would they be linked? Do you want all drugs that can cause the same side effect linked, for instance? Or all drugs that can be used for the same disease? Have you looked at any of the various drug ontologies listed in http://bioportal.bioontology.org? Which ones have you discarded so we don't suggest the same ones? What were they missing that you need?

@terdon - added some more information to the question, hopefully that answers your questions.  Apologies for the slight vagueness in some of my answers, the reason is purely because I don't know - for me it's almost just a popularity contest as to which is the best for me!

The thing is, we don't really do such broad questions on the Stack Exchange network. There's even a specific close reason for them ("too broad") so this is not the right place to ask for a general discussion of the pros and cons. We might be able to help you with a specific case, but not on "what's best" which is just a question of preference and will depend on what you want to do. What's best for you might not be best for me. So please make it as specific as you can.

@terdon - changed the wording, is that better? :)

Since this question's been asked there's been a (what looks like a first useful?) release of Drug Ontology DRON. Perhaps worth checking again?


text:
problem of "ordering in physical positions" phasing SNPs with Shapeit
I'm trying to phase my data (whole genome resequencing SNPs) with Shapeit (from ped file) and I get this error message :
[33mERROR:[0m Wrong ordering in physical positions curr_pos=54 prev_pos=65

I found that it's due to the fact that my .map file is like this (4 columns):
Chromosome_number    SNP_ID   SNP_genetic_position    SNP_physical_position
Contig100131            .              0                        65
Contig100220            .              0                        54

And Shapeit only looks at SNP_physical_position but does not check if the SNPs are on the same contig or not and so it says it is not happy because position 54 is after 65!
How could I make it understand that the positions are not in the same contig and so there is no mistake? Must I run Shapeit per contig (I have a lot of them so I hope not)?

Why don't you just sort the index file? Will that work? It might still complain if you have multiple instances of the same position, I guess, but try: sed 's/  */\t/g' foo.map | sort -nk4 > new.map and then using the new.map file. If that works, let me know and I will post an answer.

On second thought, that shouldn't work. That would mean you would have the contigs not in order and if the program doesn't take the contig name into account, that would mess everything up.

The manual suggests that each "chromosome" needs to have its own input file.

You have to split the dataset by chromosomes. ShapeIT can just phase one chromosome at a time. If you are working on GWA data, you can proceed as follows:

for chr in $(seq 1 22) ; do plink --file example --chr $chr --recode --out example_chr$chr ; done

You will obtain 44 files; example_chr1.ped, example_chr1.map, ..., example_chr22.ped, example_chr22.map that you can phase separately.

ok, indeed it's written in big in the manual ! When you don't want to see something ....  Thanks a lot to both of you for your help !

@G.Dumont If this answer solved your issue, please take a moment and accept it by clicking on the check mark to the left. That will mark the question as answered and is the way thanks are expressed on the Stack Exchange sites.

what format does the file example need to be in? i.e. what order of columns? header? my map is in order `pos    chr     cM'


text:
Comparing two genome annotations
I have one-dimentional array (human genome). Also I have two annotations for it, we can think about them as different peaks (it's nucleosome and secondary structures). How can we find correlation and/or causation between this annotations? Currently we are forming hypothesis like "one element is always before other" and statistically testing them. Can we be smarter and use known tools or at least automate this?

Do you mean with something like a Markov model?

Yeah, but how will you implement them?

I am also interested in applying such models to multiple 'omics data, I hope this question will get some useful answers.

Could you specify which kind of hypothesis do you want to test? What I have seen is centering all the positions you are interested in and observe if there is a common pattern on the other variables observed (annotations, CpG, mutation rates...) which is not noise. See here an example

Perhaps phrase this as time series analysis? Does one event precede another in a statistical sense? I think if you look for tools in this domain, you'll find some good solutions.

See answers to this plot https://bioinformatics.stackexchange.com/questions/83/how-to-correlate-two-zero-inflated-bedgraph-like-signals?rq=1

You can express the annotations of the genome with a dot plot, x and y represent the two annotations along the genome. The you can summary the result with Spearman correlation coefficient(s) based on whole set of annotations or particular gene families.

I would start with metaplots for a visual/exploratory analysis to guide your statistical testing.  Cut the genome into short segments e.g. centred on the centre of each nucleosome peak, overlay all the segments on top of each other, and plot the abundance of secondary structure features at each base count as you step away from the centres.  This should give an easy view on whether there is any clear structure in the relative positions of the different objects.
Do the plots both ways round, and try different subsets of the objects (e.g. ones overlapping genes or not, different stringency thresholds for peak calling).  This should hopefully give you some insights on what is worth testing for.


text:
Use paste function with apply families
I am trying to assign different window sizes to my SNPs dataset to identify regions under selection.
this is the head of my data
head(snp_ids)
   snp_id   chr  pos
Chr01__912    1  912
Chr01__944    1  944
Chr01__1107   1 1107
Chr01__1118   1 1118
Chr01__1146   1 1146
Chr01__1160   1 1160

class(snp_ids)
data.frame

I have chosen 4 different window sizes, win_size <- c(15000, 30000, 50000, 100000).
I have assigned each of these different window sizes to my snp_ids dataset to identify how many SNPs are distributed within each window by looping through each window size
for (i in 1:length(win_size)){
  windows <- sapply(snp_ids$pos, function(x) (ceiling(x/win_size[i])))
}

I would like to help, but do not completely follow. I am not sure what you're trying to do. Furthermore, if you want help with code, you should provide a minimal data set to work with (for us to test it).

How are you getting 1.2 for window_30000?

@Anna To help others you can post your own answer to the question. But please do not include it in the question body (unless they are unsuccessful attempts)

Not clear what we are trying to do, but here is the solution using lapply that should replicate your existing forloop.
# example input data
snp_ids <- read.table(text = "
snp_id  chr  pos
Chr01__912   1  912
Chr01__944   1  944
Chr01__1107  1 1107
Chr01__1118  1 1118
Chr01__1146  1 1146
Chr01__1160  1 1160", header = TRUE)

# window sizes to loop through
win_size <- c(15000, 30000, 50000, 100000)

res <- cbind(snp_ids,
             data.frame(
               lapply(setNames(win_size, paste("window", win_size, sep = "_")), function(w)
                 as.numeric(paste(snp_ids$chr, ceiling(snp_ids$pos/w), sep = "."))
               )))

# result
res
#        snp_id chr  pos window_15000 window_30000 window_50000 window_1e.05
# 1  Chr01__912   1  912          1.1          1.1          1.1          1.1
# 2  Chr01__944   1  944          1.1          1.1          1.1          1.1
# 3 Chr01__1107   1 1107          1.1          1.1          1.1          1.1
# 4 Chr01__1118   1 1118          1.1          1.1          1.1          1.1
# 5 Chr01__1146   1 1146          1.1          1.1          1.1          1.1
# 6 Chr01__1160   1 1160          1.1          1.1          1.1          1.1

Edit: If we really want to use forloop, then try below:
for (i in win_size){
  windows <- ceiling(snp_ids$pos/i)
  snp_ids[, paste0("window_", i)] <- as.numeric(paste(snp_ids$chr, windows, sep = "."))
}

thanks so much @zx8754. I have edited my post a bit and hope it is more understandable. Your code worked for me. But, as I am learning R programming I wanted to know is there any other way that I could fix my own code above? Thanks so much for your help.

@Anna1364 See edit, yes we could use forloop but I prefer apply family.

This explicitly doesn't use apply() or lapply(), but rather the dplyr package. This is probably easier to read and debug:
library("dplyr")
snp_ids %>% mutate(window_15000=as.numeric(paste(chr, ceiling(pos/15000), sep=".")),
                   window_30000=as.numeric(paste(chr, ceiling(pos/30000), sep=".")),
                   window_50000=as.numeric(paste(chr, ceiling(pos/50000), sep=".")),
                   window_100000=as.numeric(paste(chr, ceiling(pos/100000), sep=".")))

This is a bit manual, don't you think?

Manual and clear is easier to debug and understand 6 months later (one could also put most of that in a function to make it shorter).

Of course, personal preferences. But if there were 100 windows, that change by 500 bases for each window in 6 month time.

Manual isn’t the problem — having duplication is. If I had to write this code, five times out of ten (at least!) I’d make a typo in at least one of the numbers, leading to mismatches between the column names and their actual values. That’s why we automate. (I’m saying “I” but what I really mean is “the average, competent programmer”; this class of mistakes is incredibly common).

A couple lines of code duplication is a nonissue. If someone neede to add a hundred columns then sure, but for three then trivial intelligibility is more useful.

authors: Bob Jarvis - Слава Україні


text:
Large range parameter identification
I have a model (ODE) of an enzymatic reaction, and I'd like to estimate four of its parameters to fit some available data. To obtain the model I am just using Mass-Action Kinetics.
Some more details:

The model is of the form $S_1+2S_2+S_3+E\leftrightarrow C \leftrightarrow P_1+P_2+P_3+E$, where $S_i$ are substrates concentrations, $P_i$ products concentrations, $E$ enzyme concentration, and $C$ the complex concentration.
$S_2$ is Hydrogen and $P_3$ is Ammonium ($NH_4$), so besides the aforementioned reaction I have one reaction for the buffer and one reaction for ammonium dissociation, however, the parameters involved are known.
Substrate concentrations are in the range of $1-10mM$.
Enzyme concentrations are in the range of $.1-10\mu M$.
The parameters I want to estimate are the 4 reaction rates, say, $k_1,k_{-1},k_2,k_{-2}$.
The available data (measurements) is only pH. Thus I have a simulation of the model where the final output is $pH=log_{10}(1/S_2)$.
Then, I try to estimate the parameters from the data and the (simulated) pH.

To estimate the parameters I am using Matlab, more specifically the 'Identification Toolbox' and the 'lsqnonlin' function.
I set the initial guess "visually", that is, I move manually the parameters and whenever I assume the simulation is "close enough", I set the initial parameters guess.
The problem is that some of the parameter of he initial guess are of order $1E20$. I am thinking that because of the large value of the parameters, the algorithm cannot find any better parameter values. I have tried moving several optimisation options without any valuable result.
Q: Is there some specific software, or method that is commonly used to estimate parameters of biochemical reactions? Is the large values of kinetic rates a common problem in parameter identification for these type of systems?
I do not know any $P_i$ though, just pH as a function of time. In principle nor E neither C are constant, but I could assume quasi-steady state conditions, no problem. My feeling is that I am not using or setting up correctly the numerical tool...

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


text:
bcftools output vs. bgzip
I've got what seems like it should be a simple question, but I can't seem to figure it out from google. bcftools has 4 output options:
Output compressed BCF (b), uncompressed BCF (u), compressed VCF (z), uncompressed VCF (v).

Are any of these the same as the output from bgzip?

Yes. I haven't found it documented anywhere, but I ran a simple test:
$ bcftools view -Oz --no-version foo.vcf > bar.vcf.gz
$ md5sum bar.vcf.gz
18a2f494710cf170ba79892936015ba8  bar.vcf.gz

$ gunzip bar.vcf.gz
$ bgzip bar.vcf
$ md5sum bar.vcf.gz
18a2f494710cf170ba79892936015ba8  bar.vcf.gz

$ gunzip bar.vcf.gz
$ gzip bar.vcf
$ md5sum bar.vcf.gz
9f5b1cf654966912b78a9cb7175ffc2c  bar.vcf.gz

I passed a vcf file through bcftools view -Oz, checked the md5sum of the resulting file, then uncompressed it and recompressed with bgzip and confirmed that the resulting file had the same md5sum. Finally, I also compressed using normal gzip to be sure that the md5sum would be different, as it was.
So yes, it does indeed look like bcftools creates files that are the same as though you'd passed them through bgzip.

Nice. Can't really argue with that.

Yes the compressed outputs should be the same format as bgzip.  bcftools also has an option to index in either the CSI or TBI formats.
The index function mentions that it indexes bgzip compressed files.  So one can assume that it also creates files using the bgzip compression.  Also all the pipelines that I have built use bgzip and bcftools.
[~]# bcftools index

About:   Index bgzip compressed VCF/BCF files for random access.
Usage:   bcftools index [options] <in.bcf>|<in.vcf.gz>

Thanks! So I guess bgzip genotypes.vcf = bcftools view -Oz genotypes.vcf, yeah?


text:
What are Approximate Read Counts (Library Sizes) and Lengths (Insert Sizes) for Next-Generation DNA Sequencers?
I am performing a simulation study and am curious about the parameters of my simulated metagenome.
What are the library and insert sizes of some of the most "used" sequencers. Mainly, I am interested in the different Illumina platforms (MiSeq, HiSeq, HiSeqX), 454 Pyrosequencing, and Ion Torrent.
I do not need the long read technologies like Pac-Bio's SMRT and Oxford Nanopore since they are not typically used with metagenomics, yet.

The library size depends usually on your budget. There are some reviews where they compare the different platforms, but these papers get outdated so quickly, that it would be better to just check wiki for the specs of different platforms. Also on the websites of Illumina or Life Tech, you'll find the specs.
NB. 454 sequencing is discontinued, so that is not very relevant anymore in my opinion.

Thanks for the Wiki link. This is what I was looking for!

Both of those depend on the application of the library.  For example whole genome assembly of a mammal would require a larger library size than targeted re-sequencing of a virus.  Also insert size varies if you are doing genome assembly a larger insert size is valuable to allow for long range linking and scaffolding.
I think that if you edit your question of be more specific people will be able to provided a more concrete answer other wise it is "it depends"


text:
Why total RNA-seq usually yields low mapping rate?
Maybe this is a silly question, but I'm really wondering why we usually get low mapping rates if we map total RNA-seq, but not poly(A)-enriched (in particular, for human, mouse, and zebrafish datasets)?
Doesn't the genome fasta file contain ribosomal RNAs (which are abundantly expected in the total RNA-seq libraries) as well?
I don't have the numbers for human right now, but I recall previously having low mapping rates for human data as well. What triggered me at the moment is actually when I mapped zebrafish data. I got around 60% mapping to the genome and 45% to the transcriptome with one dataset and 36% - 35% with another dataset (note that these are from 2 different studies and both are total RNA-seq, mapping was done using STAR and Salmon, which have their own methods).

It can also happen if you map to a transcriptome sequence instead of the genome sequence because there can be a lot of intron sequence from unprocessed mRNA and from unannotated non-coding RNAs

A likely explanation is that total RNA-Seq contains a high fraction of reads from ribosomal RNAs. Ribosomal RNAs are present in multiple copies across the genome, hence many reads map to multiple genomic locations and get discarded by the aligner. For example, STAR with default parameters considers a read as unmapped if it maps to more than 10 genomic loci (this behavior can be changed with the --outFilterMultimapNmax option).
To confirm if this is actually the case, you could check the number of multi-mapping reads in the aligner's log files.

I tried to increase the number of allowed multimapping reads and indeed the number of reads mapped to multiple loci has increased, albeit slightly (6.92% to 7.50%). Most of the non-mapped are categorised into "too short" (in STAR). Do you have any ideas what are these (well, it refers back to my initial question, since this happen in Salmon as well)?

Are you trimming adapters from the reads prior to mapping? That might explain the many "too short" reads. Also, have a look at this issue on STAR's Github repo: https://github.com/alexdobin/STAR/issues/169

@kaka01 If accounting for multi-mapping doesn’t solve your problem then there may simply be something wrong with your data: on high quality data sets, mapping total RNA to a genomic reference should typically yield >80% mapped reads. Many “too short” mappings may be indicative of RNA fragments. Did you perform size selection before sequencing?

Tom: No trimming was done prior to mapping. @Konrad The libraries weren't size selected. Your explanation makes sense, aligning the unmapped reads to the collection of rRNAs yielded high mapping rate (which is a sign that those are really rRNAs). Still not sure why STAR classified them as too short, maybe I just have to play around more with the parameters.

@kaka01 “Too short” here simply means that STAR is unable to align the read with high quality: either because the initial read (after trimming) is so short that it could match the reference virtually anywhere (< ~14 bases for something like the human genome) so we have low confidence in the correct origin. Or because STAR, when running with --alignEndsType Local (which is the default) is only able to match a small part of the read (which then leads back to the first problem).

@kaka01 Since you didn't do size selection, it's possible that you have a high number of short fragments, and you end up sequencing the adapter. You could try to trim the adapter and map the reads again. This might help to mitigate the issue because in STAR the minimum mapped length if a fraction of the total read length.

@KonradRudolph in you answer you say that mixed library of short and long RNA doesn’t seem to cycle well through the whole fragment length. I wasn't aware of this issue, do you know why this is the case? (sorry for commenting here, I don't have enough reputation to comment you answer directly)

@Tom Let’s talk offline (are you at the RNA club tonight? I’ll attend.) In short, I have no confirmation but it’s what we observe, and I don’t find it implausible given the chemistry.

In particular for mouse, this can happen if the ribo-depletion wasn't terribly efficient, since there's no Rn45s sequence in the reference genome. That, combined with numerous copies of tRNAs, 5S rRNAs, etc. causing issues with multimappers (see the answer from @Tom) can heavily decrease alignment rates.
Note that the human reference genome contains a few copies of the 45S (e.g., on GL000220.1 and chr21), so these reads will all multimap (possibly too many times). I don't know about the zebrafish reference genome.

When sequencing RNA from freshly harvested tissues under proper conditions, you should generally expect > 50% mapped reads. In fact, everything < 80% would usually raise concerns.
From your description (in the question and comments) it sounds as if your samples are potentially degraded and therefore saturated with short RNA fragments, either because you’ve got old tissue, you’re looking at environmental samples, or from improper handling.
And since you’re not performing any enrichment, the sequencing data will then be highly saturated with degraded RNA, too. The reason for this is that RNA-seq is a (hopefully somewhat uniform) sampling of your sample aliquot: everything in it is proportionally represented in the RNA-seq data. That’s why enrichment steps are important, to increase the proportion of whatever fraction we’re interested in.
In particular,

Total RNA is mostly composed of genes with many copies: mostly rRNA and (much less) tRNA. Hence the importance of ribosome depletion or poly(A) selection.
RNA in general is unstable. Depending on the sample collection (non-tissue sample origin? age of the tissue before harvesting?), a large fraction of the RNA may have degraded into minuscule fragments, which are further shortened by the library preparation. Size selection would get rid of these undesired fragments.

In my experience, the second issue is exacerbated on certain sequencers (I’ve seen it on an Illumina HiSeq 1500) because a mixed library of short and long RNA doesn’t seem to cycle well through the whole fragment length. So even if you have a (small but still present) fraction of long RNA fragments, the sequencer may not be able to synthesise their whole length efficiently. The effect is that there are proportionally more short reads than there were short fragments in the sample.
You can verify this by calculating the insert size distribution of your reads after adapter trimming (e.g. using picard, or by simply tallying read lengths1). Small reads (< ~14 nt) from degradation fragments are essentially unmappable, since they are too short to have a specific identity. Even if we allowed the mapper to align them their coordinates would be essentially random.

1 This works:
awk 'NR % 4 == 2 {c[length($0)]++} END {for (i in c) print i, c[i]}' in.fastq

You can't assume that the rRNA genes are in the genome. There many copies of them and not all are placed on chromosomes. We've had this issue with human and mouse data for total RNA preps. You need to ensure you map against the whole genome not just the chromosomes.

I 'm mapping to the whole genome indeed. But anyways, assuming that you are mapping against the whole genome, did you get significantly higher mapping rate?


text:
EMT/EMT-like processes in bioinformatics to study cancer progression
First I am not sure why, but it seems EMT can be used to infer the status of cancer progression.
What kind of databases are there for the EMT process in the context of cancer?
I saw this book which mentions ArrayExpress, but I need some help to understand what kind of data it stores.
Are there similar/better resources? How can I use EMT in the context of cancer progression study?

Are there other process like EMT that are being used to study cancer progression?
Are there EMT like events in leukemia progression?

Please note that the fields of Bioinformatics and Cancer Genomics move very quickly with new technologies. A source from 2010 would be regarded by many as out of date.

@TomKelly so you don't recommend on ArrayExpress?

Not unless you need to access the data from a particular study. Large scale genomics projects have far higher sample sizes and used newer technologies. What exactly are you trying to do? It does not seem like you are familiar with this field.

EMT describe cells' change in their state from being epithelial to the mesenchymal class. So, if a cancer cell line is gaining properties that allow it to move, it might reach the metastasis phase, were a cell can move to any point of the body and start a new tumor there.  As in the wikipedia page you linked:

EMT has also been shown to occur in wound healing, in organ fibrosis and in the initiation of metastasis in cancer progression.

I am not aware of any kind of database for the EMT process, the most relevant resource might be the gene ontology, which has a term for it, which lists 185 (human) genes as related to this biological process.
ArrayExpress is a database of transcriptomics, it :

stores data from high-throughput functional genomics experiments

You can evaluate the role of EMT in the data stored in the database.

What do you mean by evaluate the role of EMT in the data stored in the DB?

@0x90 Use one of the gene set enrichment techniques to the data, from doing a PCA with those genes to perform GSEA (like in the broad institute), or look for different expressed pathways like in SIPA or ...

thanks, if you have a good reference for similar work it would be much appreciated.

Similar work to what? To use analyze a pathway in a dataset?

yes to see expressed pathways in dataset which will be adjustable for my case.

This is a hard topic, the pathways are loosely described and there isn't any canonical pathway database, but GSVA is a good software tool that let you convert your matrix of genes x samples to pathways x sample, so you can use it to evaluate if a pathway is more expressed than another

Note that EMT is not the motility of cells. It is a change in their state from being epithelial to the mesenchymal class. In development, this is important in the formation of different tissues. In cancer, it is a disruption of regulatory processes, including cellular polarity and checkpoints on growth. Essential the tumour cells are becoming more like stem cells and less like the epithelial tissue. This is not "metastasis" that is another process by which cancer cells spread but EMT can contribute towards progression of disease towards that.

The epithelial-mesenchymal is a process by which cells lose some of the structures specific to their cell type. This is normal during embryonic development or wound healing for the formation of new cells. However, it is very uncommon in healthy adults and is thereby associated with the cancers in which it occurs. EMT allows cancer cells to become more like stem cells and could facilitate metastasis or immortal cell growth, key processes for cancer progression.
If you can detect the molecular components involved in EMT, you may be able to make some inferences about the progression of the patient’s disease. This process is specific to cancers in epithelial tissues (breast, stomach, and skin cancers for example). It is just one of vastly many pathways by which cancers can grow and develop. Many of these are specific to different cancers types by tissue or molecular subtype. Cancer research is therefore very challenging and requires a deep understanding of the particular system you are focusing on.
Gene expression data sets are a great resource for this. You can compare many pathways to clinical properties from previous studies. Microarray datasets such as ArrayExpress and Gene Expression Omnibus (GEO) were very popular for this but are now largely out of date. Many databases now support higher quality RNASeq datasets with larger sample sizes and fewer batch effects. The Cancer Genome Atlas (TCGA) and International Cancer Genome Consortium (ICGC) have among the largest Cancer databases. These are a tool for researchers to investigate any genes or pathways. They will not specifically address EMT, unless you focus on the genes involved in a relevant cancer type.

What do you mean by They will not specifically address EMT, unless you focus on the genes involved in a relevant cancer type. when referring to TCGA and ICGC. How can I find EMT related stuff there?

Exactly the same way as you would with ArrayExpress data. None of them are designed for this purpose: they all simply store expression data. It is only with 3rd party databases or knowledge of the biological context that you could extract this information from expression datasets of all genes.

Pathway analysis is a broad topic, far beyond the scope of an SE question. There are entire courses dedicated to it. Please be more precise in your questions if you wish to get meaningful answers. An example of how the data would be used would be valuable.