text:
Do we design primers for single stranded DNA with same criteria as we do for double stranded DNA?
I am trying to design primers for a single stranded DNA virus.
But I am having a dilemma as whether we design the primers same way as double stranded DNA?
do we have to design forward and reverse primers both or just a single primer forward or reverse works?

If you've e.g. a parvovirus you design primers as normal, i.e. must be a pair of primers as for double stranded DNA. Albeit there's only one strand present after the first round of PCR there both strands will be present and without the reverse primer the PCR can't work.
If you want evidence of this, you've had a SARS-CoV-2 diagnosis? Well thats a forward and reverse primer PCR starting with ssRNA converted to cDNA, i.e. single copy and then a pair of primers PCR the resulting locus.

Even long before covid, starting from single-stranded DNA is a standard PCR technique, for example gene cloning from extracted cDNA. The principle of the chain reaction does not change, you always need two primers for amplification. Else, you always amplify from the original template over and over again, but this is not exponential and you'd produce an excess of forward (or reverse) strand product which has no use.

So we will just remove both the denaturation steps from PCR conditions, and the rest of the steps remains same, is it right?

No, everything is the same. Exact same. Get local supervision for lab work, this is a bioinformatics community.


text:
Remote NCBI's Blast Perl API: maximising hits returned via command line
Question and Background I'm using NCBI's Blast Perl API to interrogate Genbank for example blastp. I don't get the throughput using their point and click web stuff.
I tried Biopython's Bio.Blast import NCBIWWW but its under access restrictions from the host, notably attempting to use it during the day (works sporadically). The API however gets through and runs as fast as a web search ... that's good.
Here's an example of the command I'm running:
CMD=Put&PROGRAM=blastp&DATABASE=nr&HITLIST_SIZE=2500&ALIGNMENTS=2500&QUERY=%3Emytest%0AMSILANKDTR%0A

The above query gives default Blast output. What I'm trying to do is increase the number of hit recovered to the maximum (5000). The options for this are HITLIST_SIZE, ALIGNMENTS, DESCRIPTIONS. In the above example I'm attempting a hit_size of 2500 ... but the same output is returned as for default. More generally  I can't get any of the "number of hits/descriptions/alignments returned" options to work, the default is returned each time.
Has anybody solved this?
Typical options for the API are given below
Typical options
Parameter   Definition  Type    CMD Allowed values
QUERY       Search query    String  Put *   Accession, GI, or FASTA.
DATABASE    BLAST database  String  Put *   Database from appendix 2 or one uploaded to blastdb_custom (see appendix 4)
PROGRAM     BLAST program   String  Put *   One of blastn, blastp, blastx, tblastn, tblastx. To enable megablast, use PROGRAM=blastn&MEGABLAST=on.
FILTER      Low complexity filtering    String  Put F to disable. T or L to enable. Prepend “m” for mask at lookup (e.g., mL)
FORMAT_TYPE Report type String  Put, Get    HTML, Text, XML, XML2, JSON2, or Tabular. HTML is the default.
EXPECT      Expect value    Double  Put Number greater than zero.
HITLIST_SIZE    Number of databases sequences to keep   Integer Put,Get Integer greater than zero.
DESCRIPTIONS    Number of descriptions to print (applies to HTML and Text)  Integer Put,Get Integer greater than zero.
ALIGNMENTS      Number of alignments to print (applies to HTML and Text)    Integer Put,Get Integer greater than zero.

Just to feedback on the answer and comment

Using -remote by locally installing blast as recommended by @Scott. This is a valid answer so definitely upvoted and recommended its further upvoted. I had tried this but not mentioned it in the question (really sorry), but it times out due to traffic issues. The API overcomes traffic issues. Again, @Scott's is a valid solution (considering I'd not explained I'd tried it).

More details ...

The Perl API is here https://blast.ncbi.nlm.nih.gov/docs/web_blast.pl
Line 95 below ...
$args = "CMD=Put&PROGRAM=$program&DATABASE=$database&QUERY=" .
$encoded_query;

This is changed to,
$args = "CMD=Put&PROGRAM=$program&DATABASE=$database&HITLIST_SIZE=$maxhits&QUERY=" . $encoded_query;

Also need to add an additional line at line 73
$maxhits = shift;

From the command line,
perl ./web-blast.pl blastp nr 2500 ./example.fa

(Note, The variable order is switched in the program).
example.fa - e.g.
>rleguminosarum
WGFDGSSTQQAEGRSSDCVLKPVAIYPDPARTNGALVMCEVMMPDGVTPHASNARATILDDEDAWFGFEQ
EYFFYQNGRPLGFPEQGYPAPQGPYYTGVGYSNVGDVAREIVEEHLDLCLAAGINHEGINAEVAKGQWEF
QIFGKGSKKAADQIWMARYLLQRLTEKYGIDI

I'm proficient in Perl5, i.e. there's no error in changes to the code in the API, albeit my expertise is object oriented and high performance Python.

Limit on the request ... I'm testing that now.
First results eeek @terdon might be right ... :-(

Can you give us a minimal example of the actual code that we can test on our local machines, please? I wonder if it could be a limit in individual post requests and if you might need to set the downloading of the results in some kind of loop to get N at a time.

While I haven't used your exact method above, using the BLAST command line tools return the desired number of target sequences
blastp  -query query.fa \
        -task blastp-short \
        -db nr \
        -max_target_seqs 5000 \
        -outfmt 6 \
        -out query.out \
        -remote

query.fa
>mytest
MSILANKDTR

The resulting query.out file does contain 5000 unique hits.

I'll mark this as "accepted" because this is a good solution given the original question I asked (I tried this approach, but didn't state it in the question). The full explanation is given below and the above answer should be read within the context of my explanation for future reference. Note unless the other Mods over-rule my "accept" please keep this comment here to provide future readers with the full solution and thanks to @scot

Right here's the response from NCBI,

Please send your future blast related questions to blast-help@ncbi.nlm.nih.gov directly. The gb-admin alias is for submitters with questions related to their sequence deposition/submission.

For your issue, you should look at the following two parameters:
DESCRIPTIONS  Number of descriptions to print (applies to HTML and Text)  Integer Put,Get Integer greater than zero.
ALIGNMENTS    Number of alignments to print (applies to HTML and Text)    Integer Put,Get Integer greater than zero.

Do please note that asking for this level of hits saved/returned could be problematic for the search. They could make the search exceeds the cpu time allocation and lead to error.

I tried this ... it was the same as "HITLIST_SIZE".
I figured the answer ...
Why its happening? The API uses a default parameter for the *maximum number of hits. It's imposing a limit on the number of hits returned, but this wasn't stated in the documentation of the parameters.
Whats the limit? It's >70 but <= 100. I suspect the limit is 100 hits. Thus anything above this its set at 100.
Summary The reason why the API is working and other approaches simply time out is NCBI is imposing an upper limit on the number of hits for the API and the Blast servers are not permitting requests above this threshold at peak time. Thus because other approaches are not imposing the restrictions on the NCBI server, the API gets through, whilst the other methods (blastp directly and Biopython) are timing out.
NCBI are recommending the API is because it is automating a 100 hit cut-off and the servers are deprioritising limits beyond this
I have filed a note of advice to NCBI because it should be clear in the documentation there's an upper limit on the HITLIST_SIZE parameter. It took me at least a full half day to figure this out (which isn't cool for limited documentation).

Proof of concept (untested) Thus, using @Scott's solution
blastp  -query query.fa \
        -task blastp-short \
        -db nr \
        -max_target_seqs 100 \
        -outfmt 6 \
        -out query.out \
        -remote

This will work with the -max_target_seqs 100 and thats all the API is doing.

source: Stack Exchange

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

sort: votes


text:
How to find amino acid sequence of a given protein
Is there a way to look up the amino acid sequence of a given protein? For example, what amino acids are used to produce Amylase? I spent hours googling but to no avail. If that's impossible, how can I find the sequence of A, C, G, U of the mRNA that instructs the ribosome to produce Amylase?

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.

Please [edit] your question and tell us what protein you are looking for. Which amylase? From which species? You can find anything you want on sites like https://www.ncbi.nlm.nih.gov/ or https://ensembl.org among others. Depending on the species of interest, you can also find more focused resources, so tell us what you're after (you mention both protein and mRNA, both are possible, but which one do you want?).

Usually you go and search in a database. NCBI and Uniprot (specialized on proteins) are ways to go here.
These two offer tons of other related information, links to other databases, advanced functionalities like similarity search based on sequence (e.g. blastp). etc.
To narrow your search, you would need to specify which species you are interested. Here a quick result for amylase in humans, using NCBI and Uniprot. The nucleotide sequence is linked in both cases.

I have no knowledge of biology, thus I had no idea a same protein name (Amylase in this case) will have a different structure for each species (thus vague OP). This question arose from a Python programming homework. Do you happen to know what these letters mkfflllfti gfcwaqyspn tqqgrtsivh... mean from the origin of this link https://www.ncbi.nlm.nih.gov/protein/AAA51724.1?

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


text:
Is there a method/web server to generate Position-specific scoring matrix (PSSM) for nucleotide sequences (DNA/RNA)?
Just like, PSI-BLAST is avaialble to generate a PSSM matrix for peptides, is there a web server where we can generate PSSM for DNA/RNA? I want to find homologous sequences to a given RNA sequence.

Please [edit] your question and explain what your final objective is here. How would you use this matrix? To do what? The choice of tool and approach will depend on what you are actually trying to do.

Not a web server, but a method is available in biopython to do produce PSSMs. I've used it a number of times and it works well as far as I can tell.

thanks, but I believe the biopython method caters to DNA sequences only. But that might help as well.

@terdon sure I have edited.

Is this RNA translated? If so, it would make more sense to search using tools like tblastn.  Does it have a specific structure? If so, you can use tools that look for similar structures. Do you have multiple homologs of it? Aligned? Or are you starting from a single sequences? What does the RNA do? Depending on what it does, you would focus on different features when looking for homologs.


text:
Using average of sequence similarities when delimiting genera - how to deal with outliers?
Sometimes I am trying to see where a (bacterial) genus ends and where another one begins based on 16S rRNA phylogenies. Of course, the length and support of the branch matter a lot, but this question in particular is about sequence similarity.
As far as I understand it, in general, sequences within the same genus are expected to have roughly >94-95% sequence similarity at the 16S rRNA level. But that leaves out a lot of variation. For example, in doodle A below, there may be <90% similarity between sequence 'a' and sequence 'e' ; however, there is always >95% similarity between, say, a and b, b and c, c and d, d and e...So all would belong to the same genus.

Previously I had considered just comparing every sequence with every other sequence in my group of interest, get an average and a standard deviation, and see if the average > 94-95%. However, it occured to me that if I have a tree shaped like that in doodle B (a couple of closely-related sequences from one genus and one very distantly-related sequence from another genus), I could get the very same average (and only slight changes in the standard deviation).
It feels like there should be a simple mathematical property to distinguish between A and B based on similarity values, and/or point out outliers like sequence 'a' in B, but I do not know what that would be. Help would be appreciated!
Mock dataset A:

