text:
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

You should post this as an issue on their GitHub repository. You might get a quicker response that way.

What is the output of repophlan_microbes.log?

thanks, it is now added in question

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?

Thanks for your time and help. Actually we don't have active FTP connection on server which is required for the download. That can be a issue?

@bioinfonext Yes, that probably explains it. This would fall under the second bullet point of my answer: it keeps on trying to use an FTP connection which your server doesn't support, which fails, and then it waits and retries, which fails, and then it waits and retries... which would explain the long run-time as noted in the answer and their documentation.

@bioinfonext Also: it would be very helpful if you posted the log in answering your question, as we have requested a couple of times.

Hi Maximilian, I have added the head and tail of log file in the question above!

Hi Maximilian, I have run the script again on iMac and now I am getting an error in the log file and it keep showing it again and again! not sure what can be solution for this? is there any other way by which I can download all the bacterial and archeal genome sequences directly from NCBI?

@bioinfonext Looking at it, it seems that they are errors related to downloading specific genomes using FTP. I have updated further.


text:
How can I run the fdnadist command from the EMBOSS package?
My operating system is Mac OS X.
I would like to compute DNA distances between different sequences. For this, the dnadist package from the Phylip software would be good but hard to automate. EMBOSS seems to have an fdnadist command for this. I have installed emboss version 6.6.0 from the bioconda package but if I type fdnadist, the command is not found. How can I use dnadist or fdnadist on the command line?

To find emboss packages type wossname, then hit enter twice.  You will see all the programs.  There is one called distmat.  Then do tfm distmat, to get the fkn manual for distmat.

I found fdnadist in embassy-phylip package https://anaconda.org/bioconda/embassy-phylip

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
Virus genome: protein inside another protein
The genome sequence for the GenBank virus JX869059.2 (a human betacoronavirus) contains, among others, two proteins:

N protein (range 28566..29807) and
Orf8b (range 28762..29100).

The interval for Orf8b is contained in the interval for N.
How are these proteins translated? Is the orf8 protein created with post-translational processing after N is translated, or does the ribosome translate the two proteins separately?

https://www.ncbi.nlm.nih.gov/nuccore/JX869059.2?report=GenBank report all proteins really build from virus or they are reported from page https://www.ncbi.nlm.nih.gov/orffinder/ query "JX869059.2"

because the 28762-28566 is not divisible by 3, the title is wrong, it is not one protein inside one other, and only ribosome can compose orf8 protein in few i know...so ribosome has to compose both (N and orf8 protein, if they exist)

Viruses are really compressed in terms of information —we've all become familiar with the single long peptide that gets chopped into functional protein by the protease of SAR-COV-2. In HIV, translational slippage gets used to make a rarer peptide (as opposed to alternative splicing say). However, whereas enzyme are generally big peptides that inhibit others (such as the MERS homologue of orf8) need only 5-10 amino acids to do the trick (like a piece of paper to fix a wobbly table).

Note that the protein sequence of ORF8b is not contained within the protein sequence of protein N. That means that they are translated from different frames (though they are in the same orientation). Indeed, if we just subtract the start coordinates we do not find a multiple of 3.
Thus, they cannot be results of the same translational event. This means that the smaller protein could result from translational slippage (as suggested by Matteo Ferla), "leaky scanning" of the translation machines for start codons, or it could simply result from a different ribosome assembly site, recruited via the viral mechanism of an internal ribosome entry site, which allows a ribosome to be recruited in a mRNA cap-independent fashion (though this may be less likely).
For some mechanistic details of how one such ORF is translated, see this paper.
It looks like such things are sometimes called "subgenomic mRNAs". For a paper on such nested genes in coronavirus, including I believe Orf8b, see here.

