# Question
Title: How to achieve blast results according to the intuitive interpretation of `-max_target_seqs`?
Very recently a BLAST parameter `-max_target_seqs n` got a lot of attention. Instead of the intuitive interpretation (return the best `n` sequences) the parameters asks blast to return the first `n` sequences that pass the e-value threshold. This is affecting thousands of workflows and analyses that assumed the intuitive interpretation.

I also found that there is another parameter called `-num_alignments`. What is the difference? Related question on BioStars seems that is misinterpreting the `max_target_seqs` parameter. Is `num_alignments` also misunderstood?

Now the main question. How to actually run blast to give me the best `n` hits in the database?

Till now I thought that this example will do the job:

```
blastp -query $QUERY -db $PROTEINS -out $BLASTOUT -evalue 1e-10 -outfmt 6 -num_alignments 5

```

-- edit --

Now I understood that i can not achieve the desired output only by parameters of blast, but I still don't see an answer to how can I get the best n hits. So I wanted to check out how to get all the results and that's not intuitive either. The parameter `-max_target_seqs` is set by default to 500. Does it mean that if there are more than 500 significant hits, I don't have a guarantee for the best one? Do I have to specify `max_target_seqs` to some crazy high number to be sure that I got them all?

# Answer
The `-max_target_seqs` is also related to `-evalue` which does not work as one would think. Looking at the blast news one finds that since version BLAST 2.2.27+ released, one should use the `-max_target_seqs`:

> 4.) The output for reports without separate descriptions and alignments sections (all –outfmt greater than 4) should now use –max\_target\_seqs to control the output rather than –num\_descriptions and –num\_alignments.

So it seems that `num_alignments` is also misunderstood. I came to the conclusion that one cannot expect to get the best `n` alignments, unless you get all and filter yourself the best ones.
> 3 votes

# Answer
As a follow up to @Llopis' answer (no native blast solution), I wrote a script blast\_filter.py and I added it into my collection of useful scripts generic\_genomics. The script filters out the blast output and it works with streams. To get the `n` (5) best blast hits by e-value you can run:

```
blastp -query $QUERY -db $PROTEINS -evalue 1e-10 -outfmt 6 | blast_filter.py 5 > $BLASTOUT

```

I assume that the blast output is sorted primarily by query name and secondarily by e-value. Also, note that default `max_target_seqs` in blast is 500, if there will be a query with 500 alignments, script will drop a warning (the query has no guarantee of the best hit found).
> 1 votes

---
Tags: blast
---


title: Simulating DNA sequence evolution in R

text:
# Question
Title: Simulating DNA sequence evolution in R
I hope someone can lend their thoughts on the below code to generate DNA sequences under the Kimura-2-Parameter model of DNA substitution.

The issue is that each time the code is run and the haplotype distribution is examined, there always is a very skewed distribution.

What I would like is for a different distribution to be generated each time the code is run. that is, there should be multiple specimens sharing the same haplotype (instead of most sharing the most common haplotype).

Here is my code:

```
library(pegas)
set.seed(17)

sim.seqs <- TRUE
length.seqs <- 500
num.seqs <- 100 # number of DNA sequences
subst.model <- "K80" # nucleotide substitution model
transi.rate <- 1e-4 # transition rate
transv.rate <- transi.rate / 2 # transversion rate

if (sim.seqs == TRUE) {

  nucl <- as.DNAbin(c('a','c','g','t'))
  res <- sample(nucl, size = length.seqs, replace = TRUE, prob = rep(0.25, 4))

 if (subst.model == "K80") {
    transi.set <- list('a' = as.DNAbin('g'),
                       'c' = as.DNAbin('t'),
                       'g' = as.DNAbin('a'),
                       't' = as.DNAbin('c'))

    transv.set <- list('a' = as.DNAbin(c('c', 't')),
                       'c' = as.DNAbin(c('a', 'g')),
                       'g' = as.DNAbin(c('c', 't')),
                       't' = as.DNAbin(c('a', 'g')))

    transi <- function(res) {
      unlist(transi.set[as.character(res)])
    }

    transv <- function(res) {
      sapply(transv.set[as.character(res)], sample, 1)
    }

    duplicate.seq <- function(res) {
      num.transi <- rbinom(n = 1, size = length.seqs, prob = transi.rate) # total number of transitions

      if (num.transi > 0) {
        idx <- sample(length.seqs, size = num.transi, replace = FALSE)
        res[idx] <- transi(res[idx])
      }

      num.transv <- rbinom(n = 1, size = length.seqs, prob = transv.rate) # total number of transversions

      if (num.transv > 0) {
        idx <- sample(length.seqs, size = num.transv, replace = FALSE)
        res[idx] <- transv(res[idx])
      }

      res
      }

    }

  res <- matrix(replicate(num.seqs, duplicate.seq(res)), byrow = TRUE, nrow = num.seqs)

  class(res) <- "DNAbin"

  # write.dna(res, file = "seqs.fas", format = "fasta")

  h <- sort(haplotype(res), decreasing = TRUE, what = "frequencies")
  rownames(h) <- 1:nrow(h)

}

```

## Output

```
h

Haplotypes extracted from: res

Number of haplotypes: 5
     Sequence length: 500

Haplotype labels and frequencies:

 1  2  3  4  5
96  1  1  1  1

```

If the code is run multiple times (without the seed), the same pattern emerges: h is always skewed toward the most dominant haplotype.

For example, I want to be able to run the code (without seed) and get something like

```
h

Haplotypes extracted from: res

Number of haplotypes: 5
     Sequence length: 500

Haplotype labels and frequencies:

  1    2    3    4   5
  35  25   20   15   5

```

i think I would have to specify a distribution for the haplotypes, but the number of haplotypes generated by the mutation process is not known beforehand.

Any thoughts?

# Answer
> 1 votes
What are you expecting these lines to do?

```
num.transi <- rbinom(n = 1, size = length.seqs, prob = transi.rate)
...
num.transv <- rbinom(n = 1, size = length.seqs, prob = transv.rate)

```

Using the variables that are specified in your code, they are extremely zero-skewed. This is not too surprising given the mutation rate and sequence length.

I expect that increasing the sequence length, or increasing the transition rate, should result in more haplotypes being generated.

After running this through myself, I see that increasing the mutation rate does increase the number of haplotypes, but the skew is retained. This is because each mutation creates a new haplotype, following something like an infinite allele model.

I guess this is not particularly surprising. It's unlikely that two random mutations from a particular haplotype will result in the same haplotype, and therefore any mutations (from *any* existing haplotype) will generally create something new.

---
Tags: r, phylogenetics
---


# Question
Title: Low percentage of reads with consistent barcodes in Split-Seq run
I am following this answer to select reads with allowed barcodes from a Split-Seq run. In particular, I have a whitelist for each round of barcodes (3 rounds) and I filter out reads which have more than 1 Hamming distance from allowed barcodes.

However, after filtering, the number of reads remaining are only the 0.7% of the original number of reads before filtering (169k from 22M). Is this normal, or what may be the underlying issue here?

# Answer
> 0 votes
In my case the problem was due a wrong set of allowed round1 barcodes, therefore the filtering pipeline was removing most part of the reads. In particular, with the right barcode1 whitelist the total number of reads after filtering is 17M (about 100 times more).

---
Tags: scrnaseq, reads, splitseq
---


# Question
Title: Extracting sequences from FASTA beginning with common 5' end
I am trying to figure out the best way to extract sequences from a FASTA file which begin with a common 5' region of 43 nucleotides. Preferably, I would like to to allow for "fuzziness" in this region to allow for mutations or read overlaps. The idea is that any sequences that begin with this region, regardless of size, could be extracted into a new file.

Example of the FASTA file:

```
>1-69050-454.08
GTACGGGGAAGGACGTCAATAGTC

>2-65989-433.95
AATCTGTAGTACGACTCACTATAGC

>3-62181-408.91
GATCTGTAATACGACTCACTATAGG

>4-49959-328.53
GGGGAAGGACGTCAATAGTCACAC

```

And what I would like to get from the code is such:

```
>2-65989-433.95
AATCTGTAGTACGACTCACTATAGC

>3-62181-408.91
GATCTGTAATACGACTCACTATAGG

```

The FASTA file is quite large and I have tried using a couple grep and awk methods retrieve the sequences. Any help you could provide is much appreciated.

# Answer
> 2 votes
Here is a 'quick' `python3` solution, requiring installation of BLASTN and a `pip3 install biopython`. The results are output in FASTA format to `stdout`.

Save the following `python3` code in `hit_finder.py`

```
#!/usr/bin/env python3

import sys
from Bio import SeqIO
from Bio.Blast.Applications import NcbiblastnCommandline

cline = NcbiblastnCommandline(query=sys.argv[1],
                              subject=sys.argv[2],
                              outfmt="'6 qseqid length qstart'",
                              max_hsps=1,
                              word_size=4)

result = [j[0] for j in [i.split('\t') for i in
                         cline()[0].rstrip().split('\n')]
          if int(j[1]) >= int(sys.argv[3]) and int(j[2]) <= int(sys.argv[4])]

for record in SeqIO.parse(sys.argv[1], 'fasta'):
    if record.id in result:
        print(record.format('fasta'))

```

Make it executable: `chmod +x hit_finder.py`.

This script takes four arguments (which in future versions could be enhanced using `argparse` instead of `sys.argv`):

```
arg1 = query_fasta_file (file to search in)
arg2 = subject_fasta_file (database file to contain the sequence to match at 5')
arg3 = min_alignment_length (this means the length of the hit, allows mismatching in string with blastn and a word_size of 4)
arg4 = start_of_hit_at_5prime (where should the alignment hit start?)

```

Let's say you wanted to search for `AATCTGTAGTACGACTCACTATAGC` in your query seqeuences. Save this sequence in `subject.fa` as:

```
>subject
AATCTGTAGTACGACTCACTATAGC

```

Your query sequences are stored in `queries.fa` as:

```
>1-69050-454.08
GTACGGGGAAGGACGTCAATAGTC

>2-65989-433.95
AATCTGTAGTACGACTCACTATAGC

>3-62181-408.91
GATCTGTAATACGACTCACTATAGG

>4-49959-328.53
GGGGAAGGACGTCAATAGTCACAC

```

Now, to search for hits of at least 23 bp starting no further than 1 base from the 5' end of your sequence, you would run this search using:

```
./hit_finder.py queries.fa subject.fa 23 1

```

To search for hits of at least 30 bp no more than 10 bp from the 5' end of your queries, do:

```
./hit_finder.py queries.fa subject.fa 30 10

```

To save your results, output `stdout` to file, e.g.,:

```
./hit_finder.py queries.fa subject.fa 30 10 > results.fa

```

---
Tags: ngs, fasta, awk
---


# Question
Title: What does it mean for an amino acid to match with itself 99% of the time in a PAM1 matrix?
In this video, the lecturer says:

> You see that phenylalanine is probably still phenylalanine because only one percent of amino acids have changed so roughly 99 percent of the time, phenylalanine is still phenylalanine. Now why is it not exactly 99 percent? Because it’s one percent for all 20 amino acids so some are greater than 99 percent and some are less than 99 percent. It averages out to 99 percent.

What on earth does he mean by this? I either do not have the background knowledge or am not thinking hard enough, **or** I need this to be explained in a way that is more understandable.

