# Question
Title: What is the file .mums created by mummer?
I am doing genome sequence alignment using MUMmer, in particular I want to do a dotplot with mummerplot. So the passages that I did are:

1.create a file .mums with the following command line:

```
mummer -mum -b -c H_pylori26695_Eslice.fasta H_pyloriJ99_Eslice.fasta > mummer.mums

```

2.and then plot all the MUMs with MUMmerplot:

```
mummerplot -x "[0,275287]" -y "[0,265111]" -postscript -p mummer mummer.mums

```

The example is taken here. The plot works, but I would like to understand several things to have more consciousness of what I am doing: the file .mums has three columns with many rows. What do they represent ? I know that it is a little explained in the link given before:

> This command will find all maximal unique matches (-mum) between the reference and query on both the forward and reverse strands (-b) and report all the match positions relative to the forward strand (-c). Output is to stdout, so we will redirect it into a file named mummer.mums. This file lists all of the MUMs of the default length or greater between the two input sequences.

So the file should report all the Maximal Unique Matches between the two sequences but I do not understand what are these three columns...

Moreover, can someone explain in simple words the reason of the term "unique" in "Maximal Unique Matches" ?

Thank you in advance .

# Answer
> 2 votes
## Maximal Unique Matches

To answer the question about MUMs, there are two important definitions:

* A "match" is *maximal* if it cannot be extended in either direction and still be a match i.e, for two strings $s,r$, a match $s\_i,...,s\_{i+k}$, $r\_{j},...,r\_{j+k}$ is maximal if $s\_i,...,s\_{i+k+1}$, $r\_{j},...,r\_{j+k+1}$ is not a match and $s\_{i-1},...,s\_{i+k}$, $r\_{j-1},...,r\_{j+k}$ is also not a match.
* A "match" is *unique* if the matched substring is not present anywhere else in either genome. i.e, for two strings $s,r$, a match $s\_i,...,s\_{i+k}$, $r\_{j},...,r\_{j+k}$ is unique if the substring $s\_i,...,s\_{i+k}$ is present only once in each genome.

With those definitions in hand, we can easily define a Maximal Unique Match (MUM) as an exact match that is both maximal and unique. Note that MUMmer also provides support for Maximal Exact Matches (MEMs), which do not require uniqueness of matches.

## MUMmer output

As far as the output goes:

> For each match, the three columns list the position in the reference sequence, the position in the query sequence, and the length of the match respectively. Reverse complemented query positions are reported relative to the reverse of the query sequence unless the `-c` option was used.

For example, the row:

```
1      1120       299

```

denotes a MUM of length 299 starting at the first base in the first input file and the 1120th base in the second input file. For matches on the reverse complement strand, you'll notice a second part of the output that contains the matches on the reverse complemented strand.

Here is a sample output I have on my device:

```
>$ $  mummer -mum -b -c England1.fna Indiana-USA-1_Saudi_Arabia_2014.fna
> gi|633896549|gb|KJ813439.1|
       1      1120       299
     301      1420       174
     476      1595        65
     542      1661        81
     624      1743        43
     668      1787       264
     933      2052      1042
    1976      3095      1299
    3276      4395       330
    3607      4726        69
    3677      4796       885
    4563      5682       812
    5376      6495       792
    6169      7288       831
    7001      8120       731
    7733      8852      1963
    9697     10816      1176
   10874     11993       617
   11492     12611        41
   11534     12653       164
   11708     12827      3007
   14716     15835      1040
   15757     16876       489
   16247     17366       412
   16660     17779       337
   16998     18117       463
   17462     18581       325
   17788     18907       686
   18475     19594       599
   19075     20194       342
   19418     20537       179
   19598     20717       532
   20131     21250       152
   20284     21403       935
   21220     22339       348
   21569     22688       764
   22334     23453       351
   22686     23805        73
   22760     23879        29
   22790     23909       180
   22971     24090       532
   23504     24623       143
   23648     24767       107
   23756     24875       434
   24191     25310        59
   24251     25370        47
   24299     25418       214
   24515     25634       224
   24740     25859      1386
   26127     27246        39
   26167     27286        63
   26231     27350       114
   26346     27465       102
   26449     27568      1399
   27849     28968       207
   28057     29176       938
   28996     30115       726
   29743     30868       369
> gi|633896549|gb|KJ813439.1| Reverse

```

Now, if you compute the reverse complement of the input and rerun, you'll notice the matches appear in reverse in the second section denoted by `> gi|633896549|gb|KJ813439.1| Reverse`

# Answer
> 0 votes
found this in the manual, so I guess the first two columns are the reference and query positions and the third is the length. I am not 100% sure either, so let me know if you find out more.

http://mummer.sourceforge.net/manual/#mummerplot

> Just a simple three column list, noting the position and length of every maximal exact match. Note that for reverse complement matches (produced with the -r option), the query start positions will reference the reverse complement of the query input sequence.

---
Tags: sequence-alignment, genome, phylogenetics, file-formats, software-usage
---


# Question
Title: How do I access sequence and MAPQ from bam using HTSlib in C++?
I am using the "htslib/sam.h" header in C++. I need to access the SEQ field of bam reads and store each sequence in a vector of strings. I also need to access the MAPQ field of each read to check if the read passes my threshold quality of 10 (which I have stored as an integer). Does anyone know which htslib/sam.h functions can accomplish these tasks?