sequence1
sequence2
sequence similarity (%)

a
b
99

a
c
93

a
d
92

a
e
95

b
c
99

b
d
98

b
e
97

c
d
95

c
e
94

d
e
88

Mock dataset B:

sequence1
sequence2
sequence similarity (%)

a
b
90

a
c
90

a
d
90

a
e
90

b
c
99

b
d
99

b
e
99

c
d
99

c
e
97

d
e
97

Wouldn't the average for all comparisons to "a" in dataset B be below 95? That would indicate that a is not the same genus as the others, following your definition of >94-95 similarity.

@Niklas Thank you. My issue is that this is simple enough for one sequence, but would get more difficult to see from the table if I add more (one hundred? five hundred? a thousand?). Of course, I can still spot the outlier easily enough in the tree...most times...but I feel that adds a degree of subjectivity that I am not quite content with. I would prefer something more rigorous.

@M__ I am not quite sure how to answer. I guess either? Basically, for me A and B correspond to a group of sequences (a-e), that could be seen as taxa, and can be easily represented as trees. So more or less interchangeable. I can visually spot outliers in the trees, but I hoped for something more objective, and working with the taxa at the level of sequence similarity felt like an option. Maybe it's not the only option?

@Laura I am not really sure if that is so complicated, basically you would just need to compute the average % identity across all comparisons that contain one of your taxa and see if that is below 95%? You could also plot it to make it more visual. It depends a lot on what you expect (is this a frequent or rare thing, will most species always be from one single genus?).. I think to "look at the tree" using some metric, total branch lengths (sum of all branches) could be informative for your problem, have you tried using that?

The basic question is where do genera start and stop within phylogenetics and the 16S 5% rule? This is not a trivial question within molecular taxonomy, especially for bacteria.

The key is to define the tree structure for 16S
Define the tree structure for other genes

Background and tree stats The difference between Tree A and Tree B is referred to as 'balance'. Tree A is unbalanced, whilst Tree B is balanced and that is a formal description. I would personally use the term Tree A was "heavily paraphyletic", or "extensive paraphyly" - at a technical level that is accepted, I have used frequently. Equally I'd describe Tree B as a "singular monophyly". The statistical basis for tree balance is called the Colless value.
"Colless statistic". It's fairly easy to calculate. The reason it is not often used is because the tree needs to be phylogenetically robust (below). Not much importance is given to it because if a branch isn't robust thats a headache for this statistic (it can't tell the difference between robust and non-robust).
Difference between branches The statistical difference between taxa A and taxa B is measured by bootstrapping (I'm sure you know this), where >= 80% is bifurcation is robust.
Colless statistic might be used, but its fairly obvious, bootstrapping is needed and iqtree is best package for this operation.

The second more difficult part of the question is the delineation of a genus. The 5% threshold of 16S has always been seen as a guide. I have never seen it rigidly enforced and it creates a conflict with phylogenetic theory as your example points out. Telling what is and ain't a genus ... you are right that a key is the tree structure.

Tree structure and bootstrapping ...

A declared monophyly must be supported by a bootstrap > 80%
A declared paraphyly ... each paraphyletic branch must be supported by a bootstrap > 80%

16S versus other genes and tree structure
The tree structure needs confirming particularly as the bootstrap values 16S will often be low. The starting point is to consider the traditional MLST loci (I forget what they are but stuff like atpD, etc ...), build trees from them and :

Don't combine 16S with protein genes (they evolve at different rates), although it could be done in a partitioned model
MLST protein genes and core-house keeping genes can be combined without partitioning the data (there's a test, but it'll basically work).

Take the MLST genes repeat the tree and assess the bootstraps. Basically it will produce the same tree, but: a) confirms the tree structure; b) can provide stronger bootstrap values (which confirm the tree).

5% stuff
All the above needs doing before any 5% stuff is considered. It is my personal view that using 5% 16S as a mean value is sort of bending the numbers, particularly in a paraphyly where the divergence is e.g. between 2 to 10%.
Monophyly vs paraphyly This where the tree structure is needed. Basically, a robust monophyly would be regarded as a taxonomic singular unit. In this case if its close to the 5% threshold genus rule:

even if the pairwise distance between taxa within the monophyly is say 8% via 16S;
and this is supported by other protein loci;
and there's no suspicion this could represent a family.

In the above criteria I'd say its evidence of a genus. 10% is a bit big distance however.
If the tree structure is paraphyletic and the maximum different is say 10% between two taxa - one at the top of the paraphyly and the other as a most recent child ... then thats tricky. It could easily be two genera within the paraphyletic tree structure. For example, the vertebrates as paraphyletic, viz. amphibian, reptile, bird, edentates, marsupial, mammals and that represents 6 monophyletic Classes. Thus if that was my result I would leave open the possible interpretations: maybe this is a single genus; alternatively it could be two (or more) genera. The key to this is defining the nodes (internal branches) via bootstraps to be confident of the tree structure.
Generally a single robust paraphyletic tree structure, i.e. heavily imbalanced, where each paraphyletic branch is robustly supported would not be regarded as a single genus especially when 16S genetic differences is >> 5%. Taxonomically, this scenario would need further phylogenetic analysis (separate question) and classical bacterial taxonomists would want phenotypes before splitting a genus.

Reviewing the question against the answer - what I would do for the monophyly (diversity 10%) is split off the basal groups (at 10%) view robust bootstrap support. Then forward the in-group as the genus and the full monophyly as a putative genus. So that 10% is not the divergence unequivocally defining a single genus. So thats quite easy on the picture drawn. Thus the in-group is <8% because 10% is a bit high.
It is harder to do this with a paraphyletic tree structure and thats a separate question on the advanced techniques needed to investigate this.

I would like to note out that the terms paraphyletic and monophyletic do not refer to trees, but to clades.

Hello, all terms of the "-phyletic" group inherently describe clades, a clade can either be mono-, para-, or polyphyletic. Since a tree has no inherent clades, these terms only make sense in combination when adding clades (e.g. taxonomy) to the tree. Clades are not necessarily monophyletic (which is why they can also be para- or polyphyletic).

Clades versus monophylys. I'm not getting into that discussion because of its acrimonious history. I have used "paraphyletic clade" and singularly used "clade" as monophyly as a published definition (a paper reviewed by Science journal). The semantics are not rigid and phylogenetics is more concerned with statistical robustness than "one person says this, another says that".

I am only raising this because you refer to the trees as  "heavily paraphyletic", or "extensive paraphyly"

Thanks, personal semantics as mentioned ("imbalance" better)

This is quite advanced and it will take a bit of time to wrap my head around, but it's very informative. I guess I'd hoped that there is something simpler and more objective than what I'd already been doing (looking at the tree shape, bootstrap support etc) but I guess that even when adding extra numbers it all comes back to that...One thing I'm not convinced about is that a robust paraphyletic structure should be split if you can 'step' from one end to the other using branches with >97% similarity to each other...but that is indeed another question. Overall, thank you!

Hi @Laura, glad it helps, you can simply ask the other questions, e.g. 97% as a separate question. because it depends what you mean by 97%, if that is base tip to the "most recent" tip, or distance between the ancestral node and the other tips.


text:
Trimmomatic QC report shows drop in the reads and presence of overrepresented sequences
This question was also asked on Biostars
I am performing a de novo genome assembly using Illumina paired-end short reads, sequenced on a NovaSeq X by our collaborator at UCLA.
At present, I am in the stage of trimming the adapters.
Here, you can have a look at the basic statistics and information on the adapter content obtained from the Fast QC report, for R1.
FOR RAW READS

I used Trimmomatic for trimming the adapter.
The following is the Trimmomatic Settings
ILLUMINACLIP:~/adapters/TruSeq3-PE.fa:2:30:10 MINLEN:36
Here, you can see the basic statistics and adapter content of the Trimmed reads.

Here, the output was:
Both surviving: 566832403
Forward only surviving: 39244376
Reverse only surving: 0.00
Dropped reads: <1%
Here you can go through the full FastQC report for:
Raw Reads
Trimmed Reads
Now I have two questions:
Question 1
Can I go ahead with the assembly process, for there is zero-adapter-presence in the reads? Should I mind the loss of reads?
Question 2
I see that there are over-represented sequences, both in read 1 and read 2. I doubt if I can leave them be, or I should trim them too. I would like to know if the over-represented sequences need to be tossed, or if they can be left untouched and proceed with the assembly. If I need to remove these sequences, can Trimmomatic do it?
If I were to trim the over-represented sequences, can I use an adapter "fasta" file as follows?
The following are the over-represented sequences for R1.

The following are the over-represented sequences for R2.

Can I have my fasta file as follows, to remove the over-represented sequences:
>PrefixPE/1
-------sequence1--------------
-------sequence2--------------
-------sequence3--------------
-------sequence4--------------
>PrefixPE/2
-------sequence1--------------

java -jar trimmomatic-0.30.jar PE s_1_1_sequence.txt.gz s_1_2_sequence.txt.gz
lane1_forward_paired.fq.gz lane1_forward_unpaired.fq.gz lane1_reverse_paired.fq.gz
lane1_reverse_unpaired.fq.gz ILLUMINACLIP:overRepresentedSeqs.fa:2:30:10  MINLEN:36

NB: Sorry about the lengthy post. Since the file sizes are large, it takes over a day for the Trimmomatic to finish the task, so that it is not feasible to go with multiple rounds of trial and error check.

Given that it's a NovaSeq, these are almost certainly "no signal" reads as a result of the two-colour chemistry and should probably be excluded (or at least ignored).
A long polyG sequence is unlikely to affect mapping by much, because they're fairly rare in most genomes. Unfortunately you can't exclude them using Trimmomatic alone.
PolyG tail trimming can be done using fastp. I don't have any experience using fastp, so you'll need to read the documentation yourself.

I don't quite get how it is > 75%! For example, for R1, counting all the reads containing polyG, accounts for 0.76% of the total reads. Was there a miscalculation?

Oh, sorry. I saw percentage and thought proportion. I've deleted that bit from the answer. Thanks for noticing that.


text:
how is log2foldchange calculated in DESeq2
I am trying to calculate log2FoldChange in DESeq2 by hand because I am trying to reproduce my results and understand how DESeq2 works.
What I want to know is:

Are the counts first fitted to a negative binomial? If so, how can I get access to those values?
Regarding log2FoldChanges, I read that the coefficents are the estimates for the log2 foldchanges for each sample group. However, these estimates do not account for the large dispersion we observe with low read counts. To avoid this, the log2 fold changes calculated by the model need to be adjusted.
Why? Didn't we just fit the counts to a negative binomial, which should take into account the dispersion.
Finally, how are the log2FoldChanges calculated? It's not possible to figure this out using the raw code because most of the real calculations call C scripts.

Have you checked if the DESeq2 vignette answers your questions? https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html

