title: Genome assembly from error-prone reads

text:
# Question
Title: Genome assembly from error-prone reads
I understand how to assemble genome from error-free reads. I implemented like this:

> Construct directed overlap graph with reads as vertices and edges as maximum overlap between two vertices. Edges represent the length of overlapping maximum substring of one read suffix with other read prefix. Find the path that visits all vertices only once in overlap graph. With small graph we can use greedy approach : Starting from first pattern always choose egde of maximum weight and store shortest common superstring spelled by this path in resulting string. Resulting string spelled by path is shortest common superstring of all reads, that is assembled genome.

For example

```
AACTAG
  CTAGAT

```

Gives the `AACTAGAT` as a result

Now let's say that two error-prone reads have at most one error (mismatch):

```
AACTAG
  CTCGAT

```

As you can see second read has one substitution in third character from **'A'** to **'C'**.

So how assembled genome should look like for the error-prone reads? Do I take **A** or **C** to resulting string? Or I'm missing something?

Note: I used short reads for example because of simplicity. I work with much larger number of reads which are 100 nucleotides long and have at most one substitution of character. So that's a big overlap graph. I know I can use De Bruijn graph, but can I use some greedy approach (like in error-free example)?

# Answer
It's common for assemblers to have a coverage threshold to make sure that there are enough reads across the same region to properly assemble a contig.

You can't work out an ideal solution from a situation where the consensus sequence is ambiguous. If your overlap were representative (e.g. a heterozygous variant in a diploid genome, based on a few tens of reads of 100bp, rather than 6bp), then the resulting assembly would have a bubble in the assembly graph at the variant point.

```
        A
       / \
A-A-C-T   G-A-T
       \ /
        C

```

A conservative assembler might produce two assembled contigs corresponding to these two paths, AACTAGAT and AACTCGAT. A greedy assembler would just pick one path and use that. Subsequent assembly steps might attempt to search for bubbles and annotate the assembled contigs with variant locations.
> 1 votes

# Answer
For short reads, the typical and the most widely used solution is to correct away sequencing errors before assembly. You can correct errors with k-mer spectrum, a trie or multi-alignment. There are many papers on this topic.

Error correction alone won't fix all sequencing errors. Remaining errors may lead to bubbles and tips in the overlap/de Bruijn graph. You can identify them by inspecting the local topology of assembly graph.
> 2 votes

---
Tags: genome, assembly, reads
---


# Question
Title: Is there a database of disease prevalence?
I want to collect data on how common different diseases are in different populations (or, at least, globally). Ideally, it would give me a way of querying the database with a disease name, and would return the number of cases per N population.

My main objective is to have a general idea of how prevalent a disease is, I am not trying to track current outbreaks, so I am not too bothered if the information is a little out of date, as long as the database is still maintained and not abandoned.

I have found FIDD, but it seems to be no longer maintained. Are there any actively maintained, similar tools that have the information I am after?

# Answer
> 5 votes
For rare diseases, you can use Orphanet. You can download an xml file with prevalence from orphadata.

There is an OWL version of Orphanet called ORDO, you can browse it on the EBI Ontology Lookup Service. See for example Arrhythmogenic right ventricular cardiomyopathy and point prevalence relationships in the lower right box. There are a variety of programmatic ways to access this, including the OLS web API, SPARQL queries, or downloading the OWL and using an ontology API.

# Answer
> 5 votes
Comment as an answer, as requested.

The Incidence and Prevalance Database (IPD - subscription required).

---
Tags: database, public-databases
---


# Question
Title: Counting number of possible alignments between two DNA sequences using python
I was looking BioPython to calculate the number of possible alignments between two sequences. Let's say there are two input sequences with the length of *m* and *n*. The program will count the possible number of alignments based on the length of two sequences.

Is it possible to do that with BioPython? If not, how can I do that?

# Answer
> 10 votes
A closed-form solution is offered in *An exact formula for the number of alignments between two DNA sequences* by Torres, Cabada, and Nieto:

$$f(m,n)=\sum\_{k=0}^{min(m,n)}2^{k}\binom{m}{k}\binom{n}{k}$$

If this solution seems reasonable, you could calculate this without BioPython, but with a simple `xrange` loop, the power operator, and `scipy.special.binom`:

```
#!/usr/bin/env python

import sys
import scipy.special

def alignments(m, n):
    if m < 0 or n < 0:
        raise ValueError('m and n should be non-negative')
    s = 0
    b = min(m, n) + 1
    for k in xrange(0, b):
        s += (2 ** k) * scipy.special.binom(m, k) * scipy.special.binom(n, k)
    return s

def main():
    m = 4
    n = 2
    a = alignments(m, n)
    assert (m == 4 and n == 2 and a == 41), "something went wrong!"
    sys.stdout.write("alignments(%d, %d) = %d\n" % (m, n, a))

if __name__ == "__main__":
    main()

```

I don't know what your `m` and `n` are, but keep in mind that this solution grows quickly. Also, I did some digging and scipy approximates the calculation of the binomial coefficient, which is probably going to be an issue for large `m` and `n`.

If you're working with large values, you may want to look into using the Python `long` type and consider the accuracy of the answer. The Stack Overflow link offers alternatives to scipy that generate an exact calculation of the binomial coefficient, at the cost of a longer runtime.

---
Tags: python, sequence-alignment
---


# Question
Title: How to retrieve InChI key for KEGG compound?
I would like to retrieve the InChI representation for a given KEGG compound but I fail to find a direct solution for this (preferably doing it via bioservices).

One could do it via ChEBI like this:

```
from bioservices import *

kegg_con = KEGG()

kegg_entry = kegg_con.parse(kegg_con.get('C00047'))

chebi_con = ChEBI()

chebi_entry = chebi_con.getCompleteEntity('CHEBI:' + kegg_entry['DBLINKS']['ChEBI'])

print chebi_entry.smiles
print chebi_entry.inchi
print chebi_entry.inchiKey

```

which will print

```
NCCCC[C@H](N)C(O)=O
InChI=1S/C6H14N2O2/c7-4-2-1-3-5(8)6(9)10/h5H,1-4,7-8H2,(H,9,10)/t5-/m0/s1
KDXKERNSBIXSRK-YFKPBYRVSA-N

```

Is there a direct way which does not require to go via `ChEBI`?

# Answer
I'm not aware of a functionality in the KEGG service that would allow you to do that directly. I believe your solution (using e.g. ChEBI) is the correct one.

You may have other useful service such as unichem, chemspider (you would need a login though for that one) could help as well.

**Disclaimer**: I'm the bioservices main author
> 2 votes

---
Tags: database, python, kegg, data-retrieval
---


title: error in heatmap using R

text:
# Question
Title: error in heatmap using R
I am learning about generating heatmaps using R. For this example I am using randomly generated data to test the process, but I encountered what appears to be an error or perhaps there is something am missing. Specifically, the values of Npy2r in groups A & B are the same, yet the heat map colors them differently. I am trying to understand why. Am including the R code I am using to generate the heat map

```
heatmap.2(x, Rowv=F, Colv = F, scale="col", col=greenred(3), key=TRUE, symkey=FALSE,
      density.info="none",cexRow=1,cexCol=1,margins=c(6,11),  trace="none",srtCol=45)

```

# Answer
You have `scale="col"` in the code. What you are plotting is the z-score calculated based on the value distribution per column. Try changing it to `"row"` or `"none"`.
> 8 votes

---
Tags: r, differential-expression, microarray
---


# Question
Title: Can exons be located outside of the coding sequence?
I have a gff file like this (I edited the name):

```
scaffold_x  source  exon    2987526 2987805 .   -   .   name "foobar";transcriptId 68892
scaffold_x  source  CDS     2987526 2987805 .   -   1   name "foobar";proteinId 68892;exonNumber 5
scaffold_x  source  exon    2988610 2988745 .   -   .   name "foobar";transcriptId 68892
scaffold_x  source  CDS     2988610 2988745 .   -   2   name "foobar";proteinId 68892;exonNumber 4
scaffold_x  source  exon    2989157 2989645 .   -   .   name "foobar";transcriptId 68892
scaffold_x  source  CDS     2989157 2989645 .   -   2   name "foobar";proteinId 68892;exonNumber 3
scaffold_x  source  exon    2990903 2991541 .   -   .   name "foobar";transcriptId 68892
scaffold_x  source  CDS     2990903 2991541 .   -   2   name "foobar";proteinId 68892;exonNumber 2
scaffold_x  source  stop_codon  2990903 2990905 .   -   0   name "foobar"
scaffold_x  source  exon    2991587 2992430 .   -   .   name "foobar";transcriptId 68892
scaffold_x  source  CDS     2991587 2992430 .   -   0   name "foobar";proteinId 68892;exonNumber 1
scaffold_x  source  start_codon 2992428 2992430 .   -   0   name "foobar"

```

So the `start_codon` is `2992428 2992430` and the `stop_codon` is `2990903 2990905`; I would therefore expect all the exons within this range, however, e.g. exon 5 is found in `2987526 2987805` which is clearly outside this range.

Is there something wrong with the gff file itself or is there a plausible explanation for that?

EDIT:

In the corresponding protein fasta I do find an amino acid sequence associated with this prtein ID and name:

```
>68892|foobar
MHTGD...

```

# Answer
Yes, of course. Exons are not limited to the protein coding regions. Many UTRs are in exons. In fact, you even have various cases of UTRs being multiple exons, and being spliced.

What is strange in your file is not so much that you have exons beyond the stop codon, but that you also have them marked as CDS (coding sequence). That isn't possible, no. While there can indeed be exons in the UTRs, those are not coding and shouldn't be called CDS. That's probably just a minor nomenclature issue though and what you are looking at is a spliced UTR.

So, the exons (or the parts of the exons) that fall between the start and stop codons will be translated, while those that fall outside that regions will not be (UTRs).

For example, this is the exonic structure shown for ENST00000617185, one of the transcripts of the human P53 gene:

The boxes are exons and the lines are introns. The colored boxes are protein coding exons, while the empty ones are UTR exons. Note how they are also spliced and that this particular transcript has 3 non-coding exons on the 5' end and one on the 3' end.
> 11 votes

# Answer
Those are the untranslated regions (UTRs). All mRNAs have a 5' UTR and a 3' UTR. These give the ribosome something to grab onto and often contain important regulatory sites such as miRNA target sites.
> 8 votes

# Answer
Although the explanation that UTR (untranscribed) regions can consist of multiple exons covers most situations, I think it is good to mention stop codon readthrough.

This would result in a CDS (coding sequence) part of your mRNA (messenger RNA) *after* a stop codon.

However, following the gff3 format specification, this should probably be encoded differently, using SO (sequence ontology) terms like SO:0000697 (gene with stop codon read through) and SO:0000883 (stop codon readthrough).
> 2 votes

---
Tags: gene, gff3, molecular-genetics, exon
---


