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

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

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

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

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

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

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

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

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

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

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

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

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

This should work
import argparse

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

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

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

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

exit()

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

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

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


text:
Help with finding an unsupervised learning algorithm for RNA sequence/structure pattern identification
I have a bunch of RNA sequences (and their optimal secondary structures) and their corresponding energy values (measured by mean free energy) and I'm trying to find a way to identify features (patterns in their structures or sequences) common between samples of similar energy values. the overall goal is to try and identify sequence or structural elements in the RNA samples that tend to correlate with low mean free energy values. would be grateful if anyone could recommend me algorithms to look into using - I'm assuming this would be an unsupervised learning project and I only have experience with supervised stuff, but I'm looking into PCA right now and not sure if that'd be useful. should I be looking into something else?
Intuitively I'm imagining it as like a clustering + feature extraction problem where I have a bunch of dots, each representing an RNA sample, and then the axis represents energy so the dots could be clustered in energy value similarity and then within each energy value cluster patterns/relationships could be found between the sequences and structures of the samples. but not sure if an ML algo exists to handle this and PCA seems like not what I'm looking for because the axes would be principal components and not energy...

It depends what is meant by "features". Features can be a particular structure component in the RNA structure. In ML/DL these "features" are called instances (or observations). "Features" is also a technical word used in machine/deep learning, which can be completely different, it means the data or statistical variables associated with an instance.
Feature weight extraction (ML/DL definition) using unsupervised learning is complicated, but doable. This is where the features corresponding to each entry are weighted, i.e. quantified. Under traditional unsupervised learning this is not achievable, unless you want to include/exclude "features" and assess cluster integrity. I personally considered it a lot of work for no real interpretable results.
If the features are components of the RNA structure (what would be technically called 'instances' in ML/DL) then the general traditional unsupervised learning following PCA analysis is tSNE. However, tSNE is often used to get better resolution the PCA feature reduction output, i.e. tightens the clusters. PCA->tSNE is a cool technique and works well. PCA and tSNE can give different clusters by the way. As a rule of thumb I consider PCA > tSNE.
Multi-dimensional scaling etc ... needs a defined analytics aim.

Information on tSNE and PCA this online resource looks fine
https://www.datacamp.com/tutorial/introduction-t-sne
The above is if you want to do tSNE and compare the results to PCA this link is good.
In terms of using tSNE to sharpen PCA clusters I always used O'Reilly's "Introduction to Machine Learning with Python". There is a very clear example of Python code of the technique there, towards the end of the book and why it works.

Thanks! Do you have any resources I could look into that would help me understand how to frame this as a PCA problem as someone new to it? If not no worries

Makes sense re: edit, thank you!

Dear @joshua213p I have answered below the --------

Thanks @joshua213p

hey! reading through the t-sne blogpost and I'm a little unclear on how i'd use either of these algorithms here for my goals. i get it'd cluster my samples based on similarity in terms of first and second principal components/t-SNEs, but what could i then use identify structural or sequence patterns within each cluster? please let me know if I'm missing something - thanks!

Hi @joshua213p can you submit this as a separate question please?

source: Stack Exchange

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

sort: votes


text:
Applying PCA or tSNE on non-numerical data
Reading through this t-sne blogpost and I'm a little unclear on how i'd use either of these algorithms here for my goals. i get it'd cluster my samples based on similarity in terms of first and second principal components/t-SNEs, but what could i then use identify structural or sequence patterns within each cluster? please let me know if I'm missing something - thanks!
I'm also curious as to how I would apply PCA or t-SNE on my data considering its non-numerical nature. I was thinking of applying one-hot encoding, but is there a generally more suggested approach? I've attached a screenshot of a few of my dataset rows for context (for the sake of having the table fit, I chose rows with simpler structures/zero MFE values, but most rows have complex structures and non-zero MFE values). Thanks.

This is based on my previous post:

I have a bunch of RNA sequences (and their optimal secondary structures) and their corresponding energy values (measured by mean free energy) and I'm trying to find a way to identify features (patterns in their structures or sequences) common between samples of similar energy values. the overall goal is to try and identify sequence or structural elements in the RNA samples that tend to correlate with low mean free energy values. would be grateful if anyone could recommend me algorithms to look into using - I'm assuming this would be an unsupervised learning project and I only have experience with supervised stuff, but I'm looking into PCA right now and not sure if that'd be useful. should I be looking into something else?

Intuitively I'm imagining it as like a clustering + feature extraction problem where I have a bunch of dots, each representing an RNA sample, and then the axis represents energy so the dots could be clustered in energy value similarity and then within each energy value cluster patterns/relationships could be found between the sequences and structures of the samples. but not sure if an ML algo exists to handle this and PCA seems like not what I'm looking for because the axes would be principal components and not energy...

I'd recommend PCA, tSNE the PCA (first)+ tSNE (second).
The sequence data - one hot key encoding is a good idea.
import pandas as pd
datadict = {'UTR':[.....],'Sequence':['AGAA','ACAG' ...], 'Optimal Secondary Structure':[....]}
df = pd.DataFrame(datadict)
df_onehot = pd.get_dummies(df['Sequence'])
df = pd.concat([df, df_onehot], axis=1)     # add the one hot encoded column to master dataframe
df.drop(...) # remove the original 'Sequence' column

The other issue is vectorising the "dots". This is really a separate question, here's what I'd do ...
For optimal secondary structure, maybe simply count the dots?
count = df_series.str.count(".")

This only works on a series so you'll need to do something like df_list = df['Dots column'].to_list() and then df_series = pd.Series(df_list). Finally the count will need to be concat back onto the df and the original Dots deleted.
The thing I'm not clear about is your 'instances' that could be 'UTR_info' thats the stuff you are calling you data or "features".

Got it, but how does pca and tsne get me closer to actually identifying structural/sequence patterns present in the created clusters?

Thanks. You can now upvote BTW. This is about the fundamental design of your experiment. If post another question I'll sketch an outline about what unsupervised learning achieves. I don't exactly know your project aims (I can guess) but I know what unsupervised learning delivers.

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


text:
What are the output files of RNA-Seq from facility?
This question was also asked on Reddit
I am new in our lab and I am going to do bulk RNA-Seq. What type of files will we get from the company (Genewiz)? Will it be a bunch of Fastq files? or they give something like bam files?

Welcome to the site. Genewiz is a commercial company which you are outsourcing to and I strongly recommend contacting them directly. Essentially you are paying for this service so I think it is fair to request this information directly. Their contact details are here https://www.azenta.com/contact-us ... The implications of the question is that this site provides  technical support for a commercial entity.

Please do not cross-post https://www.biostars.org/p/9571387/ and if you do at least leave a comment with the cross-reference. It's just double effort for people in both communities, especially for these sorts of very basic questions.

OP has deleted their post on biostars.

This cannot be predicted unless you have interacted with them previously. Most commonly a sequencing facility will send raw data as gzipped fastq files, often together with md5 checksums so the user can check file integrity after transfer/download from them.
We worked with a facility before who sent data in unaligned bam format, which then for typical workflows first had to be converted back to fastq.
After all, you're a paying customer so feel free to shoot them an email and ask whatever question comes up. Ask what exactly they will provide so that you can prepare your workflows accordingly. Most likely it will be fastq.gz.


text:
Inconsistent microarray expression levels after normalizing with log2
I know it is necessary to normalize the GEO microassay data before analysing. Howerver, after using log2 to normalize the data and checking with boxplot, I found the data is still not align in the same level:

I notice the expression in some samples is lower overall, can I add a value manually to it?
Update: As @ATPoint had pointed out, log2 is not a normalization method but transformation. So it is normal to see the median number of sample data not at the same level. Actually, I need to normalize it via limma or other R package.

As I already told you here (When analysing microarray data, is it need to do normalization and standardization both?) log2 is a transformation, not a normalization so unless these data were already normalized you cannot expect equal medians.
What the method section tells is

"The scanned images were analyzed with Feature Extraction Software 9.5 (Agilent) using default parameters (protocol : GE1-v5_95_Feb07  and Grid: 014850_D_F_20060807) to obtain background subtracted and spatially detrended Processed Signal intensities. Features flagged in Feature Extraction as Feature Non-uniform outliers were excluded."

So that is probably what these values are. Hence, you might want to apply something like quantile normalization (for example via the limma package) to the log2-transformed data before applying limma for differential testing or any other downstream normalization. Read the limma manual. Note that I did not invest much time to read my self into that experiment. Treat advise with care and collaborate with local experienced people if you need hands-on guidance.

Thank you for your patience in explaining. I confused these two concepts and mistakenly thought that log2 transformation is also a form of normalization. I will try to normalize the data via limma package after log2-transformation.


text:
When analysing microarray data, is it need to do normalization and standardization both?
When we do differential expression analysis, it is common to normalize the data using a log transform:
exprs(gse) <- log2(exprs(gse))

This ensures that the expression data are more comparable and consistent.
However, I have recently learnt that there is a method called standarization which can also be used to preprocess data. So what is the difference between these two methods? Can I just do normalization for differential expression analysis?

Actually, there is three terms to introduce, which is normalization, transformation and standardization.
Normalization corrects for biases such as different baselines of array intensities, or in RNA-seq different read counts and library compositions.
What you do with log2 here is a transformation, because it brings data to a different scale, so it transforms. log2 is necessary because on normal scale the variance of genes is an almost linear transformation of the magnitude of intensity values (or counts in sequencing experiments), but what you want is not this technical variance, but an information for variability due to treatment or group effects.
Then standarization is strictly speaking another transformation that brings data from an absolute scale to a relative one. This is also called Z-scoring. Data are transformed to represent the deviation from the mean. See: https://en.wikipedia.org/wiki/Standard_score
To the question "can I do...", you should read established analysis pipelines, for example the limma documentation, as different arrays need different preprocessing, and this is all established and standard these days so don't reinvent the wheel. Generally, you always analyze arrays on log2 scale.