# Answer
The MAPQ is stored in `b->core.qual`, given a `bam1_t *b`. You can access the sequence with `bam_seqi(bam_get_seq(b), some_offset)`(see https://github.com/samtools/htslib/blob/develop/htslib/sam.h#L317). This will get you a base at a time, so you'd need to iterate over it (the length of the sequence being in `b->core.l_qseq`).
> 3 votes

---
Tags: sam, c++, htslib
---


# Question
Title: The biological meaning of the random variables and the responses in Seurat analysis
In the Seurat analysis, if we suppose that Xg and Xr denote the random variables that associate to the expression level of the gene g and of the gene r, respectively. Let Y and Yv represent the responses for Xg and Xr, respectively.

Is it true that biologically, the random variables Xg and Xr represent the distributions of the gene expression measure per cell and the gene expression measure per spot, respectively and the responses Y and Yv represent the range of the random variable Xg and Xr, respectively. The values for Y and Yv are the cell types?

Reference:

https://projecteuclid.org/journals/statistical-science/volume-18/issue-1/Multiple-Hypothesis-Testing-in-Microarray-Experiments/10.1214/ss/1056397487.full

\[Method section 2.1\]

I try to understand what is the Random variable and the response in this context biologically?

I am confused to make a biological statements about these response. Please help. Thank you.

# Answer
> 4 votes
In the linked article the authors formalize microarray analysis as the study of the joint distributions of $\overrightarrow{X}\_i$ and $Y\_i$, where $\overrightarrow{X}\_i$ is a vector of random variables, the distribution of each of which (i.e. each of $X\_{ji}$) is determined by the level of expression of gene j in sample i, and $Y$ is some response or covariate of interest.

The response/covariate is the easiest to interpret in terms of biology. It is the condition. So it might be treated or untreated if we were interested in genes that respond to drug treatment, or it might be diseased vs healthy - this would be covariates. Or it might be "time to death", "body mass index", etc. - these would be responses. The difference between a co-variate and a response is which is the dependent and which is the independent variable. In differential expression we have "gene expression ~ disease" (i.e. X ~ covariate), where `~` is pronounced "is predicted by". In the case of survival analysis "time to death ~ gene expression" (i.e. response ~ X).The malleability of the relationship between X and Y is what is meant by

> "The population and sampling mechanism will depend on the particular application"

A random variable is any variable that we cannot know the value of until it has been measured. Note, this is *not* "do not" know the value of, its "can not" know the value of. So for example, in this case of RNA-seq, X is not the gene expression, it is the number of reads. An all-knowing daemon could know the gene expression for a gene in a sample before we measured it, but they couldn't know how many reads we will get because of the randomness introduced by sampling.

Random variables don't *represent* distributions, they *have* distributions. That is, if you were to realize many trials of a random variable, the values you obtained would form a distribution.

So the number of reads you get mapping to a particular gene in a sequencing experiment is a random variable. As is the brightness of a spot on a microarray. Its an interesting question to ask if $Y$ is a random variable. The answer is "it depends" - and what it depends on is how good we are at controlling an expriment. You might argue that "drug treatment" is not a random variable because we treat the cells or not - we know, but that, say "time to death" or even "age" is a random variable, because we don't determine it ahead of the experiment. Either way, this subtlety is covered by the above quoted sentence.

So, in a traditional bulk sequencing experiment we have a matrix of random variables associated with gene expression $X\_{ji}$ \- the number of reads for gene j in sample i, and $Y\_i$ the condition for sample i. Note that $Y$ is indexed by the sample, not the gene - every gene in the sample uses the same value of $Y$. And for each gene j we test the model `X_i ~ Y_i` \- the expression in sample i is predicted by the condition of sample i.

In a single cell sequencing our matrix of random variables associated with gene expression is actaully a tensor, that is, it has three indexes, not two: $X\_{jki}$ \- that is the number of reads from gene j in cell k in sample i. In a scarily large number of single cell sequencing experiments, i only has one value (i.e. there is only one sample). And in most single cell experiments I've seen $Y\_i$ only takes one value, even if there are multiple samples. Even if $Y$ takes multiple values, its not clear what sort of model you'd fit to the structure `X_{ki} ~ Y_i` \- i.e predicting a matrix of many cells by a single vector of covariate values.

I think one analysis you can do with Seurat is to treat each cell as a separate sample, and each cluster as a condition, and thus try to find gene DE between cell types. I'm not convinced how great this sort of analysis, but here, the i in $X\_{ji}$ would represent which cell a measurement was from and $Y\_{i}$ would represent which cluster. But this sort of analysis seems prone to all sorts of date-peeking and training set leakage problems.

> if we suppose that Xg and Xr denote the random variables that associate to the expression level of the gene g and of the gene r, respectively. Let Y and Yv represent the responses for Xg and Xr, respectively.

First, $Y$ is a vector of which $Y\_v$ is a component, so $Y$ contains $Y\_v$. And in the formulation presented in the linked paper (and in almost all RNAseq analysis literature), the response is the same for two different genes, so there is only one $Y$ across both genes g and r.

> the random variables Xg and Xr represent the distributions of the gene expression measure per cell and the gene expression measure per spot, respectively

I'm not quite sure what you mean by "per spot". Spots are a microarray term, and don't really have any counterpart in RNA-seq analysis. They were, literally, the spots of DNA oligos printed onto the glass slides used in microarray analysis.

$x\_{gi}$ is the number of reads mapped to gene g in cell i, and has a distribution determined by the level of expression of gene g in the cell i.

$x\_{ri}$ is the number of reads mapped to gene r in cell i, and has a distribution determined by the level of expression of gene r in cell i.

> and the responses Y and Yv represent the range of the random variable Xg and Xr, respectively

No. The response $Y$ represents some biological dicotomy you would like to associate with gene expression. You could use cell type, but this is neither response, nor covariate, it is not measured independently of X, but is simply a phantom of the correlation structure in $X\_{ji}$, so I can see why the "response" type language would be confusing.

---
Tags: gene, seurat, differential-expression, gene-expression
---


# Question
Title: Is there a way to calculate the Ka/Ks value in R/unix?
I am investigating the rate of synonymous non synonymous substitutions and would like to find the Ka/Ks value also known as dn and ds between multiple alleles? I want to investigate multiple alleles because I am examining the MHC region which is highly polymorphic. Any suggestions would be very helpful!

# Answer
I have also used the seqinr package for Ka/Ks
> 2 votes

# Answer
A while back, I used `PAML` for this purpose for a phylogenetics project. It has Windows and macOS executables but for \*nix, you will need to compile it from source code.
> 1 votes

---
Tags: r, snp, bash
---


# Question
Title: Removing cells zero for a gene from a scRNA-seq data
I have a big single-cell RNA seq data

```
> dput(head(new.dat[,1:10]))
structure(list(cell1 = c(0.793763840992639, 0, 1.96843530982957,
0.461736429639991, 0.717968540649498, 0), cell2 = c(3.61741696702738,
0.231662370550224, 0, 0, 0, 0), cell3 = c(4.14348883366621, 0.118161316317251,
0.08074552209482, 2.27968429766934, 0.0470313356296409, 0), cell4 = c(1.34783143327084,
0.0094666040612932, 1.14392942941128, 0.652535826921119, 0.357542816432864,
0.149587369334621), cell5 = c(1.27104023273899, 1.55185229643731,
0, 0, 0, 0.0117525723115277), cell6 = c(1.92307653575663, 0,
0, 0.319156642478379, 0, 0), cell7 = c(3.9343015424917, 0.132824589520901,
0.119679885703561, 0.772516422897241, 0.0236884909844904, 0),
    cell8 = c(3.74969491678643, 0.103404975609384, 0.0354753982873036,
    0, 0, 0), cell9 = c(1.19084857532713, 3.9213265721495, 0,
    0.0341973245272891, 0.0419122921627454, 0), cell10 = c(4.1224255501566,
    0.301871669274068, 0.0633536200981225, 0.389959552469879,
    0, 0.0405296102106492)), row.names = c("PTPRC", "MHC-II",
"ITGAM", "Ly6C", "Ly6G", "EMR1"), class = "data.frame")
>

> dim(new.dat)
[1]     33 263086
>

```

How I `remove every columns` which are `zero` for one gene, let's say `PTPRC`?

# Answer
> -3 votes
```
WhichCells(seurat_object made by this matrix, slot = 'counts', expression = PTPRC > 0 )

```

By Seurat R package this is possible

# Answer
> 0 votes
This is basic R subsetting operation. Most subset operations give you a subset of rows based on a column, while you're doing the opposite here - but the logic is the same.

This gives you a subset of rows when the expression is based on a column

```
df[logical_expression_that_subsets_rows,]

```

Why not try this - it should give you a subset of columns based on a row value. (Untested, but you should be able to build on this):

```
df[, logical_expression_that_subsets_columns]

```

The logical expression for the first could be `df$cell1 < 0.5` (or the equivalent `df[,"cell1"] < 0.5`. You should be able to get a logical expression for the second based on this.

---
Tags: r, scrnaseq, 10x-genomics
---


# Question
Title: Seurat clustering Methods-resolution parameter explanation
I am learning the Seurat algorithms to cluster the scRNA-seq datasets. I found this explanation, but am confused. Can someone explain it to me, "The FindClusters function implements the procedure, and contains a resolution parameter that sets the ‘granularity’ of the downstream clustering, with increased values leading to a greater number of clusters. We find that setting this parameter between 0.6-1.2 typically returns good results for single cell datasets of around 3K cells". My questions are:

1. What is "granularity" of the downstream analysis?
2. How does the parameter "resolution" work?
3. How to know the optimum value for parameter "resolution"?

I found this paper about the resolution that might help, but not clear to me yet. https://iopscience.iop.org/article/10.1088/1742-5468/2008/10/P10008/pdf.

I think the higher the resolution, the more clusters produced meaning the smaller clusters obtained? The smaller the clusters the more specific cell types we get, is that correct?

# Answer
> What is "granularity" of the downstream analysis?

Low granularity implies basic cell types (e.g. T cells, B cells), implies few larger clusters.

High granularity implies more refined cell types (e.g. CD4+ Central Memory T, Treg, Memory B cell etc.), implies more smaller clusters

> How does the parameter "resolution" work?

A high resolution parameter value produces more smaller clusters, i.e. higher granularity. If you want to know how modularity clustering works, I recommend the Leiden algorithm paper https://www.nature.com/articles/s41598-019-41695-z

> How to know the optimum value for parameter "resolution"?

Tricky, and it depends on what questions you want to answer. Normally it's a subjective iterative process of trying a clustering resolution, checking if the clusters make sense, re-trying etc. but I also quite like the clustree approach https://cran.r-project.org/web/packages/clustree/index.html

> The smaller the clusters the more specific cell types we get, is that correct?

Yes correct
> 2 votes

---
Tags: scrnaseq, seurat, single-cell, clustering
---


title: Using txt file to query Ensembl API

text:
# Question
Title: Using txt file to query Ensembl API
I am trying to use an external file with a list of gene names to query the ENSEMBL IDs using the perl API. This script works:

```
#!/usr/local/bin/perl

use strict;
use warnings;
use Bio::Ensembl::Registry;

my $registry = 'Bio::EnsEMBL::Registry';

$registry->load_registry_from_db(
    -host => 'useastdb.ensembl.org',
    -port => 5306,
    -user => 'anonymous'
);

my $gene_adaptor = $registry->get_adaptor( 'Human', 'Core', 'Gene' );

my @genes;
@genes = ('CCL2', 'CCL5');

foreach $name(@genes){
    my $gene = $gene_adaptor -> fetch_by_display_label($name);
    print "GENE ", $gene->stable_id(), ",";
}

```

However I can't figure out how to get this to work with an external file that contains a list of genes like this:

`genes.txt`

```
CCL2
CCL5

```

Here's what I tried:

```
#!/usr/local/bin/perl

use strict;
use warnings;
use Bio::Ensembl::Registry;

my $registry = 'Bio::EnsEMBL::Registry';

$registry->load_registry_from_db(
    -host => 'useastdb.ensembl.org',
    -port => 5306,
    -user => 'anonymous'
    );

my $gene_adaptor = $registry->get_adaptor( 'Human', 'Core', 'Gene' );

open (IDs, '<', "genes.txt") or die "can't read input file";
while (<IDs>){
    my $gene = $gene_adaptor -> fetch_by_display_label($_);
    print "GENE ", $gene->stable_id(), ",";
}

```

But I'm getting this error: `Can't call method "stable_id" on an undefined value at ./get_ensembl_ids.pl line 22, <IDs> line 1.`

I was trying to follow this tutorial for reading txt files in perl, but what am I missing for getting perl to actually fetch the ENSEMBL ids with each of the .txt file lines?

# Answer
Similar to reading files in other languages, it is necessary to trim the newlines from each line. If you add `chomp;` to the first line of the `while` loop, the code runs just fine (note that I also print a newline at the end of the loop).

```
open (IDs, '<', "genes.txt") or die "can't read input file";
while (<IDs>){
    chomp;
    my $gene = $gene_adaptor -> fetch_by_display_label($_);
    print "GENE ", $gene->stable_id(), ",";
}
print '\n';

```

The new loop produces the following result:

```
>$ perl query.pl
GENE ENSG00000108691,GENE ENSG00000271503,

```

Additionally, you may need to catch any exceptions/errors thrown if you have an invalid gene in your input file:

```
open (IDs, '<', "genes.txt") or die "can't read input file";
while (<IDs>){
    chomp;
    eval {
        my $gene = $gene_adaptor -> fetch_by_display_label($_);
        print "GENE ", $gene->stable_id(), ",";
    };
    if ($@) {
        warn "Error:\t$@\n"
    }
}
print '\n';

```
> 2 votes

---
Tags: ensembl, perl
---


# Question
Title: Time complexity of the Neighbors function
I implemented the following **find\_neighbors\_with\_expected\_hamming\_distance** function that generates all k-mers of Hamming distance at most d from the given Pattern.

```
from .Hamming_Distance import find_hamming_distance

dp_map_that_stores_all_neighbors_with_expected_hamming_distance = {}

unique_nucleotides = {'A', 'C', 'G', 'T'}

def suffix(pattern):
    if len(pattern) <= 1:
        return None
    else:
        return pattern[1: (len(pattern))]

def first_symbol(pattern):
    if len(pattern) == 0:
        return None
    else:
        return pattern[0]

def find_neighbors_with_expected_hamming_distance(pattern, expected_hamming_distance):
    if expected_hamming_distance == 0:
        return {pattern}

    if len(pattern) == 1:
        return {'A', 'C', 'G', 'T'}

    if (pattern, expected_hamming_distance) in dp_map_that_stores_all_neighbors_with_expected_hamming_distance.keys():
        return dp_map_that_stores_all_neighbors_with_expected_hamming_distance[(pattern, expected_hamming_distance)]

    neighborhood_set = set()

    suffix_neighbors_set = find_neighbors_with_expected_hamming_distance(suffix(pattern), expected_hamming_distance)

    for suffix_neighbor in suffix_neighbors_set:

        if find_hamming_distance(suffix(pattern), suffix_neighbor) < expected_hamming_distance:

            for nucleotide in unique_nucleotides:
                neighborhood_set.add(nucleotide + suffix_neighbor)

        else:
            neighborhood_set.add(first_symbol(pattern) + suffix_neighbor)

    if (pattern, expected_hamming_distance) not in dp_map_that_stores_all_neighbors_with_expected_hamming_distance.keys():
        dp_map_that_stores_all_neighbors_with_expected_hamming_distance[
            (pattern, expected_hamming_distance)] = neighborhood_set

    return neighborhood_set

```

**I want your suggestion on my time complexity calculation:**

For a given pattern of length **N** and **D** as maximum expected hamming distance, I calculated time complexity as -

$\sum\_{n=1}^N$ $\sum\_{d=1}^D$ $n\choose d$ * $3^i$

**Reasoning:** For any length **n**, we will choose any **d** letters (nucleotide in this case) to be replaced with the other 3 nucleotides that is different from the nucleotide present at that particular position.

For example: Pattern = ATGCAT and **D** = 2, then suppose for one of the case, I select 2nd and 4th positions to be replaced.

Then, at the second position (current nucleotide = 'T'), I can substitute any nucleotide from the set {'A', 'C', 'G'} and

for fourth position (current nucleotide = 'C'), I can substitute any nucleotide from the set {'A', 'T', 'G'}.

**Kindly tell me if this approach is correct?**

**EDIT:**

I have provided the pseudocode for the **Neighbors** function as follows:

```
Neighbors(Pattern, d)
    if d = 0
        return {Pattern}
    if |Pattern| = 1
        return {A, C, G, T}
    Neighborhood ← an empty set
    SuffixNeighbors ← Neighbors(Suffix(Pattern), d)
    for each string Text from SuffixNeighbors
        if HammingDistance(Suffix(Pattern), Text) < d
            for each nucleotide x
                add x • Text to Neighborhood
        else
            add FirstSymbol(Pattern) • Text to Neighborhood
    return Neighborhood

```

The symbol • means concatenation of 2 strings. For example, 'S' • 'V' = 'SV'.

If we remove the first symbol of Pattern (denoted FirstSymbol(Pattern)), then we will obtain a

(k − 1)-mer that we denote by Suffix(Pattern).

**Illustration:**

Given pattern = 'ACG'

Expected Hamming distance (D) = 1

Answer = {'CCG', 'TCG', 'GCG', 'AAG', 'ATG', 'AGG', 'ACA', 'ACC', 'ACT', 'ACG'}

Note that the Answer is a set of all strings who have hamming distance of 1 with respect to the given pattern 'ACG'.

# Answer
> 2 votes
Let's first establish a lower bound on any solution:

Given a hamming distance $D$ and an input pattern $P$ of length $N$, there are $\sum\_{k=0}^D{N \choose k} 3^k$ other strings within hamming distance $D$ of $P$. Of course in the worst case, when $D=N$, there are $4^n$ possible strings. We can also get a closed bound on the number of solutions. We can show that the sum of the first $D$ binomial coefficients for a fixed $N$ is bounded by ${N \choose D} {N-(D-1) \over N-(2D-1)}$ as shown here. Therefore we have the following closed upper bound on the number of strings within hamming distance $D$ of some pattern of length $N$

$$ \sum\_{k=0}^D{N \choose k} 3^k \leq {N \choose D} {3^D (N-(D-1)) \over N-(2D-1)} \leq \left(\frac{3N}{D}\right)^D {N-(D-1) \over N-(2D-1)} $$

However, for your approach, you are going through recursively, so when the subproblem is on strings within $D$ of the suffix of length $n$, you have $\sum\_{k=0}^D{n \choose k} 3^k$ of these. Given that you have $N$ subproblems (suffixes), at first glance it does indeed look like you have roughly $\sum\_{n=1}^N\sum\_{k=0}^D{n \choose k} 3^k$ iterations of the innermost for loop. However, in that `for` loop, you are also computing the Hamming Distance function, which is $O(n)$. This means that you have roughly

$$ \sum\_{n=1}^Nn\sum\_{k=0}^D{n \choose k} 3^k $$

computations.

I would suggest finding a way to get rid of having to compute the hamming function. Since you are constructing the sequences yourself, you can keep track of each sequences hamming distance from the pattern during construction.

---
Tags: python, k-mer, algorithms
---


# Question
Title: Can HTSlib extract bam reads occurring in a specific region without iterating through the whole file?
I am using htslib/sam.h to write a C++ program. As part of this program, I must extract reads occurring on specific scaffolds in specific regions from a bam file. Essentially, I want to perform the following samtools command in C++:

```
samtools view name.bam -o name.sam scaffold:pos-pos

```

Since I must extract reads from thousands of regions, I do not want to iterate through the whole bam file each time I extract reads from a region. I am wondering if HTSlib can extract reads occurring in a specified scaffold and range without iterating through the whole bam. I found the function bam\_itr\_querys() in the htslib/sam.h documentation on GitHub, but don't know how to format the inputs for this function, or what the "idx" input is. I also saw somewhere that I would need to pair this function with bam\_itr\_next(), but don't know the format to do so.

EDIT: My bam files are sorted and indexed, if that matters.

# Answer
```
#include <stdio.h>
#include "htslib/sam.h"
int main(int argc, char *argv[])
{
  samFile *fp = sam_open(argv[1], "r");
  hts_idx_t *idx = sam_index_load(fp, argv[1]);
  bam_hdr_t *h = sam_hdr_read(fp);
  bam1_t *b = bam_init1();
  hts_itr_t *itr = sam_itr_querys(idx, h, argc < 3? "." : argv[2]);
  while (sam_itr_next(fp, itr, b) >= 0) {
    if (b->core.tid < 0) continue;
    printf("%s\t%d\t%d\n", h->target_name[b->core.tid], b->core.pos,
      b->core.pos + bam_cigar2rlen(b->core.n_cigar, bam_get_cigar(b)));
  }
  bam_destroy1(b);
  hts_itr_destroy(itr);
  bam_hdr_destroy(h);
  sam_close(fp);
}

```

The usage would then be `./a.out alignments.bam chr1:123-456`, assuming your BAM file is sorted an indexed. Importantly this is making use of the index via `sam_itr_querys()` so it's not iterating over the entire file.
> 2 votes

# Answer
Have you considered using a sorted and indexed bam file? e.g.

```
samtools sort in.bam > out.sorted.bam
samtools index out.sorted.bam

```

This makes arbitrary region lookups trivial. I am more familiar with pysam (HTSlib for python), which has the `AlignmentFile.fetch()` functionality.

Surely there is an analog in HTSlib API, though I don't know what to search to get it. It is in samtools as demonstrated here.

This answer seems to put the components together that you are interested in, among them `sam_index_load()`, `sam_itr_querys()`.
> 0 votes

---
Tags: sam, samtools, c++, htslib
---


# Question
Title: P value recalculation from multiple analysis (Transcriptomics)
I've done a differential expression analysis with transcriptomic information (microarrays and RNAseq) from different datasets by applying for some of them t-test and for the rest Wald-test. I've got many coincidences for these significant genes so I want to recalculate my p-values (or my adjusted p-values, obtained by the BH method) using all these datasets. What is the proper method to do this (if possible)? I am completely new in this field so any additional information would be appreciated.

# Answer
For cDNA-derived data have a look at the DESeq2 analysis guide.

The general process is to convert all your input data into an integer *count* matrix (not expression matrix), with genes in rows and samples in columns, and pair it with a metadata data frame that has one line per sample (ordered in the same way as the columns), and indicates the conditions for each sample (including sequencing batch). By incorporating the sequencing batch into the statistical model, it's possible to reduce any technical error associated with consistent linear changes in expression based on the batch (e.g. different runs giving different numbers of total reads).

Microarrays *can* be combined, bearing in mind that microarray results are often inconsistently presented as normalised / non-normalised, or linear / log. I've found that normalising within-batch then converting to log space usually seems to work well. There is a verbose Users Guide for Limma Voom, but I've found it to be tricky to navigate - it's less friendly to new users.

Combining microarrays and RNASeq is complicated, because the data have different dynamic ranges, and different biases. I think the best thing to do is to process the microarray data separately from the RNASeq data and generate two lists of biologically-interesting genes, then look at the set relationships (e.g. intersection, exclusion).
> 0 votes

# Answer
> What is the proper method to do this (if possible)?

Please use dedicated software for RNA-seq and microarray analysis and **not** naive tests such as the t-test. For microarray one commonly uses `limma` and for RNA-seq either `DESeq2`, `edgeR` or `limma-voom`, all are available as Bioconductor packages in R.
> 1 votes

---
Tags: statistics, microarray, p-values
---


# Question
Title: How to export summary data from Bowtie2 for MultiQC to read?
BowTie2 puts out summary info to the terminal but doesn't allow me to save each to a separate file. How can I go about saving the output of BowTie2 so that MultiQC can read it and create a website from it?

# Answer
After trying on my own this is the solution I came up with. My solution was to use a pipe of the stderr into a new file.

```
bowtie2 --local --threads 4 -x reference_genome.fa -f -U inputdata.fa -S OutputData.sam --no-hd --no-sq --no-unal -k 5 2> 'bowtie2_summary/inputdata.txt'

```
> 1 votes

---
Tags: bowtie2
---


# Question
Title: Pseudocode for gapped BLAST seed extension
I am trying to understand in detail how does the gapped BLAST seed extension DP algorithm works. Despite searching for days now, I could not find anywhere a pseudocode that would allow implementing this in python. Not being a C programmer, it would be impossible for me to grasp this from the BLAST source code.

Can anyone help me with this?

Thank you all. Best,

Daniel

# Answer
So it turns out the Gapped Blast paper doesn't actually give the extension algorithm. It gives a high level overview and refers to an article in preparation at the time:

> The central idea is to consider only cells for which the optimal local alignment score falls no more than $X\_g$ below the best alignment score yet found. Starting from a single aligned pair of residues, called the seed, the dynamic programming proceeds both forward and backward through the path graph (Zheng Zhang et al., manuscript in preparation) (Figs 3a and 4). The advantage of this approach is that the region of the path graph explored adapts to the alignment being constructed. The alignment can wander arbitrarily many diagonals away from the seed, but the number of cells expanded on each row tends to remain limited, and may even shrink to zero before a boundary of the path graph is encountered (Fig. 4). The $X\_g$ parameter serves a similar function to the band-width parameter of the earlier heuristic, but the region of the path graph it implicitly specifies be explored is in general more productively chosen.

The motivation for this approach is to be a "smarter" version of the banded alignment approach. The issue with banded alignment approaches is that they implicitly won't allow for large insertions/deletions within the alignment, even if the areas around the alignment are high scoring. Increasing the band with to allow for this reduces the speedup gained by the band constraint. This is where the idea of "Alignments Without Low-Scoring Regions" comes into play. This paper is the one cited as "in-prep" by the gapped blast paper, and it describes an early version of the gapped extension approach. However, that approach seems to have been superseded by the one described here 2 years later by the same lead author (the former paper in '98 having 67 citations and the latter in 2000 having 4700+).

For completeness, I've pasted some pseudocode from the more popular x-drop paper from 2000. Note that $X$ is a user specified parameter, with larger $X$ allowing for a larger search space but slower speed.
> 3 votes

---
Tags: blast
---


# Question
Title: Histosketch vs Count-Min sketch: preserving similarity
A preprint describing a new tool and its application to microbiome analysis was recently published in bioRxiv\[1\]. At the core of this new tool, HULK, is a new data structure called a histosketch\[2\] which is similar in spirit to CountMin sketches and MinHash sketches that have become quite popular in bioinformatics in the last several years.

One of the main advantages the histosketch authors claim over CountMin sketch is that CountMin sketches do not preserve similarity. From section II of their paper\[2\] (emphasis mine):

> It is practically difficult to compute histograms for data streams with typically unknown cardinality and which thus require and unbound amount of memory to maintain the histogram. In this context, the count sketch and the count-min sketch and other online histogram building methods were proposed to approximate the frequency table of elements (i.e., histograms) from a data stream with a fixed-size data structure. However, **the resulting sketches do not preserve similarity between different data streams**.

In what sense does the CountMin sketch not preserve similarity between streams? Does this refer to similarity at intermediate states (while data is still streaming) or does it refer to inherent characteristics of the data structure itself?

---

\[1\] Rowe WPM, Carrieri AP, Alcon-Giner C, Caim S, Shaw A, Sim K, Kroll JS, Hall L, Pyzer-Knapp EO, Winn MD (2018) Streaming histogram sketching for rapid microbiome analytics. *bioRxiv*, 408070, doi:10.1101/408070.

\[2\] Yang D, Li B, Rettig L, Cudré-Mauroux P (2017) HistoSketch: Fast similarity-preserving sketching of streaming histograms with concept drift. *Proc. IEEE Int. Conf. Data Mining*, doi:10.1109/ICDM.2017.64.

# Answer
> 1 votes
*Answer from @user172818, converted from comment:*

My understanding is that you can't compute Jaccard similarity given two count-min sketches. You can with two histosketches, though I don't know exactly how yet.

---
Tags: k-mer, algorithms, sketch
---


# Question
Title: Get genomic coordinates using GenomicFeatures by HGNC gene names
I want to get coordinates of human genes from my list (consisting of hgnc genes id) using GenomicFeatures and TxDb.Hsapiens.UCSC.hg19.knownGene R packages from Bioconductor.

```
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb=(TxDb.Hsapiens.UCSC.hg19.knownGene)

my_genes = c("INO80","NASP","INO80D","SMARCA1")

select(txdb, keys = my_genes,
   columns=c("TXCHROM","TXSTART","TXEND","TXSTRAND"),
   keytype="GENEID")

```

However, it doesn't' work because txdb doesn't take hgnc identifiers; how can it be solved? I couldn't find any appropriate keytype that will support hgnc and not sure how to match hgnc id I have and GENEID from txdb.

# Answer
*Answer from @devon-ryan, converted from comment:*

Try a UCSC gene ID instead. For those genes they'd be: `uc001zni.3` (INO80), `uc002vaz.4` (INO80D), `uc004eup.4` or `uc004eun.4` (SMARCA1) and `uc001coj.2`, `uc010olq.2`, `uc010olr.2`, `uc001coi.2` or `uc001col.2` (NASP). You might have to omit the versions (the `.4` part, for example).
> 0 votes

---
Tags: r, bioconductor, sequence-annotation
---


title: Combine VCF files

text:
# Question
Title: Combine VCF files
I have a problem with using `rbind` to combine VCF files using the library VariantAnnotation from Bioconductor.

I am reading two VCF files, when I try to combine them in a certain order with `rbind` I'm getting an error. When I combine them in a different order, it works fine.

```
vcf.full <- rbind(vcf1, vcf2) ## Doesnt work

```

> Error in as(from, to\_class, strict = FALSE) : no method or default for coercing “CompressedCharacterList” to “DNAStringSetList”

```
vcf.full <- rbind(vcf2, vcf1) ## Work

header(vcf1)
class: VCFHeader
samples(503): NA06984 NA06989 ... NA20828 NA20832
meta(2): META contig
fixed(2): FILTER ALT
info(27): CIEND CIPOS ... EX_TARGET MULTI_ALLELIC
geno(1): GT

header(vcf2)
class: VCFHeader
samples(503): NA06984 NA06989 ... NA20828 NA20832
meta(2): META contig
fixed(2): FILTER ALT
info(27): CIEND CIPOS ... EX_TARGET MULTI_ALLELIC
geno(1): GT

dim(vcf1)
[1]  30 503

dim(vcf2)
[1] 149 503

```

The two VCF files seem identical. Any idea of why it happens? Or any suggestions on a different method to combine multiple VCF files using R?

# Answer
> 1 votes
*Answer from @Kozolovska, converted from comment:*

Found the problem, some of ALT were structural variant, so the list of vcf2@info$ALT was a charterer list and not a 'DNAStringSetList' after removing that variant it worked. Thanks.

---
Tags: r, vcf
---


title: SPAdes error during assembly

text:
# Question
Title: SPAdes error during assembly
I try to perform a hybrid assembly using Unicycler. Unicycler use SPAdes to assembly the illumina sequences, and the program crashes while in SPAdes. The spades.log file says something about an invalid argument and error code 22. I have run the `spades.py --test`, and everything seems to work fine.

```
Command line: /home/user/SPAdes-3.11.1-Linux/bin/spades.py    -1  /home/user/minion/organism/Illumina/BD_S3_L001_R1_001.fastq.gz -2  /home/user/minion/organism/Illumina/BD_S3_L001_R2_001.fastq.gz -o  /home/user/minion/organism/hybrid_assembly/spades_assembly/read_correction --threads   8   --only-error-correction

System information:
  SPAdes version: 3.11.1
  Python version: 2.7.12
  OS: Linux-4.10.0-40-generic-x86_64-with-Ubuntu-16.04-xenial

Output dir: /home/user/minion/organism/hybrid_assembly/spades_assembly/read_correction
Mode: ONLY read error correction (without assembling)
Debug mode is turned OFF

Dataset parameters:
  Multi-cell mode (you should set '--sc' flag if input data was obtained with MDA (single-cell) technology or --meta flag if processing metagenomic dataset)
  Reads:
    Library number: 1, library type: paired-end
      orientation: fr
      left reads: ['/home/user/minion/organism/Illumina/BD_S3_L001_R1_001.fastq.gz']
      right reads: ['/home/user/minion/organism/Illumina/BD_S3_L001_R2_001.fastq.gz']
      interlaced reads: not specified
      single reads: not specified
Read error correction parameters:
  Iterations: 1
  PHRED offset will be auto-detected
  Corrected reads will be compressed (with gzip)
Other parameters:
  Dir for temp files: /home/user/minion/organism/hybrid_assembly/spades_assembly/read_correction/tmp
  Threads: 8
  Memory limit (in Gb): 87

======= SPAdes pipeline started. Log can be found here: /home/user/minion/organism/hybrid_assembly/spades_assembly/read_correction/spades.log

===== Read error correction started.

== Running read error correction tool: /home/user/SPAdes-3.11.1-Linux/bin/hammer /home/user/minion/organism/hybrid_assembly/spades_assembly/read_correction/corrected/configs/config.info

  0:00:00.000     4M / 4M    INFO    General                 (main.cpp                  :  83)   Starting BayesHammer, built from N/A, git revision N/A
  0:00:00.001     4M / 4M    INFO    General                 (main.cpp                  :  84)   Loading config from /home/user/minion/organism/hybrid_assembly/spades_assembly/read_correction/corrected/configs/config.info
  0:00:00.045     4M / 4M    INFO    General                 (memory_limit.hpp          :  51)   Memory limit set to 87 Gb
  0:00:00.045     4M / 4M    INFO    General                 (main.cpp                  :  93)   Trying to determine PHRED offset
  0:00:00.071     4M / 4M    INFO    General                 (main.cpp                  :  99)   Determined value is 33
  0:00:00.071     4M / 4M    INFO    General                 (hammer_tools.cpp          :  36)   Hamming graph threshold tau=1, k=21, subkmer positions = [ 0 10 ]
  0:00:00.072     4M / 4M    INFO    General                 (main.cpp                  : 120)   Size of aux. kmer data 24 bytes
     === ITERATION 0 begins ===
  0:00:00.100     4M / 4M    INFO   K-mer Index Building     (kmer_index_builder.hpp    : 289)   Building kmer index
  0:00:00.101     4M / 4M    INFO    General                 (kmer_index_builder.hpp    : 106)   Splitting kmer instances into 128 buckets using 8 threads. This might take a while.
  0:00:00.101     4M / 4M    INFO    General                 (file_limit.hpp            :  32)   Open file limit set to 1024
  0:00:00.101     4M / 4M    INFO    General                 (kmer_splitters.hpp        :  74)   Memory available for splitting buffers: 3.62484 Gb
  0:00:00.101     4M / 4M    INFO    General                 (kmer_splitters.hpp        :  82)   Using cell size of 524288
  0:00:00.102     8G / 8G    INFO   K-mer Splitting          (kmer_data.cpp             :  96)   Processing /home/user/minion/organism/Illumina/BD_S3_L001_R1_001.fastq.gz
  0:00:57.279     8G / 8G    INFO   K-mer Splitting          (kmer_data.cpp             : 106)   Processed 1338594 reads
  0:03:30.505     8G / 8G    INFO   K-mer Splitting          (kmer_data.cpp             : 106)   Processed 2152556 reads
  0:03:30.505     8G / 8G    INFO   K-mer Splitting          (kmer_data.cpp             :  96)   Processing /home/user/minion/organism/Illumina/BD_S3_L001_R2_001.fastq.gz
  0:09:23.290     8G / 8G    INFO   K-mer Splitting          (kmer_data.cpp             : 106)   Processed 3442412 reads
  0:14:08.692     8G / 8G    INFO   K-mer Splitting          (kmer_data.cpp             : 106)   Processed 4305112 reads
  0:14:08.692     8G / 8G    INFO   K-mer Splitting          (kmer_data.cpp             : 111)   Total 4305112 reads processed
  0:14:08.726    32M / 8G    INFO    General                 (kmer_index_builder.hpp    : 109)   Starting k-mer counting.
  0:20:36.408    32M / 8G    INFO    General                 (kmer_index_builder.hpp    : 115)   K-mer counting done. There are 205776690 kmers in total.
  0:20:36.409    32M / 8G    INFO    General                 (kmer_index_builder.hpp    : 121)   Merging temporary buckets.
  0:24:43.505    32M / 8G    INFO   K-mer Index Building     (kmer_index_builder.hpp    : 298)   Building perfect hash indices
  0:25:40.125   100M / 8G    INFO    General                 (kmer_index_builder.hpp    : 137)   Merging final buckets.
  1:11:46.553   100M / 8G    INFO   K-mer Index Building     (kmer_index_builder.hpp    : 320)   Index built. Total 95426752 bytes occupied (3.70991 bits per kmer).
  1:11:46.598   100M / 8G    INFO   K-mer Counting           (kmer_data.cpp             : 357)   Arranging kmers in hash map order
  1:12:55.597     3G / 8G    INFO    General                 (main.cpp                  : 155)   Clustering Hamming graph.
  1:21:35.596     3G / 8G    INFO    General                 (main.cpp                  : 162)   Extracting clusters
=== Stack Trace ===
[0x407a6a]
[0x409561]
[0x40c79b]
[0x40129a]
[0x5442b0]
[0x40620d]
Verification of expression '(intptr_t) MappedRegion != -1L' failed in function 'void MMappedWriter::reserve(size_t)'. In file '/spades/src/common/io/kmers/mmapped_writer.hpp' on line 94. Message 'mmap(2) failed. Reason: Invalid argument. Error code: 22'.
Verification of expression '(intptr_t) MappedRegion != -1L' failed in function 'void MMappedWriter::reserve(size_t)'. In file '/spades/src/common/io/kmers/mmapped_writer.hpp' on line 94. Message 'mmap(2) failed. Reason: Invalid argument. Error code: 22'.
hammer: /spades/src/common/io/kmers/mmapped_writer.hpp:94: void MMappedWriter::reserve(size_t): Assertion `(intptr_t) MappedRegion != -1L' failed.

== Error ==  system call for: "['/home/user/SPAdes-3.11.1-Linux/bin/hammer', '/home/user/minion/organism/hybrid_assembly/spades_assembly/read_correction/corrected/configs/config.info']" finished abnormally, err code: -6

In case you have troubles running SPAdes, you can write to spades.support@cab.spbu.ru
or report an issue on our GitHub repository github.com/ablab/spades
Please provide us with params.txt and spades.log files from the output directory.

```

# Answer
*Answer from @mjoppich, converted from comment:*

Late answer, but it looks like you are running out of memory, hence the intptr is invalid.
> 1 votes

---
Tags: fasta, assembly
---


# Question
Title: Estimating computing resources needed for a GWAS?
One of my dissertation papers is going to involve a GWAS. I have never actually run a GWAS before, and do not know how to estimate the computing resources I need for it. I asked someone on my committee who said "it's all guesswork" but I don't even know how to come up with a reasonable guess!

I need to put some estimate of how much cluster time I will need (XX cores for YY hours) in order to do a dissertation grant. I've tried Googling and searching literature for how to do an estimate of this, but I'm not finding much.

Details:

* GWAS, total n=350,000, 10-15% cases
* Already imputed by the time I get it, will have about 96 million variants
* Will be using Plink
* Cluster nodes have 24 cores, and jobs can use up to 4 gigs of RAM per core. If Plink can handle being spread across nodes, I can use up to 100 cores at once.
* Time is billed in core-hours, using 20 cores for one hour is equal to using one core for 20 hours
* We may want to do a few smaller GWAS where we restrict by age (eg. one of all cases and controls under 25, and one of all cases and controls who are 65+)
* Need to come up with an estimate of cluster time needed for each GWAS that we might run (XX cores for YY hours)

Can you please advise me on how to estimate something like this?

# Answer
> 0 votes
*Answer from @llrs, converted from comment:*

I found an article where they said that they could do GWAS in 6000 samples in 6 hours in a regular user computer. But I couldn't find any other explanation about how does it take with other kind of data. It wasn't with plink too and I think it was a bit outdated.

*Answer from @christopher-chang, converted from comment:*

Incidentally, you may want to use plink 2.0 rather than 1.9 for this. It automatically makes good use of all your cores; in contrast, you’d need to start ~24 copies of plink 1.9 to saturate a 24-core machine.

---
Tags: gwas
---


title: How to compare 2 whole genome sequences

text:
# Question
Title: How to compare 2 whole genome sequences
I have the files of 2 persons that got whole genome sequencing, one by nebula genomics and the other by dante labs.

1. I want to compare their genomes to see if they are father and son.
2. I want to see exactly which genetic sequences they share.
3. In the shared sequences, I want to know what mutations exist only in the son.

Which program should I use?

Each company provide fastq and vcf files to the genome. so I have each genome in fastq and vcf. which is better to use?

(I have the files on my hard drive, so is it better to use offline program?)

# Answer
I think you can't do (2) and (3) without the mother.

(1) is doable to a very high accuracy. It is better to start with fastq. You need to generate a multi-sample VCF first. Here is a high-level overview: a) align reads with bwa-mem; b) sort alignments with samtools and then index the resulting BAM; c) call variants *jointly* with freebayes (or GATK if you like challenges). After you get the two-sample VCF, you go through every site where the father is homozygous and it is a *high-quality* site in both samples (e.g. average coverage, high mapping quality, not in difficult regions, etc). Among these sites, check if the child has the paternal allele. If you have a father-son pair, the child is expected to have all the paternal alleles in theory. In practice, if \>99% of paternal homozygotes are present in the child, the two samples will almost certainly be a father-son pair. The actual threshold may depend on how you define "high-quality" sites. Errors may lead to apparent mendelian violations.

