text:
How to analyse qualitatively the penetration ability of particles in spheroids using fluorescent z-stacks?
0
I have to establish the penetration profile of particles in tumour spheroids.
For this I have spheroids composed of cells which were exposed to particles, both of which are fluorescently labeled. The spheroids were then imaged using a confocal microscope by taking z-stacks of all spheroids. For analysis, I determined the intensity profile of each focal plane using an ImageJ plugin, namely "distribution_shapes".
What should be done next?
I currently did not find any step-by-step guide for this analysis for manual assessment and any type of help would be greatly appreciated.

This looks like an interesting image analysis question, which is great because we don't get as many of those here. Could you be more specific about what you mean by "penetration profile"? Can you point to any publications that use a similar analysis? I find one possibly related paper here but they don't go into any detail unfortunately.

I am not sure if this is exactly what you're looking for, but I was able to find this ImageJ script that estimates a penetration profile in tumor spheroids.
I am not an ImageJ person but I suspect that with some minimal changes this could be used to answer your question.

This is actually what I was using to determine the intensity profile of my particles. However, I am not sure what to do next, because I have the intensity profile of a spheroid over z-axis and don't know how to put the information together or to visually represent the quantification of the penetration profile.

@Timi thanks for clarification. Can you post what your data look like (just a snippet is probably enough)? That might be something we can help with.

I added a screenshot of how to profile intensities look across a z-stack (1st image) and how the average profile intensity looks like in the last 5 focal planes in the z-stack. This is what we are aiming to get.
The papers that I found that do something similar are: https://www.sciencedirect.com/science/article/pii/S0168365921006672
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7824314/
, but they also don't go into details, and the script in Matlab is not working for me.

@Timi Can you show the original raw data (text, not plots) output of the imagej protocol, rather than plotted summaries? As in, what is the data that you used to make those plots. It seems that you want to make better or more informative plots similar to some that you are seeing in those papers (which figures? maybe include them in the question). Just a subset of the data is fine, to show what the format is, or even faked data that is in the same format if you don't feel comfortable posting raw data. I don't see the matlab script you are talking about, a direct link might help.

I don't know how to attach a table to this post. Would it be possible to send you an Excel sheet in an email or private message?

@Timi What I'd suggest doing is using table markdown or similar to embed a table in your question. There are web converters for converting spreadsheet to markdown. Alternately you could use a google drive link, pastebin, dropbox, box, etc. In each case I'd suggest embedding the link directly in your question so that anyone who needs to can look at it.

The next step depends heavily on what your question is. Some general ideas, not specific to tumor particle penetration.

Your current visualization seems OK to me, if you want to make some features more obvious you can try density estimation.

If you have a generative model (for example you believe this is a mixture of a "shallow" and "deep" Gaussian distribution, or you have reason to believe the particles follow an exponential distribution), you can fit your data to this distribution and estimate its parameters.

If you want to compare conditions without too strong assumptions on the distributions, you could compute the cumulative distribution of the particles. There are statistical approaches to compare such curves.

To compare conditions, you might even be able to come up with a simpler metric, for example "total fluorescence in the core of the spheroid" normalized by "total fluorescence in its periphery".

Ultimately it all depends on your research question (are there treatments to compare to each others, biological replicates whose consistency you want to check, known structure within the tumor spheroid, a time course, ...). My suggestions above would be easier to implement in R/Matlab/Python than pure ImageJ.

I want to compute the cumulative distribution of the particles. But I don't know how to do this with a z-stack as the distribution across z-direction will not be the same as in x and y-directions. I also don't have any knowledge of programming, that's why I was wondering if anyone knows a guide or course that can help in this regard. I did find a matlab script, but I get the same errors and I asked several people to help me with that and we didn't manage to solve it.


text:
SNP How to detect regions of erratic pairwise alignment?
I have clusters of DNA sequences that contain some mutations. All the sequences in a cluster should contain the same mutations.
As the mutations aren't known, those clusters of sequences were pairwise-aligned against the non-mutated reference with Minimap2. Then an interesting fact cropped-up: a few mutations close to each other will make the alignment go awry in that zone.

Zone with a few SNPs:
ref  GTATGCGGATACCG

sq1  GTATCGTGACCG
sq2  GTATCTCGTGACCG
sq3  GTGTCTATTCGTGACCG
sq4  GTATCTCGTGACCG

Independent pairwise alignments against the reference (I removed the insertions from the output):
ref  GTATGCGGATACCG

sq1  GTAT-CGG--ACCG
sq2  GTATCTCGTGACCG
sq3  GTATTCGG--ACCG
sq4  GTATCTCGTGACCG

So I made a MSA of this small cluster (adding the reference in it) with Muscle5; the difference is noticeable, as we can clearly "see" the mutations now:

MSA alignment that includes the reference:
ref  GTATG---CGGATACCG

sq1  GTAT-----CGTGACCG
sq2  GTATC---TCGTGACCG
sq3  GTGTCTATTCGTGACCG
sq4  GTATC---TCGTGACCG

Same MSA, but trimmed:
ref  GTATGCGGATACCG

sq1  GTAT--CGTGACCG
sq2  GTATCTCGTGACCG
sq3  GTGTCTCGTGACCG
sq4  GTATCTCGTGACCG

I would say calculating the MSA seems to be the way to go but the clusters sizes range from 2 to 1000 sequences, the sequences are 3000nt long, and there's 2,000,000 sequences in total; we cannot afford to do the MSA of each cluster on the full size of the sequences.

On the other hand, we could probably manage to do a MSA of a few small zones in the sequences.
So, given a cluster and the pairwise alignments of its sequences against the reference, how would you determine the zones that need to be re-aligned?
note: I'm using Muscle5 because there's a Python binding of it and the results seem good. I also tried FAMSA, which is very fast, but it didn't yield good alignments.

Are sequences 1-4 in your 3 example supposed to be the same sequence? It is not really clear to me what you are trying to do / which problem you are trying to solve. You seem to be dealing mainly with indels rather than SNPs? Its also not clear why the alignment is awry, it is probably the best alignment with the given penalty matrix..

@Niklas Thank you for the comment. Yes, they are the same sequences (I updated the examples to make it clearer). By "awry" I mean that I cannot  make a consensus using the pairwise alignments, as they are subjectively "wrong" (mostly because of the indels)

OK. Generally, all alignments are equally "correct" or "wrong", since they minimize penalty score given the penalty matrix, which you can edit as you prefer. You would expect a msa to yield different results from pairwise alignments since it is a different problem.. regarding your actual question: I assume you could identify "bad" regions from a long stretch by the presence of many gaps in that region. But, keep in mind that when you re-align short sequences you may get different results than for the long sequences.

Maybe something like pileup is what you are looking for? https://samtools.github.io/bcftools/howtos/variant-calling.html

@Niklas Thank you, bfctool will definitely be useful

source: Stack Exchange

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

sort: votes


text:
Generating CIGAR strings for MSA with reference
I have a list containing the results of a MSA:
msa = [
    '-G-ATG-CGGATACCG-', # reference sequence
    'CGTATCTCG--TACCGA',
    '-GTATCTCGT-GACCG-',
    'CG-AT--CGG-GACCGA',
    '---AT--CG--GACCGA',
    '-GTAT---GC-GACCG-',
    'CG-AT--CGG-GACCGA',
    'CGTAT--CG--GA CG-',
]

I would like to generate the CIGAR strings of those alignments (against the reference sequence); the expected result would be something like:
cigars = [
    '13M',
    '1I1M1I3M1I2M2D5M1I',
    '1M1I3M1I3M1D5M',
    '1I3M1D3M1D5M1I',
    '1D2M1D2M2D5M1I',
    '1M1I2M2D2M1D5M',
    '1I3M1D3M1D5M1I',
    '1I1M1I2M1D2M2D5M',
]

I don't know what would be a good approach for this. I managed to get the "aligned pairs" of each alignment but I don't think that'll help for determining the cigar strings:
# mapping = position in alignment => position in reference
mapping = {y: x for x,y in enumerate(y for y,nuc in enumerate(msa[0]) if nuc != '-')}

for aln in msa:
    aligned_pairs = [(y, mapping[z]) for y,z in enumerate(z for z,nuc in enumerate(aln) if nuc != '-') if z in mapping]
    print(aligned_pairs)