Essentially @MaximilianPress hit the nail on the head. sgRNA is a major feature of SARS-CoV-2.
This is where the 5' UTR is spliced at numerous locations in the genome, such that although it is a single stranded positive sense RNA genome it behaves like a segmented viral genome because it is translated as a series of sub-genomic RNA molecules. So influenza is a classic example of a segmented genome. SARS-CoV-2 is more slightly complicated because it is using frameshifts to achieve to achieve the different sgRNAs.
However it isn't a particularly strong example of frameshifts because virus such as HBV (hep B) are pass masters, where two different proteins are encoded within a single locus.
My personal view is that slippage isn't operating here, unlike say HIV, but that it is a function of the sgRNA construction, thus N protein occurs on one sgeRNA and ORF8b occurs on a second sgRNA and sg signal is within the N protein.
If SARS-CoV-2 is your thing then stuff like this is a big player in interpreting bioinformatic output.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
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?

It is usually helpful to do a little research with a search engine or to give more information about what exactly you need. See e.g. https://www.schrodinger.com/kb/936

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.


text:
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:

I don't think there is a principled way of choosing such cutoffs. The first and last set of filters look much better than the second set.

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.

"gringer says the Q-Q plot conform to the beta distribution" - no, I didn't say this. The p-values should conform to the linear distribution if the association is not omnigenic, but the test statistic [which is not a p-value] will not fit a linear distribution. Most GWAS programs produce beta values as their association statistic, which is why I mentioned beta.

Thanks Michael for your opinion. I am coming from the ML field were evaluation of a model with training data is more clear to me. I found it very hard to estimate GWAS filtered by "feeling" and looking at plots (also because I want to do it automatically). I am new to GWAS and all I need is a simple GWAS. I was unsure about using PCA or LMM  or I could go for a simple linear analysis and whether this looks kind of "normal".Thanks for your feedback.

Welcome to 'traditional' statistics @snowflake ... Q-Q is seen as a general approach. If it is approximately linear you good to go. In this instance that isn't quite the case. HW has meaning here (biological meaning). PCA doesn't have much meaning, its okay for data-mining, whilst GLM or GLMM is the sort of analysis where transformation is used to linearise a Q-Q. Peronally I'd stick with what you have and look at the biologial understanding. Just my opinion.

@gringer post corrected

@Michael Ok thanks, well that's a bit strange for me but I have to live with it :) But I appreciate your explanation and opinion a lot!

Hi @snowflake keep in mind the central rule of ML is that ANY equation can be replicated in ML, particularly deep learning. I'm not sure a population genecist will agree :-) but thats the theory. So you could use the classification to train an ML model and ditch the 'traditional approach' in future.

PCA are used in GWAS; they're commonly used for covariate adjustment in the statistical models.

@Michael That's the main goal and my interest. However I first need a running GWAS for comparison. I think looking at the entire SNP pattern instead of single SNPs  could be more informative (similar to  image patterns or finger prints). Let's see...

@snowflake okay I see, personally I think your HW results are fine given the caveats. However, if population genomic structure through HW isn't your thing (I would say it is a clear biological model) then first I'd do PCA and then depending on your purpose assess the populations e.g. via Fst. Personally I would just do a phylogenetic tree given these are SNPs. From a phylogenetics point of view the problem with PCA is that it is biologically difficult to interpret why the populations or different regions of the genome are separating. The thing which supercedes this dogma is ML partic. DNN

A phylogenetic tree depends on linkage disequilibrium, which is the predominant basis of your data. So you could use your HW analysis to justify a phylogeny. Obviously I don't know the species.

Follow up question https://bioinformatics.stackexchange.com/questions/14500/gwas-phenotype-data-format-and-preprocessing

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:

https://doi.org/10.1006/tpbi.2001.1543
https://doi.org/10.1038/s41467-018-05257-7

Thanks a lot for this detailed explanation, I will have a look on the beta statistics as well. plink is reporting beta and std, you are right.


text:
Error while splitting PDBQT file
I am trying to split to a pdbqt file containing 178 compounds into individual pdbqt files. I have been using the Autodock Vina_split.exe file to split the multiple entry file however I keep running into the following error. What can I do to prevent the error.

(A) Empty file? Does the file open in PyMol? Also, you can always split it there with split_states.
(B) Is the file really called output.pdbqt or is it output.pdbqt.txt and Windows is hiding the extension?