were you able to figure this out? I am trying to do something similar where I have plot showing the differential expression results for a single gene from the DESeq2 output. I plotted the marplot of the log2FoldChange and errors bars with lfcSE. In the barplot I want to include the gene expression ratios that were used in estimating the statistics. Was wondering if it is possible? Somewhat related to your question, but the difference is I am interested in the ratios of the log-counts
Any thoughts or hints will be appreciated.


text:
Enrichment analysis of DMPs | new Ilumina EPIC V2 900K
Im analyzing methylation data from the new illumina's EPIC V2 array (900K).
I would like to perform enrichment analysis on DMPs. In the past, I have always used missMethyl but since there are a lot of differences between EPIC V1 and EPIC V2 arrays I'm unsure if missMethyl can be used (Im afraid results might be biased)
Does anyone know any enrichment tool already available for the EPIC V2 array or alternatives to missMethyl?
Thanks


text:
How to pack up a Python script in the form a Dockerfile
I have a simple Python script with this header
from flask import Flask
app = Flask(__name__)

@app.route('/')
import os
import pandas as pd
from statsmodels.formula.api import ols
import statsmodels.api as sm

I have installed GitHub Desktop, Docker Desktop
I have created a folder in that I have
(.venv) (base) angel@Angels-MacBook-Pro docker % ls
Dockerfile      app.py          compose.yaml        requirements.txt

In requirements.txt I have
(.venv) (base) angel@Angels-MacBook-Pro docker % cat requirements.txt
blinker==1.6.2
click==8.1.6
colorama==0.4.6
Flask==2.3.2
itsdangerous==2.1.2
Jinja2==3.1.2
MarkupSafe==2.1.3
Werkzeug==2.3.6
gunicorn==21.2.0

This is the commands I have run so far
docker init

Then this one returns error
(.venv) (base) angel@Angels-MacBook-Pro docker % docker compose up --build
[+] Building 3.9s (14/14) FINISHED                                                                                docker:desktop-linux
 => [server internal] load .dockerignore                                                                                0.1s
 => => transferring context: 680B                                                                                0.1s
 => [server internal] load build definition from Dockerfile                                                                          0.1s
 => => transferring dockerfile: 1.63kB                                                                                0.1s
 => [server] resolve image config for docker.io/docker/dockerfile:1                                                                  2.3s
 => [server auth] docker/dockerfile:pull token for registry-1.docker.io                                                              0.0s
 => CACHED [server] docker-image://docker.io/docker/dockerfile:1@sha256:ac85f380a63b13dfcefa89046420e1781752bab202122f8f50032edf31b  0.0s
 => [server internal] load metadata for docker.io/library/python:3.9.12-slim                                                         1.3s
 => [server auth] library/python:pull token for registry-1.docker.io                                                                 0.0s
 => [server base 1/5] FROM docker.io/library/python:3.9.12-slim@sha256:364ada4cb4f943b603d330cbb3d26f689c8b54de51cb12bad1202de7dba3  0.0s
 => [server internal] load build context                                                                                0.0s
 => => transferring context: 1.39kB                                                                                0.0s
 => CACHED [server base 2/5] WORKDIR /app                                                                                0.0s
 => CACHED [server base 3/5] RUN adduser     --disabled-password     --gecos ""     --home "/nonexistent"     --shell "/sbin/nologi  0.0s
 => CACHED [server base 4/5] RUN --mount=type=cache,target=/root/.cache/pip     --mount=type=bind,source=requirements.txt,target=re  0.0s
 => [server base 5/5] COPY . .                                                                                0.0s
 => [server] exporting to image                                                                                0.0s
 => => exporting layers                                                                                0.0s
 => => writing image sha256:c207ecf929eb3dc4f172db22e56a19e7b89892fcdc31da80501cb3524e137f9d                                         0.0s
 => => naming to docker.io/library/docker-server                                                                                0.0s
[+] Running 1/1
 ✔ Container docker-server-1  Recreated                                                                                0.1s
Attaching to docker-server-1
docker-server-1  | [2023-09-22 11:09:20 +0000] [8] [INFO] Starting gunicorn 21.2.0
docker-server-1  | [2023-09-22 11:09:20 +0000] [8] [INFO] Listening at: http://0.0.0.0:8000 (8)
docker-server-1  | [2023-09-22 11:09:20 +0000] [8] [INFO] Using worker: sync
docker-server-1  | [2023-09-22 11:09:20 +0000] [9] [INFO] Booting worker with pid: 9
docker-server-1  | [2023-09-22 11:09:20 +0000] [9] [ERROR] Exception in worker process
docker-server-1  | Traceback (most recent call last):
docker-server-1  |   File "/usr/local/lib/python3.9/site-packages/gunicorn/arbiter.py", line 609, in spawn_worker
docker-server-1  |     worker.init_process()
docker-server-1  |   File "/usr/local/lib/python3.9/site-packages/gunicorn/workers/base.py", line 134, in init_process
docker-server-1  |     self.load_wsgi()
docker-server-1  |   File "/usr/local/lib/python3.9/site-packages/gunicorn/workers/base.py", line 146, in load_wsgi
docker-server-1  |     self.wsgi = self.app.wsgi()
docker-server-1  |   File "/usr/local/lib/python3.9/site-packages/gunicorn/app/base.py", line 67, in wsgi
docker-server-1  |     self.callable = self.load()
docker-server-1  |   File "/usr/local/lib/python3.9/site-packages/gunicorn/app/wsgiapp.py", line 58, in load
docker-server-1  |     return self.load_wsgiapp()
docker-server-1  |   File "/usr/local/lib/python3.9/site-packages/gunicorn/app/wsgiapp.py", line 48, in load_wsgiapp
docker-server-1  |     return util.import_app(self.app_uri)
docker-server-1  |   File "/usr/local/lib/python3.9/site-packages/gunicorn/util.py", line 371, in import_app
docker-server-1  |     mod = importlib.import_module(module)
docker-server-1  |   File "/usr/local/lib/python3.9/importlib/__init__.py", line 122, in import_module
docker-server-1  |     raise TypeError(msg.format(name))
docker-server-1  | TypeError: the 'package' argument is required to perform a relative import for '.venv.lib.python3.9.site-packages.eventlet.wsgi'
docker-server-1  | [2023-09-22 11:09:20 +0000] [9] [INFO] Worker exiting (pid: 9)
docker-server-1  | [2023-09-22 11:09:20 +0000] [8] [ERROR] Worker (pid:9) exited with code 3
docker-server-1  | [2023-09-22 11:09:20 +0000] [8] [ERROR] Shutting down: Master
docker-server-1  | [2023-09-22 11:09:20 +0000] [8] [ERROR] Reason: Worker failed to boot.
docker-server-1 exited with code 3

yaml
# Comments are provided throughout this file to help you get started.
# If you need more help, visit the Docker compose reference guide at
# https://docs.docker.com/compose/compose-file/

# Here the instructions define your application as a service called "server".
# This service is built from the Dockerfile in the current directory.
# You can add other services your application may depend on here, such as a
# database or a cache. For examples, see the Awesome Compose repository:
# https://github.com/docker/awesome-compose
services:
  server:
    build:
      context: .
    ports:
      - 8000:8000

# The commented out section below is an example of how to define a PostgreSQL
# database that your application can use. `depends_on` tells Docker Compose to
# start the database before your application. The `db-data` volume persists the
# database data between container restarts. The `db-password` secret is used
# to set the database password. You must create `db/password.txt` and add
# a password of your choosing to it before running `docker compose up`.
#     depends_on:
#       db:
#         condition: service_healthy
#   db:
#     image: postgres
#     restart: always
#     user: postgres
#     secrets:
#       - db-password
#     volumes:
#       - db-data:/var/lib/postgresql/data
#     environment:
#       - POSTGRES_DB=example
#       - POSTGRES_PASSWORD_FILE=/run/secrets/db-password
#     expose:
#       - 5432
#     healthcheck:
#       test: [ "CMD", "pg_isready" ]
#       interval: 10s
#       timeout: 5s
#       retries: 5
# volumes:
#   db-data:
# secrets:
#   db-password:
#     file: db/password.txt

Dockerfile
FROM python:3.9
WORKDIR /docker
COPY . /docker
RUN pip install -r requirements.txt
CMD ["python", "./app.py"]

Add the contents of your compose.yaml and Dockerfile here as well. The latter can be deduced from your log but it'd be nice to view the former.

Thank you, I have added these to my main post

I can't tell what to do next, my experience ends here. Hopefully someone else can help. It might be worth consulting StackOverflow as this looks like a programming problem.

I’m voting to close this question because it belongs on StackOverflow

Not sure which tutorial you're following, but I think it's overly complicated. You do not need the compose file, and the Dockerfile has some issues:
First, specify what you want to copy into the container. That is requirements.txt and app.py. Then also, tell pip where the requirements file is. This works in terms of putting a container together. I am not testing the python script itself, this is on you.
$ cat Dockerfile
FROM python:3.9
WORKDIR /docker
COPY ["app.py", "requirements.txt" ,"/docker/"]
RUN pip install -r /docker/requirements.txt
CMD ["python", "./app.py"]

$ cat app.py
from flask import Flask
app = Flask(__name__)

@app.route('/')
import os
import pandas as pd
from statsmodels.formula.api import ols
import statsmodels.api as sm

$ cat requirements.txt
blinker==1.6.2
click==8.1.6
colorama==0.4.6
Flask==2.3.2
itsdangerous==2.1.2
Jinja2==3.1.2
MarkupSafe==2.1.3
Werkzeug==2.3.6
gunicorn==21.2.0

$ docker build -t atpoint/foo:1.0 .
[+] Building 5.8s (9/9) FINISHED
 => [internal] load build definition from Dockerfile                                                                                0.0s
 => => transferring dockerfile: 38B                                                                                0.0s
 => [internal] load .dockerignore                                                                                0.0s
 => => transferring context: 2B                                                                                0.0s
 => [internal] load metadata for docker.io/library/python:3.9                                                                                5.7s
 => [1/4] FROM docker.io/library/python:3.9@sha256:031b2a713079b2436e6c2102dbffbce109c9c2582ad0867c380e3c8d7fed477f                                                                                0.0s
 => [internal] load build context                                                                                0.0s
 => => transferring context: 64B                                                                                0.0s
 => CACHED [2/4] WORKDIR /docker                                                                                0.0s
 => CACHED [3/4] COPY [app.py, requirements.txt ,/docker/]                                                                                0.0s
 => CACHED [4/4] RUN pip install -r /docker/requirements.txt                                                                                0.0s
 => exporting to image                                                                                0.0s
 => => exporting layers                                                                                0.0s
 => => writing image sha256:6a6f1b08c8343608d2eacaa8463a1638491b504a3230ab00d45a7b32a41e1a43                                                                                0.0s
 => => naming to docker.io/atpoint/foo:1.0                                                                                0.0s

Use 'docker scan' to run Snyk tests against images to find vulnerabilities and learn how to fix them

Your container is now available as atpoint/foo:1.0 which reads as username (atpoint), name of container (foo) and a version number / tag (1.0), all of which is free to choose. This is the most basic way of packing a container.