text:
How to identify index sequences for cutadapt for atac-seq
I am just starting on my atac-seq analysis (not very experienced, so my apologies for a stupid question), and I am wanting to trim the adapters of off my paired-end reads.
I am following the galaxy atac-seq protocol that says to enter in the custom 3' adapter sequence. I took the primer sequences from Buenrostro.
>Ad1_noMX
AATGATACGGCGACCACCGAGATCTACACTCGTCGGCAGCGTCAGATGTG
>Ad2.7_CTCTCTAC
CAAGCAGAAGACGGCATACGAGATGTAGAGAGGTCTCGTGGGCTCGGAGATGT

Following the ATAC protocol, I see that the primer sequences they used are the reverse complement sequences of the end of the provided sequences (ACATCTCCGAGCCCACGAGAC/CACATCTGACGCTGCCGACGA). I took this part and added the reverse complement index primer of 2.7 (GTAGAGAG) to the front.
So for my R1 my sequence is:
GTAGAGAGACATCTCCGAGCCCACGAGAC
And for my R2 my sequence is:
GTAGAGAGCACATCTGACGCTGCCGACGA
However, when I run the command I only get 1.4% read1/2 with adapter. Could anyone suggest where I am going wrong?
Thank you

The sequence to trim for ATAC-seq is CTGTCTCTTATACACATCT. It's the Nextera adapter sequence. The sequence is the reverse complement of the 3' ends of your primers, both forward and reverse.
You trim adapters which come before any index sequence so trimming adapter (and by this everything beyond it) is sufficient.


text:
Single Node Deletion in Cheng and Church Algorithm
Cheng and Church Algorithm is a biclustering algorithm or an unsupervised learning task with which both rows and columns in a data matrix can be placed in concurrent clusters.
In the original paper, authors explained that a row or column of an input matrix $A$ with the largest $d(i)$ or $d(j)$ should be removed in single node deletion phase:

Find the row $i \in I$ with largest: $d(i) = \frac{\sum_{j \in J} (a_{ij} - a_{iJ} - a_{Ij} + a_{IJ})}{\lvert J \rvert}$
Find the column $j \in J$ with largest: $d(j) = \frac{\sum{i \in I} (a_{ij} - a_{iJ} - a_{Ij} + a_{IJ})}{\lvert I \rvert}$

I found one implementation of the algorithm trying to do the following operations for single node deletion:
row_max_msr = np.argmax(row_msr) #finding the index of rows having largest d(i)
col_max_msr = np.argmax(col_msr) #finding the index of columns having largest d(j)

#Continue with removing row where $d(i) >= d(j)$ or column $d(j) > d(i)$
if row_msr[row_max_msr] >= col_msr[col_max_msr]:
    row2remove = row_indices[row_max_msr]
    rows[row2remove] = False
else:
    col2remove = col_indices[col_max_msr]
    cols[col2remove] = False

I was wondering whether anyone has a suggestion for this matter. I don't mean that the code has installation or execution problem. Actually, I want to be sure of the logic to be aligned with the algorithm.

Padilha, V. A. & Campello, R. J. G. B. (2017). A systematic comparative evaluation of biclustering techniques. BMC Bioinformatics, 18(1):55.

https://github.com/padilha/biclustlib
It needs the r-bicluster library which is available on conda https://anaconda.org/conda-forge/r-biclust
Citing the authors GitHub page for instructions
import numpy as np

from biclustlib.algorithms import ChengChurchAlgorithm
from biclustlib.datasets import load_yeast_tavazoie

# load yeast data used in the original Cheng and Church's paper
data = load_yeast_tavazoie().values

# missing value imputation suggested by Cheng and Church
missing = np.where(data < 0.0)
data[missing] = np.random.randint(low=0, high=800, size=len(missing[0]))

# creating an instance of the ChengChurchAlgorithm class and running with the parameters of the original study
cca = ChengChurchAlgorithm(num_biclusters=100, msr_threshold=300.0, multiple_node_deletion_threshold=1.2)
biclustering = cca.run(data)
print(biclustering)

Thank you for the response, however, I don't mean that the code has installation or execution problem. Actually, I want to be sure of the logic to be aligned with the algorithm.


text:
Understanding DESeq2 output for resultNames
I'm following along with some pre written code to learn how to do differential expression. I am using the same data set and have gone back through the code and copy and pasted what they used and the follow along gets the two headings rather than the single one that I am getting, so I'm a bit confused how I'm getting different results.
I have some data which I have put through DESeq but the output for resultsNames isn't what I expect. The code is
dds <- DESeqDataSetFromMatrix(countData = exprs(bottomly.eset), colData = pData(bottomly.eset), design = ~ strain)
dds <- DESeq(dds)

but this code produces this output "Intercept"  "strain_DBA.2J_vs_C57BL.6J" whereas the output that I intended was "Intercept"  "strain_DBA.2J"   "C57BL.6J"
Is there a simple way to change this within the code that I have or is this something related to the dataset I used?

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.

In your expectation, what is the DBA.2J strain tested against?

DESeq2 is a package designed for calculating differential expression, in other words for comparing expression between two states. The primary results it produces are therefore comparative values (e.g. comparing DBA.2J and C57BL.6J in your case).
Given that this does not seem to be your intended calculation, it would be helpful if you could provide more context around what you want to use DESeq2 for. If it's not for differential expression, there's likely a better tool / approach for solving your problem.

So I was following along with a coding lesson from Youtube learning how to do differential expression in R, if the data had been my own I probably would have been able to figure out what needed to be done.
The code that worked for me was
dds <- DESeqDataSetFromMatrix(countData = exprs(bottomly.eset), colData = pData(bottomly.eset), design = ~ strain)
dds <- DESeq(dds, betaPrior = TRUE)
I was then able to use this code below which was resulting in the error previously
contrast.deseq2 <- list(c("strainC57BL.6J", "strainDBA.2J"))
Thank you to everyone who commented on this question

I don't understand the difference between this and the vs offered by DESeq2 previously in your question. Given that the default test is Wald and betaPrior is TRUE by default for Wald tests, how is this result any different?

I'm not sure of the difference but this is the result I got when running the new code

Run it both ways and see if the number of genes or fold change value differs by much. I honestly don't think there should be any significant difference.


text:
How do I delete an amino acid residue from a loop and then join the loose ends?
I want to delete an amino acid residue from a loop and then join the loose ends. I intend to increase the thermal stability of the protein by reducing the size of the loop and its agitation. Pymol just deletes the residue and breaks the loop.
I intend to increase the optimal temperature of the protein, even if it costs some reduction in activity. The objective is to reduce protein agitation due to temperature increase. I will try point mutations that increase hydrogen interactions and eventually delete a loop located glycine to reduce its length and consequently its thermal agitation. All this work is aimed at preserving the catalytic site, but some loss of activity is acceptable.
As a complement to my question, I present an excerpt from the paper by Mahdie Rahban et al (2022) ; Thermal stability enhancement: Fundamental concepts of protein engineering strategies to manipulate the flexible structure.

Loop deletion or shortening by insertion of proline or deletion of
Gly, coupled with loop anchoring by increasing the number of salt
bridges, hydrogen bonds, or hydrophobicity, are good strategies to
enhance the overall rigidity of a protein conformation leading to
better compactness and a decreased number of internal cavities.

Does anyone have a suggestion how to solve this problem? I already used Modeller to insert an amino acid residue in place of a deleterious mutation, could Modeller be used to reconstruct the loop? How do I do it?

Modeller .. no thats for point mutations

https://salilab.org/modeller/9.15/manual.pdf 2.3 Loop optimization pag 26

It depends on what software you have in hand and want to use. But basically, when you delete/add a residue, you can't just stick two parts together since it affects the whole conformation.
I suggest using AlphaFold2, and input it with different sequences (including deletion mutation) to see how it affects the structure.
Here is a link to AF2 colab:
https://colab.research.google.com/github/sokrypton/ColabFold/blob/main/AlphaFold2.ipynb


text:
Unequal dbNSFP in-silico predictor results counts to alternatives in vcf file
I would love to ask a question regarding unequal counts in the dbNSFP database and data annotated by the hs38DH reference sequence.
The data was reannotated at 2021 using the GRCh38 p13 reference with snpEFF and dbNSFP. The data was a targeted gene panel for schizophrenia.
Variations were called with freeBayes v1.3.2-44-gfce9620.

Because my data analysis didn't work with multi-variation data, I decided to improve it (technically pick max value inside the cell, creating upward bias - 10,25.5,7.52 → my analysis would be: every variation as 25.5).
I used bcftools norm -m-both -o output.vcf.gz -Oz input.vcf.gz with the intention of separating every variant on separate lines to make analysis more precise. I got multiple errors like this after that command:
Error: wrong number of fields in INFO/dbNSFP_PHRED at chr1:1034085, expected 6, found 4
Error: wrong number of fields in INFO/dbNSFP_1000Gp3_AF at chr1:32695609, expected 2, found 1