# Answer
> 0 votes
It means that the 99% of conservation of amino acids is just a mean value for all the amino acids. That there can be amino acids that are less conserved and some that are more conserved.

So from the 21 basic amino acids, arginine is more conserved than lysine in those alignments. But in general they are conserved the 99% of times.

---
Tags: statistics
---


# Question
Title: Random addition method for phylogenetic tree reconstruction
I have been working on nucleotide and protein data to reconstruct maximum parsimony trees using MEGA software to get parsimony tree, via "random addition" method is used to generate the initial tree(s).

I need to know that when randomly chosen sequences are added to reconstruct an initial tree, on what basis the best tree is chosen??

I would appreciate if anyone can guide me about this method or suggest me some good literature about this in order to make me understand the concept of this method.

# Answer
MEGA manual may help to explain the criteria used in selecting the final maximum parsimony tree.

The random tree additions refers to the generation of random initial trees for a heuristic search to find MP trees. In this case, a tree is generated by randomly selecting a sequence and adding it to the growing tree on a randomly-selected branch.

Instead of using the random tree additions you can use the min-mini algorithm for the initial tree set.
> 2 votes

---
Tags: phylogenetics, phylogeny
---


# Question
Title: Why are my Chi-squared test results different from those in a published table?
I recently read the paper “A novel long non-coding RNA *linc-ZNF469-3* promotes lung metastasis through *miR-574-5p-ZEB1* axis in triple negative breast cancer”. In this I see Table1 showing correlation of Linc-ZNF469-3 with different features in TNBC and non-TNBC patients. p-values are calculated based on Pearson's Chi-sq test.

For non-TNBC patients: I calculated the p-value with the Chi-sq test just to have a look if everything is right.

Recurrence:

```
data <- data.frame(x= c(40,62), y=c(56,65))
data
chisq.test(data, correct = F)

    Pearson's Chi-squared test

data:  data
X-squared = 1.1269, df = 1, p-value = 0.2884

```

Brain:

```
data <- data.frame(x= c(6,99), y=c(6,119))
data
chisq.test(data, correct = F)

    Pearson's Chi-squared test

data:  data
X-squared = 0.096461, df = 1, p-value = 0.7561

```

You can see that the Recurrence p-value doesn't match with what I see in the table. Others, like Brain, Bone, Liver, Lung match with the p-values I calculated.

And for TNBC patients, not a single p-value matches with my calculations. I don't even get Lung and Bone as having significant p-values with my calculations.

I did it for Lung as shown here:

```
data <- data.frame(x= c(25,143), y=c(9,35))
data
chisq.test(data, correct = F)

    Pearson's Chi-squared test

data:  data
X-squared = 0.8044, df = 1, p-value = 0.3698

```

But in the paper they mentioned:

> We found that high linc-ZNF469-3 expression was significantly correlated to lung metastasis in TNBC patients (Table 1).

Can anyone tell me what could be the wrong here? Is it my calculation wrong? or in the paper?

# Answer
> 2 votes
Your calculations seem right, perhaps there was an error on their side.

I also looked at the number of samples reported, but they use the same amount of samples in each case. Because they are testing whether Linc-ZNF469-3 is different in several tissues, they should have used a multiple test correction for the p-value (although reporting the original values is usually also ok).

To answer your question in the comments. You can't correlate a categorical variable with a continuous variable. Or they did another test or they are not explaining clearly what they did.

---
Tags: r, rna-seq, statistics, gene-expression, correlation
---


# Question
Title: finding unique headers in a fasta file using linux command line
I tried to use the following command

```
uniq -u reference.fasta >> reference_uniq.fasta

```

I'd like a count of the unique headers.

# Answer
> 5 votes
If you just want the number of unique headers, you can do this:

```
grep '>' reference.fasta | sort | uniq | wc -l

```

If you want a list of the unique headers, you can do this:

```
grep '>' reference.fasta | sort | uniq

```

If you want a histogram of how many times each header occurs, you can do this:

```
grep '>' reference.fasta | sort | uniq -c | awk '{printf("%s\t%s\n", $1, $2)}'

```

# Answer
> 3 votes
The `uniq` command expects sorted input. Interestingly, the `sort` command actually has a "unique" option, `-u`, which means `uniq` is not strictly needed. For the fastest processing, you can look for the '\>' character at the start of lines with grep:

```
grep '^>' reference.fasta | sort -u > reference_headers_unique.fasta

```

For returning the *number* of unique lines, pipe through `wc -l`:

```
grep '^>' reference.fasta | sort -u | wc -l

```

For more information about regular expressions, see here.

# Answer
> 1 votes
If, for some reason, you need to keep the original order of the input file, so you only want the *first* occurrence of each header, you can do:

```
awk '/^>/ && !a[$0]++' reference.fasta > headers

```

Or, to get the number only:

```
awk '/^>/ && !a[$0]++{k++}; END{print k}' reference.fasta

```

However, these only make sense if you need to keep the order of the original file. If not, use Gringer's `grep` approach which will be much faster (`reference.fasta` is hg19):

```
$ time awk '/^>/ && !a[$0]++{k++}; END{print k}' reference.fasta
93

real    0m34.510s
user    0m31.861s
sys     0m2.634s

$ time grep '^>' reference.fasta | sort -u | wc -l
93

real    0m16.597s
user    0m0.820s
sys     0m3.725s

```

---
Tags: fasta, linux
---


# Question
Title: Seurat Merged objects tSNE - How to paint on original IDs?
I am working with single-cell RNA-seq data, using the R package "Seurat" to cluster and visual data-points.

I had two single cell datasets from which I generated two Seurat objects. I then combined the two using MergeSeurat. I did differential gene expression analysis, performed clustering, and ran a tSNE plot.

Now what I want to do is create that same tSNE plot, but color each cell according to the original object from which it came. How do I do this?

# Answer
> 3 votes
Below a few lines of code that accompany BC Wang's answer. After using MergeSeurat the sample name will be added to meta data under orig.ident. this can then be used to color the tSNE either using group.by or pt.shape. The former will show colors for each sample, the latter will color each cluster and give sample id another shape.

```
path1 <- file.path("path to 10x sample 1")
path2 <- file.path("path to 10x sample 2")

S1 <- Read10X(data.dir = path1)
S2 <- Read10X(data.dir = path2)
S1 <- CreateSeuratObject(raw.data = S1, min.cells = 3, project = "S1")
S2 <- CreateSeuratObject(raw.data = S2, min.cells = 3, project = "S2")

s.combined <- MergeSeurat(object1 = S1, object2 = S2, add.cell.id1 = "S1",
            add.cell.id2 = "S2", project = "2samples")

```

All data analysis goes here: RunTSNE etc...

## Plot the data

```
TSNEPlot(s.combined, group.by ="orig.ident")
TSNEPlot(s.combined, pt.shape="orig.ident")

```

# Answer
> 1 votes
If I understand your question correctly, you are trying to present the tSNE plot coloring and labeling the cells representations based on the original sample ID, i.e. which Seurat object it was from.

I think they offered a vignette on that. Please check this link. https://satijalab.org/seurat/merge\_vignette.html

Basically, you can add the sample name for the cells in the argument option in `MergeSeurat` function. Please check the documentation on `MergeSeurat` function as well.

When I did this, I have a lot of samples to be merged, so I simply saved the sample info for every cell in each sample to different objects. You can also use `AddMetadata` to do it as well.

Hope that helps.

# Answer
> 0 votes
It is impossible to create the same tSNE plot without knowing which seed you used. tSNE plots are not like PCAs you don't get the same results when using the same data.

But if you want another tSNE plot with a color differentiating from which object each sample comes you'll need a column with these information and color the plot by it.

---
Tags: r, rna-seq, visualization, seurat
---


# Question
Title: Tool to remove a PCR contamination in NGS data
### BACKGROUND

I work on NGS data (illumina paired ends reads) coming from a full extract of RNA (metagenomic). We are interested in the viral fraction of this extract.

I observed a contamination with a PCR amplicon. This contamination has been identified as coming from a PCR result that we used to make near the library preparation. So I know exactly the primers used in the PCR.

Here you can see a tablet visualization illustrating the contamination (2 contaminations). I obtained that image with a mapping on a known reference of my virus.

### PROBLEM

The problem I have is that these PCR products are exactly the type of sequence I want from my extract (identifying part of the virus). So I have some "false positive" and some misidentification of viruses present in the extract.

I have thought about removing all reads containing primer sequences but it leads to a gap in my data and I can't identify the viruses in my extract anymore. I am thinking that all of the read pairs that begin before or end after my primers are "true-positives" and I would really like to keep them (even having the primer sequence in them).

**Do you know any tool that can help me in that situation?** i.e. removing reads containing a sequence but only if it begins with it.

*If the tool can also handle the degeneration of my primer it would be perfect.*
An example of primer I used is

```
CNTGGGAGGGCGATCGCAA

```

**Complementary information (asked in comments)**
I have around 30 million read pairs (60 millions reads in total). Contamination reads are between 200.000 to a 1 millions reads (depending on sample). "true" reads (not from the contamination) represent around 0.5% of the contamination.

# Answer
This bash script will get you all of the reads that do not start with your primer sequence. It also handles degeneracy. It requires the tool bioawk.

It works by looking at the first bases in the read and only prints out the headers that do not have a match with the degenerate primer. Then, `comm` makes sure that the read was retained in both the R1 and R2 filtering. This allows reads to be filtered whether or not the primers start at the beginning of the R1 or R2 read.

`seqtk subseq` is where you recover the reads that don't start with the primer.

Be warned though, any off-by-one or sequencing errors will allow some of the sequences with primers at the beginning. Also, this script requires that the headers for the R1 and R2 reads are the same and aren't appended with anything like "\1", "\2", et cetera.

```
#!/bin/bash
get_primer_free_starts () {
    bioawk -cfastx 'BEGIN{seqs[0] = "CATGGGAGGGCGATCGCAA"; seqs[1] = "CCTGGGAGGGCGATCGCAA"; \
        seqs[2] = "CGTGGGAGGGCGATCGCAA"; seqs[3] = "CTTGGGAGGGCGATCGCAA"} \
        {thisseq = substr($seq, 1, length(seqs[0])); printme=1\
        for (i=0; i<=3; i++) { \
            if (thisseq != seqs[i]){ \
                printme=0 \
            }
        }; if (printme==1){ \
            print($seq)
        }}' "${INPUT} | sort > "${OUTPUT}
}

INPUT="R1.fastq.gz"
OUTPUT="R1_OK.txt"
get_primer_free_starts

INPUT="R2.fastq.gz"
OUTPUT="R2_OK.txt"
get_primer_free_starts

comm -12 R1_OK.txt R2_OK.txt | sort | uniq > passing_reads.txt

seqtk subseq R1.fastq.gz passing_reads.txt > R1.filtered.fq
seqtk subseq R2.fastq.gz passing_reads.txt > R2.filtered.fq

```
> 2 votes

---
Tags: rna-seq, ngs, covid-19, pcr
---


# Question
Title: Why are there more barcodes than GEMs in 10X chromium data?
# Question: Why are there more barcodes than GEMs in 10X chromium data?

## Introduction

