# Question
Title: How to query a network for a bacterium, specifically Streptomyces caatingaensis or Streptomyces thioluteus?
I have a list of gene/protein IDs and I could not query a network using String, BioGrid, Intact and IIS, for Streptomyces thioluteus.

The bacteria I am working with is the *Streptomyces caatingaensis*. Since it is newer than others, I thought I would find an interaction network looking for *Streptomyces thioluteus* (which shares one pathway I am analyzing in common, others do not share). As a result, I have this list of proteins, but I could not query any network from the IDs. Sometimes I could not find the bacterium in the tools and when it was found, as in IIS the result was **zero** interactions.

Is there anything I could use to have some network that could be used as insight for the analysis?

# Answer
> 1 votes
I realized that some tools work better with Gene Names instead of some type of ID.

Using IIS I was capable of creating a network using this list of proteins (*Streptomyces caatingaensis*) by querying proteins by Genes Names having as base the E. Coli.

The tool has 3 main tabs: proteins, annotation, and interactome. Search the Gene Names in proteins section, go to Interactome, set the parameters and create the network.

P.S. If you are not used to the tool you should read the research article available on the home page before using the network for research purposes.

---
Tags: database, proteins, networks
---


# Question
Title: How can I colour boxes in Gviz AnnotationTrack in R?
I'm learning the Gviz bioconductor package, I generate a plot as follows:

```
library(Gviz)
track <- AnnotationTrack(start=c(1,5,7), end=c(2,6,10), strand=c('*','*','*'), stacking="dense", showFeatureId=TRUE, id=c('red','blue', 'red'))
pdf(file=paste("test","pdf", sep="."))
plotTracks(track)
dev.off()

```

`test.pdf` looks like this:

However, what I want is for the boxes to be coloured rather than labelled with the colour word, how to do this?

Note that my real example is much more complex than this.

# Answer
Did you try the `fill` argument?

Something like this:

```
track <- AnnotationTrack(start=c(1,5,7), end=c(2,6,10), strand=c('*','*','*'), stacking="dense", showFeatureId=TRUE, id=c('red','blue', 'red'), fill=c('red','blue', 'red'))

```
> 8 votes

---
Tags: r, bioconductor, gviz
---


# Question
Title: What happens if a major bug is discovered in a bioinformatic package that has been used in published literature?
Yesterday I was debugging some things in R trying to get a popular Flow Cytometry tool to work on our data. After a few hours of digging into the package I discovered that our data was hitting an edge case, and it seems like the algorithm wouldn't work correctly under certain circumstances.

This bug is not neccessarily the most significant issue, but I'm confident that it has occurred for other users before me, perhaps in a more insidious and less obvious way.

Given that since this tool was released, literature has been published that utiltised it, **how does the profession handle the discovery of bugs in these widely used packages?**

I'm sure that given the piecemeal nature of a lot of the libraries out there, this is going to happen in a much more significant way at some point (i.e. invalidating a large number of published results)

For some context I'm a programmer who's been working across various other parts of software development, and is quite new to Bioinformatics.

# Answer
I prefer to treat software tools and computers in a similar fashion to laboratory equipment, and in some sense biology in general. Biologists are used to unexpected things happening in their experiments, and it's not uncommon for a new discovery to change the way that people look at something. Things break down, cells die off quicker on a Wednesday afternoon, results are inconsistent, and that third reviewer keeps on about doing that thing that's been done a hundred times before without anything surprising happening (just not *this* time). It's a good idea to record as much as can be thought of that might influence an experiment, and for software that includes any input data or command line options, and especially software version numbers.

In this sense, a discovered software bug can be treated as a new discovery of how the world works. If the discovery is made public, and other people consider that it's important enough, then some people might revisit old research to see if it changes things.

Of course, the nice thing about software is that bugs can be reported back to the creators of programs, and possibly fixed, resulting in an improved version of the software at a later date. If the bug itself doesn't spark interest and the program gets fixed anyway, people unknowingly use newer versions, and there might be a bit more confusion and discussion about why results don't match similar studies carried out before the software change.

If you want a bit of an idea of the biological equivalent of a major software bug, have a look at the barcode index switching issue, or the cell line contamination issue.
> 21 votes

# Answer
While I agree with gringer's answer I would first report the bug to the maintainers with a test case (if possible) in order to have the bug corrected.

The bug might be an erroneous implementation in a corner case or in a more common use case of the software. But it is a good opportunity to collaborate to increase the quality of the software, either by introducing tests (or expanding them), by re-evaluating the state of the art with other tools, or to communicate with users of other software that depends on the tool.

Depending on the bug and the time one is able to spend, do a collaboration to study the effect of such a bug in the published studies using the software.

Sometimes the bug doesn't affect much, but sometimes it can have quite an impact which can be difficult to estimate. For instance, I found that the reactome.db package didn't have the same annotations that the available online database (in version 1.58.0), due to an incorrect mapping. However, it was impossible for me to evaluate the damage of such an error on the published literature, because gene set enrichment tools (usually) don't report the data they are based on and the versions of the databases used are often under reported in the literature.
> 9 votes

# Answer
If you have a software development background you probably know about filing bug reports. Your question seems to be more about how to correct the scientific literature which has been affected by the bug.