After deleting multiple columns of dbNSFP data I found out that from 13/31 in-silico prediction tools data have inconcistencies with counts of alternatives and results of dbNSFP data. When I deleted them, it splited without problems.
=="INFO/dbNSFP_PHRED" ==  ## CADD_phred values
=="INFO/dbNSFP_1000Gp3_AF"==
=="INFO/dbNSFP_REVEL_score" ==
=="INFO/dbNSFP_VEST4_score" ==
=="INFO/dbNSFP_DANN_score"==
=="INFO/dbNSFP_Eigen-phred_coding"==
=="INFO/dbNSFP_Eigen-PC-phred_coding" ==
=="INFO/dbNSFP_LINSIGHT"==
=="INFO/dbNSFP_phastCons100way_vertebrate"==
=="INFO/dbNSFP_gnomAD_genomes_AF"==
=="INFO/dbNSFP_1000Gp3_EUR_AF"==
=="INFO/dbNSFP_ExAC_NFE_AF"==
=="INFO/dbNSFP_gnomAD_genomes_NFE_AF"==

After that i found out with vcf-validator, that this problems is not regarding only dozens variatins, but rather thousends.
So I tried to reannotate that with the dbNSFP database GRCh38.p13 v.a4.2 from SnpSift. It didn´t help; I got the same message. Then I decided to remove variations without full samples Genotypes, considering it as a possible indication of CNVs etc., but the error is still showing.

Has anyone already solved this problem where the number of predictive tool results not matching the number of alternative alleles, even though the VCF header includes a condition confirming that these numbers should match? If so, did you use the modified VCF or did you create a new one?
I'm just asking if the discrepancy in data with the header is effectively a matter of the VCF being unusable.

Thank you for any suggestions and your help

Any suggestions would be welcome.


text:
Gene Set Comparison Without Expression Data
I have been looking all over the web to find some answers to my problem but unfortunately, I was unsuccessful. I wish to determine whether an a priori defined set of genes in my case genes associated with Epidermolysis bullosa shows statistically significant, concordant differences between two biological states (e.g. phenotypes).
I have a list of genes from my own experiment and a predefined list of genes associated with Epidermolysis bullosa, and want to find out if the genes of Epidermolysis bullosa are overrepresented in my list.
Here is an article about GSEA, which is related to what I want to do:
https://www.pnas.org/doi/epdf/10.1073/pnas.0506580102

so I have the following data:
I'm researching patients with severe or mild disease; in this research, I have found a list of  200 genes after conducting a SKAT test in R (https://www.rdocumentation.org/packages/SKAT/versions/2.2.4/topics/SKAT) with a p-value <0.05 ( which I presume can cause a severe disease instead of mild  ) I want to focus on this 200 significant genes that can perhaps affect the disease level, this genes can be  -> the L list as described in the article above.
Also, the S list of predefined genes ( as is called in the article above )is associated with Epidermolysis bullosa disease. and I want to know if the genes from the Epidermolysis bullosa are overrepresented in my L list.
After reading the article and many more resources, I understood that to perform GSEA  (https://www.gsea-msigdb.org/gsea/index.jsp), I had to have gene expression data which  i don't have .
so I wish to understand if there is a tool that could do a gene set comparison without having the expression data whether it's in R or Python that uses my gene list and predefined list ( in my case associated with Epidermolysis bullosa)  ,thank you :)

You can't do a GSEA on a set alone because the differentiating statistic (or, occasionally, the rank) is used to calculate the enrichment score that is produced by the GSEA algorithm.
The other gene list type tools that I'm aware of (that my biologist colleagues use) are DAVID and PANTHER. But those are general discovery algorithms, not for comparing specific lists.

thank you , is there a way to rank the genes not based on expression data ? Or a tool that calculates the enrichment not using the   GSEA  algorithm?

@ElizaR Sounds a bit like a pathway analysis (e.g. see gprofiler). Perhaps a hypergeometric test


text:
Handling duplicates in untargeted lipidomics
I have a MS-lipidomics dataset - from human samples split into a treatment group (n = 2) and control (n = 2) - which contains several detected targets that matched to the same lipid. I don't have access to any of the raw data, only the Tier 1 and Tier 2-matched lipid names and ID numbers.
There are around 300 duplicates out of a total of 1000 detected markers, so I'm reluctant to exclude them entirely, but I couldn't find any descriptions of best-practices for what to do in this situation. I'm considering testing the correlation between duplicates and removing the ones lie below some threshold (e.g: 0.5), then averaging the remaining duplicates - but I'm not sure if this would be a valid way of approaching this.
Are there any good ways of addressing these duplicate values?

Sorry how is the threshold derived thats quite key here?

I don't know lipids, but if 1/3 of your data set comprises duplication thats a skewed distribution with a low mean and high variance (of the duplications). It is likely to follow a negative binomial distribution. You need to plot a histogram of your data and you will find that 0 (no duplicates) is massively represented, 1 duplicate is next in frequency, 2 duplicates are lower frequency ... and on until a very small proportion of the data comprises a large number of duplicates. This is a very common distribution throughout biology.
The advice is:

Do not remove the duplicates because they are a vital part of the distribution, i.e. biological information.
Construct a histogram of each data set.
Don't take "an average", i.e. don't take the mean of skewed data, thats very important if it has a high variance.

Easy solution The easiest way to deal with this data is take a median and construct upper and lower quartiles around the median (IQR), especially the frequency (i.e. the distribution duplicates). This will give an approximate comparison between data sets. This alone might work and certainly gives you "cut-offs".
If it is a negative binomial distribution then testing the significance between two negative binomial distributions beyond the IQR is an advanced area of statistics and a separate question.


text:
Cell Ranger returns too few cells!
I recently performed single cell RNA seq using core services and the web summary output from 10x that they have given shows that there are very few cells detected (100-300) cells per sample. Once this same issue happened in 1 out of 3 samples and this time it happened with 4 out of 4 samples.
What might be the issue? The core says the cell preparation was not good but the cells were sorted for viability and then used to make the library. Any help?
Full Summary snap:

I have uploaded the full html file here : https://drive.google.com/file/d/11eQduL4OC4kjC0BfhGp3Nhoci69bRAy6/view?usp=sharing

For one sample, Cell ranger 7.0.1 gave 284 cells while 6.8.2 gave 2200 cells.

Very odd, can you say something about the say something about the sample type?

Also cellranger 6.8 does not exist to my knowledge ;) can you check that

I suggest runnnjng fastqc and maybe some Samtools modules then multiqc to aid diagnosis


