text:
Alignment tools on desktop for peptide alignment sequence and how to make annotations on alignments look presentable?
I have 10 transcripts where I need to do 5 peptide sequence alignments (2 transcripts per alignment). I need to do this in order to indicate which amino acids lead to which domains.
My questions are two:

Is there a tool I can use from my laptop that aligns peptide sequences? Because results from the clustal omega peptide alignments web service only last two weeks, they're not ideal for replicability in a paper.

The ultimate aim is to use those alignments and domain classifications to a project I am writing and I am trying to make it as presentable and scientific as possible. Any suggestions on how to make peptide alignments and annotations on peptide alignments (which indicate protein domains) look presentable for a paper or if a paper has already shown peptide alignments in a presentable way that would be highly appreciated.

Could you specify what do you mean by peptide, i.e. are we talking about oligopeptides, 5-20 AA, commonly referred to as simply peptides? Or do you mean protein chains. For the former, domain makes no sense, while for the latter, many tools better than Clustal Omega exist (as Omega is for large datasets, most common being Muscle). About annotations on protein, the term features is used and there are a few feature viewers... So you may want to specify webapp, Python or R here.

Hi Matteo, thank you for your response.
I want to differentiate between classical class I and non classical class I MHC molecules in a model organism using well conserved structural features within classical MHC I molecules (eg intradomain disulfide bridges within class I a1 domain). To do this I need to align the peptide sequences of transcripts of that model organism with transcripts of classical MHC I human genes and see which genes from those transcripts have those structural features. So I am looking for a tool to align and annotate those structural features in a replicable manner

Additionally if you have any other bioinformatic ideas /computational analyses on how to find classical MHC Class I molecules on an organism which is poorly annotated that would highly appreciated!


text:
DESeq2 multiple treatments, multiple time points, multiple cell lines
Note: this question has also been asked on Bioconductor Support
I know there are a lot of questions asking similar things here on this forum and I checked the vingette and did a lot of other research but I still have a hard time wrapping my head around it.
In my position, data is unfortunately just put in front of me and I'm asked for pretty pictures. Even making it clear to the wetlab people in the lab that we need replicates is an issue.
I have multiple cell lines, multiple time points and multiple treatments:
Cell lines: CL1, CL2, CL3 Time points: 6h, 24h Treatments: T1, T2, T3, Control
For each cell line and each time point, there are 3 different treatments plus a control. 3 replicates for each sample -> 72 samples
What I want is actually quite simple; I want to measure control vs each treatment at each time point in each cell line. No testing across time points or across cell lines.
My first thought was to separate the data into 6 different data sets (CL1, 6h | CL1, 24h | CL2, 6h | CL2 24h | CL3, 6h | CL3, 24h) and simply do the DE analysis separately, but I read that this is not the way to go. In the end it would also be nice to get a normalized count matrix with all data normalized together, for PCAs and similar.
I hope I made it understandable. How do I design the DESeqDataSet?

As Mike Love alluded to on the Bioconductor pages, it's best to have specialised help for datasets that are this complex; doing analysis on anything which has multiple compounding factors is not simple. It's also a good idea to talk through the analysis before doing sequencing, so you can discuss things like having six replicates per sample (including controls).
In any case, what I did for something similar to this (eventually, after lots of discussion with biologists about possible models) was to combine all the non-tested factor columns into a single factor. Something like this:
meta.combined.df <- meta.df %>%
  mutate(lineTime = paste0(cellLine, "_", timePoint))

And then subset the data based on this statistic:
counts.sub.mat <-
   counts.sub.mat[,meta.df$lineTime == desired.lineTime]
meta.sub.df <-
   subset(meta.combined.df, lineTime == desired.lineTime)

Then use ~ Treatment as the design formula (assuming "Treatment" is a field containing either the treatment type or "control").
But... this is explicitly treating each line/time point as a separate experiment, which may not be what you want. There are complications associated with batch effects and outliers that need to be talked through, because you could end up identifying "differential expression" that is not reproducible.

When comparing one subset of samples to another subset of samples, you generally do not have to make a subset dds object.  It's preferable to leave all the samples in for dispersion estimation.  (Though if you wanted to look for linear trends across two treatments and multiple time points, you will have to either omit the unwanted treatments, or do some trickery with the model matrix).

Feel free to add this information to your answer. In my case, subsetting worked best for the data we were working with.

To compare a subset of samples to another subset of samples, make one column of column data which will distinguish them all the way you want (so probably join cell line, treatment, timepoint), and make that column your design, and then when you call results, you specify the pair you want to contrast, one pair at a time.

source: Stack Exchange

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

sort: votes


text:
How can I use statistics to compare microbial phenotypes?
Note: this question has also been asked on Biostars
I am currently trying to create a theoretical argument that a microbe's phenotype can affect gene expression in their host. I have 5 species of microbes, each with a different COG (Cluster of Orthologous Genes) profile.  I've been using python and a large database (DualSeqDB) of significant changes in gene expression to compare gene expression in a host and the COG profiles of my 5 species.
This is an example of the dataframe I'm working off of:

columns are counts of COGs present in a microbial species and rows are species
I was wondering if there is any way to compare the species and their COG profiles. My end goal is to draw parallels between how prevalent a certain phenotype is and how much that species can affect a host's gene expression.
If this seems infeasible, would a better approach be to BLAST the species and compare their genome sequences phylogenetically? (this would change the nature of my argument, so I'm hesitant)
I'm not the best with statistics so I'm sorry if there's a really obvious answer I'm not thinking of. If anyone knows of any resources I could use to read up on this stuff more on my own, I would really appreciate that as well.

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


text:
Sliding window optimisation
I would like to check transposable elements through a genome with a sliding window approach but I am not sure about which window and step size to use.
Is there a way to optimise the window size and step size when doing sliding window based operations? I mean how to choose which sizes to pick?
Thanks!

This question appears sparse on detail. What genome is this? What database of genomes? Traditionally a lot of work on transposable elements was performed on Drosophila. Also what is the algorithm for identifying transposons, e.g. via LTRs, blasting against a transposon database, or looking for the transposase protein or something else?

This is rabbit genome (Orycun2.0). And basically I would like to check the correlation between transposable elements and nuclear insertions with mitochondrial origin (numts). In the corresponding literature there are informations about the window sizes used but not about the step sizes.

I will reiterate @M__'s point: the sliding window technique refers not to a single algorithm but to a family of algorithms and procedures. It is very difficult for us to give you useful feedback if you don't provide more details on the tool(s) or approach(es) you are using.

I just would like to repeat the same experiment as they did in this great article:
https://www.mdpi.com/2075-4450/11/10/680/htm# in the 4th figure. Where the authors correlated the density distribution of numts and TEs. They provided the window size but not the step size.


text:
Statistical significance of pairwise protein sequence alignment
I have identified a conserved region in a protein I'm studying, by doing multiple sequence alignments against a dozen of its orthologs. I've then found that this region also aligns relatively well with a functional region (binding motif) in a completely unrelated protein, and I'm trying to estimate whether it's reasonable to believe the two regions might be bound by the same protein. blastp doesn't identify this match, would I then assume that this is insignificant?
Let's say I have a 15 aa long peptide sequence X, which I align against a longer 200 aa protein sequence Y. I find some best match (gaps allowed).
What statistical test could I now perform to answer the question:
Is the match I found for sequence X in sequence Y significantly better than what I would expect to find between X and a random 200 aa sequence Z?

Just formalise the answer, its a straight use of the Blast E-value and I wouldn't personally advise against anything beyond that (explained below) within molecular evolution.
So the Google gives the definition

The BLAST E-value is the number of expected hits of similar quality
(score) that could be found just by chance. E-value of 10 means that
up to 10 hits can be expected to be found just by chance, given the
same size of a random database.

Thus the smaller the E-value the better the hit.
Using a local blast server will provide a specific query verse reference E-value.
blastp -query gene.fa -subject genome.fa

The broader question concerns the whether the Blast hit is a genuine homology and that is very complicated. This is simply because amino acid loci could have arisen by convergent evolution rather than direct descendancy which are homologues and this can become controversial.
When trying to assign a function to a amino acid locus this starts becoming relevant because e.g. the selection pressures of convergence could very difference.
Identifying homologues would require a phylogenetic tree which is compatible with the established phylogenetic tree, e.g. for house-keeping genes as a starting point. It gets complex quickly beyond that. Functional domains within the Blast variant (Blast-Psi) could be of use, its field specific however.
Thus the simple answer is to use Blast E-value (bit-score) and it is a universally recognised metric, the reader can then make up their own mind what the value represents. You could simply simple state the definition of the E-value in the report.

Thanks for the recommendations, this is very clear. I was able to tweak the blastp parameters and got an expect value of 0.15, but like you point out I'd have to make up my own mind whether this is relevant or not (a least it's below 1 ;)). Also, in this specific case I'd expect the two regions to have arisen independently via convergent evolution, so no shared phylogeny expected.

Good its sorted, thanks


text:
Removing duplicate FASTA sequences based on headers with Bash
I used the following command to remove duplicate FASTA sequences based on the header sequence:
paste -d $'\t' - - <fastaFileWithNoLinebreaksInSeq | sort -t $'\t' -uk1,1 | awk 'BEGIN{FS="\t";OFS="\n"}{print $1,$2}'

After this command, the output result is shown below:
>_10|_10
AAAGCAAAACACAGAATTGACGGAAAGACTTACGTTATTAAACGTGTTAAATATAATAACGAGTTTGGCATGGATTTTAAAGAAATAGAATTAATTGGCTCAGGTGGATTTGGCCAAGTTTTC
>_10|_11
ATTATGATCCTGAGACCAGTGATGATTCTCTTGAGAGCAGTGATTATGATCCTGAGAACAGCAAAAATAGTTCAAGGTTTGGCATGGATTTTAAAGAAATAGAATTAATTGGCTCAGGTGGATTTGGCCAAGTTTTCAAAGCAAAACACA
>_10|_12
TTTGGAACTCTTTGAACAAATAACAAAAGGGGTGGATTATATACATTCAAAAAAATTAATTCATAGAGATCTTAAGGTTTGGCATGGATTTTAAAGAAATAGAATTAATTGGCTCAGGTGGATTTGGCCAAGTTTTCAAAGCAAAACACA
>_10|_13
ACTTGTAACATCTCTGAAAAATGATGGAAAGCGAACAAGGAGTAAGGGAACTTTGCGATACATGAGCCCAGAACAGGTTTGGCATGGATTTTAAAGAAATAGAATTAATTGGCTCAGGTGGATTTGGCCAAGTTTTCAAAGCAAAACACA
>_10|_14
AACGAATTTCTTCGCAAGACTATGGAAAGGAAGTGGACCTCTACGCTTTGGGGCTAATTCTTGCTGAACTTCTTCATGTTTGGCATGGATTTTAAAGAAATAGAATTAATTGGCTCAGGTGGATTTGGCCAAGTTTTCAAAGCAAAACAC