[(0, 0), (1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (6, 6), (7, 7), (8, 8), (9, 9), (10, 10), (11, 11), (12, 12)]
[(1, 0), (3, 1), (4, 2), (5, 3), (7, 4), (8, 5), (9, 8), (10, 9), (11, 10), (12, 11), (13, 12)]
[(0, 0), (2, 1), (3, 2), (4, 3), (6, 4), (7, 5), (8, 6), (9, 8), (10, 9), (11, 10), (12, 11), (13, 12)]
[(1, 0), (2, 1), (3, 2), (4, 4), (5, 5), (6, 6), (7, 8), (8, 9), (9, 10), (10, 11), (11, 12)]
[(0, 1), (1, 2), (2, 4), (3, 5), (4, 8), (5, 9), (6, 10), (7, 11), (8, 12)]
[(0, 0), (2, 1), (3, 2), (4, 5), (5, 6), (6, 8), (7, 9), (8, 10), (9, 11), (10, 12)]
[(1, 0), (2, 1), (3, 2), (4, 4), (5, 5), (6, 6), (7, 8), (8, 9), (9, 10), (10, 11), (11, 12)]
[(1, 0), (3, 1), (4, 2), (5, 4), (6, 5), (7, 8), (8, 9), (9, 10), (10, 11), (11, 12)]

How can I generate the CIGAR strings efficiently?

If you loop through the characters of two aligned sequences in parallel then you can easily determine the operations that compose the CIGAR. Then you just have to "collapse" the generated list of operations:
import itertools

ref = msa[0]

for aln in msa:
    cig = []

    for i,nuc in enumerate(aln):
        if nuc != '-':
            if ref[i] != '-':
                cig.append('M')
            else:
                cig.append('I')
        elif ref[i] != '-':
                cig.append('D')

    cig = [(len(list(it)), op) for op,it in itertools.groupby(cig)]

    print( ''.join(str(n) + op for n,op in cig) )

13M
1I1M1I3M1I2M2D5M1I
1M1I3M1I3M1D5M
1I3M1D3M1D5M1I
1D2M1D2M2D5M1I
1M1I2M2D2M1D5M
1I3M1D3M1D5M1I
1I1M1I2M1D2M2D5M

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


text:
Changing transparency of markers in Chimera X
I am using the software Chimera X to visualize the 3d composition of molecules. I am very new to Chimera and only use it for visualization. The files I'm using are cmm files (xml-formatted) basically looking like below:
...
<marker id="14" x="0.018" y="-0.180" z="0.215" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="15" x="-0.032" y="-0.191" z="0.241" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="16" x="-0.033" y="-0.117" z="0.246" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="17" x="0.020" y="-0.138" z="0.244" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="18" x="0.047" y="-0.083" z="0.280" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="19" x="0.073" y="-0.061" z="0.240" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="20" x="0.010" y="-0.068" z="0.271" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
...
<link id1="0" id2="1" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
<link id1="1" id2="2" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
<link id1="2" id2="3" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
<link id1="3" id2="4" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
<link id1="4" id2="5" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
...

i.e. just a simple polymer consisting of alot of markers (spheres) and links between them (cylinders), and with colors for each one specified in rgb (see here). This is all well and good, but now I would like to set transparency for a subset of the markers and associated links. Preferably I would set this directly in the attributes in the cmm files, but I cannot find how to do this. The manual does mention that colours in general can be set with rgba where a means alpha level, but when trying to set a like seen in the xml block above, nothing happens (I tried values close to both 0 and 1; no difference)
The manual documents a function called 'transparency' which I tried to use with for example transparency 50,a #0:10-1200
to make atoms 10-1200 50% transparent. But I just get the error: "Missing or invalid "percent" argument: Expected a number". If I only type transparency 50 everything does turn transparent, however, I wouldn't really call it transparent since I cannot really see through them to see the other molecules on the other side, they just turn bleak.
Questions:

(Would someone explain how I could simply change the actual transparency of a large set of markers and links) (preferably by loading some file with specified marker id:s)? Also, is it possible to assign this directly in the cmm files, and if so, how?

UPDATE: The command transparency #1:2492-end 90 'All' worked to set the transparency of markers e.g. 2492-end. What I would still like to know is if it is possible to set this directly in the cmm-files.

RGBA color values are an extension of RGB color values with an alpha channel - which specifies the opacity for a color. An RGBA color value is specified with: rgba(red, green, blue, alpha). The alpha parameter is a number between 0.0 (fully transparent) and 1.0 (fully opaque). , what does the  a="0"  do in your xml formatted markers and links ?

So, the 'a="0"' was a naive attempt by me to set the alpha-level directly in the xml file, which as I noted above did not work. I have since figured out how to set the transparency and will update my question.

you should try to report the bug a = 0 not working to chimerax developers. To report a problem with ChimeraX:
If possible, please use Help... Report a Bug in the ChimeraX menu. Otherwise (for example, if ChimeraX will not start), send e-mail to chimerax-bugs@cgl.ucsf.edu, making sure to include any information needed to reproduce the problem: program version, specific data, actions/steps taken, etc. The information will be publicly accessible, so do not include data or descriptions that you wish to remain private.

from https://www.cgl.ucsf.edu/chimerax/docs/contact.html

I will look in to that, however I am not certain it is a bug: Even though it seems very reasonable to be able to set it that way, that was a guess on my part. The general support of rgba colors is however otherwise documented in the manual. Regarding your suggested edit: The markers are literally referred to as 'atoms' in Chimera, and it is those which I wish to make transparent. (I am however changing it to 'marker' as a compromise =).)

A marker is a pseudoatom, typically used to mark the position of some feature within volume data; markers are implemented as atoms, so actions that apply to atoms in ChimeraX also apply to markers: they can be specified in commands, selected, etc. https://www.cgl.ucsf.edu/chimerax/docs/user/markers.html

dont know why I didnt like ChimeraX at the time , ages ago,  but haved worked only with just plain Chimera, had to dig an old VMwareVM to take these few pics

users mailing list is https://mail.cgl.ucsf.edu/mailman/mailman3/lists/chimera-users.cgl.ucsf.edu/ .. dont get why  a='50' or a='0.5' dont work with chimera too

Ok, this is not an answer, but to better visualize the problem I recreated the question/answer on Chimera that has a different syntax for the transparency  command see : https://www.cgl.ucsf.edu/chimera/docs/UsersGuide/midas/transparency.html and https://www.cgl.ucsf.edu/chimerax/docs/user/commands/transparency.html .
Input file is marker.cmm , to be opened in FILE MENU of VOLUME TRACER from VOLUME MENU :
<marker_set name="marker set 1">
<marker id="14" x="0.018" y="-0.180" z="0.215" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="15" x="-0.032" y="-0.191" z="0.241" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="16" x="-0.033" y="-0.117" z="0.246" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="17" x="0.020" y="-0.138" z="0.244" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="18" x="0.047" y="-0.083" z="0.280" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="19" x="0.073" y="-0.061" z="0.240" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="20" x="0.010" y="-0.068" z="0.271" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<link id1="14" id2="15" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
<link id1="15" id2="16" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
<link id1="17" id2="18" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
<link id1="18" id2="19" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
<link id1="19" id2="20" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
</marker_set>

Image before command :