text:
how to search and calculate occurence of certain word (like TAA) in fastafile
>003197.2:c1179912-1178219_reverse_complement
AGGCTTTTAAAAGCCTCCTGGGTTTTTAGTGAAGCTGAGTGATAGAAGCTCTGTATTTGCATAGCGTTTT
GTCGGGCGCTTTAGCCGGCGCTTTGCCCTGGCCTGAGAGAATCTGCATTCCGTTGTATAAGGTTTTTTGT
CCTTCTGATGCTGTAGATAATTCCCCCATGTCGCGCCGGGTTCTCGCAGGACAATAATTTCCGGGCGAGC
AGAACGGTTGCCGCGACGGTAAGAAGATCGCGCTGTAAGTTATAGAGGTTATGCAGCGAGTGGTTAGACG
CCGGTCCGCTTTAACTTTGGCTAACTCCATTTGGTTTGCCATTGACGTTAGAACCGGGTCTTGTTTACCC
GTGCTGCAATAAGTTCGATAAGATTTTTCTTCAATGCTTTTGCCGCGGCTTCTTCTTGTTTTGTTGCTGG
GCATCTCTAAACGCTACTGCCGCAAAGCGATGAGCTTCTTTTGCAGGTAAGCCATCCTGCTGCTGAGTGC
GTGATGCCCGTTATGCGTGAGTGTATTTTTTATGGTTTGCCAGGGCTGGTTATTAAGCTGCTTGACCTGA
CATAAGCACTGGGAAAGATATCTTTTGCGCCGATTTTCATCTCTGCTGCAGGGAGCTGCGTGTTGGTATA
ACCGTGGACATCCACAAATTATTGGCGTGATGAATATTCTTGGTATCCCAACTGCATACGCCCTTTCCCT
ATAGGGGGAAAGCACGCCATGACGTATCCCGCAAAAAAGCGTTTTATCTTTACCGTCCTCATGCACACTC
GTGCCGCAGTTAATACTTCTTTGGCTTTGTTTTCAGCGCCGACGTGACGCAGAAGCGGATCTTTTTCATG
ACGGATACCAGTTTCAGGCTTACCGCCTCGCCCGCTAAGGCTTTGTTAAGCAACTCAGGTTTACTAAAAA
CCATGCGCGCATTTGGTCCTCGACCATCGTTCCCTCTTTGCCGAAAATATTCGACGCGGTGAGTAACCCG
CCGTCTGTAGATCGCCATCTTTATTGCGGATTTTTAAATGAATCATTTTTCCCGGCTGGGTCAACGATTG
CCAAAGCCGAGCTTGAGCGCCAGCTCATTAACACCCACATTAAATGCGGCGACGTCCGGTTTTATTTTTA
TTCAGGGCGTAAATCATTGCCTAATAACTGATGTAGCGCCTCGGCATTATAGCTATCCGATGCCTTAAGG
ATGTATTGACGACCTCATAATTATCCGGGTATTGCGCCAGCCATTCGCCAACCCAGCCACCTGGTCTGGC
AGTTTATAGGGTTCGCCGCCATCTTTATGATGTTGGTTATTTTTCCATATATCCTTAATCTGGCGCGCTA
TTTGCCGCTTTTACAATTCCAGGCGGGTACCGCGTCAATTTCATGGGCTAACATGGCAAGGCGTTGTGCG
TTAACATATGGGTCTGATGTAAGGAAATGATCTCTCGCTTGATTTCTGAATCCATCATCCCTGTACGATC
TTACCGCTATTCAGTAATACTTTTTGGAAAATTTTCTGTCCACCGCTATCCGGAAGACTACCAGGCGCAC
CACCTCTGGCGATAAATTTTTCATTACTTTGTTTCCCGCCCCGCCCGTATTTTGTTTCTGAATCTCCAGG
ATGAAGAAATGCCTTTTACTGACTGCCAAATATTTTCATCCCCAACTCGTTTTTGATAGGAAAGATTGAG
TCAAGATGTGATTA
>003197.2:c3026750-3024686_reverse_complement
TTTGATCTCGGTCTGCATACCTGGCATTATGACGGGGGGCTGAGTCCTTACACTTGTAACCATTATTAAT
CCCCTGACAGCGTCGCTGTGGCACTTTCTCTGACTGCGGAAAGATTCGCCGCAAGATTCGTGGCCTGCGT
TTAAACAATGCATCCAGACTCGCCTGTTTGATCAGCGCGGGAAAATCTTCCAGTTGCGGGCCTTTAATTT
CGCTTTCTTACCGGCGACATTATTTGAATTGGTAAAAACTTCTTTTAACGCCTCAGCGTCTTTCCCGCAT
GCGTCTTCGCCTCAGGAGAATCACTGGTGGCGTTGAGCGCTGAACGAAAGAGCCCGGCAAACTCCATTAT
AAAGCTGACTTTTCTTTCAGGCCATCTTTTATGATTTGCGCGGTATACTCTGCCCCCACCTTCATCAGTA
ATCTAACAGCTGCTTTTCCAGGTTTTCCAGCTTCGCTTCCGCTTTCTTTGTTTCTGGCAGCCATGGCCCA
GCATCACCTCCTGCATAACAAGATTCGTCGATAATTTACTTAATTCACCGCCAGTGTTATTTTTGATAAT
GTAAGATTGCTGCGGGTTAACGGATCTAACGTACAGCCAAAGTTATGTTCAATGCAGCTGGCAATATAGG
TTGCGTTAATTTTTGGTGGCACATATCAAGCGCCTCAACAGCCTTCGCCGCCGCTTTGTCAACAAGGTGC
GTAGCAGCAGGGGTATCAACGTTTGCATTTCAAGGTGCCGGGCTTCCCGTCCTACGCTGGTACCCTGCTC
GGTCCAGCCGCAGGCTCCGGAATTTTAGGGTCAGGCAGTTTATCTGCAGGAATTTGTGCAAAAACATTAC
GTTACTGCTATCAATATTAATATTGATATTAATACCCGTAGGTTCTGCTTTTTTCCCACCATCAGGTACT
CGATATGGCGCTGGCTATTATCGACATGGCTTCGGCTATTGTTAATATGCTTACTGTTATCCACGCTATG
CGTTGGCTGTTATCAATAAATGTTCGGCTATTATCAATCGTCTTCCGGCTATTATCATGGTTGCTATTAT
AGTCGAATCCACACGCGAATGACTATGACTTACGTTGCTGCTGTTAGTCGAATGGTGTGACTCGCCATTG
CAATTTTACCCGCGATGCCATGGTCAATTGCCCCTGTGCTGGCGCTGTGTGCCGTCTCCGTTTGATGCGT
ACCACCTTTCCATCTTTACTTTCGGATGAAGCGTTGGTCACAGCCTCTGCCGTCGCATGAGCCGTGACGT
GCTCTTGCTGGTTACGCCATCGACTTCATCAAAAGAGGTTGTTTCACCCGTAGTGCCCTTTTCTGACGTG
CTGCCGTCTGGCTTTGCTGTTTATTATCATCAACGCCGTGCACCGTGGCCTGTACCGGCTTGCCGATAAT
TGTAAGGTATCGGTTTTAACATTCTCGACGCCAGCCAGTTGACTGGCAAGCGGCTTCACATTCACAATCT
CTTATCTGTCTGACTATTATCATTATCGGTCGTACCGGCTTTATTACCTGTAATTACCGTCGTGTCTGAT
TATTCGTAACGCTCGGTACTGTTGAGAGGAAACTATTTTGCTTAAACTTCAAACCAGAAAATGGGCCAGT
GTATCTGGCTCCAGGGCTTTACGTATCACACCCAAAAACGTCTCACGGGCGTCGAAATGCATTGAATGCA
TATCGTGTCATTGGGTAATATTTCCGCACGCAAGCCGTCAAACGCTCTGCGTACAGGGAACGGTGTGGAG
TCTCCATTTTCGCTTTCAATTCTGGGTGCTTGTCAATATCGCTACATTGCGCTTTCAACGCTGCACTCTT
CTGGCAAACTGCAGTGCAATCTCAGCCAGTTTTTCTTTTTGTGGATGATGAGTAATGACCTCTTTGAGAG
GATAAGGCCGTCCAGCACATTACTTAACACATAGTCAGTTTCCCCTTTTAGCCTGGTCAGCCCCGCCTCT
TAACCCGATCCACGCCAGGTTTATTCAGGTATGACTCGTAAGCCGGGCCAGCTCGCCAGCTACCGTCTCC
TTAACGCTGCATGTGCAAGCCATCAACGGTAGTAA
>003197.2:c3010816-3009905_reverse_complement
TATTACGAGTCGTAGGACTTAGCATACTTTTCCCTCCACATGATAGCTCCTGCACCGAAAATATCATCTT
ACAATACAGGTCAATGCTTCGGTGTTTAGCTCCCCGCTGACGTCAGGCTGAGGCGATAAACTCGCCCCCA
TTCGAGATCGGTTTCCTCGTAACTGATATGAATCCAGCTGCCATCTATAATTTCACTTTCCAGACGCTCA
CTGACATAACAAATTTAAGATTTAACTCTGGATACTTCTTATTGGCCTGTGCAACAAGAAAAGGCATCAT
TTTACGAGAAATCGGGCAGATTCAACGCCTTCCATTCTGGTCTCCTTTATAGAGGAAACAAGCTCATGGA
AATCACCGATGTCTTTCCGTCCATAACTCTAAAATCGACGACTGAAATATGGATACCTGAACTTCCCATA
CAAGAGCTGCTTTGGTCCTCAACGCCAGTAAAGCAGGTCCAAAAGCGCTACACGCTGCTGGTTCGAACAA
CATTCAGAAGAGCTTCGTTGAATGTCCAGCTCGACCATAGCAAAATAACAATCAGGCAGTTGTTCACGTT
TTCATGAATTTTTACTAAGTTCATAAATTCAAGCTGAAGTTTTTTGGCGAGCGCCAGGCTAAAAATACCG
GCAGATAGCGGTCTGCTTTATCGGACGGGAGAAAAGGTTCTTCACCACACAGACGTTCACAAATATTATC
GCCGCCGGATTGGCCTCCACATATTCATTTAATCGTTGTGCGCCTTGAGTATGTTTGTAAAAACTCACCG
TAGCATAACGGCATTGTTATCGAATCGCTCATAAAGCGTTTCATTCTTTTTGTTTACTATACTGCTTCCC
ACGGTTTAAGTAAAGACTTATATTCAGCTATCCTTTTTTTATGAGCGGATATAGAGAGTTTTTTATCGTT
TT
>003277.2:c27903-27178_reverse_complement
CCCCATCATAAGCAGCTACAATATTCAAAGGAGGGATGTTTAGATTTAGATTAGGCCTATTTATGGGCAT
ATCTTCCGCATTTGGTTGTGCAAGGAGTTGAAATTATTTTTCAGGTGCTTATTTGTAGATGGTATTTCCG
AAAGAAACCACTCTGGCGCATCCCTGAATAGTCAGGCACATCCAGCGCCTCTTTAAAGTGGGATACCGGC
GTGATTCTCGACGAGCATGGATGAAAACATCGTAACCATGGTTATTCAGCTGAAAACCTTGGCTCATAGC
GCTTGTGGCACCATATCCCTGAGCACACTGATGTGGAACTTGTCACCGGCAAATTTGCCCTGAGACTGAG
CATATCGGTCACTTTCCACTTATCTACCGGACTGTCCTCTGAAAACAGCAATCCGGACAGCGCTTGAAAT
CCTGGTCTGGTTTTATATACAACGTGAACTGAGCGCCCAGGCTAACACGGGCTTGTTGAACGACCTTCTC
AGTCTGGATTCCAGACACTCTATAAATTGCCGTGTCTTGTGGAGAAACGACGCACTGTACTGCGAATTTT
GAGATATTTCCAATTTTCAGGATGAACGTCTGACTCAGGACACTGTCCTGAAATAACCCCATTTTCGGAT
GTTCCTCACGTAAAGCCTGTCTCTGCATTTCGCCACCATCACGCCCACTTCGTAGTTCATTACGATAACT
TTACTCTGTCATCAAACGATAAAACG

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.

Welcome to the site @tonu. Please make your question clear within the question. Finally, what coding are you happy with, bash(shell), R, Python, Perl, awk etc ...?

@Tonu this has been correctly answered could you please consider upvoting/accepting?

Since you tag python, I suggest a simple Biopython solution:
from Bio import SeqIO

print(sum(s.seq.count("TAA") for s in SeqIO.parse("file.fa", "fasta")))
# 97 on your input file

If you wanted an overlapping sequence count, Biopython has a method for this too:
print(sum(s.seq.count_overlap("TAA") for s in SeqIO.parse("file.fa", "fasta")))

Here they are the same, but for a string like "TAT" they could be different:
>>> from Bio.Seq import Seq
>>> Seq("TATAT").count("TAT")
1
>>> Seq("TATAT").count_overlap("TAT")
2