This is my actual FASTA file: I don't want to lose header information from my FASTA sequence after removing duplicate sequences based on header information.
>ENST00000642616_2|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTAAAATGAAATCACAGGCGCATTACTGCCTTGTCTTGATGAGTCTCGTTTTGAAAATCTTGGAGTAAGCATGAAC
>ENST00000642616_2|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTAAAATGAAATCACAGGCGCATTACTGCCTTGTCTTGATGAGTCTCGTTTTGAAAATCTTGGAGTAAGGAAGTTG
>ENST00000642616_2|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTAAAATGAAATCACAGGCGCATTACTGCCTTGTCTTGATGAGTCTCGTTTTGAAAATCTTGGAGTAAGGGTGGGG
>ENST00000642616_2|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTAAAATGAAATCACAGGCGCATTACTGCCTTGTCTTGATGAGTCTCGTTTTGAAAATCTTGGAGTAAGGTAATTA
>ENST00000642616_2|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTAAAATGAAATCACAGGCGCATTACTGCCTTGTCTTGATGAGTCTCGTTTTGAAAATCTTGGAGTAAGGTCATGG
>ENST00000642616_2|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTAAAATGAAATCACAGGCGCATTACTGCCTTGTCTTGATGAGTCTCGTTTTGAAAATCTTGGAGTAAGTTCCTTG
>ENST00000642616_8|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTATTTGGTGCCTATCCTAAAACTTCCAGTGGGTTTTCTTGCCTTCCAGTGGCCATATAAAGGGCGTCCTGAAAAC
>ENST00000642616_7|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTCATGAACAAGGCTCAGTTATGATGTTTGAGCACCTTATTAATTCTAATGGAATTAAGCCTGTCATGGAACAATA
>ENST00000642616_1|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTGCGCCTGCGCGCGGGTCCGGCGCCGAGGTTCTTGACTGCTGTGCCGGACGCCAGGTGTAGCCATGCAGCGAGCC
>ENST00000642616_9|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTGGACTGCCATCATCTTGGAATCCAAAATAATTTTGATTACAAGCGCTTTATTAAGTTTGCCCGTGTCTGTGAAG
>ENST00000642616_5|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTGGTGGGGTATCTAGCAGGATGTCTAGTTCACGCACTGGGTGAAAAACAACCAGAGCTGCAGATAAGTGAACGAG
>ENST00000642616_4|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTGTAATTAATGATCCTATCCATGGCCACATTGAGCTCCACCCTCTCCTCGTCCGAATCATTGATACACCTCAATT
>ENST00000642616_6|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTGTCATGGGCCATTTTCTCACATGTTTGATGGACGATTTATTCCACTTGCTCGCCCGGAGGTGAAATGGACGATT
>ENST00000642616_6|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTGTCATGGGCCATTTTCTCACATGTTTGATGGACGATTTATTCCACTTGCTCGCCCGGAGGTGAAATGGACGCAT
>ENST00000642616_6|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTGTCATGGGCCATTTTCTCACATGTTTGATGGACGATTTATTCCACTTGCTCGCCCGGAGGTGAAATGGACGGAA
>ENST00000642616_3|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTTTCCTTGGGGGAGAGGAAGAAGCTGCTTAGTTATATCCAGCGATTGGTTCAAATCCACGTTGATACAATGAAGA
>ENST00000642616_3|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTTTCCTTGGGGGAGAGGAAGAAGCTGCTTAGTTATATCCAGCGATTGGTTCAAATCCACGTTGATACAATGAAGC
>ENST00000642616_3|ENST00000642616_10
AAAAAAAATCTAGCATAAAGTTAATTTAGTTGGAGGAAAAGGAGGGAAGTGGATAAAGGGGCAGATGGAGGGAAGTTTCCTTGGGGGAGAGGAAGAAGCTGCTTAGTTATATCCAGCGATTGGTTCAAATCCACGTTGATACAATGAAGG
>ENST00000462861_2|ENST00000462861_3
AAAAAACTCTTCTACAGAAATTACTCTCAAAGAAACCTGAGGATCGACCTAATTTTTCACAGACCTACGGGATGGCATCATCTCAGATATATTTGATAAAAAAG

Can you please clarify in what form/language you need the solution? I see that the solutions provided so far are not in the same form/language as your example.

I don't understand what you actually have and what you need. What output are you expecting from this example? Did the command you used work or not? I tried running your command and it seems to work fine, I don't get the output you show in the first example. What operating system are you using?

If you want to avoid using extra libraries for any reason, you can just use a simple Python script (version 3.6 and above) to do this:
fr = open("dup_test.fasta", "r")
fw = open("dup_edited.fasta", "w")
seq_dict = {}
curr_header = ''

for line in fr:
   line = line.strip()
   if line[0] == '>':
        if line not in seq_dict:
          seq_dict[line] = ''
          curr_header = line
   else:
        if seq_dict[curr_header] == '':
          seq_dict[curr_header] = line
          fw.write(f"{curr_header}\n{line}\n")

fr.close()
fw.close()

You can use seqkit for this purpose.
seqkit rmdup -n seqs.fa -o seqs_without_duplicate.fa

Very cool kit, thanks for sharing!

@Bioinfomrat you're welcome :) Your Script is also a quite nice solution.

One way using awk:
awk '/^>/ { f = !a[$0]++ } f' seqs.fa

Explanation:
The above is basically the same as:
awk '/^>/ { f = !($0 in a); a[$0]++ } f' seqs.fa

That is, for each header line (lines that start with a '>' character), set a flag, 'f' to true if the current line is not in an array, 'a'. Then add the (header) line to the array. Note that in this form, it's not necessary to increment the value for each element - the following would be sufficient:
awk '/^>/ { f = !($0 in a); a[$0] } f' seqs.fa

The default action is to print, and we print lines only when 'f' is true. This is the same as:
awk '/^>/ { f = !($0 in a); a[$0] } f==1 { print $0 }' seqs.fa

Thanks for the explanation! This is a nifty solution. Any idea on its limitations? I can see that it will work as long as the array can be managed in memory, which is quite a small limitation.

No worries @RamRS. Nothing comes to mind. It should work on most platforms. All bets are off if the input is badly formatted like in the example. It's worth mentioning also that AWK arrays are associative, so only unique keys will be stored in memory. This might be an issue if the input contains a large number of sequences and the number of duplicates is low.

Also, note that only the header lines are kept in memory, not the sequence, meaning that it is vanishingly unlikely that you will ever come across a file with enough sequences not to fit in RAM on any modern machine.


text:
How to interpret the Copy Number from pseudocount depmap portal file?
This question has also been asked on Biostars
How do I convert the Copy Number (CN) file into useful information such as HOMDEL and AMP?
This is an already normalized file with pseudo count + 1. In this file there are values from 0 to 1 for each cancer cell line, in wide format.
I have already looked at the forum post about the file, but it's not helpful for me.
This is what the file looks like:
,DDX11L1 (100287102),WASH7P (653635),FAM138A (645520),OR4F5 (79501)
ACH-000621,1.0166833616923812,1.0166833616923812,1.0166833616923812,1.0166833616923812
ACH-002512,2.9659871795507433,2.9659871795507433,2.9659871795507433,2.9659871795507433
ACH-001188,0.9881267634909964,0.9881267634909964,0.9881267634909964,0.9881267634909964
ACH-000471,0.8053837915961594,0.8053837915961594,0.8053837915961594,0.8053837915961594

Thank you for your help

I'm not familiar with this file structure, but have you looked at the forum post about relative copy number linked to in the FAQ?

Since we do not have matched normals, the output is a “copy ratio” or relative copy number. It is relative to the rest of the genome for that cell line. E.g. if the cell line is tetraploid we would not be able to see it from the relative copy number. These values are reported as log2(relative CN + 1) in the portal.

Yes I checked but I still have the same doubts


text:
How should I imagine the pre-indexing of suffix arrays in the STAR RNA Seq-aligner?
I'm using STAR for an internship (fusion-genes in ALL cancer) and I'd like to understand what exactly it is doing. The people at the internship don't know either how it works, given their biology backgrounds. They told me "its enough to know that it works". That was not an satisfactory answer for me which is why I am fighting my way trough the suplementary materials trying to wrap my head around the algo. My background is bio as well but I've spent reading a book about algos and algos in BI this year.
I am trying to understand the section Pre-indexing of suffix arrays from supplementary materials of the Dobin et al. STAR RNA seq-aligner Publication.

While suffix array search is theoretically fast owing to its binary nature, in practice it may suffer from non-locality resulting in persistent cache misses which deteriorate the performance. To alleviate this problem we developed a pre-indexing strategy. After the SA is generated, we find the locations of all possible $L-mers$ in the SA, $L<=L_{max}$, where L_{max} is user defined and is typically 12-15. Since the nucleotide alphabet contains only four letters, there are NL=22L different L-mers for which the SA locations have to be stored.

The authors also state:

Suffix Array (SA) of the whole genome is utilized to find the Maximum Mappable Prefixes (MMP).

I am confused by the use of the term $L-mer$. When googleing I only keep finding $k-mers$ in combination with suffix arrays. I did however find this GitHub page mentioning $L-mers$ in the following way:

This program computes either the suffix array or the L-mer array of a string. The suffix array gives the lexicographical order of all the suffixes of a string. The L-mer array gives the lexicographical order of all the L-mers in the string (sorts the suffixes with respect to only the first L characters).

→ Can I assume that $L-mer$ was used instead of $k-mer$ to stress the lexicographical ordering? But as far as I understood all suffix arrays contain the suffixes' index based on its lexicographic ordering as illustrated by wikipedia's page on suffix arrays:

Do I have to imagine the entire genome being a single string that is converted into one single suffix array? Could you provide sources fit for a beginner with biology background and very little idea about algorithms?

According to BI SE this falls into the category ".... However, if your motivation is “I would like others to explain ______ to me”, then you are probably OK."
From https://bioinformatics.stackexchange.com/help/dont-ask
Which is why I was hoping to find some help here.

Thanks for the clarification. It's fine to ask these types of specific questions, but context helps a lot in providing useful answers. You've mentioned "a beginner with biology background and very little idea about algorithms", which is inconsistent with the level of algorithmic detail you've provided. That's why I asked for clarification.