I have several 10X genomics chromium libraries made with the de novo assembly/genome protocol. The 10X chromium device is advertised to produce approximately 1,000,000 GEMs (gel-emulsion droplets). In theory each GEM should have a single barcode that is used to track the provenance of each read pair. If this is true, a single library should have approximately one million barcodes, whether or not it was overloaded or underloaded.

## Results

The libraries actually have more than one million barcodes, confirmed by `longranger basic`. Library 1 has 2,369,949 barcodes. Library 2 has 2,311,405 barcodes. Both libraries have approximately 150,000,000 read pairs.

## Why this is a problem

I'm using the 10X data for some non-standard analyses and want a high number of LinkedReads per molecule (LPM). Having multiple barcodes in a GEM will reduce the number of LPMs per barcode proportional to how many barcodes were present in the GEM.

## Possible explanations

* If there are more than approximately one million barcodes in a library then this leaves several possibilities:
    1. Each GEM truly has a single barcode. Therefore the Chromium device actually produced significantly more GEMs than advertised.
    2. Each GEM has a variable number of barcodes. Therefore the number of barcodes is independent of how many GEMs are created.
    3. Some combination of the above two options. Some variable number of GEMs and a variable number of GEMs per barcode.
    4. Additional barcodes are created by sequencing errors and misclassification by `longranger basic`. This would probably manifest in many barcodes that only have one or two linked reads total.

## Deeper analysis

I don't know how to test hypotheses 1, 2, & 3 bioinformatically. Hypothesis 4 is testable though by making a histogram of read counts per barcode. To do so, we will use the output of `longranger basic` for Lib2.

Longranger basic outputs an interleaved `fastq.gz` file that contains the whitelisted barcode in the header. For example (pretending the reads are shorter than 150bp):

```
@J00113:232:HFWCKBBXX:4:2224:17756:14889 BX:Z:AAACACCAGACAATAC-1
AATTACCCAAAATTAATATTTGTATTAATTTTCATATTAATAAC
+
JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ
@J00113:232:HFWCKBBXX:4:2224:17756:14889 BX:Z:AAACACCAGACAATAC-1
CACAATATCAAAAACAAATCAAAAATTAAGAAAATATAACTAAG
+
AAF<7FF-JJFJJFJJJJJJAJFJJJFJJJJJFJJJJJJJJJJJ

```

The `BX:Z:AAACACCAGACAATAC-1` tells us that the index is whitelisted (`BX:Z`), the index sequence (`AAACACCAGACAATAC`), and the library index for the analysis (`-1`).

So, we can count the number of occurrences of each barcode with this code:

```
bioawk -cfastx '{print($comment)}' barcoded.fastq.gz | sort | uniq -c > bc_counts.txt

```

The first column of `bc_counts.txt` tells us how many individual reads had that barcode. We divide this number by two to see how many read pairs had that barcode and make a histogram.

```
awk '{print($1)/2}' barcode_counts.txt | sort | uniq -c | sort -k2 -n | awk '{printf("%d\t%d\n", $2, $1)}' > pairs_per_bc.txt

```

The first column of this file has the number of read pairs in the barcode, and the second column has how many barcodes have that number of read pairs. We see that there are many barcodes that only have one or two read pairs:

```
1       593583
2       207192
3       79409
4       34174
5       17353
6       10662
7       7997

```

It looks like every barcode with fewer than 10 read pairs may be noise. How many barcodes contain at least 10 reads?

```
>awk '{if ($1 >= 10){sum += $2} } END{print(sum)}' pairs_per_bc.txt
1347686

```

Around one million.

## Conclusion

This analysis doesn't give conclusive evidence for the four hypotheses above, but it seems like there is a noise-generating process in 10X data that causes there to be many barcodes with only a few readpairs.

# Answer
> 3 votes
Your initial guess is almost certainly correct. I don't know about the linked read libraries, but in the 10X single cell sequencing protocol, separating real barcodes from noise barcodes is an important and sensitive step in the analysis pipeline.

The current analysis pipeline is to look at a plot, like the one you produced and find the "knees". The plot you produced has two knees, on at around 10 and one at around 200. I'd guess that the ones below 10 are what we call "error barcodes" generated by sequencing is PCR errors from real barcodes, and the ones between 10 and 200 are what we'd call "background" barcodes - in scRNAseq they represent things like GEMs that didn't encapsulate a cell, but did encapsulate some non-cellular DNA or RNA.

As there are more than 1,000,000 possible barcodes, not all of them are used. I believe 10X knows which sequences are used and which are not (ie the barcodes are not random), so I'm surprised that these sequences make it out of the other end of the analysis.

---
Tags: ngs, illumina, 10x-genomics, barcode
---


# Question
Title: IonTorrent variant caller failing, lack of ZM: tags in BAM files
I have a lot of IonTorrent sequencing data in BAM format that I want to call variants on, using the IonTorrent variant caller. However, for some of the BAM files (for about 12 samples out of ~150), the variant caller stops due to what it states is lack of ZM: tags in a few individual reads.

The thing is, most of the reads in the BAM file do not have ZM: tags (when viewed through samtools view), but these do not apparently cause problems.

What settings do I need to change to get these BAM files to process, or do I just filter out the problematic reads and analyse the remainder?

# Answer
It turns out, after speaking to tech support, that these are old files from a previous version of the software so won't run through variant caller.

The files themselves need reanalysing from the initial FASTQ -\> BAM with the newest version.
> 1 votes

---
Tags: variant-calling
---


# Question
Title: convert tRNA coordinates Result to Fasta format
I believe many of the experts here have gone through this issue,

I have run tRNAscan-SE on my genome to predict all the tRNAs, and have successfully gotten all the tRNA coordinates like start, end and type all arranged in a table.

My question is how can I convert this result to fasta format as long I have all the details about the tRNA ready?

Is there any tool that I can use to get the fasta sequences of the regions described in my table from the original genome?

---

```
Sequence                tRNA    Bounds  tRNA    Anti    Intron Bounds   Inf
Name            tRNA #  Begin   End     Type    Codon   Begin   End     Score   Note
--------        ------  -----   ------  ----    -----   -----   ----    ------  ------
NW_011590950.1  1       36704   36632   Ala     AGC     0       0       40.8    pseudo

```

# Answer
You can select “Output BED format” on the tRNAscan-SE server:

Given the output BED format and your input reference in FASTA format, you can easily generate a FASTA file of tRNA sequences via `bedtools getfasta`:

```
bedtools getfasta -name -s -fo ‹output.fasta› -fi ‹reference.fasta› -bed ‹tRNAscan-SE-output.bed›

```
> 3 votes

# Answer
You can do this using the tools from the `exonerate` package. I assume you have all the sequences shown in your table in a single, multifasta file. If so, first create a single file per sequence using `fastaexplode`:

```
fastaexplode reference.fa

```

This will create one file for each sequence in `reference.fa`. Now, you can read through the results table, and extract the sequences:

```
tmpFile=$(mktemp)
awk -F"\t" '{
                if(NR>3){ printf "%s %s ", $1, $2;
                    if($4>$3){
                        print $3, $4-$3
                    }
                    else{
                        print $4,$3-$4,"-"
                    }
                }
            }' file | while read name num start len rev; do
                fastasubseq "$name".fa $start $len > $num.fa
                if [[ ! -z $rev ]]; then
                    fastarevcomp $num.fa > $tmpFile &&
                    mv $tmpFile $num.fa
                fi
                done

```

I would suggest you go with Konrad's approach, however, that looks much cleaner and simpler.
> 1 votes

---
Tags: fasta, trnascan-se
---


# Question
Title: How do I rewrite a read group using pysam?
I am trying to rewrite a SAM/BAM file with altered read group entries using pysam. In this simplified version, I want to take a BAM, and rewrite the SM tags in all read groups to the same string and copy the alignments with new header into a new BAM. The contents of the read group dictionaries I am changing appear to have no effect on the header that prints or is written to file output. How can I alter the read groups?

```
import argparse
import pysam

def rewrite_read_groups(input_filename, output_filename, sample):
    try:
        bam_input = pysam.AlignmentFile(input_filename, 'rb')
        # replace header fields
        output_header = bam_input.header
        read_groups = output_header['RG']
        print(read_groups)
        for read_group in read_groups:
            read_group['SM'] = sample
            print(read_group)
        output_header.to_dict()
        print(output_header)
    finally:
        bam_input.close()

    print(output_header)

    try:
        bam_input = pysam.AlignmentFile(input_filename, 'rb')
        bam_output = pysam.AlignmentFile(output_filename, 'wb', header=output_header)

        # copy alignments
        for read in bam_input.fetch(until_eof=True):
            bam_output.write(read)

    finally:
        bam_input.close()
        bam_output.close()

if __name__ == "__main__":
    parser = argparse.ArgumentParser(description="Rewrite bam read groups")
    parser.add_argument('-s', "--sample", help="Set sample for all read groups")

    parser.add_argument('-i', "--input_file", help="Input bam", required=True)
    parser.add_argument('-o', "--output_file", help="Output bam", required=True)

    args = parser.parse_args()

    rewrite_read_groups(args.input_file, args.output_file, args.sample)

```

# Answer
> 2 votes
If the goal is to assign all alignments to a single read group (which this code seems to do), then Picard AddOrReplaceReadGroups might help: https://broadinstitute.github.io/picard/command-line-overview.html#AddOrReplaceReadGroups.

Modifying headers using PySAM might be tricky – the internal representation of headers has recently been changed and seems not to be stable across PySAM versions.

---
Tags: sam, pysam
---


# Question
Title: Tool to show DNA sequence and allowing upload of own graph data
**Background**

We want to be able to load (or request) data for a genome including the sequence and gene annotation (bacteria). Then, we want to load our own annotation which should be displayed as a line plot:

```
position score
1        5
2        6
etc.

```

Where each position corresponds to the position in the genome (i.e. FASTA). Something like this:

**Question**

Is there a tool that can load a genome sequence, gene annotation, and a position ~ score file, which can be displayed as a line plot?

# Answer
You can use UCSC genome browser, which supports custom genome sequence and annotation. Specifically, one needs to create a Track Hub to hold genomes and annotations.

Here is an example:

> http://genome.ucsc.edu/cgi-bin/hgTracks?hubUrl=http://zlab-trackhub.umassmed.edu/yfu/insects/hub.txt&db=hub\_129615\_cl\_v0.8.2

And here is the corresponding hub.txt file holding the metadata, and URLs to actual fasta, bigWig and bigBed files.
> 1 votes

---
Tags: genome, sequence-analysis, visualization
---


# Question
Title: How to generate Genbank format out of paired end fastq data?
I've unmapped cleaned paired end sequence data in fastq format of a bacterial genome. I want to get a sequence data in Genbank format in the end. What are the exact steps that I have to follow in galaxy webserver and if in some point necessary, also using linux based apps? I would appreciate a lot of details.

# Answer
The steps needed to assemble a genome in Galaxy are layed out in the Assembly training material from the Galaxy Training network. The software for this should be available on the main galaxy servers (e.g., usegalaxy.eu). Since you have a prokaryote, I suggest you follow the "Unicycler Assembly" workflow, which is also available on usegalaxy.eu.
> 1 votes

# Answer
I think the reason for your confusion is vocabulary.