text:
Extract compounds/proteins from plant genomic data
I want a process by which I can identify genes within plant genomes; i.e. a process that works for all plants, not just for a single species.
I'm wondering, if, by bioinformatic means only, it is possible to "extract" genes and, after that proteins and natural compounds from plant genome data that is available on NCBI (please see https://www.ncbi.nlm.nih.gov/datasets/genome/GCA_029618835.1/ as an example)?
Example: just by processing the data from the link above we can get to Flavoxanthin or any other substance.
What other useful information from a chemical compound perspective can be found from digging into the genome itself?
Any hints will be greatly appreciated.
Thank you!

I believe I understand the project, to obtain a biochemical/metabolic profile of a large number of poorly annotated plant genomes. There isn't a perfect solution, like Prokka for bacteria. It is clearly an area subject to recent development and interest.
I would recommend using Gene Ontology Meta Annotator for Plants (GOMAP) which is described in "Plant Methods" journal.
There are limitations:

KEGG is not used (below);
its a singularity container - and thats its only distribution.

Learning Singularlity You will need to learn running the container through Singularity and that will be limiting if you are not familiar with containers. However, this is the best of the available "out of the box" solutions I can see in plants. A Singularity HPC is a good call here, tricky if you don't have access to one. I know Docker, Singularity is similar but different.  The information on Singularity is here. I would start reading from "Executing Commands"
singularity pull docker://godlovedc/lolcow # get the container
singularity exec lolcow_latest.sif cowsay moo # run it

(the cow is a famous teaching tool for containers)
KEGG One approach that can be used is to perform COG (NCBI/Genbank code) and KEGG no. identification (e.g. eggNOG). Tools such as eggNOG are setup for metagenomics data however. The COG/KEGG no. output is only the starting point, interrogating each of these to assign gene function and metabolic pathways will require a KEGG database and custom code. Unless someone has published, e.g. KEGG database in SQL.
In my experience the Kyoto Encyclopaedia of Genes and Genomes (KEGG) is the best singular functional database available, thats why I see GOMAPs exclusion of this as a limitation. COG doesn't provide a particularly useful information functional description IMO. eggNOG, or any other method to identify COGs  and KEGG numbers, interrogated against the KEGG database works for plants very well. I know this because I'm not focused on plants but often recover plant  defence mechanisms genes from KEGG. Basically, KEGG identifies >> half of all genes with known COG nomenclature and notably higher than other functional databases.Thus I see GOMAPs exclusion of KEGG as a limitation. However, to reiterate to go identify COG/KEGG -> KEGG database interrogation and analysis is a lot of custom code.
Databases GOMAP implements Arabidopsis annotation, which may be cool and provides a starting point. It's also running UniProt, UniProt-plants, interpro, GO databases and those + Arabidopsis annotated DB could be really cool. It definitely sounds promising.
Processor time Any method of gene identification takes AGES and I have no doubt that GOMAP will be REALLY slow, looking at the pipeline it's using - there are a lot of databases its using. This is GOMAPs justification for using Singularity, it going to take ages to run therefore containerise it in a format often used by modern HPCs. To be honest it's a bit lazy, especially as it was published in 2021 (should be uploaded e.g. to conda). The authors talk about Python and R ... hmmmm, it looks to be a MySQL backend, thats where its speed will come from. The container will be large BTW hopefully is <1GB, but it will be >200MB.
Summary Its worth learning Singularity just keep to the basic commands to run GOMAP. You will need access to a large machine and the output will be okay but they're better general databases. However, I don't know how useful the Arabidopsis database it - thats about the biology, i.e. what you are looking for and whether Arabidopsis is a good model of that investigation
Note, Springer Nature refused a link image request for the container overview, but the image I wanted to share is here.

Just to mention the essential missing part the Singularity container test described singularity exec lolcow_latest.sif cowsay moo. The output is moo. Cows say "moo".

Thank you for pointing to GOMAP, seems a good fit for my plan. Basically, I'm doing a plant database (already have an online demo) as my Master's final project and I wanted to add these detectors besides aggregating other data.
I never heard about Singularity before, but, at base, I'm a software developer and I worked with docker, vagrant, virtualbox and other similar systems so it shouldn't be an issue to get it rolling.

Thats good. The project makes a lot of sense. Singularity is huge in biological science (Docker uptake was low). The only thing to add is entering the Singular container is very different to Docker's -it. From there both the code and databases are open for inspection.


text:
probability of finding a 5 amino acids in a row within a proteome
How to calculate the probability of finding two proteins that share a 5 amino acid long motif from a proteome of around 1067 proteins that have an average length of 65 residues.
The probability of a specific sequence that is completely the same is 1067/20^65

Its just significant at 2%

0.055 . 65 . 1067 = 0.021 (its actually 0.02 - see note)

Where critical probability is 0.05, i.e. 5%
There are 20 aa so 1/20 probability of random occurrence. The oligopeptide is 5 amino acids long, hence to the power 5. Obviously there are 1067 proteins and the probability is per protein, thus if these were independent its a straight multiplication because that is how independent events are measured.
Note The 65 is really 60 because it's a sliding window moving 1 amino acid at a time and the proteins are only 65 in length. The actual value is 60 rather than 65 because at the 3' terminal of the sequence the sliding window can't capture the remaining 4 amino acids, because once the oligopeptide < 5 the calculation is no longer valid. I used 65 so you could see where the value is originating.
Obviously the proteins in question are not independent so you need to state the assumption, i.e. assuming independence.

There are 22 amino acids, not 20. Almost all animals (the only exceptions I know are a few insect species) have 21, because of selenocysteine, the 21st amino acid. I know less about pyrrolysine, the 22nd, but I believe it is found in some bacteria and archaea only. Given how few occurreneces of this amino acid most species have, I doubt it will change the analysis much, but there are 21 amino acids in most non-plant eukaryotes.

Thanks @terdon, now I know.  Selenocysteine is interesting because of the importance of cysteine disulphide bonds. I assume selenocysteine doesn't form them. Anyway, the probability of random occurrence is 1/21 or 1/22 depending.

I... don't remember. Wow. I did my PhD on these things but it's been a while and I always focused on detecting them rather than on their biochemistry, but still! Anyway, Sec has selenium where Cys has sulfur, so presumably it cannot make disulphide bonds by definition, but I do seem to recall that the selenium allows for very similar bonding. We tended to treat Sec as a more reactive version of Cys, but my ignorance of biochemistry is deep so I don't really know.

This question has two aspects.
1 birthday problem
You state that

The probability of a specific sequence that is completely the same is 1067/20^65

but the probability is a bit larger. You can solve this as the birthday problem. If you have $k = 1067$ random draws out of $n=20^{65}$ proteins, then the probability of at least one double case is approximately $$p \approx \frac{0.5k^2}{d} = \frac{569244.5}{20^{65}} $$
2 Markov chain problem
For two random sequences to contain a similar subsequence the probability is approximately

$$\begin{array}{rcl}
P(n \geq 1) &\approx& 1-e^{-E[n]/(1+1\cdot E[S])} \\
E[n] &=& (y-x+1)(z-x+1)(1/k)^x \\
E[S] &=& \frac{1}{k-1}
\end{array}$$
Where $x$ is the length of the sub-sequence, $y$ and $z$ are the lengths of the random sequences among which we search for a matching subsequence of length $x$, $k$ is the number of different types of letters out of which the sequences can be composed (e.g. $k=6$ for six-sided dice rolls, $k=4$ for RNA/DNA-sequences, $k=20$ for your proteins), $n$ is the number of sub-sequence matches, $S$ is the length by which the sequence of length $x$ is extended (finding more matches).

With $k=20$, $x=5$ and $y=z=65$ this gives
$$P(n \geq 1) \approx 0.0011$$
This is described in this question on cross validated Probability of a similar sub-sequence of length X in two sequences of length Y and Z and this question on biostars Probability of finding a common sub-sequence of at least 'k' nucleotides between two sequences
Combining the two
For $m=1067$ proteins there are approximately $n(n-1)/2$ pairs and each of those have a $0.0011$ probability of having a match. This give a total probability of approximately
$$P_{total} = 1-(1-P_{individual})^{1067 \cdot 1066/2} = 1-e^{\log(1-P_{individual})*1067\cdot 1067/2} \approx 1-e^{-628} $$
comparing the formula's with simulation
We can verify the computation with a formula by testing it with a simulation for a smaller number of proteins. Say we have 30 proteins, then the probability of at least one match is 0.3815459.
... to be continued, my phone battery is running out.

Of course, the above is only true for a model where each protein is a complete random draw of a sequence of amino acids. When certain patterns of aminoacids have different probabilities, then matches may be more likely.


text:
Classical alignment with HISAT2
This question was also asked on Biostars
I have this script
#!/bin/bash

# Defining the number of threads to be used
threads=20

# Alignment of sequences to the reference genome using 'HISAT2'
echo -e "\nNow Running Classical Alignment-based with HISAT2...\n"

mkdir -p hisat2

# Download the human reference genome
wget -c ftp://ftp.ensembl.org/pub/release-100/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz -O hisat2/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
gunzip hisat2/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz

# Download the human annotation file
wget -c ftp://ftp.ensembl.org/pub/release-100/gtf/homo_sapiens/Homo_sapiens.GRCh38.100.gtf.gz -O hisat2/Homo_sapiens.GRCh38.100.gtf.gz
gunzip hisat2/Homo_sapiens.GRCh38.100.gtf.gz

# HISAT2 indexing to build the reference genome index
hisat2-build -p $threads hisat2/Homo_sapiens.GRCh38.dna.primary_assembly.fa hisat2/Homo_sapiens.GRCh38v3_hisat2.idx