# Question
Title: How to provide a GFF/GTF file to Chromonomer
I am trying to use Chromonomer. It should be possible (according to the software page and a paper) to add a GFF/GTF file to the analysis.

Alas, neither the built-in help nor some light source code browsing were enough to find a way of passing the file to the software.

How can I use a GFF or a GTF file with Chromonomer?

# Answer
I think the part of the code that should allow you to do this is in `chromomer.cc`.

Lines 2857-2983 form the function `parse_command_line`. Nowhere is there any reference to anything resembling parsing a GFF file, so I think the documentation is not in line with the code at this point.

For reference, I looked at version 1.07. It would be helpful if the code were on github.
> 1 votes

---
Tags: gff3, gtf
---


# Question
Title: Finding out proteins that have a similar domain to a specific toxin
I am trying to find out that which DNA sequences(in the below) have a similar domain to a specific toxin which is *Clostridium botulinum*

You can find the fragment of the DNA sequence that belong to that toxin gene below:

> gi|27867582 (fragment of the known Clostridium botuninum toxin gene) GTGAATCAGCACCTGGACTTTCAGATGAAAAATTAAATTTAACTATCCAAAATGATGCTT ATATACCAAAATATGATTCTAATGGAACAAGTGATATAGAACAACATGATGTTAATGAAC TTAATGTATTTTTCTATTTAGATGCACAGAAAGTGCCCGAAGGTGAAAATAATGTCAATC TCACCTCTTCAATTGATACAGCATTATTAGAACAACCTAAAATATATACATTTTTTTCAT CAGAATTTATTAATAATGTCAATAAACCTGTGCAAGCAGC

I ran a blastx job over the above sequence which translates nucleotide into protein. Here's the result.

The result tells me that. The particular DNA sequence listed above, belongs to *botulinum neurotoxin type E \[Clostridium botulinum\]*

Besides, I am given two DNA sequences:

```
>Sequence A
TCTATCAAGTAGATTATTAAATACTACTGCACAAAATGATTCTTACGTTCCAAAATATGA
TTCTAATGGTACAAGTGAAATAAAGGAATATACTGTTGATAAACTAAATGTATTTTTCTA
TTTATATGCACAAAAAGCTCCTGAAGGTGAAAGTGCTATAAGTTTAACTTCTTCAGTTAA
TACAGCATTATTAGATGCATCTAAAGTTTATACGTTTTTTTCTTCAGATTTTATTAATAC

> Sequence B
TCCTGGCTCAGGACGAACGCTGGCGGCGTGTGCTTAACACATGCAAGTCGAGCGATGAAG
CTTCCTTCGGGAAGTGGATTAGCGGCGGACGGGTGAGTAACACGTGGGTAACCTGCCTCA
AAGTGGGGGATAGCCTTCCGAAAGGAAGATTAATACCGCATAACATAAGAGAATCGCATG
ATTTTCTTATCAAAGATTTATTGCTTTGAGATGGACCCGCGGCGCATTAGCTAGTTGGTA

```

Again, I translated them into proteins. The result are listed respectively(A and B):

https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Get&RID=Z423E71C015

https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Get&RID=Z42AKVWC015

The first job tells that, Sequence A belongs to *botulinum neurotoxin type F \[Clostridium baratii\]* which means that similar domain with the above *botulinum neurotoxin*\>

Second job tells me that SequenceB belongs to *unknown protein \[Clostridium botulinum NCTC 2916\]*. So that means it is not in the same domain with the above botulinum neurotoxin right?

# Answer
A protein domain is a conserved part of a given protein sequence and (tertiary) structure that can evolve, function, and exist independently of the rest of the protein chain.

This is different from global protein similarity, so similarity to neurotoxin type F doesn't necessarily mean that sequence A has similar domain to your known toxin gene. Your blast results do include conserved domain information though.
> 1 votes

---
Tags: sequence-alignment, homework
---


title: IGV jump to rsID location

text:
# Question
Title: IGV jump to rsID location
Is it possible to load a list of rsID identifiers (e.g. rs66823261) to IGV to then be able to jump to the location of each of them when copy+pasting the rsID on the locus text box?

For human and other genomes, there is currently the feature of typing in a gene name, then jumping to the gene location.

I wonder if there is a way of achieving the same for rsIDs.

EDIT: thanks for the answers so far. I downloaded the Infinium GSA bed file from Illumina's website (https://support.illumina.com/downloads/infinium-global-screening-array-v1-0-support-files.html), loaded on an IGV session, which asked to create an index. I did a liftover of the file to hg38 and sorted with bedtools sort. Then I went to the text box, introduced one of the rs ids:

rs6569648

Got a positive result. Thank you all!

# Answer
Yes, if you load the list you are interested in as a non-indexed "BED" file then you can jump to an rsID or any other name in the BED file. You'll need a targeted list, dbSNP for example is much too large to load in this way.
> 2 votes

# Answer
you can use mysql ucsc to get the position of the SNP and then ask IGV to take some picures: see https://software.broadinstitute.org/software/igv/batch

```
$ echo -e 'rs66823261\nrs25' | awk '{printf("select concat(chrom,\":\",chromStart-10,\"-\",chromEnd+10) from  snp150 where name=\"%s\";\n",$1);}' | mysql --user=genome --host=genome-mysql.soe.ucsc.edu -A -D hg38 -N | awk '{printf("goto %s\nsnapshot\n",$1);}'

goto chr8:170681-170702
snapshot
goto chr7:11584131-11584152
snapshot

```
> 0 votes

# Answer
If you want all rsIDs (complete dbSnp database) then Pierre's answer is the best option currently. Any bed file with all the rsIDs would be too big to load. For a targeted list of, say, 100,000 or fewer snps you can use a bed file.

BTW igv.js already does the equivalent of Pierre's suggestion, this will be ported to IGV desktop in the near future.
> 0 votes

---
Tags: snp, igv
---


# Question
Title: Where are .motif files from homer knownResults?
I have been using homer's `findMotifsGenome.pl`, but with my new version (v4.9.1) of `homer` I don't get `.motif` files in the `knownResults` folder. I do get them in the `homerResults` (*de novo*) folder, though. With my previous version of `homer` I did get them and used them to find which peaks are in which motif. Which is what I would like to do

Does anyone know why there are no `.motif` files? I have only `.svg` files. I have seen that for the `homerResults` the `.motif` files are concatenated already into a file called `homerMotifs.all.motifs`, but no such file for `knownResults`.

PS. I tried to call Chuck for help, but he was not answering his phone...

PS2. I will ask this question on Biostars to hopefully get more response than here.

# Answer
> 3 votes
I have provided a solution in Biostars:

If you edit the file `findKnownMotifs.pl` and move the call to `printMotif` to come before testing `$svgFlag` you will find that `.motif` files are printed again.

In the following I have commented out the call from where it need not be, and introduced it as the first line, right before the `if` block:

```
   printMotif($_, $directory . "/$knownDir/$matrixFile"); # this should be done regardless of $svgFlag
    if ($svgFlag) {
        my $str = HomerSVGLogo::motif2svg($_,$bitFlag,$rnaFlag2);
            open SVG, ">$directory/$knownDir/$logoFile.svg";
            print SVG $str;
            close SVG;
            print MAIN "$str\n";
    } else {
  #printMotif($_, $directory . "/$knownDir/$matrixFile"); # this should be done regardless of $svgFlag
            `profile2seq.pl "$directory/$knownDir/$matrixFile" 100 > "$tmpFile"`;
            my $width = $widthFactor * $_->{'len'};
            `seqlogo -a $bits -f "$tmpFile" -F PNG -S -c -o "$directory/$knownDir/$logoFile" -h $height -w $width`;
  `seqlogo -a $bits -f "$tmpFile" -F PDF -S -c -o "$directory/$knownDir/$logoFile" -h $height -w $width 2> /dev/null`;
  print MAIN "<IMG src=\"$knownDir/$logoFile.png\"/>";
}

```

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

I have reported this to the author of this code.

# Answer
> 1 votes
The relevant help page of homer seems to be this one.

On section 13 it says:

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

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

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

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

---
Tags: chip-seq, motifs, homer, atac-seq
---


# Question
Title: How to access organism IDS, taxon nodes from UniProt?
I try to find organism IDs used by UniProt (e.g. 9606 for "Homo sapiens") and information about the 'taxonomic kingdom' this organism belongs to. The only source on UniProt I could find is this file which contains all required information. So in the Python world it would be a dictionary of the form:

```
{org_id: {'name': 'org_name', 'kingdom': 'org_kingdom'}}

```

or a dataframe.

However, the file from above is not particular nice to parse due to its inconsistent column format (not impossible to parse though), so I am wondering whether there is an alternative source to retrieve this information or whether there is maybe even a webservice available that takes care of this parsing.

# Answer
That file is actually trivial to parse. You only care about lines that have `N=$species`, so you can simply do:

```
sed -En 's/.* ([0-9]+): N=(.*)/\1\t\2/p' speclist

```

That will return a tab separated list of taxID and species:

```
$ sed -En 's/.* ([0-9]+): N=(.*)/\1\t\2/p' speclist | head
648330  Aedes albopictus densovirus (isolate Boublik/1994)
10804   Adeno-associated virus 2
648242  Adeno-associated virus 2 (isolate Srivastava/1982)
118452  Abacion magnum
72259   Abaeis nicippe
102642  Abax parallelepipedus
392897  Abalistes stellaris
75332   Abbottina rivularis
515833  Abdopus aculeatus
56673   Antarctic bacterium DS2-3R

```

So to get human, for example, you could do:

```
$ sed -En 's/.* ([0-9]+): N=(.*)/\1\t\2/p' speclist | grep sapiens
63221   Homo sapiens neanderthalensis
9606    Homo sapiens

```

And, for exact matches:

```
$ sed -En 's/.* ([0-9]+): N=(.*)/\1\t\2/p' speclist |
    awk -F'\t' '$2=="Homo sapiens"'
9606    Homo sapiens

```

If you also want to keep the taxon node (kingdom), you can use:

```
sed -En 's/.* ([A-Z])  *([0-9]+): N=(.*)/\1\t\2\t\3/p' speclist

```

Which prints:

```
$ sed -En 's/.* ([A-Z])  *([0-9]+): N=(.*)/\1\t\2\t\3/p' speclist  | head
V   648330  Aedes albopictus densovirus (isolate Boublik/1994)
V   10804   Adeno-associated virus 2
V   648242  Adeno-associated virus 2 (isolate Srivastava/1982)
E   118452  Abacion magnum
E   72259   Abaeis nicippe
E   102642  Abax parallelepipedus
E   392897  Abalistes stellaris
E   75332   Abbottina rivularis
E   515833  Abdopus aculeatus
B   56673   Antarctic bacterium DS2-3R

```

---