The keywords you search for are genome assembly (see @Devon Ryan's answer) a process that is trying to stick together reads to recreate the sequenced genomic sequence. The assembled sequence needs to be annotated if you are interested where are the genes or what they could possibly do. If you will search for tutorials using these keywords you will find plenty.

I am not sure why you need it, but Genbank format contains the genome (the product of genome assembly), metadata (who is the submitter, what species it is, etc) and genome annotation (what is in the region of genome).
> 0 votes

---
Tags: sequence-annotation, genome
---


title: Generating basic dna sequences in R

text:
# Question
Title: Generating basic dna sequences in R
Looking for some guidance.

In an introductory bioinformatics class I am tasked with doing the following in R:

1. Create the sequence ACTGACTG
2. Generate a random string of nucleotides with equal probabilities.

Are there certain functions that I should be looking at? No library or package or set of functions was provided so I am completely lost. I am assuming there must be functions in R that do these things that I need to manipulate, but I cannot seem to find any resources explaining that. I looked at seq() (thinking it might generate DNA sequences), but it only seems to generate sequences of integers. I also looked at randDNA(), but it is random, how could I make it generate the exact string in question 1).

Any direction will be greatly appreciated. Thanks.

# Answer
> 1 votes
Without giving away too much for your homework check out the `sample()` function in R. Think about how you can pass DNA bases to this function to generate random sequences.

---
Tags: genome-sequencing, homework
---


# Question
Title: tRNAscan-SE error: FATAL: Unable to find /usr/local/bin/cmsearch executable
I have downloaded tRNAscan-SE from here. After decompressing and untaring the file, I installed it using:

```
./configure
make
make install

```

When I type `tRNAscan-SE --help`, I get the help page:

```
tRNAscan-SE 2.0 (December 2017)

FATAL: No sequence file(s) specified.

Usage: tRNAscan-SE [-options] <FASTA file(s)>

  Scan a sequence file for tRNAs
          -- default: use Infernal & tRNA covariance models
             with eukaryotic sequences
             (use -B, -A, -M, -O or -G to scan other types of sequences)

Basic Options
  -E         : search for eukaryotic tRNAs (default)
  -B         : search for bacterial tRNAs
  -A         : search for archaeal tRNAs
  -M <model> : search for mitochondrial tRNAs
                 options: mammal, vert
  -O         : search for other organellar tRNAs
  -G         : use general tRNA model (cytoslic tRNAs from all 3 domains included)
  -L         : search using the legacy method (tRNAscan, EufindtRNA, and COVE)
                 use with -E, -B, -A, -O, or -G
  -I         : search using Infernal (default)
                 use with -E, -B, -A, -O, or -G
  -o <file>  : save final results in <file>
  -f <file>  : save tRNA secondary structures to <file>
  -m <file>  : save statistics summary for run in <file>
               (speed, # tRNAs found in each part of search, etc)
  -H         : show both primary and secondary structure components to
               covariance model bit scores
  -q         : quiet mode (credits & run option selections suppressed)

  -h         : print full list (long) of available options

```

I then tried to use it with `tRNAscan-SE c_elegans.PRJNA13758.WS263.genomic.fa` but it gave the following error:

`FATAL: Unable to find /usr/local/bin/cmsearch executable`

After googling cmsearch, it seems that is is part of the Infernal package. So, I installed infernal with conda using `conda install infernal`.

I pulled the help from `cmsearch` successfully:

```
cmsearch -h
# cmsearch :: search CM(s) against a sequence database
# INFERNAL 1.1.2 (July 2016)
# Copyright (C) 2016 Howard Hughes Medical Institute.
# Freely distributed under a BSD open source license.
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Usage: cmsearch [options] <cmfile> <seqdb>

Basic options:
  -h        : show brief help on version and usage
  -g        : configure CM for glocal alignment [default: local]
  -Z <x>    : set search space size in *Mb* to <x> for E-value calculations  (x>0)
  --devhelp : show list of otherwise hidden developer/expert options

Options directing output:
  -o <f>       : direct output to file <f>, not stdout
  -A <f>       : save multiple alignment of all significant hits to file <s>
  --tblout <f> : save parseable table of hits to file <s>
  --acc        : prefer accessions over names in output
  --noali      : don't output alignments, so output is smaller
  --notextw    : unlimit ASCII text output line width
  --textw <n>  : set max width of ASCII text output lines  [120]  (n>=120)
  --verbose    : report extra information; mainly useful for debugging

Options controlling reporting thresholds:
  -E <x> : report sequences <= this E-value threshold in output  [10.0]  (x>0)
  -T <x> : report sequences >= this score threshold in output

Options controlling inclusion (significance) thresholds:
  --incE <x> : consider sequences <= this E-value threshold as significant  [0.01]
  --incT <x> : consider sequences >= this score threshold as significant

Options controlling model-specific reporting thresholds:
  --cut_ga : use CM's GA gathering cutoffs as reporting thresholds
  --cut_nc : use CM's NC noise cutoffs as reporting thresholds
  --cut_tc : use CM's TC trusted cutoffs as reporting thresholds

Options controlling acceleration heuristics*:
  --max      : turn all heuristic filters off (slow)
  --nohmm    : skip all HMM filter stages, use only CM (slow)
  --mid      : skip first two HMM filter stages (SSV & Vit)
  --default  : default: run search space size-dependent pipeline  [default]
  --rfam     : set heuristic filters at Rfam-level (fast)
  --hmmonly  : use HMM only, don't use a CM at all
  --FZ <x>   : set filters to defaults used for a search space of size <x> Mb
  --Fmid <x> : with --mid, set P-value threshold for HMM stages to <x>  [0.02]

Other options*:
  --notrunc     : do not allow truncated hits at sequence termini
  --anytrunc    : allow full and truncated hits anywhere within sequences
  --nonull3     : turn off the NULL3 post hoc additional null model
  --mxsize <x>  : set max allowed alnment mx size to <x> Mb [df: autodetermined]
  --smxsize <x> : set max allowed size of search DP matrices to <x> Mb  [128.]
  --cyk         : use scanning CM CYK algorithm, not Inside in final stage
  --acyk        : align hits with CYK, not optimal accuracy
  --wcx <x>     : set W (expected max hit len) as <x> * cm->clen (model len)
  --toponly     : only search the top strand
  --bottomonly  : only search the bottom strand
  --tformat <s> : assert target <seqdb> is in format <s>: no autodetection
  --cpu <n>     : number of parallel CPU workers to use for multithreads

*Use --devhelp to show additional expert options.

```

I tried the tRNAscan-SE command again but it failed with the same error message.

So, tRNAscan-SE is looking for cmsearch in the wrong place. How can I make tRNAscan-SE use cmsearch successfully?

Cheers.

### Update

Forgot to tell you I am working on Mac OS X and that cmsearch is installed here:

`/Users/user/anaconda/envs/python3env/bin/cmsearch`

and my `python3env` anaconda environment is active.

# Answer
> 5 votes
An easy way to build tRNAscan and their dependencies "all Infernal packages" without confusion is to run

```
conda install -c bioconda trnascan-se

```

# Answer
> 8 votes
You may try to set the correct path through the `infernal_dir` variable in `tRNAscan-SE.conf.src`. You might need to do this before the installation process.

Original file starts with:

```
# tRNAscan-SE 2.0
# Configuration File

# default paths
bin_dir: @bindir@
lib_dir: @libdir@/tRNAscan-SE
infernal_dir: {bin_dir}

```

Try the following:

```
# tRNAscan-SE 2.0
# Configuration File

# default paths
bin_dir: @bindir@
lib_dir: @libdir@/tRNAscan-SE
infernal_dir: "/Users/user/anaconda/envs/python3env/bin"

```

# Answer
> 4 votes
tRNAscan-SE 2.0 requires separate installation of Infernal 1.1.2. The source code and prebuilt binaries can be obtained at http://eddylab.org/infernal/.

The installation directory of Infernal should be the same as the one where tRNAscan-SE 2.0 is installed.

When install infernal, run the configure command with the option `--prefix <your tRNAscan-SE-2.0 install path>`. If you didn't add option `--prefix` when installing tRNAscan-SE-2.0, then you should install infernal without option `--prefix`.

# Answer
> 2 votes
**A much cleaner solution is given in @bli's answer.**

---

The error message is telling you that tRNAScan-SE is looking for the `cmsearch` executable at `/usr/local/bin/cmsearch`. However, you installed it though conda which apparently placed it at `/Users/user/anaconda/envs/python3env/bin/cmsearch`.

A simple, if not very elegant solution, would be to make a link to it:

```
sudo ln -s /Users/user/anaconda/envs/python3env/bin/cmsearch /usr/local/bin/

```

---
Tags: sequence-annotation, software-installation, trna, trnascan-se
---


# Question
Title: Store contigs and genome substrings based on mismatches
I have a string, a list of substrings, and a list of positions where substrings are found in string.

```
string = 'ATGGTACACGCTACGAGCTAG'
substrings = ['TACA', 'CTAGAAAAG']
positions = ['11, 4', '']

```

Substring `'TACA'` appears at positions 11 and 4 in string (i.e. `'11, 4'`).

Substring `'CTAGAAAAG'` does not appear in string (i.e. `''`).

I need to remove `'CTAGAAAAG'` from `substrings` and I need to remove `''` from `positions`, while creating a new list with `'CTAGAAAAG'`. **I need your help here to create code which will give the following results:**

```
new_substrings = ['TACA']
new_positions = ['11, 4']
new_list = ['CTAGAAAAG']
string = 'ATGGTACACGCTACGAGCTAG' (obtained from above)

```

Once this is complete, the **given** code:

```
new_positions = [map(int, x.split(',')) for x in new_positions]
seq_pos = dict(zip(new_substrings,new_positions))
seq_ref = {k: [string[start:start+length] for start in starts] for k, starts, length in zip(new_substrings,new_positions,[len(x) for x in new_substrings])}

```

Gives:

```
seq_pos = {'TACA': [11, 4]}
seq_ref = {'TACA': ['TACG', 'TACA']}

```

Then, from seq\_ref, I need to create a list of integers which represents the number of mismatches between the key: `'TACA'` and both of its values: `'TACG'` and `'TACA'`.

```
mismatches = [[1, 0]]

```

The code to find the # of mismatches between two strings is provided:

```
def HammingDistance(p, q):
    dist = 0
    for i in range(len(p)):
        if p[i] != q[i]:
            dist += 1
    return dist

```

Finally, I need to save my data in a CSV file (i.e. Excel or Numbers).

```
new_positions = ['11, 4'] (obtained above)
new_substrings = ['TACA'] (obtained above)
actual_substrings = [['TACG', 'TACA']] (values of seq_ref)
mismatches = [[1, 0]] (obtained above)

```

Info needs to be sorted based on # of mismatches.

I understand this is a multi-step problem but any help would be appreciated.

# Answer
> 1 votes
For the first question try something like:

```
new_substrings = []
new_positions = []
new_list = []
for pos, sub in zip(positions, substrings):
   if not pos:
       new_list.append(sub)
   else:
       new_substrings.append(sub)
       new_positions.append(pos)

```

this gives:

```
>>> new_substrings
['TACA']
>>> new_positions
['11, 4']
>>> new_list
['CTAGAAAAG']

```

