# Question
Title: repophlan script to download bacterial genome
I am trying to donwload the microbial genome using the repophlan\_get\_microbes.py (https://github.com/SegataLab/repophlan) but now it is running more than 10 days and still not finished on slurm HPC.

Could anyone please suggest what is the issue with the script?

```
#!/bin/bash
#SBATCH --job-name=download.GENOME    # Job name
#SBATCH --ntasks=60
#SBATCH --nodes=3                  # Number of nodes
#SBATCH --time=2000:00:00              # Time limit hrs:min:sec
#SBATCH --partition=k2-himem
#SBATCH --mem=800G

module add apps/python3/3.7.4/gcc-4.8.5
module add libs/biopython/1.73/gcc-4.8.5+python-2.7.8+numpy-1.16.0

python repophlan_get_microbes.py --taxonomy taxonomy_reduced.txt --out_dir microbes --nproc 20 --out_summary repophlan_microbes.txt | tee repophlan_microbes.log

$ head repophlan_microbes.log
2020-09-17 20:44:59,964 ./repophlan_get_microbes.py INFO     Reading the taxonomy from taxonomy_reduced.txt...
2020-09-17 20:45:04,921 ./repophlan_get_microbes.py INFO     Done.
2020-09-17 20:45:10,716 ./repophlan_get_microbes.py WARNING  GCF_000001215.4 [Drosophila melanogaster ] excluded from download because of uninteresting phyla!
2020-09-17 20:45:10,716 ./repophlan_get_microbes.py WARNING  GCF_000001405.39 [Homo sapiens ] excluded from download because of uninteresting phyla!
2020-09-17 20:45:10,716 ./repophlan_get_microbes.py WARNING  GCF_000001635.26 [Mus musculus ] excluded from download because of uninteresting phyla!
2020-09-17 20:45:10,716 ./repophlan_get_microbes.py WARNING  GCF_000001735.4 [Arabidopsis thaliana ecotype=Columbia] excluded from download because of uninteresting phyla!
2020-09-17 20:45:10,716 ./repophlan_get_microbes.py WARNING  GCF_000001895.5 [Rattus norvegicus strain=mixed] excluded from download because of uninteresting phyla!
2020-09-17 20:45:10,716 ./repophlan_get_microbes.py WARNING  GCF_000001905.1 [Loxodonta africana ] excluded from download because of uninteresting phyla!
2020-09-17 20:45:10,716 ./repophlan_get_microbes.py WARNING  GCF_000002035.6 [Danio rerio ] excluded from download because of uninteresting phyla!
2020-09-17 20:45:10,716 ./repophlan_get_microbes.py WARNING  GCF_000002075.1 [Aplysia californica ] excluded from download because of uninteresting phyla!

$ tail repophlan_microbes.log

2020-09-21 20:42:26,871 repophlan_get_microbes.py INFO     Parsing of ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/001/635/075/GCA_001635075.1_ASM163507v1/GCA_001635075.1_ASM163507v1_genomic.gbff.gz
2020-09-21 20:42:27,198 repophlan_get_microbes.py INFO     Parsing of ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/010/319/175/GCA_010319175.1_PDT000567389.1/GCA_010319175.1_PDT000567389.1_genomic.gbff.gz
2020-09-21 20:42:27,482 repophlan_get_microbes.py INFO     Parsing of ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/003/518/025/GCA_003518025.1_ASM351802v1/GCA_003518025.1_ASM351802v1_protein.faa.gz
2020-09-21 20:42:27,608 repophlan_get_microbes.py INFO     Parsing of ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/001/788/535/GCA_001788535.1_ASM178853v1/GCA_001788535.1_ASM178853v1_protein.faa.gz
2020-09-21 20:42:27,628 repophlan_get_microbes.py WARNING  No remote file found (some ffn and faa are know to be missing remotely). Aborting the download of ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/010/997/445/GCA_010997445.1_ASM1099744v1/GCA_010997445.1_ASM1099744v1_protein.faa.gz. <urlopen error ftp error: [Errno ftp error] 550 GCA_010997445.1_ASM1099744v1_protein.faa.gz: No such file or directory>
2020-09-21 20:42:27,628 repophlan_get_microbes.py ERROR    Error in downloading of ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/010/997/445/GCA_010997445.1_ASM1099744v1/GCA_010997445.1_ASM1099744v1_protein.faa.gz <urlopen error ftp error: [Errno ftp error] 550 GCA_010997445.1_ASM1099744v1_protein.faa.gz: No such file or directory>
2020-09-21 20:42:27,628 repophlan_get_microbes.py ERROR    Error in downloading of ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/010/997/445/GCA_010997445.1_ASM1099744v1/GCA_010997445.1_ASM1099744v1_protein.faa.gz: empty fna file!
2020-09-21 20:42:28,239 repophlan_get_microbes.py WARNING  No remote file found (some ffn and faa are know to be missing remotely). Aborting the download of ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/902/778/895/GCA_902778895.1_Rumen_uncultured_genome_RUG12414/GCA_902778895.1_Rumen_uncultured_genome_RUG12414_protein.faa.gz. <urlopen error ftp error: [Errno ftp error] 550 GCA_902778895.1_Rumen_uncultured_genome_RUG12414_protein.faa.gz: No such file or directory>
2020-09-21 20:42:28,239 repophlan_get_microbes.py ERROR    Error in downloading of ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/902/778/895/GCA_902778895.1_Rumen_uncultured_genome_RUG12414/GCA_902778895.1_Rumen_uncultured_genome_RUG12414_protein.faa.gz <urlopen error ftp error: [Errno ftp error] 550 GCA_902778895.1_Rumen_uncultured_genome_RUG12414_protein.faa.gz: No such file or directory>
2020-09-21 20:42:28,239 repophlan_get_microbes.py ERROR    Error in downloading of ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/902/778/895/GCA_902778895.1_Rumen_uncultured_genome_RUG12414/GCA_902778895.1_Rumen_uncultured_genome_RUG12414_protein.faa.gz: empty fna file!

```

Many thanks

# Answer
> 1 votes
**UPDATE 2: new log posted**

**take-home message: I think it is working, it's just that a few errors happen sometimes when downloading things from big databases.**

It looks like there are a few warnings or errors associated with getting *specific* files. I think this is actually ok. For example for one of the genomes, I logged into the ftp manually to check if the file was there:

```
$ ftp ftp://ftp.ncbi.nlm.nih.gov
Trying 130.14.250.13...
Connected to ftp.wip.ncbi.nlm.nih.gov.
220-
 This warning banner provides privacy and security notices consistent with
 applicable federal laws, directives, and other federal guidance for accessing
 this Government system, which includes all devices/storage media attached to
 this system. This system is provided for Government-authorized use only.
 Unauthorized or improper use of this system is prohibited and may result in
 disciplinary action and/or civil and criminal penalties. At any time, and for
 any lawful Government purpose, the government may monitor, record, and audit
 your system usage and/or intercept, search and seize any communication or data
 transiting or stored on this system. Therefore, you have no reasonable
 expectation of privacy. Any communication or data transiting or stored on this
 system may be disclosed or used for any lawful Government purpose.
220 FTP Server ready.
331 Anonymous login ok, send your complete email address as your password
230 Anonymous access granted, restrictions apply
Remote system type is UNIX.
Using binary mode to transfer files.
200 Type set to I
ftp> get genomes/all/GCA/902/778/895/GCA_902778895.1_Rumen_uncultured_genome_RUG12414/GCA_902778895.1_Rumen_uncultured_genome_RUG12414_protein.faa.gz
local: genomes/all/GCA/902/778/895/GCA_902778895.1_Rumen_uncultured_genome_RUG12414/GCA_902778895.1_Rumen_uncultured_genome_RUG12414_protein.faa.gz remote: genomes/all/GCA/902/778/895/GCA_902778895.1_Rumen_uncultured_genome_RUG12414/GCA_902778895.1_Rumen_uncultured_genome_RUG12414_protein.faa.gz
ftp: Can't access `genomes/all/GCA/902/778/895/GCA_902778895.1_Rumen_uncultured_genome_RUG12414/GCA_902778895.1_Rumen_uncultured_genome_RUG12414_protein.faa.gz': No such file or directory

ftp> ls genomes/all/GCA/902/778/895/GCA_902778895.1_Rumen_uncultured_genome_RUG12414/

229 Entering Extended Passive Mode (|||50025|)
150 Opening ASCII mode data connection for file list
-r--r--r--   1 ftp      anonymous     2725 Feb 23  2020 GCA_902778895.1_Rumen_uncultured_genome_RUG12414_assembly_report.txt
-r--r--r--   1 ftp      anonymous     3456 Feb 23  2020 GCA_902778895.1_Rumen_uncultured_genome_RUG12414_assembly_stats.txt
dr-xr-xr-x   3 ftp      anonymous     4096 Feb 23  2020 GCA_902778895.1_Rumen_uncultured_genome_RUG12414_assembly_structure
-r--r--r--   1 ftp      anonymous      149 Apr 21 19:29 GCA_902778895.1_Rumen_uncultured_genome_RUG12414_feature_count.txt.gz
-r--r--r--   1 ftp      anonymous   949241 Feb 23  2020 GCA_902778895.1_Rumen_uncultured_genome_RUG12414_genomic.fna.gz
-r--r--r--   1 ftp      anonymous  1353865 Feb 23  2020 GCA_902778895.1_Rumen_uncultured_genome_RUG12414_genomic.gbff.gz
-r--r--r--   1 ftp      anonymous      150 Feb 23  2020 GCA_902778895.1_Rumen_uncultured_genome_RUG12414_genomic_gaps.txt.gz
-r--r--r--   1 ftp      anonymous      840 Feb 23  2020 GCA_902778895.1_Rumen_uncultured_genome_RUG12414_wgsmaster.gbff.gz
-r--r--r--   1 ftp      anonymous    44529 Sep  2 20:26 README.txt
-r--r--r--   1 ftp      anonymous      410 Feb 23  2020 annotation_hashes.txt
-r--r--r--   1 ftp      anonymous       14 Sep 21 07:11 assembly_status.txt
-r--r--r--   1 ftp      anonymous     1399 May 20 23:00 md5checksums.txt
226 Transfer complete

```

So in this particular case, that genome *exists* but is just missing a `faa.gz` file. So **the script appears to be working as expected**, it's just that some of the NCBI submissions don't have everything they are supposed to have.

(As a tangent, that Rumen genome came from a paper I'm an author on I'm pretty sure!! So I guess you could blame me for the missing file if you wanted, though I was in no way involved with annotations or submission...)

**UPDATE: now the log is posted**

It does not look like a connectivity issue based on the new information provided in the first and last lines of the log file. There it looks like there are a lot of processes finishing quickly towards the beginning and the end of the process (e.g. 48 hours apart).

Based on this new information, I am wondering if you have simply managed to get a taxonomy file that is much larger (e.g. \>10X larger to explain the time difference) than anticipated when the README for the script was written, and thus it simply takes a long time.

It would still probably be best to go through the log file rather carefully and see if there is a big hang in the process, errors, warnings, etc.

---

Here is what they say in the repo:

```
*** Download the microbial (bacteria+archaea) genomes (fna,ffn,faa,frn) with
    corresponding taxonomy
$ ./repophlan_get_microbes.py --taxonomy taxonomy_reduced.txt --out_dir /
  microbes --nproc 20 --out_summary repophlan_microbes.txt /
  | tee repophlan_microbes.log
[ ~6 hours ]. Downloads and save all available sequences for microbes. It uses
multiple parallel connections to speed up the retrieval: common exceptions
including issues with the NCBI ftp which temporarily rejects connections for
ftp load problems are handled gracefully with a 'delay and retry' policy.
IMPORTANT: specifying more than 20 processors to run in parallel causes serious
problems in terms of exceeding the allowed number of connections by NCBI causing
long delays.

```

Since your script is taking much much longer than 6 hours, it probably didn't work.

Some thoughts based only on the info in the repo:

* "IMPORTANT: specifying more than 20 processors to run in parallel causes serious problems in terms of exceeding the allowed number of connections by NCBI causing long delays." **this sounds like your problem. I think you are asking for 60 nodes ("tasks"? "nodes" is 3 which should be ok). Drop tasks to 20 or less and see what happens, maybe.**
* "...temporarily rejects connections for ftp load problems are handled gracefully with a 'delay and retry' policy." **It is possible that there is a bad request in there that it keeps trying over and over again unsuccessfully, it could keep going forever if such a problem exists and you don't stop it.**

As pointed out in the comments, you should check the log and post that here.

As also pointed out in the comments, you might consider posting an issue at the repo as they are more familiar than we are.

Other things to check:

* how much disk space do your nodes have access to? is it enough for the download?
* are you satisfying all dependencies of the script? In other words do you need to load more modules?

---
Tags: python, bash, microbial-genomics, hpc
---


# Question
Title: What programs account for structural alignment of different parts of distant homologs which have significant structural differences?
If there is a need to perform structural alignment of different parts of distant homologs, which program one should use?

Since distant homologs often have significant structural changes, meaning the assumption of a similar physico-chemical environment for aligned residues would be violated.

Any program/tools or suggested reference will be a great help.

Thanks.

# Answer
TM-align has some options to do this, see for example https://zhanglab.ccmb.med.umich.edu/TM-align/help.html.

TM-align is similar to TM-score but the alignment is purely structural.
> 1 votes

# Answer
You can use `progressiveMauve` (http://darlinglab.org/mauve/user-guide/progressivemauve.html), or `bubbz` (https://www.cell.com/iscience/fulltext/S2589-0042%2820%2930409-0)
> 0 votes

---
Tags: sequence-alignment, proteins, protein-structure, 3d-structure, structural-biology
---


# Question
Title: Interpreting GWAS results with different settings
I did a bunch of GWAS analysis (linear model without covariates) with applying different quality controls. How to choose the optimal settings when filtering for minor allele frequency (maf), Hardy-Weinberg equilibrium (HW) and missing call rates (mcr)? Are there any rules? Here are some examples:

Mahnhatten and QQ plot with settings maf=0.02, HW=1e-10, mcr = 0.2:

Manhattan and QQ plot with settings maf=0.01, HW=0.0001, mcr = 0.01:

Mahnhatten and QQ plot with settings maf=0.02, HW=1e-10, mcr = 0.05:

# Answer
The genomes under analysis are under predominant linkage disequilibrium, i.e. linkage. Chromosome 1 shows strong (or stronger - cut-off depending) HW signal as do one or two telomeres - telomeres are known to undergo recombination, so this may not be unsurprising. Nothing wrong with the HW analysis for the core conclusions, its the stuff between the red and blue lines thats difficult (Q-Q plots).

The Q-Q plots are a bit more complex and I would go with the lower graph you presented together with its HW-analysis. There is clearly part of the data you are not accounting for, @gringer says the Q-Q plot conform to the beta distribution .... they should be bang linear. Lets look in more detail,

* 0 to 2 ... predominant disequilibrium that looks fine, nice and linear
* 5.5+ .... HW signal that looks ok, nice and linear
* 2 to 5.5 ... complicated

Given you have a cut-off (red line) at around 5.5 then I'd sort of try and talk your way out of this one. You can do a full investigation, such as transformation to linearise the entire Q-Q plot. Personally I've never done this within HW analysis so I dunno, I've done it for GLM. I'd just say they between the blue and red line represents a complex HW signal which we can't be interpreted but either side of that boundary we're good to go. In other words keep it simple.
> 0 votes

# Answer
Based on the QQ plots, filter out variants with p \> 0.01, then display the beta statistic (or whatever other statistic you're using to demonstrate association). If you want to incorporate the p-value in the displayed results, show the minimum test statistic from a 95% confidence interval. Here's an attempt I made at doing that with some GWAS summary statistics relating to household income (something that wouldn't be expected to have a heritable component):

https://twitter.com/gringene\_bio/status/1207617723586371584

After looking at the test statistic, you may be able to make better decisions about thresholding with the other values. One thing you could try is to incorporate those values into a scatter plot that shows more than two dimensions. The other statistics could be represented as size, colour, or Z location on a 3D scatter plot.

P-values should not be overemphasised, used for ranking, or displayed on their own. See the American Statistical Association's report on the p-value for more information:

https://amstat.tandfonline.com/doi/full/10.1080/00031305.2016.1154108#\_i30

As far as I'm aware, the p-values as calculated by GWAS software represent uncertainty in the test statistic (i.e. t, or beta), and *not* the probability of association. If you want to calculate that association probability, carry out bootstrap sub-sampling of the individuals and then determine the proportion of sub-samples in which the tested variant appears in, say, the highest 5% of variants, ranked by the test statistic. Instead of filtering on HWE, MAF and MCR, you could try filtering on the bootstrap proportion instead, which should incorporate all statistics that introduce uncertainty into the dataset.

More details and figures can be found on this poster:

https://f1000research.com/posters/5-2190

If you're going to be generating a PCA on these data, please consider labelling the cases and controls (or low/high for quantitative traits) on the PCA to verify that the groups look similar in terms of overall population structure. Hopefully you won't find any blobs in the PCA that have only one of the two groups in them, but if you do, they should be dealt with first. Any differences in population structure (even very small, subtle differences) will dominate the GWAS results because structure-associated signal will cover a huge region of the genome in comparison to effects from a few (or a few hundred) SNPs. Here are some papers I've found insightful regarding this issue:
> 2 votes

---
Tags: python, gwas, plink
---


# Question
Title: Advantages of paired-end sequencing compared to single end
One of the advantages of paired end sequencing over single end is that it doubles the amount of data. Another **supposed** advantage is that it leads to more accurate reads because if say Read 1 (see picture below) maps to two different regions of the genome, Read 2 can be used to help determine which one of the two regions makes more sense. However, I don't understand how this second advantage pertains only to paired end sequencing. With single-end sequencing, couldn't you just use the read that comes right after read 1 to essentially get a longer net read?

# Answer
Like you mentioned, having read pairs can help with alignment. Your FASTQ for a single-end sequencing run will look something like this:

```
@SEQ_ID_0001
AGCTAGCGCGGTTGGCTTAGCGACT
+
!''*((((***+))%%%++)(%%%%
@SEQ_ID_0002
AGGTGGTTGTAGGGAAAAAAGTCTC
+
!''*((((***+))%%%++)(%%%%
...

```

There is not necessarily a relationship between `@SEQ_ID_0001` and `@SEQ_ID_0002`. They could come from neigbouring regions in the genome or they could come from different species if that's what your experiment involves. So aligning `@SEQ_ID_0001` doesn't exactly give you information about where `SEQ_ID_0002` goes.

`@SEQ_ID_0001` could equally map to multiple locations in the genome with the same probability. How do we decide where it "should" go? We can randomly distribute multi-mapped reads or do some other fancy stuff, but that usually relies on facts about the data as a whole, and not something specific to the original DNA molecule that sequencing read came from.

If we have paired-end sequencing, then like your diagram shows, you know that read 1 (`R1`) and read 2 (`R2`; these have the same `@SEQ_ID`) will relate to each other in some way. `@SEQ_ID_0001 R1` may be multi-mapped, but if `@SEQ_ID_0001 R2` isn't, then we know that `R1` should probably be relatively close to `R2`. Because of how you prepare your DNA samples before sequencing, you'll have some information about the insert length distribution. You can use this information with sequence aligners to better identify where `R1` should be placed.

You can also do this if both `R1` and `R2` are multi-mapped; certain pairwise alignments will be more likely than others, even if the marginal alignment probabilities for each read are the same.

Finally, there are other unexpected benefits that you can get from paired-end sequencing. All of the above applies to just the mapping step, but applications like structural variant detection and chromosome conformation capture really benefit from paired-end sequencing.

If you have single-end reads for structural variant detection, then you will only detect breakpoints if the read overlaps a breakpoint. But if you have paired-end, you not only have twice the amount of DNA sequenced to look for breakpoints, you can also estimate if there is a breakpoint between the two reads (in the inner distance from your diagram) if the two reads align really far apart from each other, or on separate chromosomes. The idea is similar for chromosome conformation capture and the ligation junction sites.
> 2 votes

# Answer
In paired-end sequencing, you have an expected fragment length. So if Read 1 maps to two different places, whereas Read 2 only maps to one place. Then you can determine which one is more likely for Read 1 to map to. It cannot be too far away from Read 2 given the expected fragment length.

It is not necessarily true that paired-end sequencing gets twice the amount of data. The total amount of data can be determined when running sequencing. If single-end has 20M reads, I can just stop sequencing at 20M as well for paired-end.
> 1 votes

---
Tags: rna-seq
---


# Question
Title: Add segregating sites to branches of a tree
I'm trying to figure out how I would plot a tree with number of segregating sites on display on the branches

I'm using acctran from phangorn to plot a parsimony tree (using phydat object from fasta files)

I have been trying to use seg.sites from the ape package to somehow merge these objects together but there has to be some easy way to plot a tree with the information that it was made from

# Answer
> 1 votes
The *segregating sites* in parsimony are formally called *steps*. Any number of programs will do this but I don't really use R, particularly for phylogeny, You can of course simply look at the raw treefile, all the data is there (described in the last paragraph). If look at the options for *genetic distance*, in parsimony this will be an integer and again the technical name is *steps* and is the *segreating sites* you want.

Python's `ETE3` package will work (but you'd need to pipe it through Biopython), as well as specific standalone packages `Figtree`, the steps or 'genetic distance' on the branches on the branches option is in the check box. The heavy weight program for tracking step associated phylogeny using parsimony is `Mesquite` and its famous grandparent is called `MacClade`.

I suspect its a nomenclature thing, the options probably talk about *genetic distance* because few people do parsimony phylogeny now. Thats the option you want and you will know because it will be \>\>1, whilst distance and model-based methods will be between 0 to 1. From a tree viewing it doesn't matter if the genetic distance is 0 to 1 or 1 and beyond because its all proportional. Its just an option thing you are missing. Thus if there is an option, e.g. to include *bootstrap support* on the tree, there'll be an option for *genetic distance*.

Just make sure you don't confuse the number of steps per branch with a consensus statistic, because there're completely different but might give similar looking outputs. Again its easily checked against the raw treefile because the "steps" will be immediately adajacent to the taxa name for the "twigs" (terminal branches) in the tree.

---
Tags: r, fasta, phylogenetics, linux
---


# Question
Title: predict the IC50 after using Schrodinger suite Glide tool for docking
How to predict the half-maximal inhibitory concentration (IC50) of my ligands.

I'm currently using Schrodinger suite software and I don't know how to predict the IC50?

# Answer
I don't think it is possible just from the docking scores. Docking scores are designed to give you an estimate of the binding affinity (entalphy), but they do not take entropic contributions into account very well. The best you can do are molecular dynamics simulations to calculate the free energy of binding. Which can be computationally quite costly e.g. with GROMACS or CHARMM.

You can estimate the IC50 with molecular dynamics simulations as well, but that is not straightforward. E.g. with Umbrella simulations. An example for a potassium channel and a very simple ligand can be found here.

Another approach would be to build a predictive model. If you have a library of compounds at hand you could use the chemical structures and e.g. the docking scores to estimate the IC50 based on statistics or machine learning. But biology is very messy in general, though if you get a representative dataset with very similar compounds to the one that you have in your dataset then you might be successful. You have to be sure that the binding site and pose are correct though.
> 2 votes

---
Tags: docking
---


# Question
Title: What should be the numbers in the fastSTRUCTURE file?
The structure file for fastSTRUCTURE seem to be a table with rows and columns but the cells are all numbers, while my genotypes are simply bases. How could I know what I should insert there?

# Answer
> 4 votes
Numbers, usually in the form of 0/1/2, is a common way many formats represent bases.

For example, if a SNPs has bases which are A/T, then these are encoded as 0/1, with the reference allele/base taking the 0 and the alternate allele/base taking 1.

However, if you read the documentation, it says:

> Specifically, rows in the data file correspond to samples, with two rows per sample (note that only diploids are handled by this software), and columns correspond to SNPs. The first 6 columns of the file will be ignored; these typically would include IDs, metadata, etc. This software only handles bi-allelic loci. **The two alleles at each locus can be encoded as desired**; however, missing data should be encoded as \`-9'.

So you should be able to use your base letters here and it will work OK as long as the columns and rows are formatted correctly.

---
Tags: file-formats
---


# Question
Title: How to get the uniparc ids for all members of uniref cluster in batch request?
I currently use the 'UniProt website REST API' (example) to collect the UniParc IDs of all the members of the given cluster. Fetching each entry individually is very slow though, so I wanted to ask if there is a way that I can do this request in batches.

**What I have already tried:** There are pages mentioning a batch retrieval of entries. As I understood it, all the examples only adress the ID mapping service. Here I have tried adjusting the given example to the following

```
import urllib.parse
import urllib.request

url = 'https://www.uniprot.org/uploadlists/'

params = {
'from': 'NF100',
'to': 'UPARC',
'format': 'tab',
'query': 'UniRef100_Q8RKD8'
}

data = urllib.parse.urlencode(params)
data = data.encode('utf-8')
req = urllib.request.Request(url, data)
with urllib.request.urlopen(req) as f:
   response = f.read()
print(response.decode('utf-8'))

```

which resulted in the output

```
From    To
UniRef100_Q8RKD8        null

```

I also checked the interface of this service on the website and there I can't seem to map NF100 to UPARC, so I'm assuming that is why I'm not getting a result via the API.

I also know that the UniRef Id contains either a UniParc ID or the UniProtKB accession, so I could use the UniParc ID directly or fetch it for the UniProtKB via the proteins api with the uniparc endpoint, but from what I can see I only get the UniParc ID/entry of one of the members.

# Answer
**Update 17.10.22**

It seems that the api changed and the previous solution doesn't work anymore. I updated my solution by using the Uniprot id mapping service. All the following instructions are copied from the documentation page but added here because link answers are discouraged. The documentation page also includes an implementation in python. The basic example is as follows (using commandline tools)

```
$ curl
    --request POST 'https://rest.uniprot.org/idmapping/run'
    --form 'ids="P99999,P31946"'
    --form 'from="UniProtKB_AC-ID"'
    --form 'to="UniRef100"'

```

This returns a job id as a response

```
{"jobId":"27a020f6334184c4eb382111fbcad0e848f40300"}

```

You then have to poll the api for the status of the job

```
$ curl -i 'https://rest.uniprot.org/idmapping/status/{jobId}'

```

Which returns you a status of FINISHED and the url where you can fetch the result in the HTTP Header. Note that since the job might take some time to finish you need to poll this url until you get the status FINISHED.

```
HTTP/1.1 303
Server: nginx/1.17.7
Location: https://rest.uniprot.org/idmapping/uniref/results/27a020f6334184c4eb382111fbcad0e848f40300
Content-Type: application/json
Access-Control-Allow-Origin: *
...
{"jobStatus":"FINISHED"}

```

Access the result url with curl

```
$ curl -s "https://rest.uniprot.org/idmapping/uniref/results/27a020f6334184c4eb382111fbcad0e848f40300"

```

Which gives you a json dictionary

```
{ "results": [...], "failedIds": []}

```

Each entry in the result list is a dictionary similar to the following example

```
{
    'from': 'P99999'\n\t
    , 'to': {
        'id': 'UniRef100_P99999'
        , 'memberCount': 17
        , 'representativeMember': {
            'uniparcId': 'UPI000013EAA0'
            , ...
        }
        , 'memberIdTypes': [
            'UniProtKB Unreviewed (TrEMBL)'
            , 'UniProtKB Reviewed (Swiss-Prot)'
        ]
        , 'members': [
            'P99999,0'
            , 'Q6WUX8,0'
            , 'P99998,0'
            , 'Q5RFH4,0'
            , 'G4XXL9,1'
            , 'A0A2J8VDX6,1'
            , 'K7CAI8,1'
            , 'A0A663D840,1'
            , 'G4XXM0,1'
            , 'G4XXM2,1'
            , 'G4XXM4,1'
            , 'G1RY04,1'
            , 'G4XXM3,1'
            , 'G4XXM6,1'
            , 'C9JFR7,1'
            , 'A0A2J8VE06,1'
            , 'A0A2J8KZC0,1'
        ]
        , ...
    }
}

```

The uniparc id of the cluster is given in the representative Member.

If you want the group members directly they are given as a list of strings that contain the accession and the memberIdType index separated with a comma. You would use the index to look up the type of the accession in the memberIdTypes array.

The result also contains information such as the sequence of the representative Member, or the number of organisms in the cluster.

**Note:** From what I can tell in the UniRef100 Clusters there will only be one Uniparc Id which you will find with the representative Member. If you want to map to UniRef90 or UniRef50 there will be other Uniparc ids in the member list, however the representative Uniparc Id will still only be in the representative Member, so you it's probably safest to parse both.

---

**Old solution**

I found a way to batch the mapping from uniref cluster ids to the ids of all members. The url `https://www.uniprot.org/uniref/?query=UniRef100_A0A003&format=tab&columns=id,members` can be used to get the link from one uniref Cluster to it's member ids. Multiple uniref cluster can be fetched at once by connecting the ids with a `+OR+` like so `https://www.uniprot.org/uniref/?query=UniRef100_P99999+OR+UniRef100_P31946&format=tab&columns=id,members`

Because the uniref cluster ids are not stable I ended up using the following syntax, where instead of querying the UniRef cluster by it's id, I query by one of it's members, namely the one that was used to create the current uniref cluster id: `https://www.uniprot.org/uniref/?query=*&format=tab&limit=10&columns=id,members&fil=identity:1.0+AND+member:(P99999+OR+P31946)`

In the notes section of the 'Retrieve / ID mapping' documentation of the uniprot they talk about fetching chunks up to 20.000 identifiers at a time. I am not sure if I use the service they are talking about though, so I'm keeping my requests at 100 identifiers per request and that seems to work.

One thing I noticed was that fetching results via the uniref endpoint is faster for small sets of identifiers (\<20) since the fetch time of the batch request is similar for smaller and larger (\<100) sets. So batching only makes sense for more than 20 identifiers at a time. (Disclaimer: I have only done very basic tests on this)

Once I have all the members of a cluster I can use the proteins/uniparc endpoint of the proteins api to fetch the uniparc data. For uniparc accessions (starting with the string 'UPI') I use the `GET /uniparc/upi/{upi}` endpoint. For the other accession, which should all be Uniprot accessions, I use the `GET /uniparc/accession/{accession}` endpoint.

This doesn't answer the question of a direct mapping between the uniref ids and uniparc member ids, but it is a workable solution.
> 2 votes

---
Tags: python, uniprot, embl
---


title: median normalization for proteomics

text:
# Question
Title: median normalization for proteomics
I am using the data from a proteomics study were the data was log2 transformed and then a median normalization was applied. The data was normalized by groups of conditions (normal, mutant), not for all the arrays at the same time, so the median for the normal group is slightly higher than for the mutant group. I am wondering if this is correct or should I normalized the data en general, I mean taking into account all the arrays and normalizing all at the same time so the 2 groups have the same median.

Regards

# Answer
> 2 votes
The kind of normalization depends on what you what to explore. So, there is no absolute answer to this. Different normalizations highlight different aspects of your dataset. There is a nice paper for a higher throughput quantitative study that I would recommend reading:

Nusinow, David P., John Szpyt, Mahmoud Ghandi, Christopher M. Rose, E. Robert McDonald 3rd, Marian Kalocsay, Judit Jané-Valbuena, et al. 2020. “Quantitative Proteomics of the Cancer Cell Line Encyclopedia.” Cell 180 (2): 387–402.e16.

What you want to aim for is a normalization that removes batch effects and that clusters together your quality control samples that you know are very similar to each other.

---
Tags: normalization, proteomics
---


# Question
Title: Targeted hybrid capture bam files showing reads outside target regions in my bed file
Sorry if this is a stupid question I am new to analysing targeted hybrid capture data.

I am analysing my HTS data from a targeted hybrid capture enrichment method and subsequent sequencing on a NextSeq.

Looking at my bam files, I see there are reads that are overlapping into my non target regions.

For example - my .bed file has two target regions in chr17:

However I still have reads overlapping between 7578739-7578918 (following example is position 7578850):

Quite a lot of these reads don't hit any part of the targeted regions in my bed file. Are these mistakes? How have they been pulled down for sequencing when they aren't in the target regions of the baits?

# Answer
> 3 votes
It strongly depends on exactly what your protocol is. Note that there are a good number of tools for QCing capture data (here, here).

Off-target reads are a standard observation in capture experiments; recall that it is target **enrichment**, not target **perfect purification**, as commenter swbarnes2 suggests.

The real question is whether the proportion of those off-target reads is too high for your application, or otherwise suggests problems. Those depend on details of the technique that you will likely have to figure out yourself.

Here are some tables from the first of those QC papers (supplement), showing on-target proportions for various capture samples (sorry, it would have been a bit of effort to get text of the tables rather than images):

These data suggest that on-target reads are 60%-80% of mapped reads. Note that they also lose a lot of reads in mapping, so we cannot exclude the possibility of mis-mapped garbage reads contributing to off-target mapping.

---
Tags: sam, hts
---


# Question
Title: How to get the genomic sequences from a blat result?
Assume we have a `query.fa` file that contains sequences and we run:

`blat -stepSize=5 -repMatch=2253 -minScore=20 -minIdentity=0 -out=pslx /genomes/mm10.fa.qz query.fa output.pslx`

the output `output.pslx` file looks like this:

```
match   mis-    rep.    N's     Q gap   Q gap   T gap   T gap   strand  Q               Q       Q       Q       T               T       T       T       block       blockSizes      qStarts  tStarts
        match   match           count   bases   count   bases           name            size    start   end     name            size    start   end     count
---------------------------------------------------------------------------------------------------------------------------------------------------------------
20      0       0       0       0       0       0       0       +       seq     20      0       20      chr9    124595110       44046930        44046950   20,      0,      44046930,       aaaagtatcagtgtgtatag,   aaaagtatcagtgtgtatag,
20      0       0       0       0       0       0       0       +       seq     20      0       20      chr9    124595110       44046930        44046950   20,      0,      44046930,       aaaagtatcagtgtgtatag,   aaaagtatcagtgtgtatag,

```

What would be a reasonable way to get the genomic contexts (5bp upsteam and 5bp downstream) for each aligned sequence.

For example, assume that `blat` found that the seq: `AAATTGGGGAAAA` aligns to `chr2:100-113`, so the question is how to get `chr2:95-118` easily without reinventing the wheel.

---

I couldn't make it work with bedtools, because my genome's index file is corrupted, but this should work for others who have successfully used `bwa` or `samtools` to index their reference genome:

```
blat -stepSize=5 -repMatch=2253 -minScore=20 -minIdentity=0 -out=pslx /genomes/mm10.fa.qz query.fa output.pslx
awk 'NR>5 {print $14 "\t" $16-10"\t" $17+10}' output.pslx > regions.bed
bedtools getfasta -fi /genomes/mm10.fa.gz -bed regions.bed

```

# Answer
> 1 votes
Via BEDOPS `convert2bed` (`psl2bed`) and `bedops` operations:

```
$ psl2bed < hits.psl | bedops --range 5 --everything - > answer.bed

```

The file `answer.bed` will contain target intervals from the PSL (BLAT) input, padded up- and downstream by five bases.

This BED file can be run through `samtools faidx` or similar to get sequence data.

References:

# Answer
> 0 votes
I couldn't get it to work with varied: `pcl2bed`, `samtools`, `bedtools` commands. I think this biopython should work:

```
import gzip
from Bio import SeqIO
n=5

exact_matches = {}
blat_qresult = SearchIO.read('blat_output.pslx', 'blat-psl', pslx=True)
for hit in blat_qresult:
    for f in hit.fragments:
        if f.hit_span == f.query_span:
            chrom = hit.id
            if chrom not in exact_matches:
                exact_matches[chrom] = []
            exact_matches[chrom].append(f.hit_range)

with gzip.open('/genomes/mm10.fa.gz', "rt") as handle:
    for record in SeqIO.parse(handle, "fasta"):
        if record.id in exact_matches:
            for subseq in exact_matches[record.id]:
                print(record.seq[subseq[0]-n:subseq[1]+n])

```

---
Tags: python, bedtools, data-preprocessing, blat
---


# Question
Title: Running GSEA or comparable analysis with multiple variable interaction?
I am working with an RNA-seq dataset generated from a complex experimental design. Our main question is focused on the interaction between two variables, one (INJ) with two conditions (L,S) and the other (SOC) with four conditions (A,B,C,D). We also have a third variable for tissue type (TIS) from 3 tissues (H,N,T), although we are not interested in this variable for this particular question.

We were able extract DE genes underlying the interaction of interest in DESeq2 using an LRT test. We also ran WGCNA to get gene clusters and looked at trait correlations. Using these WGNCA clusters and VST normalized count data, I would like to run GSEA or a similar analysis to find WGCNCA clusters that differ based on the interaction between the two main variables.

What I've tried so far was generating a pseudo-variable that is a combination of INJ and SOC with 8 separate conditions for running in GSEA, but this doesn't quite get at the interaction between variables since we really want to see if the effect of INJ differs between SOC conditions.

One thing I was thinking was running a GSEA-like analysis on log2 ratios between INJ states for each SOC condition, but these would have to be summary statistics for groups of samples and I would lose the variance of expression within groups (maybe that doesn't matter though for exploratory analysis?). This would also reduce my "sample" count to 1 per condition, which isn't ideal for this type of analysis and I would probably have to look at differences between gene sets instead of phenotypes for more robust stats.

Any thoughts on ways to do this? I don't need to use GSEA if there is another analysis I could run on gene sets that allows for more complex experimental models. Thanks in advance.

# Answer
> 2 votes
I'm not sure I understand "I would like to run GSEA or a similar analysis to find WGCNCA clusters that differ based on the interaction between the two main variables". I would run an association analysis (regression) for module eigengenes and the appropriate interaction term. The (most) significant modules are your candidates. This step is simply regression (in the simplest case a linear model should be sufficient), no GSEA needed. You can then subject the identified modules to GSEA or similar to get an idea of their biological function.

---
Tags: rna-seq, wgcna, gsea, interactions
---


title: How to run Jupyter script on Slurm HPC

text:
# Question
Title: How to run Jupyter script on Slurm HPC
Now jupyter installed on the server and I am using below code to plot the rarefaction plot but I am still getting some error. could you please suggest how I can get rid of it?

```
    #Run it in jupyter lab or notebook

import seaborn as sns
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import os

#%matplotlib inline

qza = 'alpha-rarefaction-20000.qzv' # path to rarefaction file from qiime2
os.system("unzip qza")

digest = a[1].split('/')[0].replace('  inflating: ','')
inf = digest + '/data/observed_otus.csv'
otus = pd.read_csv(inf,index_col=0, sep=',')
os.system('rm -r $digest')

cols = [col for col in otus.columns if 'iter' not in col]
mean,data = otus[cols],pd.DataFrame(columns=cols)
depths = [col.split('_')[0] for col in otus.columns if 'depth' in col]
otus = otus.drop(cols,axis=1)
otus.columns = depths
for depth in depths:
    mean['ASV'] = otus[depth].mean(axis=1)
    mean['depth']= depth.split('-')[-1]
    data = pd.concat([data,mean])
# here provide colors for each item that will be plotted
pal={'Irrigated':'brown','Rainfed':'red'}
fig,ax = plt.subplots(figsize=(2.5,2),dpi=600,tight_layout=True)
sns.set(style='ticks',rc={"lines.linewidth":0.7,"axes.linewidth":0.5})
# use your column name to plot here
sns.lineplot(x='depth',y='ASV',data=data,palette=pal,hue='Season',sort=False,err_style='bars',\
             dashes=True,style='Season',ci=67)
ax.set_xlabel('Sequencing depth',fontsize=8)
ax.set_ylabel('Observed ASVs',fontsize=8)
ax.tick_params(axis='x', labelrotation=90)
ax.tick_params(axis='both',which='major',length=2,pad=0.5,labelsize=6)
handles, labels = ax.get_legend_handles_labels()
ax.legend(handles=handles[1:],labels=labels[1:],fontsize=5,frameon=False,numpoints=4,borderaxespad=0,handletextpad=0.2,loc=2,)
plt.savefig('Infected_observed_ASVs.png', bbox_inches='tight')

```

Error;

```
 $ python3 myscript.rarefaction.py
Traceback (most recent call last):
  File "myscript.rarefaction.py", line 4, in <module>
    import seaborn as sns
  File "/opt/apps/seaborn/0.11.0/python3-3.7.4/lib/python3.7/site-packages/seaborn/__init__.py", line 2, in <module>
    from .rcmod import *  # noqa: F401,F403
  File "/opt/apps/seaborn/0.11.0/python3-3.7.4/lib/python3.7/site-packages/seaborn/rcmod.py", line 7, in <module>
    from . import palettes
  File "/opt/apps/seaborn/0.11.0/python3-3.7.4/lib/python3.7/site-packages/seaborn/palettes.py", line 9, in <module>
    from .utils import desaturate, get_color_cycle
  File "/opt/apps/seaborn/0.11.0/python3-3.7.4/lib/python3.7/site-packages/seaborn/utils.py", line 11, in <module>
    import pandas as pd
  File "/opt/apps/seaborn/0.11.0/python3-3.7.4/lib/python3.7/site-packages/pandas/__init__.py", line 11, in <module>
    __import__(dependency)
  File "<frozen importlib._bootstrap>", line 983, in _find_and_load
  File "<frozen importlib._bootstrap>", line 967, in _find_and_load_unlocked
  File "<frozen importlib._bootstrap>", line 668, in _load_unlocked
  File "<frozen importlib._bootstrap>", line 638, in _load_backward_compatible
  File "/opt/gridware/depots/54e7fb3c/el7/pkg/libs/matplotlib_python34/1.5.1/gcc-4.8.5+python3-3.4.3+numpy_python34-1.11.3/python/lib/python3.4/site-packages/pytz-2019.1-py3.4.egg/pytz/__init__.py", line 1101, in <module>
  File "/opt/gridware/depots/54e7fb3c/el7/pkg/libs/matplotlib_python34/1.5.1/gcc-4.8.5+python3-3.4.3+numpy_python34-1.11.3/python/lib/python3.4/site-packages/pytz-2019.1-py3.4.egg/pytz/lazy.py", line 101, in _lazy
  File "/opt/gridware/depots/54e7fb3c/el7/pkg/libs/matplotlib_python34/1.5.1/gcc-4.8.5+python3-3.4.3+numpy_python34-1.11.3/python/lib/python3.4/site-packages/pytz-2019.1-py3.4.egg/pytz/__init__.py", line 1098, in <genexpr>
  File "/opt/gridware/depots/54e7fb3c/el7/pkg/libs/matplotlib_python34/1.5.1/gcc-4.8.5+python3-3.4.3+numpy_python34-1.11.3/python/lib/python3.4/site-packages/pytz-2019.1-py3.4.egg/pytz/__init__.py", line 114, in resource_exists
  File "/opt/gridware/depots/54e7fb3c/el7/pkg/libs/matplotlib_python34/1.5.1/gcc-4.8.5+python3-3.4.3+numpy_python34-1.11.3/python/lib/python3.4/site-packages/pytz-2019.1-py3.4.egg/pytz/__init__.py", line 102, in open_resource
  File "<frozen importlib._bootstrap>", line 983, in _find_and_load
  File "<frozen importlib._bootstrap>", line 967, in _find_and_load_unlocked
  File "<frozen importlib._bootstrap>", line 668, in _load_unlocked
  File "<frozen importlib._bootstrap>", line 638, in _load_backward_compatible
  File "/opt/gridware/depots/54e7fb3c/el7/pkg/apps/setuptools/24.0.1/python-2.7.8/python/lib/python2.7/site-packages/setuptools-24.0.1-py2.7.egg/pkg_resources/__init__.py", line 2927, in <module>
  File "/opt/gridware/depots/54e7fb3c/el7/pkg/apps/setuptools/24.0.1/python-2.7.8/python/lib/python2.7/site-packages/setuptools-24.0.1-py2.7.egg/pkg_resources/__init__.py", line 2913, in _call_aside
  File "/opt/gridware/depots/54e7fb3c/el7/pkg/apps/setuptools/24.0.1/python-2.7.8/python/lib/python2.7/site-packages/setuptools-24.0.1-py2.7.egg/pkg_resources/__init__.py", line 2952, in _initialize_master_working_set
  File "/opt/gridware/depots/54e7fb3c/el7/pkg/apps/setuptools/24.0.1/python-2.7.8/python/lib/python2.7/site-packages/setuptools-24.0.1-py2.7.egg/pkg_resources/__init__.py", line 956, in subscribe
  File "/opt/gridware/depots/54e7fb3c/el7/pkg/apps/setuptools/24.0.1/python-2.7.8/python/lib/python2.7/site-packages/setuptools-24.0.1-py2.7.egg/pkg_resources/__init__.py", line 2952, in <lambda>
  File "/opt/gridware/depots/54e7fb3c/el7/pkg/apps/setuptools/24.0.1/python-2.7.8/python/lib/python2.7/site-packages/setuptools-24.0.1-py2.7.egg/pkg_resources/__init__.py", line 2515, in activate
  File "/opt/gridware/depots/54e7fb3c/el7/pkg/apps/setuptools/24.0.1/python-2.7.8/python/lib/python2.7/site-packages/setuptools-24.0.1-py2.7.egg/pkg_resources/__init__.py", line 2097, in declare_namespace
  File "/opt/gridware/depots/54e7fb3c/el7/pkg/apps/setuptools/24.0.1/python-2.7.8/python/lib/python2.7/site-packages/setuptools-24.0.1-py2.7.egg/pkg_resources/__init__.py", line 2047, in _handle_ns
  File "/opt/gridware/depots/54e7fb3c/el7/pkg/apps/setuptools/24.0.1/python-2.7.8/python/lib/python2.7/site-packages/setuptools-24.0.1-py2.7.egg/pkg_resources/__init__.py", line 2066, in _rebuild_mod_path
AttributeError: '_NamespacePath' object has no attribute 'sort'

    metadata file ;
    #SampleID       Villages        Region  Treatment       Batch   Season
    PN0086D.1.S1    Kholabari       Mymensingh      T1      batch1  Irrigated
    PN0086D.2.S2    Bogajan Mymensingh      T1      batch1  Irrigated
    PN0086D.3.S3    Bonkua  Mymensingh      T1      batch1  Irrigated

```

# Answer
> 2 votes
This is a bash script (#!/bin/bash) and you try to run python in it. You should save your python script (e.g. in "myscript.py") and call it within your bash script. Maybe first try to see if you are able to run it in your terminal and then run it on your server with SLURM.

EDIT1: Also I don't think you need 3 nodes to run such script and this will probably not be even possible here. I think you need some specific modules to parallelize a python script with MPI (see https://stackoverflow.com/questions/25792357/how-to-distribute-a-python-script-across-multiple-nodes), but again I think this would be an overkill. Try first with `#SBATCH -c 1` instead.

Edit2: Finally it might be possible to run jupyter-notebook on a node allocated by SLRUM instead of putting this into a python script by using `jupyter-notebook --no-browser --port=${port} --ip=${node}` command in your bash slurm script. Maybe you'll find this doc interesting: https://docs.ycrc.yale.edu/clusters-at-yale/guides/jupyter/ .

# Answer
> 0 votes
I'll make a contribution and keep editing it rather than use 'comments', the error is about

```
import numpy as np

```

At the moment this isn't about Jupyter, its about getting Python working on your server. Is that a good assessment? With Jupyter notebook you write a little bit of code, check it, write a bit more. That's how it is designed to operate. The other way in is to use Python in your shebang statement otherwise it will default to `sh` or `bash`.

The code you supplied doesn't explicitly use `numpy`. Pandas does use numpy and `np` is the standard shorthand, `Multiarray` is probably something pandas would use. Did you get this code working on your local machine?

I am wondering whether the numpy module on your server is up to date? Did you creat a specific Python environment, use `pip` or `conda`? Checking the Python version on your server and whether `numpy` needs updating is perhaps where I might start.

It is not pretty getting a HPC to run with your stuff and can take real work.

---

I read your new post, I would simply reinstall Python from `Anaconda3` and run the script for now from shell pointing the shebang into Anaconda3 (I can't remember if you need to do this, I can check if needed). This a more professional way of using a HPC than using the resident Python. In good old days of Perl, when everything was nice and simple, `perlbrew` was essentially the Perl equivalent and no-one would use the system Perl.

The reason why I'm saying this is I get this feeling it will be an endless cycle of import errors. Anaconda3 will give you control of your Python system without admin privelidge.

Just to explain if you are not calling `np` explicitly (and looking at your code I don't think you are) then explicitly stating `import numpy` isn't needed. The fact you wrote it and it then jumped onto a different error is a bit worrying. I would just use your own version of Python. In Jupyter I would guess, but don't know they use their own Python ... but from a data science point of view Jupyter is the front end to AWS and GCS because those guys want some control of what you are doing, not running riot. In DNN jupyter is big. If you want an absolute diagnosis then Stackoverflow might help, but I think you'll just another bug and another etc... One thing is that `pandas` changes alot and the version you are using might be ahead of the curve, thats my guess

Final point is if you are just plotting the `ggplot2` is really good, but that is `R`. `seaborne` is really good, `matplotlib` but hmmmm, but nothing is as clear or explicit as `ggplot2`. Python is good at lots of things, pandas is fabulous, but plotting isn't on par with `ggplot2`

---

What is 100% clear is that you are loading the HPC version of Python and this is out of date (3.4) should be 3.7 (or even 3.8?). The thing I don't understand is why is it loads Python 2 modules?? This is a complete mystery to me, I thought Python 2 and Python 3 were incompatible, this is suggesting backwards compatibility or else the reason it is falling over. I have no idea, but I wouldn't use the resident Python.

---
Tags: python, bash, hpc, jupyter
---


# Question
Title: NCBI blast for exact match of a short sequence
I'm trying to Blast for exact matches to the sequence: 'ATTGNNNNGCAAACCA' in the human transcriptome using NCBI Blast on its 'refseq\_rna' database. However, when I do a basic query I get "No significant similarity found."

Note the N's in the middle of the sequence, if I get rid of those N's, yielding the sequence 'AGCGGATTGCAAAGCAAACCA', I get a match in the 3' UTR of the human MeCP2 gene (which is correct). I don't understand why adding the N's here make it not work.

I looked at the help section and I believe this may be due to shorter sequences being less "statistically significant". Statistical significance of output is completely irrelevant to me, so I set 'Expect Value' to 100000000000000. However, I still get the same result.

I also tried to modify the HTML of the page to allow me to submit a word size of 4 to the form. This went through, but I still got "No significant similarity found".

Could someone help me out with this? I feel like searching for a short sequence shouldn't be this hard.

# Answer
**UPDATE:** while procrastinating, I worked up the below throwaway code into a more generally applicable (but still bare-bones) script.

---

**ORIGINAL:**

I feel like you might want to use something other than BLAST to do this.

BLAST is not really designed to solve this problem, I think it's just easy to use. Also, your N-containing sequence (16nt) is shorter than the version without Ns (21 nt), which can't be correct if you care about "exact" matches.

As a side note, the idea of exact matches is a bit tricky with Ns involved. The whole point of ambiguity codes is that things won't be exact.

And all of this is on top of the short-string problem- most alignment programs won't even return results that don't meet some minimum alignment score, which is exacerbated a bit by the Ns. For how BLAST handles ambiguity codes, you can see section 7.2 here.

Why not use regular expressions instead to search your sequence against the transcriptome?

```
# python
import re
from Bio import SeqIO

fwd = r'ATTG[ACGT]{4}GCAAACCA'
rev = r'TGGTTTGC[ACGT]{4}CAAT'  # reverse complement

# reference file
seqfile = SeqIO.parse("your_reference_file.fa", "fasta")  # a transcriptome reference fasta

for seq in seqfile:
    # maybe use re.match() to get coords but easier to show it with findall
    a = re.findall(fwd, str(seq.seq))
    b = re.findall(rev, str(seq.seq))
    if len(a) > 0:
        print(seq.id, a)
    if len(b) > 0:
        print(seq.id, b)

```

When I put this in a file and run it against a test case with `ATTGAAAAGCAAACCA`, it works:

```
$ python exact_match.py
the_contig ['ATTGAAAAGCAAACCA']

```

Alternately, you could just blast the 16 non-N versions of that sequence and combine them, as you've shown.
> 2 votes

---
Tags: blast
---


# Question
Title: Truble to run a multiprocessing kmer count script
Hi there I have this code:

```
from functools import partial
import multiprocessing

def kmer_count(sequence, alphabet, k):
    """Returns a dictionary with kmers and it counts."""
    seq = sequence.upper()
    seq_len = len(seq)
    kmers = [seq[i:i+k] for i in range(0, seq_len - k + 1)]
    filterd_kmers = [kmer for kmer in kmers if
                     all(base in set(alphabet) for base in kmer)]
    return Counter(filterd_kmers)

pool = multiprocessing.Pool(4)

args = [''.join(alphabet.iupac_dna), 4] # in the real code it came from argparse
f = partial(kmer_count,seq)
pool.starmap(f, *args)

```

Result:

```
---------------------------------------------------------------------------
RemoteTraceback                           Traceback (most recent call last)
RemoteTraceback:
"""
Traceback (most recent call last):
  File "/opt/anaconda3/lib/python3.7/multiprocessing/pool.py", line 121, in worker
    result = (True, func(*args, **kwds))
  File "/opt/anaconda3/lib/python3.7/multiprocessing/pool.py", line 47, in starmapstar
    return list(itertools.starmap(args[0], args[1]))
TypeError: kmer_count() missing 1 required positional argument: 'k'
"""

The above exception was the direct cause of the following exception:

TypeError                                 Traceback (most recent call last)
<ipython-input-79-0f5ba8381a51> in <module>
      1 f = partial(kmer_count,seq)
----> 2 pool.starmap(f, *args)

/opt/anaconda3/lib/python3.7/multiprocessing/pool.py in starmap(self, func, iterable, chunksize)
    274         `func` and (a, b) becomes func(a, b).
    275         '''
--> 276         return self._map_async(func, iterable, starmapstar, chunksize).get()
    277
    278     def starmap_async(self, func, iterable, chunksize=None, callback=None,

/opt/anaconda3/lib/python3.7/multiprocessing/pool.py in get(self, timeout)
    655             return self._value
    656         else:
--> 657             raise self._value
    658
    659     def _set(self, i, obj):

TypeError: kmer_count() missing 1 required positional argument: 'k'

```

If I run like this:

```
f(*args)

```

I got this:

```
Counter({'GTTC': 9780,
         'TTCG': 11291,
         'TCGC': 14803,
         'CGCC': 14428,
         'GCCA': 19254,
         'CCAG': 11512,
         'CAGA': 11641,
         'AGAG': 10387,
         'GAGC': 13164,
         'AGCG': 16025,
         'GCGG': 13067,
         'CGGT': 12329,
         'GGTT': 14217,
         'GTTT': 18733,
         'TTTT': 25488,
         'TTTG': 22006,
         'TTGA': 20004,
         'TGAC': 11351,
         'GACT': 7262,
         'ACTA': 5663,
         'CTAG': 4190,
         'TAGC': 8543,
         'AGCT': 12898,
         'GCTT': 19114,
         'CTTG': 16557,
         'TGAA': 17703,
         'GAAC': 9334,
         'AACA': 14050,
         'ACAC': 8812,
         'CACA': 11050,
         'ACAT': 9705,
         'CATC': 17357,
         'ATCC': 11157,
         'TCCC': 5754,
         'CCCG': 5926,
         'CCGT': 8601,
         'CGTC': 7401,
         'GTCC': 3811,
         'TCCT': 5890,
         'CCTG': 7985,
         'CTGC': 15324,
         'TGCG': 14856,
         'TTTA': 15607,
         'TTAG': 8366,
         'TAGG': 5393,
         'AGGC': 9580,
         'GGCA': 13939,
         'GCAC': 11686,
         'CACC': 14496,
         'ACCA': 16950,
         'CCAC': 13028,
         'CACT': 13194,
         'ACTT': 12585,
         'TTGC': 18877,
         'TGCA': 13770,
         'GCAT': 13307,
         'CATT': 14102,
         'ATTA': 9776,
         'TTAA': 12585,
         'TAAG': 9048,
         'AAGC': 18834,
         'AGCC': 12734,
         'CCAT': 13799,
         'CATA': 8188,
         'ATAA': 12278,
         'GCGT': 12546,
         'CGTT': 12599,
         'GTTG': 15136,
         'TTGG': 18903,
         'TGGC': 19391,
         'GGCC': 7685,
         'GCCC': 7404,
         'CCCA': 9752,
         'CCAA': 18755,
         'CAAC': 14852,
         'AACT': 11886,
         'ACTC': 10028,
         'CTCA': 13518,
         'TCAC': 15385,
         'ACCG': 12027,
         'CCGC': 13105,
         'CGCA': 14853,
         'ATAG': 7066,
         'AGGG': 5273,
         'GGGG': 4927,
         'GGGT': 7972,
         'GAAA': 16511,
         'AAAC': 18215,
         'ACTG': 12639,
         'TGCT': 16663,
         'CTTT': 18228,
         'TTTC': 17497,
         'TCGG': 10116,
         'TGAT': 20593,
         'GATG': 17301,
         'ATGA': 14218,
         'GATC': 14479,
         'ATCG': 16604,
         'CGCG': 13405,
         'CGTA': 8086,
         'GTAG': 6829,
         'TAGA': 6656,
         'AGAT': 11833,
         'ATCA': 20344,
         'TCAG': 13868,
         'ATCT': 11892,
         'TCTT': 14706,
         'CTTC': 13897,
         'TTCC': 9522,
         'TCCG': 7099,
         'CCGG': 5935,
         'CGGG': 5786,
         'GGGA': 5615,
         'GGAT': 11054,
         'GATA': 12085,
         'ATAT': 9205,
         'TATC': 11941,
         'TCGT': 8802,
         'GTAT': 8661,
         'TATT': 12620,
         'ATTT': 18149,
         'TGAG': 13150,
         'GAGA': 8658,
         'GAGG': 5644,
         'GCCT': 9673,
         'CCTT': 11377,
         'CTTA': 9100,
         'TTAC': 11418,
         'TACG': 8250,
         'ACGC': 12774,
         'CACG': 10171,
         'ACGA': 8951,
         'CGAT': 16617,
         'ATTG': 17372,
         'GGCG': 14257,
         'TTGT': 13287,
         'TGTA': 7808,
         'ACGG': 8529,
         'CGGA': 7109,
         'ATGT': 9975,
         'TGTC': 8286,
         'GTCG': 9088,
         'TCGA': 11986,
         'GATT': 16232,
         'GACC': 7608,
         'CATG': 11552,
         'ATGG': 13703,
         'TGGG': 9728,
         'CGAC': 9020,
         'GACG': 7394,
         'ACGT': 8085,
         'GTAA': 11360,
         'TAAT': 9895,
         'AATA': 12342,
         'GGCT': 12713,
         'AAAT': 17834,
         'AATC': 16220,
         'TTCA': 18526,
         'TCAA': 20304,
         'CAAA': 21752,
         'AATG': 14200,
         'TGGT': 16952,
         'GTTA': 10179,
         'GCGA': 14624,
         'AGGT': 9338,
         'GGTA': 10326,
         'CGGC': 12423,
         'TAAC': 9956,
         'AACG': 12524,
         'GCGC': 18306,
         'CTGG': 11344,
         'TGGA': 10805,
         'GCTG': 17455,
         'CTGA': 13532,
         'CGAG': 8649,
         'GAGT': 10014,
         'AGTA': 8668,
         'GTAC': 8038,
         'TACC': 10217,
         'ACCT': 9204,
         'ATAC': 8527,
         'TACA': 7453,
         'ACAA': 12591,
         'AAAA': 24650,
         'AAAG': 17848,
         'AAGG': 11240,
         'AGTT': 12431,
         'CAAT': 17456,
         'GTCT': 6235,
         'AGGA': 5937,
         'GGAG': 5294,
         'AGAA': 12516,
         'AAGA': 13980,
         'CGTG': 9948,
         'GTGC': 11762,
         'TCAT': 14388,
         'AGTC': 7451,
         'GTCA': 11557,
         'CAGT': 13149,
         'GGTC': 7553,
         'GTGA': 15043,
         'TAAA': 15634,
         'GGGC': 7200,
         'GGTG': 14495,
         'TTCT': 13230,
         'TCTA': 6621,
         'CTAC': 6510,
         'TACT': 8573,
         'CTAT': 6856,
         'TATG': 7993,
         'ATTC': 11395,
         'GAAG': 13449,
         'AAGT': 12578,
         'GCTA': 8346,
         'CTAA': 8310,
         'TAGT': 5859,
         'GACA': 8060,
         'TTAT': 12293,
         'TCCA': 11138,
         'CAGC': 17389,
         'GCTC': 13661,
         'CTCG': 8724,
         'CCCC': 5036,
         'TGCC': 13824,
         'GCCG': 12340,
         'GGAA': 9364,
         'CTGT': 9651,
         'TGTT': 14580,
         'ATGC': 13150,
         'AGTG': 13050,
         'GTGT': 9052,
         'TCTG': 11772,
         'AATT': 13738,
         'CGAA': 11026,
         'GTGG': 12927,
         'TCTC': 8866,
         'ACAG': 9505,
         'CCTC': 5686,
         'CTCC': 5391,
         'CCTA': 5236,
         'CGCT': 16302,
         'AACC': 13897,
         'ACCC': 8037,
         'GCAA': 18697,
         'CAAG': 16287,
         'AGCA': 16273,
         'GCAG': 15145,
         'CCGA': 9751,
         'CAGG': 7851,
         'TATA': 4461,
         'GGAC': 3754,
         'CTCT': 10608,
         'TGTG': 11291,
         'AGAC': 6199,
         'CCCT': 5517,
         'GAAT': 11315}

```

What I am doing wrong with the multiprocessing step. I want to speed up my analysis, because I have a lot of genomes in a directory to process. I just need to figure out how can I implement the pool, because probably all other functions a similar.

Thanks for your time and attention

# Answer
This code help me to get the answer I was looking for.

```
r = pool.apply_async(f, a)
print(r.get())

```
> 2 votes

---
Tags: python, genome, k-mer
---


# Question
Title: Identifying nested/overlapping reading frames in nucleotide sequences
I'll start by saying I'm a total stranger to this field, so I'd love to be corrected if I misunderstood something along the way.

I'm writing a python code to identify all the reading frames in a sequence of nucleotides. To my understanding, a reading frame is defined by a start codon and a stop codon, which made total sense to me until I learned you can have nested/overlapping reading frames.

My question is:

Assuming the start codon "AUG", and the stop codon "UAA":

**AUG** AAA GGG **AUG** CCC UUU **UAA** AGA CUC **UAA**

How do you identify the reading frames? Is each stop codon referring to a single start codon? If so, how can you tell to which one?

I tried looking this up and could only come up with articles explaining the effects of overlapping reading frames, rather than how are they interpreted. I would greatly appreciate even a nudge in the right direction.

Please let me know if I was not clear enough or if I got something wrong.

Thanks!

EDIT: After reading the answers here and researching further I learned that ORF identification on DNA sequences isn't as simple as finding a single start and stop codon, and is complicated beyond the scope of my question. I'll probably default to using online tools to do this from now on.

# Answer
There are 6 reading frames. Three forward, three in the reverse direction.

An Open Reading Frame, an ORF, is a sequence which begins with a start codon, and ends with a stop **in the same frame**. You find them by looking for sequences that fit that pattern.

> Is each stop codon referring to a single start codon? If so, how can you tell to which one?

You could have two different transcripts for the same gene with different start codons but the same stop codon. But in general, the protein translation apparatus will not accidently read through a stop codon.

And yeah, I hope you are doing this for fun, because obviously people have already made publicly available tools to find ORFs.
> 0 votes

# Answer
An RNA sequence has six possible reading frames depending on the offset (0, 1, 2) and the direction (forward, reverse complement). All frames can potentially be used for protein translation, but only one frame per protein (i.e. there is no shifting going on when translation is happening, excluding shifts introduced by splicing).

If you need to find all the possible protein sequences within an RNA sequence, there are plenty of programs already available (e.g. getorf and transseq from the EMBOSS tool kit).
> 0 votes

---
Tags: rna-seq
---


# Question
Title: Is there an established method for comparing expression of groups of genes (gene sets)?
I have a figure in a paper where I show the log2(fold-change) (obtained with DESeq2) between two groups (based on genotype) for genes within a specific hallmark gene set, and a reviewer is asking whether the differences are statistically significant.

I have not been able to find any methods for comparing expression of groups of genes like this, and i'm sure there is a method or a reason why it isn't done that I'm not aware of.

To be extra clear, what i have is two sets of log2FC values for a set of genes (one for one genotype (drug-treated vs. untreated) and the same for a second genotype), and I'm looking for a statistical comparison of these two groups of Log2FC measurements.

Thanks

# Answer
You could apply the GSVA transformation to later compare them. You would need to give a step back and start from the expression values not log2FC, but I think it would be better than using the log2FC directly.

If you need to use them you could compare the gene set enrichment value of these two sets on the same comparison. I would use fgsea with the log2FC and the two sets of genes. If the normalized enrichment value (NES) is different then you could explain that these sets expression are indeed different.
> 2 votes

---
Tags: differential-expression, deseq2, gsea
---


# Question
Title: Single cell RNAseq cell cluster (true cluster or sub cluster)
I am trying to run seurat on ~5000 sinngle cells. I am expecting minimum 15 cell types to be present in the data. I tried to runn it with multiple conditions;I can see there is 27 clusters; I believe many of them might be sub clusters. Is there a way to know which clusters are actually a sub-clusters.

As suggested by the members I have used clustree to check the cluster stability also changes the clustering parameter.

# Answer
> 2 votes
I wouldn't be surprised if the 15 cell types that you are expecting were mostly characterized via "protein-based" methods; antibody staining, flow cytometry using fluorophore coupled antibodies, ... The data that scRNA-seq is RNA based and should not be expected to provide the same information as protein based assays. Moreover, the scRNA-seq data, especially those from microfluidic systems like 10x, are sparse and might not be enough to resolve closely related cell types/states such as those of T cells.

On top of the aforementioned, evaluating clusters or clustering stability is a difficult task and in my experience, quite some cluster evaluation metrics cannot be easily applied to single cell data simply because of the shear size of the data. Having asked your question and looked for answers (and I am still doing that), I am doing the following:

i) Try to use transcriptomic markers as much as possible. Even then a lack of a marker would not mean much in terms of scRNA-seq, it might very well not be detected simply due to sensitivity

ii) Use clustree do select the most plausible number of clusters (resolution parameter in Seurat). The package is compatible with Seurat and some other scRNA-seq packages.

iii) Use silhouette width as a metric for clustering. This is computationally expensive, however, should be fine for your 5000 cells. For larger cell numbers, I use another package that approximates this (will add a link but first I have go through my scripts for the name).

iv) Check if batch effects are accounting for what you call "sub-clusters". For example clusters 0, 14, 24 and 25 above might correspond to different samples processed in different days.

# Answer
> 1 votes
You should browse the web for methods to assess cluster stability. The idea is to permutate clusters by randomly adding and subtracting cells of a given cluster and see whether it collapses or not. I would bet money that many of your clusters are not stable, e.g. 17 and 2. Try to change clustering parameters to decrease resolution in order to get larger clusters. I am not a Seurat user but typically there are resolution and nearest-neighbor parameters, check the documentation. See for example section 10.6 in https://osca.bioconductor.org/clustering.html as well as https://www.biorxiv.org/content/10.1101/2020.05.26.116640v1 for inspiration towards cluster stability assessment.

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


# Question
Title: Mapping statistics from bam file using bbtools and sambamba
Below are the statistics for RNA-seq mapped and unmapped paired-end reads to rice genome using reformat.sh from bbtools on bam files. It gives 77% mapped and 5% unmapped, what about the remaining 18% reads? How I can get information for all reads?

```
Command;
reformat.sh in=Leaf_T1_F_R10_S1_L001.bam out=Leaf_T1_F_R10_S1_L001.mapped.bam, mappedonly

Unspecified format for output Leaf_T1_F_R10_S1_L001.mapped.bam,; defaulting to fastq.
Found sambamba.
Input is being processed as unpaired
Input:                   67471075 reads           9312399845 bases
Output:                 52020538 reads (77.10%) 7208830835 bases (77.41%)

Time:                         171.212 seconds.
Reads Processed:      67471k 394.08k reads/sec
Bases Processed:       9312m 54.39m bases/sec

For unmapped read;
Command;
reformat.sh in=Leaf_T1_F_R10_S1_L001.bam out=Leaf_T1_F_R10_S1_L001.unmapped.bam, unmappedonly

Unspecified format for output Leaf_T1_F_R10_S1_L001.unmapped.bam,; defaulting to fastq.
Found sambamba.
Input is being processed as unpaired
Input:                   67471075 reads           9312399845 bases
Output:                 3435326 reads (5.09%) 465456267 bases (5.00%)

Time:                         142.978 seconds.
Reads Processed:      67471k 471.90k reads/sec
Bases Processed:       9312m 65.13m bases/sec

```

With sambamba; Can I say 94.91% reads are mapped and 5.09% are unmapped?

```
67471075 + 0 in total (QC-passed reads + QC-failed reads)
12015211 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
64035749 + 0 mapped (94.91%:N/A)
55455864 + 0 paired in sequencing
27727932 + 0 read1
27727932 + 0 read2
51006046 + 0 properly paired (91.98%:N/A)
51419222 + 0 with itself and mate mapped
601316 + 0 singletons (1.08%:N/A)
263816 + 0 with mate mapped to a different chr
237111 + 0 with mate mapped to a different chr (mapQ>=5)

```

# Answer
> 1 votes
Try `samtools flagstat` or `sambamba flagstat` to get similar type of information about your bam file. It may be easier to interpret. Both of these tools can be installed using `conda` if needed.

**EDIT:**

Yes, the results of `sambamba flagstat` you showed in your edited post indicate exactly what you said: `94.91%` mapped and `100% - 94.91% = 5.09%` unmapped:

```
64035749 + 0 mapped (94.91%:N/A)

```

# Answer
> 0 votes
More of a comment than an answer... In my opinion the terminology is a bit confusing since the term *reads* actually refers to a more generic *records*.

Look at this test file:

```
samtools view test.sam | cut -f 1-4
HSQ9103:404:C6F0VANXX:1:2208:4363:50381 0   chr7    5529094
HSQ9103:404:C6F0VANXX:1:2208:4363:50381 0   chr7    5529194
HSQ9103:404:C6F0VANXX:1:2208:4363:50381 0   chr7    5529194
HSQ9103:404:C6F0VANXX:1:2208:4363:50381 0   chr7    5529194
HSQ9103:403:C6F0HANXX:8:2215:4792:20065 4   chr7    5529871

```

I would say that here there are 5 records (of which 4 are alignments) but only 2 **reads**.

`samtools flagstat` reports 80% mapping rate but maybe I would say this is 50%.

```
samtools flagstat test.sam
5 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
4 + 0 mapped (80.00% : N/A)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (N/A : N/A)
0 + 0 with itself and mate mapped
0 + 0 singletons (N/A : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

```

---
Tags: rna-seq, samtools, bash, java, sambamba
---


# Question
Title: e-utils api: get chromosome position for each snp on a gene
So I have the following e-utils api-link:

> https://eutils.ncbi.nlm.nih.gov/entrez/eutils/elink.fcgi?dbfrom=gene&db=snp&id=5726&retmode=json

it returns all the SNP-ids for the gene with the id 5726 as an example. Additionally to the SNP-ids I also need the SNP position and the genomic distance. Also I would need the chromosome but I could query that from the gene as well rather than for each SNP. Do you know how I can add all these informations to the resulting json/xml?

The first one in the resulting list is the SNP with the id 1490381980. Further information about it can be found on dbSNP: https://www.ncbi.nlm.nih.gov/snp/?term=1490381980%5BUID%5D

I really want to avoid creating a request for each single SNP and would prefer to query everything at once. Due to the unavoidable performance-issues. The Tas2R38 doesn't have that many SNPs but there are other genes with way more.

# Answer
If your dataset is limited to human SNPs, and you know the list of SNPs, I would advise taking a look at the dedicated NCBI Variation Services. This would be one request for each SNP, but if you fire off 1-2 per second, I think you'll be ok.

For example, for rs1490381980, you could try something like:

```
curl -X GET 'https://api.ncbi.nlm.nih.gov/variation/v0/refsnp/1490381980' \
   | jq '.primary_snapshot_data.placements_with_allele[] | "\(.seq_id) \(.alleles[0].allele.spdi.position)"'

```

It would be better to query for the SNP and parse the JSON in python, but `jq` helps provide an idea of what is possible.

Alternatively, you could download the JSON or VCF data from the NCBI dbSNP FTP site, and process out the information you need.
> 2 votes

---
Tags: snp, phylogenetics, api
---


# Question
Title: Should one be especially conscious of removing low counts, when applying batch correction?
When (as a mistake) I did **not** remove low counts at all (beside those that equal zero for all samples), I got the following ma plot (using Glimma) :

On the right you see the individual counts. The gene that is displayed **is** signficant. However, the reality seems that the gene is not expressed both in treatment and in control. The batch I have corrected for made it a significantly DE however.

The batch correction algorithm of DESeq2 doesn't care whether the counts of the gene are almost all zeroes, and does not take this as a special case. In reality, it seems to me that it **is** a special case - if it is questionable whether the gene is expressed in the first place, it makes no sense to apply batch correction to zero counts.

The solution to this seems to be to make sure to remove low counts when applying batch correction. Usually it's only a recommendation; but when applying batch, it seems a necessity.

Would you agree?

I even wonder if one should remove counts that are low per-condition. Suppose that a gene is not expressed in treatment (and only in treatment), does it really make sense to apply batch correction to this gene?

# Answer
> 1 votes
I don't think this has anything to do with batch correction.

The averages of that gene are distinctly different between groups, so it's correct that the gene is statistically significant. But with one sample driving all the difference, it's likely to be an artifact, or otherwise not fruitful to follow up on.

As geek\_y suggested, a filter like "omit genes unless 2 samples have \> 10 counts" would drop that gene.

---
Tags: batch-effects
---


# Question
Title: if I know the number of sequencing circles can I give this information to DESeq2?
I am trying to understand library normalization in DESeq2. I would like to ask the following: I know that some samples have been run 15 cycles and some others 20, can I give this information to DESeq2, would it be useful? I mean, I know that DESEQ2 uses an algorithm to normalize library size, but since I have this piece of information, I wonder if it would be useful to provide it somehow to DESEQ2. If it is not, would it make sense to incorporate this information in the design matrix, something like: design(~cycle+condition)? Thank you a lot!

# Answer
> 2 votes
You can include that info in the design, though you'd have to decide whether to include it as a number of just a factor. I've never heard of anyone including PCR cycles. I've heard of people including RIN scores, though I've never done it.

If the PCA looks alright, I wouldn't worry about. I've got tons of examples of the lab people saying "I did a few more cycles on these few samples" and then you can't tell any difference in the PCA. I think needing to do more cycles likely has more to do with variance in the prep than a real biological difference.

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


# Question
Title: Comparing groups of 16S rRNA sequences - how to?
I have several groups of 16S rRNA sequences associated with taxonomic groups (let's call them A, B and C). The sequences within each group have a common ancestor and are on average more closely related to sequences within the same group. The groups are also related, but less.

Ideally I'd like to have a table like this (where percentages represent the average similarity):

```
           Group A     Group B     Group C
Group A      99%
Group B      97%         98%
Group C      94%         95%         99%

```

Is there software that can calculate these percentages for me? I'm thinking that I could set up a many-against-many BLAST search and then try to use pandas in python in order to parse the resulting table...But if something like this already exists, I'd rather not reinvent the wheel. Thank you!

# Answer
If you want to make an all vs all with blast this can be helpful. You can create a blast database with your sequences and then search against itself.

```
makeblastdb -in 16S_sequences.fasta -dbtype nucl -out my_16S_sequences_db

blastn -db my_16S_sequences_db -query 16S_sequences.fasta -outfmt 6 \
       -out allvsall_16S.tsv -num_threads "$(nproc)"

```

If you want a percent identity matrix for a pairwise alignment you can try with this approach: https://www.biostars.org/p/220154/
> 1 votes

---
Tags: phylogenetics, taxonomy, sequence-homology
---


# Question
Title: What is a quick way to find the reverse complement in bash
I have a DNA sequence for which I would like to quickly find the reverse complement. Is there a quick way of doing this on the bash command line using only GNU tools?

# Answer
> 32 votes
Thanks to Manu Tamminen for this solution:

```
echo ACCTTGAAA | tr ACGTacgt TGCAtgca | rev

```

# Answer
> 7 votes
Reverse complement FASTA/Q:

```
seqtk seq -r in.fa > out.fa

```

https://github.com/lh3/seqtk

# Answer
> 2 votes
You may either install the emboss package (http://emboss.sourceforge.net/download/) and use the revseq program. It contains many useful command-line programs such as the pairwise alignment tool needle.

---
Tags: sequence-analysis, phylogenetics, bash
---


# Question
Title: How to query NCBI (Nucleotide database) by a feature qualifier?
I cannot seem to find a way to query nucleotide database (https://www.ncbi.nlm.nih.gov/nuccore) by a specific qualifier of the feature. For example, by /host under the feature "source". For example, below is a part of one record, where the qualifier /host contains "Bemisia tabaci B". How could I search for all sequences that contain the word "Bemisia" in the qualifier \host?

```
FEATURES             Location/Qualifiers
 source          1..351658
                 /organism="Candidatus Portiera aleyrodidarum BT-B-HRs"
                 /mol_type="genomic DNA"
                 /strain="BT-B"
                 /host="Bemisia tabaci B"
                 /db_xref="taxon:1206109"
 gene            1..1926
                 /gene="dnaK"
                 /locus_tag="C530_RS00005"
                 /old_locus_tag="C530_001"
 CDS             1..1926
                 /gene="dnaK"
                 /locus_tag="C530_RS00005"
                 /old_locus_tag="C530_001"
                 /inference="COORDINATES: similar to AA
                 sequence:RefSeq:WP_014894920.1"
                 /note="Derived by automated computational analysis using
                 gene prediction method: Protein Homology."
                 /codon_start=1
                 /transl_table=11
                 /product="molecular chaperone DnaK"
                 /protein_id="WP_014894920.1"
                 /translation="MIRIIGIDLGTTNSCVAVLENGKSRVIENSEGARTTPSIIAYTD
                 DNEILIGQTAKRQAVTNPSNTLYAIKRLIGRKFNDDIVQKDIKMVPYKIVSAENGDAW
                 VQVKNRKLAPPQISAEILKKMKKTAEDYLGEKITEAVITVPAYFNDSQRQATKDAGRI
                 AGLNVKRIINEPTAAALAYGMDKNRGDKKIAVYDLGGGTFDISIIEIADVDGETQFEV
                 MATNGDTFLGGEDFDLKLINYIVKQFKLESGVDLSGDSLAMQRLKEAAEKGKIELSSS
                 QQTDINLPYITADKTGPKHLTLKITRAKLESLVEDLIKRSLLPCKKAIDDAKLTTKEI
                 DDIILVGGQTRMPLVQKKVAEFFGKEARKDVNPDEAVAMGAAIQGGVLGGDVKDVLLL
                 DVTPLTLGIETMGGIMTPLIEKNTTIPTKKTQIFSTAEDNQTAVTIHALQGERKKALQ
                 NKSLGRFDLTDIPPAPRGTPQIEVAFDLDANGILNITAKDKATGKEQSIIIKSSGGLS
                 EEEINKMIQDAETNKEEDKKFEELVQTRNQADGMIHTIRKTLKEIKNEDEKNKLESLI
                 NNLEQTLKSDDIKTIKEDLNKLTEEYSTIYQKIYSEKKPQNHNSEPAQKNAEKSSSKK
                 EGDVVDAEYEEVNNKTK"
                 ...

```

Best, Milan

# Answer
> 1 votes
For this particular query, try host=Bemisia tabaci B\[WORD\]. However, you may also be interested in related queries like host=Bemisia tabaci B biotype\[WORD\].

But, I am not sure a general solution exists, as I do not believe all qualifiers are indexed.

---

Also, you can see how your query is interpreted on the right side under "Search details". For the original attempt (`[host]"Bemisia tabaci B"`), it is interpreted as:

```
host[All Fields] AND "Bemisia tabaci B"[All Fields]

```

while the goal is to produce

```
host=Bemisia tabaci B biotype[WORD]

```

---
Tags: database, phylogenetics, entrez
---


# Question
Title: Does rRNA depletion protocol give higher number of mapped reads in Intronic regions?
Recently, I have downloaded a publicly available dataset, which are 350 tumor samples. I see the following information from the published paper.

They used Ribo Zero Gold and rRNA was depleted. Strand specific data. After aligning the data I did some alignment quality check with Qualimap RNA-Seq QC tool. I visualised the bam files in IGV. Alignment is good. For all samples 90% alignment rate was seen. I observed that in all samples Higher percentage of mapped reads were originating in Intronic regions. Followed by Exonic and intergenic regions.

I have seen a post here Reads mapped to exonic, intronic and intergenic regions where they say high intronic reads could be because of contamination. I googled about higher reads in intronic regions and found some papers Evaluation of two main RNA-seq approaches for gene quantification in clinical RNA sequencing: polyA+ selection versus rRNA depletion and some other links RIBO-DEPLETION IN RNA-SEQ – WHICH RIBOSOMAL RNA DEPLETION METHOD WORKS BEST? in which they said Greater intronic reads were with rRNA depletion protocol.

And even in this RNA-SEQ tutorial, it is mentioned that - A higher intronic mapping rate is expected for rRNA removal compared to polyA selection.

So, my question:

I am working with lncRNAs. So, I'm using the samples prepared with rRNA depletion protocol. Is this higher intronic rate is common in rRNA depleted dataset or do I have to check anything else to proceed further with these samples?

# Answer
It's not so much that you have "intronic contamination" or "genomic contamination", rather you're not selecting explicitly for full-length mature transcripts with rRNA depletion. That is the most common cause for higher intronic read rates. There's nothing you can do about this post-hoc, just continue along.

BTW, many lncRNA's are polyadenylated, so you'll keep them with poly-A selection.
> 1 votes

# Answer
I can't possibly see how intron contamination is linked to removal of rRNA depletion.

The only reason it would appear to have increased the number of contaminants is because post-rRNA removal the proporation of intron contaminants has increased against the total remaining RNA content. However, the actual total number of contaminants remains exactly the same pre- and post-rRNA removal. By the same token the RNA of interest lncRNAs will also have proportionally increased, so you get a better depth of its predominance and diversity.

Thats just life and perhaps just filter this bioinformatically.
> 0 votes

---
Tags: rna-seq, reads, rna-alignment, ribosomal, strandedness
---


# Question
Title: How does sequence alignment help with protein folding prediction?
I've been seeing using sequence alignment to help with prediction of protein folding, but am unable to understand why. Google returns other research papers for these search terms, hence I'm hoping the stackoverflow community can help thank you.

# Answer
Protein sequences from the same family form a sequence alignment. Some positions in the alignment will be conserved, i.e. the same in all sequences. These are often associated with function, e.g. a catalytic residue at position 50.

However there is a more subtle signal in sequence alignments called *covariation*. Two positions, say residue 30 and residue 60, may appear to change together across the alignment. We infer that they are close in 3D space in the structure:

Initially statistical techniques (DCA, PSICOV etc.) were used to separate direct from indirect (chained) couplings and extract predicted contacts from sequence alignments. Now deep neural networks have improved performance and even allowed the prediction of residue-residue distances and orientations.

These constraints on the structure allow model building, for example using the CNS software. A good set of constraints effectively guarantees the correct fold. The depth and quality of the sequence alignment is the major factor in the accuracy of the predicted constraints.

This approach to protein structure prediction is the current state of the art. Use of covariation from sequence alignments is the only approach that has reliably given accurate models for larger proteins without templates. See this paper and this paper for the best current methods, and tune into the CASP conference this December to see if anyone can improve on them.
> 4 votes

---
Tags: sequence-alignment, protein-structure
---


# Question
Title: Pathway enrichment analysis with or without custom background
I am doing a pathway enrichment analysis of a sets of proteins. Specific cell types (cardiospheres) were used to perform the experiment and for the quantification of proteins was used mass spectrometry. I am using gProfiler to perform pathway analysis, should I use custom background of all of those proteins identified from the experiment? or should I used without custom background?

Thank you!

UPDATE: Also, I am doing an enrichment analysis of the interactome of the set of proteins of interest. In that case should I used a custom background for the enrichment analysis or not?

Thanks!!

# Answer
> 1 votes
I would defninitely go with a custom background which should be all the proteins that were actually measured in your experiment since you should only take into account what was measured.

---
Tags: pathway, go-enrichment
---


# Question
Title: Motivation behind the neighbor-joining distance matrix recomputation
In each iteration of the neighbour-joining method of phylogenetic trees construction, after joining the nearest neighbours the (additive) distance matrix $D\_{m \times m}$ is recomputed as

$$Q\_{ij} = (m-2)d\_{ij} - \Bigg( \sum\_{k \neq i} d\_{ik} \Bigg) - \Bigg( \sum\_{k \neq j} d\_{jk} \Bigg)$$

As is stated in numerous sourced that I found, it prevents the joining of non-direct neighbours, but am looking for a more detailed explanation and motivation.

Sources: Wikipedia

# Answer
> 2 votes
Let's look at the components of the equation one by one first since it is easier to describe each part on its own. We can do this because they are added together in computing $Q\_{ij}$. Component1, $(m-2)d\_{ij}$ is saying that we want a quantity which would be the total penalization, if every distance for $m$ sequences was $d\_{ij}$ excluding the 2 distances between themselves (it's a direct proportionality so that the $Q\_{ij} \propto (m-2)d\_{ij}$). Component2, $\left(\sum\_{k\neq i}d\_{ik}\right)$, says that we look at the sum of the distances sequence $i$ has against all other sequences excluding itself (it's an inverse proportionality so that the $Q\_{ij} \propto \left(\sum\_{k\neq i}d\_{ik}\right)^{-1}$), and $\left(\sum\_{k\neq j}d\_{jk}\right)$ as component3 is analogous.

*non-direct neighbours* means that on the phylogenetic tree, 2 extant sequences (sequences in the alignment) have another sequence set that can link them together. If two sequences are definitely related to each other then on the phylogenetic tree they will produce a bifurcation and the smallest distance between themselves as being *directly* related, $(m-2)d\_{ij} \< \left(\sum\_{k\neq i}d\_{ik} + \sum\_{k\neq j}d\_{jk}\right)$, eg. the $d\_{ij}=1$ and the rest are all greater than 1 gives us a low number. If $f\_{max\_{i,j}}(d\_{ij})$ is used we'd be using the maximum distance for component1 so the overall measure would be positive and as large as possible.

In summary it provides a measure for placing the sequences together with the closest distance. You may say that this is 'overkill' as we can just search the $m^2/2-m$ set of distances and choose the minimum. In the goal for producing the tree, things get more interesting. Let's say look at 3 sequences $i,j,k$ of a larger set $m$, and see that sequences together have the same distance between them $d\_{ij}=d\_{ik}$ in respect to $i$. $d\_{ij}=d\_{j\*}$ would mean that $j$ is equally distance to all other sequences as it is to $i$ and let's say $d\_{kz}=max \land d\_{ky,y\neq z}=d\_{ij}$ (the sequence $k$ is equidistance to all sequences except to a particular sequence $z$ and can happen in some region in the alignment). Which do we join with $i$? According to the formula we choose $k$ over $j$ because $k$ is more distant from other sequence(s) as it would provide a more negative number (lower values cause of the 2nd and 3rd components although component1 is the same in both cases). Effectively we want the minimum distance to $i$ and simultaneously the maximal distance from other sequences (penalize sequences which could be closely related to the other sequences helps produce an overall better tree).

*Anders Gorm Pedersen (of DTU)* in his coursera videos (also on YT) provides some numerical examples to see how the numbers are produced on an example tree.

# Answer
> 2 votes
NJ as with all clustering methods there isn't an explanation, It is a hierarchical clustering method but there are many algorithms to cluster a distance matrix. NJ is closely associated with the concept behind UPGMA (unweighted pair group mean algorithm)

The key difference with UPGMA is that after joining two taxa by NJ the reconstructed genetic distance between the node and one taxa versus the other taxa is not assumed to be equal, i.e. the rate of mutation is heterogeneous. All other clustering algorithms I know assume a "global molecular clock" if that makes any sense, the mutation is constant against time. NJ removes that assumption and thats why we like it.

---
Tags: phylogenetics, sequence-alignment, phylogeny, sequence-analysis, matrix
---


# Question
Title: How to get rows with similar values in two different columns using command line?
I have following example data from a large file.

```
cont100 1128 1125
cont1005 3642 3642
cont1006 103 19
cont1037 3146 3146
cont104 895 890
cont1056 934 934
cont1059 1750 1750
cont1072 2577 2577
cont1078 43 42
cont1081 727 2

```

I need to get only the rows with similar values in column 2 and column 3. The output should be like this.

```
cont1005 3642 3642
cont1037 3146 3146
cont1056 934 934
cont1059 1750 1750
cont1072 2577 2577

```

Also I need to get the inverse of that. That means I need to get the rows with values that are not similar in the columns. The output should be like this.

```
cont100 1128 1125
cont1006 103 19
cont104 895 890
cont1078 43 42
cont1081 727 2

```

# Answer
> 8 votes
I think it is easy to do with awk, e.g. see here.

For your example, with your data in `file`:

```
%  awk '$2 != $3 {print $0}' file
# outputs
cont100 1128 1125
cont1006 103 19
cont104 895 890
cont1078 43 42
cont1081 727 2

%  awk '$2 == $3 {print $0}' file
# outputs:
cont1005 3642 3642
cont1037 3146 3146
cont1056 934 934
cont1059 1750 1750
cont1072 2577 2577

```

Note that where you say "similar" I am interpreting that to mean "exactly the same", because that is the case in your example.

# Answer
> 2 votes
The Perl solution is similar to the awk solution from *MaximilianPress*:

Print lines with identical values in columns 2 and 3:

```
perl -lane 'print if $F[1] == $F[2];' in.txt

# Prints:
cont1005 3642 3642
cont1037 3146 3146
cont1056 934 934
cont1059 1750 1750
cont1072 2577 2577

```

Print lines with different values in columns 2 and 3:

```
perl -lane 'print if $F[1] != $F[2];' in.txt

# Prints:
cont100 1128 1125
cont1006 103 19
cont104 895 890
cont1078 43 42
cont1081 727 2

```

The Perl one-liner uses these command line flags:
`-e` : Tells Perl to look for code in-line, instead of in a file.
`-n` : Loop over the input one line at a time, assigning it to `$_` by default.
`-l` : Strip the input line separator (`"\n"` on \*NIX by default) before executing the code in-line, and append it when printing.
`-a` : Split `$_` into array `@F` on whitespace or on the regex specified in `-F` option. Arrays in Perl are 0-indexed.

**SEE ALSO:**
`perldoc perlrun`: how to execute the Perl interpreter: command line switches

---
Tags: bash
---