text:
Different results of spearman correlation between TPM and FPKM
TPM and FPKM of RNA-Seq data form GDC TCGA calculated based STAR were retrieved, respectively. The correlation between a specific gene, e.g. HIF1A, and other genes were calculated based on TPM and FPKM, respectively. And the significant genes were filtered based on P value < 0.05. Here is the intersection of these two significant gene set shown in venn diagram. IMO, the two gene sets should not be obiviously different. Ang advice for this discrepance would be appreciated!

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.

This shouldn't be surprising that you see different correlations between gene expression data when expressed in different units.
To see why, let's look at how these units are defined.
Let's denote the absolute number of fragments we sequence from an RNA-seq measurement as $X_{gi}$ for gene $g$ in sample $i$, and $l_g$ as the length of gene $g$ in kilobases.
Then our different expression values will be:
$$
f_{gi} = \frac{ X_{gi} \cdot 10^6 }{ l_g \cdot \sum_h X_{hi} } \text{ in FPKM}
$$
or
$$
t_{gi} = \frac{ X_{gi} \cdot 10^6 }{ l_g \cdot \sum_h \frac{ X_{hi} }{ l_h } } \text{ in TPM}
$$
You convert from FPKM to TPM units like so:
$$
t_{gi} = 10^6 \cdot \frac{ f_{gi} }{ \sum_h f_{hi} }
$$
As you can see, within any single sample, $i$, this last denominator is fixed.
But it is not fixed across samples (i.e. $\sum_h f_{hi} \ne \sum_h f_{hj}$ if $i \ne j$).
This means that the ordering for a single gene's expression across samples is not preserved when you switch between units (i.e. if $f_{gi} > f_{gj}$, it is not guaranteed that $t_{gi} > t_{gj}$).
Because Spearman correlation is based on the order of each gene across samples in your dataset, and the order of genes in FPKM will be somewhat different from the order of genes in TPM, you will find different correlations between HIF1A and other genes when using the different units.
Here's a simplified, concrete example.
Imagine that your transcriptome is made of up 3 genes (A, B, and C) with lengths 1 kb, 2 kb, and 3 kb, respectively.
And imagine you're doing RNA-seq in 4 samples.
If your absolute count matrix is this
$$
X = 10^6 \cdot \begin{bmatrix}
  1 & 2 & 9 & 7 \\
  2 & 2 & 6 & 5 \\
  3 & 2 & 2 & 5
\end{bmatrix}
$$
Then the expression matrix in FPKM is
$$
F = \begin{bmatrix}
  166\,666.7 & 333\,333.3 & 529\,411.76 & 411\,764 \\
  166\,666.7 & 166\,666.7 & 176\,470.59 & 147\,058.82 \\
  166\,666.7 & 111\,111.1 & 39\,215.69 & 98\,039.22
\end{bmatrix}
$$
and in TPM its
$$
T = \begin{bmatrix}
  333\,333.3 & 545\,454.5 & 710\,526.32 & 626\,865.7 \\
  333\,333.3 & 272\,727.3 & 236\,842.11 & 223\,880.6 \\
  333\,333.3 & 181\,818.2 & 52\,631.58 & 149\,253.7
\end{bmatrix}
$$
If you calculate the Spearman correlation for genes B and C against A, then you get these values.

Gene
Correlation with A (FPKM)
Correlation with A (TPM)

B
0.316
-0.8

C
-1
-1

You can see that gene C is anti-correlated with A, but B actually switches signs and goes from moderately correlated in FPKM to strongly anti-correlated in TPM.
So your observation that the same comparison on the different units gives pretty different results is expected.
If you made a scatter plot with Spearman's $\rho$ from TPM along the x-axis and Spearman's $\rho$ from FPKM along the y-axis for all the genes that had a correlation p-value < 0.05, you'll probably have a lot of genes that also flip signs like I showed with this toy example.
This paper by Wagner et al. is probably a good paper to review. It can give some context about RPKM and how to think about count-based sequencing data. Correlation-based analyses on data like this are typically used as a first-pass approach to get an idea about how things are behaving. But it's not very precise, as you've experienced. If you're doing exploratory work, I'd be inclined to pick TPM units, look at the results you have, then try to come up with a differential analysis question that can be handled in a more statistically robust way.

Got that! Thanks so much!


text:
rRNA genes missing from metagenomic bins; is there a way to recruit them?
I have a couple of metagenomic bins which are okay in quality, but often missing rRNAs (16S, 23S...). I assume this is due to high population variability, combined with the high conservation of rRNAs, and reads of at most 300 bp. My binning software of choice is vamb.
Could this be fixed? Once a metagenomic bin is built, is there a way to see which 16S (or multiple 16S possible sequences) are likely to be associated with it? I considered using tools which explore the assembly graph, either manually (Bandage) or automatically (spacegraphcats), but I'm not sure if this is the right approach. After all, if the rRNA genes were linked in the assembly graph to genes that now clearly belong to a bin, they would have been added to that bin...It feels like this is a problem others may have encountered before, so I was wondering if anyone has an idea : )

I am not super familiar with vamb, but it doesn't appear to use the assembly graph: "Vamb is a metagenomic binner which feeds sequence composition information from a contig catalogue and co-abundance information from BAM files into a variational autoencoder and clusters the latent representation." Most binners I'm familiar with don't pay attention to the assembly graph. I'd expect that there is some assembly graph signal, possibly the issue is that there is a really nasty knot making it hard to navigate the graph properly.

@MaximilianPress Sorry, I was confusing. I meant using the assembly graph from my assembly program (MEGAHIT) and then trying to correlate the contigs with the bins. If you'd like to expand on the knot it would still be helpful! ('here are the reasons this would be difficult' is a valid answer : )

Ah, I mostly meant that "After all, if the rRNA genes were linked in the assembly graph to genes that now clearly belong to a bin, they would have been added to that bin..." isn't necessarily true if your binner doesn't pay attention to the graph. In conclusion: yes please do use spacegraphcats or a similar tool.

I considered mentioning this but considered it implausible. However, apparently people are actually able to use PE shotgun data to link rRNA genes to MAGs.
This is something that is quite trivial with linked-read or Hi-C data, but I am assuming that you don't have these data (and don't have the desire to go back and generate them). However you presumably have PE shotgun data that you used to generate your MAGs.

Really impressed that you linked something so recent, thank you!

source: Stack Exchange

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

sort: votes


text:
Generate a single column from two different columns in the same format
I have multiple data frames in R in the following format:
Dataframe1:
Var1    Var2    cor p   fdr
Bacteroides_galacturonicus  Lactobacillus_rogosae   0.967374392060265   0   0
Bacteroides_galacturonicus  Lachnospira_pectinoschiza   0.526779994497013   0   0
Lactobacillus_rogosae   Lachnospira_pectinoschiza   0.29419393866076    0   0
Bacteroides_galacturonicus  Bacteroides_salyersiae  0.374581964109831   0.548   0.829267117726658
Lactobacillus_rogosae   Bacteroides_salyersiae  0.45736458243079    0.002   0.03622024291498
Lachnospira_pectinoschiza   Bacteroides_salyersiae  0.566636776534  0.21    0.56910941475827
Bacteroides_galacturonicus  Holdemanella_biformis   0.342443352 0.094   0.394447279549719
Lactobacillus_rogosae   Holdemanella_biformis   0.15592768883205    0.174   0.527328455284553

Dataframe2:
Var1    Var2    cor p   fdr
Lactobacillus_rogosae   Bacteroides_galacturonicus  0.467374392060265   0   0
Lachnospira_pectinoschiza   Bacteroides_galacturonicus  0.426779994497013   0   0
Lachnospira_pectinoschiza   Lactobacillus_rogosae   0.69419393866076    0   0
Bacteroides_galacturonicus  Bacteroides_salyersiae  0.074581964109831   0.548   0.829267117726658
Lactobacillus_rogosae   Bacteroides_salyersiae  0.320636458243079   0.002   0.03622024291498
Lachnospira_pectinoschiza   Bacteroides_salyersiae  0.132857736776534   0.21    0.56910941475827
Bacteroides_galacturonicus  Holdemanella_biformis   0.19490130543352    0.094   0.394447279549719
Lactobacillus_rogosae   Holdemanella_biformis   0.14472768883205    0.174   0.527328455284553

These are two correlation matrices that I want to bind through rbind(). The problem is that the exact correlation is represented in two ways in those two dataframes.
A bacteria is found in Var1 of Dataframe:1 while its pair is found in Var2 of Dataframe:2 (see first correlation).
I want to make a column so that each pair of correlations can be represented in the same format. e.g.:
Var1:Var2   cor p   fdr
Lactobacillus_rogosae:Bacteroides_galacturonicus    0.967374392060265   0   0
Lactobacillus_rogosae:Bacteroides_galacturonicus    0.467374392060265   0   0

Can someone help me achieve this, please?

Could please kindly consider @user438383's answer for upvoting/accepting. It is a good answer.

I think this should be fairly straightforward:
# working on Dataframe2
Dataframe2$tmp <- Dataframe2$Var2   # Create a temporary column of Var2
Dataframe2$Var2 <- Dataframe2$Var1  # Replace Var2 values by Var1
Dataframe2$Var1 <- Dataframe2$tmp   # Re-attribute tmp values to Var1
Dataframe2$tmp <- NULL              # delete tmp column

# Binding the two dataframes
df <- rbind(Dataframe1, Dataframe2)

# Concatenating the two columns
df$Var1 <- paste(df$Var1, df$Var2, sep = ":")
df$Var2 <- NULL                     # delete Var2 column

PS: I haven't tested this code so let me know if it works for you.
EDIT1:
InvertedRows <- Dataframe1[!(Dataframe1$Var1 %in% Dataframe2$Var2),]
InvertedRows$tmp <- InvertedRows$Var2   # Create a temporary column of Var2
InvertedRows$Var2 <- InvertedRows$Var1  # Replace Var2 values by Var1
InvertedRows$Var1 <- InvertedRows$tmp   # Re-attribute tmp values to Var1
InvertedRows$tmp <- NULL              # delete tmp column
Dataframe2<- rbind(InvertedRows, Dataframe1[Dataframe1$Var1 %in% Dataframe2$Var2,])
df <- rbind(Dataframe1, Dataframe2)

Thanks, @user324810. The problem is not all the combination is just the opposite.Some are the same in both dataframe. In this case, all of them will just flip.

I am not sure I understand you. Do you mean that not all the rows are inverted?

Yes, that's correct. Not all rows are inverted. We can just frame the question in another way: like if we first bind the rows, then format the inverted combinations in identical way. Is it possible?

I understand better. I am on the move right now. If nobody answers in an hour or so, I will come back to resolve this issue once I am back on my PC.

Meanwhile, you can see my comment, I wrote it but didn't test.

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


text:
How does split reads look like in sam files?
I used bwa mem to align the DNA with the reference genome. If there are split reads (chimeric reads, come from two different parts of the genome), will it be split into two lines rather than one line?

Split reads as in spliced or split reads as in a read with a supplemental alignment? "Split reads" has many possible meanings.

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.

split reads here refer to chimeric reads which come from two different parts of the genome. For example, the first 1-70 bp is from a region in chromosome 1, 71-150bp is from a region in chromosome 3.