There could be many reasons. Maybe lack of permissions or the space in the path. Can you open the file with a different program? Are you on Linux or Windows etc.


text:
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

The question is fine. Are you fixed on using R or would you use external packages to do this analysis? Do you use Macs?

Fixed on Linux, checked out Beast but had problems installing. Semifixed on using R but open for suggestions (maybe py(thon)curious)

Okay thanks ... really this is an R coding issue. Biopython is a whole new script.

In the end I will say "whatever works", but I feel like I'm looking for my glasses that are on my nose

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.

I went away for some days. You may be right, but now I'm having problems using both ETE and Mesquite. I will try some more and hopefully figure something out. Manuals and cookbooks look promising though.

It is really easy to do what you want to do. ETE3 needs Biopython to import the tree. I could post code but you should be able to do this via R easily. As a last resort simply export your R tree and open it in Figtree and just click the option to 'view branch lengths'. I presume you want to automate the process so at a guess 'ape' will have the option. Mesquite is for identifying each individual alignment position that contributes to a given branch, my understanding is you just want number of sites not alignment positions.

Great. ETE had bugs, Mesquite did not like my xml format. Figtree smiled, patted me on my back, took my nexus file from R and showed me the branch length.

Great! I am very happy for you! Congratulations to FigTree. ETE3 is just complicated, I like it because I'm use to it (as a Python coder I'm a bit protective). I do use ETE3 to print bootstraps, so it should just another option therein, but admittedly unless your performing a ''big data' operation its not worth it.


text:
summarising read group information from a .bam file
I have merged together 2 different .bam files in order to simulate sample contamination. So the reads can come from one of two samples, as shown by the read group info:
@RG ID:0    PL:ILLUMINA SM:LP4100018-DNA_C11_Proband    PU:HGY3WDSXX:1:none
@RG ID:1    PL:ILLUMINA SM:LP4100018-DNA_C11_Proband    PU:HGY3WDSXX:2:none
@RG ID:2    PL:ILLUMINA SM:LP4100018-DNA_C11_Proband    PU:HGY3WDSXX:3:none
@RG ID:3    PL:ILLUMINA SM:LP4100018-DNA_C11_Proband    PU:HGY3WDSXX:4:none
@RG ID:0-11EFC00B   PL:ILLUMINA SM:LP4100018-DNA_E11_Proband    PU:HGY3WDSXX:1:none
@RG ID:1-B8A1099    PL:ILLUMINA SM:LP4100018-DNA_E11_Proband    PU:HGY3WDSXX:2:none
@RG ID:2-330086F    PL:ILLUMINA SM:LP4100018-DNA_E11_Proband    PU:HGY3WDSXX:3:none
@RG ID:3-7681F092   PL:ILLUMINA SM:LP4100018-DNA_E11_Proband    PU:HGY3WDSXX:4:none

I'd like to check that the correct proportion of read groups originate from each sample.
Currently I am using:
samtools view example.bam | rev | cut -f 1 | rev > output.txt

, but this is not very elegant and only works because the RG field is last in the .bam.
Is there a quick way to tabulate the number of reads groups with different IDs? E.g. produce an output like:
ID:0 1000
ID:1 2000
ID:2 3000
...

A solution in samtools would be ideal, along the lines of the output produced in samtools stats.

Does this answer your question: https://bioinformatics.stackexchange.com/q/10559/6251 ?

Also, try sambamba.

@TimurShtatland, yes, using samtools view -R read_groups.txt example.bam | wc -l will work. thanks.

Would you like to self-answer the question then? So it's easier for others to find the solution.


text:
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?

What format are your genotypes in at the moment? There is likely some tool someone has written to convert it to the correct format.

Just table of 'A/C' bases

'Table' isn't a format. What are the rows? What are the columns? Is it genotypes or alleles - i.e. is there one or two rows per marker? the more info you give, the easier it will be for someone to help you. Could you add a sample of your genotype format to the main answer so others can see it.

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.