text:
Is a classification tree appropriate method to use for my project?
I have about ~3000 genomes of organisms and presence (recorded as 1)/absence (recorded as 0) data of ~300 genes in them. Among 3000 genomes, we have experimental data on 20 organisms, which were represented by 2 different phenotype groups (a, b). According to a previous study, the 300 genes I mentioned are believed to be associated with phenotypes 'a' and 'b'.
To predict the phenotypes of 2980 genomes with no phenotype information, I did hierarchical clustering of 3000 genomes based on the presence/absence of the 300 genes. I cut the resulting dendrogram at a certain height that allowed distinction of clusters that included either group 'a' or 'b', but not both. This gave me 120 clusters, including clusters with no group 'a' or 'b' organisms, clusters with at least one group 'a' organisms, and clusters with at least one group 'b' organisms. Among these 120 clusters, there are 15 clusters that had either group 'a' or 'b' organisms.
I am questioning which collections of genes are responsible for causing phenotypes 'a' or 'b'. To address this, I am thinking about assigning different names for the 15 clusters that contain group 'a' or 'b' organisms, and make a classification tree using these 15 clusters. Does this sound like an appropriate application of classification tree approach? If not, could there be a more sound way to address this question?
Thanks!

This isn't a cool strategy because dendrograms are very poor representations of the relationships between genomes, for example they use UPGMA style metrics fail to account for mutation rate variation between genomes and prone to long-branch attraction artefacts.
I would recommend doing whole genome alignments using e.g. muscle5 -super option and then performing a formal phylogenetic tree using e.g. iqtree or even just MEGAX (distance methods) and using bootstrapping to identify the robustness of clades comprising a and b.
If two clades are supported by bootstraps > 80% then your hypothesis looks fine. At present your downstream analysis could be an artefact of the dendrogram phylogeny.


text:
efetch multiple records at once
If I want to efetch two records, I can do the following:
efetch -db nuccore -id AB610939,AB610940 -format fasta

I have hundreds of records to fetch. Can I put all the IDs in a file and efetch. If so, how?

Yes, efetch can be given a file with IDs. See efetch -h:
$ efetch -h | grep -m1 file
  -input         Read identifier(s) from file instead of stdin

So you can pass your file of IDs (e.g. ids) directly:
$ cat ids
AB610939
AB610940

$ efetch -db nuccore -input ids -format fasta
>AB610939.1 Synthetic construct RNA, spike-in microarray control RNA CRM 6204-a, clone: 6204-a-500-1
CTCGACTAGTTAATACGGTACAGGATAACCGATCGGCTTGCAACATAACGGCGTTAAGAATGCGGGAGTG
CAGTTTCCGATTCTCACATCAATCGCCAATAAGGCCTTGTCGCAATATAGACTCAACGGTTCTAGTAGCT
GATCGGTATTACGTGACGCAACCGATTAGACATGCACAATTCCTTGGTCGCTATACTACGGAAATCGTCA
GGTACTATAACCCGTCGCAGGCCTAATACGTGTCGTCACATCGCCAACCTATCGTCAGTCGGAAAGACGT
TGCTGTCTACCATCGAAACTATTTACCGCTCCGAGATTCACGAGTACGAACTCACGAGGAAGTTGCCCTA
TGTAAGGTATCACTCCAGGTACTGCGCCGATAGTACCAGGTGATCAAACGGTTGCAAGAAGGCCACGACG
TATCGGGCTCTTTAGACGTACGCTCGAGATTAAACGCGCACTGATTCACTTTAGCCCGGAATGTCTCGGT
GCGATGTAGA
>AB610940.1 Synthetic construct RNA, spike-in microarray control RNA CRM 6204-a, clone: CRM6204-a-500-2
AGACTAAATCTCGGCGTCGGTTCATACGCGCGATCGTTTGCTGTCAGGGCATACTCGAATCCGGACTCCG
ACAATTATAGGCCATCCTGAATAGCCGATCATGCGAGTCACGATAAGGCAGGCTCTGCGATATCCCGATA
TACTGGAGAAGCTGAATCCCACCTAGAGCGAACTGTCAGAGGATCGACCTCAGGCTCGCTATCATCATAA
CGGCGGACGACCTGTGTCACATTCCGAACGCTACGTGACGATATTATCTGTCGAAAGGCATAGAACGCCG
GTCAATATCCTGCGGCATTCTCTTTATCACCGGCTATAACTACTAGGTTCCGCAGATATAGACTGCGCAC
GGAACATGTAGATAGATCGAGTAGGGTAGCGATTTAACGACTCGACTTACAGACAGAGACGTAGAACGTC
AGACGAGTGGTATGCCCACCAGAGGCGATACAGGCTGTACCTGCGTAGCACTAGAGTCGTGCGTCATGCG
GACCCTATCT

If that were not possible, another option would be to put all your IDs as a comma separated string. If your IDs are in a file named ids, one per line, you can convert that to a comma-separated string like this:
$ tr '\n' , < ids | sed 's/,$/\n/'
AB610939,AB610940

That assumes GNU sed. If your sed doesn't handle \n, you can use this instead:
printf "$(tr '\n' , < ids | sed 's/,$//')" "$(cat ids)"

Armed with this, you can use command substitution to pass the comma-separated list to efetch:
efetch -db nuccore -id $(tr '\n' , < ids | sed 's/,$/\n/') -format fasta

If there are limits to te number of IDs you can fetch at once, or as may well be the case, if this results in a command that is too long, then you can use a loop as suggested in another answer.

However, be aware that while looping seems to work for efetch, for reasons that are not clear to me, the esearch developers have chosen to make it consume stdin. This means that the obvious, straightforward approach fails with esearch. See https://unix.stackexchange.com/q/682748/22222 for details.

If you're happy looping through each id in turn and each id in the single file (idfile.txt) is separate by a newline:
while read -r x; do
  efetch -db nuccore -id  "$x" -format fasta
done < idfile.txt > fetchout.fa

Speculation If the ids are not separated by a newline but a comma on a single line that might give the simultaneous efetch you're wanting. I strongly suspect that each id would need to be surround in quotation marks in such a file. Again this is just speculation and the above solution is preferred IMO.

Sorry, my bad. That issue seems to be specific to esearch, I cannot reproduce it with efetch. I deleted my comment.

@terdon world-first from the BioInfo 'nix tools police ;-)


text:
Linking GenBank records to biosamples (and vice versa) using edirect
Assume that I wish to find all complete human mitochondrial genome records on GenBank (or rather, NCBI nuccore) that also have an entry in NCBI's Biosample database.
MYQUERY="(mitochondrion[TITLE] OR mitochondrial[TITLE]) \
    AND complete genome[TITLE] \
    AND (human[TITLE] or homo sapiens[TITLE])"

I know that 61,719 such records exist on NCBI nuccore as of Sep-23-2023.
esearch -db nuccore -query "$MYQUERY" | \
    xtract -pattern ENTREZ_DIRECT -element Count

I also know that 452 such records exist on NCBI Biosample as of Sep-23-2023.
esearch -db biosample -query "$MYQUERY" | \
    xtract -pattern ENTREZ_DIRECT -element Count

However, these record sets either do not have an intersection or are not correctly linked (or I am using the edirect tools incorrectly), as 0 hits are identified.
esearch -db nuccore -query "$MYQUERY" | \
    elink -target biosample | \
    xtract -pattern ENTREZ_DIRECT -element Count

On the other hand, first esearching in the biosample database and then linking the hits to the nuccore database does seem to give 271 hits.
esearch -db biosample -query "$MYQUERY" | \
    elink -target nuccore | \
    xtract -pattern ENTREZ_DIRECT -element Count

What is the reason for this seemingly unpredictable behavior?

It's because the elink is messed up and not consistent. It's quite a lot of code to sort it out, but it's doable (I've done it). In an ideal world the elink would be the key across the relational database - it ain't always the case but its the best thing there is.

Not related to your issue, just a couple of points about good practices in shell scripts: i) avoid using CAPS for shell variable names. Since, by convention, global environment variables are capitalized, it is bad practice to also use caps for your own, local variables because this can lead to naming collisions and hard to find bugs. ii) There is no need for \  after a | to break a line. | is a list terminator, so it's fine as a line break character. See


text:
How come bigWigSummary shows coverage greater than 1 for some region of the genome?
I am running this command, using the UCSC utilities:
bigWigSummary -type=coverage ./galGal6.phastCons77way.bw chr4 25526708 25528708 1

to get phastCons77way coverage from 25526708 to 25528708 on chromosome 4 of the chicken genome.
The file I am using can be downloaded from here. I assumed "coverage" would be the ratio of loci with phastCons scores (so if a region of 2,000 bp has phastCons for every bp, the coverage should be 1). I think the phastCons was calculated vis multi-alignment, so it's not like each of the 77 species would have a score for that locus on the chicken genome, but rather they together has one score calculated for each bp when available.
However, I got 1.1345 as output. I am confused, 'cause I thought coverage cannot be more than 1, aka 100%.


text:
Would DESeq2 be appropriate or should I use another tool for differential gene expression?
Hay, this question is sort of an extension of previous post. We have sequenced extra cellular RNA from culture media. Since the number of cells were limited, we did not have enough material for biological or technical replicates, so each sample has only been sequenced once, and there are essentially three groups and each group has 2 subgroups, main groups are control, group A and group B. I think the table below should help -

group
subgroup
Number of unique samples

Control
Day 3
~30

Control
Day 5
~5

Treatment_1
Day 3
~30

Treatment_1
Day 5
~50

Treatment_2
Day 3
~40

Treatment_2
Day 5
~20

We want to compare the gene expression (Or rather abundance of RNA) between samples based on the Day subgroup.
Since the same sample / culture batch has not been sequenced with replicates, and there is variance in number of samples per group, I was wondering if DESeq2 would be a valid option, as the documentation does mention that without replicates, the output from DESeq2 should be considered only for exploratory purpose.
If DESeq2 is not appropriate, which other tool would be better?
EDIT - I came across NOISeq which I will try once I get any output from DESEq2, if anyone has any other tools / papers which can potentially help me out, do let me know!
Thanks!

If I am understanding your question right, you are after Day 3 vs Day 5. In that sense you should be able to test this with a formula like ~ group + number-of-samples + subgroup. This will enable comparing "days" while controlling for the number of samples and treatment. You might also consider adding an interaction term group:subgroup if you think treatment (or lack of it) has differential effects on day 3 and day 5.

We are after same day comparison between groups, sorry if I had not worded it correctly.

When you are lacking replicates, most statistical tools would either complain or even refuse to run the test you would like to perform. Even when a stats tool gives you an output without replicates, you will not be able to generalize your results anyway. In your case, a gene that you get from such an analysis might be the result of biology or just some technical artefact. I don't think your questions can be answered with the data that you have in hand.

I am aware of the caveats you've mentioned, I just want to do the best I can with the data we currently have. I came across NOISeq, but it doesn't seem to be in development any longer, but can work without replicates.