The traditional way that this has been done is by writing an paper that explains the bug and gives examples of how it may impact published results (see gringer's examples of index switching and HELA contamination). This way the community will be aware that older results using the implicated method need to be verified before building upon them. Unfortunately, awareness of these issues doesn't always diffuse through the community as quickly as it should, but it's probably still the best tool we have to ensure that science is self-correcting.
> 7 votes

# Answer
In extreme cases, it can destroy careers:

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

Most cases are not that extreme. In fact, realistically, I would guess that *most* papers have some data processing errors of some kinds, yet still we muddle on towards truth. And this isn't unique to software – the whole idea of independent replication in science is that we're always probably messing things up in ways we don't realize, so we need to re-do things and hope that the errors are different this time. But software is unusual in how tiny errors can cause very large effects, and often the effects are hard to detect. If you drop a test tube on the floor, you'll probably notice. If a piece of new software gives incorrect results... well, we didn't know what the results were supposed to look like, that's what makes it science! There are things you can do to mitigate this somewhat, but to some extent it's just inherent to the scientific process.

I still like to show that blog post to new grad students to scare them into writing tests. I actually think it's a good habit to, whenever you start using a new software package, start by writing some tests for it – it's a good way to check your assumptions and understanding about what it actually does, and even if it's a popular package by eminent researchers, you'll still sometimes be shocked by what you find.
> 7 votes

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

---
Tags: r, software-quality
---


# Question
Title: Filter BAM file for read pairs where one or both of the reads starts with a given sequence pattern
What's the most straightforward way to filter a BAM file for read pairs where one or both of the reads starts with a given sequence pattern? Would a tool that deals with UMIs work for this? Which one would be best?

# Answer
**First option**

My solution uses the following steps:

* use picard sortsam to sort the records on query-name (not samtools sort because the order is not the same between java and C )
* use jjs (java scripting engine) and the htsjdk library to build a `buffer`of reads having the same name. If any read starts with a pattern, print the whole `buffer`. The htsjdk should be already contained in the picard.jar.
* sort back to coordinate order using samtools.

All in one:

```
java -jar /path/to/picard.jar SortSam I=in.bam O=/dev/stdout  SO=queryname |\
jjs -cp /path/to/picard.jar  script.js |\
samtools sort -T TMP -o out.bam -

```

with script.js

```
var out=java.lang.System.out;
var File = Java.type("java.io.File");
var SamReaderFactory = Java.type("htsjdk.samtools.SamReaderFactory");
var SamInputResource = Java.type("htsjdk.samtools.SamInputResource");
var SAMFileWriterFactory = Java.type("htsjdk.samtools.SAMFileWriterFactory");
var samReader =  SamReaderFactory.makeDefault().open(SamInputResource.of(java.lang.System.in));
var samWriter =  new SAMFileWriterFactory().makeBAMWriter(samReader.getFileHeader(),true,java.lang.System.out);
var iter=samReader.iterator();
var buffer=[];
for(;;)
        {
        var rec=null;
        if(iter.hasNext()) rec=iter.next();
        if(rec==null || (buffer.length>0 && !rec.getReadName().equals(buffer[0].getReadName())))
            {
            var i=0;
            for(i=0;i< buffer.length;++i)
                {
                if(buffer[i].getReadString().match(/^AAAA/)) break;
                }
            if(i!=buffer.length)
                {
                for(i=0;i< buffer.length ;++i)
                    {
                    samWriter.addAlignment(buffer[i]);
                    }
                }
            if(rec==null) break;
            buffer=[];
            }
        buffer.push(rec);
        }

samReader.close();
samWriter.close();

```

**A second option** is to use samtools view and awk to export the read names:

```
 samtools view  input.bam | awk -F '\t' '{if($10 ~ /^AAAA/) print $1;}' | sort | uniq > names.txt

```

and to use Picard FilterSamReads and READ\_LIST\_FILE to filter the reads.

**UPDATE: 3rd solution:**

I've updated my tool samjdk: there is now a new option `--pair` to work with query-name-sorted BAM files.

The command would be:

```
java -jar picard.jar SortSam I=S1.bam O=/dev/stdout  SO=queryname |\
java -jar dist/samjdk.jar --samoutputformat BAM --pair \
     -e 'return records.stream().anyMatch(R->R.getReadString().startsWith("AAAAAAG"));' |\
samtools sort -T TMP -o out.bam -

```
> 3 votes

# Answer
For at least **one in a pair** matching the sequence, if the read names of pairs are identical, some grep magic could do the job:

1. Print the sequence at the beginning of the line and read the name in the second column
2. Find all lines starting with the sequence
3. Extract only read names

And then just grep all found reads and save into the BAM file:

```
SEQ=<my_sequence>
samtools view unfiltered.bam | \
   awk '{print $10 "\t" $1}' | \
   grep "^$SEQ" | \
   cut -f 2 > list_of_matching_names
samtools view unfiltered.bam | \
   grep -f list_of_matching_names | \
   samtools view -b > filtered.bam

```

Two small notes:

1. You might want to consider also reverse complementary sequences.
2. You also might want to keep header of BAM file (you can extract it at the beginning and add it in the end).
> 3 votes

# Answer
I am one of the authors of the `UMI-tools` pacakge, and although you probably could hack together a way to do this with our package, its not what its designed for and the solutions above are probably better.

The best way to do this would be to filter the fastq before mapping if thats a possibility, if it has to be the BAM, then here is a python implementation similar to Pierres java solution above. This one fiddles about with strands and cigar operations a bit more.

First sort the reads by name with `samtools sort -n`

then

```
import pysam
inbam = pysam.AlignmentFile("mybam.bam")
outbam = pysam.AlignmentFile("mybam.filtered.bam", "wb", template=inbam)
seq_to_find = "AAGG"
rc_seq_to_find = "CCTT"

def test_and_out_buffer(read_buffer, outbam, seq_to_find, rc_seq_to_find):
    '''This function takes a list of reads and check if any of them
    have the prespecific sequences at their start. If the sequence
    is detected the whole buffer is output. The start of the read
    is the end of the alignment if the read is on the '-' strand.
    Also checks that the read doesn't start with hard-clipped sequence'''

    do_output=False

    for segment in read_buffer:

        if segment.is_reverse:
            # the last operation is stored in the first part of the
            #last of the cigar tuples. Operation 5 is Hard-clip
            output_condition = (not segment.cigartuples[-1][0] == 5) and \
              segment.query_sequence.endswith(rc_seq_to_find)
        else:
            output_condition = (not segment.cigartuples[0][0] == 5) and \
              segment.query_sequence.startswith(seq_to_find)

        if output_condition:
            do_output=True
            break

    if do_output:
       for segment in read_buffer:
           outbam.write(segment)

last_name=""
read_buffer=[]

# Loop over reads in the file
for read in inbam.fetch(until_eof=True):

    # if the name hasn't changed, the name is part of the
    # same "pair"
    if read.query_name==last_name:
        buffer.append(read)
        continue

    # if we got here then the a new pair is being examined
    # test the old pair and output if neccesary.
    test_and_output_buffer(read_buffer, outbam, seq_to_find, rc_seq_to_find)

    # make a new buffer with just the single new read.
    last_name=read.query_name
    read_buffer=[read]

# the last set of reds will not have been output.
test_and_output_buffer(read_buffer, outbam, seq_to_find, rc_seq_to_find)

```
> 0 votes

---
Tags: sam, sequence-analysis, filtering
---


# Question
Title: How to colour multiple residues in Pymol?
We can use

```
color yellow, resi 3
color yellow, resi 3-8

```

to colour one residue or a range of residues. But what if I want to colour residue 3-8 and 10-20, how to write it?

Also, my protein has two chains, light chain (L) with 100 residues and heavy chain (H) with 100 residues. How to colour the residue 10-20, 30-40 of light chain, and the residue 5-10, 20-30 of heavy chain in one go?

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

```
select toBecolored, resi 3-8
color yellow, toBecolored

```

Regarding the L/H chain part, if they are effectively annotated as L and H, something like the following should work:

```
select toBeColored, ((i. 10-20 or i. 30-40) and c. L ) or ((i. 5-10+20-30) and c. H)
color yellow, toBeColored

```

Combining both commands is probably possible as:

```
color yellow, ((i. 10-20 or i. 30-40) and c. L ) or ((i. 5-10+20-30) and c. H)

```

Pymol allows synonyms (eg. resi and i. or chain and c.).

The following pages may be also helpful:

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

---
Tags: pymol
---


# Question
Title: Least present k-mers in the human genome
What are the least present *k*-mers in the human genome at different sizes?

Starting with k=4 and going up in size until k=10, what are the *k*-mers least seen (or not at all) in the human genome? I am only interested in the reference human genome, so I am discounting SNPs/Indels in the population.

If this is not pre-calculated somewhere, what tool is recommended to use that starts with the *GRCh38* reference as input?

# Answer
> 12 votes
You can use the Jellyfish software to calculate the k-mer profiles up to length 31.

From the instructions in the user guide:

The basic command to count all k-mers is as follows:

```
jellyfish count -m 21 -s 100M -t 10 -C reads.fasta

```

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

```
jellyfish histo mer_counts.jf

```

To query the counts of a particular k-mer, use the query subcommand (see section 3.3):

```
jellyfish query mer_counts.jf AACGTTG

```

To output all the counts for all the k-mers in the le, use the dump subcommand (see section 3.2):

```
jellyfish dump mer_counts.jf > mer_counts_dumps.fa

```

# Answer
> 10 votes
You can also use R. I give you an example of only chr1 and only kmer=4.

```
library(BSgenome.Hsapiens.UCSC.hg38)
library(Biostrings)

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

m <- as.matrix(kmers)

m[order(m),]

```

# Answer
> 4 votes
Grabbing chromosomes for `hg38`:

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

```

Via `kmer-counter` and Python, here's how to search for kmers of length 7 from chromosome `chrY`:

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

import sys
import subprocess
import itertools

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

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

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

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

```

This script will print a two-column tab-delimited text file to standard output, where the first column is the 7mer and the second column is the count of that 7mer and its reverse complement over `chrY` (including zero-counts):

```
CGACGCG 20
CGTCGCG 20
CGCGATA 23
TACGCGC 25
...
AAATAAA 33521
TTTCTTT 34014
GAATGGA 35361
AATGGAA 36906
TTTTTTT 103093

```

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

---
Tags: human-genome, sequence-analysis
---


title: IMPUTE2 missing filter option

text:
# Question
Title: IMPUTE2 missing filter option
I went through the manuals and can't find the meaning of the below flag:

```
-filt_rules_l 'filter==0'

```

Here is the link to manual for this flag: filt\_rules\_l

Can anyone explain?

IMPUTE version:

```
impute_v2.3.1_x86_64_dynamic

```

---

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

# Answer
> 3 votes
I have never used this program, but reading the manual in the link it has to mean the following.

The manual says:

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

So, `-filt_rules_l 'filter==0'` would mean that the program removes any variants with `filter` values of 0 from the reference panel. In the legend file there should be a column named `filter`, which are filtered out when the value is 0.

Here is the example legend file, in your case there must be a custom additional column named *"filter"*:

```
$ head ALL_1000G_phase1integrated_v3_chr1_impute.legend

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

```

---
Tags: impute2
---


# Question
Title: How do PCR duplicates arise and why is it important to remove them for NGS analysis?
I am trying to understand PCR duplicates in NGS analyses (actually whole-genome). I searched, and the best answer I found is in this blog.

However I don't understand if I understood how PCR duplicates arise correctly because I cannot see the problem of having them in the downstream analysis - aside from computational problems, i.e. unnecessary redundancy.

If I understood correctly, PCR duplicates arise during library preparation when PCR amplifies the fragments with adapters. In this case, if you have duplicate fragments, you'll end up amplifying some of them twice or more.

However, when you are doing the mapping, they should fall on the same region, most likely decreasing the quality of the mapping (since they increase the consensus on a specific sequence, which could have been subject to sequencing errors). But aside from that, you should have the same mapping with respect to the mapping where duplicates are removed, although with decreased quality.

Is the quality problem the actual reason for removing PCR duplicates or is there anything that I'm missing?

# Answer
> 13 votes
In any scenario where depth of coverage is an important factor, PCR duplicates erroneously inflate the coverage and, if not removed, can give the illusion of high confidence when it is not really there.

For example, consider the following hypothetical scenario.

```
                                         *
TTTCATACTAACTAGCCTGCGGTCTGTGTTTCCCGACTTCTGAGTCATGGGGTTTCAATGCCTATAGATTC
               ..........................C.
                ............................
                      .............T..............
                       ..................C.........
                           ............................
                               ............................
                                   ............................
                                    .....C......................
                                       ..C.........................

```

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

Now consider the following hypothetical scenario.

```
                                         *
TTTCATACTAACTAGCCTGCGGTCTGTGTTTCCCGACTTCTGAGTCATGGGGTTTCAATGCCTATAGATTC
                ............................
                      .............T..............
                       ..................C.........
                       ..................C.........
                       ..................C.........
                       ..................C.........
                               ............................
                               ............................
                               ............................

```

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

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

---
Tags: genome, ngs, genome-sequencing
---


# Question
Title: How to determine a protein's cellular location based on its sequence?
I am wondering about the appropriate workflow to determine a protein's cellular location based on its sequence.

Let's say I have a sequence like this from a fasta file

```
MTIKEHKVVYEAHNVKALKAPQHFYNSQPGKGYVTDMQHYQEMYQQSINEPEKFFDKMAK
EYLHWDAPYTKVQSGSLNNGDVAWFLNGKLNASYNCVDRHAFANPDKPALIYEADDESDN
KIITFGELLRKVSQIAGVLKSWGVKKGDTVAIYLPMIPEAVIAMLAVARIGAIHSVVFAG
FSAGSLKDRVVDANSKVVITCDEGKRGGKTINTKKIVDEGLNGVDLVSRILVFQRTGTEG
IPMKAGRDYWWHEEAAKQRTYLPPVSCDAEDPLFLLYTSGSTGSPKGVVHTTGGYLLGAA
LTTRYVFDIHPEDVLFTAGDVGWITGHTYALYGPLTLGTASIIFESTPAYPDYGRYWRII
QRHKATHFYVAPTALRLIKRVGEAEIAKYDTSSLRVLGSVGEPISPDLWEWYHEKVGNKN
CVICDTMWQTESGSHLIAPLAGAVPTKPGSATVPFFGINACIIDPVTGVELEGNDVEGVL
AVKSPWPSMARSVWNHHDRYMDTYLKPYPGHYFTGDGAGRDHDGYYWIRGRVDDVVNVSG
HRLSTSEIEASISNHENVSEAAVVGIPDELTGQTVVAYVSLKDGYLQNNATEGDAEHITP
DNLRRELILQVRGEIGPFASPKTIILVRDLPRTRSGKIMRRVLRKVASNEAEQLGDLTTL
ANPEVVPAIISAVENQFFSQKKK

```

I can feed it to NCBI's blastp which then gives me a list of genes. I can check e.g. the first entry which also provides a link to SwissProt. There I finally see that this protein is located in the nucleus and also the cytosol.

My questions are:

1) Is there a more straightforward way of achieving this information?