The `sed` command might appear a little daunting if you don't know `sed`, but it's actually quite simple. The `-E` flag enables extended regular expressions which give us a simpler syntax. The `-n` suppresses normal output, so that `sed` won't actually print anything unless explicitly told to.

The `s/from/to/` is the substitution operator and will replace `from` with `to`. Here, we are replacing zero or more characters (`.*`) until a space followed by one or more numbers (`[0-9]+`), a `:`, another space, then the string `N=` and then everything else (`.*`). The parentheses let us *capture* the matched strings, so that we can then refer to the 1st captured string as `\1`, the second as `\2` etc. We are therefore replacing the entire line (because the regex matches the entire line) with the 1st pattern captured (the taxID, `\1`), a tab (`\t`) and the species name (`\2`). Finally, the `p` at the end means "print the resulting line if the substitution was successful", so only lines matching the desired pattern are actually printed.
> 1 votes

# Answer
These are NCBI taxonomy ID's rather than anything specific to UniProt. One way to parse them would be using BioPython via `Bio.Entrez`, for example:

```
from Bio import Entrez

def get_taxon_name_from_NCBI_id(id_):
    Entrez.email = 'name@provider.com'
    handle = Entrez.efetch(db='taxonomy', id= id_, retmode='xml')
    records = Entrez.read(handle)[0]
    return records['ScientificName'], records['TaxId']

print(get_taxon_name_from_NCBI_id(9606))
# ('Homo sapiens', '9606')

```

Of course you could query Entrez directly: https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=taxonomy&id=9606

Or even more basically do a manual search of NCBI taxonomy: https://www.ncbi.nlm.nih.gov/taxonomy/?term=9606
> 4 votes

---
Tags: database, data-retrieval, uniprot, webservice
---


title: Find Patterns in Cluster

text:
# Question
Title: Find Patterns in Cluster
I have a heatmap and I would like to find some rectangles.

I have already used clustermap. But here, I can not calculate these rectangles. The order of the data should not be changed.

This Code is not working:

```
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
tips2 = pd.read_csv('heatmapExport.csv')
sns.set(color_codes=True)
g = sns.clustermap(tips2, cbar="false")
plt.show()

```

Does somebody have an idea?

Thanks a lot

Here is an example of what I want:

# Answer
> 3 votes
Based on your description I think you should have a look at a technique called 'biclustering'.

The example on this page defines the goal of this technique as 'Finding subgroups of rows and columns which are as similar as possible to each other and as different as possible to the rest.'

Since your examples are python-based, you could check out scikit-learn's implementations of biclustering.

---
Tags: python, clustering
---


title: What is deep sequencing?

text:
# Question
Title: What is deep sequencing?
People talk about deep sequencing. Is there any way to calculate how deep the sequencing is ? What should be the optimum depth to get reliable data ?

I am doing whole genome sequencing of a virus genome which is 10 kb long. I got 80000 reads from Illumina sequencing. Is there any way to tell about my sequencing depth ?

# Answer
I found a post useful for this topic.

It explains the difference of coverage and depth. It also has a useful explanation on how to calculate coverage and depth.

Here is a copy of what the link says just encase the post is removed:

> **Depth of coverage**
>
> *How strong is a genome "covered" by sequenced fragments (short reads)?*
>
> Per-base coverage is the average number of times a base of a genome is sequenced. The coverage depth of a genome is calculated as the number of bases of all short reads that match a genome divided by the length of this genome. It is often expressed as 1X, 2X, 3X,... (1, 2, or, 3 times coverage).
>
> **Breadth of coverage**
>
> *How much of a genome is "covered" by short reads? Are there regions that are not covered, even not by a single read?*
>
> Breadth of coverage is the percentage of bases of a reference genome that are covered with a certain depth. For example: 90% of a genome is covered at 1X depth; and still 70% is covered at 5X depth.

# SAMtools: get breadth of coverage

Get breadth of coverage of a reference genome.

How many bases of my reference genome are covered by my short reads sample at a level of 5X or higher?

```
 # set shell variables to your file-names

 REF_GENOME_FILE=my_ref_genome.fna
 WGS_SAMPLE=my_short_read_sample.fastq
 MIN_COVERAGE_DEPTH=5

```
1. Map short reads against reference genome (bowtie2)

```
# create bowtie2 index database (database name: refgenome)

bowtie2-build ${REF_GENOME_FILE} refgenome
ls
refgenome.1.bt2
refgenome.2.bt2
refgenome.3.bt2
refgenome.4.bt2
refgenome.rev.1.bt2
refgenome.rev.2.bt2

# map reads and sort bam file
bowtie2 -x refgenome --no-unal -U ${WGS_SAMPLE} -S - -p 12 | \
samtools view -bS - | \
samtools sort -m 5G - mapping_result_sorted.bam

-U ${WGS_SAMPLE} consider all reads as unpaired reads
-S -             bowtie output in SAM format, write to standard out
-x refgenome     use reference genome bowtie2-database
--no-unal        suppress SAM records for reads that failed to align
-p 12            use 12 processors

```
2. Get breadth of coverage (SAMtools)

```
# create samtools index
samtools index mapping_result_sorted.bam

# get total number of bases covered at MIN_COVERAGE_DEPTH or higher
samtools mpileup mapping_result_sorted.bam | awk -v
X="${MIN_COVERAGE_DEPTH}" '$4>=X' | wc -l
32876

# get length of reference genome
bowtie2-inspect -s refgenome | awk '{ FS = "\t" } ; BEGIN{L=0}; {L=L+$3};
END{print L}'
45678

```

Result: breadth of reference genome coverage total number of covered bases: 32876 (with \>= 5X coverage depth) Depth of coverage (average per-base coverage): 0.719 X (32876 ÷ 45678) (total number of covered bases divided by reference genome length) percent: 71.9% (0.719 × 100) (coverage breadth in percent)
> 11 votes

# Answer
There are several questions in your post I'll try to answer each one:

> Is there any way to calculate how deep the sequencing is ?

See gringer's answer. TLDR: The depth of the sequencing is how many times each position has been sequenced.

> What should be the optimum depth to get reliable data ?

The optimum depth depends on what you want to do with that data. Usually, 30X or 50X is done, but in targeted sequencing you can sequence a lot more reaching 1000X (which would be really deep).

> Is there any way to tell about my sequencing depth ?

If you know the length of the region sequenced and how many bases you have in the reads, you can calculate the general sequencing depth.
> 9 votes

# Answer
Sequencing depth is typically calculated as the *number of total bases sequenced* divided by the *number of bases in the target genome*. An Illumina sequencing run with 2x125 bp reads and 500 million read pairs sequenced would be a sequencing depth of about 40X (assuming my calculations are correct for a 3 billion base genome).

The sequencing depth depends on how long each read is, as well as whether it is a paired-end sequencing run. 80,000 read pairs with each having 125 bp would be a coverage of 2000X for a 10 kbp viral genome (assuming the sample were entirely composed of that virus).

The optimum depth changes depending on the application. A genotyping-by-sequencing analysis, combined with imputation, could work well with 96 samples at a 1X coverage (or possibly less). On the other end of the scale, a cell-free DNA analysis to identify structural variation in metastasing cancer may need a few thousand times coverage of the host DNA.
> 8 votes

---
Tags: ngs, covid-19, sars-cov-2
---


title: Adding entries to bigwig file

text:
# Question
Title: Adding entries to bigwig file
I generate bigwig files using a shell script based on `bedtools genomecov` (to generate a bedgraph from a bam file) `bedmap` (to compute means across 10 bp bins) and `bedGraphToBigWig` to convert the binned bedgraph into bigwig.

Sometimes, the bam file has no data on some chromosomes, this results in a bedgraph file with no entries for the empty chromosomes.

I would like to add the missing entries. How could I do this?

I tried to see what I could do with pyBigWig, opening the bedgraph file in "a" mode, but don't really know how to proceed:

```
# I have the chromosome sizes for my genome of interest already in a list:
chrom_sizes
#[('I', 15072434),
# ('II', 15279421),
# ('III', 13783801),
# ('IV', 17493829),
# ('MtDNA', 13794),
# ('V', 20924180),
# ('X', 17718942)]
bw_test = pyBigWig.open("path/to/bigwig", "a")
bw_test.chroms()
# {'IV': 17493829, 'V': 20924180, 'X': 17718942}
bw_test.addHeader(chrom_sizes)
# RuntimeError: Received an error in bwCreateHdr

```

# Answer
> 3 votes
Presuming you still have your BAM file sitting around somewhere:

```
#!/usr/bin/env python
import pysam
import pyBigWig

bam = pysam.AlignmentFile("alignments.bam")
ifile = pyBigWig.open("your current file.bw")
ofile = pyBigWig.open("your new file.bw", "w")

# Create/Add the header
hdr = []
for chrom in bam.header['SQ']:
    hdr.append((chrom['SN'], chrom['LN']))
bam.close()
ofile.addHeader(hdr)

# Write the bigWig
for _ in hdr:
    try:
        ints = ifile.intervals(_[0])
        starts = [x[0] for x in ints]
        ends = [x[1] for x in ints]
        vals = [x[2] for x in ints]
        ofile.addEntries([_[0]] * len(starts), starts, ends=ends, values=vals)
    except:
        # intervals() will throw an error if the chromosome isn't present
        pass
ofile.close()
ifile.close()

```

The reason you can't append to a bigWig file is that it has a various indices and zoom levels typically built in. You'd have to add to them without overwriting neighboring blocks of data. In theory this might sometimes be theoretically possible, but not typically and I'm not aware of tools that support this.

---
Tags: python, bigwig
---


# Question
Title: Does picard markduplicate toggle PCR duplicate samflag
I have a RNA-seq bam file and there are few reads that are puzzling me.

According to the bam header, this bam file is sorted by coordinate, created using tophat and markduplicate step is not done. But some reads are marked for being duplicate in the samflag. What is worse is when I run picard markduplicate, these reads' pcr duplicate flag is toggled, marking them not a duplicate. Also I manually found the duplicate of this read (identical reads with same start positions and mates start position) so initial marking looks true.

So my questions are:

Any idea why would this happen?