It's very simple. You can fiddle with these unreplicated data and get unreliable results. That takes time and in the end will be unsatisfying. Or you repeat this experiment twice for n=3 and then have a triplicate. That's what I'd do. These rank-based analysis and custom playing with unreplicated data just never yield anything substantial.

Again, I agree with you, but what I'm looking for now is a tool / methodology for rank-based analysis so that we can make the best of the data we already have. I'm looking at NOISeq for now, if you can suggest any tools, please let me know.

Definitely a fan of NOIseq :-) (long story) ... there's criticism it lacks power. What you are looking at is the significance of the observation against Type 1 error. The general criticism of non-parametrics is that they have low power, so few results. Personally, I'd use it and given how widely used DSeq2 is perform the same analysis via Dseq2. I'm not a fan of parametrics in this context, but in this field they dominate and DSeq2 is the best. Thus I'd use a parametric (DSeq2) and non-parametric NOIseq. The justification is you are using one of each. Personally, I do neither, I use machine learning to direct the thresholds of a parametric - but thats outside the public domain.
The time series is more complicated. Its essentially a 2-way problem: within a time point and between time points. Binary log difference is a good analysis between two time points in this situation.
I assume there was treatment between day 3 and day 5 for the treatment samples but not the control.
The issue of replicates. Personally what I would attempt to take a short cut with replicates: if Treatment 1 and 2 have a shared impact on the sample, then treating them as replicates and compare them against the control would be the first stage of analytics. This way if you go to the expense of producing the replicates described here you can be sure you've got a result before conducting further wet-lab work.
Clearly I don't know what Treatment 1 nor Treatment 2 are, but they could be drugs attacking the same metabolic pathway. Basically to do this you need a biological rationale for combining Treatments 1 and 2 to be considered replicates. I've attempted to re-describe this approach to biological analytics below the "-----"
The post-hoc analysis is more complex. The post-hoc analysis is a database operation assessing the relative binary log fold difference between time points for treatment versus control.

@haci likely represent the orthodox approach (e.g. at review level). What I would attempt is to provide biological justification that Treatment 1 and 2 were replicates. That depends on precise justification, i.e. this time series models the impact of generic interruption of this generic class of drug on pathway X, i.e. something that biochemically combines treatment 1 and 2.
The other analysis is a within time point analysis and if you did that between Treatment 1 and Treatment 2 for day 3 and day 5 separately (2 analyses), this would be used as justification if there ain't much different within a time point. In addition if treatment is post day 3, then Treatment 1, 2, and control are replicates, i.e. a triplicate (that's unequivocal). That provides a basis for moving forward.
However, The interaction @haci is describing is likely a formal model because this is a two-way problem. I rarely disagree but you've got to be careful here because its a time-series which aren't  independent events and most "interaction" based models require independence. Thus day 5 is a result of day 3, therefore they are not independent. I've encountered this as an editor and it can go very wrong. Summary analysis between time series - those are independent, analysis between data points within a time-series are not independent and interaction mixes those up.
Note
Time series experiments are always complex and performing wet-lab replicates would keep things simple (but obviously more expensive). So, if you've good analytics support or you're an analyst then you probably do it without further wet-lab, otherwise send it back to the wet-lab and then follow the recipe.

Thanks for your input, I wouldn't mind listening to the NOISeq long story (Once I'm done with my analysis, but you're welcome to PM). The time series does complicate things and inherently Day 3 will be different than Day 5. There are a good number of samples for each group, but lets see what the analysis show.


text:
Mapping List of Probes/Primers/Short Oligos to a Reference Fasta/q
Could you help me bioinformatics SE people.
So I have been duckduckgo-ing up a storm and even consulted the AI overlords and haven't come up with a solid way to do this, but it is so fundamental to bioinformatics that it has to already exists and I don't need to be re-inventing the wheel...
So imagine you have a fasta/q file with a "reference" sequence:
AATTTTAAAAAAACCCCAAAAAAAAAAAGGGGGAAACCACAAAA
And a list of short oligo sequences:
TTTT, CCCC, GGGG
How would be the best way to map (if that is the right terminology) the short oligos to the reference seq so that a configurable mismatch amount is allowed:
  TTTT       CCCC           GGGG    CCCC
AATTTTAAAAAAACCCCAAAAAAAAAAAGGGGAAAACCACAAAA

Blast, BWA and Bowtie2 from my understanding aren't ideal for this because they aren't built for really short sequences, something about them using seeded reference lookup?
Using Smith-Watterman or Needleman-Wunsch technically gets the local aligning right but the mismach allowance is based on score not X total allowed. Plus, iterating over hundreds of oligos and aligning against dozens of large references seems inefficient.
The closest thing I have found to what I'm trying to do is Sequence Manipulation Suites' Primer Map. This tool is almost perfect except a) it is a web tool, not a module I could import into python or incorporate into a CLI and b) does not appear to offer any tolerance for mismatches.
Am I just searching with the wrong vocabulary? This has to be an already solved bioinformatics process right?

A reference FASTQ is not a thing: A reference sequence is curated in some way, and a curated sequence does not come with quality scores.

Here is the source code of Primer Map You are not restricted to use it as a web tool only. Moreover you can modify it according to your needs.

@RamRS this is why I put reference in quotes lol. Basically the bare sequences are parsed out of a fasta/q file and are used as targets to align all the shorter oligos to.

@M__ I'm trying to determine coverage of a pool of primers against arbitrary target(s)

@haci Thanks! I cloned a copy and am working on dissecting it and figuring out if it has the ability to be modified for mismatching. Still feel like this is something that already exists and I just don't have the right search terms to find it...

@RPINerd I don't quite understand what you mean by "coverage of primers" but you might want to check out primer3. I used this in the past to design/test a large number of primers. You can customize and run via CLI.

@haci taking the example I gave in the OP, the "coverage" would be (4 bp * 4 oligos)/ 44 bp target sequence, so about 36%. And taking start stop positions into account I could also make some statistics on the distance between aligned sites

Both of @haci's solutions are very cool:

Primer Map
Primer3 CLI

I wasn't aware that Primer3 was available as a CLI (command line interface).
The solution I had in mind was in silico PCR one example being Kelpie and the associated paper is here, however I am not exactly sure about the application of the "primer mapping", but most "primer" applications concern PCR.