2) Let's say someone annotated this protein to be a mitochondrial one. How could I check that i.e. (dis)prove that?

# Answer
> 6 votes
You could use Blast2go to find the predicted relevant gene ontologies of your sequences. If you look for the cellular compartment sub ontology you should find the location of such protein. I have never used this tool, but here is the paper describing how it works.

I wouldn't restrict to the first predicted match of blastp see if some more genes with a good fit are also located at the same sites. But otherwise your workflow seems right to me.

But to prove this is in X or Y you, or your team or collaborators, should always check with histochemistry if it really allocates on the predicted place. All the softwares predictions must be validated with some biological experiments to prove them.

# Answer
> 2 votes
You can run the sequence through tools like SignalP to look for cell surface signals and so on, which will give you a clue about where I is bound for.

None of these tools will give you an absolutely definitive answer though, often because we simply don’t have the experimental data.

---
Tags: proteins, blastp
---


title: TruSeq Exome Kit bed file old version

text:
# Question
Title: TruSeq Exome Kit bed file old version
I have to run an analysis using Whole Exome Sequencing (WXS) data prepared with the TruSeq Exome Kit (2012).

I need the `.bed` file for the kit and I found it here.

However, this file was posted on Nov. 16, 2015, while the kit used in my case is from 2012. I searched but I didn't find the bedfile for that specific version online.

Will I run into problems if I use this version, and in this case, where can I retrieve the bedfile for 2012's version?

# Answer
> 0 votes
Usually it is more problematic to work with older versions of data with newer protocols than the other way around.

In this case, I think that the features annotated in the (newer) bed file should be more accurate than the older one. Unless you specifically want to compare your experiment with a previous experiment performed with a bed file from 2012 you'll be fine using the newer one.

However make sure that you are correctly using the right file and that you report (or record) each version of protocol and data used in case some doubts arise.

---
Tags: bed, exome
---


# Question
Title: TCGA bam file duplicates (aligned reads) and related last part of TCGA barcode
I am running an analysis on TCGA bam files (GCRh37/hg19) of whole exome sequencing data.

I usually find one bam file for each sample. However, at least one sample has two bam files, in which the last part of the name (before the extension) has the suffix `.1` or `.2` or `.3`...

For example, `C500.TCGA-BL-A3JM-01A-12D-A21A-08.1.bam` and `C500.TCGA-BL-A3JM-01A-12D-A21A-08.2.bam`.

I didn't find an explanation for this last suffix. Are they duplicates, or.. ? Is there some wiki about the last part of the filename?

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

---
Tags: sam, sequence-alignment, human-genome
---


# Question
Title: Why do I obtain very strange results with DESeq2?
I am using DESeq2 to perform a differential expression analysis, but I obtained very strange results for some genes. For some genes, I have very high log2FoldChange with very low p-value as displayed in the following table.

When I look the normalized expression values for these genes. I can see that these genes are expressed in only one single sample (condition1= samples from 1 to 7, condition2 = samples 8 to 13).

Could you please tell me how to avoid DESeq2 to detect these genes as differentially expressed ?

I am using to following commands:

```
dds                      <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory = directory, design= ~ condition)
dds                      <- DESeq(dds)
res                      <- results(dds, contrast=c("condition",condition1,condition2))

```

# Answer
> 2 votes
I'm guessing here, but I suspect this is at least in part due to the fact that all your samples in condition 1 are not expressed at all.

Formally you have an infinitely large fold-change between condition1 and condition2. Most tools solve this by adding a small pseudocount to unexpressed samples, but this will still result in a very high fold change. So far for the big fold-change.

If I had to guess at the reason why this is significant I would think in terms of dispersion estimation. This Bioconductor post has some discussion with the authors of DESeq2 on how this works.

Finally: calculating differential expression on these genes in these samples will never be really informative, as they really are only expressed in one sample. The way I see it you have two options:

1. Exclude the sample from your analysis. Straightforward, but does not solve this problem for other samples.
2. Filter your genes before calculating differential expression, i.e. remove all genes that are not expressed in a large number of samples. This approach is standard in for instance sleuth. A good starting point would be to remove genes that are not expressed in 75% of your samples.

---
Tags: differential-expression, statistics, deseq2
---


# Question
Title: How to align genomic sequence with corresponding amino acid sequence
Does anyone know of a program that can align a genomic sequence with introns with the corresponding amino acid sequence?

I have both the genomic sequence and the correct amino acid sequence but no information on the genemodel in e.g. gff or genbank format.

This is not what one would call a 'normal' alignment, but I would like it to look like this:

### Hypothetical example:

nucleotide sequence with intron in italic:

ATGCACGATACGACTGACGTACGT*ACGTACGTACGT*ACGTACGTACGGAGACGTAGACTC

corresponding amino acid sequence:

MHDTTDVRTYVRRRRL

The intron does not encode for amino acid which should thus result in a gap in the 'alignment'. Desired 'alignment':

```
ATGCACGATACGACTGACGTACGTACGTACGTACGTACGTACGTACGGAGACGTAGACTC
  M  H  D  T  T  D  V  R              T  Y  V  R  R  R  R L

```

# Answer
> 3 votes
Have a look at exonerate's protein2genome! From the documentation:

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

# Answer
> 3 votes
The classic tools for this are `genewise` and the already mentioned `exonerate`. Both are capable fo aligning protein sequences to genomic DNA (or cDNA) while modelling splice sites to give you a gene model.

For example, this is the output1 of `genewise` on your sequences:

```
EMBOSS_001         1 MHDTTDVRTYVR----RRRL
                     MHDTTDVRTYVR    RRRL
                     MHDTTDVRTYVRTYVRRRRL
EMBOSS_001         1 acgaaggcatgcatgcacac
                     taaccatgcatgcatggggt
                     gctgtcatgcatgcagatac

```

The output also includes a GFF of your sequence:

```
EMBOSS_001  GeneWise    match   1   60  22.88   +   .   EMBOSS_001-genewise-prediction-1
EMBOSS_001  GeneWise    cds 1   60  0.00    +   0   EMBOSS_001-genewise-prediction-1

```

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

---

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

---
Tags: sequence-alignment
---


# Question
Title: problem of "ordering in physical positions" phasing SNPs with Shapeit
I'm trying to phase my data (whole genome resequencing SNPs) with Shapeit (from ped file) and I get this error message :

```
[33mERROR:[0m Wrong ordering in physical positions curr_pos=54 prev_pos=65

```

I found that it's due to the fact that my .map file is like this (4 columns):

```
Chromosome_number    SNP_ID   SNP_genetic_position    SNP_physical_position
Contig100131            .              0                        65
Contig100220            .              0                        54

```

And Shapeit only looks at `SNP_physical_position` but does not check if the SNPs are on the same contig or not and so it says it is not happy because position 54 is after 65!

How could I make it understand that the positions are not in the same contig and so there is no mistake? Must I run Shapeit per contig (I have a lot of them so I hope not)?

# Answer
> 5 votes
The manual suggests that each "chromosome" needs to have its own input file.

> You have to split the dataset by chromosomes. ShapeIT can just phase one chromosome at a time. If you are working on GWA data, you can proceed as follows:
>
> > `for chr in $(seq 1 22) ; do plink --file example --chr $chr --recode --out example_chr$chr ; done`
>
> You will obtain 44 files; `example_chr1.ped, example_chr1.map, ..., example_chr22.ped, example_chr22.map` that you can phase separately.