text:
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?

Are you asking why your pull down experiment was not 100.00000% perfect?  You should read the literature or ask your co-workers who have done pulldowns before what % is typical for samples like yours.

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.


text:
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

Just ask your cluster admin / manager if it's possible to run notebooks via the queueing system / workload manager, if it is not, you won't get it running anyhow.

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/ .

thanks I am running it as python script after loading jupyter and python module but still shows an error? I am going to update the script and the error.

%matplotlib inline is a jupyter notebook command. Try to comment this line.

Again error. it is updated now.

Those are again some jupyter notebook writings to use bash commands within a notebook. You should try with subprocess module instead or build a bash script that does those steps and then call your python script.
I think you should think by yourself now and googling until you are really stuck, this is a good training to write real python script ;)

Hi @bioinfonext please don't forget to upvote this post and if you feel it answered your question to mark as 'resolved' if resolved.

You look to be running Jupyter notebook interactively and if run from PC/laptop/Mac via a jupyter-notebook command it should connect interactively. Pointing this onto your server I'm not certain however. Normally jupyter- notebook is run on one of the big cloud HPC platforms such as AWS or GCS.

I think now it is a bit confusing because the post has been edited a lot. But I think @bioinfonext wrote a jupyter notebook to plot his/here data and then wanted to do the same task but by running it on a cluster using SLURM job scheduling manager.