I'm assuming you mean supplementary/chimeric alignment? The SAM Format Specification has a really detailed explanation as well as its Optional Fields Specification.
In general though, you will have an additional line in the SAM file to represent a chimeric alignment yes. For a single read, one line will represent the primary alignment; this is arbitrarily chosen. The other line will represent the split/supplementary alignment and have its supplementary alignment flag set. If there are multiple chimeric alignments, each one will have its own line.
You check if a primary alignment has any supplementary alignments associated to it by checking the SA tag which will list the location(s) of any supplementary alignment(s).


text:
Can I Incorporate svaseq() into GSEA/GSVA analysis?
I understand GSEA/GSVA can take microarray-like expression matrix such as the output of voom() or vst(). However, I have a question on how we can also use svaseq() variables, correct the output from voom() or vst(), and apply GSVA on it.
I assume this is possible since I have read papers showing PCA after correction by SVA (such as this one and this one in their supplementary figure), but my statistical and coding knowledge does not allow me to recreate what they do. Is there a known pipeline that can do this
Edit: to clarify, I would like to generate a normalised matrix (e.g. rlog, or out put of voom from limma, or vst from DESeq2) that is adjusted for batch effect. ComBat, for example, can provide what I need. However, since I am doing my DESeq2 DGE analysis with svaseq(), I want to adjust it with svaseq variables. I am wondering if there are pipelines that can use svaseq vairables (or any varaibles) and adjust the value in the normalised expression matrix.

Use removeBatchEffects from limma. The input counts should be on log scale, so vst and rlog are ok. The covariate argument can take the surrogate matrix from svaseq.
https://rdrr.io/bioc/limma/man/removeBatchEffect.html


text:
What is the contributing factor in single measurements of nanopore sequencing? One base, a k-mer or difference between the base left and entered?
My question is about nanopore sequencing and specifically about the current that is measured by the device in each measurement. The question is:
In each measurement in nanopore sequencing, the change in the current between two sides of the membrane at the sensing region is recorded by the sequencer. Is it a single nucleotide or a k-mer or the difference between the first and last base in the kmer that contributes to one single measurement?
I am citing two papers here that I looked at:

Nanopore sequencing technology, bioinformatics and applications
(Nature 2021)
Three decades of nanopore sequencing (Nature 2016)

Based on my finding in these papers, I realized that all the nucleotide molecules in the k-mer in the sensing region contribute to the measurement, and not only one base, yet I am not sure if it is a good assumption.
I know some models like HMM and NN sequential models can computationally infer the individual bases (for example some basecallers). But my question is more about what is physically measured?
In the second cited paper, it is stated that, in more recent nanopore designs, DNA strand can be fixed such that one individual nucleotide is positioned in the sensing region. Still, the diameter of the sensing region is not exactly one base.
Moreover, since the current change is recorded not the current itself, the contributor could be the difference between the base which just passed from the sensing region and the one base that just entered?
So I restate my question; what is the contributing factor in single measurements of nanopore sequencing?
In case it is the k-mer that is measured, how could we find k? Is it based on the properties of the device and the kit information?

The sequencing channels measure electric current, i.e. the flow of charged particles through the nanopore. The circuits are not measuring the bases directly; they're measuring the way the flow of ions is changed when different things are introduced into the pores.
Older R9 pores had a length that could fit about five bases in the pore at any one time, the newer R10 pores have a double-size pore that fits about ten bases. Any of these bases can potentially influence the flow of current, and that degree of influence is increased at the narrower pinch points of the Nanopore. Inconsistent unwinding of double-stranded DNA (or hybrid DNA/RNA) can also affect the current profile over time, so there may also be some effects on the current that transfer from where the unwinding occurs, 20 or so bases prior to the pore introduction.
The 5/6 kmer models are (or were, with R9) a reasonable sequencing model, because that's the number of bases that would fit in the pore, so those are the bases that have the greatest influence on the current.
I'll explain using a bigger physical model of hydroelectricity:
The mental model that I have in my head of a flow cell is a dam with a lake / reservoir, and holes in the dam that slowly drain the lake, spinning turbines that produce electricity (i.e. current). The interesting bits in the lake are the things going through the holes that aren’t water, but the water flow is the only thing that can be measured: there are no cameras attached to the holes, only current meters attached to the turbines.
It just so happens that those non-water things have a strange enough shape that they interrupt the water flow in large and predictable ways, which means that observing the water flow can be used as a good proxy for observing the things themselves.
If the holes get blocked, nothing can move through the holes, so there is no current flow: sequencing cannot happen.
Eventually that water reservoir will get depleted. If the reservoir is fully drained [ionic balance of the carrier ions in the flow cell is equalised], then it doesn’t matter how clean the holes are; the dam will produce no more electricity [sequencing cannot happen].
As an aside, this model helps me to understand why keeping the pores occupied with long DNA prolongs the sequencing life of the flow cell. The flow of carrier ions is reduced when a pore is actively sequencing, so the carrier ion reservoir drains slower, and can keep sequencing for longer.