---
Tags: ngs, software-recommendation, errors
---


# Question
Title: Visualization tools for 3C/Hi-C long-range interaction data?
I have identified a genomic region (CpG island) that is enriched for a characteristic of interest in the vast majority of my samples.

After viewing it in the UCSC genome browser, I have confirmed that it does not overlap any annotated gene promoters, gene bodies, or retrotransposons, but am convinced that it is functionally interacting with something long range.

Is there a publicly available visualization tool where I can input my genomic region and get a list of topologically associated domains from 3C or other chromatin-chromatin interaction data?

# Answer
I would also recommend two very recent Hi-C visualization frameworks (with some public data available in both): HiGlass and JuiceBox.
> 6 votes

# Answer
If you have more than a dozen or even hundreds of regions that you want to compare you could check a companion tool of HiGlass called **HiPiler**. HiPiler lets you arrange, cluster, and pile up to thousands of genomic regions. Take a look at this 5min video introduction to see if the tool could help you to answer your questions.

Disclaimer: I am the creator of HiPiler.
> 3 votes

# Answer
Depending on the genome you want, you may be able to use the chorogeneome navigator, which is from some of my colleagues. This uses a number of public HiC datasets and allows you to view interactions.
> 2 votes

---
Tags: software-recommendation, visualization
---


title: Use paste function with apply families

text:
# Question
Title: Use paste function with apply families
I am trying to assign different window sizes to my SNPs dataset to identify regions under selection.

this is the head of my data

```
head(snp_ids)
   snp_id   chr  pos
Chr01__912    1  912
Chr01__944    1  944
Chr01__1107   1 1107
Chr01__1118   1 1118
Chr01__1146   1 1146
Chr01__1160   1 1160

class(snp_ids)
data.frame

```

I have chosen 4 different window sizes, `win_size <- c(15000, 30000, 50000, 100000)`. I have assigned each of these different window sizes to my snp\_ids dataset to identify how many SNPs are distributed within each window by looping through each window size

```
for (i in 1:length(win_size)){
  windows <- sapply(snp_ids$pos, function(x) (ceiling(x/win_size[i])))
}

```

# Answer
Not clear what we are trying to do, but here is the solution using `lapply` that should replicate your existing `forloop`.

```
# example input data
snp_ids <- read.table(text = "
snp_id  chr  pos
Chr01__912   1  912
Chr01__944   1  944
Chr01__1107  1 1107
Chr01__1118  1 1118
Chr01__1146  1 1146
Chr01__1160  1 1160", header = TRUE)

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

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

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

```

---

**Edit:** If we really want to use `forloop`, then try below:

```
for (i in win_size){
  windows <- ceiling(snp_ids$pos/i)
  snp_ids[, paste0("window_", i)] <- as.numeric(paste(snp_ids$chr, windows, sep = "."))
}

```
> 4 votes

# Answer
This explicitly doesn't use `apply()` or `lapply()`, but rather the `dplyr` package. This is probably easier to read and debug:

```
library("dplyr")
snp_ids %>% mutate(window_15000=as.numeric(paste(chr, ceiling(pos/15000), sep=".")),
                   window_30000=as.numeric(paste(chr, ceiling(pos/30000), sep=".")),
                   window_50000=as.numeric(paste(chr, ceiling(pos/50000), sep=".")),
                   window_100000=as.numeric(paste(chr, ceiling(pos/100000), sep=".")))

```
> 2 votes

---
Tags: r
---


title: bcftools output vs. bgzip

text:
# Question
Title: bcftools output vs. bgzip
I've got what seems like it should be a simple question, but I can't seem to figure it out from google. bcftools has 4 output options:

```
Output compressed BCF (b), uncompressed BCF (u), compressed VCF (z), uncompressed VCF (v).

```

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

# Answer
> 5 votes
Yes. I haven't found it documented anywhere, but I ran a simple test:

```
$ bcftools view -Oz --no-version foo.vcf > bar.vcf.gz
$ md5sum bar.vcf.gz
18a2f494710cf170ba79892936015ba8  bar.vcf.gz

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

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

```

I passed a vcf file through `bcftools view -Oz`, checked the `md5sum` of the resulting file, then uncompressed it and recompressed with `bgzip` and confirmed that the resulting file had the same `md5sum`. Finally, I also compressed using normal `gzip` to be sure that the `md5sum` would be different, as it was.

So yes, it does indeed look like `bcftools` creates files that are the same as though you'd passed them through `bgzip`.

# Answer
> 2 votes
Yes the compressed outputs should be the same format as bgzip. bcftools also has an option to index in either the CSI or TBI formats.

The index function mentions that it indexes bgzip compressed files. So one can assume that it also creates files using the bgzip compression. Also all the pipelines that I have built use bgzip and bcftools.

```
[~]# bcftools index

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

```

---
Tags: file-formats, vcf, htslib
---


# Question
Title: What are Approximate Read Counts (Library Sizes) and Lengths (Insert Sizes) for Next-Generation DNA Sequencers?
I am performing a simulation study and am curious about the parameters of my simulated metagenome.

What are the library and insert sizes of some of the most "used" sequencers. Mainly, I am interested in the different Illumina platforms (MiSeq, HiSeq, HiSeqX), 454 Pyrosequencing, and Ion Torrent.

I do not need the long read technologies like Pac-Bio's SMRT and Oxford Nanopore since they are not typically used with metagenomics, yet.

# Answer
The library size depends usually on your budget. There are some reviews where they compare the different platforms, but these papers get outdated so quickly, that it would be better to just check wiki for the specs of different platforms. Also on the websites of Illumina or Life Tech, you'll find the specs.

NB. 454 sequencing is discontinued, so that is not very relevant anymore in my opinion.
> 0 votes

# Answer
Both of those depend on the application of the library. For example whole genome assembly of a mammal would require a larger library size than targeted re-sequencing of a virus. Also insert size varies if you are doing genome assembly a larger insert size is valuable to allow for long range linking and scaffolding.

I think that if you edit your question of be more specific people will be able to provided a more concrete answer other wise it is "it depends"
> 1 votes

---
Tags: ngs, metagenome, reads, illumina
---


# Question
Title: Using PLINK to find SNPs in LD (Linkage Disequilibrium) with another set of SNPs
I have a list of SNPs (index SNPs). As my tagging SNPs, I want to use SNP data from here. The data has the following column headers: Chromosome, Position, MarkerName, Effect\_allele, Non\_Effect\_allele, Beta, SE, Pvalue.

My ultimate aim is to identify SNPs in the tagging set that are in LD with each of index SNPs. I need to calculate r2 between index SNPs and tagging SNPs within 1MB window and get a table of index vs tagging SNPs (r2\>0.8).

How can I use PLINK to carry out this work? I'm new to PLINK and LD clumping

# Answer
> 1 votes
You can't generate LD statistics for summary data from SNPs, as the LD calculation relies on the comparison of variants at the genotype level. Plink expects data to be in a genotype matrix with either each marker for each line (probably better in your case), or each individual for each line. Plink v1.9, in particular, will accept VCF-formatted input files, when populated with 'GT' information that describes the variants for each individual.

If you're okay with using the 1000 genome population for calculating LD, plink v1.9-compatible datasets should be accessible from here.

If not, you'll need to discuss your project with IGAP to see if you can have access to genotype-level data for LD testing. That's unlikely, because there are personal identification issues associated with non-summary data.

There's comprehensive documentation about the LD command(s) on the Plink 1.9 website here.

---
Tags: snp, plink
---


# Question
Title: Integration of different microarray dataset to run GSEA
I'm planning to run a GSE Analysis on some microarray datasets. Until now what I have done is to pre-processed them using frma algorithm.

The reason I chose frma instead of gcrma or rma is that I wanted to merge these different dataset into one.

Here is how I'm planing to proceed:

i.e

```
dataset1 : 3x Controls , 3x Treatment1
dataset2 : 2x Controls , 2x Treatment2
dataset3 : 2x Controls , 2x Treatment3

```

After merging

```
7x Controls , 3x Treatment1 , 2x Treatment2 , 2x Treatment3

```

Then I will execute the GSEA for controls , Treatmen1 , Treatment2 and Treatment3 separably.

So the question is if it's more reasonable to follow the path that I described or to run the GSEA on the separate datasets (for example for each control ) and then get the intersect of the gene sets of these three repetitions.

# Answer
It makes more sense to evaluate by separate the pathway or gene set you want and see if in the three datasets result in a coherent message than to merge these datasets, as you will mix batch effects and treatment effects.

When in several datasets one looks for the same statistic it is usually considered meta analysis, not integration.
> 2 votes

---
Tags: microarray, gse
---


# Question
Title: Why total RNA-seq usually yields low mapping rate?
Maybe this is a silly question, but I'm really wondering why we usually get low mapping rates if we map total RNA-seq, but not poly(A)-enriched (in particular, for human, mouse, and zebrafish datasets)?

Doesn't the genome fasta file contain ribosomal RNAs (which are abundantly expected in the total RNA-seq libraries) as well?