For the hamming distance you just need to iterate over the values in `seq_ref` and append the results of the hamming distance function in `mismatches`. Based on the given code and the example above you should be able to create that block. Below is an untested quick example to get you on track.

```
mismatches = []
for seq in seq_ref.keys():
   for value in seq_ref[seq]:
       mismatches.append(HammingDistance(seq, value))

```

As far as the CSV and sorting is concerned I suggest that you look into pandas or numpy. Both of those packages provided functions to help with this.

---
Tags: python
---


# Question
Title: Machine learning using protein-sequences
I'm participating in a bioinformatics machine-learning seminar at my university. The main task is predicting binary classification of protein-protein interactions using sequence data as input.

One of the subtasks is familiarization with the dataset and presenting the dataset. Now I'm wondering which additional information I can get out of the sequence data.

I want to start out with binary classification for the protein interaction with the labels “Interact” and “Non-Interact. The dataset that I got provided consists of two fasta files. One containing the sequence header with the species and the corresponding amino-acid sequences ~300 entries. The other containing the exact same species header and the label “Interact” or "Non-interact”. Species are completely mixed.

Is there something interesting that I can additionally include in my analysis? I'm used to exploratory data analysis though this is different due to the nature of the protein sequences.

# Answer
There are two types of features which can be extracted from protein sequences. It is not possible to known which, if any, of these features would be useful in classification tasks. It may not be possible to build a classifier at all, or it may be very straight forward. To know this, a feature selection technique must be used, such as e.g. forward feature selection.

**Sequence based features**

Sequence based features are derived from sequence-only information.

The most common used features are the PseAAC (pseudo amino acid composition) style features, including:

* Average Position of each amino acid - 20 features
* Average Distance between each amino acid - 20 features
* Motifs - 20 features, 20\*20 features, and/or 20\*20\*20 features (depending on motif length)
* OAAC (Overall Amino Acid Composition) - 20 features
* SAAC (Split Amino Acid Composition) - 60+ features
* AAIndex Composition - 1 to 567+ features
* PCC (Physicochemical Composition) - ~13 features
* PSSM encoded as 1, 20, 200, 400 or 600 features (i.e. Psi-Blast output using NR or swissprot)
* Average predicted residue burial/exposure level
* Average predicted sequence entropy

**Structure based features (predicted)**

Using the sequence, you can either predict the 3d structure using homology, ab initio, or directly use a homologous structure. You can then extract predicted 3d features, such as:

* 3d size (i.e. max of x,y,z)
* Calculated ASA/SASA
* 3d contacts (type, frequency, etc)
* any other local/global environmental information
* cath/scop classification
* dssp/stride composition
> 8 votes

---
Tags: machine-learning, protein-protein-interaction
---


title: What is a site in a phylogenetic tree?

text:
# Question
Title: What is a site in a phylogenetic tree?
When I read about a term related to phylogenetics, for example, $K\_a$/$K\_s$ ratio, I constantly see the word *site*, whose meaning isn't explained and it puzzles me.

For example substitutions per synonymous site, substitutions per non-synonymous site, substitutions at a site etc.

Could you define a phylogenetic *site*?

So there were 0.6 substitutions per the blue-violet site, what does it mean? What is the length of the site? Is the length of sites constant? Is it the length of a sequence read? Are introns the part of the site?

# Answer
A site is simply an individual discrete position, normally a single nucleotide (or amino acid). For example, consider a toy alignment:

```
sequence1 ATGC
sequence2 ATCC
#           ^ arrow points to the third site

```

Here there is 1 mutation in 4 sites, so 0.25 substitutions per nucleotide site.

I'm not sure where the term originated, but probably one of the earlier definitions is from renowned population geneticist Motoo Kimura, who simply stated in 1968:

> In this paper, “site” refers to a single nucleotide pair

A pair in this context, means either a homozygous or heterozygous site.

Synonymous sites are positions where nucleotide changes lead to no amino acid change, these are typically at the 3rd position of the codon where the nucleotide is often degenerate, for example, `TAT`, `TAC` `TAG`, and `TAA` all encode the Serine amino acid. By contrast non-synonymous sites are where a single nucleotide change does lead to an amino acid change, and therefore possibly a more important functional change to the resulting protein.

Thus, the ratio of Ka/Ks (or Dn/Ds) is thus often used as a crude proxy of selection, with Ka/Ks \< 0, Ka/Ks = 0 and Ka/Ks \> 0 implying purifying (negative) selection, neutral evolution and positive selection respectively. Note that this assumes the neutrality of synonymous sites, which is increasingly questioned due codon usage bias and regulatory elements overlapping coding regions (and also ignores mutational biases and biased gene conversion, although these can be accounted for in an evolutionary model, and also assumes no sequencing errors and perfect alignments, which can impact a lot on the results too).
> 4 votes

---
Tags: phylogenetics, phylogeny, terminology
---


title: How to calculate AUC in coverage graph

text:
# Question
Title: How to calculate AUC in coverage graph
Is there a way to calculate the area under a curve in a coverage graph?

Thank you in advance.

EDIT:

I'd like to calculate the AUC per strand per gene in a coverage graph. I have two wig files (sense and antisense) generated from a single bam (RNA-seq) file (using the strand-coverage software).

# Answer
> 4 votes
Yes, just add up the coverage at each location; that's integration in a nutshell.

Or, alternatively / equivalently, count the total number of reads within the region.

---
Tags: rna-seq, coverage
---


# Question
Title: Reference genome for allele specific expression
We are trying to sort out a pipeline for doing allele specific expression. Our plan is to call SNPs from RNA-seq data and combine with known SNP annotations. A well known problem in ASE is reference bias, where reads are more like to map if they have the reference sequence than the alt sequence.

1. We have come accross two possible ways to deal with this. In many papers (e.g. Zhuo et al, Sigurdsson et al) a new reference is generated, substitute N bases where there are SNPs using FastaAlternativeReferenceMaker from GATK and reads are remapped to that.
2. HISAT2 offers the possibility to provide SNP calls to the index builder, and build an index that accounts for these SNPs.

Does anyone have any advice on which of these might be better?

# Answer
> 3 votes
As you said, this is a well-known problem for ASE. In Degner et al. 2009, they state that, "Perhaps surprisingly, we found that masking known SNPs does little to eliminate inherent biases in read-mapping."

As such, substituting N bases is not likely to completely solve your issue. Depending on the aligner, these N bases may also negatively affect your ability to map reads.

Of your two suggestions, I believe HISAT2 will be the better approach. There are also more intensive methods you can use such as empirical Bayesian methods or more traditional variant phasing (ex. AlleleSeq).

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


# Question
Title: Does read sequence identifier need to be unique for proper alignment?
I am preprocessing a set of scTHSSeq reads in which each read sequence identifier has been replaced with the cell barcode, e.g.:

```
@CGAGGCTG_TATCCTCT_AGTCCA
...
...
...
@CGAGGCTG_TATCCTCT_GGTACG
...

```

The read sequence identifier has been replaced with the cell barcode to be able to demultiplex alignments after bwa into a different bamfile for each cell.

I ran the alignment using bwa mem and I received no errors. However I would like to know if there are caveats or reasons why this approach may lead to misleading results.

I will then proceed to remove PCR duplicates with `samtools rmdup` executed in each cell's bamfile and then I'd run peak calling using spp.

# Answer
What tool did you use to replace the read names? Can you modify that to split up the reads into separate files or encode the cell barcode in the read description instead of the read name?

Typically is it assumed by many tools that the FASTQ read names are unique. you could have issues with BWA and split alignments but I don't know if that pertains to your data. I think the best you can do is the sanity check the results at each step of the process.

