# Question
Title: What is the pvalue from anova(glm1, glm2, test = “Chisq”) telling me in R when the glms have a negative binomial distribution?
Having constructed two nested negbin glms in R, one has a dropped factor, We can do an `anova` with `test = "Chisq"` to see if there is a significant difference between the models. But what does the Pvalue as a result of this actually tell me?

```
anova(mnegbin1, mnegbin2, test = "Chisq")
Likelihood ratio tests of Negative Binomial Models

Response: count
                        Model    theta Resid. df    2 x log-lik.   Test    df      LR stat. Pr(Chi)
1 origin + substrate + sample 18.46502      1232       -2459.886
2          substrate + sample 18.46502      1232       -2459.886 1 vs 2     0 -1.500666e-11       1

```

# Answer
> 2 votes
When you fit a model, you need to check whether it's overly determined, meaning some terms cannot be determined. We can look at one of your samples, for example only for RB5 and Nocardia\_nova\_SH22a\_RB20\_RB56:

```
df = structure(list(X = c(226L, 227L, 228L, 229L, 230L, 231L, 232L,
233L, 234L, 235L, 236L, 237L, 238L, 239L, 240L, 241L, 242L, 243L,
244L, 245L, 246L, 247L, 248L, 249L, 250L, 251L, 252L, 253L, 254L,
255L, 256L, 257L, 258L, 259L, 260L, 261L, 262L, 263L, 264L, 265L,
266L, 267L, 268L, 269L, 270L, 1126L, 1127L, 1128L, 1129L, 1130L,
1131L, 1132L, 1133L, 1134L, 1135L, 1136L, 1137L, 1138L, 1139L,
1140L, 1141L, 1142L, 1143L, 1144L, 1145L, 1146L, 1147L, 1148L,
1149L, 1150L, 1151L, 1152L, 1153L, 1154L, 1155L, 1156L, 1157L,
1158L, 1159L, 1160L, 1161L, 1162L, 1163L, 1164L, 1165L, 1166L,
1167L, 1168L, 1169L, 1170L), sample = structure(c(1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Nocardia_nova_SH22a_RB20_RB56",
"RB56"), class = "factor"), substrate = structure(c(10L, 37L,
29L, 11L, 24L, 9L, 39L, 30L, 42L, 14L, 21L, 40L, 38L, 45L, 28L,
5L, 22L, 3L, 41L, 16L, 32L, 27L, 18L, 26L, 20L, 34L, 7L, 12L,
23L, 17L, 44L, 35L, 4L, 19L, 25L, 13L, 1L, 33L, 31L, 43L, 36L,
6L, 2L, 8L, 15L, 10L, 37L, 29L, 11L, 24L, 9L, 39L, 30L, 42L,
14L, 21L, 40L, 38L, 45L, 28L, 5L, 22L, 3L, 41L, 16L, 32L, 27L,
18L, 26L, 20L, 34L, 7L, 12L, 23L, 17L, 44L, 35L, 4L, 19L, 25L,
13L, 1L, 33L, 31L, 43L, 36L, 6L, 2L, 8L, 15L), .Label = c("acylsphingosine",
"agarose", "alcohol", "alginate", "alpha-glucan", "amylaceous polysaccharides",
"arabinogalactans", "aribinose", "beta-glucan", "cellulose",
"chitin", "cutin", "cyclomaltodextrin", "D-galactose", "D-glucan",
"D-glucose", "dermatan sulfate", "dextran", "fructans", "fructose",
"galactan", "galactose", "galacturonide", "glucan", "glucose",
"glucuronoside", "glycan", "glycans", "glycogen", "hemicellulose",
"hyaluronate", "L-rhamnose", "lacto-N-tetraose", "levan", "lignin",
"maltose", "mannan", "mannose", "pectin", "peptidoglycan", "pullulan",
"trehalose", "Vanillyl-alcohol", "xanthan", "xylose"), class = "factor"),
    count = c(2L, 0L, 2L, 2L, 0L, 1L, 2L, 2L, 3L, 0L, 0L, 1L,
    0L, 0L, 0L, 2L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
    0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
    0L, 0L, 0L, 3L, 1L, 3L, 3L, 0L, 3L, 3L, 5L, 2L, 0L, 0L, 4L,
    0L, 0L, 1L, 2L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,
    1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L,
    0L, 0L, 0L), origin = structure(c(1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L,
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("free", "isolate"
    ), class = "factor")), row.names = c(226L, 227L, 228L, 229L,
230L, 231L, 232L, 233L, 234L, 235L, 236L, 237L, 238L, 239L, 240L,
241L, 242L, 243L, 244L, 245L, 246L, 247L, 248L, 249L, 250L, 251L,
252L, 253L, 254L, 255L, 256L, 257L, 258L, 259L, 260L, 261L, 262L,
263L, 264L, 265L, 266L, 267L, 268L, 269L, 270L, 1126L, 1127L,
1128L, 1129L, 1130L, 1131L, 1132L, 1133L, 1134L, 1135L, 1136L,
1137L, 1138L, 1139L, 1140L, 1141L, 1142L, 1143L, 1144L, 1145L,
1146L, 1147L, 1148L, 1149L, 1150L, 1151L, 1152L, 1153L, 1154L,
1155L, 1156L, 1157L, 1158L, 1159L, 1160L, 1161L, 1162L, 1163L,
1164L, 1165L, 1166L, 1167L, 1168L, 1169L, 1170L), class = "data.frame")

```

Now we look at whether how it is distributed:

```
df$sample = droplevels(df$sample)
table(df$origin,df$sample)

          Nocardia_nova_SH22a_RB20_RB56 RB56
  free                               45    0
  isolate                             0   45

```

You can see, the samples are either free or isolate, and never both. So you cannot discern the effect of origin from samples. This means your model is over-determined and hence you see the error above.

To know the effect of origin, you have to fit, i used a poisson for this example because if i use only the subset above, the counts are not overdispersed hence glm.nb throws some error:

```
fit = glm(count~ origin+substrate,data=df,family=poisson)
summary(fit)

 head(coefficients(summary(fit)))
                           Estimate   Std. Error       z value   Pr(>|z|)
(Intercept)           -2.266968e+01 2.930556e+04 -7.735624e-04 0.99938279
originisolate          6.931472e-01 2.970443e-01  2.333481e+00 0.01962291
substrateagarose       2.282092e-06 4.144430e+04  5.506408e-11 1.00000000
substratealcohol       2.280854e-06 4.144430e+04  5.503420e-11 1.00000000
substratealginate      2.278865e-06 4.144430e+04  5.498621e-11 1.00000000
substratealpha-glucan  2.295736e+01 2.930556e+04  7.833790e-04 0.99937495

```

You can see that in this model, free is the reference, and isolate is compared to free in the "originisolate". The coefficient 6.931472e-01 means isolate samples are higher than free samples, by exp(6.931472e-01) times.

---
Tags: statistics, logistic-regression
---


# Question
Title: How do I get repeating regions from fasta file or directly from internet
If I go to ncbi site, apply `customize view -> Customize -> All features -> Update view` and then `ctrl + f : Alu` I can get all alu repeaing regions coordinates from this chromosome.

It looks like:

```
 repeat_region   complement(37490..37800)
                 /note="AluY"
                 /rpt_family="SINE/Alu"
                 /rpt_type=DISPERSED

```

How could I download all this coordinates and names of repeating regions from the internet to .txt file? Or some tool

# Answer
Click on the 'GenBank' dropdown menu at the top left corner of this page to see a graphical view of the sequence. You should see a 'Repeat region' track in the Sequence Viewer. Click on the 'Download' button in the top right corner of the Sequence Viewer, choose 'Download Track Data' to launch a dialog box where you can enter a range, select the 'Repeat region' track and download data in BED format. Unfortunately, the names for the features are not great.

If you want the names to be the same ones in the flat file, you can download the feature table and parse it. To download the feature table, click on the 'Send To' link at the top right corner of this page, choose 'Complete Record', 'File' and 'Feature Table' as 'Format'. After you download the feature table, you can parse it using `awk` as shown below to create a 4-column table with the data.

```
cat seq.ft \
  | awk 'BEGIN{FS="\t";OFS="\t"}{
    if(NF==3) {print s,e,f,q ; s=$1; e=$2; f=$3; q=""}
    else if(NF==2) {print s,e,f,q ; s=$1; e=$2; q=""}
    else if(NF==5) {q=q"|"$4"="$5} }' \
  > feat_table.tsv

```
> 1 votes

---
Tags: fasta
---


# Question
Title: Too many requests error when querying NCBI
I keep getting a "too many requests error" when querying the NCBI SRA database, even though I'm running less than 10 requests per second, and I have an API key, which supposedly should allow me to run 10 per second.

Here is my (Python) code:

```
import subprocess
import concurrent.futures
import time

bioprojects = [
    "PRJNA644722",
    "PRJNA644892"
]

all_metadata_generators = []

def fetch_metadata(bioproject):
    return subprocess.run(f"esearch -db sra -query '{bioproject}[bioproject]' | efetch -format runinfo", shell=True, capture_output=True, text=True)

with concurrent.futures.ThreadPoolExecutor() as executor:
    for i in range(0, len(bioprojects), 10):
        all_metadata_generators.append(executor.map(fetch_metadata, bioprojects[i:i+10]))
        time.sleep(1)

for metadata_generator in all_metadata_generators:
    for metadata in metadata_generator:
        print(metadata)

```

# Answer
As @vkkodali pointed out, the problem turned out to be that you shouldn't use multithreading.
> 0 votes

# Answer
The execution of your EDirect commands may be limited unless you are using an API Key, which you can create from your MyNCBI account. Your API Key can be used in one of several different ways.

**In an argument to an EDirect command:**

```
esearch -db nuccore -query 'some query' -api_key 12345

```

**In a UNIX environment:**

You can either run this just before you are running some edirect commands or add the following line to your `.bashrc` file.

```
export NCBI_API_KEY=12345

```

**In the environment, from a Perl script:**

```
$ENV{NCBI_API_KEY} = "12345"

```
> 2 votes

---
Tags: python, phylogenetics
---


# Question
Title: Bacterial DNA at the tail of transcriptome reads. What does that mean?
I am assembling a transcriptome obtained from the Internet. The transcriptome was extracted from a human cancer tissue that had been previously grafted into a mouse. I have detected that many assembled RNA sequences do not match precisely to their reference but have additional sequences at one of their tails containing bacterial DNA. It's not that there are reads pertaining to that bacterium (always the same one) but that the RNA sequences of the cell have incorporated fragments of bacterial sequences (the reads corroborate these insertions with high coverage), therefore I can't just filter out bacterial reads as there are not reads exclusively pertaining to bacteria but hybrids. My question is how can it be possible? What possible mechanisms could have brought these bacterial fragments to the RNA of the cell. I have read in some papers that there is speculation that bacteria could be involved in the genesis of cancer cells. That would be an explanation if it wasn't because there are also mouse RNA for which the same circumstance applies: they also have bacterial sequences at the tails. I have also read in papers that there is speculation with the possibility that cancer cells are more susceptible to get DNA from bacteria transferred to their DNA. I am speculating that is the reason in this case (and somehow that weakness was transferred to contiguous mouse cells) but I would like to know if there is another possibility I am missing. Thanks.

# Answer
> 1 votes
From time to time these situations are encountered, notably when PCR is being deployed in rare DNA samples.

In this instance the wiki states,

> Ralstonia solanacearum is an aerobic non-spore-forming, Gram-negative, plant pathogenic bacterium.".

I'd leave it. The simple explanation is there has been a bad mix up. The best explanation is the people loading the Illumina machine, designating the adaptors etc ... In this scenario it has nothing to do with bacterial genome integration, plant bacteria don't touch mice and aerobic bacteria wouldn't grow in cancer cells.

The second issue is what is a soil bacteria doing in a mouse cancer biological facility? For the staff this wouldn't be good news.