I don't have the numbers for human right now, but I recall previously having low mapping rates for human data as well. What triggered me at the moment is actually when I mapped zebrafish data. I got around 60% mapping to the genome and 45% to the transcriptome with one dataset and 36% - 35% with another dataset (note that these are from 2 different studies and both are total RNA-seq, mapping was done using STAR and Salmon, which have their own methods).

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

# Answer
> 3 votes
In particular for mouse, this can happen if the ribo-depletion wasn't terribly efficient, since there's no Rn45s sequence in the reference genome. That, combined with numerous copies of tRNAs, 5S rRNAs, etc. causing issues with multimappers (see the answer from @Tom) can heavily decrease alignment rates.

Note that the human reference genome contains a few copies of the 45S (e.g., on GL000220.1 and chr21), so these reads will all multimap (possibly too many times). I don't know about the zebrafish reference genome.

# Answer
> 2 votes
When sequencing RNA from freshly harvested tissues under proper conditions, you should generally expect \> 50% mapped reads. In fact, everything \< 80% would usually raise concerns.

From your description (in the question and comments) it sounds as if your samples are potentially degraded and therefore saturated with short RNA fragments, either because you’ve got old tissue, you’re looking at environmental samples, or from improper handling.

And since you’re not performing any enrichment, the sequencing data will then be highly saturated with degraded RNA, too. The reason for this is that RNA-seq is a (hopefully somewhat uniform) sampling of your sample aliquot: everything in it is proportionally represented in the RNA-seq data. That’s why enrichment steps are important, to increase the proportion of whatever fraction we’re interested in.

In particular,

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

In my experience, the second issue is exacerbated on certain sequencers (I’ve seen it on an Illumina HiSeq 1500) because a mixed library of short and long RNA doesn’t seem to cycle well through the whole fragment length. So even if you have a (small but still present) fraction of long RNA fragments, the sequencer may not be able to synthesise their whole length efficiently. The effect is that there are proportionally more short reads than there were short fragments in the sample.

You can verify this by calculating the insert size distribution of your reads after adapter trimming (e.g. using picard, or by simply tallying read lengths1). Small reads (\< ~14 nt) from degradation fragments are essentially unmappable, since they are too short to have a specific identity. Even if we allowed the mapper to align them their coordinates would be essentially random.

---

1 This works:

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

```

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

---
Tags: rna-seq, read-mapping, mappability
---


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

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

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

# Answer
Here’s a function that I keep around for calculating the individual lengths of fasta sequences

```
falens(){
awk '/^>/ {if (seqlen){print seqlen}; print ;seqlen=0;next; } { seqlen += length($0)}END{print seqlen}' $1
}

```

And one I use for summing a column of numbers

```
sum(){
awk '{ sum += $1 } END { print sum}'
}

```

You can combine these by piping the output to do the sum for you:

```
falens myseqs.fa | grep -v “>” | awk '{ sum += $1 } END { print sum}'

```

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

---
Tags: fasta
---


# Question
Title: How can I find out the code DRIMSeq uses to make dmplotProportions objects?
I'm using DRIMSeq to do some analysis but I'd like to make my own figures from the results ideally using what they've done as a base.

I'm trying to find how they make their ggplot objects in the function `plotProportoins` but when I use `findMethods('plotProportoins')` it ends up with a method `dmPlotProportoins`.

It says in the manual that it's a *ggplot* object but I can't seem to find how they make it. Any ideas?

# Answer
I was curious how they have done it as well, so here expanded @Liopis comment :

```
wget http://bioconductor.org/packages/release/bioc/src/contrib/DRIMSeq_1.6.0.tar.gz
tar -zxf DRIMSeq_1.6.0.tar.gz
# if you need to find it, or use grep
# find DRIMSeq -name "dm_plotProportions.R"
less DRIMSeq/R/dm_plotProportions.R

```
> 4 votes

---
Tags: rna-seq, r, rna-splicing, ggplot2
---


# Question
Title: Automating sequence download from Wormbase
I'm working on the *C. elegans* model organism, for which there exists an online biological information database: Wormbase.

I'm able to manually search for a transcript (say "K06C4.12") and this gets me the following URL: http://www.wormbase.org/species/c\_elegans/cds/K06C4.12

Say I have this transcript ID: "F08G2.2". I guess that I can form a valid similar URL, as follows: http://www.wormbase.org/species/c\_elegans/cds/F08G2.2

Indeed, it works as expected.

Now, I would like to go further and be able to download in an automated way what can be obtained by clicking on one of the "download" buttons in the "Sequence" section of the above pages. More specifically, I'm interested in the "spliced + UTR" and "conceptual translation" fasta sequences.

If I use the DOM inspector in Firefox (Shift-Ctrl-C), I can locate things like the following:

```
<div class="button ui-corner-all generate-file-download sequence-download">
    <span class="ui-icon ui-icon-arrowthickstop-1-s"></span>
    <span>download</span>
    <div id="content" class="ui-helper-hidden">
        >K06C4.12 spliced + UTR
        ATGGCCCCACCAAAGCCGTCTGCTAAAGGAGCAAAGAAAGCCGCCAAGACCGTCACCAAGCCAAAGGACGGAAAGAAGAG
        ACGTCATGCCCGTAAAGAATCCTACTCCGTGTACATCTACCGTGTGCTCAAGCAGGTTCATCCAGACACTGGAGTCTCTT
        CCAAGGCTATGTCCATCATGAACTCTTTCGTCAACGACGTCTTCGAGCGTATCGCTGCTGAGGCATCTCGTCTTGCTCAC
        TACAACAAGAGATCCACCATCTCTTCTCGCGAAATTCAAACTGCCGTCCGTTTGATTCTTCCAGGAGAGCTTGCCAAGCA
        CGCCGTCTCTGAGGGAACCAAGGCTGTCACCAAGTACACCTCTAGCAAGTAA
    </div>
    <div id="filename" class="ui-helper-hidden">K06C4.12conceptual translation.fasta</div>
    <div id="content" class="ui-helper-hidden">
        >conceptual translation
        MAPPKPSAKGAKKAAKTVTKPKDGKKRRHARKESYSVYIYRVLKQVHPDTGVSSKAMSIMNSFVNDVFERIAAEASRLAH
        YNKRSTISSREIQTAVRLILPGELAKHAVSEGTKAVTKYTSSK*
    </div>
</div>
<div class="toggle ">
    <span style="float:left" class="ui-icon ui-icon-triangle-1-s"></span>
    view conceptual translation (123 aa)
</div>
<div class="returned sequence-container" style="display:block">
    >conceptual translation
    <br>
    MAPPKPSAKGAKKAAKTVTKPKDGKKRRHARKESYSVYIYRVLKQVHPDTGVSSKAMSIMNSFVNDVFERIAAEASRLAH
    <br>
    YNKRSTISSREIQTAVRLILPGELAKHAVSEGTKAVTKYTSSK*
</div>

```

I don't find these if I simply download the page with `wget` and search in the source code, so they must be dynamically generated.

**How could I proceed to get the sequences?**

The documentation about Wormbases's REST API is not very detailed, so I don't understand how to form an URL that would give me the desired sequences.

# Answer
Further inspecting the source code of the downloaded page, I found the following:

```
<li id="nav-sequences"
    class="module-load transcript"
    wname="sequences"
    href="/rest/widget/transcript/K06C4.12/sequences">
    <span class="ui-icon ui-icon-close module-close" wname="sequences"
    </span>
    Sequences
</li>

```

This suggested me to have a look at http://www.wormbase.org/rest/widget/transcript/K06C4.12/sequences and I obtained a page containing the elements found in the DOM inspector.

This can be parsed, for instance using the `bs4` python module:

```
# To get the page
import requests
# To parse the page content
from bs4 import BeautifulSoup
# To get this into Biopython
from io import StringIO
from Bio import SeqIO

rest_template = "http://www.wormbase.org/rest/widget/transcript/{accession}/sequences"
rest = rest_template.format(accession="K06C4.12")
# This did not work with "html.parser"
soup = BeautifulSoup(requests.get(rest).text, "html5lib")

def generate_fasta(soup):
    name = None
    for div in soup.find_all("div"):
        if "ui-helper-hidden" in div["class"]:
            if div["id"] == "content":
                fasta = SeqIO.read(StringIO(div.text), "fasta")
                description_fields = fasta.description.split()
                if len(description_fields) == 4:
                    fasta.id = "_".join(description_fields[:2])
                    if name is None:
                        name = fasta.name
                elif len(description_fields) == 2:
                    # This assumes that the "conceptual_translation"
                    # is after the "unspliced + UTR" or "spliced + UTR"
                    # Otherwise, one would need to pass
                    # the accession as function parameter
                    fasta.id = f"{name}_translation"
                else:
                    fasta.id = "_".join(description_fields)
                yield (fasta.id, fasta)

records_dict = dict(generate_fasta(soup))
print(records_dict['K06C4.12_spliced'])
print(records_dict['K06C4.12_translation'])

```

Output:

```
ID: K06C4.12_spliced
Name: K06C4.12
Description: K06C4.12 spliced + UTR
Number of features: 0
Seq('ATGGCCCCACCAAAGCCGTCTGCTAAAGGAGCAAAGAAAGCCGCCAAGACCGTC...TAA', SingleLetterAlphabet())
ID: K06C4.12_translation
Name: conceptual
Description: conceptual translation
Number of features: 0
Seq('MAPPKPSAKGAKKAAKTVTKPKDGKKRRHARKESYSVYIYRVLKQVHPDTGVSS...SK*', SingleLetterAlphabet())

