text:
How to shade text in a table
I have currently made a tool that takes input sequences and maps them to a reference sequence highlighting it. You can then save the word document and it creates a table of the sequence name and the sequences.
What I want is to

Highlight the sequences in specific colours.
Highlight those sequences in the reference sequence.

Below is the create a word function.
def create_word_document(output_text, reference_sequence):
    document = Document()
    reverse_complemented_sequences = {}

    # Add reference sequence
    document.add_heading("Reference Sequence", level=1)
    document.add_paragraph(reference_sequence, style='BodyText')

    # Add a table for sequences
    table = document.add_table(rows=1, cols=2, style='Table Grid')
    table.autofit = True
    table.cell(0, 0).text = "Sequence Name"
    table.cell(0, 1).text = "Highlighted Sequence"

    # Colors for highlighting sequences
    highlight_colors = {
        "F3": "FF0000",  # Red for F3
        "B3": "0000FF",  # Blue for B3
        "F2": "FFC0CB", # Pink for F2
        "F1c": "0000FF", # Blue for F1c
        "B2": "FFFF00", # Yellow for B2
        "B1c": "00008B", # Dark blue for B1c
        "LF": "FFA500", # Orange for LF
        "LB": "800080", # Purple for LB
        # Change colors as needed
    }

    for seq_name, highlighted_sequence in output_text.highlighted_sequences.items():
        # Use the reverse complemented sequence for this section
        reversed_sequence = reverse_complement_if_needed(seq_name, highlighted_sequence)
        reverse_complemented_sequences[seq_name] = reversed_sequence

        # Add sequence information to the table
        row_cells = table.add_row().cells
        row_cells[0].text = f"{seq_name} Sequence"
        row_cells[1].text = reversed_sequence

        # Add shading to the cell containing the sequence text if available
        if seq_name in highlight_colors:
            color_code = highlight_colors[seq_name]
            apply_shading_to_cell(row_cells[1], reversed_sequence, color_code)

        # Add FIP and BIP sequences after B3
        if seq_name == "B3":
            f1c_sequence = reverse_complement_if_needed("F1c", output_text.highlighted_sequences["F1c"])
            f2_sequence = output_text.highlighted_sequences["F2"]
            bip_sequence = reverse_complement_if_needed("B1c", output_text.highlighted_sequences["B1c"])
            b2_sequence = output_text.highlighted_sequences["B2"]

            # Add a row for FIP sequence
            row_cells_fip = table.add_row().cells
            row_cells_fip[0].text = "FIP Sequence"
            row_cells_fip[1].text = f"{f1c_sequence}-{f2_sequence}"
            # Apply shading to FIP sequence if available
            if "F1c" in highlight_colors:
                apply_shading_to_cell(row_cells_fip[1], row_cells_fip[1].text, highlight_colors["F1c"])

            # Add a row for BIP sequence
            row_cells_bip = table.add_row().cells
            row_cells_bip[0].text = "BIP Sequence"
            row_cells_bip[1].text = f"{bip_sequence}-{b2_sequence}"
            # Apply shading to BIP sequence if available
            if "B1c" in highlight_colors:
                apply_shading_to_cell(row_cells_bip[1], row_cells_bip[1].text, highlight_colors["B1c"])

    # Save the document
    document.save("LAMP Primer mapping.docx")

def apply_shading_to_cell(cell, text, color_code):
    # Clear existing text
    cell.text = ''
    # Add the text to the cell
    paragraph = cell.paragraphs[0]
    run = paragraph.add_run(text)
    # Create shading element
    shading_elm = OxmlElement('w:shd')
    shading_elm.set(qn('w:val'), 'clear')
    shading_elm.set(qn('w:color'), 'auto')
    shading_elm.set(qn('w:fill'), color_code)  # Set fill color
    # Append shading element to run properties
    run._r.append(shading_elm)
    # Set font size for the entire paragraph
    paragraph.style.font.size = Pt(12)

Full code:
import tkinter as tk
from Bio.Seq import Seq
from Bio import pairwise2
from docx import Document
from docx.shared import Inches
from docx.shared import Pt
from docx.oxml.shared import OxmlElement
from docx.oxml.ns import qn

class HighlightedOutputText(tk.Text):
    def __init__(self, *args, **kwargs):
        tk.Text.__init__(self, *args, **kwargs)
        self.highlighted_sequences = {}

    def map_and_highlight(self, seq_name, reference_sequence, input_sequence):
        best_alignment = self.find_best_alignment(reference_sequence, input_sequence)

        if best_alignment is not None:
            reference_aligned, input_aligned, start_position, end_position = best_alignment

            # Configure highlight tag based on sequence name
            highlight_tag = f"{seq_name.lower()}_highlight"
            self.tag_configure(highlight_tag, background=get_highlight_color(seq_name))

            # Apply the color tag to the aligned region in the reference sequence
            self.insert(tk.END, f"\n{seq_name} Input sequence:\n{input_sequence}\n\n")
            self.insert(tk.END, "\nHighlighted sequence:\n")
            self.insert(tk.END, reference_sequence[:start_position], "default")
            self.insert(tk.END, reference_sequence[start_position:end_position], highlight_tag)
            self.insert(tk.END, reference_sequence[end_position:], "default")

            # Save highlighted sequence for later use
            self.highlighted_sequences[seq_name] = reference_sequence[start_position:end_position]
        else:
            self.insert(tk.END, f"\n{seq_name} Input sequence:\n{input_sequence}\n\n")
            self.insert(tk.END, "\nNo matching region found in the reference sequence.\n", "default")

    def find_best_alignment(self, reference_sequence, input_sequence):
        alignments_forward = pairwise2.align.localxx(reference_sequence, input_sequence)
        alignments_reverse = pairwise2.align.localxx(reference_sequence, str(Seq(input_sequence).reverse_complement()))

        if not alignments_forward and not alignments_reverse:
            print(f"No alignment found for input sequence: {input_sequence}")
            return None

        # Choose the best alignment based on the score
        best_alignment = max(alignments_forward + alignments_reverse, key=lambda a: a.score)

        # Extract aligned sequences
        reference_aligned = best_alignment[0]
        input_aligned = best_alignment[1]

        # Find the start and end positions of the input sequence in the reference
        start_position = reference_aligned.find(input_aligned.replace('-', '').lstrip())
        end_position = start_position + len(input_aligned.replace('-', '').rstrip())

        return reference_aligned, input_aligned, start_position, end_position

    def clear_output(self):
        self.delete(1.0, tk.END)
        self.highlighted_sequences = {}

def reverse_complement_if_needed(seq_name, input_sequence):
    # Check if the sequence name indicates reverse complementing is needed
    if seq_name in ["B3", "F1c", "B2", "LF"]:
        return str(Seq(input_sequence).reverse_complement())
    else:
        return input_sequence

def get_highlight_color(seq_name):
    # Return RGB values based on the sequence name
    if seq_name == "F3":
        return "#808080"  # Gray
    elif seq_name == "B3":
        return "#404040"  # Dark Gray
    elif seq_name == "LF":
        return "#008000"  # Green
    elif seq_name == "LB":
        return "#FF0000"  # Red
    elif seq_name == "F2":
        return "#40E0D0"  # Turquoise
    elif seq_name == "F1c":
        return "#FFFF00"  # Yellow
    elif seq_name == "B2":
        return "#FFC0CB"  # Pink
    elif seq_name == "B1c":
        return "#00008B"  # Dark Blue
    else:
        return "#FFA500"  # Default color for other sequences

def create_word_document(output_text, reference_sequence):
    document = Document()
    reverse_complemented_sequences = {}

    # Add reference sequence
    document.add_heading("Reference Sequence", level=1)
    document.add_paragraph(reference_sequence, style='BodyText')

    # Add a table for sequences
    table = document.add_table(rows=1, cols=2, style='Table Grid')
    table.autofit = True
    table.cell(0, 0).text = "Sequence Name"
    table.cell(0, 1).text = "Highlighted Sequence"

    # Colors for highlighting sequences
    highlight_colors = {
        "F3": "FF0000",  # Red for F3
        "B3": "0000FF",  # Blue for B3
        "F2": "FFC0CB", # Pink for F2
        "F1c": "0000FF", # Blue for F1c
        "B2": "FFFF00", # Yellow for B2
        "B1c": "00008B", # Dark blue for B1c
        "LF": "FFA500", # Orange for LF
        "LB": "800080", # Purple for LB
        # Change colors as needed
    }

    for seq_name, highlighted_sequence in output_text.highlighted_sequences.items():
        # Use the reverse complemented sequence for this section
        reversed_sequence = reverse_complement_if_needed(seq_name, highlighted_sequence)
        reverse_complemented_sequences[seq_name] = reversed_sequence

        # Add sequence information to the table
        row_cells = table.add_row().cells
        row_cells[0].text = f"{seq_name} Sequence"
        row_cells[1].text = reversed_sequence

        # Add shading to the cell containing the sequence text if available
        if seq_name in highlight_colors:
            color_code = highlight_colors[seq_name]
            apply_shading_to_cell(row_cells[1], reversed_sequence, color_code)

        # Add FIP and BIP sequences after B3
        if seq_name == "B3":
            f1c_sequence = reverse_complement_if_needed("F1c", output_text.highlighted_sequences["F1c"])
            f2_sequence = output_text.highlighted_sequences["F2"]
            bip_sequence = reverse_complement_if_needed("B1c", output_text.highlighted_sequences["B1c"])
            b2_sequence = output_text.highlighted_sequences["B2"]

            # Add a row for FIP sequence
            row_cells_fip = table.add_row().cells
            row_cells_fip[0].text = "FIP Sequence"
            row_cells_fip[1].text = f"{f1c_sequence}-{f2_sequence}"
            # Apply shading to FIP sequence if available
            if "F1c" in highlight_colors:
                apply_shading_to_cell(row_cells_fip[1], row_cells_fip[1].text, highlight_colors["F1c"])

            # Add a row for BIP sequence
            row_cells_bip = table.add_row().cells
            row_cells_bip[0].text = "BIP Sequence"
            row_cells_bip[1].text = f"{bip_sequence}-{b2_sequence}"
            # Apply shading to BIP sequence if available
            if "B1c" in highlight_colors:
                apply_shading_to_cell(row_cells_bip[1], row_cells_bip[1].text, highlight_colors["B1c"])

    # Save the document
    document.save("LAMP Primer mapping.docx")