In principle, it is possible to do a similar analysis with the two provided VCFs, but you need to be very careful about various technical artifacts.
> 2 votes

# Answer
Assuming you've got genotyped VCF files from the same type of data, and they are similarly normalised, I expect that the VCF files will be easier to work with initially.

That's a lot of assumptions, so until there's more information about the input data, I'll just point you to bcftools merge, which I think would be a reasonable first step given those assumptions:

```
bcftools merge -o merged.vcf.gz A.vcf.gz B.vcf.gz

```

This will produce a merged VCF file that contains genotypes from both files on each line. To look for similarities and/or differences, you would then inspect each line to check whether the genotypes were the same or different (e.g. `same: 1/1 1/1`; `different: 0/1 1/1`).
> 0 votes

---
Tags: wgs, human-genome
---


title: What do the FASTQ file names mean here?

text:
# Question
Title: What do the FASTQ file names mean here?
I would like to run an analysis on Nanopore data. I'm trying to download sample data from https://github.com/nanopore-wgs-consortium/chm13. I would like to run minimap2 on the FASTQ files, resulting in a single BAM file.

Q: What do the file names mean? What is "S13"? What is "L002\_l1"? What is "L002\_R1"? Am I supposed to download them and combine them into a single FASTQ file, and then align them?

# Answer
*What is "S13"?*