```
> 4 votes

---
Tags: python, public-databases, data-retrieval, api
---


# Question
Title: Estimating BAM file from compressed fastq file size
Is there a way to estimate the size of a BAM file will have after mapping with BWA?

The input file are two mates fastq files, compressed with gzip, each one about 70G.

# Answer
> 4 votes
A BAM file should have a roughly similar size to a compressed FASTQ file, because they're both compressed files that contain similar information. BAM will be slightly larger, because it contains information about mapping and differences from the reference.

If you are concerned about space and have a reference sequence available, you can store mapped reads as CRAM format, which will save a lot more space.

---
Tags: sam, fastq, bwa
---


title: Identifying relevant SNPs from a list

text:
# Question
Title: Identifying relevant SNPs from a list
I have a list of all SNPs (a list of RsIDs) that fall in the coding sequence of several thousand human genes, and I'm looking for a convenient method to come up with a list of known SNPs that alter an amino acid in those genes.

I'll be glad to hear your thoughts.

# Answer
What you are looking for is SNP annotation. If you have the chromosome:position reference and alternate alleles for your SNPs of interest, it can be as simple as uploading them to the variant effect predictor.

http://grch37.ensembl.org/Homo\_sapiens/Tools/VEP

This will give you the predicted protein change and novelty of the variant with respect to known variation in large scale sequencing projects like the 1000 genomes.

If you cannot afford to upload data on the web, there are downloadable command line software like SNPeff, Variant effect predictor etc but it would require some more bioinformatics knowledge of converting your data to a VCF format and obtaining a human reference genome.
> 5 votes

# Answer
ANNOVAR is another tool that will help you functionally annotate genetic variants. It will tell you if the variant is known in COSMIC or Clinvar databases, which AA change will occur, if it's a frameshift, synonymous, non-synonymous, etc...
> 3 votes

---
Tags: database, gene, snp, sequence-analysis
---


# Question
Title: dividing genome into non-overlapping windows using R
I have nearly 10 million SNPs located on 10 chromosomes. I want to divide the genome into non-overlapping windows of 15, 20 and 30 kb. Here is part of my SNP table:

```
head (sap_ids)
       snp_id    chr  pos
    Chr01__15043   1 15043
    Chr01__15079   1 15079
    Chr01__15139   1 15139
    Chr01__15165   1 15165
    ...
    ...
    ...
    Chr17__214708424  17 214708424
    Chr17__214708451  17 214708451
    Chr17__214708484  17 214708484
    Chr17__214708508  17 214708508
    Chr17__214708574  17 214708574

```

I have been using this code, but it gives me wrong output because each chromosome starts from position 1

```
win_size<-c(15000,30000,50000,100000)
res<- cbind(snp_ids,
       data.frame(lapply(setNames(win_size, paste("window",win_size, sep = "_")),
       function(x)as.numeric(ceiling(snp_ids$pos/x)))))

```

For instance, as I you see in the example below in window 4 I get SNPs from chromosomes 1, 3 and 5:

```
snp_id     chrom poistion   window
Chr01__58332   1 58332            4
Chr01__58335   1 58335            4
Chr01__58341   1 58341            4
Chr01__58450   1 58450            4
Chr01__58471   1 58471            4
Chr01__58530   1 58530            4
Chr01__58542   1 58542            4
Chr01__58641   1 58641            4
Chr03__45457   3 45457            4
Chr03__45604   3 45604            4
Chr04__56873   4 56873            4
Chr04__57387   4 57387            4
Chr04__57399   4 57399            4
Chr04__57528   4 57528            4
Chr04__58419   4 58419            4
Chr04__59670   4 59670            4
Chr04__59704   4 59704            4

```

I guess the solution would be to run the code above by looping through each chromosome I tried something like this:

```
for (i in 1:10){
       res <- cbind(cats,
             data.frame(
               lapply(setNames(win_size, paste("window", win_size, sep = "_")), function(w)
                 paste(cats$chr[cats$chr == i], ceiling(cats$pos/w), sep = "_"))
               ))
       }

```

But still not working appropriately (column 4 should be 1\_1, for the chr 1)!

```
Chr01__912    1  912          2_1
Chr01__944    1  944          2_1
Chr01__1107   1  1107         2_1
Chr01__1118   1  1118         2_1

```

Can anyone help me to figure out how can I modify the code to create these non-overlapping windows properly?

# Answer
The question is a bit confusing to me, at it's core I understand you want non-overlapping windows across chromosomes.

One way to achieve is to use `GenomicRanges::tileGenome` function, which needs the chromosome lengths as input, e.g.

```
library(Biostrings)
mygenome <- readDNAStringSet(list.files(mypath,"mygenome.fa$",full=TRUE))
chrSizes <- width(mygenome)
names(chrSizes) <- names(mygenome)
print(chrSizes)
#     Chr1     Chr2     Chr3     Chr4     Chr5
# 18585056 19698289 23459830 26975502 30427671

```

This can than be fed into tileGenomes:

```
library(GenomicRanges)
bins   <- tileGenome(chrSizes, tilewidth=5e5, cut.last.tile.in.chrom=T)

```

This creates a GRange object which gives you a whole repertoire of post processing options, i.e. overlapping with SNPs or NGS reads (I use the latter for genome wide coverage plots):

```
print(bins)
# GRanges object with 240 ranges and 0 metadata columns:
#         seqnames               ranges strand
#            <Rle>            <IRanges>  <Rle>
#     [1]     Chr1   [      1,  500000]      *
#     [2]     Chr1   [ 500001, 1000000]      *
#     [3]     Chr1   [1000001, 1500000]      *
#     [4]     Chr1   [1500001, 2000000]      *
#     [5]     Chr1   [2000001, 2500000]      *
#     ...      ...                  ...    ...
#   [236]     Chr5 [28000001, 28500000]      *
#   [237]     Chr5 [28500001, 29000000]      *
#   [238]     Chr5 [29000001, 29500000]      *
#   [239]     Chr5 [29500001, 30000000]      *
#   [240]     Chr5 [30000001, 30427671]      *
#   -------
#   seqinfo: 5 sequences from an unspecified genome

```
> 7 votes

---
Tags: r, genome
---


# Question
Title: How can I download genePred file from UCSC?
I need a gtf annotation file. UCSC doesn't give us a proper gtf file with distinct gene\_id and transcript\_id. It asked us to get a genePred file to convert to gtf.

But where can we get genePred files?

# Answer
It turned out that using table browser on UCSC we could download an output with all fields. That output actually contains a GenePred file but has one extra column on the left end.

The genePred file format definition is here.
> 1 votes

---
Tags: ucsc
---


# Question
Title: EMT/EMT-like processes in bioinformatics to study cancer progression
First I am not sure why, but it seems EMT can be used to infer the status of cancer progression.

What kind of databases are there for the EMT process in the context of cancer?

I saw this book which mentions ArrayExpress, but I need some help to understand what kind of data it stores.

Are there similar/better resources? How can I use EMT in the context of cancer progression study?

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

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

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

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

ArrayExpress is a database of transcriptomics, it :

> stores data from high-throughput functional genomics experiments

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

# Answer
> 1 votes
The epithelial-mesenchymal is a process by which cells lose some of the structures specific to their cell type. This is normal during embryonic development or wound healing for the formation of new cells. However, it is very uncommon in healthy adults and is thereby associated with the cancers in which it occurs. EMT allows cancer cells to become more like stem cells and could facilitate metastasis or immortal cell growth, key processes for cancer progression.

If you can detect the molecular components involved in EMT, you may be able to make some inferences about the progression of the patient’s disease. This process is specific to cancers in epithelial tissues (breast, stomach, and skin cancers for example). It is just one of vastly many pathways by which cancers can grow and develop. Many of these are specific to different cancers types by tissue or molecular subtype. Cancer research is therefore very challenging and requires a deep understanding of the particular system you are focusing on.

Gene expression data sets are a great resource for this. You can compare many pathways to clinical properties from previous studies. Microarray datasets such as ArrayExpress and Gene Expression Omnibus (GEO) were very popular for this but are now largely out of date. Many databases now support higher quality RNASeq datasets with larger sample sizes and fewer batch effects. The Cancer Genome Atlas (TCGA) and International Cancer Genome Consortium (ICGC) have among the largest Cancer databases. These are a tool for researchers to investigate any genes or pathways. They will not specifically address EMT, unless you focus on the genes involved in a relevant cancer type.

---
Tags: cancer, pathway
---


# Question
Title: Normalizing microarray data for clustering heat map
I wanted to generate a clustering heat map for the microarray data. This is the first time I'm working on Microarray data. I read some tutorials but have few doubts.

I'm using microarray (Affymetrix SNP 6.0 data) gene expression data. For example the data looks like following:

```
ProbeID S1  S2  S3  S4  S5  S6  S7
10008   131.4311    369.4926    222.0441    687.4181    176.8892    258.1233    316.5573
10010   78.73022    83.97501    81.56039    86.11443    78.09758    81.88231    84.17101
10014   90.02816    95.07267    101.1761    93.35585    81.96468    94.3553 93.89527
10017   79.86837    81.63064    88.19524    79.47265    76.3437 101.6351    93.71674
10019   99.03493    109.1835    104.97  102.7423    108.3677    98.93459    101.4052
10020   79.58075    84.28915    90.53562    74.47786    75.96112    96.39649    95.8828
10021   121.5373    149.9351    146.5956    122.8523    110.5759    132.4268    130.4409
10025   616.5994    1326.735    1358.187    2315.851    1068.745    3229.759    4435.021
10035   70.44073    69.56772    68.25446    68.35857    70.86771    74.3843 67.93569