Does Tophat mark reads for being duplicate? ( I don't think so) And does `picard markduplicate` toggle the duplicate flag if the reads are already marked as being duplicates?

Here is how the read looks before and after mark duplicate step.

Before:

```
C0RTF 1187 17 7579880 255 61M10754N40M = 7579927 10902 CTC...
0UNP1 163 17 7579880 255 61M10754N40M = 7579927 10902 CTC...

```

After Markduplicate:

```
C0RTF 163 17 7579880 255 61M10754N40M = 7579927 10902 CTC...
0UNP1 163 17 7579880 255 61M10754N40M = 7579927 10902 CTC...

```

# Answer
Yes, if MarkDuplicates encounters a pair that's marked as a duplicate that it considers (for whatever reason) to not be a duplicate then it will unset the duplicate mark. You can test this yourself by making a small BAM file either with or without duplicate entries but some marked regardless. In the case of actual duplicates, the duplicate flag will be maintained. In the case where there are no actual duplicates then the flag will be removed. At least that's the case in version 2.14.
> 6 votes

---
Tags: rna-seq, samtools, markduplicate
---


# Question
Title: Differential Expression With 2 Treatments
I have a expression data from a small cohort of samples taken at baseline and after 2 independent treatments. I can do differential expression contrasting T1 and T2 or I can contrast T1 vs baseline and T2 vs baseline and look at the differences. What is the difference between these two analyses? Are they both valid? If so what inferences can be drawn from each?

# Answer
If you're interested in looking at the differences between two treatments then you'll end up wanting to do both a direct contrast as well as the individual comparisons to baseline.

The direct contrast will give you the genes actually differentially expressed between the two conditions. In practice, you may want to filter this a bit so you only have genes differentially expressed vs. baseline in at least one condition (e.g., to get rid of genes only slightly higher due to T1 and slightly lower in T2, but not different enough in either case to be DE). Use a pretty lax p-value threshold (e.g., 0.1 or 0.2) for this.
> 4 votes

# Answer
Depending upon the type of treatment used the set of DEGs will change.If the treatments have similar kind of effect you will get a small list (less variable genes will have higher p-val) using a cutoff of p\<0.05. So, it's better to start with control vs treated comparisons then T1 vs T2.

Compare the lists form CTRL vs T1 and CTRL vs T2 you will get the genes that are expressed in both conditions as well as unique to individual treatments.
> 3 votes

---
Tags: differential-expression
---


# Question
Title: Coverage calculation: long reads (RNA-seq)
Say your aim is to calculate the coverage of an RNA-seq experiment generated with long-read sequencing (so, uneven read length).

Up to now, I relied on the Lander/Waterman equation:

$$C = L\*N / G$$

where:

$C$ = final coverage
$G$ = haploid genome (or transcriptome) length
$L$ = read length
$N$ = number of reads

I have two major conceptual issues with this formula:

* $L$ is uneven in this case. Should we use the median/mean read length instead? If so, side question: if the RNA-seq experiment was executed following the size-fractionated protocol (PacBio), is it advisable to merge all the fractions prior to coverage calculation or should it be calculated separately for each fraction?
* Should $N$ consist of the mapped reads to the reference or simply the number of reads as found in the subreads whole dataset?

# Answer
The read length is irrelevant when calculating the mean coverage statistic. It's simply the *total number of bases sequenced* divided by the *target Xome length*.

In the example provided in the question, $L\*N$ is equivalently expressed as $\sum^{N} l\_x$, or the total length of sequenced reads. For equal-length reads, this is easier to calculate as *number of reads* $\*$ *length of reads*, but the full sum is necessary if read lengths are not equal.

For what it's worth, this value, $L\*N$, is also equivalent to the average read length multiplied by the number of reads, i.e.

$$\frac{\sum^{N} l\_x}{N} * N$$

... which makes it a bit more obvious why the total length of sequenced bases is *very slightly* easier to calculate.

In answer to the second part of the question, $N$ should ideally correspond to only the sequences that are from the reference. For an assembly, only the reads that were used to create the assembly should count when considering coverage. When there's no likelihood of contamination (and especially for RNASeq with variable coverage per transcript), a statistic using the total number of reads should be good enough for most purposes that I can think of. High-throughput sequencing is such a fuzzy process that a few tens of millions of bases here or there shouldn't make a difference.
> 6 votes

---
Tags: rna-seq, transcriptome, long-reads, pacbio, coverage
---


# Question
Title: What is this 5-column sequencing format from 2009?
There is high throughput sequencing data here, and I don't know what format it is in.

It was submitted in 2009, and the description says the following:

* Library strategy: ncRNA-Seq
* Library source: transcriptomic
* Library selection: size fractionation
* Instrument model: Illumina Genome Analyzer II
* Description: CIPPNK, tar file of Illumina \*_seq.txt files provided as supplementary file

I got the archive here:

ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM455nnn/GSM455387/suppl/GSM455387%5FWT%5FCIPPNK%5Fseq%5Fs1%2Etar%2Egz

Inside, there are 330 files from `s_1_0001_seq.txt` to `s_1_0330_seq.txt` which are tab-delimited text files where the first column is always `1`, the second has the number found in the file name, then 2 mysterious integers, and then what looks like a read of length 36, with sometimes a dot instead of a letter:

```
$ head s_1_0330_seq.txt
1   330 690 785 TTCCTACATTGTTCCCCCATGCTGTTGGCACCATCA
1   330 44  145 TTTTTATCACGAGTTTTAAATCTGTAGTCACCATCA
1   330 53  141 AATAATGCATAACAAAACGGAATCTGTAGAA.AAA.
1   330 784 461 TAATTGTAGTGATTGATCAATCTGTAGGCACCATCA
1   330 588 634 TATTATGCACATTTTCTAGTTCACTGTAGGCACCAT
1   330 718 678 TTACATGTTTCGGGTAGGAGCCTGTAGGCACCATCA
1   330 635 834 TGTGATCATTAGTTCAAAGCCCCCTGTCGGCACCCT
1   330 494 523 TGAAAATCAAAAATGCTGAACTGTAGGCACCATCAA
1   330 393 783 TTTTTTTTTAAATTTAAAAAAACTGTAGGCACCATC
1   330 48  148 GTTTAACCGTGTAGACGTTGGTTTCTGTAGGCACCA

```

This seems to be some "old" high-throughput sequencing format. I think someone in a 2008 message in sequanswers was dealing with this type of file:

http://seqanswers.com/forums/showpost.php?p=1841&postcount=8

What is this format that seemed so standard back then, that the authors did not give more info than describing the files as "Illumina \*_seq.txt files"? I do not dare ask them for such a trivial question (the given contact is a Nobel prize laureate and is probably too busy to answer random bioinformatics questions).

**In particular, what are the columns 3 and 4, and what do the dots mean?**

# Answer
This is an early Solexa/Illumina sequencer format. The columns are the identifier to location on the flowcell. I believe the "." was the original placeholder for an unread base, which has been replaced by and "N" in current Illumina sequencing output.

From on http://www.crg.eu/en/content/processing-and-analysis-illumina-sequencing-data.

> seq.txt (Gerald)
>
> 4 1 23 1566 ACCGCTCTCGTGCTCGTCGCTGCGTTGAGGCTTGCG
>
> The first 4 columns are an ID that provides location details of the cluster in the flow cell, followed by the read sequence. Seq.txt file size typically ranges between 250 Megabytes and 1 Gigabyte. One file is generated per lane.
> 7 votes

---
Tags: file-formats, hts
---


# Question
Title: Retrieve abstract/summary for NCBI bookshelf via Entrez
I want to retrieve abstracts/summaries for NCBI bookshelf entries, eg: "NBK1440"

Docs say the dbname is "books" and efetch guide says a rettype of "docsum" works for all databases. However when I make the query:

```
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=books&id=NBK1440&rettype=docsum

```

returns:

```
<ERROR>UID=1440: cannot get document summary</ERROR>

```

The default query is useless (just echos back the ID):

```
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=books&id=NBK1440&retmode=text

```

returns

```
1440

```

Picking another 3:

```
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=books&id=NBK1391,NBK410087,NBK1391&rettype=docsum

```

Fails the same way:

```
<ERROR>UID=1391: cannot get document summary</ERROR>
<ERROR>UID=410087: cannot get document summary</ERROR>
<ERROR>UID=1391: cannot get document summary</ERROR>

```

The same goes for esummary:

```
https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=books&id=1440

```

Dies on NBK1440 and removing NBK gives:

```
<ERROR>UID=1440: cannot get document summary</ERROR>

```

The values come from the (9Mb) ClinVar citations file: https://ftp.ncbi.nlm.nih.gov/pub/clinvar/tab\_delimited/var\_citations.txt \- the first is from tihs line:

```
15048   9   1800562     NCBIBookShelf   NBK1440

```

**ANSWER**

Thanks to arupgsh - yes the ClinVar citation\_id was a RID for bookshelf not an ID like it is for Pubmed and PMC. FYI my solution was:

```
def get_summary_for_bookshelf_rid(bookshelf_rid):
    handle = Entrez.esearch(db="books", term=bookshelf_rid)
    search_results = Entrez.read(handle)
    handle = Entrez.esummary(db="books", id=','.join(search_results['IdList']))
    results = Entrez.read(handle)

    for r in results:
        if r["RID"] == bookshelf_rid:
            return r

    return None

get_summary_for_bookshelf_rid("NBK1440")

```

returns:

```
{'PubDate': '2000/04/03 00:00', 'Title': 'HFE-Associated Hereditary Hemochromatosis', 'Text': '', 'RType': 'chapter', u'Id': '1475938', u'Item': [], 'Book': 'gene', 'Parents': '', 'RID': 'NBK1440', 'ID': 'gene/chapter/hemochromatosis/PMC'}

```

# Answer
Entrez requires the unique identifiers(UID) for fetching related info. The id you are using in the query is RID. I guess that why you are getting `<ERROR>UID=1440: cannot get document summary</ERROR>`.

List of example uids using the term "ncbi+blast":

https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=books&term=ncbi+blast

Summary of the first item:

https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=books&id=4459608

Using the Entrez python package:

```
from Bio import Entrez
Entrez.email = "youremail@domain.com"
handle = Entrez.esummary(db="books", id="4459608,4456634", retmode="xml")
records = Entrez.parse(handle)
for record in records:
     print(record['Title'])
handle.close()

```

I'm not sure about the availability of book abstracts.
> 5 votes

# Answer
As arupgsh says, you need to use esearch to get a list of unique identifiers before using efetch to retrieve info about each result. I think the easiest way to do this is to use Entrez Direct, which allows you to simply pipe esearch output to efetch:

```
esearch -db books -query NBK1440 | efetch -format docsum

```

or

```
esearch -db books -query NBK1440 | esummary

```
> 3 votes

---
Tags: api, data-retrieval
---


title: Biopython code refactoring

text:
# Question
Title: Biopython code refactoring
I want to make function get\_abstract, but i don't know what this can look like if the function must return idlist.

```
from Bio import Entrez, Medline, SeqIO

list_of_ids = []
Entrez.email = 'ski89@gmail.com'

handle_one = Entrez.esearch(db="pubmed", term='e coli',
                           retmode="text", retmax=100000)
record = Entrez.read(handle_one)
idlist = record["IdList"]
for ID in idlist:
    list_of_ids.append(ID)

    if len(list_of_ids) == 100000:
        handle_one = Entrez.esearch(db="pubmed", term='e coli',
                           retmode="text", retstart=100000, retmax=100000)
        record = Entrez.read(handle_one)
        idlist = record["IdList"]
        for ID in idlist:
            list_of_ids.append(ID)

            if len(list_of_ids) == 200000:
                handle_one = Entrez.esearch(db="pubmed", term='e coli',
                           retmode="text", retstart=200000, retmax=100000)
                record = Entrez.read(handle_one)
                idlist = record["IdList"]
                for ID in idlist:
                    list_of_ids.append(ID)

```

so the args will be the start and TERM:

```
get_abstracts(start,TERM):
handle_one = Entrez.esearch(db="pubmed", term=TERM,
                               retmode="text", retstart=start, retmax=100000)
                    record = Entrez.read(handle_one)
                    idlist = record["IdList"]
                    for ID in idlist:
                        list_of_ids.append(ID)

```

Can someone help refactor this? cheers

# Answer
Here is a sample of code that should get you started.

I created a separate `get_abstracts` function. Then that function get called repetitively until less than 100K ids are returned. If you wanted to get fancy this whole process could be called as a recursive function but that seemed overly complex.

```
BUFFER_SIZE = 100000

def get_abstracts(start, term):
    """Get N number of records from Entrez
    :returns list: List of ids returned from search
    """
    handle_one = Entrez.esearch(db="pubmed", term=term,
                       retmode="text", retstart=start, retmax=BUFFER_SIZE)
    record = Entrez.read(handle_one)
    idlist = record["IdList"]
    return idlist

list_of_ids = []
start_at = 0

# Iterate the search returning 100000
while true:
    ids = get_abstracts(start_at, 'e coli')

    list_of_ids.extend(ids)

    # Break if less than 100K results were returned
    if len(ids) < BUFFER_SIZE:
        break
    start_at += BUFFER_SIZE

# Do something with list_of_ids

```
> 4 votes

---
Tags: python
---


title: Simple Affy Differential Expression

text:
# Question
Title: Simple Affy Differential Expression
I think I have a fairly easy-to-solve problem.

I have an expression table generated from an Affy set. It has one row of gene names (already symbols) and then 180 sample rows. 181 rows total.

So I made another excel file with one row with three different phenotypes I want to compare. That row has a label ("classes") and 180 names. 181 rows total.

In R, I can import each of these

```
dat <- read.csv("expression.csv")
dat <- data.matrix(dat)
pheno <- read.csv("pheno.csv")
pheno <- data.matrix(pheno)

```

I can then fit them

```
fit <- lmFit(dat,pheno)

```

And ask for an eBayes from limma

```
fit <- eBayes(fit)

```

Now, when I view the topTable,

```
table <- topTable(fit)

```

My first column has numbers (like "834","1142"), when I'd expect to have gene names there (like "ACTR", "CHRNA1").

What am I doing wrong? Or, rather, am I doing anything right?

Thanks!

# Answer
> 1 votes
What do the first few lines of your lmFit input variables look like (i.e. `dat` and `pheno`)? Are the rows numbered, or do they use the actual symbols?

By default, I don't think `read.csv` expects a file to have row numbers (even though `write.csv` puts them in by default). Try adding a `row.names` parameter to the expression matrix read call:

```
dat <- read.csv("expression.csv", row.names=1)

```

---
Tags: rna-seq, bioconductor, differential-expression, microarray, rna
---


# Question
Title: sum of products of mapq and mapped bases for each read in a from a BAM file
Given a BAM file, I'd like to calculate the sum, over all reads, of the mapping quality and the number of mapped bases (i.e. number of M's in the CIGAR string).