S# is most likely the sample number. If so, you do not want to combine these.

*What is "L002\_l1"?*

L00# is the lane number (of the 8 Illumina flow cell lanes). 'I1\_' (not 'L') are index files and are not needed (see also https://www.biostars.org/p/374581/).

*What is "L002\_R1"?*

R1 is one end of paired-end sequencing, R2 is the other end. These filepairs must have the same number of lines. Use an aligner that utilizes both.
> 5 votes

# Answer
Peter's answer is a good one and to the point. I'll be a bit more verbose.

The breakdown of FASTQ file names that come directly from the sequencer typically have the following format:

`{sample_name}_S{sample number}_L{lane number}_{R/I}{read or index number}_001.fastq.gz`

`sample_name` is the sample name provided by you (or whoever sequenced the data) to the sequencer.

To save time and money, samples are often "multiplexed" on sequencing runs. This multiplexing means you aggregate samples ("S") on the same run and split samples on different lanes ("L") to limit some technical noise. Lanes are essentially technical replicates for the sequencing data on the same sequencing run, and these numbers tell you which sample and which lane that FASTQ comes from.

Sequencing can be done in single-end or paired-end. Single-end sequencing generates a single FASTQ file (`_R1`), whereas paired-end sequencing produces 2 FASTQ files (`_R1` and `_R2`).