The singular biological possibility is its a Pseudomonas and some species are serious opportunistic pathogens of immune-weakened humans, but there are a lot of Pseudomonas species and this one doesn't ring any bells, BUT its mice so who knows? One species in particular is associated with hospital outbreaks, but a biomedical unit that isn't quite the same thing.

---

The pathogen in humans is called Pseudomonas aeruginosa. If hospitals get this infection it is BAD NEWS, for example burns victims are really susceptible and without intervention it kills them. If you have evidence that the bacterial equivalent in mice is the bacterial pathogen you have found that is a substantial step. Biological Safety Facilities (BSFs) DEFINITELY get infected. They hate it when this happens and can be very picky for fear of a contamination, but it certainly happens. If you have evidence of a major contaminent pathogen in the long run it would be of HUGE help to BSFs globally. No-one really bothers with the disease pathology of mice - because their mice and spend their entire lives going 'squeek' - but when they are models for cancer their status changes beyond recognition. So I'd erase my initial thoughts and suggest, what you have found is publishable and the application is BSF contamination - don't say that though say this could be the mouse equivalent of P. aeruginosa. Its certainly a publication and goodluck.

# Answer
> 1 votes
What exactly is the sequence? Are you sure it's not an artifact of the library prep?

5 seconds of googling turns up this:

https://support.illumina.com/content/dam/illumina-support/documents/documentation/chemistry\_documentation/experiment-design/illumina-adapter-sequences-1000000002694-13.pdf

This sequence:

ACCACCGAGATCTACA

is on the adapters.

So is this

GACGGCATACGAGAT

Bacterial contamination would not put bacterial sequence at the end of host transcripts! But adapter contamination sure could be.

Here's the index 2 adapter

AATGATACGGCGAC**CACCGAGATCTACAC\[i5\]TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG**

Here's your sequence

**CACCACCGAGATCTACAC**CTAAGCCT**TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG**GTATCAACGCAGAGTACGGG

Here's it's blast against nr

AC**CACCGAGATCTACACCTAAGCCTTCGTCGGCAGCGTCAGATGTGTATAAGAGACAG**GTATC

It's all adapter.

Whoever did the transcriptome assembly didn't trim for adapters first. That's all.

---
Tags: ngs, assembly, rna, cancer, bacteria
---


# Question
Title: Display a score or graph to visualize the degree of conservation of each residue from alignment data
I would like to display a score or graph to visualize the degree of conservation of each residue from the amino acid alignment data. If possible, I'd like to extract the parts of the game that scored above a certain level.

Is there a tool or source code (python) to do that?

Sorry for my bad English ...

# Answer
I would suggest a logo plot, e.g.:

I would take @Devon Ryan's suggestion to create a multiple sequence alignment (MSA), and plug that into that logo plot tool (or others, there are many like it).
> 1 votes

---
Tags: sequence-alignment
---


# Question
Title: Calculating Isoelectric point from a multifasta file
I need to calculate Isoelectric point from a multifasta file, is there any python code or web tool that allows me to do that?

# Answer
> 4 votes
You can use EMBOSS Pepstats for this. It takes a multi-fasta as input and produces a table that contains various statistics related to the protein, including isoelectric point.

---
Tags: fasta, proteins
---


# Question
Title: What is the best way to get a large number of RNA seq data from SRA in Python without being denied access
I have the following code to download data from the SRA using multithreading in Python. After running this a few times now (for testing purposes), I keep getting denied access to the data. Not sure how to fix this. In particular, `cat`ing the output files gives:

```
<?xml version="1.0" encoding="UTF-8"?><Error><Code>AccessDenied</Code><Message>Access Denied</Message><RequestId>cb0f0e98-cafb-1dd7-9b7b-d8c49756ec52</RequestId><HostId>QeDGVwBXYp61J0B4_OUTn7UsEsiQEec0n18DAeR0kaE</HostId></Error>

```

Here is my code:

```
import threading
import requests
import time

start = time.perf_counter()

class MyThread(threading.Thread):
    def __init__(self, url):
        threading.Thread.__init__(self)
        self.url = url
        self.result = None
        self.filename = url.split('/')[-1]
    def run(self):
        res = requests.get(self.url)
        with open(self.filename, 'wb') as f:
            f.write(res.content)
urls = [
        'https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos1/sra-pub-run-1/SRR000001/SRR000001.1',
        'https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos1/sra-pub-run-1/SRR000001/SRR000001.2',
        'https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos1/sra-pub-run-1/SRR000002/SRR000002.1',
        'https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos1/sra-pub-run-1/SRR000002/SRR000002.2']

threads = [MyThread(url, ) for url in urls]
for thread in threads:
    thread.start()
for thread in threads:
    thread.join()

finish = time.perf_counter()

print(f'Finished in {round(finish-start, 2)} second(s)')

```