For example, given two reads like this:

```
read1: mapq=40, mapped_bases=10
read2: mapq=20, mapped_bases=50

```

The number I'm looking for is:

```
(40*10) + (20*50) = 1400

```

One option here is to stream the file and get the two quantities line-by-line. But this seems likely to be slow. Thoughts welcome.

# Answer
BAM files are stored line by line, so your only option is to process it that way. In python with `pysam`:

```
import pysam

bam = pysam.AlignmentFile("something.bam")
tot = 0
for read in bam.fetch():
    if read.is_unmapped:
        continue
    tl = 0
    for op, l in read.cigartuples:
        if op in [0, 7, 8]:
            tl += l
    tot += tl * read.mapping_quality
bam.close()
print("the total is {}".format(tot))

```
> 3 votes

# Answer
One possibility uses bamutils (http://ngsutils.org/modules/bamutils/removeclipping/) to remove the clipped bases, then awk to get the sum of the products, something like:

```
bamutils removeclipping val.bam val_rc.bam
samtools view val_rc.bam | awk '{print length($10)*$5}' | awk '{sum+=$0} END { print "mapq*len sum =",sum}'

```

the bamutils step is rather slow, but it will do the job.
> 0 votes

---
Tags: sam, samtools, mapq
---


# Question
Title: How to release an R package with several bed files?
I'm currently creating an R package, and my scripts require before any analysis is done that 1-2 bed files be loaded. Normally, I would run the scripts with the following:

```
library(data.table)
session_bed = read.table("/path/name/to/bedfile.bed")

```

and then refer to `session_bed` throughout the analysis.

However, in order to release this as an R package, I believe that these files must be downloaded at the time on the package installation.

(1) Normally, how are bed files distributed in R packages? Are these gzipped? Where are they placed?

According to http://r-pkgs.had.co.nz/data.html, all .bed files should be `inst/extdata` I guess. However, it's not clear to me how to "link" to this data within the R package.

How have other packages dealt with this? (Is there a size limit?)

(2) In this case, should the function loading the bed file be executed about loading the R library, `library(package_name)`? Is there more appropriate option?

# Answer
> 7 votes
You can store it in `data/` subfolder. If one of your functions need this bed file, you can include importing it in your function.

Here are some examples how to get the data from `data/` subfolder.

There are probably more ways to do this, this is just one way.

# Answer
> 6 votes
As benn said you can store them in `data/` subfolder of your package. There is no size limit (AFAIK) but if you plan to publish your package in some repositories they might have limits or other restrictions. Bioconductor requires to divide the data required for a software package at certain threshold, but for small amount of data it is fine.

You can find them via `system.file` like: `read.csv(system.file("data/mydata.csv", package = "mypackage"))`. Of course your bed files will need the apropiate function to load them.

You can load when your library is called by modifying the `local` variable, or using `.First` (I am not so sure of the last one), but this pollutes the environment and sometimes if preferably to use `data` function. You might need to convert your files to `.Rdata` file.

---
Tags: r, bed
---


# Question
Title: How to resolve in snakemake error : "Target rules may not contain wildcards."
I would like to do easily reproducible analysis using publicly available data from NCBI, so I have chosen a `snakemake`.

I would like to write a single rule, that would be able to download any genome given a species code name and separated table of species and their NCBI IDs. So I wrote a script `scripts/download.sh` that takes a species code and download the genome to `data/<species_code>/genome.fa.gz`. The script internally reads the table `tables/download_table.tsv`, where are corresponding species code names and NCBI IDs.

So I tried to do a `snakemake` like this :

```
species='Cbir Avag Fcan Lcla Dcor Dpac Pdav Psp62 Psp79 Minc1 Minc2 Mjav Mare Mflo Mhap Pant'

rule download:
    input :
        "tables/download_table.tsv"
    output :
        "data/{sp}/genome.fa.gz"
    shell :
        "scripts/download.sh {sp}"

```

However, `snakemake` returned an error message I do not really understand :

```
Target rules may not contain wildcards. Please specify concrete files or a rule without wildcards.`.

```

Is there a way to write a single rule for downloading all the genomes?

# Answer
> 9 votes
The problem is that you need a master rule that requires all of your desired outputs as inputs, in your case it would be :

```
rule all:
    input:
        expand("data/{sp}/genome.fa.gz", sp=species.split(' '))

```

You'll also need separate download link inputs for each species. You could make a separate `download_table.tsv` for each species, but it would probably be easier to make a config file with this information, and add a params keyword to your rule. Something like:

```
rule download:
   params:
       url=config['locations']['sp']
   output :
       "data/{sp}/genome.fa.gz"
   shell :
       "scripts/download.sh {params.url}"

```

# Answer
> 3 votes
Given the Makefile you provided and which someone then deleted from the question, you should probably just add this line (this should be the first rule):

`rule all: input: ["data/{}/genome.fa.gz".format(x) for x in species.split()]`

This rule just specifies a list of expected output files and corresponds to the following lines from the original `Makefile`:

`GENOMES=$(patsubst %, data/%/genome.fa.gz, $(SPECIES))
all : $(GENOMES)`

---
Tags: snakemake
---


# Question
Title: Extract LINE-1 methylation status from WGBS data
I've been trying to figure out the correct way of extracting methylation status of particular elements from WGBS data, and found several possible approaches, but would like to discuss which one could be the best:

1. using software dedicated to the extraction of these elements on WGS (like TraFiC or MELT and then quantify methylation (like Bismark)
2. using bedtools to scan bam file with a filtered version of UCSC RepeatMasker table in bed format) and then quantify methylation
3. I've seen the use of custom perl scripts on fastq files like in this paper, in which the Supplemental Methods Section states: *"A custom Perl script was written to assay methylation of CpGs within a previously published motif of LINE-1 elements (Yang et al., 2004). Methylation levels for each of the four CpGs within the consensus LINE-1 element were determined globally for each WGBS sample prior to alignment using Line1\_FASTQ.pl. Dup15q and control were compared and tested for significance using an unpaired t-test."*

Any thoughts on this?

# Answer
> 0 votes
Have you checked out WGBS Tools on github? All of the custom scripts used in the paper that you referenced are in that repository. Each tool has a brief description of how to use it.

# Answer
> 1 votes
After a more extensive search, I found a useful answer that I'll share in case anyone else wonders how to do this. I'm using MethylKit package for differential methylation expression and it is quite easy to integrate with Genomic Ranges and rtracklayer in R to extract particular element methylation differences:

```
element.of.interest <- import.bed("FILENAME.bed") #download custom track from Genome Browser
element.meth <- diffmeth[as(diffmeth, "GRanges") %over% element.of.interest]

```

where "diffmeth" is the resulting MethylKit object from the differential methylation analysis.

# Answer
> 1 votes
Via BEDOPS:

```
$ bedmap —-echo —-echo-map-id-uniq <(bam2bed < reads.bam) <(rmsk2bed < repeats.rmsk) > answer.bed

```

---
Tags: methylation, transposable-elements
---


# Question
Title: How can I specify to DEseq2 to only perform comparisons on which I am interested with?
I am currently performing a large RNA-seq analysis from mice PBMCs. The dataset contains around 6,000 transcriptomic profiles and I would like to use DESeq2 to identify the sets of differentially expressed genes in the different conditions. In total, I have 100 biological stimulations, and for each stimulation I have 30 control samples and 30 samples treated with a molecule of interests (so is have: (30 controls+30 treated) * 100 stimulations = 6,000 samples).

I want to identify, for each stimulation, the sets of differentially expressed genes between control samples and treated samples. I do not want to compare samples from the different stimulations. In total, I would like to have thus 100 lists of differentially expressed genes.

I have started to use deseq2 to identify these lists but deseq2 is spending lot of time to perform comparisons on which I am not interested with (comparisons between the biological stimulations).

For now, I have a table sampleTable which looks like that:

I am using DEseq2 using the following command:

```
DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory = directory, design= ~ condition)