@gringer Thanks, I see (-:  trust me, algorithms are still a big challenge and mystery to me.

You can take a look at slide 40 of the lectures slides I use when I teach suffix arrays (reproduced here)

An L-mer is a K-mer — STAR builds an explicit lookup table from all length 14 (by default) strings that maps the length 14 string to interval of the suffix array that contains all suffixes that begin with this length 14 prefix.  This means that, to find a MMP, instead of having to do a binary search on the entire suffix array, you only have to do a binary search on the sub-interval that begins with the length 14 prefix that your query string starts with.  If all substrings suffixes began with length 14 prefixes with equal frequency, then this would reduce the size of the interval you have to look at by a factor of $268435456 = 4^{14}$.  Of course, in biological sequence, all prefixes are not equally likely, so the speedup is more variable in practice.  Nonetheless, it makes lookup in the suffix array much, much faster.
Of course, in addition to the core algorithm, STAR contains many heuristics not all of which are likely described in the paper (especially since the tool has been and continues to be under active development and improvement since the paper's publication).  Nonetheless, this is the idea behind the suffix array pre-indexing.  Since STAR still has the whole suffix array, if there is no length 14 exact prefix match, it is still possible to find a smaller match using a more standard binary search in the suffix array.  However, I do not know (and am not sure where it is documented) exactly when/if STAR will fall back to a more naive binary search if a prefix L-mer is not found.  If you are truly curious about this, I'd recommend asking Alex Dobin (the STAR author and maintainer) over on the STAR GitHub repository.  He is generally very responsive to user queries there.


text:
How to draw a 3d structure using 3d coordinates from pdb file?
I am trying to create a 3d structure from 3d coordinated from a pdb file.
I am trying to draw 3d coordinates from pdb file "ATOM 1 N VAL A 14 22.692 -27.591 -39.402" and using webgl javascript.
I've managed to get the 3d coordinates from a pdb file but i don't know the exact logic that is used to create 3d structure from the coordinates found in pdb file.
I tried those coordinates and managed to create a cloud of 3d points but it is not anywhere near the structure that is found on pdb site.
So wondering how the structure is created on the pdb website?

can you please add the code you tried so far?

JSmol sounds like what you want - it'll let you display 3D structures from PDB files onto webpages using JavaScript. If you just want to view the structure rather than display them on webpages, there's many programs that exist already like Swiss-PdbViewer and Pymol. There's even web-interfaces that let you upload a PDB file and they display it for you like Jena3D Viewer.

If you want to share a PDB file, you could try Michelanglo. It allows you to upload a PDB file (among other things) edit a description panel (which can feature special links that control the protein view and representation) and share the link with whomever —without needing any installations on either side.
Here is the page of shared protein: gallery
(disclaimer: I am the admin).
If you are interested in implementing a molecular viewer on your site, I'd suggest NGL. The PDB website uses NGL too, but is recently switching to MOL* ("molstar"), which is a new development, so less documented, but fixes some corner cases and is marginally faster and is a collaboration with a PDBe. NGL does not depend on any technology that will be depracated mind you, so it actually is a very solid choice for a JS molecular viewer.
The PDB format is a working format, not a deposition format, as is mmCIF. As a result NGL and Mol* have some issues when it comes to PDB files, such as lack of secondary structure records (HELIX/SHEET) results in anaemic sheets, fancier linkages (LINK/SSBOND reconds) are not handled only CONECT records.

disclaimer: I am the admin ..... I though you were the creator (!!) of Michelangelo ?

@pippo1980 Ops. Yes, I am: I was trying to be subtle, but at the same time abiding by the transparency rules/guidelines/custom of SO.

I had a look at Michelanglo and as a third  party I can say that it rocks !!!! Very nice tool indeed

Thank you —that is always nice to hear!

As others said, there already exist javascript libraries to visualize a protein. The most complete and recent one would probably be Mol*, which is what is used on RCSB PDB website. Mol* actually works using NGL, a WebGL engine. But probably the oldest library used on RCSB PDB to visualize proteins would be JSMol. You can actually switch between these viewers on RCSB PDB 3D view using "Select a different viewer" in the bottom right (try it out on 1UBQ)
Now to address the first exercise: there are several things going on after getting the atom coordinates. A PDB file may or may not contain a CONECT part, where bonds between atoms will be clearly stated (contrarily to a CIF file which I believe always contain the bonds). Henceforth these packages may or may not have to guess atom bonds. Given that a PDB file also contains the residues of each atom it's fairly easy to do, less so when dealing with small molecules.
However what you see most of the time is a cartoon representation of proteins, where each secondary structure is represented differently. The library will compute the secondary structure of each residue based on the dihedral angles of the backbone, then represent them. As of why they look like this, check the work of Irving Geis and Jane Richardson.
So if you want to have a website showing protein structures I would strongly advise to use an already existing library, although the exercise of coding your own can always be interesting. If you want to visualize structures on your own computer without using a web interface I can only recommend PyMOL, VMD and ChimeraX.

I am not good as you are and needed to use a parser for the pdb file
see https://mmcif.wwpdb.org/docs/software-resources.html:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Nov  7 17:40:03 2020

@author: Pietro

"""

import matplotlib.pyplot as plt  #matplotlib

from pdbx.reader.PdbxReader import PdbxReader #PDB file parser

import numpy as np #numpy

openz = open('./1msc.cif')

print('opened')

# Initialize a PdbxReader object with the input file handle

pRd = PdbxReader(openz)

# Initialize a list to be propagated with DataContainer (and/or DefinitionContainer) objects (of the DataContainer class, which inherits from ContainerBase) parsed from the CIF file, where data blocks map to DataContainer objects

data = []

print(data)

# Call the read(self, containerList) method with your list

pRd.read(data)

print ('data ###################################################################\n' , data)

print('data lenght' , len(data))

# Your list is now propagated with one or more DataContainer objects, which represent data blocks. To get the first data block, just use list notation:

block = data[0]

print(' block ##############################################################\n',block)

pippo=block.getObjNameList()

print(pippo)

# To retrieve a category object, use the getObj(self, name) method

atomsite = block.getObj("atom_site")

print('atom site : ', atomsite)

atom = atomsite.getValue("group_PDB",0)

# print('atom sites id :', atom)

atomlist=[]

i=0
while True:
    atom =  atomsite.getValue("group_PDB",i)
    atomid = atomsite.getValue('label_atom_id',i)
    if atom =='ATOM':
        # print(atomsite.getValue("group_PDB",i))
        if atomid == 'CA':                         #select only carbon alpha (CA) to work with

            atomlist.append([float(atomsite.getValue('Cartn_x',i)),\
                           float(atomsite.getValue('Cartn_y',i)),float(atomsite.getValue('Cartn_z',i))])

        i +=1
    else:
        print('loop ends')
        False
        break

atomlist=np.array(atomlist)

print('end')

print(atomlist)

print('lunghezza : ', len(atomlist))

fig = plt.figure()
ax = plt.axes(projection="3d")

for p in atomlist:       #print CA as dots
    ax.scatter3D(p[0], p[1], p[2])
    print(p)

for i in range(0,len(atomlist)-1):   #connect with a line CA n with n+1
    x,y,z = [atomlist[i+1][0],atomlist[i][0]],[atomlist[i+1][1],atomlist[i][1]],[atomlist[i+1][2],atomlist[i][2]]
    ax.plot(x,y,z, color='red')

plt.show()

my input is https://files.rcsb.org/download/1MSC.cif
Result is just the CAs (carbon alpha) shown as dots and connected to form the main chain:

Didn't manage how to show the dimer though


text:
Fastq: how can I check if they are from DNA or RNAseq data?
I have (gave me) Illumina fastq files, which I want to use for variant calling, and I do not know if they are DNAseq or RNAseq data. How can I check this? I do not have any report or who to ask.
Many thanks

Honestly, ask who made them.  If the reads are from an organism without splicing, it will not be trivial to figure it out with 100% certainty.
If you are calling variants, I don't think it matters much, except that higher organisms have such a small % of their genome transcribed, it should be obvious very quickly if you have the desired coverage across the entire genome, or just the transcripts.


text:
Find species from FASTA files
In a school project we got the task of finding the species from two FASTA files. We got some hints of what to look for, like it is a unicellular eukaryote and that we can use blast. We don't know if the two files are the same species.
So I wonder if there is a standard way of finding species from FASTA files using the terminal?
I have a background in informatics, so I'm not well known with bioinformatics tools.

Are you not allowed to use the web BLAST tool? That's what I use if I need to quickly find out the likely origin of a DNA sequence.
There are command-line ways to do the same thing, but I only use them if I need to search hundreds or thousands of sequences at once.


text:
blast p-error with making directory
I tried running this command python blastp.py CmydPro tlr TLR.txt to extract the query sequence for the target species I'm looking at. However, whenever I try running it always comes up as a blast p error and says there's no directory. I downloaded blastp.py as a file made from my lab's GitHub but I have no idea how to troubleshoot the error since it's not a required directory. I made sure all my directories were correct and that all my files were intact. Any help with this would be greatly appreciated.
This is the error I keep getting:
Traceback (most recent call last):
File "blastp.py", line 53, in
blastp_scout = subprocess.Popen(['blastp', '-db', organism_db, '-query', 'temp.fa', '-out', 'blast', '-num_descriptions', '5', '-num_alignments', '0'])
File "/home/dylan_w190/miniconda3/lib/python3.8/subprocess.py", line 854, in __init__
self._execute_child(args, executable, preexec_fn, close_fds,
File "/home/dylan_w190/miniconda3/lib/python3.8/subprocess.py", line 1702, in _execute_child
raise child_exception_type(errno_num, err_msg, err_filename)
FileNotFoundError: [Errno 2] No such file or directory: 'blastp'

These are all the files I have stored in my directory:
'$RECYCLE.BIN' 'Genome Files' TLR.txt
'ASM 28112' 'Genome Files (1)' Tham
'ASM 28964' 'Genome Files (2)' Vim
'ASM 3597' 'Genome Files (3)' Xenopus
'ASM 45574' 'Genome Files (4)' blast-2.11.0+
AnoCar Git blastp.py
Chrysemys GopAga get-pip.py
Cliv GopAga2.0-Genome human
CmydPro Halia ls
CmydPro_results.txt LICENCE mouse
Config.Msi Lib ncbi-blast-2.6.0+
CroPor Miniconda3-latest-Linux-x86_64.sh ncbi-blast-2.6.0+-win64.exe
'FASTAs (2)' OpenSSH-Win64 pip
Format Ornitho pvi
GCF_015237465.1_rCheMyd1.pri_genomic.fna PelSin pyzo
GFFs Python setup.py
'GFFs (2)' PythonX temp.fa
GRCh38 README.txt tlr_CmydPro.fa
GRCm38 Rnor tlr_ref.fasta
Gallus 'Rnor - Copy' vs_buildtools__576366542.1609221376.exe
GavGan Scripts vs_community__576366542.1609221376.exe
Gekko 'System Volume Information'

There is a windows version. It seems that I already downloaded blast but I need to do this now:
Command 'blastp' not found, but can be installed with: sudo apt install ncbi-blast+
Do you have any idea how to might download this?
I already downloaded blast as a tar file and the 2.11.0+ version. I'm just not sure how to download just ncbi-blast+ on its own.

That looks like you are using Windows. have you asked the developers if the script you are trying to run is even supposed to work on Windows? Is there a Windows version of blast? That can be launched from the command line?

Could you please tell us what operating system you're running (e.g. Windows, Mac, Linux)? If Linux, could you please show the output of the command uname -a?

You don't have blastp in your path. Blast is either not installed on your system or it's just not been added to path right.


text:
multi--platform GEO analysis based in expression of single gene
This question has also been asked on Biostars
I have a challenge in analyzing GSE34941( It has 4 platforms) according to the expression of a specific gene in R. First, which is the best way to group my patients according to the expression of a gene to high- and low-expression? Second, which is the best way to remove the batch effect of multi-platform analysis here?
I really appreciate any help you can provide.

Maybe can try the ComBat to remove the batch effect


text:
How to select the best ligands in a virtual screening matrix
I have the results of a virtual screening experiment using docking simulations with Autodock Vina. The result is a matrix of 7 (proteins) by 28000 (ligands) with the calculated binding energies for each protein/ligand pair:

Ligands - Proteins
Protein1
Protein2
Protein3
Protein4
Protein5
Protein6
Protein7

Ligand1
-5.4
-5.9
-5.0
-5.2
-5.8
-5.7
-4.5

Ligand2
-9.5
-8.6
-8.9
-9.0
-9.2
-7.1
-6.4

Ligand3
-9.5
-8.3
-8.2
-9.9
nan
-6.8
nan

Ligand4
-9.4
-7.8
-8.0
-6.9
-9.5
-1.3
-7.8

Ligand5
-9.3
-6.9
-8.5
-7.9
-9.2
-6.8
-6.9

---

Ligand28000
-6.3
-7.9
-5.4
-2.7
-10.1
-9.4
-3.9

nan results are either failed dockings, or positive values.
The scores are not only affected by how well a ligand docks with each protein, but also bigger (higher molecular weight) ligands tend to show better (more negative) scores. However, this is in part an artifact and it does not imply that bigger molecules are better and smaller ones are worse.

Typically, for a virtual screening with just one protein, you would select the top ranking ligands and that would be it (e.g. keep the top 1% and test them). However, this results in some ligands selected for most proteins because they invariably score better due to their size.
In this case, I could also take into account the variability of each ligand because I have observations for different unrelated proteins (i.e. see if a particular ligand systematically binds with better scores across proteins).
Is there any way to take into consideration both factors when selecting the best ligands?
I thought of calculating Z-scores 1) across rows and 2) across columns, and select the ligands for each protein that have less than -1.65 standard deviations (p=0.05) and that also fill the same criteria across proteins. However, there are only 7 observations for each ligand to calculate the Z-scores, and it does not feel correct to give it the same weight to the Z-scores calculated from the 28000 observations for each protein. Plus, in the end, protein docking effect is more important than the ligand effect.
Any suggestions would be greatly appreciated!
Ligand efficiency looks like a promising suggestion. However, upon trying (calculating ∆G/MW or ∆G/n.heavy_atoms) resulted in very small ligands with mediocre or bad energies going to the top (e.g. urea with -3.1 became a Z-score of -8.5 SDs).

Urea?! —what dataset are you docking and what's the application? I had assumed it was Enamine REAL or Wuxi biologics dataset for purchasing for drug discovery.

If it is to find the native ligand of an enzyme, the task is completely different, but the monetary aspect is worse —say the native ligand is something like phosphoribosyl pyrophosphate and you have 10 candidates that are close, that will cost as much as the GDP of a small state. Therefore ancillary checks are needed like operon guilty-by-association etc.

I seem to have glitched the system by adding an extra comment and my previous comment has shrunk. If it's for drug discovery I'll re type it... :shrung:

Sorry Matteo I don't know where my comment went either! I used the ZINC database of natural compounds here, with the aim of drug discovery. It is expected to find things like urea I guess, and of course any result like this would be manually curated

My comment was about cost and success rate: smaller compounds are cheaper and more likely to be true so LE + a lower cutoff of 150–200 Da depending on the assay or co-crystal soak. So generally smaller compound hits are elaborated upon until one has a lead and jumping to the final product is a waste of cash. ZINC natural products is not make-on-demand & may require a chemist to source the plant (read: import paperwork) and extract the product (read: low yield). Plus they may not be not "drug-like" (cf. Lipinski's Ro5 vs. nature's proclivity towards PO4, CoA or glycosyl groups as a handle).

What I did was to keep any ZINC "natural-products" entry also described as "for-sale", which (in theory...) should be available for buying from distributors, and also discard any that violated 2 or more Lipinski's rules. But yeah, still it left with around 28000 compounds. I will try to see if I can use LE without boosting the "very" small compounds.

Ah. 28k is more than I thought. I find natural compounds awesome to stare at to see their biosynthesis, but I'd worry the routes for expansions may be more isolated than the modular building block approach —say you have a cyclised sesquiterpene (farnesol (3x isoprene) is 222 Da) and you would like to change a substituent it would need to exist in nature or on sale or easily modifiable (Codexis & co. sell plates of exciting enzymes but AFIK that is for production optimisation not exploration).


text:
Generate simulated bam for certain snps
I want to benchmark my DNA sequencing pipeline. In order to do that I want some gold-standard files in every step(.vcf,.bam,.fastq). I want to generate/simulate a bam file of reads for a given set of vcf. How to do that?
The GIAB project has gold standard vcf, bam, and fastq. But they are not perfect for pipeline testing. Their size is quite big. Also, they are real data. My question was how to generate simulated Bam from vcf.

Have you tried making a modified fasta containing the variants and then simulated reads using that. Normally then you'd mix those together with reads simulated against the reference sequence at different ratios.

Actually, I am planning to do that now. The thing is, I couldn't determine which one is the most convenient way. I am planning to do that using Python. Can you write it in the answer and suggest a convenient way to do that?

Hint, you want to make a fasta from a VCF (e.g., https://www.biostars.org/p/360900/) and then use any random read simulator. No need to program anything, that'd be a waste of time.

Sorry, I didn't read your comment earlier properly. Actually, I used many simulators like NEAT, ART, MASON, etc. But, unfortunately, they are too much buggy. If you go to their issues, you will see their problems. That's why I wanted to write my own code in the first place. However, among the simulators, NEAT probably is the least buggy one.


text:
BioPython - Retrieve sequence records from pubmed database
I've been assigned a task to fetch sequence records from "pubmed" database:

fetch 30 record of type "fasta" from (pubmed database, with term hemoglobin AND alpha) as a sequence records

It doesn't make any sense to me!
I know the PubMed database is for publications and citations, but how can I retrieve those records from it? is that even possible?
I tried to search for the IDs with a given term in it and then retrieve the records from another database (like "nuccore"), but I don't know if this is right.
Any ideas?

The PubMed database will not recover any sequence data - that's a database of publications. I think the question was Genbank

Answer from @devon-ryan, converted from comment:
They meant the nucleotide database.


text:
Query on htseq count
This question has also been asked on Biostars
I am trying to run htseq-count for carrying out rna-seq analysis for solanum tuberosum and i used the following command:
htseq-count --format bam --order pos -s no -a 10 -t exon -i gene_name --idattr gene_id SO_8612_L1.bam GCF_000226075.1_SolTub_3.0_genomic.gff > L1_htseq_count.tsv
and im getting an error message saying:
[E::idx_find_and_load] Could not retrieve index file for 'SO_8612_L1.bam' Error processing GFF file (line 12 of file GCF_000226075.1_SolTub_3.0_genomic.gff): Feature exon-XM_015312074.1-1 does not contain a 'gene_id' attribute [Exception type: ValueError, raised in features.py:326]

How to resolve this? and also can gff be used instead of gtf file. what is the standard gtf file used for Solanum tuberosum?


text:
Metagenome simulation from a concatenated FASTA file
I am trying to simulate metagenome sequencing. So I will start with a file with a lot of concatenated genomes. From there, I would like to randomly extract 10,000 sequences of length 200bps.
I got help subtracting a random 200 sequence from a fasta file like this:
>JQ086372.1 Enterobacteria phage HK446, complete genome
ATGGCAGATAAACGAATCCGTTCCGACAGTTCGGCGGCAGCGGTTCAGGCCATGAAAAATGCAGCTGTGGCTTCTCCCGCTTCAGTCTGCCAATATGGATGTGAAACTCGTCGGAAAAAAAGTGGTTTATCGCTATCAACGCAACAGCGAATACGCCGACTTTTCGCAGAAAGAGATTTTTCACCTTAAAGGCTTCGGATTCACCGGGCTTGTAGGCCTGTCAATTTTTTCCACTTATTTAAATGTTTGTAAATTAAGTAAATAACGTTTTTAATCTAAAAATTCCACAATTTCCACAATAACTACTACTACTAAATTTTTTTAATTATTTTATTTATTTATATTAGAAAAGAGAACAATATAAGTTATTGAATTTAATTTACGAAAAAAGTATATTGTGTATGTGTTAAAAACTTCCAAATTCCGATGGAAATTATTTGTAATCAAAATGAGTTAAACAATGCTATACAACTAGTAAGTAAAGCAGTTGCTTCACGACCAACTCATCCAATTCTTGCAAATATACTTTTAACAGCTGACCAAGGAACTAATAAAATTAGTTTGACAGGATTTGATCTTAATTTAGGAATTCAAACTTCTTTTGATGGAACTGTTAAAAATAGCGGTGCTATTACGATACCTTCAAAACTTTTATCTGAAATAGTAAACAAATTAC...(this still goes for a lot more A or C or G or T)

using this:
tail -n +2 file.fasta | tr -d '\n' > newfile
n=$(stat -c "%s" newfile)
r=$(shuf -i1-"$((n-200+1))" -n1)
newfile tail -c+"$r" | head -c200

And as an output I got: GTTAAGGGCGCCGTTCTGCTTGTCTTTGCTGACATGTTTGAACACCGGACGGCACAAAGCGAAGTGCAGCTTTATGAGAATGCAGCCGCCGAACGCATGATGTTCATTCATCGCAACTGGCGCGGTAAATCTGAACCTGAGGAGGGCTCCTGATGGAACCTGGACGATTCAGGCACAGGGTAAAAATTCTCACCTTCACG
But now I want to extract instead of 1 fragment of 200 characters; 10,0000 random fragments of 200 characters using the strategy above. For this, I am looking to exclude all those lines that begin with ">" (p.e. >JQ086372.1 Enterobacteria phage HK446, complete genome).The input file looks like this:
>JQ086372.1 Enterobacteria phage HK446, complete genome
ATGGCAGATAAACGAATCCGTTCCGACAGTTCGGCGGCAGCGGTTCAGGCCATGAAAAATGCAGCTGTGGCTTCTCCCGCTTCAGTCTGCCAATATGGATGTGAAACTCGTCGGAAAAAAAGTGGTTTATCGCTATCAACGCAACAGCGAATACGCCGACTTTTCGCAGAAAGAGATTTTTCACCTTAAAGGCTTCGGATTCACCGGGCTTGTAGGCCTGTCAATTTTTTCCACTTATTTAAATGTTTGTAAATTAAGTAAATAACGTTTTTAATCTAAAAATTCCACAATTTCCACAATAACTACTACTACTAAATTTTTTTAATTATTTTATTTATTTATATTAGAAAAGAGAACAATATAAGTTATTGAATTTAATTTACGAAAAAAGTATATTGTGTATGTGTTAAAAACTTCCAAATTCCGATGGAAATTATTTGTAATCAAAATGAGTTAAACAATGCTATACAACTAGTAAGTAAAGCAGTTGCTTCACGACCAACTCATCCAATTCTTGCAAATATACTTTTAACAGCTGACCAAGGAACTAATAAAATTAGTTTGACAGGATTTGATCTTAATTTAGGAATTCAAACTTCTTTTGATGGAACTGTTAAAAATAGCGGTGCTATTACGATACCTTCAAAACTTTTATCTGAAATAGTAAACAAATTAC...(this still goes for a lot more A or C or G or T)

>NC_010391.1 Salmonella phage Fels-1, complete genome
TGATTGCATAACCACTTAACATCTTGTTTTATCTAAATAAAATTAAGCATGTTATCTTTTTGGGGCACTCCTGGGGCAGTAGATGCCAGTTGTTGATTCAGTATATCTACTTGTGACTGGTTATTATCCCGATTTTTTTAGTTTTAAGGTGTTGACATAGCCATCCATGCTCCATATACTGTATAGACCATCTGAGCGTTTGAGTGCCCCATTTGTGAAGCAATAAAGTTCGGGTTCGCGCCAGCGGCAAGCGCCCAGCATGCACCGATTTTTTTAGTTTTAAGGTGTTGACATTAGGTATGTCGGGACTGGTATGCTTTCCTGTGTCGCAGCCCGGCGCGTCTCAATGCAGATTCCCATATCCTGTTCAGGGAAGGTGCAGCATAATGTGCTTTGGATTTACCTGAGCGATCGGTAATTTGCGGATTGAAGACAAAGGTAATTTTTTCCACTTATTTAAATGTTTGTAAATTAAGTAAATAACGTTTTTAATCTAAAAATTCCACAATTTCCACAATAACTACTACTACTAAATTTTTTTAATTATTTTATTTATTTATATTAGAAAAGAGAACAATATAAGTTATTGAATTTAATTTACGAAAAAAGTATATTGTGTATGTGTTAAAAACTTCCAAATTCCGATGGAAATTATTTGTAATCAAAATGAGTTAAACAATGCTATACAACTAGTAAGTAAAGCAGTTGCTTCACGACCAACTCATCCAATTCTTGCAAATATACTTTTAACAGCTGACCAAGGAACTAATAAAATTAGTTTGACAGGATTTGATCTTAATTTAGGAATTCAAACTTCTTTTGATGGAACTGTTAAAAATAGCGGTGCTATTACGATACCTTCAAAACTTTTATCTGAAATAGTAAACAAATTAC...(this still goes for a lot more A or C or G or T)

>NC_007577.1 Prochlorococcus marinus str. MIT 9312, complete genome
AATTTTTTCCACTTATTTAAATGTTTGTAAATTAAGTAAATAACGTTTTTAATCTAAAAATTCCACAATTTCCACAATAACTACTACTACTAAATTTTTTTAATTATTTTATTTATTTATATTAGAAAAGAGAACAATATAAGTTATTGAATTTAATTTACGAAAAAAGTATATTGTGTATGTGTTAAAAACTTCCAAATTCCGATGGAAATTATTTGTAATCAAAATGAGTTAAACAATGCTATACAACTAGTAAGTAAAGCAGTTGCTTCACGACCAACTCATCCAATTCTTGCAAATATACTTTTAACAGCTGACCAAGGAACTAATAAAATTAGTTTGACAGGATTTGATCTTAATTTAGGAATTCAAACTTCTTTTGATGGAACTGTTAAAAATAGCGGTGCTATTACGATACCTTCAAAACTTTTATCTGAAATAGTAAACAAATTACAATTTTTTCCACTTATTTAAATGTTTGTAAATTAAGTAAATAACGTTTTTAATCTAAAAATTCCACAATTTCCACAATAACTACTACTACTAAATTTTTTTAATTATTTTATTTATTTATATTAGAAAAGAGAACAATATAAGTTATTGAATTTAATTTACGAAAAAAGTATATTGTGTATGTGTTAAAAACTTCCAAATTCCGATGGAAATTATTTGTAATCAAAATGAGTTAAACAATGCTATACAACTAGTAAGTAAAGCAGTTGCTTCACGACCAACTCATCCAATTCTTGCAAATATACTTTTAACAGCTGACCAAGGAACTAATAAAATTAGTTTGACAGGATTTGATCTTAATTTAGGAATTCAAACTTCTTTTGATGGAACTGTTAAAAATAGCGGTGCTATTACGATACCTTCAAAACTTTTATCTGAAATAGTAAACAAATTAC

So for my new output, I will end up with a file with 10,000 subsets of 200 characters that would look like this:
>1
ATGGCAGATAAACGAATCCGTTCCGACAGTTCGGCGGCAGCGGTTCAGGCCATGAAAAATGCAGCTGTGGCTTCTCCCGCTTCAGTCTGCCAATATGGATGTGAAACTCGTCGGAAAAAAAGTGGTTTATCGCTATCAACGCAACAGCGAATACGCCGACTTTTTCAGTCTGCCAATATGCGTTCCGACAGCAGATAAAC

>2
TGATTGCATAACCACTTAACATCTTGTTTTATCTAAATAAAATTAAGCATGTTATCTTTTTGGGGCACTCCTGGGGCAGTAGATGCCAGTTGTTGATTCAGTATATCTACTTGTGACTGGTTATTATCCCGATTTTTTTAGTTTTAAGGTGTTGACATAGCCATCCATGCTCCATATACTGTATAGACCATCTGCAGATA

>3
AATTTTTTCCACTTATTTAAATGTTTGTAAATTAAGTAAATAACGTTTTTAATCTAAAAATTCCACAATTTCCACAATAACTACTACTACTAAATTTTTTTAATTATTTTATTTATTTATATTAGAAAAGAGAACAATATAAGTTATTGAATTTAATTTACGAAAAAAGTATATTGTGTATGTGTTAAAAACTTCCAAAT

Makes sense? Any help?
Thanks in advance

This is a great question for discussion. @gringer  is 100% correct, the question is which and why

It might be better to formulate your question, or append the summation of the goals of your study to match them against a simulator

Given that you want to simulate short-read sequencing from reference genomes, there are read simulator programs that already exist. For example, see this answer.

Its a fantastic discussion


text:
Multiple Alignment cost application
Is there any biological application/case where someone would be interested in the estimated total cost of the alignment between a set of sequences (genes or amino acids) without the aligned sequence itself?
Edit: to clarify by cost I mean the cost function such as the edit distance. I have some ideas for heuristics on the multiple sequences alignment problem (which is intractable) but these would not retrieve the aligned sequence, therefore, I was wondering whether such an algorithm would be useful.

UPDATE
In response to comments, it seems that the question is more along the lines of:
"what are some applications where computing summary statistics associated with a specific MSA result could be useful, even if we can't return the sequences?"
I have 2 responses to this:

in keeping with my original answer, I would say that it would still be useful to obtain some measure of identity, similarity, or pairwise distances between sequences (or any summary statistic about the relationship of the sequences). There are lots of applications for this level of data. I don't know whether this is something that can also be extracted from the computation, I'm a little unclear on what exactly the proposed computation is doing if it's getting a cost but not sequences.

Michael says below (and I agree) that an alignment cost is likely to strongly depend on particularities of the underlying alignment model as well as the data, and things get nasty when you try to account for gaps or weirdness. Nonetheless, it's possible that a cost could be used as a sort of guide heuristic for a separate follow-up alignment procedure, to indicate a portion of parameter space which is likely to yield good results. I don't know what this looks like, but it could be simple as resource provisioning for the computation (e.g. CPUs/RAM/disk). Alternately, maybe the cost result can be a little richer and actually seed some parameters in the next computation.

Again, without more detail I'm not sure what this looks like (and I'm not an alignment expert so my suggestions wouldn't be very good anyways).

The simplest application might be to come up with a quick "check" of the sequences to ask whether they are alignable in the first place. Like, as a QC step (did you forget to reverse complement?) maybe. (this assumes that the alignment in question is pretty onerous and you don't want to have to try several times.)

In summary, a lot depends on how you get that cost and what else you could get in addition to the cost. I wonder if rather than focusing on cost you can output some other richer summary of the alignment (such as e.g. sequence similarities, or gapping stats, etc.).
ORIGINAL ANSWER
The short answer is "yes".
Quite frequently we are interested in measuring "identity" or "similarity" of two sequence strings. As one totally arbitrary example, in metagenomics, quite frequently two rDNA amplicons that show 97% identity are considered to be in the same species, 95% same genus, 98% same strain, etc. So these measures are frequently useful as heuristics.
Identity itself is just a length-normalized edit distance across a sequence, after all.
But more generally, such measures of identity or similarity can themselves be understood as cost measures, in the sense that e.g. a dynamic programming alignment algorithm is trying to find an optimal set of matches between two sequences, which necessarily will maximize identity. Generally what it is actually maximizing is similarity (according to some cost matrix such as e.g. BLOSUM or Jukes-Cantor or whatever).

Thanks for the answer do you have any application in mind where the number of strings under analysis is >2. To be clearer I have some ideas for heuristics on the multiple sequences alignment problem (which is intractable) but these would not retrieve the aligned sequence, therefore, I was wondering whether such algorithm would be useful

@gcorso i see. i share Michael's skepticism below. it seems to me like it will be highly dependent on the specifics of parameterization and the actual sequences, such that it will be low generalizability. possibly it could be useful as an intermediate heuristic computation that is an input to another procedure, but I don't think it would be easy to interpret on its own.

Thank you for the answer, do you have any specific procedure in mind where it could be useful as input?

i'll just edit answer with response.

I'm not a fan personally, but it is context specific. The problem is that it is heading towards optimisation, particularly algorithmic optimisation. If the alignment is the final calculation - that is cool, but the situations are select. It depends how far you want to take it.
The rationale of my skepticism is simple - if an alignment was decided on a 'maximization' criteria, be it likelihood or some other metric - the subsequent calculations based on the alignment, such as phylogenetics Jukes Cantor, Kimura 2 and 3-parameter etc ... we enter circular mathematics. The critiera that defined the alignment is the criteria we are using to perform unrelated post hoc analysis. In truth there is a degree of circular mathematics in any algorithm used to produce an alignment, but we overcome this by optical verification precisely avoid this.
In this sense the alignment is simply to establish homologous amino acid/nucleotide positions for downstream analysis. Thus optical optimisation, whilst cautioned against in every other area of bioinformatics, is considered useful in alignments because of the risk that a logical error in the alignment is amplified in downstream analysis, potentially leading to an erroneous conclusion.
Algorithms are a problem where indels are concerned and aligning across an indel. It is very easy to make errors here and 'optical analysis' would simple introduce a blockwise deletion to avoid the uncertainty, i.e. if you are not sure just chop it out - better to lose data than risk artefacts on an equivocal judgement call.
I'm not saying don't try, I am saying don't look to overstate the applicability of the metric.

Thank you for your answer. What do you mean when you refer to a logistical error? Would it be more useful if it could produce some sort of heat map of the regions in the sequences that are conserved (again without retrieving the discrete aligned sequence)?

Sorry! I meant a logical error. An error that doesn't give a bug but could give an erroneous conclusion

Yes a heatmap is a very good idea. Highlighting areas of mutational saturation would be very useful. The theory is a bit complex, but it is one area - particularly in phylogeny  but is applicable to any analysis - where we unknowingly introduce saturated signal


text:
Gene names from components in NMF analysis
This question was also asked on Biostars
I am using Rcppml package in R for my NMF analysis, I have a matrix from single-cell analysis. I have cells scores and I wanted to know how to extract genes instead of the components.
i used the following codes.
t1.cpm<- LogNormalize(data = as.matrix(t1@assays$RNA@counts)) # Log converting the data form seurat obj

t1.model <- RcppML::nmf(as.matrix(t1.cpm), 30, verbose = F, seed = 1234)

t1.w <- t1.model$w
t1.d <- t1.model$d
t1.h <- t1.model$h

dim(t1.w)
rownames(t1.w)<- rownames(t1.cpm)
colnames(t1.w)<- paste0("component", 1:30)
rownames(t1.h)<- paste0("component", 1:30)

and my output is as follows
        T_1_TACCTATCAGATGGCA T_1_CCTATTATCAGGCAAG
component1             0.000000000            0.0000000000
component2             0.000000000            0.0023590277
component3             0.001512931            0.0000000000
component4             0.000000000            0.0021574559
component5             0.001175799            0.0020222144
component6             0.000000000            0.0004463346

It would be helpful if you can let me know how to extract Gene names instead of components.
Thanks,
Dave.

Rcppml is a machine learning library.
https://github.com/zdebruine/RcppML
NMF is Non-negative matrix factorization which is used in bioinformatics and AI, which can examine the latent relationships in experimental data sets.

NMF of a gene expression matrix into $k$ components works to decompose the expression matrix into an activity matrix ($H$) and latent factor matrix ($W$).
$W$ describes the $k$ latent factors (often called components) by providing a score for how much each gene contributes to each component. $H$ describes how much each component contributes to the expression of each cell.
To extract the activity and component matrices from your NMF model fit using RcppML, use the following code:
# cells is a Seurat object
model <- RcppML::nmf(cells@assays$RNA@data, k = 5, verbose = F, seed = 1234)
W <- model$w # Components Matrix
H <- model$h # Activity Matrix

# Name genes/components for latent factor matrix
colnames(W) <- paste0("component", 1:k)
rownames(W) <- rownames(cells) # Gene names

# Name components/cells for activity matrix
colnames(H) <- colnames(cells) # Cell IDs
row names(H) <- paste0("component", 1:k)

> head(W)
              component1   component2   component3 component4   component5
MIR1302-2HG 0.000000e+00 0.000000e+00 0.000000e+00          0 0.000000e+00
FAM138A     0.000000e+00 0.000000e+00 0.000000e+00          0 0.000000e+00
OR4F5       0.000000e+00 0.000000e+00 0.000000e+00          0 0.000000e+00
AL627309.1  2.927886e-07 9.248074e-07 7.961105e-07          0 2.802758e-06
AL627309.3  3.486284e-08 0.000000e+00 1.916799e-07          0 7.938089e-08
AL627309.2  0.000000e+00 0.000000e+00 0.000000e+00          0 0.000000e+00

> head(H)
                 cell_1       cell_2       cell_3       cell_4       cell_5
component1 0.0001370465 0.0002195795 0.0001628898 0.0002239623 0.0001484157
component2 0.0003604048 0.0002239471 0.0005688147 0.0005242678 0.0004821531
component3 0.0003948121 0.0001485325 0.0002885584 0.0002242272 0.0002865587
component4 0.0001391838 0.0002418495 0.0001361317 0.0001382889 0.0001558255
component5 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000

People will often use the top $N$ genes to describe a component. The following code gets the top $N$ genes for each component and stores them in a list
# This code assumes the above code block has been already run
require(dplyr) #used for %>% piping
N <- 100
top_genes <- list()
for (i in 1:dim(W)[2]) {
  name <- colnames(W)[i]
  top_genes[[name]] <- sort(W[, i], decreasing = T) %>% head(N) %>% names()
}


text:
Is ./1 a valid GT value to add to VCF?
I have a VCF file coming from a structural variant detection tool that doesn't add GT value. I can add GT values, for example by using this tool: https://vatools.readthedocs.io/en/latest/vcf_genotype_annotator.html. I know there exist the alternate allele (1), but I don't know if it's homozygous or heterozygous.
My collaborator wants a VCF with genotype. If an SV is called but not genotyped, I was wondering if ./1 could mean unknown allele and alternate allele for a diploid genome (working with human DNA). Although, if it is not a standard format, it may be a problem.
Is ./1 a valid value to add or does it have to be 0/1 or 1/1? Thanks.

No, ./1 isn't valid. You can find the VCF specifications at https://samtools.github.io/hts-specs/. The current version is VCFv4.4 and this is the genotype section which hasn't really changed much from previous versions as far as I recall:

GT (String): Genotype, encoded as allele value preceded by either of /
or | depending on whether that allele is considered phased. The first
phasing indicator may be omitted and is implicitly defined as / if any
phasing indicators are / and | otherwise. The allele values are 0 for
the reference allele (what is in the REF field), 1 for the first
allele listed in ALT, 2 for the second allele list in ALT and so on.
For diploid calls examples could be 0/1, 1 | 0, /0/1, or 1/2, etc.
Haploid calls, e.g. on Y, male non-pseudoautosomal X, or mitochondria,
should be indicated by having only one allele value. A triploid call
might look like 0/0/1, and a partially phased triploid call could be
|0/1/2 to indicate that the first allele is phased with another
variant in the VCF. If a call cannot be made for a sample at a given
locus, ‘.’ must be specified for each missing allele in the GT field
(for example ‘./.’ for a diploid genotype and ‘.’ for haploid
genotype). The meanings of the phasing indicators are as follows (see
the PS and PSL fields below for more details on incorporating phasing
information into the genotypes):

/ : allele is unphased
| : allele is phased (according to the phase-set indicated in PS or PSL)

For symbolic structural variant alleles, GT=0 indicates the absence of any of the ALT symbolic structural
variants defined in the record. Implementer should note that merging a VCF record containing only symbolic
structural variant ALT alleles with a record containing other alleles will result a change of the meaning of the
GT=0 haplotypes from the record containing only symbolic SVs.

This is the bit you are most interested in:

If a call cannot be made for a sample at a given locus, ‘.’ must be specified for each missing allele in the GT field (for example ‘./.’ for a diploid genotype and ‘.’ for haploid genotype).

So, ./1 doesn't mean anything. You could use ./. which is the standard for "not genotyped", or you could simply not have a GT field at all since that isn't actually required by the standard, so a VCF file with no genotype field at all can be a perfectly valid and standard compliant VCF.
In other words, just think about what you really need. Do you absolutely need a GT field? If not, just leave it out and you're done. If you do, you can use ./. if you want to say "I don't know what the genotype here is" (but be aware that some tools interpret ./. as "this might not be a real call"), but don't use ./1 since that won't be understood by anyone.

./1 is valid - it's contained a lot in the output from decompose operations.

@RamRS I don't see why that makes it valid. It isn't mentioned in the specs as far as I can tell, and I don't see anything in the specs that would support that, so why do you say it is valid? There are loads of tools, including commonly used popular ones, that produce invalid vcf files. To be fair, it is an extremely complex format, and it isn't easy to get it right, but the fact that one tool does something is really no indication that that thing is valid. Do you know what they use ./1 to mean?

I assumed that it does not break the format even if it does not adhere to specs 100%, so it's not really invalid. I think they mean that it's pseudo-haploid. When a multiallelic 1/2 GT is broken into 2 biallelic sites, it does not really yield 0/1 and 0/2 but ./1 and ./2 as each entry is kinda sorta haploid.

@RamRS as someone whose day job involves convincing clients that just because their file has a vcf extension, that don't make it no valid vcf, the standard is all I can rely on. If it doesn't comply with the standard, it isn't a valid vcf file and whether or not software can read it correctly is anyone's guess. That said, bcftools, which is famously strict, doesn't complain about 1/. but other tools might. And since it isn't in the specs, we can only say it isn't valid.


text:
Export multiple variables from shell script into a Nextflow output channel
Within a Nextflow workflow, some process runs a shell script that creates multiple variables, e.g.
var_a=$(echo "1")
var_b=$(echo "2")
var_c=$(echo "3")
var_d=$(echo "4")

These variables contain statistics (e.g. number of input reads, quality scores, ...) that I want to keep until the end of the pipeline, to then save them and perform QC. So I'd like to export them out of this process and carry them to the end of the workflow. To be clear, the number of variables is not unknown, it's just too many for convenience if I need to list them all every time.
I understand that I can use the env qualifier or an intermediate file to export the variables individually:
process get_vars {
  output:
    tuple env(var_a), env(var_b), env(var_c), env(var_d)

  shell:
  '''
  var_a=$(echo "1")
  var_b=$(echo "2")
  var_c=$(echo "3")
  var_d=$(echo "4")
  '''
}

workflow{
  get_vars | view
}

with output:
[ea/a145b4] process > get_vars [100%] 1 of 1 ✔
[1, 2, 3, 4]

However that becomes very impractical if I have a lot of such variables (that I do want to continue passing along to later processes, along with a much longer tuple of other relevant variables).
My question: is there a convenient way to "encapsulate" all these variables in the output channel, for example with a Groovy map similar to the meta map in nf-core? So I could directly write something like:
process get_vars {
  output:
    val(my_vars)

  script:
  '''
  var_a=$(echo "1")
  var_b=$(echo "2")
  var_c=$(echo "3")
  var_d=$(echo "4")
  '''

  my_vars = [:]
  my_vars.a = var_a
  my_vars.b = var_b
  my_vars.c = var_c
  my_vars.d = var_d
}

workflow{
  get_vars | view
}

which obviously fails because var_a is not defined outside of the ''' block.

Can't you simple write all vars to a plain text file and then have some groovy/nf code to loop over the content and pass the it to any data structure that is appropriate, be it a map, list, anything?

That does sound like a solution, at least I would only have to pass along a path. But that means parsing again a file of variables I have already parsed, which doesn't feel like the most direct approach. If you post it as an answer I would accept it.

I think you can use output for this, ouput: env varname, as is mentioned here https://www.nextflow.io/docs/latest/process.html#outputs

I am not sure if i get what you mean by encapsulate, you can give your tuple a name using emit?
process get_vars {
  output:
    tuple env(var_a), env(var_b), env(var_c), env(var_d), emit: vars

  script:
  '''
  var_a=$(echo "1")
  var_b=$(echo "2")
  var_c=$(echo "3")
  var_d=$(echo "4")
  '''
}

workflow {
  get_vars()
  get_vars.out.vars.view()
}

That gives:
 [1, 2, 3, 4]

If I'm not mistaken, using emit (described in the DSL2 docs) is to create a separate (named) channel, so I still have to ensure that it keeps the order in sync with other channels, right?

Here is an illustration keeping in sync with a main channel:
`process get_vars {
input:
val(main_chan)
output:
tuple env(var_a), env(var_b), env(var_c), env(var_d), emit: vars
val(main_chan), emit: main_chan
shell:
'''
var_a=!{main_chan}
var_b=$((!{main_chan} + 1))
var_c=$((!{main_chan} + 2))
var_d=$((!{main_chan} + 3))
'''
}`
`process next_process {
input:
val(vars)
val(main_chan)
output:
tuple(val(vars), val(main_chan))
script:
'''
'''
}`

workflow {   main_chan = Channel.of( 1, 2, 3 )   get_vars(main_chan)   next_process(get_vars.out.vars, get_vars.out.main_chan)   next_process.out.view() }

The key here is that next_process(get_vars.out.vars, get_vars.out.main_chan) needs to receive both channels from the previous process at once. Feel free to include this in your answer if you think it would be relevant

I see, you can use a id (like sample-name or meta.id in nfcore) that you include in both channels and then use to join them together into a single channel to make sure that information for one sample moves though the pipeline together.


text:
How to parallelize BEAST
I'm running BEAST v1.10.4 with BEAGLE v4.0.0. My dataset has one large partition and I have access to a machine with many CPUs so I would like to optimize performance by parallelizing the process. Following the documentation and various (quite old and probably outdated) help Q&A pages I tried
beast -beagle -beagle_CPU -beagle_instances X -beagle_thread_count Y -threads Z mydata.xml

with different values for XYZ but BEAST is ever only using a single CPU regardless.
How do I make BEAST use on more than one CPU?
We have a 'server' in-house, that is essentially a big computer with 144 CPUs that is shared by all users. That's why it would be important to limit the max amount of threads with the -threads option. No queuing system, no GPUs. I want to do some tip date sampling which afaik is not available in Beast2.

The machine is important because of the introduction of M1 type chips which will run GPU.
From historical use of beast1 the following will work ...
beast -overwrite -beagle_auto -seed 2001 mydata.xml > myfile.runtime.log

This will parallelise happily within the number of processors available. I would recommend it. This is because it may not parallelise across 144 CPUs. It depends on the machine architecture. Also from memory (mine not RAM) it dynamically readjusts, so even if it was 144 CPUs parallelised when another job lands on a CPU that drops to 143 CPUs for beast1.
Anyway, I use the Cloud so I no longer keep up with computer architecture.
If this floods 144 processors and there's no dynamic readjustment I would suggest,
 beast -beagle -beagle_SSE -beagle_CPU -threads 10 -overwrite mydata.xml > myfile.runtime.log

I would limit the number of threads to get it working before considering e.g. 40.
However, shifting to GPU is strongly recommended, this will solve any beast2 issues and its very fast. However, I don't know about beast1 for GPU, I've never used it for that.

If none of the above work and you've not get any sort of M1 or higher chip (to move to GPU) then you need to look at the error log, especially if you see this error.

Error 122 parsing the xml input file
Cannot create class: ThreadedTreeLikelihood.

Note although you've got a nice CPU heavy machine, simple fairly modern Mac will do the same operation probably quicker via GPU, at least for beast2, beast1 I dunno

Just to note I'm acting in solidarity with an SE strike until an internal matter is resolved about ChatGPT and requested a meeting with staff to explain why. I've formally notified SE. I'll be answering very select questions and doing only essential admin  (breaches of CoC and spam). Hopefully all is back to normal shortly.

Thanks! Unfortunately I couldn't get either to work (the first suggestion doesn't run at all), however, I followed your general advice to use more/different resources to speed up beast worked for me.


text:
I have really skewed RNA-seq data, what's the best way to normalise it? Preferably in R!
My data frame compares the RNA-seq reads from many genes in different tissues.
The reads look as follows.

I tried using log to make it better but still looks pretty skewed.

Are there any better transformations/functions for this in R?
This is part of an assignment I have, I have used in the past edge software and agree it will do a better job but in this work I'm not meant to use any of the pipelines.

I assume that the eventual goal is differential expression or something along those line, in which case you're much better off using one of the many existing Bioconductor packages.

I like using DESeq2. There's a great document written by the developers about how to process gene count tables and look at differential expression:
https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
The document also mentions other tools that have similar approaches:

Other Bioconductor packages with similar aims are edgeR, limma, DSS, EBSeq, and baySeq.

Another approach, that does not replace the methods already mentions, but could be useful if you don't have raw counts to start with, or if you want to screen hundreds of genes, and you want to focus on the most interesting.

either log transform or zscale your expression values
do a boxplot to see how your gene is expressed in the different tissues
run a non-parametric test to compare the distributions. I usually use a kruskal.test with multiple comparisons and a Mann-Whitney U test (wilcox.test) if I only have two groups.

This should give you an idea if the expression of a gene is different between tissue types.

Thanks for all your answers.
Although the pipeline packages would have been useful I was trying to normalise the data before analysing it later and so I used the bestnormalize package which looked like it did a pretty good job.

What is your goal with this analysis? Your RNA-Seq data looked relatively typical at a glance and I would really recommend not using a general normalization as you did. A lot of research has gone into normalizing RNA-Seq and you are at best hurting your analysis, at worst doing something outright wrong. Use - as recommended before - edgeR, limma or DESeq2 unless you have a good reason not to.

See e.g. Soneson & Delorenzi, 2013

This is part of an assignment I have, I have used in the past edge software and agree it will do a better job but in this work I'm not meant to use any of the pipelines

@Connor you should edit your question and include that information, so that people trying to help don't offer you a solution you can't use.


text:
How to download ligands for PDB structure
I have some PDB IDs and for every structure, I need all its ligands, so I want to automate the process.
The ligands are different for every chain:
$ cat ./7rdx_N_1N7.smi
C1[C@@]2(C)[C@@H](C[C@@H](C1)O)[C][C]([C]1[C@H]3[C@]([C@H](C[C]21)O)([C@H](CC3)[C]([C])[C][C])C)O   1N7
$ cat ./7rdx_U_1N7.smi
C1[C@@]2(C)[C@@H](C[C@@H](C1)O)C[C@H]([C@H]1[C@H]3[C@]([C@H](C[C@H]21)O)([C@H](CC3)[C@H](C)CCC(=O)NC[C][C][N+]([C])([C])[C])C)O 1N7

So, I cannot just download the ligands by its name (e.g. https://files.rcsb.org/ligands/download/ZN_ideal.sdf).
I've tried to find the structures in .pdb file for a given structure, but I see only their names:

Also, as you could see above, I'll convert the .sdf / .mol2 files to SMARTS which seems to be a simpler representation - maybe this information'll make the task easier.

I want to make a script for the process in Python, and the only way to access the RSCB PDB with a library I've found was using biopython but I found there an option to download just a .pdb file (Bio.PDB.PDBList.retrieve_pdb_file), not ligands.
I've also tried to make my own script, but still I don't know how to access ligands files automatically for every structure. The .pdb files I can access simply by ftp (e.g. https://ftp.wwpdb.org/pub/pdb/data/biounit/PDB/divided/rd/).

How to save each ligand from a PDB file separately with Bio.PDB?  ---> https://stackoverflow.com/questions/61390035/how-to-save-each-ligand-from-a-pdb-file-separately-with-bio-pdb

have a look at ModelServer API too https://models.rcsb.org/#/ and for Pyhton https://github.com/OpenAPITools/openapi-generator

there is a (new for me) https://data.rcsb.org/graphql/index.html GraphQL API : The GraphQL interface offers more flexible data retrieval, essentially making it possible to grab any piece of data from any level of the hierarchy in a single query. To use it programmatically POST your GraphQL queries under the data.rcsb.org/graphql endpoint. No idea how it works but should be the one you are looking for to download all Ligands ??

dont understand what does:  " /v1/{id}/ligand
Coordinates of the first group satisfying the given criteria."  means for the ModelServerAPI https://models.rcsb.org/#/

Thank you for your response. It seems they meant first ligand for the given structure, but not unique one. I tired it on 7RDZ and got ASP which is not an unique ligand for the structure.
Still trying to get all unique ligands.

Played a little bit on modelserver with no luck, graphQL should be easier because the test page has autocompletion  but wont work with OpenAPI generator for Python bindings. Ligands used to belong to a Chain in the old PDBs standard now they are non polymer entity or chemical compound so buffer components get mixed up with 'proper' ligands. Need to figure out how to get these entity out from API just sending entry ID and then recover them one by one (for all the Asymmetric units !?) ; not sure it'll be a,less hassle then parsing the full PDBx/mmCIF file

ADDENDUM https://www.apollographql.com/blog/graphql/python/complete-api-guide/ : Using GraphQL with Python – A Complete Guide

Thanks, I'll try to use a query I've made earlier in GraphiQL. It seems to provide the ligands in .smi format right away.
Still, I'm curious if it's possible to get the files from the webpage (.sdf, .mol2) with any GrapgQL query.

Found this on rcsb pdb : Ligand of Interest (LOI)
A small molecule is termed a Ligand of Interest (LOI) when it is the subject of the author’s research rather than a byproduct of the experimental conditions. For example, in entry 1OBV, the small molecule palmitic acid. (PLM) is designated as an LOI as it is bound inside the cavity of the the fatty-acid-binding protein EgFABP1 and is integral to the study. Now need to figure out where LOI is written

Pdb 101 has a typo example pdb is 1O8B not 1OBV

More, I am learning mmcif :  residues/ligands do not have a residue number in mmCIF

my memoruy is failing please find this post I tried to answer about 1 year ago https://bioinformatics.stackexchange.com/questions/18076/help-me-in-understanding-the-pdb-file

Here can use GraphiQL API provided by RCSB.
Here is an example
import requests
from rdkit import Chem
pdbid = '7RDX'
ls_smi = []

query = f'''{{
  entry(entry_id: "{pdbid}") {{
    nonpolymer_entities {{
      nonpolymer_comp {{
        rcsb_id
        pdbx_chem_comp_descriptor {{
          type
          comp_id
          program
          descriptor
        }}
      }}
    }}
  }}
}}'''
url = "https://data.rcsb.org/graphql"
r = requests.post(url, json={'query': query})
for i in r.json()['data']['entry']['nonpolymer_entities']:
    if i['nonpolymer_comp']['pdbx_chem_comp_descriptor'][0]['type'] == 'SMILES':
        ls_smi.append(Chem.CanonSmiles(i['nonpolymer_comp']['pdbx_chem_comp_descriptor'][0]['descriptor']))

Thanks @AustinApple and welcome to the site. A good answer

OK, seems ligands and PDBs are not really a trending topic here.... (or my posts are really poor).
Nevertheless the question is literally eating my brain so I went into a second direction: "the API" already skimmed here Help me in understanding the PDB file.
So inspecting and copying directly from the PDB entry 7RDX webpage the download value: https://models.rcsb.org/v1/7rdx/ligand?auth_seq_id=1005&label_asym_id=M&encoding=sdf&filename=7rdx_M_1N7.sdfI got this bit of code:
import requests

def main():
    risposta = requests.get(f"https://models.rcsb.org/v1/7rdx/ligand?auth_seq_id=1005&label_asym_id=M&encoding=sdf",
                      allow_redirects=True)
    open(filename, 'wb').write(risposta.content)
    print(risposta.status_code)

if __name__ == '__main__':
   filename = '7rdx_M_1N7.sdf'
   main()

That is equivalent to feeding this:
{
  "atom_site": [
    {
      "auth_seq_id": 1005,
      "label_asym_id": "M"
    }
  ]
}

to the 'https://models.rcsb.org/#/General/ligand-post' test webpage Request body.
now need to figure out how to get the auth_seq_id and the label_asym_id  parameters for all the
'non_polymer_entity_ids' of my targed PDB 7RDX from the https://data.rcsb.org/rest/v1/ DATA API.
Went till here:
import requests
import json

def main(entry_id: str):
    risposta = requests.get(f"https://data.rcsb.org/rest/v1/core/entry/{entry_id}/")
    parsed = json.loads(risposta.content.decode())
    print(risposta.status_code)
    print((parsed))
    print(parsed["rcsb_entry_container_identifiers"]["non_polymer_entity_ids"])
    print('\n', type(parsed))

    parsed_lig = parsed["rcsb_entry_container_identifiers"]["non_polymer_entity_ids"]
    print('\n\n',parsed_lig, type(parsed_lig))

if __name__ == '__main__':
   main('7rdx')

output :
 ['7', '8', '9', '10', '11'] <class 'list'>

The list values are the ids of the 'non_polymer_entity_ids' of my targed PDB 7RDX.
Right now I am trying to recover the the auth_seq_id and the label_asym_id  for every instance of the 'non_polymer_entity_ids' inside my mmCIF/PDBx using the DATA API.
Not sure this is the fastest and right way using the RCSB PDB API to accomplish the OP question (download all ligands of the model structure (not the ligands models)) from the PDB database. Please feel free to add more or point me towards the best solution.
EDITED:
It took a while, I am sure this is not the right way to make good use of all the databases, there should be a way to navigate the API schemas better and faster, and to get and manage the intermediate values needed to get the two keys for the download from ModelServer "Post /v1/{id}/ligand Coordinates of the first group satisfying the given criteria", anyway:
import requests
import json

parsed_lig_dict = {}

def main(entry_id: str):
    risposta = requests.get(f"https://data.rcsb.org/rest/v1/core/entry/{entry_id}/")
    parsed = json.loads(risposta.content.decode())
    # print(risposta.status_code)
    # print((parsed))
    # print(parsed["rcsb_entry_container_identifiers"]["non_polymer_entity_ids"])
    # print('\n', type(parsed))

    parsed_lig = parsed["rcsb_entry_container_identifiers"]["non_polymer_entity_ids"]
    # print('\n\n',parsed_lig, type(parsed_lig))

    # for i in parsed:
    #     print('\n', i  , parsed[i])

    for lig in parsed_lig:
        risposta = requests.get(f"https://data.rcsb.org/rest/v1/core/nonpolymer_entity/{entry_id}/{lig}")
        parsed = json.loads(risposta.content.decode())
        # print('\n',risposta.status_code, risposta.content)

        # print(parsed["rcsb_nonpolymer_entity_container_identifiers"]["asym_ids"],
              # parsed["pdbx_entity_nonpoly"]["comp_id"])

        parsed_lig_dict[lig] = [parsed["pdbx_entity_nonpoly"]["comp_id"] ,
                                dict.fromkeys(parsed["rcsb_nonpolymer_entity_container_identifiers"]["asym_ids"])]

    for lig in parsed_lig_dict:
        for chain in parsed_lig_dict[lig][1]:
            risposta = requests.get(f"https://data.rcsb.org/rest/v1/core/nonpolymer_entity_instance/{entry_id}/{chain}")
            # print(risposta.status_code)

            parsed = json.loads(risposta.content.decode())
            # print('\n',risposta.status_code, risposta.content)
            # print(chain , '--> ',parsed["rcsb_nonpolymer_entity_instance_container_identifiers"]["auth_seq_id"])

            parsed_lig_dict[lig][1][chain] = parsed["rcsb_nonpolymer_entity_instance_container_identifiers"]["auth_seq_id"]

    # print(parsed_lig_dict)

    # for i in parsed_lig_dict:
    #     print(parsed_lig_dict[i])

    # for i in parsed_lig_dict:
    #     for y in parsed_lig_dict[i][1]:
    #         print(i,parsed_lig_dict[i][0], y , parsed_lig_dict[i][1][y])

    cnt = 1
    for i in parsed_lig_dict:
        for y in parsed_lig_dict[i][1]:
            seq_id = parsed_lig_dict[i][1][y]
            chain = y
            comp_id = parsed_lig_dict[i][0]
            print(seq_id , y)

            risposta = requests.get(f"https://models.rcsb.org/v1/7rdx/ligand?auth_seq_id={seq_id}&label_asym_id={chain}&encoding=sdf" ,
                              allow_redirects=True)
            print(risposta.status_code)

            filename = f"lig_{comp_id}_{chain}_{seq_id}_{cnt}.sdf"
            open(filename, 'wb').write(risposta.content)
            cnt += 1

if __name__ == '__main__':
   main('7rdx')

This approach is very slow on my machine, the other two answers that parse the downloaded full mmCIF file via Biopython Structure Object are faster even adding the translation .sdf files to SMILES via Openbabel. I think this is because I translate the DATA API JSON formatted schema into a dictionary to extract the data I need, there should be another way to make use of that data. This way doesn't need installed Biopython though.

Ok, cannot try the API approach, no resources during vacation. Got the download and parse/split approach, added conversion to SMILES by Openbabel. I believe it is exactly the same question of : https://stackoverflow.com/questions/61390035/how-to-save-each-ligand-from-a-pdb-file-separately-with-bio-pdb from where I stole the juicy bit of the code:
## stealing from https://stackoverflow.com/questions/61390035/how-to-save-each-ligand-from-a-pdb-file-separately-with-bio-pdb

from Bio import __version__

print('BIOPYTHON Version : ' , __version__)

from Bio.PDB import MMCIFParser, PDBIO, Select, PDBList

import subprocess

filename = '7rdx.cif'

pdblist = PDBList()

pdblist.retrieve_pdb_file(filename.upper().split('.')[0] , pdir = '.')

pars = MMCIFParser(QUIET = True)

struct = pars.get_structure('7rdx',  filename)

print(struct, type(struct))

def is_het(residue):
    res = residue.id[0]

    if res != " " and res != "W":

        return True #1

    else:

        return False #0

class ResidueSelect(Select):
    def __init__(self, chain, residue):
        self.chain = chain
        self.residue = residue

    def accept_chain(self, chain):

        if chain.id == self.chain.id:

            return True #1

        else:

            return False #0

    def accept_residue(self, residue):
        """ Recognition of heteroatoms - Remove water molecules """

        if residue == self.residue and is_het(residue):

            return True #1

        else:

            return False #0

for model in struct :

    dizio = model.child_dict

    print('\n'+ str(dizio))

for model in struct:

    i = 1

    print(model, model.id)

    for chain in model:

        print(chain, chain.id)

        for res in chain:

            # print('\n'+'res : ', res, 'res.resname : ', res.resname)

            if res.id[0] != " " :

                print('\n'+'res : ', res, res.id, res.resname)

                io = PDBIO()

                io.set_structure(struct)

                io.save(f"lig_{chain.id}_{res.resname}_{i}.pdb", ResidueSelect(chain, res))

                try:
                    subprocess.run(["openbabel.obabel", "-ipdb", f"lig_{chain.id}_{res.resname}_{i}.pdb", "-osmi" ,"-O",

                                    f"lig_{chain.id}_{res.resname}_{i}.smi"])

                except Exception as excp:

                    print('\n ', excp, '\n')

                i += 1

Problem is Biopython parses mmCIFs files; that are now the official standard for RCSB PDB Submission of PDBx/mmCIF format files for crystallographic depositions to the PDB will be mandatory from July 1st 2019 onward. PDB format files will no longer be accepted for deposition of structures solved by MX techniques. A nice review about the long history of the 'new' format :[ PDBx/mmCIF Ecosystem: Foundational Semantic Tools for Structural Biology.
Unfortunately the Bio.PDB.MMCIFParser module return a a structure that is:
<Structure id=7rdx> <class 'Bio.PDB.Structure.Structure'> so all the ligands are related to PDB chains classification and not PDBx entity.
I could be wrong here, I am a novice in both Biopython and Structural Bioinformatics but this is what I've noticed, please feel free to elaborate on this and add a solution using the RCSB APIes ... before curiosity kills the cat
PS:
NEEDs to be tested with sample containing water

adding a different way to solve the question , still without using RCSB PDB database, but using Biopython.PDB.MMCIF2Dict module to start to get a feeling on how info are stored in .cif format files. Ligands are selected amongs the ones in the structure for their Mw.
Ligands are then saved using the same approach of first answer or a new one that takes advantage of Biopython.PDB.Dice module to extract and save portion of PDB. Again Openbabel is used to translate PDBs to SMILES.
Code :
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Jan  5 11:43:07 2023

@author: bob

https://bioinformatics.stackexchange.com/questions/20301/how-to-download-ligands-for-pdb-structure?noredirect=1#comment29972_20301

NEEDs to be tested with sample containing water

"""

## stealing from https://stackoverflow.com/questions/61390035/how-to-save-each-ligand-from-a-pdb-file-separately-with-bio-pdb

from Bio import __version__

print('BIOPYTHON Version : ' , __version__)

from Bio.PDB import (
                    MMCIFParser,
                    PDBIO,
                    Select,
                    PDBList,
                    MMCIF2Dict,
                    Dice
                    )

import subprocess

filename = '7rdx.cif'

# filename = '7RDX.cif'  # this will  work , but biopython saves in lower case

# pdblist = PDBList()

# pdblist.retrieve_pdb_file(filename.split('.')[0] , pdir = '.')

pars = MMCIFParser(QUIET = True)

struct = pars.get_structure('7rdx',  filename)

# print(struct, type(struct))

mmcifdict = MMCIF2Dict.MMCIF2Dict(filename)

# print(mmcifdict)

# for i in mmcifdict:

#     print('\n', i ,'   ---> ', mmcifdict[i])

print('##############################################################')

# print(mmcifdict['_entity'])

print(mmcifdict['_pdbx_nonpoly_scheme.entity_id'])
print(mmcifdict['_pdbx_nonpoly_scheme.asym_id'])
print(mmcifdict['_pdbx_nonpoly_scheme.entity_id'])
print(mmcifdict['_pdbx_nonpoly_scheme.mon_id'])

ligands_map_mmcif = list(zip(mmcifdict['_pdbx_nonpoly_scheme.entity_id'],

                             mmcifdict['_pdbx_nonpoly_scheme.asym_id'],

                             mmcifdict['_pdbx_nonpoly_scheme.entity_id'],

                             mmcifdict['_pdbx_nonpoly_scheme.mon_id']))

print('#########################################')

print(ligands_map_mmcif)

print(mmcifdict['_pdbx_nonpoly_scheme.entity_id'])
print(mmcifdict['_pdbx_nonpoly_scheme.pdb_seq_num'])
print(mmcifdict['_pdbx_nonpoly_scheme.pdb_strand_id'])
print(mmcifdict['_pdbx_nonpoly_scheme.pdb_mon_id'])

ligands_map_pdb = list(zip(mmcifdict['_pdbx_nonpoly_scheme.entity_id'],

                           mmcifdict['_pdbx_nonpoly_scheme.pdb_seq_num'],

                           mmcifdict['_pdbx_nonpoly_scheme.pdb_strand_id'],

                           mmcifdict['_pdbx_nonpoly_scheme.pdb_mon_id']))

print('#########################################')

print(ligands_map_pdb)

print('WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWw')

print(mmcifdict['_entity.id'])
print(mmcifdict['_entity.formula_weight'])

entity_Mw = dict(zip(mmcifdict['_entity.id'],

                     mmcifdict['_entity.formula_weight']))

print(entity_Mw)

print('-------------------------------------------------------------')

ligands_map_pdb_red = [i for i in ligands_map_pdb if float(entity_Mw[i[0]]) >= 200 ]

print(ligands_map_pdb_red)

#########################################################   first approach

# class ResidueSelect(Select):
#     def __init__(self, chain, residue):
#         self.chain = chain
#         self.residue = residue

#     def accept_chain(self, chain):

#       if chain.id == self.chain.id:

#           return True #1

#       else:

#           return False #0

#     def accept_residue(self, residue):
#         """ Recognition of heteroatoms - Remove water molecules """

#         if residue == self.residue :

#             return True #1

#         else:

#             return False #0

# for model in struct:

#     print(model)

#     for chain in model:

#         if chain.id == 'A':

#             for res in chain:

#                 print(res, res.id, res.get_resname(), type(res))

# i = 1

# for res in ligands_map_pdb_red :

#     print('\n'+'res : ', res[3], res[0], res[1], type(res[1]), res[2])

#     print(struct[0][res[2]] , type(struct[0][res[2]]) )

#     print(struct[0][res[2]][('H_'+res[3], int(res[1]), ' ')] , type(struct[0][res[2]][('H_'+res[3], int(res[1]), ' ')]) )

#     io = PDBIO()

#     io.set_structure(struct)

#     io.save(f"lig_{res[3]}_{res[2]}_{str(i).zfill(5)}.pdb", ResidueSelect(struct[0][res[2]],

#                                                     struct[0][res[2]][('H_'+res[3], int(res[1]), ' ')]))

#     try:
#         subprocess.run(["openbabel.obabel", "-ipdb", f"lig_{res[3]}_{res[2]}_{str(i).zfill(5)}.pdb", "-osmi" ,"-O",

#         f"lig_{res[3]}_{res[2]}_{str(i).zfill(5)}.smi"])

#     except Exception as excp:

#         print('\n ', excp, '\n')

#     i += 1

########################################################   second approach

import warnings

from Bio import BiopythonWarning

def accept_residue_mod(self, residue):
        """Verify if a residue sequence is between the start and end sequence."""
        # residue - between start and end
        hetatm_flag, resseq, icode = residue.get_id()
        # if hetatm_flag != " ":
        #     # skip HETATMS
        #     return 0
        if icode != " ":
            warnings.warn(
                f"WARNING: Icode {icode} at position {resseq}", BiopythonWarning
            )
        if self.start <= resseq <= self.end:
            return 1
        return 0

Dice.ChainSelector.accept_residue = accept_residue_mod

i = 1

for res in ligands_map_pdb_red :

    print(struct, res[2], int(res[1]), int(res[1]))

    Dice.extract(struct, res[2], int(res[1]), int(res[1]), f"lig_{res[3]}_{res[2]}_{str(i).zfill(5)}.pdb")

    try:
        subprocess.run(["openbabel.obabel", "-ipdb", f"lig_{res[3]}_{res[2]}_{str(i).zfill(5)}.pdb", "-osmi" ,"-O",

        f"lig_{res[3]}_{res[2]}_{str(i).zfill(5)}.smi"])

    except Exception as excp:

        print('\n ', excp, '\n')

    i += 1

```


text:
Why are conversion output SMARTS strings are not the same for different file types of the same structure?
I've just used Open Babel to convert .sdf and .mol2 files to SMARTS, and I'm a bit surprised the output can be different for the same structure. It's confusing me that one structure, one piece of information can be different depending on the file format in which it was saved.
I converted a ligand 1N7 in two available formats (.sdf and .mol2)  of the 7RDX structure from RSCB PDB for chain N, and also for chain U.
After the conversion, the SMARTS strings for the ligand in both chains are not the same for both input files types (i.e. the difference between the same ligand, from the same chain, just in two different formats available to download from RSCB PDB):
example no. 1:
[additional @H added in places, some square brackets added]
$ cat ./7rdx_N_1N7_from_mol2.smi

C1[C@@]2(C)[C@@H](C[C@@H](C1)O)[C][C]([C]1[C@H]3[C@]([C@H](C[C]21)O)([C@H](CC3)[C]([C])[C][C])C)O   1N7

$ cat ./7rdx_N_1N7_from_sdf.smi

C1[C@@]2(C)[C@@H](C[C@@H](C1)O)C[C@H]([C@H]1[C@H]3[C@]([C@H](C[C@H]21)O)([C@H](CC3)[C@H](C)CC)C)O   1N7

example no. 2:
[some square brackets added]
$ cat ./7rdx_U_1N7_from_mol2.smi

C1[C@@]2(C)[C@@H](C[C@@H](C1)O)C[C@H]([C@H]1[C@H]3[C@]([C@H](C[C@H]21)O)([C@H](CC3)[C@H](C)CCC(=O)NC[C][C][N+]([C])([C])[C])C)O 1N7

$ cat ./7rdx_U_1N7_from_sdf.smi

C1[C@@]2(C)[C@@H](C[C@@H](C1)O)C[C@H]([C@H]1[C@H]3[C@]([C@H](C[C@H]21)O)([C@H](CC3)[C@H](C)CCC(=O)NCCC[N+](C)(C)C)C)O   1N7

In both examples, the structures look quite similar to each other, but the difference is a little worrisome, because there is no chance to check every output by hand.
Is it a normal behaviour during the conversion?
Is it related to a different notation in the input files?
Which input format is better?

Mmmhhh be sure that each of the 3 modeled ligands correspond to the full one, could be that the density is poor and pdb shows only what is really ((??) there and not what they soaked the crystal with, actually 1N7 is CHAPSO I believe a detergent

Prior to grid freezing, 3-([3-cholamidopropyl]dimethylammonio)-2-hydroxy-1-propanesulfonate (CHAPSO, Anatrace) was added to the sample (8 mM final), resulting in a final complex concentration of 8 μM. The final buffer condition for the cryo-EM sample was 20 mM HEPES pH 8.0, 150 mM K-Acetate,10 mM MgCl2, 2 mM DTT, 1 mM ADP, 1 mM AlF3, 8 mM CHAPSO.

to me your 2 strings seem different, no N in first molecule.
Here pictures, I used a PDB version of the original mmCIF file
but dindt use the RCSB PDB available one, imported .cif into
Biopython and parsed the resulting Struct object [so there could be
mistakes].
so chain N --> A and chain U --> E

as per in my comment:

Mmmhhh be sure that each of the 3 modeled ligands correspond to the full one, could be that the density is poor and pdb shows only what is really ((??) there and not what they soaked the crystal with, actually 1N7 is CHAPSO I believe a detergent

... I could be wrong though, this mmCIF format is driving me crazy