The data you highlight isn't actually Nanopore:

> Approximately 50x of data was generated on a NovaSeq instrument.

10x includes "indexes" ("I") in addition to the sample and lane information to keep track of things like single cells or other unique molecular identifiers, and that information is listed after the `I`. `_R1`, `_R2`, and `_I1` all go together for a unique sample and lane combination
> 1 votes

# Answer
The file names mean this isn't Nanopore data. Read the bit just above these file names:

> Approximately 50x of data was generated on a NovaSeq instrument. Based on the summary output of Supernova, there are 1.2 billion reads with 41x effective coverage.

That linked github page has Nanopore, Illumina/10x, BioNano, HiC, and PacBio data. For the nanopore files, you should be looking at this section.
> 0 votes

---
Tags: fastq, genomics
---


# Question
Title: How to align Illumina DNA seq reads to CDS?
I have the CDS sequences for my genome. I want to align Illumina DNA seq reads from related fungal strains to these. This is to identify SNPs in just the coding sequences. I am not sure if I should use a DNA seq aligner like BWA and Bowtie, or to use a split aware aligner like STAR.

I have the full genome. But I wanted the SNPs for each of the CDS sequences, so I can create consensus sequences for these CDS for the related strains.

Can you provide examples or links to get alignment from CDS region or using VCFtools to get SNPs and INDELs from the CDS region? I have tried looking but haven't found anything useful. Thank you.

# Answer
*Answer from @fabio-marroni, converted from comment:*

I suggest that you align to the full genome and after this, you only extract the SNPs in the CDS. This procedure is commonly used in all exome sequencing projects, so you can find a lot of instructions software. You can choose between extracting the regions from the alignment (using samtools), or keeping everything, calling the SNPs and then only select those in CDS using vcftools.

For looking at the CDS region, you could try bedtools, as suggested here.
> 1 votes

---
Tags: sequence-alignment
---


# Question
Title: How should I deal with segmental duplications when aligning NGS reads to a reference genome?
This is a follow-up of my other question. I have been having trouble calling variants in the human SMN1 and SMN2 genes, because the human genome has a large segmental duplication there and these two genes (and their surrounding sequences) are essentially identical :

Click here to see the original image on the UCSC browser.

Since variants in the SMN1 gene have been linked to disease, it is important for me to be able to call them in my NGS data. However, because of this segmental duplication, aligning reads to this region correctly is very hard since the presence of the duplication makes each read align almost perfectly to separate locations in the genome, leading to low alignment scores. The low scores will, in turn, block variant callers from calling in those regions.

I thought that the only real solution for this would be to mask the SMN2 gene and its associated duplication in my reference genome which would allow the aligner to align all reads to the region of SMN1. I would then miss any real variants in SMN2, but since the primary gene of interest for spinal muscular atrophy is SMN1 and not SMN2, I would be OK with that.

My questions are:

1. Is this a reasonable approach? I was told that "generally, masking segdups is not recommended" by someone who most certainly knows what he's talking about. So is this not a good approach for my stated goal?
2. If not, then what else can I do that will let me call variants in the SMN1 gene?

I would like something similar to the approach described here: mask one of the regions to rescue mutations in the other.

# Answer
Short answer for short reads: you can't. Not reliably, at least.

There's a known error mode with aligning short reads to reference genomes where it's quite good at mis-calling reads within repetitive regions. See, for example, this paper, which compares Illumina vs Nanopore for detecting bacterial strains.

If you know that a region is likely to contain a segmental duplication and/or deletion (e.g. due to unexpected changes in read coverage), then it'd be a good idea to double-check any discovered variants (or unexpected non-variants) using a different technology that does not have the same coverage issues. I recommend a long-read sequencing technology that can handle regions at least as large as the segmental duplication of interest (e.g. PacBio up to ~15kb, Nanopore for anything longer), so that you have a high chance of mapping reads uniquely.
> 1 votes

---
Tags: sequence-alignment, variant-calling, reference-genome
---


# Question
Title: What is protein secondary structure assignment?
I am absolutely new to Structural Bioinformatics (I only started last week). I am working on the secondary structure assignment/prediction (actually I am not sure) of proteins using ML.

When it comes to protein analysis, most of the books or research papers talk about protein structure **Prediction**:

https://bip.weizmann.ac.il/course/structbioinfo/pdf/SecondaryStructurePrediction\_WIS.pdf

Only a handful of books or articles talk about protein structure **Assignments**:

https://www.biorxiv.org/content/10.1101/2021.02.02.429329v1

Can anyone tell me what is the basic difference between Assignment and Prediction?

Kindly explain with examples.

# Answer
To make it crystal clear (and make than pun):

* `Assignment`: a 3D structure is known and the residues are *assigned* a secondary structure
* `Prediction`: an algorithm predicts from linear (primary) sequence what the secondary structure may be —often incorrectly. This can be homology based, covariance, deep learning etc. But most commonly just the vector norm of window residue composition, e.g. proline is never in a helix and helices and sheets have different periodicity (two for a sheet and three and a bit for a helix).
> 7 votes

# Answer
*Answer from @jan-wilamowski, converted from comment:*

That second paper seems to distinguish between the known data that is used for ML training (i.e. amino acid sequences are assigned some secondary structure) and the attempt to replicate this natural assignment by an "artificial" prediction.
> 1 votes

---
Tags: proteins, protein-structure, secondary-structure
---


# Question
Title: Precisely clipping bam file to bed coordinates
I have a bed file of coordinates and would like to clip anything outside these regions in my aligned bam files. However, trying programs such as bedtools intersect and samtools ampliconclip (on complementary bed file), they always let some reads slip through outside the ranges. I tried using hard-clip and both-ends options for ampliconclip aswell.

Example (top is clipped bam, middle is unclipped bam and bottom is bed file with regions I would like to keep):

Although most reads outside the designed amplicons are removed, there are still reads outside the regions in the Designed\_pool1 file. The problem is that the bed file contains amplicons and outside the amplicons are primers where I do not want to call variants. I cannot get the exact primer sequences because I don't know their length and the manufacturer does not state them, so I cannot use programs that require primer sequences.

Does anyone have some experience with this and know why the programs behave this way?

EDIT:

The exact commands I have used are:

```
bedtools intersect -a $BamFile -b $BedFile > $BamOut

```

and

```
samtools ampliconclip --hard-clip --both-ends -b $ComplBedFile $BamFile > $BamOut

```

as well as without --hard-clip and --both-ends options. All methods let some sequences stay out of the regions in bed files.

I am not completely familiar with the bedpe format, but I am not using paired end data. I cannot get primer start and stop positions, only amplicon positions.

# Answer
I am curious whether your pileup visualization tool is respecting the clipping. It is possible that there is a subtlety of samtools ampliconclip that your visualization is not understanding, with regard to how it treats e.g. supplementary alignments or something.

Perhaps more importantly, **I notice that the reads that are not being properly clipped in your visualization tend to have one end anchored in each of your BED regions.**

Looking at the samtools ampliconclip docs, it seems to be assuming that **reads will only be clipped from the ends**. They will not be clipped if the region outside the BED is in the middle, if I read it correctly. I think that you may need to include an additional step in which you either trim reads (so that they are not long enough to bridge your regions) or you somehow postprocess your alignments to explicitly clip in the middle.

If true, this seems like a significant enough oversight of samtools that it may be worth logging an issue with the samtools maintenance team.