```

Could you please help me with that? How can I specify to deseq2 to not perform the comparisons between the biological stimulations, but rather between control and treated sample within each stimulation ?

Thank you and best,

# Answer
You can specify the exact comparisons you want in the `results()` function. So:

```
dds = DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory = directory, design= ~ condition)
dds = DESeq(dds)
res = results(dds, contrast("condition", "treatment1", "control1"))

```

Note that `condition` should not be, "`stim001-control1`, `stim001-control2`, etc.", but instead, "`stim001control`, `stim001control`, `stim001control`, `stim001treatment`, etc.". Don't put minus signs in your levels and don't number your samples in them (they're not in different groups.

The last command would be repeated for each of the comparisons. Note that this will still be quite slow due to the number of samples and the number of groups. With so many samples, you might just make a separate `sampleTable` for each of the comparisons you want to make. Each would then only contain the samples relevant for that comparison, so 60 total in each.

Alternatively, if you want to fit all of the groups at once then it's likely that limma/voom will prove to have better performance. It uses different math that happens to be quicker with large models.

BTW, if the slowness is coming from `DESeqDataSetFromHTSeqCount()` then you can mere everything on the command line into a single matrix. A couple lines of python or even the `join` command can do that quickly enough.
> 3 votes

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


# Question
Title: How to proceed after contigs.fa file obtained from VELVET assembly
I want to assemble a linear genome and I have a reference genome too.

I ran the following commands for de novo velvet assembly:

```
velveth auto 31,45,2 -fastq -shortPaired1 -separate 2-CA-SVA-2_S2_L001_R1_001.fastq 2-CA-SVA-2_S2_L001_R2_001.fastq

```

I chose auto\_33 for further analysis

```
velvetg auto_33 -exp_cov auto

```

Now my output files inside auto\_33 folder are:

```
contigs.fa
Graph2
LastGraph
Log
PreGraph
Roadmaps
Sequences
stats.txt

```

But I have no clue how to move ahead. Can I do anything else with velvet or do I have to switch to another program?

# Answer
> 3 votes
The contigs.fa file is a FASTA formatted file of the final contigs generated by velvet. I would start assessing the quality of the assembly that was generated. Here are some questions to get you started. Really the next steps depend on how you plan to use the assembly.

What is the distribution of contig lengths? Average and N50 contig length

How does the assembly compare to the reference?

* How much of the reference does the assembly cover?

When you align the original reads back to the contigs how well are they covered?

* Do any of the read alignments not make sense?
* Are there sections of the contigs that are not covered?

---
Tags: ngs
---


# Question
Title: How to confirm exon shuffling in a gene?
*Note: this question has also been asked on reddit*

I'm trying to confirm that the sequence of a novel gene is derived by exon shuffling between several different genes. I have the promoter sequence, gene sequence, and mRNA (with defined exon/intro boundaries). I've tried performing several database searches with the sequences, but each search produces a set of different unrelated hits. How would one go about using this info to confirm the exon shuffling hypothesis? Any tips are appreciated. Thanks!

# Answer
> 1 votes
These suggestions are based on your comment that you have the gene sequence.

For a statistical approach, you could formulate a null model and show that the sequence similarity you observe is much higher than expected under the null model, but reasonable under your exon shuffling model. Ideally, you would show this individually for each exon, to convince the reader about the exons' disparate origins. I'll propose a reasonable null model, though suggestions are gladly accepted here. Pick a random gene, assume it's built using exon shuffling, and search the rest of the transcriptome for sequences that best substantiate that argument.

For an evolutionary approach, you could search for orthologs (homologues? I never remember the difference) of your gene in other species. Hopefully, the results would show that your gene formed more recently than the exons it's built out of. If your gene actually has some other origin, maybe you'll find out it's older than some of its putative components. This can't exactly "confirm" your hypothesis, but it will test your hypothesis and just be generally a cool thing to try.

---
Tags: ngs, sequence-analysis, molecular-genetics
---


# Question
Title: Download data from the Human Microbiome Project via ascp
I have asked this question in biostars, but I am trying here as well as people working with "omics" data might be able to help. I think my issue relates understanding how large data storage on online server works. I am trying to download large data from the Human Microbiome Project. These are available on online server and can be downloaded via aspera (ascp).

I can't figure out the right path to get the data. On their website, they give an example of command to retreive an example dataset:

```
 ascp -l 100M -QT user@aspera.hmpdacc.org:/ReferenceGenomes/test_genome.tar /path/to/local/directory

```

I executed this, using my own username, from this I am asked my password, and I can authenticate. But then I get the following error:

```
(Error: Server aborted session: No such file or directory)

```

Then, I downloaded the manifest file via the HMP portal. In this manifest, they provide the path to which these files are supposed to be available. I tried to execute again the command with one of those path, provided by the manifest file:

```
ascp -l 100M -QT user@aspera.hmpdacc.org:/files/5a950f27980b5d93e4c16da124000378  /path/to/local/directory

```

I again get the same error of "No such file or directory".

The full url provided by the manifest

I tried putting different chunks of this address into to ascp command but nothing works. I think there is something I don't understand about how to access the server via aspera and what is expected to go in this path. When I copy paste this url in google, this directly download the file I want. So I'm close, but I need to do it via command line as I need to download hundreds of these files and they are huge (6-10GB each).

How can I download the data of the human microbiome project using ascp?

# Answer
> -1 votes
I can't test for ascp because they do not longer support my current setup but you can follow the FTP instructions. With FileZilla:

url: public-ftp.ihmpdcc.org

---
Tags: data-download
---


# Question
Title: Mark BWA-SW split alignments in output for long reads
Is there a way to remove split alignments of single-end long reads from BWA-SW output? BWA-MEM has an option to include or flag split alignments in the result, but it seems BWA-SW method include them in the result unless the split occur in paired-end reads.

# Answer
> 1 votes
From BWA manual, it is advised to use BWA-MEM in your longs read case (unless you have a good reason not to):

> BWA-MEM and BWA-SW share similar features such as long-read support and split alignment, but BWA-MEM, which is the latest, is generally recommended for high-quality queries as it is faster and more accurate.

I guess that BWA-SW returns several lines per transcript (based on local alignments). In that case, the issue has already been reported. The user asks:

> How can I differentiate between a read that has multiple hits in the reference and a read that has a split alignment?

The author replies:

> In principle, it is hard to tell. Bwa-mem makes an arbitrary choice. Mappings flagged with 0x100 are multiple hits. Mappings flagged with 0x800 are from a chimeric alignment.

So you could actually filter by flag. But I am afraid the output could include something you do not necessarily want to (since there is no specific flag for the split reads).

Maybe a better solution: you could switch to some mapper that will report the proper split in the hit of the read map, such as minimap2 from the same author. It is suitable for long reads and comes with a `splice -x` feature.

Then, the laziest solution I see here would be to remove the split alignments exploiting the sam format and in particular the CIGAR string. Meaning, to remove the sam fields containing a `N` in their CIGAR (aka, skipped region from the reference). Something like this would work:

```
cat alignment.sam | awk '$6 !~ /N/ {print $0}'

```

Hope this can be somewhat helpful.

---
Tags: sequence-alignment, bwa, sam, long-reads
---


title: 5'UTR and 3'UTR annotation in yeast

text:
# Question
Title: 5'UTR and 3'UTR annotation in yeast
I am working on a project in which I need to compute several parameters (such GC content and length) of 5'UTR and 3'UTR sequences of *Saccharomyces cerevisiae* yeast genes.

The problem is finding a proper annotation for these regions in yeast. I have tried with BiomaRt and these sequences are not available for yeast. I have also tried with UCSC table browser to obtain a BED file with 5'UTR or 3'UTR sequences but it is returning only coordinates of tRNAs, pseudogenes, ncRNAs.

In summary, **is there a straight-forward way to obtain an annotation for 5'UTR and 3'UTR regions in *Saccharomyces cerevisiae* yeast**?

# Answer
I am unaware of any "official" or gold-standard UTR annotations in *S. cerevisiae*.

One option is to use the annotations from the TIF-Seq publication (**Pelechano et al. 2013**).

The `GSE39128_tsedall.txt.gz` file contains the major isoforms identified. It would be up to you to computational associate each transcript with a given gene. It is also up to you to determine which major isoforms are present in your dataset. Furthermore based on my own experience there is often a significant number of discrepancies between what I see in an NGS dataset (e.g. RNA-Seq) and what is called as a possible isoform in these annotations. However, a lot of this could be related to the inherent technological limitations of the method.

If you look at the SGD genome browser they provide some additional UTR annotations. Derived from **Nagalakshmi et al. 2008**.

You can find the relevant annotations for those UTR's in their supplemental data section. The file is `1158441_tables_s2_to_s6.zip` and inside you will find *TableS4* which seems to contain what you seek.

There might be more resources but these are the two main ones I'm aware of.
> 9 votes

# Answer
I found the following two files in https://downloads.yeastgenome.org/sequence/S288C\_reference/:

According to the README files in the same directory, these are (the README for the 5' file is equivalent):

> SGD\_all\_ORFs\_3prime\_UTRs.README
>
> Information about the SGD\_all\_ORFs\_3prime\_UTRs.fsa file.
>
> This file contains sequences for predicted 3' UTR sequences for all characterized, uncharacterized and dubious ORFs in S288c.
>
> 3' UTR coordinates were retrieved from YeastMine (yeastmine.yeastgenome.org) on July 3, 2017 using the stop coordinates of both the predicted transcript and gene sequences. A BED file with these coordinates was made. BED files were fixed on Sept 22, 2017 to account for strandedness. Sequences were batch downloaded again in FASTA format from USCS Genome Table Browser on Sept 22, 2017 (https://genome.ucsc.edu/cgi-bin/hgTables?hgsid=604915183\_aUc9FZ1vwyAccGti1uI0wteYEM7x)
>
> File created: 2017-08-25

Those therefore seem to be the 3' and 5' UTR sequences of all yeast genes in FASTA format and should be exactly what you need.
> 7 votes

---
Tags: sequence-annotation, bed, gtf
---


title: Somatic tumor only variant calling?

text:
# Question
Title: Somatic tumor only variant calling?
I'm evaluating possibilites for somatic tumor variant calling without paired-normal samples. I'm aware of the consequences without a normal sample.

All the popular variant calling tools such as Strelka, VarScan etc require a normal sample.

What are the some of the possibilites I have? I'm looking for good sensitivity in exchange for higher false positives.

# Answer
There is a recent paper that attempts to do this:

> ISOWN: accurate somatic mutation identification in the absence of normal tissue controls.
>
> In this work, we describe the development, implementation, and validation of ISOWN, an accurate algorithm for predicting somatic mutations in cancer tissues in the absence of matching normal tissues. ISOWN is available as Open Source under Apache License 2.0 from https://github.com/ikalatskaya/ISOWN
> 5 votes

# Answer
I've been using the LoFreq\* caller for exactly this. It is designed to find variants with very low frequency, so is well suited for this type of analysis.

> LoFreq\* (i.e. LoFreq version 2) is a fast and sensitive variant-caller for inferring SNVs and indels from next-generation sequencing data. It makes full use of base-call qualities and other sources of errors inherent in sequencing (e.g. mapping or base/indel alignment uncertainty), which are usually ignored by other methods or only used for filtering.
>
> LoFreq\* is very sensitive; most notably, it is able to predict variants below the average base-call quality (i.e. sequencing error rate). Each variant call is assigned a p-value which allows for rigorous false positive control. Even though it uses no approximations or heuristics, it is very efficient due to several runtime optimizations and also provides a (pseudo-)parallel implementation.

(In my experience it is not fast at all, but that's what they claim).

I have used it many times to call somatic variants using the germline caller (and not its somatic one which requires a tumor/normal pair). Because the tool is so sensitive, its germline caller actually deals quite well with somatic samples.

You can find instructions for its use here.
> 5 votes

# Answer
> All the popular variant calling tools such as Strelka, VarScan etc require a normal sample.

Strelka and VarScan require a normal sample in somatic mode, but they both have a germline mode for unpaired analysis.

One of the most popular somatic variant callers is MuTect and that does not require a normal.
> 2 votes

---
Tags: variant-calling, cancer, genomics
---


title: BiomaRt error: Error in martCheck

text:
# Question
Title: BiomaRt error: Error in martCheck
I am experiencing an error while using getSequence() in biomaRt.

The code is given below:

```
library("biomaRt")
ensemblp <- useMart(biomart="plants_mart",host="plants.ensembl.org")
osj = useDataset("osativa_eg_gene", mart = ensemblp)
genes <- read.csv("osj.csv")
osj_id <- genes$ensembl_transcript_id

