# Question
Title: How can longest isoforms (per gene) be extracted from a FASTA file?
Is there a convenient way to extract the longest isoforms from a transcriptome fasta file? I had found some scripts on biostars but none are functional and I'm having difficulty getting them to work.

I'm aware that the longest isoforms aren't necessarily 'the best' but it will suit my purposes.

The fasta was generated via Augustus. Here is what the fasta file looks like currently (sequence shortened to save space)

```
>Doug_NoIndex_L005_R1_001_contig_2.g7.t1
atggggcataacatagagactggtgaacgtgctgaaattctacttcaaagtctacctgattcgtatgatcaactcatca
ttaatataaccaaaaacctagaaattctagccttcgatgatgttgcagctgcggttcttgaagaagaaagtcggcgcaagaacaaagaagatagaccg
>Doug_NoIndex_L005_R1_001_contig_2.g7.t2
atggggcataacatagagactggtgaacgtgctgaaattctacttcaaagtctacctgattcgtatgatcaactcatca

```

The format is as such:

```
Gene 1 isoform 1
Gene 1 isoform 2
Gene 2 isoform 1
Gene 2 isoform 2

```

and so forth. There are several genes that have more than one pair of isoforms (up to 3 or 4). There are roughly 80,000 total transcripts, probably 25,000 genes. I would like to extract the single longest isoform for each gene.

# Answer
> 6 votes
While the solution from https://bioinformatics.stackexchange.com/users/96/daniel-standage should work (after adjusting for possible python3 incompatibility), the following is a shorter and less memory demanding method that uses biopython:

```
#!/usr/bin/env python
from Bio import SeqIO
import sys

lastGene = None
longest = (None, None)
for rec in SeqIO.parse(sys.argv[1], "fasta"):
    gene = ".".join(rec.id.split(".")[:-1])
    l = len(rec)
    if lastGene is not None:
        if gene == lastGene:
            if longest[0] < l:
                longest = (l, rec)
        else:
            lastGene = gene
            SeqIO.write(longest[1], sys.stdout, "fasta")
            longest = (l, rec)
    else:
        lastGene = gene
        longest = (l, rec)
SeqIO.write(longest[1], sys.stdout, "fasta")

```

If you saved this as `filter.py`, then `filter.py original.fa > subset.fa` would be the command to use.

# Answer
> 3 votes
Here is a solution in R. Could get really slow with big files. Works for the example you posted.

```
library(Biostrings)

## read your fasta in as Biostrings object
fasta.s <- readDNAStringSet("sample.fa")

## get the read names (in your case it has the isoform info)
names.fasta <- names(fasta.s)

## extract only the relevant gene and isoform id (split name by the period symbol)
gene.iso <- sapply(names.fasta,function(j) cbind(unlist(strsplit(j,'\\.'))[2:3]))

## convert to good data.frame = transpose result from previous step and add relevant column names
gene.iso.df <- data.frame(t(gene.iso))
colnames(gene.iso.df) <- c('gene','isoform')

## and length of isoforms
gene.iso.df$width <- width(fasta.s)

## split data.frame into list with entry for each gene
gene.iso.df.split <- split(gene.iso.df,gene.iso.df$gene)

## optional to keep all the information but really you just need indices
##gene.iso.df.split.best <- lapply(gene.iso.df.split,function(x) x[order(x$width)[1],])

## pull out the longest isoform ID for each gene (in case of a tie just take the first one)
best.id <- sapply(gene.iso.df.split,function(x) row.names(x)[order(x$width)[1]])

## subset your original reads with the subset
fasta.s.best <- fasta.s[best.id]

## export new fastafile containing longest isoform per gene
writeXStringSet(fasta.s, filepath='sample_best_isoform.fasta')

```

# Answer
> 3 votes
Late to the party here, but I like to try to avoid writing scripts when some command line magic will do. It's good practice to index your `FASTA` so use it.

```
samtools faidx <myfasta.fa>

```

