# Question
Title: How to access the patient-derived xenografts (PDXs) repository?
I saw this nature news item, it sounds that Cellminer is obsolete, is it right?

What are the new tools to analyze the "new cell lines"? Where's the PDX repository?

What about the European initiative that the article refers to?

I would be happy to see correlation between drugs and genes and similar capabilities as we can see on Cellminer.

# Answer
I already see for links when I try to search. Take a look at them below. Also if you can be precise about kind of data are you trying to find, it would be better

https://pdmr.cancer.gov/

http://cdt.northwestern.edu/news/patient-derived-xenograft-repository-now-available-researchers

http://www.epo-berlin.com/epo-tumor-models-xenografts.html

http://www.proxe.org/

http://data-analysis.charite.de/care/

http://www.oasis-genomics.org/

http://www.cbioportal.org/

Hope one of these above will serve the purpose for you.

Update: added some more links that might serve. As for the Cellminer the last update was 2012.
> 2 votes

---
Tags: cancer, cell-line
---


# Question
Title: How can I build a protein network pathway from a gene expression quantification file?
Assume I have found the top 0.01% most frequent genes from a gene expression file.

Let's say, these are 10 genes and I want to study the protein protein interactions, the protein network and pathway.

I thought to use string-db or interactome, but I am not sure what would be a plausible way to approach this problem and to build the protein network etc. Are there other more suitable databases?

How can I build a mathematical graph or network for these data?

# Answer
There are multiple ways to do this, and multiple protein interaction databases besides the ones you mentioned, such as BioGRID or IntAct. Interaction databases are different in how interactions are defined, sometimes it can be experimental evidence of interaction, sometimes coexpression, orthology-based predictions, etc.

There is no single solution to your problem. For String-DB you can use their R package STRINGdb.

```
# string-db
library(STRINGdb)

# retrieve full graph for H. sapiens
string_db <- STRINGdb$new(version="10", species=9606,
                          score_threshold=400, input_directory="" )

# define genes of interest
genes.of.interest <- ...

# get their neighbors
neighbors <- string_db$get_neighbors(genes.of.interest)

# get the subgraph of interest (your genes + their neighbors)
my.subgraph <- string_db$get_subnetwork(c(genes.of.interest, neighbors))

# look how many genes and interactions you have
my.subgraph

```

If you prefer GUI, you might consider using Cytoscape. They have multiple addons for interaction databases, e.g. for STRINGdb. There you can just provide your gene list and get a network in a few clicks.
> 9 votes

# Answer
Many interaction databases now work with PSI format files. Most of the main databases can do this and the EBI has set up PSICQUIC View, a very useful page where you can query multiple databases at once.

Note that it is very important to limit the results according to the detection method. There is a lot of noise in protein interaction databases. Depending on what you want to do you could limit to only experimentally verified interactions or to only direct, binary interactions (so exclude the results of, for example, ChIP analyses which can also find complexes) etc.

That said, here's a simple example script that will query the APID database using its PSICQUIC service:

```
#!/usr/bin/perl
use strict;
use warnings;
use LWP::Simple;

my @proteins;
## Read a list of target proteins, one per line (this expects UniProt names)
open(my $fh, "$ARGV[0]") or die "Need a list of proteins as the 1st argument: $!\n";
while (<$fh>) {
    chomp;
    push @proteins, $_;
}
close($fh);
## Get human interactions only
my $species="9606";

## Get the interactions for each target protein
foreach my $protein (@proteins) {
    my $queryUrl= "http://cicblade.dep.usal.es/psicquic-ws/webservices/current/search/query/$protein";
    $queryUrl  .= "?taxidA:$species%20ANDtaxidB:$species";

    my $tries=1;
    my $content = get $queryUrl;
    while ($tries<=10) {
        if (defined($content)) {
            $tries=11;
        } else {
            print STDERR "Could not retrieve $queryUrl, retrying($tries)...\n";
            $content = get $queryUrl;
        }
        $tries++;
    }

    # Now list all interactions
    my @lines = split(/\n/, $content);
    my $LINES= @lines;
    my $count = 0;
    for my $line (@lines) {
        $count++;
        my @flds = split(/\t/, $line); # split tab delimited lines
        # split fields of a PSIMITAB 2.5 line
        my ($idA, $idB, $altIdA, $altIdB, $aliasA, $aliasB, $detMethod, $author, $pub, $orgA, $orgB, $intType, $sourceDb, $intID, $conf) = @flds;
        ## Here you can add logic to limit the interactions by specific detection method codes ($detMethod)
        ## or database of origin etc.

        ## Print
        print "$line\n"
    }
}

```

I tested it by using a file with the UniProt ID for human TP53 (P04637) and it returns a list of 3988 interactions (showing 10 randomly selected results below):

```
$ foo.pl names.txt | shuf -n 10
uniprotkb:Q9BWC9    uniprotkb:P04637    -   -   uniprotkb:CCDC106(gene_name)    uniprotkb:TP53(gene_name)   psi-mi:"MI:0007"(anti tag coimmunoprecipitation)    Zhou, J. et al.(2010)   pubmed:20159018 taxid:9606(Homo sapiens)    taxid:9606(Homo sapiens)    -   psi-mi:"MI:0469"(intact)    intact:EBI-7812926  -
uniprotkb:P04637    uniprotkb:O15126    -   -   uniprotkb:TP53(gene_name)   uniprotkb:SCAMP1(gene_name) psi-mi:"MI:0018"(two hybrid)    Lim, J. et al.(2006)    pubmed:16713569 taxid:9606(Homo sapiens)    taxid:9606(Homo sapiens)    -   psi-mi:"MI:0463"(biogrid)   biogrid:720622  -
uniprotkb:P63165    uniprotkb:P04637    -   -   uniprotkb:SUMO1(gene_name)  uniprotkb:TP53(gene_name)   psi-mi:"MI:0018"(two hybrid)    Minty, A. et al.(2000)  pubmed:10961991 taxid:9606(Homo sapiens)    taxid:9606(Homo sapiens)    -   psi-mi:"MI:0463"(biogrid)   biogrid:262339  -
uniprotkb:P04637    uniprotkb:P31350    -   -   uniprotkb:TP53(gene_name)   uniprotkb:RRM2(gene_name)   psi-mi:"MI:0416"(fluorescence microscopy)   Xue, L. et al.(2003)    pubmed:12615712 taxid:9606(Homo sapiens)    taxid:9606(Homo sapiens)    -   psi-mi:"MI:0465"(dip)   dip:DIP-40167E  -
uniprotkb:P04637    uniprotkb:Q00987    -   -   uniprotkb:TP53(gene_name)   uniprotkb:MDM2(gene_name)   psi-mi:"MI:0004"(affinity chromatography technology)    Dai, MS. et al.(2004)   pubmed:15308643 taxid:9606(Homo sapiens)    taxid:9606(Homo sapiens)    -   psi-mi:"MI:0463"(biogrid)   biogrid:478073  -
uniprotkb:Q99576    uniprotkb:P04637    -   -   uniprotkb:TSC22D3(gene_name)    uniprotkb:TP53(gene_name)   psi-mi:"MI:0428"(imaging technique) Ayroldi, E. et al.(2015)    pubmed:25168242 taxid:9606(Homo sapiens)    taxid:9606(Homo sapiens)    -   psi-mi:"MI:0463"(biogrid)   biogrid:1255896 -
uniprotkb:P04637    uniprotkb:Q00987    -   -   uniprotkb:TP53(gene_name)   uniprotkb:MDM2(gene_name)   psi-mi:"MI:0415"(enzymatic study)   Lui, K. et al.(2013)    pubmed:23572512 taxid:9606(Homo sapiens)    taxid:9606(Homo sapiens)    -   psi-mi:"MI:0463"(biogrid)   biogrid:859223  -
uniprotkb:P25685    uniprotkb:P04637    -   -   uniprotkb:DNAJB1(gene_name) uniprotkb:TP53(gene_name)   psi-mi:"MI:0004"(affinity chromatography technology)    Qi, M. et al.(2014) pubmed:24361594 taxid:9606(Homo sapiens)    taxid:9606(Homo sapiens)    -   psi-mi:"MI:0463"(biogrid)   biogrid:938952  -
uniprotkb:P04637    uniprotkb:Q8IW41    -   -   uniprotkb:TP53(gene_name)   uniprotkb:MAPKAPK5(gene_name)   psi-mi:"MI:0424"(protein kinase assay)  Sun, P. et al.(2007)    pubmed:17254968 taxid:9606(Homo sapiens)    taxid:9606(Homo sapiens)    -   psi-mi:"MI:0469"(intact)    intact:EBI-1202077  -
uniprotkb:Q13526    uniprotkb:P04637    -   -   uniprotkb:PIN1(gene_name)   uniprotkb:TP53(gene_name)   psi-mi:"MI:0096"(pull down) Mantovani, F. et al.(2007)  pubmed:17906639 taxid:9606(Homo sapiens)    taxid:9606(Homo sapiens)    -   psi-mi:"MI:0469"(intact)    intact:EBI-6112688  -

```
> 6 votes

# Answer
If you run this Rscript using the gene name as an argument, you'll get a file with the pathway written to the directory.

```
#!/usr/bin/Rscript
args = commandArgs(trailingOnly=TRUE)
library(paxtoolsr)
id <- args[1]
write.table(graphPc(source=id,kind='neighborhood',
            format='BINARY_SIF',verbose=TRUE),
            file=paste(id,'pathway',sep='_'), sep='\t',quote=FALSE)

```
> 2 votes

# Answer
Why not just use string-db's online tool? There you can adjust all parameters of interest such as interaction confidence as well as pull in extra protein interactions.

If you make a user you can also upload a custom background gene-set set which means the analysis of go-terms and interconnectedness that sting provides is usable.
> 1 votes

---
Tags: proteins, networks, gene, pathway
---


# Question
Title: 10x Genomics Chromium single-cell RNA-seq data analysis options?
Provide an overview of 10x data analysis packages.