# HISAT2 Alignment
for sample in $(cat sample_ids.txt)
do
    hisat2 -p $threads -x hisat2/Homo_sapiens.GRCh38v3_hisat2.idx \
        -1 trimmed_data/${sample}_1_val_1.fq -2 trimmed_data/${sample}_2_val_2.fq \
        -S hisat2/${sample}_hisat2.sam

    # Compress the SAM files to binary format, sort them, and index them.
    samtools view -Sb hisat2/${sample}_hisat2.sam | samtools sort > hisat2/${sample}_hisat2_sorted.bam
    samtools index hisat2/${sample}_hisat2_sorted.bam
rm hisat2/${sample}_hisat2.sam  # Remove the .sam files for storage purposes
done

# Quantifying Aligned reads using the Subread package's script 'featureCounts'
for sample in $(cat sample_ids.txt)
do
featureCounts -T $threads -a hisat2/Homo_sapiens.GRCh38.100.gtf \
        -o hisat2/${sample}_hisat2_counts.txt \
    hisat2/${sample}_hisat2_sorted.bam
done

# Generating analysis report with MultiQC
echo "Running MultiQC on HISAT2 alignment results ........."
multiqc -o multiQCreports hisat2
echo "MultiQC on HISAT2 alignment results completed successfully."

# Display message to show that the script is complete
echo "Script completed!"

However, I have receieved an error below
Load annotation file Homo_sapiens.GRCh38.100.gtf ...                       ||
||    Features : 1378275                                                      ||
||    Meta-features : 60683                                                   ||
||    Chromosomes/contigs : 47                                                ||
||                                                                            ||
|| Process BAM file SRR19170514_hisat2_sorted.bam...                          ||
ERROR: Paired-end reads were detected in single-end read library : hisat2/SRR19170514_hisat2_sorted.bam
ERROR: Paired-end reads were detected in single-end read library : hisat2/SRR19170514_hisat2_sorted.bam

Question: What can I do to solve this?

I believe the error is from featureCounts and the solution would be to add the -p flag to the featureCounts call, stating that you are providing paired-end data.
EDIT: You might want to check this r/bioinformatics (reddit) thread out.

Well, thanks. It worked when I added the -p flag to the featureCounts call

@EmmanuelAroma Accept haci's answer to mark the post as answered.

thank you so much comrades, I've got a wonderful answer


text:
Longstitch error make: command: Command not found *** No rule to make target
I installed Longstitch and ran the test script with no issues. The output files matched the expected output files. But when I am now trying to run Longstitch on my own data I am getting this error.
(medaka) operations@LabPC1:~/medaka/longstitch-1.0.5$ ./longstitch run draft=~/medaka/~medaka/medaka_results/consensus.fa reads=~/medaka/combined.fastq G=4.5e8
make: command: Command not found
make: command: Command not found
make: *** No rule to make target'/home/operations/medaka/~medaka/medaka_results/consensus.fa.fa', needed by'/home/operations/medaka/~medaka/medaka_results/consensus.fa.cut250.tigmint.fa'.  Stop.

I changed the file name from .fasta to .fa for my draft assembly because it requires a .fa file, but now I am getting ".fa.fa" What is wrong here?

Please provide a link and/or additional information about what Longstitch is. It would also be useful if you could provide a few lines of input and output (e.g. from both the example data and your own data).

