text:
Does using a microarray chip that matches the chip used in training lead to higher explained variance of a polygenic score?
I'm trying to replicate an existing polygenic score (i.e. test the accuracy in a new sample), and want to know if matching the original study's microarray chip will improve the accuracy (that is, variance explained in a new test set) of the polygenic score.
The labs I have access to have different multiple microarray chips to choose from (e.g., GSA, GDA, UK Biobank Axiom Array, UK BiLEVE Axiom array, GeneChip 2.0 array, Infinium Core-24 Kit, etc.). Will using a chip that is different than the chip used in the study that produced the polygenic score impact the explained variance (accuracy) of the polygenic score? If so, by how much?
Example: a study on heart disease runs a GWAS on 900,000 individuals who are genotyped with the UK BiLEVE Axiom array. This study uses these summary statistics to produce weights for a polygenic score. Lab A, B, and C are conducting polygenic score validation/replication studies, and they each have a random sample of 10,000 individuals drawn from the general population, and they all use the same summary statistics from the original study. After imputing the genotype data of the 10,000 people, they calculate the polygenic score for each person, and then evaluate how much variance in the phenotype they are able to explain. The only difference is that Lab A genotypes their 10,000 people with the UK BiLEVE Axiom array, whereas lab B uses the Illumina GSA, and lab C uses the Illumina GDA.
Do we have any reason to suspect that there will be a substantial difference in variance explained (r^2) by the polygenic score between labs? Are there any studies which have analyzed something that could be informative to this question?
I've been unable to find research on this question. This should not be dependent on which polygenic score is being studied; there is no reason to specify a particular phenotype.
Note: it is standard to impute both in the original study and the replication.

Genotyping can be validated by whole-genome sequencing, and this has been done on the International Genome Sample Resource. You should be less concerned about the reproducibility of genotyping results, and more concerned about the correctness or validity of genotyping results. Companies will often revisit old SNPchips to mark invalid spots, update rs numbers, and add new genotyping probes that are better on a population-wide scale. In that sense, choosing a SNPchip platform that has an established track record of QC and information about version updates would be a better idea than trying to match a previous study as precisely as possible.
However, with regards to replicating existing studies, a bigger concern for you should be the variation in sampled populations, which plays a big role in the outcome of a polygenic test.
This variation is so extreme that directions of associations can flip for the same variant, especially for low-frequency variants with extremely low p-values (i.e. the ones that often show up as manhattan plot peaks).
The issues around population and study replication are explained well in the following paper from 2001:
https://doi.org/10.1006/tpbi.2001.1543
This is not something that can be fixed by attempting to replicate studies in more populations. As Cecile Janssens says, larger sample size only affects random error, not systematic error.
https://twitter.com/cecilejanssens/status/1119235476567801856


text:
Need to make circos plot from Enrichment data
I have a list of genes and their associated pathways. The data was generated using Partek. The data looks like this -

Gene_in_list
Gene_set
Description
Enrichment score
P.value

Gene1
hsa
Disease C
69
0.005

Gene1
hsa
Disease A
45
0.05

Gene2
hsa
Disease B
23
0.007

Gene2
hsa
Disease Y
40
0.009

Gene2
hsa
Disease C
69
0.0001

I want to 'split open' the Description column (Which I think means going long to wide in tidyverse terminology?), and transform it to the below table -

Gene_in_list
Disease A
Disease B
Disease C
Disease Y

Gene 1
1
0
1
0

Gene 2
0
1
1
1

So that I can use the online circos plot to essentially show, which genes are involved in the pathways.
I was also debating if Enrichment score / p.value could be used as the thickness of connections in circos, but I think it might look too thin / uneven.
Any help appreciated, thanks!

source: Stack Exchange

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

sort: votes


text:
Adding step size to cross validation in Python package lifelines
I'm using the Python package lifelines to build a CoxPH model where I reduce the step size to 0.1, like so:
cph = CoxPHFitter()
cph.fit(df, 'time', 'event', fit_options = {"step_size":0.1})

I used k_fold_cross_validation to build the model, but I don't know how to add in step size. I expected to add it to fitter_kwargs, like this, but this does not work:
k_fold_cross_validation(cph, df, duration_col='time', event_col='event', k=5, fitter_kwargs={fit_options:{"step_size":0.1}})

have you tried simply cph.fit(df, 'time', 'event', step_size=0.1} ?

@JRodrigoF I wouldn't be able to do cross-validation then

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


text:
retrieving past BLAST requests in bioPython

is there any way to retrieve a BLAST request made using the Biopython library using the request's query_id?
Alternatively is it possible to retrieve a query_id before the blast result is complete, and checking on it later?
(bit off-topic) I seem to only retrieve 60 results even at higher  hitlist_size although i get much more results when blasting the same entry on uniprot

I'm working on a script which processes a BLAST result and having to re-request the search every time is quite temporally cumbersome since i'm still learning the library. currently i'm using:    blastResult = NCBIWWW.qblast("blastp", "swissprot", targetAminoAcids,format_type="Tabular", hitlist_size=200)

Welcome to the site. This sounds like you just want to add some cache, is that right? Thats fairly easy via Python. Setting up a cache database is more complicated however - but absolutely doable - details about when you are running the different requests and how would help here.

Please [edit] your question and give us i) a script that we can run locally to reproduce your issue and ii) an example input so we can test it.

1 and 2 —> Yes, but you need to have the query ID in order to do so: having the description would be a bonus, but not required.  If you go to 8.1.1  Query Result of the Biopython manual, it will tell you how to get a query id, if you don’t already have it!  Go to section 7.1 as a reminder of how to run  BLAST via the internet.  Go to section 7.2 to. Run BLAST locally.   Weigh your options out with running via Internet or locally.  Running BLAST locally is faster usually and gives you more user freedom to create your own database with your own parameters, etcetera., but certain tools will need to be downloaded via the command line/terminal.  Follow this link for the BLAST command line user manual:
https://www.ncbi.nlm.nih.gov/books/NBK279690/
I will paste the link to the Biopython user manual if you do not already have it:
http://biopython.org/DIST/docs/tutorial/Tutorial.html#sec123
3 —> Try the Protein Data Bank (PDB) instead of UniProt or you can combine both UniProt and PDB (I recommend this route as it’s a combination of the 2 databases).  You can use these two databases to process a BLAST result!  I have found in my experience, combining data from more than 1 database has been beneficial in my case!


text:
Getting ALL Protein entries from Entrez for a particular taxonomic id
I'm trying to get ALL Proteins for a Taxonomic ID from Entrez using Esearch/Efetch and the following command:
esearch -db "protein" -query "txid374666[Organism]" | efetch -format fasta > txid_374666.fasta

Can somebody tell if this command will get ALL proteins or not? Besides, this appears to be a bit slow.
Can somebody point out any Biopython workflow, particularly how to get ALL Proteins
Thanks

That would work - except for the precise search term used. The issue is the NCBI server may time out, certainly if this was Blast it would do. Doing a single command line would prevent checking. Personally I would use a script to check what was downloaded with the esearch output (which will certainly work). If I get chance I will provide a full answer, e.g. search term stuff.

Thanks for your response. I'll await your answer

I've to download a huge lot of sequences for hundreds of bacterial genus (not species)  for a metaproteomics study. Which method will work best for me? esearch/efetch, python scripting , anything else?

If I've time I will answer in detail about how I'd do this. In general further details would help (even if not me answering), i.e. how many different species are needed? This sounds like a script more than a one liner.

Around 700 genus are there, and under each jenus on an average 3 species. I'm trying to download using esearch/efetch . I belong to University, where I've good internet speed , yet, the link is breaking very often

Hi there, not sure if it can help but I asked a similar question some times ago here : https://stackoverflow.com/questions/74727918/handle-500-server-error-type-during-api-request-using-rentrez/75055203#75055203. I am using R however. General idea was that it was taking too long so I had to deal with server errors using a try/catch/sleep function.

As suggested by somebody, the datasets tools from NCBI  is blazing fast:
datasets download genome taxon 187492  --include protein

NCBI Datasets Tool

Thanks for sharing this, its really useful.

If you do get chance @TheAugust please mark your answer as "accepted" (helps with site stats). I recommend upvoting the answer because its efficiency allows local blast (and other methods) to interrogate data when remote NCBI blast is inefficient. NCBI servers do have a restricted access policy for some of their services.


text:
Blastn through NCBIWWW not giving any results, as opposed to the online BLAST
I am trying to create a simple pipeline that takes a file containing various short sequences, and run them through blast. I have an example sequence;
>some_random_sequence
AAGGTTCGGTCCAAATTGAA

which returns a hit for FUT3 (using the RefSeq_Gene DB and using grch38 as an entrez_query). However, when I try to run this using NCBIWWW (Biopython) I don't get any results and I would like to figure out why.
A snippet that contains the call that I am using, is listed here;
from Bio.Blast import NCBIWWW
from io import StringIO

sequence=['>some_random_sequence\n','AAGGTTCGGTCCAAATTGAA']

lines = "".join(sequence)
file = StringIO(lines)
results = NCBIWWW.qblast(
    program="blastn",
    database="RefSeq_Gene",
    entrez_query="grch38",
    sequence=file.read(),
    # Manually setting the parameters because short_query checks the length
    # of the sequence, but for some reason includes the name (resulting in
    # a length of 42, and thus defeating the entire purpose of using
    # short_query (which does a if len(sequence) < 31).
    expect=1000,
    word_size=7,
    nucl_reward=1,
    filter=None,
    lcase_mask=None,
    # short_query=True,
)
with open('test_output.xml', 'w') as save_file:
    blast_results = results.read()
    save_file.write(blast_results)

This seems like it might be an XY problem. Could you please add some more context around the question that you're asking? What are you trying to do with BLAST? What do you plan to do with the BLAST results once you see them?

This is easy. If you set
expect=1000,

to
expect=100,

It will work. NCBI servers impose a tighter restriction on the CLI whilst the web restriction is 5000.
If you require a 1000 restriction then out of (US) hours blasting is required, e.g. weekends.

The Perl Blast API is here
https://blast.ncbi.nlm.nih.gov/docs/web_blast.pl
The documentation is here https://ncbi.github.io/blast-cloud/dev/api.html
Normally you need to modify it, that would be a separate question. This will work in default mode.

I will try this straight away, I am curious if that was documented anywhere though? as I didn't recall any mention of that in the doc of NCBIWWW's function (could have missed it).

This would also mean that the short_query parameter can only be run on the weekends, as it sets expect to a 1000 ponder.

I have raised the issue. If this doesn't work there is a direct Perl API that definitely works for 100.

I am still waiting for a response of the snippet, for well over 10 minutes, as such I can't say if it is working, yet.

@BasJansen the Perl Api is at the web address posted above.

Much appreciated, the expect was indeed the issue it seems.

Many thanks, glad its worked


text:
New to DNABERT confused about kmers input
I'm trying to use DNABERT to predict promoter regions. Right now, I have a df with 500k rows, 1st column is the DNA sequence (all different lengths), 2nd column says "promoter" vs. "non-promoter". I'm trying to understand the GitHub, read the paper, as well as other background reading, but I'm still confused.
If I want to get started, they said it's optional to pre-train (so you can skip to step 3). This is where I got tripped up: "Note that the sequences are in kmer format, so you will need to convert your sequences into that." From what I understand, you need to do this so that all of the sequences are the same length? So kmer=6 means all of the sequences are length 6? Someone suggested that I take the first nucleotide in the promoter and grab 3 nucleotides before and 3 nucleotides after (+/-3 bases). I don't think that's how the kmer thing works though? I tried replicating how I think it works down below (I got confused on the last row of the 'after' df). Please correct me if I'm wrong!
Also, if there is a good package that does this for me, please let me know!
Before:
seq                          label
AAAGAGGT                      1 [promoter]