```

I created expression matrix and did log2 conversion. Now the data looks like following.

```
            S1  S2  S3  S4  S5  S6  S7
Gene1   6.429276339 6.338451158 6.333760753 6.419191996 6.503471181 6.329103499 6.211373601
Gene2   6.379471993 6.296018518 6.237316465 6.2696332   6.329489132 6.359770303 6.240070336
Gene3   12.84498365 12.00265682 13.92741965 12.553162   13.39001307 13.8933423  12.58695704
Gene4   8.661860382 6.723004202 6.300975176 7.661012019 7.905219709 6.957578023 6.70945883
Gene5   6.945382967 6.814979733 6.575916303 6.63460198  6.627380524 6.733926424 6.280618235
Gene6   9.280222581 8.81560969  9.683073561 9.480038673 8.801438707 9.190943942 7.743705471
Gene7   6.593564871 6.63502488  6.389962535 6.511360029 6.694572404 6.6492763   6.527199544
Gene8   6.431615372 6.309515078 6.248153876 6.329288965 6.46768078  6.355268547 6.384284754
Gene9   7.513349406 7.234654595 7.490935892 6.801368215 7.323811386 7.018733196 7.055932044

```

Before clustering for normalisation:

I applied a function for normalisation

```
data = t(apply(data, 1, function(x) {
  q10 = quantile(x, 0.1)
  q90 = quantile(x, 0.9)
  x[x < q10] = q10
  x[x > q90] = q90
  scale(x)
}))

```

When I apply the above function I see there are values which are positive and also negative.

And when I used "rma" for background correction and "normalise" function from "oligo" package I don't see any negative values.

I don't understand.

Which one should I use to check the clustering among samples? Any help is appreciated.

# Answer
> 5 votes
You see negative values with your function because you're setting the average of each row to 0 and its standard deviation to 1.

In general, I would trust a standard normalization method (rma in this case) more than some random "truncate and then scale the rows" method. Your method isn't even doing any between-array normalization, which is the benefit of rma.

# Answer
> 1 votes
Negative values are to be expected as you have log-transformed the data. Any data value below 2 will have a negative log2.

There is a difference between normalising for analysis and truncating for plotting. The RMA normalisation or similar is intended for differential expression analysis. I would recommend using this normalisation first.

When plotting, the normalised data may still contain outliers which skew the plotting parameters (such as the colour scale). Due to the large number of genes, there may be some far beyond 3 standard deviations and plotting these at their true value will make it hard to discern trends in the middle. Truncating the data for a heatmap is common, typically at |z| \> 3 but quantiles should work much the same. This is a matter of visualisation, not clustering analysis. It would not be needed unless you were plotting the data.

---
Tags: r, bioconductor, normalization, microarray
---


# Question
Title: Gene function annotation - bacterial genome
I'm trying to annotate a genome to find all genes with a specific function. I have a FASTA and the read FASTQs - I'd like to assign the functional group of the identified proteins (e.g. Kegg orthology) automatically.

For more context, I have whole genome sequences of bacteria - I'm trying to find all genes involved in membrane transport. I'd either like a tool to retrieve gene ontologies from a list of gene names (input: genA output: Membrane Transport), or a tool that outputs the function as part of the annotation pipeline (e.g. genbank output with GO function defined).

Or posed in a different way, how would you identify all genes with a specific function from prokaryotic whole genome sequence data?

Are there any scripts that can do this?

# Answer
The PATRICBrc.org is a free resource you can annotate your genome and compare it with publicly available genomes.

Refer this tutorial for more detail.

Once you have annotated your genome in PATRIC, you will have PATRIC Local families (k-mer based protein family), which will help you identify all genes with a specific function.
> 4 votes

# Answer
Im not 100% on what you are looking for - you seem to jump from annotation to detection of genes in your question.

But if you have a database of genes you want to detect (presence/ absence) you could try srst2 which basically takes your raw reads (fastq) and maps them to the database. it will give you a report on any hits, and any snps/ indels in the hit. The only problem is making you own database is difficult as you must adhere to the fasta header format.
> 1 votes

---
Tags: bacteria, sequence-annotation
---


# Question
Title: How to get Nanopore MinION fast5 from SRA
I found some Nanopore MinION data on SRA, which I would like to investigate. I use `sratoolkit` for Illumina data all the time, but I am not sure how to get the `fast5` file from the `sra` file.

What I have tried is to `prefecth` it, but then I don't know which command to use for `.sra` to `.fast5` conversion. Not with fastq-dump, I assume?

```
prefetch -v SRR5286963 --max-size 50000000

```

# Answer
When I was trying to upload MinION data to SRA, they wouldn't allow me to upload only the raw signal data (i.e. FAST5 files). Every file needed to contain a FASTQ sequence, and this sequence was extracted from the file and stored in the archive. Given this, I don't expect it will be possible to retrieve signal or event information from SRA. I also don't expect this will be possible in the future either, given that the SRA database expects every record to have a FASTQ-like sequence attached.

However, EBI/ENA *do* allow the upload of raw FAST5 files to the database (e.g. here), but there are no FAST5 reads for the specific single-cell B1a study you have indicated. If you're interested specifically in that study, then I recommend that you ask Miten Jain directly for access to those reads; he should have another method to give you access to them.

ENA is "ready" for ONT reads, but researchers need to upload their reads to ENA in order for that database to be useful.
> 3 votes

# Answer
I would suggest trying an alternative way, as the FAST5 is a fairly new entry in SRA - 6 public datasets as of today.

You can search the ERX/ERR entry you got from SRA in ENA and either (1) direct download it via ftp or (2) redirect to Galaxy.

Here's a tutorial that might be interesting for you: basically, you first get the archive from ENA, then find the FAST5 files in the top level of the downloaded archive, all to be processed in R using `Fast5Summary` objects from IONiseR package.
> 2 votes

---
Tags: nanopore, fast5, sratoolkit
---


title: Scaffolding a genome with hybrid data

text:
# Question
Title: Scaffolding a genome with hybrid data
I am assembling a ~500MB genome, and have ~150x long reads and ~200x 150bp PE short reads, with a ~400bp insert size.

I've done a lot of work with minimap+miniasm, and have what I think is a good set of unitigs and a .gfa file from miniasm. Depending on the parameters, I have a few hundred to ~1K unitigs. I can easily polish these unitigs with pilon, racon, and nanopolish, so that they are as error free as possible.

The question I have is: given access to

1. Polished (or not) unitigs
2. 150x long reads (long enough to span almost all repeats in the genome)
3. 200x 150bp PE short reads (~400bp insert size)

What would you do to scaffold a genome? And what are the pros and cons of your suggested approach?

What I need here is a draft genome. I would prefer a conservative approach that is more fragmented and has fewer misassemblies. So of course one option is just to keep the unitigs and make no attempt at scaffolding. But it does strike me there may be suitably conservative approaches I haven't considered. E.g. perhaps I could use ABYSS to scaffold my polished unitigs using my short-read data, thereby leveraging whatever information I can from the PE short reads?

# Answer
> 2 votes
Try LR\_Gapcloser.

I've used L\_RNA\_Scaffolder for trying to scaffold a genome (which turned out to be more complex than I had expected). It looks like LR\_Gapcloser (written by the same people) is similar, but designed specifically for scaffolding using long-reads.

That page also suggests PBJelly and GMCloser as competing tools.

---
Tags: assembly, long-reads, illumina, scaffold
---


# Question
Title: differential expression multiple studies
I was wondering if we can perform differential expression on samples from different library and read lengths.

What we intend to do is,

1. n = 6 tumors (recently sequenced using Hiseq 2500 read length 150, paired end)
2. n = 5 normals (previously sequenced using Hiseq read length 100, paired end, these are not matched normal)
3. n = 100 ( normals, from tcga read length 100, paired)

What we would like to do is compare n=6 Vs 100 + 5 normals. What steps should we be considering before we move ahead with this analysis?

# Answer
> 2 votes
Adding experiments 1 & 2 should be fine without modification on the computer side of things, as long as your sample preparation is consistent. To be ultra-cautious, you could add in the experiment as an additional covariate for the differential expression model.

I'm less confident about adding in TCGA data; it's probably going to look better than your own samples, which might cause problems in calculating differential expression from raw reads. You can do a separate normalisation/correction of the TCGA data, but that will limit the inferences that can be made about differential expression -- without raw count information, the shot noise associated with read sampling cannot be included in the statistical model. I notice that TCGA has already done their own filtering / expression, and that would be preferable to use (mRNA sequencing, data subtype 'Gene'):

https://tcga-data.nci.nih.gov/docs/publications/tcga/datatype.html

However, even with data that has been properly normalised and filtered to exclude batch effects, there will still be population differences that influence results. My recommendation would be to do some subsampling to reduce the population effects, something like at least the following comparisons:

* 6 Vs 100 + 5 normals
* 6 Vs 100 normals
* 6 Vs 50 + 3 normals \[randomly sampled from the two normal groups\]

And only consider differential expression that is consistently observed in all three of these comparisons. Note that the third "subsampled" dataset can be generated more than once, taking different subgroups each time.

If you are going to use separately normalised data, then DESeq2 may not be appropiate (because it expects count-level data). Something like Limma could work for DE analysis, see chapter 15 of the user's guide for RNASeq analysis. Unfortunately I'm not very familiar with using Limma, so can't provide any advice on things like what to watch out for in the input data and/or results.

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


# Question
Title: Is there a more elegant solution to the bwa-mem: paired reads have different names error?
I'm currently trying to run bwa-mem on Influenza substrains using the following command:

```
~/bwa mem h5n1_1_cons.fa h5n1_1_read1.fq h5n1_1_read2.fq