Image after transparency  50,a,b #0:16-end
(#0 being 0 model as the only loaded model on my Chimera instance) :

Better transparency visualizations :

Dont get why a='50' or a='0.5' in marker,link definition dont work with chimera too

Together with my edited question I think this all suffice as "good enough" for now at least. I am satisfied for the moment just using something like transparency #1:2492-end 90 'All', but I guess when I get more complicated regions it would be good to be able to load the marker spec, i.e. e.g. #1:2492-2500, 3000-3500, ... from a file if that is possible. But I guess one can run command-files right? and I suppose there's no limit on the length on the spec, so perhaps I could just create one very long command in a file. (Many thanks for your efforts!!! unfortunately I need more rep to upvote)


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

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

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

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

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


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

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

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

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


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

CO1 mithocondrial DNA
16S rRNA
18S rRNA

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

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

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

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

Name
SeqID

Oncorhynchus mykiss
EU186789.1

Xenopus laevis
AB278691.1

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

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

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

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

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

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

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

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

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

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

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

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

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

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


text:
What does the different sequences represent?
I am using this package nsdpy to download genome sequences from NCBI nucleotide database.
Specifically I am interested in the whole mitochondrial genome of different species, here I will use a subset of species just as an example.
My Query is:
"{mitochondrion[Title]" -i -t -L subset_species_list.txt
Where :

-i Adds the taxonomic information to the information lines of the sequences written in the output files (not relevant for this example)
-t option produces a tsv file with the following fields: Taxon name, SequenceID, TaxID, Taxonomic lineage, Sequence length, Sequence.
subset_species_list.txt is a list of species, below you can see its content

Xenopus laevis
Gasterosteus aculeatus
Chironomus riparius
Rasbora heteromorpha
Planktolyngbya limnetica

The output file is quite long and includes the actual sequences. Here I will post the same file where I dropped the column for the sequence and the lineage to keep it compact

My aim is to obtain phylogenetic distances between species.
My questions are:

What do the sequences for the same species represent?
Which sequences are best suited for an alignment and how do I know that?
Should I also look into the orientation of the sequences or are they all on the forward strand?
Is it better if I filter for Coding sequences (CDS) instead of whole genomes

Also if you noticed any rookie mistake in my process feel free to point that out. It is the first time for me working with DNA sequences and I really want to learn to correct way to do that.

Welcome to the site. This question is too broad could you limit the number of questions please?

Welcome! I have some quick basic answers for the questions you're asking so hopefully this helps:

Sequences for the same species represent variation/diversity/evolutionary changes. Think of humans; how do different people's DNA sequences represent the sequence of a human? This is especially important for phylogenetic studies.
For the most part, all sequences are best suited for alignment. The main thing to think about is the amount of data you have (if it's a lot of GB the alignment will take a lot of computing power), the type of sequnece data you have (simple rule of thumb is to not mix short read with long read sequences), and the tool you're using for alignment.
Always a good idea to think about orientation but most tools will account for orientation when aligning sequences.
This depends on what you're looking for but more often than not using the whole genomes is best as it lets you look at the whole picture instead of just a small portion but that's definitely a lot of computing power. If you're just looking at the mitochondria just focus on that region

Hope this helps!


text:
How do I select a subset of genes for functional enrichment(GO/KEGG) analysis from WGCNA results?
I am working on using WGCNA for a bulk RNA-seq experiment. I have three experimental conditions (control, treatment 1, and treatment 2) and have a total sample size of 18 (6 per group). Using the online WGCNA tutorials, I could analyze my datasets after variance stabilizing transformation (vst). I got the summary of network analysis results as a CSV with Gene significance (GS), p.GS, module membership(MM), and associated p-value (p.MM).
I want to extract a set of significant genes from the module with the highest positive correlation to my treatment conditions. As I am just getting acquainted with WGCNA interpretation, which of these parameters in the CSV should I use for filtering further (GS, p.GS, and module membership and associated p-value (p.MM))?
Also, how would I choose the limits of the criteria? Should it be using the scatter plot (MM vs GS)?
As I am pretty new to WGCNA, any help is appreciated. Thanks in advance!

Consider filtering on every parameter: Gene Significance (GS), p-value for Gene Significance (p.GS), Module Membership (MM), and p-value for Module Membership (p.MM).
This doesn't mean you need to necessarily filter on every parameter, but it depends on what exactly you want to do, and you might reasonably choose any/all.

For GS, you may want to ensure you're capturing genes strongly associated with your experimental conditions.

For p.GS, you may want to ensure statistical relevance, (so filter out genes with high p-values). Similar for p.MM.

For MM, you may want it to be high to make sure the genes you select are strongly co-expressed and are relevant members of the module.

Plot MM vs GS to see if you can see an obvious filter threshold.


text:
How to make a POST "vep/:species/region" request that is able to accomodate multiple alternate alleles for a single variant?
There is an R code example on Ensembl's website on how to execute a POST "vep/:species/region" request for variants that only have one alternate allele. The part that holds the variant data in the POST request is structured as follows:
..., body = '{ "variants" : ["21  26960070  rs116645811 G A . . .", "21  26965148  rs1135638 G A . . ." ] }'...

Where, for the first variant, 21 is the chromosome number, 26960070 is the start point coordinate on the chromosome, rs116645811 is the variant identifier, "G" is the reference allele and "A" is the alternate allele.
Alternatively, you can substitute the variant identifiers - rs116645811 and rs1135638, in the example above - for the endpoint coordinate of the variant. I use this option myself and can personally confirm that it works.
As mentioned before, the example above shows a POST request of variants that have only single alternate alleles.
Does anybody, on the other hand, know how (one should adjust the structure in the example above) to execute a POST request that accommodates multiple alternate alleles for a single variant?
If you explain your answer with an example like the one above, please use an example with the endpoint coordinate (and not the variant ID).
Thanks!

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

Thanks for your question. The question has received a "closed" request. The site convention is for the OP to present an example for others - to  help them answer - here you as OP are  requesting the person answering to present the answer. This is can lead to misunderstandings about what is required. There are some very skilled R coders on the site who could probably answer this.


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

Thank you all in advance!

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

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

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

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

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

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


text:
Convert genotyping results to vcf or other plink formats
I used a Genotyping software called T1K to generate KIR genotyping, the output file contains genes and alleles for each sample, how can i convert the results to a vcf file or other formats which can be used in plink?
The file i have for now looks like this:
Gene_name Allele_1 Allele_2
KIR2DL1 KIR2DL1*004 KIR2DL1*042
KIR2DL4 KIR2DL4*005 KIR2DL4*005
KIR3DL1 KIR3DL1*005 KIR3DL1*005
KIR3DP1 KIR3DP1*006 KIR3DP1*010

My thought is to first convert the above file to a vcf file and use plink to generate bedbimfam from it, a standard vcf is like this:
##fileformat=VCFv4.3
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=file:///seq/references/1000GenomesPilot-NCBI36.fasta
##contig=<ID=20,length=62435964,assembly=B36,md5=f126cdf8a6e0c7f379d618ff66beb2da,species="Homo sapiens",taxonomy=x>
##phasing=partial
##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples With Data">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
##INFO=<ID=AA,Number=1,Type=String,Description="Ancestral Allele">
##INFO=<ID=DB,Number=0,Type=Flag,Description="dbSNP membership, build 129">
##INFO=<ID=H2,Number=0,Type=Flag,Description="HapMap2 membership">
##FILTER=<ID=q10,Description="Quality below 10">
##FILTER=<ID=s50,Description="Less than 50% of samples have data">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
##FORMAT=<ID=HQ,Number=2,Type=Integer,Description="Haplotype Quality">
#CHROM POS      ID         REF   ALT    QUAL  FILTER   INFO
FORMAT       NA00001         NA00002          NA00003
20     14370    rs6054257  G     A      29    PASS    NS=3;DP=14;AF=0.5;DB;H2           GT:GQ:DP:HQ  0|0:48:1:51,51  1|0:48:8:51,51   1/1:43:5:.,.
20     17330    .          T     A      3     q10     NS=3;DP=11;AF=0.017               GT:GQ:DP:HQ  0|0:49:3:58,50  0|1:3:5:65,3     0/0:41:3
20     1110696  rs6040355  A     G,T    67    PASS    NS=2;DP=10;AF=0.333,0.667;AA=T;DB GT:GQ:DP:HQ  1|2:21:6:23,27  2|1:2:0:18,2     2/2:35:4
20     1230237  .          T     .      47    PASS    NS=3;DP=13;AA=T                   GT:GQ:DP:HQ  0|0:54:7:56,60  0|0:48:4:51,51   0/0:61:2
20     1234567  microsat1  GTC   G,GTCT 50    PASS    NS=3;DP=9;AA=G                    GT:GQ:DP     0/1:35:4        0/2:17:2         1/1:40:3

Any tools can be used to realize it? For now i only have Gene names and allele names for each sample, How to obtain other genetic informtation a vcf file needed?

Hi! I am not so familiar with this output format. Could you explain what we are looking at in the output sample?

@winni2k Thank you for your reply. I have added more details to my question as you suggested. Could you please take a look at it again?

I am not at all familiar with KIR genotyping, but is there a way of converting things like KIR2DL1*042 into chromosome, position, ref and alt? Is there some database that lists what these KIR notations actually mean?

I would also recommend searching dbsnp for the kir mutations.

your output file does not seem to match what is described in T1K docs. Are there other files output? T1K is supposed to output a VCF even.

@MaximilianPress  Thank you for reply. The output i posted is based on the genotype ouptput, The vcf output only contains noval SNP, i dont know how to obtain those info for other varients. Any methods to solve this? Thanks~

Genotype output according to the docs looks like this: gene_name num_diff_alleles allele_1 abundance_1 quality_1 allele_2 abundance_2 quality_2. That doesn't match what you show. But regardless, if they really don't want to give you non-novel alleles there is likely a reason, I'd suggest contacting the authors on how you might do that. I glanced at the code but was not sure how to get non-novel, though this line might be relevant.

@MaximilianPress Thank you, i really appreciate ur advice.


text:
Get a certain gene sequence from bam/vcf and reference
I need to get a fasta sequence of a certain gene for a certain worm strain that is different from reference. I have a reference genome, BAM for the strain of interest, and coordinates of the gene. I know that vcftools can convert bam to fasta, but I do not see if I can do it only for a certain coordinates (like, for example, chr1:10100500-10200500). I converted the whole bam to fasta and then I tried using bedtools getfasta, but my system ran out of space. Is there a more rational way to do it?

The first approach that comes to mind is to cut the bam file and then convert the already cut bam to fasta:
samtools view -hb chr1:10100500-10200500 > small.bam
samtools fasta small.bam > small.fa

(note that I am not certain about the samtools fasta command, I have never used it and I don't know if it will produce many sequences, one for each haplotype supported by reads; please test this.)

https://www.htslib.org/doc/samtools-fasta.html -- the documentation for samtools fastait should work as @terdon wrote it but just in case you need more information


text:
How long should it take to cluster Uniref?
I have the Uniref90 dataset containing around 76M sequences. I have made modifications to all sequences based on a reduced amino acid alphabet.
I am trying to cluster it again at 90% identity. I am using  mmseq cluster on a 64 core, 500GB RAM machine. It has been more than 12h now.
Starting prefiltering scores calculation (step 1 of 1)
Query db start 1 to 74390789
Target db start 1 to 74390789
[=

How long does it normally take to cluster such a database? I know I can use linclust instead to cluster in linear time, but it might not be what I need in this case.
I am using mmseq cluster as follows:
mmseqs createdb "uniref90-$1.fasta" "uniref90-DB-$1"
mmseqs cluster "uniref90-DB-$1"  "uniref90-DB-$1_clu" tmp/ --cov-mode 2 -c 0.8 --min-seq-id 0.9

Thanks. Hmmm ... when you use htop what do you see when its running?

CPU is used at the maximum and 20% of RAM. Anyway, I have looked into the mmseqs2 in more detail and found that it should take at least a few days.

Thanks and thank you for your question. Thats very clear and answered below.

What does "CPU is used at the maximum" mean? Do you see 100% in the CPU% column? Or do you see 6400%?

6400 %CPU, 12.8 %MEM

When I want to estimate the time taken to do something big, I subset my input data to something that takes a few minutes or less, and calculate forward from there. If there's alignment involved, it's usually a good idea to work it out for multiple dataset sizes, as there could be a linear relationship between input size and time taken, square relationship, or worse.
Have you looked at how long it takes for, say, 1000 sequences, 2000 sequences, and 5000 sequences?

The mmseqs2 code should be specifically compiled and optimised for your system. mmseqs2 optimisation is automated and might be available via Spack* (I'm not sure its actually done for mmseqs2), it is definitely available for EasyBuild. I wouldn't advise attempting a manual optimisation unless you've prior experience.
How much time this will save is more difficult to assess but 50% faster is a conservative assessment. This approach is good practice for a job like this size and its good for the associated report, its standard HPC practice. If you're doing cloud computation (i.e. $$$), HPC optimisation is essential. AWS make serious money from suboptimal HPC practice (in fairness to them they also do a lot of development work to minimise costs).
Finally, if you're an expert in optimisation you can save serious computational time, but this is a major investment of your personal development time and technique. The rationale of how and why is complicated.
*, its very likely available on AWS.

You might simply consider its been running for a few days already so HPC technique isn't needed for this specific job. However, you'll know for the future.

Is the tool multithreaded? The OP mentioned "CPU is used at the maximum" in htop which makes me think that it might be 100% of a single core.

@terdon. Yes it is. Whether it's default - dunno - I was thinking that too. Which is one reason why I requested the htop, but OPs response strongly suggests the parallelisation is default. It's easy to code parallelisation as default (at least in Python), normally its Total CPU - 1. What is a bit unusual I'll admit ... one CPU should be free (standard coding practice) but the OP didn't report that. However, we don't know how the original coder did this. To be honest mmseqs coder is hot, so they'd know the convention. If mmseqs2 is intended for HPC Total CPU - 1 is dropped


text:
Generating force-fields for ANY residue
Is this possible to generate force fields for any residue in my .pdb file? I am trying to generate conformations for the proteins with arbitrary residues binded to it from .pdb files using OpenMM, but there are always errors like 'missing residue VAL, missing residue HEM, etc'. I tried to use pdbfixer from openmm library, skip\replace missing\nonstandard residues using GROMACS and generate .top .gro, but even after that it shows me errors or generates nonsence if I replace residues. I read some articles and looks like generating force fields is a complicated task. But still wondering is there any open-source framework\library\neural network that can generate such force fields?
For example I follow this guide
http://docs.openmm.org/latest/userguide/application/02_running_sims.html#a-first-example
Whenever I try to run the following code
from openmm.app import *
from openmm import *
from openmm.unit import *
from sys import stdout

pdb = PDBFile('input.pdb')
forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME,
        nonbondedCutoff=1*nanometer, constraints=HBonds)
integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
simulation = Simulation(pdb.topology, system, integrator)
simulation.context.setPositions(pdb.positions)
simulation.minimizeEnergy()
simulation.reporters.append(PDBReporter('output.pdb', 1000))
simulation.reporters.append(StateDataReporter(stdout, 1000, step=True,
        potentialEnergy=True, temperature=True))
simulation.step(10000)

It shows
No template found for residue (number) (residue_name)

for every .pdb file that I am trying to run and every ForceField parameter permutation.
As an example I tried to run simulations on this protein
https://www.rcsb.org/structure/1a3n

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

when you say "force field", can you be more specific? it seems like a force field is a step in your workflow rather than your goal. "generate conformations for the proteins with arbitrary residues binded to it": what does this mean? What exactly are you modeling?  what exactly are you doing to your pdb file? can you give examples?

Do you mean "electron density clouds"?


text:
What is the troubleshoot for this error: conversion of .SRA to FASTA file on command prompt?
I am getting this error message after using the following code:
C:\sratoolkit.3.0.7-win64\sratoolkit.3.0.7-win64\bin>fastq-dump --fasta SRR1658345

Error:
2023-12-11T06:08:04 fastq-dump.3.0.7 err: timeout exhausted while waiting condition within process system module - failed SRR1658345

=============================================================

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

=============================================================

2023-12-11T06:09:46 fastq-dump.3.0.7 err: connection failed while opening file within cryptographic module - Cannot KStblHttpFileTimedReadChunked: to=0

What should be done?
Link to ncbi error report
Edit:
I am on Windows, using Command Prompt.

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

Please [edit] your question and explain what you are doing. First of all , tel us your operating system. Are you running Windows or some reason? If yes, what exactly are you running? Is that WSL? Or are you running it from a Windows shell? Also, please upload that error report to something like https://pastebin.com instead of making it a pdf. We need to be able to copy text from it and that site doesn't allow copying.

Isn't the syntax    --fasta <line-width> where <line-width> can be                                      "default" or "0" for no wrapping"

@Sam I don't think so. I just tried running fastq-dump --fasta SRR1658345 on my Arch Linux machine and it brought me a fasta file just fine.


text:
Uploading manifest from CCDI to CGC project
Thank you in advance for your help!
I am trying to download data from the CCDI (Children Cancer Data Initiative) and I have a manifest csv file that I downloaded from the CCDI. I'm trying to upload this data to my CGC project, but there is not option to import data from manifest file from the CCDI specifically, and when I do try to upload it, only the manifest csv is added and not the files. What is the process for uploading manifest that are not from one of these databases?

Thank you so much!

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

Please [edit] your question and give us (a lot) more context. What is this screenshot from? Where could we go to find it? Why are you using "import", is there no "upload" option? Do you have a manifest file or a csv file? Can we see an example? What is CGC? Do you mean https://cgc.umn.edu/? Or maybe https://cancergenomics.org/? What would uploading these files actually do? What kind of data do you have?


text:
What method uses "diverge probability" in RNA-Seq differential expression analysis?
I am replicating differential expression analysis from a paper (here) on RNA-Seq data. I extracted the regulated genes using DESeq2::results with a threshold of 1 for log2 fold change and 0.05 for the p-value. The results are very different from that of the paper. The paper does not mention much information about their statistical methods except for this paragraph:

Finally, the libraries were sequenced using the BGISEQ-500 platform. Differentially expressed genes were identified based on fold change ≥ 2 and diverge probability ≥ 0.8. The functional categories of the differently expressed genes and the pathway were obtained using MapMan as the classification source.

Now, I have no idea what they mean by "diverge probability" but it looks like another method to get the regulated genes. My questions are:
Question 1: What could be the method they are using in their analysis?
Question 2: Are the results from the "diverge probability" method different from the ones we get using DESeq2::results?
Question 3: If the answer to Question 2 is yes, how do we get from one to the other?
I appreciate any insights or experiences you can share regarding this approach. Thank you!

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

I'm pretty informed about stats but haven't a clue what "diverge probabilities" are. Are you referring to Kullback–Leibler divergence by any chance? Whilst I understand the question, at present its unanswerable but the metric doesn't appear to correspond to a known critical value approach.

@M__ I tried to add more info and edit my questions to be more clear. Hope that helps, and let me know if you want me to add more info.


text:
Can I interpret 15x coverage area as a part of one strand is missing?
This is a BAM file from a 30x coverage Nebula Genomics test opened with IGV.
I found an area (1:30169255-30169600) where the coverage is less than 15x. Can I interpret it as a part of one strand of Chromosome 1 is missing?

UPDATE:
Seeing a different question where similar vocabulary is used: it would not be one strand that would be lost. single stranded DNA is not very happy in the nucleus, and would either lead to something seriously bad or be gap-filled to make it dsDNA again.
It would instead be that one of two chromosomes has a deletion at that position, which would be a classic CNV.
ORIGINAL ANSWER:
No.
There are a variety of copy number variant (CNV) calling tools that put a lot of effort into statistically robust ways of asking this question, all of which are a lot more complex than just looking at raw coverage. I'd recommend looking more into them.
In terms of explanations for lower coverage in some regions that aren't CNVs: sequenceability of DNA is not even, replication of DNA is not even, mappability of reads is not even, etc. You need to reject those hypotheses before you go to variants.


text:
Gene list from eclip-seq .bed files
I have some eclipseq data, three IP and three input datasets, and have the alignments for all three. I used Bedtools "Intersect Interval" with the genome gtf file (hg19) I used to align them with, and that creates a bed file with gene names for each peak. I then wrote an R script to sum the clipper score for each gene within a sample, then do a simple student's t test across groups. I now have a gene list and p values for differential expression, but not sure if this is correct.
I've tried looking at a lot of sources but cant really find any canonical way to generate a gene list from bed files.
I've also tried IDR and it doesn't work on the clipper output, unfortunately.
Thanks!
Edit:
eCLIP-sequencing, it's essentially a variant of iCLIP-sequencing and data analysis methods are transferrable, to my knowledge.
I have alignments (bam files) made with Star aligner against Hg19. After deduplication and sorting, I called peaks using Clipper (a peak-caller made by the creators of eCLIP). It outputs bed files with scores for each sample.
My question is how exactly do I turn this bed file with scores into a gene list? Is my method above correct?

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

What is "eclipseq data"? Can you please provide some example input files? Anything you can add to clarify the most important problem you're having will be useful for people who are answering this.

Added more information. Thanks!

As someone who works with eCLIP data from ENCODE, I think I can help out.
You need to think of this more similar to ATAC or RNA-Seq, whereby you actually get a counts matrix, normalize that, and then apply math.
Another good idea would be to look through this eCLIP analysis document provided by Gene Yeo of ENCODE for more specific ideas.
https://www.encodeproject.org/documents/1f171ac6-a36a-41ac-b632-741aeb47aad2/@@download/attachment/eCLIP_analysisSOP_v2.3.pdf

As others have said it's not very clear what you're asking here, but I'm assuming you want to know how to determine relevant genes that were enriched in your IP.
As per the group behind eCLIP, Van Nostrand et al. 2020
: "standard peak analysis used the set of peaks identified as irreproducible discovery rate (IDR) reproducible and meeting fold-enrichment (≥ 8-fold) and significance (p value ≤10−3) in immunoprecipitation versus paired size-matched input."
If those cutoffs are too stringent you could change them to fit the context of your project. And what do you mean by differential expression? It sounds like you're addressing binding.


text:
LD local clump-GRCh38 bfile
I want to do LD local clumping.
The tutorial (https://mrcieu.github.io/ieugwasr/articles/local_ld.html) teaches how to, but the downloaded bfile is GRch37.
I followed this question to download GRch38.
LD clump GRCh38 GWAS results
I downloaded and it is "all_hg38.pgen". How can I transform it into bfile? Thank you.

I have the same problem when facing the analysis of eQTL data for clumping and I do not find any answer to solve it. So I try to find a method to transform 1000 Genomes phase 3 ailgned by hg38 data which is pfile format (pgen, pvar, psam). I found a method at: https://cran.r-project.org/web/packages/snpsettest/vignettes/reference_1000Genomes.html. This is a method about how to extract the EUR population and transfrom to bfile format using PLINK 2.0. After this transfrom we can use PLINK 1.9 to clump the GWAS data. I think this is the method which we want to use and I will try it. Hope this can help you.


text:
Featurecount output .txt file from bam file
After running this command:
featureCounts -p -s 2 -a $genome -o $dir"/"$specie"/count_table/"$value".txt" $output_loc -T 8

on the BAM output of HISAT (output_loc), I got this output:

Now, does the last column of this file show the gene count value?
And if I merge the same column of different samples, how should I go through the normalization step?
Thanks a lot for your help!

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

The first 6 column in standard featureCounts output represent what is in the column names. All columns after than (starting at 7) represent the counts for the sample(s). If you use a single bam file as input then it's one column, if you use many bams as input then it is one column per bam. If the counts are gene level, exon transcript or anything else depends on how you run the tool. Typically, with default settings, it is usually the raw gene level counts.
Generally, you can sum up technical replicates (=lane/sequencing replicates to increase depth) prior to running normalization.

Thank you so much I appreciate your help!
My other question is that since I didn't know that I could give multiple bam files to the feature count I instead ran the feature count on each bam file and then merged column 7 of the output text files and made my final count table in which rows are genes and columns are samples and cells are the raw gene counts. Now, as I know I should do normalization at this step, would TPM be a good method? if it's so off the table can you please let me know the sources I should use to learn how to proceed?
Thanks a lot!

You simple list all bam files to be processed.
featureCounts (options...) one.bam two.bam / or just *.bam
Read manuals of edgeR or DESeq2, use these normalizations. TPM is poor. Many posts on this at biostars.org why that is.

Thanks, ATpoint!


text:
How to prove correlation between gene expression and functions using omics data in bioinformatics?
I am new to bioinformatics, coming from a background in mathematical biology. I am currently participating in a project where the Principal Investigator (PI), who is a biologist, has tasked me with proving a correlation between the expression and functions of two or more elements. Specifically, one of them is a kinase enzyme and the other is mTOR. The PI has asked me to use omics data or databases for this purpose.
As someone with limited experience in bioinformatics, I am struggling to understand how a biological pathway or correlation can be definitively proven, beyond basic statistical methods. Additionally, the ultimate goal of this research is to develop treatments using substances that our lab has discovered, which have shown promise in killing specific cancer cells in murine and cell lines (although our data on this is limited).
My question goes as followes- is it possible to demonstrate a connection between these elements using only databases and omics methods? Why might a dynamic modeling approach, such as ordinary differential equations (ODEs) not be enough?
I am feeling quite confused and would greatly appreciate any guidance or direction on how to approach this problem effectively.
Thank you.

In general terms data mining is needed, then supervised learning and finally ODE, because you are performing DE against enzyme kinetics and thats doable, but you'll need a lot of biochemical knowledge.
I would first take the 'omics' data sets and perform a PCA analysis on it and see where your mTOR and kinase genes are located. PCA is known as unsupervised learning and can be accessed by Seurat or better using Scanpy packages. That will help establish the overall behaviour of your key enzymes pre- and post- treatment ... against the control.
Basically you want:

Data sets of cancer mice - pre-treatement
Data set of cancer mice under treatment (drugs targeting kinases)
Healthy controls

See what unsupervised learning looks like between the 3 groups and where your enzymes are. The PCA plots then compared for the 3 groups and you want to see the targets you are interested in changing groups, especially the treatment group shifting between the control and the diseased unsupervised plots. This would provide a basis for exploring supervised learning (outlined below).

Supervised learning What you will do is look at gene expression and train it against your phenotype of interest (cancer and recovery from cancer). The training 'target' is the drug the mice are being treated with for each given data set, obviously they'll be a "untreated" and "control" categories. If your model trains correctly >0.85 to 0.9 accuracy ... post-back and I'll explain how to extract the answer. This approach will extract the biochemical pathway for how the cancer treatment works, but requires more work than described here and that provides the framework for the DE.
Again, a lot of serious knowledge of biochemistry is needed to deploy ODE, unless you've got a brilliant a priori - then you could dive straight in. To do that you'd need to know the exact mechanism of how these enzymes are involved in cancer and then parse that data from your data sets (totally different approach to what is described above). The problem with this approach is if someone knows both the maths and genomics (like a reviewer or examiner) and you've made an error, that will that will be a problem. If you've been through an orthodox approach (above), its a lot stronger justification for the ODE, thus difficult to criticise.
Note The key is what data set types comprise 'omics', I'm assuming this is RNA-based (thats normal), but it might not be.


text:
How to build hisatgenotype extra index files (.link, .locus, etc)
I'm trying to run hisatgenotype on a custom genome, originally specified in fasta format. I'm able to use hisat2-build to build the .1.ht2, .....8.ht2 files, however I've not been able to figure out how to get the other files that hisatgenotype requires. I've looked into hisatgenotype_toolkit build-genome but it errors out due to problems with argument parsing. Is this even possible, or do I have to use the genomes/indices downloaded by the setup.sh script?
edit:
The error I get is the same as this issue: https://github.com/DaehwanKimLab/hisat-genotype/issues/69
% hisatgenotype_toolkit build-genome --base hla
Traceback (most recent call last):
  File "/Users/ted/miniforge-x86/envs/hisat/hisat-genotype/hisatgenotype_tools/hisatgenotype_build_genome.py", line 546, in <module>
    typing_common.clone_hisatgenotype_database()
  File "/Users/ted/miniforge-x86/envs/hisat/hisat-genotype/hisatgenotype_modules/hisatgenotype_typing_common.py", line 52, in wrapper
    func(*args, **kwargs)
TypeError: clone_hisatgenotype_database() missing 1 required positional argument: 'destination'
Script exited with error 1

If you try to run with an extra positional argument it fails too
% hisatgenotype_toolkit build-genome --base hla indices
usage: hisatgenotype_build_genome.py [-h] [--base BASE_FNAME]
                                     [--locus-list LOCUS_LIST]
                                     [--inter-gap INTER_GAP]
                                     [--intra-gap INTRA_GAP]
                                     [--aligner ALIGNER] [--linear-index]
                                     [--commonvar] [--clinvar] [-p THREADS]
                                     [-v]
hisatgenotype_build_genome.py: error: unrecognized arguments: indices
Script exited with error 2

And similar problems if you add --destination indices

"but it errors out due to problems with argument passing" - show us your exact command and the error being displayed.


text:
Loop and TAD calling from HiC reads using HOMER
I have HiC reads from an HiC experiment and I am using HOMER to call loops and TADs. The numbers I am getting are like:
Total TADs: 5024; Total Loops : 254
As multiple loops are expected in a TAD so numbers of loops should be higher than the number of TADs. That's why I am not sure how to justify the numbers I am getting.
I made sure I am not doing any mistake in the HOMER pipeline of loop/TAD calling and I followed all the steps in the tutorial here: http://homer.ucsd.edu/homer/interactions2/index.html Unfortunately I don't have any other loop calling program installed right now and HOMER is easy to install and run.
If someone working with HiC data analysis and TAD/loop calling from HiC reads can please comment on this?
Thanks
Following is the set of commands I used (based on HOMER tutorial here http://homer.ucsd.edu/homer/interactions2/index.html) to get loop/TAD calls from HiC reads.
#Trimming/processing of HiC fastq files:
homerTools trim -3 GATC -mis 0 -matchStart 20 -min 20 Hu_concat_1.fastq
homerTools trim -3 GATC -mis 0 -matchStart 20 -min 20 Hu_concat_2.fastq
#Aligning the two paired-end fastq files to get sam:
bowtie2 -p 20 -x hg38.fa -U Hu_concat_1.fastq.trimmed > Hu_concat_1.sam
bowtie2 -p 20 -x hg38.fa -U Hu_concat_2.fastq.trimmed > Hu_concat_2.sam
#Making tag directory (the HOMER makeTagDirectory program is used to process the alignment file into a HOMER-style tag directory which is used for subsequent analysis steps)
makeTagDirectory Hu_Hic_TagDir/ Hu_concat_1.sam,Hu_concat_2.sam -tbp 1
#Calling loops and TADs
findTADsAndLoops.pl find Hu_Hic_TagDir/ -cpu 10 -res 10000 -window 15000 -genome hg38.fa
The whole process ran smoothly without errors. It is just that in the output, the number of loops I am getting are far fewer than no. of TADs (Total TADs: 5024; Total Loops : 254) which I think is not biologically correct. If someone using HOMER for calling loops and TADs from HiC reads can justify this output, that would be really helpful.

I haven't looked into the scripts of Homer on what loop size they're setting to detect loops. Have you considered looking into loop size? Recently, I used Chromosight to detect loops and I was able to find more loops when I increased loop size.

Without being terribly familiar with your data or the tool or the tutorial, I think that there is a simple technical explanation: your statistical power to detect TADs will be much higher than your power to detect loops. I'm not sure how many reads you're working with, but in general, TADs are very big strong features and loops are smaller features. Bigger features will have more reads due to just being bigger pieces linear chromosome, and stronger features with more robust chromatin structure will show up better in Hi-C data than the more labile loop structures.
See for example this Figure from a review showing qualitatively the expected strength of each feature. You just have a lot more data supporting TADs than you do loops!
That said, I have not been able to track down a citation explicitly saying so- but intuitively it makes sense to me that TADs would be easier to define than loops.

You can evaluate the power of your dataset here for detecting loops. If you find that it has low power, that is probably the source of your small number of loops detected.
Note also that different tools may have very different power thresholds! Trying a few other loop detection tools may lead you to the conclusion that Homer just isn't very sensitive. There is indeed some reason to believe that Homer is not sensitive in comparisons of different loop detection tools.
You don't mention how much data you're working with or what kind of hi-c it is, so knowing that may partially answer your question.

Generally, loops are harder to identify than TADs because of the length scales and data involved.
TADs are often identified by detecting changepoints in Hi-C contact matrices averaged over Mbp along the diagonal, whereas loops are kbp scale increases in contact matrix values at some unknown distance away from the diagonal.
In practice, this makes TAD detection often more sensitive than loop detection, regardless of the underlying biology.
Theoretically, there should be multiple loops within each TAD, but Hi-C data isn't exactly optimized for detecting loops.
Maximilian's answer also raises a good point that HOMER isn't that sensitive of a loop caller.
Other methods like Fit-Hi-C (1 or 2), Peakachu, and Mustache are largely better at loop calling, so you could try those and see if you get better results.
Loop detection is also very dependent on the library size for each sample, whereas TAD detection seems to be a bit more robust.
I'm not sure how many reads you have per sample, but that may be something to consider.
Overall, loop detection is pretty unreliable for representing the underlying biology of cis-regulatory interactions.
It's one of the best genome-wide methods we've got, but is still unreliable.
So if you want to draw inferences from the number of loops called in different experimental conditions, I would be careful, here.
There are lots of caveats and uncertainties that just come from the detection / measurement part of your experiment.


text:
Finding homologous regions in multiple whole genomes
Right now I have six genomes that I want to compare and identify homologous regions in the genomes. I have run nucmer, show-coords and obtained the output files. An example is shown below with Genome 1 vs Genome 2. More files go through Genome 1 vs Genome 3 until all are compared.
I was trying to use bedtools merge after concatenating the multiple output files from show-coords command to find the homologous regions. However, I am not quite sure on the command to use. Any help is appreciated.
Edit: added commands. I just want to know what regions are homologous across the six genomes (99/100% IDY). I don't know what I am supposed to merge to find this out. I was reading and I believe bedtools should be but I don't know. I am usually a wetlab but I try to do a bit of bioinformatics.
 # Loop through all pairs of genomes
for i in {1..6}; do
    for j in $(seq $((i + 1)) 6); do
        genome1="genome${i}.fasta"
        genome2="genome${j}.fasta"

        output_prefix="output_${i}_${j}"

        # Run NUCmer
        nucmer -p "${output_prefix}" "$genome1" "$genome2"

        # Show coordinates
        show-coords -r -l -c -T "${output_prefix}.delta" > "${output_prefix}_coords.txt"
    done
done

After using the above I then concatenate the files with
cat output_*_coords.txt > combined_coords.txt

This gives me something that looks like this below.
[S1]    [E1]    [S2]    [E2]    [LEN 1] [LEN 2] [% IDY] [LEN R] [LEN Q] [COV R] [COV Q] [TAGS]
1       14347   1       14347   14347   14347   99.02   5064019 5004926 0.28    0.29    Genome1      Genome2
14474   32175   14350   32051   17702   17702   99.33   5064019 5004926 0.35    0.35    Genome1      Genome2
16541   22042   4698474 4692975 5502    5500    99.64   5064019 5004926 0.11    0.11    Genome1      Genome2
16548   22042   3728461 3722968 5495    5494    99.84   5064019 5004926 0.11    0.11    Genome1      Genome2
16583   22041   3904238 3898781 5459    5458    99.65   5064019 5004926 0.11    0.11    Genome1      Genome2
18727   22097   4426876 4423506 3371    3371    99.73   5064019 5004926 0.07    0.07    Genome1      Genome2

Please [edit] your question and show us the exact commands you used. All of them. Including how you concatenated the files. We need to know exactly what each of the columns is and then we need to see or at least understand your expected output. You only say you want to use bedtools merge, but you don't tell us what you are merging or what you want to do. When should a line be merged? When it is 100% identical? More than N% identical? Any size limits? How much identity across genomes? Please explain what you need in more detail.

According to show cords the format you are getting contains columns which you can use to obtain BED format with relevant information.
You need following columns:
TAGS_genome1
[S1]    Start of the alignment region in the reference sequence.
[E1]    End of the alignment region in the reference sequence.
TAGS_genome2
[% IDY]

I would start with not merging all coordinates but getting bed files for each pairwise genome_a__vs__genome_b.bed etc. and start from there, taking a subset of 3 genomes first (say a, b,c).
You can find the regions shared by all three bedtools intersect:
bedtools intersect -a a_vs_b.bed -b a_vs_c.bed

The issue is that you will have many to many relationships, since region_1 from b can have several similar regions in b, and each one of these (== regions in b) also can be similar to a number of regions in a. Depending on what you are looking for and how dis-similar the genomes are you may end up with reporting repetitive elements or if you filter hard missing i.e. exons.

Thanks @darked89. That's exactly what I needed, those 5 columns. I didn't consider that there are too many relationships.
I will then continue to do comparisons. I think will give me the answer that I want.
bedtools intersect -a a_vs_b.bed -b a_vs_c.bed > a_vs_b_vs_c.bed
bedtools intersect -a a_vs_b_vs_c.bed vs a_vs_d_vs_e.bed > all_compared.bed


text:
split fastq file containing a sequence block at different locations
I have some fastq files (obtained from nanopore sequencing) that contain reads that can be of either of these 5 forms:

a known CDS with 3'UTR:

CDS----------------------- (seq1: original sequence)

the same cds with a block of 50bp represented by ||||| (that substitutes the 3'UTR at 4 different locations)

CDS|||||------------------ (sub1: substitution 1)
CDS-----|||||------------- (sub2: substitution 2)
CDS----------|||||-------- (sub3: substitution 3)
CDS---------------|||||--- (sub4: substitution 4)

I need to split my fastq files in 5 smaller files clustering the different types of sequences (seq1, sub1, sub2, sub3 and sub4). The problem is that those sequences all share some segments and I don't know which method I should apply. I tried finding the best match by building a BLAST database with the 5 sequences but as it is local alignments, I get HSPs for each sequences and cannot discriminate between them.
Does anyone has an idea how I could achieve that?
Any help would be amazing.
edit
Here are 5 sequences as I expect them (i.e., fasta reference)
>seq1 (this one does not have the 50bp block)
ACCTGCCAGTTCCTGAAGGGTGTACTGATCCTGTGGCTGAAAACTTTGATCCAACGGCTAGAAGTGACGATGGAACCTGTGTCTACAACTTTTGAGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCCTATTGTCTCTTTTGGTTTATTTTTCTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTTTTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTCATAAA
>sub1
_GTAAGTGAGGGAGTCTCGGACTCAATTAGAGGCTCTTCTTTCACA_ TTGATCCAACGGCTAGAAGTGACGATGGAACCTGTGTCTACAACTTTTGAGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCCTATTGTCTCTTTTGGTTTATTTTTCTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTTTTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTCATAAA
>sub2
ACCTGCCAGTTCCTGAAGGGTGTACTGATCCTGTGGCTGAAAACT _GTAAGTGAGGGAGTCTCGGACTCAATTAGAGGCTCTTCTTTCACA_ TTTGAGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCCTATTGTCTCTTTTGGTTTATTTTTCTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTTTTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTCATAAA
>sub3
ACCTGCCAGTTCCTGAAGGGTGTACTGATCCTGTGGCTGAAAACTTTGATCCAACGGCTAGAAGTGACGATGGAACCTGTGTCTACAACT _GTAAGTGAGGGAGTCTCGGACTCAATTAGAGGCTCTTCTTTCACA_ TTGTCTCTTTTGGTTTATTTTTCTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTTTTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTCATAAA
>sub4
ACCTGCCAGTTCCTGAAGGGTGTACTGATCCTGTGGCTGAAAACTTTGATCCAACGGCTAGAAGTGACGATGGAACCTGTGTCTACAACTTTTGAGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCCTA _GTAAGTGAGGGAGTCTCGGACTCAATTAGAGGCTCTTCTTTCACA_ GATCTTGTCCTCTTTTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTCATAAA

examples (50bp block if present is in lowercase here)

this read should be associated with sub3 because of the flanking sequences (...CTACAACT in 5' and TTGTCTCTTTT... in 3' of the 50 bp block)

@e213fdd5-bbfb-4dfe-9b03-29eea256da0d
TGTGTACTTCGTTCAGTTACGTATTACTAGTTATTGAGTGTCTTTGTGTTTCTGTTGGTGCGTCTTCGCACAAGGCTAATCTTACTCTCCTCTCCAATGTACTGCTTCTATGTTCATCTCAGCCGAATTCAATAAGGAGAAGAACTTTCACTAGGTTTCCCCATTCTGTTCGATTAATTGGTGATGTTAGTAAATTTCAATTTTCTGTCGGTGGAAGGTGAAGGTGATTGCTGCTGGCTTTTCATTGTTTACTTTTGGCTTCTGTTTTGTAGTAACGCATCACTACTTTCTCTTATAGTGTTCAATGCTTTCCAAGTCTGGAATCATGAAGCAGCGGCTTCTTCAAGAGCGCCATGCCTGAGGGATACGTGCAGGAGAGGACCATCTTCTTCAAGGACGACGGGAACTACAAGACGTGCTGAAGTCAAGTTTGGGGAGACACCCTCGTCAATGAGATCGAGCTTAAGGGAATCGATTTCAAGGAGGACGGAAACATCCTCGGCCACAAGTTGGAATACAACTACAACTCCCACAACGTATACATCATAGCCGACAAAACCCAAAAGACGGCATCAAGAAACCAACTTCAGACCTTACCATCAGGCGGCGGCCATCACACAGTCATTGTCAACAAATGTAATTGGCGGTGGCTATTCTTTTTACCCAGACAACCATTACCTTGTCCACACAGTCACACTGTGAAAAGATCCCAACTTTGAAAGAGAGGCCACTGGTCTTCTTGAGTTTGTAACAACTGCTGGGATTACACATGGCATGGATGAACTATACAAACATGACAGACTCTAAAGCCTGCAGTTCCTGAAGAGATTGTACTGATCCTGTGGCTGAAAACTTTGATCCAACGGCTAGAAGTGACGATGGAACCTGTGTCTACAACTgtaagtgagggagtctcggactcaattagaggctcttctttcacaTTGTCTCTTTTGGTTTATTTTCTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTTTTGTTTTCCTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTCATGCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACTGTAGGCACCATCAATGAAGATAGAGCGACAGGCAAGTCACAAAAACACCGACAACTTTCTTGTCACGGTAGGCGATAGCAATACGTAACTG
+
$&&&%234441124--585545;>B99968650-++*()'&'(;D540488?=888,++*+*)*$$$##&%$$$%)*****('%$$'1/+++-,('))'%$$%$$$&)&%%%&(*'&),/'''&&'06872..*&%)),,3/49:;<?A773)('()/066+'&&)('(+&%$&''$$&&()*--888)('()'')')'''++*+)&')&$$$%()'')*''&&))&'&$&(&$#$%&%&&&%%$%&..&%%%%&&%$%&'&$$'-)(%%%$&(*,+****&'''('+17<{..+*('((*+)*&&&'**.*%$&%&&%$$$%&&'&&(+'(&&&%%&&($$%&&'+55536988;;74632+(*)''%6<=(((&*''(68998===?ABBBB=B?<;;96667>>?>>?@?A=-+.7556<@>A<0//002(()-989::::;<<>>;:5++++->;::;<?F:7<()988;;<?>9;9)36:5-963/002B@A@@?554543,,9<=<@?@>777461166.../B88<;<@=<:888;;>42(''%%(%%+*('&%%&&**)(()*+020*'('&$$$'('''&%'((%(&$##$&)*()&%%%$$#$%'&(&'%$$%%&'%%%$%%%$&('$##$&&&&$$$%$%$$%&%%&040(((&'.023266,,,++0)++-.-*,,0+&&'$&''''%%'&-+,0((+))*+*(%&)+00.3)''%'&'''(*(%%)++++,--+.....0'(()7:<=;=<<<<<0.<;9:;<@;;;:<<=;2+)++,1212611)())*:1++*'&&%((&&056/,454))*''''(3346;<;;;;;000.0/043546??=;<*)+64335663455665333788:78779<>9755446::<><198662*+7899722248867>?8++4))+6779<<BA><9;<;:;;{<?BB-,,-13444/154///933{3==@B:9989::::::;;<<=:3335=EBB;:*)..'-:69?@CKKKDC<<?>:)()((012.,,+,55;=?>@>:98:78544344*)%$$$&)))))))))(((((''''&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')+.1866677978:;;:::;;;99::988889653456:::9889;::8899.,+0668:;:::8436735//:632/-,+**+013459;9332/+($

this read should be associated with sub4 because of the flanking sequences (...TTTTACTCCTA in 5' even if a C is missing, and GATCTTGTCCT... in 3')

@107dfd6d-2101-4fd8-aeb5-2c432e176d72
ATGTACTTCGTTCAGTTACGTATTGCTATCGCCTACCGTGACAAGAAAGTTGTCGGTGTCTTTGTGTTTCTGTTGGTGCTGATATTGCATGAAGACTAATCTTTTCTCTTTCTCATCTTTTCACTTCTCCTATCATTATCCTCGACCGAATTCAGTAAAGGAGAAGAACTTTTCACTGGGTTGTCCCAATTCTTGTTGAATTAGATGGTGATGTTAATGGGTACAAATTTTCTGTCAGTGGAAGAGGGTGAAGGTGATGCAACATACGGAAAACTTACCCTTAAATTTATTTGCACTACTGGAAAACTACCTGTTCCATGGCCAACACTTGTCACTACTTTCTCTTATGGTGTTCAATGCTTTCAAGATACCAGATCATATGAAGCGGCACGACTTCTTCAAGAGCGCCATGCCTGAGGGATACGTGCGGGAGAGGACCATCTTCTTCAAGGACGACGGGAACTACAAGACACGTGCTGAAGTCAAGTTTGAGGAGACACCCTCGTCAACAGGATCGAGCTTAAGGGAATCGATTTCAAGGAGGACGGAAACATCCTCGGCCACAAGTTGGAATACAACCGCAACTCCCACAACGTATACATCATGGCCGACAAGCAAAAGAACGGCATCAAAGCCAACTTCAAGACCCGCCACAACATCGAAGACGGCGGCGTGCAACTCGCTGATCATTATCAACAAAATACTCCAATTGGCGATGGCCTGTCCTTTTACCAGACAACCATTACCTGTCCACACAATCTGCCCTTTCGAAAGATCCCAACGAAAAATTTTTGCATGGTCCTTCTTGAGTTTGTAACGGCTGCTGGGATTACATGGCATGGATGAAGAACTGCAACACAAACACTTGACGAACTCTAAACCTGCCAGTTCCTGAAGAAATTGTACTGATCCTGTGGCTGAAAACTTTGATCCAACGGCTAGAAAGTGACGATGGAACCTGTGTCACCAACTTTTGAGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCTAgtaagtgagggagtctcggactcaattagaggctcttctttcacaGATCTTGTCCTCTTTTGTTCCCTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTCATAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACTGTAGGCACCATCAATGAAGATAGAGCGACAGGCAAGTCACAAAGACACCGACAACTTTCTTGTCACGGTAGGCGATAGCAATACGTAACT
+
-0,+))),:::<??===>=<<<?@>=<<=?;:889:<;<<=;;;;CD@?:9:A>9999=999DBAA@<BE=3220111<B@@AA@A?>>>=>>;:9999:9;<E9A<;?????C8788;==@?;;;==9989@A6677?=7788+**+368>B87::==@<99;<<<<<<??>@=:9/.))*698.*,/8A=<==>>=?BAB=::::;;;<?@>@B@@?:888:>?@CCCEDB>??>===::.+)()2157=7:99<B?>???@??=;::=DGC=<6?;84568>AB@CCBCBAA====AAACEI{AA=<::;<?@=<<>=>@AAA@>>>>?<?====A?==<=@>??@>?==??CA@AAA@AB=757:965(*.;<;<::?AAA@>=>::::<=?@BCC5336899=?>>?@A??42220;9:89.-***+77;97<;:;<<??BC{878<:9:<:2112<<>===@=<;;<>???877/,.-.-///0,++,)299:=;:;;=@@<;=7778?=::;;<;=<<C>AB<BBBCEDCDC?A>;;9:655:DAA@?>>>=?=><;<<@=>><=B=@=:870++,,;=;>><<==CBA@?>=>>=667<;;7778B>B?B<<==<99:::4:659;6642212:.-<=??=6556;??@@>::;:==;8887899==;2143235<<<?>==?CEBD?>>>ADB=720/.02334<00//06568>:888:711*)-776'''&06:<=><><866;{C@A?>8897=>EK558568552-179664)(()),,--,()))3:;:<BAAA>=>>>><<>>*))+/245999;>>=96;?B9888<<77-,,,,&&%%%&&'+6;A;2111--,,.:;><<ACCA@@?=9;;;>=>==@930../+)(())07>>8<><;:;9238<>=>998878:7700016545**1236434689<9<<;.--0-,,,&'():78862237@???><;;;;<<=@>@??@@===;<=>>;:8;)(67511223/.---.2001023:;21<?=**6223;;>?AACGA=73366:;<=:766;ACC=<;86000?8=@DHEC@@=;;<9999<<==9><==8889>A>=;30266777:972+***)))'&%%%%%$$$$$$######"""""""""""""""""""""""""""""""""""""""""#&')+-1122/...46688<<=0('%%&258899;92136:999;:::9:;;--2-.343668=<?@EF????==:7668;6667;><<;;>;;8851.*

@terdon, I added two examples. I cannot provide all cases, as I actually have 16 of them and i'm having a hard time doing it by hand. I still hope that it is clearer now.

You can try seqkit locate with say start or end of CDS sequence. This will give you a tabulated output with locations of your CDS query in each read.
Since you expect that the CDS sequence should be located at one of the five different locations inside your reads you can figure out to which one of the five a given read belongs.
You may also use clumpify.sh from BBMap to cluster reads prior to searching for the CDS locations. The idea is that maybe the five groups will cluster in fairly well separated blocks, but this is fairly speculative.
edit
Since the alignment of your sequences looks like this:
CLUSTAL 2.1 multiple sequence alignment

sub1            NNNNNNNNNNGTAAGTGAGGGAGTCTCGGACTCAAT---TAGAGGCTCTTCTTTCACANN
sub2            --ACCTGCCAGTTCCTGAAGGGTGTACTGATCCTGTGGCTGAAAACTNNNNNNNNNNGTA
seq1            --ACCTGCCAGTTCCTGAAGGGTGTACTGATCCTGTGGCTGAAAACTTTGATCCAACGG-
sub3            --ACCTGCCAGTTCCTGAAGGGTGTACTGATCCTGTGGCTGAAAACTTTGATCCAACGGC
sub4            --ACCTGCCAGTTCCTGAAGGGTGTACTGATCCTGTGGCTGAAAACTTTGATCCAACGG-
                          **   *** **     * **  *  *   *  *  **

sub1            NNNNNNNNTTGATCCAACGGCTAGAAGTGACGATGGAACCTG---TGTCTACAACTTTTG
sub2            AGTGAGGGAGTCTCGGACT-CAATTAGAGGCTCTTCTTTCAC---ANNNNNNNNNNTTTG
seq1            --------------------CTAGAAGTGACGATGGAACCTG---TGTCTACAACTTTTG
sub3            TAGAAGTGACGATGGAACCTGTGTCTACAACTNNNNNNNNNNGTAAGTGAGGGAGTCTCG
sub4            --------------------CTAGAAGTGACGATGGAACCTG---TGTCTACAACTTTTG
                                              *                          * *

sub1            AGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCCTATTGTCTCTTT-----TGGT
sub2            AGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCCTATTGTCTCTTT-----TGGT
seq1            AGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCCTATTGTCTCTTT-----TGGT
sub3            GACTCAATTAGA--GGCTCTTCTTTCACANNNNNNNNNN--TTGTCTCTTT-----TGGT
sub4            AGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCCTANNNNNNNNNNGTAAGTGAG
                  *   ****    *  * **  **  *                            **

sub1            TTATTTTT---------------CTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTT
sub2            TTATTTTT---------------CTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTT
seq1            TTATTTTT---------------CTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTT
sub3            TTATTTTT---------------CTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTT
sub4            GGAGTCTCGGACTCAATTAGAGGCTCTTCTTTCACANNNNNNNNNNGATCTTGTCCTCTT
                  * * *                *** * * *   *          **************

sub1            TTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTC
sub2            TTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTC
seq1            TTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTC
sub3            TTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTC
sub4            TTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTC
                ************************************************************

sub1            ATAAA
sub2            ATAAA
seq1            ATAAA
sub3            ATAAA
sub4            ATAAA
                *****

instead of CDS end/start you have to select regions around your deletions with say 15-20bp flanks on each side. Put these in fasta and run locate.
caveat
running with zero mismatches will miss some reads.
edit 2
The command:
zcat input.fq.gz | seqkit locate --threads 8 \
--max-mismatch 0 \
--pattern-file fasta_with_5_frags_extracted_from_algnt.fa > result.tsv

Some trial and error is required to select good fragments matching exactly one of the five groups plus increasing max-mismatch if due to sequencing errors a lot of reads were dropped.
To speed up IO one can compress result.tsv on the fly, i.e. using pigz or zstd.
In the next step we need to create five separate files listing read names containing a given fragment:
grep frag_name_1 result.tsv | cut -f 1 > frag_name_1.read_names.txt

seqkit grep --pattern-file frag_name_1.read_names.txt input.fq.gz -o frag_name_1.reads.fq.gz

# repeat five times with different pattern file and changing the output

Thanks for all the tips. I actually think that this method can work.
I have tried to use seqkit locate to find the 50bp block and ran it a second time to locate the end of the CDS. Then, if I compute the distance between the two, I get distances roughly multiple of 50, which allows me to distinguish the sequences group I am after. I am going to try and make something quite clean and then I will post here my solution!


text:
How to align two FASTAs and extract the aligned part?
I have two fastas and I'd like to align them to extract the aligned part. I've just found examples (Several of them) of alignment of FASTQs using FASTA references, so I'm not sure if what I want really makes sense.
Just to give a context, I want to find the alignment between the FASTAs, then extract the aligned part to feed a classifier that will use this data to estimate the DNA age based on the age range labels of the fastas I have in the database. I'm currently having an error larger than the expected, and I'm trying to find out if it is caused by the differences between the sequences, so I'd like to use just the aligned part in the classifier.
Is it possible/Does it make sense to align the FASTA sequences, and which would be the best tool for that?
EDIT: All the sequences are from human mitochondrial DNA. We have some complete modern samples and some ancient non-complete, damaged samples. They are expected to be similar, but we are expecting that the age range estimations can be derived from the DNA modifications that happens with time. We have several samples for validation, we expect most of the validation samples to be partially aligned, since we have several ancient damaged samples. So we discard the parts with poor alignment and hopefully we'll have a more reliable material for the classifier.

I've added the information you recommended. Well, maybe the aligners have some threshold so we discard the parts based on the quality? Does it exist?

Given that you're aligning human mitochondrial DNA sequences, consider aligning each sequence to the Reconstructed Sapiens Reference Sequence (RSRS) as a base for alignment:
http://phylotree.org/resources/RSRS_annotated.htm
Aligning to a well-curated outgroup should give you a better estimate of DNA age, and make it easier to compare variants.
My preference when aligning large sequences is to use LAST, because it produces good base-level matching, and can be trained and/or tweaked as needed to improve alignment accuracy in unusual situations. LAST will produce alignments including [only] the aligned part as its default MAF output:
https://gitlab.com/mcfrith/last
Information on how to use it can be found in the cookbook. It generally follows the process for most aligners:

Generate an index for the reference sequence
Align the query sequences to the generated reference

In the case of LAST, an additional training step can be added after index generation; the process is as follows:
$ lastdb RSRS.fa RSRS.fa
$ last-train RSRS.fa query_sequences.fa > trained.mat
$ lastal -p trained.mat RSRS.fa query_sequences.fa > query_aligned.maf

When a phylogenetics calculation goes wrong the bug is almost certainly the alignment - so that bit of the question I understand. This suggests a machine learning algorithm is being used. In that case it's the vectorisation of the alignment. If thats the program then the alignment needs to absolutely as occurred in the original training data (if its ML). That could be the problem.
The key to the question is understanding the "classifier" because its understanding the input into that algorithm.
Anyway, any aligner will do MAFFT or Muscle 5 -superfast option. Simply input sequences in fasta format and click or activate the "align" function. Including an example of the original sequence the ML (you mentioned it was a "classifier") was based on is very good idea.
Automated editing is via trimal using -nogaps option. However, you should manually inspect your alignment and edit it prior using trimal and Sequtron (Macs) is cool. Manually editing the alignment is standard practice.

Do you recommend any aligner in specific? Some that I could use in command line to add it to the pipeline in my python code. I've tried to use bowtie2, but it gives me an error: "Error: reads file does not look like a FASTQ file"

This is an "assembler" not an "aligner" in my definitions, bwa-mem2 is a highly memory efficient assembler. bowtie2 hasn't been recommended for some time.

bowtie2 and bwa-mem are aligners, not assemblers

Not in my field. An aligner is MAFFT, Muscle, Clustal.


text:
Is there some tool so find run accessions metadata?
I have a project in which I'll use some FASTQs to feed a neural network and, as label, I'll use the estimated age range. I've got the data from AMTDB, and there they have the metadata.txt, in which we are supposed to have the metadata, including this specific data I need.
Nevertheless, most of the available FASTQs are repeated in every project id (i.e., I have project id I363 and I250 and both have the FASTQs with run access from ERR001 to ERR100).
That happens for the ~300 project ids available. That happens although the metadata of each project id indicates a different range, so I infer that the age range corresponds to a specific subset of the FASTQs.
What I'd like to know is if there's some tool available to get the metadata specific to a run access (Just like with fastq-dump it's possible to download them by RA).
Thank you in advance for any insight.

What is the AMTDB? Can you link us to it? Or at least to this metadata.txt? What age range are you referring to? Is this the age of the individual who was sequenced? The age of the file? The sequencing date? We really need to see the data to understand what you ask. Why don't you just give random names?


text:
Is it possible to extract deamination patterns from VCF files?
I have some FASTQs from ancient DNA, and I can run them on mapDamage to get the damage profile and use the parameters generated as input to my ML model that try to infer samples' ages.
I also have access to VCF files from ancient DNA. My doubt is if the variant calls from ancient DNA samples in a VCF file (Or other equivalent variant calling formats) reflect the deamination pattern (C to T substitution) that is related to the sample age. If yes, I could try to use these, along with their positions on the chromossomes, to try to infer the age of the samples based on these features.
This is the data I have (Currently transformed to VCF format locally): https://reich.hms.harvard.edu/allen-ancient-dna-resource-aadr-downloadable-genotypes-present-day-and-ancient-dna-data
It is the largest dataset I've found, so I'd really like to work around it as much as possible to extract the most of it! Any tip is welcome!


text:
How to find a database for proteins that bind to different molecules and have different binding sites?
I'm doing a study in which I'm trying to analyze the effect of molecules (ligands) on binding sites of proteins.
Therefore, I'm trying to find a database in which the proteins bind to multiple different molecules, and as a result have different binding sites. The database that I'm working with currently is from UniProt, and the problem with it is that it has proteins that bind to only a small class of molecules, and hence the binding site is similar for each protein-ligand pair.
For reference, here is the database that I'm using right now: https://bit.ly/ProteinLigandBindingSite
Can someone help me find a database which contains proteins that bind to different molecules, and hence have different binding sites?
Thank you!

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


text:
Why is my phylogenetic tree always splitting into one leaf and one node?
I used RAxML software on about 300 species mtDNA region (COI gene).
Unfortunately I used a highly variable region so my tree only has highly support values for branches close to the leaves.
But I also noticed that all the nodes split into one leaf and one other node... Why is this happening?
Isn't it common that a node will simply split into 2 other nodes?
I am using the GTRGAMMA function in RAxML (I don't know if that may be the cause).
Do you have any suggestions on this topic? Should I be suspicious of my tree if I don't see any node splitting into 2 other nodes?
More info

Alignment algorithm: MUSCLE
Graphical example:
My tree looks like the one on the left. Every bifurcation links to a tip/leaf (terminal node) of the tree and a non-terminal node of the tree.  It sounds suspicious to me that I never observe a bifurcation that splits into 2 non-terminal nodes of the tree (figure on the right)

What do you mean by "But I also noticed that all the nodes split into one leaf and one other node... Why is this happening?" ?

Anyway if you can provide a clearer explanation I can probably answer your question @Mirk

What I suggest is the description sounds like a polytomy (star tree). In a phylogenetic tree, a polytomy is a node that has more than two immediate descending branches.
A polytomy can have many branches originate from one node. In contrast, a bifurcating tree only two branches stem from one node.
A polygomy simply means that the evolutionary relationships can not be fully resolved to dichotomies, i.e. bifurcations. There's nothing wrong with this result and it makes reading the tree much easier. The model you describe looks fine.
The only thing to watch if you're not looking at a dendrogram, but are looking at the phylogram. I doubt thats the explanation in any case.

From the update. The figure you presented is an asymmetry phylogeny and is perfectly normal, for example this is observed in molluscs, described in Palmer 1995.  If a species is moving directionally e.g. east to west across a massive landmass and at uniform geographic intervals establishing an ecological niche and forming a new species this will result in an asymmetric tree and might form a "cline".

Thanks but my tree is fully resolved so i don't have explicit polytomies. I tried to create an example figure to clearly show what i mean.


text:
Hybrid assembly versus polishing for hifi and illumina reads
I will have to carry out a project of assembly using hifi reads for which I have already illumina reads and I am wondering which of the hybrid assembly or polishing would be the best option for this kind of data? I wanted to use hifiasm which seems the gold standard and then nextpolish2 but I am wondering if using for example spades to do an hybrid assembly would lead to something better or worse?
This question is only intended to see if someone did such benchmarking before or not.
Best regards,

Usually, if you have HiFi reads, there is no way you can improve the assembly with short reads. Not just that they don't have higher quality, but they also map a lot worse. People sometimes use short reads for genome profiling and QCing the whole thing using Merqury.
Not sure if there is benchmarking published though.


text:
Does the value of RMSD have any impact on my vina results? Are results with RMSD greater than 2 unreliable?
This is part of my result.I saw someone saying that Vina docking results with rmsd greater than 2 are unreliable, but I don’t know how to judge.Help,plz！

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

Please paste your text-based table as text into your question. Images are harder to read, reformat and interpret.


text:
Protein symmetry classification from PDB data
Given the PDB or mmCIF information of a homomeric protein assembly whose global symmetry is unknown, how might one design an algorithm to find which point group symmetry class it belongs to (cyclic $C_n$, dihedral $D_n$, or neither) without requiring manual inspection?
One idea I had for detecting cyclic symmetry was to calculate the centroid of each subunit in the homomer. If there indeed is some cyclic symmetry, the axis must pass through the global centroid, and it must be normal to the plane formed by the individual subunit centroids. To test if this axis is indeed valid, one can rotate the entire protein about the axis by 2π / n, and check that the atomic coordinates align closely with those of the original (using RMSD for example).
But I'm having trouble thinking about dihedral symmetry, specifically for simple homomers. Any help on this would be appreciated.

https://core.ac.uk/download/pdf/190349304.pdf PROSHADE for a coordinate input file , calculates the theoretical map neverthless; there should bee a reason for that approach ...

what if you have cyclic symmetry and only two subunits :  it must be normal to the plane formed by the individual subunit centroids , my geometry is really bad but how many plane passes for 2 points ?

copying from PROSHADE : Following from the definition of the chiral dihedral symmetry, a simple al-
gorithm can be devised to detect them from the already available list of detected
cyclic symmetries. Specifically, the list of the cyclic symmetries detected by the
algorithm described in the section 4.3 can be searched for any pairs of symmetries
with symmetry axes perpendicular to each other. The resulting list of symme-
try pairs then, by definition, is the list of all detected chiral dihedral symmetries
present in the structure.