os_genes <- getSequence(id = osj_id, type = "ensembl_transcript_id", seqType = "peptide", mart = ensemblp)

```

The error message I am getting is :

```
Error in martCheck(mart, c("ensembl", "ENSEMBL_MART_ENSEMBL")) :
This function only works when used with the ensembl BioMart.This
function only works when used with the ENSEMBL_MART_ENSEMBL BioMart.

```

Did I make any mistake in code? Please help me fix the error. I tried asking this in Bioconductor support, but no one responded yet.

# Answer
getSequence has only been enabled for the main Ensembl (vertebrates) biomaRt. This is all at the Bioconductor end. You'll need to use a getBM instead, for example:

```
os_genes <- getBM(attributes =c("ensembl_transcript_id", "peptide"), filters = "ensembl_transcript_id", values = osj_id, mart = osj)

```
> 9 votes

# Answer
Biomart seems to be down most of the time so i think you can parse the ensemble ID and names from the annotation file which i did as i was having similar problem .

You can use this to make a file with ensemble ID and gene names for you use

```
awk '{
    for (i = 1; i <= NF; i++) {
        if ($i ~ /gene_id|gene_name/) {
            printf "%s ", $(i+1)
        }
    }
    print ""
}' Homo_sapiens.GRCh37.70.gtf | sed -e 's/"//g' -e 's/;//g' -e 's/ /\t/' | sort -k1,1 | uniq > Homo_sapiens.GRCh37.70.txt

```
> 1 votes

---
Tags: r, biomart
---


# Question
Title: bigWig to bed for regions above/below threshold
I'd like to threshold a bigWig file on a value and get all regions above or below that value in bed format. Ideally without too many file conversions. Are there any tools to do that?

# Answer
You can do that with a bit of python:

```
#!/usr/bin/env python
import pyBigWig
threshold = 10  # Change me
bw = pyBigWig.open("file.bw")  # Change me
of = open("regions.bed", "w")  # Change me

for chrom, len in bw.chroms().items():
    intervals = bw.intervals(chrom)
    for interval in interals:
        if abs(interval[2]) > threshold:
            of.write("{}\t{}\t{}\n".format(chrom, interval[0], interval[1]))
bw.close()
of.close()

```

You'll need to install pyBigWig (it's available via `pip` or `conda install`) and change the lines with `Change me` in them.
> 1 votes

# Answer
Convert bigWig to WIG via `bigWigToWig`:

```
$ bigWigToWig signal.bw signal.wig

```

Convert WIG to BED via BEDOPS `wig2bed`:

```
$ wig2bed < signal.wig > signal.bed

```

Filter BED by its score column with `awk`:

```
$ awk '($5 < 1)' signal.bed > signal.filtered.bed

```

The test `$5 < 1` filters elements from `signal.bed` that have a score value less than `1`. Change this to the threshold test you would want to apply.
> 1 votes

---
Tags: bed, bigwig
---


title: What does "motif 2 bp or longer" mean?

text:
# Question
Title: What does "motif 2 bp or longer" mean?
I'm reading a preprint on a new somatic variant discovery tool, and the first sentence of the Discussion section has me confused.

> *Across the four datasets analyzed in this study, we discovered that the major source of disagreement between callers originates from somatic variants called within STRs, in particular if the motif is two base pairs or longer.*

The "*two base pairs or longer*" bit is what confuses me. I interpreted "motif" as referring to the short tandem repeats (STRs), which are by definition 2 bp or longer. A motif of \< 2 bp is...a homopolymer, not really a motif in the common sense.

Am I missing something here? Or is this simply an awkward word choice on the part of the authors?

# Answer
What's your reference for the definition of STRs? I think it is still ambiguous among the community. Wikipedia states motifs of \>=2 bp. However other references include homopolymers as well:
Source 1, source 2 and source 3.
> 4 votes

---
Tags: variant-calling, genomics, repeat-elements
---


# Question
Title: Interpreting large z scores from microarray data
I have microarray z scores for gene set enriched (GSE) pathways that range approximately from -3 to 14. I understand that an absolute value of 1.96 has a p-value = .05 and significant for our purposes. My question is whether I can consider extreme z scores as more significant?

# Answer
This is a statistical question. What a higher z value means is that it is more extreme (if you assume a normal distribution), thus a more extreme value is less likely to happen by chance. Which is translated to a (lower) p-value.

What it is significant and not is arbitrary. Prior to the experiment you must decide which is your threshold, how much chance you allow. Apparently you decided to select 0.05 as your threshold, which means that one of 20 results might be due to chance. But your extreme z score is not **more** significant, it is **less expected** (under the null hypothesis, those p-values might change if you change your null hypothesis, ie you think your genes expression don't follow a normal distribution) but **as significant as** your 0.0499999 p-values.
> 2 votes

---
Tags: microarray
---


# Question
Title: combining the use of workdir and option --jobscript in snakemake
It seems that in snakemake the script specified after the `--jobscript` cannot be used properly if a `workdir:` is specified in the snakefile. The path of the script specified becomes relative to the workdir defined in the snakefile instead of being relative to the current working directory. I do not know if this is a feature or a bug !!

Nevertheless, it becomes quite painful to workaround this problem. The only way I found was to parse the command line in the snakefile, for example: `snakemake -c "qsub" -j 30 --js ./sge.sh --latency-wait 30 -rp` in order to copy `./sge.sh` in the output directory defined by the keyword `workdir` in the snakefile.

As a side effect, if you specify option for qsub on the command line, for example: `snakemake -c "qsub -e ./logs/ -o ./logs/" -j 30 --js ./sge.sh --latency-wait 30 -rp` the `logs` folder must be created under the workdir directory.

Is there something I don't understand with the `--jobscript` option? Am I not using it the right way? Am I going against the best practice?

# Answer
> 2 votes
To quote from the snakemake documentation

> All paths in the snakefile are interpreted relative to the directory snakemake is executed in. This behaviour can be overridden by specifying a workdir in the snakefile:
>
> `workdir: "path/to/workdir"`
>
> Usually, it is preferred to only set the working directory via the command line, because above directive limits the portability of Snakemake workflows.

It's not explicitly stated, but it's sort of implied that all relative paths then become relative to the working directory. I would expect that specifying an absolute path would get around this.

As an aside, in my mind setting the working directory is usually only needed on clusters without a shared file system (presumably shared between the worker nodes, but not with the head node), since there you can't `cd $WORKDIR` before running snakemake. This is then normally done with your scheduler, in such cases.

# Answer
> 1 votes
It's just a guess, but I think you might want to use `--keep-target-files` flag to not to adjust paths of targets. From snakemake manual :

> `–keep-target-files`
>
> Do not adjust the paths of given target files relative to the working directory.
>
> Default: False

So if I understand right Devon's answer, combination of this flag and setting of `workdir` could force compute target in different directory that the snakemake is executed in.

---
Tags: snakemake
---


# Question
Title: What is local realignment and what is the problem it solves?
I am trying to understand local realignment but I could not get a clear idea of what is the problem solved by it.

For example, reading Homer and Nelson (2010):

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

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

# Answer
SNPs are likely to be created and InDels are likely to be missed. Suppose you have a read, `ACTGACTGACTGTAC` and you align it to a reference sequence `ACTGACTGACTGTTAAGAACGACTACGAC`. If you aligned that, you would either get:

```
ACTGACTGACTGTac    (lower case denotes soft-clipping)
ACTGACTGACTGTTAAGAACGACTACGAC

```

or

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

```

But what happens if you look at a bunch of reads? See below:

```
ACTGACTGACTGTAC
  TGACTGACTGTACTA
    ACTGACTGTACTAAG
ACTGACTGACTGT--TAAGAACGACTACGAC

```

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

Raw alignments in IGV:

Alignments after realignment:

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

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


title: Highlight CpG sites with Texshade

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

# Answer
The best solution I could find was to use the `\shaderegion` function on all the sequences so it ended up looking something like this:

```
\shaderegion{1}{CG}{Red}{OrangeRed}

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

...

```

Then I cast the rest of shading similarity to grey so the CG pairs were more easily distinguishable.

```
\shadingcolors{grays}

```

Two examples one with `\fingerprints` and one regular.
> 4 votes

---
Tags: sequence-alignment, methylation, visualization
---


# Question
Title: Verify a predicted protein in one genome in a different genome of the same species
I have two genome assemblies of the same non-model species, call them Assembly 1 (generated from Illumina data) and Assembly 2 (generated from PacBio data).

For Assembly 1, I also have predicted proteome data, generated with EVM. Say there is a protein, call it Protein X, for which I have the fasta sequence from predictions, but suspect an error in Assembly 1 (because the Protein X has a long segment of auto-identity). For this reason, I would like to verify the assembly of the locus coding for Protein X in Assembly 2 and/or verify the presence of Protein X in the Assembly 2.