The main issue I can see is the filenames, and should be easily solvable.
Simply put, LongStitch is very particular about filenames. It expects the input draft to have the extension .fa, and the input reads to have any from a range of extensions (your .fastq is perfectly fine, though).
The main issue is that LongStitch expects the extensions to be missing from the command. So if your draft is draft is called consensus.fa and your reads are  combined.fastq, then the command would include something like draft=consensus reads=combined, not draft=consensus.fa reads=combined.fastq.
(unsure if it accepts input from different directories, or whether you'll need to move your input to the same directory as well)
More generally, when having issues with a particular tool which is available via github, I would encourage you to have a look at the 'Issues' or 'Help' section. What I normally do is to delete the 'is:issue is:open' terms in the search bar, which allows me to search for my bug among older, potentially solved issues (e.g. https://github.com/bcgsc/LongStitch/issues/54), and also generally to check that the software is still supported.
Good luck!
(*that said, M_'s suggestion of using mamba is generally a good one, and has helped me with another LongStitch bug)

Could you download it from conda and try again? Its here, https://docs.conda.io/en/main/miniconda.html
Then,
conda create --name myenv
conda activate myenv
conda install -c bioconda longstitch
longstitch ...

Conda avoids all the build stuff.

Notes
If you find that "conda" is slow and doesn't ever load.
There is an alternative to conda called mamba, this is fast and designed to address this exact issue. I have never needed to use it so the following is a rough guide only.
Download and install via
conda create --name myenv
conda activate myenv
conda install -c conda-forge mamba

The documentation is here, https://mamba.readthedocs.io/en/latest/user_guide/mamba.html
I've never used it but the idea is ...
mamba create -n mambaenv <list of packages>
mamba activate mambaenv

The list of packages could be left blank for now, then,
mamba install "Python>=3.10.0" numpy pandas pathlib
mamba install longstitch

Can also do
mamba install -n mambaenv longstitch

NOTE longstitch exists strictly in -c bioconda. You do not seem to search and install via channels in mamba
You search for a package and its dependencies via
mamba repoquery search longstitch
mamba repoquery depends -a longstitch

It is important to place a virtual environment between the system Python and package Python (don't want a system corruption).

Finally ...
For fa.fa ... hmmmm, try ls *.fa and check for .fa.fa files.


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

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

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

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

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

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

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


text:
Mitochondrial genes - TPM calculation bulk RNA-Seq
This question was also asked on Biostars
I was wondering if any of you have encountered a situation for bulk RNA-Seq where, possibly due to low sample quality or presence of dead cells, mitochondrial genes are expressed to a very large degree relative to other genes, thus skewing TPM values of all nuclear genes (by effectively scaling them down).
In such a situation, what could I potentially do to alleviate this beyond preparing fresh samples? Could I, for example, exclude 'MT-' genes from samples and then recalculate TPM based on this filtered set of genes?
My aim was to get some idea of within-samples expression levels, hence why going for TPMs as I know that counts normalised by DESeq2, edgeR and the like do not correct for gene length and are inappropriate for within-samples comparisons.
Many thanks in advance for any insight

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.

Yes, this is a common situation. Not just for mitochondrial genes, but for any gene that may be an outlier (e.g. the B-cell receptor gene in a contaminant/unexpected B-cell population). The way to get around this is to not use TPM for differential expression analysis.
DESeq2 by default uses a quantile-based normalisation method for bulk differential expression calculations, which works well when the majority of genes are not differentially expressed. When this assumption is broken (e.g. for a targeted panel), a control gene list can be supplied, or custom size factors can be applied to each sample.
More information about using DESeq2 for differential expression can be found in the DESeq2 manual:
https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
For comparing expression from different genes, I normalise variance-stabilised data using the longest annotated gene length from transcripts.
See this answer and this paper for more details.

Hi, thanks very much for your response. I should have explained this at the outset, but my aim was to get some idea of within-samples expression levels, hence why going for TPMs as I know that counts normalised by DESeq2, edgeR and the like do not correct for gene length and are inappropriate for within-samples comparisons.

Added a little more detail about how I do it, using the VST from DESeq2.

Thank you, appreciate the help I will take a look at those links


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

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

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

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


text:
How to create a phylogenetic tree from diverse mitochondrial genomes
I would like to create a phylogenetic tree for the most species in my dataset.
I'm starting with around 1200 species, but since it's not good practice to align short and long sequences I tried filtering only for the species with the [15k-18k bp] range in length, but this leaves me with around 350 species in total.
Those sequences vary a lot in size, here's an histogram of the number of species I have data for

As you can see I have a lot of sequences with length [12k-15k bp] and [15k-18k bp] (values on the x axis are scaled by a factor of 3000 and then rounded).
I would like maximise the species genetic information.
To do so I would like to extract the part of genome that is shared the most among all the species.
Is this achievable by using a regex?
I am not sure on the length and type of the query to be sure of not extracting non-homologous regions as it would create a big problem later on for the phylogenetic tree creation.
I can't make an alignment for all species because of the lack of computational power. That is why I'm trying to subset a piece of DNA that is the most shared (prior to the alignment, but I don't know if that is doable or logic)
My question is: how would you proceed to get the most phylogenetic data out of this dataset?
If you would use a regex, what type of query would you use?
More Info
Alignment
When I try to align 10 sequences drawn from the [15k-18k bp] range I get this result (the numbers on the right are the TaxID of the species).

We can't see much here so I used an online MSA viewer
[https://alignmentviewer.org/]

Data
My data is one Fasta file containing mitochondrial genomes for all the species.
To do this preliminary analysis I only used 10 randomly chosen species with a sequence of length in the range [15k-18k bp].

Have you considered using only genes for the alignment instead of the entire sequence?

Is there a way for you to find an existing vertebrate phylogeny (this seems to be your group of interest) based on a larger set of sequences and to prune it to the set you care about? e.g. idk here, see .newick and .tre dowloads. iTOL makes interactive tree exploration/pruning pretty easy.

Your overall approach is cool, given caveats (below) to adapt the method. I can fairly authoritatively say your specific methodological approach will not work regarding its phylogenetic strategy without changing it. This is due to a phenomenon where the number of mutations becomes so great it is impossible to recover the ancestral state, basically it randomises the phylogenetic signal, this especially occurs at the 3rd codon nucleotide positions, i.e. 1 in every 3 nucleotides, 33.3% of the entire data set. This violates a central tenant of phylogenetic theory, therefore I am afraid that you cannot build a tree using this approach, but again you can adapt it to work (below).
If you want only to use nucleotides then 16S is the only target available. However, 16S has been a good target for amphibians, therefore just aligning that gene and building a tree is cool. Its perfectly solid given prior work in vertebrates ...

Chan, K.O., Hertwig, S.T., Neokleous, D.N. et al. Widely used, short 16S rRNA mitochondrial gene fragments yield poor and erratic results in phylogenetic estimation and species delimitation of amphibians. BMC Ecol Evo 22, 37 (2022).

If you wish to use protein coding genes you need to switch amino acids and align those. In terms of conservation there alot of established tools rather than "reinvent the wheel". One program you could consider is AAcon, available here. This would give you can empirical basis for deducing which gene is

the most conserved
the most represented, i.e. "read-depth", or species representation.

... after its aligned, thereby yielding a genetic marker for phylogenetics. It uses a large number of metrics therefore most scenarios would be accommodated.
COX1 is good target and likely gene to emerge from AAcon analysis.
Building protein phylogenetic trees is a separate question.

Lacking computational power to align. That is solvable, even on small computers, and you can place that as a separate question stating your computer capacity (RAM and CPUs). It's not difficult to align 1000s even on a small machine, it's simply the strategy to do this. Algorithm performance and accuracy change rapidly.

It looks like these are random mitochondrial sequences pulled from NCBI, which is probably why the alignments aren't working well.
As a starting point, I recommend you make sure that all the mitochondrial DNA sequences are rotated and reverse-complemented as necessary to make sure that the same gene matches position 1 of the sequence. I prefer doing this for COX1, because it's a well-conserved gene. After doing that, it'll be much easier to look at the alignment of that gene alone, which will give you an idea of the complexity of extending the task to the full length of the mitochondrial genome.
I don't expect you'll find a long exact-match sequence that is consistently conserved across all species, even within COX1. It wouldn't surprise me if the longest conserved sequence were less than 16bp, with a few repeats of that sequence throughout the genome. If that is the case, you need to think about what to do with that information. Given the complexity of this intermediate problem you've asked about, it would be a good idea to think first about what you would do after you discovered there were no such sequence.


text:
time*treatment series with repeated (i.e: NOT independent) sampling of replicates
I am performing the analysis of cell cultures in suspension, untreated(U) and treatment A and B, at t0, t1 and t2, 4 replicates per treatment.
The experiment started with 12 cultures, 4U, 4A and 4B, and cultures were sampled at the three times. I fathom this is different from starting with 36 independent cultures, each contributing to a single library, but I have trouble figuring out how to set my "design" object to use it in edgeR.
Having a samples.df with columns: sample(36 unique identifiers), time(3), treatment(3), treatment_time(9) and replicate(4),  I built my DGEList in edgeR using group as a factor of treatment_time, since I am primarily interested in exploring differences between treatments A and B using U as baseline at any time point, although the two treatments might show similarities in the mechanism of action, so for example At1-Ut1 might be also compared to Bt2-Ut2 and so on.
The replicates are not independent though, hence I tried set my design as model:
matrix(~0 + treatment_time + replicate)

Is this accurate and am I good for any downstream analysis as per the simple edgeR workflow as described in the manual?
Replicates originate from splitting, after treatment, the same cell culture into four separate vessels, but this is not my main concern. What is not independent is the time points of any given treatment for any given replicate. Imagine replicate 2 treatment A is, say, contamined, then every time point for that culture is contamined. Now, contamination is expected to have a huge effect on the library profiles, and I am using this example just to emphasize the concept. What I wonder is, how do I treat the data? Setting contrast matrix is all it takes?
Any suggestions, even deviating from my chosen route, are welcome!

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.

Partial answer: because your replicates involved a physical splitting of the input, I would consider them to be sufficiently distinct to be considered biological replicates. EdgeR should consider distinct samples to be replicates if all modeled fields are identical, so you should not include the replicate field in the model design.


text:
GWAS in python?
A question similar to this has also been asked on Biostars
I understand GWAS is done with R, but are any written in Python?
I'm used to python and java, acclimating to R will take a bit of time and my advisor wants some data analyzed. I plan on using R but given the deadline, I want to finish the analysis in Python.

Please [edit] your question and give us (a lot, a huge lot) more context. Tell us what input data you have, what kind of analysis you want to do (do you have the results of the GWAS study or are you trying to run the study itself?), what species and anything else that could be relevant. Of course you can do it in python, you can do anything in any programming language, so you really need to narrow this down to something very specific that can be answered in a few lines.

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

I would adapt PLINK. Whilst Biopython pop gen could be used it is better to use a standalone package.

You can use Python to wrap around command line tools

It doesn't (really) matter what language a piece of GWAS software is written in. SAIGE is written in a mix of R and C++. BOLT is written in C++. You shouldn't be writing a GWAS analysis from scratch, as it ends up being a very complex piece of software. You will almost certainly be interacting with the command line, when running a GWAS, so my advice is to find the software which is most suitable for your piece of analysis, not one which is written in a programming language you like.


text:
1.What size PCR product will be generated using the primers (bold and underlined) in the sequence below?
I would like to know the steps of this question and how to determine the forward and the reverse primers and also the length of the sequence, without using NCBI. If you could explain and give an answer step by step I would appreciate it. Thank you

Wouldnt it be a question for chemistry or biology stack ?

You are the second person to have raised this @pippo1980 exact point, Biology SE was suggested. If you flag this formally it helps to decide about migration.

I think it should be on SE Biology. Trying to flag with this.

This question tells me that you do not understand how PCR works. You should figure it out. Use Google search.

Just to note, the question has now received 2 closed requests.

Without giving too much of the "what they want you to do yourself" bits of this away, here are a few explanations of the image that may help you answer this question:

As mentioned in the text, the forward and reverse primers are underlined in the image.
The - and + symbols below the forward-direction sequence help determine the positions of the bases above and below. A + appears every 10 bases.
The first base in the sequence is at position 1441 (as indicated by the numbers) to the left of the sequence. The first + to the right of that first base represents position 1450.
Primers are usually represented as a sequence in a "5-prime to 3-prime direction", which is almost always the left to right direction of the top sequence in an image like this, and the right to left direction of the bottom sequence.
The length of a PCR product includes the primer sequence.

For help beyond that explanation... this looks like a canned question. Unless you give more context around what bioinformatics project you're working on and how this question relates to your project (i.e. explain why it's a real-world problem rather than a text-book question), I'm going to leave answers to the specific questions asked as an exercise to the reader.

So in this case what would be the length?

I'm happy to help you understand the question, but I don't want to help you answer this question beyond that understanding. Other people are welcome to do that if they desire.

Thank you for that! But as the questions says, I need to determine the size of the PCR product. Is there a formula for that?

@Teo you should be able to figure out the formula from what gringer has told you. We are not a homework writing service, our job is to help you understand so you can figure it out. So please read through what gringer has written, and see if you can figure it out. Then, if you still need help, update your question and add what you think the solution is and how you found it and then we can help you take it from there. If the information given here isn't enough for you, then update your question and add what you understand and what you feel is missing.

something about 1.5MDa in weight ?

@Teo Indeed there is a formula: b - a, where a is the start and b is the end point :)

@Teo if I start reading a book on page 20, and then close the book after page 30 and go have a cup of coffee, how many pages have I read?


text:
Producing just network files from BiG-SCAPE
I was wondering if anyone knows off-hand if there is a way to run BiG-SCAPE without producing the svg or html files? All I want is the networking files (which I can play with in CytoScape).
There isn't anything obvious from reading the parameters section in the wiki but perhaps I just missed it? Maybe like a minimal mode or something similar?
Normally I would simply delete them, but I'm running around 20 000 BGCs through it. And I thought it would save time and disk space if those files are just never written out. I'm simultaneously trying out bigslice but the documentation is surprisingly tricky to follow.


text:
BLAST+ not working on my Command Line
As you know, ncbi has several databases, including nt, nr, etc.
What I am going to do is to blast a sequence of approximately 30 base pairs to the nr database and extract all similar sequences.
For this purpose, with the help of YouTube, I downloaded and installed Command Line Blast on my Windows PC.
Then, from the database section I downloaded all the nr databases, which were 83 zip files and 200 GB in size.  Then, I unzipped it and stored it to a folder on my computer.
Because I saw that the size of the first two zip files that came out of the zip increased from 25 GB to 50 GB, I only unzipped two files to see if my command line works or not.
So my database became a folder of almost 50 GB.
Then I ran the simplest command in the cmd of my computer.
blastn -query D:\query.fa -out D:\results.txt -db D:\nrDataBase

It's been almost two weeks now that I'm stuck day and night with this simple line of code and it doesn't work.
It says once that the query file is not accessible. Once it says that no suitable files were found in the database, etc.
I tested every solution that came to my mind, I used the text format for the query, but it didn't work
I downloaded one or two files from the nt database, I said maybe there is a problem with the nr database, but it still didn't work.
I made a small fasta file, introduced it as a database, but it still didn't accept.
I tried also on the Macbook, it still didn't work, either the query got an error or the database.
The blast has been well installed, When using -h it shows the full info of the tool. It doesn’t work means that I encounter different errors like: The query is not acceptable Or The alias or index file is found in the database I encounter these errors while e. g., the query file is a normal fasta file. I mean, the errors don’t seem to be logical, it looks like something else is going wrong.
I'm really upset, I told you to send me a message to see if you can help me.
Thank you very much and sorry for the long message.

"it doesn't work" is not enough information. What do you see on the command prompt when you just run blastn -h?

@M__ the question was migrated; the OP does not have a profile on bioinfo SE yet.

Thanks @RamRS of course I should have spotted.

You state that you are not using the entire database, you're using just a subset of the resource files. That could be causing the error. Just because you're using a subset of the resource files does not mean you're picking a functional subset of the database.

Please provide more context around what you want to do. It is not necessary to download the entirety of the nr database in order to extract only the sequences that match a specific subsequence. There's probably a better way to achieve your goal without massive downloads.

try running makeblastdb on your downloaded database files, then re-run your blast command

Are you working in the correct directory? Meaning, are you in the directory of the script you are trying to run? What did the error message look like? Please run again and provide that detail so it gives us more to work with!:)


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:
ZIP the fasterq-dump files in loop
I want to zip the files generated by fasterq-dump simultaneously as they are generated. The code that I tried is this.
cat SRR_Acc_List.txt | while read i
do  fasterq-dump ${i} --split-spot -v --skip-technical  -e 15| gzip -k  ${i}> ${i}.gz
done

Haven't tried, but I believe you may add -Z to fasterq-dump and remove -k ${i} from gzip. In addition, I would recommend the following:
cat SRR_Acc_List.txt | xargs -i echo fasterq-dump {} --split-spot -Z -v --skip-technical -e 15 \| gzip > {}.gz | sh

It is simpler and more flexible. For example, replacing sh with parallel would allow you to download in parallel.

I tried this and it worked
cat SRR_Acc_List.txt | while read i ; do  fasterq-dump ${i} --split-spot -v --skip-technical  -e 20 ; for j in i ; do   gzip ${i}.fastq  ; done; done


text:
How to convert FASTA to BED
I have a FASTA file:
> Sequence_1
GCAATGCAAGGAAGTGATGGCGGAAATAGCGTTAGATGTATGTGTAGCGGTCCC...
> Sequence_2
GCAATGCAAGGAAGTGATGGCGGAAATAGCGTTAGATGTATGTGTAGCGGTCCC....
....

I want to generate a BED file for each sequence like:
Sequence_1 0 1500
Sequence_2 0 1700

The BED regions will simply be the size of the sequences.
Q: I did that before with a one-line command. I don't remember what that is, it was on Biostars. I can't find the post now. What's the simplest way to do the conversion?

Do you mean this biostars post? https://www.biostars.org/p/15476/

@Greg The post doesn't give me one-line command to do what I want. I'm 100% sure I saw it somewhere on biostars a while ago, but I don't remember the name of the program that can do it.

What about this one?  https://www.biostars.org/p/15476/#189089

@Greg It was a Python program. I don't remember the name and thus no idea how to search.

Would we need to create a new tag, maybe something like "format-conversion", "file-conversion"? As I expect there will be many more convert from "X to Y format" type of questions.

We have many excellent answers! This will be an excellent reference for future users.
I found what exactly what I was asking in my question:

https://www.biostars.org/p/191052/

$ pip install pyfaidx
$ faidx --transform bed test.fasta > test.bed

This is the one-line command I was asking. The other answers also work, but I want to accept my own answer.

Developer of pyfaidx here. I was just about to write this answer and lament that I didn't pick a more pronounceable or memorable package name.

You can do this easily with bioawk, which is a version of awk with added features facilitating bioinformatics:
bioawk -c fastx '{print $name"\t0\t"length($seq)}' test.fa

-c fastx tells the program that the data should be parsed as fasta or fastq format. This makes the $name and $seq variables available in the awk commands.

It's good practice to have your FASTA indexed, so you can leverage the .fai you are likely to already have. If not, you can just generate the index with samtools and use some awk to make your BED:
samtools faidx $fasta
awk 'BEGIN {FS="\t"}; {print $1 FS "0" FS $2}' $fasta.fai > $fasta.bed

This will maintain tab separation but you can drop the BEGIN statement to use spaces. The BED spec only requires "whitespace" for the simple BED format.

You could adapt this awk one-liner. Note that it assumes that sequence IDs are not longer than 100 characters and that there is no description following the sequence ID on the header line.
cat myseqs.fasta | awk '$0 ~ ">" {print c; c=0;printf substr($0,2,100) "\t0\t"; } $0 !~ ">" {c+=length($0);} END { print c; }'

Otherwise, any Bio* library (Perl, Python, Ruby) provides FASTA format parsers which will extract sequence IDs and lengths.
I'd point out that whilst this resembles BED it is not, strictly-speaking, since BED maps to coordinates on a chromosome or some longer sequence object.

Here is an approach with BioPython. The with statement ensures both the input and output file handles are closed and a lazy approach is taken so that only a single fasta record is held in memory at a time, rather than reading the whole file into memory, which is a bad idea for large input files. The solution makes no assumptions about the sequence ID lengths or the number of lines that the sequences are spread across:
from Bio import SeqIO

with open('sequences.fasta') as in_f, open('sequences.bed','w') as out_f:
    for record in SeqIO.parse(in_f, 'fasta'):
        out_f.write('{}\t0\t{}\n'.format(record.id, len(record)))

Inspired by this answer to a related question on read length distributions, you could do this with Biopython:
from Bio.SeqIO import parse
with open("regions.bed", "w") as bed:
    for record in parse("regions.fasta", "fasta"):
        print(record.id, 0, len(record.seq), sep="\t", file=bed)

I love Python. I prefer this answer simply because it is in Python and that is what we should be using in 2017. But why does this look so much like Perl? Readability counts. I wanted to propose a cleaner snippet, but it looks like comments do not allow code. Would it be considered good tone to edit your answer?

Point taken! I've given it a go, but you're welcome to improve on it. At present there is nobody on the site with enough reputation to accept an edit. I'm happy for you to suggest an edit through comments, or edit the post and wait until someone hits 350 rep. We can delete the meta-comments after the edit is done.

OK! I edited the answer; seems like a StackExchange-y way of doing things after all.

I too love Python, but I feel like you shouldn't need to write a script to do something as trivial as this - there are lots of command line tools that can accomplish this more quickly.

@SamStudio8 On the surface level, I agree. At the same time, if a person has a Python REPL constantly open (which is a lot of people), this has minimal overhead, and, in my humble opinion, quick algorithmic thinking should be promoted at all times in this profession. Moreover, if this conversion is part of a Snakemake pipeline, then it may be even more natural than calling shell() or writing a shell block.

@KirillG I usually have a REPL open, but I think I'd still favour using awk (or such). Years ago was a different story, but learning how to use the basics of awk has been one of the most productivity improving things I have done since switching to use Linux. But sure, I definitely agree if this is part of a larger script or pipeline, I wouldn't dare call shell. Though I'd probably use pysam over biopython, but that's just personal preference :)

@KirillG I wrote a BioPython solution that is arguably more readable/pythonic and doesn't put the whole input file into memory https://bioinformatics.stackexchange.com/a/106/104

@Chris_Rands' answer is much better than mine, please upvote it!

Thanks Scott, although with @KirillG's edit your answer is more pythonic too! However, only more recent versions of BioPython allow a string as input for SeqIO.parse (rather than a file handle) and I don't think BioPython ever explicitly closes the input file handle in this case?

Hmm. Good point about Biopython version - I have only tested the most recent version on Python 2.7.11 and 3.5.1. I was a little nervous about closing file handles too, but SeqIO.parse returns a generator, not a file, so wouldn't it automatically close the handle when we reach StopIteration? (disclaimer - I have no proof of this!)

If the FASTA sequences are all on a single line, then the following perl one-liner should work:
cat myseqs.fasta | perl -ne 'if(/^>([^ ]+)/){print $1} else {print " 0 ",length,"\n"}'

Explanation:

If the line begins with a '>', then print everything up to the first space (but don't put a line break at the end)
Otherwise, print " 0 ", followed by the length of the line, followed by a line break

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


text:
Somatic tumor only variant calling?
I'm evaluating possibilites for somatic tumor variant calling without paired-normal samples. I'm aware of the consequences without a normal sample.
All the popular variant calling tools such as Strelka, VarScan etc require a normal sample.
What are the some of the possibilites I have? I'm looking for good sensitivity in exchange for higher false positives.

There is a recent paper that attempts to do this:

ISOWN: accurate somatic mutation identification in the absence of
  normal tissue controls.
In this work, we describe the development, implementation, and
  validation of ISOWN, an accurate algorithm for predicting somatic
  mutations in cancer tissues in the absence of matching normal tissues.
  ISOWN is available as Open Source under Apache License 2.0 from
  https://github.com/ikalatskaya/ISOWN

I've been using the LoFreq* caller for exactly this. It is designed to find variants with very low frequency, so is well suited for this type of analysis.

LoFreq* (i.e. LoFreq version 2) is a fast and sensitive variant-caller for inferring SNVs and indels from next-generation sequencing data. It makes full use of base-call qualities and other sources of errors inherent in sequencing (e.g. mapping or base/indel alignment uncertainty), which are usually ignored by other methods or only used for filtering.
LoFreq* is very sensitive; most notably, it is able to predict variants below the average base-call quality (i.e. sequencing error rate). Each variant call is assigned a p-value which allows for rigorous false positive control. Even though it uses no approximations or heuristics, it is very efficient due to several runtime optimizations and also provides a (pseudo-)parallel implementation.

(In my experience it is not fast at all, but that's what they claim).
I have used it many times to call somatic variants using the germline caller (and not its somatic one which requires a tumor/normal pair). Because the tool is so sensitive, its germline caller actually deals quite well with somatic samples.
You can find instructions for its use here.

I have used LoFreq as well and think it does a good job.  It is very sensitive and almost always requires some post filtering to get a reasonable number of calls.  One hazard I have found with its model is that if the variant is actually the majority call at the position it can have trouble reporting it.  This can happen with samples that have a large tumor cellularity.

@Bioathlete oh? Thanks, that's good to know!

All the popular variant calling tools such as Strelka, VarScan etc
  require a normal sample.

Strelka and VarScan require a normal sample in somatic mode, but they both have a germline mode for unpaired analysis.
One of the most popular somatic variant callers is MuTect and that does not require a normal.

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