should I use below line of code instead "a = !unzip $qza" ;  import subprocess
a = subprocess.Popen(["unzip'', qza])

unzip is your command, qza is a python variable, if you look for "give variable to subprocess call" you find: https://stackoverflow.com/questions/4795190/passing-variables-to-a-subprocess-call

Thank you so much. Now I am trying to run this on server terminal after loading the modules and with some changes, but still getting error? I have updated the code and error.

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.

script and error updated after importing numpy.

@bioinfonext commented below

Yes, thanks, this is all very confusing, R is much easier to adopt for potting.But I could not able to find any good R code for rarefaction plot that can plot it based on metadata column.

@bioinfonext answered above

@bioinfonext answered above


text:
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.

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.


text:
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.

You could read about differential co-expression. But that would be feasible depending on the number of samples you have for each condition...

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.

Thanks for the input and sorry for the late response. I've been looking into association analyses for modules/eigengenes and I'm thinking of running LRT analyses using GLMs with and without the interaction term of interest. Would you suggest I run separate analyses using each module's eigengene values individually (vs phenotypes) then correcting for multiple hypothesis, or can you recommend a more elegant way to test all modules at once? We would also like to run 3 separate analyses for each of 20 modules based on tissue sample subsets, so running 60 different analyses is not ideal. Thanks.

For module eigengenes, a linear model such as ones implemented in limma should be appropriate, no reason to run GLMs. With a linear model, you can simply test the significance of the interaction coefficient which is output by all linear model functions that I know (as long as the full and reduced model only differ by one coefficient). As far as I am aware, the only way to run such association analysis is fitting the models to the eigengenes one by one (perhaps with some empirical Bayes modifications as in limma) and then correcting for multiple comparisons.


text:
Using dssp after chain extraction
I have a list of PDB IDs with realtive chains IDs that have to be extracted, and then run on dssp.
For the single chain extraction I tried several methods, such as:
import pymol2
for entry in list:
    with pymol2.PyMOL() as pymol:
        pymol.cmd.fetch(entry)
        pymol.cmd.save(entry + '.cif')

Where, for example, if I want to extract chain AC for 4YBB, entry = 4YBBAC.
This works well in saving the cif files (pymol correctly visualizes them), but makes them not suitable for mkdssp. Indeed, if i run
mkdssp 4YBBAC

I get
DSSP could not be created due to an error:
Is this an mmCIF file?

I tried also using Biopython:
from Bio.PDB import MMCIFParser
from Bio.PDB.mmcifio import MMCIFIO
parser = MMCIFParser()
io = MMCIFIO()
structure = parser.get_structure("4YBB", "4YBB.cif")
model = structure[0]
chain = model['AC']
io.set_structure(chain)
io.save("4YBBAC.cif")

But again, when running mkdssp, I get:
DSSP could not be created due to an error:
bad lexical cast: source type value could not be interpreted as target

On a closer inspection, both methods drop the columns _atom_site.pdbx_formal_charge  _atom_site.auth_comp_id _atom_site.auth_atom_id from the file, while retaining all the others.
This is the first line of the coordinates from the biopython/pymol file:
ATOM   1    N N   . GLY A ? 1   ? -125.713 46.378 -19.108 1.0 77.15  2   AC 1

While this is the original file:
ATOM   34684  N  N     . GLY C   3  1    ? -125.713 46.378   -19.108  1.00 77.15  ?  2    GLY AC N     1

I think that renaming the AC chain to a single-character code would in principle resolve the issue, but it's more of a workaround.
mkdssp does support multicharacter chain ids. I've proved it first it by cutting and pasting the chain into the file and then with an ad-hoc parser. I think the issue is with how biopython and pymol extract the chains.
Is this a bug? Am I using the tools incorrectly? Can my purpose be achieved with other (possibly hassle free) methods?
MORE ON THIS TOPIC ON Biopython GitHub :
https://github.com/biopython/biopython/issues/3439

Could you post the PDBx/mmCIF loop_ before the coordinate lines of the pre and post structures/chains to better show what dropping the column mean

Biopython 1.81 development https://biopython.org/docs/dev/api/Bio.PDB.DSSP.html has a module for DSSP

to me looks like https://github.com/biopython/biopython/blob/master/Bio/PDB/MMCIFParser.py is missing 'pdbx_formal_charge' while  https://github.com/biopython/biopython/blob/master/Bio/PDB/mmcifio.py has it.

If you are you using PyMOL, why not use its command dss?
PyMOL is annoying for removing most entries bar for ATOM/HETATM data from mmCIF and PDB filenames. So the SS data would not be saved, but it's easy to append to the file (I have done it here for example as NGL.js requires SHEET and HELIX)

@saidem why do you think ; "I think that renaming the AC chain to a single-character code would in principle resolve the issue, but it's more of a workaround." , have you tried it ?

not sure if its relevant didi you checked laast verision of parser ?? https://biopython.org/docs/1.80/api/Bio.PDB.MMCIFParser.html it has : auth_chains - True by default. If true, use the author chain IDs. If false, use the re-assigned mmCIF chain IDs. AND auth_residues - True by default. If true, use the author residue numbering. If false, use the mmCIF “label” residue numbering, which has no insertion codes, and strictly increments residue numbers. NOTE: Non-polymers such as water don’t have a “label” residue number, and will be skipped.


text:
Truncated Chemical Names
[Note: this question has also been posted on KNIME]
Good Morning,
I am currently researching treatment options for the Corona Virus from the CAS COVID-19 Antiviral Candidate Compounds Dataset available on the www.CAS.org. One issue we have run into is that the Chemical Names of the Antiviral Compounds in the Dataset are truncated as depicted in the example below. We will be needing the entire Chemical Names for Downstream Analysis. Is there anyway we can obtain the complete/ untruncated Chemical Names of the Antiviral Compounds?
Example of an Antiviral Agent with a truncated Chemical Name
Acetamide 2-amino-N-(1-tricyclo[3.3.1.137]dec-1-ylethyl)-
  Mrv1921 09212018412D

 17 19  0  0  0  0            999 V2000
    0.0000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.0000    0.0000    0.0000 N   0  0  0  0  0  0  0  0  0  0  0  0
    0.0000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.0000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.0000    0.0000    0.0000 N   0  0  0  0  0  0  0  0  0  0  0  0
    0.0000    0.0000    0.0000 O   0  0  0  0  0  0  0  0  0  0  0  0
    0.0000    0.0000    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.8517   -1.6585    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.1163   -2.0322    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.5657   -1.2419    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    1.3012   -1.6156    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    2.0368   -1.2419    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    1.5873   -2.0322    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    2.0368   -0.4945    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    1.3012   -0.1208    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.8517   -0.9111    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.5657   -0.4945    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
  1  2  1  0  0  0  0
  2  3  1  0  0  0  0
  3  4  1  0  0  0  0
  4  5  1  0  0  0  0
  3  6  2  0  0  0  0
  1  7  1  0  0  0  0
  1  8  1  0  0  0  0
  8  9  1  0  0  0  0
  9 10  1  0  0  0  0
 10 11  1  0  0  0  0
 11 12  1  0  0  0  0
 12 13  1  0  0  0  0
  8 13  1  0  0  0  0
 15 17  1  0  0  0  0
 14 15  1  0  0  0  0
 15 16  1  0  0  0  0
  8 16  1  0  0  0  0
 12 14  1  0  0  0  0
 10 17  1  0  0  0  0
M  END

Thank you.

You should be using SMILES or Inkeys not human readable names regardless of whether they are IUPAC valid.


text:
install bowtie2 from sources cannot find -ltbb
I am trying to install bowtie2 aligner from sources without the root access.
Bowtie2 needs tbb package to be installed and it is recommended to install oneTBB.
# install one TBB
mkdir -p SHOME/soft
cd $HOME/soft
git clone https://github.com/oneapi-src/oneTBB.git
cd oneTBB
make

Then I add C_INCLUDE_PATH to .bashrc:
echo "export C_INCLUDE_PATH=$HOME/soft/oneTBB/include" >> ~./bashrc

and continue to install bowtie2:
# install bowtie2
mkdir -p $HOME/soft/bowtie/source
cd $HOME/soft/bowtie/source
wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.4.1/bowtie2-2.4.1-source.zip
unzip bowtie2-2.4.1-source.zip -d ..
cd ../bowtie2-2.4.1
make

TBB library is recognized but in the end I get this error:
/usr/bin/ld: cannot find -ltbb
/usr/bin/ld: cannot find -ltbbmalloc_proxy
collect2: error: ld returned 1 exit status

ltbb and ltbballoc_proxy are not in TBB but libtbb and libtbballoc_proxy instead.
I also tryed to add variables to .bashrc:
TBB_LIB_PATH=$HOME/soft/oneTBB
TBB_INCLUDE_PATH=$HOME/soft/oneTBB/include
TBB_RELEASE_DIR=$HOME/soft/oneTBB/build/home/suvar/apps/oneTBB/build/linux_intel64_gcc_cc9.3.0_libc2.31_kernel5.4.0_release

but seems it does not help.
I need help to complete the installation of bowtie2 on our server.

Did you remember to source your bashrc before compiling?

sure: source ~./bashrc

Do you really need to compile Bowtie2? There are compiled binaries for Linux and MaxOsX available:
https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.4.1/bowtie2-2.4.1-linux-x86_64.zip
https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.4.1/bowtie2-2.4.1-macos-x86_64.zip

If you are using Ubuntu you can install tbb with sudo apt install libtbb-dev, in MacOS you can install it with Homebrew brew install tbb.

You need to give the path to the source files for the compiler and the full path to the compiled libraries (all parent directories included) for the linker.  You also need to use either CPATH or CPLUS_INCLUDE_PATH instead of C_INCLUDE_PATH to point to the oneTBB source, since bowtie2 is C++ and C_INCLUDE_PATH will only have an effect for C (see Environment Variables Affecting GCC in the GCC documentation for more info.)
For example like this, inside bowtie2-2.4.1, without making any lasting changes to the shell's setup in this case:
$ ONETBB=$HOME/soft/oneTBB
$ export LIBRARY_PATH=$ONETBB/build/linux_intel64_gcc_cc5.4.0_libc2.23_kernel4.4.0_release
$ export LD_LIBRARY_PATH=$LIBRARY_PATH
$ export CPATH=$ONETBB/include
$ make

After it's built you'll still need the library path to oneTBB defined (via LD_LIBRARY_PATH) so bowtie2 can find it at run-time.
But if you don't have a specific reason to compile from source I'd recommend avoiding it.  You could go with what h.mon suggested and use an existing binary, or use a package manger to install it like zorbax suggested.  Conda has a package for bowtie2 in the bioconda channel, up through and beyond version 2.4.1, and you can use that as a package manager without root access.


text:
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)

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)