# Answer
I once needed to download more than 1000 files from SRA, and I found a blog post that is very useful (https://reneshbedre.github.io/blog/fqutil.html). The idea is to use aspera to fetch SRA files. This is fast becasue of aspera. Then use fasterq-dump to convert SRA file to fastq. fasterq-dump is faster than fastq-dump partly because it does not compress the output file. I think compressing is the rate-limiting step here. After this, you can use pigz, which is a multi-threading version of gzip to compress the fastq files. Note that if you want to download files to a computing cluster, you can submit jobs for fasterq-dump and pigz. This makes things much faster. You cannot submit jobs for fetching SRA, because you need to connect to internet.
> 0 votes

---
Tags: python
---


title: DiffBind, diffferentially binding site

text:
# Question
Title: DiffBind, diffferentially binding site
I have data for 3 histone marks (2 for silencing and 1 for activation) each mark has three replicates.

when I run the `diffBind` package I have three contrast:

```
Silence1-Activation
Silence1-Silence2
Silence2-Activation

```

When I run dba.plotBox for any of contrasts I will have the same distribution of reads.like below:  I even tried to start with 1000 peaks of each one that showed higher fold change but again got the same pattern for plotbox. How can I see which peaks belong to which histone marks?

# Answer
> 1 votes
I got the answer from Rory Stark. the samples should have a unique ID. (https://support.bioconductor.org/p/132502/#132514)

---
Tags: sequence-alignment, bioconductor, sequence-analysis, phylogenetics, chip-seq
---


# Question
Title: How to fragment genomes into non-overlapping sequences of differing sizes?
I have bacterial chromosomes/plasmids that I have downloaded off refseq in a multi-fasta that I would like to fragment into non-overlapping fragments of between 1 kb and 15 kb.

I have tried googling for related tools but to no avail. Does anyone know how I would achieve this? Something very similar is done in this paper: https://link.springer.com/article/10.1007/s40484-019-0187-4

Thanks in advance!

# Answer
> 3 votes
You could use `numpy.random.randint` to sample integers from a uniform distribution, from 1000 to 15000 nt, extracting subsequences of that length from a linearized multi-FASTA file.

Here's some code:

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

import sys
import numpy as np

low = int(sys.argv[1])
high = int(sys.argv[2])

def process_record(header, sequence, low, high):
    next = 0
    fragment = 1
    while next < len(sequence):
        new_next = np.random.randint(low, high) + next
        if new_next > len(sequence) - high:
            new_next = len(sequence)
        sys.stdout.write('{}:{}\n{}\n'.format(header, fragment, sequence[next:new_next]))
        next = new_next
        fragment += 1

header = ""
sequence = ""
for line in sys.stdin:
    if line.startswith('>'):
        if len(sequence) > 0:
            process_record(header, sequence, low, high)
            sequence = ""
        header = line.rstrip()
    else:
        sequence += line.rstrip()

if len(sequence) > 0:
    process_record(header, sequence, low, high)

```

To use it:

```
$ ./13879.py 1000 15000 < in.fa > out.fa

```

Fragments will be disjoint (non-overlapping, contiguous) subsequences.

I add a check to determine if the last fragment of a FASTA record would be shorter than `low` (1kb) in length, and, if so, append that fragment to the final subsequence. It is possible for that last fragment to be up to 16 kb - 1 nt in length, given your parameters:

```
if new_next > len(sequence) - high:
    new_next = len(sequence)

```

You could change this logic to discard the last fragment (i.e. `break`), or simply allow fragments less than 1kb in size, depending on what you need.

**EDIT**

Here's a slight modification that should allow to you preserve fragment sizes within your bounds.

A cartoon may help illustrate the problem. Let's set up a FASTA record, which has been broken up into `n` fragments:

```
   1       2           n-2      n-1      n
|-----|---------|...|-------|---------|-----|

```

We may get unlucky and sample a random integer for fragment `n-1` (`frag(n-1)`) which makes the length of `frag(n)` less than `low` (1 kb in your example).

In my first suggested solution, `frag(n-1)` and `frag(n)` are appended to one another. This is why you get some fragments longer than `high` (15 kb).

Another option is to discard `frag(n)`, which is of length less than `low` (1 kb). Here's some code:

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

import sys
import numpy as np
import timeit

low = int(sys.argv[1])
high = int(sys.argv[2])
in_fn = sys.argv[3]

"""
Globals
"""
fragment_lengths = []

"""
Process a record
"""
def process_record(header, sequence, low, high):
    next = 0
    last = 0
    fragment = 1
    fragment_sequence = ""
    while last < len(sequence) - high:
        next = np.random.randint(low, high + 1) + last
        fragment_sequence += sequence[last:next]
        # Comment the sys.stdout.write and uncomment the fragment_lengths line, if measuring fragment length distributions
        sys.stdout.write('{}:{}-{}:{}\n{}\n'.format(header, last, next, fragment, sequence[last:next]))
        #fragment_lengths.append(next - last)
        last = next
        fragment += 1
    """
    Uncomment to compare input sequence with fragments
    """
    #sys.stderr.write('In > {}\nOut > {}\n'.format(sequence, fragment_sequence))

"""
Parse records
"""
def parse_records(in_fn, low, high):
    header = ""
    sequence = ""
    with open(in_fn, 'r') as ifh:
        for line in ifh:
            if line.startswith('>'):
                if len(sequence) > 0:
                    process_record(header, sequence, low, high)
                    sequence = ""
                header = line.rstrip()
            else:
                sequence += line.rstrip()
        if len(sequence) > 0:
            process_record(header, sequence, low, high)

def count_elements(seq):
    hist = {}
    for i in seq:
        hist[i] = hist.get(i, 0) + 1
    return hist

def wrapper(func, *args, **kwargs):
    def wrapped():
        return func(*args, **kwargs)
    return wrapped

def main(in_fn, low, high):
    """
    Comment below if testing the distribution of fragments
    """
    parse_records(in_fn, low, high)
    """
    Uncomment to test the distribution of fragments from
    running fns on same input 100k times
    """
    #wrapped = wrapper(parse_records, in_fn, low, high)
    #timeit.timeit(wrapped, number=100000)
    #print(count_elements(fragment_lengths))

if __name__ == "__main__":
    main(in_fn, low, high)

```

If `frag(n)` is needed, and you need to preserve the uniform distribution (or whatever fragment length distribution you want to model), then rejection sampling might be an option.

With rejection sampling, you keep fragmenting your input sequence, until one of the distributions of lengths gets you a set of fragments that entirely cover the input sequence.

If you need that, follow up and I'll take another look.

# Answer
> 1 votes
I've previously written my own script to fragment DNA sequences (input as either fastq or fasta files) into identical-length fragments, with a given overlap. I used this because there was a particular assembler that demanded equal-length sequences, and I was working with nanopore reads, but it's since been useful for other purposes as well.

Because it is most concerned about making lengths identical, there may be some forced overlap when setting the overlap to zero. The overlap will only be exactly zero when the sequence length is an exact multiple of the fragment size. Here's an example use:

## zero overlap

```
normalise_seqlengths.pl -fraglength 25 -overlap 0 oneRead.fq
@A00187:180:HLGNGDSXX:3:1101:8486:1031#0 [150 bp, 6 fragments, 0 overlap] 1:N:0:CGCAAGAA+ACTCCTAA
ATTGCGCAATGACAGCCAGGGGAGT
+
FFFFFFFFFFFFFFFFF:,FFFFFF
@A00187:180:HLGNGDSXX:3:1101:8486:1031#1:25-50  1:N:0:CGCAAGAA+ACTCCTAA
GGGAGCTGTCTGGGGTAACCTGGCT
+
:FFFF:F:FFFFF::FFFFF,::FF
@A00187:180:HLGNGDSXX:3:1101:8486:1031#2:50-75  1:N:0:CGCAAGAA+ACTCCTAA
CTAGGACCAGGTGGAAAAAGCACCT
+
F:FFFFFFFFFFFFFFFFFFFFFFF
@A00187:180:HLGNGDSXX:3:1101:8486:1031#3:75-100  1:N:0:CGCAAGAA+ACTCCTAA
TGGCTGGCCCAGCCCTGTTCTTCTC
+
:FF,FFFFF,:,FFF::F:FFF,:F
@A00187:180:HLGNGDSXX:3:1101:8486:1031#4:100-125  1:N:0:CGCAAGAA+ACTCCTAA
CTGGGGAGGGAAGGGCAGGTCCCCT
+
F:FFFFFFFFF,FFFFF:FFF,F,,
@A00187:180:HLGNGDSXX:3:1101:8486:1031#5:125-150  1:N:0:CGCAAGAA+ACTCCTAA
AATGCCTGAACCACATGGTGGGCCA
+
FFF,F:FFFFFFF:F,,,FFF,F,,

```

## partial overlap

```
normalise_seqlengths.pl -fraglength 63 -overlap 0 oneRead.fq
@A00187:180:HLGNGDSXX:3:1101:8486:1031#0 [150 bp, 3 fragments, 20 overlap] 1:N:0:CGCAAGAA+ACTCCTAA
ATTGCGCAATGACAGCCAGGGGAGTGGGAGCTGTCTGGGGTAACCTGGCTCTAGGACCAGGTG
+
FFFFFFFFFFFFFFFFF:,FFFFFF:FFFF:F:FFFFF::FFFFF,::FFF:FFFFFFFFFFF
@A00187:180:HLGNGDSXX:3:1101:8486:1031#1:43-106  1:N:0:CGCAAGAA+ACTCCTAA
CCTGGCTCTAGGACCAGGTGGAAAAAGCACCTTGGCTGGCCCAGCCCTGTTCTTCTCCTGGGG
+
FF,::FFF:FFFFFFFFFFFFFFFFFFFFFFF:FF,FFFFF,:,FFF::F:FFF,:FF:FFFF
@A00187:180:HLGNGDSXX:3:1101:8486:1031#2:87-150  1:N:0:CGCAAGAA+ACTCCTAA
CCCTGTTCTTCTCCTGGGGAGGGAAGGGCAGGTCCCCTAATGCCTGAACCACATGGTGGGCCA
+
FFF::F:FFF,:FF:FFFFFFFFF,FFFFF:FFF,F,,FFF,F:FFFFFFF:F,,,FFF,F,,

```

---
Tags: assembly, bacteria
---


# Question
Title: Chromatin accessibility level on promoters with different CpG ratio calculation
I want to categorize promoters based upon high, intermediate and low-CpG content (high-CpG-density promoters (HCPs), intermediate-CpG-density promoters(ICPs), and low-CpG-density promoters (LCPs)).

So the data I have is for promoter which i have annotated 1000+/- around TSS and taken them as something less than 1kb is promoter region and beyond distal.

In term of tool i have used chipseeker to annotate and do the above step.

So now if i have to find the CpG density as i have mentioned above how do i get that?

# Answer
Below I import the peaks call in narrowbed format from macs2 into a GRanges object. If you have a bed file, you can just use rtracklayer for it:

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

peaks_narrowbed = "https://www.encodeproject.org/files/ENCFF846JAO/@@download/ENCFF846JAO.bed.gz"

extraCols <- c(signalValue = "numeric", pValue = "numeric",
                          qValue = "numeric", peak = "integer")
gr <- import(peaks_narrowbed, format = "BED",extraCols = extraCols)

head(gr)

GRanges object with 6 ranges and 6 metadata columns:
      seqnames        ranges strand |
         <Rle>     <IRanges>  <Rle> |
  [1]     chr1   10049-10326      * |
  [2]     chr1 180686-181040      * |
  [3]     chr1 186648-186903      * |
  [4]     chr1 191293-191604      * |
  [5]     chr1 267872-268146      * |
  [6]     chr1 586110-586324      * |
                                                   name     score signalValue
                                            <character> <numeric>   <numeric>
  [1] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_1       135     7.10907
  [2] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_2       309    11.65888
  [3] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_3        81     5.40289
  [4] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_4        56     4.54981
  [5] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_5       145     7.39343
  [6] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_6       116     6.54035

```

Then we extract the sequences that are within these regions, you can of course redefine the region by extending it from the peak:

```
library(Biostrings)
SEQ = getSeq(Hsapiens,gr)
head(SEQ)

A DNAStringSet instance of length 6
    width seq
[1]   278 TAACCCTAACCCTAACCCTAACCCTAACCCTAAC...CCCCAACCCCAACCCCAACCCCAACCCCAACCC
[2]   355 CTAACCCCTAACCCCTAACCCTAACCCTACCCTA...CTCAGCCGGCCCGCCCGCCCGGGTCTGACCTGA
[3]   256 GCAACCACCTGAGCGCGGGCATCCTGTGTGCAGA...CCCCTTCTTTCCATTGGTTTAATTAGGAACGGG
[4]   312 CAATCAGCAGGGACCGTGCACTCTCTTGGAGCCA...CCCTGCCCCTGTCTCTTCCGTGCAGGAGGAGCA
[5]   275 TGGATGTGTGCATTTTCCTGAGAGGAAAGCTTTC...TTGTTTAGTTTGCGTTGTGTTTCTCCAACTTTG
[6]   215 GAGGCTGAAGGAGACTGATGTGGTTTCTCCTCAG...CCAACTTTGTGCCTCATCAGGAAAAGCTTTGGA

```

This is the CG, and you can decide how you want to bin your regions:

```
freq = dinucleotideFrequency(SEQ,as.prob=TRUE)
hist(freq[,"CG"],br=100)

```
> 2 votes

---
Tags: atac-seq
---


# Question
Title: Error given while trying to index a BAM file with Samtools Index - NO COOR?
I am currently working on my own Metagenomic pipeline, utilizing Bowtie 2 to map. Bowtie 2 outputs a SAM file, which I convert to a .BAM and sort it using Samtools. When I try to utilize Samtools to index my .BAM file it gives me this error:

```
[E::hts_idx_push] NO_COOR reads not in a single block at the end 1652 -1.

```

Here is my code:

```
source /opt/asn/etc/asn-bash-profiles-special/modules.sh
module load samtools/1.3.1

cd ~/gz_files/soapdenovo
r=20
###SAM to BAM and sort

samtools faidx soapdenovo.fa
samtools view -b -T soapdenovo.fa -o sample${r}bowtiemapping.bam sample${r}bowtiemapping.sam
samtools sort -n -T ~/gz_files/soapdenovo/temp -o sample${r}bowtiemapping_sorte$
samtools index sample${r}bowtiemapping_sorted.bam

################LETS MAKE THE ABUNDANCE TABLES###########################$
samtools view -c -f 4  sample${r}bowtiemapping_sorted.bam > sample${r}alignment$

samtools idxstats sample${r}bowtiemapping_sorted.bam >  sample${r}idxstats.txt

#git clone https://github.com/metajinomics/mapping_tools.git
#python mapping_tools/get_count_table.py *.idxstats.txt > contig_counts.tsv
##less contig_counts.tsv

```

Any help would be appreciated!

---

*Notes on pipeline*

* `samtools import` was replaced with `samtools view`
* `faidx`corrected
* Instructions used are here on github
* Suggestion for future trouble shooting is given here by @winni2k

# Answer
By now this will be of little interest to the original poster, but for future reference:

```
$ samtools sort -n -T ~/[etc] -o sample${r}bowtiemapping_sorted.bam [etc]
$ samtools index sample${r}bowtiemapping_sorted.bam
[E::hts_idx_push] NO_COOR reads not in a single block at the end 1652 -1.

```

This error message indicates (admittedly somewhat obscurely) that the reads are not sorted by coordinate — in particular, there are mapped reads appearing after the unplaced unmapped reads (those with RNAME `*`) that should be at the very end of the file.

This is because `samtools sort -n` has been used to sort the reads by name instead. Remove `-n` to sort by position, which is what is needed to prepare a BAM file for indexing with `samtools index`.
> 1 votes

---
Tags: samtools, bash, bugs
---


# Question
Title: What pitfalls exist with running FastQC on a bam file?
I was analysing a bad sequencing run of some RNA data using FastQC, I supplied it RNA-STAR-aligned bam files (and the human reference). The output of MultiQC had many more "unique reads" in the sequence counts table than there were reads in the FastQ files that were fed to RNA-STAR.

At first, I realized that all the FastQC stats were about "mappings", not reads (as they were labeled), but "unique mappings" doesn't seem to make sense either, because even if I consider the possibility that segments of reads were mapping in different places, the numbers don't make sense. Here is the worst case:

* Max read length after trimming: 47 bases
* Number of reads in the original FastQ file that went into STAR: 59M

Bam Stat Table Numbers:

* Total records: 391 million
* Unique mappings: 3.1 million
* Non-primary Hits: 335 million

FastQC Sequence Counts Bar graph:

* "Unique reads": 185 million
* "Duplicate reads": 206 million

While I found documentation stating that you can run FastQC on bam files, there were no caveats I could find about it needing to be unmapped reads.

Can anyone explain precisely why there can be 185 million "unique reads" when bam stats says there were 3 million unique mappings and the number of reads is 59M?

# Answer
> 1 votes
I don't see why you are trying to compare two totally different measurements.

When FASTQC says a read is unique, it means that no other read shares its sequence.

When STAR says a read is uniquely **mapped**, that means it maps to one and only place. That does not mean that there aren't a hundred more identical reads that map nicely to the exact same place.

You could have a very high level of PCR duplication of reads that nicely align to one and only one place. Or you could have a whole bunch of reads that differ slightly from each other who all map to 5 locations each.

FASTQC was designed for DNA sequencing, not RNASeq. Some of its warnings make no sense in the context of RNASeq. If I was doing RNASeq of a filtered set of pancreatic alpha cells, I would expect a hell of a lot of reads for glucagon, and a lot of them would be natural duplicates, and that would be biologically normal, and not a problem.

You could take the same batch of reads, and align them to one genome, and align them to a mix of two closely related genomes. The second will have far more non-uniquely mapped reads, but the FASTQC stats will not have changed.

That said, if you doing human or mouse RNASeq, 1% unique alignment means something is wrong. Find out where all those reads are aligning, you might have rRNA contamination. If you are doing an organism with a high level of duplication in its genome that might be normal.

---
Tags: sam, rna-alignment, fastqc
---


# Question
Title: downloaded software in Conda environment does not show in ls command
I installed gzip package in Conda environment, then to make sure the software was there, I used "ls" command, it was not in the list. Tried to install again, system says all the packages have already installed, Why? How do I check the list of software installed in the environment? Thanks!

# Answer
> 1 votes
Just type

```
conda list gzip

```

OR

```
conda list | grep gzip

```

OR if its installed in a given environment

```
conda activate myenv

```

or else just `activate myenv` system depending followed by either of the first two commands.

It is better to use `conda-navigate` a very nice and under used GUI, and would work in this particular instance. The details are found in an earlier post I made on esentially the same question here together with a pretty picture.

`ls` will only work if you type `ls ~/Anaconda3/bin | grep gzip` , but depends on the installation, e.g. you might have installed miniconda.

Having said all that `conda` is very cool, but whatever you do NEVER use `pip` or `pip3` from this moment onwards EVER, it messes up your conda and you'll need to `env` your way out of it.

---
Tags: conda
---


# Question
Title: Searching for start and stop codons for protein sequencing of contigs
I need to convert contigs into their respective protein sequences given a reference genome (i.e. I need to take a substring, whose position is already known on the string, and I need to locate the nearest start and stop codons).

This is tricky because sometimes the first position of triplet codon within the contig will not be a multiple of three and could be located in the intergenic region (i.e. between genes). The goal is to separate both coding and non-coding DNA.

This is my code so far:

```
from Bio.Seq import Seq
from Bio.Alphabet import generic_dna, generic_protein

start_codons = ['ATG']
stop_codons = ['TAG', 'TAA', 'TGA']

string = 'GG*TAG*CCAATT*ATG*AACGAA*TAG*GAC' #remove '*', just for visual
contigs = ['CCAA', 'TGAAC', 'GAA', 'GGAC']
positions = [5, 12, 17, 22] #position indices for each contig on string

extended_contigs = []
extended_position_contigs = []
intergenic_contigs = []
intergenic_position_contigs = []
for i in contigs:
    extended_contigs.append(#some code)
    extended_positions_contigs.append(#some code)
    intergenic_contigs.append(#some code)
    intergenic_positions_contigs.append(#some code)

```

I should get `extended_contigs = ['ATGAAC', 'ATGAACGAA']` and `extended_positions_contigs = [12, 17]`. These are the contigs that are located within genes. In order to code them into peptides, I need to reach back in the string until I find the start codon and extend the initial contig (e.g. `TGAAC` -\> `ATGAAC` and `GAA` -\> `ATGAACGAA`)

I should also get `intergenic_contigs = ['CCAA', 'GGAC']` and `intergenic_positions_contigs = [5, 22]`. When the missing code is run, the computer searches to the left of the string and finds a stop codon (e.g. `TAG`) before the start codon. Thus, the contig is located in between two genes and nothing needs to be added. These intergenic contigs are just stored in a new list.

My code continues:

```
prot_contigs = []
for i in extended_contigs:
    my_dna = Seq(i, generic_dna)
    my_prot = my_dna.translate()
    prot_contigs.append(str(my_prot))

```

Here, no new code needs to be added. After the above is run, `prot_contigs = ['MN', 'MNE']`.

The last step of code (which I need help with) will transform `prot_contigs` into `new_prot_contigs = ['MN', 'E']`.

How? If for any contig (e.g. `'TGAAC'`), the beginning or end is part of another codon (not perfect multiple of 3), the extra codon on either end will be maintained (e.g. `'MN'` stays `'MN'`). Otherwise, if the contig (e.g. `'GAA'`) perfectly matches a codon, anything that was added onto it will be removed (e.g. `'MNE'` becomes `'E'`).

I would try to solve the 2 portions of code myself, but I am unsure on how to take a position on a string (i.e. beginning points of contig) and look along the string to find the nearest start/stop codons, so I could determine the DNA's function and accurately sequence the protein coding contigs into peptides.

Any help would be appreciated!

# Answer
If your string is 'GG*TAG*CCAATT*ATG*AACGAA*TAG*GAC' then I would assume that you could separate this into coding and non-coding regions such as:

```
GG*TAG*CCAATT*    (non-coding)
ATG*AACGAA*TAG    (coding region)
*GAC              (non-coding)

```

It starts with a start codon ATG and ends with an end codon TAG.

Your first event position 12 does not make sense to me. What does * mean? It seems of that your sequence starts with GG which is not a triplet? Perhaps did you mean GGG?

I'm the author of node based DNA visualisation `AminoSee`. It has a bug in that it only renders half right: because DNA can be read backwards, and my program only does it *forwards*. I actually came to BioInformatics today to try and get my head around reading backwards and how that might work.

AminoSee only transcodes 3 base pairs at a time, and discards triplets that aren't, well, 3 base pairs longs. I'm self taught and so I guess unaware if it's possible to have any sequence that is not a multiple of 3. Perhaps it's possible. I assumed not. My program is more of an art piece that a science tool, but I plan to use it one day to perform a "visual diff" between two species (Humans and Chimp). See AminoSee DNA Visualisation over at aminosee.funk.nz
> 1 votes

---
Tags: proteins, python, phylogenetics, genomics
---


# Question
Title: How to install/import Biopython into Python 3.8/ PyCharm IDE
I'm new to Python/Biopython programming. I'd like to install Biopython into my PyCharm IDE but I've tried all instructions on this, but it's not confirming it. I have the available biopython file already extracted. There's a SyntaxError for 'install' in: pip install biopython What's the proper way of installing and importing biopython? And how does one confirm the installation.

# Answer
> 0 votes
Install Anaconda3 the instructions are here. Then making sure you are connected to the internet type

```
conda install -c anaconda biopython

```

Once Anaconda is installed don't attempt to use `pip`, if you do then you need to set an 'environment' in Anaconda and work from that environment.

I suspect you used `pip` or `pip3` via

```
pip install biopython

```

---
Tags: python, software-installation
---


title: Running MaxQuant on Linux

text:
# Question
Title: Running MaxQuant on Linux
I'm trying to run MaxQuant on a linux laptop, hower I'm constantly running into problems with MaxQuant crashing at different stages

I tried to use

* mono 6.8
* mono 5.4.1

on

* centos
* ubuntu

each with MaxQuant

* 1.6.3
* 1.6.7
* 1.6.11

(everything using docker), but I had no success so far.

I also tried to use the ahead of time compilation of mono, also without sucess. MQ mostly stalls at the "preparing searches" stage.

Has anyone managed to get MQ running properly on linux and can tell me, which linux distro/version, mono version and MQ version you are using?

Any other help is of course also welcome

# Answer
> 4 votes
There is a conda package for maxquant that should work on most platforms (I've used it on CentOS and I think it's been used on Ubuntu as well). There's also a docker container version of that that we (Bioconda) have posted to the biocontainers project that should also work.

Please note that large parts of a maxquant analysis are single threaded, regardless of what you specify. We've also had single analysis runs lasting over a month (even when using 90 threads), so the program can be very slow if your search space is large.

# Answer
> 2 votes
I have been running MaxQuant quite frequently on Linux. I would use conda to create a dedicated environment and install `mono`.

```
# create the environment
conda create -n maxquant -c conda-forge mono
# activate the environment
conda activate maxquant
# run any maxquant version
mono /path/to/maxuant/MaxQuantCmd.exe mqpar.xml

```

This installs a conda environment with `mono` which is a `.NET` framework for Linux. I had problems in the past running certain MaxQuant versions that way.

The `mqpar.xml` that is required can be created with the MaxQuant GUI. And then the path for input files needs to be contained correctly. The path to the *fasta* file and the *raw* file. I would create a template file and then parse these filenames into the specific locations.

```
# mqpar.xml
...
<fastaFiles>
  <FastaFileInfo>
     <fastaFilePath>__FASTA__</fastaFilePath>
     <identifierParseRule>>([^\s]*)</identifierParseRule>
     <descriptionParseRule>>(.*)</descriptionParseRule>
     <taxonomyParseRule></taxonomyParseRule>
     <variationParseRule></variationParseRule>
     <modificationParseRule></modificationParseRule>
     <taxonomyId></taxonomyId>
  </FastaFileInfo>
</fastaFiles>
...
<filePaths>
  <string>__RAW__</string>
</filePaths>
...

```

The combination of `maxquant=1.6.10.43` and `mono=5.14` worked and a few other versions. Some versions between 1.6 and 2.x did not work at all.

Specific versions of MaxQuant are also available directly via `conda` and there is no need to download the zip file. Instead, you can simply install it with:

```
conda install -c conda-forge -c bioconda maxquant mono=5.14

```

Recently MaxQuant 2 has been added, and not all previous versions are available and run maxquant with:

```
maxquant mqpar.xml

```

---
Tags: linux, proteomics, mass-spectrometry, maxquant
---


# Question
Title: Removing data from OrthoFinder output prefix
I am using OrthoFinder to find orthogroups between species. I am using the latest version unpacked through tar as per their instructions.

Orthofinder creates an output folder according to the current date, ie, Results\_Dec25. However, I want to integrate orthofinder in snakemake, so having a stable output directory is necessary. The -o option allows specification of an output prefix, but the dated folder is then created within.

Is it possible to force a given name to the output folder? Or is there a way in snakemake to allow ambiguous output directory names?

I know this would be possible via editing the source code of OrthoFinder, but this workflow needs to be reproducible without faff.

Many thanks

This question has already been asked here but remains unanswered.

Is there an option in Orthofinder to choose the output directory name?

# Answer
> 2 votes
The latest version solves the issue.

Source : https://github.com/davidemms/OrthoFinder/issues/403
Latest Version: https://github.com/davidemms/OrthoFinder/releases/tag/2.4.0

---
Tags: snakemake, orthofinder
---


# Question
Title: Human Herpes Virus, fusion motif in protein vs. SARS-CoV-2?
Is it possible that the fusion with cell propriety of Human Herpes virus, I see X14112.1 in GenBank is due to protein UL36 (tegument of virus) starting in linear position 2685 ???

```
H R R A R S L A

```

there are similarities with SarsCov2 virus protein S in name MN908947.3 in GenBank position 681.

```
P R R A R S V A

```

It is strange but the now (21-12-2020) the new Sars2 Great Britain strain is more near to that UL36 piece, it seems already:

```
H R R A R S V A

```

Herpes virus protein UL36 and SarsCov2 virus protein S\[1200 position\], have too

```
L Q E L G K

```

i find in internet they call it "Glutaminase kidney isoform, mitochondrial"

# Answer
I can answer that ... I'm afraid no. What you are asking is whether there is any similarity in their cellular receptor binding or "fusion loops" because the motifs look a bit similar.

These are very different viruses and you need to appreciate just how different they are,

*Human herpes* is a double stranded DNA virus resulting in a chronic infection. It has a vast genome of \>150kb, which includes cytokine mimics and proteins that we have no idea what they do. The protein you are referring to is likely a surface antigen in a virus covered in different surface antigens. However, this protein is involved in blocking the beta inferon response is considered a protease - not fusion loop and appears to have some sort of complex replication function involving the capsid. Finally, there is also the question of which human herpes 1 to 8. Here you are referring to HHV1 but for example HHV3 is chickenpox (okay its an acute not chronic infection) and HHV5 is cytomegliovirus (congenitcal infection).

*COVID-19* is a simple virus of a 30Kb genome of a single stranded positive sense RNA, causing an acute infection. I agree the Spike protein is much bigger than most RNA viruses, but the receptor-binding is tightly defined with ACE-2. It forms a classic RNA icosahedral virus, singly covered in Spike protein. I have no doubt it will go through conformational change like other RNA viruses in receptor binding but the protein involved is singularly involved in cell entry - thats its job.

**Summary** There are large differences between HHV1 and COVID-19 for,

1. Viruses; DNA double stranded vs RNA single stranded
2. Infections; chronic vs. acute
3. Transmission; vertical/integration? direct contact vs. respiratory
4. Ultimately very different proteins
5. Finally very different functions

At best this would be heavily scrutinised as a virological assumption.

*Bioformatics* If you really wanted to follow this up you could compare the structures around these motifs, because COVID-19 is close to SARS and assess this possibility ... but you would need to be careful about the conformational changes Spike goes through in receptor-binding.

**Advice** before investing large amounts of time in a bioinformatics question, it is well-worth understanding the biological rationale for the question. In this instance that rationale appears to be lacking. If you found a comparison for example with gamma-coronaviruses, which was absent from many beta-coronaviruses that would be interesting. Here there is simply too much possibility of a chance similarity and given the size HHV1 that cannot be easily overlooked. However it is good that you are thinking and questioning, if you apply the same techniques, e.g. psi-blast, protein-blast into different questions ... that can work.
> 2 votes

---
Tags: proteins, covid-19, sars-cov-2, motifs, coronavirus
---


# Question
Title: How can I restrict the search space when querying GEO?
I am using the following command to retrieve GEO Datasets entries for daf-2 c. elegans that have as their experiment type High Throughpout sequencing:

```
esearch -db gds -query "daf-2" | efilter -query "expression profiling by high throughput sequencing [DataSet Type] AND Caenorhabditis elegans [organism]"

```

However, among the results that I got, I noticed that I got this:

```
23. Impaired Insulin-/IGF1-Signaling Extends Life Span by Promoting Mitochondrial L-Proline Catabolism to Induce a Transient ROS-Signal
(Submitter supplied) Transcriptome profiling of three models with impaired insulin/IGF1 signaling. 1. Deep sequencing of endogenous mRNA from Caenorhabditis elegans N2 var. Bristol (wildtype) and daf-2(e1370) mutant; 2. Deep sequencing  of endogenous mRNA from murine embryonic fibroblasts (MEF)  wildtype and irs1-/- knockout; 3. Deep sequencing of endogenous mRNA from murine embryoinic fibroblast (MEF) insr+/- -lox and insr+/- knockout  Jena Centre for Systems Biology of Ageing - JenAge (www.jenage.de)
Organism:   Mus musculus; Caenorhabditis elegans
Type:       Expression profiling by high throughput sequencing
Platforms: GPL11002 GPL13776 14 Samples
FTP download: GEO (CSV) ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE36nnn/GSE36041/
SRA Run Selector: https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA151765
Series      Accession: GSE36041 ID: 200036041

```

As you can see, in the `Organism` field, there are both `Mus musculus` and `Caenorhabditis elegans`. How can I restrict my search so that I only get entries that exclusively work on `Caenorhabditis elegans`? One naive way to do this would be to include `NOT Mus Musculus` in the query, but of course that would mean that I could get something else too.

Another way I can think of solving this is just to write a script that would do the extra filtering using regexes, but I was wondering if there's a simpler solution using the e-utilities functionalities

# Answer
You can use the following query to first search GEO to get a list of GEO Series of interest, then find all linked SRA runs and perform a second filter step to keep only the SRA runs that satisfy another set of criteria as shown below.

```
esearch \
  -db gds \
  -query "daf-2[all fields] \
    AND expression profiling by high throughput sequencing [DataSet Type] \
    AND Caenorhabditis elegans [organism]" \
  | elink \
  -db gds -target sra \
  | esearch \
  -query "(#3) \
    AND Caenorhabditis elegans [organism]" \
  | efetch -format runinfo

```

Here, I am combining two independent queries using the `(#3)` term in the second query. More information about this is in the "Combining Independent Queries" subsection of "Searching and Filtering" section of the Entrez Direct documentation.
> 0 votes

---
Tags: geoquery
---


# Question
Title: Merge / Reconciliate several de novo transcriptome assemblies with different kmers
I am building a *De Novo* transcriptome reference assembly for an eukaryotic organism for which I have a genome.

I've created several assemblies with rnaSpades using different kmer sizes (19 to 69 with step 10). Now I would like to merge them into one final transcriptome.

How could I do that?

Is using a genome assembly reconciliation tool such as metassembler a good idea?

# Answer
> 2 votes
You may combine all your transcriptomes into a single file and then apply a clustering method to group very similar transcripts into a single one. For that purpose, you may try CD-HIT-EST or MMseqs2. For each identity threshold you are going to test, you may assess the final quality with BUSCO or by blasting against reference sequences.

# Answer
> 2 votes
A simple option would be to simply supply SPAdes with some contigs that you like with the `--trusted-contigs` option.

There exist tools specifically addressing the problem of merging transcriptome assemblies. I am unsure whether this is notably different from genome assemblies sufficiently that they are better than the genome assembly mergers. Here is a partial list:

1. transfuse
2. this paper
3. this paper
4. this paper
5. DRAP

This paper seems to have some more information on comparing tools, though it is not focused specifically on merging.

See also the SeqAnswers on this topic, and this discussion too.

---
Tags: assembly, transcriptome
---


# Question
Title: Adding treatment groups via metadata to Seurat object?
I'd like to add metadata to 6 individual Seurat objects so that after I merge the objects into one, I can later label or split by using these identifiers. For example, I'd like to append an age group and then interval across these 6 objects. For adding an interval, I've tried using the below:

> AddMetaData(control, metadata = 1hr, col.name = Interval)

However, this error will occur:

```
Error in .AddMetaData(object = object, metadata = metadata, col.name = col.name) :
  object '1hr' not found

```

Appreciate any help, thank you!

# Answer
> 0 votes
You should have fed `"1h"` to your function call but you fed `1h` instead. The latter is looking for an object with that name and as the error shows, there is no such object.

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


# Question
Title: How to get a list of genes corresponding to the list of SNPs (rs ids)?
Is there a way for me to get a list of genes given a list of SNP rs ids? I found several questions asked with a similar goal years ago, and the answers are always about using multiple online tools, and/or R programming.

I wonder if there is a simpler solution recently? (and hopefully, I can solve it within the scope of python. )

The relevant question I found (that are asked and answered more than five years ago.)

# Answer
> I wonder if there is a simpler solution recently? (and hopefully, I can solve it within the scope of python. )

A *simpler* solution, I don't know... but this is at least *one* Python solution using Biopython's ELink method via NCBI's Entrez E-utils.

The Biopython library is flexible enough (they have an in-depth tutorial worth reading) to modify the code below to to fit project-specific needs. For example, the mapping between snp id and gene id is lost with this bulk request method.

```
from Bio import Entrez
Entrez.email = "youremail@gmail.com"

snp_ids = ["rs8179065",
           "rs12899845",
           "rs8179066",
           "rs80357332",
           "rs80356920"]

record = Entrez.read(Entrez.elink(dbfrom="snp",
                                  id=",".join(snp_ids).replace('rs', ''),
                                  db="gene"))

```

This returns a `record` object

```
[{'LinkSetDb': [{'DbTo': 'gene', 'LinkName': 'snp_gene', 'Link': [{'Id': '374654'}, {'Id': '90381'}, {'Id': '672'}]}], 'ERROR': [], 'LinkSetDbHistory': [], 'DbFrom': 'snp', 'IdList': ['8179065', '8179066', '12899845', '80356920', '80357332']}]

```

which contains several iterable fields -- 'Link' has a list of the Gene **UID**'s that can be converted to names using `ESummary`.

```
for gene_id in record[0]['LinkSetDb'][0]['Link']:
    handle = Entrez.esummary(db="gene", id=gene_id['Id'])
    uid_record = Entrez.read(handle)
    handle.close()

    uid_summary = uid_record["DocumentSummarySet"]['DocumentSummary'][0]
    print (gene_id['Id'], uid_summary['Name'])
    # ...or something more interesting than print

374654 KIF7
90381 TICRR
672 BRCA1

```

An option to preserve discrete mapping, is to loop over each snp id and extract the gene(s):

*although I don't endorse the method below, especially for long lists of snps, you may flood the server with too many requests, maybe someone more familiar can comment on request limits or best practices for Entrez*

```
for snp_id in snp_ids:
    record = Entrez.read(Entrez.elink(dbfrom="snp",
                                  id=snp_id.replace('rs',''),
                                  db="gene"))
    results = record[0]['LinkSetDb'][0]['Link']
    for result in results:
        uid = result['Id']
        handle = Entrez.esummary(db="gene", id=uid)
        uid_record = Entrez.read(handle)
        handle.close()
        uid_summary = uid_record["DocumentSummarySet"]['DocumentSummary'][0]
        gene_name = uid_summary['Name']

        print(snp_id, uid, gene_name)

rs8179065 374654 KIF7
rs12899845 374654 KIF7
rs12899845 90381 TICRR
rs8179066 374654 KIF7
rs80357332 672 BRCA1
rs80356920 672 BRCA1

```
> 8 votes

# Answer
I can show you a simple way in R, using biomaRt. Let's say you have two snp ids you want to exmine.

```
SNPids <- c("rs431905509", "rs431905511")

library("biomaRt")

# To show which marts are available
listMarts()

# You need the SNP mart
mart <- useMart("ENSEMBL_MART_SNP")

# Find homo sapiens
listDatasets(mart)

# This will be the dataset we want to use
dataset <- useDataset("hsapiens_snp", mart=mart)

# Show available filters
listFilters(dataset)

# Now list all available attributes
listAttributes(dataset)

# To get the ensembl gene id belonging to the SNPs
getBM(attributes=c('refsnp_id', 'ensembl_gene_stable_id'),
      filters = 'snp_filter',
      values = SNPids,
      mart = dataset)

    refsnp_id ensembl_gene_stable_id
1 rs431905509        ENSG00000070371
2 rs431905509        ENSG00000100075
3 rs431905509        ENSG00000260924
4 rs431905511        ENSG00000145335

```
> 7 votes

# Answer
If you want a fast (binary) search with a locally-stored dataset, where you know exactly what SNPs and gene annotations are being queried, here's a way to set up files to do your own querying, which should scale up nicely if you have a lot of SNPs.

1) Get annotations of interest and write them into a sorted BED file.

For example, here's a way to get Gencode v26 gene annotations:

```
$ wget -qO- ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_26/gencode.v26.basic.annotation.gff3.gz \
    | gunzip -c \
    | convert2bed --input=gff - \
    > gencode.v26.bed

```

2) Get SNPs and write them into a text file sorted by SNP ID.

For `hg38`, for instance:

```
$ LC_ALL=C
$ wget -qO- http://hgdownload.cse.ucsc.edu/goldenpath/hg38/database/snp147.txt.gz \
   | gunzip -c \
   | awk -v OFS="\t" '{ print $5,$2,$3,($3+1) }' \
   | sort -k1,1 \
   > hg38.snp147.sortedByName.txt

```

3) Install pts-line-bisect:

```
$ git clone https://github.com/pts/pts-line-bisect.git
$ cd pts-line-bisect && make && cd ..

```

4) Run a query:

```
$ rs_of_interest=rs2814778
$ ./pts-line-bisect/pts_lbsearch -p hg38.snp147.sortedByName.txt ${rs_of_interest} \
   | head -1 \
   | cut -f2- \
   | bedmap --echo --echo-map-id-uniq --range 1000 - gencode.v26.bed \
   > answer.bed

```

The file `answer.bed` contains the genomic interval of the SNP of interest, along with the IDs of any Gencode v26 gene annotations within 1000 nt upstream and downstream of the SNP.

To repeat step 4 within a Python script, you can use `subprocess.check_output()`:

```
...
pts_lbsearch = os.path.join(project_dir, "pts-line-bisect", "pts_lbsearch")
if not os.path.exists(pts_lbsearch):
  print "Content-type:application/json\r\n\r\n"
  print json.dumps({'error' : 'could not locate pts_lbsearch [%s]' % (pts_lbsearch)})
  sys.exit(errno.ENOENT)

if term.startswith('rs'):
  query_fn = os.path.join(annotations_dir, "hg38.snp147.sortedByName.txt")
  cmd = "%s -p %s %s" % (pts_lbsearch, query_fn, term)
  try:
    res = subprocess.check_output(cmd, shell=True)
    if res and len(res) > 0:
      (chr, start, stop) = res.strip().split('\n')[0].split('\t')
    else:
      pass
  except subprocess.CalledProcessError as err:
    sys.stderr.write("Query failed!\n")
    sys.exit(errno.EIO)

#
# map 'chr', 'start', and 'stop' variables to annotations with
# another subprocess call (not shown -- use example above)
#
...

```
> 2 votes

# Answer
Have you tried using SnpEff and SnpShift (http://snpeff.sourceforge.net) for annotation of the SNPs?

Even if your desired organism available in their annotation database you can create them easily with the gene bank files. Other than annotation the tools features a bunch of other amazing functionalities.
> 1 votes

# Answer
using ncbi-efetch and **XSLT** with the following stylesheet:

```
<?xml version='1.0' encoding="UTF-8"?>
<xsl:stylesheet
    xmlns:xsl='http://www.w3.org/1999/XSL/Transform'
    xmlns:r="https://www.ncbi.nlm.nih.gov/SNP/docsum"
    version='1.1'
    >
<xsl:output method="text"  encoding="UTF-8"/>

<xsl:template match="/">
<xsl:apply-templates select="r:ExchangeSet"/>
</xsl:template>
<xsl:template match="r:ExchangeSet">
<xsl:apply-templates select="r:Rs"/>
</xsl:template>
<xsl:template match="r:Rs">
<xsl:apply-templates select="r:Assembly/r:Component/r:MapLoc/r:FxnSet"/>
</xsl:template>

<xsl:template match="r:FxnSet">
<xsl:text>rs</xsl:text>
<xsl:value-of select="../../../../@rsId"/>
<xsl:text>  ncbi:</xsl:text>
<xsl:value-of select="@geneId"/>
<xsl:text>  </xsl:text>
<xsl:value-of select="@symbol"/>
<xsl:text>
</xsl:text>
</xsl:template>

</xsl:stylesheet>

```

usage (group the rs# by 3):

```
$ echo -e "rs8179065\nrs12899845\nrs8179066\nrs80357332\nrs80356920\n" | xargs  -L 3 echo  | tr " " "," | while read R; do curl -sL "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=snp&id=${R}&retmode=xml" | xsltproc --novalid  transform.xsl -; done | uniq
rs8179065   ncbi:374654 KIF7
rs12899845  ncbi:90381  TICRR
rs8179066   ncbi:374654 KIF7
rs80357332  ncbi:672    BRCA1
rs80356920  ncbi:672    BRCA1

```
> 0 votes

# Answer
```
BiocManager::install('biomaRt')
require(biomaRt)
ensembl <- useMart("ENSEMBL_MART_SNP", dataset = "hsapiens_snp")
getBM(attributes=c("refsnp_id", "chr_name", "chrom_start", "chrom_end",
                   "allele", "mapweight", "validated", "allele_1", "minor_allele",
                   "minor_allele_freq", "minor_allele_count", "clinical_significance",
                   "synonym_name", "ensembl_gene_stable_id"),
      filters="snp_filter", values=c("rs6025","rs431905509"),
      mart=ensembl, uniqueRows=TRUE)

```
> 0 votes

---
Tags: database, gene, snp
---


# Question
Title: How to add cluster name to metadata in Seurat?
I'm working on a Seurat object and want to name the clusters according to 2 values alone (yes/no). So I want to add a new column to metadata and annotate the clusters (UMAP) with it.

```
head(B@meta.data)
        orig.ident nCount_RNA nFeature_RNA percent.mito RNA_snn_res.0.6 seurat_clusters
BC01_02       BC01   999789.2         9225 0.2385090196               6               6
BC01_03       BC01   999776.9         9568 0.2925109851               6               6
BC01_04       BC01   999595.1         7482 0.0009888404               6               6
BC01_05       BC01   999992.0         8301 0.3801925053               6               6

```

I want to add a new column imported from a csv file that contains the row.names same as the meta.data with the status column as yes and no.

```
head(status)
               status
BC01_02           yes
BC01_03           yes
BC01_04           yes
BC01_05           yes

```

After adding the column, I want to name the clusters with yes or no, how is it possible?

# Answer
To add the column `status` can try the following:

```
B@meta.data$status <- status[row.names(B@meta.data),]$status

```

To rename the clusters you can try the following:

```
B@meta.data$seurat_clusters <- B@meta.data$status

```

However, I think the plotting functions in Seurat do not use the cluster information from `B@meta.data`, rather they use information from `B@active.idents`. So, I guess you have to do the following instead:

```
B@active.idents <-  B@meta.data$status

```

This should work, but, I will verify and update this answer later.
> 1 votes

---
Tags: rna-seq, scrnaseq, seurat, sequence-annotation, clustering
---


# Question
Title: Can Maven (mass spectrometry software) read mzML files?
I did not find a documentation. Can Maven (mass spectrometry software) read mzML files?

# Answer
> 0 votes
No, Maven still can read only mzXML files.

---
Tags: mass-spectrometry
---


# Question
Title: Principle of TMT Tags in Multiplex Proteomics
I am new to proteomics research and analyzing mass-spec data. I am trying to learn more about the theory/principle of TMT tagging (a type of isobaric labeling) coupled with LC-MS/MS spectrometry. I understand that each tag has a

1) protein binding group (amine reactive ONLY?) 2) reporter group (so we know the tag came from which sample?) 3) normalisation group (to ensure that each tag has same mass?) 4) cleavage group (so tag can be cleaved from the peptide?)

Did I get the purpose of all the four parts in the tag?

Also, what does it mean by 4-plex, 10-plex, n-plex tag? And why do the tags have different isotopes?

If I could have some clarity on the workflow of TMT tagging and how it leads to analysis, it would be really helpful.

# Answer
TMT stands for tandem-mass-tag. The idea is to tag different samples with one of the tags then pool all samples together and run them through the mass spectrometer together. Then the reporter groups are split off and quantified (ideally) for each peptide individually. There are different methods with different numbers of isobaric labels. Isobaric means all the labels have the same mass therefore they behave exactly the same in the mass spectrometer at first. E.g. TMT11 means you have 11 different tags with identical masses, so you can pool and process 11 samples at a time.

The protein binding group attaches the label to the peptide. Then the tagged peptides from all samples are run through the chromatography and the spectrometer. At first, the MS instruments simply collects the tagged peptides with different m/z ratios at each point in time for a while. When enough peptides are collected reporter groups are cut off and the peptides then quantified. This cutting-off happens by increasing the internal energy of the collected peptides. Because the cleavage group has the weakest chemical bond it brakes first and releases the reporter group. Again, this step separates the reporter group from the rest. Finally, the cut-off reporter groups are quantified.

The different reporter groups don't have the same mass, otherwise you would not be able to differentiate between them in the MS instrument. But you need same masses before quantification to ensure that the same peptides are all collected together. Therefore, you need to balance out the different masses of the reporter groups. That is the job of the normalization group. So, the normalization group ensures that overall the tags have the same mass and the peptides with different tags are not treated differently in the MS instrument before quantification. Then for quantification the reporter groups (with different masses between each other) are cut off and the MS instrument can therefore detect and quantify them individually. There is a time dependency of the whole process which happens in multiple phases. The whole process makes it possible to compare the amounts of peptides coming from the different samples.

The chemistry in the chromatography column is quite complicated which makes it difficult to compare samples from different experiments. The presence of one peptide might enhance or suppress the signal of another peptide or impact the chromatographic retention time. Therefore, you should only compare within the same sample. More tags TMT11, TMT16 etc. allow you to run more samples in one run.

You can then use one (or more) channels (mass-tags), to tag a standard sample and then compare everything to this standard that is assumed to be identical over multiple experiments. This way you can compare even more samples than available tags. e.g. more than 11.
> 2 votes

---
Tags: proteins, proteomics, mass-spectrometry
---


# Question
Title: How does picard's MarkDuplicate handle unmapped reads?
Our BAM files are created according to a "lossless" alignment procedure \[\] from the Broad Institute `GATK` documenation and involves re-adding the unaligned/unmapped reads into an aligned BAM, using Picard's `MergeBamAlignment`.

The BAM files are produced in the end contain both the mapped and the unmapped reads. These files are then sorted with SortSam \[\]- so that the sort order in the header becomes:

`@HD VN:1.6 SO:coordinate`

How does `MarkDuplicates` handle the unmapped reads of a BAM file containing both unmapped and mapped?

**Note** `MarkDuplicates` seems to normally take the BAM's ordering into account, namely, it accepts arguments such as `--ASSUME_SORT_ORDER X`. However it's not specified whether reads without a position are ignored, or have to be compared with all other possible reads.

---

*Disclaimer*: I initially posted this question on the GATK forum \[\], but I'm reaching out to hopefully a broader audience.

Citations:

# Answer
From the Picard documentation:

> DUPLICATION METRICS: Metrics that are calculated during the process of marking duplicates within a stream of SAMRecords.

> UNMAPPED\_READS The total number of unmapped reads examined. (Primary, non-supplemental)

It won't alter the flags on these reads, but it will count them in the summary report it generates. You should be able to test this yourself with a small set of mapped + unmapped reads
> 2 votes

---
Tags: picard, markduplicate
---


# Question
Title: Standard for when to collapse phylogenetically uninformative nodes
I have been creating phylogenetic trees using RAxML (using Boostrap + ML) and I end up with a bestTree file, outputted in newick format. When I open up the tree in FigTree, I notice that it collapses many of the internal nodes with small values (1 x 10-6). I am working on visualizing the trees using ETE, which when I render the tree includes these small breakpoints. I was wondering if there is a standard for when to collapse nodes that are considered phylogenetically uninformative that I can apply. Thank you for any help.

# Answer
> 1 votes
If 1E6 is genetic distance, this means one nucleotide difference per million base pairs, which is of course incredibly small and bootstrapping at this level doesn't really work. So what I suspect is happening is that the bootstrap is so low it is not present on the tree. You can assign a cut-off bootstrap but anything below 75% is not interpretable in any case. What might be happening is you are producing the bootstrap "cladogram" which is showing a default branch length, but when \<1E6 genetic distance it is showing a trifurcation. This simply means the bootstrap of the consensus tree is so low it is collapsing the branch and this is a function of the very small genetic distance at this node.

---

Looking at your treefile I don't see any collapsed clades.

```
(Group_9:0.00000100000050002909,(Group_50:0.00507659947127924935, (.

```

The clade above has not been colapsed it is just a very small distance in a paraphyletic group. It looks fine to me, i.e. Group\_9 is 1E6. It could look 'collapsed' simply because it is so small it is difficult to see. I hope thats ok?

---
Tags: phylogenetics, phylogeny
---


# Question
Title: How to calculate the number of peaks that are upstream/downstream of some other peaks
I have 3 histone marks,I have used Macs2 for peak calling and diffBind to analyze the peaks. I was wondering if you know any way to calculate the peak numbers of one specific histone mark that occur for example within 2kb of peaks belong to another histone mark.
Much appreciated

# Answer
You can use windowBed from bedtools with `-c` if you just want counts. Or use `-wa` `-wb` and then use groupBy for any other operations.
> 1 votes

---
Tags: sequence-analysis, chip-seq, bedtools, peak-calling, deeptools
---


# Question
Title: How identifiable are human omics data and how to mitigate their identifying features?
Say a database were to store human omics datasets. The human subjects are known and the sample size is rather small in size initially (n=500). The database contains genomics, transcriptomics, proteomics, gut metabolomics, and epigenomics. Since the sample size is rather small initially, it would be important that certain identifying features are mitigated so that individual subjects cannot be identified. This leads me to two questions:

1. What identifiable features could be present in these different types of raw data? (Such as race, sex, age, hair color, eye color, height, where someone has lived, and anything else I may not even be considering)? Which of these omics types would be the most and least dangerous for identifying individuals? How can the identifiable nature of these data be mitigated?
2. If environmental metagenomics is also collected by these human subjects, is it possible to identify the human subjects by contamination? (i.e. some of the reads from the metagenomics data inadvertently contain human reads?) How can the identifiable nature of these data be mitigated?

I think this subject may be a bit futuristic, but I am very uninformed. If there are any references that provide additional thinking about these topics, please kindly share. Thank you for sharing your thoughts.

# Answer
This is one of the major problems with genomic information in todays research. This was highlighted some years ago with the police using publicly available genomic data bases to idenitfy unkown murder suspects.

The full extent of this issue was exemplified in a paper by Yaniv Erlich from 2018 that should be a good starting point for you.
They claim that giving the current (2018) amount genetic information more than half the searches for a person of european descent will result in a thrid cousin or better match. And these odds will only inprove with time.

Therefore, the genetic information not only can give you a hint at inherited variables like the approximate height of a person but can even yield their name and address.

About your second question, I would think that the amount of host (i.e human) material collected alongside the microbiome sample is rather small. This probably is not sufficient to obtain enough coverage for a profile. But I could be mistaken.
> 2 votes

# Answer
PPK provides a great answer, but for question 2 I can provide a different perspective. For shotgun metagenome sequencing (without any enrichment/depletion protocols) it is common for \>90% of reads to map to the human genome. There is variation across body sites, for example the gut microbiome is rich and so the % of human reads from a stool sample will be less compared to say a blood sample.

Checkout this study, where they mined human microbiome datasets for host (human) reads and were able to re-construct draft full human genomes with sometimes 20X coverage. See also part of their discussion:

> We show here that it is possible to reconstruct complete host genomes using metagenomic sequence data, which is potentially identifiable. However, this was possible due to the unique study design of the HMP, whereby multiple body sites from each individual were sequenced at a high depth, allowing us to pool data across body sites and reach a 10x mean coverage per host genome. Common metagenomic shotgun sequencing studies, which usually include an order of magnitude less sequence data, are unlikely to enable such an analysis. Moreover, the majority of studies sequence stool samples, which include many fewer host-derived reads. Nevertheless, we anticipate that future shotgun metagenomics sequencing studies would consider these potential privacy concerns.
> 1 votes

---
Tags: genomics, metagenome, proteomics
---


# Question
Title: How to run MaxQuant in command line mode?
MaxQuant is a software package for mass spectroscopy and proteomics. There is a windows version and a linux version. To run on linux you have to use a program that is called `mono`. I think, it is developed by Microsoft, which I find quite nice from them. I installed mono on Linux and started MaxQuant.

```
mono MaxQuant.exe

```

The documentation claims it is for running MaxQuant in commandline mode, though a GUI is started anyway. How can I start the program in commandline mode to perform batch analysis? I think, there is not enough information given on the web-page.

# Answer
> 2 votes
Give

```
mono MaxQuantCmd.exe [OPTIONS]

```

a try.

University of Oslo description of how to use the tool in their infrastructure.

# Answer
> 0 votes
You can now install MaxQuant simply with Conda:

```
conda install -c bioconda maxquant

```

and to run MaxQuant type:

```
maxquant mqpar.xml

```

Reference

---
Tags: mass-spectrometry
---


# Question
Title: How to extract metadata from NCBI's experiment?
I want to extract metadata from experiment SRX1596422. I have sratoolkit and Entrex direct but not sure what line of code to use. Here's the link for the experiment: experiment SRX1596422

# Answer
> 3 votes
You'll want something like:

```
esearch -db sra -query SRX1596422 | efetch -format runinfo

```

This will produce a CSV output to the screen with columns containing the meta information available in SRA:

```
Run,ReleaseDate,LoadDate,spots,bases,spots_with_mates,avgLength,size_MB,AssemblyName,download_path,Experiment,LibraryName,LibraryStrategy,LibrarySelection,LibrarySource,LibraryLayout,InsertSize,InsertDev,Platform,Model,SRAStudy,BioProject,Study_Pubmed_id,ProjectID,Sample,BioSample,SampleType,TaxID,ScientificName,SampleName,g1k_pop_code,source,g1k_analysis_group,Subject_ID,Sex,Disease,Tumor,Affection_Status,Analyte_Type,Histological_Type,Body_Site,CenterName,Submission,dbgap_study_accession,Consent,RunHash,ReadHash
SRR3184375,2016-02-24 00:00:00,2016-02-23 23:14:45,556018223,111203644600,556018223,200,57997,GCA_000001635.1,https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos2/sra-pub-run-7/SRR3184375/SRR3184375.1,SRX1596422,,WGS,size fractionation,GENOMIC,PAIRED,0,0,ILLUMINA,Illumina HiSeq 2000,SRP070623,PRJNA312530,,312530,SRS1306521,SAMN04492320,simple,10090,Mus musculus,MCNT_E1,,,,,female,,no,,,,,THE SCRIPPS RESEARCH INSTITUTE,SRA356784,,public,52A6F836C13F2D15755EDE3E52A18D78,240E7130CFE36405818DA5EE2A10E24C

```

The more relevant fields are: Subject\_ID, Sex, Disease, Tumor, Affection\_Status, Analyte\_Type, Histological\_Type, Body\_Site. Note that not all of these are available for that particular sample.

---
Tags: phylogenetics, sratoolkit, entrez
---


# Question
Title: Separating peaks of chip-seq with specific length
My data contain several chip-seq results. I have the peaks called by MACS2.I wanted to only look at those peaks that their size is e.g 500bp to 1000bp. How can I separate those peaks efficiently?

I know one way is to take the peaks.xls file sort them based on the length column and take out the peaks less than 1kb. But I was looking for a more efficient way to do that. Any suggestion?

Many Thanks

# Answer
You can do it by awk.

```
awk -v OFS="\t" '{ if ( ($3-$2) >= 500 && ($3-$2) <= 1000) { print } }' in.bed > out.bed

```

Assuming there is no header.
> 1 votes

# Answer
You can use `rtracklayer` from `R bioconductor` to import it, filter and do much more stuff, i use the code from this blog for importing the narrow bed file. Below I use an example file since you did not provide yours:

```
library(rtracklayer)
fl = "zas1_default_peaks.narrowPeak"

extraCols_narrowPeak <- c(signalValue = "numeric", pValue = "numeric",
                          qValue = "numeric", peak = "integer")
gr <- import(fl, format = "BED",
                        extraCols = extraCols_narrowPeak)

head(gr)

GRanges object with 6 ranges and 6 metadata columns:
      seqnames        ranges strand |                name     score signalValue
         <Rle>     <IRanges>  <Rle> |         <character> <numeric>   <numeric>
  [1] AB325691     8982-9219      * | zas1_default_peak_1        18     1.31697
  [2] AB325691   16685-18050      * | zas1_default_peak_2        59     1.51283
  [3]        I 101971-102188      * | zas1_default_peak_3        17     1.23717
  [4]        I 102390-103538      * | zas1_default_peak_4        19     1.24967
  [5]        I 152286-152496      * | zas1_default_peak_5        23     1.37631
  [6]        I 175025-175577      * | zas1_default_peak_6        69     1.62063
         pValue    qValue      peak
      <numeric> <numeric> <integer>
  [1]   3.37358   1.89101       146
  [2]   8.12841   5.92232       419
  [3]   3.26012   1.79545        64
  [4]   3.47875   1.97673       708
  [5]   3.89921   2.32948       136
  [6]   9.34429    6.9608       374

```

The narrowBed and peaks.xls contains the same information. If you really really need to read in the xls, you can do:

```
library(GenomicRanges)
peaks = read.table("zas1_default_peaks.xls",comment="#",header=TRUE)
gr = makeGRangesFromDataFrame(peaks,keep.extra.columns=TRUE)

```

Once you have the genomic ranges object, filtering is the same:

```
 gr[width(gr)>500 & width(gr)<1000]
GRanges object with 208 ranges and 7 metadata columns:
        seqnames          ranges strand |    length abs_summit    pileup
           <Rle>       <IRanges>  <Rle> | <integer>  <integer> <numeric>
    [1]        I   175025-175577      * |       553     175399    186.99
    [2]        I   186807-187489      * |       683     187009    208.66
    [3]        I   339070-339700      * |       631     339258    222.91
    [4]        I   431557-432088      * |       532     431734    195.66

```
> 1 votes

---
Tags: sequence-analysis, differential-expression, chip-seq, peak-calling, macs2
---


title: Extract sequences from partial Header

text:
# Question
Title: Extract sequences from partial Header
I have a library I would like to split them based on their headers. Some are related to RNA and DNA.

The header contains much information, but the most important the presence of DNA or RNA/LTR.., these partial words could be in between or at the beginning.

**The point is the knowledge of how to extract sequences from a partial header occur in between the ID**

My question is how to use grep or awk to grep the header that has one of these words along with sequences? Note that The sequences are more than one line.

Or perhaps grep the specific word and ignore what before and after?

```
>Tigger16a#DNA/TcMar-Tigger DF0000028 TcMar-Tigger **DNA** transposon
>rnd-4_family-38#SINE/MIR ( Recon Family Size = 20, Final Multiple Alignmen
>rnd-6_family-31751#LTR/Gypsy ( Recon Family Size = 26, Final Multiple Alignment Size = 22 )
>RNA2558#LTR/ERVL
>NonDNA1#LINE/I-Jockey
>DNA5#DNA/TcMar-Tc1

```

I have tried to create a list of the Required IDs using grep creating the list and extracting the sequences but for some reasons, the output has more sequences than those specified in the DNAID.txt list.

```
grep -A1000 -w -f DNAID.txt.fa MyLibrary > DNA_Sequence.fa

```

# Answer
You could use Biopython. If its a fasta file, it will be complicated to write fasta output (especially multiline fasta) with `grep` or `awk`. Simple solution is to use biopython, so that you can even match for any complex patterns in the fasta header.

```
from Bio import SeqIO

rna_records  = []
dna_records = []

for seq in SeqIO.parse("in.fa","fasta"):

    if "RNA" in seq.id:
        rna_records.append(seq)

    elif "DNA" in seq.id:
        dna_records.append(seq)

SeqIO.write(rna_records, "RNA_out.fa","fasta")
SeqIO.write(dna_records, "DNA_out.fa","fasta")

```

**This code needs to be optimised if you are dealing with large sequences like human genome fasta as it holds everything in memory.**
> 2 votes

# Answer
Here is a solution with `awk`:

This is the example fasta file which I used:

```
>Test DNA 1
ACGTAAGGATATAATAC
ACGTA
AGGATA
>Test 2 RNA
ACGGT
GGTA
CGGA
>DNA - TEST 3
ACG
GATA
AGGT

```

You can create a file called `fasta_search.awk` and insert following code:

```
#!/usr/bin/awk -f
BEGIN {
    hit=0
}
{
    if ($0 ~ /^>/)
    {
        if ($0 ~ search)
        {
            hit=1;
            print $0;
        }
        else
        {
            hit=0;
        }
     }
    else if(hit==1)
    {
        print $0
    }
}

```

Now you can run: `àwk -v search="DNA" -f fasta_search.awk my_fastas.fa`

`search` is therefore a awk variable which could be any String you are looking for.

As result I get this:

```
>Test DNA 1
ACGTAAGGATATAATAC
ACGTA
AGGATA
>DNA - TEST 3
ACG
GATA
AGGT

```

To get all your keywords, you could run:

```
for i in DNA RNA LTR; do awk -v search="$i" -f fasta_search.awk my_fastas.fa > lib_$i.fasta; done

```

This results in `head lib_*`:

```
==> lib_DNA.fasta <==
>Test DNA 1
ACGTAAGGATATAATAC
ACGTA
AGGATA
>DNA - TEST 3
ACG
GATA
AGGT

==> lib_LTR.fasta <==

==> lib_RNA.fasta <==
>Test 2 RNA
ACGGT
GGTA
CGGA

```
> 1 votes

---
Tags: awk, filtering, grep
---


# Question
Title: Subsetting from seurat object based on orig.ident?
I am pretty new to Seurat. I want to subset from my original seurat object (BC3) meta.data based on orig.ident. however, when i use subset(), it returns with Error.

```
ER_HER_P <- subset(BC3, idents = c("BC03"))
Error: No cells found

table(Idents(BC3))
BC01   BC02   BC03 BC03LN   BC04   BC05   BC06   BC07 BC07LN   BC08   BC09   BC10   BC11
21     50     33     51     42     74     15     48     52     21     55     13     11

```

There are 33 cells under the identity. I checked the active.ident to make sure the identity has not shifted to any other column, but still I am getting the error?

```
head(BC3@active.ident)
[1] BC01 BC01 BC01 BC01 BC01 BC01
Levels: BC01 BC02 BC03 BC03LN BC04 BC05 BC06 BC07 BC07LN BC08 BC09 BC10 BC11

```

# Answer
> 1 votes
If you are going to use idents like that, make sure that you have told the software what your default ident category is.

This works for me, with the metadata column being called "group", and "endo" being one possible group there.

```
Idents(combined.all) <- "group"
endo_subset <- subset(combined.all, idents = c("endo"))

```

I think this is basically what you did, but I think this looks a little nicer

# Answer
> 0 votes
A bit too long to type as a comment.

To use subset on a Seurat object, (see `?subset.Seurat`) , you have to provide:

```
Subset a Seurat object

Description:

     Subset a Seurat object

Usage:

     ## S3 method for class 'Seurat'
     x[i, j, ...]

     ## S3 method for class 'Seurat'
     subset(x, subset, cells = NULL, features = NULL, idents = NULL, ...)

Arguments:

       x: Seurat object to be subsetted

i, features: A vector of features to keep

j, cells: A vector of cells to keep

```

So you either use the matrix to subset:

```
library(Seurat)
data(pbmc_small)
Idents(pbmc_small) = paste0("BC",Idents(pbmc_small))
table(Idents(pbmc_small))

BC0 BC2 BC1
 36  19  25

test = pbmc_small[,Idents(pbmc_small)=="BC0"]
table(Idents(test))

BC0
 36

```

Or you provide the cells:

```
subset(pbmc_small,cells=colnames(pbmc_small)[Idents(pbmc_small)=="BC0"])
An object of class Seurat
230 features across 36 samples within 1 assay
Active assay: RNA (230 features, 20 variable features)
 2 dimensional reductions calculated: pca, tsne

```

What you have should work, but try calling the actual function (in case there are packages that clash):

```
Seurat:::subset.Seurat(pbmc_small,idents="BC0")

An object of class Seurat
230 features across 36 samples within 1 assay
Active assay: RNA (230 features, 20 variable features)
 2 dimensional reductions calculated: pca, tsne

```

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


# Question
Title: maftools tcgaCompare plot: x-axis labels getting clipped
I have a *{maftools}* plot:

whatever I try to save that completely I see the down part of plot is trimmed like

I have tried both save as pdf and save as image options which did not differ

Any suggestion please?

Plot comes by this code

```
tcgaCompare(maf = laml,cohortName = "Cohort",rm_zero = FALSE,decreasing = TRUE,logscale = TRUE,primarySite = TRUE,col = c("grey30", "red"),medianCol = "black")

```

# Answer
One way is to increase the outer margins of plot.

```
par(oma=c(7,0,0,0))
tcgaCompare(maf = laml, primarySite = TRUE)

```
> 2 votes

---
Tags: r, maftools
---


title: Comparing aligned amino acids to codon

text:
# Question
Title: Comparing aligned amino acids to codon
I have a set of 4 amino acids which i have aligned,but i want to compare then with respective nucleotide sequence in terms of change codon level.

I have the alignment file where i want to take only for the C-terminal domain region region only the CTD starts with amino acids `“NITNLC”` till it ends with `“HAPATV”` and rest of the sequence i dont want to take.

In short i want to take the part sequence which starts from "NITNLC" to "HAPATV" then i would like to compare what are the changes in the nucleotide sequence at codon level.

1. To start with I have used the seqinr library to generate nucleotide alignment from aligned protein sequence.

reverse.align {seqinr}
2. Then the next step that i want to compare each of my amino acid sequence to my reference sequence which in my case is `P0DTC2`to rest of the sequence for example I want to compare `P0DTC2` to `K9N5Q8`to the above mentioned CTD region which starts from `NITNLC` and ends with `HAPATV`and find what are the codon level changes if there is any I would like to report both in amino acid level and codon level.

First part I guess my approach is correct. Second part I'm not sure how to proceed i guess it's much more than simple parsing I guess!!.

Any help or suggestion would be really appreciated ,and if there exist a R based solution that would be really welcome.

Files Amino acid file nucleotide file Aligned amino acid file Reverse alignment nucleotide file using seqinr

# Answer
> 2 votes
To set the ball rolling, I tried to implement it manually with your files. Overall, there are three steps:

1. Locate `NITNLC` (or `HAPATV`) in protein `P0DTC2`.
2. Locate mismatched amino acids between proteins `P0DTC2` and `K9N5Q8` within the range from step 1.
3. Print the amino acid and the DNA codon of both proteins from step 2.

It works, but only for the first 60 amino acids. I wonder if it's due to the amino sequences in `AMINOOO_seq_removed.fasta` repeat itself every 60 acids. But why?

```
#For example, the first three lines of protein P0DTC2
>P0DTC2
MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDLFLPFFS MFVFLVLLPLVSSQCVNLTTRTQLPPAYTNSFTRGVYYPDKVFRSSVLHSTQDLFLPFFS
NVTWFHAIHVSGTNGTKRFDNPVLPFNDGVYFASTEKSNIIRGWIFGTTLDSKTQSLLIV NVTWFHAIHVSGTNGTKRFDNPVLPFNDGVYFASTEKSNIIRGWIFGTTLDSKTQSLLIV
NNATNVVIKVCEFQFCNDPFLGVYYHKNNKSWMESEFRVYSSANNCTFEYVSQPFLMDLE NNATNVVIKVCEFQFCNDPFLGVYYHKNNKSWMESEFRVYSSANNCTFEYVSQPFLMDLE
...

```

Step 0: Read the files.

```
library('seqinr')

#align file containing protein sequences
count_added <- read.alignment('count_added_.clustal_num', format='clustal')
names(count_added$seq) <- count_added$nam

#DNA sequences
rev3.aln <- read.alignment('rev3.aln', format='fasta')
names(rev3.aln$seq) <- rev3.aln$nam

```

Step 1: Locate `NITNLC` (or `HAPATV`) in protein `P0DTC2`. Again, there are two repeats of `NITNLC` 60 acids apart (at 807 and 870).

```
locate <- function(seq, find)
  {address <- gregexpr(paste(strsplit(find, '')[[1]], collapse='[^a-z]*'), seq)
  #substr(seq, address[[1]][1], address[[1]][1]+attr(address[[1]], 'match.length')[1]-1)

  return(list(start=as.numeric(address[[1]]),
             end=as.numeric(address[[1]] + attr(address[[1]], 'match.length') - 1)))
  }

locate(seq=count_added$seq[['P0DTC2']], find='nitnlc')
#start
#807 870
#end
#812 875
locate(seq=count_added$seq[['P0DTC2']], find='hapatv')
#start
#1212 1279
#end
#1217 1285

```

Step 2: Locate mismatched amino acids between proteins `P0DTC2` and `K9N5Q8` within the range from step 1. Instead of using the range 812 - 1279, I chose range 1 - 20 for demonstration purposes.

```
compare <- function(seq1, seq2, after=0, before=100000)
  {seq1_ = strsplit(seq1, '')[[1]]
  seq2_ = strsplit(seq2, '')[[1]]
  ind = which(seq1_ != seq2_ & grepl('[a-z]',seq1_) & grepl('[a-z]',seq2_))
  ind = ind[ind>after & ind<before]
  #seq1_[ind[1]]
  #seq2_[ind[1]]

  return(ind)
  }
compare(seq1=count_added$seq[['P0DTC2']], seq2=count_added$seq[['K9N5Q8']], after=1, before=20)
# [1] 5  7  8  9 10 13 14

#protein comparison
#K9N5Q8      MIHSVFLLMFLLTPTESYVD
#P0DTC2      ----MFVFLVLLPL------
#<mismatch>      5 7890  34

```

Step 3: Print the amino acid and the DNA codon of both proteins from step 2. Note that `ind` is based on the align file.

```
print_amino_codon <- function(ind, seq, seq_gene)
  {locate_amino <- gregexpr('[a-z]', seq)[[1]]
  if (!ind %in% locate_amino) return(NA)
  ind2 = match(ind, locate_amino)

  return(c(amino=substr(seq, ind, ind), codon=substr(seq_gene, ind2*3-2, ind2*3)))
  }

codon(ind=5, seq=count_added$seq[['K9N5Q8']], seq_gene=rev3.aln$seq[['K9N5Q8']])
#amino codon
#  "v" "atg"
codon(ind=6, vseq=count_added$seq[['K9N5Q8']], seq_gene=rev3.aln$seq[['K9N5Q8']])
#amino codon
#  "f" "ttc"
codon(ind=7, seq=count_added$seq[['K9N5Q8']], seq_gene=rev3.aln$seq[['K9N5Q8']])
#amino codon
#  "l" "ttg"

codon(ind=5, seq=count_added$seq[['P0DTC2']], seq_gene=rev3.aln$seq[['P0DTC2']])
#amino codon
#  "m" "atg"
codon(ind=6, seq=count_added$seq[['P0DTC2']], seq_gene=rev3.aln$seq[['P0DTC2']])
#amino codon
#  "f" "ttt"
codon(ind=7, seq=count_added$seq[['P0DTC2']], seq_gene=rev3.aln$seq[['P0DTC2']])
#amino codon
#  "v" "gtt"

#protein comparison
#K9N5Q8      MIHSVFLLMFLLTPTESYVD
#P0DTC2      ----MFVFLVLLPL------
#<mismatch>      5 7890  34
#<print>         ^^^

#K9N5Q8 gene codon
#gtg ttt cta ctg atg ttc ttg tta aca
#                ^^^ ^^^ ^^^
#P0DTC2 gene codon
#atg ttt gtt ttt ctt
#^^^ ^^^ ^^^

```

---
Tags: r, sequence-alignment
---


# Question
Title: parsimony and maximum likelihood tree comparison in R
Maximum Likelihood analysis compares the tips or the species based on their sequence similarity while Parsimony analysis compares the characteristic features among the species.

So far I did comparison of both the tree using dendextend. Is there any other methods or library to make the comparison more quantitative.Due to my limited understanding I have this in mind. Can i even compare these two methods?

# Answer
I don't use R to build trees, I use established subprocessing stripts in Python and ETE3. So I've not run this code, but I do understand the critical parts of the code.

The hard bit is bootstrapping a maximum likelihood tree, first is the ML model 1.

```
fitGTR <-optim.pml(fitGTR, model="GTR", optInv=TRUE, optGamma=TRUE, + rearrangement = "NNI", control = pml.control(trace = 0))

```

The above model is perfect, general time reversible model with a gamma correction and an estimate of the invariant sites. This is the most important part of tree estimation. The NNI for tree rearrangements is a good one. So THIS is the model you want. You can emperical test this, but these days we don't bother.

The information given for bootstrapping,

2.

```
BStree <- bootstrap.pml(x, bs = 500, trees = TRUE, multicore = FALSE,
  mc.cores = NULL, ...)

```

The critical parameter is bs=500. By default it uses bs=100 which I don't think is enough. I would set multicore = TRUE to speed the calculation (it will speed it up alot) and I assume there must be a setting for the number of cores (?).

I assume you can then stack all the parameters from 1/ into code 2. So in other words the model you used in part 1. (GTR PINVAR and Gammma) is same model used in the bootstrapping.

This will produce a bootstrap file, probably of 500 bootstrapped ML trees, so you need to produce a consensus and represent the values on the optimal ML tree.

Thus you need to overlay the bootstrap onto the tree that is output in 1. This says fitGTR and I suspect that will be a tree file with a likelihood, you might need to set a tree output and that would be 'tree' in the code below. 3.

```
plotBS(tree, BStrees, type = "myroot", bs.col = "black",
  bs.adj = NULL, p = 50, frame = "none", ...)

```

The type = myroot is your defined outgroup for the tree. This should be your nicely bootstrapped tree.

There was also this code in the help manual and I'm not sure what this is doing. I suspect it is just producing bootstrapped datasets, so its not needed here.

```
bootstrap.phyDat(x, FUN, bs = 100, multicore = FALSE,
      mc.cores = NULL, jumble = TRUE, ...)

```

To do the same with parsimony then just replace pml for parsimony. Comparing 2 bootstrap trees is a node is robust when it scores \>80%
> 1 votes

# Answer
> Can i even compare these two methods?

See if this paper helps: https://www.nature.com/articles/nature02917
> 2 votes

---
Tags: r, phylogenetics
---