text:
Shell script to validate fastq issue
I have two GSM ID as test case where both of them are having data when I checked through sra explorer tool but when I try to fetch through a script only one of them returns output where the other one fails.
GSM ID
GSM3603268
GSM2683458
Output
For GSM3603268
$ cat GSM3603268.test.fastq
@GSM3603268.1 1 length=100
TTTTTTAAAGAAAACTTGAGCTTTTGGAGCATGGCAACCTAGCCTGCAGACACCGTATCCCCTTGTCCACTTCCCCCTGCAAACCATAAGTCCATTCCTA
+GSM3603268.1 1 length=100
=@@DDDDA++2CDABCEG@GEHHEHGA9EDFHGHD3::B??<DBFFBFB===)5;CA;(7=?;?CE;7)7;>C((583;<9>(,89@@@@CCCCDD:BC:
@GSM3603268.1 1 length=100
AGTTCCACGAGTTTTCTTTTTTTTAAGTGGTAGGAATGGACTTATGGGTTGCAGGGGGAAGTGGACAAGGGGCTACGGTGTCTGCAGGCTTGGTTGCCAT
+GSM3603268.1 1 length=100
:?=BD:D?C)<CDGH4+<CCCGHI)??D<D?00?D*?39C<CG4BFE28;4(.7AB############################################

$ cat GSM3603268.test.log
2023-09-26T21:18:46 fastq-dump.3.0.1 int: string unexpected while executing query within virtual file system module - multiple response SRR URLs for the same service 's3'
2023-09-26T21:18:46 fastq-dump.3.0.1 int: string unexpected while executing query within virtual file system module - multiple response SRR URLs for the same service 's3'
2023-09-26T21:18:46 fastq-dump.3.0.1 int: string unexpected while executing query within virtual file system module - multiple response SRR URLs for the same service 's3'
Read 1 spots for GSM3603268
Written 1 spots for GSM3603268

Output for GSM2683458
$ cat GSM2683458.test.fastq
    bash-5.1#

$ cat GSM2683458.test.log
2023-09-26T21:16:55 fastq-dump.3.0.1 int: string unexpected while executing query within virtual file system module - multiple response SRR URLs for the same service 's3'
2023-09-26T21:16:55 fastq-dump.3.0.1 int: string unexpected while executing query within virtual file system module - multiple response SRR URLs for the same service 's3'
2023-09-26T21:16:55 fastq-dump.3.0.1 int: string unexpected while executing query within virtual file system module - multiple response SRR URLs for the same service 's3'
2023-09-26T21:17:11 fastq-dump.3.0.1 warn: directory not found while opening manager within virtual file system module - can't open NC_000076.5 as a RefSeq or as a WGS
2023-09-26T21:17:11 fastq-dump.3.0.1 err: directory not found while opening manager within virtual file system module - failed GSM2683458

=============================================================
An error occurred during processing.
A report was generated into the file '/root/ncbi_error_report.txt'.
If the problem persists, you may consider sending the file
to 'sra-tools@ncbi.nlm.nih.gov' for assistance.
=============================================================

fastq-dump quit with error code 3

The script
#!/bin/sh

set -x
PS4='[\\d \\t] '

# Check parameter for error
check=0
# Print fastq-dump executable path
echo $(which fastq-dump)

# Function to download FastQ with retries
download_with_retry() {
    local id="$1"
local retries=5
for attempt in $(seq "$retries"); do
    # Download start of fastq
    fastq-dump $(get_ngc) -X 1 -Z --split-spot "$id" > "${id}.test.fastq" 2> "${id}.test.log"
    numLines=$(cat "${id}.test.fastq" | wc -l)
    if [ $numLines -gt 0 ]; then
            echo "${id} has data... OK"
        return 0
    else
        echo "${id} does not have data on attempt $attempt... Retrying in 5 seconds..."
        sleep 5
    fi
done
echo "${id} could not be downloaded after $retries attempts... ERROR"
    check=1
    return 1
}

# Loop through all parameters to check validity
for file in "$@"; do
cp "${file}" .
    # Extract filename for sampleID
    file_basename=$(basename "${file}")
    id="${file_basename%".id"}"
# Start validation with retries
echo "Checking ${id}..."
    download_with_retry "$id"
done

# Exit with error if some FastQs are not accessible
if [ $check -gt 0 ]; then
    echo "ERROR: One or more samples have inaccessible FastQs.. exiting"
    exit 1
fi

What exactly is going wrong in the case of the other ID? No matter how many retries I give, it fails.
Any suggestion or help would be really appreciated.
The error log:
$ cat ncbi_error_report.txt
<Report>
 <Run>
  <Date>
   <Start value="Tue Sep 26 2023 9:16:55 PM"/>
   <End value="Tue Sep 26 2023 9:17:11 PM"/>
  </Date>
  <Home name="HOME" value="/root"/>
  <Cwd>/tmp</Cwd>
  <CommandLine argc="6">
   <Arg index="0" value="fastq-dump"/>
   <Arg index="1" value="-X"/>
   <Arg index="2" value="1"/>
   <Arg index="3" value="-Z"/>
   <Arg index="4" value="--split-spot"/>
   <Arg index="5" value="GSM2683458"/>
  </CommandLine>
  <Result rc="RC(rcVFS,rcMgr,rcOpening,rcDirectory,rcNotFound)"/>
  <User admin="true"/>
 </Run>
 <Configuration>
  <Files count="2">
   <File name="/etc/ncbi/settings.kfg"/>
   <File name="/root/.ncbi/user-settings.mkfg"/>
  </Files>
  <refseq state="not found"/>
  <krypto state="pwfile: not found"/>
  <sra>
   <quality_type>raw_scores</quality_type>
  </sra>
  <Config>
  <ConfigurationFiles>
/etc/ncbi/settings.kfg
/root/.ncbi/user-settings.mkfg
    </ConfigurationFiles>
    <APPNAME>"fastq-dump"</APPNAME>
    <APPPATH>"/tmp/"</APPPATH>
    <BUILD>"RELEASE"</BUILD>
    <HOME>"/root"</HOME>
    <HOST></HOST>
    <LIBS>
      <GUID>"119c217a-7b81-47e8-91d6-56d19c8c9f15"</GUID>
      <IMAGE_GUID>"119c217a-7b81-47e8-91d6-62229c64ee59"</IMAGE_GUID>
    </LIBS>
    <NCBI_HOME>"/root/.ncbi"</NCBI_HOME>
    <NCBI_SETTINGS>"/root/.ncbi/user-settings.mkfg"</NCBI_SETTINGS>
    <OS>"linux"</OS>
    <PWD>"/tmp"</PWD>
    <USER></USER>
    <VDB_CONFIG></VDB_CONFIG>
    <VDB_ROOT></VDB_ROOT>
    <kfg>
      <arch>
        <bits>"64"</bits>
        <name>"56d19c8c9f15"</name>
      </arch>
      <dir>"/root/.ncbi"</dir>
      <name>"user-settings.mkfg"</name>
    </kfg>
    <libs>
      <cloud>
        <report_instance_identity>"true"</report_instance_identity>
      </cloud>
    </libs>
    <repository>
      <user>
        <ad>
          <public>
            <apps>
              <file>
                <volumes>
                  <flat></flat>
                  <flatAd>"."</flatAd>
                </volumes>
              </file>
              <refseq>
                <volumes>
                  <refseqAd>"."</refseqAd>
                </volumes>
              </refseq>
              <sra>
                <volumes>
                  <sraAd>"."</sraAd>
                </volumes>
              </sra>
              <sraPileup>
                <volumes>
                  <ad>"."</ad>
                </volumes>
              </sraPileup>
              <sraRealign>
                <volumes>
                  <ad>"."</ad>
                </volumes>
              </sraRealign>
              <wgs>
                <volumes>
                  <wgsAd>"."</wgsAd>
                </volumes>
              </wgs>
            </apps>
            <root>"."</root>
          </public>
        </ad>
      </user>
    </repository>
    <sra>
      <quality_type>"raw_scores"</quality_type>
    </sra>
    <vdb>
      <lib>
        <paths>
          <kfg>"/usr/local/bin"</kfg>
        </paths>
      </lib>
    </vdb>
  </Config>
  <RemoteAccess available="false"/>
  <CurrentProtectedRepository found="false"/>
 </Configuration>
 <Object path="https://sra-pub-run-odp.s3.amazonaws.com/sra/SRR5755657/SRR5755657" type="database" fs_type="unexpected">
  <Dependencies>
   <List count="22" missing="22">
    <Dependency index="0" seq_id="NC_000067.5" local="false" path=""/>
    <Dependency index="1" seq_id="NC_000068.6" local="false" path=""/>
    <Dependency index="2" seq_id="NC_000069.5" local="false" path=""/>
    <Dependency index="3" seq_id="NC_000070.5" local="false" path=""/>
    <Dependency index="4" seq_id="NC_000071.5" local="false" path=""/>
    <Dependency index="5" seq_id="NC_000072.5" local="false" path=""/>
    <Dependency index="6" seq_id="NC_000073.5" local="false" path=""/>
    <Dependency index="7" seq_id="NC_000074.5" local="false" path=""/>
    <Dependency index="8" seq_id="NC_000075.5" local="false" path=""/>
    <Dependency index="9" seq_id="NC_000076.5" local="false" path=""/>
    <Dependency index="10" seq_id="NC_000077.5" local="false" path=""/>
    <Dependency index="11" seq_id="NC_000078.5" local="false" path=""/>
    <Dependency index="12" seq_id="NC_000079.5" local="false" path=""/>
    <Dependency index="13" seq_id="NC_000080.5" local="false" path=""/>
    <Dependency index="14" seq_id="NC_000081.5" local="false" path=""/>
    <Dependency index="15" seq_id="NC_000082.5" local="false" path=""/>
    <Dependency index="16" seq_id="NC_000083.5" local="false" path=""/>
    <Dependency index="17" seq_id="NC_000084.5" local="false" path=""/>
    <Dependency index="18" seq_id="NC_000085.5" local="false" path=""/>
    <Dependency index="19" seq_id="NC_000086.6" local="false" path=""/>
    <Dependency index="20" seq_id="NC_000087.6" local="false" path=""/>
    <Dependency index="21" seq_id="NC_005089.1" local="false" path=""/>
   </List>
  </Dependencies>
 </Object>
 <SOFTWARE>
  <VDBLibrary vers="2.7.47"/>
  <Build static="true">
   <Module name=""/>
  </Build>
  <Tool date="Nov 18 2022" name="fastq-dump" vers="3.0.1">
   <Binary path="/usr/local/bin/fastq-dump" type="alias" md5="c461c39bfa514aff3c4f7c0416ced617">
    <Alias resolved="fastq-dump.3">
     <Alias resolved="fastq-dump.3.0.1">
      <Alias resolved="sratools.3.0.1"/>
     </Alias>
    </Alias>
   </Binary>
  </Tool>
 </SOFTWARE>
 <Env>
 </Env>
</Report>

How are you running the script? What arguments are you giving it? What is get_ngc? What are the files you are copying (cp "${file}" .), what do they contain? And why are you running this as root!? Never run things as root unless it is absolutely necessary. Any minor error or bug could result in destroying your system if you run the wrong command as root!

to test it locally I i first start the container and then I copy the script inside the container and run using a single GSM id as input here the issue is for this particular series of GSM id none of them works

'and why are you running this as root!? Never run things as root unless it is absolutely necessary' im running it in a container

Use the SRR accessions, not the sample-IDs you are using.
With entrez e-utils, I found the SRR of GSM2683458:
esearch -db sra -query GSM2683458 | efetch -format xml | xtract -pattern IDENTIFIERS -if PRIMARY_ID -contains SRR -element PRIMARY_ID is SRR5755657.
For GSM3603268, using esearch -db sra -query GSM3603268 | efetch -format xml | xtract -pattern IDENTIFIERS -if PRIMARY_ID -contains SRR -element PRIMARY_ID, the SRR is: SRR8571942
Then just do fasterq-dump SRR5755657; fasterq-dump SRR8571942.

well the pipeline is already made to use GSE id so i have no option but to use it


text:
tensorQTL interaction issue
I hope that the question has its place on this board.
I am trying to run the mapping of cis eQTL using an interaction term with tensorQTL. I am running into the following index error when running the software:
>>> cis.map_nominal(genotype_df, variant_df,
...                 phenotype_df,
...                 phenotype_pos_df,
...                 prefix = prefix,
...                 covariates_df=covariates_df,
...                 interaction_df=interaction_s,
...                 maf_threshold_interaction=0.05,
...                 run_eigenmt=True,
...                 output_dir=output,
...                 write_top=True,
...                 write_stats=True)
cis-QTL mapping: nominal associations for all variant-phenotype pairs
  * 569 samples
  * 24665 phenotypes
  * 68 covariates
  * 69763935 variants
  * including 1 interaction term(s)
    * using 0.05 MAF threshold
  * cis-window: ±1,000,000
  * checking phenotypes: 24665/24665
  * Computing associations
    Mapping chromosome 1
    processing phenotype 137/24665Traceback (most recent call last):
  File "/QRISdata/Q5776/software/tensorqtl/tensorqtl/core.py", line 197, in calculate_interaction_nominal
    Xinv = torch.matmul(torch.transpose(X_t, 1, 2), X_t).inverse() # ng x (1+2*ni) x (1+2*ni)
torch._C._LinAlgError: linalg.inv: (Batch element 2704): The diagonal element 3 is zero, the inversion could not be completed because the input matrix is singular.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/QRISdata/Q5776/software/tensorqtl/tensorqtl/cis.py", line 294, in map_nominal
    res = calculate_interaction_nominal(genotypes_t, phenotype_t.unsqueeze(0), interaction_t,
  File "/QRISdata/Q5776/software/tensorqtl/tensorqtl/core.py", line 200, in calculate_interaction_nominal
    i = int(re.findall('For batch (\d+)', str(e))[0])
IndexError: list index out of range

However, when I run the same command but only on 50 genes, the code runs properly.
>>> cis.map_nominal(genotype_df, variant_df,
...                 phenotype_df.iloc[1:50],
...                 phenotype_pos_df.iloc[1:50],
...                 prefix = prefix,
...                 covariates_df=covariates_df,
...                 interaction_df=interaction_s,
...                 maf_threshold_interaction=0.05,
...                 run_eigenmt=True,
...                 output_dir=output,
...                 write_top=True,
...                 write_stats=True)
cis-QTL mapping: nominal associations for all variant-phenotype pairs
  * 569 samples
  * 49 phenotypes
  * 68 covariates
  * 69763935 variants
  * including 1 interaction term(s)
    * using 0.05 MAF threshold
  * cis-window: ±1,000,000
  * checking phenotypes: 49/49
  * Computing associations
    Mapping chromosome 1
    processing phenotype 49/49
    time elapsed: 0.44 min
    * writing output
done.

The interaction factor is a categorical variable (0/1), I am not sure why the issue is being raised and how to fix
I would appreciate any help/direction to go forward.
Thanks

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

Writing a solution for anyone interested. The issue stems from a dosage error, some variant contains only heterozygous individuals and no homozygous ones needed for QTL calculation.
I used plink2 to calculate genotypes for all variant using the following command:
plink2  --pfile PATH \
     --geno-counts  \
     --out OUT

I then removed variant that contains less than 10 Homozygous invidivuals (HOM_REF_CT < 10) using plink then was able to run tensorQTL interaction mode.
Hope it helps.


text:
Align two different structures with respect to their active sites or heme cofactor
I have two enzymes, with heme cofactors. but their overall structure is very different.
I use Pymol regularly to visualize protein pdb structures.
I want to align both structures, but the structures don't align well due to very low homology. However, their active sites/substrate entrance channels may be more similar than the rest of the structure.
so,

I want to align only the substrate channels from both enzymes and compare their structure.
or
align both proteins only using the heme as a reference point. So both hemes from two enzymes align and then I can manually do the rest.

I am more interested in how the substrate would be accessible in both enzymes differently.
Thanks in advance.
Edit: Sorry for being a bit unclear. We can pick two proteins like cytochrome P450 BM3 with pdb id 6H1L or 6H1O. The second protein can be a similar peroxygenase UPO, e.g. 7pn4.
Thanks.

@terdon , do the same rules of SO apply to Bioinformatics , the question seems too wide , at least irfan could have named two pdbs code as inputs, toghrther with the aa defining the channels  to have us to try to get a result with

If you feel the question is too broad, then ask for clarifications and/or vote to close as "needs more focus". That said, I don't see why this is too broad: the main question here is "given two pdb structures, how can I align and compare only specific parts of those structures instead?" which seems like it should be answerable by those working with this kind of data.

Please add more context to this question; be as specific as possible about the problem you're struggling with. Example proteins, or other input files, are very helpful.

OK, trying to answer question 2, think I need to be sure that both enzymes have same identical heme cofactor and then find a rotation/translation matrix that applied to heme of protein A moved exactly on top of heme of protein B , than I have to move the entire protein backbone of protein A using that matrix, is it reasonable ? Am I getting it wrong ? With PyMOL I tried to use Pair/Pair_fit command with two heme proteins :
2d0t (https://www.rcsb.org/structure/2d0t and 4zv8 (https://www.rcsb.org/structure/4ZV8) :
pair hem_2d0t , hem_4zv8,  same as pair_fit hem_2d0t , hem_4zv8, where hem_2d0t,hem_4zv8  are the heme of the two proteins [select model 2d0t and resn HEM and set_name sele, hem_2d0t or select hem_2d0t, model 2d0t and resn HEM].
See results, superposed heme:

Superposed Proteins :

As for question 1 I tried to answer using the tools available to me.
I picked the two same proteins : 2d0t (https://www.rcsb.org/structure/2d0t and 4zv8 (https://www.rcsb.org/structure/4ZV8). Using parKVFinder (https://github.com/LBC-LNBio/parKVFinder , https://www.softxjournal.com/article/S2352-7110(20)30319-8/fulltext) pymol plugin on both proteins stripped of all their HETATM, I calculated their cavities. As an example , here a png of 2d0t  cavities :

Then I extracted the residues surrounding the cavity using :
select model 2d0t_merged near_to 5.0 of resn KAF ; [named sele_2d0t]
select model 4zv8_merged near_to 5.0 of resn KAB ; [named sele_4zv8]
being KAF and KAB cavities described as non_bounded atom in pdb format as output of parKVFinder .
I believe parKVFinder has is own output for the residues surrounding cavities, but wanted to try PyMOL
Example for 2d0t :

Then I used
cealign sele_4zv8 , sele_2d0t (https://pymolwiki.org/index.php/Cealign) to align the residues of the two cavities. (cealign uses just CA of the selections provided).
Results , only cavities and surrounding residues :

the two heme protein cavities are not alignrd even if the HIS (2dt0, red)  and CYS (4zv8, yellow) are on the same side :

NOTE:
2dt0 in magenta
4zv8 in cyan
I choose 2 different heme protein (could be more extreme though) , just because the question was open to broad attempts, could be interesting too see for example results on not too distant cytochromes

Thank you very much. I will try to follow the process and report back.


text:
For loop in Python goes wrong
I have mutation status (0,1) for ten genes in ten cancer 3D models
I have also gene expression in these 3D models, like
expr1 <- structure(list(GeneA_KO = c(0.0756383119639074, -0.162841965682723, -0.110703100204517, -0.457154616239399, -0.505815193885352, 0.366072775033448, 0.288650458935083, -0.0287646498640309, -0.690732824311905, 0.157366122620995, -1.29451546737684, 0.124167425460344, -0.168045313861751, 0.0186504943306465, -0.495880827400706, -0.206668627339351, -0.148609490335447, 0.066844529089804), GeneB_KO = c(-0.870157285953621, -0.35900872537626, -0.990406888128838, -0.120656981579292, -1.41427631791049, -0.866421095836781, -0.816936112146099, -0.894597236248507, -1.00084584207823, -0.850754488133456, -1.05314239761653, -1.07580675428372, -1.24708324665688, -1.04098460460899, -0.232388675985086, -1.48513703581865, -1.48373348870346, -1.76484058374354)), class = "data.frame", row.names = c("ModelM", "ModelO", "ModelS", "ModelA", "ModelI", "ModelC", "ModelT", "ModelH", "ModelE", "ModelR", "ModelA2", "ModelP", "ModelE2", "ModelU", "ModelT2", "ModelC2", "ModelS2", "ModelX"))
mut <- structure(list(ModelM = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelO = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelS = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelA = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelI = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelC = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelT = c(0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 0L), ModelH = c(1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 0L), ModelE = c(1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L), ModelR = c(1L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 0L), ModelA2 = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelP = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L), ModelE2 = c(1L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L), ModelU = c(1L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L), ModelT2 = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelI2 = c(1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelC2 = c(1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L),     ModelS2 = c(0L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 1L)), class = "data.frame", row.names = c("Gene1_mut", "Gene2_mut", "Gene3_mut", "Gene4_mut", "Gene5_mut", "Gene6_mut", "Gene7_mut", "Gene8_mut", "Gene9_mut", "Gene10_mut"))

mut.csv and expr1.csv exported for use in Python:
write.csv(mut, "mut.csv", quote = FALSE, row.names = FALSE)
write.csv(expr1[1:10,], "expr1.csv", quote = FALSE, row.names = FALSE)

I want to do a one-way ANOVA in Python so that calculate P-VALUE for each pair of mutation-expression in each pair of gene-model
This code returns error
import os
import pandas as pd
from statsmodels.formula.api import ols
import statsmodels.api as sm
m = pd.read_csv("mut.csv", sep=",")
m=pd.DataFrame(m)
m.info()
g = pd.read_csv("expr1.csv", sep=",")
g=pd.DataFrame(g)

lst = list(g.columns)
m_lst = list(m.columns)

dfs = pd.DataFrame({})
for Mut in m_lst:
    for Model in lst:
        df = pd.concat([m[[Mut]].reset_index(drop=True), g[[Model]].reset_index(drop=True)],axis=1)
        df = df.rename(columns={Mut: "mut", Model: "g"})
        model = ols('g ~ C(mut)', data=df).fit()
ValueError: endog has evaluated to an array with multiple columns that has shape (10, 18)

        aov_table = sm.stats.anova_lm(model, typ=2)
        aov_table[['Model']] = Model
        aov_table[['Mut']] = Mut
        dfs = dfs.append(aov_table)

The ideal output is something like
with pd.option_context('display.max_rows', None, 'display.max_columns', None):
    print(dfs)
            sum_sq   df         F    PR(>F)     Model      Mut
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO   ModelM
Residual  0.978877  8.0       NaN       NaN  GeneA_KO   ModelM
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO   ModelM
Residual  0.586718  8.0       NaN       NaN  GeneB_KO   ModelM
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO   ModelO
Residual  0.978877  8.0       NaN       NaN  GeneA_KO   ModelO

I can see what you're wanting to do, I can see bugs but the key bit of output is missing (below). It seems a bit weird to me the nested loop in context, but thats fine if you're happy.
The critical step is this
df = pd.concat([m[[Mut]].reset_index(drop=True), g[[Model]].reset_index(drop=True)],axis=1)

What I'm concerned about is the pd.concat. Are you sure pd.concat is the right command here? I thought m is one data table, e.g. mutation and g is another e.g. expression. pd.concat isn't going to work here, are you sure this isn't a merge or join command?
Anyway, the fitting is complaining the .shape() is wrong. Needs reshaping? Thats the error. Printing df at that point would be needed.

ValueError: endog has evaluated to an array with multiple columns that has shape (10, 18)

To summarise normally the error you have received is about re-shaping the data. I've not put forward that solution here because,

I am not sure what the data represents
If I was doing this I wouldn't use pd.concat

It seems a lot of columns for a 1-way ANOVA, but could be right, I've not done 1-way ANOVA across 10 columns.
Bug 1
m = pd.read_csv("mut.csv", sep=",")
m=pd.DataFrame(m)
m.info()
g = pd.read_csv("expr1.csv", sep=",")
g=pd.DataFrame(g)

This should simply be:
m = pd.read_csv("mut.csv")
g = pd.read_csv("expr1.csv")

Bug 2?
dfs = pd.DataFrame({})
...
dfs = dfs.append(aov_table)

This will work given aov_table output is a pandas dataframe (I'm not at all sure it is). You need aov_table.info() to be sure. You can always place it in a dataframe however.


text:
How do you download metadata for TCGA pathology slides
I downloaded the TCGA pathology slides for the SKCM project and the files are all a name such as "TCGA-3N-A9WB-01Z-00-DX1.A9950ED4-9480-455C-AE0D-8E076D4DA432.svs".
How do you get the clinical metadata for this dataset? When I click "clinical" for the project, the .tsv file has no case IDs that match to my file names.
https://portal.gdc.cancer.gov/exploration?filters=%7B%22content%22%3A%5B%7B%22content%22%3A%7B%22field%22%3A%22cases.case_id%22%2C%22value%22%3A%5B%22set_id%3AYavV3YoBvzAafV-N5zuW%22%5D%7D%2C%22op%22%3A%22IN%22%7D%5D%2C%22op%22%3A%22AND%22%7D

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

The patient ID is the first 10 letters, e.g. "TCGA-3N-A9WB-01Z-00-DX1.A9950ED4-9480-455C-AE0D-8E076D4DA432.svs" relates to the patient "TCGA-3N-A9WB". If you download the clinical information for SKCM, there is clinical info in the tsv for patient "TCGA-3N-A9WB". Does that make sense?

Thank you! That makes sense!


text:
How to obtain the KEGG mudule based on KEGG KO? (KEGG mudule enrichment)
With HUMAnN3, I got the result of the KEGG KO annotation (by humann_regroup_table:humann_regroup_table --input hmp_subset_genefamilies.tsv --groups uniref90_ko --output hmp_subset_genefamilies_uniref90_ko.tsv).
I need to quantify/transform the KEGG KO annotation to KEGG module definitions.
However, I did not find the relevant guidance/instruction. Are there any relevant tutorials available for reference?
Just like the packages ("clusterProfiler") using R for transform to KEGG pathway (KEGG pathway enrichment analysis).
Thanks a lot!
Any help will get me out of trouble!

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

The problem with the question is it's not clear what "quantify/transform" KEGG really means. Nor is it clear what your background in coding or SQL is because this is a fairly simple operation. Sure I've done loads of this, but whether that transferable here its difficult to assess.

@M__  Thank you for replying. "quantify/transform" KEGG means that the result of the KEGG KO annotation (genome information) is mapped to KEGG modules for the functional annotation like "KEGG Pathway Enrichment Analysis" (by R package: clusterProfiler). I'm just starting to learn these. Can you give me some advice? Thanks a lot.


text:
Merge fragment reads
I'm writing a function (python) that merges overlapping paired reads of a fragment, and a corresponding reference sequence. There's quite a lot of edge cases (e.g. when the reads disagree on indels) so I wonder if there's some standard library that does that?
Context - the input is pysam AlignedSegment DNA reads. The output is fragment, reference - both strings corresponding to the merged fragment and its corresponding reference.
e.g.
read1:           ACGTACGTACGT
read2:                CGTACGTACGT
expected output: ACGTACGTACGTACGT

read1:           ACGTANGTACGT
read2:               ACGTACGTACGT
expected output: ACGTACGTACGTACGT

this is just an illustration. The idea is to merge the overlapping DNA read sequences into the fragment and its corresponding reference using the information in the AlignedSegment (i.e. query_sequence, reference_start, query_qualities, cigartuples, and get_reference_sequence())
Thanks.

Thanks. But why is this not an assembler? This is pretty much exactly what an assembler does, it assembles reads into contigs. Or, if you already have the reference genome, what a read aligner does.

well, I mean it's not an aligner in the sense that's not aligning reads to a reference but to each other.

OK, so it's an assembler. But then what do you mean by the word "reference" in your question? Are you referring to the consensus sequence you can build by assembling the reads, the inferred genomic source sequence? Also, do you absolutely have to write this from scratch and in a scripting language? Can't you use one of the many existing genome assemblers?

I believe the OP wants to work at the "read-pair level" to obtain contigs of mate pairs and not all of the available reads. If that is the case, I would check make.contigs from the mother suite.

I would use existing high-performance software for this, for example BBmerge from BBmap suite, or FLASH, probably many others exist as well. I see no point in re-writing new software when it already exists, unless this is a learning exercise, or you think you can do better algorithm-wise, but then, with respect, you would probably not be asking here.

I'd suggest looking at similar tools that already exist, e.g. here is one that I hadn't heard of that lives on github. Just glancing at the code, it seems to use a scoring heuristic to decide how to derive a merge point, so it's using that to get around the absence of unambiguous points, probably at some rate of error.
BioPython has some functions that will create consensus sequences for you as well based on a pairwise alignment. (Bio.Align.AlignInfo.SummaryInfo.dumb_consensus, linked in that post.)
As a command line tool pandaseq does the same thing, I'd suggest looking at how they handle edge cases.
There are similar utilities in QIIME, NGMerge, BBMerge, [illumina-utils]- this is a pretty well-trodden path, I think all of those will be illuminating to look at in terms of reusing or writing original code.
I'm guessing it will be easier to merge FASTQs and then align, rather than align and then merge. (Note that these tools don't seem to operate on alignments, because that is trickier).
Alternately, if your real goal is to create a consensus genome sequence from your reference based on your alignments, there are a large number of tools that do that, notably samtools (see e.g. https://www.biostars.org/p/419927/). No new code needed, you just input your BAM file and reference sequence.


text:
How can I compute energy from a protein PDB file?
Suppose I have a protein PDB file and want to compute energy using the following formula from the coordinates in the file.

How can I do that?
Sometimes people suggest using the CONECT field for that purpose. However, some PDB files do not have this field.
Example: 4OSK.pdb (Crystal structure of TAL effector reveals the recognition between asparagine and guanine)

https://mattermodeling.stackexchange.com/questions/11678/how-can-i-retrieve-bond-information-froma-protein-pdb-file

https://www.ebi.ac.uk/pdbe-srv/pdbechem/bond/list/CYS Molecule : CYS
Bonds between two atoms of the molecule , then of course in each protein chain you'll have a Nterminal and Cterminal residues and all the ones in between having a peptide bond CO-N between Cter and Nter

http://ligand-expo.rcsb.org/pyapps/ldHandler.py?formid=cc-index-search&target=cys&operation=ccid here with pics http://ligand-expo.rcsb.org/reports/C/CYS/CYS_D3L1.gif

think you'll need to create a dictionary for each residue

see here for better explanation https://chemistry.stackexchange.com/a/162689/118509


text:
SARS-CoV-2 sequence used in the AlphaSeq Antibody Datasets to predict binding affinity
Currently we are building a sequence based deep-learning model to predict binding affinity between antibody and antigen. For this we are training a sequence based model with AlphaSeq Antibody dataset (https://github.com/mit-ll/AlphaSeq_Antibody_Dataset). But in this dataset they haven't provided the sequence for SARS-CoV-2 virus.
So I read the paper and there they have mentioned they have tested against peptide in the HR2 region of the SARS-CoV-2 spike protein and they have given the sequence as "PDVDLGDISGINAS". (https://www.nature.com/articles/s41597-022-01779-4#Sec3)
But this sequence is very small. I'd like a larger protein sequence that includes non-epitope portions so that the prediction is more meaningful (e.g. predicts binding to the actual epitope when provided a sequence which includes sequence beyond the epitope). I searched online and only found a covid genome sequence which isn't very helpful.
Searching the RCSB PDB website for Spike protein HR2 gave two SARS-CoV-2 results (RCSB PDB website).
One of them had sequence like the following:
>8CZI_1|Chains A, B, C|Scaffolded Spike protein S2' HR1|Nostoc punctiforme (strain ATCC 29133 / PCC 73102) (63737)
MSHHHHHHGSQTLLRNFGNVYDNPVLLDRSVTAPVTEGFNVVLASFQALYLQYQKHHFVVEGSEFYSLHEFFNESYNQVQDHIHEIGERLDGLGGVPVATFSKLAELTCFEQESEGVYSSRQMVENDLAAEQAIIGVIRRQAAQAESLGDRGTRYLYEKILLKTEERAYHLSHFLAKDSLTLGFAYENQKLIANQFNSAIGKIQDSLSSTASALGKLQDVVNQNAQALNTLVKQLSSNFGAISSVLNDILSRLDKVE

>8CZI_2|Chains D, E, F|Spike protein S2' HR2|Severe acute respiratory syndrome coronavirus 2 (2697049)
KNHTSPDVDLGDISGINASVVNIQKEIDRLNEVAKNLNESLIDLQ

So I am not sure which one is the correct sequence for our task. Can someone provide some insights regarding this matter?
Also why are there so many sequences popping up when I search for SARS-CoV-2? Is there a version that can be considered as the root from which other versions are derived?

Can you please replace the image with the actual sequence? The image is hard to read, impossible to copy from and just makes understanding your question harder. You can use the {} button to format it as code.

Is there a reason you don't just use the peptide sequence that they indicate as an antigen? if you're trying to replicate their results then why do you need a different sequence?  Note that your RefSeq COVID link is to the genome (DNA) not protein sequence.

@MaximilianPress Yeah that's a valid question. Actually I wanted to input a bigger sequence to our model rather than the specific part of the whole sequence to which it actually binds, so that our machine learning model learns to guess the correct region from the sequence as well as the correct binding affinity.

@MaximilianPress thanks for indicating that's a DNA

Thanks for explanation, I've edited your question to include this.

PDB has a whole page talking about various slices of the covid spike protein. I'd strongly recommend reading that for a better understanding of why people choose different pieces to represent. In that page they include several different solved structures of the protein that may be helpful for your benchmarking.
In terms of the "canonical" sequence, you could look at a reference genome sequence. Here for example is the EMBL reference genome for covid in GenBank format, which includes the protein sequences. I believe in that the spike protein has ID QHD43416.
You can additionally look at NCBI, which has a slightly different version of the reference (though they should be functionally identical, with just slightly different identifiers and annotations).
Assuming that you want the ancestral initial isolate of covid, you can find a similar GenBank file here from NCBI.


text:
Biopython PhyloTree to PyVis Network
We are trying to convert a Phylogenetic tree to networkx previously calculated with
from Bio.Phylo.TreeConstruction import DistanceMatrix
from Bio.Phylo.TreeConstruction import DistanceTreeConstructor
from Bio.Phylo import draw

lower_triangle = []

for i in range(0, len(distance_matrix)):
    lower_triangle.append(list(distance_matrix[i][:i+1]))

lower_triangle

# Make a list of strings from '0' to '9'
labels = list(map(str, range(0, len(lower_triangle))))

distance_matrix = DistanceMatrix(names=labels, matrix=lower_triangle)

constructor = DistanceTreeConstructor()
upgma_tree = constructor.upgma(distance_matrix)

and convert with
from pyvis.network import Network
from Bio import Phylo

g=Network(height=500, width=500)
g.toggle_physics(True)
network_x_graph = Phylo.to_networkx(upgma_tree)
g.from_nx(network_x_graph)

It fails with
assert isinstance(n_id, str) or isinstance(n_id, int)

here
https://github.com/WestHealth/pyvis/blob/ccb7ce745ee4159ce45eac70b9848ab965fc0906/pyvis/network.py#L236
so it seems the nodes are of wrong type? How can I check / adjust it?
print(upgma_tree)

Tree(rooted=True)
    Clade(branch_length=0, name='Inner9')
        Clade(name='2')
        Clade(name='Inner8')
            Clade(name='7')
            Clade(name='Inner7')
                Clade(name='5')
                Clade(name='Inner6')
                    Clade(name='6')
                    Clade(name='Inner5')
                        Clade(name='4')
                        Clade(name='Inner4')
                            Clade(name='Inner1')
                                Clade(name='9')
                                Clade(name='1')
                            Clade(name='Inner3')
                                Clade(name='3')
                                Clade(name='Inner2')
                                    Clade(name='8')
                                    Clade(name='0')


text:
How many ML tree should I build?
I created a phylogenetic tree of 328 species using RAxML. I run the ML algorithm 20 times and I chose the best tree (the one with the highest likelihood).
The choice of 20 iterations was based on a tutorial that was carrying out the same analysis but on fewer species.
Now I have a doubt that my search is actually not optimal, and could have got better by increasing the number of iterations.
Do you know if there is a way to estimate the optimal number of iterations to find the best robust ML Tree?

Simply use the iterations option (e.g. -n 1000) of IQTREE.
It's a lot quicker (several logs) and is strongly of identical performance to RAxML.
iqtree -s alignment.phy -n 1000 -m GTR+I+G -nt 4

It takes phylip, fasta, nexus  ... and no-one questions 1000 iterations, so you're safe.
I personally I'd then go on to bootstrap the tree (IQTREE is extremely fast). If you were bootstrapping with iterations RAxML is unworkable, you'd need e.g. 10 iterations per bootstrap (1000 bootstraps) and would be time consuming. IQTREE could easily perform 10 iterations per bootstrap and would accommodate 100 iterations (1000 would be a bit excessive for bootstrapping).
RAxML is an old program now, IQTREE is its successor and massively refactored for maximum likelihood phylogeny.

If you don't understand the command line, drop a note and I'll explain the model it's using.