A little bit of `awk` and `sort` can determine the largest isoform of each gene (and they don't even have to be sorted by name in the source `FASTA`).

```
awk -F'[\t.]' '{print $1,$2,$3,$4}' <myfasta.fa>.fai | sort -k4nr,4 | sort -uk1,2 | cut -f1-3 -d' '| tr ' ' '.' > selection.ls

```
* The `-F[\t.]` specifies both the `fai` tab delimiters, and the full stops in your contig names (`Doug_NoIndex_L005_R1_001_contig_2.g7.t1`) as delimiters for `awk`.
* For each line we print four columns. `$1` to `$3` are the three components of your contig name (after splitting on the `.`), and column $4 is the column in the index for sequence length.
* We use sort to organise the lines by `$4` (the contig length). Reverse sort with `-r` so the longest genes appear first.
* We then subsort by columns `$1` and `$2` (the contig and gene name (`g1`...)). `-u` then drops any repeat of column `$1` and `$2` pairs (that is, all lines after the first pair -- the longest length).
* Finally, we reassemble `$1` to `$3` with `cut` and `tr` to give you a list of sequence names.

A loop with our old friend `samtools faidx` can put your chosen sequences out to a new `FASTA`.

```
while read contig;
    do samtools faidx <myfasta.fa> $contig >> selection.fa;
done < selection.ls

```

Don't forget to index it, they are useful!

```
samtools faidx selection.fa

```

# Answer
> 1 votes
An ex-coworker (Josep Avril) has written a couple of very useful little scripts that convert fasta to tbl (`seqID<TAB>Sequence`) and back again. These are extremely handy for this sort of thing (and are included at the end of this answer). Using them, you can convert your fasta to a one sequence per line format, keep the longest sequence with a simple awk script and convert back to fasta.

I am assuming that the part after the last `.` in the ID line should be removed (so that `Doug_NoIndex_L005_R1_001_contig_2.g7.t1` and `Doug_NoIndex_L005_R1_001_contig_2.g7.t2` both map to `Doug_NoIndex_L005_R1_001_contig_2.g7`). If that is a correct assumption, this should work for you:

```
$ FastaToTbl file.fa | sed 's/\.[^.]*\t/\t/' |
    awk -F"\t" '{
                    if(length($2)>length(a[$1])){
                        a[$1]=$0
                    }
                }
                END{
                    for(i in a){
                        print a[i]
                    }
                }' | TblToFasta
>Doug_NoIndex_L005_R1_001_contig_2.g7
atggggcataacatagagactggtgaacgtgctgaaattctacttcaaagtctacctgat
tcgtatgatcaactcatcattaatataaccaaaaacctagaaattctagccttcgatgat
gttgcagctgcggttcttgaagaagaaagtcggcgcaagaacaaagaagatagaccg

```

One possible issue with that approach is that it needs to keep one sequence per ID in memory. If you have a huge fasta file and not much memory, that might be a problem. If that's the case, you can first sort the file to ensure that similar IDs always appear together and then print each line as soon as you find the next ID:

```
FastaToTbl file.fa | sed 's/\.[^.]*\t/\t/' | LC_ALL=C sort -t$'\t' -k1,1 |
    awk -F"\t" '{
                    if(NR>1 && prev!=$1){
                        print a[$1];
                        prev=$1
                    }
                    if(length($2)>length(a[$1])){
                        a[$1]=$0
                    }
                }
                END{
                    print a[$1]
                }' | TblToFasta

```

---

* FastaToTbl

```
#!/usr/bin/awk -f
{
        if (substr($1,1,1)==">")
        if (NR>1)
                    printf "\n%s\t", substr($0,2,length($0)-1)
        else
            printf "%s\t", substr($0,2,length($0)-1)
        else
                printf "%s", $0
}END{printf "\n"}

```
* TblToFasta

```
#! /usr/bin/awk -f
{
  sequence=$NF

  ls = length(sequence)
  is = 1
  fld  = 1
  while (fld < NF)
  {
     if (fld == 1){printf ">"}
     printf "%s " , $fld

     if (fld == NF-1)
      {
        printf "\n"
      }
      fld = fld+1
  }
  while (is <= ls)
  {
    printf "%s\n", substr(sequence,is,60)
    is=is+60
  }
}

```

# Answer
> 1 votes
This solution is will work for your example and probably all Augustus-derived Fastas, but mileage will vary beyond that.

```
#!/usr/bin/env python
from __future__ import print_function
import sys

def parse_fasta(data):
    """Stolen shamelessly from http://stackoverflow.com/a/7655072/459780."""
    name, seq = None, []
    for line in data:
        line = line.rstrip()
        if line.startswith('>'):
            if name:
                yield (name, ''.join(seq))
            name, seq = line, []
        else:
            seq.append(line)
    if name:
        yield (name, ''.join(seq))

isoforms = dict()
for defline, sequence in parse_fasta(sys.stdin):
    geneid = '.'.join(defline[1:].split('.')[:-1])
    if geneid in isoforms:
        otherdefline, othersequence = isoforms[geneid]
        if len(sequence) > len(othersequence):
            isoforms[geneid] = (defline, sequence)
    else:
        isoforms[geneid] = (defline, sequence)

for defline, sequence in isoforms.values():
    print(defline, sequence, sep='\n')

```

Save as `longest.py`, and invoke on the command line like so.

```
python longest.py < intput.fasta > output.fasta

```

---

The line `geneid = '.'.join(defline[1:].split('.')[:-1])` is key. Let me break it down.

* `defline[1:]`: ignore the first character of the defline (the `>` symbol)
* `defline[1:].split('.')`: split the string on `.` symbol
* `defline[1:].split('.')[:-1]`: ignore the last value after the split (the isoform name)
* `'.'.join(defline[1:].split('.')[:-1])`: join the split values again

This value will be the same for all isoforms for the same gene, so we use it in the script to keep track of the longest sequence associated with each gene ID.

---
Tags: fasta, isoform, filtering
---


# Question
Title: Why are there missing calls in a VCF file from exome sequencing?
My data is a VCF file generated from an exome sequencing variant call pipeline. I'm not very familiar with the sequencing and variant calling process. I noticed that there are some missing genotypes, which are recorded as "./." at the GT field. From googling I learned that they're not homozygous reference genotype ("0/0"), but missing calls due to some sequencing failure.

In order to make my data "cleaner", I thought it would be better to filter out the loci with missing calls if there're not too many of them. I also checked the corresponding "DP" of the loci with and without missing calls. For example:

```
chr1:123 GT:DP 0/0:2 1/1:4 ./.
chr1:234 GT:DP 0/0:10 1/1:11 1/1:20
chr1:345 GT:DP 0/1:40 1/1:37 0/0:78
chr1:456 GT:DP 0/1:7 0/0:23 ./.
chr1:567 GT:DP 0/1:34 1/1:39 0/0:58

```

In the above toy example, there are 3 people and 5 loci. I checked the mean DP of all 5 loci and I found that the mean DP of loci with missing calls (1st and 4th loci) to be significantly lower than the mean DP of loci without missing calls (2nd, 3rd and 5th loci). This happens to my real data with a lot more loci and samples than this toy example. Is that a coincidence? Or are there any specific reason that the loci with missing calls have lower coverage than the "normal" loci? Thanks!

# Answer
> 4 votes
Missing variant calls due to lack of coverage shouldn't happen in the targeted capture region and I'd think most of these would come from off-target regions where some samples had reads mapped. I'd filter out the VCF to only include on-target loci before proceeding with further analyses.

If you did your variant calling on samples separately and then merged them together, you'll have missing calls at loci called variant in one sample and homozygous reference in the other.

You can get over this by asking the variant caller to emit all sites within the target region, which would lead to a large VCF or do multi-sample calling, which is more computationally intensive.

One other way is to do variant calling on all samples separately, generating a BED file of all variant regions in the cohort and then asking the variant caller to genotype those loci across all samples before merging them together.

# Answer
> 3 votes
I've just been generating data like this, so can tell you about why/how missing calls are created in my dataset. There are two main reasons:

## 1. Sequencing failure

When reads don't map across the variant region, then it's impossible to accurately determine a genotype for that region. This will commonly happen just outside the borders of the selected regions for exome sequencing, but can also happen through natural random sampling of the genome, or through repeated sequences coupled with systematic error in either the assembled genome or the sequencer. In this case, it would be expected that called variant regions including missing data would be lower coverage. It may be helpful to plot the coverage within a particular region to work out if this is likely to be happening.

## 2. Merged non-variant data

Variant calling can be processed in two different modes, one which takes into account *all* samples when doing calling, and another which does the calling one sample at a time, followed by subsequent merging of the results. Multi-sample variant calling is more accurate, but computational limits (or experimental design) might mean that it is more appropriate for samples to be called one-by-one. It's fairly common (for the purpose of preserving space and computational power) for calling algorithms to only report variant information so if a single sample is identical to the reference, and that sample is the only one being called, then the variant and coverage information for that sample will be lost. When samples are merged after being separately called, the declared call for the "same as reference" samples will be set to missing. In this case, it would be generally expected that the other called samples would have high-coverage within the variant region.

---
Tags: variant-calling, vcf, ngs, exome, reads
---


# Question
Title: Get the mapping statistics of a single read from a BAM file
I have a BAM file, and I have a read ID. What is the simplest way to get mapping statistics of that read in human-readable format?

E.g. I might want: % identity of aligned bases; number of insertions and deletions; number of bases aligned to the reference; etc.

Extra points for one-liners that give detailed summary statistics, and extra-extra points if you can extend the answer to output summaries for all reads longer than some predetermined length.

# Answer
> 7 votes
`samtools stats` seems to be able to do most of this, excluding the CIGAR-string parsing stuff (i.e. INDELs):

```
$ samtools view -h mapped.bam | grep -e '^@' -e 'readName' |
    samtools stats | grep '^SN' | cut -f 2-

raw total sequences:    2
filtered sequences: 0
sequences:  2
is sorted:  1
1st fragments:  2
last fragments: 0
reads mapped:   2
reads mapped and paired:    0   # paired-end technology bit set + both mates mapped
reads unmapped: 0
reads properly paired:  0   # proper-pair bit set
reads paired:   0   # paired-end technology bit set
reads duplicated:   0   # PCR or optical duplicate bit set
reads MQ0:  0   # mapped and MQ=0
reads QC failed:    0
non-primary alignments: 0
total length:   13632   # ignores clipping
bases mapped:   13632   # ignores clipping
bases mapped (cigar):   13502   # more accurate
bases trimmed:  0
bases duplicated:   0
mismatches: 3670    # from NM fields
error rate: 2.718116e-01    # mismatches / bases mapped (cigar)
average length: 6816
maximum length: 6816
average quality:    8.9
insert size average:    0.0
insert size standard deviation: 0.0
inward oriented pairs:  0
outward oriented pairs: 0
pairs with other orientation:   0
pairs on different chromosomes: 0

```

For a length filter, an awk length check can be slotted in instead of the grep:

```
$ samtools view -h mapped.bam |
    awk -F'\t' '{if((/^@/) || (length($10)>500)){print $0}}' |
    samtools stats | grep '^SN' | cut -f 2-

raw total sequences:    131
filtered sequences: 0
sequences:  131
is sorted:  1
1st fragments:  131
last fragments: 0
reads mapped:   131
reads mapped and paired:    0   # paired-end technology bit set + both mates mapped
reads unmapped: 0
reads properly paired:  0   # proper-pair bit set
reads paired:   0   # paired-end technology bit set
reads duplicated:   0   # PCR or optical duplicate bit set
reads MQ0:  0   # mapped and MQ=0
reads QC failed:    0
non-primary alignments: 0
total length:   569425  # ignores clipping
bases mapped:   569425  # ignores clipping
bases mapped (cigar):   453884  # more accurate
bases trimmed:  0
bases duplicated:   0
mismatches: 113755  # from NM fields
error rate: 2.506257e-01    # mismatches / bases mapped (cigar)
average length: 4346
maximum length: 16909
average quality:    9.1
insert size average:    0.0
insert size standard deviation: 0.0
inward oriented pairs:  0
outward oriented pairs: 0
pairs with other orientation:   0
pairs on different chromosomes: 0

```

# Answer
> 2 votes
Use htsbox:

```
htsbox samview -p in.bam | less -S
htsbox samview -pS in.sam | less -S

```

It outputs mapping positions in the PAF format, which looks something like:

```
read1  4983  774  4982  +  chr18  80373285  26911072  26915544  3835  4631  60 \
       mm:i:214  io:i:119  in:i:159  do:i:339  dn:i:423

```

This is one long line in terminal. I folded the line for display purposes. For the first 12 fixed fields, see the table in PAF page. They give mapping positions and identity. The optional fields tell you #mismatches (mm), #insOpens (io), #insertions (in), #delOpens (do) and #deletions (dn). Note that the alignment MUST have an "NM" tag; otherwise the number of matching bases (col. 11) is overestimated and "mm" will always be zero.

You can easily compute summary statistics from PAF if you are familiar with command lines. For example:

```
htsbox samview -p in.bam | awk '{x+=$10;y+=$11}END{print x/y}'  # identity
htsbox samview -p in.bam | awk '$2>1000{x+=$10;y+=$11}END{print x/y}'  # identity for >1000bp
htsbox samview -p in.bam \
  | perl -ane '{$g+=$1 while /[id]n:i:(\d+)/g;$y+=$F[10]}END{print $g/$y,"\n"}' # gap rate

```

If you don't like long one-liners, feel free to write a proper script to collect all statistics in one go.

---
Tags: sam
---


title: Count genomic ranges

text:
# Question
Title: Count genomic ranges
I have a set of genomic ranges that are potentially overlapping. I want to count the amount of ranges at certain positions using R.

I'm Pretty sure there are good solutions, but I seem to be unable to find them.

Solutions like cut or findIntervals don't achieve what I want as they only count on one vector or accumulate by `all values <= break`.

Also `countMatches {GenomicRanges}` doesn't seem to cover it.

Probably one could use Bedtools, but I don't want to leave R.

I could only come up with a hilariously slow solution

```
# generate test data
testdata <- data.frame(chrom = rep(seq(1,10),10),
                       starts = abs(rnorm(100, mean = 1, sd = 1)) * 1000,
                       ends = abs(rnorm(100, mean = 2, sd = 1)) * 2000)

# make sure that all end coordinates are bigger than start
# this is a requirement of the original data
testdata <- testdata[testdata$ends - testdata$starts > 0,]

# count overlapping ranges on certain positions
count.data <- lapply(unique(testdata$chrom), function(chromosome){
    tmp.inner <- lapply(seq(1,10000, by = 120), function(i){
        sum(testdata$chrom == chromosome & testdata$starts <= i & testdata$ends >= i)
    })
    return(unlist(tmp.inner))
})

# generate a data.frame containing all data
df.count.data <- ldply(count.data, rbind)

# ideally the chromosome will be columns and not rows
t(df.count.data)

```

# Answer
> 6 votes
`GenomicRanges::countOverlaps` seems to be what you’re after:

```
position_range = GRanges(position$chrom, IRanges(position, position, width = 1))
ranges_at_position = countOverlaps(position_ranges, granges)

```

---
Tags: r
---


# Question
Title: For what bioinformatics tasks is Biopython more adapted than Bioperl?
Are there any advantages to learning Biopython instead of learning Bioperl?

Ideally, we would learn both, but someone starting out in bioinformatics may have to choose what to learn first depending on the kind of problems actually encountered.

**Are there problems for which Biopython is better than Bioperl (or vice-versa)?**

# Answer
> 5 votes
Regading the perl vs python discussion, there is no final answer which language is better, but I have some advice for you:

**Learn the language your colleagues or your advisor use**. This way you are able to discuss your code with them and also get help if you run into problems.

# Answer
> 5 votes
Currently you could use either but a major question is which platform will others be using in the future. AFAIK Perl is only superior to Python for regex.

Based on the trend I see for new programmers and new software being released: Perl is on the way out and Python is still growing.

https://trends.google.com/trends/explore?date=all&q=bioperl,biopython

# Answer
> 4 votes
This usually comes down to religious issues, so let me try and steer it back to more objective grounds:

* What language do you know (better)? Use the library for that one.
* If you know neither and will be learning a language to use the library, the majority opinion would be that Python is easier to learn.
* However, some people say that they "click" with Perl better
* Python is what the majority of bioinformaticians are using at the moment (about 60%). And there is a virtue and aid in using what everyone else is.
* What are your colleagues / collaborators going to be using?
* Broadly, newer libraries tend to be written in Python.
* Conversely, many sing the praises of regular expressions and string handling in Perl, citing speed and ease of use.
* I think the numerics support in Python may be markedly better than that in Perl.

# Answer
> 1 votes
I agree with @story and @agapow here. As it stands right now BioPerl is still much more mature than BioPython. If there's a tool or method that only exists in the BioPerl library then you should probably use BioPerl to accomplish what you need to accomplish. But it's undeniable that BioPython is progressing very quickly and will soon overtake BioPerl. That being said there are still a lot of plus side regarding Perl in the Bioinformatics realm and it will definitely help to be somewhat competent in it.

So in short my answer is to use both! I find myself mostly using Python and BioPython but if I need to accomplish something using a tool that isn't available in the BioPython library or can be accomplished more quickly using Perl I will look towards BioPerl/Perl.

---
Tags: python, perl
---


# Question
Title: How can we distinguish between true zero and dropout-zero counts in single-cell RNA-seq?
In single-cell RNA-seq data we have an inflated number of 0 (or near-zero) counts due to low mRNA capture rate and other inefficiencies.

How can we decide which genes are 0 due to gene dropout (lack of measurement sensitivity), and which are genuinely not expressed in the cell?

Deeper sequencing does not solve this problem as shown on the below saturation curve of 10x Chromium data:

Also see Hicks et al. (2017) for a discussion of the problem:

> Zero can arise in two ways:
>
> 1. the gene was not expressing any RNA (referred to as *structural zeros*) or
> 2. the RNA in the cell was not detected due to limitations of current experimental protocols (referred to as *dropouts*)

# Answer
> 15 votes
Actually this is one of the main problems you have when analyzing scRNA-seq data, and there is no established method for dealing with this. Different (dedicated) algorithms deal with it in different ways, but mostly you rely on how good the error modelling of your software is (a great read is the review by Wagner, Regev & Yosef, esp. the section on "False negatives and overamplification"). There are a couple of options:

* You can **impute values**, i.e. fill in the gaps on technical zeros. CIDR and scImpute do it directly. MAGIC and ZIFA project cells into a lower-dimensional space and use their similarity there to decide how to fill in the blanks.
* Some people straight up **exclude genes** that are expressed in very low numbers. I can't give you citations off the top of my head, but many trajectory inference algorithms like monocle2 and SLICER have heuristics to choose informative genes for their analysis.
* If the method you use for analysis doesn't model gene expression explicitly but **uses some other distance method** to quantify similarity between cells (like cosine distance, euclidean distance, correlation), then the noise introduced by dropout can be covered by the signal of genes that are highly expressed. Note that this is dangerous, as genes that are highly expressed are not necessarily informative.
* **ERCC spike** ins can help you reduce technical noise, but I am not familiar with the Chromium protocol so maybe it doesn't apply there (?)

since we are speaking about noise, you might consider using a protocol with unique molecular identifiers. They remove the amplification errors almost completely, at least for the transcripts that you capture...

EDIT: Also, I would highly recommend using something more advanced than PCA to do the analysis. Software like the above-mentioned Monocle or destiny is easy to operate and increases the power of your analysis considerably.

# Answer
> 1 votes
Some people use imputation to differentiate between true zeros and dropout in single-cell data. Some approaches you can look into:

---
Tags: rna-seq, ngs
---


# Question
Title: How should the SAM MD tag match the CIGAR string?
I am trying to understand how the MD:Z tag is used. The following is from the SAM Optional Fields Specification, which gives an example but is not thorough.

> The MD field aims to achieve SNP/indel calling without looking at the reference. For example, a string ‘10A5^AC6’ means from the leftmost reference base in the alignment, there are 10 matches followed by an A on the reference which is different from the aligned read base; the next 5 reference bases are matches followed by a 2bp deletion from the reference; the deleted sequence is AC; the last 6 bases are matches. The MD field ought to match the CIGAR string.

Suppose I have a read that I want to soft clip at both ends. If the read starts with CIGAR `100M` and MD `50G49`, and I want to change the CIGAR to `7S86M7S`, what should the MD field become?

Is there a more complete explanation of what appears in the MD tag field?

# Answer
> 5 votes
The MD string doesn't apply to soft or hard clipped regions,. so your example read becomes `43G42`. Since variant calling using something simplistic like this is only ever going to use the aligned portion, I guess there was never much reason to bother with adding in soft-clipping information. This also means that you need to parse the CIGAR string if you're using the MD string.

---
Tags: file-formats, sam, cigar
---


title: Customizing bigWig file

text:
# Question
Title: Customizing bigWig file
I generate bigWig files using `bamCoverage` from deeptools, in part for my colleagues to visualize their mapped libraries in the IGV viewer.

A problem is that the displayed track name is apparently the file name, which is not convenient because some files for different libraries have the same name. The reason is that the files are organized with one directory per library. I would prefer to avoid changing the way I name and organize the files.

**Is there a way to edit my bigWig files in order to have control of the track name independently from the file name?**

As far as I understand, bigWig format derives from wiggle format, which have a track definition line that can contain a `name` and a `definition` attribute: https://genome.ucsc.edu/goldenpath/help/customTrack.html#TRACK.

Maybe there is a way to convert a bigWig file back to wiggle, edit the track definition, and then rebuild an updated bigWig file.

# Answer
> 4 votes
There's no equivalent to the wiggle header in bigWig (or bigBed) files, which is why UCSC uses the file name. This is actually the reason for the track line stuff that you linked to, since you can then specify a name and just point to where the bigWig (or other format) file is on the internet.

BTW, you can certainly convert your bigWig to wiggle, add the track line, and convert it back with UCSC tools (e.g., `bigWigToWig` and `wigToBigWig`). It shouldn't change the track name when you display it in UCSC (though if it does, then let me know, since I'll modify pyBigWig and libBigWig to support that then).