After:
seq                          label
AAAGAG                         1
AAGAGG                         1
AGAGGT                         1
??[GAGGT]??                    1

Edit the post and add the exact command you used.

I didn't use any command to get the 'after' rows. I just did this by hand.

The next line in the manual speaks of a python function they provide that can get k-mers from your sequence, so you don't need to worry about doing this a different way. The function slices the sequence starting from the first base with a sliding window of k to obtain k-sized pieces of the sequence and then space-separates these pieces into a single string.
seq="ATGCCGTA"
k=2
[seq[x:x+k] for x in range(len(seq)+1-k)]
['AT', 'TG', 'GC', 'CC', 'CG', 'GT', 'TA']

" ".join([seq[x:x+k] for x in range(len(seq)+1-k)])
'AT TG GC CC CG GT TA'

k=3
" ".join([seq[x:x+k] for x in range(len(seq)+1-k)])
'ATG TGC GCC CCG CGT GTA'


text:
HISAT2: RNA strandedness
My library is unstranded and the code I'm trying to use is this:
# Run alignment with unstranded method
hisat2 -q --rna-strandness none -x HISAT2/grch38/genome_index_grch38 -U hipsci/Kolf-1_R1_001.fastq | samtools sort -o HISAT2/Kolf-1_R1_001.bam

echo "HISAT2 finished running!"

However, I get this error:
(base) ~ % RNASeq_pipline/data/script/RNASeqpipline.sh

Error: should be one of F, R, FR, or RF
Error: Encountered internal HISAT2 exception (#1)
Command: /usr/local/bin/hisat2-align-s --wrapper basic-0 -q --rna-strandness none -x
HISAT2/grch38/genome_index_grch38 -U hipsci/Kolf-1_R1_001.fastq
(ERR): hisat2-align exited with value 1
[W::hts_set_opt] Cannot change block size for this format
samtools sort: failed to read header from "-"
HISAT2 finished running!
0 minutes and 0 seconds elapsed.

Looking at the HISAT2 Manual:
For single-end reads:

Use 'F' if the read corresponds to a transcript.
Use 'R' if the read corresponds to the reverse complemented counterpart of a transcript.

For paired-end reads:

Use 'FR' if the first read in the pair corresponds to a transcript on the forward strand, and the second read corresponds to the reverse strand.
Use 'RF' if the first read in the pair corresponds to a transcript on the reverse strand, and the second read corresponds to the forward strand.
When you use the --rna-strandness option with either 'FR' or 'RF' for paired-end reads, HISAT2 will assign an XS attribute tag to each read alignment, indicating whether the read belongs to a transcript on the '+' (plus) or '-' (minus) strand of the genome.

It doesn't say what to input in the case of an unstranded library. Any help would be greatly appreciated!

According to the help produced when I type hisat2 --help, the default strandedness when running HISAT2 is unstranded:
  --rna-strandness <string>          specify strand-specific information (unstranded)

So taking out the rna-strandedness option should map the reads as if they were unstranded:
# Run alignment with unstranded method
hisat2 -q -x HISAT2/grch38/genome_index_grch38 -U hipsci/Kolf-1_R1_001.fastq | samtools sort -o HISAT2/Kolf-1_R1_001.bam

If this is a problem for you, I recommend you create an issue about it on the HISAT2 github page:
https://github.com/DaehwanKimLab/hisat2/issues/new

Kind of a bad move to not have a way of explicitly recreating the default behavior. There should be an option you can pass to rna-strandness that replicates not specifying that parameter at all. The usage reads like specifying --rna-strandness unstranded will work but it exits if the arg is not one of those 4 allowed values.


text:
How can I limit blast results to exclude a specific taxonomic ID?
I would like to perform blast+, using the command line, here is the code:
-task blastn -db nr -outfmt 7 -query myQuery.fas -out result.txt -remote.

I have to exclude SARS-CoV-2 hits from my results, so I must include -negative_taxids mentioning SARS-CoV-2 NCBI taxid number. However, the problem is, -negative_taxids is NOT compatible with -remote command and my command line encounters an error.
I have to use -remote as I cannot download entire +400 GB NR database to run it locally.
Could you please kindly help me out with this contradiction in my work?
My output txt file for this command line is:
# Fields: query acc.ver, subject acc.ver, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
# 600483 hits found
AstV-VA1    MN464146.1  100.000 32  0   0   1   32  6362    6393    7.14e-06    59.0
AstV-VA1    MN148428.1  100.000 32  0   0   1   32  6685    6716    7.14e-06    59.0
AstV-VA1    MK671313.1  100.000 32  0   0   1   32  6682    6713    7.14e-06    59.0
AstV-VA1    MK671312.1  100.000 32  0   0   1   32  6682    6713    7.14e-06    59.0
AstV-VA1    MK671311.1  100.000 32  0   0   1   32  6682    6713    7.14e-06    59.0

Have you considered running your query against all of nr and then removing the SARS hits from the results instead?

Thanks terdon for your constant and positive contribution to my questions.
Actually, I have already done it and I got +600,000 hits but the point is the txt result file contains only the following information:
"query", "acc.ver", "% identity", "alignment length", "mismatches", "gap opens", "q. start", "q. end", "s. start", "s. end", "Evalue", "bit score"
so, there are no organism and taxid columns in the result. One possible solution can be using R to download the annotation to add the results automatically, but it has two disadvantages:

downloading annotation of +600,000 hits will take around 200 hours (more than a week)
not all hits automatically yield annotations.
Therefore, I think filtering the results from the beginning can be a better option.

No, no, there is no reason to download entire annotations! All you need to do is extract the IDs of the target sequences and then get the taxids for those sequences. If you [edit] your question and show us a few examples of your output, we should be able to help, but it will just be a simple entrez query whose input is the list of sequence identifiers and output will be the associated taxids.

What do you want to do with these BLAST hits?

This isn't particularly elegant, but it should work. First, extract the accessions into a file. Next, get the taxid for each of the accessions. Once you have that, you can identify the accessions that are for SARS-CoV-2 sequences, and remove them from your output. Something like this:

Extract the accessions into a new file
awk '!/^#/{print $2}' blast.out > names

Map those accessions to their taxids.
while read -r accession; do
   printf '%s\t%s\n' "$accession"  \
       $(curl -s "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=nuccore&id=${accession}&rettype=summary&retmode=text" |
          awk -F'[<>]' '/Name="TaxId"/{print $3}');
done < names > taxids

With your example input, that produces:
$ while read -r accession; do
>    printf '%s\t%s\n' "$accession"  \
>        $(curl -s "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=nuccore&id=${accession}&rettype=summary&retmode=text" |
>           awk -F'[<>]' '/Name="TaxId"/{print $3}');
> done < names
MN464146.1   1435485
MN148428.1   2661902
MK671313.1   1239566
MK671312.1   1239566
MK671311.1   1239566

Now, this will probably be way too slow, doing them one at a time. So if you have many IDs, try and combine multiple queries into one. Like this:
perl -ne 'chomp; print $. % 50 == 0 ? "$_\n" : "$_," ' names |
 | sed 's/,$/\n/' |
    while read accession; do
      curl -s "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=nuccore&id=${accession}&rettype=summary&retmode=text" |
        awk -F'[<>]' '/Name="Caption"/{ acc=$3 }
                   /Name="TaxId"/{printf "%s\t%s\n", acc, $3}'
    done  < names > taxids

Here, the first Perl command will read through the list of accessions and print the accession followed by a comma and, once every 50 lines ($. % 50 == 0), it will instead print it with a newline. The sed command will replace a trailing comma with a newline for the very last line of output where, unless the number of accessions is an exact multiple of 50, you will have a line ending with a comma instead of a newline.
The result is a list of 50 comma-separated identifiers per line (fewer, on the last line). This is then fed to the curl command and then a naive awk that parses it and prints out the accession and the taxid.
Important: this approach loses the version information. So MK671311.1 becomes MK671311. This shouldn't be a problem, and we will take it into account later, but it should be understood. This is what the output looks like:
MN464146    1435485
MN148428    2661902
MK671313    1239566
MK671312    1239566
MK671311    1239566

Remove the SARS-CoV-2 sequences.
The taxid for SARS-CoV-2 is 2697049, so you can get the list of its accessions with:
awk -F'\t' '$2=="2697049"{print $1}' taxids

You can now filter your file to remove them using:
awk '(NR==FNR){
       if($2=="2697049"){
         sars[$1]
       }
       next
      }
      {
        acc=$2;
        sub(/\.[0-9]*$/, "", acc);
        if(acc in sars){ next }
          print
      }' taxids blast.out > blast.filtered.out