Hi Timur, Thanks, I have added the results from sambamba as well, I think it gives the accurate result. Could you please respond to query related to sambamba result?

@bioinfonext Edited the answer.

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)


text:
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

I would guess that the kmer counting will be much faster with tools like KMC. Is there any particular reason you want your own counter in python?

I know what you mean, but I am learning to write code and I just want to use python to do that! There are many tools out there like khmer, jellyfish, kat etc... Thank you

I guess that's fair. The question is related to this one about kmer counting in python. Although it does not specificly talk about multiprocessing, I thought it's still relevant.

counting the mers it is not the problem. My problem is  just the multiprocesing step thats giving me a hard time

This code help me to get the answer I was looking for.
r = pool.apply_async(f, a)
print(r.get())


text:
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

They want you to find genes whose response to the drug is different based on the genotype? If so, you can add an interaction term to your DESeq2 regression formula ~ condition + genotype + condition:genotype. Refer to the DESeq2 documentation and help("results", package="DESeq2") for more information on multifactor design and interaction terms.

Okay yeah that makes sense for how to obtain differences in the response, but then how might I statistically compare these differences in response for a group of genes (I.e a gene set)? Would you just do something simple like a t test / Wilcox test?

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.


text:
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.

Hi @JayLeeman better to post the code. Simple answer is to translate from each start point and assess length of the exon against the number of stop codons, i.e. maximising the length of the exon. There are 6 possible reading frames, 3 forward, 3 reverse. If you have an idea of the structure of the protein, e.g. surface antigen you can assess GPI anchors and hydrophobic alpha helices as the end of the gene