---
Tags: file-formats, format-conversion, bigwig, visualization, genome-browser
---


# Question
Title: Why do BLASTn and prokka not seem to be searching the whole fasta file?
When I use blastn and prokka (I will detail exactly how I did so below) on a 2.8 million bp fasta file I get output start/end numbers that do not seem to cover the entire genome.

Starting with a .fna genome such as genome.fna I ...

1 - blastn

Searched the genome for 3135 different 28-mers using BLASTn.

```
makeblastdb -in genome.fna -dbtype nucl -parse_seqids -out ./output/genome

```

*next command in python*

```
blast_tsv = NcbiblastnCommandline(query=Q, db=DB, perc_identity=100, outfmt=6, out=(OUT))
    stdout, stderr = blast_tsv()

```

*Q is this list of k-mers. DB is the database created*

This outputs a list of search results but the important thing is that none of the start/end (columns 7 and 8) integers were greater than 100,000 yet the entire genome is 2.8 million base pairs long. I can provide the link to this file in the comments. Does this have to do with blastn stopping searching after finding one match? And if so how can I tell blastn to search for every match for each k-mer? (I'm open to using other blastn programs other than biopython)

2 - prokka

I used prokka to create gff files for the genome.

```
prokka-1.12/bin/prokka —setupdb
prokka-1.12/bin/prokka -kingdom Bacteria -rfam -outdir ./prokka_database/genome -force -prefix "${genome/.fna/}" genome.fna

```

Th gff file produced by this command only described genes with start and end positions less than 200k. I can provide the link to this file as well.

# Answer
With a k-mer size of 28 it shouldn't be finding that many matches. And the prokka results are suspicious as well. Maybe you have multiple contigs (none larger than 100kb) in that file? What is the result of

```
grep ^'>' fasta_file | wc -l

```

? This would show how many contigs you have in the file.
> 6 votes

---
Tags: sequence-annotation, genome, k-mer, python, blast
---


# Question
Title: What methods are available to find a cutoff value for non-expressed genes in RNA-seq?
I have a gene expression count matrix produced from bulk RNA-seq data. I'd like to find genes that were **not** expressed in a group of samples and were expressed in another group.

The problem of course is that not all *effectively* non-expressed genes will have 0 counts due to sequencing errors, or because they were expressed in a small subset of cells.

I'm interested in solutions using `R`.

# Answer
> 13 votes
> I'd like to find genes that were not expressed in a group of samples and were expressed in another group.

This is, fundamentally, a differential expression analysis, with a twist. To solve this, you’d first use a differential expression library of your choice (e.g. DESeq2) and perform a one-tailed test of differential expression.

Briefly, you’d perform the normal setup and then use

```
results(dds, altHypothesis = 'greater')

```

To perform a one-tailed test. This will give you only those genes that are significantly upregulated in one group. Check chapter 3.9 of the vignette for details.

Of course this won’t tell you that the genes are *unexpressed* in the other group. Unfortunately I don’t know of a good value to threshold the results; I would start by plotting a histogram of the (variance stabilised) expression values in your first group, and then visually choose an expression threshold that cleanly separates genes that are clearly expressed from zeros:

```
vst_counts = assay(vst(dds))
dens = density(vst_counts[, replicate])
plot(dens, log = 'y')

```

(This merges the replicates in the group, which should be fine.)

Counts follow a multimodal distribution, with one mode for unexpressed and one or more for expressed genes. The expression threshold can be set somewhere between the clearly unexpressed and expressed peaks:

Here I used `identify(dens)` to identify the threshold interactively but you could also use an analytical method:

```
threshold = identify(dens)
quantile = sum(dens$x < dens$x[threshold]) / length(dens$x)

# Using just one replicate here; more robust would be to use a mean value.
nonzero_counts = counts(dds, normalized = TRUE)[, replicates[1]]
nonzero_counts = nonzero_counts[nonzero_counts > 0]

(expression_threshold = quantile(nonzero_counts, probs = quantile))

```

```
26.5625%
4.112033

```

# Answer
> 10 votes
A common method is to use zFPKMs, which you can find implemented in R here.

Having said that, there's an inherent problem in declaring a difference between two things on either side of a given threshold. Given that, I would encourage you to use more than "expressed in one and not in another", likely adding at least a "with a minimal difference of X between them" metric. You may also find the tau metric useful, which is implemented in R here. This is meant to measure tissue-specificity, which is more akin to what you're probably interested in doing.

# Answer
> 5 votes
Depending on how much effort you wish to put into this, here is one suggestion I have see used before (uses more than just R). Steps 1-4 come from this paper (see supplimentary material section "Calculation of per gene Local FDR").

1. *Take your gene models and collapse the introns down to maybe 100bp*
Thus if we have a gene with two exons and a 1kb intron eg. exons (1000, 1100) and (2000,2100), we reduce the intron size so that the exons are (1000,1100) and (1200,1300). This is because we need to find sufficient gene free space to fit the null gene in.
2. *Shift each gene into the nearest genomic space at least 5 kp from a genomic annotation and free of ESTs.*
3. *Calculate the FPKM distribution of this shifted null set*
4. *From these two distributions it should be possible to calculate a local FDR for expression.*
You will need to use a "mixing proportion" to do this. The above reference used Qvality to do this.
5. *Carry out differential expression analysis as described in Konrad's answer*
6. *Subset the differentially expressed genes to only consider those that are not expressed in the control condition*
7. *You might want to also recalculate the FDRs for differential expression as you are only considering a subset of the tests you might have.*

# Answer
> 0 votes
Unfortunatly that is not doable - there is now way of distinguishing between features not expressed and features not expressed for technical reasons. Furthermore a lot of genes might seem lowly expressed to to technical noise.

This also means that you will have to choose a cutoff and simply use that. For choosing a cutoff looking at the replicate distribution of log10(counts) and log10(FPKM) are probably the best options.

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


# Question
Title: Spliced vs. unspliced ratios for transcripts in RNA-seq data
Is there a computational tool for measuring what percentage of RNA is spliced in an RNAseq experiment?

I'm not particularly interested in complicated analyses that give ratios for all possible alternative splicing variations. I'd rather have a binary classification. For example: (RNA1 = 35% spliced, 65% unspliced)

**Edit: (Additional information)**

We're very early in the exploratory stage, so a specific hypothesis is still forming. But we're interested in the sequence content of lncRNAs, many of which are poorly spliced. So we want to get a sense for how prevalent unspliced version of these transcript are.

Because this is currently exploratory, estimations and noise are fairly acceptable. A very basic approach might be to calculate the ratio of reads that map to exon-to-exon junctions, verses those that map from exon-to-intro junctions.

A more advanced approach I'm also considering is Salmon.

# Answer
I don't know about tools, but i've used the following python code to calculate the ratio of reads that overlap the 5' or 3' ends of introns or that are spliced. We sum these across all introns in a gene set (we actaully use this for iCLIP analysis to see if RNA binding proteins bind pre-mRNA or spliced RNA).

```
import pysam
from collections import Counter
from CGAT import GTF, IOTools

def calculateSplicingIndex(bamfile, gtffile, outfile):

    bamfile = pysam.AlignmentFile(bamfile)

    counts = Counter()

    for transcript in GTF.transcript_iterator(
            GTF.iterator(IOTools.openFile(gtffile))):

        introns = GTF.toIntronIntervals(transcript)

        for intron in introns:
            reads = bamfile.fetch(
                reference=transcript[0].contig,
                start=intron[0], end=intron[1])

            for read in reads:
                if 'N' in read.cigarstring:
                    blocks = read.get_blocks()
                    starts, ends = zip(*blocks)
                    if intron[0] in ends and intron[1] in starts:
                        counts["Exon_Exon"] += 1
                    else:
                        counts["spliced_uncounted"] += 1
                elif (read.reference_start <= intron[0] - 3
                      and read.reference_end >= intron[0] + 3):
                    if transcript[0].strand == "+":
                        counts["Exon_Intron"] += 1
                    else:
                        counts["Intron_Exon"] += 1
                elif (read.reference_start <= intron[1] - 3
                      and read.reference_end >= intron[1] + 3):
                    if transcript[0].strand == "+":
                        counts["Intron_Exon"] += 1
                    else:
                        counts["Exon_Intron"] += 1
                else:
                    counts["unspliced_uncounted"] += 1

    header = ["Exon_Exon",
              "Exon_Intron",
              "Intron_Exon",
              "spliced_uncounted",
              "unspliced_uncounted"]

    with IOTools.openFile(outfile, "w") as outf:

        outf.write("\t".join(header)+"\n")
        outf.write("\t".join(map(str, [counts[col] for col in header]))
                   + "\n")

```

Unfortunately this uses a bunch of libraries you may or may not have, including CGAT (for the GTF parser and IOTools package) and pysam.

Once you've got these statistics you can calculate the "splicing index" as the log2 ratio of 2 times the number of spliced reads divided by the number of reads overlapping the 3' and 5' ends of introns.
> 3 votes

# Answer
Take this first comment with a grain of salt, since this isn't an area I've worked in much, but: is binary classification possible? If a gene has 3 introns, and 2 are spliced out but 1 is retained, is this "spliced" or "unspliced". My first impression is that an analysis would be a bit more nuanced than a binary classification.

That said, I'm aware of a tool called Keep Me Around for intron retention analysis (preprint, code). I've never used this software, but it looks like it's actively maintained and was created by a research group with a lot of experience and influence in this area, for whatever that's worth. :-)