def apply_shading_to_cell(cell, text, color_code):
    # Clear existing text
    cell.text = ''
    # Add the text to the cell
    paragraph = cell.paragraphs[0]
    run = paragraph.add_run(text)
    # Create shading element
    shading_elm = OxmlElement('w:shd')
    shading_elm.set(qn('w:val'), 'clear')
    shading_elm.set(qn('w:color'), 'auto')
    shading_elm.set(qn('w:fill'), color_code)  # Set fill color
    # Append shading element to run properties
    run._r.append(shading_elm)
    # Set font size for the entire paragraph
    paragraph.style.font.size = Pt(12)

# GUI setup
root = tk.Tk()
root.title("Sequence Alignment")

# Entry widgets for reference sequence
label_reference = tk.Label(root, text="Reference Sequence:")
entry_reference = tk.Entry(root, width=30)

# Entry widgets for multiple sequences
entries = []
output_text = HighlightedOutputText(root, height=30, width=40)
output_text.configure(font=("TkDefaultFont", 10))

for seq_name in ["F3", "B3", "F2", "F1c", "B2", "B1c", "LF", "LB"]:
    label = tk.Label(root, text=f"{seq_name} Sequence:")
    entry = tk.Entry(root, width=30)
    entries.append((seq_name, entry))
    label.grid(row=len(entries), column=0, padx=10, pady=5)
    entry.grid(row=len(entries), column=1, padx=10, pady=5)

def on_align_button_click():
    reference_sequence = entry_reference.get()

    output_text.clear_output()  # Clear previous output

    for seq_name, entry in entries:
        input_sequence = entry.get()
        input_sequence = reverse_complement_if_needed(seq_name, input_sequence)
        output_text.map_and_highlight(seq_name, reference_sequence, input_sequence)

    # Display all highlighted sequences in the output window
    for seq_name, highlighted_sequence in output_text.highlighted_sequences.items():
        output_text.insert(tk.END, f"\n{seq_name} Input sequence:\n{highlighted_sequence}\n")

# Button to trigger alignment
align_button = tk.Button(root, text="Align", command=on_align_button_click)

# Button to create a Word document
document_button = tk.Button(root, text="Create Word Document", command=lambda: create_word_document(output_text,entry_reference.get()))

# Grid layout
label_reference.grid(row=0, column=0, padx=10, pady=5)
entry_reference.grid(row=0, column=1, padx=10, pady=5)

align_button.grid(row=len(entries) + 1, column=0, columnspan=2, pady=10)
document_button.grid(row=len(entries) + 2, column=0, columnspan=2, pady=10)
output_text.grid(row=len(entries) + 3, column=0, columnspan=2, padx=10, pady=5)