There is no reason at all for an exon to have either a valid stop or start codon.  If introns are in the sequence the whole enterprise is a waste of time.

@swbarnes2 This actually clears a lot out of me, thanks! Just clarifying: If I'm reviewing a DNA sequence, it still has introns. If there are introns, how can I locate the ORF's without being able to splice the introns "manually"?

You cannot programicly know where introns are.  You can only make good guesses.  Or, you have to look up what introns people have empirically determined in your sequence.

I'll look into that, sounds extremely interesting. Thank you!

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.

Yes I am doing this fun, and looking at orf finders results at the same time to compare and figure out how it does what it does. I'm still having some trouble understanding: Given a single frame (out of the six) - I may encounter a frame with 100 start codons and 300 stop codons, and rest coding for amino acids. Can we actually locate the ORF's in a case like this? If so, these are the "rules" I'm looking for. Thanks for the reply!

If you want to know how ORF finders work, why don't you...look at how ORF finders work?

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).


text:
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.

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.

@haci Thanks for such a wonderful explanation. I have just added the clustree plot. Its dont look that great; but would you please provide your comment? Thanks

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.


text:
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?

Usually removing genes that are not expressed in atleast n samples might help. E.g If you have  3 replicates per condition, keeping genes that are expressed in atleast 2 samples of one of the conditions.

DESeq2 doesn't have any built-in way to remove batch effects. Did you use limma::removeBatchEffect to do that, or are you referring to modelling the batch effect in the design formula, e.q. ~batch+condition? The only reason to filter genes before running DESeq2 is to speed up the process. The package already performs independant filtering when you run the deseq() function unless you explicitly turn that off. Can you share your data/pipeline here? This gene should not show up with that much of a fold change, so there is something else going on.

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.


text:
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!!

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.


text:
Finding the relevant annotations for biomarkers
I am looking to find a way to map a list of biomarkers of interest, say genes or bacterial functions, to existing annotation databases. Usually this is done in reverse, as such that based on metagenomics profiling you can see which pathways or genes are present. However, I want to be able to determine which genes or pathways to extract beforehand i.e. as I think it will be quite difficult to condense pathways or genes of interest from annotation data.
Example:
Biomarker gene for better sleep - maps to annotation database at ...
Biomarker gene for glucose metabolism - maps to annotation database at..
.. which in turn allows to directly quantify those biomarkers.
Thanks