**UPDATE**: It looks like the original post was modified since I submitted my answer, but I'll leave the text about binary classification around for posterity. :-)
> 2 votes

---
Tags: rna-seq, rna-splicing
---


# Question
Title: What motif finding software is available for multiple sequences ~10Kb?
I have around ~3,000 short sequences of approximately ~10Kb long. What are the best ways to find the motifs among all of these sequences? Is there a certain software/method recommended?

There are several ways to do this. My goal would be to:

(1) Check for motifs repeated within individual sequences

(2) Check for motifs shared among all sequences

(3) Check for the presence of "expected" or known motifs

With respect to #3, I'm also curious if I find e.g. trinucleotide sequences, how does one check the context around these regions?

Thank you for the recommendations/help!

# Answer
For (3), this page has a lot of links to pattern/motif finding tools. Following through the YMF link on that page, I came across the University of Washington Motif Discovery section. Of these projection seemed to be the only downloadable tool. I find it interesting how old all these tools are; maybe the introduction of microarrays and NGS has made them all redundant.

Your sub-problem (2) seems similar to the problem I'm having with *Nippostrongylus brasiliensis* genome sequences, where I'd like to find regions of very high homology (length 500bp to 20kb or more, 95-99% similar) that are repeated throughout the genome. These sequences are killing the assembly.

The main way I can find these regions is by looking at a coverage plot of long nanopore reads mapped to the assembled genome (using GraphMap or BWA). Any regions with substantially higher than median coverage are likely to be shared repeats.

I've played around in the past with chopping up the reads to smaller sizes, which works better for hitting smaller repeated regions that are such a small proportion of most reads that they are never mapped to all the repeated locations. I wrote my own script a while back to chop up reads (for a different purpose), which produces a FASTA/FASTQ file where all reads are exactly the same length. For some unknown reason I took the time to document that script "properly" using POD, so here's a short summary:

> Converts all sequences in the input FASTA file to the same length. Sequences shorter than the target length are dropped, and sequences longer than the target length are split into overlapping subsequences covering the entire range. This prepares the sequences for use in an overlap-consensus assembler requiring constant-length sequences (such as edena).

And here's the syntax:

```
$ ./normalise_seqlengths.pl -h
Usage:
    ./normalise_seqlengths.pl <reads.fa> [options]

  Options:
    -help
      Only display this help message

    -fraglength
      Target fragment length (in base-pairs, default 2000)

    -overlap
      Minimum overlap length (in base-pairs, default 200)

    -short
      Keep short sequences (shorter than fraglength)

```
> 1 votes

# Answer
The MEME Suite web site contains a collection of tools for motif analysis (I'm one of the maintainers). It contains two de novo motif discovery tools: MEME and DREME. Public web applications are provided, but you can also download and build command line tools for a local installation.

For your first goal you could use MEME and select the "Any number of repetitions model" (ANR). For your second goal, you'd use MEME with the "Zero or One Occurrences Per Sequence" (ZOOPS) model. For your third goal you could use FIMO (Find Individual Motif Occurrences), and one or more of the motif databases provided on the software and database download page.

It sounds like your sequence data is about 30Mb. The MEME web application is limited to 60kb of sequence data, so you'd have to install a local copy of the MEME Suite. MEME would take a long time to analyze a 30Mb sequence database unless you have MPI configured, and lots of cores available. You might want to consider analyzing a randomly selected subset of your sequences. The running time of MEME grows as the cube of the number of sequences.

For short motifs, you may want to use DREME rather than MEME. DREME is better than MEME at identifying short motifs, but is limited to motifs \<= 8 positions wide.
> 3 votes

# Answer
Check out HOMER. "Software for motif discovery and next generation sequencing analysis", it's what my lab uses currently for finding eRNA motifs.

Edit: For @ShanZhengYang "HOMER was designed as a de novo motif discovery algorithm..." HOMER De Novo Motif
> 2 votes

# Answer
Most tools I know of looks for enrichment of specific motifs - but that requires that you have a set of sequences which are of special interest and a background set to test against.

Is that your case?

Update after comments 12th June 2017.

You could try the meme suite more specifically the motif finder
> 1 votes

# Answer
it is under development, but maybe BaMMmotif! is something for you? Its main selling point is that it can look for motifs enriched in a set of sequences of equal length *de novo*. If you can't/don't want to supply a negative set it learns one from the positive sequences. There is a wealth of options to choose from if you have more information about your sequences: there are different models for "zero or one", "one" and "multiple" occurrences of the motif.

You can also use it to look for known motifs, if you encode them as an XXmotif PWM. If you have a file with motifs (like binding sites) you can use that as initialization as well.

While I have not used the software myself, the authors are very responsive on git and the installation instructions seem pretty straightforward.

EDIT: Apparently the software was developed with ChIP experiments as the usual use case scenario, it might misbehave for larger sequences or take super long to run.
> 1 votes

---
Tags: motifs, software-recommendation
---


# Question
Title: Time and memory efficient processing of many intervals for a bigWig file
I am trying to calculate the mappability adjusted length of introns as described by Boutz et al.

Briefly, for each intron I wish to calculate the length minus the number of bases that are non-uniquely mappable. Mappability tracks can be downloaded as bigWigs from here. The the score at each base is 1/number of mapping positions, so 1 indicates a read can be uniquely mapped at this location.

I am struggling to do this computation in either a reasonable amount of time, or a reasonable amount of memory.

First I tried importing the whole bigwig:

```
library(rtracklayer)
mappability <- import(BigWigFile("wgEncodeCrgMapabilityAlign50mer.bigWig"),
                      as="NumericList",
                      selection=BigWigSelection(intron_ranges))

```

where `intron_ranges` is a `GRanges` object with the introns (about 800,000 of them).

This used **way** too much memory, and soon caused by machine to fall over.

Second I tried processing one intron at a time:

```
mappability_file = BigWigFile("wgEncodeCrgMapabilityAlign50mer.bigWig")
effective_length <- function (gr) {
  intron_selection <- BigWigSelection(gr)
  scores <- import(mappability_file, as="NumericList", selection=intron_selection)
  non_unique <- sum(scores[[1]] < 1.0)
  eff_len = width(gr)[1] - non_unique
  return(eff_len)
}

mappability <- sapply(intron_ranges, effective_length)

```

This has been running for hours and shows no sign of finishing.

Is there a way to do this that uses less than, say 4GB of RAM, but finishes in say less than 30 minutes? It feels like there should be.

I'm not tied to R, just what I've tried so far. Happy with answers using binary packages, common shell tools, python or R.

# Answer
Assuming you can make a BED file of introns, you can then use the `pyBigWig` module in python:

```
import pyBigWig
import numpy

bw = pyBigWig.open("some file.bw")
bed = open("introns.bed")
for line in bed:
    cols = line.strip().split("\t")
    vals = bw.values(cols[0], int(cols[1]), int(cols[2]), numpy=True)
    effLen = cols[2] - cols[1] - (vals < 1.0).sum()
    # Do something with this.

```

One can make that parallel (with the deeptools API) should it prove too slow. This will only use an appreciable amount of memory for very large introns, since there you need to store the value at each base. One could make this much more memory efficient by using `bw.intervals()` instead, but that'd require quite a bit more code (you'd need to determine the number of bases overlap per interval and then assume that either each base you want to query has an overlapping interval or, if not, keep track of that fact).
> 2 votes

---
Tags: bigwig
---


# Question
Title: validating identified sub-populations of cells in scRNA-seq
In the analyses of single-cell RNA-seq data there are different unsupervised approaches to identify putative subpopulations (e.g. as available with Suerat or SCDE packages).

Is there a good way of computationally validating the cluster solutions? Different methods may results in slightly different clustering results. How to know which one is the best i.e. representative of biological sub-populations?

# Answer
> 6 votes
A SC3, single-cell consensus clustering, approach could be helpful here. It aims at achieving "high accuracy and robustness by combining multiple clustering solutions through a consensus approach" https://www.nature.com/nmeth/journal/v14/n5/full/nmeth.4236.html

# Answer
> 3 votes
While better methods of evaluating your clusters would be to use an external dataset or a dataset with known truth, there are a variety of internal validation metrics that can be used to compare clustering solutions without another dataset.

Here are a few metrics:

* Davies-Bouldin Index
* Calinski-Harabasz Index
* Root-Mean-Square Standard Deviation

Many more can be found in this clustering review: http://stke.sciencemag.org/content/9/432/re6

These internal validation metrics grade your clustering solution based on three measures: compactness, connectedness, and separation. When using these metrics to compare clustering solutions, be sure to consider which metric is appropriate for your results as some algorithms work by optimizing certain measures.

---
Tags: scrnaseq, clustering
---


# Question
Title: Is the optional SAM NM field strictly computable from the MD and CIGAR?
From SAM Optional Fields Specification the NM field is

> Edit distance to the reference, including ambiguous bases but excluding clipping

Assuming both the MD and CIGAR are present, is the edit distance simply the number of characters `[A-Z]` appearing in the MD field plus the number of bases inserted (`xI`, if any) from the CIGAR string? Are there any other complications?

# Answer
> 4 votes
Assuming both the MD and CIGAR are present **and correct**, then yes, you can parse both to get the edit distance (`NM` auxiliary tag). One big caveat to this is that there's a reason that the `samtools calmd` command exists, since it's historically been the case that not all aligners have output correct MD strings. It's rare for the CIGAR string to be wrong and that'd be more of a catastrophic error on the part of an aligner. For what it's worth, if the `NM` auxiliary is absent on a given alignment but present on others produced by the same aligner then it's fair to assume `NM:i:0` for a given alignment by default (many aligners only produce `NM:i:XXX` if the edit distance is at least 1).