In the fastq format (https://en.wikipedia.org/wiki/FASTQ\_format) on the first line denoted by the '@' there can be two values separated by a space. The first is the name and the second is the description.

So an example would be:

```
@Original_read_name CGAGGCTG_TATCCTCT_AGTCCA
ATCGATCATGATGACA
+
@&%*@(%@%*($!##!

```
> 2 votes

---
Tags: read-mapping, reads, identifiers
---


# Question
Title: How to calculate Gene Ontology terms in python
I am testing different `Protein-Protein Interaction networks` alignment tools. I observed that all alignment tools show some `Gene Ontology` terms similarity. How can I calculate `Functional similarity` / `average functional similarity` for `Gene ontology`'s terms `Biological processes`, `Cellular components`, `Molecular functionality` of any two proteins networks

# Answer
> 2 votes
You can calculate gene ontology similarity with the GOSemSim R package (paper)(disclaimer I'm a contributor to this package). You have several similarity scores implemented, some of them are across the three subontologies.

In python you have the GOATOOLS (paper). I think there was another Python package but now I don't remember its name.

---
Tags: sequence-alignment, proteins, networks, protein-protein-interaction
---


# Question
Title: What are "split reads" and "intron clusters?"
I'm working through the example data set for LeafCutter and the documentation mentions "split reads":

> This will cluster together the introns fond \[*sic*\] in the junc files listed in test\_juncfiles.txt, requiring 50 split reads supporting each cluster and allowing introns of up to 500kb.

Searching "split reads" on Google yielded few results, and the results I did find were hard for me understand,

> Chimeric reads occur when one sequencing read aligns to two distinct portions of the genome with little or no overlap. Chimeric reads are indicative of structural variation. Chimeric reads are also called split reads.

One sequencing read aligns to two distinct portions of the genome, but then with little or no overlap? Isn't that a self-contradiction?

I can more or less deduce what "intron clusters" are but I figured I might as well try to solicit a more formal definition. I know what introns are and I can figure what it means to cluster them together. I guess I just want to know what split reads and intron clusters mean in the context of working out the example data for LeafCutter.

# Answer
**split reads** \- These are read that have two or more alignments to the reference from unique region of the read. In this example a 150bp read sequenced from RNA could have base 1-75 aligning to the 3' end of exon2 and bases 76-150 aligning to the 5' end of exon3. This would be a split read because it have two alignments (exon2 and exon3) and those alignment are from unique (non overlapping) parts of the read, 1-75 and 76-150.

Chimeric reads are an example of split reads. Many people use "chimeric reads" to describe a type of noise that occur during PCR or sequencing where the bases from two different molecules are reported as one read and the alignments of that read become split when it aligns back to the location of the two original molecules.

**"intron clusters"** \- In this case I think that the LeafCutter software is describing stacks of split reads to cluster to support the detection of an intron.
> 5 votes

---
Tags: sequence-alignment, read-mapping, rna-alignment, rna-splicing
---


# Question
Title: Reject reads with low quality bases from a Bam file through pysam
I have a code below:

```
def FilterReads(in_file, out_file):

def read_ok(read):
    """
    read_ok - reject reads with a low quality (<5) base call
    read - a PySam AlignedRead object
    returns: True if the read is ok
    """
    if any([ord(c)-33 < _BASE_QUAL_CUTOFF for c in list(read.qual)]):
        return False
    else:
        return True

_BASE_QUAL_CUTOFF = 30

bam_in = pysam.Samfile(in_file, 'rb')
bam_out = pysam.Samfile(out_file, 'wb', template=bam_in)

for read in bam_in.fetch():
    if read_ok(read):
        bam_out.write(read)

```

This code works fine by

* Rejecting reads having a base with phred quality score below five
* But it first takes a BAM file
* And then creates a filtered BAM file for further analysis but this takes a lot of time.

So is there a way to reject these reads using pileup in pysam so that I may not have to create a file then sort it and again read it? Or can I modify this code to perform the same function

# Answer
> 1 votes
For can do this by accessing the basecall qualities from `PileupRead.alignment`. For example:

```
_BASE_QUAL_CUTOFF = 30
with pysam.Samfile(in_file, 'rb') as bam_in:
    for column in bam_in.pileup():
        for read in column:
            qual = read.alignment.query_qualities
            if all([ord(c)-33 >= _BASE_QUAL_CUTOFF for c in qual]):
                # do some stuff

```

This involves checking every read at every position. You could memoise this task:

```
def check_quality(read):
    try:
        return passed_quality[read.query_name]
    except KeyError:
        qual = read.alignment.query_qualities
        passed = all([ord(c)-33 >= _BASE_QUAL_CUTOFF for c in qual])
        passed_quality[read.query_name] = passed
        return passed_quality[read.query_name]

_BASE_QUAL_CUTOFF = 30
passed_quality = dict()
with pysam.Samfile(in_file, 'rb') as bam_in:
    for column in bam_in.pileup():
        for read in column:
            if check_quality(read):
                # do some stuff

```

Or, using `functools.lru_cache` (thanks @winni2k):

```
from functools import lru_cache

@lru_cache(maxsize=1000)
def check_quality(read):
    qual = read.alignment.query_qualities
    return all([ord(c)-33 >= _BASE_QUAL_CUTOFF for c in qual])

_BASE_QUAL_CUTOFF = 30
passed_quality = dict()
with pysam.Samfile(in_file, 'rb') as bam_in:
    for column in bam_in.pileup():
        for read in column:
            if check_quality(read):
                # do some stuff

```

---
Tags: sam, ngs, sequence-alignment, python, pysam
---


# Question
Title: Order of batch effects removal, data imputation and library size normalization in scRNA-seq data
I am preprocessing scRNA-seq data. What is the best practice in use to run both ComBat for batch effects removal, data imputation (to mitigate dropout) and library size normalization?

I thought that library size should be run first, since it is per-cell normalization, then ComBat batch effects removal. On the original paper - Johnson et al. (2007) \- it is stated that:

> We assume that the data have been normalized and expression values have been estimated for all genes and samples.

However, I want to apply it to scRNA-seq data. Does this statement still hold? Additionally, I plan to apply imputation (e.g. with MAGIC) in the end. Is there any problem you can spot?

### Update

I attach the PCA regarding an example Mus Musculus dataset in which different colors represent different mice. It seems clear to me that the first two principal components are affected by batches (mouse id).

### Update 2

I rerun the PCA on raw counts data (the first PCA was on log-transformed data) and I obtain a different description of the dataset, in which batch effects seem not to be prevalent.

# Answer
> 3 votes
MAGIC assumes input data has been both library-size normalized, and either log or sqrt transformed prior to imputation (see also: MAGIC tutorial). Additionally, any graph-based methods (MAGIC, PHATE, t-SNE, UMAP, spectral clustering, Louvain, etc etc) will give flawed results if your data contains a batch effect, since the neighbourhood graph would reflect that structure of your batch effect, and worse, imputation would further reinforce this batch effect.

Thus I would recommend the following pipeline:

* Library-size normalization
* Square root (or log) transform
* Batch effect removal
* Imputation

Regarding your **update**, the reason you don't see the batch effect in the raw counts data is simply that the batch effect is not visible in the most highly expressed genes. Prior to transformation, the principal source of variation in your data is simply the expression of the most highly expressed genes - this is essentially masking the source of the batch effect, not removing it. I recommend never working with raw molecule counts in scRNAseq, as the raw counts data hides much of the heterogenity in your dataset, which is precisely what you are looking in doing single-cell RNA-seq.

---
Tags: normalization, scrnaseq, imputation, combat, batch-effects
---


# Question
Title: Fitting a principal curve into a diffusion map
I returned a Seurat object with a diffusion map by

```
 seurat_object <- RunDiffusion(seurat_object,genes.use = seurat_object@var.genes)

```

as below

Now, how I could know the order of a cell's projection onto this principal curve representing its predicted progression through the development by fitting a principal curve on my diffusion map?

I know Seurat slot for DM1 and DM2 is inside the seurat@dr

by

```
fit=principal_curve(as.matrix(seurat_object@dr$dm@cell.embeddings))
plot(fit)
line(fit)
points(fit)

```

I draw this graph of my cells

But how I know the first and last cells on curve? I mean how I know cells from each cluster is being placed where on the curve? By `fit$ord` I will obtain orders of cells on the line but how I know the first and end of the line?????

# Answer
If `fit$ord` gives you the order of cells on the curve, you should be able to obtain the first and last cells by doing the following:

```
first_cell_index <- which.min(fit$ord)
last_cell_index <- which.max(fit$ord)

```

However, I'm not sure this is so informative - if you're looking to get representative gene expression for the start and end of the curve, you should really be looking at the cluster of cells towards the start and end of the curve - for example, by comparing cells in clusters 1 and 3 that you showed in your first plot. A single cell will suffer from dropout and not tell you a lot about the cluster.
> 1 votes

---
Tags: r, scrnaseq, seurat
---


# Question
Title: Compare clusters from different datasets
I have two `scRNAseq` datasets, and I select a cluster from one of them, and another cluster from the other. I want to find differentially expressed genes between the two clusters. What would be the right way of doing that? Right now I am using `DESeq2` in the following way. I select only the selected clusters' cells, merge them together forming a matrix and then feed them into `DESeq2`:

```
mtx <- merge(data_1,data_2,by="row.names",all=TRUE)
mtx[is.na(mtx)] <- 0
mtx_conv <- matrix(as.numeric(unlist(mtx[,2:ncol(mtx)])),nrow=nrow(mtx))
mtx_int <- round(mtx_conv, 0)
rownames(mtx_int) <- as.character(mtx$Row.names)

cluster <- as.factor(c(rep(0, dim(data_1)[2]), rep(1, dim(data_2)[2])))
colData <- data.frame(
  Cluster = cluster
)

dds <- DESeqDataSetFromMatrix(countData = mtx_int,
                          colData = colData,
                          design = ~ Cluster)
dds <- DESeq(dds)

```

I am not sure though that it is the correct way to go because I am not using the full original datasets. The reason is that I want the computations to run faster. Another way of doing that is probably combining the whole two datasets, defining 3 clusters: two groups I am interested in, and all other cells, and then running the analysis. Would it differ significantly? Am I on the right track at all?

I have another idea: to use `Seurat` package. There, I would run `CCA` algorithm to align two full datasets, and then run `FindMarkers` function between the two clusters. I feel like it may be wrong, because the two datasets may need to be re-normalized together but Seurat does not seem to be doing that:

https://satijalab.org/seurat/Seurat\_AlignmentTutorial.html

Any suggestions would be greatly appreciated.

# Answer
I've been working with the multi-CCA alignment workflow in the Seurat package and am very impressed with functionality. As usual, the CCA tutorial is excellent for an alignment workflow and will get you off the ground quickly. You can even use the DESeq2 algorithm for determining differential expression from within Seurat.

When you use the function `RunCCA()` you can pass `rescale.groups=TRUE` if you wish to re-scale.

**Abstracted workflow**

1. Data preprocessing and gene selection
2. Define a shared correlation space with canonical correlation analysis
3. Remove dataset specific cells
4. Align correlated subspaces using dynamic time warping
5. Integrated analysis across datasets (clustering, trajectory building, etc ...)

---

Butler, A., and Satija, R. (2017). Integrated analysis of single cell transcriptomic data across conditions, technologies, and species.
> 1 votes

---
Tags: scrnaseq, differential-expression, deseq2, seurat
---


title: chromPlot R package: remove scales

text:
# Question
Title: chromPlot R package: remove scales
I am using chromPlot R package to plot some SNPs and their associated score, however, I cannot find a way to remove the scale of such score (circled in the attached picture). Does anybody know how to do it?

Thank you so much.

# Answer
> 1 votes
Setting the `legChrom` argument to `NA` will omit plotting a legend.

Found in the manual for chromPlot on Bioconductor for release 3.7.

---
Tags: snp, chromplot
---


# Question
Title: Is it possible to create a DESeqDataSet with a user-provided design matrix?
I'm trying to run a differential gene expression analysis using DESeq2, with counts coming from kallisto. I have imported them using tximport and I'm creating the `DESeqDataSet` (dds) using the `DESeqDataSetFromMatrix` function.

```
> dds <- DESeqDataSetFromMatrix(counts,
                                samples,
                                ~ strain + batch)

```

And I get the following error, expected given my experimental design:

```
Error in checkFullRank(modelMatrix) :
  the model matrix is not full rank, so the model cannot be fit as specified.
  One or more variables or interaction terms in the design formula are linear
  combinations of the others and must be removed.

```

Now, I know that I can just remove one column from the design matrix to make it work, but is there a way to supply my own design matrix to DESeq2? The following code raises an error:

```
> design <- model.matrix(~strain+batch, s2c)
> design = design[, -9] # removing a column to avoid full-rank
> dds <- DESeqDataSetFromMatrix(counts, s2c, design=design)
converting counts to integer mode
Error: $ operator is invalid for atomic vectors

```

Is there a way to provide my own `model.matrix`?

p.s. the modified model works in sleuth, but I would like to use DESeq2 for this particular analysis.

# Answer
> 4 votes
Provide rank sufficient design to `DESeqDataSetFromMatrix` and then use your custom model matrix in `DESeq`. In essence:

```
dds = DESeqDataSetFromMatric(counts, s2c, design=~batch)
design <- model.matrix(~strain+batch, s2c)
design = design[, -9]
DESeq(dds, full=design)

```

See this thread on the bioconductor site for details.

# Answer
> 1 votes
If your strain and batch are confounded, (like you did all of strain one on one day, and all of strain two on another), then there is no clever way to fix that. Just use one column, and inform whomever designed the experiment that you can't separate batch effect from strain difference.

---
Tags: rna-seq, deseq2
---


# Question
Title: Can I create a CRAM file with a relative reference path?
I’m trying to create a CRAM file that stores its path to the FASTA reference as a relative path, rather than an absolute path, so that I can move the files around. Unfortunately I can’t get this to work; I was expecting the following to work:

```
⟩⟩⟩ samtools view -C -T ../reference/ref.fa -o output.cram input.bam

```

However, the resulting file contains an absolute path in its header:

```
⟩⟩⟩ samtools view -H output.cram
…
@SQ     SN:1    LN:249250621    M5:hash     UR:/absolute/path/to/data/mapped/../reference/ref.fa
…

```

As a result, I am unable to open the file via a different path mount that results in different absolute paths, and I can’t move the file (+ its reference) around, or to different machines.

I know that I could set the `REF_PATH` environment variable or specify `-T` when reading the file but I would like to avoid this (the result file needs to be readable by IGV, launched by users who don’t know how to set environment variables).

Is there a way of creating a CRAM file that stores a *relative path* to its reference?

# Answer
> 1 votes
In your CRAM file, make all of the reference sequences present in the `./` directory. Then make a bash script to make relative links to your reference for whatever mount you are using at the time - ie `ln -s TARGET LINK_NAME`.

Now your CRAM is portable and you just have to tailor the local working dir to it.

---
Tags: sam, samtools, cram
---


# Question
Title: featureCounts segmentation fault on Arch Linux
I am encountering a segmentation fault when attempting to run featureCounts from `subread-1.6.3` on even small test data. I installed featureCounts from SourceForge using the precompiled Linux binary.

Most existing help threads on this topic (, , ) show some splash screen output before the error, but I get a segmentation fault before any error is shown.

```
$ featureCounts -T 1 -p -t exon -g gene_id -a /data/sai/annots/musmus_cdna_all_ensembl92.gtf -o test.txt test.bam test2.bam
Segmentation fault (core dumped)

```

# Answer
> 3 votes
I resolved this problem by installing from source; I assume the precompiled binary was not compatible with the Arch Linux kernel.

---
Tags: rna-seq, sam, subread, featurecounts
---


title: about dotplot legend meaning

text:
# Question
Title: about dotplot legend meaning
Here is the code I used:

```
DotPlot(object = cluster, genes.plot = c("Hcn4","Shox2","Nkx2-5","Nppa","Tnnt2","Gja5","Tbx18","Tbx3","Isl1","Scn5a","Cacna1d"),  plot.legend=TRUE)

```

I was confused about the legend with the Negative number shown for `avg.exp.scale`. Does it mean that the cluster expressed the genes at low level or that gene was down regulated? (like the Shox2 in the "0" compared with "2").

# Answer
The expression values for each gene are scaled / standardized by subtracting the genes mean expression and dividing by its standard deviation. A value of -1 would imply it's one standard deviation below the mean expression for that gene.
> 3 votes

---
Tags: r, seurat, single-cell
---


# Question
Title: What are some of the major questions pursued by the field bioinformatics?
I am an occasional user of tools in bioinformatics and I am fascinated by the possibility of uncovering deeper patterns in the biological realm using the help of computer algorithms. Yet I don't know much about the current state of this science. As I am considering a post-graduate degree in bioinformatics, I am interested to dig deeper into this subject and check whether the subject indeed aligns with my interest. So I start with the basic question stated in the title: What are some of the major questions pursued by the field bioinformatics?

Thanks in advance.

# Answer
Welcome! Bioinformatics is just a term that describes an approach to doing biology. People who use ( DNA sequencing/genomics/protein informatics + programming ) to answer biological questions are after the same answers as the the people who use ( pipettes/microscopes/molecular biology + a lab notebook ).

My best suggestion is to learn about some biological concepts that are interesting to you, see how people have applied bioinformatics to that branch of biology, and think if there is something you might be interested in doing to work on that question. It is critical to have a strong understanding of the underlying biology to get to the core of new findings using bioinformatics.

A lot of the hardcore algorithm development going on today is in genome assembly with long + noisy DNA sequencing reads, aligning long+noisy reads to genomes and trying to figure out the differences between the maternal and paternal haplotypes, working on evolutionary questions with population-level genetic variation data, and figuring out how things are related to each other with phylogenetics, to name a few.

As far as a machine-learning type approach which you seem to be hinting at in your question, one very hot topic at the moment is single-cell sequencing. This is a set of technologies that allows us to get DNA, RNA, or protein sequencing from a population of single cells to look at heterogeneity in tissues, timepoints, et cetera.

Good luck! If you have some free time on your hand this website, ROSALIND, is a fun project Euler-type website where you can learn some biology while programming bioinformatics solutions.
> 2 votes

---
Tags: algorithms
---


# Question
Title: vcftools: indel size histogram command returns empty file
I would like to get some summary statistics on a vcf file from one individual, which has over a million variant calls. I've tried to make a histogram of indel sizes with this command,

```
vcftools --vcf sample.vcf --out sample --hist-indel-len

```

but it returns an empty file. I see thousands of deletions in the vcf, although strangely no insertions. A typical line in the vcf file I have looks like

```
NC_025968.3 168363 191 N <DEL> . . SVTYPE=DEL;STRANDS=+-:19;SVLEN=-1026;END=169389;CIPOS=-9,330;CIEND=-152,9;CIPOS95=-3,61;CIEND95=-32,3;IMPRECISE;SU=19;PE=19;SR=0 GT:SU:PE:SR ./.:19:19:0

```

Update: I aligned reads to my reference with BWA MEM, then used samblaster to extract discordant pairs and split-reads. I sorted and indexed the three bam files, then used lumpyexpress to generate the vcf.

```
samblaster -i sample.sam --excludeDups --addMateTags --maxSplitCount 2 --minNonOverlap 20 | samtools view -S -b - > sample.bam

samtools view -h sample.bam | samblaster -a -e -d sample.disc.sam -s sample.split.sam -o /dev/null

lumpyexpress -B ../sample.sorted.bam -S ../sample.split.sorted.bam -D ../sample.disc.sorted.bam -o sample.vcf

```

I would like to see the size distribution of deletions, and if possible, figure out why there are no insertions.

# Answer
> 1 votes
It looks like `lumpy` looks for breakpoints in your reference. Presumably places where the is a haplotype switch detectable with the data on-hand?

My only thought is that since `lumpy` is being used to produce the VCF file that it is inherently not looking for insertions. I would suggest running the same BAM file through the GATK best practices pipeline (easier said than done) and vcftools pipeline. If those produce insertions but `lumpy` does not, it must be due to how `lumpy` produces the calls. Good luck!

---
Tags: vcf, ngs, vcftools
---


# Question
Title: cell type specific genes - heatmap using rank based approach
I was wondering if some one can help me with pseudo r code for my approach here based on this paper to plot a heatmap

# Figure 3

**Figure description** The top 40 enriched genes percelltypeare shown in a heat map.Onlyhighlyexpressed genes with FPKM ⬎20 areincludedin this analysis. Fold enrichment iscalculatedasFPKM of one cell type divided by the average FPKM of all other cell types. The majority of these genes showed specific expression by only one cell type, with the exception that some are expressed during more than one maturation stage in the oligodendrocyte lineage.

**Similar description in their methods:** Differential expression was calculated as the FPKM of a given cell type divided by the average FPKM of all other cell types. Genes were ranked by their fold enrichment in each cell type, and top enriched genes for each cell type were identified

In my case I have 9 cell types(columns) and ~24000 genes(rows) from mouse and the values are TPMs.

I know to order celltype only based on one column like this

```
expr.order = expr[order(expr[,1],decreasing =T),]

```

Which clearly not the right way and would need some help with code in r

# Answer
The figure caption says "the top 40 genes per celltype" - so we can simply use the code you provided for the top 40 genes for a single celltype, repeat this for all celltypes, and then remove any duplicate genes. Here's an example:

```
> expr <- matrix(runif(9*24000), nrow=24000, ncol=9)
> top_genes_each_celltype <- sapply(1:ncol(expr), function(i) order(expr[,i],decreasing =T)[1:40])
> dim(top_genes_each_celltype)
[1] 40 9
> top_genes_all_celltypes <- unique(as.vector(top_genes_each_celltype))
> length(top_genes_all_celltypes)
[1] 358

```

Then you can build a biclustered heatmap with `expr[top_genes_all_celltypes,]`, for example, using `heatmap2`.
> 2 votes

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


# Question
Title: How to query the Human Microbiome Project (HMP) to find all subjects with both 16s and WGS workups?
I am looking for a query to run on the HMP database that will return all subjects who have had BOTH 16s *and* whole genomes sequence (WGS) workups. I am currently using this query...

```
file.matrix_type = 16s_community OR file.matrix_type = wgs_community AND (file.format = FASTA OR file.format = FASTQ)

```

and sifting through the results by adding them to my shopping cart, browsing the file types until I have added all of the samples from one subject, and then I check to see if I have at least one 16s and one WGS. So far, I have not found one person who has had both 16s and WGS where both files were publicly accessible, and most do not have WGS at all. There must be a way to get subjects with both workups done to make this faster.

# Answer
After trying a few different queries, I don't believe there is a more specific one than this and I found that by searching and clicking through all of the results given by this query that eventually I found patients that fit the criteria of having both WGS and 16s.
> 0 votes

---
Tags: public-databases, genomics, data-download, metagenome
---


# Question
Title: Source for whole genome comparisons at NCBI Genomes
The NCBI Genomes database has these dendrograms for (presumably) whole genome comparisons for certain species, e.g. Pseudonomas aeruginosa or Escherichia coli.

How were these comparisons done? Someone knows the source / paper?

# Answer
I inquired into the details of the dendrograms after becoming frustrated with the lack of information. As with Ensembl, I'm sure that the folks at NCBI have a standardized pipeline that they run the sequences through to generate these dendrograms as no specific source deserving attribution was mentioned:

> The tree is based on a pairwise, BLAST comparison of chromosome sequences (an assembly is used for genomes with 2 or more chromosomes). BLAST identity is normalized on sequence length. The distance map is generated by a Neighbor Joining algorithm, where the distance is the blast score.
>
> *Dr. Wayne Matten, NCBI*
> 1 votes

---
Tags: phylogenetics, genbank
---


# Question
Title: Using external list of PCs for clustering
I am going to use principal components (PCs) comes from calcPCA function in URD program for clustering my cells in Seurat; So

```
seurat@dr$pca=URD_object@pca.load

> seurat <- FindClusters(
+     object = seurat,
+     reduction.type = "pca",
+     dims.use = 1:8,
+     resolution = 1.0,
+     print.output = 0,
+     save.SNN = TRUE
+ )
Error in GetDimReduction(object = object, reduction.type = reduction.type,  :
  cell.embeddings slot doesn't exist

```

There is `seuratToURD` function but there is not URD to seurat function I guess so I could use URD object for clustering in seurat How I can overcome this error please? I want to see when I am using PCA from URD what would be my clustering by SNN algorithm in seurat

Thanks a lot for any help

# Answer
> 3 votes
There is a way to do this, and even better--there is documentation for how to do it! No surprise coming from the Satija Lab. In the vignette they perform multidimensional scaling, but the idea is the same. `cmdscale()` returns the cell embeddings. `SetDimReduction()` is the Seurat function you are looking for. No manual editing using `@` required.

The authors use add mds embeddings, but if you want to add your own pca embeddings swap out `"mds"` for `"pca"` and `mds` for `URD_object@pca.load`.

```
# Before running MDS, we first calculate a distance matrix between all pairs
# of cells.  Here we use a simple euclidean distance metric on all genes,
# using object@scale.data as input
d <- dist(x = t(x = pbmc@scale.data))
# Run the MDS procedure, k determines the number of dimensions
mds <- cmdscale(d = d, k = 2)
# cmdscale returns the cell embeddings, we first label the columns to ensure
# downstream consistency
colnames(x = mds) <- paste0("MDS", 1:2)
# We will now store this as a new dimensional reduction called 'mds'
pbmc <- SetDimReduction(object = pbmc, reduction.type = "mds", slot = "cell.embeddings",
    new.data = mds)
pbmc <- SetDimReduction(object = pbmc, reduction.type = "mds", slot = "key",
    new.data = "MDS")

```

See the dimensional reduction vignette for more.

---
Tags: seurat, clustering, pca
---


# Question
Title: How to interpret duplication from MultiQC/FastQC?
FastQC of my sample files, aggregated into a single plot by MultiQC. Blue represents unique reads. Black represents duplicate reads. The x-axis is the number of reads.

I see I have more duplicate reads than unique reads, how do I interpret this plot? This is my attempt:

> FastQC took the first 100,000 sequences in each sample (details). Generally, higher duplicated read percentage is welcome as it represents better sequencing depth and lower sequencing error contamination. This dsitribution in this plot looks reasonable and nothing unusual.

# Answer
I would agree with your statement. Typically the ratio of unique to duplicated reads is a function of input library complexity (number of input molecules) verse sequencing depth (the amount of molecules sequenced).

In the example you gave you have about a 2:1 ratio of duplicates to unique reads which suggests that these libraries have been sequenced to their fullest and additional sequencing would yield few new unique reads.
> 0 votes

---
Tags: fastqc
---


# Question
Title: How to deal with duplicate genes having different expression values?
I have RNA-Seq data which is `FPKM`. In the dataframe `df` first column is `gene_name` and the other 100 columns are `samples`.

Usually if it is `counts` data I do like following:

```
df2 <- aggregate(. ~ gene_name, data = df, max)

```

I'm not sure what do with the `FPKM` data if there are duplicate genes with different `FPKM` value for the same sample.

Lets say it looks like this:

```
gene_name     sample1        sample2        sample3
5S_rRNA      0.3206147    0.3327312      0.377578
5S_rRNA      0.3342000    0.0000000      0.1305166

```

Any suggestions please.

# Answer
> 5 votes
I assume you're familiar with the various issues surrounding FPKMs, so I'll not expound upon them.

As a general rule, you should be using gene IDs rather than gene names, since the former are unique while the latter are not. If you only have access to data quantified on gene names, then the appropriate way to merge RPKMs is with a weighted sum:

$FPKM\_{gene} = \frac{FPKM\_{copy1} * Length\_{copy1} + FPKM\_{copy2} * Length\_{copy2}}{Length\_{copy1} + Length\_{copy2}}$

As an aside, rRNA expression levels will tend to be wrong, since they're normally excluded by either poly-A enrichment or the use of ribo-zero.

# Answer
> 3 votes
The best way to deal with this is to use unique gene IDs, for example ensembl accession numbers. So use the ensemble gtf annotation when quantifying the read counts and not the gene symbols. Just to illustrate, when I look for "5S\_rRNA" in ensembl's annotation, i see 18 different "genes" with that gene symbol. But which 2 you have is unclear now.

```
grep "5S_rRNA" ensembl_symbol.txt
"ENSG00000252830"       "5S_rRNA"
"ENSG00000276442"       "5S_rRNA"
"ENSG00000274408"       "5S_rRNA"
"ENSG00000274059"       "5S_rRNA"
"ENSG00000276861"       "5S_rRNA"
"ENSG00000274759"       "5S_rRNA"
"ENSG00000280646"       "5S_rRNA"
"ENSG00000277411"       "5S_rRNA"
"ENSG00000201285"       "5S_rRNA"
"ENSG00000212595"       "5S_rRNA"
"ENSG00000277418"       "5S_rRNA"
"ENSG00000277049"       "5S_rRNA"
"ENSG00000274097"       "5S_rRNA"
"ENSG00000277488"       "5S_rRNA"
"ENSG00000274663"       "5S_rRNA"
"ENSG00000283433"       "5S_rRNA"
"ENSG00000275305"       "5S_rRNA"
"ENSG00000278457"       "5S_rRNA"

```

---
Tags: r, rna-seq, gene-expression, fpkm
---


# Question
Title: In calculating the retention index, why do we use the character state with the lowest frequency?
To calculate the retention index for a phylogenetic tree, we use the following formula:

$$\frac{\text{maximum number of steps on tree - number of steps on the tree}} {\text{maximum number of steps on tree - minimum number of steps in the data}}$$

To calculate the maximum number of steps, we look at the number of observed character states for a character and take the lowest number. For each character, we take the lowest number of observed traits and then sum across all characters.

Take the following dataset (columns are characters, taxa are rows):

```
          123456
Outgroup  000000
Taxa A    100001
Taxa B    110001
Taxa C    111001
Taxa D    111101
Taxa E    111111
Taxa F    111110

```

The maximum number of steps for the first character is 1, since the character state 0 is attested once and character state 1 occurs six times.

What I don't understand is why we do not use 6 for the maximum number of changes for the first character.

# Answer
> 1 votes
The maximum number of steps (or changes) is the number of taxa with state 1 or 0, whichever is smaller. Conceptually, this is the number of character changes if each taxon evolved its state independently of the other taxa.

If for a given character the character state 0 occurred twice and character state 1 occurred five times, 2 would be the maximum number of changes.

---
Tags: phylogenetics, phylogeny, taxonomy, quantification
---


# Question
Title: Error in quantile.default(x, p = p) in EgdeR calcNormFactors
I am trying to run `TMM` normalization using `rpy2` and when I run `calcNormFactors()` function:

```
dge_list = edgeRLib.DGEList(counts=rawcounts)
dge_list = edgeRLib.calcNormFactors(dge_list, method="TMM")

```

I am getting the error:

> RRuntimeError: Error in quantile.default(x, p = p) : missing values and NaN's not allowed if 'na.rm' is FALSE

So, I thought that maybe my `rawcounts` `2D` `numpy` array contains `nan` values but it does not. I checked it with both `pd.isnull()` and `np.isnan()` functions. So, I suspect now that the `2D` array itself can not be normalized for some reason and I do not know why. The error is most probably not related to `rpy2` set up and `R` installation because it was working before for similar `numpy` arrays but failed this time. Any suggestions would be greatly appreciated.

# Answer
Ok, it turned out that some rows and/or columns contained all values as zeros and that was the reason for the error. I wrote a function to fix that:

```
def fix_rawcounts(rawcounts, cell_ids, gene_symbols):
    '''The function checks whether there are
    any rows or columns with all zeros
    and removes them'''
    print('Checking whether rows or columns contain all zeros')

    # Checking columns
    count_arr = np.apply_along_axis(np.count_nonzero, axis=0, arr=rawcounts )
    rawcounts = rawcounts[:,count_arr != 0]

    # Checking rows
    count_arr = np.apply_along_axis(np.count_nonzero, axis=1, arr=rawcounts )
    rawcounts = rawcounts[count_arr != 0,:]

    print('finished checking zero rows and columns')
    return rawcounts

```
> 4 votes

---
Tags: edger, tmm
---


# Question
Title: Find a cutoff value for genes that are expressed in single cell RNA-seq?
I want to find a cutoff value for each gene, above which we can consider a gene expressed.

The problem is that not all *effectively* non-expressed genes will have 0 counts due to sequencing errors for example.

My question is related to this one, but is concerned with single-cell data. The technical differences between bulk and single cell RNA-seq warrant a separate question for this problem. Also similar to this one, but focuses on gene expression and not differences between zero counts.

# Answer
Not sure if this has anything to do with your question but I think the way how Monocle approach this question sort of makes sense to me. The normalization methods for Seurat and Monocle are somehow similar (correct if I am wrong, this is simply my personal experience), as @kohlkopf mentioned.

For monocle specifically, it has this function where you can display dispersion v.s. expression of the genes to decide the cut off. Here is an example of this plot that I generated when processing the PBMC data downloaded from the Seurat tutorial.

You can assign the cut-off manually but I doubt I answered your question on how to decide the cut-off value. Maybe you need to try a set the expression cut-off from a range to experiment what value may potentially answer a biologically relevant question in your context.

If you are interested in how monocle approach this, I think the Monocle introduction and Dave Tang's Blog are excellent in explain how it works.

Again, you are right that 0 could be sequencing error. But I would also suggest that, if you are looking at genes at very small expression level, scRNA-seq, at least droplet-based technique, could be really difficult in such context. I happened to hear someone trying to find lncRNA expression talking about how a "0" sometimes may be a true "0" due the low expression nature of lncRNA in scRNA-seq, etc.

Anyhow, not sure if I answer the question but I do hope it helps you somehow. I am not able to make comments but this is the best I can do.
> 1 votes

# Answer
Much of the Seurat workflow incorporates normalized, scaled, and mean-centered expression values.

Raw gene expression counts for a cell are normalized by total expression, multiplied by 10,000, and log transformed. Expression level for each gene is scaled by dividing the centered gene expression levels by their standard deviations, resulting in z-scores. It's not an ideal solution, but I find a threshold of 1 (in the realm of z-scores) to make sense. A z-score of 1 is one standard deviation from the mean, or what may be interpreted as baseline expression.

This is my 2¢. I would love to hear what other people are doing.
> 0 votes

---
Tags: rna-seq, scrnaseq, gene-expression
---


# Question
Title: How to transform a DNAStringSet from the Bioconductor package Biostrings to a data frame?
I am working on Mac OS X. I am using R version 3.5.1.

I have important a FASTA file into R using `Biostrings::readDNAStringSet`. This creates a DNAStringSet object in R.

I would like to transform it into a data frame so that I can then plot the widths with ggplot2.

`as.data.frame` on my DNAStringSet object returns a data frame with only one column. Is there a function that returns a data frame with columns `width`, `seq` and `names` present in the DNAStringSet object?

# Answer
> 8 votes
According to the documentation (`?Biostrings::DNAStringSet`):

> `width(x)`: A vector of non-negative integers containing the number of letters for each element in `x`. Note that `width(x)` is also defined for a character vector with no NAs and is equivalent to `nchar(x, type="bytes")`.
>
> `names(x)`: `NULL` or a character vector of the same length as `x` containing a short user-provided description or comment for each element in `x`. These are the only data in an XStringSet object that can safely be changed by the user. All the other data are immutable! As a general recommendation, the user should never try to modify an object by accessing its slots directly.
>
> `as.character(x, use.names=TRUE)`: Converts `x` to a character vector of the same length as `x`. The `use.names` argument controls whether or not `names(x)` should be propagated to the names of the returned vector.

This should be enough for you to create a `data.frame`:

`dss2df <- function(dss) data.frame(width=width(dss), seq=as.character(dss), names=names(dss))`

Also, if you have a huge `DNAStringSet`, perhaps you do not want to convert all the sequences to characters; storing sequences as characters supposedly is less efficient. For plotting it might be sufficient to have just widths.

`qplot(width(dss), geom='histogram'))`

---
Tags: r, bioconductor, format-conversion, biostrings
---


# Question
Title: Can I stop my nanopore sequencing run if there are no more reads being produced?
I am sequencing a whole genome on MinION. I have used this flow cell for 6 hours for the phage lambda control. That gave me around 300,000 reads. I started a 48h sequencing run on the same flow cell after washing it. The sequencing has been going on for 1 day and 21 hours. It is not producing any reads at the moment (Strand = 0). The number of single pores is also 0. Can I "Stop acquisition" without negatively affecting the results?

The image below is a screenshot taken from MinKNOW, the software used to control the MinION.

# Answer
Yes, you can click "stop acquisition" and your run won't be negatively affected. All of the reads are saved as they are generated. I am not sure how this will impact live basecalling though if that is something that you do.
> 6 votes

# Answer
Yes, you can stop the run and it should not affect the base calling. We normally do the runs and then do the base calling after the run has completed.
> 1 votes

---
Tags: nanopore, ngs, minion, long-reads, minknow
---