```

h5n1\_1\_cons.fa is the consensus sequence for substrain h5n1\_1, and the fq files are paired end reads 1 and 2. Running the program gives me this error:

```
[mem_sam_pe] paired reads have different names: "JQ714243-920-1", "JQ714243-920-2"

```

Visual inspection shows that the read names are identical until the last character. I have found a package (bbmap) with a utility, bbrename.sh, that renames the sequences, producing:

```
./bbrename.sh in=h5n1_1_read1.fq out=fixed prefix=h5n1_1

```

which yields the following output

```
-bash-4.2$ head fixh5n1_1_read1.fq
@h5n1_0
gaaaatgtctcaggagttgatgcatggtgtttaatctttggtttgcgctg
+
@4>C5?7FA?9@@@;1,-4B+E??4E9,DA>6CB)*2.)?7/.'65+75/
@h5n1_1
catcagttggattttgcctaaggatgtccaccatccttatcccaccaatt
+
7(9?)=C3AEF+@@F=C-E>EE57CD:;'?.*?7;>:3>2.;29..)=2'
@h5n1_2
catgcatttgttggaataatgttctcacaaatcaactgtattgaccgctg

```

This fixes the problem and allows BWA-MEM to work. However, header information is discarded. Is there another way to ensure the header names are the same while retaining all header information?

# Answer
> 7 votes
Assuming the read name lines only contain things like `@JQ714243-920-1` and nothing else, then `sed -i '1~4 s/-[12]$//' file.fastq` (N.B., example updated a bit for safety, thanks to @terdon) will remove the `-1` and `-2` suffixes.

However, you should try to figure out **why** the read names were munged into that form to begin with (at least Illumina machines won't produce that, I don't know about others).

---
Tags: bwa, errors
---


title: Joining scaffolds to get full genome?

text:
# Question
Title: Joining scaffolds to get full genome?
I am using Ray to assemble a about 7.2 kb picornavirus genome. I am using k-mer of 55. Following is the exact code that I used.

```
mpiexec -n 1 Ray -k 55 -s 2-CA-SVA-2_S2_L001_R1_001.fastq_q30_trim_LR_15.fasta  -o rayd55tr

```

This worked quite well as I got a `scaffolds.fasta` file with two scaffolds one with 5730 bp and another with 1548 bp.

Now, I am confused how I can join those scaffolds to get one linear genome. I joined them manually and blasted just to check and I got something like below :.

Is there any tool to join these two contigs together to have one complete linear genome? I have reference genome too if needed.

# Answer
Ragout worked well for me. It is a reference-based scaffolder.

```
ragout.py --outdir output/ input/ragout.rcp > message.txt 2> error.txt --threads 8 --overwrite

```

The ragout.rcp file contains a 'recipe' with the reference sequence and the contigs in FASTA format.:

```
.references = ensembl
.target = assembly

ensembl.fasta = ref.fa
assembly.fasta = contigs.fa

```

In my case, it reduced the number of contigs 10-fold and increased the N50 10-fold. It added N's between the contigs.
> 2 votes

# Answer
In some cases you might consider reference-based scaffolding, as per charlesdarwin's answer. I would argue that in most cases this will only provide limited information, as this will only really provide you with the order of the contigs *as it is in your reference*, as opposed to how it is in your sampled organism. In that sense it really only provides pseudo-scaffolding that makes your assembly look artificially good.

In all cases you will have to check wether your raw data supports your scaffolding solution.

Example: you will never be able to call structural variants across contig boundaries if you use reference based scaffolding.

I would suggest you try different k-mer sized for your assembly algorithm, and look into dedicated scaffolding tools. Without knowing your type of sequencing data it is hard to suggest anything specific.
> 0 votes

---
Tags: ngs
---


# Question
Title: Transferring genomic features on new coordinates
I have a eukaryotic genome for which an updated sequence for a chromosome was recently obtained. I want to map RNAseq reads on the genome (and perform other downstream analyses) and would like to use the most up-to-date information possible (so the 'new' sequence of that chromosome).

However, before performing the read mapping, I'd like to update the annotation of my genome (GFF3) to match with the 'new' coordinates of this chromosome.

Basically, I would like to transfer information from the 'old-coordinates' GFF3 to obtain a 'new-coordinates' GFF3, and if possible conserve all the information/hierarchy of the file (gene, mRNA, exon, etc.).

**Example:**

```
chrXX   source  gene    222 5942    .   -   .   ID=gene_1;Name=gene_1;length=5720
chrXX   source  mRNA    222 5942    .   -   .   ID=gene_1.1;Parent=gene_1;Name=gene_1.1;length=5720
chrXX   source  exon    222 5794    .   -   .   ID=gene_1.1.2;Parent=gene_1.1
chrXX   source  exon    5889    5942    .   -   .   ID=gene_1.1.1;Parent=gene_1.1
chrXX   source  CDS 222 5794    .   -   1   ID=CDS:gene_1.1.2;Parent=gene_1.1;Name=gene_1.1
chrXX   source  CDS 5889    5942    .   -   0   ID=CDS:gene_1.1.1;Parent=gene_1.1;Name=gene_1.1

```

... should be updated to ...

```
chrXX   source  gene    333 6053    .   -   .   ID=gene_1;Name=gene_1;length=5720
chrXX   source  mRNA    333 6053    .   -   .   ID=gene_1.1;Parent=gene_1;Name=gene_1.1;length=5720
chrXX   source  exon    333 5905    .   -   .   ID=gene_1.1.2;Parent=gene_1.1
chrXX   source  exon    6000    6053    .   -   .   ID=gene_1.1.1;Parent=gene_1.1
chrXX   source  CDS 333 5905    .   -   1   ID=CDS:gene_1.1.2;Parent=gene_1.1;Name=gene_1.1
chrXX   source  CDS 6000    6053    .   -   0   ID=CDS:gene_1.1.1;Parent=gene_1.1;Name=gene_1.1

```

**The approach I tried** (but I am not sure it is the way to do it):

1. Extract sequence of all the features that are going to be updated to fasta (`bedtools getfasta`)
2. Map these sequences to the genome with the updated chromosome (`gmap` with `--nosplicing` since sequences we are mapping correspond to genomic regions).
3. Create the updated GFF3 file. When doing that, I apply one supplementary rule: if a feature mapped on another chromosome, but there was also another alignment on the updated chromosome, prioritize the one that correspond to the same (updated) chromosome.

What would be the proper way to do such a thing? Any suggestion of methods/tools is welcome!

# Answer
I think the standard way of doing this is to make a chain file, then use it to do a liftOver of the annotations:

```
mkdir psl
for i in ../ci3/rm/masked/*.masked; do blat ../ci2.2bit $i -tileSize=12 -fastMap -minIdentity=98 psl/`basename $i .fa.masked`.psl -noHead -minScore=100; done

```

Translate psl files to chains in the directory chain:

```
mkdir chain
for i in psl/*.psl; do axtChain -linearGap=medium -psl $i ../ci2.2bit ../ci3/ci3.2bit chain/`basename $i .psl`.chain; done

```

Merge short chains into longer ones into the directory chainMerge:

```
mkdir chainMerge
chainMergeSort chain/*.chain | chainSplit chainMerge stdin -lump=50

```

concat and sort the chains:

```
cat chainMerge/*.chain > all.chain
chainSort all.chain all.sorted.chain

```

Need info about chromosome sizes for netting:

```
twoBitInfo ../ci3/ci3.2bit ci3.chromInfo
twoBitInfo ../ci2.2bit ci2.chromInfo

```

Netting: identify alignable regions from chains:

```
mkdir net
chainNet all.sorted.chain ci2.chromInfo ci3.chromInfo net/all.net /dev/null

```

Finally, select the right alignable regions using the nets, creating a "liftOver" file:

```
netChainSubset net/all.net all.chain ci2ToCi3.liftOver

```

Run liftOver:

```
CrossMapy.py bed ci2ToCi3.liftOver test.hg18.bed

```
> 4 votes

# Answer
In addition to the answer by heathobrien you might take a look at RATT: Rapid Annotation Transfer Tool. The approach is essentially the same, only using more modern approaches to for instance whole genome alignments. It uses nucmer, which when using the latest version (4.0) allows you to do whole genome alignments using parallel computing.
> 2 votes

---
Tags: gff3, gmap
---