---
Tags: file-formats, sam, cigar
---


# Question
Title: Given a transcription factor, what genes does it regulate?
I have a list of transcription factors and I am interested in finding out which which genes might be transcribed as a result of the formation of transcription factor complexes with that transcription factor.

Any ideas on good databases? I've seen that ENCODE hosts data from CHIPseq experiments, but I'm not sure if I can find conclusions regarding interactions from this site. I'm looking for a database with known or putative transcription/gene interactions.

Thanks in advance.

# Answer
MSigDB has a collection (C3:TFT) of gene sets corresponding to transcription factor targets.

Harmonizome has functional terms for genes extracted from over a hundred publicly available resources.
> 5 votes

# Answer
iRegulon takes a sequence-based approach to finding transcription factor targets. There's a Cytoscape app that you can use to find the regulators of a given gene list, or the targets of a particular transcription factor.

Transcription factor binding sites are predicted using a collection of position weight matrices (PWMs) from a number of sources, including Jaspar and TRANSFAC.

iRegulon: From a Gene List to a Gene Regulatory Network Using Large Motif and Track Collections PLoS Comput Biol. 2014 Jul 24;10(7):e1003731. doi: 10.1371/journal.pcbi.1003731. eCollection 2014 Jul.
> 4 votes

# Answer
In general, there is two options to identify targets for transcription factors: experimental (ChIP-seq) and sequence-based predictions.

### TF binding from experimental data

The are multiple projects that produce binding data of transcription factors and quantify their peaks across the genome. The advantage here is that you know that binding actually occurs, as opposed to motif predictions. But you need the corresponding experiments.

You already mentioned ENCODE, which is probably the biggest producer of such data. NIH's Roadmap Epigenomics is another one, but that focusses less on TFs.

If you have a gene list and would like to know which transcription factor was likely involved, you can do an enrichment test of your gene set in known TF targets. The ChEA (ChIP Enrichment Analysis) database does this.

### Prediction using motifs

Another possibility is to look at binding motifs and see whether your gene list is enriched in those. These will, however, be inactive in a given tissue or cell type if the chromatin is packed or the methylation state of the promotor is unfavorable.

Examples for those are the JASPAR and TRANSFAC databases, or the MSigDB motif gene set (as burger mentioned). You can also query those features using the Ensembl genome database (BioMart, REST).

### Calculating enrichment in gene sets

Most likely, an analysis you will want to perform is to calculate the enrichment in your gene list, e.g. that you got from differential expression in two conditions.

The most convenient way is to use the Enrichr platform, which is a web page that accepts a gene list and will compute enrichment in ChEA, JASPAR, TRANSFAC, etc. You can also download their gene sets.
> 4 votes

# Answer
TRUSST is a manually curated database with around 790 TFs and their target genes. TRUSST Website
> 0 votes

---
Tags: chip-seq, encode, multi-omics
---


# Question
Title: How do I find the most similar sequences from a large set of short sequences?
I've been asked to evaluate whether the MinION will be sufficient to distinguish between 20bp CRISPR guide RNAs from the GeCKO v2 set. I know that this set has some exactly identical sequences, which wouldn't be distinguishable regardless of the sequencing method used, so I would like to know whether any two non-indentical sequences differ by less than 15% of the 20bp (i.e. only 3 bases different).

There are around 130k 20bp sequences in the dataset, and I'm pretty sure that an all-vs-all approach is not feasible, but perhaps there's some trick that I'm missing.

What I thought of doing was looking at dimer and trimer/codon counts to hunt down sequences that were most similar. I found about 100 sequences that shared the same trimer count signature (but were not identical), and within this group found \[only\] one pair with 85% identity; all the remainder were 70% or less.

... but I'm having trouble ignoring the doubt I have that these two sequences are the most similar. It might be possible that there are sequences in the dataset that don't have identical trimer count signatures, but differ by less than three bases. I'll give a very obvious example (which isn't in the dataset):

```
1: AAAAA AAAAA AAAAA AAAAA
2: AAAAC AAAAA AAAAA AAAAA

```

The trimer counts for these two sequences are different (Sequence 1 -- AAA: 18; Sequence 2 -- AAA: 15; AAC: 1; ACA: 1; CAA: 1), and yet the sequences are only 1 base different (which might be a problem for the MinION). How can I work out if any such sequence pairs exist in the dataset?

# Answer
130k * 20bp is a small data set. At this scale, SSE2 Smith-Waterman may work well:

```
git clone https://github.com/attractivechaos/klib
cd klib && gcc -O2 -D_KSW_MAIN ksw.c -o ksw -lz
./ksw -a1 -b1 -q1 -r1 -t14 20bp.fa 20bp.fa > out.tsv

```

I simulated 130k 20bp reads from E. coli K-12. The last command line takes about ~2 CPU hours based on the partial output.

It is possible to make this much faster by using 5bp seeds. You can also try ssearch and fasta from the fasta aligner package. They are slower, probably because they are not optimized for such input. Computing Hamming distance with SSE2 will be much faster, but is probably of little use as it does not allow gaps. There are also Gene Myers' edit-distance based algorithms, which can be faster than SSE2-SW.

Bowtie2 (as well as most short-read mappers) won't work well at the default setting. Bowtie2 uses long exact seeds. It will miss many 1-mismatch hits, let alone 3-mismatch ones. The number you get from bowtie2 is an underestimate. You might be able to tune bowtie2, but to find 3-mismatch hits, you have to make it a lot slower.
> 3 votes

# Answer
For this application, you could probably also do something like calculate the Hamming distances between all of the strings in an all-vs-all approach (it should not take too long or too much overhead). You could use something like the Hamming distance tools in Julia. Here is an example of what I mean (using Julia):

```
using StringDistances

k = ["AATTGGCC", "AATTGGCA", "AATTCCGG"]

for s in k
    for y in k
        print(s, " compared to ", y, ":  ", compare(Hamming(), s, y), "\n")
    end
end

```

The output looks like this:

```
AATTGGCC compared to AATTGGCC:  1.0
AATTGGCC compared to AATTGGCA:  0.875
AATTGGCC compared to AATTCCGG:  0.5
AATTGGCA compared to AATTGGCC:  0.875
AATTGGCA compared to AATTGGCA:  1.0
AATTGGCA compared to AATTCCGG:  0.5
AATTCCGG compared to AATTGGCC:  0.5
AATTCCGG compared to AATTGGCA:  0.5
AATTCCGG compared to AATTCCGG:  1.0

```

This is just a simple small example but the core is there. Should not be too bad scaling it up. Hope this helps and provides another helpful option. :)
> 2 votes

# Answer
> I'm pretty sure that an all-vs-all approach is not feasible, but perhaps there's some trick that I'm missing.

The "trick" is that I hadn't actually tried this, and Bowtie2 happens to be quite good at doing this. When there are other similar matches, the MAPQ score that Bowtie2 produces are reduced, so all that is needed is to identify those sequences with less than a MAPQ of 42. I've found 2810 such sequences that share some similarity to others (2416 if I exclude reverse-complement matches):

```
$ bowtie2-build Mouse_GeCKOv2_Library.fasta Mouse_GeCKOv2_Library.fasta
$ bowtie2 -f -x Mouse_GeCKOv2_Library.fasta -U Mouse_GeCKOv2_Library.fasta |
   samtools view | awk '{if($5 != 42){print ">"$1"\n"$10}}' > similar_seqs.fasta
125795 reads; of these:
  125795 (100.00%) were unpaired; of these:
    1 (0.00%) aligned 0 times
    122982 (97.76%) aligned exactly 1 time
    2812 (2.24%) aligned >1 times
100.00% overall alignment rate

```
> 1 votes

---
Tags: minion, crispr
---


title: Aligning many long sequences

text:
# Question
Title: Aligning many long sequences
I'm faced with having to align many (some 100s) bacterial genomes, where the genome length is in the millions. Obviously, this is beyond normal alignment techniques and it's unclear to me what the best practice is for such circumstances:

* conventional alignment on a very powerful computer with lots of memory
* break up the genome into smaller fragments and align them individually
* some exotic different procedure

What possible avenues of attack are there?