10x provides Cell Ranger which prepares a count matrix from the bcl sequencer output files and other files (see bottom of page https://support.10xgenomics.com/docs/license for the programs it uses).

What can we do with the output files?

# Answer
# Data preparation

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

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

---

# Data analysis

**Loupe Cell Browser visualization of** `cloupe.cloupe` **files**

**Count table** `matrix.mtx` **analysis options:**

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

---
Tags: rna-seq, r, bioconductor, scrnaseq, 10x-genomics
---


# Question
Title: Searching for gene expression data by cell line
I have two cancer cell lines (OCI-Ly18 & riva) that I want to find gene expression data for, but I'm not aware of many gene expression databases that allow searching by cell-line without searching by gene.

I tried Genevestigator on the recommendation of this thread but found no data on either cell line.

What would be the best way to find gene expression data for these cell lines?

# Answer
Try the Gene Expression Omnibus \- it looks like they have some datasets.
> 5 votes

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

---
Tags: database, cancer, differential-expression, cell-line
---


# Question
Title: Is there a software or web server to calculate thermodynamic parameter enthalpy change of cellulase?
I want to calculate thermodynamic paremeter enthalpy change of cellulase, but I don't know what software can do this work. Please tell me, thank you very much !

# Answer
You could use PISA assuming you have a PDB file of cellulase? Sorry it's a short answer and should be a comment, I don't have the reputation, however.
> 2 votes

---
Tags: protein-structure, thermodynamics
---


# Question
Title: Dataset: Locations of regulatory sequences in the human genome?
I am looking for the positions of annotated regulatory sequences (promoters, enhancers and suppressors) in the human genome.

I looked at Ensembl regulatory Build and PAZAR but I am not used to look for datasets and I failed to find what I was looking for.

**Can you help to find annotated regulatory sequences (start and end positions) in the human genome?**

# Answer
> 5 votes
This can be done quite easily using Ensebl's BioMart.

1. Choose the Ensembl Regulation database:
2. Select the "Human Regulatory Features" dataset:
3. That's basically it right there, just click on "Results":
4. Export to file and click "Go":

This will download a file called `mart_export.txt` which looks like this (I chose TSV for tab separated values):

```
$ head mart_export.txt
Chromosome/scaffold name    Start (bp)  End (bp)    Feature type
18  76429380    76430144    Open chromatin
8   66405962    66406502    Open chromatin
4   61184401    61184600    CTCF Binding Site
X   40733600    40737000    Promoter
5   97407001    97407200    CTCF Binding Site
X   73946201    73946600    Promoter Flanking Region
15  19948201    19949200    CTCF Binding Site
5   11302601    11303143    Open chromatin
2   208407801   208408000   CTCF Binding Site

```

As you can see the format is pretty straightforward. The first column has the chromosome or scaffold name (it should always be a chromosome for human), the second and third are the star and end positions and the last field is the type of region.

You can also choose to limit by region type by choosing a filter (click the "Filters" link) in BioMart before downloading or by simply parsing the file once you've downloaded it.

# Answer
> 0 votes
Promoters and enhancers based on CAGE-seq from FANTOM5: http://fantom.gsc.riken.jp/5/data/

Just choose the category and download bed-file, e.g. CAGE-peaks would be promoters. There are READMEs in each folder that explain the data.

---
Tags: database, human-genome
---


# Question
Title: GRCh38 vcf file with common cancer mutations
Is there a vcf file on the GRCh38 assembly with common cancer mutations I can download somewhere? Maybe from one of the big international cancer genomics consortia?

By common, I mean whichever mutations have been found recurrent in different types of cancer.

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

1. Get the VCF file describing all variants in Clinvar from NCBI:

```
wget ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/clinvar.vcf.gz

```
2. Extract any variants whose VCF line contains the word "cancer":

```
zgrep -iE '^#|CLNDBN=[^;]*cancer' clinvar.vcf.gz > cancer.vcf

```

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

# Answer
If you are looking for cancer mutations, the primary resource is COSMIC and they provide GRCh38 VCFs. The download page is here: http://cancer.sanger.ac.uk/cosmic/download

It'll be up to you how you define "common", but the VCF includes a lot of information you can use.
> 8 votes

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

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

---
Tags: vcf, cancer
---


# Question
Title: How exactly is "effective length" used in FPKM calculated?
According to this famous blog post, the effective transcript length is:

$\tilde{l}\_i = l\_i - \mu$

where $l\_i$ is the length of transcript and $\mu$ is the average fragment length. However, typically fragment length is about 300bp. What if when the transcript $l\_i$ is smaller than 300? How do you compute the effective length in this case?

A related question: when computing the FPKM of a gene, how to choose a transcript? Do we choose a "canonical" transcript (how?) or combine the signals from all transcripts to a gene-level FPKM?

# Answer
> 12 votes
I have a blog post that describes the effective length (as well as these different relative abundance units). The short explanation is that what people refer to as the "effective length" is actually the *expected* effective length (i.e., the expectation, in a statistical sense, of the effective length). The notion of effective length is actually a property of a transcript, fragment pair, and is equal to the number of potential starting locations for a fragment of this length on the given transcript. If you take the average, over all fragments mapping to a transcript (potentially weighted by the conditional probability of this mapping), this quantity is the expected effective length of the transcript. This is often approximated as simply $l\_i - \mu$, or $l\_i - \mu\_{l\_i}$ --- where $\mu\_{l\_i}$ is the mean of the *conditional* fragment length distribution (conditioned on the fragment length being \< $l\_i$ to account for exactly the issue that you raise).

# Answer
> 15 votes
The effective length is $\tilde{l}\_i = l\_i - \mu + 1$ (note the R code at the bottom of Harold's blog post), which in the case of $\mu \< l\_i$ should be 1. Ideally, you'd use the mean fragment length mapped to the particular feature, rather than a global $\mu$, but that's a lot more work for likely 0 benefit.

Regarding choosing a particular transcript, ideally one would use a method like salmon or kallisto (or RSEM if you have time to kill). Otherwise, your options are (A) choose the major isoform (if it's known in your tissue and condition) or (B) use a "union gene model" (sum the non-redundant exon lengths) or (C) take the median transcript length. None of those three options make much of a difference if you're comparing between samples, though they're all inferior to a salmon/kallisto/etc. metric.

Why are salmon et al. better methods? They don't use arbitrary metrics that will be the same across samples to determine the feature length. Instead, they use expectation maximization (or similarish, since at least salmon doesn't actually use EM) to quantify individual isoform usage. The effective gene length in a sample is then the average of the transcript lengths after weighting for their relative expression (yes, one should remove $\mu$ in there). This can then vary between samples, which is quite useful if you have isoform switching between samples/groups in such a way that methods A-C above would miss (think of cases where the switch is to a smaller transcript with higher coverage over it...resulting in the coverage/length in methods A-C to be tamped down).

# Answer
> 4 votes
For the effective length part please see to Devons answer. I just have a small addition: Kallisto/Salmon/RSEM incorporate all bias estimates into the effective length meaning the effective length not only represent the length bias if you take the values from those tools (given that they were run with the bias algorithms enabled naturally).

With regards to getting gene level estimates you should not choose a specific transcript. Instead you should extract/calculate the RPKM/FPKM/TxPM (transcript per million that Kallisto/Salmon/RSEM outputs) for each transcript and sum them up to get the gene level estimate.

---
Tags: rna-seq
---


# Question
Title: How to estimate whether a long-read is meaningful sequence?
**The setup**

Imagine that I work on an organism without a reference genome, and that the closest reference genome I can get is quite diverged. E.g. ~10% diverged in terms of SNVs when measured with short reads, and also has a lot of structural variants too.

Now imagine I get a 1 million base-pair long-read (e.g. from Nanopore data) for my organism. The question is this:

How can I estimate the proportion of the read that is meaningful sequence vs garbage?

**Some things that probably won't work**

Most standard approaches won't work here. E.g. I could try mapping the read to the reference, but *even if the read was perfectly good* I wouldn't expect most of it to map thanks to true structural variations between the read and the reference. The same goes for standard alignment or BLAST.

**Some things that might work**

The best naive method here seems to be to cut the read up into smaller pieces (either overlapping or not) and use standard approaches to map/align each of these.

So, what have people tried for this? And what tools have you used and why?

# Answer
As a first pass, you could check if the read is chimeric. Porechop searches for known nanopore adaptors both on the ends and through the middle of the read. This won't resolve issues around blocked or empty pores, but it will at least check if you have found two long-ish reads lumped into the same file.

By default, Porechop splits chimeric reads into two (or more, I suppose!) non-chimeric reads, but the option `--discard-middle` would be a quick and easy way to check - run it on a fasta file containing only your long read, and if the output is empty, the read is chimeric.
> 4 votes

# Answer
If you're looking at a single organism, in the absence of a reference genome you can map other reads to suspect reads and look at coverage. Looking at the actual sequence can also be useful: real DNA usually doesn't have an abundance of two different bases.

Nanopore reads also give another way to see if the read looks weird by having a look at the raw signal. Our chimeric reads paper gives a few examples of what DNA sequence should look like under normal circumstances. Here's the first raw signal figure from that paper:

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

---
Tags: sequence-alignment, nanopore, read-mapping
---


title: Finding homologs of a protein sequence

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

# Answer
> 5 votes
It depends on what you're looking for. If you're just looking for sequence homology, then you can simply pick the best hits from a blast search. If, however, you are referring to *functional* homology, if you are looking for the protein which has the same functions as your query, then it's more complicated.

Sequence homology is not enough to infer functional homology. For example, you can have cases of gene duplication and subsequent functional divergence. Such *paralogs* are still *homologs* (paralogs are a subset of homologs), but they don't necessarily have the same function. It is also often the case that the homolog (be it orthologous or paralogous) of a protein in species B has a completely different function than its homolog in species A despite a high level of sequence similarity. This is usually very hard to determine *in silico*.

To find the functionally homologous protein(s), you would ideally need to identify the essential residues that allow your protein to perform its function. This could be done using something like PFam which will identify protein domains. You can then check whether the homologs you find also have this domain.

This is essentially what PSI-blast does. Although it doesn't take domains into account, each successive iteration is used to build a model of your proptein. The model is built under the reasonable assumption that highly conserved residues are important. So it will consider more diverged sequences as homologous if those residues are conserved.

If you know how your protein works and what residues are important, you can use that knowledge to refine the results of your PSI-blast. If you don't, you'll have to use only "good" hits to make the model. One way to do this, for well studied proteins, is to only add proteins that are already annotated as homologs of what you are looking for to build your model, then use that model to search in un-annotated species.

If you don't know, you could try looking for recognizable protein domains in your query protein (use PFam) and then use the HMM (hidden markov model) of the domain to identify important residues. For example, this is the HMM logo for the zf-A20 zinc finger domain:

The huge cysteine (C) residues are shown that size because they are very conserved across proteins carrying this domain and, presumably, are functionally important for the domain. So, if you pass your protein through PFam and identify domains, find the important residues and make sure all your homologs have those conserved. If using PSI-blast, only include sequences where those residues are conserved in the results you keep.

Finally, another useful tool that works in the same way is HMMER. This takes a protein alignment as input, like PSI-blast builds an HMM model from it and then can use this model to query a protein database for more hits. Methods like HMMER and PSI-blast are far better than simple sequence similarity approaches when looking for homologs.

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

> SmartBLAST is a new and experimental NCBI tool that makes it easier to complete common sequence analysis tasks, such as finding a candidate protein name for a sequence, locating regions of high sequence conservation, or identifying regions covered by database sequences but missing from the query.
>
> To do this, SmartBLAST performs the following tasks in much less time than it takes to run a typical BLASTp search:
>
> * a BLASTp comparison of the query with the closest matching sequences available;
> * a parallel BLASTp search to find the closest matches to high quality sequences from model organisms;
> * a multiple alignment between the query and five of the closest matching sequences (usually including two high quality sequences);
> * an analysis that produces a phylogenetic tree from the multiple sequence alignment.

\[from NCBI Insights\]

---
Tags: sequence-homology
---


title: How to extract fasta from a blastdb

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

# Answer
> 9 votes
You can extract fasta sequence from a blastdb constructed from a fasta file using blastdbcmd which should be installed when you install blast/makeblastdb.

```
blastdbcmd -entry all -db <database label> -out <outfile>

```

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

* `my_database.nhr`
* `my_database.nsq`
* `my_database.nin`

and you wanted your fasta output file to be called `reference.fasta` you would run the following:

```
blastdbcmd -entry all -db my_database -out reference.fasta

```

---
Tags: fasta, blast
---


# Question
Title: Getting a "system is computationally singular" error in sleuth
I am analysing 142 samples belonging to 6 batches. Additionally, those samples belong to 72 strains, which means that for most of the strains there are two samples.

I could fit simple models (for strain and batches for instance), but when I get to the "full" model (~batch+strain), I get the following error:

```
so <- sleuth_fit(so, ~strain+batch, 'full')
Error in solve.default(t(X) %*% X) :
  system is computationally singular: reciprocal condition number = 5.2412e-19

```

I should point out that of the 72 strains, only 15 have samples in distinct batches. This means that most strains (57) have both samples in the same batch.

Is the error due to an unknown bug or rather to the experimental design? Does it mean that the information on batches cannot be used?

Thanks

EDIT I've posted the experimental design in a gist

```
batch   strain  replica
batch_1     strain_41   1
batch_4     strain_41   2
batch_1     strain_28   1
batch_4     strain_28   2
batch_1     strain_26   1
[...]

```

# Answer
You should be able to remove any one of the following strains to end up with a rank-sufficient model matrix: 5, 10, 12, 13, 14, 15, 19, 26, 28, 3, 30, 32, 36, 39, 41, 45, 46, 49, 5, 50, 52, 53, 58, 59, 60, 69, 8. As an aside you can figure this sort of thing out as follows (I read your dataframe in a the `d` object):

```
> m = model.matrix(~batch+strain, d)
> dim(m) # 142 row, 77 columns, so minimum rank is 77
> qr(m)$rank # 76, so just barely rank insufficient
> #see if we can remove a single column and still get rank 76
> colnames(m)[which(sapply(1:77, function(x) qr(m[,-x])$rank) == 76)]

```

You obviously don't want to remove the batch columns or the intercept. The normal tricks that you can sometimes use to get around this issue with case-control studies don't appear to help here, which is why I would just drop a strain and call it done. Keep in mind that your power is still likely terrible. I generally recommend at least 6 replicates per group (scale down the number of groups to fit your budget).

EDIT:

Once the desired strain is removed from the model, it can be fit directly into the `sleuth_fit` function to obtain the full model:

```
> m = m[, -9] # whatever column to drop to get the appropriate rank
> so <- sleuth_fit(so, m, 'full')

```
> 2 votes

# Answer
This happen when the variables (strain +batch) create a design matrix like this:

```
batch strain
1 1 #
1 1 #
1 2
2 2
3 3
4 3
...
16 72

```

Which means that some of the covariates are not linearly independent (ie batch 1 and strain 1), all the strain 1 is in batch 1.

You can correct for batch effects, but not in this design of the linear model (if you want to take into account the strain). You could do one batch more with those 15 strains that are in a single batch (if they are in different batch between them) that way you would get an independent design.

Many strains in a single batch are in the same batch. You need to increase the number of samples (recommended anyway due to the low number of samples per strain) to avoid this problem.

---

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

---
Tags: rna-seq, r, differential-expression
---


# Question
Title: Visualisation of long read RNA-Seq splicing
I have a dataset of Oxford Nanopore cDNA reads. Many of my reads are full-length or close to full-length transcripts, and I and am interested in examining alternative splicing. For this, I would like to begin by visualising my reads and comparing variants qualitatively.

I have tried visualising my data in both IGV and SeqMonk, but neither has given a satisfying result.

**IGV**

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

**SeqMonk**

SeqMonk gives a much cleaner visualisation, but unfortunately shows no links between aligned exons, and worse, I cannot find the read IDs. This means there is functionally (as far as I can see) no way to tell which exons come from the same read.

**EDIT: Sashimi Plot**

As recommended in this answer, I have tried using IGV's Sashimi Plot. Unfortunately, I believe the high error rate and relatively low coverage is causing this to create somewhat confusing output. The Sashimi Plot shows mostly nearly every junction to have just a single read supporting it.

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

# Answer
IGB gives the intended result without much hassle. It doesn't have the bells and whistles of IGV, but presents a clean and intuitive view of the individual reads.

In IGB, each read is presented as a single, contiguous line.

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

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

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

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

# Answer
DNASTAR's software is for purchase, but high quality. GenVision Pro does genomic visualization, including Sashimi plots.

Edit: not sure why this answer is being downvoted, unless it's because the software isn't free. OP has tried IGV and SeqMonk, I mentioned an alternative he might not have heard of.

Here is a video demonstrating the use of Sashimi plots in GenVision Pro:

http://www.dnastar.com/t-support-videos.aspx?video=YJvcERoSIsg
> 5 votes

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

For example:
> 3 votes

---
Tags: rna-seq, nanopore, read-mapping, visualization
---


# Question
Title: What is the distribution of indel sizes in a healthy human genome? of insertion:deletion ratios?
My understanding is that indels are from 1bp to 10Kb, and a healthy genome has ~400K-500K Indels. Surely most of these are small.

What is the distribution of insertion sizes in a healthy human genome? What is the distribution of deletion sizes? What is an average ratio of insertions:deletions?

I have not found a definitive reference.

# Answer
> 6 votes
One of the 2015 papers from the 1000 genomes project has a nice figure (figure 1) showing the size distribution of medium to large sized insertions and deletions:

From another 2015 1000 genomes paper, one can see that the absolute number of smaller indels is much larger, though an exact size range isn't given (as far as I saw). If you really want to know that, just download the most recent 1000 genomes VCF file(s) and compute the exact median size and/or distribution that you want.

# Answer
> 5 votes
Genome-In-A-Bottle (GIAB; version 3.3.2) contains 3.21M SNPs on auto+X chromosomes and 0.51M INDELs in 2.58Gb confident regions. The ins:del ratio is 0.92. On the CHM1-CHM13 pacbio assembly (European ancestry), there are 3.57M auto+X SNPs and 0.58M INDELs in 2.71Gb confident regions \[reference\] with an ins:del ratio 0.99. These give you an idea of relative short INDELs.

For indels longer than 50bp, I would highly recommend to read this paper by Chaisson et al. This call set is constructed from the whole-genome assembly of the CHM1 genome. It is far more comprehensive and probably more accurate than all the other call sets. The following is Table 1 from the paper:

I have to say that the high ins-to-del ratio is surprising to me, but this is what the data tells us. I have done a similar analysis and reached a similar ratio.

PS: These numbers are essentially from one haplotype. They will be higher for a diploid sample.

EDIT: some numbers in the original post are wrong.

---
Tags: variation, indel, literature-search
---


# Question
Title: How does Cellminer's "Cross-correlations of transcripts, drugs, and microRNAs" work
I have seen one of Cellminer tools. I am not sure how do they calculate the cross correlation of the genes, what does it actually mean? Based on what databases?

For example if I take their example input (HUGO format):

```
abcb1
BRCA2
CNBP

```

I get the following cross correlation matrix (based on what signals).

```
Identifier      abcb1   BRCA2   CNBP
abcb1           1      -0.142   0.069
BRCA2          -0.142   1       0.176
CNBP           0.069    0.176   1

```

# Answer
> 2 votes
In the article describing Cellminer says about the cross-correlations:

> ... a cross-correlation table of the resultant z-scores can be generated ...

On page 3503

And the z-scores are defined previously as:

> The tool output includes relative transcript intensity presented as z-scores...

---
Tags: statistics, gene
---


title: Duplicate long hits from PSI-BLAST

text:
# Question
Title: Duplicate long hits from PSI-BLAST
I had a protein Refseq ID and I PSI-BLASTed this sequence against Refseq database. We all know that the Refseq is a Reference sequence database and it shouldn't have redundancy. After BLASTing my sequence, at first iteration I got 1000 hits and among them there were a lot of redundant sequences!

My sequence had 241 amino acids and I found a lot of sequences with 100% identity, 100% cover and 0 E-value exactly the same as my sequence but with different IDs. All of these IDs were from RefSeq! In other iterations and after adjusting format options, I got this redundancy with other sequences from other species. My sequence is related to a chain of a multichain protein (E.coli fumarate reductase)

For example, when using NP\_418578 as a query, I found WP\_078165098.1, WP\_064226696.1, WP\_062863447.1, WP\_001401474.1 and other that were identical.

I want to know what is wrong with Refseq. Is it really a Reference sequence database? If it is, where does this redundancy come from and why?

# Answer
This is what is known as a feature, not a bug. Note that your identical proteins all have accessions starting with `WP_`. These are special "non-redundant proteins". Many sequences—particularly bacterial sequences—are identical between various different species so having a separate RefSeq entry for each of them would be inefficient. Therefore, RefSeq combines multiple such proteins into a single `WP_` record. This is documented here (emphasis mine):

> A new type of RefSeq protein record which represents non-redundant protein sequences was introduced in mid-2013. This record type was introduced to address a growing issue with redundancy in the Prokaryotic RefSeq protein dataset that coincided with a significant increase in bacterial genome submissions from individual isolates and closely related bacterial strains. For example, a large number of high-quality bacterial genomes may be submitted during a disease outbreak. The submitted sequences may reflect pathogen evolution during the course of the outbreak but the majority of the encoded proteins from these genomes may be identical to each other. As RefSeq includes these genomes, per community requests, this resulted in increased redundancy. **By representing identical proteins using a single non-redundant protein accession number (with the prefix 'WP\_'), redundancy in the database is significantly reduced.**
>
> \[ . . . \]
>
> Because a non-redundant protein sequence may be found in RefSeq genomes from multiple species, **the organism information provided on the protein record reflects the lowest-common taxonomic node ranging from the genus species level to super-kingdom**. A non-redundant protein record that provides organism information at the level of a genus, family, or even super-kingdom does not mean that the protein is found in all RefSeq genomes below that taxonomic classification. It only indicates that the protein is found in more than one genome of different species for which the genus, family, or super-kingdom classification is the lowest common taxonomic node.

So, your query was NP\_418578.1, anaerobic fumarate reductase catalytic and NAD/flavoprotein subunit from *E. coli* strain K-12, substrain MG1655. The first thing to notice is how specific that is. This is the protein found from one specific substrain of one specific strain of one specific bacterial species. It is reasonable to expect that there will be identical sequences from many, many closely related species. Both from, most probably, all other strains and substrains of *E. coli* and from other, related bacteria.

Now, the specific sequences you mention are actually *slightly* different and not 100% identical. Below is a multiple alignment of NP\_418578.1 and the 4 WP\_ sequences you mentioned. Note that each of the 5 entries is slightly different. Each has one residue that differs from the rest. Look for the `:` in the identity line, there are 4 `:` and all others are `*` (I am only showing the relevant alignment blocks here and have removed those where all 4 sequences were identical):

```
WP_001401474.1      MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAVAQDH
WP_062863447.1      MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAVAQDH
WP_064226696.1      MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAVAQDH
NP_418578.1         MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAVAQDH
WP_078165098.1      MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAIAQDH
                    *******************************************************:****

[ . . . ]

WP_001401474.1      KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHVRGLVAMNMMEGT
WP_062863447.1      KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHVRGLVAMNMMEGT
WP_064226696.1      KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHIRGLVAMNMMEGT
NP_418578.1         KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHVRGLVAMNMMEGT
WP_078165098.1      KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHVRGLVAMNMMEGT
                    ***********************************************:************

[ . . . ]

WP_001401474.1      GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKVSQAFWHEWR
WP_062863447.1      GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKISQAFWHEWR
WP_064226696.1      GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKVSQAFWHEWR
NP_418578.1         GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKVSQAFWHEWR
WP_078165098.1      GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKVSQAFWHEWR
                    **************************************************:*********

WP_001401474.1      KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGIDPVKEPIPVRPTAHYTMG
WP_062863447.1      KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGVDPVKEPIPVRPTAHYTMG
WP_064226696.1      KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGVDPVKEPIPVRPTAHYTMG
NP_418578.1         KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGVDPVKEPIPVRPTAHYTMG
WP_078165098.1      KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGVDPVKEPIPVRPTAHYTMG
                    *****************************************:******************

[ . . . ]

```

Your sequence (NP\_418578.1) is only identical to one `WP_*` multi-species sequence, WP\_001192973:

```
WP_001192973.1      MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAVAQDH
NP_418578.1         MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAVAQDH
                    ************************************************************

WP_001192973.1      DSFEYHFHDTVAGGDWLCEQDVVDYFVHHCPTEMTQLELWGCPWSRRPDGSVNVRRFGGM
NP_418578.1         DSFEYHFHDTVAGGDWLCEQDVVDYFVHHCPTEMTQLELWGCPWSRRPDGSVNVRRFGGM
                    ************************************************************

WP_001192973.1      KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHVRGLVAMNMMEGT
NP_418578.1         KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHVRGLVAMNMMEGT
                    ************************************************************

WP_001192973.1      LVQIRANAVVMATGGAGRVYRYNTNGGIVTGDGMGMALSHGVPLRDMEFVQYHPTGLPGS
NP_418578.1         LVQIRANAVVMATGGAGRVYRYNTNGGIVTGDGMGMALSHGVPLRDMEFVQYHPTGLPGS
                    ************************************************************

WP_001192973.1      GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKVSQAFWHEWR
NP_418578.1         GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKVSQAFWHEWR
                    ************************************************************

WP_001192973.1      KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGVDPVKEPIPVRPTAHYTMG
NP_418578.1         KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGVDPVKEPIPVRPTAHYTMG
                    ************************************************************

WP_001192973.1      GIETDQNCETRIKGLFAVGECSSVGLHGANRLGSNSLAELVVFGRLAGEQATERAATAGN
NP_418578.1         GIETDQNCETRIKGLFAVGECSSVGLHGANRLGSNSLAELVVFGRLAGEQATERAATAGN
                    ************************************************************

WP_001192973.1      GNEAAIEAQAAGVEQRLKDLVNQDGGENWAKIRDEMGLAMEEGCGIYRTPELMQKTIDKL
NP_418578.1         GNEAAIEAQAAGVEQRLKDLVNQDGGENWAKIRDEMGLAMEEGCGIYRTPELMQKTIDKL
                    ************************************************************

WP_001192973.1      AELQERFKRVRITDTSSVFNTDLLYTIELGHGLNVAECMAHSAMARKESRGAHQRLDEGC
NP_418578.1         AELQERFKRVRITDTSSVFNTDLLYTIELGHGLNVAECMAHSAMARKESRGAHQRLDEGC
                    ************************************************************

WP_001192973.1      TERDDVNFLKHTLAFRDADGTTRLEYSDVKITTLPPAKRVYGGEADAADKAEAANKKEKA
NP_418578.1         TERDDVNFLKHTLAFRDADGTTRLEYSDVKITTLPPAKRVYGGEADAADKAEAANKKEKA
                    ************************************************************

WP_001192973.1      NG
NP_418578.1         NG
                    **

```

So, in summary, RefSeq will combine multiple identical sequences into a single `WP_*` multi-species accession. You should therefore expect to find one 100% identical `WP_*` sequence for your query and multiple, *almost* identical `WP_*` entries. And that's precisely what you see here.
> 2 votes

---
Tags: blast, sequence-homology, refseq
---


title: QC measures for NGS sequencing

text:
# Question
Title: QC measures for NGS sequencing
What are good means for performing quality control (QC) or NGS reads?

I'm aware of:

* FastQC
* NGS Screen
* Kraken (e.g., for screening against contaminants)

What are other popular means for such QC?

# Answer
> 12 votes
MultiQC can merge all your different reports into a single one. Which could be useful once you manage to know which QC tools to use.

# Answer
> 7 votes
We routinely run both FastQC and FastQ Screen on all of our raw sequencing reads. FastQ Screen is a tool for detecting cross-species contamination. MGA is another similar tool.

There are then lots of QC tools specific to different types of data, most of which run after alignment. For example RSeQC (RNA data), Qualimap and many many others. Without specifying what kind of data you have this is a bit difficult to make recommendations for though.

Phil

# Answer
> 5 votes
The quality control of ngs reads is heavily dependent on type of the project.

**For genome assembly** projects based on short reads, beside already covered checking quality of sequencing, you would like to look at the kmer spectra to find out, if your reads are going to make sense when they will be translated to De Brujin graph. You will get also a clues about genome coverage, genome size, repetitive content and small clue about heterozygocity. A lot of useful info about interpretation you can find in the README of GenomeScope.

A list tools I used:

* jellyfish \- to count k-mer frequencies
* GenomeScope \- a package for analysis and visualisation of k-mer frequencies (they recommend to use jellyfish for counting k-mer frequencies)
* kmergenie \- for prediction of optimal kmer for assembly

Using these tools can save you a lot of frustration if you accidentally sequence a contaminated sample or if your colleagues / a sequencing company have sent you a wrong file!

# Answer
> 2 votes
The kind of QC you do routinely depends on what your lab's focus is. We do a lot of low-quality, multiplexed DNA and RNA. If you routinely do fresh frozen whole genomes, your QC will be different.

Weighing in from the resequencing side of things (i.e. sequencing an organism that has a good reference), there are several things we are testing with quality control:

* Did the library preparation succeed? Was the capture successful? Is our input material of high enough quality?
* Did the sequencing work? For Illumina, did enough clusters pass filter? What's our Q30? Did we sequence into adapters?
* Can we proceed with analysis? Are the indices on multiplexed samples correct? Do we have the correct reference sequence? Do we have the right targets?

You may use a combination of tools:

* bcl2fastq : if you multiplexed, check the undetermined indices file after converting to fastq to make sure you caught all of the indexed reads
* FastQC: General sequencing quality, sequencing into adapters, handy flagging for obvious problems
* Illumina interop/SAV: Cluster pass filter, Q30
* Alignment (with BWA, bowtie, etc), bedtools, samtools, Picard: check for contamination, make sure we're using the right reference, overlap with target regions, depth of coverage, soft-clipping for low-quality sequence

# Answer
> 0 votes
If you are looking for NGS QC for your fastq, bam, bed and vcf files I would suggest a commercial tool called omnomicsQ.

It automates coverage analysis and sample pass, warn fail according according to your defined SOPs. Unlike fastqc it also comes with a database so that you can compare protocols, samples, and your performance to peer organisations. It also charts performance over time, exceptions,and correlations between metrics

---
Tags: ngs, quality-control
---


# Question
Title: After artificially creating events in a FASTA file, how do I keep track of the old coordinates?
I'm beginning with the reference genome in FASTA format, hg19. I am reading the sequence into a Python dictionary with `BioPython`:

```
genome_dictionary = {}
for seq_record.id in SeqIO.parse(input_fasta_file, "fasta"):
    genome_dictionary[seq_record.id] = seq_record

```

This creates a Python dictionary with keys as chromosome names and values as strings from the FASTA.

My goal is to take the FASTA, and simulate artificial SNPs and InDels along the FASTA to benchmark various algorithms. If I was only creating SNPs in the FASTA, this would be a trivial problem. The reference FASTA has the original coordinates, and the SNPs would also be in a coordinate system which is the same.

However, when I randomly insert an insertion (e.g. an insertion of 12 bp), I move all events into a new coordinate system such that

```
new coordinate system = old coordinate system + 15

```

This becomes very difficult to keep track off throughout the entire genome after several interations. I therefore cannot keep track of the new and old coordinate system without running an algorithm and creating a VCF, which may have errors. This defeats the purpose of the simulation.

The output would be either a BED file or VCF which has keep track of the changes.

What method/data structure is used to solve this problem? Would LiftOver work? Maybe CrossMap?

http://genome.sph.umich.edu/wiki/LiftOver http://crossmap.sourceforge.net/

# Answer
I've written a handful of programs from scratch to simulate mutations and variations in real or simulated sequences.

The trick has always been to sort the variants by genomic coordinate, **apply the variant with the largest coordinate first**, then apply the variant with the second largest coordinate, all the way down to the variant with the smallest coordinate.

Indels and other variants that affect sequence structure only affect *subsequent* coordinates. So going in reverse order ensures variants at the beginning of the sequence don't throw off variants at the end of the sequence.
> 5 votes

# Answer
If you look for a program which would randomly introduce SNPs + short indels and then would save everything into a VCF file, DWGsim or Mason Variator could be a good choice. Then you can create a corresponding Chain file using `bcftools consensus -c` and transform various formats between these two coordinate systems using CrossMap.
> 4 votes

---
Tags: python, bed, structural-variation, indel, simulated-data
---


title: How does blastn sort output?

text:
# Question
Title: How does blastn sort output?
I used blastn to search one a genome database for the sequences in a file.

```
blastn -out output.txt -outfmt 6 -query sequence_list -db genome_database -perc_identity 100

```

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

I also checked if it was sorted by which contig the hit was in and this is not the case.

# Answer
> 1 votes
**Edit** I stand corrected. @Daniel Standage correctly pointed out that the below is only true for a single query, or for the sorting of hits within those found by a given query. Otherwise they're sorted in the order they were given.

BLAST results are sorted in descending order by the "Bit score" (column 12). Sorting by E. value would not be quite the same.

i.e. column 11 (E value) will be presented in ascending order (an E. value of 0 is as good as it gets), and as it gets higher, your hit is getting worse. Likewise, bitscores are highest at the top and decrease down the list. A hit with E. value = 0 should have a high score. How high exactly depends on the hit length, so it's common to see many E values with 0.0, but the scores differ. Here's a snippet of one of my blast results:

```
PVCcif_ATCC43949  PAU_06042014  100.000  24180  0     0   1      24180  2233012  2257191  0.0  43606
PVCcif_ATCC43949  PAU_06042014  78.098   5593   1164  16  2765   8323   3219956  3225521  0.0  4547
PVCcif_ATCC43949  PAU_06042014  77.549   5394   1139  16  2948   8323   2376163  2370824  0.0  4260
PVCcif_ATCC43949  PAU_06042014  74.140   5669   1342  34  8490   14100  3225688  3231290  0.0  3573
PVCcif_ATCC43949  PAU_06042014  76.835   3665   825   9   16507  20150  3234191  3237852  0.0  2764
PVCcif_ATCC43949  PAU_06042014  76.871   2646   571   11  16625  19248  2362597  2359971  0.0  2001
PVCcif_ATCC43949  PAU_06042014  70.733   4131   1005  51  9047   13063  2368869  2364829  0.0  1954
PVCcif_ATCC43949  PAU_06042014  75.040   1891   448   7   8500   10383  2370771  2368898  0.0  1272
PVCcif_ATCC43949  PAU_06042014  78.986   1361   269   7   22766  24119  890051   891401   0.0  1150
PVCcif_ATCC43949  PAU_06042014  77.819   1082   217   5   536    1610   3217983  3219048  0.0  868

```

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

# Answer
> 2 votes
BLASTN should definitely report alignments in the order that the query sequences were provided.

Based on your description and my personal experience, I think the most likely explanation is some kind of mix up. It happens to all of us, especially when we're in a rush to get an answer, but even sometimes when we're being careful and disciplined.

If you can go back and recreate the result, I would be very interested to see the two pairs of input and output files. If not, I wouldn't spend too much more time worrying about it. :-)

**UPDATE** Alignments are reported in the order that the query sequences are provided, even in tabular mode. The order of each query's hits is sorted by bitscore or e-value, but that is a secondary sorting. I've included an example below that has two queries, `seq1` and `seq2` with `-outfmt 6`.

```
seq1    scaffold2       100.000 280     0       0       1       280     11131   11410   7.73e-146       518
seq1    scaffold294     88.693  283     25      4       1       279     12160   12439   6.46e-92        339
seq1    scaffold1525    85.315  286     29      9       1       277     27142   26861   3.07e-75        283
seq1    scaffold1874    83.505  291     32      10      1       280     8558    8273    1.86e-67        257
seq1    scaffold147     82.653  294     34      9       1       279     10797   10506   1.45e-63        244
seq1    scaffold2478    82.374  278     38      8       11      280     7622    7896    1.13e-59        231
seq1    scaffold3405    87.940  199     17      5       1       196     14390   14584   1.46e-58        228
seq1    scaffold395     83.060  183     23      4       106     280     12021   11839   5.40e-38        159
seq1    scaffold2853    85.965  114     13      3       1       113     20864   20753   9.16e-26        119
seq2    scaffold9       100.000 140     0       0       1       140     1121    1260    2.34e-68        259
seq2    scaffold2950    90.909  110     6       3       1       106     110     1       6.84e-34        145
seq2    scaffold3416    95.385  65      3       0       1       65      11882   11818   1.16e-21        104
seq2    scaffold3103    95.385  65      3       0       1       65      9736    9672    1.16e-21        104
seq2    scaffold5297    88.608  79      7       1       64      140     150     228     6.99e-19        95.3

```

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

---
Tags: sequence-alignment, blast, sequence-analysis
---


# Question
Title: What is a good pipeline for using public domain exomes as controls?
I'm currently attempting association analysis with an extremely small set of patient exomes (n=10), with no control or parental exomes available. Downloading the ExAC VCF of variant sites (http://exac.broadinstitute.org/downloads) or the 1000G integrated call sets (http://ftp.1000genomes.ebi.ac.uk/) and combining this with our pooled patient VCFs has not been successful (I suspect the approach of attempting to merge such large VCFs generated from different pipelines is rather naive).

Looking at the primary literature, I have gathered it should be possible to use these resources to help increase statistical power for our analysis. My question is how do I take these large .vcfs with many samples and successfully merge them to our patient .vcfs, such that the combined VCF can be used downstream to run analysis packages? (PODKAT, PLINK, etc.)

# Answer
I gather that you want to use the background SNP frequency as a prior for input to your SNP calling algorithm?

I'm not sure of a canned algorithm for doing this, but a quick google shows up some promising links:
> 1 votes

---
Tags: ngs, public-databases, variation, exome
---


# Question
Title: what is the difference between local identity and global identity for homolog finding?
We believe that if after running blast, the global identity between a resulting sequence and our query is at least 30%, we can say that two sequences are homologs. what is the difference between local identity and global identity and how can we calculate them? the following file is the global alignment result of one of my PSI-BLAST hits against query:

```
Aligned_sequences: 2
# 1: NP_418577.1
# 2: WP_036312822.1
# Matrix: EBLOSUM62
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 264
# Identity:      85/264 (32.2%)
# Similarity:   125/264 (47.3%)
# Gaps:          29/264 (11.0%)
# Score: 376.5
#
#
#=======================================

NP_418577.1        1 -----------MAEMKNLKIE-------VVRYNPEVDTAPHSAFYEVPYD     32
                                ..|..:..|:       :.|::||||..|....|:|...
WP_036312822.      1 MTIVDSGAPADTQEANDSGIQSYLVTFIIRRFDPEVDAEPRWVDYDVEMY     50

NP_418577.1       33 ATTSLLDALGYIKDNLAPDLSYRWSCRMAICGSCGMMVNNVPKLACKTFL     82
                     .|..:||||..||.::...||:|.||...||||..|.:|...:|||||.:
WP_036312822.     51 PTDRVLDALHRIKWDVDGTLSFRRSCAHGICGSDAMRINGRNRLACKTLI    100

NP_418577.1       83 R--DYTDGMKVEALANFPIERDLVVDMTHFIESLEAIKPYIIGNSRTADQ    130
                     :  |.:..:.|||:...|:|:||:|||..|.||...::|::...|.....
WP_036312822.    101 KDLDISKPIYVEAIKGLPLEKDLIVDMDPFFESFRDVQPFLQPKSAPEPG    150

NP_418577.1      131 GTNIQTPAQMAKYHQFSGCINCGLCYAACPQFGLNPEFIGPAAITLAHRY    180
                     ....|:....|.|...:.||.|..|.::||.|..:.::.|||||..|||:
WP_036312822.    151 KERFQSIKDRAVYDDTTKCILCAACTSSCPVFWTDGQYFGPAAIVNAHRF    200

NP_418577.1      181 NEDSRDHGKKERMAQLNSQNGVWSCTFVGYCSEVCPKHVDPAAAIQQGKV    230
                     ..||||.....|:..||.:.|||.|.....|:|.||:.::...||.:.|.
WP_036312822.    201 IFDSRDDAADVRLDILNDKEGVWRCRTTFNCTEACPRGIEITKAIAEVKQ    250

NP_418577.1      231 ESSKDFLIATLKPR    244
                     ...:.
WP_036312822.    251 AVLRG---------    255

#---------------------------------------
#---------------------------------------
the following file is the local alignment of that hit against query:
Aligned_sequences: 2
# 1: NP_418577.1
# 2: WP_036312822.1
# Matrix: EBLOSUM62
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 219
# Identity:      83/219 (37.9%)
# Similarity:   119/219 (54.3%)
# Gaps:           2/219 ( 0.9%)
# Score: 390.5
#
#
#=======================================

NP_418577.1       13 RYNPEVDTAPHSAFYEVPYDATTSLLDALGYIKDNLAPDLSYRWSCRMAI     62
                     |::||||..|....|:|....|..:||||..||.::...||:|.||...|
WP_036312822.     31 RFDPEVDAEPRWVDYDVEMYPTDRVLDALHRIKWDVDGTLSFRRSCAHGI     80

NP_418577.1       63 CGSCGMMVNNVPKLACKTFLR--DYTDGMKVEALANFPIERDLVVDMTHF    110
                     |||..|.:|...:|||||.::  |.:..:.|||:...|:|:||:|||..|
WP_036312822.     81 CGSDAMRINGRNRLACKTLIKDLDISKPIYVEAIKGLPLEKDLIVDMDPF    130

NP_418577.1      111 IESLEAIKPYIIGNSRTADQGTNIQTPAQMAKYHQFSGCINCGLCYAACP    160
                     .||...::|::...|.........|:....|.|...:.||.|..|.::||
WP_036312822.    131 FESFRDVQPFLQPKSAPEPGKERFQSIKDRAVYDDTTKCILCAACTSSCP    180

NP_418577.1      161 QFGLNPEFIGPAAITLAHRYNEDSRDHGKKERMAQLNSQNGVWSCTFVGY    210
                     .|..:.::.|||||..|||:..||||.....|:..||.:.|||.|.....
WP_036312822.    181 VFWTDGQYFGPAAIVNAHRFIFDSRDDAADVRLDILNDKEGVWRCRTTFN    230

NP_418577.1      211 CSEVCPKHVDPAAAIQQGK    229
                     |:|.||:.::...||.:.|
WP_036312822.    231 CTEACPRGIEITKAIAEVK    249

#---------------------------------------
#---------------------------------------

```

# Answer
> 3 votes
If you want a purely alignment based approach to begin to *infer* homology (note I say *infer*, as we stated in the comments, this isn't sufficient to say 2 sequences are homologs), I would use a Multiple Sequence Alignment (MSA) instead of a pairwise like you've done so far, and include as many similar proteins as you can from BLAST.

The choice of algorithm for the MSA probably won't matter too much, as you will be looking for conserved regions anyway.

Once you have your MSA, you can easily spot conserved domains. If anything is known about this protein structurally, it's quite trivial to work out whether a conserved region is a functional domain (active site, that kind of thing).

With this information, you would be somewhat justified in calling these proteins structural homologs at least if they share domains, but, to reiterate, without doing phylogenetic analyses, you *cannot* call them true homologs.

# Answer
> 1 votes
Defining homology isn't as simple as an identity cutoff, however, you're right that homology seems likely given a certain identity. This is where e-value comes into play. See: https://en.wikipedia.org/wiki/Expected\_value

To answer the question, global-global requires that the alignment account for every base in the query and the target sequence. local-global requires a the full length of the query on a part of the target. local-local allows only the highest scoring part of the alignment between the query and target to be reported. i.e. one functional domain against another. these hits are often called HSPs.

---
Tags: sequence-homology
---


# Question
Title: How to represent trans-spliced genes in GTF?
For example, see this gene (nad1) in ENA: http://www.ebi.ac.uk/ena/data/view/ABI60879

If you look at the XML for that gene you see the following:

```
join(
             DQ984518.1: 324706 .. 325091 ,
  complement(DQ984518.1:  24417 ..  24498),
  complement(DQ984518.1:  22828 ..  23019),
             DQ984518.1:   3484 ..   3542 ,
  complement(DQ984518.1: 153702 .. 153960)
)

```

Which shows 5 exons joined out of phase and out of order. Is there a valid GTF representation of this?

How to dump a 'non-canonically spliced' gene into GTF? i.e. what's the recommendation?

Cross-posted on biostars

# Answer
To my knowledge there's no defined way to deal with that in GTF. GFF3 handles trans-splicing (you'll have to scroll down to "trans-spliced transcript") by giving an individual transcript multiple parents (e.g., `ID=some_transspliced_gene;Parent=gene1,gene2`). You could use the same methodology with GTF files, but just note that it'll break most downstream programs.
> 3 votes

---
Tags: file-formats, gtf, rna-splicing
---


# Question
Title: Find gene at position from gff or gbk file
I have a VCF file with SNPs from a bacterial genome and want to find if the SNPs are located inside genes, is there some CLI-tool where you can pass a VCF file and a gff or gbk file and it returns the name of the genes?

# Answer
Via BEDOPS:

```
$ gff2bed < annotations.gff > annotations.bed
$ vcf2bed < snps.vcf > snps.bed
$ bedmap --echo --echo-map-id-uniq snps.bed annotations.bed > answer.bed

```

This can be reduced to a one-liner if you're using `bash`:

```
$ bedmap --echo --echo-map-id-uniq <(vcf2bed < snps.vcf) <(gff2bed < annotations.gff) > answer.bed

```
> 3 votes

---
Tags: snp
---


# Question
Title: Correlating gene expression with qualitative variables
I have a gene expression dataset that I want to investigate. Particularly, I would like to understand whether there is any correlation between each gene's expression and some quantitative or qualtitative data (say, correlation between gene 'XPTO' , body mass index, and race).

One possible way to test this would be through logistic regression, but is this a good approach or are there caveats that I should know about using such a statistic?

My question is the following: which methods would you advise to measure such correlations, and why?

(this post was crossposted on Biostars)

# Answer
Logistic regression would generally be a bad choice for non-binary outcomes. In such cases, linear regression (or a GLM more generally) still works fine. You can already do that in the standard R packages for RNAseq (DESeq2, edgeR, and limma), where the fold-change is in whatever units you're measuring your quantitative trait in.
> 3 votes

---
Tags: rna-seq, r
---


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

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

I've searched in PDB and UNIPROT.

# Answer
Tried looking for an explicit database? i.e.

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

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

---
Tags: database, protein-structure
---


# Question
Title: qPCR: Why is fold change and standard deviation calculated after transformation?
I am analyzing data from a quantitative polymerase chain reaction (qPCR) using R. After cleaning the raw data, it looks something like this:

```
> dput(x)
structure(list(Reporter = c("FAM", "FAM", "FAM", "FAM", "FAM",
"FAM", "FAM", "FAM", "FAM", "FAM", "FAM", "FAM", "VIC", "VIC",
"VIC", "VIC", "VIC", "VIC", "VIC", "VIC", "VIC", "VIC", "VIC",
"VIC"), Number = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L,
12L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L), A = c(22.19,
22.24, 22.5, 22.54, 22.6, 22.59, 23.07, 23.46, 22.43, 22.74,
24.09, 23.91, 24.52, 25.03, 25.25, 25.82, 25.13, 24.71, 25.34,
25.85, 25.25, 26.15, 25.81, 25.29), B = c(21.72, 21.78, 22.86,
22.73, 19.88, 20.07, 21.06, 21.06, 20.96, 21.11, 19.46, 19.43,
24.75, 24.69, 25.64, 25.19, 23.76, 23.69, 24.35, 25.05, 24.1,
23.81, 22.81, 23.13), C = c(21.37, 21.56, 20.07, 20.01, 21.17,
21.08, 20.54, 20.36, 33, NA, NA, NA, 23.91, 24.31, 23.61, 23.88,
24.33, 24.31, 23.37, 23.53, 33, NA, NA, NA), E = c(26.26, 27.33,
25.93, 26.56, 25.76, 23.03, 24.72, 25.27, 24.43, 24.31, 26.98,
23.33, 24.04, 25.02, 25.1, 25.1, 24.68, 25.48, 25.87, 26.22,
25.35, 25.36, 25.11, 25.98), F = c(25.81, 26.9, 25.58, 26.61,
25.06, 21.85, 23.59, 24.04, 23.19, 23.19, 25.17, 20.8, 24.12,
24.26, 25.32, 25.25, 24, 23.78, 24.7, 24.48, 23.52, 23.87, 23.05,
23.05), G = c(26.12, 27.02, 24.08, 25.15, 25.99, 23.18, 24.2,
24.05, 33, NA, NA, NA, 23.47, 23.45, 23.7, 23.74, 24.46, 24.19,
23.56, 23.53, 33, NA, NA, NA)), .Names = c("Reporter", "Number",
"A", "B", "C", "E", "F", "G"), class = c("tbl_df", "tbl", "data.frame"
), row.names = c(NA, 24L))

```

I think that each well measures two genes: target which is 12 different genes (color FAM), and the reference or housekeeping gene, GAPDH (color VIC). Also, control is triplicate A, B, C (3 x 12 wells), and treatment is E, F, G (3 x 12 wells).

I have created a function `ddCt_` for analyzing the data with the ddCt algorithm (Livak & Schmittgen, 2001). It takes one argument `x` which is a data.frame of the form exemplified above.

```
ddCt_ <- function(x) {
  # Subset x by control/treatment and target/reference
  # Then, calculate Ct averages for each triplicate
  TC <- x %>%
    select(Reporter, Number, A, B, C) %>%
    filter(Reporter == "FAM") %>%
    rowwise() %>%
    mutate(Ct_Avg = mean(c(A, B, C), na.rm = TRUE))

  RC <- x %>%
    select(Reporter, Number, A, B, C) %>%
    filter(Reporter == "VIC") %>%
    rowwise() %>%
    mutate(Ct_Avg = mean(c(A, B, C), na.rm = TRUE))

  TT <- x %>%
    select(Reporter, Number, E, F, G) %>%
    filter(Reporter == "FAM") %>%
    rowwise() %>%
    mutate(Ct_Avg = mean(c(E, F, G), na.rm = TRUE))

  RT <- x %>%
    select(Reporter, Number, E, F, G) %>%
    filter(Reporter == "VIC") %>%
    rowwise() %>%
    mutate(Ct_Avg = mean(c(E, F, G), na.rm = TRUE))

  # Normalize Ct of the target gene to the Ct of the reference gene
  dCt_control   <- TC$Ct_Avg - RC$Ct_Avg
  dCt_treatment <- TT$Ct_Avg - RT$Ct_Avg

  # Normalize dCt of the treatment group to the dCt of the control group
  ddCt <- dCt_treatment - dCt_control

  # Calculate fold change
  fc <- 2^(-ddCt)

  # Calculate avg and sd
  dCt_control_avg   <- mean(dCt_control)
  dCt_control_sd    <- sd(dCt_control)
  dCt_treatment_avg <- mean(dCt_treatment)
  dCt_treatment_sd  <- sd(dCt_treatment)

  # Create output
  df <- data_frame(
    Sample = 1:12, "dCt Control" = dCt_control, "dCt Treatment" = dCt_treatment,
    "ddCt" = ddCt, "Fold Change" = fc, "dCt Control Avg" = dCt_control_avg,
    "dCt Control SD" = dCt_control_sd, "dCt Treatment Avg" = dCt_treatment_avg,
    "dCt Treatment SD" = dCt_treatment_sd
    ) %>% round(2)
  write.csv(x = df, file = "result.csv")
  saveRDS(object = df, file = "result.rds")
  df
}

```

However, I am not a biochemist or molecular biologist, so I am unsure of the basic concepts involved and the overall approach. So, my questions are:

1. Is the function correct?
2. Why is it important to calculate fold change and standard deviation after normalization?

# Answer
1. The functions look correct, but calculate a few by hand and ensure they match. One thing I should note is that the subtraction of the Ct values usually happens before an average is made, since generally both are in the same well. In other word, subtract the reference Ct from the gene of interest Ct from the same well and then average the dCts. This order is important since one usually does a multiplexed qPCR.
2. You have to normalize to a reference gene to control for how much cDNA was used, since that will alter the Ct values. If you calculated the fold-changes without normalization then they could be purely due to using more/less cDNA in the reaction (i.e., the output would be meaningless).
> 4 votes

---
Tags: r, normalization
---


# Question
Title: How to dump genes from GenBank in GFF3 format?
*This question has also been asked on BioStars*

If I look at this record in GenBank I see about 6k genes: https://www.ncbi.nlm.nih.gov/nuccore/CM000760?report=gbwithparts

I'd really like to be able to dump those genes in GFF3 format, but I'm guessing I can't do that in a single step, so...

1. How to dump all genes from a GenBank record in 'whatever' format?
2. How to convert 'whatever' format to GFF3?

# Answer
> 3 votes
Turns out you can just grab the GFF3 from the NCBI's FTP site!

ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/003/195/GCA\_000003195.3\_Sorghum\_bicolor\_NCBIv3/

See:

Or "Access the data" on the right here:

# Answer
> 1 votes
When accessing all of the annotated genes for a reference genome, downloading a GFF3 file directly from the Genbank or RefSeq FTP sites is definitely the way to go.

But when submitting a search query into NCBI Entrez, there is no easy way to download the results in GFF3 format. Probably the best route is to download the results in ASN.1 format, and then use `annotwriter` from the NCBI C++ toolkit to convert to GFF3. But getting the latter to compile and run properly isn't trivial in my experience.

---
Tags: gff3, data-download
---


# Question
Title: How can I find the relevant pathway map from gene-gene or protein-protein interaction list?
I have some difficulties to understand/interpret the pathway map and how a gene-gene interaction list or DNA sequencing can map into pathways.

In addition what's the difference between MARK/ERK pathway and the output of KEGG

# Answer
See a project like Plant reactome, where they infer the pathways of all sequenced plant genomes using orthology to well annotated species. http://plantreactome.gramene.org/
> 2 votes

# Answer
As you alluded to in your question, KEGG does provide a curated set of enzyme and metabolite information. This information can be parsed and used to create a network that you can analyze to look at how gene products could be working together. Additional steps could include mapping transcriptome data to the map.

See this paper as an example of how this can be done.

I hope this helps. :)
> 2 votes

---
Tags: gene, ngs, phylogenetics, pathway, kegg
---


# Question
Title: Intersection of two genomic ranges to keep metadata
I am trying to find intersection of two genomic ranges (gr1 and gr2) and keep metadata from one of them

```
gr1
chrI [1, 100] * | 0.1
chrI [101, 200] * | 0.2

gr2
chrI [50, 150] + |

```

`intersect(gr1, gr2)` will remove metadata, where `subsetByOverlaps(gr1, gr2)` return:

```
chrI [1, 100] * | 0.1
chrI [101, 200] * | 0.2

```

where metadata is saved, but intersection is bigger than real (the same as in initial gr1 start/end)

How can I get an intersection range and keep metadata? to get the range:

```
chrI [50, 100] * | 0.1
chrI [101, 150] * | 0.2

```

# Answer
Presuming you don't have cases of multiple ovelaps

```
> library(GenomicRanges)
> gr1 = GRanges(c("chrI", "chrI"), IRanges(start=c(1, 101), width=c(100, 100)), mcols=data.frame(foo=c(0.1, 0.2)))
> gr2 = GRanges(c("chrI"), IRanges(start=c(50), width=c(101)))
> o = findOverlaps(gr1, gr2)
> grl1 = split(gr1[queryHits(o)], 1:length(o)) # You can't mendoapply on a GRanges object
> grl2 = split(gr2[subjectHits(o)], 1:length(o))
> foo = function(x, y) {
+     rv = x
+     start(rv) = max(start(x), start(y))
+     end(rv) = min(end(x), end(y))
+     return(rv)
+ }
> unlist(mendoapply(foo, grl1, y=grl2))

```

The output is:

```
GRanges object with 2 ranges and 1 metadata column:
    seqnames     ranges strand | mcols.foo
       <Rle>  <IRanges>  <Rle> | <numeric>
  1     chrI [ 50, 100]      * |       0.1
  2     chrI [101, 150]      * |       0.2
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

```

It'd be more memory efficient to iterate over `o`, but this gives you the idea.

In cases of multiple overlaps, one could use `intersect()`, then `findOverlaps()` with that and then continue on with the procedure demonstrated above (or iterate over the overlaps).
> 4 votes

---
Tags: r, bioconductor
---


# Question
Title: Publicly available, free, complete database for antibiotics names and classes?
This is a tough one I think: is there a publicly available, up-to-date, free, complete database for antibiotics names and classes?

I am specifically looking for information like, e.g., `cefpirome (is a) cephalosporin`.

I've looked in a couple of different places, but nothing really comes close. E.g.:

1. CARD comes closest, but is not free (companies need to pay a license) and not complete/up-to-date (e.g., cefpirome is missing)
2. @Karel mentioned KEGG in his answer, but this also needs a license (companies for sure, academics if they download data via FTP).
3. Novel Antiobiotics Database has no class ... and stops in 2003
4. The Antimicrobial Index is web only, and terms of use apparently prohibit any kind of web-scraping.
5. A suggestion by @Pierre on Twitter was to try Antimicrobial stubs on Wikipedia. Good idea, but the class is not stored in the ChemBox and must be free-text parsed. Or writing a script which evaluates the Wikipedia category tags.
6. I drew blank at different government agencies (EMA, FDA, EFSA, NIH, etc.). The closest I found was from the WHO a list of critically important antibiotics. But not complete, and in a PDF.

Any pointer appreciated.

# Answer
**ChEBI**

It should be possible to get this information from the ChEBI database, see the exported tables. You could download the ontologies (in OWL / OBO), parse them using some ad-hoc parser or using a dedicated library (e.g., Pronto), build a directed acyclic graph based on the `is_a` relationship, and extract the antibiotic subtree.

**KEGG**

You could also use KEGG: Kyoto Encyclopedia of Genes and Genomes, but a subscription is probably required. The exported drug table should be easy to parse (e.g., in Python).

For a quick verification of the completeness, you can download the table, filter out everything except antibacterial agents and check the number of records:

```
curl "ftp://ftp.genome.jp/pub/kegg/medicus/drug/drug" > keg_all.txt

cat keg_all.txt \
   | perl -pe 's/\n/BREAK/g' \
   | perl -pe 's@///BREAK@///\n@g' \
   | grep 'CLASS       Antibacterial' \
   | perl -pe 's/BREAK/\n/g' \
   > keg_antibacterial.txt

cat keg_antibacterial.txt | grep "///" | wc -l

```

There are 646 records. However, some of them are almost identical (e.g., only 1 atom difference).

If you add `| grep -i "antibiotic" \` to the pipeline, then you get only 445 records.
> 6 votes

---
Tags: database
---


# Question
Title: How to maximize fastq parsing with FastqGeneralIterator (Bio.SeqIO.QualityIO)
I'm contributing to a python-based project that uses `Biopython` to analyze fastq files. It currently uses `SeqIO.parse`, which populates various structures with all of the fastq information (including converting quality scores). There is apparently a faster (lighter-weight) parser called FastqGeneralIterator that doesn't populate all of these items.

I'd like to use `FastqGeneralIterator`, but be able to perform some of the functions `SeqIO.parse` offers on the fly (not for *every* sequence). For example, I'd like to convert base quality to PHRED scores for specific sequences, but I don't see that function available.

Is there an easy way to take a tuple output by `FastqGeneralIterator` and convert it into a proper `Bio.SeqIO` `SeqRecord`?

# Answer
I had not come across `FastqGeneralIterator` before, but I will start using it for some of my work!

One answer is to replicate the code in the `FastqPhredIterator`, which does exactly what you are looking to do.

However, My personal preference would be to use the biopython infrastructure as much as possible, especially when it comes to the tricky aspects of record parsing. If the records you want extra info on are uncommon, then using `FastqPhredIterator` for phred quality extraction might work for you:

```
import io
from Bio.SeqIO.QualityIO import FastqGeneralIterator, FastqPhredIterator

with open("input.fastq", "rU") as handle:
    for title, sequence, quality in FastqGeneralIterator(handle):
        # do something with title, sequence, and quality
        # Set special_case to True when you want more info on the record
        special_case = False
        if special_case:
            record_stream = io.StringIO("\n".join([title, sequence, "+", quality]))
            record = next(FastqPhredIterator(record_stream))
            # do something special with record

```

What we are doing here is combining the record strings into one and dumping them into a `StringIO` object, which implements the file handle duck-type. We then hand the in-memory handle to `FastqPhredIterator` to parse the FASTQ record again. Make sure to set `special_case` to `True` only when you want extra info on a record, as the `special_case` branch makes at least two extra copies of the FASTQ record.
> 3 votes

---
Tags: fastq, python
---


# Question
Title: Which quality score encoding does PacBio use?
Do you know which quality score encoding PacBio uses now? I know some of their file formats have changed in the past year or two, but I haven't found much on their quality score encoding.

The most recent answer I found is from 2012, where one user states it's "phred style", but doesn't state whether it's 33 or 64: DeNovo assembly using pacBio data

# Answer
> 12 votes
PacBio does use PHRED 33, *but it turns out the question may be irrelevant* for the newer PacBio Sequel Sequencer, because it reports all base qualities as PHRED 0 (ASCII `!`). The RS-II reports PHRED 33 quality scores.

The Sequel provides data in the form of *subreads*, which are the circular consensus sequences (CCS) from a single zero-mode waveguide (ZMW). Essentially, a single circularized DNA molecule enters a ZMW well, and is sequenced as many times as possible. (Sometimes multiple molecules get into a single well, but that's not the normal case.) The number of times this single molecule can be sequenced depends on several variables, including molecule length. The consensus sequence for this single molecule is then reported out as a subread.

As I started playing with the quality scores in my data, I realized the quality score for every base is `!`, which equates to a PHRED quality of 0. I contacted someone at PacBio who informed me that they learned the base quality values were not any more useful than the bases themselves for generating a high-quality circular consensus sequence. Thus, rather than spend the computational effort to calculate the scores, they skipped it, and reported `!` for every base. They also report an overall read quality of 0.8, but this is simply a placeholder. They are working to make this information more readily available.

---
Tags: ngs, pacbio
---


# Question
Title: NaN values after ComBat analysis on TCGA COAD RNAseq
I have FPKM-UQ data from COAD-TCGA.

I generated an expression set of this data using:

```
> edata = log2(data + 1)
> edata[1:2,1:2]
      X01240896.3f3f.4bf9.9799.55c87bfacf36
1                                  8.540967
10                                13.528968
100                               16.296422
1000                              13.658199
10000                             15.143788
      X01ad5016.f691.4bca.82a0.910429d8d25b
1                                  9.886537
10                                16.719682
100                               17.212312
1000                              10.317842
10000                             13.166767

```

Row name is entrezid, and column names are case\_ids.

The batch data provided by TCGA which indicates a total batch num = 42

```
> omics[1:5,]
                                      snf    boo order batch
X01240896.3f3f.4bf9.9799.55c87bfacf36   3 192807     1    22
X01ad5016.f691.4bca.82a0.910429d8d25b   2  82213     2     2
X01f493d4.229d.47a6.baa8.32a342c65d01   1 A08907     3    29
X022f39e9.57ee.4b2b.8b3a.8929e3d69a37   2 A16S13     4    32
X02f9668c.71e6.485f.88b1.b37dc8bdd2ab   3 102113     5     7
> batch = omics$batch
> batch
  [1] 22  2 29 32  7  1  8  8  7 41 26 14 26 26  1  5 22  5 15 18 28 12 12  8 19
 [26] 13 38  5 27  7 41  1 22 26 28 22 33 14 13  8 10 13 11 17  5 38  5  7  7 31
 [51] 22 42 11  8  7 17 15 17 12 17  2 19  2 12 18 14 22 11 21 15 17 17 33 15 19
 [76] 15  6  1 10 11  2 28 36 17  9 15 12 17 15 10  6  8 13 25 11 15  1 13 24  8
[101] 13 15 17  8  2 22 17 11 13 37 19 38 13 19  5 16  5  5 22 35 35 19 13 17 22
[126]  7  8  8 41 39 14  6 19 18 38 22 14 30 19 24 11 14 13 16  7 13  8 22 12  8
[151] 22 17  2 40 13 35 15  5 24 19 41 22 22 19 17  8 13 33 10 17 29 12  9 11 24
[176] 27  5 14  7 25 13 35  9 18 16 13 31 25 18 10 35 18 14 33  5 19 17  8 11 11
[201] 42  7 10  8 39 19  8  7 11 19 31 35 28 36 13 37 10 38  3 18 13 34 39 13  3
[226] 40 31 17 19 26  7 19 12 22 19 17 12 17 13 12 28 33 37 17 35 18 40 13 20 35
[251] 17  6 15 14 12 27  5 41 31 41 18 21 19 17 40  8 41  2  5 15  5 15 15  5 39
[276]  1 22 33 13 17  7 11 15  3  7 38 15  7 14 14 22 31 14 18  3 19  7  9  6  5
[301] 22  8 11  1 13 10 19  8 14 15  6 32  1  2 19  4  5 13 18 42 11 17  8 36 19
[326] 15 22 28 15  8 11  8 15 38 19 39  7 11 42 23  5  1 22 17 35 13 40 25  7 41
[351] 38  5  8 19 31  6 38  8 36  4 15 15  7 12 22 38 29 35 15  2 10  2 19 15 25
[376]  5 39 30 19 18 15 10 10 33  1 19 32 19 19 15 20 15  7 22 35 22 30 14 12 16
[401] 21 13 15 18  7  5  6 31 12  7 15 32 33 18 14 15 15 26 31  1 14 19 17 18  1
[426]  5 10 32

```

Unfortunately, when I run the ComBat analysis it generates all NANs, and no values:

```
> modcombat = model.matrix(~1, data=omics)
> combat_edata = ComBat(dat=edata, batch=batch, mod=modcombat, par.prior=TRUE, prior.plots=FALSE)
> combat_edata[1:5,1:2]
      X01240896.3f3f.4bf9.9799.55c87bfacf36
1                                       NaN
10                                      NaN
100                                     NaN
1000                                    NaN
10000                                   NaN
      X01ad5016.f691.4bca.82a0.910429d8d25b
1                                       NaN
10                                      NaN
100                                     NaN
1000                                    NaN
10000                                   NaN

```

# Answer
I think I might know what the problem is.

I plotted the log2 corrected values, and the distribution may not be normal because of the 0 values present in RNASeq data. If I understand correctly, this data is actually missing, whereas with my current analysis it is viewed by R as an actual value of 0. So either the missing data needs to be imputed, or the data needs to be scaled and normalized to mean of 0, then subsequently log2 transformed.

```
> data <- standardNormalization(data)
> data <- log2(data + 1)
> combat_data = ComBat(dat=data, batch=batch, mod=modcombat, par.prior=TRUE, prior.plots=FALSE)
> combat_data[1:5,1:2]
          X01240896.3f3f.4bf9.9799.55c87bfacf36
100507661                           -0.13857428
57103                               -0.06700506
22838                                0.16295096
55567                               -0.13579991
6147                                 1.64990772
          X01ad5016.f691.4bca.82a0.910429d8d25b
100507661                           -0.13789534
57103                               -0.06210232
22838                                0.01546224
55567                               -0.13443289
6147                                 1.48473072

```
> 2 votes

---
Tags: r, combat
---


# Question
Title: how to set database other than nr for remote blast+ search
I am attempting to run a BLAST search remotely using BLAST+.

I can get search to work correctly at the command line with the following commands:

```
blastp -query proteins.fasta -remote -db nr -out proteins_nr.txt -outfmt 6 -evalue 1e-30

```

However, I would like to us the remote database titled "Microbial proteins from nr" which is what would be used for the microbial blast search accessing BLAST from the NCBI site.

I'm struggling to find the correct code to access this database instead of just "nr".

Is there a list with the codes of different databases?

Thanks in advance

# Answer
I believe you're looking for `env_nr`? It's listed as such, under `Metagenomic proteins` in the blastp webpage. It appears that the word within brackets should be supplied alongside the `-db` parameter. A quick test with a dummy amino acid fasta file does turn up a result to a valid NCBI protein accession.

Edit: I followed up on this with a little more digging around, and it appears the database you're looking for is called `Microbial_proteins`. I ended up quote-searching `Microbial proteins from nr` on Google and exactly 2 hits turn up. The first hit, I think, is what the OP posted a screenshot of in a comment to this answer. The second is a list of databases along with keywords from a Spain-based bioinformatics firm: http://data.biobam.com/ncbi\_blast\_dbs\_protein.pdf. I retried the dummy example with the `Microbial_proteins` database, and it seems to work. Again, I am not sure why I can't easily find this on an NCBI resource. Also, I am curious as to how different this database is from `env_nr`.
> 9 votes

---
Tags: database, blast
---


# Question
Title: Read length distribution from FASTA file
I have a single ~10GB FASTA file generated from an Oxford Nanopore Technologies' MinION run, with \>1M reads of mean length ~8Kb. How can I quickly and efficiently calculate the distribution of read lengths?

A naive approach would be to read the FASTA file in Biopython, check the length of each sequence, store the lengths in a numpy array and plot the results using matplotlib, but this seems like reinventing the wheel.

Many solutions that work for short reads are inadequate for long reads. If I'm hey output a single (text) line per 1/10 bases, which would lead to a text output of upwards of 10,000 lines (and potentially more than 10x that) for a long read fasta.

# Answer
> 20 votes
If you want something quick and dirty you could rapidly index the FASTA with `samtools faidx` and then put the lengths column through R (other languages are available) on the command line.

```
samtools faidx $fasta
cut -f2 $fasta.fai | Rscript -e 'data <- as.numeric (readLines ("stdin")); summary(data); hist(data)'

```

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

# Answer
> 12 votes
Statistics for nanopore reads are tricky because of the huge range of read lengths that can be present in a single run. I have found that the best way to display lengths is by using a log scale on both the x axis (length) and the y axis (sequenced bases, or counts, depending on preference).

I have written my own scripts for doing this: one for generating the read lengths, and another for plotting the length distribution in various ways. The script that generates read lengths also spits out basic length summary statistics to standard error:

```
$ ~/scripts/fastx-length.pl > lengths_mtDNA_called.txt
Total sequences: 2110
Total length: 5.106649 Mb
Longest sequence: 107.414 kb
Shortest sequence: 219 b
Mean Length: 2.42 kb
Median Length: 1.504 kb
N50: 336 sequences; L50: 3.644 kb
N90: 1359 sequences; L90: 1.103 kb

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

```

Here are a couple of of the produced graphs:

The scripts to generate these can be found here:

# Answer
> 7 votes
Using Biopython and matplotlib would seem like the way to go, indeed. It really just boils down to three lines of code to get that graph:

```
import Bio, pandas
lengths = map(len, Bio.SeqIO.parse('/path/to/the/seqs.fasta', 'fasta'))
pandas.Series(lengths).hist(color='gray', bins=1000)

```

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

```
#!/usr/bin/env python2

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

Written by Lucas Sinclair.
Kopimi.

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

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

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

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

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

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

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

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

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

# Data #
values = pandas.Series(lengths)

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

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

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

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

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

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

```

EDIT - an example output:

# Answer
> 5 votes
There are several potential approaches. For example:

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

# Answer
> 5 votes
***You specifically asked about FASTA files, but it's important to always consider read length and quality jointly when assessing high-error long-read data. FASTA files do not provide the quality. This information will help you determine how successful the run was, how many reads were 'high quality', etc.***

I originally posted a full answer here, suggesting pauvre, but I decided it was a little off topic since you seem to only have the FASTA files.

I recommend generating FASTQ files, but I'm not sure whether you have the original base-called fast5 files. If so, generate FASTQs using `poretools` as follows (poretools doc for generating FASTQ files):

```
poretools fastq fast5/

```

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

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

# Answer
> 4 votes
bioawk could be reasonably efficient for this kind of task.

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

```

The `-c fastx` tells the program to parse the data as fastq or fasta. This gives access to the different parts of the records as `$name`, `$seq` (and `$qual` in the case of fastq format) in the awk code (bioawk is based on awk, so you can use whatever language features you want from awk).

Between the single quotes come a series of `<condition> {<action>}` blocks.

The first one has no `<condition>` part, which mean it is executed for each record. Here, it updates the lengths counts in a table which I named "histo". `length` is a predefined function in awk.

In the second block, the `END` condition means we want it to be executed after all the input has been processed. The action part consists in looping over the recorded length values and print them together with the associated count.

The output is piped to `sort -n` in order to sort the results numerically.

On my workstation, the above code took 20 seconds to execute for a 1.2G fasta file.

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

---
Tags: fasta, nanopore
---