text:
How to remove redundancy from a gtf file?
I have an annotation file. I would like to remove redundancy, as shown in the example (in the real file, I have a lot of these redundant cases).
I would like to consider only one of the following genes (the longest could be a good choice). In the case below, the 2 genes have several common exons. As you can see here in the picture, there are two identified genes in the same locus (one identified through Augustus, the other one through GeneMark. They have the same start codon, and the same initial codons, but different final codons and stop codons. Hence, I'd like to retain only one of the two genes:

How can I do that? I have already tried with bedtools, agat_sp_keep_longest_isoform.pl but with no result.
Thanks for the support!

can you show your attempts with those two tools?


text:
How to align integer (non-DNA/protein) sequences?
I am looking for an algorithm to find the "best" alignment between two sequences of integers similar to how one aligns nucleic acids or amino acids for homology comparisons. For example, the best alignment for the two sequences below is:
                10 06 13 91 22 16 06
                |  |  |  |  |  |  |
       01 01 02 10 04 11 92 22 17 05

or an example with gaps:
                10 _  _  91 22 16 06 22 12
                |  |  |  |  |  |  |
          11 02 10 04 11 92 22 17 05

The alignment would be based on the distance between each number similar to aligning amino acid sequences where one would prefer to align a hydrophobic amino acid with another hydrophobic amino acid because they are more "similar". I hope this explanation makes sense.
Another way to think about this problem is if you took the 20 amino acids and replaced them with a number from 1-20 instead of their single letter abbreviation and then tried to align 2 sequenced based on how close the numbers where, how would someone align them? In this example the alignment is based on how close two integers are from one another whereas when aligning amino acids the alignment is based on the chemical similarity of the amino acids, ie leucine is very similar to isoleucine whereas valine is very dissimilar to tryptophan.
I've been experimenting with using the biopython pairwise library and trying to utilize it for my problem. There is a similar approach to what I am looking for here:
https://stackoverflow.com/a/45512865/1255817
This solution solves my problem, but only for the single character digits [0-9] and my numbers are more than one character.
Does anyone know of an approach (algorithm, package, library) to align sequences of integers? Ideally, I am looking for something in Python with a minimal working example.

The tricky part about your request is the size of the alphabet.  If you are happy with alphabets of up to about 250 "characters", then the programs in the FASTA package (github.com/wrpearson/fasta36) can be modified to do what you need.  The are designed to work with unsigned char strings, and could work with matrices of up to 250 x 250 or so.
There is a non-vectorized version of Smith-Waterman that can calculate optimal local alignments with alphabets extended to 250 characters (with the appropriate scoring matrix).
However, I suspect that pairwise could do the same thing (with modifications), again using an expanded scoring matrix.  The affine-gap Smith-Waterman algorithm is pretty simple, the trick is to get your multi-digit "characters" into a single byte (or unsigned char).

To parse data into pairwise2 so you can follow the SO post here fairly easy by simply converting your output values into a higher base.
from numpy import base_repr
base36 = []
for count in range(1,20):
    base36.append( base_repr(count, 36))
myseq36i = str(''.join(base36))

Output

123456789ABCDEFGHIJ

This represents 1 to 20 as single digits (of course there's 36 single digit integers) and then implement the pairwise2 solution. The advantage of this solution is you're not stuck on 20 amino acids. str may not be needed.
The custom matrix described in that post you would supply can be based around the exact base 10 differences under the scoring matrix - to avoid subtracting in the higher base.

Please do upvote the answer on SO as it is relevant to your problem. The first solution they presented is cool and am sure there's broader application.


text:
How to convert 10x matrix.mtx.gz files to hdf5 format?
I have a bunch of folders containing barcodes.tsv.gz, features.tsv.gz, and matrix.mtx.gz from the Cellranger Count output for a single-cell dataset that was sent to me from another lab. I need the .h5 version of these count matrices, but I do not have access to the original Cellranger Count output which contains the generated .h5 files.
Is there an R or Python function that can reliably convert these files into the proper hdf5 format? I haven't found anything online yet.
(x-post: https://www.reddit.com/r/bioinformatics/comments/veku5r/how_to_convert_10x_matrix_tsv_files_to_hdf5_format/)

Conversion into HDF5 from Python3 is straight forward. Could you please forward your CellRanger output example? The key is ... the key and I'd need to see the stuff you want to write and look at the keys in a CellRanger file (keys are everything in HDF5)

Answered on biostars (https://www.biostars.org/p/9527636/#9527638) as well as on reddit: (https://www.reddit.com/r/bioinformatics/comments/veku5r/how_to_convert_10x_matrix_tsv_files_to_hdf5_format/icqsfwn/)
The biostars answer is better as it needs just one package (DropletUtils) as opposed to the 2-3 required by the Seurat-allied packages suggested on reddit. Plus, Seurat can take forever to load a dataset.

Conversion into HDF5 from Python3 is straight forward, if this is a simple key: value dictionary format then it should be easy to mimic the CellRanger format. If it is a more complex data structure that would be tricky.  The key to HDF5 is 'the HDF5 key' which you need to find (you would call it the directory)
To find the 'key' one approach is do is download https://www.hdfgroup.org/downloads/hdfview/ and then forward a screenshot of the CellRanger HDF5 file here. The directories are what HDF5 call 'keys' and you need these to write to and call from HDF5. There is a CLI for this operation.
An example of the data you are writing to would be needed.If per chance the data is in a Python dataframe you're in business, because these easily IO to HDF5
All you do is put the data into a simple dictionary and assign the dictionary to a key which is written to HDF5. Parsing the data from e.g. an AlignIO/SeqIO object to the dictionary is more complicated. When the HDF5 is read, the key has to be supplied then the dictionary is extracted and assembled back into an object.


text:
the following errors are found during the installation of an open source software, UNCeqR on z shell
I am trying to do somatic mutation detection via the open source software, UNCeqR
I got following error was detected on mac z shell during installation.
BEGIN failed--compilation aborted at test.pl line 11.
not ok 1
make: *** [test_dynamic] Error 2
  CALLAHAN/Math-CDF-0.1.tar.gz
  /usr/bin/make test -- NOT OK
//hint// to see the cpan-testers results for installing this module, try:
  reports CALLAHAN/Math-CDF-0.1.tar.gz

I would like to ask for any advice to proceed.

Did you look at the cpan-testers results as suggested? Is this the entirety of the output? What command did you run?

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.

It is a Perl compilation error:
(Perl libraries) Math::CDF-0.1 -- Generate probabilities and quantiles from several
statistical probability functions

conda install -c bioconda perl-math-cdf

cpan INSTALL Math::CDF

This module provides cumulative probabilities [P(X <= x)] and
quantiles [given p, returns x such that P(X <= x) = p] from several
statistical probability distributions.
This package provides a perl interface to DCDFLIB
https://metacpan.org/release/CALLAHAN/Math-CDF-0.1/source/README
make and make test of Perl libraries Math::CDF-0.1 in the installation are not okay

https://drive.google.com/file/d/1Bd8TklsUPMqDRDn7zImGEXI9Nk84vkfU/view?usp=sharing


text:
How to remove the unpaired reads in sam/bam files?
I have sam and bam files for the chimeric reads, which come from two different parts of the genome (For example, the first half of the read from part of Chromosome 1 and the second half of the read from part of Chromosome 3). I have removed the low mapping quality reads (MAPQ <= 30) and reads with high mismatches. This might result in some unpaired reads which need to be removed before further analysis. May I ask whether there is any way to remove the unpaired reads in sam/bam files? Thanks!

The -f option of samtools view is for flags and can be used to filter reads in bam/sam file matching certain criteria such as properly paired reads (0x2) :
samtools view -f 0x2 -b in.bam > out.bam

Thanks. If I remove the low mapping quality reads and reads with high mismatches after the alignment, will the flags be updated? If they are not updated, there will possibly be unpaired reads even if I use flags for filtering in my point of view.

That is a good question and there is a good change those flags won't be updated, it would need some testing. Here is a post on seqanswers that might help you.


text:
BWA-mem and sambamba read group line error
this question has been asked [and answered] on Stack Overflow
This is a two-part question:

help interpreting an error;
help with coding.

I'm trying to run bwa-mem and sambamba to align raw reads to a reference genome and to sort by position. These are the commands I'm using:

bwa mem  \
  -K 100000000 -v 3 -t 6 -Y \
  -R '\@S200031047L1C001R002\S*[1-2]' \
  /path/to/reference/GCF_009858895.2_ASM985889v3_genomic.fna \
  /path/to/raw-fastq/'S\S[^_]*_L01_[0-9]+-[0-9]+'_1.fq.gz \
  /path/to/raw-fastq/'S\S[^_]*_L01_[0-9]+-[0-9]+'_2.fq.gz | \
/path/to/genomics/sambamba-0.8.2 view -S -f bam \
  /dev/stdin | \
/path/to/genomics/sambamba-0.8.2 sort  \
  /dev/stdin \
  --out host_removal/${SAMPLE}/${SAMPLE}.hybrid.sorted.bam

This is the error message I'm getting: [E::bwa_set_rg] the read group line is not started with @RG.
My sequences were generated with an MGI sequencer and the readgroups are identified like this @S200031047L1C001R0020000243/1, i.e., they don't begin with an @RG. How can I specify to sambamba that my readgroups start with @S and not @RG?

The commands written above are a published pipeline I'm modifying for my own research. However, among several changes, I'm not confident on how to define sample id as such stated in the last line of the code: --out host_removal/${SAMPLE}/${SAMPLE}.hybrid.sorted.bam (I'm referring to ${SAMPLE}). Any insights?

Thank you very much!

Read groups must start with RG, that is a hardcoded requirement of the SAM format, nothing you can do about that.
See https://samtools.github.io/hts-specs/SAMv1.pdf

My sequences were generated with an MGI sequencer and the readgroups are identified like this @S200031047L1C001R0020000243/1

Those are MGI read names, not read groups. (And the @ is the FASTQ syntax that introduces a read name in a FASTQ file, not part of the read name or qname itself.)
Read groups more or less correspond to the sample that was sequenced. If a sample was sequenced across multiple sequencing machine runs and the resulting data merged into one file, there might be several read groups corresponding to the individual runs all pertaining to the same sample — though with increasing data volumes from individual sequencing runs this is perhaps less common than it used to be.
You get to make up the ids yourself, usually corresponding to what the sample is known as in your workflows. So your -R option should be something like
-R '@RG\tID:S123\tSM:S123\tDS:Sample 123\tPL:DNBSEQ'

(BWA's -R syntax is unfortunately not very easy to use. See addreplacerg for a similar -r option that is rather more user-friendly.)
The name you use for ${SAMPLE} would usually be the same as the value of the SM: field within that -R option — the identifier that is used in your LIMS and workflows to identify that particular sample.


text:
deseq2 makeContrasts function
This article talked about various deseq2 designs etc. One of the designs I would like to use is explained as this:
Control versus treatment average
makeContrasts((treatmentI+treatmentII+treatmentIII)/3-treatmentCTL,
  levels=colnames(design))

I'm not sure how to incorporate this in deseq2 run which step?
Is it after running this
dds <- DESeq(dds)

Can I call the makeContrasts function?

The makeContrasts function is from edgeR so it does not help here. If you want to compare one group vs the average of several others here is a suggestion of how this could be done. It uses a design without intercept so all groups are easily accessable as coefficients in resultsNames(dds):
suppressMessages(library(DESeq2))

#/ Example with four groups
set.seed(1)
dds <- makeExampleDESeqDataSet(m=12)
dds$condition <- factor(as.character(rep(LETTERS[1:4], each=3)))

#/ Use a design without intercept so all four groups are accessable as coefficients
design(dds) <- ~0+condition
dds <- DESeq(dds)
#> estimating size factors
#> estimating dispersions
#> gene-wise dispersion estimates
#> mean-dispersion relationship
#> final dispersion estimates
#> fitting model and testing

#/ results using listValues to get A vs the average ov B,C,D
res <- results(dds, contrast=list(c("conditionA"), c("conditionB", "conditionC", "conditionD")), listValues=c(1, -1/4))

Created on 2022-06-21 by the reprex package (v2.0.1)

Determining whether or not it is appropriate to do this in a particular situation is up to you. You have to know what question you want to answer for your work. I cannot help with that.


text:
deseq2 model design : Different gene output
Here I have two designs:
Design 1
 dds <- DESeqDataSetFromMatrix(countData=rpkm_ordered,
           colData=coldata, design= ~ Sex + Age + TMB + WBC + BM_percentage + FAB)

To see the results I did this
resultsNames(dds)
 [1] "Intercept"          "Sex_Male_vs_Female" "Age"                "TMB"                "WBC"                "BM_percentage"      "FAB_M1_vs_M0"
 [8] "FAB_M2_vs_M0"       "FAB_M3_vs_M0"       "FAB_M4_vs_M0"       "FAB_M5_vs_M0"

I extracted only the Age from the first design
Design 2
I again used the design using Age variable :
dds <- DESeqDataSetFromMatrix(countData=rpkm_ordered,
          colData=coldata, design= ~ Age)

As an example I take this gene this output from the first design
  gene            baseMean log2FoldChange  lfcSE  stat   pvalue       padj
  <chr>              <dbl>          <dbl>  <dbl> <dbl>    <dbl>      <dbl>
1 ENSG00000196188     97.4         -0.742  0.113 -6.58 4.57e-11 0.00000125

This from the second design:
  gene            baseMean log2FoldChange  lfcSE  stat   pvalue       padj
  <chr>              <dbl>          <dbl>  <dbl> <dbl>    <dbl>      <dbl>
  ENSG00000196188     97.4         -0.757  0.119 -6.34 2.28e-10   3.43e-07

Another example this gene from the first design; the result is not significant:
           gene  baseMean log2FoldChange  lfcSE    stat pvalue   padj
ENSG00000173826  74.54467     0.01283834  0.589  0.0218  0.983  0.993

The same gene from Design 2 is significant:
           gene  baseMean log2FoldChange  lfcSE    stat   pvalue     padj
ENSG00000173826 74.54467       -2.004743  0.190   -10.6 4.76e-26 1.35e-21

Is this drastic outcome is due to the maths behind it or am I doing something wrong?
So would my interpretation be as such when i consider Age as my only variable I can consider ENSG00000173826significant ? But it is not significant when my main effect which i want to test is FAB which is my design 1.
I would like to know how to interpret and report the same.
Any suggestion or help would be really really appreciated

It's not surprising to me that the adjusted fold change statistic changes in different ways for different genes when covariates are added or removed. That's sort of the point of using covariates, as long as it makes sense to include them (which is an experimental decision, not a mathematical or statistical decision).
One thing I'd recommend checking is to make sure the results you're looking at are what you think you are looking at. You say that you're interested in testing FAB, but that is not included in the second (Age-only) design, which means that it won't be possible to extract any FAB comparisons from the results. It looks like your first design might be looking at Sex (based on the resultsNames ordering) or one of the FAB comparisons (based on reverse ordering of variables in the model), whereas the second might be looking at Age. This should be explicitly stated when displaying the results, but you haven't included that information in your output.
If you want FAB comparisons in the results, that variable needs to be in the specified model, e.g. "design = ~ Age + FAB". The Design 2 model that you have in your question only includes Age.
A good way to do this without getting the wrong contrast is to fully specify the variables when fetching results, as explained in the DESeq2 vignette:

In any case, the contrast argument of the function results takes a character vector of length three: the name of the variable, the name of the factor level for the numerator of the log2 ratio, and the name of the factor level for the denominator. The contrast argument can also take other forms, as described in the help page for results and below

res.age <- results(dds, name="Age")
res.FAB.M3.M0 <- results(dds, contrast=c("FAB","M3","M0"))

design= ~ Sex + Age + TMB + WBC + BM_percentage + FAB)  I used this one this resource https://github.com/hbctraining/DGE_workshop_salmon_online/blob/master/homework/DGE_assignment_2_answer_key.R  "which means that it won't be possible to extract any FAB comparisons from the results" i could extract the result for FAB

this biostar post https://www.biostars.org/p/278684/ explains i hope im getting it right

We can account for the different types of sequencing, and get a clearer picture of the differences attributable to the treatment. As condition is the variable of interest, we put it at the end of the formula. I followed this ..

What gringer means is you can't get FAB results from design2, because it doesn't include FAB at all.  If you specify the contrast when getting results, the order of the elements in the design does not matter at all.   But since you won't say how you called results, no one really knows what you are doing.

Also, unless you have a huge number of samples, I am skeptical that you can really model so many different variables at once.  I would keep only those elements which correlate to the first few PCs.

@swbarnes2 yes I tried the second suggestion so if i want to achieve that 80% variance threshold it would be like 40 components may I need to filter or group the gene in better way...


text:
Construct phylogeny from a Fasta file
I have a set 189 taste receptor protein sequences (not aligned) in a fasta file. I like to get the phylogenetic tree in newick format. I was using earlier https://ngphylogeny.fr/. Unfortunately now it is not working due to some bugs. Necessarily, I just need newick information for the input fasta file.
Any help would be highly appreciated.
Here is one newick example:
((raccoon:19.19959,bear:6.80041):0.84600,
  ((sea_lion:11.99700, seal:12.00300):7.52973,
     ((monkey:100.85930,cat:47.14069):20.59201,
       weasel:18.87953):2.09460):3.87382,
  dog:25.46154);

Please [edit] your question and add it there. Are these sequences of homologous genes? Entire genomes? Are they DNA? Protein? Are they already aligned to each other?

You will first need to align your sequences, and then feed the alignment to a program that can generate trees. I am very much out of date in this field and I haven't done anything like that in at least ten years, but back when I used to do it, my workflow was:

Align the sequences using MAFFT or t-coffee.
Open the alignment in the Jalview alignment editor/viewer which has options to generate the tree:

then just save the tree in Newick format:

The screenshots were taken using the default alignment provided by JalView itself, and this is the tree I got when I saved as Newick:
$ cat tree
((FER1_MESCR:24.484375,(((O80429_MAIZE:23.699999999999996,FER1_MAIZE:17.300000000000004):21.240625,(FER2_ARATH:13.067708333333336,((Q93Z60_ARATH:8.0,FER1_ARATH:8.0):8.401785714285715,(FER_BRANA:14.796875,FER3_RAPSA:8.203125):1.5982142857142847):0.9322916666666643):18.259375):9.75390625,((Q7XA98_TRIPR:24.8515625,FER1_PEA:27.1484375):6.643229166666671,((Q93XJ9_SOLTU:10.125,FER1_SOLLC:8.875):5.61363636363636,(FER_CAPAN:5.0,FER_CAPAA:0.0):19.38636363636364):20.48177083333333):8.97265625):9.52734375):20.2578125,FER1_SPIOL:20.2578125);

Here's the result from FigTree. I thought there was a bug in the newick file, but it was fine.
I strongly recommend the answer by @terdon because there are other issues here, for example 111.16 value for the root and the 100.86 which occurs twice isn't cool.
I should have placed monkey as the outgroup, anyway just please ensure you do this for the report.


text:
postgreSQL error: "ERROR: invalid input syntax for type numeric: NULL"
I'm new to SQL. I have a table with a column of numeric  type that include NULL values. I built the table structure in postgresql, when I try to import my table I get the following error:
ERROR: invalid input syntax for type numeric: NULL
I tried both
UPDATE MY_TABLE  SET MY_COLUMN = CAST(NULL AS NUMERIC)
and setting manually the default value of the column as NULL, but none of them are working.
Is there another way to handle this problem?

Hi and welcome to the site! We only deal with questions about bionformatics here and this looks like a very standard SQL question with no connection to bioinformatics so I am afraid it is off topic. That said, I suspect you just want UPDATE MY_TABLE  SET MY_COLUMN = NULL.

I’m voting to close this question because it is about SQL and not related to bioinformatics.

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:
What is the function of the heteroatoms
I am trying to make a prediction application using the sequence of proteins.
If we think about the 6COU protein, I'm running the code below
pdb = PDBParser().get_structure("6c0u", "6c0u.pdb")

for residue in pdb.get_residues():
    print(residue)

Some of the output looks like below when looking at residues. What is the function of the heteroatoms here, and follow up question;
When reading a paper on Binding Site Detection of Proteins, there is exactly one sentence in the paper.
"We first clean the input structure by removing all heteroatoms and solvent molecules from the protein structure using the Biopython47 library."
What is the meaning of this sentence if you are working on protein sequences, heteroatoms are not already there, what does it mean to delete them?
If you want to check the paper, here

The text you cited mentions "input structure", so they work with structures, not just "protein sequences".
Heteroatoms = atoms marked as HETATM in the PDB and mmCIF formats.
In these formats each atom is listed as either ATOM or HETATM. These are not particularly meaningful categories, but they are there, so they are easy to use. It's a common practice to remove heteroatoms with grep -v ^HETATM.
Note that in the files from the Protein Data Bank not only solvent and ligands are marked as HETATM, but also "non-standard" amino-acids in a protein. An amino acid is considered standard if it can be genetically encoded. So in the files from the PDB selenocysteine is ATOM and selenomethionine is HETATM. Some programs don't follow these rules and output all protein residues as ATOM.
Removing heteroatoms works well for the vast majority of structures that have no modified residues. But in some cases it will remove also some protein residues, which may not intended.
Update: if you are wondering why there are ATOM and HETATM records, I have never came across the rationale. The PDB format
specification has had both types of records since 1976.
Different versions of the format spec use different wording,
but the meaning is always similar:

ATOM records give atomic coordinates for "standard" groups.
HETATM records give atomic coordinates for "non-standard" groups.

(Although standard residues that not in a polymer are also HETATM.)
In the mmCIF format the same record is preserved in _atom_site.group_PDB with the description:

This data item is provided for compatibility with the original Protein Data Bank format, and only for that purpose.

Most of the structures have only standard residues in protein/DNA/RNA polymers, so removing HETATM records removes only ligands and solvent.

What I understood from the paper's goal is that deleting residues where heteroatoms are not empty. The value het appears blank (het=" ") next to 95% residue of protein. When you get rid of these, the protein sequence will shorten. Could they delete the residues here because they sees it as unimportant for predicting pocket (binding) area of protein? I'm also trying to estimate pocket area using protein sequence should I delete it too?

@drorhun In the image in your question the last residue in the protein chain is PHE (resseq=350). The next Residue objects are ligands and waters.


text:
Filtering paired-end reads with sambamba: avoid discarding reads on the minus strand
I have a BAM file (DNA, shallow whole genome sequencing at ~1X) where I want to filter reads (using sambamba) to keep only those which have a template length > 20 and mapping quality > 20, discarding unmapped reads.
The filter I have used is this:
sambamba view -h \
   -t 5 \
   -F "paired and template_length > 0 and mapping_quality > 20 and not (unmapped or mate_is_unmapped)" \
   -f bam 142-ready.bam \
   -o 142-filter.bam

However, inspection of the reads in IGV led me to think that perhaps reads on the minus strand are being somehow filtered out.
Pre filter:

Post filter:

The only way I don't get such an aggressive filtering is by setting only not (unmapped or mate_is_unmapped). But that, obviously, defeats the purpose as it doesn't include what I need.
I can't, however, see what exactly is wrong with the filter. What is the correct way to filter for template length while not accidentally discarding reads on the minus strand?

The reason is that sambamba does not use the absolute value of the template length, so, if the read maps so the minus strand, the template length is actually negative. Setting a template length > 0 means that all reads mapping to the minus strand are removed.
In fact, you can't still filter for negative template lengths. See this GitHub PR.


text:
Extracting base sequences from ABI/AB1 sanger sequencing chromatogram
I am trying to understand the sanger sequencing ABI/AB1 file format better, and extract base calls from given signal intensities over time.
As I understand, reading in a raw AB1/ABI file into python, I can access the channel corrected values for the different bases in ['DATA9'] to ['DATA12'], and understand which bases belong to a specific channel through ['FWO_1'] i.e either ordered GATC or ATGC etc.
However, the file I am looking at has identified 412 bases as the output from SnapGene, yet has 4950 data records for each channel in the raw ['DATA9'] to ['DATA12']. How do you convert those 4950 data records into the bases correctly? Is there some form of normalisation - e.g every 10 raw records gives data for 1 base position, or peak calling signal detection algorithms. And if so, do you take the average over 10 records, or the highest peak in the 10 records from the channels and the channel with the highest peak/average is the correct base? Do you start at the beginning or do you have some form of cutoff? Does it make sense to convert this figure into a quality score similar to NGS, and do you use the highest peaks in the 10 record interval or the average?
I've tried using scipy's signal.find_peaks() function (python) which can get me various peaks with a minimum height and a minimum distance to other peaks, currently using a minimum peak height of 150 and distance to other peaks of 5, but still getting around 20-30 more peaks than the bases identified.
Could anyone suggest any logic, methods, algorithms or resources that could help with translating the raw signal intensities for the 4 channels into a sequence of accurate base calls?
I know software exists to identify the bases, but that's not what I am looking for. I intend to do this from scratch to understand the workings better so I can write some custom processing algorithms.
I hope this question makes sense.
Thanks

Please consider the answer below because it appears to answer the question and if you are happy to upvote and/or accept it

You can access the peak positions of Sanger Sequencing chromatograms from ["PLOC1"].
The position values in ["PLOC1"] correspond to the index values in ["DATA9"], ["DATA10"], ["DATA11"], ["DATA12"], so you can extract the signal values of each base at the positions.


text:
VFDB search output
I want to check virulence factor of bacteria protein sequence by using Virulence Factor Database and it shows No hit in VFDB BLAST search output. Will someone please explain what is meant by it?

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.

Is your question about why blast doesn't get a hit, or about why your protein isn't in VFDB? BLAST and VFDB both have documentation about how they work, can you try to investigate there?

if VFDB result show No Hit for a protein, so is it means it is non-virulent?


text:
deseq2 single factor design output interpretation
This is regarding the single factor design For example if I have Age or other continuous numerical variable how to provide that into the design formula.
For this post do i need to 'You could dichotomise your continuous variables into meaningful groups' or it can go without it?
My simple design

dds <- DESeqDataSetFromMatrix(countData=rpkm_ordered, colData=coldata, design= ~ Age)

dds <- estimateSizeFactors(dds)

dds <- estimateDispersions(dds)

dds  <- DESeq(dds, parallel = TRUE)

Here in case of metadata/coldata I m giving a single numerical value which is Blast percentage. In case of sample I have like 5 sub-types from M0 to  M5.
Now for the interpretation part How do I interpret the result?
Would it be as such
the expression differences between the my subtypes due to 'Age'
I'm bit confused since in my coldata I'm not providing any information regarding my subtypes. So when I run this
resultsNames(dds)  I get this

[1] "Intercept"     "Age"

So if I would like to know if there is a difference between which subtype due to this Age variable how do i know or my question is violating the statistical design and rationale here?
To know the differences in subtype I have do which is proving the FAB which are basically my sub-types information where I have tested pairwise.
I would like to know if I give any numerical variable to my design how do I interpret output?
Any suggestion or help would be really appreciated.

Your design is ~ Age so you can only test how Age affects each gene and it is assumed that the effect is linear. If you would like to test the effect of sub type on gene expression you'll have to introduce it to the design e.g. ~ Age + Sub.Type (or whatever the column name is). This will not test how the effect of Age is different between between the subtypes, if you would like to test this you will have to include the interaction term in the design formula: ~ Age + Sub.Type + Age:Sub.Type, this design have a lot of parameters and you will need a lot of samples to have a reliable test.
Hope it helped

thank you after reading other similar queries i had a bit of clarity, one more question if I m testing only the age what goes into the intercept in case of factor or category I normally set a reference level  in case of numerical variable what is turned into intercept is it the lowest age is turned into intercept?

Age is treated as linear so Age=0 goes to the intercept and the Age coefficient times the age is the Age effect for each individual.

Age=0 this although i don't have it in my data but what I understand is the lowest age is turned into intercept?

The intercept is the "background" effect without the effect of the covariates in the formula, it's the beta0 in the formula y = beta0 + Age*beta1

"background" I would request what would you mean by that in biological term can you explain ?I understand the beta0 when age is 0 but what it means biologically that would be really helpful

That the gene expression level at birth is the intercept value. Is it biologically true? I doubt it but that's the model.

thank you for making it easier to understand .


text:
How to normalise the histogram height in matplotlib
I am trying to create a histogram plot with matplotlib. I am using bmi column from four different data frames (d1, d2, d2, d4).
The issue is that these 4 datasets have different sample sizes. So there is a disparity between the heights of the bar.
Is there any additional option in the matplotlib to normalize the heights of the bar?
import matplotlib.pyplot as plt

def bmiplot(graphWidth, graphHeight):
    f = plt.figure(figsize=(graphWidth/110.0, graphHeight/150.0), dpi=100)
    axes = f.add_subplot(111)

    d1 = dm_group1
    d2 = dm_group2
    d3 = dm_group3
    d4 = dm_group4

    plt.hist([d1.bmi.values.flatten(), d2.bmi.values.flatten(), d3.bmi.values.flatten(), d4.bmi.values.flatten()],
    label=['bmi_dm1', 'bmi_dm2', 'bmi_dm3', 'bmi_dm4'])
    plt.legend()
    plt.rc('xtick')
    plt.rc('ytick')

graphWidth = 700
graphHeight = 800
bmiplot(graphWidth, graphHeight)

```

One approach is density, traditionally the code would be ...
kwargs = dict(histtype='stepfilled', alpha=0.4, normed=True, bins=50)

plt.hist(d1, **kwargs)
plt.hist(d2, **kwargs)
plt.hist(d3, **kwargs);

alpha is transparency allowing overlapping distributions.
Note normed is deprecated, the blurb here describes density now replacing it.
The blurb is here:

matplotlib.pyplot.hist(x, bins=None, range=None, density=False,
weights=None, ... **kwargs)

The definition of density is here

densitybool, default: False If True, draw and return a probability
density: each bin will display the bin's raw count divided by the
total number of counts and the bin width (density = counts /
(sum(counts) * np.diff(bins))), so that the area under the histogram
integrates to 1 (np.sum(density * np.diff(bins)) == 1). If stacked is
also True, the sum of the histograms is normalized to 1.

A different approach is relative frequency, via the weights command.
 plt.hist(mydata, weights= *)

*, relative counts per data type / total counts over all data types


text:
Converting .bw files to .fq (fastq)
I am able to convert .bw files to .fq(fastq) manually in Shell, but I would like to automate the process coz I have hundreds of .bw file that I need to convert. Till now, I could think of following script, but unable to solve the problem. The .bw files are in H3K27ac folder.
import subprocess

path_bw = "~/Documents/H3K27ac"
import os

files = os.path.expanduser(path_bw)
for file in files:
    if file.endswith('.gz'):
        files.remove(file)
files = [x[:-3] if x.endswith('.bw') else x for x in files]

# for .bw files

for file in files:
    bigwigtowig = "bigWigToWig " + path_bw + file + ".bw " + path_bw + file +".wig"
    subprocess.call(bigwigtowig, shell=True)

    wigtobed = "wig2bed < " + path_bw + file + ".wig > " + path_bw + file + ".bed"
    subprocess.call(wigtobed, shell=True)

    getfasta = "bedtools getfasta -fi " + path_bw + "hg19.fa -bed " + path_bw + file + ".bed -fo " + path_bw + file + ".fa"
    subprocess.call(getfasta, shell=True)

    perl = "perl " + path_bw + "fasta_to_fastq.pl" + path_bw + file + ".fa > " + path_bw + file + ".fq"
    #os.system(perl)
    subprocess.run(perl, shell=True)

Manually, I used the following to convert the .bw files:

# bigWigToWig file1.bw fileout.wig
# wig2bed <fileout.wig> fileout.bed
# bedtools getfasta -fi hg19.fa -bed fileout.bed -fo fileout.fa
# perl fasta_to_fastq.pl fileout.fa > file1.fq

All the above packages were preinstalled via shell and works perfectly. hg19 is the genom that i am using. fasta_to_fastq.pl is openly available in github and helps to convert .fa to .fq.
MAin thing i require is to convert all the files from .bw to .fq and this is the only way in my research that i could find. I need all the names to be the same as output.  For example file1.bw shoud be file1.fq , file2.bw --> file2.fq and so on.
Do you have any ideas if this is possible or are there easier ways to convert it? It would be awesome to know your view point.

Can you please explain why you want to do this? It's a very unusual situation, and it's likely that there are better ways to approach what you want to do. It is often better to ask questions about your actual goal, rather than intermediate steps (see here for more information).

What you are doing makes no sense at all. A bigwig file contains intervals and then indicates the coverage of that interval. That coverage again comes from the alignments covering that region. There is no direct connection between the interval and the fastq file other than that the fastq file is the basis for the coverage calculation. Converting the interval to a fasta and then forcing it into a fastq file is technically possible, but has absolutely no practical meaning. You cannot recover sequencing reads from a bigwig, that information is lost once you convert your alignments (usually a BAM file) into a bigwig. You should open a new question describing the problem that this is actually based on.

I'll add slightly to this excellent answer and mention that if this is public data (highly likely), there is almost certainly a link somewhere to download the underlying fastq files.

Excellent answer that misuses the word "convert" like most people on these forums do. One does not convert BAM to (big)wig, one extracts/computes wig information from a BAM file. "Convert" kind of implies information content preservation (convert sam to bam, etc.), which as you say in your answer is not true when getting WIG from BAM.

The code is verbose and I ain't at all sure some of it will work, at least in two places.
The first point is os.path.expanduser ... admittedly I prefer pathlib, however that is definitely not how I understand that command. My understanding of expanduser is if ~ is being used. You definitely want to look at globb-ing for the opening part of the code.
from pathlib import Path
import os

mypath = Path(os.path.expanduser('~/Documents/H3K27ac'))
files = mypath.glob('*.bw') # or use ('*' + '.bw')

I'm pretty certain that is what you want and its a lot easier than lines 1 to 10 in your code. The files list will then go straight into the loop.
Also, please check
files = [x[:-3] if x.endswith('.bw') else x for x in files]

Comprehensions are never logical in Python but convention is always to place the if statement after the list of interest.
Thus the following is much closer:
files = [x for x in files if not x.endswith('.bw') else x[:-3]]

Note not tested, so this may throw a bug, but its better. However, using glob bing you don't need this statement.
The other issue is the sub.process call ... Normally they are written:
os.chdir(path_bw)
file1 = 'file1.bw'
file2 = 'fileout.wig'
subprocess.call(["bigWigToWig",  file1, file2], shell=True)

The way you've written is:

Really hard to read
I don't think it would work because of the absence of the [], certainly if you are using popen it will not work

Also check the following variant:
subprocess.call(["bigWigToWig",  file1, file2], text=True)

I know its seems a bit weird, but it very likely work (don't ask why)
I agree subprocess.call is a good idea in a Python pipeline. Issuing a hard path within subprocess call I'm not a fan hence using os to get Python to move into the cwd is a better approach.
A further point is the way you've used 'file' throughout your pipeline line, hmmmm... okay, but its difficult to read. Its not generally how a pipeline is written, but just saying.
nextflow is an increasingly popular pipeline language and you may consider learning that rather than going down a Python route (Python is very good for data science), of course this particularly question could be very easily solved in bash (but for some reason it ain't trendy).


text:
How can I generate ASV file from nanopore sequencing data?
I am converting the fast5 file to fastq by using guppy basecaller after that by using kraken2 classified the sequence.
Now I am trying to generate ASV file.
Is this possible to generate ASV file from kraken output?
Or is this any other way to generate ASV file from nanopore sequencing?
I am new to nanopore sequencing data. I need help or suggestions to generate ASV file.
Thank you

What is an ASV file? You should add information about that file format to your question. The first information I get talks about MATLAB, which doesn't seem correct.

ASV: amplicon sequence variant

Dear @twinkle please do check the answer below and do consider upvoting/accepting. Its a good answer and it does help with site stats.

I am not sure that kraken/kraken2 does anything similar to this.
Usually to get OTUs/ESVs I think that you want to use a tool that clusters your sequences into groupings that represent the same sequence, allowing for technical error.
Examples might include dada2 or dadasnake, or qiime. This is after the technical steps for adapter/trimming/denoising steps.
kraken2 is more for classifying into known taxonomic categories based on a reference sequence, which is a related but distinct task.
I don't think that I understand what the ESV file format is supposed to look like. This resource seems to have some info on processing data into ESVs, but it is not for nanopore. Possibly it is still useful for the final step of creating that file format. A linked repository has a ton of one-off scripts that may be helpful but are not very well documented. For example, here is a script that they use to generate TSV files which may be the "ASV"/"ESV" format:
#!/bin/bash
source activate qiime2-2018.8

for item in `ls */*table.qza`
    do
    name=`echo $item | cut -d\/ -f1`
    qiime tools export --input-path $item --output-path EXPORTS-GTDB-tax/$name
    mv EXPORTS-GTDB-tax/$name/feature-table.biom EXPORTS-GTDB-tax/$name.biom
    rmdir EXPORTS-GTDB-tax/$name
done

qiime tools export --input-path 11_ANT-28-5-SPA_seqs_classified_against_GTDB_r86.1_807F-1050R_primer_sliced/classification.qza --output-path EXPORTS-GTDB-tax/

sed -i '1c#OTUID    taxonomy    confidence' EXPORTS-GTDB-tax/taxonomy.tsv

#Tutorial here: https://forum.qiime2.org/t/exporting-and-modifying-biom-tables-e-g-adding-taxonomy-annotations/3630
for item in `ls EXPORTS-GTDB-tax/*biom`
    do
    filestem=`basename $item .biom`
    biom add-metadata -i $item -o EXPORTS-GTDB-tax/$filestem.with-tax.biom --observation-metadata-fp EXPORTS-GTDB-tax/taxonomy.tsv --sc-separated taxonomy
done

for item in `ls EXPORTS-GTDB-tax/*.with-tax.biom`
    do
    filestem=`basename $item .biom`
    biom convert -i $item -o EXPORTS-GTDB-tax/$filestem.tsv --to-tsv --header-key taxonomy
done


text:
How can I calculate the radius of a residue?
I want to read a PDB file and display/render the atoms using a graphics engine.
If I consider each residue in the protein as a solid sphere, what would be the radius of the sphere?
How can I calculate it? Or, is there any database for that?

You may want to check how existing molecular viewers (PyMOL, Jmol, Chimera, web-based viewers on RCSB, PDBe and PDBj websites) render it.

@marcin, I am developing my own software for learning ourposes.

Van der Waals radii are used to construct a special kind of molecular model called a space-filling model. These models are constructed by drawing each atom as a van der Waals sphere with the atom’s nucleus at the center of the sphere. taken from https://www.reed.edu/chemistry/ROCO/Geometry/space_filling.html not sure how to expand it to a full residue

@pippo1980, thanks! good information. :)

try this paper not sure why 5A radius is choosen https://www.researchgate.net/figure/Each-protein-is-modeled-as-a-union-of-spheres-of-radius-5-A-one-sphere-per-amino-acid_fig2_8614866 the picture seems OK to me tough

ok full article is for free; it states: "....However, it has been suggested that the 5 Å radius is appropriate to model the hydration shell around the protein
[19]. The problem of computing
r for such a model reduces
to computation of the distance to the surface of each sphere,
and taking the minimum over all spheres.... " Worth to notice article 19 its about globular protein, no idead if 5A rules could be extended to other kind of protein too

The protein is
represented as a union of spheres, each of radius 5 Å, and
centered on the C
a atom of each amino acid residue


text:
What do these characters mean in a PDB file?
The 33rd atom information in 3nir.pdb is written as follows:
ATOM     33  N  ATHR A   2       4.935 -11.632  15.046  0.74  2.95           N
How should I interpret N  A in the above line?
I mean, what does it mean? Is it Sodium or Nitrogen?
If it is Nitrogen, what does A mean then?

have a look here : https://www.wwpdb.org/documentation/file-format

from https://www.wwpdb.org/documentation/file-format

Protein Data Bank Contents Guide:
Atomic Coordinate Entry Format Description
Version 3.30
Document Published by the wwPDB
This format complies with the PDB Exchange Dictionary (PDBx)
http://mmcif.pdb.org/dictionaries/mmcif_pdbx.dic/Index/index.html.
I believe PDB format is kind of outdated, but not sure..
At page 177 it says about altLoc:

AltLoc is the place holder to indicate alternate conformation. The alternate conformation can be
in the entire polymer chain, or several residues or partial residue (several atoms within one
residue). If an atom is provided in more than one position, then a non-blank alternate location
indicator must be used for each of the atomic positions. Within a residue, all atoms that are
associated with each other in a given conformation are assigned the same alternate position
indicator. There are two ways of representing alternate conformation- either at atom level or at
residue level (see examples).
For atoms that are in alternate sites indicated by the alternate site indicator, sorting of atoms in
the ATOM/HETATM list uses the following general rules:
• In the simple case that involves a few atoms or a few residues with alternate sites, the
coordinates occur one after the other in the entry.
• In the case of a large heterogen groups which are disordered, the atoms for each
conformer are listed together.

your question could be more useful if you post the entire description for residue 2:
ATOM     33  N  ATHR A   2       4.935 -11.632  15.046  0.74  2.95           N
ATOM     34  N  BTHR A   2       5.078 -11.406  15.180  0.31  2.78           N
ATOM     35  CA ATHR A   2       5.757 -11.521  13.850  0.81  3.02           C
ATOM     36  CA BTHR A   2       5.773 -11.153  13.921  0.20  2.67           C
ATOM     37  C  ATHR A   2       7.070 -10.839  14.210  0.74  2.82           C
ATOM     38  C  BTHR A   2       7.155 -10.559  14.193  0.29  1.80           C
ATOM     39  O  ATHR A   2       7.152  -9.941  15.050  0.80  3.31           O
ATOM     40  O  BTHR A   2       7.214  -9.641  15.012  0.25  2.41           O
ATOM     41  CB ATHR A   2       4.976 -10.693  12.813  0.87  5.53           C
ATOM     42  CB BTHR A   2       4.896 -10.354  12.941  0.25 12.07           C
ATOM     43  OG1ATHR A   2       4.611  -9.432  13.388  1.00  6.88           O
ATOM     44  OG1BTHR A   2       3.743 -11.083  12.501  0.25  9.57           O
ATOM     45  CG2ATHR A   2       3.858 -11.584  12.293  0.75 10.03           C
ATOM     46  CG2BTHR A   2       5.683  -9.885  11.726  0.27  5.90           C
ATOM     47  H  ATHR A   2       4.547 -10.814  15.527  0.75  3.44           H
ATOM     48  H  BTHR A   2       5.510 -10.211  15.754  0.25  2.90           H
ATOM     49  HA ATHR A   2       5.962 -12.339  13.548  0.75  3.32           H
ATOM     50  HA BTHR A   2       4.036  -9.929  13.477  0.25  2.86           H
ATOM     51  HB ATHR A   2       5.648 -10.589  11.938  0.75  5.43           H
ATOM     52  HB BTHR A   2       4.644  -9.326  13.574  0.25  5.67           H
ATOM     53  HG1ATHR A   2       5.030  -9.344  14.216  0.75  8.74           H
ATOM     54  HG1BTHR A   2       3.236 -11.198  13.399  0.25 10.21           H
ATOM     55 HG21ATHR A   2       4.096 -12.441  11.924  0.75 10.92           H
ATOM     56 HG21BTHR A   2       6.542  -9.278  12.024  0.25  9.66           H
ATOM     57 HG22ATHR A   2       3.222 -10.974  11.650  0.75 10.92           H
ATOM     58 HG22BTHR A   2       5.039  -9.142  11.179  0.25  9.66           H
ATOM     59 HG23ATHR A   2       3.163 -11.738  13.200  0.75 10.92           H
ATOM     60 HG23BTHR A   2       5.904 -10.639  11.169  0.25  9.66           H


text:
How can I parse alternative atom information?
I am trying to parse PDB files.
Say, a PDB file has the following data:
ATOM     33  N  ATHR A   2       4.935 -11.632  15.046  0.74  2.95           N
ATOM     34  N  BTHR A   2       5.078 -11.406  15.180  0.31  2.78           N
ATOM     35  CA ATHR A   2       5.757 -11.521  13.850  0.81  3.02           C
ATOM     36  CA BTHR A   2       5.773 -11.153  13.921  0.20  2.67           C
ATOM     37  C  ATHR A   2       7.070 -10.839  14.210  0.74  2.82           C
ATOM     38  C  BTHR A   2       7.155 -10.559  14.193  0.29  1.80           C
ATOM     39  O  ATHR A   2       7.152  -9.941  15.050  0.80  3.31           O
ATOM     40  O  BTHR A   2       7.214  -9.641  15.012  0.25  2.41           O
ATOM     41  CB ATHR A   2       4.976 -10.693  12.813  0.87  5.53           C
ATOM     42  CB BTHR A   2       4.896 -10.354  12.941  0.25 12.07           C
ATOM     43  OG1ATHR A   2       4.611  -9.432  13.388  1.00  6.88           O
ATOM     44  OG1BTHR A   2       3.743 -11.083  12.501  0.25  9.57           O
ATOM     45  CG2ATHR A   2       3.858 -11.584  12.293  0.75 10.03           C
ATOM     46  CG2BTHR A   2       5.683  -9.885  11.726  0.27  5.90           C
ATOM     47  H  ATHR A   2       4.547 -10.814  15.527  0.75  3.44           H
ATOM     48  H  BTHR A   2       5.510 -10.211  15.754  0.25  2.90           H
ATOM     49  HA ATHR A   2       5.962 -12.339  13.548  0.75  3.32           H
ATOM     50  HA BTHR A   2       4.036  -9.929  13.477  0.25  2.86           H
ATOM     51  HB ATHR A   2       5.648 -10.589  11.938  0.75  5.43           H
ATOM     52  HB BTHR A   2       4.644  -9.326  13.574  0.25  5.67           H
ATOM     53  HG1ATHR A   2       5.030  -9.344  14.216  0.75  8.74           H
ATOM     54  HG1BTHR A   2       3.236 -11.198  13.399  0.25 10.21           H
ATOM     55 HG21ATHR A   2       4.096 -12.441  11.924  0.75 10.92           H
ATOM     56 HG21BTHR A   2       6.542  -9.278  12.024  0.25  9.66           H
ATOM     57 HG22ATHR A   2       3.222 -10.974  11.650  0.75 10.92           H
ATOM     58 HG22BTHR A   2       5.039  -9.142  11.179  0.25  9.66           H
ATOM     59 HG23ATHR A   2       3.163 -11.738  13.200  0.75 10.92           H
ATOM     60 HG23BTHR A   2       5.904 -10.639  11.169  0.25  9.66           H

We see that there are many alternative atoms in the 2nd residue.
How should I choose the atoms?
Should I just randomly choose one atom from the alternatives? How can my parser learn to differentiate between those two atoms (say, 33 and 34)? Coz, from the parsing point of view, there is no indication that they are related.
How can I parse alternative atom information in a PDB file?

Generally alternate location are used in structures with partial occupancy.
In https://www.wwpdb.org/documentation/file-format-content/format33/sect9.html you will see that the fourth numerical column is occupancy. So you could take the higher occupancy version.
The tool PDBfixer can do it for you: github.com/openmm/pdbfixer

Do not assume the sum of occupancy is 1 —it should but errors do creep up especially for close alternate location that overlap.
The anisotropic and isotropic B-factors will be lower than were the structure refined with a single conformation (like normal atoms).
Same letters in different residues should not be assumed to be paired up. Models as seen for NMR structures are a solution for that.

Lastly... Your structure has protons, which is unusual (exceptions exist). It is fine if it is a neutron diffraction structure or NMR. If you are not expecting protons (X-ray structure), skip it out of safety.

not sure your answer is the best one, maybe the question was not precise enough, lets says it was: " Is there any PDB parser out there that could differentiate between alternate conformation ? " What would be your answer ?

On a similar note is there any PDB viewer that can easily be told to show one of the possible altLoc conformation ? Easy enough to duplicate by hands structure where you have alternate loop conformation, but what about single aa in differnet places, so that you could have protein AAA, BBB, ABB, AAB and so on ... I tried to google pymol command related to this  but seem overly complicate to me a like you can select either A or B conformation for the full sequence , but maybe I am wrong and YOu can point to nice examples

3nir.pdb Crambin     Method: X-RAY DIFFRACTION
Resolution: 0.48 Å
R-Value Work: 0.127
R-Value Observed: 0.127

PDBfixer can do it: https://github.com/openmm/pdbfixer

Cool, spent hours on google .... guess biomformatics SE is a better expert source, keep it up


text:
Caching and parallelization
Currently algorithmically developing a calculation for "co-dependent" mutation (quotation appropriate). The algorithm uses chunked biopython alignment objects which are parallelized.
Aim
One optization is caching a database within parallelization, but I ain't sure whether the parallelization duplicates the cache.
Method
Each alignment chunk (biopython object) is passed onto a separate CPU and refers back to the original database to obtain its metadata in the final stage of the calculation (dragging metadata around this algorithm makes no sense).
The idea is to cache the database to prevent pulling it in from the harddrive:
from multiprocessor import Pool
from myOOPlibrary import mystuff

Facade design, the parallelised index to the alignment chunks enter into the myOOPlibrary:
from cachetools import cached, TTLCache

cache = TTLCache(maxsize=1, ttl=500)

Class mystuff (stuff2):
... # init omitted

    def mymethod(self, ..):
    # to be updated

Class stuff2 ():
... # init omitted

    @cached(cache)
    def getDB(self, ...):
.... #  pull it in (harddrive) and parse to a dictionary
    return metadataDB

Question
I'm wondering if caching is the right approach in parallelization, specifically whether this will end up with multiplying the database held in RAM...
Alternatives
The alternativees to caching are:

IO to JSON (not cool)
IO via pickle (binary IO): a JSON will readilypickle

The possible advantage of the alternativees are the database is written once and accessed by all cores for certain.
Underlying question

When caching is the database  cached for each core? If so the RAM overhead becomes an issue because the database is duplicated total CPUs - 1 times
Alternatively is the cache shared between CPUs centrally and each CPU refers to this location?

As currently written, this might be more appropriate for Stack Overflow

Possibly but parallelization of biopython objects must be fairly common in algorithm/pipeline development. Alot work with biopython here

What are you using for the parallelization part?  If it's just the built-in threading module, I don't think anything will be duplicated there, you just have to make sure you don't step on your own toes with any read/write aspects (e.g. use something like Queue objects if needed).  But in that case, why using a caching library at all, versus just a single shared DB object you load once?  (More of the rest of the code might help clarify.)

Updated, @Jesse comment is cool and closes in on an easy solution.

If python code is the bottleneck, rewriting in C/C++ may get you a 50-fold speedup.

Thanks, full Cython coding is the next development (C within Python and compiling into C). You're 100% right @user172818 what I've not described is there's extensive looping and that will be speeded up via C. For numpy there's no gain in Python (which I'm sure you know). Coding in raw C ... lots core-dumps and not quite so enthusiastic.

Summarising

@Jesse's comment I think is the right answer, open the database outside parallelisation. What I didn't thoroughly describe (and what was essentially requested) is the architecture, the alignments entering into an entire OOP library and the whole library is parallelised, i.e. every new alignment entering library is on its own CPU. The answer is obvious: instantiate the OOP library with the entire database and thats a key purpose of the design.
The technical answer is joblib. This is built to cache across parallelised archetechture as the joblib team point out here. The criticism is that there's parallelisation within parallelisation. The discussion* is as much in the comments as in the response.
@user172818 response is correct too, the algorithm would benefit from C. Its a lot of code to rewrite in C, but sidestepping via Cython is a future development (its a hybrid language connecting C and Python).

*, https://stackoverflow.com/questions/25033631/multiple-processes-sharing-a-single-joblib-cache

There was a comment hinting this is Stack Overflow, that's difficult. Ultimately the technical answer (which is Stack Overflow) isn't the correct answer as far as I understand OOP.

Feedback
Just to wrap this, point 1. worked fine, thank you. It means a database sits on each processor in parallelisation, but is imported  once. In truth, it worked nicely and provided cleaner code. Re-writing to instantiate once would affect other stuff so would be more costly than a cost-saving.
Python can be very fast (e.g. numpy), but also really slow (hence it has critiques), extensively parallelising helps bridge the gap.
Python parallelisation has critiques which I strongly, strongly believe are unfounded: I parallelise loads via multiprocessor library no problems.
Caching here, i.e. parallelising across a parallelisation could stretch things too far.
Last point
Point 3, re-writing in Cython/C, is a greater priority than the testing required for parallelised caching: can't parallelise all data and that code is better as C.


text:
Retrieve protein sequence from Mgnify given only accession code
I only have the accession codes of several proteins from the MGnify database (https://www.ebi.ac.uk/metagenomics/). I would like to retrieve the full amino acid sequence data from the database, but I have not been able to find a way to do this using the web server or available API.
The website provides sequence search with the HMMER tool (https://www.ebi.ac.uk/metagenomics/sequence-search/search/phmmer). However, I need the sequence data to search with HMMER.
How can I get the protein sequence only with the accession code?
P.S: An example of the accession code is MGYP001242958860.


text:
What tool I can use to map short-reads sequences to reference genome and get specific mapped size
I have about 300 90-bp sequences which I would like to map to a reference genome to make it longer to a 300bp sequence, wherein the 90bp is at the middle.
Anyone knows what bioinformatics tool I can use for this? Thank you!
Actual scenario:
The SNP sequences are from a SNP panel data (90bp sequence, SNP is at the middle). We would like to extend the sequence to 300bp, retaining the SNP at the middle so we can make a sequence capture panel of the SNP.

You mean the target capture panel of SNP 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.

@envs_h_gang_5 I'm actually not very familiar, but something like 'genotyping by sequencing' kind of panel.

If you have a SNP panel then you have the coordinates on the genome. You can use tools such as bedtools slop to expand the coordinates to the desired length and then extract those regions from the target genome using bedtools getfasta.

Thank you! That's very helpful. The dilemma, however, is we'd like to extract the 300bp sequences from a different genome, and not from the reference genome of the SNP panel, therefore the coordinates on the SNP panel will not be exactly on our genome.

Then you need liftover https://genome.ucsc.edu/cgi-bin/hgLiftOver
I don't know if the genomes you're working on have mapping files though.


text:
Find enriched processes in modules obtained by WGCNA
I am doing a gene enrichment of different clusters obtained by WGCNA, but there are some of them in which I don't get enriched processes, is there any way to find processes to these modules, for example to associate them to orthologues of other organisms. If there are different ways, could you send me links on how to do it?

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

Hi María José. What kind of gene enrichment are you using? Which databases are you using to infer functionality? What organism are you analyzing? If it is human I doubt there are other organisms more annotated than that.

Well, not all modules will have significant enrichment in gene sets such as GO terms or KEGG pathways. In my experience, the smaller modules are often not significantly enriched. They could be the result of some technical variation, they could contain genes localized to a particular genomic region (do you also calculate enrichment in genomic regions?) or they could be a particular gene family. WGCNA modules often are enriched in known pathways and processes, but that does not mean that all will be.

thank you so much for your clarifications @Peter Langfelder , I liked the idea of calculating the enrichment by gene regions, which packages would you recommend for this.


text:
multi-sequence alignment of samples with multiple contigs each
I want to perform a multi-sequence alignment on 12 samples that clustered based on cgMLST. Ultimate goal is to find out whether they differ by the presence of a few genes.
I have performed multi-sequence alignment in the past using MAFFT, the main difference being that, back then, I was aligning protein sequences, and now I want to align multiple denovo assembled genomes, with multiple contigs each
I first had a look at MAFFT, the program that I used in the past for protein multi-sequence alignment, though I noticed that it considers each contig in a denovo assembly as independent sequence, which is not what I want
Then I had a look at Mauve, but was discouraged by installation trouble-shooting. Also I read on their webpage:

When an individual file contains several sequence entries, they will
be concatenated and the whole concatenated sequence will be aligned to
the sequences in the other files. This behavior allows
multi-chromosomal genomes to be aligned by including all chromosomes
in a single sequence file. Similarly, incomplete genomes that consist
of several sequence contigs can be aligned, though beware that
incorrectly ordered contigs will appear as genome rearrangements in
the Mauve viewer.

So I guess concatenating all contigs in a single fasta file is not an option, as artefacts (genome rearrangements) might appear
It turns out that MAFFT has also an experimental feature to multi-align fastq reads, which sounds attractive to me. Therefore I would like to ask if anyone can recommend this new/ experimental feature. Also if this option is available on the command-line tool, as the mafft help is quite limited (it does not specify the format of input and output files):
  MAFFT v7.505 (2022/Apr/10)
  https://mafft.cbrc.jp/alignment/software/
  MBE 30:772-780 (2013), NAR 30:3059-3066 (2002)
------------------------------------------------------------------------------
High speed:
  % mafft in > out
  % mafft --retree 1 in > out (fast)

High accuracy (for <~200 sequences x <~2,000 aa/nt):
  % mafft --maxiterate 1000 --localpair  in > out (% linsi in > out is also ok)
  % mafft --maxiterate 1000 --genafpair  in > out (% einsi in > out)
  % mafft --maxiterate 1000 --globalpair in > out (% ginsi in > out)

If unsure which option to use:
  % mafft --auto in > out

--op # :         Gap opening penalty, default: 1.53
--ep # :         Offset (works like gap extension penalty), default: 0.0
--maxiterate # : Maximum number of iterative refinement, default: 0
--clustalout :   Output: clustal format, default: fasta
--reorder :      Outorder: aligned, default: input order
--quiet :        Do not report progress
--thread # :     Number of threads (if unsure, --thread -1)
--dash :         Add structural information (Rozewicki et al, submitted)

Any other recommendation on another strategy to align multiple denovo assembled genomes with multiple contigs each are welcome

Thanks its a very good question, for alignments MAFFT in very recent history was the optimal aligner on 2 criteria.

Were you able to resolve your issue? And which tool did you use? I have a similar problem, though the number of sequences I have to align are in the thousands.

@Raghuram in the end I just loaded the .bam files into artemis and inspected the regions of interest. Looking at the depth, I could inspect the presence/ absence of the genes I wanted to detect.

MUMmer4 is a versatile alignment tool for DNA and protein sequences. It supports one reference genome and up to 32 query genomes. MUMmer4 will align every contig in each genome to the reference genome.

do you concatenate all query genomes in one fasta file? If so, how does MUMmer identify which contig belongs to which genome? Based on a header identifier, for example?

@BCArg I did not concatenate all query genomes in one fasta file.

@ForrestVigor does it support anything like nucmer -p <prefix> ref.fa qry1.fa qry2.fa ... qryn.fa ? on the help page it reads  specifies the multi-FastA sequence file that contains the query sequences, to be aligned with the references

@BCArg Yes. You can use mummer ref.fa qry1.fa qry2.fa ...

If you have a reference genome (or are willing to designate one of your de novo assemblies a reference), you may find QUAST helpful.
QUAST will perform an alignment of all genomes against the reference genome and use that as the basis for comparisons. I believe that under the hood it is running minimap2 for this purpose.
It includes a large-scale alignment viewer ("Icarus") to investigate possible differences within contigs, computes high-level statistics on genome completeness/etc, and has other handy features.
You can see the most recent QUAST paper or the manual for more information.

thanks for the comment, though I am looking for a program to perform the multi-sequence alignment, not to visualise the alignment, for which purpose QUAST is indeed excellent

@BCArg QUAST will generate the alignment. Updating answer accordingly.

@Maximillian Press you are right, and obviously it makes sense that Quast needs to perform the alignment to evaluate/ compare two assemblies. Will explore it further

in the end, I just loaded the .bam files into artemis and, by inspecting the depth (heat-map), I could check which samples had the genes I was looking for:

It won't work on more than 3 genomes at a time, but SynMap3D from the Comparative Genomics webserver (https://genomevolution.org/coge/) should do the trick. I'm mentioning it here because the advantage, for me (in contrast to, say, Mauve), was user-friendliness.
Small note: Be mindful of what you name your uploaded samples. You can set your sequence not to be publicly available, but if you give it a new name, as opposed to a species name already in the database, that will be visible.

thanks, though I am looking for a MS aligner that supports > 3 sequences


text:
How good does the assembly of an NCBI prokaryotic genome have to be in order to argue gene loss?
NCBI has several labels for assembly completeness - Complete, Scaffold, Chromosome and Contig. Complete would be a circularized genome (or linear, rarely)
For a Complete genome it's fairly straightforward - if not present, a gene can safely assumed to be missing (or possibly present on plasmids, if those aren't available in the database).
But what about the others? With Contig-level genomes, I'd feel uncertain about making any statements (especially since some of these projects are split into thousands of contigs). But what about Scaffold, or Chromosome? Are they expected to be complete in terms of genes, but just missing repetitive sequences?

I would suggest that in the absence of "Complete" no very confident assertion can be made about gene presence/absence in a genome assembly. You need to actually  have the chromosome sequence to make strong statements.
I think that most people (me included) would reluctantly accept that genomes with very good BUSCO or CheckM numbers (depending on your organism's taxonomy) are probably more or less accurate in terms of presence/absence. For metagenomics, I would extend that to the MIMAG criteria, which are pretty stringent.
These criteria are of course independent of the NCBI categories, which are pretty coarse-grained.
That said, I think that you can make stronger arguments about gene presence/absence from the sequencing reads from which the genome assembly is derived. Find all the reads for your genome, map them to the genome assembly and also to the closest related ortholog of the gene in question. How well is the gene covered? How do the alignments look? Harder for metagenomics, but there are still probably tricks one could use.

Thank you! Regarding your last paragraph, are you saying that it is probable that genes would be lost during assembly / binning (therefore not fully trusting a non-complete assembly), but not so much during sequencing (therefore trusting the raw data)? I imagine that might not hold for low sequencing depths, though. I was wondering, how high do you thing the sequencing depth have to be in order for the raw data to be reliable in terms of gene absence?

@Laura yes, I believe that's accurate. Many things are lost/misassembled during assembly that are obviously present in inspecting raw data. You would be looking for a similar coverage to the rest of the genome, unless you expect some funny business in that regard. The issues that arise are in the case of repetitive or highly homologous sequences, which is more of an issue in metagenomics as you might have separately binned closely related or homologous genomes. For example if you have abundant genome A, that will tend to suck up a lot of reads for homologous sequences in genome B.


text:
Example of a small but moderately complex PDB file?
I am developing a bioinformatics framework only for learning purposes. I want to calculate --

dihedral angles
detect H, E, and C components
detect Hydrogen bonds
etc.

Therefore, I want a PDB file that represents a protein with a small chain but is complex enough for me to experiment with.
Can you suggest such a PDB file?

detect H, E, and C components ? what are they ?

@pippo1980, Helix, Sheet, Coil.

was missing the E. What about Turns  and Loops

@pippo1980, Firstly, I am only talking about secondary structures (so, turns are not discussed). Secondly, Loops and Coils are mostly synonymous.

Turns may not be SSEs (although they are related) but it may be worthwhile detecting them anyway, as it is fairly simple to check for a Calpha-Calpha distance less than 4Å

Although Insulin is a great, well studied example of a small structure it does have the slight downside of not having any sheet structure in the monomer. There is a small cross-chain sheet between residues 24-26 in both the B and D chains (from https://www.ebi.ac.uk/pdbe/entry/pdb/1zni/ ).
Another possibility might be crambin (http://www.cathdb.info/version/latest/domain/1ejgA00) which does have a small antiparallel sheet, although it is tiny at only 2 residues per strand.
A more general approach might be to use something like the top2018 structure set from the Richardson lab (https://zenodo.org/record/5777651#.YtQB59LMJH4) . I used the top500 back in the day, and they are a high-quality set of structures.

Insulin, perhaps?
https://www.rcsb.org/structure/1trz
From here:

Insulin is a perfect molecule for exploring protein structure. It is small enough that you can display all of the atoms and still have a picture that is not too confusing. Human insulin is pictured here, using entry 1trz. This structure, like many other insulin structures, includes a hexameric complex of the hormone, and only one monomer (composed of one A-chain and one B-chain) is shown here. In the structure, you can see many of the key features that stabilize protein structure. Carbon-rich amino acids, like leucine and isoleucine, cluster in the middle of the molecule, forming a hydrophobic core, and the surface is covered with charged amino acids like arginine and glutamate that interact favorably with the surrounding water. Also notice the three disulfide bridges between cysteine amino acids, which stabilize this tiny protein.

Although an excellent, canonical example ... it does have the slight downside of not having any sheet structure. It does have an bit of extended chain (shown as an isolated strand in the RCSB link). Maybe crambin (http://www.cathdb.info/version/latest/domain/1ejgA00) ?

Can you please add this comment as a separate answer?


text:
Which atoms are not found in protein PDB files?
I am developing an educational Bioinformatics framework, I need to know

Which atoms are absolutely not found in any PDB files?

in the following list -

Atomic radii of the elements (data page)

PDB file readability
It is a question of implementation.
For example:
from rdkit import Chem

mol = Chem.MolFromSmiles('[Og]')
pdbblock = Chem.MolToPDBBlock(mol)

import pymol2

with pymol2.PyMOL() as pymol:
    pymol.cmd.read_pdbstr(pdbblock, 'Oganesson')

shows that both RDKit and PyMol can deal with oganesson.
However, this is just handling the atoms in coordinate space. But as soon as we touch it in molecular mechanics way things blow up:
from rdkit.Chem import AllChem

AllChem.EmbedMolecule(mol)  # UFFTYPER: Unrecognized atom type: Og (0)

But generally most molecular mechanics tools struggle with coordinated iron or similar. Take the mangled case of https://www.rcsb.org/ligand/SIR which ought to have been planar...
PDB DB
In the PDB polymeric units and ligands are called chemical components.
They have a 3 letter code. The list is easily parsable.
The problem is that some lists are incomplete: https://ftp.ebi.ac.uk/pub/databases/pdb/data/component-models/complete/chem_comp_model.cif.gz
when read
from pdbecif.mmcif_tools import MMCIF2Dict
models: dict = MMCIF2Dict().parse('chem_comp_model.cif')

has a large number of entries missing because of <reasons>.
So using a more crude approach and simply querying component by component:
with open('chem_comp.list') as fh:  # this is from somewhere else in the ftp
    chem_resns = *map(str.strip, fh),

import requests
import re
elements = set()
for i in range(len(chem_resns)//100):
    response = requests.post('https://www.ebi.ac.uk/pdbe/api/pdb/compound/summary/', ','.join(list(mia)[100*i:100*(i+1)]))
    response.raise_for_status()
    for resn in response.json().values():
        for entry in resn:
            if entry['formula'] is None: # the unknown ligand/ions etc.
                continue
            elements.update(re.findall('[^\d ]+', entry['formula']))

Gives the following elements that are present:
{'Ag',
 'Al',
 'Am',
 'Ar',
 'As',
 'Au',
 'B',
 'Ba',
 'Be',
 'Bi',
 'Br',
 'C',
 'Ca',
 'Cd',
 'Ce',
 'Cf',
 'Cl',
 'Cm',
 'Co',
 'Cr',
 'Cs',
 'Cu',
 'D',
 'Dy',
 'Er',
 'Eu',
 'F',
 'Fe',
 'Ga',
 'Gd',
 'H',
 'Hf',
 'Hg',
 'Ho',
 'I',
 'In',
 'Ir',
 'K',
 'Kr',
 'La',
 'Li',
 'Lu',
 'Mg',
 'Mn',
 'Mo',
 'N',
 'Na',
 'Nd',
 'Ni',
 'O',
 'Os',
 'P',
 'Pb',
 'Pd',
 'Pr',
 'Pt',
 'Pu',
 'Rb',
 'Re',
 'Rh',
 'Ru',
 'S',
 'Sb',
 'Sc',
 'Se',
 'Si',
 'Sm',
 'Sn',
 'Sr',
 'Ta',
 'Tb',
 'Te',
 'Th',
 'Ti',
 'Tl',
 'U',
 'V',
 'W',
 'X',
 'Xe',
 'Y',
 'Yb',
 'Zn',
 'Zr'}

While those that are absent are:
pt = Chem.GetPeriodicTable()
full = set(map(pt.GetElementSymbol, range(1, 118+1)))
missing = full - elements
for m in sorted(missing, key=pt.GetAtomicNumber):
    print(pt.GetAtomicNumber(m), m)

2 He
10 Ne
32 Ge
41 Nb
43 Tc
61 Pm
69 Tm
84 Po
85 At
86 Rn
87 Fr
88 Ra
89 Ac
91 Pa
93 Np
97 Bk
99 Es
100 Fm
101 Md
102 No
103 Lr
104 Rf
105 Db
106 Sg
107 Bh
108 Hs
109 Mt
110 Ds
111 Rg
112 Cn
113 Nh
114 Fl
115 Mc
116 Lv
117 Ts
118 Og


text:
Program for SOV calculation
Are there any programs (web server or installable) out there to calculate segment overlap score (SOV) for protein secondary structures?
Thanks

https://github.com/nanjiangshu/calSOV

Thanks a lot, that works!

Just googled it not sure it calculates the SOV in use today, seems SOV99 to me


text:
BiocManager::install("DiffBind") mergeOne.c:341:3: note: use option -std=c99 or -std=gnu99 to compile your code
I am using
BiocManager::install("DiffBind")

to install the package "DiffBind" on CentOS_7,
I got the error:
gcc -I"/opt/common/CentOS_7/R/R-4.0.5/lib64/R/include" -DNDEBUG -D_FILE_OFFSET_BITS=64 -I'/lila/home/lid/R/x86_64-pc-linux-gnu-library/4.0/Rhtslib/include' -I'/admin/opt/common/CentOS_7/R/R-4.0.5/lib64/R/library/Rcpp/include' -I/usr/local/include   -fpic  -g -O2  -c mergeOne.c -o mergeOne.o
mergeOne.c: In function ‘mungeTarget’:
mergeOne.c:341:3: error: ‘for’ loop initial declarations are only allowed in C99 mode
   for (int i=l;i<l+r;i++) {
   ^
mergeOne.c:341:3: note: use option -std=c99 or -std=gnu99 to compile your code
make: *** [mergeOne.o] Error 1
ERROR: compilation failed for package ‘DiffBind’

Can anyone tell me what should I do?
Thanks

You can set a specific compiler in ~/.R/Makevars file, create the folder and write it down in Makevers:
CC = gcc -std=c99

or:
CC = gcc -std=gnu99

The previous answer is correct, in that the compiler wants to be able to use C99 to compile your library.
However, the solution of putting CC = gcc -std=c99 into your ~/.R/Makevars is a little bit blunt, since everything that you now compile not directly from source is going to use the -std=c99 flag. This is a recipe for some confusing and hard to track down error messages. It's a bit like installing everything in your base conda environment and waiting for everything to break.
It's a safer option to use withr, which allows you to specify compiler flags for individual installations, keeping everything nice and neat. It works a bit like the with (or within) function in R which constructs an environment to execute some command in.
library(withr)

with_makevars(
    c(PKG_CFLAGS = "-std=c99"),
    BiocManager::install("DiffBind"),
    assignment = "+="
)

It can also be performed directly with BiocManager:
BiocManager::install(
    "DiffBind",
    configure.vars = c("CC" = "gcc", "-std" = "c99")
)

Couldn't that be done via the ... argument on BiocManager::install that it is passed to install.packages? BiocManager::install("DiffBind", configure.vars = c("CC" = "gcc", "-std" = "c99") never used it but I think this is the "standard" way.

@llrs ah, didn't realise that. Thanks, will add to answer!

I might be wrong. how the variables are passed but I think there is the option to do that via install.packages and by extension BiocManager::install


text:
Pulling Multiple Files via rsync with WDL
I'm trying to run a simple task to pull multiple files via rsync with WDL. The task should be quite simple but how to pull the multiple files subsequently is giving me issues. Below is a snippet of my task and how it's being called:
task MultiFiles {
  input {
    String wallclock
    String memory
    Int tmpspace_gb
    Int cpu
    Array[String] remote_files #files we want to pull
    Array[File] local_files
    String data_host
    Boolean? checksum_flag = false

  }
    #Optionality for WDL
    String checksum = if checksum_flag == true then "--checksum" else ""

  output {
    Array[File] output_files = "${local_files}"
  }

  command <<<

    set -e

    #rsync command
    rsync -avP
      --log-file=rsync.log \
      ~{checksum} \
      "~{data_host}":"~{remote_files}" \
      "~{local_files}"

  >>>
  runtime {
    cpu: 4
    memory: "2GB"
  }
}

Call
    call MultiFiles {
      input:
        #Task Resources
        remote_files=remote_files,
        data_host=data_host,
        local_files=local_files
}