Data in my possession:

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

Things I was thinking of:

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

Any ideas on how I could do that?

# Answer
> 3 votes
This should be very easy to do. Here are some options:

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

Of the two, using the first approach is much better since it will try to build a valid gene model for you instead of simply finding the regions of high homology as blast will do.

Finally, just for completeness, you can go old school and run a full de-novo gene prediction tool on your assembly. Something like GeneID for instance. But that's certainly overkill if you're just looking for a single gene.

---
Tags: proteins, assembly, validation
---


title: Matlab scatter diagram

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

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

# Answer
> 4 votes
A reverse image search on google leads to the exact code used to produce the second image:

```
[Xcoarse, Ycoarse] = meshgrid([0 1 2 3], [0 1 2 3]);
[Xfine, Yfine] = meshgrid(linspace(0,3,3000), linspace(0,3,3000));
DataCoarse = [ 1 2 4 1; ...
           6 3 5 2; ...
           4 2 1 5; ...
           5 4 2 3];
DataBicubicFine  = interp2(Xcoarse, Ycoarse, DataCoarse, Xfine, Yfine, 'bicubic');
figure
surf(Xfine, Yfine, DataBicubicFine); shading flat; colormap(jet); view(0, 90)
 hold on; plot3(Xcoarse, Ycoarse, 10*ones(size(Xcoarse)), 'k.', 'MarkerSize', 20)
colorbar
print -dpng -r300 BicubicInterpolationExample.png

```

---
Tags: software-recommendation, statistics
---


# Question
Title: Filtering out all seqs with mutations of list2 from list1
I have 2 lists, `list1` and `list2`, of protein sequences of the same given gene in different strains. In `list1` are reported amino-acid sequences from resistant strains, and in `list2` are reported amino-acid sequences from susceptible strains. Each item of each list is an amino-acid sequence.

I want to filter out (in Python 3) from `list1` all of the variants that are also present in `list2`. The final aim is to have a list of the sequences that have unique mutations that are not in the 2nd list. The mutations can be of any type (e.g., indels, mismatches).

All the sequences are of at least 90% identity to some reference sequence (filtered out the low matches using seq alignment).

Any enlightenment / package / tool for this job will be appreciated!

Toy example:

`a`, `b`, `c` : amino-acids

For this input:

1. Reference seq: `aaabbb`
2. Resistant strains: `list1` = \[aa**c**bb**c**, a**b**abbb\]
3. Susceptible strains: `list2` = \[aa**c**bbb, aaabb**c**\]

The expected output would be: \[a**b**abbb\], since the mutations in `list2` cancels out the item aa**c**bb**c**.

Of course this is a simple example and it is possible to have more complex scenarios (e.g., longer indels and mismatches).

# Answer
> 1 votes
Samtools/VCF is overkill for comparing amino acid sequences from a single protein. I would suggest making pairwise alignments between each sequence in list2 and the reference and making a list of mutations (2c and 6c in your example). Then repeat the procedure for list1, discarding sequences if they contain mutations in your stored list. This can all be done using biopython. Here's how to make the alignments to get you started:

```
from Bio import pairwise2
from Bio import SeqIO
seq1 = SeqIO.read("alpha.faa", "fasta")
seq2 = SeqIO.read("beta.faa", "fasta")
alignments = pairwise2.align.globalxx(seq1.seq, seq2.seq)

```

---
Tags: python, proteins, sequence-alignment, filtering
---


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

---

I have a strange situation. I have two BLAST databases (`20170821` and `20170824`) that both contain a certain subject sequence `LC074724`. When I query them both with a given query, `Q`, in one case (`20170824`) there is a **very strong match** and in the other (`20170821`) **no match at all**.

Why is there not even a lower quality match found although the databases are nearly the same?

```
BLAST version

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

```

Both databases contain an identical subject sequence `LC074724`:

```
$ blastdbcmd -entry LC074724 -db 20170824 | md5
fcfbc167c7f81dfd75aad1bf3db6220b

$ blastdbcmd -entry LC074724 -db 20170821 | md5
fcfbc167c7f81dfd75aad1bf3db6220b

```

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

# Answer
I'm trying to figure out why. This will be a longer read, **Tl;dr at the end**:

## Doing a match

There is a good match of `Q` against `LC074724` in `20170824`:

```
$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170824 -query Q.fasta \
  -task blastn -max_hsps 1 | grep LC074724
Q   LC074724.1  2581

```

Note that the bitscore of the `LC074724` match is `2581`.

But there is **no match** of `LC074724` **at all in** `20170821`:

```
$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170821 -query Q.fasta \
-task blastn -max_hsps 1 | grep LC074724

```

The option `-outfmt '6 qaccver saccver bitscore'` tells BLAST to print the query accession number and version, the subject accession and version, and the bitscore of the match.

## What are the best hits for 20170821

So what are the bitscores of the Q matches in 20170821?

```
$ blastn -outfmt '6 bitscore' -db 20170821 -query Q.fasta -task blastn -max_hsps 1 | sort -nr | head
3301
2616
2590
2590
2563
2554
2547
2544
2540
2540

```

The bitscore of `Q` against `LC074724` in `20170824` (2581) would have placed it **in 5th position in the top 10 matches of Q**, but `LC074724` is **not matched at all** in `20170821`! The lowest bitscores matched in `20170821` is 2354.

## What's the difference between the databases?

Get the ids of the subjects in each database

```
$ blastdbcmd -entry all -db 20170824 | egrep '^>' | cut -c2- | sort > 20170824.ids
$ blastdbcmd -entry all -db 20170821 | egrep '^>' | cut -c2- | sort > 20170821.ids

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

```

`20170821` has 643 sequences that are not in `20170824`:

```
$ comm -23 20170821.ids 20170824.ids | wc -l
     643

```

`20170824` has no sequences that are not in `20170821`:

```
$ comm -13 20170821.ids 20170824.ids | wc -l
       0

```

and the two have 9111 sequences in common:

```
$ comm -12 20170821.ids 20170824.ids | wc -l
    9111

```

These numbers match with the numbers from `wc -l` above.

So `20170821` is just `20170824` with an extra 643 sequences.

## Maybe LC074724 matches, but with a low bit score?

BLAST only shows the top 500 matched subjects by default. Let's ask it to show more:

```
$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170821 -query Q.fasta \
  -task blastn -max_hsps 1 -max_target_seqs 10000 > 20170821-10000-matches.txt

```

And....

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

```

There it is! And **it's in 5th position, just as it should be** (based on the bitscores discussed above).

Just to confirm, taking away the `-max_target_seqs` option gets no match:

```
$ blastn -outfmt '6 qaccver saccver bitscore' -db 20170821 -query Q.fasta \
  -task blastn -max_hsps 1 | grep LC074724
$

```

## Three confusingly named BLAST options

I've always found the names and descriptions of the following three BLAST command-line options very confusing.

From the Formatting options section of blastn -help:

> -num\_descriptions =0\> Number of database sequences to show one-line descriptions for Not applicable for outfmt \> 4
> Default = \`500' * Incompatible with: max\_target\_seqs
>
> -num\_alignments =0\> Number of database sequences to show alignments for Default = \`250' * Incompatible with: max\_target\_seqs

And from the Restrict search or results section:

> -max\_target\_seqs =1\> Maximum number of aligned sequences to keep Not applicable for outfmt \<= 4 Default = \`500' * Incompatible with: num\_descriptions, num\_alignments

From the section names, it seems like only the latter option (`-max_target_seqs`) would have any effect on the search, and that the former two are **just about what is displayed**.

But... using either just `-num_descriptions` or `-num_alignments` with a big value also gets a match:

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

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

```

So BLAST is apparently finding the `LC074724` match even without the high `-max_target_seqs` option, but it's not displaying it unless there is a high `-num_alignments` or `-num_descriptions` value. (**Keep in mind, this is the 5th best hit! Not somewhere \>500**)

Just to confirm, when none of the 3 options are given and no -outfmt option is either, the match is not found:

```
$ blastn -db 20170821 -query Q.fasta -task blastn -max_hsps 1 | grep LC074724
$

```

---

**Tl;dr:**

In summary, the options

```
 -evalue
 -max_target_seqs
 -num_alignments
 -num_descriptions

```

are **NOT** just filters applied to the results, but instead directly influence the results **found**.

This means you can End up with a **bacterial** sequence as **BEST hit** when using a value `x`, and with no bacterial sequences **at all** in your top20 when using value `y`.

To quote the NCBI user service:

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

And

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

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

---
Tags: blast, sequence-alignment, filtering
---


# Question
Title: Quantification of EST/cDNA based retained intron transcripts
*Note: this question has also been asked on Biostars*

We are seeing a retained intron transcript event for some RNAseq samples, and we want to assess at the sequence level which intron retention events are actually occurring so that we can do protein domain predictions. For example, we have transcript level expression support for DROSHA-203, which according to Ensembl is a retained intron transcript that potential gives rise to multiple splicing variants via the cDNA sequence information:

In theory, the orange sites indicate where a mutation might occur which could obliterate a splicing recognition site, thereby causing an retained intron event. So in total there are 6 potential mutational events which might give rise to 4 possible retention events. Given that 4 intronic regions might be retained, the number of unique retention events is a combinatorics question equal to 2^4, e.g., intron1 retained only, vs intron1 + 2 retained, all possible unique combination, etc, correct?

Is there a better way to more granularly quantify these possible transcripts? Or do we have to PCR validate all of them to figure out which one's are actually occurring? Is there a way that I can look at the reads mapping to the corresponding intronic regions to make even a probabilistic assessment of which transcripts are present and to what extent they are expressed? The libraries are prepped for total RNA with a ribozero depletion, so I should still be able to quantify any enrichment to these intronic regions, correct?

# Answer
> 1 votes
You can do full-length whole-transcriptome sequencing using either the direct RNA kit or a cDNA kit (with or without PCR) on a MinION. This will allow you to quantify isoforms directly, rather than estimating by looking at read pairs.

However, it'd be a bit expensive to do it that way if you're looking at a single transcript (and aren't barcoding to multiplex samples on the same run).

If you've already done Illumina whole-transcriptome sequencing (e.g. strand-specific paired-end sequencing), tools like Kallisto, RSEM, or Salmon can be used to estimate isoform expression based on the mapped reads. These tools use the information about transcript placement and/or similarity to alter the observed read counts in an attempt to better represent reality.

Unfortunately, if you want these tools to quantify intronic sequence (i.e. the $2^4$ possibilities you have suggested), then those sequences need to be included in the sequences provided. They won't estimate frequencies for non-existent sequence.

As an alternative to including every possible gene model, something like a SuperTranscript model might work. I'm not familiar with Alicia Oshlack's method, but it at least gives the appearance that it might reduce the complexity of searching for multiple unannotated isoforms within a gene.

---
Tags: rna-seq
---