# Question
Title: reverse translation from amino acid string to DNA strings
what is the opposite of .translate() function calls ?

I mean let's say I am given an amino acid string `CYCLIC`, how do I obtain **all the possible combinations of DNA strings** before the translation process ?

Note: I was told that degeneracy has an important role in this particular aspect, but I am not entirely sure about the degeneracy calculation in the following screenshot highlighted in yellow color.

The degenerate code is as follows:

```
table = {
        'ATA':'I', 'ATC':'I', 'ATT':'I', 'ATG':'M',
        'ACA':'T', 'ACC':'T', 'ACG':'T', 'ACT':'T',
        'AAC':'N', 'AAT':'N', 'AAA':'K', 'AAG':'K',
        'AGC':'S', 'AGT':'S', 'AGA':'R', 'AGG':'R',
        'CTA':'L', 'CTC':'L', 'CTG':'L', 'CTT':'L',
        'CCA':'P', 'CCC':'P', 'CCG':'P', 'CCT':'P',
        'CAC':'H', 'CAT':'H', 'CAA':'Q', 'CAG':'Q',
        'CGA':'R', 'CGC':'R', 'CGG':'R', 'CGT':'R',
        'GTA':'V', 'GTC':'V', 'GTG':'V', 'GTT':'V',
        'GCA':'A', 'GCC':'A', 'GCG':'A', 'GCT':'A',
        'GAC':'D', 'GAT':'D', 'GAA':'E', 'GAG':'E',
        'GGA':'G', 'GGC':'G', 'GGG':'G', 'GGT':'G',
        'TCA':'S', 'TCC':'S', 'TCG':'S', 'TCT':'S',
        'TTC':'F', 'TTT':'F', 'TTA':'L', 'TTG':'L',
        'TAC':'Y', 'TAT':'Y', 'TAA':'_', 'TAG':'_',
        'TGC':'C', 'TGT':'C', 'TGA':'_', 'TGG':'W',
    }

```

# Answer
Given the table above the solution for CYCLE is:

```
from itertools import product

prot = 'CYCLE'
degenDict = dict()
for v,k in table.items():
    if k in degenDict:
        degenDict[k].append(v)
    else:
        degenDict[k]=[]
        degenDict[k].append(v)
nucs = [degenDict[resid] for resid in prot]
for degenNuc in product(*nucs):
    print(''.join(degenNuc))

```

> TGCTACTGCCTAGAA TGCTACTGCCTAGAG TGCTACTGCCTCGAA TGCTACTGCCTCGAG TGCTACTGCCTGGAA TGCTACTGCCTGGAG TGCTACTGCCTTGAA TGCTACTGCCTTGAG TGCTACTGCTTAGAA TGCTACTGCTTAGAG TGCTACTGCTTGGAA TGCTACTGCTTGGAG TGCTACTGTCTAGAA TGCTACTGTCTAGAG TGCTACTGTCTCGAA TGCTACTGTCTCGAG TGCTACTGTCTGGAA TGCTACTGTCTGGAG TGCTACTGTCTTGAA TGCTACTGTCTTGAG TGCTACTGTTTAGAA TGCTACTGTTTAGAG TGCTACTGTTTGGAA TGCTACTGTTTGGAG TGCTATTGCCTAGAA TGCTATTGCCTAGAG TGCTATTGCCTCGAA TGCTATTGCCTCGAG TGCTATTGCCTGGAA TGCTATTGCCTGGAG TGCTATTGCCTTGAA TGCTATTGCCTTGAG TGCTATTGCTTAGAA TGCTATTGCTTAGAG TGCTATTGCTTGGAA TGCTATTGCTTGGAG TGCTATTGTCTAGAA TGCTATTGTCTAGAG TGCTATTGTCTCGAA TGCTATTGTCTCGAG TGCTATTGTCTGGAA TGCTATTGTCTGGAG TGCTATTGTCTTGAA TGCTATTGTCTTGAG TGCTATTGTTTAGAA TGCTATTGTTTAGAG TGCTATTGTTTGGAA TGCTATTGTTTGGAG TGTTACTGCCTAGAA TGTTACTGCCTAGAG TGTTACTGCCTCGAA TGTTACTGCCTCGAG TGTTACTGCCTGGAA TGTTACTGCCTGGAG TGTTACTGCCTTGAA TGTTACTGCCTTGAG TGTTACTGCTTAGAA TGTTACTGCTTAGAG TGTTACTGCTTGGAA TGTTACTGCTTGGAG TGTTACTGTCTAGAA TGTTACTGTCTAGAG TGTTACTGTCTCGAA TGTTACTGTCTCGAG TGTTACTGTCTGGAA TGTTACTGTCTGGAG TGTTACTGTCTTGAA TGTTACTGTCTTGAG TGTTACTGTTTAGAA TGTTACTGTTTAGAG TGTTACTGTTTGGAA TGTTACTGTTTGGAG TGTTATTGCCTAGAA TGTTATTGCCTAGAG TGTTATTGCCTCGAA TGTTATTGCCTCGAG TGTTATTGCCTGGAA TGTTATTGCCTGGAG TGTTATTGCCTTGAA TGTTATTGCCTTGAG TGTTATTGCTTAGAA TGTTATTGCTTAGAG TGTTATTGCTTGGAA TGTTATTGCTTGGAG TGTTATTGTCTAGAA TGTTATTGTCTAGAG TGTTATTGTCTCGAA TGTTATTGTCTCGAG TGTTATTGTCTGGAA TGTTATTGTCTGGAG TGTTATTGTCTTGAA TGTTATTGTCTTGAG TGTTATTGTTTAGAA TGTTATTGTTTAGAG TGTTATTGTTTGGAA TGTTATTGTTTGGAG
> 3 votes

---
Tags: python, phylogenetics, phylogeny, sequence-analysis
---


# Question
Title: mass spectrum matching with cyclic peptide and linear peptide
I was doing a simple quiz which I encountered the following two questions related to mass spectrum matching with cyclic peptide and linear peptide.

I have done some preliminary google search, but I am not sure how to use either https://codereview.stackexchange.com/questions/193859/find-a-cyclic-peptide-with-a-given-non-ideal-spectrum or https://github.com/rochoa85/pepfun in this case.

> Which of the following cyclic peptides could have generated the theoretical spectrum 0 71 101 113 131 184 202 214 232 285 303 315 345 416? (Select all that apply.)

```
TAIM

TMIA

MAIT

MTAI

TMLA

TLAM

```

> Which of the following linear peptides is consistent with Spectrum = {0 71 99 101 103 128 129 199 200 204 227 230 231 298 303 328 330 332 333}? (Select all that apply.)

```
AVQ

TCQ

TVQ

QCV

TCE

ETC

```

# Answer
You are way overthinking this.

As this is homework I won't give any code, but I will throw back some questions:

* What is the definition of m/z?
* Were the only tool available Wikipedia, could you simply read the infobox of each residue?
* If so, what values would you need and what elements do you need to subtract?
* Does residue order matter for linear and backbone cyclised peptides? (Let's ignore isoaspartate etc.)
* Have you been introduced to RDKit in your course? E.g. what would this give you?

```
from rdkit import Chem
from rdkit.Chem import Descriptors
mol: Chem.Mol = Chem.MolFromFASTA('ELVIS')
print( Descriptors.ExactMolWt(mol) )

```
> 3 votes

---
Tags: proteins
---


title: Parsing text within a taxonomy.csv file

text:
# Question
Title: Parsing text within a taxonomy.csv file
I have the following file in tab separated values (TSV).

| ID | Taxonomy |
| --- | --- |
| 76980f6d906c6baaef3f96559fc7ba1b | d\__Eukaryota; p\__Cercozoa; c\__Vampyrellidae; o\__Vampyrellidae; f\__Vampyrellidae; g\__uncultured |
| 4f808d3adec894395304661391037369 | d\__Eukaryota; p\__Ochrophyta; c\__Xanthophyceae; o\__uncultured; f\__uncultured; g\__uncultured |
| 777584465097568f5b30bc9e2b6584c1 | d\__Eukaryota; p\__Labyrinthulomycetes; c\__uncultured; o\__uncultured; f\__uncultured; g\__uncultured |
| 5a99dafabe78847d8f5fed0d37f8672a | d\__Eukaryota; p\__Labyrinthulomycetes; c\__Labyrinthulomycetes; o\__Labyrinthulomycetes; f\__Thraustochytriaceae; g\__uncultured |

I would like to modify the "g\__uncultured" with the last known rank possibly using either shell or R scripts.

The new table would then read:

| ID | Taxonomy |
| --- | --- |
| 76980f6d906c6baaef3f96559fc7ba1b | d\__Eukaryota; p\__Cercozoa; c\__Vampyrellidae; o\__Vampyrellidae; f\__Vampyrellidae; g\__f\__Vampyrellidae |
| 4f808d3adec894395304661391037369 | d\__Eukaryota; p\__Ochrophyta; c\__Xanthophyceae; o\__uncultured; f\__uncultured; g\__c\__Xanthophyceae |
| 777584465097568f5b30bc9e2b6584c1 | d\__Eukaryota; p\__Labyrinthulomycetes; c\__uncultured; o\__uncultured; f\__uncultured; g\__p\__Labyrinthulomycetes |
| 5a99dafabe78847d8f5fed0d37f8672a | d\__Eukaryota; p\__Labyrinthulomycetes; c\__Labyrinthulomycetes; o\__Labyrinthulomycetes; f\__Thraustochytriaceae; g\__f\__Thraustochytriaceae |

# Answer
> 1 votes
The `pandas` solution is here:

```
#!/usr/bin/env python3
import pandas as pd
from pathlib import Path

path = '/Whatever'
infile = 'myfile.csv'
outfile = 'outfile.csv'
pd.set_option('display.max_colwidth', None)
df = pd.read_csv(Path(path,infile), sep='\t')
taxaChange = {'76980f6d906c6baaef3f96559fc7ba1b':'g__f__Vampyrellidae', '4f808d3adec894395304661391037369':'g__c__Xanthophyceae'}
df2 = df.set_index('ID')
df3 = df2
for x,y in taxaChange.items():
    df3.loc[x, 'Taxonomy'] = df2.loc[x, 'Taxonomy'].replace('g__uncultured', y)
print (df3)

df.to_csv(Path(path,outfile), sep='\t')

```

Output

```
ID                                       Taxonomy
76980f6d906c6baaef3f96559fc7ba1b         d__Eukaryota; p__Cercozoa; c__Vampyrellidae; o__Vampyrellidae; f__Vampyrellidae; g__f__Vampyrellidae
4f808d3adec894395304661391037369         d__Eukaryota; p__Ochrophyta; c__Xanthophyceae; o__uncultured; f__uncultured; g__c__Xanthophyceae
777584465097568f5b30bc9e2b6584c1         d__Eukaryota; p__Labyrinthulomycetes; c__uncultured; o__uncultured; f__uncultured; g__uncultured
5a99dafabe78847d8f5fed0d37f8672a         d__Eukaryota; p__Labyrinthulomycetes; c__Labyrinthulomycetes; o__Labyrinthulomycetes; f__Thraustochytriaceae; g__uncultured

```

Both `g__f__Vampyrellidae` and `g__c__Xanthophyceae` are there, i.e. the code was tested and worked fine.

The point for the `pandas` approach is you've a natural table format so it makes sense to take it straight to a dataframe (table -\> dataframe). regex-ing using `\t` and `\n` is of course the alternative and you need to formulate the table as a hash (Perl) or dictionary (Python).

**Notes** The `df` and `df2` can reduced a single dataframe by updating the `df` (`inplace` was omitted and will work for an index), rather than creating a new `df`. However, the `df3` is needed (need to copy the dataframe) because 'inplace' didn't work. `inplace` is a command to force the dataframe to update. The `set_option` is just so I can print the display without truncation.

*Finally* You substitute the infile and outfile to correct path and name values, save it as a `.py` file and just `python3 myscript.py`.

---
Tags: r, phylogenetics, taxonomy, parsing
---


# Question
Title: Improving list speed: glycosylation example
I discussed a question with @gaspanic Python/Biopython - Replace amino acid residue on MSA with "z" from a list of unaligned positions . The issue emerged was speeding up lists in Python. The solutions we discussed were

1. Raw code
2. Predefining a list size
3. `deque`

| 1. @gaspanic | my code | 2. @gaspanic | my code | 3. @gaspanic | my code |
| --- | --- | --- | --- | --- | --- |
| 3.45-3.48 | 7.17-7.2 | 4.38-4.85 | 7.09-7.23 | 3.44-3.53 | 7.21-7.39 |

Run timings are seconds (3 replicated) above with @gaspanic's code being faster.

What improvements would speed this up?

Here I extended the data set

```
from itertools import repeat
def iterit(x):
    seq = 'M----KFLAFLCLLGFANAQ-------------D-GKCG----TLSNKSPS--------------------K'
    return ''.join(list(repeat(seq, x)))

if __name__ == '__main__':
    seqE3 = iterit(2800) # length is 201600
    seqE4 = iterit(28000) # length is 2016000
    seqE5 = iterit(280000) # length is 20160000

```

@gaspanic code

```
for seqbit in [seqE3,seqE4,seqE5]:
    seq = seqbit
    pos = [2, 5, 10, 11, 14]
    new_seq = []
    count = 0
    for s in seq:
        if s == "-":
            new_seq.append(s)
        else:
            count += 1
            if count in pos:
                new_seq.append('Z')
            else:
                new_seq.append(s)
   # print(''.join(new_seq))

```

My code:

```
# same variables and loop structure
newseq = list()
count = 0
for n,s in enumerate(prot):
    if str(s) == '-':
        count += 1
        s = '-'
    else:
        if (n + 1 - count) in pos:
            s = 'Z'
    newseq.append(s)
# print(''.join(newseq))

```

Predefining lists code (additional counter needed for @gaspanic's code:

```
# same variables and loop structure
new_seq = [] * len(seq)
count = 0
for n,s in enumerate(prot):
    if str(s) == '-':
        count += 1
        s = '-'
    else:
        if (n + 1 - count) in pos:
            s = 'Z'
    newseq[n] = s

```

Deque code simply any of the above substituting

```
newseq = list()

```

with

```
from collections import deque
newseq = deque()

```

Suggestions for efficiency improvements welcome.

A solution in `Cython` code would be cool. `PyPI` is not sought here.

# Answer
I'd start by saying that Python is inherently slow and this feels like a case of perhaps unnecessary optimisation. With that in mind here are my suggested (untested) optimisations of @gaspanic code:

```
pos = {2, 5, 10, 11, 14} # make pos a set for O(1) member testing, and don't define inside the loop
for seq in (seqE3,seqE4,seqE5): # Define seq here rather than inside the loop, and use an immutable tuple instead of a mutable list
    new_seq = []
    count = 0
    for s in seq:
        if s == "-":
            new_seq.append(s)
        else:
            count += 1
            new_seq.append('Z' if count in pos else s) # use ternary operator

```
> 2 votes

# Answer
Huge thanks to @Chris\_Rands, this solution increases runtime by ~ 0.5 seconds uniformly across all coded solutions, so this code has a runtime of just over 3 seconds runtime on the benchmark machine. I incorporated my solution with @Chris\_Rands solution and the final runtime I got was:

@gaspanic

> 1.7 - 1.8 seconds

This being twice the fastest speed in the question and four times the speed of my raw code.

My code

> 4.4 - 4.5 seconds

However a unique combination is needed:

1. `tuple` (Chris Rands)
2. `set` \- important (Chris Rands)
3. an `iterator`
4. a `list comprehension`

In particular 3. and 4. are needed in combination. This is explained below.

```
def iterit(x):
    seq = 'M----KFLAFLCLLGFANAQ-------------D-GKCG----TLSNKSPS--------------------K'
    return tuple(''.join(list(repeat(seq, x))))

def mIterit(seqE3,seqE4,seqE5):
    pos = {2, 5, 10, 11, 14}
    for seqbit in (seqE3,seqE4,seqE5):
        seq = seqbit
        count = 0
        # loc = 0
        for s in seq:
            if s == "-":
                yield '-'
            else:
                count += 1
                if count in pos:
                    yield 'Z'
                else:
                    yield s

f __name__ == '__main__':

    seqE3 = iterit(2800) # length is 201600
    seqE4 = iterit(28000) # length is 2016000
    seqE5 = iterit(280000) # length is 20160000

    listComp = [i for i in mIterit(seqE3,seqE4,seqE5)]

```

The list comprehension (variable `listComp`) is the very last line and is crucial to the solution otherwise there is little increase in speed to @Chris\_Rand's code. This alternative is to replace the last line with:

```
newseq = list()
for i in mIterit(seqE3,seqE4,seqE5):
    newseq.append(i)

```

A simple:

```
for i in mIterit(seqE3,seqE4,seqE5):
    pass

```

Is the same speed as the list comprehension, but its of no value because the output isn't saved to anything.

I was impressed by the interaction between the `iterator` and `list comprehension`. The other notable influence on speed is the `set` described by @Chris\_Rand.

How does it work? `tuples` do increase speed by removing some functional of a `list` (immutability). Here the gain is small. Placing the `pos` variable outside the loop speaks for itself (0.1 seconds gain). `set` is interesting.

The `iterator` spews out a stream of results via `yield` which is where it attains its speed. The problem is capturing that stream outside the `iterator` is expensive. However `list comprehensions` are faster than standard loops and in this scenerio carry little/no detectable overhead. I was impressed. Up to this point I thought `list comprehensions` were a bit faster and a lot 'cooler' than loops. I hadn't understood there was a synergy with an `iterator`.

Is it all worth it? `Python` can easily achieve faster run times than `Java` using Cython but I will place the discussion in the *chat* forum next week.
> 1 votes

---
Tags: python, phylogenetics, phylogeny, algorithms, benchmarking
---


title: Count reads at specific gene features

text:
# Question
Title: Count reads at specific gene features
I have BAM files from an RNA-Seq experiment and for all genes (or a subset) I want to get the number of reads in regions around the TSS (e.g. 2kbp) and the TES (e.g. 2 kbp) and calculate the ratio between them. I was thinking of first creating these regions for each gene and store them in a BED file and then use `bedtools intersect` with the BAM files to get the reads in these regions. My questions is how would I make such a BED file that contains only the specified regions downstream of the TSS and specified regions upstream of the TSS?

# Answer
> 3 votes
First of all, you need gene coordinates. You can retrieve them in GTF format from GENCODE or other sources.

Next, extract TSS and TES +/- 2kb windows from each gene, taking in mind that genes can be in the forward or reverse strand. Something like this:

```
awk -F"\t|\"" -vOFS="\t" '$3=="gene" {
  TSS=($7=="+" ? $4-2000"\t"$4+2000 : $5-2000"\t"$5+2000)
  TES=($7=="+" ? $5-2000"\t"$5+2000 : $4-2000"\t"$4+2000)
  print $1,TSS,$10"_TSS_2kb",0,$7
  print $1,TES,$10"_TES_2kb",0,$7
}' annotation.gtf > tss_tes.bed

```

Here I used GNU Awk from Bash, but you can write your own solution using Python, R or whatever you like.

Finally, count reads falling in these intervals with your tool of choice (you mentioned bedtools, which is fine).

---
Tags: sam, ngs, bedtools, bed
---


# Question
Title: Tool for cross referencing a newly found MOTIF to database of known MOTIFS
I need a tool or a function I can use in my code (R, or Python) that I can cross reference a MOTIF against known MOTIFS, a function that will take as input a MOTIF (a probability weight matrix, PWM which is a 4 by m matrix for an m-length motif) and spit out the closest -n- known motif(s). Is there such a function/tool or something close to it? Thanks much.

# Answer
> 2 votes
A popular tool for doing this is tomtom, part of the meme package.

You can try it via interactive web form here.

Once installed locally it can be run via command line

Additionally, it is available in the Bioconductor package `memes` as explained in Motif Comparison using TomTom

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

Hi @ZebraFish Motifs package Biopython

https://biopython.org/docs/1.76/api/Bio.motifs.html

---
Tags: r, rna-seq, python, motifs
---


# Question
Title: Permutations of 4 amino acid sequence oligopeptide
I have an assignment for a project to build a tetrapeptides library to screen against some target enzymes *in silicon*.

I am looking for a script to construct about 160000 tetrapeptide from the 20Amino acids permutations.

I am quite familiar with Bash, but am struggling to get a function to working.

# Answer
The Python solution is here. IMO Bash definitely isn't a good solution for this calculation.

```
from itertools import permutations, combinations

one2all ={'A': ('A', 'ALA', 'alanine'),
          'R': ('R', 'ARG', 'arginine'),
          'N': ('N', 'ASN', 'asparagine'),
          'D': ('D', 'ASP', 'aspartic acid'),
          'C': ('C', 'CYS', 'cysteine'),
          'Q': ('Q', 'GLN', 'glutamine'),
          'E': ('E', 'GLU', 'glutamic acid'),
          'G': ('G', 'GLY', 'glycine'),
          'H': ('H', 'HIS', 'histidine'),
          'I': ('I', 'ILE', 'isoleucine'),
          'L': ('L', 'LEU', 'leucine'),
          'K': ('K', 'LYS', 'lysine'),
          'M': ('M', 'MET', 'methionine'),
          'F': ('F', 'PHE', 'phenylalanine'),
          'P': ('P', 'PRO', 'proline'),
          'S': ('S', 'SER', 'serine'),
          'T': ('T', 'THR', 'threonine'),
          'W': ('W', 'TRP', 'tryptophan'),
          'Y': ('Y', 'TYR', 'tyrosine'),
          'V': ('V', 'VAL', 'valine'),
          'X': ('X', 'GLX', 'glutaminx'),
          'Z': ('Z', 'GLI', 'glycine'),
          'J': ('J', 'NLE', 'norleucine'),
          'U': ('U', 'CYC', 'cysteinc')}
myprotein = [x for x in one2all.keys()]
tetrapeps = combinations(myprotein,4)
peptide = []
count = 0
for x in list(tetrapeps):
    peptide.append(''.join(x))
for tetra in peptide:
    permpeptides = list(permutations (tetra))
    for i in permpeptides:
        count +=1
        print (''.join(i))

```

I get 255024 permutations so more than the estimate of 160000, however the code is correct (please see **Validation** below).

Truncated output

> .... JUXV UVXJ UVJX UXVJ UXJV UJVX UJXV VZJU VZUJ VJZU VJUZ VUZJ VUJZ ZVJU ZVUJ ZJVU ZJUV ZUVJ ZUJV JVZU JVUZ JZVU JZUV JUVZ JUZV UVZJ UVJZ UZVJ UZJV UJVZ UJZV XZJU XZUJ XJZU XJUZ XUZJ XUJZ ZXJU ZXUJ ZJXU ZJUX ZUXJ ZUJX JXZU JXUZ JZXU JZUX JUXZ JUZX UXZJ UXJZ UZXJ UZJX UJXZ UJZX

I cut and paste the amino dictionary from here (I wasn't going to write it by hand), `AminoAcid.py` by Pasteur.

---

**Validation**

The theoretical number of permutations from 4 samples of 20 amino acids is

> **116280**

However, the amino acid chart I used comprised **24** amino acids (I just copied and pasted it). The theoretical number for 24 amino acids is therefore:

> **255024**

This was the number observed in the output of the above code and validated the calculation.

A permutations calutation can be performed online here

I am not sure how the expected output of 160000 was derived but suspect it was a guess.

---

The permutation equation is as follows, the output of this equation and the code above coincide,

$$ P\_{n, r}=\frac{n !}{(n-r) !}$$

n = no. amino acids

r = size of oligopeptide

---

#### Permutation by repetition

If the permutation by repetition equation is sought (comment below). The code is:

```
from itertools import product
myprotein = [x for x in one2all.keys()]
print(len([p for p in product(myprotein, repeat=4)]))

```

for 20 amino acids this is:

> 160000

for 24 amino acids:

> 331776
> 2 votes

# Answer
Here's a bash solution using brace expansion and the single-letter codes for the 20 common amino acids:

```
aalist=D,E,K,R,H,N,Q,S,T,P,G,A,V,I,L,M,C,F,Y,W
eval printf '"%s\n"' {"$aalist"}{"$aalist"}{"$aalist"}{"$aalist"}

```

This prints 160,000 lines, one for each tetrapeptide.
> 1 votes

---
Tags: python, phylogenetics, phylogeny, protein-structure, bash
---


# Question
Title: Refactoring a script mapping ID to sequence
#

I wrote simple Python script which has an excessive run time, notably when using large data sets. For example a data set of 1 000 000 sequences.

I am seeking assistance refactoring the code and would be very grateful for any assistance.

```
from Bio import SeqIO
import re

def move_item(old, new):
    k, v = old.popitem()
    new[k] = v

def map_ids_to_aa_sequence_list(id_):
    id_ = id_.strip()
    for s in aa:
        if re.search(id_,s, re.I):
            return s

file = open(r'list_of_ids', 'r').readlines()
names = [for name.strip() in file if len(name) > 1]

handle = open('sequences.fasta','r')
sequences = list(SeqIO.parse(handle, 'fasta'))

seq_dict = {}

for s in sequences:
    seq_dict[str(s.description)] = str(s.seq)

aa = []

for k, v in seq_dict.items():
    sequence = '>' + k +'\n' + v + '\n'
    aa.append(sequence)

separated_aminoacid_sequences_set = set(map(map_ids_to_aa_sequence_list, names))

separated_aminoacid_sequences_list = list(separated_aminoacid_sequences_set)

with open('separated_sequences.fasta','w') as fs:
    for s in separated_aminoacid_sequences_list:
        fs.write(str(s))

```

# Answer
This is beyond a single answer. 'Refactoring' (code efficiency, high performance code) is a huge topic. Like its own project.

The first thing I would personally absolutely target first is `re` here:

```
    if re.search(id_,s, re.I):

```

Python is not Perl and the `re` is EXTREMELY slow and in particular using `search` is especially slow, most notably when uncompiled (see notes and point 3 below).

### How to avoid `re.search`

Could you revise this to a direct `C` style operation? This is the best solution if you get it working. For example,

1. Using `==`:

e.g.

```
if str(id_) == str(s):

```

if the case is important (like they're different) then...

```
if str(id_).lower() == str(s).lower():

```

You could dump the `str()` it ain't doing much its just bug insurance.

The other approach if you want to stick with `re` is to switch to `match`, that will be quicker.

2. Using `match`

e.g.

```
if re.match(id_, s, re.I):

```

`match` anchors the search at the first character, i.e. id\_\[0\], `search` goes everywhere. `match` is faster.

3. Sticking with `search` but combining with compile:

e.g.

```
id_ = id_.strip()
x = re.compile (id_)
if x.search(s, re.I): # don't use re. x replaces it

```

or

```
id_ = id_.strip()
x = re.compile (id_)
if x.match(s, re.I):

```

I have not put these in order of preference. My preference is for 1. but it will risk throwing logical errors (missing stuff).

*If you are using `re` you must compile the search term*. Its a pain but whilst compiled `re` is slow, uncompiled `re` is really, really slow. Disclaimer, I'm from a Perl background and from this perspective `re` is the pits because previously relied on the efficient regex (Perl).

The `==` and `match` could fall over if `>` is in one variable and not in the other, but the `search` will return the correct 'match' (i.e. not the output of 'match') at the cost of algorithm efficiency.

This isn't an exhaustive analysis of `re` but thats where I would start.

Please note this is a starting point from prior experience to pick a 'big issues'. It is useful for future, particularly for questions, to place time checks across the code to identify the bottlenecks.

---

**Notes**

* The thing we don't know is what the fasta headers look like.
* For every query in `search` there are $1.10^6$ fasta headers. So if you've $1.10^5$ queries thats $1.10^{11}$ searches for `search` thats why it can be taking up a lot of time.
* Rewriting the code using Biopython objects throughout is good practice. `SeqIO` in particular are cool. `AlignIO` are algorithmically expensive however.
> 2 votes

---
Tags: python, algorithms, benchmarking
---


# Question
Title: Find a genomic coordinates for a protein aminoacid position
Is there a function that can map genomic position (hg19) back from a protein position? I can have name of a particular transcript, and exon number.

For example, I have KRAS gene for which I would like to know the exact genomic position of G12, and I expect to have an output saying:

chr12, from 25398283 to 25398285.

transvar could be an excellent solution if the current built was working (broken links for fasta annotations). The web search is working but not allowing files heavier than 5kB which is my case.

# Answer
You can do this manually by looking up the protein variant on https://varsome.com. For your example, that would mean searching for `KRAS G12V` (change `V` for anything else, the idea is to give an HGVS protein notation for a variant, so any single letter amino acid code would work). That would take you to this page https://varsome.com/variant/hg19/KRAS%20G12V?annotation-mode=germline, which includes:

Of course, I imagine you want to automate this. VarSome can help there too, through our API, but that is a paid service. If you have the option of paying for this, the VarSome API can be used to (among other things) look up arbitrary variants and translate between different notations, so you could easily adapt this to get the information you need.

***Disclosure***: I work for the company behind VarSome.
> 1 votes

---
Tags: sequence-annotation, codon, genomicranges
---


# Question
Title: How do I read mmCIF files from Alphafold DB as coordinates to use in replicating their visualization for Unity VR?
https://alphafold.com/entry/Q9Y499

Starting at line 918 of the downloadable mmCIF file, it looks like there might be coordinates in x, y,z:

```
ATOM 1    N N   . TRP A 1 1   ? 12.549  -2.565  -21.811 1.0 96.06 ? 1   TRP A N   1 Q9Y499 UNP 1   W
ATOM 2    C CA  . TRP A 1 1   ? 12.477  -2.290  -20.365 1.0 96.06 ? 1   TRP A CA  1 Q9Y499 UNP 1   W
ATOM 3    C C   . TRP A 1 1   ? 11.221  -2.928  -19.812 1.0 96.06 ? 1   TRP A C   1 Q9Y499 UNP 1   W
ATOM 4    C CB  . TRP A 1 1   ? 12.469  -0.785  -20.103 1.0 96.06 ? 1   TRP A CB  1 Q9Y499 UNP 1   W
ATOM 5    O O   . TRP A 1 1   ? 10.206  -2.922  -20.501 1.0 96.06 ? 1   TRP A O   1 Q9Y499 UNP 1   W
ATOM 6    C CG  . TRP A 1 1   ? 13.795  -0.114  -20.270 1.0 96.06 ? 1   TRP A CG  1 Q9Y499 UNP 1   W
ATOM 7    C CD1 . TRP A 1 1   ? 14.286  0.508   -21.369 1.0 96.06 ? 1   TRP A CD1 1 Q9Y499 UNP 1   W
ATOM 8    C CD2 . TRP A 1 1   ? 14.791  0.098   -19.228 1.0 96.06 ? 1   TRP A CD2 1 Q9Y499 UNP 1   W

```

I am a bioinformatics noob (took an intro class in my undergrad CS degree at a Canadian university) but know some Unity game engine with which I want to build a standalone visualization tool for users with VR headsets.

# Answer
## Don't

A protein is not only list of atomic coordinates. Bar for a proof of concept (i.e. learn to code in Unity VR), it is an endless spiral of hurt in a field where there is a lot of competition and were you to push a product you'd need to invest more manhours in publicity than in coding...

## Connectivity

Here is the short of it. Chemistry is classically a graph network: nodes=atoms, edges=bonds. The connectivity of the atoms is dictated by their residue's 3-letter name, which is either

* assumed —software know the 20 AAs (and several modifications, eg. `SEP`, 4x2 bases and some ligands, e.g. `HOH`— or
* a special entry in the PDB or mmCIF file —in PDB it's a CONECT entry. E.g. https://www.rcsb.org/ligand/CFF

AlphaFold2 does not have ligands, AlphaFill does. But scientists (not PR specialists) use PDB entries or run their own modelling. AlphaFold2 from EBI is rarely used: oligomers with ligands is where it's at. Especially if doing a drug discovery campaign (and red biotech can and does afford VR headsets).

## Format

mmCIF is a deposition standard, PDB is the workhorse standard. mmTF is a cool idea that is not catching on (cf XKCD comic strip about standard proliferation...). Using PDB is easier: https://www.wwpdb.org/documentation/file-format-content/format33/v3.3.html

## Representation

Problems arise when offering representations —everyone has their favourite. See NGL gallery for an array: https://nglviewer.org/ngl/gallery/index.html Surface display has additional issues —PyMOL simply and openly reuses the Advance Poisson-Boltzmann Solvation tool as it's tricky.

## Dragons in the minefield

In a PDB entry not from AF2, you can have disulfides (`SSBOND` and/or `CIZ` entry), isopeptide bonds (`LINK`), missing atoms, the `UNK`, `UNX` and `UNL` residues, alphatraces, insertion sequences, alt. occupancy, gaps, non-standard usages of all of these etc etc. and my favourite implied proximity bonding.
> 3 votes

---
Tags: proteins
---


# Question
Title: Problems in conducting batch effect correction
I am now trying to conduct batch effect correction for my data matrix, whose samples (human) have different gender and come from different sequencing platforms. I have tried `Combat`, `Combat-seq`, and `limma`. However, the outcomes are bad. It seems that `Combat-seq` is based on Negative binomial regression models, `Combat` is based on Gaussian distribution, and `limma` is based on the linear model. I wonder if my data violates their assumption, so I check my data and get this outcome. It seems that my data follows Gamma distribution. I am wondering what I could do in this situation.

---

The data I am working on is DNAse-seq data, which is a matrix with 16 samples. As the samples come from different platforms, I would like to correct the batch effect. I am also considering gender information as there are male and female samples. To sum up, my questions are:

1. Do we need to consider the distribution of the data when we apply the batch correction (limma, combat, and combat-seq)?
2. If we need to consider the distribution, which method do I need to correct the data which follows gamma distribution?
3. Are there methods designed for read counts? (I know there are many methods designed for surat objects, but do we have others that are more general?)

Thank you very much for all kinds of help.

# Answer
The ComBat-seq method requires raw counts which can be modelled as a negative binomial. The other methods are typically used with normalized counts transformed to log2 scale. Obviously for any method that tries to regress a factor out of a dataset you have to make sure that none of the experimental variables are confounded with the batch variables that you want to regress. More details can be given if you edit the post to say what the data are, how the experimental setup is and what the factors to regress are.
> 1 votes

---
Tags: limma, batch-effects, combat
---


# Question
Title: Why is library(Seurat) causing my R session to have a fatal error?
I recently tried to load library(Seurat) after installing Seurat package. It caused my R studio to have a fatal error and abort my session. I have successfully downloaded it before. I have uninstalled R and R studio and reinstalled, restarted my computer. I am new to R too, so I'm not sure what I'm doing.

I am following the directions provided on the satijalab.org website on how to install Seurat.

Here's a screenshot of the error message.

I ran install.packages("Seurat") in R. R is also crashing as soon as I run library(Seurat). Here's a screenshot right before I hit run, and when it crashes, R application just closes.

# Answer
> 1 votes
If you haven't already done so I would make a complete installation of Anaconda, particularly good with bash. The I would either install using Anaconda Navigator, or ...

```
# Create an R environment and activate it
conda create -n r_env r-essentials r-base python=3.7
conda activate r_env

# install Reurat
conda install -c bioconda r-seurat

# installing standard packages
conda install -c r r-randomforest --yes

```

If you have a Mac the instructions are here. In some systems you need to just `activate r_env` .

Please note if you've installed e.g. `python3` with **`pip`** a conflict occurs with the Anaconda package.

# Answer
> 0 votes
I just saw this question. The solution to your problem is here

---
Tags: r, seurat, rstudio
---


# Question
Title: I am trying to run nab to create the pdb file but it isn`t working (AMBER)
Here is the input file for nab:

```
molecule m;
m = arna( "gcuucuucuucuucgc" );
putpdb( "lr16.pdb", m, "-wwpdb -nocid -tr");

```

I am trying to run nab to create the pdb file with this command :

```
$AMBERHOME/bin/nab lr16.nab\./a.out

```

But it comes out saying :

```
gcc: error: lr16.nab./a.out: No such file or directory
cc failed!

```

I have loaded the appropriate modules , please could someone help

# Answer
I have no idea about `nab` ... however

```
“$AMBERHOME/bin/nab lr14.nab./a.out”

```

Looking at your error `lr14.nab./a.out:`

It suggests you didn't put a space character between `lr16.nab` and `a.out`. If you run

```
$AMBERHOME/bin/nab lr14.nab ./a.out

```

I suspect it will work because at present it is looking for a file named lr14.nab./a.out which could be interpreted as a directory of `lr14.nab.` and a input file of `a.out` or the error suggests it is reading `/` as part of the file name. `a.out` of course doesn't exist its your specified output file name

---

output of `ls -l`

```
-rw-r--r-- 1 phys 89 Oct 6 12:58 lr14.nab – Peter 5 hours ago

```
> 1 votes

# Answer
Try this command:

```
$AMBERHOME/bin/nab -v lr14.nab -o ./a.out,

```

then execute a.out separately using

```
./a.out.

```

You will get the `pdb` file of nab extension.
> 1 votes

---
Tags: rna-seq, pdb
---


# Question
Title: statistical test and visualization of genes present in an experiment
I have a list of genes (n = 120) that are involved in breast cancer. And I have a list of differentially expressed genes (n = 70) that are present in the breast cancer gene list. My question is: Are differentially expressed genes significantly more likely to be breast cancer related vs non-differentially expressed genes? which statistical test will be suitable to address this question and how I can visualize its results in `R`.

# Answer
> 1 votes
This question can be addressed by gene set enrichment analysis. You should use Gene Ontology and KEGG Pathways to find out which genes from your gene set are enriched in GO terms and KEGG pathway terms. You can use clusterProfiler package in R for this anlysis and for visualization also. Here is the link to the guide: https://yulab-smu.top/biomedical-knowledge-mining-book/enrichment-overview.html

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


# Question
Title: How to use the contrast function in ebayes for making different comparisons
The design matrix i have is this I would like to know how to use they way its done in deseq2 where we can use the contrast function to make particular comparison

The code im running to make differential testing is this

```
tcgaGlm = lmFit(as.matrix(tcgaExpr), design = tcgaDesign)
tcgaGlm = eBayes(tcgaGlm)

```

Now I will show the part of the data frame and the question Im trying to ask

```
head(tcgaDesign)
             Offset CEBPA DNMT3A FAM5C FLT3_ITD FLT3_TKD IDH1 IDH2 KIT KRAS NPM1 NRAS PHF6 PTPN11 RAD21 RUNX1 SMC1A SMC3 STAG2 TET2 TP53 U2AF1 WT1
TCGA-AB-2803      1     0      0     0        0        0    0    0   0    0    0    0    0      0     0     0     0    0     0    0    0     0   0
TCGA-AB-2805      1     0      0     0        0        0    0    1   0    0    0    0    0      0     0     1     0    0     0    0    0     0   0
TCGA-AB-2806      1     0      0     0        0        0    0    0   0    0    0    0    0      0     0     0     0    0     0    0    0     0   0
TCGA-AB-2807      1     0      0     0        0        0    0    1   0    0    0    0    0      0     0     1     0    0     0    0    0     0   0
TCGA-AB-2808      1     1      0     0        0        0    0    0   0    0    0    1    0      0     0     0     0    0     0    0    0     0   0
TCGA-AB-2810      1     0      0     0        0        0    0    1   0    0    1    0    0      0     0     0     0    0     0    0    0     0   0
               complex  minus5_5q   minus7q  plus8_8q     plus21    t_15_17     t_8_21 inv16_t16_16
TCGA-AB-2803 0.0000000 0.00000000 0.0000000 0.0000000 0.00000000 1.00000000 0.00000000   0.00000000
TCGA-AB-2805 0.0000000 0.00000000 0.0000000 0.0000000 0.00000000 0.00000000 0.00000000   0.00000000
TCGA-AB-2806 0.0000000 0.00000000 0.0000000 0.0000000 0.00000000 0.00000000 1.00000000   0.00000000
TCGA-AB-2807 0.1410256 0.08974359 0.1217949 0.1153846 1.00000000 0.09615385 0.04487179   0.04487179
TCGA-AB-2808 0.0000000 0.00000000 0.0000000 0.0000000 0.00000000 0.00000000 0.00000000   0.00000000
TCGA-AB-2810 0.1410256 0.08974359 0.1217949 0.1153846 0.05063291 0.09615385 0.04487179   0.04487179

```

Here if I have to make a comparison between CEBPA vs DNMT3A mutation in samples which carries CEBPA mutation and other which carries DNMT3A what or how should I do it.

To begin with I followed this manual I m still not sure once you run the model do i need to form a design matrix which comprises of sample which contain mutation in CEBPA and DNMT3A ? or I can go ahead with this experimental design

```
dput(tcgaDesign)
structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0,
1, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0,
1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1,
1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0,
1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0,
0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0,
0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0,
1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1,
0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0,
1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0,
0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1,
0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0,
1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0,
0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1,
0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1,
1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0,
1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0,
0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0.141025641025641,
0, 0.141025641025641, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0.141025641025641, 0, 0, 0, 0.141025641025641, 0, 0,
0, 1, 0, 0, 0, 0.141025641025641, 0.141025641025641, 0.141025641025641,
0, 0, 0, 0, 0.141025641025641, 0, 0, 0.141025641025641, 1, 0,
1, 0, 0, 1, 0, 0, 0, 0.141025641025641, 0, 0, 1, 0, 0, 0, 1,
0, 0, 0, 0.141025641025641, 0, 0, 0, 0, 0, 1, 0, 0, 0.141025641025641,
0, 0, 1, 0.141025641025641, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0,
1, 0, 0, 0, 1, 0.141025641025641, 1, 0, 0, 1, 0, 0, 0, 0, 0,
1, 0, 0.141025641025641, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1,
0.141025641025641, 0.141025641025641, 0, 1, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0.0897435897435897, 0, 0.0897435897435897, 0, 0,
1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0897435897435897,
0, 0, 0, 0.0897435897435897, 0, 0, 0, 1, 0, 0, 0, 0.0897435897435897,
0.0897435897435897, 0.0897435897435897, 0, 0, 0, 0, 0.0897435897435897,
0, 0, 0.0897435897435897, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0.0897435897435897,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0.0897435897435897, 0, 0, 0, 0,
0, 1, 0, 0, 0.0897435897435897, 0, 0, 0, 0.0897435897435897,
0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0.0897435897435897,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0897435897435897, 0, 0, 0,
1, 0, 0, 0, 1, 0, 1, 0, 1, 0.0897435897435897, 0.0897435897435897,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.121794871794872, 0, 0.121794871794872,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.121794871794872,
0, 0, 0, 0.121794871794872, 0, 0, 0, 1, 0, 0, 0, 0.121794871794872,
0.121794871794872, 0.121794871794872, 0, 0, 1, 0, 0.121794871794872,
0, 0, 0.121794871794872, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0.121794871794872,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0.121794871794872, 0, 0, 0, 1,
0, 1, 0, 1, 0.121794871794872, 0, 0, 0, 0.121794871794872, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.121794871794872,
1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0.121794871794872, 0, 0, 0,
1, 0, 0, 0, 0, 0, 1, 0, 0, 0.121794871794872, 0.121794871794872,
0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0.115384615384615, 0, 0.115384615384615,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0.115384615384615,
0, 0, 0, 0.115384615384615, 0, 0, 0, 1, 0, 0, 0, 0.115384615384615,
0.115384615384615, 0.115384615384615, 0, 0, 0, 0, 0.115384615384615,
0, 0, 0.115384615384615, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0.115384615384615,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0.115384615384615, 0, 0, 0, 0,
0, 0, 0, 0, 0.115384615384615, 0, 0, 0, 0.115384615384615, 0,
0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0.115384615384615,
0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0.115384615384615, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0.115384615384615, 0.115384615384615,
0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0.0506329113924051,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0506329113924051,
0, 0, 0, 0.0506329113924051, 0, 0, 0, 1, 0, 0, 0, 0.0506329113924051,
0.0506329113924051, 0.0506329113924051, 0, 0, 0, 0, 0.0506329113924051,
0, 0, 0.0506329113924051, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
1, 0, 0, 0, 0, 0, 0, 0, 0.0506329113924051, 0, 0, 0, 0, 0, 0,
0, 0, 0.0506329113924051, 0, 0, 0, 0.0506329113924051, 0, 0,
0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0.0506329113924051,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0506329113924051, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 0, 0, 0.0506329113924051, 0.0506329113924051,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0.0961538461538462, 0, 0.0961538461538462,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0.0961538461538462,
0, 0, 0, 0.0961538461538462, 0, 0, 0, 0, 0, 1, 1, 0.0961538461538462,
0.0961538461538462, 0.0961538461538462, 0, 0, 0, 0, 0.0961538461538462,
0, 0, 0.0961538461538462, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0961538461538462,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0.0961538461538462, 0, 0, 0, 0,
0, 0, 0, 0, 0.0961538461538462, 0, 0, 0, 0.0961538461538462,
0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0961538461538462,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0961538461538462, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0961538461538462, 0.0961538461538462,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0,
1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0.0448717948717949, 0, 0.0448717948717949,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0.0448717948717949,
0, 0, 0, 0.0448717948717949, 0, 0, 0, 0, 0, 0, 0, 0.0448717948717949,
0.0448717948717949, 0.0448717948717949, 0, 0, 0, 0, 0.0448717948717949,
0, 0, 0.0448717948717949, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0.0448717948717949,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0.0448717948717949, 0, 0, 0, 0,
0, 0, 1, 0, 0.0448717948717949, 0, 0, 0, 0.0448717948717949,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0448717948717949,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0448717948717949, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0.0448717948717949, 0.0448717948717949,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0448717948717949, 0, 0.0448717948717949,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0448717948717949,
0, 0, 0, 0.0448717948717949, 0, 0, 0, 0, 0, 0, 0, 0.0448717948717949,
0.0448717948717949, 0.0448717948717949, 0, 1, 0, 0, 0.0448717948717949,
0, 0, 0.0448717948717949, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0448717948717949,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0.0448717948717949, 0, 0, 1, 0,
0, 0, 0, 0, 0.0448717948717949, 1, 0, 0, 0.0448717948717949,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0.0448717948717949,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0448717948717949, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0.0448717948717949, 0.0448717948717949,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0), .Dim = c(173L, 31L), .Dimnames = list(
    c("TCGA-AB-2803", "TCGA-AB-2805", "TCGA-AB-2806", "TCGA-AB-2807",
    "TCGA-AB-2808", "TCGA-AB-2810", "TCGA-AB-2811", "TCGA-AB-2812",
    "TCGA-AB-2813", "TCGA-AB-2814", "TCGA-AB-2815", "TCGA-AB-2816",
    "TCGA-AB-2817", "TCGA-AB-2818", "TCGA-AB-2819", "TCGA-AB-2820",
    "TCGA-AB-2821", "TCGA-AB-2822", "TCGA-AB-2823", "TCGA-AB-2824",
    "TCGA-AB-2825", "TCGA-AB-2826", "TCGA-AB-2828", "TCGA-AB-2830",
    "TCGA-AB-2832", "TCGA-AB-2833", "TCGA-AB-2834", "TCGA-AB-2835",
    "TCGA-AB-2836", "TCGA-AB-2837", "TCGA-AB-2838", "TCGA-AB-2839",
    "TCGA-AB-2840", "TCGA-AB-2841", "TCGA-AB-2842", "TCGA-AB-2843",
    "TCGA-AB-2844", "TCGA-AB-2845", "TCGA-AB-2846", "TCGA-AB-2847",
    "TCGA-AB-2848", "TCGA-AB-2849", "TCGA-AB-2851", "TCGA-AB-2853",
    "TCGA-AB-2854", "TCGA-AB-2855", "TCGA-AB-2856", "TCGA-AB-2857",
    "TCGA-AB-2858", "TCGA-AB-2859", "TCGA-AB-2860", "TCGA-AB-2861",
    "TCGA-AB-2862", "TCGA-AB-2863", "TCGA-AB-2865", "TCGA-AB-2866",
    "TCGA-AB-2867", "TCGA-AB-2868", "TCGA-AB-2869", "TCGA-AB-2870",
    "TCGA-AB-2871", "TCGA-AB-2872", "TCGA-AB-2873", "TCGA-AB-2874",
    "TCGA-AB-2875", "TCGA-AB-2877", "TCGA-AB-2879", "TCGA-AB-2880",
    "TCGA-AB-2881", "TCGA-AB-2882", "TCGA-AB-2884", "TCGA-AB-2885",
    "TCGA-AB-2886", "TCGA-AB-2887", "TCGA-AB-2888", "TCGA-AB-2889",
    "TCGA-AB-2890", "TCGA-AB-2891", "TCGA-AB-2895", "TCGA-AB-2896",
    "TCGA-AB-2897", "TCGA-AB-2898", "TCGA-AB-2899", "TCGA-AB-2900",
    "TCGA-AB-2901", "TCGA-AB-2903", "TCGA-AB-2904", "TCGA-AB-2908",
    "TCGA-AB-2909", "TCGA-AB-2910", "TCGA-AB-2911", "TCGA-AB-2912",
    "TCGA-AB-2913", "TCGA-AB-2914", "TCGA-AB-2915", "TCGA-AB-2916",
    "TCGA-AB-2917", "TCGA-AB-2918", "TCGA-AB-2919", "TCGA-AB-2920",
    "TCGA-AB-2921", "TCGA-AB-2924", "TCGA-AB-2925", "TCGA-AB-2927",
    "TCGA-AB-2928", "TCGA-AB-2929", "TCGA-AB-2930", "TCGA-AB-2931",
    "TCGA-AB-2932", "TCGA-AB-2933", "TCGA-AB-2934", "TCGA-AB-2935",
    "TCGA-AB-2936", "TCGA-AB-2937", "TCGA-AB-2938", "TCGA-AB-2939",
    "TCGA-AB-2940", "TCGA-AB-2941", "TCGA-AB-2942", "TCGA-AB-2943",
    "TCGA-AB-2944", "TCGA-AB-2946", "TCGA-AB-2948", "TCGA-AB-2949",
    "TCGA-AB-2950", "TCGA-AB-2952", "TCGA-AB-2954", "TCGA-AB-2955",
    "TCGA-AB-2956", "TCGA-AB-2959", "TCGA-AB-2963", "TCGA-AB-2964",
    "TCGA-AB-2965", "TCGA-AB-2966", "TCGA-AB-2967", "TCGA-AB-2969",
    "TCGA-AB-2970", "TCGA-AB-2971", "TCGA-AB-2972", "TCGA-AB-2973",
    "TCGA-AB-2975", "TCGA-AB-2976", "TCGA-AB-2977", "TCGA-AB-2978",
    "TCGA-AB-2979", "TCGA-AB-2980", "TCGA-AB-2981", "TCGA-AB-2982",
    "TCGA-AB-2983", "TCGA-AB-2984", "TCGA-AB-2985", "TCGA-AB-2986",
    "TCGA-AB-2987", "TCGA-AB-2988", "TCGA-AB-2990", "TCGA-AB-2991",
    "TCGA-AB-2992", "TCGA-AB-2993", "TCGA-AB-2994", "TCGA-AB-2995",
    "TCGA-AB-2996", "TCGA-AB-2998", "TCGA-AB-2999", "TCGA-AB-3000",
    "TCGA-AB-3001", "TCGA-AB-3002", "TCGA-AB-3005", "TCGA-AB-3006",
    "TCGA-AB-3007", "TCGA-AB-3008", "TCGA-AB-3009", "TCGA-AB-3011",
    "TCGA-AB-3012"), c("Offset", "CEBPA", "DNMT3A", "FAM5C",
    "FLT3_ITD", "FLT3_TKD", "IDH1", "IDH2", "KIT", "KRAS", "NPM1",
    "NRAS", "PHF6", "PTPN11", "RAD21", "RUNX1", "SMC1A", "SMC3",
    "STAG2", "TET2", "TP53", "U2AF1", "WT1", "complex", "minus5_5q",
    "minus7q", "plus8_8q", "plus21", "t_15_17", "t_8_21", "inv16_t16_16"
    )))

```

# Answer
A design matrix is a matrix of values of the grouping variable.

ANOVA needs such a matrix to know which samples belong to which group. Since `limma` performs an ANOVA it needs this design matrix.

1. You can create it using the `model.matrix()` method.
2. Afterwards you can tell `limma` which groups you want to compare.

For point 2 you define a contrast matrix defining the comparisons of interest by using the `makeContrasts()` method.

It is well explained in this manual's section 9.2
> 2 votes

---
Tags: rna-seq, limma
---


title: STAR vs Bowtie2

text:
# Question
Title: STAR vs Bowtie2
What is the fundamental difference between STAR and Bowtie(2). Specifically, what is the difference in their final output (regardless of run-time differences, speed, memory usage etc.). Both seem to work on the same kind of input \*.fastq files without discriminating what they are associated with and generate an output version of the \*.fastq files. Would my downstream analysis be incorrect if I use Bowtie2 where STAR was supposed to be used, and vice versa? Thank you.

# Answer
> 6 votes
> What is the fundamental difference between STAR and Bowtie(2).

STAR for mapping spliced (i.e. with introns) short RNA-seq reads against a genome. Bowtie2 for mapping short reads without splicing.

> Would my downstream analysis be incorrect if I use Bowtie2 where STAR was supposed to be used (?)

Yes. If you use bowtie2 to map spliced short reads against a genome, STAR-based downstream analysis is likely to be wrong.

> ... and vice versa?

If you map genomic reads with STAR, the result will be suboptimal. Use the right tool for the right task.

# Answer
> 3 votes
The seeds are clustered together by proximity to a selected set of ‘anchor’ seeds in the 2nd phase of the STAR algorithm, but this •clustering/stitching/scoring process• is not present in Bowtie2.

mapping by performing a two-step process

---
Tags: sequence-alignment, bowtie2, star
---


title: nested samples in design DEseq2

text:
# Question
Title: nested samples in design DEseq2
I'm running a DEG analysis with Deseq2 by specifing the following design that includes 12 samples from 6 pts in 4 different conditions:

```
coldata <- data.frame("sample"=c(
  "T58", "T86", "T87",
  "D58", "D86", "D87",
  "C62", "C80", "C84",
  "C62", "C80", "C84"
))
coldata$Pts <- c(
  "m58", "m86", "m87",
  "m58", "m86", "m87",
  "c62", "c80", "c84",
  "c62", "c80", "c84"
)
coldata$class <- c(
  "Treated", "Treated", "Treated",
  "Distal", "Distal", "Distal",
  "CLEFT", "CLEFT", "CLEFT",
  "CRIGHT", "CRIGHT", "CRIGHT"
)

```

I would model the differences among the conditions by considering that the same individuals are nested in different conditions

I tried the following code

```
ddsMat <- DESeqDataSetFromMatrix(
  countData = as.matrix(matrix),
  colData = coldata,
  design = ~ class+Pts
)

```

but this produces the famous error

> "the model matrix is not full rank..."

If I reduce the previous dataset to

```
coldata2 <- as.data.frame(c(
  "T58", "T86", "T87",
  "D58", "D86", "D87"
))
coldata2$Pts <- c(
  "m58", "m86", "m87",
  "m58", "m86", "m87"
)
coldata2$class <- c(
  "Treated", "Treated", "Treated",
  "Distal", "Distal", "Distal"
)

```

then everything goes well and I'm able to create a `DESeq` object with the following design:

```
ddsMat <- DESeqDataSetFromMatrix(
  countData = as.matrix(matrix[,1:6]),
  colData = coldata2,
  design = ~ class+Pts
)

```

Please can you provide to me some suggestion on how to solve the problem?

Thank you

# Answer
> 2 votes
It is actually simple. The first two conditions are nested with the Pt's starting with m and the other two conditions are nested with those PT's starting with c. That means you cannot run this design. The only solution is to subset the data into two chunks and analyze separately.

# Answer
> 1 votes
The DESeq2 documentation explains this error as follows:

> While most experimental designs run easily using design formula, some design formulas can cause problems and result in the DESeq function returning an error with the text: “the model matrix is not full rank, so the model cannot be fit as specified.” There are two main reasons for this problem: either one or more columns in the model matrix are linear combinations of other columns, or there are levels of factors or combinations of levels of multiple factors which are missing samples.

In your case, I think the error is popping up because there are too many defined values in the "Pts" variable. DESeq2 doesn't require you to separately define replicates as a variable; it assumes that different data columns are taken from different biological replicates.

Getting the models right is tricky, and needs a deeper discussion about what the goal of analysis should be. My guess is that the following should work, but only if you're not interested in comparing expression results from different patients:

```
coldata$Pts <- c(
  "m", "m", "m",
  "m", "m", "m",
  "c", "c", "c",
  "c", "c", "c"
)
coldata$class <- c(
  "Treated", "Treated", "Treated",
  "Distal",  "Distal",  "Distal",
  "CLEFT",   "CLEFT",   "CLEFT",
  "CRIGHT",  "CRIGHT",  "CRIGHT"
)

```

This may not be what you want to do; other ways to deal with this are explained in the DESeq2 documentation on the Bioconductor website:

https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#model-matrix-not-full-rank

# Answer
> 1 votes
check if `Pts` and `class` are of `factor` class or not by using `is.factor(coldata$Pts)` and `is.factor(coldata$class)`

if it is not then convert it into `factor` by

`coldata$Pts = factor(coldata$Pts)`

`coldata$class = factor(coldata$class)`

Then try to design your matrix (A design matrix is a matrix of values of the grouping variable). I think this would work as it works with

`design <- model.matrix(~0+coldata$Pts+coldata$class)` in `limma` package.

You can also check in this RNA-seq analysis workflow of `DESeq2` where they also used nested samples for their design matrix.

# Answer
> 0 votes
Your original was,

```
ddsMat <- DESeqDataSetFromMatrix(countData = as.matrix(matrix), # colData = coldata, design = ~ class+Pts)

```

What happens when you use

```
ddsMat <- DESeqDataSetFromMatrix(countData = as.matrix(matrix), colData = coldata, design = ~ class+Pts)

```

Thus `#` is removed. The suspicion is `#` is causing the code to fall over. `#` is used for comments in `R` so the presence of this within a line of code will cause the following code `colData = coldata, design = ~ class+Pts` to not be interpreted because it can't be seen.

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


# Question
Title: How to reduce the occupied RAM when you are dealing with a very sparse matrix in a single-cell Experiment in R?
I'm dealing with a very large and sparse dataset and the first issues I met occurred when I tried to use quickCluster that reported me this error:

```
                'cannot allocate vector of size 156.6 Mb'

```

So, given that I cannot wait to change the RAM of my computer and I can't afford to use a cluster, I want to rely on some other strategies like some package that would allow me to handle sparse matrices. I'm thinking about sparseM but given that I don't know well this package I'd like to know how to shrink the ram allocation for these kind of matrices. Any suggestion will be very appreciated!

Link to the dataset

# Answer
> 5 votes
Ah, looks like I can't even procrastinate on StackExchange anymore without seeing work-related stuff. Oh well.

Anyway, the other answers and comments are way off. **scran** has supported sparse matrices for years, ever since we switched over to the `SingleCellExperiment` class as our basic data structure. `quickCluster` does no coercion to dense format unless you tell it so explicitly, e.g., with `use.ranks=TRUE` (in which case you're asking for ranks, so there's little choice *but* to collapse to a dense matrix).

You don't provide an MWE or your session information, but this is how it rolls for me:

```
# Using the raw counts in the linked dataset. Despite being
# called a CSV, it's actually space delimited... typical.
library(scater)
mat <- readSparseCounts("GBM_raw_gene_counts.csv", sep=" ")

# Making an SCE just for fun. Not strictly necessary for
# this example, but you'll find it useful later.
sce <- SingleCellExperiment(list(counts=mat))

library(scran)
system.time(clust <- quickCluster(sce))
##   user  system elapsed
##  3.170   0.174   3.411

```

This is running on my laptop - 16 GB RAM but I'm definitely not using all of it. I only go full throttle when I'm working on some real data, e.g., the 300k HCA bone marrow dataset. Check out the book for more details.

Session info below, I don't know quite enough SO-fu to know collapse it.

```
R version 4.0.0 Patched (2020-04-27 r78316)
Platform: x86_64-apple-darwin17.7.0 (64-bit)
Running under: macOS High Sierra 10.13.6

Matrix products: default
BLAS:   /Users/luna/Software/R/R-4-0-branch/lib/libRblas.dylib
LAPACK: /Users/luna/Software/R/R-4-0-branch/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets
[8] methods   base

other attached packages:
 [1] scran_1.16.0                scater_1.16.1
 [3] ggplot2_3.3.2               SingleCellExperiment_1.10.1
 [5] SummarizedExperiment_1.18.1 DelayedArray_0.14.0
 [7] matrixStats_0.56.0          Biobase_2.48.0
 [9] GenomicRanges_1.40.0        GenomeInfoDb_1.24.2
[11] IRanges_2.22.2              S4Vectors_0.26.1
[13] BiocGenerics_0.34.0

loaded via a namespace (and not attached):
 [1] beeswarm_0.2.3            statmod_1.4.34
 [3] tidyselect_1.1.0          locfit_1.5-9.4
 [5] purrr_0.3.4               BiocSingular_1.4.0
 [7] lattice_0.20-41           colorspace_1.4-1
 [9] vctrs_0.3.1               generics_0.0.2
[11] viridisLite_0.3.0         rlang_0.4.6
[13] pillar_1.4.4              glue_1.4.1
[15] withr_2.2.0               BiocParallel_1.22.0
[17] dqrng_0.2.1               GenomeInfoDbData_1.2.3
[19] lifecycle_0.2.0           zlibbioc_1.34.0
[21] munsell_0.5.0             gtable_0.3.0
[23] rsvd_1.0.3                vipor_0.4.5
[25] irlba_2.3.3               BiocNeighbors_1.6.0
[27] Rcpp_1.0.4.6              edgeR_3.30.3
[29] scales_1.1.1              limma_3.44.3
[31] XVector_0.28.0            gridExtra_2.3
[33] dplyr_1.0.0               grid_4.0.0
[35] tools_4.0.0               bitops_1.0-6
[37] magrittr_1.5              RCurl_1.98-1.2
[39] tibble_3.0.1              crayon_1.3.4
[41] pkgconfig_2.0.3           ellipsis_0.3.1
[43] Matrix_1.2-18             DelayedMatrixStats_1.10.0
[45] ggbeeswarm_0.6.0          viridis_0.5.1
[47] R6_2.4.1                  igraph_1.2.5
[49] compiler_4.0.0

```

# Answer
> 3 votes
Not a direct solution but some workarounds:

* As far as I know, `Seurat` can work with sparse matrices.
* The particular function of `scran` that you are using eats up quite some memory. I believe it is needed for the "normalization" step (that is how I used it anyway). Whereas the scaling normalization performed by this function is superior to crude "log normalization", you can give a try with the latter, which is far less computationally-intensive (does not do clustering). Seurat, once again, can help with this.
* You can downsample your data to the extent that it fits your into your RAM.
* You can give a try with Python. More and more single cell packages are written in Python, to some extent because of the problem you have experienced. For example `Scanpy` output is comparable to that of `Seurat`, although I am not sure if you can use scaling normalization with `Scanpy`.

# Answer
> 1 votes
Essentially you've hit a RAM bottleneck and the calculation will slow to zero, or in this instance refuse to go forward. The way to do this normally is to parallelize the calculation across the cores of your machine. This will likely remove the RAM bottlenack, dont ask me the computer archetectural reasons for why it works - but its works.

However, my knowledge of R is minimal. I wouldn't know how to parallelize an R calculation. It is certainly doable in Perl and Python, but the calculation needs to be written to ensure parallelisation.

The other way in is to reconfigure your calculation to remove sparse matrices OR find someone doing NGS where they've confgured their machine around heavy RAM.

Looking at your calculation I don't quite understand why you need to use a specific package, it looks like unsupervised machine learning, PCA - tSNE that sort of thing and you don't need a given package to do that you just need to vectorise the inputs. If you worked out the statistical components of `Scran` then there are a few extremely strong R statisticians/bioinformaticians on site who would have zero problem in replicating this within a few lines of code. Its not hard in Python's Sci-kit learn either. At a guess they perform PCA and resolve it via tSNE and this give nice clear clusters.

GCP is free for 3 months so it will cost you zero in context to a single calculation.

# Answer
> 0 votes
Not a direct solution to your problem because it does not use SCE. But if you (and others who ended on this post) are facing the RAM bottleneck when analysing single-cell data, then **Scarf** might be useful.

Here is the Github link: https://github.com/parashardhapola/scarf

Documentation: https://scarf.readthedocs.io/

Paper: https://www.nature.com/articles/s41467-022-32097-3

---
Tags: rna-seq, scrnaseq, bioconductor, single-cell, scran
---


# Question
Title: How to manage memory contraints when analyzing a large number of gene count matrices? I keep running out of RAM with my current pipeline
I have several hundred scRNA-seq count matrices, each from a different sample. For my other dataset containg a few dozen samples, I simply merged everything together into one Seurat object, but that won’t work here as far as I can tell. When I try to merge them in the typical way, I get the error:

```
Error in .cbind2Csp(x, y) :
  Cholmod error 'problem too large' at file ../Core/cholmod_sparse.c, line 89

```

This basically means I ran out of RAM. I even ran it on the cluster with 1 TB of RAM, and I got the same error. This is unsurprising, as memory requirements increase exponentially in relation to matrix size. It seems that this happens when the *merge()* function attempts to convert all the sparse matrices to dense matrices. Here is my code:

```
seurat_object_list <- "path/to/matrices" |> list.dirs(recursive = FALSE) |> lapply(\(matrix_folder) {
    matrix_folder |> Seurat::Read10X() |> Seurat::CreateSeuratObject()
})
seurat <- merge(seurat_object_list[[1]], seurat_object_list[2 : length(seurat_object_list)])

```

The Seurat object list is created with no issue, but it always fails on the merge. I tested this code with a smaller subset of the matrices, and it works just fine. It seems that I just have too many cells to process in this manner.

At this point, It looks like I have a few options.

1. I could hack together some way to force the matrices into one Seurat object. But I don't know if that would even be useful, since I would likely keep bumping my head on the same memory contraints whenever I try to do anything with the combined object, such as running *SCTransform()* or *FindMarkers()*.
2. I could divide the samples into groups that are small enough to merge, and process them in the normal way. But this imposes limitations on the types of comparisons/analyses I could perform.
3. I could keep the matrices separate and analyze them iteratively. This would solve the memory problems, but I haven't found any info online about anyone analyzing count matrices this way. I imagine that I would lose a lot of out-of-the-box functionality that comes with the standard single-cell R libraries. For example, I'm not sure how I would normalize all the gene counts against each other as separate objects. This option seems like it might require a lot of re-inventing the wheel and a lot of code.

I like option 3 the best, but I'm not sure where to begin. I am still pretty new to bioninformatics, so I would love to hear some input from those who are more experienced. How do you deal with memory constraints when dealing with a large number of cells/samples? Are there any standard ways to deal with this problem? Is there some tutorial/doc online that I failed to find? Any advice you could offer is greatly appreciated.

# Answer
> 1 votes
Scarf has been developed to specifically solve this issue. Scarf handles large-scale scRNA-Seq/scATAC-Seq/CITE-Seq datasets that will otherwise not fit into memory. Unfortunately, this is not available in R yet.

Here is the Github link: https://github.com/parashardhapola/scarf

Documentation: https://scarf.readthedocs.io/

Paper: https://www.nature.com/articles/s41467-022-32097-3

Please feel free to create an issue on Github for you use case (if you cant already find a solution in the documentation).

/PD

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


# Question
Title: How to extract clusters from a phylogenetic tree?
I am using Biopython Phylo and RDkit Cluster to obtain a UPGMA tree from a distance matrix of 5k x 5k entries. How can I extract the taxa names within a given clade(s)?

For example, Bio.Phylo can return plain-text representation of a tree, so I think I can parse it by the number of tabs:

```
Phylogeny()
     Clade()
         Clade(branch_length=0.06)
             Clade(branch_length=0.102, name='A')
             Clade(branch_length=0.23, name='B')
         Clade(branch_length=0.4, name='C')
...

```

or, ASCII format:

```
          _____________ A
  _______|
_|       |_______________________________ B
 |
 |_______________________________________________________ C

```

Which, I think is more difficult to parse but possible. Either way, those are not established methods for this kind of task, so I would like a hint.

# Answer
The example given a 3 taxa tree, isn't an ideal example. What is required is a raw treefile in a common format (nexus, philip, newark). The example below uses ETE3 for a 5 taxa tree A,B,C,D,E. Its ROOTED - very important to know the direction.

ETE3 is providing a means of traversing down the tree from the root, however the coder needs to interrogate the path it has taken.

1. Tree file imported
2. Ladderizing ... so you know the horizontal direction of the tree (see \*)
3. Ask the ETE3 to search down the nodes ... each child is a node from the root.
4. Parse the output and identify the clade.

The two outputs from point 3 are either it encounters a node - in this case it returns all the children below it (taxa names). The alternative is that it encounters a terminal node or a leaf, like a species name or isolate. In this case the taxa name is returned.

The caveat is the terminal node (taxa label or leaf) is always bifurcation and the taxa names are split into \[1\] and \[0\] on the ETE3 object (see code) ... you have to interrogate both.

Thats what the OP is looking for and this part of the code is marked **'Part 1'**

Parsing **Part 2** is more complicated, because usually the coder wants the clade a particular taxa (or leaf) belongs to. Essentially, the coder needs to ask what the leaf name is, if the target leaf name is found then the clade is then identified, i.e. the other stuff next to it.

The code is below. I agree if you don't think trees, this is a rushed explanation.

```
from ete3 import Tree

# Part 1
tree = Tree(intree, format=1)
tree.ladderize(direction=1)
clade1_ = tree.children[0].children[0].children[0].name
leaf1_ = tree.children[0].children[0].children[0].children[0].name
leaf2_ = tree.children[0].children[0].children[0].children[1].name

# Part 2
if leaf1 == 'A' or leaf2 == 'A':
    CladeA = clade1_
if leaf2 == 'B' or leaf2 == 'B':
    CladeA = clade1_

```

\*, ladderizing that can be important for classifying \[0\] and \[1\] in ETE3 ... its a weird one in trees 'cause it shouldn't matter
> 2 votes

---
Tags: phylogenetics, python, phylogeny, clustering
---


# Question
Title: Can I use a tool for highly variable gene selection for selecting highly variable "bins" from an ATACSeq dataset?
I'm working with ATACSeq data from multiple tissue/cell types. The data is binned in 1 megabase bins. I'd like to identify the bins that are "highly variable" across the different tissue types. I can't seem to find a tool to do this, like there exists for identifying highly variables genes in a gene expression setting (e.g Scanpy's `scanpy.pp.highly_variable_genes`). I was wondering if I could use the gene selection tools for my data, and if there are any other suggestions?

# Answer
I don't think is a good idea to apply the gene selection method on the ATAC-seq genomic regions because of the high sparsity of ATAC-seq data.

A common way is to select the most frequent regions (like the Signac does), or to discard the house-keeping regions (like the SnapATAC does).

And you can also use the episcanpy, a epigenomics version of scanpy.
> 0 votes

---
Tags: atac-seq
---


title: bash script for interproscan REST

text:
# Question
Title: bash script for interproscan REST
hey i have the problem with bash script. It should run interproscan REST Im green in bash scripting. I found this script on the web:

```
#!/bin/bash
i=1
waitevery=30

mkdir -p out

for j in $(find `pwd` -type f -name "*.fa")
do
    echo "Iteration: $i; File: $j"
    filename=$(basename "$j")
    python3 iprscan5.py \
        --goterms \
        --pathways \
        --email=blabla@protonmail.com \
        --outfile=out/${filename} \
        --outformat=tsv \
        --quiet \
        $j & (( i++%waitevery==0 )) && wait
done

```

out:

```
Iteration: 1; File: /home/marcin/interproscan-5.57-90.0-64-bit/interproscan-5.57-90.0/test/xtr.fa
interpro_batch_submit.sh: 18: i++%waitevery==0: not found
Iteration: 1; File: /home/marcin/interproscan-5.57-90.0-64-bit/interproscan-5.57-90.0/test/xts.fa
interpro_batch_submit.sh: 18: i++%waitevery==0: not found
Iteration: 1; File: /home/marcin/interproscan-5.57-90.0-64-bit/interproscan-5.57-90.0/test/xtt.fa
interpro_batch_submit.sh: 18: i++%waitevery==0: not found
Iteration: 1; File: /home/marcin/interproscan-5.57-90.0-64-bit/interproscan-5.57-90.0/test/xtu.fa
interpro_batch_submit.sh: 18: i++%waitevery==0: not found
Iteration: 1; File: /home/marcin/interproscan-5.57-90.0-64-bit/interproscan-5.57-90.0/test/xtv.fa
interpro_batch_submit.sh: 18: i++%waitevery==0: not found
Iteration: 1; File: /home/marcin/interproscan-5.57-90.0-64-bit/interproscan-5.57-90.0/test/xtw.fa

```

Edit:

I think I forgot to add `--sequence`, but something was still wrong after it.

# Answer
First of all *never* do `for foo in $(find ...)`. That is very bad practice as it can't handle even slightly odd file names like those with a space. See bash pitfall #1.

Next, you don't need `pwd`, just run `find -type f -name "*.fa"` which will search in the current directory by default and unlike `find `pwd`...` , won't break if your path contains whitespace.

That's just general best practice, the real issue here is that you are probably running the script like this:

```
sh /path/to/script.sh

```

Right? And you are using a Debian-based system like Ubuntu, right? On such systems, `/bin/sh` is actually not bash but a very basic POSIX compatible shell called `dash` and `dash` doesn't understand the `(( i++%waitevery==0 ))` syntax.

Since you have a shebang line (`#!/bin/bash`), you can simply make your script executable (`chmod a+x /path/to/script.sh`) and then run the script directly like this:

```
/path/to/script.sh

```

Or, alternatively, run it explicitly with `bash` instead of `sh`:

```
bash /path/to/script.sh

```

Either should work and if your file names and paths are sane, the issues I mentioned above won't be relevant. But to keep things nice and clean, I would also change the script to (note that this assumes you are using Linux and GNU `find`):

```
#!/bin/bash
i=1
waitevery=30

mkdir -p out

find -type f -name "*.fa" -print0 |
  while IFS= read -r -d '' j; do
    echo "Iteration: $i; File: $j"
    filename=$(basename "$j")
    python3 iprscan5.py \
        --goterms \
        --pathways \
        --email="blabla@protonmail.com" \
        --outfile=out/"$filename" \
    --outformat=tsv \
    --quiet \
    "$j" & (( i++%waitevery==0 )) && wait
done

```
> 1 votes

# Answer
I'd guess that the specific error that's coming up is because there are spaces missing from the wait line. This seems to be a problem in the source script, i.e.:

```
$j & (( i++%waitevery==0 )) && wait

```

should probably be:

```
$j & (( i++ % waitevery == 0 )) && wait

```

but for me (following @terdon's reproduction advice via `dash -c "i=9; ((++i%5==0)) && echo 'worked!'"`) that only changes the thing that is not found from `i++%waitevery==0` to `i++`, suggesting that the post-increment operator is not supported in the shell. Post-increment *is* supported in bash, suggesting that this script is being run (as @RamRS commented) in an unexpected way.

Try running the script directly:

```
./interpro_batch_submit.sh

```

Or explicitly using bash:

```
bash ./interpro_batch_submit.sh

```

... but see @terdon's answer for a better approach to doing this.
> 1 votes

---
Tags: bash
---


# Question
Title: Help with performing PCA analysis on data from 3 different datasets of normalized counts data for RNA-seq experiment
So, I have limited knowledge of R but I need to do a PCA analysis of 3 different datasets of gene expression as a result of combined growth or mono-culture growth.

The 3 different datasets I performed DESeq2 analysis on are as follows:

Dataset 1: This was sequenced by Ion Torrent and it was single-end reads. Samples that were compared were biological quadruplicates of yeast species A grown as pure culture versus species A grown in direct contact with yeast species B. This was all done for 1 time point @ 24h.

Dataset 2: This was sequenced by Illumina MiSeq and paired-end reads were generated, but this dataset had two parts to it.

Samples that were compared were biological quadruplicates of yeast species A grown as pure culture versus species A grown in direct contact with yeast species B. This dataset contains 2 time points @ 2h and 24h.

For the other part of this, samples that were compared were biological quadruplicates of yeast species A grown as pure culture versus species A grown via indirect (metabolic exchange possible) contact with yeast species B. This dataset contains 2 time points @ 2h and 24h.

Dataset 3: This was sequenced by Illumina MiSeq and paired-end reads were generated, but this dataset had two parts to it.

Here, two biological repeats of samples of yeast species A grown as pure culture versus two biological repeats of species A grown in direct contact with yeast species B, under AEROBIC conditions were performed. This dataset contains 1 time point @ 48h.

Then one biological repeat of sample of yeast species A grown as pure culture versus two biological repeats of species A grown in direct contact with yeast species B, under ANAEROBIC conditions were performed. This dataset contains 1 time point @ 48h.

I wrote the results out of the normalized counts per sample combined with the DESeq2 results as follows:

```
res <- results(dds)
table(res$padj<0.05)
res <- res[order(res$padj), ]
resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)), by="row.names", sort=FALSE)
names(resdata)[1] <- "Gene"
head(resdata)
write.csv(resdata, file="diffexpr-results.csv")

```

So the issue I am having is, because these are different datasets and conditions differ across them there are different numbers of genes expressed in each dataset that are not present across all datasets. So I need to combine all the datasets to only compare the genes found in all the datasets mentioned and then perform a PCA analysis on those genes only.

An example of the setup of the excel files I have, without the other info that you get with DESeq2 analysis e.g. columns `baseMean`, `log2FoldChange`, `lfcSE`, `stat`, `pvalue` and `padj`, is as follows:

| Gene\_ID | Sample1\_PC | Sample2\_PC | Sample3\_PC | Sample4\_PC | Sample5\_Mixed | Sample6\_Mixed | Sample7\_Mixed |
| --- | --- | --- | --- | --- | --- | --- | --- |
| YL09980 | 1240,016 | 759,683 | 985,322 | 1025,855 | 5002,654 | 4200,963 | 4002,855 |
| YLW9980 | 589,456 | 452,236 | 630,228 | 523,821 | 2325,556 | 2566,987 | 1800,633 |

So these tables are just the normalized counts of the genes per sample.

I was told that the data in the table has to be transposed basically, so that the genes are the variables and the samples are the rows. I am just not sure how to do this in R.

I was also told there is a good chance that the noise resulting from the use of different platforms would likely cause these datasets to separate from each other based on that, making this analysis meaningless or not allowing for comparison of gene expression that is conserved across the datasets as a result of mixing the two species during growth.

However, I just need to do this analysis and find that out as soon as possible, as this is a bottleneck to my progress at this point. Please any assistance would be greatly appreciated. I am currently teaching myself more about "tidyverse" in an effort to try and do this, but I am on a deadline with this.

Edit: So I did the PCA plot using the following function:

```
data_total.pca <- prcomp(data_total[,c(2:32)], center = TRUE, scale. = TRUE)

```

and this was the output I got:

```
Importance of components:
                          PC1    PC2    PC3    PC4     PC5     PC6     PC7     PC8     PC9    PC10    PC11    PC12    PC13    PC14
Standard deviation     4.2180 2.5400 1.9188 1.2866 0.65393 0.47799 0.40826 0.38879 0.35585 0.34915 0.27036 0.18549 0.14335 0.11830
Proportion of Variance 0.5739 0.2081 0.1188 0.0534 0.01379 0.00737 0.00538 0.00488 0.00408 0.00393 0.00236 0.00111 0.00066 0.00045
Cumulative Proportion  0.5739 0.7820 0.9008 0.9542 0.96800 0.97537 0.98074 0.98562 0.98971 0.99364 0.99600 0.99711 0.99777 0.99822
                          PC15    PC16    PC17    PC18    PC19    PC20    PC21    PC22    PC23    PC24    PC25    PC26    PC27    PC28
Standard deviation     0.10153 0.09239 0.08407 0.07765 0.07536 0.05978 0.04883 0.04556 0.04406 0.03960 0.03764 0.03378 0.03167 0.02897
Proportion of Variance 0.00033 0.00028 0.00023 0.00019 0.00018 0.00012 0.00008 0.00007 0.00006 0.00005 0.00005 0.00004 0.00003 0.00003
Cumulative Proportion  0.99855 0.99883 0.99906 0.99925 0.99943 0.99955 0.99963 0.99969 0.99976 0.99981 0.99985 0.99989 0.99992 0.99995
                          PC29    PC30    PC31
Standard deviation     0.02695 0.02126 0.02085
Proportion of Variance 0.00002 0.00001 0.00001
Cumulative Proportion  0.99997 0.99999 1.00000

```

So according to the tutorial I am following a general rule of thumb is that the first 3 PCs should explain at least 70% of the variance. The issue is my PCA plot looks weird all of the points are squished up at the x = 0; y = 0 axes of my PCA plot. So I think I need to transpose the data so that the genes are the columns (variables) and the samples are the rows (objects). So I tried this:

```
data_total_t <- transpose(data_total)
rownames(data_total_t) <- colnames(data_total)
colnames(data_total_t) <- rownames(data_total)

```

And this worked, but now I need to exclude row 1 from the function I used to plot the PCA as the PCA cannot have character (gene name) info in it.

Any help on how I can do this please?

# Answer
I think you are making this harder than you have to. Why can't you follow the code used in the vignette for making PCA plots? It doesn't even involve getting all the way to the DE genes.
> 1 votes

# Answer
> So I need to combine all the datasets to only compare the genes found in all the datasets mentioned and then perform a PCA analysis on those genes only.

You can merge all your datasets for common genes in all datasets and perform `PCA`.

> I was told that the data in the table has to be transposed basically, so that the genes are the variables and the samples are the rows. I am just not sure how to do this in R.

Your choice of transposing data will depend on whether you want to show variability between genes or samples. However you can get results for both individuals (genes) as well as variables (samples), using the `prcomp()`, `factoMineR` and `factoexta` packages in 'R'. Usually I like to use `$rotation` values to plot the variability between samples using `ggplot2` package. You can use this PCA Tutorial to learn more.

> I was also told there is a good chance that the noise resulting from the use of different platforms would likely cause these datasets to separate from each other based on that, making this analysis meaningless or not allowing for comparison of gene expression that is conserved across the datasets as a result of mixing the two species during growth.

This is true as the datasets are from different platforms, there may be some noise or batch effects. However, you can plot the 'PCA' and decide for yourself whether it is good or not.
> 1 votes

---
Tags: r, rna-seq, deseq2, pca, meta-analysis
---


# Question
Title: Beginner with DESeq2 having issue with analysis
I am trying to run DESeq2 analysis on my featureCounts output counts table.

I ran the following code:

```
library("DESeq2")

samplefile<-"C:/Users/RGB King/Documents/R/~/R/DESeq2_analysis/sample_table.txt"

coldata<-read.table(file=samplefile, sep="\t", header=FALSE, col.names = as.factor(c("","treatment")), row.names=1)

countfile<-"C:/Users/RGB King/Documents/R/~/R/DESeq2_analysis/feature_counts_test_AvABC.txt"
countdata<-read.table(file=countfile, sep="\t", header=FALSE, row.names = 1)

countdata<-countdata[-c(1:5)]

dds <- DESeqDataSetFromMatrix(
    countData = countdata,
    colData = coldata,
    design = ~ treatment
)

```

I then get the following error:

```
Warning message:

In DESeqDataSet(se, design = design, ignoreRank) :

  some variables in design formula are characters, converting to factors

```

The countdata table looks as follows:

```
|------|Variable1|Variable2|Variable3|Variable4|Variable5|Variable6|Variable7|Variable8|
|GeneID|A_1.bam--|A_2.bam--|A_3.bam--|A_4.bam--|AB_1.bam-|AB_2.bam-|AB_3.bam-|AB_4.bam-|
|THI3--|:12:-----|:34:-----|:8:------|:15:-----|:6:------|:22:-----|:41:-----|:33:-----|
|THI7--|:17:-----|:78:-----|:9:------|:23:-----|:11:-----|:54:-----|:19:-----|:7:------|

```

The coldata looks like this:

```
|--------|treatment|
|A_1.bam-|:A:------|
|A_2.bam-|:A:------|
|A_3.bam-|:A:------|
|A_4.bam-|:A:------|
|AB_1.bam|:B:------|
|AB_2.bam|:B:------|
|AB_3.bam|:B:------|
|AB_4.bam|:B:------|

```

So after I got the error I found posts on a forum that stated that DESeq2 only recognises "integer" class input values for the analyses. Therefore, I ran the following code:

```
sapply(countdata,class)

```

This then showed me that the data in the countdata table was in "character" class, so to fix this I ran the following:

```
chars <- sapply(countdata, is.character)
countdata[ , chars] <- as.data.frame(apply(countdata[ , chars], 2, as.integer))

```

But this then caused the values in the GeneID row to become "NA". So my questions are how can I change the class of just the GeneID row back to what it was? Will this then allow the analysis to run correctly? And in my code for the countdata I do specify "row.names =1" so why does this not keep my values in the GeneID row separate from the rest of the table?

# Answer
There is nothing wrong with your colData. It is encoded as character rather than factor, and DESeq2 is turning it into factors, that is all that the warning (warning != error) tells. I never needed to convert my matrix to load featureCounts output into DESeq2. Just load the output with the header, use `rownames()` to move the gene column to rownames, remove all the non-relevant columns (that is usually 1-7) and then feed the remaining one into DESeq2. I have no featureCounts output at hand right now. If you need code then please provide suitable example data via the `dput` function for easy copy/pasting.
> 1 votes

# Answer
As the warning message is suggesting, you should convert the `class` of your `treatment` column in `coldata` from `character` to `factor`.

Use this:
`coldata$treatment = factor(coldata$treatment)`
> 1 votes

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


# Question
Title: What is the difference between AnnotationDbi/org.Mm.eg.db and biomaRt/Mus.musculus for converting to gene symbols?
I am interested in the differences between AnnotationDbi/org.Mm.eg.db and biomaRt/Mus.musculus. They yield the same results as seen in the code below.

```
library(edgeR)
library(biomaRt)
library(Mus.musculus)
library(AnnotationDbi)
library(org.Mm.eg.db)

url <- "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file"
utils::download.file(url, destfile="GSE63310_RAW.tar", mode="wb")
utils::untar("GSE63310_RAW.tar", exdir = ".")
files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt", "GSM1545538_purep53.txt",
           "GSM1545539_JMS8-2.txt", "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
           "GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt", "GSM1545545_JMS9-P8c.txt")
for(i in paste(files, ".gz", sep=""))
  R.utils::gunzip(i, overwrite=TRUE)

files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt",
           "GSM1545538_purep53.txt", "GSM1545539_JMS8-2.txt",
           "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
           "GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt",
           "GSM1545545_JMS9-P8c.txt")

x <- readDGE(files, columns=c(1,3))
geneid <- rownames(x$counts)

genes <- AnnotationDbi::select(org.Mm.eg.db, keys=geneid, columns="SYMBOL", keytype="ENTREZID")
dim(genes)
genes2 <- biomaRt::select(Mus.musculus, keys = geneid, columns = c("SYMBOL"), keytype = "ENTREZID")
dim(genes2)

```

**Do they use different databases? Is one preferred or recommended over the other?**

# Answer
`AnnotationDBi` uses `AnnotationDb` objects which is the virtual base class for all annotation packages. It contain a database connection and is meant to be the parent for a set of classes in the Bioconductor annotation packages. Examples of classes derived from the `AnnotationDb` object include (but are not limited to): `ChipDb`, `OrgDb`, `GODb`, `OrthologyDb` and `ReactomeDb`. The database from one of these classes is used to extract annotations. `org.Mm.eg.db` belongs to `OrgDb` class.

`Mus.musculus` represents a collection of annotation packages that can be used as a single object. It combines different information about same organisms (Mus musculus) into one place. It also uses `org.Mm.eg.db` as a dependency. So, both uses same database to extract `Gene symbols`. That is the reason they are producing same results. There is helpful documentation available on Bioconductor.
> 2 votes

---
Tags: r, differential-expression, edger
---