root.mainloop()
```

this is a pretty good question in terms of detail and complete code etc. However i think it is more of a stackoverflow question about writing word docs in python.

I agree but I have done this (years ago) - there's validity to it being here. However, reading that amount of code for customised colouring for something I can hardly remember now isn't my idea of fun.

I dunno about Word, but you can definitely do this kind of thing in Excel via Python with openpyxl, if that helps.


text:
Looking for guidance on Microbial Association Networks (MANs) Inference with samples from different metagenomics/metatranscriptomics studies
I'm Sebastian, working on my MSc thesis, aiming to create a knowledge graph of the wastewater treatment microbiome. We're using Netcomi for constructing Microbial Association Networks (MANs) from MGnify database metagenomics and metatranscriptomics data. Specifically, we are using the CCLasso method. However, we are facing these challenges and have some questions:

Insufficient Samples: Many MGnify studies lack enough samples for MANs inference, leading us to combine data from various studies. Is combining samples from different studies a viable strategy for MANs inference in this context?
Data Heterogeneity: The data exhibits variability in sequencing technologies (e.g., Illumina Hiseq, Nextseq, Nanopore, etc), experiment types (metagenomics, assembly, or metatranscriptomics), sampling countries, and other variables. How best can we account for the data's heterogeneity in our analyses using Netcomi?

Your expertise and advice would be incredibly valuable. Thanks for your time and assistance.
Sebastian

Can you give more details on what you are looking at in MGnify? I've not used it, I use SRS directly, so why not just download the FASTQ files and process them? Also what specific question are you seeking to address, i.e. whats your hypothesis? Basically something isn't right when you say "combined data" (it's a bad idea). What do you mean by that exactly? What are you combining? To address the downvote, the question is badly explained, but it's a valid investigation and I can explain the methodology and the rationale in the answer.

Just to reiterate, I can probably answer this, but I'd need to know more detail and will be free to respond last thing Monday/early Tuesday.

source: Stack Exchange

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

sort: votes


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.

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


text:
Is there any tools to predict the potential RNA candidate that can bind to a protein?
I have a list of RNA binding proteins from HyDRa prediction. I have a long list of mRNA, I'm interested in.
For the list of RNA binding proteins, I want to find the potential RNA binding partners. Lets say each protein can bind to 5-20 RNAs, I want to find this list.
Now for each RNA-protein binding, I want to calculate the binding affinity or any equivalent measurement.
Is there any tools, I can use to do this? or any related insilico approach I can try?

Almost certainly I would look at published machine learning studies


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:
Comparing sc-seq combination and bulk-sequencing
I have expression matrices for different cell types, representing the expression of individual cells of that type. They were learned through a generative model, so I am confident they represent the approximate expression patterns of specific cell types. Now, I want to implement a bulk sequencing deconvolution using the aforementioned expression patterns. The pseudo-bulk I used is the summation of a large number of single-cell sequencing data (sc-seq). That's the background.
My first approach is to design an optimization process to optimize a series of weights, so that the product of the weights and the cell type expression approximates the pseudo-bulk. I was advised to use Poisson loss as the loss function because it aligns with the biological characteristics of RNA-seq. However, I couldn't add non-negativity constraints to the weights during optimization, resulting in negative values in the optimization results, which is meaningless. Then I found the optimiza.nnls method in the scipy package, which implements non-negativity constraints, but it uses Euclidean distance to compare the pseudo-bulk and the sc-seq combination. I obtained some good results using this method, but I have the following questions.

Can I use Euclidean distance to compare the differences between two sequencing methods? To me, this problem seems to become a linear regression problem, i.e., combining sc-seq to approximate pseudo-bulk. At this stage, there don't seem to be any biological distribution assumptions, so I guess it's feasible.
If the answer to the previous question is no, what biological assumptions does using Poisson loss follow, and what am I ignoring when using Euclidean distance for comparison?
If I want to continue using Poisson loss to optimize weights, how should I set the non-negativity constraints on the weights? I have tried methods such as ReLU and softmax in machine learning, but the results are not good.

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

I do not know the specific answers to your questions but I know there are several published papers out there about single cell deconvolution that might help.


text:
Best reference book for basic concepts when working with bioinformatics
I'm starting to study Bioinformatics, I already have basic knowledge of Biochemistry and Genetics, but I want to buy a reference book, to use when I forget a concept (I'm a Computer Science student, not a Biotechnologist).
I have already had contact with the following two books (for basic concepts):

Molecular Biology of the Cell, Alberts;
Principles of Biochemistry, Lehninger.

Which is better for my use or do you have another recommendation?
I prefer a generalist book, which covers the basic concepts, and not a Bioinformatics book, which covers basic concepts in one chapter.

check Voet's Principles of Biochemistry By Donald Voet, Charlotte W. Pratt, Judith G. Voet · 2018 .... from google  Voet & Voet is definitely the more comprehensive and it's the one I use most frequently as a reference (it's sitting above my desk in lab whereas Lehninger is back in a box in my apartment).  However, Lehninger probably does a better job of presenting the material and might be better for a novice at biochemistry looking to self-study.

and of course https://chemistry.stackexchange.com/questions/tagged/biochemistry?tab=Newest

For a reputable up to date text book on biochemistry I would only consider

'Biochemistry" Tenth Edition ©2023 by Jeremy Berg; Gregory Gatto Jr.; Justin Hines; John L. Tymoczko; Lubert Stryer published by MacMillan.

Stryer has been in the text book biochemistry business for a long time. I agree it's not a cheap option and I stopped reading biochemistry text books some time ago (see below), but this should be  as good as it gets.
Don't use Lehninger although admittedly the 7th Edition is 2017, but biochemistry is a very fast moving field and you need to be as up to date as possible if you are relying on text books as a primary source of information. To be relying on text books they need to have a "Grade A" reputation and it to be extremely to date. Few text books stand the test of time.
I've taught text book stuff from courses that relied on "out of date editions", i.e. the authors have not updated the edition for a while. It can be a nightmare because essentially you need to navigate material that you know is incorrect and insubstantial, then accommodate and moderate these deficiencies in the end of course exam. Thats okay at a certain level, at research grade it gets unstuck.


text:
missing additional visits in GSE data
When downloading GSE data, often I encounter that the authors refer to additional visits in the data description, yet, they are not found in the data.
For example, the following study. They mention additional visits, however inspecting the phenotype data extracted using getGEO function does not show them (only baseline). The same happens with different studies as well, for example this study.
Example to how I extract the phenotype data:
gsm <- getGEO('GSE88885', GSEMatrix = TRUE)
gsm <- gsm[[1]]
pheno_dat <- gsm@phenoData@data %>% as.data.frame()

Am I missing some additional data?

Can you please explain what you mean by "visits"?

If it's not in there then you can either email the authors or just live with it. It is fully on the authors what you upload to NCBI. Nobody at NCBI checks the actual content of the uploads in terms of whether it indeed contains what methods texts or the paper pretend that it contains. Nor do (usually) the reviewers. Published data are great in terms of having a large repository of resources and knowledge available, downside is that metadata are often a mess or absent. There is no general and easy answer to this.


text:
Is it normal to have low mapping to bacterial genome in a total RNA sample of plant root-colonized bacterial cells?
I had isolated the total RNA from a sample of plant roots with colonized bacteria by excising the roots (harboring the bacterial cells) and crushing them into a powder using liquid nitrogen. The RNA isolation procedure was done following the kit's instructions specific for bacterial RNA.
The RNA was sequenced and the raw reads obtained were checked for mapping percentage with the bacterial genome's GTF file using bowtie2, as I need to study the expression of the bacterial genes. While the mapping was close to 80% for samples coming from the planktonic state, it was quite less - 0.4 to 0.56%, for all the samples coming from root colonized cells, obviously because a major part of it is comprised of plant RNA. Firstly, is this acceptable? Can I use these reads for further analysis?
Secondly, I use DESeq2 to analyze the data for upregulated and downregulated genes, comparing the colonized cells to the planktonic ones. Since the package includes a normalization step, is it safe to assume that although the mapping is less it does not influence the analysis significantly?


text:
Practical challenges for microarray analysis and single cell
currently, I am learning single-cell sequencing and microarray analysis. I am searching for somewhere (or any alternative way)to practice and get hands-on with these techniques. please help me with that.
Thank you in advance.

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

Please be more specific on which level of analysis you are planning to do. Are you looking for point-and-click graphical software or do you plan to construct your own workflows, and if yes, in which language?


text:
How to design a synthetic wastewater FASTQ file that could be flagged as "engineered"?
I am an apprentice-level python user who occasionally works in medical laboratories and mainly does business work. I am presented with the unenviable task of orchestrating a proficiency test of some sequencing teams who have the almost equally unenviable task of evaluating metagenomic FASTQ files of wastewater (sewage) data where the question is "could this have been engineered ?"
Trick: how do you make some FASTQ files that look like wastewater (plenty of those online), definitely has some pathogen (SARS-CoV-2, Ebola, whatever), oh, and the pathogen has to look like it could have been engineered? And make an equal number that do not have the pathogen but otherwise look the same?
Here is my initial thoughts:

World's simplest version: cat a virus FASTQ onto the end of a
wastewater FASTQ.

change all the identifiers in the virus file so
they look like a continuation of the wastewater, chop the virus file
with visual mode in Vim so the reads are the same length as the
wastewater.

instead of cat, do some python to build a new
FASTQ by randomly selecting lines from the virus and wastewater
"feedstock" FASTQs.

Very fancy: grab a assembled virus FASTA, splice in an interesting
bit of pangolin, raccoon, pig, etc, then randomly select and subtly
mutate substrings as reads of appropriate length, and assemble a de
novo synthetic virus FASTQ. Mix that into the wastewater FASTQ as
above.

Surely I am missing many, many steps and "very fancy" is probably not that fancy (or maybe the sun will grow cold before it succeeds, I do not know). Help and suggestions welcome, it could be fun.

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 pathogen has to look like it could have been engineered  ---> after CRISP editing it could be very difficult to prove that

The question needs breaking up into separate questions. I would suggest 3 questions. No-one is going to spend time considering an entire project design in a single question. The whole issue about "is engineered" needs explaining. In reality the project looks relatively easy - but again you are posting an entire project as a single question.

The law considers genetically engineered different from genome edited ... please take that onyo your question

Use existing software to generate simulated fastq files. The input for these tools is generally a set of reference organisms that you'd expect to find in your sample. Simply add your 'engineered' viral genome to the list of references. These tools have error models to simulate all the errors that real sequencing instruments make.
How you'd go about making a viral sequence recognisable as genetically engineered, I have no clue. I'm sure the person who gave you this assignment has an idea on how to engineer viruses and what recognisable features that would result in, maybe consult with them.


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.


text:
How to download the same 18S rRNA gene region for multiple species?
I need 18S rRNA gene sequences for a wide variety of eukaryotic species.
I know that Ensembl has this gene sequence, but I don't understand what is the correct protocol to download the sequences. I already know  the basics of Ensembl's Rest API, but there are still a few things that I don't understand.
I only know the gene name (18S rRNA) and the list of taxa I want to select it from. But I do not know the region of the gene.
For example: I tried looking for '18srrna' on the Ensembl browser and I noticed that for the Fruit fly there are 2 sequences corresponding to this gene (Ensembl IDs: FBgn0085813 and FBgn0267521).

How do I know which sequence I should use to make an alignment with other species?
How can I batch download the same region for all the species of interest?
Do you know any protocol/guide that goes into detail of how to download and select sequences to make phylogenies?

I don't have a solution for question (1) but for question (2) I am planning on requesting from the API the output of the endpoint GET sequence/region/:species/:region. But then again the problem is that I would need the position of the region for each species, and it sounds unfeasible to do for more than a dozen of species.

Hi @Mirk. I'll answer over the weekend if no-one else does. You wrote, "eucariotic" do you mean eukaryotes?

Yes, sorry I did a spelling mistake but now I corrected it. Thanks for pointing it out.

How do I know which sequence I should use to make an alignment with other species?

Use both. For fruit flies there should thousands of 18S sequences.

How can I batch download the same region for all the species of interest?

Using Blast (below)

Do you know any protocol/guide that goes into detail of how to download and select sequences to make phylogenies?

Sure - Blast ... align via muscle -super5 and build via IQ-TREE
EMBL I don't know - using REST API via a CLI is possible. The problem is that you could potentially 1000s of 18S genes but a simple search term can ignore everything that isn't part of a genome explicitly labelled "18S". Thats where a Blast style approach is useful.
However, for NCBI Blast they intentionally block access to CLI with >100 hits requested. This is a complete pain and really, really limits using Blast to extract out a gene remotely using a remote request script.
Basically the way to do it is to ditch the API and run this via the Blast NCBI web interface, or download the entire eukaryotic database and perform local Blast, MMseqs2 or Minimap2.
It's a bit bizarre I know, NCBI just doesn't like bots running its Blast service, but it gives more permission to use its web-service. The central advantage is that Blast will nicely remove the 18S gene from any given genome which can be downloaded as a single locus, without needing to download loads of genome stuff thats not needed. NCBI still does impose a limit on how much can be blasted via the web which sadly now is only 5000 (I'm sure it was much more like 50000) and its really not cool.
What you can do is use the Blast search "exclude/include" parameters to target the taxa of interest within the 5000 limit and using a tree of relevant eukaryotic evolution navigate the data sets. I assume you are looking for around 100000 sequences. Basically, once you've blasted you have the option of downloading just the region blasted and nothing else ... or the whole genome (avoid this option)... you will then have 5000 nicely 18S genes excised in silico
You can also set a lot other blast parameters, in particular only considering full length hits.

I can suggest two ways to batch-download aligned sequences for the species of interest. Both ways are integral to the Ensembl database and toolsets.
But before that, I suggest FBgn0085802 (which is rRNA) instead of FBgn0085813 and FBgn0267521 (which are pseudogenes). You may check FlyBase database.

Since you are interested in Fruit fly 18S rRNA homologs, Metazoa Ensembl is preferable. You query Ensembl ID in the Ensembl search box (at www.metazoa.ensebml.org) and click the popping link. 18SrRNA:CR41548 comes for FBgn0085802. On the left-hand side, Gene-based displays, select Sequence. under the Marked-up sequence, pick "BLAST this sequence".

Opening the BLAST/BLAT search automatically sets a DNA alignment query. You may search against the species of interest and add "18S" as Description. Later, run without changing any other parameter. It may take several minutes. Then, view results for hits found. You can investigate Sequence/s with the best alignment and closest genomic locations, then export data in FASTA or other formats.
While this assures you of the alignment and lets you investigate in detail the homologs, you may not batch-download particularly in this way.

As mentioned, BioMart Restful API is a way to batch-download sequences that you want. The website does only allow one search at a moment but editing the Perl script, which is generated according to your query, you may get sequences for multiple datasets. For setup, you should follow BioMart Perl API guideline cautiosly.

For your example query, set proper datasets, filters and attributes. Again for fruit fly, you may refer to:
https://www.metazoa.ensembl.org/info/data/biomart/biomart_perl_api.html
My example query seems like below in popping Perl script:
$query->setDataset("dyakuba_eg_gene");
$query->addFilter("with_dmelanogaster_eg_homolog", ["Only"]);
$query->addFilter("biotype", ["rRNA"]);
$query->addFilter("transcript_biotype", ["rRNA"]);
$query->addAttribute("ensembl_gene_id");
$query->addAttribute("ensembl_transcript_id");
$query->addAttribute("gene_biotype");
$query->addAttribute("description");
$query->addAttribute("transcript_biotype");


text:
Abnormal representation of differentially expressed genes in volcano plots located at the extremities
This question was also asked on reddit

Does anyone know how to correct the abnormal genes that seemed to be located on the curves near the extremities?
The analysis was done using R package edgeR on samples measured once. However, we couldn't figure out what was causing the anomaly or how to fix it.
Does anyone have any advice on this matter?
The code used for analysis is from the following website:  https://blog.csdn.net/u012110870/article/details/102804557

On samples measured once as in no replicates?

@RamRS Yes. I know that this makes the data unreliable, but long story short, I am trying to make the best of this situation.

In "normal" (that is replicated) experiments such structures usually are due to poor/insufficient/no prefiltering when there are unreliable count outliers (like two samples have zero counts for a gene and one has a lot).
As you cannot really meaningfully do any prefiltering other than requiring that either one or both samples have a certain read count you should probably just see these "DE results" (it's not differential really as it is unreplicated) as a ranking. Fold changes and p-values are very unreliable at 1 vs 1 so just rank by something like signed p-value and go along with that. Don't overinterpret these sorts of plots which are mostly meaningless in such a setup.
That having said, for prefilter I would probably require that at least one sample has a good number of reads, like 50 to at least dampen unreliably large fold changes a little bit in this setup.

Duly noted!
As for the filter cutoff, how do you choose the number? Is there someway of determining this or does this come with experience and intuition? Also if I have multiple comparisons with the same cell type (just different culturing conditions) would you keep this cutoff constant or does it depended on the data?

Might as well transfer my reddit answer here...
If your p-value is smaller than the reciprocal of the number of atoms in the universe, the model used for calculating that differential expression is almost certainly wrong.
Do an MA plot; that will give you a better idea of what's going on.

X axis - mean of (log expression for all samples); should be in the range of 0~15
Y axis - log2 fold change
colour - p-value (linear scale)

I suspect that the curves you're seeing are likely to be low count genes, in other words, genes on the left hand side of the graph in an MA plot.
I also suspect, given the super-astronomical p-values, that you've got an offset L2FC (i.e. the "spine" of the MA plot doesn't sit at L2FC = 0), suggesting that a covariate hasn't been appropriately accounted for. This is common when processing data from multiple samples from multiple different methods.


text:
Method for reproducing minimal, cartoon-style protein structure visualisation
I recently came across a preprint ( DOI, PDF link ) that has some figures of protein structures in a style that I have not seen before.  To me it appears easier on the eyes and less noisy than the more traditional 3D models with lighting that "reflects" on a glossy surface.  I am interested in creating similar plots of my own protein structures.
Does anyone recognise this style / know how to generate these plots from a PDB file or similar?
(Extended Data Figure 2, page 20)

see https://ccsb.scripps.edu/illustrate/

Thanks. Any chance you could put the link with a small explanation as an answer @pippo1980?

https://ccsb.scripps.edu/illustrate/documentation/

https://github.com/ccsb-scripps/Illustrate

I've never used it, so cant really help you out , sorry


text:
How to use CIRCOS plot pairwise syteny blocks for multi species
Hi, I want to use Circos to draw syntenic relationships among multiple species, like the Fig.1b in this article. How to define the ideograms for different species?
Thank you.
https://doi.org/10.1038/s41588-024-01683-

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:
How to create a dataset of contact maps of homologous proteins?
Disclaimer: I am a machine learning researcher working on network science and want to use proteins to test an algorithm of mine using real-world data. My bioinformatics knowledge is minimal.
I need to create a dataset of contact maps. Namely, for a given set of proteins, I want to find the set of binary matrices that represent each of their contact maps. I understand that this binary matrix is obtained by thresholding the distances between atoms with some criteria. Thus, it would also be helpful to find the set of distance matrices for these proteins (I can just threshold them afterward).
For computational reasons, I would prefer proteins that are as small as possible. Namely, proteins that have the fewest amount of amino acids. Also, I need the proteins to be homologous since I want to be able to do a Multiple Sequence Alignment (MSA) to them and analyze the matrix returned by the MSA with my machine learning algorithm.
So my questions is: given a reference protein (let's say YES_HUMAN for example), where can I find a list of homologous proteins and the contact maps for all of them?

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

Finding homologs is easy enough, there are many databases of homologous sequences, you can get this from Ensembl or similar places. I would instead focus this question on getting the contact maps since that is something I personally have no idea where to find. Once you've figured out how to get a set of contact maps, then selecting the homologous sequences won't be hard.

@terdon As far as I understand, contact maps can be computed (or estimated) by thresholding the distances between the amino acids that form the sequence. Even if you don't know where to obtain the contact maps directly, would you know where I can find a dataset of proteins with the positions of the amino acids (i.e., of their $C_\alpha$ atoms)?

@Tendero The Protein Data Bank has both experimentally determined, and computationally predicted 3D models of proteins https://www.rcsb.org/

according to wikipedia The average size of a protein increases from Archaea to Bacteria to Eukaryote (283, 311, 438 residues and 31, 34, 49 kDa respectively so your yes_human with 543 is not that small , or is it ?

try https://ngdc.cncb.ac.cn/hgd/gene and browse for yes_human

https://ngdc.cncb.ac.cn/hgd/gene?type=&keyword=P07947%20

@pippo1980 In Table 1 from this paper it says that YES_HUMAN has a length of 48. What does that mean then?

https://www.uniprot.org/uniprotkb/P07947/entry#sequences

@Tendero Crystallographic structure of the SH3 domain of the human c-Yes tyrosine kinase: Loop flexibility and amyloid aggregation. https://www.rcsb.org/structure/2HDA

https://en.m.wikipedia.org/wiki/SH3_domain


text:
How to add corresponding "chr" before location (start, end) of bed alike file?
I have some bed alike files but don't have chr term in the fisrt column, that I want to add chr and its corresponding chromosome number, I'm wondering how can I do that in bash!
head dat.tsv
9901    10000   0.186130156985236
10001   10100   244.668091357092
10101   10200   367.420929888855
10201   10300   179.801731647738
10301   10400   141.551984387272
10401   10500   92.5066880216621
10501   10600   92.3205578646769
10601   10700   127.96448292735
10701   10800   74.6381929510795
10801   10900   62.6327978255318

Thanks for including an input file. Please show what you would like the output to look like after the chr term has been correctly added.

I think the biggest problem is that it's not clear how you/we would know what chromosome numbers go with each row.

You could use the following perl one-liner, however if genome you are using contains more than one chromosome you need to figure out which lines get which chromosome chr1, chr2. If it's single chrosome like E.coli, usually name of chromosome is "chr".
perl -npe "s/(.*)/chr\t\1/" dat.tsv

I appreciate your responses.
Using the R dplyr package after using bedtools, I was able to deal with this problem.
Using the makewindows function from bedtool, I first obtain a 100 bp window.
The full_join function from the dplyr R package then allows me to obtain a satisfactory result.
I was looking for a script that only uses bash.


text:
How should multi allelic sites be handled in calculation per site Fst?
I'd like to calculate per site Fst(fixation index) to investigate genomic differentiation between two populations(haploid).
I did SNP calling and got vcf file.
In the vcf file, multi allelic sites exist:
#CHROM  POS ID  REF ALT
1   154 .   G   "A, T"

For example this site is triallelic sites.

Could I calculate Fst in this site by normalizing vcf file?

Normalization means dividing into 2 lows.
REF G ALT A(T would be convert G)
REF G ALT T(A would be convert G)
In calculation using Popgenome package, we can specify missing sites(Missing sites would not be included in the calculation).
In this case,　which is better: replacing with 'REF' or replacing with missing sites?

Why most studies calculate only biallelic sites?

Before analyzing, removing variants based on MAF is make sence?
For example, one sites include A(56%), T(44.5%), C(0.5%).
In this case, could I remove variant C(replace missing site(NA)) before calculation of Fst?

Best regards,

Please copy-paste tables of text as text into your question posts. This makes them easier to read on some devices, and makes it more easy to modify the tables to improve the question.

Thank you for your advice.

Hmmmm .... I personally would not use a single VCF to calculate Fst in the way you are describing. It will also miss one of the nucleotides. I would aggregate across VCFs for a homologous position and use those frequencies. Its designed for heterozygotes, you are describing a haploid genome, like mtDNA (I think). Fst needs a genuine allelic frequency, whereas you've shoehorned it (possibly).
Bi-allelic because that is a synonymous variation for most amino acids, i.e. without triggering an amino acid change.
In terms of removing a variant ... unless you believe this represents the in vitro error rates (which you can define from an invariant site) ... then okay. If it is above the in vitro error rate then you can't remove it. I know the read-depth here, so I can't comment.

Just to clarify Fst can be used for prokaryotic genomes (haploids), but unless I've misunderstood the question it seems to being derived from a single VCF ... single Metagenome (?). Thats not how it works. I can describe the what the calculation is measuring will artificially inflate Fst and you can't correct for it.


text:
WGCNA soft threshold selection
Here I am conducting WGCNA analysis. When calculating the soft threshold, it is basically below 0.8. I choose Power=6, which only has 0 and 7. Is it feasible?

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 copy-paste text as text into question posts. This makes the text easier to read on some devices, and allows it to be more easily modified when required.

I don't know the context (how many samples, what tissue, experimental design etc) but the mean and median connectivities are pretty high at power 6. I would first check the data for strong drivers of correlation and possibly adjust for them if they do not seem biologically relevant to the question you study. Assuming you find no issues, I would go with a higher power here, maybe 12.


text:
How can I filter variants in a VCF file based on size?
I have a filtered VCF file with indels (no SNPs). Does anyone know of a easy way to use bcftools/vcftools to filter out small or large indels (I'm interested in indels between 50-100bp in size)? 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.

What is the output to bcftools view -H <file.vcf> | head -n 5 where <file.vcf> is the actual VCF file?

Vcftools allows you to filter to include or not include a list of sites.
https://vcftools.github.io/man_latest.html#SITE%20FILTERING%20OPTIONS
If you can first grab the locations of indels of this size (for example you could write a python script  that goes through your vcf and grabs the location of anything with ref or alt of a length of at least 50 characters), then you can provide vcftools with a list of variants to filter on.
If you are not comfortable writing your own script to look through the file, you could delete the vcf header lines, change the file extension to tsv, then open it in excel. Then you could use excel to filter for those sites and get their positions.


text:
How can I check if a PDB file is incomplete?
One of the requirements of simulating the protein chain is to check for missing atoms or residues. I.e., using software tools to identify (and, if necessary, model) any missing regions.
How can I check if a PDB file is complete or incomplete using Python?

Please [edit] your question and add more detail. Ideally, give us an incomplete file and a complete file so we can test any solutions. I know next to nothing about PDB files, so if the information is not in the file itself, you will also need to provide a way to compare it to something, so we would need to know what proteins you are aligning and from what species, so we can try and find the full sequences to check for completion.

OK, lets say you need that information for 2PE4.pdb, if you look at its webpage on RCSB PDB repository: 2PE4 and display its .pdb format file at https://files.rcsb.org/view/2PE4.pdb, you'll notice at REMARK 465:
REMARK 465
REMARK 465 MISSING RESIDUES
REMARK 465 THE FOLLOWING RESIDUES WERE NOT LOCATED IN THE
REMARK 465 EXPERIMENT. (M=MODEL NUMBER; RES=RESIDUE NAME; C=CHAIN
REMARK 465 IDENTIFIER; SSSEQ=SEQUENCE NUMBER; I=INSERTION CODE.)
REMARK 465
REMARK 465   M RES C SSSEQI
REMARK 465     GLY A   437
REMARK 465     HIS A   438
REMARK 465     HIS A   439
REMARK 465     HIS A   440
REMARK 465     HIS A   441
REMARK 465     HIS A   442
REMARK 465     HIS A   443

this should be the missing residues from the full peptide sequence of the protein that is listed as SEQRES:
SEQRES   1 A  424  ARG SER PHE ARG GLY PRO LEU LEU PRO ASN ARG PRO PHE
SEQRES   2 A  424  THR THR VAL TRP ASN ALA ASN THR GLN TRP CYS LEU GLU
SEQRES   3 A  424  ARG HIS GLY VAL ASP VAL ASP VAL SER VAL PHE ASP VAL
SEQRES   4 A  424  VAL ALA ASN PRO GLY GLN THR PHE ARG GLY PRO ASP MET
SEQRES   5 A  424  THR ILE PHE TYR SER SER GLN LEU GLY THR TYR PRO TYR
SEQRES   6 A  424  TYR THR PRO THR GLY GLU PRO VAL PHE GLY GLY LEU PRO
SEQRES   7 A  424  GLN ASN ALA SER LEU ILE ALA HIS LEU ALA ARG THR PHE
SEQRES   8 A  424  GLN ASP ILE LEU ALA ALA ILE PRO ALA PRO ASP PHE SER
SEQRES   9 A  424  GLY LEU ALA VAL ILE ASP TRP GLU ALA TRP ARG PRO ARG
SEQRES  10 A  424  TRP ALA PHE ASN TRP ASP THR LYS ASP ILE TYR ARG GLN
SEQRES  11 A  424  ARG SER ARG ALA LEU VAL GLN ALA GLN HIS PRO ASP TRP
SEQRES  12 A  424  PRO ALA PRO GLN VAL GLU ALA VAL ALA GLN ASP GLN PHE
SEQRES  13 A  424  GLN GLY ALA ALA ARG ALA TRP MET ALA GLY THR LEU GLN
SEQRES  14 A  424  LEU GLY ARG ALA LEU ARG PRO ARG GLY LEU TRP GLY PHE
SEQRES  15 A  424  TYR GLY PHE PRO ASP CYS TYR ASN TYR ASP PHE LEU SER
SEQRES  16 A  424  PRO ASN TYR THR GLY GLN CYS PRO SER GLY ILE ARG ALA
SEQRES  17 A  424  GLN ASN ASP GLN LEU GLY TRP LEU TRP GLY GLN SER ARG
SEQRES  18 A  424  ALA LEU TYR PRO SER ILE TYR MET PRO ALA VAL LEU GLU
SEQRES  19 A  424  GLY THR GLY LYS SER GLN MET TYR VAL GLN HIS ARG VAL
SEQRES  20 A  424  ALA GLU ALA PHE ARG VAL ALA VAL ALA ALA GLY ASP PRO
SEQRES  21 A  424  ASN LEU PRO VAL LEU PRO TYR VAL GLN ILE PHE TYR ASP
SEQRES  22 A  424  THR THR ASN HIS PHE LEU PRO LEU ASP GLU LEU GLU HIS
SEQRES  23 A  424  SER LEU GLY GLU SER ALA ALA GLN GLY ALA ALA GLY VAL
SEQRES  24 A  424  VAL LEU TRP VAL SER TRP GLU ASN THR ARG THR LYS GLU
SEQRES  25 A  424  SER CYS GLN ALA ILE LYS GLU TYR MET ASP THR THR LEU
SEQRES  26 A  424  GLY PRO PHE ILE LEU ASN VAL THR SER GLY ALA LEU LEU
SEQRES  27 A  424  CYS SER GLN ALA LEU CYS SER GLY HIS GLY ARG CYS VAL
SEQRES  28 A  424  ARG ARG THR SER HIS PRO LYS ALA LEU LEU LEU LEU ASN
SEQRES  29 A  424  PRO ALA SER PHE SER ILE GLN LEU THR PRO GLY GLY GLY
SEQRES  30 A  424  PRO LEU SER LEU ARG GLY ALA LEU SER LEU GLU ASP GLN
SEQRES  31 A  424  ALA GLN MET ALA VAL GLU PHE LYS CYS ARG CYS TYR PRO
SEQRES  32 A  424  GLY TRP GLN ALA PRO TRP CYS GLU ARG LYS SER MET TRP
SEQRES  33 A  424  THR GLY HIS HIS HIS HIS HIS HIS

Pay attention that the .pdb starts at residue position 20; that is a cloning artifact, as shown by:
SEQADV 2PE4 ARG A   20  UNP  Q12794              CLONING ARTIFACT
SEQADV 2PE4 SER A   21  UNP  Q12794              CLONING ARTIFACT
SEQADV 2PE4 THR A  436  UNP  Q12794              CLONING ARTIFACT
SEQADV 2PE4 GLY A  437  UNP  Q12794              CLONING ARTIFACT
SEQADV 2PE4 HIS A  438  UNP  Q12794              EXPRESSION TAG
SEQADV 2PE4 HIS A  439  UNP  Q12794              EXPRESSION TAG
SEQADV 2PE4 HIS A  440  UNP  Q12794              EXPRESSION TAG
SEQADV 2PE4 HIS A  441  UNP  Q12794              EXPRESSION TAG
SEQADV 2PE4 HIS A  442  UNP  Q12794              EXPRESSION TAG
SEQADV 2PE4 HIS A  443  UNP  Q12794              EXPRESSION TAG

See PDB residue numbering
Using Biopython lib you could try using 2pe4.cif instead of .pdb file (https://files.rcsb.org/view/2PE4.cif) as input :
# First import MMCIF2Dict function
from Bio.PDB.MMCIF2Dict import MMCIF2Dict

# read in the CIF file.
pdb_info = MMCIF2Dict("2pe4.cif")

#
print(pdb_info["_pdbx_unobs_or_zero_occ_residues.label_seq_id"])
print(pdb_info["_pdbx_unobs_or_zero_occ_residues.label_comp_id"])
print(list(zip(pdb_info["_pdbx_unobs_or_zero_occ_residues.label_seq_id"],
pdb_info["_pdbx_unobs_or_zero_occ_residues.label_comp_id"])))

Result:
['418', '419', '420', '421', '422', '423', '424']
['GLY', 'HIS', 'HIS', 'HIS', 'HIS', 'HIS', 'HIS']
[('418', 'GLY'), ('419', 'HIS'), ('420', 'HIS'), ('421', 'HIS'), ('422', 'HIS'), ('423', 'HIS'), ('424', 'HIS')]

but for this kind of stuff it's probably better suited
PDBx Python Parser


text:
Any other (good) long read aligners apart kma?
I wonder if there are any other more efficient long read aligners apart kma? I have heard there is some hisat2 long read application but it's not really efficient. Any suggestions?

"Apps"? You mean aligners?

Yes, aligners. I shouldn't call them apps, I just meant some module that can be attached to a script, like bwa, minimap or that kma for ex.

"module that can be attached to a script" is not a thing. Modulefiles are a way of packaging binaries and related dependencies, so you're looking for an aligner like ATPoint says, not an app or a module. If you need a binary as a module, consult your HPC sysadmin.

Please [edit] your question and clarify your needs. If you need an aligner, then please say so (apps is meanigless). If you instead mean a variant caller, say that. Then, explain what more efficient means for you (faster? better resource utilization? less RAM specifically? something else?), and whether you are tied to a specific operating system. Perhaps also tell us what kind of sequencing you will be doing (targeted? WES? WGS?).


text:
how to annotate the effect of a genomic variant to predict what the resulting protein is?
Spoiler: This is probably a theoretical question. But, of course, if you have any tool that does this I will be extremely happy.
So, I have WGS data for an organism O, both wild type (WT-O) and modified (M-O). I was able to assemble the reads using the O reference genome. Then, I was able to compare WT-O against M-O. I finally have the variants files (vcf), which I annotated with SnpEff and VEP.
Now what I need to do is to evaluate the effect of the variants at protein level. What I know from the vcf files is which nucleotides were affected, their position, which chromosome, gene, and whatnot (in most cases, although this is not always true for all the variants).
Question is: how to translate these variants from nucleotides to amino acids (AAs) and then, eventually, compare (BLAST?) the AAs sequences containing the variants to check whether they match with any known protein?
EDIT:
so the final goal would be to know whether the variants are causing some weird protein to be produced which may have a deleterious effect. It can be something like a toxin, or an allergen that are normally not present but that can appear "by random chance". I need to make sure that nothing like that happened and to be sure i need to know which proteins may result from the variants.

Doesn't VEP provide the effect at the protein level? At the very least it should give the reference and modified amino acids, right? Can you please show us a few lines of your VCF so we know what you have to work with? Also, I assume you don't have a budget for paid tools here, right? The company I work for has an API that can translate variants to HGVS protein notation and from that you can write a script that can build the protein sequence. But it's a paid service and you would still need to write your own script to get the full protein sequence.

Can you also give us some background on the biological question you are trying to answer? Why would you want to compare the proteins? It seem like a strange thing to do, I can't really think of why that would be useful. What is the final objective here?
What information do you want to get from this comparison? And what is this organism, does it have splicing or could you just translate the DNA directly?

@terdon, thank your for the link. i updated the question and i hope the final goal is more clear, now. also...i am not sure i can provide you with a sample from the VEP annotated file but i'll look into it.

@terdon, maybe i have what i need in this VEP file, but it is not very clear to me how it works. i am quite new to this field of variant calling. i am looking into it.

Just open the vcf file you have annotated with VEP, and show us a couple of variant lines with the annotations. Feel free to change the position/chromosome if you fear you cannot share it (although a couple of variants cannot be sensitive information). I want to see what kind of protein-level information you have. But the rest of your experiment doesn't really make sense to me. It sounds like you want to investigate the pathogenicity of the variants. There is no reason to try and re-create the protein and search for it, especially since the variants will often only change one amino acid.


text:
multiplatt from indicspecies package not producing results for all groups
Based on presence absence data, I'm attempting to find indicator species for different sample sites. I have 25 different groups or "clusters" as I've labelled them. Running the R code below, only results for 13 of the 25 groups was returned.
Is the output truncated, i.e. is there a max output parameter in R I need to fix? Or does this imply that there simply aren't any statistically relevant indicator species in the groups not included in the output?
Code used:
library(indicspecies)

data = read.csv("ISAI.txt", header= TRUE, sep= "\t") #Load up data from file

abund = data[,3:ncol(data)] #Collect only numerical data, starts in 3rd column
cluster = data$Cluster #Groups are defined in the Cluster column of the imported data

inv = multipatt(abund, cluster, func = "r.g", control = how(nperm=1000),duleg = TRUE) #Run analysis.

summary(inv) #print out results

Notes:
For the func parameter, "r" was used due to data being in presence/absence format. The "g" corrects for groups with different numbers of members. The duleg = True prevents group site combinations from being assessed as I'm not interested in those.

Got it! :)
I switched the summary(inv) to summary(inv, alpha=1) to get all possible results and all my groups showed up.
So, the answer to my own question is that the groups are missing from the output because they don't have any statistically significant OTUs correlated with those groups.


text:
Gene set enrichment ontology analysis issue
I have been dealing with this problem for the past 48 hrs. I do not know what is it, which prevents the function from assessing the parameters it gets passed.
params <- GSEAGOHyperGParams(name = "My Custom GSEA based annot Params",
                              geneSetCollection = new_df,
                              geneIds = new_df$GOIds,
                              universeGeneIds = gene_universe,
                              ontology = "BP",
                              pvalueCutoff = 0.05,
                              conditional = FALSE,
                              testDirection = "over",
                              organism = "human")
Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘geneIdType’ for signature ‘"character"’

In an example below (edited), you could see that it is character, but either way if I change it do integer, it does not work. Maybe character vector ?
> str(new_df)
'data.frame':   68 obs. of  4 variables:
 $ GOIds   : chr  "GO:0001725" "GO:0001725" "GO:0003779" "GO:0005515" ...
 $ evi1    : chr  "IBA" "IDA" "IEA" "IPI" ...
 $ genes1  : chr  "23002" "23002" "23002" "23002" ...
 $ ontology: chr  "CC" "CC" "MF" "MF" ...

Same goes, when I try to cast it into an integer type by using as.integer().
Any suggestions for this function from an R package clusterProfiler?
I believe this is from this package, but there could be one or two other packages, where this function maybe present.

What does ensembl_gene_ids have to do with your function call?

Hi! What package are you using in R? Thanks

Edited the query and answered to both ! Thanks !

According to the documentation, geneSetCollection should be a GeneSetCollection Object : https://www.rdocumentation.org/packages/Category/versions/2.38.0/topics/GSEAGOHyperGParams

Okay, thank you every1. With a duplicate answer from Stack Overflow, I was able to succeed a bit only to stumble upon a different issue, for some reason the function OverMF <- hyperGTest(params2). Error in getUniverseHelper(probes, datPkg, entrezIds) :    After filtering, there are no valid IDs that can be used as the Gene universe.   Check input values to confirm they are the same type as the central ID used by your annotation package.   For chip packages, this will still mean the central GENE identifier used by the package (NOT the probe IDs).


text:
How to do post hoc comparisons after a repeated measures ANOVA in R
I have a data set of several samples with their expression of proteins, in response to four different doses of a drug and two genotypes. I have been able to generate the two-way ANOVA using ezANOVA but I am having difficulty in generating multiple comparisons results. In GraphPad, I am able to compare means of each group to the control, the 0 dose, within each genotype and apply a Dunnett's correction. However, I am trying to implement this in R and my significant values are not lining up with GraphPad. I am using, emmeans and lme4 package to try generate these multiple comparisons. Any ideas or other readings to try? Thanks.
The code for how I'm generating the anova results below
library(ez)
ezANOVA(data = data_long,
                         dv = .(Expression),
                         wid = .(Subject),
                         within = .(Dose),
                         between = .(Genotype),
                         type = 3,
                         detailed = TRUE,
                         return_aov = TRUE)

Then the code I used to generate the Dunnett's Test
library(DescTools)

DunnettTest(data_long$Expression, data_long$Genotype)

However the output of the Dunnett's test only compares from the Genotype and not the dosages within each Genotype at each dose.
The output in R

The output in GraphPad

I've tried switching from Genotype to Dose, but it doesn't separate the pairs by Genotype.


text:
Fisher's exact test for COG categories in pan, core genome analysis
I have come across several papers that do a Fisher's exact test to show over represented genes in COG categories specifically in a pan-core-accessory-species specific genome analysis. The test is done to show over represented genes in the core, accessory and species specific datasets as compared to the pan genome. Here is one such paper - Vibrio cincinnatiensis - comparative genome analysis
Usually, the way the pan and core genome analysis is represented in papers is by showing the proportion of genes in each of the COG categories in both. For example, the proportion of C (Energy production and conversion) could be 6.07 in the pan but 6.8 in the core genome. The Fisher's exact test is used to show that this particular category is over represented in the core genome as compared to the pan.
I already have the csv files from eggNOG from which I filtered out the genes in major COG categories for the pan, core, accessory and species-specific genomes into a single Excel sheet for 8 closely related genomes. Could someone tell me how to perform this test on this dataset?
My null hypothesis is that the proportion of genes in the pan, core, accessory and species specific genomes is the same. But, obviously this analysis is done to show how the proportions vary in some COG categories between the above mentioned 4 genome datasets.
I did try using python and R as well, but the test requires 2 rows and columns and here I am just comparing two columns (pan, core) with respect to 1 row - the 'C' category. It seems impossible to apply the test to this problem and I am sure I am missing something and that's why I posted it here.
I tried doing the Fisher's test on the proportions I mentioned previously (6.07 and 6.8) but I am getting an error : 'Fisher's test requires two rows and two columns'. Should I compare the absolute values, like '3745' genes in the pan and '4195' genes in the core?

Please [edit] your question and clarify what you want to test.

Sorry I've got my head around your question. For Fisher's Exact you need 2 "independent" observations and for SciPy its a 2x2 table.
So in your case ...

Genome
Species A
Species B

Pan
6.07
6.10

Core
6.08
6.12

If you think about it, thats a better way of doing it because thats comparing pan vs core between species. A single observation isn't sufficient. To really do this you need 2 x 10 type contingency tables but SciPy will not do that, but you can do it in RStudio. For example using 10 species rather than just 2 ... and thats a real test of pan/core. My advice is to have your test species and in the other column have your control species ... like E. coli that sort of thing.
The SciPy solution is here. Place your data in a list of lists here its mylist. Next place into a numpy array and place that in a scipy function, like this...
import numpy as np
from scipy.stats import fisher_exact
mylist = [[6.07, 6.08], [6.10, 6.12]]
mytable = np.array(mylist)
res = fisher_exact(mytable, alternative='two-sided')
print (res.pvalue)

It should print the p value where > 0.05 is not significant.
To start doing 2 x >2 contingency tables in Python r2py is needed and from my experience (albeit about 5 years ago) I wouldn't advise it. r2py is a disaster when I tried to use it and just did the calculations in RStudio.

Hi @M__ thank you so much for your input. Will try doing it with Rstudio.

Thanks @K_081 it should work fine.

Hi @M__ , I have added a table to my question for clarity, could you possibly go through that? I want to compare the proportion of "C" or other COG categories between the pan, core and the others, not between species as all the species I need for the analysis are included in the pan, core. Sorry for all the trouble and thanks a lot, if you can look at this again.

@K_081 could you kindly ask this as a separate question please and link this question within it? It doesn't look easy but I know what COGs are.

Done @M__ : https://bioinformatics.stackexchange.com/questions/22210/calculating-fishers-exact-test-for-cog-categories


text:
Gene ensemble vs Hugo id
Hi
I have been recently struggling with the gene annotations and changing gene ID to ensemble ID. While I was doing that there were some genes for the same species that had the same gene ID but different ensemble ID. Can someone let me know how this happens? for example, for these two ensemble IDs ENSCAFG00845015896, and ENSCAFG00000017519 which are both for canines there is the same symbol of "TNMD".
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.

If you search for those gene IDs in Ensembl, you will see that ENSCAFG00845015896 is in the primary assembly, but ENSCAFG00000017519 is deprecated, and apparently has never been in the primary assembly but originates from the genome annotation of one individual dog. I did not align the protein sequences, but on a quick look they seem identical, which is why both are connected to the same gene symbol. If your identifiers are not all from the same assembly (and sometimes even if they are), you may see several ids for the same gene symbol.

Dear Cloudberry, thanks for your answer!
I have another question if you please answer, you mentioned that 17519 might originated from the genome annotation of one individual dog. I thought in genome annotation for any individual, after genome sequencing and gene assembly, the gene prediction part we name regions using the refrence genome. How is that, an individual dog gene name differs from its reference genome naming? is it possible that ENSCAFG00845015896 is the updated name for ENSCAFG00000017519?

@ToTheMoon I'm not familiar with how the dog genome has been built, but I would imagine there have first been de novo assemblies from individual dogs, from which the primary assembly has been built. ENSCAFG00000017519 doesn't have a location in the primary assembly, so it wasn't mapped to the reference genome (because the reference didn't exist yet). When the 'offficial' assembly was built, genes were given new IDs to avoid confusion. So you're right, it's an updated ID for that gene.

thanks a lot for your answer.


text:
What is it about gene names starting with "LOC"?
I was struggling to use AnnotationDbi to change my ensemble ID to gene name, for datasets of three different species (human, canine, mouse). Among gene names for all three species there are genes with names starting with "LOC" like ENSG00000228037:
LOC100996583, ENSCAFG00845012021: LOC613013 and ENSMUSG00000095523
LOC100038995.
What type of genes are they? and is the naming unique through species, I mean if I have LOC100996583 in humans, is it possible to have it in canines or mice? if yes, are they doing the same thing?
Thanks in advance!

These are "genes" that don't have an official name. As explained in https://www.ncbi.nlm.nih.gov/books/NBK3840/#genefaq.Nomenclature:

Symbols beginning with LOC. When a published symbol is not available, and orthologs have not yet been determined, Gene will provide a symbol that is constructed as 'LOC' + the GeneID. This is not retained when a replacement symbol has been identified, although queries by the LOC term are still supported. In other words, a record with the symbol LOC12345 is equivalent to GeneID = 12345. So if the symbol changes, the record can still be retrieved on the web using LOC12345 as a query, or from any file using GeneID = 12345.

In other words, when a gene has no known function or homologs in other species, when all we know about it is that it is a locus that seems to be actively transcribed, then it gets the name of LOC plus a numerical gene ID. So these genes, the LOCNNNN are essentially lesser known and studied genes. Some of them may graduate to a full gene name, others might not.
As for uniqueness, no absolutely not. GenBank gene IDs are unique, they don't imply homology. Only one gene can have a given ID. For example, gene ID 1 corresponds to human A1BG, alpha-1-B glycoprotein. The mouse homolog of human A1BG, mouse A1bg, alpha-1-B glycoprotein, has the gene ID 117586. The numerical value of the ID tells you nothing at all. It is just a unique identifier for a specific gene in a specific species.
Finally, a more general note, homology is rarely as straightforward as your question suggests you might think it is. Even actually homologous genes often have different functions with specific functions having been gained or lost in a given lineage. So it is very hard to know and never safe to assume that two homologous genes in different species "do the same thing".

Thanks a lot for your well-described answer.


text:
What does the absence of a variant in a VCF file mean?
I have a individual.snp.vcf.gz file of an individual genome and the referencegenome.snp.vcf.gz file of the reference genome.
When I run the following code on the individual genome
gunzip -c individual.snp.vcf.gz | grep -v '^##' | awk 'BEGIN {OFS="\t"} {print $1, $2, $3, $4, $5, $6, $7}' | awk 'BEGIN{OFS="\t"} {print $0 , "SNP"}' | sed '1s/SNP/Type/'

I get this:
#CHROM  POS ID  REF ALT QUAL    FILTER  Type
chr1    10067   .   T   A   0   LowGQX  SNP
chr1    10070   .   C   A   0   LowGQX  SNP
chr1    10073   .   T   C   0   LowGQX  SNP
chr1    10105   .   A   C   0   LowGQX  SNP
chr1    10108   rs62651026  C   T   0   LowGQX  SNP
chr1    10157   .   T   C   0   LowGQX  SNP
chr1    10177   rs201752861 A   C   17  LowGQX  SNP
chr1    10180   rs201694901 T   C   0   LowGQX  SNP
chr1    10181   rs1246412344    A   T   0   LowGQX  SNP
chr1    10250   rs199706086 A   C   0   LowGQX  SNP
chr1    10257   rs111200574 A   C   0   LowGQX  SNP
chr1    10285   rs866375379 T   C   67  LowGQX  SNP
chr1    10321   rs1002315756    C   T   0   LowGQX  SNP
chr1    10327   rs112750067 T   C   36  LowGQX  SNP
chr1    10332   rs1175748383    C   A   0   LowGQX  SNP

When I run the equivalent code on the reference genome
gunzip -c referencegenome.snp.vcf.gz | grep -v '^##' | awk 'BEGIN {OFS="\t"} {print $1, $2, $3, $4, $5, $6, $7}' | head -n100 | awk 'BEGIN{OFS="\t"} {print $0 , "SNP"}' | sed '1s/SNP/Type/'

I get this:
#CHROM  POS ID  REF ALT QUAL    FILTER  Type
1   10001   rs1570391677    T   A   .   .   SNP
1   10002   rs1570391692    A   C   .   .   SNP
1   10003   rs1570391694    A   C   .   .   SNP
1   10008   rs1570391698    A   G   .   .   SNP
1   10009   rs1570391702    A   G   .   .   SNP
1   10015   rs1570391706    A   G   .   .   SNP
1   10020   rs1570391708    A   C   .   .   SNP
1   10021   rs1570391710    A   G   .   .   SNP
1   10026   rs1570391712    A   C   .   .   SNP
1   10027   rs1570391716    A   C,G .   .   SNP
1   10032   rs1570391720    A   C   .   .   SNP
1   10033   rs1570391722    A   G   .   .   SNP
1   10039   rs978760828 A   C   .   .   SNP
1   10043   rs1008829651    T   A   .   .   SNP
1   10045   rs1570391729    A   C,G .   .   SNP
1   10051   rs1052373574    A   C,G .   .   SNP
1   10055   rs892501864 T   A   .   .   SNP
1   10056   rs1570391738    A   C   .   .   SNP

Now I want to get the genotype of the individual genome at all the sites listed in the reference genome table, so my question is: When I don't find a particular variant ID in the individual genome table (for example, rs1570391716), can I assume that the genotype of the individual is equal to the homozygous REF (AA, for rs1570391716) in the reference genome table? How certain can I be (assuming a low sequence reading error)?

What exactly is the referencegenome.snp.vcf.gz file? Is that a list of all known, reported, SNPs in this particular species? And how did you obtain the individual file? Note how all of the variants you show have a QUAL of 0 and are marked with LowGQX, it is extremely unlikely that those are real calls.

Just FYI, there is no need for gunzip, you can just use zgrep -v '^##' vcf.gz | awk ....

Yes, the referencegenome.snp.vcf.gz file is the list of all known, reported, SNPs in this particular species. And the individual file was obtained from the sequencing company (I don't know how exactly they obtained it) but I believe it was obtained through some standard variant calling methods.
Regarding the QUAL of 0, let's ignore that for now.

The absence of a variant means the variant was not found. However, that does not mean you can assume the position in the sequenced individual was homozygous for what is in the reference. Perhaps the position was simply not sequenced. Or the coverage was too low. So your first step will have to be to focus only on the regions covered by whatever sequencing assay was used to sequence your data, and only those regions with sufficient coverage for a variant to be called.
There is no magic number here, the threshold will depend on how the data were generated. For instance, a minimum coverage of 8 is a reasonable lower bound, but these days, we often have coverage values of 60x or more, and in the case of amplicon assays, we can have coverage in the thousands. So the threshold is something you will need to decide on.
Next, you really need to take the quality and filters into account. You haven't shown us the full VCF line, but my bet is that all the variants you show in your individual file are most likely not genotyped at all. These are not actually present. They all have a score of 0 and are marked with LowGQX which means:

The genotyping quality locus (GQX) is less than 30 or not present.

And GQX is:

The minimum of genotype quality assuming variant position and genotype quality assuming nonvariant position.

Variants with a QUAL score of 0 and low genotype quality are almost certainly not variants at all and you really shouldn't consider them.
So, since you apparently don't have the raw data, neither fastq nor bam, that means you cannot easily filter by coverage in the various positions. However, at the very least, you need to find out what sequencing assay and method was used to sequence your sample. Then, you can get a bed file with the regions of interest, and then limit your analysis to variants in those regions.
Also, I urge you to have a look at tools like bedtools and bcftools that will help you filter your data.
Finally, while your command will work, it is a bit inefficient, and you could simplify
gunzip -c individual.snp.vcf.gz | grep -v '^##' |
  awk 'BEGIN {OFS="\t"} {print $1, $2, $3, $4, $5, $6, $7}' |
    awk 'BEGIN{OFS="\t"} {print $0 , "SNP"}' | sed '1s/SNP/Type/'

to
zgrep -v '^##'  individual.snp.vcf.gz |
  awk 'BEGIN {OFS="\t"} {print $1, $2, $3, $4, $5, $6, $7, NR==1 ? "Type" : "SNP"}'

Thank you for the answer!
I do have the fastq1 and fastq2 files. How should I proceed to filter by coverage in the various positions?

@williantafsilva you would have to align them to a genome to get a bam file and then you can filter the bam. If you need help with that, read gatk best practices and if you still need help ask a new question and we will try and help out.

Thank you! Very helpful!


text:
Continuous predictor variable to identify responder genes using DESeq2
I have samples from 45 individuals at 2 time points and these individuals were given a specific diet. The BMI of these individuals were recorded at both time points. The responders to the diet are the ones that show decrease in BMI from T1 to T2. I want to use continuous variable (BMI) to identify genes that are predictors of response to diet. I also have other covariates that need to be controlled like age & gender. The metadata table looks like following

I am not sure the design that I should go for. The one that I think is appropriate is
dds <- DESeqDataSetFromMatrix(countData=cts,colData=meta,design=~0+Time+Gender+AGE+Time:BMI)
dds <- DESeq(dds,test="LRT", reduced=~Gender+AGE)
res=results(dds, contrast=list("T1.BMI","T2.BMI"))

Is this the correct way to go ahead and to identify responder genes to the diet?
All suggestions are very welcome.
Thank you


text:
SHAP correlation patterns for model prediction of protein are confusing
I have developed an ANN model to predict protein post-translational modification patterns from Elisa data. For simplicities sake, the how, feasibility and parameters aren't important for my question and some details were omitted.
For a given protein, I will call it protein X, it has ubiquitin as a modification but no phosphorylation patterns.
I trained the ANN with a variety of post-translational modification patterns but a KEY PIECE OF INFO: my training data did not contain any ubiquitin patterns (just assume there's a reason)
So when I performed an ELISA with a set of antibodies and antibody a is specific for ubiquitin patterns and antibody b is specific for phosphorylated patterns. When I ran the prediction of protein x modification patterns using antibody a, antibody b (and others) we got pretty good accuracy.
To explain how the model worked, I ran a SHAP analysis and a bipartite graph to show the feature importances (antibodies) and the modifications but got confounding results

For antibody a, there were positive and some negative SHAP values to other modification patterns except for ubiquitin, which it is specific for

For antibody b, we also found positive and negative SHAP values to modification patterns except phosphorylation, which protein x doesn't really have.

So how do I explain why SHAP produced this pattern: 1)antibody a did not have any SHAP correlation to its target, ubiquitin but it did to others, 2) antibody b also did not have SHAP correlations to its target but to others instead.
Again this is confounding because I expected it antibody a to have SHAP correlations to ubiquitin but not the others, and then antibody b should not have any SHAP correlations since it targets phosphorylation but protein x has no phosphorylation.
Im not too convinced or can't relate some of the limitations with SHAP that it shows the hidden patterns/relationships for the model, not what we expect in "real life"
Can someone provide a more nuanced insight to this observed SHAP analysis?


text:
Confusing naming in MANE GFF file
I'm working with the MANE v1.2 annotation (with RefSeq keys) and am puzzled about the ID naming convention being used. Specifically, there are several transcripts (18 by my count) that have hyphenated numbers attached to their IDs, which are otherwise just RefSeq IDs. For example, the entry NM_000937.5 (for POLR2A) is listed as "NM_000937.5" in FASTA, but as "NM_000937.5-2" in the GFF. In fact, everywhere in the GFF any reference to POLR2A is listed as "POLR2A-2". This perhaps could make sense if there were multiple transcripts included for the gene, but it is not the case - there are no plain "POLR2A" entries nor any "POLR2A-1" entries.
What is going on here? Why do these entries have additional numbers attached?

For reference, the 18 transcripts with hyphenated names are:
MANE ID                 RefSeq ID
rna-NM_001405666.3-2    NM_001405666.3
rna-NM_001291281.3-2    NM_001291281.3
rna-NM_001137667.2-2    NM_001137667.2
rna-NM_017680.6-2       NM_017680.6
rna-NM_002457.5-2       NM_002457.5
rna-NM_001164277.2-2    NM_001164277.2
rna-NM_178167.5-2       NM_178167.5
rna-NM_000937.5-2       NM_000937.5
rna-NM_001324073.3-2    NM_001324073.3
rna-NM_001097639.3-2    NM_001097639.3
rna-NM_152503.8-2       NM_152503.8
rna-NM_001372044.2-2    NM_001372044.2
rna-NM_001417890.1-2    NM_001417890.1
rna-NM_022555.4-2       NM_022555.4
rna-NM_021983.5-2       NM_021983.5
rna-NM_014513.3-3       NM_014513.3
rna-NM_020535.3-5       NM_020535.3
rna-NM_014512.1-9       NM_014512.1


text:
How to Calculate Delta Phi Psi between aligned pair of residue between Reference and Target Protein Structure in Pymol
I aligned multiple protein structures from the same protein family to a reference structure in pymol and now I want to calculate the Difference in Phi psi values for each 'aligned' residue pair between Reference Structure and each Target Structure, I have extracted the phi psi values using phi_psi command and I have the sequence alignment but I don't know how to extract Delta Phi Psi values for every aligned residue between Reference and Target.

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


text:
What does this BLASTN mean: Error pre-fetching sequence data?
I am trying to run BLASTN on a 14GB Pseudomonas database using the command: blastn -db ~/pseudomonas-db/pseudomonas -query ~/desktop/inputfiles/poz176.fasta -out poz176-out -num_threads 8
I get the following error:
Error: NCBI C++ Exception:
    T0 "c++/include/corelib/ncbidiag.hpp", line 980: Critical: (CCoreException::eNullPtr) ncbi::CObject::ThrowNullPointerException() - Attempt to access NULL pointer.
    T0 "c++/include/corelib/ncbidiag.hpp", line 789: Critical: (CCoreException::eNullPtr) ncbi::BlastFormatter_PreFetchSequenceData() - Error pre-fetching sequence data
     Stack trace:
      /usr/lib/ncbi-blast+/libxncbi.so ???:0 ncbi::CObject::ThrowNullPointerException() offset=0xA6 addr=0x7fa1ea902f76
      /usr/lib/ncbi-blast+/libxobjmgr.so ???:0  offset=0x209609 addr=0x7fa1eb653609
      /usr/lib/ncbi-blast+/libxobjmgr.so ???:0 ncbi::objects::CBioseq_Handle::IsSetInst_Length() const offset=0x9 addr=0x7fa1eb653a09
      /usr/lib/ncbi-blast+/libxobjmgr.so ???:0 ncbi::objects::CBioseq_Handle::GetBioseqLength() const offset=0xE addr=0x7fa1eb653dce
      /usr/lib/ncbi-blast+/libxblast.so ???:0 ncbi::blast::LoadSequencesToScope(std::vector<ncbi::objects::CSeq_id_Handle, std::allocator<ncbi::objects::CSeq_id_Handle> >&, std::vector<ncbi::CRange<unsigned int>, std::allocator<ncbi::CRange<unsigned int> > >&, ncbi::CRef<ncbi::objects::CScope, ncbi::CObjectCounterLocker>&) offset=0x329 addr=0x7fa1ebd90e79
      blastn ???:0 ncbi::BlastFormatter_PreFetchSequenceData(ncbi::blast::CSearchResultSet const&, ncbi::CRef<ncbi::objects::CScope, ncbi::CObjectCounterLocker>, ncbi::blast::CFormattingArgs::EOutputFormat) offset=0x5FE addr=0x555af7fb31ee
      blastn ???:0 CBlastnApp::Run() offset=0xFDC addr=0x555af7fa99ac
      /usr/lib/ncbi-blast+/libxncbi.so ???:0 ncbi::CNcbiApplicationAPI::x_TryMain(ncbi::EAppDiagStream, char const*, int*, bool*) offset=0x193 addr=0x7fa1ea84edd3
      /usr/lib/ncbi-blast+/libxncbi.so ???:0 ncbi::CNcbiApplicationAPI::AppMain(int, char const* const*, char const* const*, ncbi::EAppDiagStream, char const*, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&) offset=0x61D addr=0x7fa1ea85032d
      blastn ???:0 main offset=0x6B addr=0x555af7fa7e5b
      /lib/x86_64-linux-gnu/libc.so.6 ???:0 __libc_start_main offset=0xEA addr=0x7fa1ea38ed0a
      blastn ???:0 _start offset=0x2A addr=0x555af7fa876a

I tried changing memory allocation but that didn't seem to do anything. I'm not sure what this error is telling me. How do I get around this?
(Had posted this on stackoverflow and was asked to post here)

Your blast run command looks fine. I'd look at the answer below. You could try running it on a bigger machine (more CPU), but I think you'll get the same answer.

Based on other posts elsewhere it seems likely that there is something funky about your database. How did you generate it?
Could you try to re-generate the db without the -parse-seqids flag?

I agree, I think this is the error output Attempt to access NULL pointer. is a bad error and suggests there's a database problem. This doesn't appear to be a memory allocation issue. 14GB is big, but there's nothing in the error message to raise concern about that


text:
Variants from multiple tools: normalization before or after annotation with VEP?
We want to compare some variant calling tools and their calls on Whole Exome Sequencing data. We will have to normalize the variants that are called (.vcf format) before comparison using for instance vt normalize or hap.py. We are comparing with a gold standard too and use hap.py for this now to get some benchmark information.
Edit: with normalization I mean representing the variants in such a way, we are able to compare them: to ensure that the representation of a variant is “parsimonious” and “left-aligned.”
I was just wondering, and I cannot find any real guidelines (yet), when we annotate using VEP, do we have to normalize before or after annotation to keep the annotation correct as well. I know VEP has some normalization functionality too see here. So might we even not need an additional normalization tool before and just let VEP do it for us?
I am just wondering about the guidelines and what what others do.
Thanks in advance,

What kind of normalization are you going to do? VEP talks about normalizing variant represenations, variant names and left-aligning indels. However, all variant callers should be producing a vcf, so there is no need to normalize the names. Indels perhaps.  You might want to normalize other things like variant quality but we need to understand what you mean before we can help. Can you please [edit] your question and clarify?

If you're referring to left alignment + parsimonious representation, do it before any annotation. Most if not all annotation happens by comparing CHROM, POS and optionally REF and ALT and since vt normalize affects POS, you should do it before annotation.
I'm not sure how smart VEP is, but most annotators are simple lookup tools at least for third party databases, so they're not built to normalize while looking up.

Yes, that makes sense. And VEP has some form of normalization before annotation: here, but it seems to be more basic than a normalization tool like vt normalize.

Yes, vt normalize and bcftools norm both preserve a record of the pre-normalization entry so data loss is avoided. I'm not sure if VEP does that or even logs it.

Thinking a bit more about your comments ^, why would you not always normalize before annotation to be sure you would get the same results over different studies (hypothtically)? Would "the same variants" called with different callers, cause discrepencies in position and therefore annotation if you do not normalize them? Or am I missing something?

You're asking the right question. Without normalization, certain annotations could be missed especially if POS changes post normalization. Most callers should call a variant at the same location based on an alignment, but I'm not sure if that can be taken for granted. I always decompose and normalize all my VCFs pre annotation.

I guess for SNPs this might be (more) the case that they would have the same location between callers, INDELs are a different story likely

SNVs/MNVs that occur in a repeat region can confuse callers a little. One also needs to keep in mind that certain variants were named/annotated way back when and used to be called using the wrong coordinate for quite a while. I've seen this happen in old medical records, not sure how prevalent this is these days.


text:
tabix to subset each VCF by the coordinates given
I have a vcf file like below
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  2072184
chr1    10109   .   AACCCT  A   21.51   PASS    AC=2;AF=1.000;AN=2;DP=10;FS=8.451;MQ=13.01;MQRankSum=0.000;QD=2.15;ReadPosRankSum=-1.300;SOR=3.236;FractionInformativeReads=0.700;R2_5P_bias=0.000  GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    ./.:3,4:0.571:7:1,1:2,3:6:30,8,0:2.1511e+01,5.8560e+00,1.3473e+00:0.00,7.00,10.01:1,2,4,0:3,0,1,3
chr1    10230   .   AC  A   29.28   PASS    AC=1;AF=0.500;AN=2;DP=21;FS=12.601;MQ=20.22;MQRankSum=-0.715;QD=1.39;ReadPosRankSum=1.109;SOR=4.407;FractionInformativeReads=0.667;R2_5P_bias=-3.357    GT:AD:AF:DP:F1R2:F2R1:GQ:PL:GP:PRI:SB:MB    0/1:6,8:0.571:14:1,4:5,4:22:59,0,19:2.9283e+01,2.9815e-02,2.2470e+01:0.00,30.00,33.01:3,3,0,8:2,4,4,4

I am going to subset this VCF by the coordinates given in a bed file using tabix
This is how my bed file looks
chr7    127471196  127472363
chr7    127472363  127473530
chr7    127473530  127474697
chr7    127474697  127475864

This is my command
tabix -h filevcf.gz -R XXX.bed > file.vcf.gz_extracted.vcf.gz

But output file is empty after headers I can not find variants in the output
I look at output by bcftools view -H file.vcf.gz_extracted.vcf.gz > output.txt which is totally empty
Where I am making a mistake here?
I got results when I used some of positions in input vcf

Do you have a typo in the name of the file? it should be file.vcf.gz (it's missing a period).

Also: is the file block-gzipped and indexed? (you can check that it's block gzipped by running "file" on it....if the results is "Blocked GNU Zip Format (BGZF; gzip compatible)" )

Could you clarify what you mean by: "got results when I used some of positions in input vcf"? used how?

From what you show us here, we see a vcf that lacks headers and a bed that describes positions not present in the vcf. So empty results are expected. Please show us an example that you expect to produce output and then we can actually test.


text:
Can anyone share insights or provide the complete mRNA sequence for tdtomato-WPRE for scRNA seq analysis?
I'm currently immersed in research involving the Rosa26-tdtomato (Ai9) mouse line. We performed scRNA seq and analysis focusing on tdtomato-expressing cells but faced a challenge with a lower-than-expected number of positive cells.
Our alignment process, using the tdtomato sequence (Genbank accession: AY678268), seems insufficient. Fellow researchers suggest success in similar cases by incorporating the full mRNA sequence, including the 3´UTR, along with tdtomato-WPRE mRNA.
I searched everywhere for the full mRNA sequence but couldn't find it. Can anyone share insights or provide the complete mRNA sequence for tdtomato-WPRE? This would greatly enhance our alignment process and improve the accuracy of tdtomato positive cell detection.
Appreciate your contributions to advancing our research!

I'd suggest looking at the mouse itself, which appears to be here on the JAX website.
This link includes RNA-seq experiments, so presumably you should be able to compare your pipeline and data to previous work to see whether there's anything weird there.
Looking at the original paper, you can see the structure of the transgene which might be enough info to work from, there is also information in the supplement. You should in principle be able to figure out the sequence of each of those components and include them.
If you are really having issues, you could do a transcriptome assembly to see whether you can recover the locus.
You can always write to the senior author asking about this line, or try to figure out why your transcript may truly not be expressed.

The Genbank sequence you've specified appears to be for RFP, not tdTomato.
Have you tried one of the other tdTomato sequences, for example a lentiviral sequence, or a sequence that is called tdTomato?
https://www.ncbi.nlm.nih.gov/nuccore/KJ081243.1