Hi @alexander and welcome. How do you propose to 'map', are you looking for a specific software tool/website or else custom code? This question is best described as metabolomics.The question needs more clarity.

Thanks for your comment. For example, I have a metagenomic dataset with annotated genes, based on this, I want to quantify presence of pathways or related genes of interest. However, finding the relevant genes within the annotated file seems to be labour-intensive. As such, I would want to see which of my genes and pathways of interest are present in the annotated file(s). As such, my question is how to 'find' my genes of interest within the annotated files or used database.


text:
Convert MeSH term to drugbank ID
I have set of annotated term of drugs extracted from PubTator those drugs are annotated in form of MeSH terms. I would like to map them into drugbank ID. I do search the the MeSH term i find the reference to the drug and link to the drugbank. But not any database as such. How do i map the MeSH term to the drugbank ID.?
My small dataset.
Terms         MESH
1        ara-C MESH:D003561
2 mitoxantrone MESH:D008942
3    etoposide MESH:D005047
4  fludarabine MESH:C024352
5   gemtuzumab MESH:C406061
6   ozogamicin            -

As an example of one of the drug I would give this drug "gemtuzumab"
DB00056 Gemtuzumab ozogamicin   Myeloid cell surface antigen CD33   target

So i want to map the above drugs into drugbank ID and add the information such as the class which is here target or it can be a transporter or enzyme etc.
The issue is there is no common ID between the drug MEsH and the drug-bank ID.
Any suggestion or help would be highly appreciated if anyone can cite me some resources.

PubTrator description from the website
PubTator Central (PTC) is a Web-based system providing automatic annotations of biomedical concepts such as genes and mutations in PubMed abstracts and PMC full-text articles.


text:
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!

What does "circle" mean in this context?  Do you mean some number of PCR cycles?

I am sorry, I meant cycles (I edited my question)!

what is dds??? is it possible to ask 1 question

Ok, I have edited my question, only one question..I am sorry..

dds is the name name of the DESeq2 object in pretty much every vignette and tutorial.

You can include it in the design matrix, but you need to be careful. If you use it as a continuous variable, then you are trying to find a linear relationship between counts and the cycle number. This may not be the case and in same cases you actually introduce more noise into the model

I think you are trying to make a model more complicated than it should be

Ok, so you recommend against..if I use it as a factor? Thank you!

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.

Thank you a lot for your reply!The problem is that PCA does not look good and I have observed that the samples with 15 cycles tend to cluster. I have tried to use removeBatchEffect from limma (I made a design matrix with cycle as a factor).then the pca got better. In the downstream analysis I used the new design matrix with cycle like in the question, and then the differentially expressed genes got dramatically reduced...(from 750 to 100).Any thoughts?

It is possible to have some real DE genes even if your PCA looks discouraging.  100 is not a lot, but maybe that's realistic if your condition doesn't split nicely along your first few PCs. You might play around with making cycles a factor vs a number.  But if it's right that the number of cycles is really driving some difference, or representative of some real difference, accounting for it in the design will attenuate those differences.

Thank you so much for your reply, it really helps. One more question, I have observed that when I use vsd instead of rlog for PCA pot, PCA plot looks better. However, my samples are less than 30 (in the tutorial, rld is recommended especially for samples less than 30). Any ideas? Thank you in advance.

I think the author recommends vst for everything now.


text:
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!

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/


text:
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.

Where did you see using sequence alignment to help with prediction of protein folding? Perhaps I just too far from the topic to understand the question, but I find hard to understand where the question comes from and I wonder if there is a one correct non-opinion-based 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.

Verified top answer @pohje please mark as "resolved". The technical theory is called "convolution" regarding direct and indirect couplings - you can't just use raw covariancce because you get artefacts, but simplicity you can think of it as just covariance.