This awk script will read the first input file (NR is the current input line number while FNR the line number of the file being read; the two are only equal while reading the first file) and if its second field is the desired taxid, it saves the first field, the accession, in the array sars. The next ensures we don't do anything else while reading the first file.
When we move to the second file, we make a copy of the second field, the accession, and save it as acc, and then we remove the version information ( sub(/\.[0-9]*$/, "", acc). If the resulting accession is found in the sars array, we skip the line with next and if not, we print.
This assumes you are using the faster, batch retrieval I gave above, the one that removes the versioning information. If you instead use the very slow first one, you can filter with
 awk '(NR==FNR){ if($2=="2697049"){ sars[$1]} next }
      { if($2 in sars){next} print}' taxids blast.out > blast.filtered.out

That should give you a filtered blast output file that no loner has sequences from taxid 2697049. However, please make sure to test this, count the number of input and output lines and make sure everything makes sense. I don't have access to your actual data and while this does work with the example you provided, I cannot guarantee it will work out of the box for your real data. As always, you need to be careful!

@FarzadBeikpour kindly consider upvoting. This is a good answer and a lot of work. If you can't get -negative_taxids working I'd use this approach

I would personally look at is the version number of the blast you are using.
Bugs around blast+ versions >= 2.9.0 (April 2019) for -taxids and I would assume -negative_taxids flag should have been fixed as described in the the Blast+ release notes here.

Fix taxID filtering combined with mask-based alias BLAST databases.

Bug fix descriptions are usually short, hence the lack of further detail.
Blast taxids can be implemented via lists here https://www.ncbi.nlm.nih.gov/books/NBK569846/


text:
DNA methylation mediator analysis for environmental exposure and health outcome
We measured blood DNA methylation using epic array in 100 no-diabetic and 100 type 1 diatbetic patients matched for age, sex and smoking and also measured arsenic exposure using ICP-MS. Could you please suggest best way to find CpG sites associated with arsenic exposure and diabetes. We know that toxic element leads to health risk but we want to study this by using DNA methylation as mediator for disease risk.
Could you please suggest any R package for doing these types of analysis. I have used previously limma for toxic element exposure and DNA methylation association analysis but not sure how make connections with health risk phenotype?


text:
What are some good miRNA databases?
I am writing a project proposal for my boss and it is on codon usage bias in the dystrophin protein and the spectrin-like repeats found in a hmmscan file using the longest dystrophin isoform protein sequence which I got from a GenBank file!
There is a lot of research on how mRNA and tRNA impact codon usage bias and codon frequency.  Long story short, I work at a gene therapy lab as a bioinformatician doing research for muscular dystrophy.  Our lab is very interested in miRNA, and I have found a certain amount of articles that show miRNA may impact codon usage bias since miRNA can either cleave or translationally repress mRNA.
I found a miRNA database called www.tools4mirs.org. It seems like a potentially good database, but I don't know if it's the right database geared for this project.  We saw some unusual lack of protein expression when there should be a lot of protein expression and our lab is really focused on an interested in miRNAs.
Are there any other miRNA databases/tools that you can recommend?
I code in Python and know some AI/ML and would like to make this project a AI/ML bioinformatics project, but I do not know what some good miRNA databases are, since I want to pull large amounts of miRNA data for my project that is relative to muscular dystrophy.
So to sum up my question, what are some really heavily and well-curated miRNA databases besides tools4mirs.org ? I am kind of new to this field, and would like to hear any suggestion since I am still learning the inner workings of the field of bioinformatics and the possibilities and directions it can go in.

Welcome to the site. I would first use Ncbi's SRA. Every miRNA study will need to upload there and I'd use the sratoolkit.

Good databases for what? What kind of data do you need? Sequences? Functional annotations? Variants? Population data? What species? A specific one? Any? What does protein expression have to do with it?

miRNA databases with functional annotations and structural annotations for humans is what I am looking for.  I am looking at the concept of translation/reverse translation and miRNA regulation of mRNA which impacts tRNA/anti-codons


text:
Downloading a list of all RCSB-PDB entries
I'm trying to download a list of all of the non-redundant protein chain IDs in the RSCB-PDB (not the actual FASTA entries, just the IDs), filtered by some % identity.
Does anyone know how to do this?  I haven't been able to figure it out with the advanced search function .

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.

% identity to what? To each other? To a specific query sequence?

Please add some additional context to your question. Why dod you want to download this list?

@terdon I believe he needs the ones not identical AKA different ones ' non-redundant' as he says

So the ones that do not have more than X% identity to any others in the database, @pippo1980? OK, so what identity? And won't all of them be non-redundant? I mean, if they are identical, they wouldn't have a different ID, would they? (I know, my deep, deep ignorance of the protein structure field is probably evident here)

@terdon I dont know a lot about it either but I guess you are right the OP should specify a number and whath the number stands for it coukd be both sequence identity or structure identity but probably for the latter we are talking about mainchain Root Mean Square Deviation (RMSD),

but maybe is just this : https://www.rcsb.org/news/feature/61b8e5e869569d045d35a9a3 Remove Redundancy in Search Results

@terdon I'm a computational biologist with a background in biology, not computers. I did my PhD work on gene prediction and comparative evolution, my post-doc on protein-protein interaction networks and I have been working in NGS analysis, specifically human variant calling and annotation, for the past several years. Jack of many trades and master of none. .. I know, my deep, deep ignorance of the protein structure field is probably evident here ... YES SURE .... ;-)

Did you try RCSB sequence similarity search?
These two announcement highlight the options you ask about:

Improved Searching with New Sequence Identity Cutoff (2017)
Remove Redundancy in Search Results (2015)

Alternatively, if you'd like to calculate identity on your own, you can download sequences of all the PDB entries from:
https://files.wwpdb.org/pub/pdb/derived_data/pdb_seqres.txt
(270MB).

Thank you!  The https://files.wwpdb.org/pub/pdb/derived_data/pdb_seqres.txt link was what I needed.  Any way to filter them by identity from the text file?

I suppose you'd need to do pairwise sequence alignment and calculate identity.


text:
Annotation result for phylogenetic analysis shows no common evolutionary gene in contigs, can i change the contigs?
I'm currently doing my thesis with the topic of phylogenetic analysis and is taking references from the previous person in my university (who have done the same topics but different species).
They used the largest contigs comprising well known markers such as matK, rbcL, and ITS. However my largest contigs doesn't have them.
Is it possible for me to use the 2nd/3rd largest contigs that have those markers?

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.

Please add some examples of these contigs; pictures would be good.

ITS: what is needed is the aim of the phylogenetic investigation. ITS evolves very quickly so it's good for a few closely related species and sub-species. Understanding the behaviour of a locus is as important as the volume of the locus analysed. Thus, the answer is yes.
I have a feeling these are insects/arthropods simply because of the fragmented nature of the genome: they're hard to assemble and genetically diverse so there's often no reference genome. The issue of genome vs locus was addressed by David Penny years ago: genome trees can be misleading and that would be justification for using the contigs to derive select loci. The reason why genome trees can be misleading is complex and a separate question.
Moreover, if you can map the loci to a chromosome and demonstrate the contigs are on different chromosomes you can use differences/similarities in trees between chromosomes as the first stage of disequilibrium analysis, i.e. selection acting via Mendelian exchange.
To reiterate there is good phylogenetic justification for switching contigs because the loci are more important than the size of the contig (you'll need supporting references however).


text:
installing and using differential methylation DSS software
I am attempting to install and use the Bioconductor package DSS on Ubuntu 22.04.
Following the instructions here: https://bioconductor.org/packages/release/bioc/html/DSS.html
I get to this error:
hfile_libcurl.c:45:10: fatal error: curl/curl.h: No such file or directory
   45 | #include <curl/curl.h>
      |          ^~~~~~~~~~~~~
compilation terminated.
make[1]: *** [Makefile.Rhtslib:128: hfile_libcurl.o] Error 1
make[1]: Leaving directory '/tmp/RtmpNry6uq/R.INSTALL10673560786c68/Rhtslib/src/htslib-1.7'
make: *** [Makevars.common:24: htslib] Error 2

In Ubuntu 22.04, install:
sudo apt install libcurl4-openssl-dev

Then reattempt:
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("DSS")

Then upon completion:
library(DSS)
Loading required package: Biobase
Loading required package: BiocGenerics

Attaching package: ‘BiocGenerics’

The following objects are masked from ‘package:stats’:

    IQR, mad, sd, var, xtabs
...


text:
Racon: empty overlap set!
Not sure if anyone can help me with a Racon issue as the developers are not very responsive on GitHub. I am trying to assemble a genome without a reference. I first did basecalling using dorado, then converted from bam to fastq using samtools, then trimmed using nanofilt from there I used flye for assembly and then ran minimap2 to get my sam file. For some reason when I then run Racon, I get the following error:
[racon::Polisher::initialize] loaded target sequences 0.000226 s
[racon::Polisher::initialize] loaded sequences 0.034854 s
[racon::Polisher::initialize] error: empty overlap set!

I used the following command to run the assembly:
flye --nano-raw sample_id.trimmed.fastq -o assembly -I 3 -t 4

I used the following command to run minimap2:
minimap2 -ax map-ont -t 4 sample_id.trimmed1.fastq assembly/assembly.fasta > sample_id.sam

I then used the following command to run racon:
racon -m 8 -x -8 -g -6 -w 500 -t 4 sample_id.trimmed1.fastq sample_id.sam assembly/assembly.fasta > Racon_polished.fasta

Could I please get some assistance on this issue as I am not sure what I am doing wrong. Thank you
First lines of the assembled fasta file:
>contig_1
TAACAGTGAAATTCAATGGTCTGGTATCCATCTAAACATTGAATGAGCAAGAAAAATACA

First lines of the trimmed fastq file:
@bf7e12fd-b52d-498e-9e61-39bfe2eacc80
TGTACCCTAGTGGAGAATTCACCAACGCAGCAATACCTCCGGTGTTCTGTCTACAGCGTTACTGTGTCT

First lines of the sam file:
contig_1    0   08ef4187-c2ab-4756-a269-275f6dfe99af    1   60  93179S96M1I81M1D29M1I10M2I420M1I118M1I107M1D63M1I18M4I660M2I242M CTAACTGTATTGCTGGGAAACAAATCAA


text:
VEP annotation INFO field Ensembl IDs and locations
I have a vcf file that I annoteted with VEP, for human data. I have run VEP to annotate my files with some additional parameters (as shown below in the ##VEP-command-line).
However, my output is rather strange (mainly the INFO column).
##VEP="v108" time="2023-04-27 15:13:08" cache="workflow/resources/variants/cache_vep/homo_sapiens/108_GRCh38" ensembl-funcgen=108.56bb136 ensembl-variation=108.a885ada ensembl-io=108.58d13c1 ensembl=108.d8a9c80 1000genomes="phase3" COSMIC="96" ClinVar="202205" HGMD-PUBLIC="20204" assembly="GRCh38.p13" dbSNP="154" gencode="GENCODE 42" genebuild="2014-07" gnomADe="r2.1.1" gnomADg="v3.1.2" polyphen="2.2.2" regbuild="1.0" sift="sift5.2.2"
##INFO=<ID=CSQ,Number=.,Type=String,Description="Consequence annotations from Ensembl VEP. Format: Allele|Consequence|IMPACT|SYMBOL|Gene|Feature_type|Feature|BIOTYPE|EXON|INTRON|HGVSc|HGVSp|cDNA_position|CDS_position|Protein_position|Amino_acids|Codons|Existing_variation|DISTANCE|STRAND|FLAGS|SYMBOL_SOURCE|HGNC_ID|GENE_PHENO|SIFT|PolyPhen|AF|AFR_AF|AMR_AF|EAS_AF|EUR_AF|SAS_AF|gnomADe_AF|gnomADe_AFR_AF|gnomADe_AMR_AF|gnomADe_ASJ_AF|gnomADe_EAS_AF|gnomADe_FIN_AF|gnomADe_NFE_AF|gnomADe_OTH_AF|gnomADe_SAS_AF|CLIN_SIG|SOMATIC|PHENO">
##VEP-command-line='vep --dir_cache workflow/resources/variants/cache_vep --vcf -i WESP_run/results/variant_calling/VP_DBS_25_7_S7.vcf -o WESP_run/results/variant_annotation/VP_DBS_25_7_S7.annotated.vcf.gz --af --af_1kg --af_gnomade --sift b --polyphen b --gene_phenotype --fork 24 --offline --compress_output gzip'
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  VP_DBS_25_7_S7
chr1    13273   .   G   C   25.3    PASS    CSQ=C|non_coding_transcript_exon_variant|MODIFIER|DDX11L1|ENSG00000223972|Transcript|ENST00000450305|transcribed_unprocessed_pseudogene|5/6||||313|||||rs531730856||1||HGNC|HGNC:37102||||0.0950|0.0204|0.1455|0.0625|0.1471|0.1401||||||||||||,C|non_coding_transcript_exon_variant|MODIFIER|DDX11L2|ENSG00000290825|Transcript|ENST00000456328|lncRNA|3/3||||521|||||rs531730856||1||EntrezGene|||||0.0950|0.0204|0.1455|0.0625|0.1471|0.1401||||||||||||,C|downstream_gene_variant|MODIFIER|WASH7P|ENSG00000227232|Transcript|ENST00000488147|unprocessed_pseudogene||||||||||rs531730856|1131|-1||HGNC|HGNC:38034||||0.0950|0.0204|0.1455|0.0625|0.1471|0.1401||||||||||||,C|downstream_gene_variant|MODIFIER|MIR6859-1|ENSG00000278267|Transcript|ENST00000619216|miRNA||||||||||rs531730856|4096|-1||HGNC|HGNC:50039||||0.0950|0.0204|0.1455|0.0625|0.1471|0.1401|||||||||||| GT:GQ:DP:AD:VAF:PL  0/1:23:55:28,27:0.490909:25,0,27

For instance, at the first location/variant I look at the Ensembl gene id and see multiple genes on the same locations and strand. This is very strange to me, as I would understand different transcripts, but different genes on the same strand? Also, apart from that, the corresponding locations/regions of genes on the Ensembl website do not agree with the location in the vcf file.
For instance, the location for this SNP is 13273, but the annotations are (corresponding ensembl archive version 108):

Agrees: http://oct2022.archive.ensembl.org/Homo_sapiens/Gene/Summary?db=core;g=ENSG00000223972;r=1:12010-13670;t=ENST00000450305
Agrees: http://oct2022.archive.ensembl.org/Homo_sapiens/Gene/Summary?db=core;g=ENSG00000290825;r=1:11869-14409;t=ENST00000456328
Does not agree: http://oct2022.archive.ensembl.org/Homo_sapiens/Gene/Summary?db=core;g=ENSG00000227232;r=1:14404-29570;t=ENST00000488147
Does not agree: http://oct2022.archive.ensembl.org/Homo_sapiens/Gene/Summary?db=core;g=ENSG00000278267;r=1:17369-17436;t=ENST00000619216

What am I missing here?
Kind regards,
Eva

You need to set your expectations clearly instead of just saying "Does not agree". I'm taking a guess here - your statement is that the 2/4 of the gene co-ordinates do not overlap the variant location.
This is because you're viewing just part of the picture - the co-ordinates - and not seeing another important part, which is the kind of effect. For the latter two genes, VEP clearly says the effect is downstream_gene_variant, which Sequence Ontology defines as "A sequence variant located 3' of a gene."
Putting those two pieces together, the entire picture makes sense. Variants affect different transcripts (of possibly different genes) in different manners. Looking at just 1 of those 2 components will of course make things seem odd.

You are right, it is just overlapping genes. And VEP annotates these, while other variant annotation tools like ANNOVAR do not (just outputs the most deleterious) ref.

That’s not true. Both tools can annotate all or most deleterious effects depending on the parameters you use.

Just saw the reference post you’d added - that post is wrong. See ANNOVAR documentation- it says “ If the users want to have all functional consequences printed out (rather than just the most important one defined by the precedence above), the --separate argument should be used. In this case, several output lines may be present for each variant, representing several possible functional consequences.”

I guess the difference would be the default behavior maybe, but you are right. Sorry for the confusion and thank you for your patience.


text:
How to run a GATK Docker Image with local files?
I'm trying to use the HaplotypeCaller from the GATK toolkit but I keep getting an error. I pulled GATK through Docker and am using this command:
docker run -v /Users/rimo/ -it broadinstitute/gatk:latest gatk HaplotypeCaller -R /Users/rimo/reference.fasta -I /Users/rimo/sample1.bam -O /Users/rimo/sample1.g.vcf.gz -ERC GVCF

/Users/rimo is my home directory it's where the files are currently located on my computer and it's where I want GATK to look for files to use in the command.
but then I get this error:
A USER ERROR has occurred: The specified fasta file (file:///Users/rimo/reference.fasta) does not exist.

I don't fully understand how Docker commands work but I found a similar problem listed with this issue. I tried to use their solution but it's still giving me the same error as listed above.

Read the manual please: https://docs.docker.com/storage/bind-mounts/

-v or --volume: Consists of three fields, separated by colon characters (:). The fields must be in the correct order, and the meaning of each field is not immediately obvious.

In the case of bind mounts, the first field is the path to the file or directory on the host machine.
The second field is the path where the file or directory is mounted in the container.
The third field is optional, and is a comma-separated list of options, such as ro, z, and Z. These options are discussed below.

Try
docker run -v /Users/rimo/:/myhome ...

Also, why are you combining a -it run with an actual command? I thought the -it is to start an interactive shell, which is the opposite of running a single command.

Thanks for your response! Didn't see it until I posted my answer below but it definitely makes sense. I'm still trying to figure out how Docker works so your comment about -it and what it actually means helped!

I went back to the documentation listed on the GATK website instead of Docker because the Docker documentation was a bit confusing. And was able to figure it out with this command:
docker run -v ~/data:/gatk/my_data -it broadinstitute/gatk:latest

And for those who might use this as a resource because they don't understand how to use local files on a Docker image the breakdown is as follows:

~/data = path where the data is located on the local machine
/gatk/my_data = path where the data will be located on the Docker image (folder will be created when the command is executed)
this link saved me: https://gatk.broadinstitute.org/hc/en-us/articles/360035889991

Upvoting @RamRS would be a cool thing to do.

The main issue here was that I didn't understand the breakdown of a docker run command and the documentation on the website is pretty confusing if you aren't familiar with all of the verbiage. Going back and breaking down the command helped me understand how the local files on my computer get onto the docker image and where to find them in order to run a separate command. I got there in the end but thanks for your contribution!


text:
Mendelian randomization
I would be thankful to you if you can help with how I can use MR from my case-control study to link DNA methylation (epic array) with toxic element exposure (arsenic) and health outcome (Cardiovascular diseases).
I do have data from 125 Cardiovascular diseases and 125 non-Cardiovascular diseases individuals matched with age, sex, and smoking.I do have phenotypes files in a csv format having sex, age, smoking, disease, samples id, iAs (arsenic) values in log format and then a separate files for methylation data. we also have SNP data for these samples.
Many thanks,

Please add more detail to this question. What format are your datasets in? Can you provide any examples of that format?

added other informations, hope this is ok now.

csv and log formats are not well specified. Please show some examples of inputs and expected outputs.

thanks gringer, I don't have whole data yet but thinking how to do this analysis step by step.

I understand linkage disequilibrium (LD) versus panmixia, which the technical term for Mendelian randomisation and for humans PLINK is a population package. The best outcome is a correlation between LD and cardiovascular disease. It is however an unlikely outcome.
This really is a machine learning question and you would need to understand "feature selection" and will require very high accuracy scores >0.95 - however you could triple the size of our control group quite easily (recommended). This is a separate question and I'm happy to answer this there.

Just to address the comments. Inheritance probability, if that is the intended statistics, I personally don't think it would work because a 125 sample is low and the patients would need to be in cohorts of familial generations. LD is unlikely the explanation with humans but other non-humans species that might work (selection advantage).
My personal view is this is a straight machine learning calculation, but it's a separate question. If your sample is heavily stratified across each familial generation then Hardy-Weinberg might work, but my experience of patient cohorts is that is really unusual and there could still be confounders. Supervised learning (subset of ML) is 100% the way forward.
.... possibly use GWAS, but you've multiple features/training targets. GWAS is usually when there's one phenotype, ML would account for interaction in a way GWAS could not (it couldn't account for non-sequence features, i.e. mixed data).

Thanks, I would like to know what step by step analysis I need to perform here; like first I can look for arsenic exposure and CpG association using linear regression and then I can also look for CpG and health outcome association and then for common CpG site for arsenic exposure and health outcome I need to take for identify with the genetics variants(IVs) that affect the methylation at these site. And finally not sure if this is the correct  approach and then what I need as input for MR.

Thanks @bioinfonext Upvotes and accepts are more appreciated than thanks. I've answered above. An alternative approach e.g. via ML is a separate question.

Thanks again, still I am not clear which R packages should I use for these analysis.

GWAS: statgenGWAS package; machine learning its caret. LD is about PLINK and other autonomous packages. If you want a single package for an easy answer then shoehorn GWAS

thanks, so first I need to find CpG sites associated with arsenic exposure and then I should use these CpG site for GWAS and should also perform GWAS with disease phenotypes (CVD/noCVD). and then use these GWAS for MR?

Or do you think here I should use structural equation modelling to find CpG sites associated with arsenic exposure and disease condition as I have very small size samples and you already point out that it may not useful for MR analysis.

@bioinfonext I think its a separate question. Small sample size changes things a lot because that means no ML and no HW. I've never use SEM; I know what it does. Do consider upvoting however.


text:
How many Ns and ns in GRCh37 / GRCh38 per 'canonical' chromosome?
This is kind of pedantic, but I'm not sure where to look...
For GRCh38 (and a lot of work...) I have the following...
Chr          Length         Ns   ns
chr1    248,956,422 18,475,229  181
chr2    242,193,529   1645,291   10
chr3    198,295,559    195,420    4
chr4    190,214,555    461,888    0
chr5    181,538,259    272,881    0
chr6    170,805,979    727,255  202
chr7    159,345,973    375,841    1
chr8    145,138,636    370,500    0
chr9    138,394,717 16,604,127   40
chr10   133,797,422    534,287  173
chr11   135,086,622    552,880    0
chr12   133,275,309    137,270  223
chr13   114,364,328 163,81,200    3
chr14   107,043,718 164,75,569    0
chr15   101,991,189 173,49,841   23
chr16    90,338,345   8532,401    1
chr17    83,257,441    336,367  870
chr18    80,373,285    283,627   53
chr19    58,617,616    176,858    0
chr20    64,444,167    499,586  324
chr21    46,709,983  6,621,300   64
chr22    50,818,468 11,658,691    0
chrX    156,040,895  1,147,765  101
chrY     57,227,415 30,812,232  140
chrM         16,569          1    0

(Where) could I have looked this up?
What are the counts for GRCh37?

I don't know where you could look this up and of course this will depend on the versions of the relevant genomes and on which exact assemblies you are working with. I doubt you would have found this specific number, how many Ns anywhere, but it's really not very hard to get them. For instance, I have hg19 and hg38 versions locally:
$ for genome in hg19 hg38 ; do ls $genome/$genome.fa; done
hg19/hg19.fa
hg38/hg38.fa

This means I can get the Ns for both using:
for genome in hg19 hg38 ; do
  FastaToTbl "$genome"/$genome.fa |
    perl -ne 'print if /^(chr)?[0-9MXY]*\s/' |
      awk -F"\t" '{ print $1, split($2,a,/[nN]/)}' > $genome.ns
done

And then combine into a nice table with:
$ printf 'chr\thg19\thg38\n-----------------------\n';
  join -j1 -t ' ' <( sed 's/chr//' hg19 | sort -gk1,1 ) \
                  <( sed 's/chr//'  hg38 | sort -gk1,1) |
                   column -t | sed 's/^/chr/'
chr hg19    hg38
-----------------------
chrM   1         2
chrX   4170001   1147867
chrY   33720001  30812373
chr1   23970001  18475411
chr2   4994856   1645302
chr3   3225296   195425
chr4   3492601   461889
chr5   3220001   272882
chr6   3720002   727458
chr7   3785001   375843
chr8   3475101   370501
chr9   21070001  16604168
chr10  4220010   534461
chr11  3877001   552881
chr12  3370503   137494
chr13  19580001  16381204
chr14  19060001  16475570
chr15  20836627  17349865
chr16  11470001  8532403
chr17  3400001   337238
chr18  3420020   283681
chr19  3320001   176859
chr20  3520001   499911
chr21  13023254  6621365
chr22  16410022  11658692

The fancy sed stuff is just to ensure that sorting works correctly by removing and then re-adding the chr string. Note that this is counting both N and n. The FastaToTbl script can be found here.


text:
Filtering columns from large file based on a column list
I have list of samples as a txt file and want to only keep these samples from this large file; column name is samples name in the large file.
library(readr)
library(data.table)

df1 <- read_csv("control.csv", col_types = cols(Chromosome = col_skip(), Start = col_skip(), End = col_skip(), Strand = col_skip()))

df1[1:5,1:5]
# A tibble: 5 × 5
  ID          DN54 DC146 DDN44 DC179
  <chr>      <dbl> <dbl> <dbl> <dbl>
1 cg26928153 0.848 0.879 0.872 0.897
2 cg16269199 0.731 0.732 0.759 0.775
3 cg13869341 0.842 0.863 0.819 0.896
4 cg24669183 0.824 0.866 0.822 0.910
5 cg26679879 0.425 0.363 0.420 0.358

What have you tried? This is a straightforward R operation.

it is resolved now, I have updated my code above.

"I have updated my code above" - Don't do that, add an answer and accept it instead.

Thanks Ram, I will add correct code as answer soon.

You can save time and memory by preprocessing your file before loading it into R.
Use miller toolkit to manipulate csv or tab-delimited files using column names.
Below is a one-liner that reads csv file, and writes columns specified in option -f, in the order (-o) listed on the command line:
mlr --csv cut -o -f col3,col1 in.csv > out.csv

To read the list of columns from a file, 1 column name per line:
mlr --csv cut -o -f $(paste -s -d, columns.txt) in.csv > out.csv

See also:
Note that miller can be easily installed, for example, using conda, specifically miniconda, like so:
conda create --channel conda-forge --name miller miller

It's really just count[,col$V1]. Please learn about basic subsetting of columns and rows, it's simple yet powerful.

thanks but it is not giving me CpG ID in first column.

There you see that you're not investing time into understanding an answer. What this here does is to provide names to the object to be filtered using the inbuilt R filtering operator. If you need the names plus something else then add this. Like c(cpg_name, col$V1). Show some effort.

Thanks @ATpoint, it is sorted now and with your command, I was also not getting column names.

Below code works perfectly of me.
library(readr)
library(data.table)

df1 <- read_csv("control.csv", col_types = cols(Chromosome = col_skip(), Start = col_skip(), End = col_skip(), Strand = col_skip()))

df1[1:5,1:5]
# A tibble: 5 × 5
  ID          DN54 DC146 DDN44 DC179
  <chr>      <dbl> <dbl> <dbl> <dbl>
1 cg26928153 0.848 0.879 0.872 0.897
2 cg16269199 0.731 0.732 0.759 0.775
3 cg13869341 0.842 0.863 0.819 0.896
4 cg24669183 0.824 0.866 0.822 0.910
5 cg26679879 0.425 0.363 0.420 0.358

col <- unlist(fread("samples.csv", header=F))
col <- c("ID", col)
extracted <- df1[,col, with=F]

write.table(extracted, file="control.beta.txt",sep=",", quote=F, col.names=TRUE,row.names=FALSE)

df <- read_csv("control.beta.txt")
df[1:5,1:5]
# A tibble: 5 × 5
  ID         DC103R DC110R DC115 DC120R
  <chr>       <dbl>  <dbl> <dbl>  <dbl>
1 cg26928153  0.929  0.885 0.862  0.877
2 cg16269199  0.836  0.826 0.604  0.739
3 cg13869341  0.826  0.903 0.830  0.841
4 cg24669183  0.767  0.719 0.808  0.812
5 cg26679879  0.388  0.366 0.428  0.360


text:
Metabolic Modeling
I am trying to identify the flux of IgG antibody in CHO ovary cells with optimised biomass and glucose input ranging from 0 to 10, but I am not sure when I run the FBA analysis what should I define my lower and upper bounds for glucose.

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

Hi, this question is not detailed enough to get a good answer. Please include information about the datasets that you have. What format are they in? Can you show a few lines of example input, and ideal output?


text:
Unable to access error logs of failed cromwell execution (WDL) in local machine running in ubuntu environment
I was trying to run bases2fastq software downloaded from the element bioscience platform. I was running this software through WDL workflow named list_fastq. For some reason, maybe syntax, my cromwell execution failed. However, I haven't been able to read the logs as they are empty. Below is the content of the 'stderr.background' file.
mkfifo: cannot create fifo '/mnt/d/element-bioscience/cromwell-executions/list_fastq/5e9a780d-c010-41f7-bba0-bcb42c870891/call-list_samples/tmp.05c20794/out.78994': Operation not supported
mkfifo: cannot create fifo '/mnt/d/element-bioscience/cromwell-executions/list_fastq/5e9a780d-c010-41f7-bba0-bcb42c870891/call-list_samples/tmp.05c20794/err.78994': Operation not supported
/mnt/d/element-bioscience/cromwell-executions/list_fastq/5e9a780d-c010-41f7-bba0-bcb42c870891/call-list_samples/execution/script: line 17: /mnt/d/element-bioscience/cromwell-executions/list_fastq/5e9a780d-c010-41f7-bba0-bcb42c870891/call-list_samples/tmp.05c20794/out.78994: No such file or directory
/mnt/d/element-bioscience/cromwell-executions/list_fastq/5e9a780d-c010-41f7-bba0-bcb42c870891/call-list_samples/execution/script: line 18: /mnt/d/element-bioscience/cromwell-executions/list_fastq/5e9a780d-c010-41f7-bba0-bcb42c870891/call-list_samples/tmp.05c20794/err.78994: No such file or directory

Please guide me to the steps on being able to prevent cromwell from deleting the contents of log files, if that's what's happening.


text:
PyRanges: drop specific record from an object
Given a specific PyRanges object, such as
import pyranges as pr
gr = pr.from_dict({"Chromosome": [1, 1, 1, 1], "Start": [1, 25000, 30000, 40000],
                   "End": [24000, 29500, 37000, 45000], "ratio": [0.5, 0.3, 0.5, 0.4],
                   "ID": (0,1,2,3), "Strand": "+"})

gr

+--------------+-----------+-----------+-------------+-----------+--------------+
|   Chromosome |     Start |       End |       ratio |        ID | Strand       |
|   (category) |   (int64) |   (int64) |   (float64) |   (int64) | (category)   |
|--------------+-----------+-----------+-------------+-----------+--------------|
|            1 |         1 |     24000 |         0.5 |         0 | +            |
|            1 |     25000 |     29500 |         0.3 |         1 | +            |
|            1 |     30000 |     37000 |         0.5 |         2 | +            |
|            1 |     40000 |     45000 |         0.4 |         3 | +            |
+--------------+-----------+-----------+-------------+-----------+--------------+

I want to get rid of one interval here (not necessarily the first or the last: for this exercise, for example the third one, without relying on the presence of the ID column.  So I would need to drop the record with Chromosome 1, Start 30000 and End 37000. The use case is that I need to replace it with a different record.
Replacement of the record is somewhat expensive but not a problem (pr.concat followed by a sort to put everything in order), but I can't figure, after reading docstrings and searching around, how to remove a record.
I can only think of making the record I want to drop a second PyRanges object then using overlap() with invert=True to remove the record. Is there a more elegant and less convoluted way to do this?


text:
What is the best Query to retrieve DNA from NCBI?
I want to retrieve a sequence for many species from the Nucleotide database in NCBI.
I'm using a command line approach and I have to figure out what is the best query that will return exclusively the region of DNA I am interested in and filter out unwanted noisy sequences.
I am using a wide range of species (that I have stored as a list of TaxIDs of about 2000 species): it includes small crustaceans, invertebrates, some algae and small vertebrates (reptiles, amphibians and few mammals).
My final goal is to obtain a phylogenetic tree for all or most of the species.
I have been suggested to use these genes to create such phylogeny

CO1 mithocondrial DNA
16S rRNA
18S rRNA

I want to formulate a specific query that will return exclusively those sequences.
I'm using the GenBank query builder to visually check for the accuracy of my search and when I find a good query I will use it in the API.
So far I came up with the following queries:

(COX1[Title] OR CO1[Title]) AND complete[Title] --> 63/2000 Species
16S[title] AND complete[title] AND rRNA[title] NOT partial[title] --> 9/2000 species
18S[title] AND complete[title] AND rRNA[title] NOT partial[title] --> 15/2000 species

As you can see the number of species that I get sequence for is very low compared to the initial 2000 species. I doubt that we have so little available sequences (especially for COX1 that is used for barcoding)
Can you help me understanding whether my queries are good or not? And if possible suggest a better alternative
More Info
A subset of 10 of my species of interest is
Rasbora heteromorpha
Elasmopus rapax
Gasterosteus aculeatus
Palaemonetes pugio
Catostomus commersoni
Daphnia magna
Oryzias latipes
Xenopus laevis
Tigriopus japonicus
Oncorhynchus mykiss

Out of these 10 only 2 have sequences available using the COX1 query

Name
SeqID

Oncorhynchus mykiss
EU186789.1

Xenopus laevis
AB278691.1

But Daphnia magna is one of the most commonly used organisms in lab and I found this paper regarding the complete mitogenome of a specific strain of D. magna thus implying that the complete mitogenome of the species is already known. That means that there must be a way to retrieve all the other species that have full mitogenome but no COX1 gene mentioned in their GenBank title.

What database are you searching? I get 1995 hits for 16S[title] AND complete[title] AND rRNA[title] NOT partial[title], across 1430 eukaryotes, 287 bacteria and 178 archaea when using "nucleotide".

Yes sorry maybe i was not very clear: the numbers on the right are the species that are both in my database and also in the results from Nucleotide. So even if you get 1400 eukariotes only 9 of these are overlapping with the species i work with. But they are mostly lab species and not rare, That is why I find it strange that i don't get many sequences from my search

You seem to be assuming gene names are shared across organisms. It's a reasonable assumption but, sadly, it is absolutely false. Gene names are not shared across species. They can be, sometimes, but not always or even often. For example, there is no CO1 or COX1 in human. "Cytochrome c oxidase assembly factor 1" has the symbol "COA1".

For anyone having the same type of problem as me (batch COI download) i recommend reading this article as it has proven very useful for me DOI: 10.1111/1755-0998.13756

The search term you need is [gene name] rather than [title]. For example,
(Gasterosteus aculeatus[Organism] AND cox1[Gene Name]) OR (Gasterosteus aculeatus[Organism] AND coi[Gene Name])

Results in 266 hits which all contain COX1 COI
Using Daphnia magna ...
(Daphnia magna[Organism] AND cox1[Gene Name]) OR (Daphnia magna[Organism] AND coi[Gene Name])

This results in 775 hits.
Daphina magna comprises 61 mtDNA genomes, you perform the above search and sort by sequence length and set the sequences per page to 200(easy).
Hope that resolves the issue.

There are two main problems I can see in your approach. First, gene names are not standardized across species and the "same" gene (what that means is a whole different discussion) can have different names. For example, human  "Cytochrome c oxidase assembly factor 1" has the symbol "COA1".
Next, I don't quite understand why you expect the word "complete" to be in the title. That is actually why you don't find the D. magna sequence: it isn't complete, the start codon hasn't been correctly annotated as you can see on the relevant NCBI gene page:

I would instead try something like this:
((chrM[Chromosome]) OR (chromosome M[Chromosome]) OR (MT[Chromosome]) OR (chromosome MT[Chromosome])) AND ((CO1[Gene Name]) OR (COX1[Gene Name]) OR ((cytochrome[Gene Name]) AND (oxidase[Gene Name]))

or
((16S[Gene/Protein Name]) and (ribosomal[Gene/Protein Name]))

All that said, trying to rely on gene annotations is rarely the best approach. There just isn't enough standardization. I would instead look at databases that provide lists of homologous genes. For instance, NCBI's Clusters of Orthologous Genes (COGs), which provides grouped homologous sequences.

Thanks for the tips they are very relevant. But the COG is only a viable option for archea and bacteria so it doesn't suit for me.

Ah, crud, sorry @Mirko, I didn't know it was limited to microorganisms. There is also KOG and Ortho MCL or you can use BioMart to extract orthologs etc. My main point here is that trying to use a gene name is unlikely to work well and you should look for dedicated resources instead.


text:
How to add charges to a PDB file?
I did some calculations using a PDB file as input. As a result, I got a file with the atomic partial charges (only the values in one column, no other information) and the output PDB file structure without the charge information.
My question: Is there a tool to add the charges to the output PDB file using the separate file?
My PDB file:
    TITLE     sistema05
    REMARK   4      COMPLIES WITH FORMAT V. 3.0, 1-DEC-2006
    REMARK 888
    REMARK 888 WRITTEN BY MAESTRO (A PRODUCT OF SCHRODINGER, LLC)
    HETATM    1  C1  BAS A   1      10.151  -3.071  -1.881  1.00  0.00           C
    HETATM    2  C2  BAS A   1       9.389  -3.683  -0.874  1.00  0.00           C
    HETATM    3  C3  BAS A   1       9.976  -3.936   0.382  1.00  0.00           C
    HETATM    4  C4  BAS A   1       9.194  -3.742   1.530  1.00  0.00           C
    HETATM    5  C5  BAS A   1       9.781  -3.201   2.695  1.00  0.00           C
    HETATM    6  C6  BAS A   1       9.029  -2.357   3.510  1.00  0.00           C
    HETATM    7  C7  BAS A   1       9.651  -1.313   4.236  1.00  0.00           C
    HETATM    8  C8  BAS A   1       8.948  -0.139   4.488  1.00  0.00           C
    HETATM    9  C9  BAS A   1       9.635   1.100   4.502  1.00  0.00           C
    HETATM   10  C10 BAS A   1       9.001   2.259   4.094  1.00  0.00           C
    HETATM   11  C11 BAS A   1       9.724   3.177   3.329  1.00  0.00           C
    HETATM   12  C12 BAS A   1       9.148   3.815   2.228  1.00  0.00           C
    HETATM   13  C13 BAS A   1       9.909   3.945   1.081  1.00  0.00           C
    HETATM   14  C14 BAS A   1       9.332   3.754  -0.198  1.00  0.00           C
    HETATM   15  C15 BAS A   1      10.098   3.266  -1.245  1.00  0.00           C
    HETATM   16  C16 BAS A   1       9.501   2.445  -2.232  1.00  0.00           C
    HETATM   17  C17 BAS A   1      10.234   1.431  -2.847  1.00  0.00           C
    HETATM   18  C18 BAS A   1       9.575   0.238  -3.215  1.00  0.00           C
    HETATM   19  C19 BAS A   1      10.249  -0.986  -3.121  1.00  0.00           C
    HETATM   20  C20 BAS A   1       9.529  -2.124  -2.713  1.00  0.00           C
...
...
...
    CONECT    1  212    2   20
    CONECT    1   20
    CONECT    2    1    3   21
    CONECT    2   21
    CONECT    3    2  213    4
    CONECT    3    4
    CONECT    4    3    5   23
    CONECT    4    3
    CONECT    5    4  214    6
    CONECT    5    6
    CONECT    6    5    7   25
    CONECT    6    5
    CONECT    7    6  215    8
    CONECT    7    8
    CONECT    8    7    9   27
    CONECT    8    7
    CONECT    9    8  216   10
    CONECT   10    9   11   29
    CONECT   11   10  217   12
    CONECT   12   11   13   31
    CONECT   13   12  218   14
    CONECT   14   13   15   33
    CONECT   14   15
    CONECT   15   14  219   16
    CONECT   15   14
    CONECT   16   15   17   35
    CONECT   16   17
    CONECT   17   16  220   18
    CONECT   17   16
    CONECT   18   17   19   37
    CONECT   18   19
    CONECT   19   18  221   20
    CONECT   19   18
    CONECT   20    1   19   39
    CONECT   20    1
    CONECT   21    2   22   40
    END

My CHARGE file:
   -0.02474906
    0.02291960
   -0.02789302
    0.02227110
   -0.02831805
    0.02230034
   -0.02876149
    0.02225527
   -0.02774076
    0.02216971
   -0.02074528
    0.02275309
   -0.02662750
    0.02177531
   -0.03007329
    0.02283631
   -0.02779055
    0.02312780
   -0.02312680
    0.02315866

You could use the b factor or occupancy value column to store your charge values, dont remember how many digits are available for each of the 2 proprierties. Parse a pdb file with biopython and change it with your charge files, I assume the pdb doesnt have any header or remark just atoms from lines 1 onwards that relate to the list of charges one per line in the second file ?? A small snippet of both files would be of help for anybody interested in helping you out

See https://stackoverflow.com/questions/53938256/how-to-modify-columns-contents-in-a-protein-pdb-file-by-keys-and-values-stored-i

Sorry just realized pdb has column for that https://www.cgl.ucsf.edu/chimera/docs/UsersGuide/tutorials/pdbintro.html , column 79-80

Could use biopython https://biopython.org/docs/1.79/api/Bio.PDB.Atom.html and set the charge for each atom   atom.set_charge( CHARGE) should do

Thanks @pippo1980. I edited the post, adding part of the files. I will go through your suggestions.

Mmmh your charges are a lot of digits... . You'll need a custom made pdb thats not a pdb format anymore to store all that and a parser to be able to reload that pdb, are you sure there are not other more suitable formats

OK given , input_connect.pdb , using Biopython library
TITLE     sistema05
REMARK   4      COMPLIES WITH FORMAT V. 3.0, 1-DEC-2006
REMARK 888
REMARK 888 WRITTEN BY MAESTRO (A PRODUCT OF SCHRODINGER, LLC)
HETATM    1  C1  BAS A   1      10.151  -3.071  -1.881  1.00  0.00           C
HETATM    2  C2  BAS A   1       9.389  -3.683  -0.874  1.00  0.00           C
HETATM    3  C3  BAS A   1       9.976  -3.936   0.382  1.00  0.00           C
HETATM    4  C4  BAS A   1       9.194  -3.742   1.530  1.00  0.00           C
HETATM    5  C5  BAS A   1       9.781  -3.201   2.695  1.00  0.00           C
HETATM    6  C6  BAS A   1       9.029  -2.357   3.510  1.00  0.00           C
HETATM    7  C7  BAS A   1       9.651  -1.313   4.236  1.00  0.00           C
HETATM    8  C8  BAS A   1       8.948  -0.139   4.488  1.00  0.00           C
HETATM    9  C9  BAS A   1       9.635   1.100   4.502  1.00  0.00           C
HETATM   10  C10 BAS A   1       9.001   2.259   4.094  1.00  0.00           C
HETATM   11  C11 BAS A   1       9.724   3.177   3.329  1.00  0.00           C
HETATM   12  C12 BAS A   1       9.148   3.815   2.228  1.00  0.00           C
HETATM   13  C13 BAS A   1       9.909   3.945   1.081  1.00  0.00           C
HETATM   14  C14 BAS A   1       9.332   3.754  -0.198  1.00  0.00           C
HETATM   15  C15 BAS A   1      10.098   3.266  -1.245  1.00  0.00           C
HETATM   16  C16 BAS A   1       9.501   2.445  -2.232  1.00  0.00           C
HETATM   17  C17 BAS A   1      10.234   1.431  -2.847  1.00  0.00           C
HETATM   18  C18 BAS A   1       9.575   0.238  -3.215  1.00  0.00           C
HETATM   19  C19 BAS A   1      10.249  -0.986  -3.121  1.00  0.00           C
HETATM   20  C20 BAS A   1       9.529  -2.124  -2.713  1.00  0.00           C
CONECT    1  212    2   20
CONECT    1   20
CONECT    2    1    3   21
CONECT    2   21
CONECT    3    2  213    4
CONECT    3    4
CONECT    4    3    5   23
CONECT    4    3
CONECT    5    4  214    6
CONECT    5    6
CONECT    6    5    7   25
CONECT    6    5
CONECT    7    6  215    8
CONECT    7    8
CONECT    8    7    9   27
CONECT    8    7
CONECT    9    8  216   10
CONECT   10    9   11   29
CONECT   11   10  217   12
CONECT   12   11   13   31
CONECT   13   12  218   14
CONECT   14   13   15   33
CONECT   14   15
CONECT   15   14  219   16
CONECT   15   14
CONECT   16   15   17   35
CONECT   16   17
CONECT   17   16  220   18
CONECT   17   16
CONECT   18   17   19   37
CONECT   18   19
CONECT   19   18  221   20
CONECT   19   18
CONECT   20    1   19   39
CONECT   21    2   22   40

and charges list input_charge.txt :
-0.02474906
 0.02291960
-0.02789302
 0.02227110
-0.02831805
 0.02230034
-0.02876149
 0.02225527
-0.02774076
 0.02216971
-0.02074528
 0.02275309
-0.02662750
 0.02177531
-0.03007329
 0.02283631
-0.02779055
 0.02312780
-0.02312680
 0.02315866

with code :
import Bio
print('########## Biopython  VERSION ##########################################')
print('\n\n BIOPYTHON_VERSION : ', Bio.__version__)
print('###################################################################')

from Bio.PDB import (
                    PDBParser,
                    Atom,

                    )

                    # PDBIO)

import PDBIOcopy

# from Bio.PDB import  PDBIO

pdb = 'input'

charges = 'input_charge.txt'

charges_list = []

with open(charges , 'r') as handle :

    for i in handle.readlines() :

        charges_list.append(i.strip())

print(charges_list)

parser = PDBParser(PERMISSIVE=True , QUIET = False)

structure = parser.get_structure(pdb, pdb+'.pdb')

atoms = structure.get_atoms()

idx = 0

for atom in atoms :

    print(atom)

    atom.set_charge(charges_list[idx])

    idx += 1

atoms = structure.get_atoms()

for atom in atoms :

    print(atom , atom.get_charge())

struct_to_save = PDBIOcopy.PDBIO()

# struct_to_save = PDBIO()

struct_to_save.set_structure(structure)

struct_to_save.save('result_input.pdb')

and modifying Bio.PDB.PDBIO as :
line 22 :
_ATOM_FORMAT_STRING = (
    "%s%5i %-4s%c%3s %c%4i%c   %8.3f%8.3f%8.3f%s%6.2f      %4s%2s%12s\n"
)

line 164 :
# charge="  ", #deleted
pqr_charge     # insertes

line 213:
charge = atom.pqr_charge  #inserted

line 389:
atom.pqr_charge,  #inserted

output result_input.pdb:
HETATM    1  C1  BAS A   1      10.151  -3.071  -1.881  1.00  0.00           C -0.02474906
HETATM    2  C2  BAS A   1       9.389  -3.683  -0.874  1.00  0.00           C  0.02291960
HETATM    3  C3  BAS A   1       9.976  -3.936   0.382  1.00  0.00           C -0.02789302
HETATM    4  C4  BAS A   1       9.194  -3.742   1.530  1.00  0.00           C  0.02227110
HETATM    5  C5  BAS A   1       9.781  -3.201   2.695  1.00  0.00           C -0.02831805
HETATM    6  C6  BAS A   1       9.029  -2.357   3.510  1.00  0.00           C  0.02230034
HETATM    7  C7  BAS A   1       9.651  -1.313   4.236  1.00  0.00           C -0.02876149
HETATM    8  C8  BAS A   1       8.948  -0.139   4.488  1.00  0.00           C  0.02225527
HETATM    9  C9  BAS A   1       9.635   1.100   4.502  1.00  0.00           C -0.02774076
HETATM   10  C10 BAS A   1       9.001   2.259   4.094  1.00  0.00           C  0.02216971
HETATM   11  C11 BAS A   1       9.724   3.177   3.329  1.00  0.00           C -0.02074528
HETATM   12  C12 BAS A   1       9.148   3.815   2.228  1.00  0.00           C  0.02275309
HETATM   13  C13 BAS A   1       9.909   3.945   1.081  1.00  0.00           C -0.02662750
HETATM   14  C14 BAS A   1       9.332   3.754  -0.198  1.00  0.00           C  0.02177531
HETATM   15  C15 BAS A   1      10.098   3.266  -1.245  1.00  0.00           C -0.03007329
HETATM   16  C16 BAS A   1       9.501   2.445  -2.232  1.00  0.00           C  0.02283631
HETATM   17  C17 BAS A   1      10.234   1.431  -2.847  1.00  0.00           C -0.02779055
HETATM   18  C18 BAS A   1       9.575   0.238  -3.215  1.00  0.00           C  0.02312780
HETATM   19  C19 BAS A   1      10.249  -0.986  -3.121  1.00  0.00           C -0.02312680
HETATM   20  C20 BAS A   1       9.529  -2.124  -2.713  1.00  0.00           C  0.02315866
TER      21      BAS A   1
END

That unfortunately is not a pdb file format anymore:
PDB Format uses for RECORD ATOM CHARGE columns 79-80 , see Protein Data Bank Contents Guide: Atomic Coordinate Entry Format Description
At the end of all of this I realized that using python and just adding the charge values at the end of the pdb file read line by line would be a lot easier & faster but anyway I had the chance to explore PDBIO
If you want to add the charges in the bfactor or occupancy columns just use :atom.set_bfactor(charges_list[idx]) or atom.set_occupancy(charges_list[idx]) be sure to modify PDBIO to make space for your values
There are other examples on SO about same problem too :
How to modify columns contents in a protein PDB file by keys and values stored in a text file using only bash tools
Just Python code :
pdb = 'input_connect.pdb'

charges = 'input_charge.txt'

charges_list = []

with open(charges , 'r') as handle :

    for i in handle.readlines() :

        charges_list.append(i.strip())

print(charges_list)

with open('input_connect.pdb' , 'r') as handle :

    pdb_file = handle.readlines()

print(pdb_file)

idx = 0

with open('result_python.pdb' , 'w') as handle :

    for i in pdb_file :

        i = i.strip()

        if i[0:6] == 'HETATM' :

            i = i +'  '+  charges_list[idx] +'\n'

            idx += 1

            handle.writelines(i)

        else :

            i = i + '\n'

            handle.writelines(i)

Am sure there is a bash or perl one liner for the code above.
... How to merge two files consistently line by line

The existing answer works I'm sure, however there are existing tools that do this or something quite similar.

OpenBabel is able to add charges using the option --partialcharge, though the docs are not impressive.
Autodock vina has a utility script that adds charges and extends the pdb to pdbqt, which is admittedly not exactly what you want.

From the comments: @pippo1980 suggests this page which shows how to use the openbabel Python API to set charges
Openbabel models (note that it includes fromfile meaning you can supply your own):
obabel -L charges
eem    Assign Electronegativity Equilization Method (EEM) atomic partial charges. Bultinck B3LYP/6-31G*/MPA
eem2015ba    Assign Electronegativity Equilization Method (EEM) atomic partial charges. Cheminf B3LYP/6-311G/AIM
eem2015bm    Assign Electronegativity Equilization Method (EEM) atomic partial charges. Cheminf B3LYP/6-311G/MPA
eem2015bn    Assign Electronegativity Equilization Method (EEM) atomic partial charges. Cheminf B3LYP/6-311G/NPA
eem2015ha    Assign Electronegativity Equilization Method (EEM) atomic partial charges. Cheminf HF/6-311G/AIM
eem2015hm    Assign Electronegativity Equilization Method (EEM) atomic partial charges. Cheminf HF/6-311G/MPA
eem2015hn    Assign Electronegativity Equilization Method (EEM) atomic partial charges. Cheminf HF/6-311G/NPA
eqeq    Assign EQEq (charge equilibration) partial charges.
fromfile    Assign charges from file containing {'atom-name', charge} pairs
gasteiger    Assign Gasteiger-Marsili sigma partial charges
mmff94       Assign MMFF94 partial charges
none    Clear all partial charges
qeq    Assign QEq (charge equilibration) partial charges (Rappe and Goddard, 1991)
qtpie    Assign QTPIE (charge transfer, polarization and equilibration) partial charges (Chen and Martinez, 2007)

The issue here is that the charges were already calculated by other software (using a semi-empirical approach that is much better than force-field).

@Camps the indicated fromfile approach using openbabel should accommodate inputs of the sort that you indicate in your question. Maybe I am misunderstanding that you have separately computed some estimates that need to be added to a PDB file.

nice use of obabel here too https://leela-dodda-yale.medium.com/atomic-charges-from-openbabel-8943394f2d62


text:
how to use CIBERSORT
I have seen in several articles, the use of CIBERSORT to calculate the proportion of immune system cells using gene expression data, I would like to try it but I have not found how to install it for R, I have tried to register on the website but it doesn't allow me to enter after introduing my e-mail and password.

Could you clarify the issue ?


text:
Intepreting and applying ordinal logistic regression coefficients to calculate probabilities?
Can someone help hint me how I can interpret ordinal logistic regression coefficients and how I can use the .L, .Q and .C terms to calculate probabilities? I am analysing a dataset, where people complain of chest or not pain when they receive external supplements for a vitamin X or not. so I performed a logistic regression, using categorized and ordered (so it is an ordinal categorical variable) values of vitamin X from these observations. So I have grouped these vitamin X values within my data as very low, sub-optimal, optimal and high values and ordered that way.
The data are as such:
chest_pain_status: binary, 0 means not having chest pain, 1 means having chest pain
suppl = binary, 0 means not supplementing and 1 means supplementing
vitx_levels: very_low < suboptimal < optimal< high
My approach is as follows:
mod <- glm(chest_pain_status ~ suppl*vitx_levels, dataset = data,family = binomial())

summary(mod)
                        Estimate       std.Error     z value          Pr( <|z|)
(intercept)             -2.7816        0.2993        -9.2927          0.000
suppl1                   0.9822        0.3209         3.0803          0.0021
vitx_levels.L           -0.1751        0.0634        -2.7391          0.0062
vitx_levels.Q           -0.1537        0.0552        -2.7689          0.0053
vitx_levels.C            0.0164        0.0460         0.3579          0.7203
suppl1:vitx_levels.L     0.16821       0.0739         2.7472          0.0229
suppl1:vitx_levels.Q     0.2433        0.0641         3.7932          0.00015
suppl1:vitx_levels.C    -0.0096        0.0537        -0.1787          0.8582

Question1: Please how do I interpret coefficients of vitx_levels.L, vitx_levels.Q and vitx_levels.C? From my googling I can understand that R tries to test linear, quadratic and cubic relationships between vitx_levels and the log odds ratios of developing chest pain. It therefore means that there is a significant linear and negative relationship between vitx_levels and chest_pain. Same for the quadratic term. I hope I interpreted it right? Any additional ideas?
Question2: What values could I substitute in to vitx_levels.L, vitx_levels.Q and vitx_levels.C If I wanted to create a formula like the one below where I will later transform y to get probabilities:
y <- -2.7816 + 0.9822 * suppl1 - 0.1751 * vitx_levels.L - 0.1537 * vitx_levels.Q + 0.0164 * vitx_levels.C +  0.16821 * suppl1:vitx_levels.L + 0.2433 * suppl1:vitx_levels.Q - 0.0096 * suppl1:vitx_levels.C

So for the function above, what values could I fit into vitx_levels.L, vitx_levels.Q and vitx_levels.C to be able to compute the value of y?

Which way around are the coefficients? Is the control first or the diseased patients first? Where's the coefficient for chest pain - it's in your equation but not in the output? I assume "intercept" is the residual is that right?

@M__. Thanks for your attention. Chest pain is the outcome and would not have a coefficient. Intercept is not the residual.

Yes of course - thanks. Where is the residual (every GLM should have one)?

Its a while since I've done GLM because ML is the method of choice.
Question1 Yes multi-linear regression is an approximate description of GLM.
Generally the weights simply mean whether a feature is positively or negatively associated with the phenomena. It's important to be sure it's the right way around (negatives and positives can flip, if that makes sense).
suppl1 is significantly and relatively strongly associated with symptoms, whilst L and Q vitamin levels appear to have some protective effects.
Suppl1 has a dose dependent effect with L and Q having a lower correlation with the observed symptoms.
The dose dependency result therefore confirm suppl1 as associated with the symptoms. In all cases levels C are not significant, I assume this is a low level supplement.
Summary vitx_levels appear to have some protective effect.
I do stress it is important to ensure the "negative versus positive" is the right way round. The residual is important and I suspect you've omitted this in the GLM. The lower the residual the better the fit of the model.
BTW LD is different from GLM although they are related (title versus body of the question).
Question 2 Yes that is correct 0.9822 * 1 I assume you mean 0.9822 * suppl1 ... However all C weights are not significant, so you could omit them.

residual versus intercept ... I'm pretty certain that these are the same thing (its the c variable) and the intercept here appears quite large.

I am very grateful for your very elaborate response. I wish to ask what are the possible values vitx_levels.L or vitx_levels.Q could take in the equation below            y <- -2.7816 + 0.9822 * suppl1 - 0.1751 * vitx_levels.L - 0.1537 * vitx_levels.Q + 0.0164 * vitx_levels.C +  0.16821 * suppl1:vitx_levels.L + 0.2433 * suppl1:vitx_levels.Q - 0.0096 * suppl1:vitx_levels.C. For example, let us say I have a x=2 could I say for the  vitx_levels.Q  term in the equation above that  vitx_levels.Q  = x^2 since  vitx_levels.Q  is the quadratic term?

@Charles Yes its predictive model


text:
Confirming Identity of clone in Phylogeny
This might be something very trivial to ask, but I am new to this and was confused as in how to identify or differentiate between very closely related isolates from clonal isolate from the whole genome/core-genome phylogeny?
For instace, I am drawing a core-genome phylogeny of Staphylococcus aureus genomes taken from NCBI as well as from my own sequenced ones. Now, my isolates have grouped with several isolates from NCBI (Bottom four greens with four red ones). I want to know, weather these isolates are actually closely related to NCBI ones or they are just those strains themselves?

Thank you all in advance!

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

Please [edit] and tell us what you are doing, what kind of data you are working with etc.

Can you provide a tree that represents your data? If you do that I'll explain how to interpret it.

Thank you very much for the reply! I have tried to add details to the best of my understanding!

Your isolates are both robustly part of the NCBI lineage of 4, i.e. now making a lineage of 8 taxa and genetically very close to it, but slightly distinct.
The robustness is given by the bootstrap values (100) and the genetic proximity is given by the genetic distances for each branch. What you want to do is look at the epidemiological information for your strains and those in Genbank, there should be some sort of recognised association, e.g. country, possibly locality, anything else?
The association has epidemiological significance, thus it's not random, if that is your question. There will definitely something connecting this lineage of 8.

Yes, they indeed are from same country. Thank you very much!!


text:
Identification of metabolites through spectral library search
I am very new to single-cell metabolomics and today I want to jump in this area. I have been reading many review papers related to identification of metabolites at single-cell resolution from spectral data through spectral library search and I have a question that cannot find out a solution.
The common workflow to identify mebolites I have discovered so far is that starting from, let say, 1000 live cells, we extract the content of each cell, put them to, let say, GC-TOF-MS technology to have spectral graphs (x-axis presents m/z values and y-axis presents intensities), and finally compare the spectral graph of each cell against a certain spectral library to infer metabolites.
The thing is we know that a single cell has multiple metabolites regulating cellular functions of that cell. However, the data I see show that each cell has only one spectrum (the figure below) => we can only find out 1 metabolite for one cell??? Is that true? Or something I am missing?

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.


text:
Need help with binding DB API
import pandas as pd
import requests
import xml.etree.ElementTree as ET
df = pd.read_excel('file.xlsx')
smiles = df['SMILES'].to_list()
metabolite =  df['Plant_metabolite'].to_list()
def downloader(smile):
    url = "https://bindingdb.org/axis2/services/BDBService/getTargetByCompound?smiles={SMILES}&cutoff={similarity_cutoff}"
    if type(smile) != str:
        return None
    else:
        similarity_cutoff = "0.85"

        url = url.replace("{SMILES}", smile)
        url = url.replace("{similarity_cutoff}", similarity_cutoff)

        response = requests.get(url)

        if response.status_code == 200:
            response = response.text

        else:
            return None
        return response
for i in range(0,len(smiles)):
    resp = downloader(smiles[i])
    if resp == None:
        pass
    else:
        tree = ET.fromstring(resp)
        dictionary = {}
        for j in range(3,len(tree)):
            for x in tree[j]:
                if x.tag[29:] not in dictionary.keys():
                    dictionary[x.tag[29:]] = []
                dictionary[x.tag[29:]].append(x.text)
        df = pd.DataFrame(dictionary)
        if len(df.columns) > 0:
            df = df.loc[df['tanimoto'] > "0.85"]
            df = df.drop_duplicates(subset='smiles',keep = 'first')
            df.replace({'na': pd.NA}, inplace=True)
            df = df.dropna()
            name = "Valeriana jatamansi/{}.csv".format(metabolite[i])
            df.to_csv(name,index = False)
        else:
            pass

This my code which I am using to download targets for my compound, but there is a difference between the output returned by the API and in the online database? Like the names of the targets and other stuff...
For example: For alpha-cadanol (SMILES: CC1=CC2C(CC1)C@@(O)CC[C@@H]2C(C)C)
This is the downloaded version:

This is the online version:

Is there something wrong in the code, or is something else the problem here? Is there any other method to do the same thing?
The similarity cutoff is 0.85 both in scripts and online, and affinity value for online portal is set up as 10000, no affinity filter is provided by the API.

It would help if you showed an example where the code and the online database show different results.

Edited, please check the updated question now and I hope it's clear now. It's a partial screenshot of the online version, the are a few more targets which are left out.


text:
How to visualize FIMO outputs?
I used FIMO to scan several of my promoter sequences for TF binding motifs. The output .tsv file gives me the motif name/sequence plus the start and stop location, how can I visualize these outputs to get a graph that shows each motif as differing length rectangles along each promoter sequence, along with a gradient that shows the p-value?

Potential workaround: use MAST.
As far as I understand, 'under the hood', so to speak, FIMO and MAST work pretty much the same. The difference might be just that FIMO focuses on scoring motifs, while MAST is scoring sequences (https://groups.google.com/g/meme-suite/c/qMSj6AdFxoM)
Example of FIMO output: https://web.mit.edu/meme_v4.11.4/share/doc/examples/fimo_example_output_files/fimo.html
Example of MAST output: https://web.mit.edu/meme_v4.11.4/share/doc/examples/mast_example_output_files/mast.html
There might be small differences in the output of the tools, but then it's up to you to see whether this is relevant for your data.

FIMO and MAST have entirely different purposes. FIMO identifies statistically significant matches to motifs in a collection of sequences. MAST identifies the sequences that have the best overall match to a collection of motifs. The output of MAST does include a cartoon showing the motif matches for a particular sequence, it only does this for the sequences that have the best overall match to ALL the motifs in the input. It won't necessarily show the  best matches to individual motifs. The typical use case for MAST is identifying promoter sequences regulated by multiple transcription factors.

(cont) While internally MAST does perform a step equivalent to FIMO the final MAST output is not equivalent to the FIMO output.

If the sequences you are scanning with FIMO are supported by the UCSC  or Ensembl genome browsers you can upload the GFF file generated by FIMO as a custom track: UCSC Custom Track, ENSEMBL Custom Track. Note that for the Ensembl website the FIMO GFF file is in GFF3 format.


text:
Remove sequences from a fasta file with IDs from a text file using Python
a python beginner here.
I have a fasta file with 2500+ sequences, and after doing some analysis I want to remove around 200+ sequences based on the matching IDs. Now, I have one fasta file (as sample.fa) and a text file with a list of IDs for the sequences that I want to remove from this fasta file. Below is an example of the fasta sequences. For simplicity, I have shortened the sequences and their respective IDs. The ID in the text file (for the sequences to remove, and in the fast file will be exactly same). The sequences were originally downloaded form NCBI and had very long headers so I have them shortened as at this point I don't need all that info, usually just the sequence ID in the headers is enough.
Sample fasta sequences in fasta file:
>WP_017417394.1-Bacillus_amyloliquefaciens_gp.
MRGKKVWISLLFALALIFTMAFGSTSPAQAAGKSNGEKKYIVGFKQTMSTMSAAKKKDVISEKGGKVQKQFKYVDAASAT
>WP_268291751.1-Bacillus_inaquosorum
MRRKKLWISLLFALTLIFTMAFSNMSAQAAGKSSTEKKYIVGFKETMSAMSSAKKKDVISEKGGKVQKQFKYVNAATATL
>WP_003327717.1-Bacillus_atrophaeus
MRSRKLWIGLLFAFTLVFTMAFSSMSPAQAAGKSSEEKKYIVGFKQTMGAMSAAAKKDVIAENGGKLQKQFKYVDAATAT
>WP_206702175.1-Bacillus_vallismortis
MAFSNMSAQAAGKSSTEKKYIVGFKQTMSVMSSAKKKDVISEKGGKVQKQFKYVNAAAATLDAKAVKELKQDPSVAYVEE
>KAF1680298.1-Bacillus_sp.
MKSKKLWISLLFALTLVFTMAFSNMSAQASGKNSTEKKYIVGFKQTMSTAKKKDIISEKGGKVQKQFKYVNAASATLDQK
>WP_172293960.1-Bacillus_sp.
MGKKRLWYSLLTVLLLVFSMAFSHPANAVKSSAKDVKKDYIVGFKSSLKTTSAKKDVIKENGGKVDKHFRSINAAKVTLD
>WP_082194748.1-Bacillus_xiamenensis
MCVKKKNVMTSVLLAVPLLFSAGFGGSMANAETASKSESEKSYIVGFKASATTNSSKKQAVTQNGGTLEKQYRLINAAQV
>WP_048353666.1-Bacillus_glycinifermentans
MKKRSLWLSVLTALLLVFTMAFSNPASAEQPAKDVEKDYIVGFKSSVKTASVKKDVIKESGGKVDKQFKIINAAKATLDQ
>WP_206702858.1-Bacillus_sonorensis
MTALLLVLSTAFSSPASAAQPAKDVEKDYIVGFKSSVKTAAVKKDVIKENGGKVDKQFKIINAAKATLDQEEVKALKKD

The list of IDs in a separate text file that need to be removed from the above fasta file:
>WP_017417394.1-Bacillus_amyloliquefaciens_gp.
>WP_003327717.1-Bacillus_atrophaeus
>KAF1680298.1-Bacillus_sp.
>WP_048353666.1-Bacillus_glycinifermentans
>WP_206702858.1-Bacillus_sonorensis

I searched and tried a few scripts shared elsewhere (e.g. from here), but I couldn't make it work with my files. So all the help will be really appreciated. Thanks in advance.
PS. I am using using python usually with its own IDLE on windows 10 PC.
I did one thing with this code. I removed the '>' character from the headers file and then this code works but it keeps those headers and their respective sequences. I think it works as intended in this code, but I want to remove those from the sequence file not to keep.
Update: as asked below, I copied this code from another post.
from Bio import SeqIO
with open('headers.txt', 'r') as fin:
    to_keep = set([])
    for line in fin:
        to_keep.add(line.strip().split(' ')[0])
with open('output.fa', 'w') as fout:
    for record in SeqIO.parse('sample.fa','fasta'):
        if record.id in to_keep:
            SeqIO.write(record, fout, 'fasta')

And how did they not work? I was about to write an answer but it is essentially what  had already written as an answer to the other question. Please [edit] and clarify why those other solutions didn't work, since right now, this looks like a duplicate question.

Does this answer your question? Filtering Sequences (entries) by headers ID from a Fasta file database

This is not a good FASTA file format. IDs and description should be separate by space. I suggest to take a look at the official documentation. Good practice will make your life easier here an example or

yes, I agree. but I do that on purpose b/c after some manipulation, I make a phylogenetic tree and visualize (usually on ITOL) which messes up the formatting, sometimes it is very long and I am interested to see only the accession ID and the bacterial name. So I delete the rest and make it as one continuous string so everything I need stays in and still within a short limit to be visible.

After modifying the code above in the original post, this one works the best for me. I got it working after putting 'not' in the second last line instead of "if record.id in to keep" to "if record.id not in to keep":
from Bio import SeqIO
with open('headers.txt', 'r') as fin:
    to_remove = set([])
    for line in fin:
        to_remove.add(line.strip().split(' ')[0])
with open('output.fa', 'w') as fout:
    for record in SeqIO.parse('sample.fa','fasta'):
        if record.id in to_remove:
            pass
        else:
            SeqIO.write(record, fout, 'fasta')

You can use something like seqkit: https://bioinf.shenwei.me/seqkit/
seqkit has a grep sub-command that is one of the most flexible I've ever seen. Your use case is quite simple, you can use seqkit grep -v like so:
seqkit grep -v -n -f list_of_IDs_to_exclude.txt in.fasta > out.fasta

This should work
import argparse

ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.ArgumentDefaultsHelpFormatter)
ap.add_argument('--fasta', required=True, type=str, help='Fasta file to extract from')
ap.add_argument('--headers', required=True, type=str, help='Headers to remove, one per line.')

conf = ap.parse_args()
f = conf.fasta
h = conf.headers
header_set = set()

with open(h) as header_file:
    for line in header_file:
        line = line.rstrip()
        header_set.add(line)

with open(f) as fasta_file:
    remove_sequence = False
    for line in fasta_file:
        if '>' in line:
            header_id = line.split()[0].replace(">", "")
            if header_id in header_set:
                remove_sequence = True
            else:
                remove_sequence = False
                print(line.rstrip())
        elif not remove_sequence:
            print(line.rstrip())

exit()

Save it as python script and then run
python remove_fasta.py --fasta sample.fa --headers listofheaders.txt > cleaned.fasta

Thanks. I tried this first thing in the morning. First it ran but didn't remove anything, then I removed '>' character from the headers list and it gave the correct expected cleaned file.
What if I want to run it directly form the python IDLE and not from cmd?

I'm glad it worked. As you can see, there are other options that will get the job done, but I understood that you want to learn how to do this in Python. I don't personally use the Python IDLE on Windows, but I believe you can save your code, open it in IDLE, and click on "run module". Make sure to modify the file locations as needed. I hope this helps, and good luck with your coding

I would argue writing a whole python script for this is overkill. Assuming you are using a*nix based operating system like Linux or macOS, you can use the two scripts I have posted here, FastaToTbl and TblToFasta, and just grep out the ones you do not want:
FastaToTbl sample.fa | grep -v -xFf listofheaders.txt | TblToFasta > cleaned.fa

Hi terdon, I am only using windows 10 PC with python. I am trying the above suggestions as well. Writing my own code is a big opportunity to learn as well.

Oh. Windows? OK. But then please add that to your question. Windows isn't really used by bioinformaticians at all so nobody would guess you are trying this in windows. If you want to learn bioinformatics, I would urge you to use Linux or macOS or basically anything other than Windows.

I agree with you. I am already on to it. ")