(I've attempted to use Mafft and Clustal with little success)

# Answer
> 9 votes
Whole genome aliment can be done using Progressive Mauve, LAST or Mummer. For bacteria I used Mauve since it has also very nice visualisation engine. A very new tool is Minimap2, a super fast mapper that supposed beside read mapping be able to handle reference vs reference. However, I do not know how performance of it compares to the tools mentioned above.

If you are interested in a rough idea of the shared genome regions, you can use bevel. Bevel is not really an aligner, it is more like a dot-plot, but it is super fast (even for mammalian sized genomes).

---
Tags: sequence-alignment, genome
---


# Question
Title: Preparing binary matrix data for Scikit classification algorithms
I made this post in regular stack overflow but I was told about this awesome feature by @nbryans.

I am a researcher (my programming knowledge is small) conducting analysis on a set of antibiotic (methicillin) resistant and a set of antibiotic (methicillin) susceptible. The sole difference between these two sets is assumed to be the resistance (R) vs susceptibility (S). I want to find if a genetic element in the resistant or susceptible genomes is correlated to resistance (positively or negatively).

I have done a complete fisher exact test already and will detail how I did so below. However I would like to figure out how to prepare my data for a more complex learning algorithm such as those offered by sci-kit.

(many of my techniques were learned from this article)

First I used KMC to split each fasta file for both R and S genomes into k-mers (I don't have enough "reputation" so I can't use \>2 links but k-mers are essentially genetic sequences of k length). The output for KMC was a dump of k-mers present in a genome.

Example: kmc dump

```
AAAAAAAAAACGAATGTACACAATCGAA    1
AAAAAAAAACGAATGTACACAATCGAAC    1
AAAAAAAAATCAAATCCTGACTATTTAG    1
AAAAAAAAATCGGTCAATTCATTAAAAG    1
AAAAAAAACAAACATGAAGACCTTGTTA    1
AAAAAAAACGAATGTACACAATCGAACA    1
AAAAAAAACGTGTTAAAGTGAATCACAC    1

```

I then merged the data into one binary matrix which had columns as genomes and rows as k-mers.

Example: 1 resistant genome and 1 susceptible genome

```
kmer    R_1    S_1
AAAAAAAAAACGAATGTACACAATCGAA    1   0
AAAAAAAAACGAATGTACACAATCGAAC    1   0
AAAAAAAAATCAAATCCTGACTATTTAG    1   0
AAAAAAAAATCGGTCAATTCATTAAAAG    1   0
AAAAAAAAATTCCCTTCTAATCTTGAAT    0   1
AAAAAAAACAAAAATTATATAAAGCGAA    0   1
AAAAAAAACAAACATGAAGACCTTGTTA    1   1
AAAAAAAACAACCACCCATACATTGAGT    0   1
AAAAAAAACCCTTACAACAAATATGTAA    0   1

```

This is the data I have been working with and would like to analyze further with sci kit or other classification algorithms. Essentially each k-mer is a feature whose correlation with resistance is based on its presence in R vs S genomes. This is where I would like help transitioning this data into a dataset format suitable for sci kit.

I conducted a Fisher Exact test on each k-mer using R. For each k-mer a 2x2 matrix was created with the first column as resistant and the second column as susceptible. The two rows were present or not. So four numbers (# times k-mer was present in R, # in susceptible, # not in R, # not in S).

```
      R    S
yes   #    #
no    #    #

```

I used the following code in R:

```
phenotype = as.numeric(grepl('^R_', colnames(raw_kmer_table)[2:ncol(raw_kmer_table)]))
fe_results = apply(kmer_pres_abs_matrix, 1,
               FUN = function(row) {
                 fe_mat = matrix(0, ncol = 2, nrow = 2)
                 fe_mat[1,1] = sum(row == 1 & phenotype == 1)
                 fe_mat[1,2] = sum(row == 1 & phenotype == 0)
                 fe_mat[2,1] = sum(row == 0 & phenotype == 1)
                 fe_mat[2,2] = sum(row == 0 & phenotype == 0)

                 fe = fisher.test(fe_mat)
                 return(fe$p.value)
               }
)

```

Now that you have read how my data is currently formatted I would like to hear suggestions on how I can fit this data into a more complex test using sci kit or another resource. Algorithms I am interested in using are any you see fit and previously field-tested algorithms such as adaboost and forest learning algorithms. Also please ask me any clarifications or things I may have left out of this post!

# Answer
Since you are using R, you probably don't want to use scikit-learn, which is for Python. However, there is a similar R library mlr ("R package to make machine learning in R easy") that provides a unified interface to all popular machine learning methods. Check their tutorial on how to get started.
> 2 votes

---
Tags: k-mer, python, machine-learning
---


# Question
Title: Should the cell sorting marker genes be excluded during clustering?
We sort different populations of blood cells using a number of fluorescent flow cytometry markers and then sequence RNA. We want to see what the transcriptome tells us about the similarity and relation between these cells. In my experience on bulk RNA-seq data, there is a very good agreement between flow cytometry and mRNA expression for the markers. It's good to remember we sort using a few markers only, while there are hundreds if not thousands of cell surface proteins.

Should we exclude genes of these sorting (CD) marker proteins when we perform PCA or other types of clustering?

The argument is that after sorting, these genes may dominate clustering results, even if the rest of the transcriptome would tell otherwise, and thus falsely confirm similarity relations inferred from sorting.

Should we at least check the influence of these genes on clustering results?

# Answer
I do not think there is a simple "yes" or "no" answer here.

A good starting point would be, as you suggest, use all the genes and assess the results in the light of the marker genes and expected results. This could both serve as as good quality control as well as give you overview of all the processes happening in the cells.

Depending on the marker genes effect and biological question you may then want to remove the marker genes with potential ordering effect, or even other genes, e.g. by GO terms or pathways.

You most likely want to account for the cell cycle phases as well. And check for other technical factors. I recommend Bioconductor pipeline to get inspired https://www.bioconductor.org/help/workflows/simpleSingleCell/ when it comes to scRNA-seq analyses.
> 3 votes

---
Tags: rna-seq, cell-line, clustering
---


# Question
Title: How to select high quality structures from the Protein Data Bank?
Models of structures deposited in the Protein Data Bank vary in the quality, depending both on the data quality and expertise and patience of the person who built the model. Is there a well-accepted subset of the PDB entries that has only "high quality" structures? Ideally these structures would be representative for classes of proteins in the whole PDB.

based on a real question from biology.SE

# Answer
There is a very nice database, pdbcull (also known as the PISCES server in the literature). It filters the PDB for high resolution and reduced sequence identity. It also seems to be updated regularly. Depending on the cut-offs, you get between 3000 and 35000 structures.

If you are specifically interested in rotamers, you may want to look at top8000 instead, where they have checked for high resolution, and good MolProbity scores. They also provide a rotamer database.

PDB also provides their own clustering. They first cluster the sequences, and then extract a representative structure for each one, based on the quality factor (`1/resolution - R_value`). This has the advantage of being comprehensive, but you will have bad structures when no good ones were ever obtained.
> 9 votes

# Answer
If you choose to perform your own culling of the PDB, resolution is probably the first thing you'll want to look at, which as Davidmh mentions is the main selection criteria for PISCES. High quality structures will also have better R-factor values. You can also give preference based on experimental technique, in descending order of quality:

Neutron diffraction, X-ray diffraction, solution/solid state NMR, electron microscopy/crystallography, fiber diffraction, solution scattering.
> 5 votes

---
Tags: protein-structure, pdb
---


# Question
Title: Get RefSeq accession numbers with versions
Google searching for **NM\_002084** gives the following result:

NM\_002084.4

This, I assume, is the latest version v4, hence the `.4` suffix.

Searching for previous versions I get the following results, along with notes saying it was updated or removed.

NM\_002084.3

> This sequence has been updated. See current version.

NM\_002084.2

NM\_002084.1

> Record removed. This record was replaced or removed.

Using *biomaRt* I can get the latest(**?**) version as follows:

```
library("biomaRt")

# define db
ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")

# get refseqs
getBM(attributes = c('refseq_mrna',
                     'chromosome_name',
                     'transcript_start',
                     'transcript_end',
                     'strand'),
      filters = c('refseq_mrna'),
      values = list(refseq_mrna = "NM_002084"),
      mart = ensembl)
#  refseq_mrna chromosome_name transcript_start transcript_end strand
#1   NM_002084               5        151020438      151028992      1

```

But querying for specific versions gives nothing:

```
getBM(attributes = c('refseq_mrna',
                     'chromosome_name',
                     'transcript_start',
                     'transcript_end',
                     'strand'),
      filters = c('refseq_mrna'),
      values = list(refseq_mrna = c("NM_002084.1", "NM_002084.2", "NM_002084.3", "NM_002084.4")),
      mart = ensembl)
# [1] refseq_mrna      chromosome_name  transcript_start transcript_end   strand
# <0 rows> (or 0-length row.names)

```

**Question:** How can I get all versions (preferably using R)?

# Answer
I don't believe this is possible using `biomaRt`, nor using `AnnotationHub`.

I have two suggestions, neither of them very satisfactory. First, you can specify an Ensembl archive for `biomaRt`, for example:

```
mart72.hs <- useMart("ENSEMBL_MART_ENSEMBL", "hsapiens_gene_ensembl",
                      host = "jun2013.archive.ensembl.org")

```

Of course, that requires that you have some idea of the date for each accession version and that the archives span that date - so not especially useful.

The other option is to access EUtils using *e.g.* `rentrez`, which does allow search by version number:

```
library(rentrez)
es <- entrez_search("nuccore", "NM_002084.1 NM_002084.2")
es$ids

[1] "4504104" "6006000"

```

So knowing the accession, you could simply append 1, 2, 3... to it, run the search and see if UIDs come back, then get them using `entrez_fetch`.

EDIT: here's a quick and dirty function which takes an accession as query, appends version = 1, fetches the ID, then increments the version and repeats until no more results are returned. It is not well-tested!

```
getVersions <- function(accession) {
  require(rentrez)
  ids <- character()
  version <- 1
  repeat({
    es <- entrez_search("nuccore", paste0(accession, ".", version))
    if(length(es$ids) == 0) {
      break
    }
    ids[version] <- es$ids
    version <- version + 1
  })
  ids
}

```

Example:

```
getVersions("NM_002084")
[1] "4504104"    "6006000"    "89903006"   "1048339180"

```
> 6 votes

---
Tags: r, bioconductor, biomart, refseq
---


title: What is the etymology of "Entrez ID"?

text:
# Question
Title: What is the etymology of "Entrez ID"?
Since I have seen NCBI gene names were called "Entrez ID" for the first time, I am wondering where that comes from. Such a weird name!

Does anybody know where that originates?

My hypothesis is: in French, "entrez" can be translated to "please enter" (it is `to enter` in the imperative mood), so I'm guessing that might come from a past website/database in French where `Entrez ID` as in "Please enter the ID" was written in the search bar...

But I could not find any reference... I cannot be the only one to ask himself this!

# Answer
## Update

As pointed out by neilfws:

According to Wikipedia and Practical bioinformatics by Michael Agostino: The name "Entrez" means "Come in!" in French (pronounced “on-tray”); it was chosen to reflect the spirit of welcoming the public to search the content.

## Old (probably incorrect) theory

The first version of Entrez database was distributed by NCBI in 1991 (on CD-ROM). At that time, it consisted of nucleotide sequences, protein sequences and associated citations and abstracts from MEDLINE.

So even at the early stage it contained objects (entries) of very different types (sequences and abstracts). It is hard to find a word which will accurately and briefly describe all the types of data.

While I don't know the correct answer and this is purely a speculation, my guess is that Entrez is an intentionally misspelled word "entries"; which accurately describes all the datatypes in the database.

Using a misspelled word as a name is not uncommon, e.g. Google comes from misspelled "googol" and Flickr is a misspelled version of "flicker".
> 3 votes

---
Tags: gene
---


title: Minion channel ID's from Albacore

text:
# Question
Title: Minion channel ID's from Albacore
The latest version of Albacore from Oxford Nanopore Technologies calls bases from raw fast5 files. A useful piece of output is the `sequence_summary.txt`, which is a big tab-delimited file with information on each read.

One of the columns in there is the `channel`, which references the channel number from 1 to 512 on the flowcell.

These channels obviously have a physical layout on the flowcell. What is the layout? I.e. how many channels are there per row and per column, and how do the channel numbers map to the rows and columns?

My best guess is that it's a 16\*32 array, but I don't know for sure, and even if I did I wouldn't know whether the `channels` map across the shorter or the longer axis in order (or maybe they are not in order at all, but that would be silly).

# Answer
> 3 votes
There are actually 2048 usable sequencing wells, hexagonally packed with four wells connected to the same sequencing sensor/channel via a multiplex (mux) selector. The combination of the mux and the channel number determines the physical location of the well. Unfortunately, the association between channel number and physical location is not obvious, and within each channel the muxes are also not in an obvious order, i.e. `[3,4,1,2,2,1,4,3]` for two adjacent channels.

Here's an "official" text description of the layout as it was for R7; it hasn't changed for the most recent flow cells (R9.5):

* Channels 1-64 occur at the top of the chip (the other low channel numbers are at the bottom)
* Channels order down the chip: 1-64, 449-512, 385-448, 321-384, 257-320, 193-256, 129-192, 65-128
* Muxes run from left to right in the order: 3, 4, 1, 2, 2, 1, 4, 3

To confirm this layout, I pulled the channel/position lookup table definitions out of MinKNOW, and it seems to be the same as this:

Here's the code to produce this image (I haven't adjusted it to have a hexagonal structure):

```
#!/usr/bin/Rscript
library(rjson);

layout.fname <- "/opt/ONT/MinKNOW/Client/resources/app/resources/channels.json";

data.list <- fromJSON(file=layout.fname);

png("MinION_FC_Layout.png", width=3072, height=1536, pointsize=24);
par(mar=c(4,4,0.5,0.5), lwd=3);
data.mat <- matrix(unlist(data.list), ncol=2, byrow=TRUE);
plot(data.mat, pch=21,
     bg=(colorRampPalette(c("white","grey"))(2048))[1:2048],
     col=c("red","yellow","green","blue"),
     xlab="locX", ylab="locY", cex=2,
     xlim=c(7,76), ylim=c(0,33));
text(x=data.mat, labels=rep(1:512, each=4), cex=0.4);
legend("top",legend=c(1,2,3,4), fill=c("red","yellow","green","blue"),
       inset=0.05, title="MUX");
text(x=6, y=data.mat[data.mat[,1] == 8,2],
     labels = ceiling(which(data.mat[,1] == 8) / 4));
text(x=77, y=data.mat[data.mat[,1] == 75,2],
     labels = ceiling(which(data.mat[,1] == 75) / 4));
arrows(x0=c(8,75), x1=c(39,44), y0=33, length=0.2);
invisible(dev.off());

```

# Answer
> 1 votes
poRe has a `show.layout()` function which shows you the 32\*16 grid on which the channels are arranged. EDIT: THIS APPEARS NOT TO BE THE CORRECT PHYSICAL LAYOUT.

```
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
 [1,]  125  121  117  113  109  105  101   97   93    89    85    81    77
 [2,]  126  122  118  114  110  106  102   98   94    90    86    82    78
 [3,]  127  123  119  115  111  107  103   99   95    91    87    83    79
 [4,]  128  124  120  116  112  108  104  100   96    92    88    84    80
 [5,]  253  249  245  241  237  233  229  225  221   217   213   209   205
 [6,]  254  250  246  242  238  234  230  226  222   218   214   210   206
 [7,]  255  251  247  243  239  235  231  227  223   219   215   211   207
 [8,]  256  252  248  244  240  236  232  228  224   220   216   212   208
 [9,]  381  377  373  369  365  361  357  353  349   345   341   337   333
[10,]  382  378  374  370  366  362  358  354  350   346   342   338   334
[11,]  383  379  375  371  367  363  359  355  351   347   343   339   335
[12,]  384  380  376  372  368  364  360  356  352   348   344   340   336
[13,]  509  505  501  497  493  489  485  481  477   473   469   465   461
[14,]  510  506  502  498  494  490  486  482  478   474   470   466   462
[15,]  511  507  503  499  495  491  487  483  479   475   471   467   463
[16,]  512  508  504  500  496  492  488  484  480   476   472   468   464
      [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
 [1,]    73    69    65    61    57    53    49    45    41    37    33    29
 [2,]    74    70    66    62    58    54    50    46    42    38    34    30
 [3,]    75    71    67    63    59    55    51    47    43    39    35    31
 [4,]    76    72    68    64    60    56    52    48    44    40    36    32
 [5,]   201   197   193   189   185   181   177   173   169   165   161   157
 [6,]   202   198   194   190   186   182   178   174   170   166   162   158
 [7,]   203   199   195   191   187   183   179   175   171   167   163   159
 [8,]   204   200   196   192   188   184   180   176   172   168   164   160
 [9,]   329   325   321   317   313   309   305   301   297   293   289   285
[10,]   330   326   322   318   314   310   306   302   298   294   290   286
[11,]   331   327   323   319   315   311   307   303   299   295   291   287
[12,]   332   328   324   320   316   312   308   304   300   296   292   288
[13,]   457   453   449   445   441   437   433   429   425   421   417   413
[14,]   458   454   450   446   442   438   434   430   426   422   418   414
[15,]   459   455   451   447   443   439   435   431   427   423   419   415
[16,]   460   456   452   448   444   440   436   432   428   424   420   416
      [,26] [,27] [,28] [,29] [,30] [,31] [,32]
 [1,]    25    21    17    13     9     5     1
 [2,]    26    22    18    14    10     6     2
 [3,]    27    23    19    15    11     7     3
 [4,]    28    24    20    16    12     8     4
 [5,]   153   149   145   141   137   133   129
 [6,]   154   150   146   142   138   134   130
 [7,]   155   151   147   143   139   135   131
 [8,]   156   152   148   144   140   136   132
 [9,]   281   277   273   269   265   261   257
[10,]   282   278   274   270   266   262   258
[11,]   283   279   275   271   267   263   259
[12,]   284   280   276   272   268   264   260
[13,]   409   405   401   397   393   389   385
[14,]   410   406   402   398   394   390   386
[15,]   411   407   403   399   395   391   387
[16,]   412   408   404   400   396   392   388

```

---
Tags: nanopore, minion, albacore
---


# Question
Title: How to select a cutoff for interaction confidence in STRINGdb?
I have a list of 100 genes that are called as hits in a genetic screening. I want to have a network of the interactions between the proteins of these 100 genes. I am using both STRINGdb web and its R API.

In both situations you can select the threshold of confidence of interaction (default 0.4) for displaying the network. Changing this parameter drastically changes the network that I am generating. I have read the FAQ section of STRINGdb and they recommend to choose some arbitrary number based on the number of interactions you need for you analysis. If I use the default threshold, I have few interactions and I don't know if lowering the threshold would be correct.

Therefore, my question is, is there any established threshold based on habitual user's experience used for that, beyond the recommendation of setting an arbitrary one?

# Answer
> 5 votes
I have used STRING pretty heavily, and have compared it to various other databases of protein interactions and signaling pathways. I do feel like it has a lot of quality interaction annotations, but you have to sift through a lot of noise to get to them. The simplest method I have found for doing this is to look at the individual scores for each interaction, and accept it if it passes one of the following tests:

* Experiment Score \> 0.4
* Database Score \> 0.9

Anything that passes one of these thresholds we consider at least an interaction of acceptable quality. Those interactions with Experiment Scores \> 0.9 are high-quality, and have been experimentally validated. The other scores represent inaccurate methods for determining signaling events, and should be ignored.

You are not going to catch every actual protein signaling event this way, but you will at least be spared a lot of false positives. The best way to construct an actual network of signaling events is to combine interaction records from multiple databases.

# Answer
> 5 votes
This depends on what you are trying to do and whether you value specificity over sensitivity. We can't tell you since it is entirely dependent on the biological question you want to answer.

However, I would recommend two things:

1. Don't use STRING. The creators of STRING made the choice to value sensitivity over all else, so they include *any* interaction they can get their hands on. That is sometimes fine, depending on what you want to do, but is more often a problem. Personally, I tend to avoid STRING as much as possible and instead use more curated databases like APID or IntAct. You can find a list of many PPI databases in the EBI's Psiqcuick View page. You might also be interested in my answer here which gives an example script using APID which can easily be modified to query other DBs.
2. Whatever you use, I recommend you filter by detection method. There are many different interaction detection methods ranging from those that only identify direct, binary interactions between two proteins (e.g. Yeast two-hybrid), those that also find proteins that are in the same complex even if there is no *direct* interaction between them (e.g. ChIP), to various non-experimental methods which are used to *infer* interactions.

Luckily, all of these have been collected into a controlled vocabulary which you can see in the Interaction Detection Method Ontology Lookup Service page. Go there, collect the interaction detection methods you consider good enough for whatever you want to do and filter your interactions using those.

For example, for a high quality human interactome I built for a project I used to work on, I was only interested in direct binary interactions, so I only kept interactions detected by the following methods:

```
    "MI:0008" => "array technology",
    "MI:0009" => "bacterial display",
    "MI:0010" => "beta galactosidase complementation",
    "MI:0011" => "beta lactamase complementation",
    "MI:0012" => "bioluminescence resonance energy transfer",
    "MI:0013" => "biophysical",
    "MI:0014" => "adenylate cyclase complementation",
    "MI:0016" => "circular dichroism",
    "MI:0017" => "classical fluorescence spectroscopy",
    "MI:0018" => "two hybrid",
    "MI:0020" => "transmission electron microscopy",
    "MI:0030" => "cross-linking study",
    "MI:0031" => "protein cross-linking with a bifunctional reagent",
    "MI:0034" => "display technology",
    "MI:0040" => "electron microscopy",
    "MI:0041" => "electron nuclear double resonance",
    "MI:0042" => "electron paramagnetic resonance",
    "MI:0043" => "electron resonance",
    "MI:0047" => "far western blotting",
    "MI:0048" => "filamentous phage display",
    "MI:0049" => "filter binding",
    "MI:0051" => "fluorescence technology",
    "MI:0052" => "fluorescence correlation spectroscopy",
    "MI:0053" => "fluorescence polarization spectroscopy",
    "MI:0055" => "fluorescent resonance energy transfer",
    "MI:0065" => "isothermal titration calorimetry",
    "MI:0066" => "lambda phage display",
    "MI:0073" => "mrna display",
    "MI:0081" => "peptide array",
    "MI:0084" => "phage display",
    "MI:0089" => "protein array",
    "MI:0090" => "protein complementation assay",
    "MI:0091" => "chromatography technology",
    "MI:0092" => "protein in situ array",
    "MI:0095" => "proteinchip(r) on a surface-enhanced laser desorption/ionization",
    "MI:0097" => "reverse ras recruitment system",
    "MI:0098" => "ribosome display",
    "MI:0099" => "scintillation proximity assay",
    "MI:0107" => "surface plasmon resonance",
    "MI:0108" => "t7 phage display",
    "MI:0111" => "dihydrofolate reductase reconstruction",
    "MI:0112" => "ubiquitin reconstruction",
    "MI:0114" => "x-ray crystallography",
    "MI:0115" => "yeast display",
    "MI:0226" => "ion exchange chromatography",
    "MI:0227" => "reverse phase chromatography",
    "MI:0231" => "mammalian protein protein interaction trap",
    "MI:0232" => "transcriptional complementation assay",
    "MI:0255" => "post transcriptional interference",
    "MI:0369" => "lex-a dimerization assay",
    "MI:0370" => "tox-r dimerization assay",
    "MI:0397" => "two hybrid array",
    "MI:0398" => "two hybrid pooling approach",
    "MI:0399" => "two hybrid fragment pooling approach",
    "MI:0400" => "affinity technology",
    "MI:0401" => "biochemical",
    "MI:0405" => "competition binding",
    "MI:0406" => "deacetylase assay",
    "MI:0410" => "electron tomography",
    "MI:0411" => "enzyme linked immunosorbent assay",
    "MI:0415" => "enzymatic study",
    "MI:0416" => "fluorescence microscopy",
    "MI:0419" => "gtpase assay",
    "MI:0420" => "kinase homogeneous time resolved fluorescence",
    "MI:0423" => "in-gel kinase assay",
    "MI:0424" => "protein kinase assay",
    "MI:0425" => "kinase scintillation proximity assay",
    "MI:0426" => "light microscopy",
    "MI:0428" => "imaging technique",
    "MI:0432" => "one hybrid",
    "MI:0434" => "phosphatase assay",
    "MI:0435" => "protease assay",
    "MI:0437" => "protein three hybrid",
    "MI:0440" => "saturation binding",
    "MI:0508" => "deacetylase radiometric assay",
    "MI:0509" => "phosphatase homogeneous time resolved fluorescence",
    "MI:0510" => "homogeneous time resolved fluorescence",
    "MI:0511" => "protease homogeneous time resolved fluorescence",
    "MI:0512" => "zymography",
    "MI:0513" => "collagen film assay",
    "MI:0514" => "in gel phosphatase assay",
    "MI:0515" => "methyltransferase assay",
    "MI:0516" => "methyltransferase radiometric assay",
    "MI:0655" => "lambda repressor two hybrid",
    "MI:0657" => "systematic evolution of ligands by exponential enrichment",
    "MI:0678" => "antibody array",
    "MI:0695" => "sandwich immunoassay",
    "MI:0696" => "polymerase assay",
    "MI:0726" => "reverse two hybrid",
    "MI:0727" => "lexa b52 complementation",
    "MI:0728" => "gal4 vp16 complementation",
    "MI:0809" => "bimolecular fluorescence complementation",
    "MI:0813" => "proximity enzyme linked immunosorbent assay",
    "MI:0824" => "x-ray powder diffraction",
    "MI:0825" => "x-ray fiber diffraction",
    "MI:0827" => "x-ray tomography",
    "MI:0841" => "phosphotransferase assay",
    "MI:0870" => "demethylase assay",
    "MI:0872" => "atomic force microscopy",
    "MI:0879" => "nucleoside triphosphatase assay",
    "MI:0880" => "atpase assay",
    "MI:0887" => "histone acetylase assay",
    "MI:0889" => "acetylase assay",
    "MI:0892" => "solid phase assay",
    "MI:0894" => "electron diffraction",
    "MI:0895" => "protein kinase A complementation",
    "MI:0899" => "p3 filamentous phage display",
    "MI:0900" => "p8 filamentous phage display",
    "MI:0905" => "amplified luminescent proximity homogeneous assay",
    "MI:0916" => "lexa vp16 complementation",
    "MI:0920" => "ribonuclease assay",
    "MI:0921" => "surface plasmon resonance array",
    "MI:0946" => "ping",
    "MI:0947" => "bead aggregation assay",
    "MI:0949" => "gdp/gtp exchange assay",
    "MI:0953" => "polymerization",
    "MI:0968" => "biosensor",
    "MI:0969" => "bio-layer interferometry",
    "MI:0972" => "phosphopantetheinylase assay",
    "MI:0976" => "total internal reflection fluorescence spectroscopy",
    "MI:0979" => "oxidoreductase assay",
    "MI:0984" => "deaminase assay",
    "MI:0989" => "amidase assay",
    "MI:0990" => "cleavage assay",
    "MI:0991" => "lipid cleavage assay",
    "MI:0992" => "defarnesylase assay",
    "MI:0993" => "degeranylase assay",
    "MI:0994" => "demyristoylase assay",
    "MI:0995" => "depalmitoylase assay",
    "MI:0996" => "deformylase assay",
    "MI:0997" => "ubiquitinase assay",
    "MI:0998" => "deubiquitinase assay",
    "MI:0999" => "formylase assay",
    "MI:1000" => "hydroxylase assay",
    "MI:1001" => "lipidase assay",
    "MI:1002" => "myristoylase assay",
    "MI:1003" => "geranylgeranylase assay",
    "MI:1004" => "palmitoylase assay",
    "MI:1005" => "adp ribosylase assay",
    "MI:1006" => "deglycosylase assay",
    "MI:1007" => "glycosylase assay",
    "MI:1008" => "sumoylase assay",
    "MI:1009" => "desumoylase assay",
    "MI:1010" => "neddylase assay",
    "MI:1011" => "deneddylase assay",
    "MI:1016" => "fluorescence recovery after photobleaching",
    "MI:1019" => "protein phosphatase assay",
    "MI:1024" => "scanning electron microscopy",
    "MI:1026" => "diphtamidase assay",
    "MI:1030" => "excimer fluorescence",
    "MI:1031" => "protein folding/unfolding",
    "MI:1036" => "nucleotide exchange assay",
    "MI:1037" => "Split renilla luciferase complementation",
    "MI:1038" => "silicon nanowire field-effect transistor",
    "MI:1087" => "monoclonal antibody blockade",
    "MI:1088" => "phenotype-based detection assay",
    "MI:1089" => "nuclear translocation assay",
    "MI:1111" => "two hybrid bait or prey pooling approach",
    "MI:1112" => "two hybrid prey pooling approach",
    "MI:1113" => "two hybrid bait and prey pooling approach",
    "MI:1137" => "carboxylation assay",
    "MI:1138" => "decarboxylation assay",
    "MI:1142" => "aminoacylation assay",
    "MI:1145" => "phospholipase assay",
    "MI:1147" => "ampylation assay"

```

If you also want those that might not be direct binary but simply in the same complex, include:

```
      "MI:0019" => "coimmunoprecipitation",
      "MI:0006" => "anti bait coimmunoprecipitation",
      "MI:0007" => "anti tag coimmunoprecipitation",
      "MI:0858" => "immunodepleted coimmunoprecipitation",
      "MI:0096" => "pull down",
      "MI:0963" => "interactome parallel affinity capture",
      "MI:0676" => "tandem affinity purification"

```

---
Tags: proteins, public-databases, networks
---


# Question
Title: What are the ways to keep track of branches in the analysis?
I'm going through an RNA-seq pipeline in `R/Bioconductor` and want to try multiple parameters at subsequent steps, for example, running clustering with different settings, running RegressOut or not on unwanted effects etc. That's a lot of "versions", even if I don't do combinations of these steps.

How can I keep track of this, and my conclusions? Not necessarily want to save the results.

* Save the different scripts with git (seems overkill)
* Make notes in the script itself

# Answer
The main purpose of git is to version code, which usually means sequential improvement of the codebase. While it is possible to use branches for multiple variants of the software, permanent branches are traditionally used for gradual integration of new features (i.e. dev/testing/master branches). Supporting multiple independent branches requires some investment, i.e. distributing common changes among branches via merge or cherry-pick. This is hard to manage when you have more than two-three branches.

If you compare different methods of analysis, you probably want to compare the results between methods. Having them on separate branches makes it hard.

In my opinion you should integrate all methods of analysis into the master branch. To avoid copy & paste, it is better to put common code in a library or an independent script. You can also specify a method as a run-time parameter of your pipeline, and create a meta-script which will execute all methods of interests.

Once you performed benchmarking, you shouldn't remove unused methods from you master branch. Having them is important for reproducible research, and your scripts could be used in the future for new datasets.
> 5 votes

# Answer
> Save the different scripts with git (seems overkill)

Whoa. I did an actual double take when reading this:1 it’s the *opposite* of overkill.

Version controlling your scripts (using Git or something similar) is the absolute minimum, and should become completely automatic. For every new project I begin, one of the very first steps is to issue the `git init` command, and to set up a remote repository (on Github).

To keep track of *different* analyses, I use a combination of the following approaches:

1. Write reusable functions/scripts and parametrise. The parameters are kept either inside the script itself (that then calls the relevant function repeatedly), or in a declarative workflow.2
2. Document the alternative analysis approaches; once again, this could be a declarative workflow with different rules for alternative analyses, or a set of notebooks (via R Markdown).
3. Have different Git branches for mutually exclusive approaches. At the end of the analysis one of these branches gets merged into the main branch, and published. If I want to publish several analysis approaches, I merge all these branches into the main branch, and use approaches (1) or (2) enable them simultaneously.

In fact, I recommend creating a declarative workflow script for *every* analysis; I have found that this is the most practical, self-documenting way to run an analysis. It most closely resembles a wet-lab lab notebook. The advantage over a single R Markdown document is that rerunning just parts of the analysis can be completely automated, and dependencies in the workflow are apparent from the dependencies of the workflow rules. This is much harder in R Markdown.

Some time ago I create an example analysis workflow to show how this can be structured.

Regarding your other point:

> Make notes in the script itself

“Notes” are a dangerous beasts: documentation is important, but experience shows that it’s sometimes very hard to keep documentation synchronised with the code. There is no mechanism that ensures that documentation and code actually agree. Differences between presumed analysis and actually executed analysis can become very problematic.

People therefore prefer using self-documenting code as much as possible; that is: writing code so that its meaning becomes immediately clear, without comments, even to somebody who hasn’t worked on the code before. Doing this well is hard and takes practice, but improves overall code quality. Once again, using a declarative workflow helps here because the dependencies between the rules are self-documenting the kind of analysis that was performed.

Jeff Atwood has written two seminal essays on this subject:

They are two of the best pieces of advice on programming that I could give.

---

1 To emphasise: take a look at the edit history of this answer.

2 I recommend ‘targets’ (for R), Nextflow or Snakemake. I used to use (and recommend) GNU make, but I no longer believe this to be an adequate tool.
> 9 votes

# Answer
I quite agree with this answer by @Konrad Rudolph.

I would like to emphasis that using parametrization for your scripts is what will help you avoid multiplying branches in git. So yes, use git, but you don't necessarily have to go "overkill" creating lots of branches.

Then, I would command these scripts from a workflow managment tool that will somehow take care of generating the branches of your analyses. If you use Snakemake, the various options taken along the path from your data to your results will be represented by the wildcards system, and this will be visible in the structure of your folders and file names, due to the fact that Snakemake works by inferring what should be done to produce a file based on its name.

This of course is not an excuse for not using other documentation approaches: Comments in the Snakefile and in the scripts, README files explaining how the workflow was run, using what configuration file. Put your scripts, the Snakefile, its configuration files and the README files under version control and document again using commit messages.
> 1 votes

---
Tags: r, bioconductor, git, workflow-management, best-practice
---


# Question
Title: Techniques for analyzing and quantifying sample bleed through in genotyping with Illumina
I am looking at the presence of viral genotypes within individual samples within an assay. Often times there is a sample whose read counts are firing off the charts and this sample tends to "bleed through" to the other samples.

I have recently tried confronting this problem with chi square analysis of conditional probability. For instance let's say I observe with 10 sample viral genotypes A and B. I then calculate the expected number of samples having both genotypes A and B from the number of A positive samples and number of B positive samples.

So far this is all I have. I do not know if anyone else has any ideas on how to determine and quantify bleed through in viral genotyping.

Slightly more information about this situation: these samples were run on Illumina HiSeq, we are planning on doing another HiSeq run with even more samples shortly so would like to take proper precautions with that, additionally the current pipeline utilizes Novobarcode as a demultiplexer so if anyone has suggestions on what I can implement during processing to minimize the effect of cross contamination that'd be helpful.

# Answer
Are you using dual-index barcodes with different barcodes at each end? There is a known phenomena of "index switching" that occurs in Illumina reads. One way to control for this is to add in unamplified water samples that have unused barcode combinations where each of the indexes is shared with another sample.

More information can be found on this SeqAnswers thread. Here are some useful associated links:

* Signal-spreading preprint -- detailed investigation of the phenomenom, including methods, graphs, and associated data
* Crossblock -- tool by Brian Bushnell that attempts to remove cross contamination
* Summary by James Hadfield -- mentions a few fixes suggested by Illumina, many centred around not using barcode combinations that could conflict
* Illumina's quick-fire response -- they've known about it for a long time, and point out that it shouldn't be a problem for most purposes

I have suspicions that at least some of the index switching is being caused by ligation of two separate barcoded fragments during the sample prep (which definitely happens during sample prep for nanopore), but I haven't done any experimentation to work out if this is also the case for Illumina (due to money/time/goal constraints). The experimental design would be to take a prepared Illumina library just prior to sequencing, ligate on ONT adapters, then sequence on both the ONT and Illumina machines to see if chimeric reads discovered via nanopore sufficiently matches the index switching observed via short-read sequencing.
> 1 votes

---
Tags: genotyping
---