Also, not clear if it's relevant but reading the samtools ampliconclip docs:

> adjustments to the left most mapping position (POS) will mean that **coordinate sorted files will need resorting. In such cases the sorting order in the header is set to unknown. Clipping of reads results in template length (TLEN) being incorrect.** This can be corrected by samtools fixmates. Any *MD and NM aux tags will also be incorrect*\*, which can be fixed by samtools calmd.
> 1 votes

# Answer
For single-end reads, generically:

```
$ bam2bed < reads.bam > reads.bed
$ bedmap --echo --fraction-ref 1 reads.bed roi.bed > answer.bed

```

The BED file `answer.bed` contains filtered reads which entirely overlap intervals in `roi.bed`. It contains fields required for backconversion to SAM/BAM.
> 1 votes

---
Tags: base-clipping
---


# Question
Title: How to run GSEA analysis on R studio using DEG file list generated from scRNAseq analysis
I want to draw the standard plot of GSEA on Rstudio. I have a data frame that consists of a list of DEGs as below.

Actually, I put the DEG list directly to PANTHER that resulted in the GO or pathway analysis data. However, I want to use GSEA.

Meanwhile, I found nice online site and follow their codes for the GSEA analysis(https://bioinformaticsbreakdown.com/how-to-gsea/).

And, here is the my code.

```
tmp = read.csv("/data/directory/DEGs_HCvsJSS_DC.csv", header = TRUE) %>% filter(Gene.Symbol !="")

GO_file = "/blue/zhangw/my.kim/HCvsJSS/c5.all.v7.4.symbols.gmt" # got this from http://www.gsea-msigdb.org/gsea/msigdb/collections.jsp

gene_list = tmp$avg_log2FC

names(gene_list) = tmp$Gene.Symbol

gene_list = sort(gene_list, decreasing = TRUE)

gene_list = gene_list[!duplicated(names(gene_list))]

head(gene_list)

res = GSEA(gene_list, GO_file, pval = 0.05) # this caused the error as follows. Error in build_Anno(TERM2GENE, TERM2NAME) :  argument "TERM2GENE" is missing, with no default

```

It's really dissapointing because I found this relatively easy protocol for GSEA, but found the error, lol.

Do you guys have any suggests for this troubles? or any recommendations for GSEA analysis that can be done with DEG list generated from seruat scRNAseq data?

# Answer
I'd probably just use fgsea. Your code won't have to change much, it also just uses a ranked list.
> 1 votes

---
Tags: r, scrnaseq, gsea
---


# Question
Title: How to demultiplex MinION fastQ files without Porechop?
I recently got fastQ files from a MinION sequencer. I would like to demultiplex these files.

Is there any other software to demultiplex (no porechop) in or outside Artic protocol? I would like to learn options different from porechop.

Thank you!

# Answer
The Guppy basecaller itself now has the ability to internally perform multiplexing. That is currently the preferred option.

If for some reason, you still want a stand-alone application to demultiplex already basecalled reads, then `qcat` is an option: https://github.com/nanoporetech/qcat
> 2 votes

# Answer
I have written up my own semi-manual method for demultiplexing Nanopore reads. As distinct from the other methods that are commonly used, it searches the *entire* read for potential barcode sequences, and splits out identified chimeric reads (containing multiple barcode sequences) into a separate file:

https://dx.doi.org/10.17504/protocols.io.7vmhn46
> 1 votes

---
Tags: r, minion
---


# Question
Title: splitCsv map with text, join with another channel, nextflow
I have a csv file which describes number of chunks per chromosome, like;

```
chrom,chunks
chr1,"CHUNKS=([1]=""0e6 3e6"" [2]=....)"
chr2,"CHUNKS=([1]=""0e6 3e6"" [2]=....)"
chr3,"CHUNKS=([1]=""0e6 3e6"" [2]=....)"
.....

```

I could crate a channel as follows,

```
channel
    .fromPath('chunks3.csv')
    .splitCsv(header:true)
    .map{row -> tuple(row.chrom,row.chunks)}
    .set{chunk_ch}

```

It gives me chromosome with respective chunk definition.

But when I try to join this channel with another channel, it does not work.

For example I have a channel "phased". It is a output channel created from a process.

```
channel
    .fromFilePairs('/some/path/genetic_map_chr*_combined_b37.txt',size:1)
    .map {group_key, file_list -> tuple(group_key.replaceFirst(/^*genetic_map_/,""),file_list.first())}
    .set {genetic_map}

channel
    .fromFilePairs('/some/path/*hrc.r1-1.ega.grch37.chr*.{hap.gz,legend.gz,samples}',size:3)
    .map {group_key,file_list -> tuple(group_key.replaceAll(/^.*.hrc\.r1-1\.ega\.grch37\.|\.haplotypes\.noibd|\.haplotypes/,""),file_list)}
    .set {ref_panel}

qc1_chr
    .map{chrom,bfiles -> tuple(chrom,bfiles.findAll{it.baseName.endsWith('step10')})}
    .join(ref_panel).join(genetic_map)
    .into {shape_in;imp_in}

process phasing {
    publishDir params.out, mode:'copy'

    input:
    tuple val(chrom), path(input_bed), path(input_ref), path(input_map) from shape_in

    output:
    tuple val("${chrom}"),path("${chrom}.*") into phased

    script:
    def (bed, bim, fam) = input_bed
    def (haplotype, legend, sample) = input_ref
    """
    shapeit --thread "${task.cpus}" --input-bed "${bed}" "${bim}" "${fam}" --input-map "${input_map}" --input-ref "${haplotype}" "${legend}" "${sample}" --duohmm --
    output-max "${chrom}.step10" --seed 54321
    """
}

```

Which emits something like:

```
phased.view()
[chr1,[file1]
chr2,[file2]
chr3,[file3]....]

```

If I try `chunk_ch.join(phased)` it gives nothing.

Any help?

# Answer
Below is a reproducible example. I tried to simulate the same channels using `chr1` as an example.

```
#!/usr/bin/env nextflow

// Simulate the chunk channel
channel
    .fromPath('chunk.csv')
    .splitCsv(header:true)
    .map{row -> tuple(row.chrom,row.chunks)}
    .set{chunk_ch}

// Simulate chr and chr_path as a single entry list
channel.from(["chr1", "chr1.txt"]).toList().map{tuple(it[0], it[1])}.set{shape_in}

// uncomment to view the content
// shape_in.view()

// process of phased_ch. I modified it to output the chr and chr_path as given in your example
process phased_ch {
  echo true
  input:
  tuple val(chr), val(chr_path) from shape_in
  output:
  tuple val("${chr}"), path("${chr_path}") into phased
  script:
  """
  echo ${chr}
  touch ${chr_path}
  """
}

// joining the two channels to produce the correct join
chunk_ch.join(phased).view()

// will produce
// [chr1, "CHUNKS=([1]=""0e6 3e6"" [2]=....)", /media/Documents/Nextflow/work/8e/2fcd9a00b04f0f5d44711068d68a0f/chr1.txt]

```

`chunk.csv` contains

```
chrom,chunks
chr1,"CHUNKS=([1]=""0e6 3e6"" [2]=....)"
chr2,"CHUNKS=([1]=""0e6 3e6"" [2]=....)"
chr3,"CHUNKS=([1]=""0e6 3e6"" [2]=....)"

```

In your case, you should verify what are the contents of `chunk_ch` and `phased` on output. The first key of `join()` should be similar otherwise it won't give anything. From the documentation of join:

> The join operator creates a channel that joins together the items emitted by two channels for which exists a matching key. The key is defined, by default, as the first element in each item emitted.
> 1 votes

---
Tags: nextflow
---


# Question
Title: How to calculate RMSD and binding energy from the already docked ligand-receptor
I have the following PDB file that looks like this. It is downloadable **here**. This PDB already contains ligand and receptor together.

The red molecule is the ligand and the yellow is the receptor. They're already docked with **ClusPro**. My question is how to calculate the RMSD and binding energy from it? Any existing tools (with or without command line) that can help?

I don't have any data to try it out at hand but have a look at the documentation I linked. There you can find ways of constructing an Ensemble which can then be superposed and RMSDs calculated, even pairwise.

Reference is the red one (receptor). I have hundreds of other complexed ligands. I want to calculate how good are those complexes. What's the reasonable measure to do that?

Example code using Prody? please :-)

# Answer
> 1 votes
> NB. This is a long comment as opposed to an answer which requires a colossal working chunk of code.

## Control

Science is all about controls failing. A great control for docking is the barnase-barnstar complex (PDB:1BRS): both small, but form a strong interactions. If this is done, the RMSD against the original structure is very meaningful and useful.

## RMSF

A RMSD is a metric of how much two structures differ. In an ensemble, the metric of how much a structure and the ensemble average is called a RMSF (root-mean square fluctuation). Here *the RMSF is sought* as the RMSD from first/top pose or worse the manually thrown together pre-docking pose does not mean much.

## Clusters

The decoys (PBDs trails) need to be clustered into different ensembles first, before a RMSF can be calculated.

## ∆∆G and water

For the calculation of difference in Gibbs free energy between the docked minus the docked poses, the structures need to be energy minimised with the forcefield scorefunction that is to be used for scoring. Water is very important for protein-protein interactions, but adds a lot of complication. Say Gromacs is used, `gmx energy` will return the single state (conformer) Gibbs free energy, but it requires a water box. For the difference, the water box of undocked needs to be the same size and preferably with the same amount of molecules. A much simpler solution is using an implicit solvent model, e.g. Rosetta. As a bonus, the energy minimisation (`relax`) is relative fast. And specific explicit waters can be added. However, calculating ∆∆G for every decoy will take hundreds of CPU-hours and is a bit nonsensical.

A better approach is to look at residue conservation (consurf), cross-linking mass spec data if available, and if human post-translational modifications (phosphosite) and healthy variants (gnomAD). But ultimately wet lab validation by ITC of different mutants/labelling based upon the docking results *etc.* is required.

---
Tags: proteins, protein-structure, pdb, docking
---


# Question
Title: How to reconstruct a target string from a large set of overlapping reads
I got this question as a homework. Does anyone know how the it can be solved?

> Consider a short string $r$ as a read from a string $S$ if $r$ is converted from a substring of $S$ by at most $k$ swaps at random locations, where a swap is an operation that switch two adjacent characters, i.e., string $s\_1s\_2...s\_i s\_{i+1}...s\_n$ is one swap away from the string $s\_1s\_2...s\_{i+1} s\_{i}...s\_n$. Two reads are called overlapping if the substrings to which they are converted are overlapping. Given the input of $k \<\< |r|$, devise an algorithm to reconstruct a target string $S$ from a large set of overlapping reads of $S$.

The question doesn't explicitly mention which of the possible superstrings of reads should be returned, so I think it's asking for *any* superstring.

# Answer
> 2 votes
Unless it's not obvious, the problem has no unique solution. An obvious example is a string which uses entirely one character; the problem is then to determine its *length*, which you can only do statistically.

As you've correctly surmised, this is more or less the short-read *de novo* assembly problem, with one modification: the "read errors" are all swaps.

The usual way that we solve this problem is:

* Pick a length (which is customarily called $k$ in bioinformatics, but the problem confusingly already used that, so we'll use $\kappa$ instead), and extract all $\kappa$-length substrings of the reads. These substrings will, from here on, be referred to as $\kappa$-mers.
* Construct a subgraph of the De Bruijn graph, where the $\kappa$-mers are the nodes, and there is an edge between two nodes if two adjacent $\kappa$-mers are connected in any read. Label the edges with the number of times this happens in all reads.
* Prune the graph of any edges that have low frequency. In your example, the kinds of errors are known, so you could in fact "correct" the pruned edges, or even correct the reads, based on the higher-consensus path.
* An Eulerian superpath through the graph corresponds to the desired string.

I'm not saying that this is the best approach for your homework question, but this is how we do it in practice.

To understand the method in depth, here are some suitable references. Some of these are tool papers, and others explain the data structures and algorithms well.

---
Tags: algorithms, homework
---


# Question
Title: How to show that one docking is better than another with Pymol visually
I have the following docking results:

`Dock_1` has better energy score than `Dock_2`. What I want to do is to show with Pymol that `Dock_1` is better visually.

Here I use hydrogen bond, highlighted with Yellow (`Dock_1`) and Gray (`Dock_2`). But as you can see, with that I can't show that `Dock_1` is clearly better. That's why I'm looking for better alternative than hydrogen bond.

# Answer
> 1 votes
maybe you can use the "**per-residue interaction scores**" (wich is available via Schrödinger).

On the other hand,to show if a Docking Pose is really "better" you could show specific interactions rather than the complete big molecule.. e.g. an important salt bridge is missing. (hydrogen bonds are not that important, due to the fact that they can be bridged with water molecules)

Sorry, this is not really a solution, but I am not allowed to comment jet -.-

---
Tags: protein-structure, pdb, docking, pymol
---


title: color by clusters and sampled in Seurat

text:
# Question
Title: color by clusters and sampled in Seurat
I have a Seurat object

I want to have both cluster numbers and coloured cells by sample names like this figure (from a Nature paper)

I have tried `group.by` argument in Seurat but only samples come not cluster numbers like below

```
> DimPlot(pbmc, reduction = "umap",group.by = "samples")
>

```

\[!\[enter image description here\]\[2\]\]\[2\]

Any help that to get something like the first figure? Having cluster numbers and coloured by sample names please

# Answer
> 2 votes
Here is a solution that makes use of `LabelClusters()` from `Seurat`:

```
# creating a plot and assigning it to the "plot" variable
# coloring by sample id and turning labels off
plot <-
DimPlot(
  pbmc,
  reduction = "umap",
  group.by = "some_random_sample_id",
  label = FALSE,
)

# here is the trick, adding cluster labels to the "data" variable of the "plot"
# without this, the "plot" object would not have the clustering information but only embeddings and sample id
plot$data$seurat_clusters<- pbmc@meta.data$seurat_clusters

# and here is a sanity check to see if cluster labels are transferred to the "right places"
all(rownames(plot$data) == rownames(pbmc@meta.data))

LabelClusters(plot, id = "seurat_clusters")

```

---
Tags: scrnaseq, seurat, 10x-genomics
---


title: Python module bamtags

text:
# Question
Title: Python module bamtags
Does anyone know this python module, bamtags? I'm using a pipeline (Dovetail Genomics, HiRise) that imports it, but it seems it doesn't exist.

It's impossible to install it and there is no information on the internet. Does anyone know it?

Additional information:

Hirise is a 3 years old pipeline that uses Chromossome Capture (HiC) data to identify misassembled regions in a genome assembly and then rescaffold them. By integrating 2 different sequencing methods, the confidence of the genome assembly is thus increased.

The workflow is free (we only paid for the data generation). Unfortunately I have no idea what the module does, as I find no information on internet and the scripts are poorly documented. f python-reading is correct, it looks like that is an internally defined module (class) which is in here. So basically, it's there. I am not a python expert so maybe others more experience with python development will confirm or deny this.

The only function that is used is the following: `BamTags.mate_mapq(aln)`, which is inside a loop in aln: `for aln in sam.fetch(reference=c1)`

Thanks, Ricardo

# Answer
> 0 votes
*Answer from @fridaymeetssunday, converted from comments:*

If python-reading is correct, it looks like that is an internally defined module (class) which is in here. So basically, it's there.

*Response from poster:*

Thank you, guys, so helpful! Yes, that's it. I wasn't aware of the simplicity of import function. I thought it only worked with builtin or pip/conda installed modules. I just had to add the scripts folder to my PATH variable or have the functions in the working directory.

---
Tags: python, sam, genomics, hi-c
---


title: Adding metadata to Seurat object

text:
# Question
Title: Adding metadata to Seurat object
I have a Seurat object of 8 patients. I want to add metadata to that so that I have origin of each cell. At the moment UMAP just shows a bunch of cells while I want to color clusters by sample

This is my metadata

```
> str(head(metadata))
'data.frame':   6 obs. of  3 variables:
 $ barcodes: chr  "AAACAGCCAAACTGCC-5" "AAACAGCCAACAGGAT-1" "AAACAGCCAAGATTCT-1" "AAACAGCCAAGGTATA-1" ...
 $ libcodes: Factor w/ 8 levels "1","2","3","4",..: 5 1 1 1 5 7
 $ samples : chr  "PN0252_0005" "PN0252_0001" "PN0252_0001" "PN0252_0001" ...
>

> pbmc
An object of class Seurat
36601 features across 26473 samples within 1 assay
Active assay: RNA (36601 features, 0 variable features)
>

> dim(metadata)
[1] 26473     3
>

```

# Answer
> 2 votes
If you have single-dimension per-cell metadata, and it's arranged identically to the cell order in the Seurat object, I find it easier to use the double bracket notation to add metadata to a Seurat object. For example:

```
metadata$barcodes -> pbmc[["barcodes"]]
metadata$libcodes -> pbmc[["libcodes"]]
metadata$samples -> pbmc[["samples"]]

```

# Answer
> 2 votes
*Answer from @user438383, converted from comment:*

Did you try the AddMetaData function? This post may also be useful.

```
cluster_letters <- LETTERS[Idents(object = pbmc_small)]
names(cluster_letters) <- colnames(x = pbmc_small)
pbmc_small <- AddMetaData(
  object = pbmc_small,
  metadata = cluster_letters,
  col.name = 'letter.idents'
)

```

---
Tags: scrnaseq, seurat, 10x-genomics
---


# Question
Title: Getting all variants in a region from ensembl
I want to get all the variants (SNPs) in a given region of the human genome. Is there a way on how to do this with the ensembl-API? It should be a representation similar to this:

There I added dbSNP to the region details. Having this data in a JSON or XML-Format would be really helpful! Here's the link to the image.

Edit: I know that the export-function on ensembl is a thing. But then I only get the rs-Numbers. To get all the info I need I would have to run (potentially thousands) all the rs-Numbers through some additional API. Getting that info in bulk right away would be a huge time-saver.

# Answer
> 1 votes
REST API overlap region endpoint. For example.

---
Tags: ensembl, variation, api
---


title: How to convert GFF3 to GTF2

text:
# Question
Title: How to convert GFF3 to GTF2
I would like to convert a file in `gff3` format to a `gtf2.2` format.

The reason why I would like to do this is: I have a set of transcripts assembled by a bunch of different software (and using RNA-seq data from different sequencing technologies) and I would like to compare each of those sets with the reference transcriptome annotation for that species (*D. melanogaster*). I already asked for the community advice about how to proceed on that, but to run the suggested software (TACO and cuffmerge) I need to have a GTF file containing the experimental transcripts to compare to the reference transcripts (GTF).

Up to now, I **unsuccessfully** tried:

1. `gffread` utility in the Cufflinks package (`gffread input.gff3 -T -o output.gtf`): this results in an empty `output.gtf` file and an empty `log` file (used `Cufflinks` v.2.2.1) - I contacted the authors via their Google group but haven't heard of them yet
2. `gff3_to_gtf` utility in the GenomeTools (`gt`) package (`gt gff3_to_gtf input.gff3 -o output.gtf`): the output is not created and the log file is not informative - I contacted the authors via their mailing list but haven't heard of them yet
3. `GFF3_to_GTF` utility in the FML package (`./gff3_to_gtf_converter.pl input.gff3 output.gtf`): the output just contains a header (`##gff-version 2.5`) and the `log` is empty

The `gff3` file was created as output of GMAP, and contains the transcripts as found by alignment to the reference (specifying option `-f gff3_match_cdna`).

\[Edit: what I found out a poteriori, is that this format is not a standard gff3, thus the conversion is not trivial...\]

Here is the beginning of the gff3 file I tried to convert:

```
$ head r9_gmap_match-cdna.gff
##gff-version   3
# Generated by GMAP version 2017-04-24 using call:  /home/aechchik/software/gmap-2017-04-24/bin/gmap.sse42 -d gmapidx -D /scratch/beegfs/monthly/aechchik/isoforms/ref/chromosomes/ -f gff3_match_cdna -n 0 -t 20 r9_2d.fasta
2L    gmapidx    cDNA_match    18442664    18443024    79    -    . ID=ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1;Name=ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d;Target=ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d 141 489;Gap=M13 I10 M8 D3 M8 D3 M32 D1 M4 D2 M5 D3 M34 I3 M6 D1 M7 D1 M5 D1 M30 D2 M16 I1 M8 D3 M10 D1 M5 D1 M6 I1 M8 I1 M8 D1 M33 D3 M33 I1 M28 D3 M25
###
3R    gmapidx    cDNA_match    15853880    15855465    96    +    . ID=dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1;Name=dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d;Target=dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d 80 1645;Gap=M46 D2 M12 D2 M157 D3 M4 D1 M50 I2 M3 I1 M12 I1 M68 I1 M66 I1 M53 D2 M47 D1 M16 D1 M35 D1 M155 D1 M28 D1 M166 D2 M47 D1 M8 D4 M69 D1 M28 D1 M5 D1 M12 D1 M202 I1 M115 D1 M61 I2 M7 D1 M7 I1 M36 D2 M41
3R    gmapidx    cDNA_match    15855529    15855742    97    +    . ID=dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1;Name=dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d;Target=dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d 1646 1856;Gap=M21 D1 M68 D1 M85 D1 M37
###
X    gmapidx    cDNA_match    14837810    14838142    93    -    . ID=960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1;Name=960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d;Target=960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d 74 406;Gap=M13 D1 M182 I1 M14 I2 M56 I1 M30 D2 M21 D1 M13
X    gmapidx    cDNA_match    14837470    14837753    92    -    . ID=960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1;Name=960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d;Target=960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d 407 688;Gap=I2 M5 I1 M64 I1 M44 D1 M9 D5 M31 D3 M23 I1 M19 I1 M25 I1 M26 D1 M19 I1 M9
###

```

# Answer
> 5 votes
The following bit of python code should work:

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

lastTranscript = [None, None, None, []]  # ID, chrom, strand, [(start, end, score), ...]

def getID(s):
    """Parse out the ID attribute"""
    s = s.split(";")
    for k in s:
        if k.startswith("ID="):
            return k[3:]
    return None

def dumpLastTranscript():
    """Print the last transcript"""
    bounds = sorted(lastTranscript[3])
    print("{}\tgmapidx\tgene\t{}\t{}\t.\t{}\t.\tgene_id \"{}\"; transcript_id \"{}\";".format(lastTranscript[1], bounds[0][0], bounds[-1][1], lastTranscript[2], lastTranscript[0], lastTranscript[0]))
    print("{}\tgmapidx\ttranscript\t{}\t{}\t.\t{}\t.\tgene_id \"{}\"; transcript_id \"{}\";".format(lastTranscript[1], bounds[0][0], bounds[-1][1], lastTranscript[2], lastTranscript[0], lastTranscript[0]))
    for start, end, score in bounds:
        print("{}\tgmapidx\texon\t{}\t{}\t{}\t{}\t.\tgene_id \"{}\"; transcript_id \"{}\";".format(lastTranscript[1], start, end, score, lastTranscript[2], lastTranscript[0], lastTranscript[0]))

def handleLine(cols):
    """Handle a single line, appending the exon bounds to the previous if relevant"""
    ID = getID(cols[8])
    assert(ID is not None)
    if lastTranscript[0] is not None and lastTranscript[0] != ID:
        dumpLastTranscript()
        lastTranscript[3] = []
    lastTranscript[0] = ID
    lastTranscript[1] = cols[0]
    lastTranscript[2] = cols[6]
    lastTranscript[3].append((int(cols[3]), int(cols[4]), cols[5]))

f = open(sys.argv[1])
for line in f:
    if line.startswith("#"):
        continue
    cols = line.strip().split("\t")
    handleLine(cols)

dumpLastTranscript()
f.close()

```

If you saved that as `gff2gtf.py` and made it executable, the usage would be `./gff2gtf.py foo.gff3 > foo.gtf`. With the example that you provided in your post, the result is:

```
2L  gmapidx gene    18442664    18443024    .   -   .   gene_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1"; transcript_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1";
2L  gmapidx transcript  18442664    18443024    .   -   .   gene_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1"; transcript_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1";
2L  gmapidx exon    18442664    18443024    79  -   .   gene_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1"; transcript_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1";
3R  gmapidx gene    15853880    15855742    .   +   .   gene_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1"; transcript_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1";
3R  gmapidx transcript  15853880    15855742    .   +   .   gene_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1"; transcript_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1";
3R  gmapidx exon    15853880    15855465    96  +   .   gene_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1"; transcript_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1";
3R  gmapidx exon    15855529    15855742    97  +   .   gene_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1"; transcript_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1";
X   gmapidx gene    14837470    14838142    .   -   .   gene_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1"; transcript_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1";
X   gmapidx transcript  14837470    14838142    .   -   .   gene_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1"; transcript_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1";
X   gmapidx exon    14837470    14837753    92  -   .   gene_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1"; transcript_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1";
X   gmapidx exon    14837810    14838142    93  -   .   gene_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1"; transcript_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1";

```

That at least looks like correct gtf 2.2 format to me.

# Answer
> 1 votes
The GenomeTools package has lots of small applications that let you massage GFF3. It may have a converter, or it may let you tweak your existing GFF3 so that gffread works

# Answer
> 0 votes
I'm quite sure it will work with `agat_convert_sp_gff2gtf.pl` from AGAT.
See the review of the different solution I made here: https://agat.readthedocs.io/en/latest/gff\_to\_gtf.html

AGAT it the only one you can decide wich version of GTF you require as output.

---
Tags: software-recommendation, format-conversion, gff3, gtf
---


title: Genome annotation

text:
# Question
Title: Genome annotation
I have helped a lab sequence a mitochondrial genome. I used then MITOS to annotate the genome (warning them that I did it just for curiosity as I am not experienced). They submitted the sequence and annotation to genbank but it was refused for different reasons.

It's the genome of a species for which there is no reference. It's a jumping spider and the genome is almost identical to the one below. I used MITOS web server and I got an annotation almost identical to this one (http://mitos.bioinf.uni-leipzig.de/result.py?hash=bJuyJw4b) which corresponds to this sequence (https://www.ncbi.nlm.nih.gov/nuccore/MH891570 ). But the submitted annotation to Genbank is slightly different so I guess they manually edited it or they used another software. I am not sure I can share more specific information without the permission of the lab.

You can download a GFF file or the raw sequence at the page in the link (for the similar sequence).

My questions are as follows:

* after applying MITOS, is there any manual or automatic review that should be applied?
* Is it common that the annotation returned by MITOS is not good enough?
* Or could that indicate that maybe there is a problem in the assembled genome? (looks unlikely except for a few SNPs as three assemblers return almost the same sequence)
* Should we edit the annotation ourselves or there is a software that does it itself?
* How often using MITOS is enough or how often manually editing is needed?

The sequence corresponds to an arachnid. The reasons given by Genbank are:

> Some or all of the sequences contain incorrect annotation and/or contain little similarity to other sequences in the database based on BLAST searches. Please also review your nucleotide sequences for issues such as low quality data or frameshifts.
>
> Incorrect annotation is:
>
> * coding regions lack stop codons
> * CDS features do not have the correct nucleotide spans or location for example, nad5 is longer at the amino terminus compared to similar proteins in the database and extends into the flanking tRNA.

The first part of their reply is apparently generic as I have seen the same answer in another submission in which I didn't participate. I am not sure it's even verified as there are quite good BLAST matches. I have no idea about the other part. I just applied MITOS.

I used the invertebrate mitochondrial translation table option in MITOS. I guess it's the correct one as the result quite makes sense except for small adjustments.

# Answer
Any regions annotated as protein sequences must have a complete protein sequence, including start and stop codons. It's somewhat common for mitochondrial genomes to require polyA extension for the stop codons to be completed; I expect that NCBI is asking that the submissions make this explicit. See, for example, our lab's Nippostrongylus brasiliensis mtDNA submission:

https://www.ncbi.nlm.nih.gov/nuccore/NC\_033886.1

As an example in this genome annotation, the ND5 gene requires extension, as mentioned in the note:

```
3858..5439
/gene="ND5"
/locus_tag="B2K98_mgp10"
/note="TAA stop codon is completed by the addition of 3' A
residues to the mRNA"
/codon_start=1
/transl_except=(pos:5439,aa:TERM)
/transl_table=5
/product="NADH dehydrogenase subunit 5"
/protein_id="YP_009347353.1"
/db_xref="GeneID:31081492"

```

It doesn't surprise me that an automated tool is not perfect. The mitochondrial genome is small enough that manual verification can be feasibly done. You should at least look at the annotated information to make sure it makes sense; the response from NCBI suggests that it doesn't (e.g. nad5 extends into an annotated tRNA sequence).

As a direct / explicit response to your questions, just in case my answers are useful:

> after applying MITOS, is there any manual or automatic review that should be applied?

I don't know. You should at least look at the annotated information to make sure it makes sense.

> Is it common that the annotation returned by MITOS is not good enough?

I don't know; I have never used MITOS.

> Or could that indicate that maybe there is a problem in the assembled genome? (looks unlikely except for a few SNPs as three assemblers return almost the same sequence)

I don't know; you haven't provided the assembled gene sequence or annotation. The response from NCBI suggests that there is a problem with either the sequence or annotation (e.g. nad5 extends into an annotated tRNA sequence).

> Should we edit the annotation ourselves or there is a software that does it itself?

If there is an error in the annotation through automated methods, you should fix it manually.

> How often using MITOS is enough or how often manually editing is needed?

I don't know; I have never used MITOS.
> 2 votes

# Answer
You might be luckier by submitting to INSDC by the ENA gate instead of NCBI.

To do so you first have to convert your GFF into the EMBL format used by ENA with the help of EMBLmyGFF3. During the conversion you can use the `--tranlate` option that can help to see if there is anything wrong in the ORFs (e.g. frameshifts).

Incomplete genes is supposed to be supported by INSDC databases, so I don't get why NCBI complains about it. Either their pipeline does not take care of this, or it is something specific to Mitochondrion annotations. Using `EMBLmyGFF3` incomplete CDS will be handled properly and the location field in the resulting EMBL flat file will have the information that it is incomplete (see `3.4.3 Location examples` in the documentation http://ftp.ebi.ac.uk/pub/databases/embl/doc//FT\_current.html#3.4)

Then once you have your EMBL flat file you can check by yourself if there is anything not accepted by the DB using the ENA webin-cli (https://github.com/enasequence/webin-cli/releases/tag/v3.7.0) It there are errors they will be listed with details and you will have to fix your file accordingly. When everything is fixed you can submit to ENA. Once accepted by ENA it will be within 24 hours within NCBI.
> 2 votes

---
Tags: assembly, sequence-annotation
---


# Question
Title: What is a good RNA seq normalization method that allows for across sample comparisons and between transcripts
What is a good RNA seq normalization method that allows for across sample comparisons, **and** allows between transcripts comparisons as well? I read that TMM for example allows across sample comparisons but not between transcripts. On the other hand, I read that TPM is not good enough for between sample comparisons

# Answer
> 0 votes
I found a recent technique that allows for both, called geTMM.

# Answer
> 0 votes
I'm having a bit of trouble understanding the question, but we've used DESeq2 for doing both between-sample and between-transcript comparisons.

When comparing transcripts for short reads I use the variant-stabilised transformation, adjusted for transcript length. This is good for a visual overview of results (e.g. in a heatmap).

When comparing samples I use the variance-adjusted DESeq2 test, producing normalised $\textrm{log}\_2FC$ values. This works well for quantitative results.

See this paper, where we use both of these approaches (e.g. Figure 2 shows both in a side-by-side view).

# Answer
> 0 votes
As mentioned by @gringer DESeq2 is a popular choice for normalization of counts between samples. However, because each gene is adjusted individually it is not well suited to compare within a sample.
For this a global normalization like `TMM` from `edgeR` would be more appropriate.

You can find more information in this tutorial from the hbc at Harvard.

---
Tags: rna-seq, normalization
---


title: ensembldb equivalent in python

text:
# Question
Title: ensembldb equivalent in python
Is there a python equivalent to ensembldb? I want to get genome coordinates for a transcript (like the function `transcriptToGenome`) but need to do it in python. I wrote a R script that I'm calling from python but I was hoping there was a more direct way to do this.

Thanks

# Answer
> 4 votes
`pyensembl` should do the job:

```
>>> from pyensembl import EnsemblRelease
>>> data = EnsemblRelease(76)
>>> data.transcript_by_id("ENST00000506751")
Transcript(..., contig='5', start=140861224, end=140863521, strand='+', genome='GRCh38')

```

---
Tags: bioconductor, python
---


# Question
Title: How can I make Snakemake ignore all previously produced/cached output files?
I want to rerun a Snakemake workflow, ignoring all previously produced intermediary outputs.

Snakemake has about 100 different flags. I just want to make it easier to find an answer to this. What's the syntax/CLI flag I need to add?

# Answer
> 1 votes
If you want Snakemake to ignore the previous output of all rules, add the flag `--forceall`. Otherwise, if you just want to ignore the output of a specific rule to rerun that rule, add `--force`:

```
snakemake --forceall -c1

```

The command line interface documentation is helpful as a reference: https://snakemake.readthedocs.io/en/stable/executing/cli.html#all-options

---
Tags: snakemake
---


# Question
Title: Copying files from multiple directories to a single
I have multiple .GFF files, each in a single directory per annotated whole genome, that I wish to copy into a single directory so that I can perform a core genome alignment.

However the following code does not seem to work:

```
for pathname in "~/ST_0075/annotations/*/*.gff"; do cp "$pathname" ~/ST_0075/alignment; done

```

The directory structure is the following:

```
Home
-ST_0075
---Annotations
-----Sample_001
----------001.gff

```

The goal is to copy from multiple Sample\_XXX subdirectories within `Annotations` the XXX.gff files into another directory named `alignment`

# Answer
Tinkered a bit and realized the solution did not neccistate a for loop.

```
mkdir ~/ST_0075/alignment; cp ~/ST_0075/annotations/*/*.gff ~/ST_0075/alignment

```
> 3 votes

---
Tags: linux, data-management
---